Learning-Induced Dynamical Transition in Recurrent Neural Networks

# Learning-Induced Dynamical Transition in Recurrent Neural Networks
Source: [https://arxiv.org/html/2609.19288](https://arxiv.org/html/2609.19288)
Varun VaidyaAffiliation:Department of Physics, University of South Dakota, Vermillion 57069, USAEmail:[Varun\.Vaidya@usd\.edu](mailto:[email protected])

###### Abstract

Learning in recurrent neural networks can fundamentally reshape their underlying dynamics, transforming initially chaotic activity into stable task\-dependent behavior\. We develop a non\-equilibrium dynamical mean\-field theory\(DMFT\) to describe this transition during learning\. We show that a slow feedback\-driven learning process generates an evolving effective feedback strength that drives the network through a transition from chaotic to stable dynamics defined by a bifurcation of the DMFT solution\. By deriving the two\-time correlation function throughout learning, we identify a critical feedback strength and a corresponding learning rate dependent critical time separating these regimes\. The transition arises from the progressive deformation of an effective dynamical landscape by the growing learned feedback structure\. Starting from the untrained state, the theory predicts the time evolution of the network output during training and shows quantitative agreement with numerical simulations\.

## 1Introduction

Recurrent neural networks provide a powerful framework for studying how complex dynamical systems can generate flexible and stable computational behavior\. Beyond their role as models of artificial intelligence, recurrent networks have also become important theoretical models for understanding the principles governing high\-dimensional biological neural circuits\. A central question in this context is how learning modifies the intrinsic dynamics of recurrent networks, transforming initially irregular activity into structured dynamics capable of performing specific computational tasks\.

Random recurrent networks exhibit a rich range of dynamical behaviors, including chaotic activity arising from the interaction between recurrent connectivity and nonlinear neuronal responses\. Seminal work using dynamical mean\-field theory\[[1](https://arxiv.org/html/2609.19288#bib.bib1)\]established a framework for understanding these collective dynamical regimes and characterized transitions between chaotic and ordered states as a function of network parameters\. Subsequent studies have shown that learning can substantially reshape these dynamics, either through changes in recurrent connectivity or through learned feedback pathways, leading to stable trajectories, attractor states, and task\-dependent dynamical regimes\[[2](https://arxiv.org/html/2609.19288#bib.bib2),[3](https://arxiv.org/html/2609.19288#bib.bib3),[4](https://arxiv.org/html/2609.19288#bib.bib4),[5](https://arxiv.org/html/2609.19288#bib.bib5),[6](https://arxiv.org/html/2609.19288#bib.bib6)\]\.

Recent theoretical studies have investigated how learned structure controls the dynamical state of trained recurrent networks\. In particular, Clark et al\.\[[7](https://arxiv.org/html/2609.19288#bib.bib7)\]showed that sufficiently strong task\-dependent recurrent restructuring can drive a transition between chaotic and stable regimes\. Such approaches provide a phase\-space description of learned networks, identifying the conditions under which stable computation becomes possible\.

However, an important question remains: how does learning itself dynamically drive a network from one dynamical regime to another? During training, the network is not a stationary system\. Its synaptic parameters evolve continuously, modifying the effective recurrent dynamics and, consequently, the statistical properties of neural activity\. Understanding this process therefore requires following the full non\-equilibrium learning trajectory, rather than characterizing only the final trained state\. In this work, we develop a dynamical mean\-field theory for recurrent neural networks undergoing slow, feedback\-driven learning\. Rather than introducing an externally controlled parameter to tune between dynamical regimes, we show that learning itself generates an evolving effective control parameter through the gradual formation of structured feedback\. This evolving feedback dynamically deforms the effective landscape governing network fluctuations, driving the system through a transition from chaotic to stable dynamics\.

A central object in our analysis is the full two\-time correlation function between neuronal outputs, which captures the non\-stationary nature of the learning process\. Building on the two\-time correlation and response formalism of dynamical mean\-field theory\[[1](https://arxiv.org/html/2609.19288#bib.bib1)\], we use this description to track the emergence of long\-time correlations and identify a critical feedback and a corresponding learning rate dependent critical time at which the dynamical regime changes\. We show that this transition corresponds to a qualitative change in the effective dynamical landscape, where the evolution of learned feedback drives the system across a critical point\. Quantitatively this corresponds to a bifurcation of the DMFT solution\. We show that during learning, the two time correlation function exhibits a separation between rapid chaotic relaxation and slow evolution of the correlation plateau, allowing the learning trajectory to be described as a slow deformation of the underlying dynamical state\.

Beyond characterizing the transition itself, our theory predicts the evolution of the network output throughout learning\. The predictions obtained from the dynamical mean\-field equations show quantitative agreement with numerical simulations, demonstrating that the theory captures both the transient learning dynamics and the emergence of the final stable computational state\.

Our results provide a framework for understanding learning as a dynamical process that can reorganize the dynamical regime of a recurrent system\. Rather than viewing stable computation as a property imposed by a fixed learned structure, we show how stability can emerge dynamically through slow plastic adaptation, providing a mechanism by which recurrent networks transition from internally generated chaotic activity to reliable task\-dependent dynamics\.

The paper is organized as follows\. In Section[2](https://arxiv.org/html/2609.19288#S2), we introduce the slow feedback\-driven learning model for our recurrent neural network\. We develop the DMFT analysis for this system in Sections[3](https://arxiv.org/html/2609.19288#S3)and[4](https://arxiv.org/html/2609.19288#S4)\. The time evolution of the system during training is derived in Section[5](https://arxiv.org/html/2609.19288#S5)followed by a comparison with simulation in Section[6](https://arxiv.org/html/2609.19288#S6)\. The conclusions along with an outlook are presented in Section[7](https://arxiv.org/html/2609.19288#S7)\.

## 2The learning model for RNN

We consider a recurrent neural network consisting of \(N\) interacting units whose dynamics are governed by

τc​x˙i​\(t\)=−xi​\(t\)\+g2​∑j=1NJi​j​ϕ​\(xj​\(t\)\)\+∑μWif​b​y​\(t\),\\tau\_\{c\}\\dot\{x\}\_\{i\}\(t\)=\-x\_\{i\}\(t\)\+g^\{2\}\\sum\_\{j=1\}^\{N\}J\_\{ij\}\\phi\(x\_\{j\}\(t\)\)\+\\sum\_\{\\mu\}W^\{fb\}\_\{i\}y\(t\),\(1\)
wherexi​\(t\)x\_\{i\}\(t\)denotes the activity of neuron \(i\),ϕ⁡\(x\)\\phi\(x\)is a nonlinear activation function, andJi​jJ\_\{ij\}represents the recurrent connectivity\.τc\\tau\_\{c\}is the characteristic time scale for the reservoir dynamics\. Throughout this paper, we use the standard activation functionϕ⁡\(x\)=tanh⁡x\\phi\(x\)=\\tanh\{x\}\. The recurrent couplings are drawn from a Gaussian distribution with zero mean and a variance1/N1/\\sqrt\{N\}, and we analyze the statistics of this ensemble\.

The network produces task\-relevant outputs through a low\-dimensional readout,

y⁡\(t\)=∑iwio​u​t​\(t\)​ϕ​\(xi​\(t\)\)\.y\(t\)=\\sum\_\{i\}w^\{out\}\_\{i\}\(t\)\\phi\(x\_\{i\}\(t\)\)\.\(2\)
These outputs are fed back into the recurrent network through feedback projectionsWif​bW^\{fb\}\_\{i\}which are also drawn from a Gaussian distribution with mean zero and varianceσf​b/N\\sigma\_\{fb\}/\\sqrt\{N\}\. Feedback from the readout to the recurrent reservoir is a standard architecture in closed\-loop reservoir computing and FORCE\-based recurrent networks\[[8](https://arxiv.org/html/2609.19288#bib.bib8),[4](https://arxiv.org/html/2609.19288#bib.bib4),[9](https://arxiv.org/html/2609.19288#bib.bib9)\], where it enables recurrent networks to autonomously generate learned trajectories and complex temporal dynamics\. Here we adopt this architecture as a biologically motivated mechanism through which task\-related information can progressively reshape recurrent dynamics during learning\. From a biological perspective, the feedback can be interpreted more generally as a top\-down signal that modulates ongoing neural activity and conveys task\-relevant information\[[10](https://arxiv.org/html/2609.19288#bib.bib10),[11](https://arxiv.org/html/2609.19288#bib.bib11),[12](https://arxiv.org/html/2609.19288#bib.bib12)\]\.

The central feature of the model is that the readout weights evolve slowly according to a learning rule\. We consider a feedback\-driven learning process in which the readout is modified according to the difference between the desired target signal \(y∗​\(t\)y^\{\*\}\(t\)\) and the generated output,

e⁡\(t\)=y∗​\(t\)−y⁡\(t\),e\(t\)=y^\{\*\}\(t\)\-y\(t\),\(3\)
with synaptic adaptation governed by

w˙io​u​t​\(t\)=1N​α​e​\(t\)​ϕ​\(xi​\(t\)\)\.\\dot\{w\}^\{out\}\_\{i\}\(t\)=\\frac\{1\}\{N\}\\alpha\\ e\(t\)\\phi\(x\_\{i\}\(t\)\)\.\(4\)Here \(α≪1/τc\\alpha\\ll 1/\\tau\_\{c\}\) sets the slow learning timescale relative to the intrinsic neural dynamics\. This separation of timescales allows the recurrent activity to relax rapidly for a given value of the slowly evolving synaptic parameters, while learning gradually reshapes the effective network dynamics\.

The purpose of this work is not only to determine the final learned output, but to understand how the dynamical state of the network evolves during learning\. Initially, the random recurrent network exhibits chaotic fluctuations\. As the feedback pathway becomes increasingly structured through learning, it modifies the effective dynamics of the recurrent system\. We therefore study the non\-equilibrium learning trajectory generated by the coupled evolution of neural activity and synaptic parameters\.

The main object characterizing this evolution is the two\-time correlation function,

C⁡\(t,s\)=1N​∑iϕ⁡\(xi​\(t\)\)​ϕ​\(xi​\(s\)\),C\(t,s\)=\\frac\{1\}\{N\}\\sum\_\{i\}\\phi\(x\_\{i\}\(t\)\)\\phi\(x\_\{i\}\(s\)\),\(5\)which captures the non\-stationary dynamics of the learning process\. Unlike equilibrium or fully trained analyses, where correlations depend only on the time difference \(t\-s\), learning generates a system whose statistical properties evolve with time\. This two\-time structure allows us to follow the emergence of long\-time correlations starting from a chaotic state and the transition between dynamical regimes during training\.

## 3Disorder averaging and DMFT analysis

The MSRJD path integral framework\[[13](https://arxiv.org/html/2609.19288#bib.bib13),[14](https://arxiv.org/html/2609.19288#bib.bib14),[15](https://arxiv.org/html/2609.19288#bib.bib15)\]applied to random neural networks\[[16](https://arxiv.org/html/2609.19288#bib.bib16)\]introduces response fieldsx~i​\(t\)\\tilde\{x\}\_\{i\}\(t\)to encode the dynamics through a generating functional written as a path integral:

Z=∫∏i𝒟​xi​𝒟​x~i​exp⁡\{i​∑i∫d​t​x~i​\(t\)​\[τ​x˙i​\(t\)\+xi​\(t\)−∑j\(g2​Ji​j\+Wif​b​wjo​u​t\)​ϕ​\(xj​\(t\)\)\]\}\.Z=\\int\\prod\_\{i\}\\mathcal\{D\}x\_\{i\}\\mathcal\{D\}\\tilde\{x\}\_\{i\}\\;\\exp\\Bigg\\\{i\\sum\_\{i\}\\int dt\\,\\tilde\{x\}\_\{i\}\(t\)\\Big\[\\tau\\dot\{x\}\_\{i\}\(t\)\+x\_\{i\}\(t\)\-\\sum\_\{j\}\\left\(g^\{2\}J\_\{ij\}\+W\_\{i\}^\{fb\}w\_\{j\}^\{out\}\\right\)\\phi\(x\_\{j\}\(t\)\)\\Big\]\\Bigg\\\}\.\(6\)The response fieldsx~i​\(t\)\\tilde\{x\}\_\{i\}\(t\)enforce the microscopic equations of motion within the generating functional\. We can perform an ensemble average over the network parametersJi​j,Wif​bJ\_\{ij\},W\_\{i\}^\{fb\}, using the distributions

P⁡\(Ji​j\)=12​π​N​exp​\{−N2​g2​Ji​j2\},P⁡\(Wif​b\)=σf​b2​π​N​exp​\{−N2​σf​b2​\(Wif​b\)2\}\\displaystyle P\(J\_\{ij\}\)=\\frac\{1\}\{\\sqrt\{2\\pi N\}\}\\exp\{\\\{\-\\frac\{N\}\{2g^\{2\}\}J\_\{ij\}^\{2\}\\\}\},\\ \\ \\ \\ P\(W\_\{i\}^\{fb\}\)=\\frac\{\\sigma\_\{fb\}\}\{\\sqrt\{2\\pi N\}\}\\exp\{\\Big\\\{\-\\frac\{N\}\{2\\sigma\_\{fb\}^\{2\}\}\(W^\{fb\}\_\{i\}\)^\{2\}\\Big\\\}\}\(7\)Following details presented in Appendix[A](https://arxiv.org/html/2609.19288#A1), this yields an effective one neuron equation

τc​x˙i​\(t\)\+xi​\(t\)=ηi​\(t\)\\displaystyle\\tau\_\{c\}\\dot\{x\}\_\{i\}\(t\)\+x\_\{i\}\(t\)=\\eta\_\{i\}\(t\)\(8\)
whereηi​\(t\)\\eta\_\{i\}\(t\)is the effective noise term that obeys Gaussian statistics encoded in the effective action

Seff=−12∫dtds∑iηi\(t\)C¯−1\(t,s\)ηi\(s\)\\displaystyle S\_\{\\text\{eff\}\}=\-\\frac\{1\}\{2\}\\int dtds\\sum\_\{i\}\\eta\_\{i\}\(t\)\\bar\{C\}^\{\-1\}\(t,s\)\\eta\_\{i\}\(s\)\(9\)where

C¯​\(t,s\)\\displaystyle\\bar\{C\}\(t,s\)=\\displaystyle=g2N​∑i⟨ϕ⁡\(xi​\(t\)\)​ϕ​\(xi​\(s\)\)⟩\+σf​b2N​⟨y⁡\(t\)​y​\(s\)⟩\\displaystyle\\frac\{g^\{2\}\}\{N\}\\sum\_\{i\}\\langle\\phi\(x\_\{i\}\(t\)\)\\phi\(x\_\{i\}\(s\)\)\\rangle\+\\frac\{\\sigma\_\{fb\}^\{2\}\}\{N\}\\langle y\(t\)y\(s\)\\rangle\(10\)≡\\displaystyle\\equivg2​C​\(t,s\)\+σf​b2N​⟨y⁡\(t\)​y​\(s\)⟩\\displaystyle g^\{2\}C\(t,s\)\+\\frac\{\\sigma\_\{fb\}^\{2\}\}\{N\}\\langle y\(t\)y\(s\)\\rangleis a correlator which is determined self consistently using this action\.We can further simplify this equation in the large N limit by noting that

⟨y⁡\(t\)​y​\(s\)⟩=∑i,j⟨wio​u​t​\(t\)​wjo​u​t​\(s\)​ϕ​\(xi​\(t\)\)​ϕ​\(xj​\(s\)\)⟩\\displaystyle\\langle y\(t\)y\(s\)\\rangle=\\sum\_\{i,j\}\\langle w\_\{i\}^\{out\}\(t\)w\_\{j\}^\{out\}\(s\)\\phi\(x\_\{i\}\(t\)\)\\phi\(x\_\{j\}\(s\)\)\\rangle\(11\)wio​u​tw\_\{i\}^\{out\}itself evolves through Eq\.[4](https://arxiv.org/html/2609.19288#S2.E4), which is proportional to the activity of theit​hi^\{th\}neuron\. Since cross correlations between distinct neurons are suppressed in the large N limit, the dominant contribution for evolution is⟨y⁡\(t\)⟩​⟨y⁡\(s\)⟩\\langle y\(t\)\\rangle\\langle y\(s\)\\rangle\. The object of interest is then⟨y⁡\(t\)⟩\\langle y\(t\)\\rangle\. Given Eqns\.[3](https://arxiv.org/html/2609.19288#S2.E3),[4](https://arxiv.org/html/2609.19288#S2.E4), we can write

wio​u​t​\(t\)=αN​∫0td​s​\(y∗​\(s\)−y⁡\(s\)\)​ϕ​\(xi​\(s\)\)\\displaystyle w\_\{i\}^\{out\}\(t\)=\\frac\{\\alpha\}\{N\}\\int\_\{0\}^\{t\}ds\(y^\{\*\}\(s\)\-y\(s\)\)\\phi\(x\_\{i\}\(s\)\)\(12\)so that

⟨y⁡\(t\)⟩=αN​∫0td​s​∑i⟨\(y∗​\(s\)−y⁡\(s\)\)​ϕ​\(xi​\(s\)\)​ϕ​\(xi​\(t\)\)⟩\\displaystyle\\langle y\(t\)\\rangle=\\frac\{\\alpha\}\{N\}\\int\_\{0\}^\{t\}ds\\sum\_\{i\}\\langle\(y^\{\*\}\(s\)\-y\(s\)\)\\phi\(x\_\{i\}\(s\)\)\\phi\(x\_\{i\}\(t\)\)\\rangle\(13\)which again, in the large N limit reduces to

⟨y⁡\(t\)⟩=α​∫0td​s​\(y∗​\(s\)−⟨y⁡\(s\)⟩\)​C​\(t,s\)\\displaystyle\\langle y\(t\)\\rangle=\\alpha\\int\_\{0\}^\{t\}ds\(y^\{\*\}\(s\)\-\\langle y\(s\)\\rangle\)C\(t,s\)\(14\)Clearly, to obtain a closed system of equations, we need an evolution equation forC⁡\(t,s\)C\(t,s\)which can be obtained given its definition and the fact thatηi​\(t\)\\eta\_\{i\}\(t\)evolves with the action Eq\.[9](https://arxiv.org/html/2609.19288#S3.E9)\. We setτc=1\\tau\_\{c\}=1henceforth so all time scales are measured in units ofτc\\tau\_\{c\}\. Formally Eq\.[8](https://arxiv.org/html/2609.19288#S3.E8)can be solved to write

xi​\(t\)=∫td​t′​exp⁡\{t′−t\}​ηi​\(t\)\\displaystyle x\_\{i\}\(t\)=\\int^\{t\}dt^\{\\prime\}\\exp\{\\\{t^\{\\prime\}\-t\\\}\}\\eta\_\{i\}\(t\)\(15\)Using this result, we can define the autocorrelation functionΔ⁡\(t,s\)\\Delta\(t,s\)of the neuronal activity

Δ⁡\(t,s\)≡⟨xi​\(t\)​xi​\(s\)⟩=∫td​t′​exp⁡\{t′−t\}​∫sd​s′​exp⁡\{s′−s\}​⟨ηi​\(t\)​ηi​\(s\)⟩\\displaystyle\\Delta\(t,s\)\\equiv\\langle x\_\{i\}\(t\)x\_\{i\}\(s\)\\rangle=\\int^\{t\}dt^\{\\prime\}\\exp\{\\\{t^\{\\prime\}\-t\\\}\}\\int^\{s\}ds^\{\\prime\}\\exp\{\\\{s^\{\\prime\}\-s\\\}\}\\langle\\eta\_\{i\}\(t\)\\eta\_\{i\}\(s\)\\rangle\(16\)which given Eq\.[9](https://arxiv.org/html/2609.19288#S3.E9)obeys the two\-time evolution equation

\(1\+∂s\)\(1\+∂t\)Δ\(t,s\)=g2C\(t,s\)\+σfb2N⟨y\(t\)⟩⟨y\(s\)⟩\\displaystyle\(1\+\\partial\_\{s\}\)\(1\+\\partial\_\{t\}\)\\Delta\(t,s\)=g^\{2\}C\(t,s\)\+\\frac\{\\sigma\_\{\\text\{fb\}\}^\{2\}\}\{N\}\\langle y\(t\)\\rangle\\langle y\(s\)\\rangle\(17\)Sincexi​\(t\)x\_\{i\}\(t\)evolves with Gaussian statistics,C⁡\(t,s\)C\(t,s\)can be expressed as a function of the autocorrelation functionΔ⁡\(t,s\)\\Delta\(t,s\)\[[1](https://arxiv.org/html/2609.19288#bib.bib1)\]

C⁡\(t,s\)=∫D​z​D​y​D​x​ϕ​\(Δ⁡\(t,t\)−Δ⁡\(t,s\)​x\+z​Δ⁡\(t,s\)\)​ϕ​\(Δ⁡\(s,s\)−Δ⁡\(t,s\)​y\+z​Δ⁡\(t,s\)\)\\displaystyle C\(t,s\)=\\int DzDyDx\\phi\(\\sqrt\{\\Delta\(t,t\)\-\\Delta\(t,s\)\}x\+z\\sqrt\{\\Delta\(t,s\)\}\)\\phi\(\\sqrt\{\\Delta\(s,s\)\-\\Delta\(t,s\)\}y\+z\\sqrt\{\\Delta\(t,s\)\}\)\(18\)where∫Dz≡∫dz/2​πexp\{−z2/2\}\\int Dz\\equiv\\int dz/\\sqrt\{2\\pi\}\\exp\{\\\{\-z^\{2\}/2\\\}\}is the integral over a normal distribution\. Eqns\.[14](https://arxiv.org/html/2609.19288#S3.E14),[17](https://arxiv.org/html/2609.19288#S3.E17)and[18](https://arxiv.org/html/2609.19288#S3.E18)form a closed system of equations that can, in principle be solved numerically to obtain the the complete evolution of the system in time\. However, we see that they are highly non linear and due to the triple integral in the definition forC⁡\(t,s\)C\(t,s\), a brute force solution would be expensive\. Likewise, a direct numerical approach also would offer little intuition about the underlying physics\. Therefore, in the next section we will explore the underlying physics dictated by these equations by exploiting the separation in time scales, which will give us some insight at the cost of quantitative precision\. To gain insight and simplify the analysis, we will the set the target to a constant valuey∗​\(t\)=A∼O⁡\(N\)y^\{\*\}\(t\)=A\\sim O\(\\sqrt\{N\}\)\.The scaling is chosen so that the feedback term, when the network is fully trained, is of the same order as the internal chaotic dynamics and is always relevant\. We can consider two other regimes:1\.A≫NA\\gg\\sqrt\{N\}in which case we can expect the system to quickly align with the output since it will overwhelm any internal fluctuations\. The case whereA≪NA\\ll\\sqrt\{N\}is also interesting and we will comment on this in Section[7](https://arxiv.org/html/2609.19288#S7)\.

## 4Learning as a dynamical deformation of an effective potential

One of the central quantities in the set of equations derived in the previous section is the equal time auto\-correlation functionΔ⁡\(t,t\)\\Delta\(t,t\)\. Given the separation of time scales1/α≫τc1/\\alpha\\gg\\tau\_\{c\}, if we restrictt−s≪1/αt\-s\\ll 1/\\alpha, this implies looking atΔ⁡\(t,s\)\\Delta\(t,s\)matrix near the diagonal\. For this part of the matrix, we can adopt an SCS\-type analysis defining a central timeT=\(t\+s\)/2T=\(t\+s\)/2andτ=t−s\\tau=t\-s\. Since evolution over central time is slow, we can ignore any derivatives with respect toTTwhich scale asα≪1\\alpha\\ll 1\. Likewise the output evolves over a time scale1/α≫11/\\alpha\\gg 1, so that it is effectively constant acrossτ\\tau\. With these approximations, Eq\.[17](https://arxiv.org/html/2609.19288#S3.E17)simplifies to

\(1−∂τ2\)Δ\(T,τ\)=g2∫Dz\[Dxϕ\(Δ⁡\(T\)−Δ⁡\(T,τ\)x\+zΔ⁡\(T,τ\)\)\]2\+1N⟨y\(T\)⟩2\\displaystyle\(1\-\\partial^\{2\}\_\{\\tau\}\)\\Delta\(T,\\tau\)=g^\{2\}\\int Dz\\Big\[Dx\\phi\(\\sqrt\{\\Delta\(T\)\-\\Delta\(T,\\tau\)\}x\+z\\sqrt\{\\Delta\(T,\\tau\)\}\)\\Big\]^\{2\}\+\\frac\{1\}\{N\}\\langle y\(T\)\\rangle^\{2\}\(19\)We can define a normalized outputy^​\(T\)=1/N​⟨y​\(T\)⟩\\hat\{y\}\(T\)=1/\\sqrt\{N\}\\langle y\(T\)\\rangle, choosing for simplicityσf​b=1\\sigma\_\{fb\}=1\. Since we are assuming a constant targety∗​\(t\)=Ay^\{\*\}\(t\)=A, the output of the networky^​\(T\)\\hat\{y\}\(T\), whenever learning is successful, is a monotonic function of time\. Hence, we can considerΔ⁡\(T,τ\)\\Delta\(T,\\tau\)as a functionΔ⁡\(y^,τ\)\\Delta\(\\hat\{y\},\\tau\),

\(1−∂τ2\)Δ\(y^,τ\)=g2∫Dz\[Dxϕ\(Δ⁡\(y^,0\)−Δ⁡\(y^,τ\)x\+zΔ⁡\(y^,τ\)\)\]2\+y^2\\displaystyle\(1\-\\partial^\{2\}\_\{\\tau\}\)\\Delta\(\\hat\{y\},\\tau\)=g^\{2\}\\int Dz\\Big\[Dx\\phi\(\\sqrt\{\\Delta\(\\hat\{y\},0\)\-\\Delta\(\\hat\{y\},\\tau\)\}x\+z\\sqrt\{\\Delta\(\\hat\{y\},\\tau\)\}\)\\Big\]^\{2\}\+\\hat\{y\}^\{2\}\(20\)Therefore, at least near the diagonal, the problem reduces to solving the SCS equation with a deformation\. Following\[[1](https://arxiv.org/html/2609.19288#bib.bib1)\], we can adopt the picture of a particle moving in an effective potential which slowly evolves with time but can be considered as stationary over the time scale of the reservoir\. In this case,y^\\hat\{y\}evolves monotonically from 0 toA/NA/\\sqrt\{N\}\. We note again that this picture is only valid whent−s≪1/αt\-s\\ll 1/\\alphawhich is a narrow band around the diagonal of the matrixΔ⁡\(t,s\)\\Delta\(t,s\)\. We call this solutionΔfast​\(y^,τ\)\\Delta\_\{\\text\{fast\}\}\(\\hat\{y\},\\tau\)\. Definingu=Δfast​\(y^,τ\)u=\\Delta\_\{\\text\{fast\}\}\(\\hat\{y\},\\tau\), the equation reduces to

u¨=−d​V​\(u,y^\)d​uwhere​V​\(u,y^\)=−12​u2\+g2​∫0ud​v​∫D​z​\[D​x​ϕ​\(Δ⁡\(y^,0\)−u​x\+z​u\)\]2\+y^2​u\\displaystyle\\ddot\{u\}=\-\\frac\{dV\(u,\\hat\{y\}\)\}\{du\}\\ \\ \\text\{where\}\\ V\(u,\\hat\{y\}\)=\-\\frac\{1\}\{2\}u^\{2\}\+g^\{2\}\\int\_\{0\}^\{u\}dv\\int Dz\\Big\[Dx\\phi\(\\sqrt\{\\Delta\(\\hat\{y\},0\)\-u\}x\+z\\sqrt\{u\}\)\\Big\]^\{2\}\+\\hat\{y\}^\{2\}u
We know that for the original SCS equation \(y^=0\\hat\{y\}=0\), the stable solution is a monotonic decay from the g\-dependent initial valueΔ⁡\(0,0\)\\Delta\(0,0\)to zero forτ≫τc\\tau\\gg\\tau\_\{c\}\. Within the mechanical analogy of the dynamical mean\-field equation, this corresponds to the motion of a particle through the effective potential, evolving fromu=Δ⁡\(0,0\)u=\\Delta\(0,0\)to the stationary point at u=0, as illustrated in Fig\.[1](https://arxiv.org/html/2609.19288#S4.F1)\(a\)\.

Asy^\\hat\{y\}increases with time, the potential becomes progressively shallower while the width of the well decreases\. At the same time, the long\-time stationary pointΔ⁡\(y^,∞\)\\Delta\(\\hat\{y\},\\infty\)moves continuously to a nonzero u, reflecting the emergence of a non\-zero plateau in the correlation function\. Most importantly, at a critical value ofy^\\hat\{y\}, the local minimum and adjacent maximum merge, causing the potential well to collapse and the system moves through the stability landscape\. Beyond this critical timetcrt\_\{\\text\{cr\}\}, the fast\-time fluctuations can no longer be sustained, and the system enters the stable regime\. This sequence is illustrated in Fig\.[1](https://arxiv.org/html/2609.19288#S4.F1)for a representative value of g=1\.3\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/potpanel.png)Figure 1:The deformation of the effective one dimensional potential as a function of the feedback parametery^\\hat\{y\}for g=1\.3\. Starting from the SCS picture in \(a\), increasing feedback with time deforms the potential to a narrower and shallower well with the simultaneous emergence of a non\- zero plateau\. Aty^=0\.2\\hat\{y\}=0\.2, subplot \(f\), we see the collapse of the well signaling a qualitative change\. Beyond this, the solution is frozen at the maxima of the potential and fast dynamics disappear\.Similar transitions between the chaotic and stable regimes have recently been obtained by varying an externally imposed feedback parameterγ\\gammaand studying the asymptotic state of the network\[[7](https://arxiv.org/html/2609.19288#bib.bib7)\]\. In contrast, here no external control parameter is introduced\. Instead, the learning dynamics itself continuously reshapes the effective potential, driving the network across the bifurcation at a finite critical feedback and hence a critical learning time\.

This interpretation immediately predicts that the disappearance of the potential well should coincide with the collapse of the fast\-time component of the two\-time correlation function, a prediction that is confirmed by the numerical solution presented below\.

The condition that the minima and maxima merge at the critical feedback means thatΔ⁡\(y^,0\)=Δ⁡\(y^,∞\)=uc\\Delta\(\\hat\{y\},0\)=\\Delta\(\\hat\{y\},\\infty\)=u\_\{c\}obeys

d2d​u2​V​\(u,y^\)=0⟹1=g2​dd​u​∫D​z​\[∫D​x​ϕ​\(uc−u​x\+z​u\)\]2\|u=uc\\displaystyle\\frac\{d^\{2\}\}\{du^\{2\}\}V\(u,\\hat\{y\}\)=0\\implies 1=g^\{2\}\\frac\{d\}\{du\}\\int Dz\\Big\[\\int Dx\\phi\(\\sqrt\{u\_\{c\}\-u\}x\+z\\sqrt\{u\}\)\\Big\]^\{2\}\\Big\|\_\{u=u\_\{c\}\}\(22\)
Expanding aboutu=ucu=u\_\{c\}, and simplifying we arrive at

1=g2​∫D​z​ϕ′2​\(uc​z\)\\displaystyle 1=g^\{2\}\\int Dz\\phi^\{\\prime 2\}\(\\sqrt\{u\_\{c\}\}z\)\(23\)Interestingly, this is identical to the marginal stability condition obtained in the SCS analysis\[[1](https://arxiv.org/html/2609.19288#bib.bib1)\]of random recurrent networks\. In the original SCS setting, this condition determines the boundary between chaotic and fixed\-point dynamics as the recurrent coupling strength g is varied\. Here, however, g is fixed and the transition is driven dynamically by learning\. The evolving feedback changes the stationary solutionuc​\(y^\)u\_\{c\}\(\\hat\{y\}\)through the self\-consistency condition

dd​u​V​\(u,y^\)=0⟹uc=g2​∫D​z​ϕ2​\(uc\)\+y^c2\\displaystyle\\frac\{d\}\{du\}V\(u,\\hat\{y\}\)=0\\implies u\_\{c\}=g^\{2\}\\int Dz\\phi^\{2\}\(\\sqrt\{u\_\{c\}\}\)\+\\hat\{y\}\_\{c\}^\{2\}\(24\)causing the system to move through the stability landscape until it reaches the marginal pointucu\_\{c\}\. This occurs when the output of the network reaches a specific valuey^\\hat\{y\}\. As shown in Fig\.[1](https://arxiv.org/html/2609.19288#S4.F1)\(f\), for g=1\.3, this number isy^=0\.2\\hat\{y\}=0\.2and in general will increase with increasing value of g\. We notes that this number is independent of the learning rateα\\alpha\. On the other hand, the specific time at which this happens during the learning trajectory defines a critical timetcrt\_\{\\text\{cr\}\}which will be a function of g andα\\alpha\. This requires us to solve explicitly for the time dependence ofy^​\(t\)\\hat\{y\}\(t\), which we do in Section[5](https://arxiv.org/html/2609.19288#S5)\.

Thus, the freezing transition does not arise from a change in the local stability criterion itself, but from the learning\-induced evolution of the correlation structure that brings the system to the SCS marginal state\. Solving Eq\.[20](https://arxiv.org/html/2609.19288#S4.E20)numerically allows us to computeΔ⁡\(y^,τ\)\\Delta\(\\hat\{y\},\\tau\)and consequentlyC⁡\(y^,τ\)C\(\\hat\{y\},\\tau\)forτ≪1/α\\tau\\ll 1/\\alpha\. The numerical results will be presented later in Section[6](https://arxiv.org/html/2609.19288#S6)after we finish the analysis for the far off diagonal elements and the time evolution ofy^\\hat\{y\}\.

Now we consider the off\-diagonal elements ofΔ⁡\(t,s\)\\Delta\(t,s\)in the regiont−s≥1/α≫τct\-s\\geq 1/\\alpha\\gg\\tau\_\{c\}\. The fast correlations have decayed away so all that remains is a slow evolution of the plateau\. Hence all the time derivatives in Eq\.[17](https://arxiv.org/html/2609.19288#S3.E17)scale asα≪1\\alpha\\ll 1and can be ignored\. In this case the two\-time evolution equation reduces to a self consistent equation

Δ⁡\(p,q\)\\displaystyle\\Delta\(p,q\)=\\displaystyle=g2​∫D​z​D​y​D​x​ϕ​\(Δ⁡\(p,p\)−Δ⁡\(p,q\)​x\+z​Δ⁡\(p,q\)\)​ϕ​\(Δ⁡\(q,q\)−Δ⁡\(p,q\)​y\+z​Δ⁡\(p,q\)\)\\displaystyle g^\{2\}\\int DzDyDx\\phi\(\\sqrt\{\\Delta\(p,p\)\-\\Delta\(p,q\)\}x\+z\\sqrt\{\\Delta\(p,q\)\}\)\\phi\(\\sqrt\{\\Delta\(q,q\)\-\\Delta\(p,q\)\}y\+z\\sqrt\{\\Delta\(p,q\)\}\)\(25\)\+\\displaystyle\+p​q\\displaystyle pqwherep≡y^​\(t\)p\\equiv\\hat\{y\}\(t\),q≡y^​\(s\)q\\equiv\\hat\{y\}\(s\)\. We have already solved for the diagonal valueΔ⁡\(p,p\),Δ⁡\(q,q\)\\Delta\(p,p\),\\Delta\(q,q\)which then allows us to solve for the far off diagonal correlation function which we callΔslow​\(p,q\)\\Delta\_\{\\text\{slow\}\}\(p,q\)\. The corresponding two time correlation functionCslow​\(p,q\)C\_\{\\text\{slow\}\}\(p,q\)can then be computed as well\.

Given the separation of time scales, we now approximate the fullΔ⁡\(t,s\)\\Delta\(t,s\)matrix as follows\. Starting from the diagonal value, we have a fast decay to the plateau dictated by the potential as shown in Fig\.[1](https://arxiv.org/html/2609.19288#S4.F1)\. This gives usΔfast​\(t,τ\)\\Delta\_\{\\text\{fast\}\}\(t,\\tau\)over a time scalet−s≪1/αt\-s\\ll 1/\\alpha\. We then match this to a slowly evolving plateauΔslow​\(t,s\)\\Delta\_\{\\text\{slow\}\}\(t,s\), over time scalest−s≥1/αt\-s\\geq 1/\\alphadictated by Eq\.[25](https://arxiv.org/html/2609.19288#S4.E25)\. The matrixC⁡\(t,s\)C\(t,s\)is approximated in exactly the same manner allowing us to write

C⁡\(t,s\)≈\(Cfast​\(t,τ\)−Cfast​\(t,τ=∞\)\)\+Cslow​\(t,s\)\\displaystyle C\(t,s\)\\approx\\left\(C\_\{\\text\{fast\}\}\(t,\\tau\)\-C\_\{\\text\{fast\}\}\(t,\\tau=\\infty\)\\right\)\+C\_\{\\text\{slow\}\}\(t,s\)\(26\)
So far, we have this matrix as a function ofτ\\tauandy^\\hat\{y\}for short time scales around the diagonal and as function ofy^​\(t\),y^​\(s\)\\hat\{y\}\(t\),\\hat\{y\}\(s\)for longer times\. We still need to solve fory^\\hat\{y\}as a function of time to obtain the full time dependence which we turn to in the next section\.

## 5Time evolution of the system

In this section, we explicitly solve for the learning trajectory as a function of time\. There are four time scales in the problem, two of which are the initial parameters of the system namelyτc\\tau\_\{c\}and1/α1/\\alphawhich govern the fast dynamics and the slow feedback learning respectively\. We also have two emergent scales,tdecayt\_\{\\text\{decay\}\}, the time scale over which correlators undergo a fast decay to a plateau\. Finally we havetcrt\_\{\\text\{cr\}\}when the system transitions from chaotic to stable dynamics\. This naturally allows us to divide the time evolution into three regimes as follows\.

### 5\.1Evolution at early time

At early times,t≤tdecay≪1/αt\\leq t\_\{\\text\{decay\}\}\\ll 1/\\alphawhen the output is small, the feedback is not strong enough to overcome the chaotic fluctuations completely\. Hence the dynamics is sensitive to the fast fluctuations over the time scaleτc\\tau\_\{c\}\. The evolution of the output is governed by Eq\.[14](https://arxiv.org/html/2609.19288#S3.E14)\. Since the integral over s is capped by t, in this regime, the off diagonal correlation matrix C\(t,s\) is relevant only over the intervalt−s≪1/αt\-s\\ll 1/\\alpha\. Similarly the diagonal valueΔ⁡\(t,t\)≈Δ⁡\(0,0\)\\Delta\(t,t\)\\approx\\Delta\(0,0\)is almost a constant over this small time interval since no appreciable learning has happened yet\. Therefore the system is essentially in an SCS chaotic state\. In that case we can write

⟨y⁡\(t\)⟩≈α​A​∫0td​s​Cfast​\(t=0,t−s\)\\displaystyle\\langle y\(t\)\\rangle\\approx\\alpha A\\int\_\{0\}^\{t\}dsC\_\{\\text\{fast\}\}\(t=0,t\-s\)\(27\)where we have also assumed that⟨y⁡\(t\)⟩\\langle y\(t\)\\rangleis small and so the error is approximately given by the target A\.

⟨y⁡\(t\)⟩=A​α​χ​\(t\)whereχ⁡\(t\)=∫0td​s​Cfast​\(t=0,t−s\)\\displaystyle\\langle y\(t\)\\rangle=A\\alpha\\chi\(t\)\\ \\ \\text\{where\}\\ \\ \\chi\(t\)=\\int\_\{0\}^\{t\}dsC\_\{\\text\{fast\}\}\(t=0,t\-s\)\(28\)To proceed further we make explicit choicesA=NA=\\sqrt\{N\}andα=0\.007\\alpha=0\.007\. This ensures that the feedback to the reservoir when it is fully trained is O\(1\) andα\\alphais sufficiently small for our separation of time scales to hold\. Our solution for the early time evolution is theny^​\(t\)=α​χ​\(t\)\\hat\{y\}\(t\)=\\alpha\\chi\(t\)and is shown in Fig\.[2](https://arxiv.org/html/2609.19288#S5.F2)\(b\)\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/EarlyEvol.png)Figure 2:Evolution at early timet≪1/αt\\ll 1/\\alpha\. \(a\) shows the fast decay of the SCS correlatorC⁡\(0,τ\)C\(0,\\tau\)\. \(b\) shows the time evolution of the network output at early times dictated by this fast decay\.We see in Fig\.[2](https://arxiv.org/html/2609.19288#S5.F2)\(a\) that the SCS fast decay happens over a times scale oftdecay∼t\_\{\\text\{decay\}\}\\sim20\-30τc\\tau\_\{c\}\. This is much smaller than1/α≈140​τc1/\\alpha\\approx 140\\tau\_\{c\}\. We want to use this approximate solution fory^\\hat\{y\}upto a time where the SCS potential has not deformed appreciably\. So we choosey^​\(25\)=\.033\\hat\{y\}\(25\)=\.033as the boundary for this early time evolution where the fast correlator has decayed to the plateau\. We see that over the time interval\{0,25\}​τc\\\{0,25\\\}\\tau\_\{c\}, wheny^\\hat\{y\}rises from 0 to0\.0330\.033, the potential is still very close to the SCS potential as shown in Fig\.[1](https://arxiv.org/html/2609.19288#S4.F1)\.Δ⁡\(t,τ=0\)≈0\.4\\Delta\(t,\\tau=0\)\\approx 0\.4, whileΔ⁡\(t,∞\)\\Delta\(t,\\infty\)rises by a very small value to0\.050\.05\.

### 5\.2Evolution at intermediate time

We now consider the intermediate\-time regime which corresponds to the intervaltcr≥t≥tdecayt\_\{\\text\{cr\}\}\\geq t\\geq t\_\{\\text\{decay\}\}\. In this regime, both fast and slow dynamics are important\. As explained in the last section, given the separation of scales, we can approximate the correlation matrix as a matched solution of the fast and slow matrix Eq\.[26](https://arxiv.org/html/2609.19288#S4.E26)\. This enables us to write the output using Eq\.[14](https://arxiv.org/html/2609.19288#S3.E14),

⟨y⁡\(t\)⟩=α​∫0td​s​\(A−y⁡\(s\)\)​Cslow​\(t,s\)\+\(A−y⁡\(t\)\)​α​∫0td​s​\(Cfast−Cfast​\(t,τ=∞\)\)\\displaystyle\\langle y\(t\)\\rangle=\\alpha\\int\_\{0\}^\{t\}ds\\left\(A\-y\(s\)\\right\)C\_\{\\text\{slow\}\}\(t,s\)\+\\left\(A\-y\(t\)\\right\)\\alpha\\int\_\{0\}^\{t\}ds\(C\_\{\\text\{fast\}\}\-C\_\{\\text\{fast\}\}\(t,\\tau=\\infty\)\)\(29\)Since the second term, governed by fast dynamics, only has support over time scalet−s≪1/αt\-s\\ll 1/\\alpha, the output essentially remains constant aty⁡\(t\)y\(t\)and hence can be pulled out of the integral\. This also allows us to set the upper limit of integration for this second term to∞\\inftysince in this regime we are only looking att≥tdecayt\\geq t\_\{\\text\{decay\}\}\. Since⟨y⁡\(t\)⟩\\langle y\(t\)\\rangleis a monotonic function of time, we can make a change of variables top=⟨y⁡\(t\)⟩p=\\langle y\(t\)\\ranglein the first term which allows us to write

p=α​∫0pd​v​\(A−q\)​f′​\(q\)​Cslow​\(p,q\)\+α⁡\(A−p\)​K​\(p\)\\displaystyle p=\\alpha\\int\_\{0\}^\{p\}dv\(A\-q\)f^\{\\prime\}\(q\)C\_\{\\text\{slow\}\}\(p,q\)\+\\alpha\(A\-p\)K\(p\)\(30\)wheref⁡\(p\)≡y−1​\(p\)=tf\(p\)\\equiv y^\{\-1\}\(p\)=tand the kernelK⁡\(p\)K\(p\)is known completely from the fast dynamics

K⁡\(p\)=∫0∞d​τ​\(Cfast​\(p,τ\)−Cfast​\(p,τ=∞\)\)\\displaystyle K\(p\)=\\int\_\{0\}^\{\\infty\}d\\tau\\left\(C\_\{\\text\{fast\}\}\(p,\\tau\)\-C\_\{\\text\{fast\}\}\(p,\\tau=\\infty\)\\right\)\(31\)We again note that since the integrand in the definition of K\(p\) goes to zero over a time scaletdecayt\_\{\\text\{decay\}\}, we can safely take the limit of integration to∞\\inftywhen looking at evolution at time t larger than the typical fast dynamics decay time\. We have then reduced the problem to solving for the functionf⁡\(p\)=y−1​\(p\)=tf\(p\)=y^\{\-1\}\(p\)=t, which when inverted gives us the solution for the learning trajectory\. This can be done by discretizing the integral and implementing a finite difference for the derivatives

pn=α​∑i=1n\(A−pi\)​\(fi\+1−fi\)​Cslow​\(pn,pi\)\+α⁡\(A−pn\)​K​\(pn\)\\displaystyle p\_\{n\}=\\alpha\\sum\_\{i=1\}^\{n\}\(A\-p\_\{i\}\)\(f\_\{i\+1\}\-f\_\{i\}\)C\_\{\\text\{slow\}\}\(p\_\{n\},p\_\{i\}\)\+\\alpha\(A\-p\_\{n\}\)K\(p\_\{n\}\)\(32\)Rearranging, we can solve forfn\+1f\_\{n\+1\}given all values tillpnp\_\{n\}andfnf\_\{n\}through the equation

fn\+1=fn\+un−α⁡\(A−pn\)​K​\(pn\)−α​∑i=1n−1\(A−pi\)​\(fi\+1−fi\)​Cslow​\(pn,pi\)α⁡\(A−pn\)​Cslow​\(pn,pn\)\\displaystyle f\_\{n\+1\}=f\_\{n\}\+\\frac\{u\_\{n\}\-\\alpha\(A\-p\_\{n\}\)K\(p\_\{n\}\)\-\\alpha\\sum\_\{i=1\}^\{n\-1\}\(A\-p\_\{i\}\)\(f\_\{i\+1\}\-f\_\{i\}\)C\_\{\\text\{slow\}\}\(p\_\{n\},p\_\{i\}\)\}\{\\alpha\(A\-p\_\{n\}\)C\_\{\\text\{slow\}\}\(p\_\{n\},p\_\{n\}\)\}\(33\)Since we knowp=⟨y⁡\(t\)⟩p=\\langle y\(t\)\\rangleat early times, upto the decay of fast dynamics \(∼t=25​τc\\sim t=25\\tau\_\{c\}forg=1\.3g=1\.3\), we already know the functionf⁡\(p\)=tf\(p\)=ttill that time\. This is the initial condition for Eq\.[33](https://arxiv.org/html/2609.19288#S5.E33)which we can then march forward to solve for all subsequent values offif\_\{i\}\.

### 5\.3Evolution at late time

As the output of the network increases, chaotic fluctuations are increasingly suppressed, until they completely disappear beyondt=tcrt=t\_\{\\text\{cr\}\}\. Hence, beyond this time, the kernelK⁡\(u\)K\(u\)goes to zero and we only have contribution from slow evolution\. So we can continue using Eq\.[33](https://arxiv.org/html/2609.19288#S5.E33)settingK⁡\(u\)K\(u\)to 0 beyondu=ucu=u\_\{c\}\(Eqns\.[23](https://arxiv.org/html/2609.19288#S4.E23)and[24](https://arxiv.org/html/2609.19288#S4.E24)\) , which is 0\.2 for g=1\.3\.

fn\+1=fn\+un−α​∑i=1n−1\(A−ui\)​\(fi\+1−fi\)​Cslow​\(un,ui\)α⁡\(A−un\)​Cslow​\(un,un\)\\displaystyle f\_\{n\+1\}=f\_\{n\}\+\\frac\{u\_\{n\}\-\\alpha\\sum\_\{i=1\}^\{n\-1\}\(A\-u\_\{i\}\)\(f\_\{i\+1\}\-f\_\{i\}\)C\_\{\\text\{slow\}\}\(u\_\{n\},u\_\{i\}\)\}\{\\alpha\(A\-u\_\{n\}\)C\_\{\\text\{slow\}\}\(u\_\{n\},u\_\{n\}\)\}\(34\)
This now allows to solve for the function f throughout the learning process, allowing us to access the temporal evolution of the network output as well as the two\-time correlation function\.

## 6Simulation and comparison with DMFT

Now we present the numerical results comparing simulation with the DMFT prediction\. To validate our theory, we choose a representative value of g=1\.3 simulated for a network of N= 5000 neurons and a constant targetA=NA=\\sqrt\{N\}\. The network was initially prepared in an SCS state by evolving it without any feedback\. Fig\.[3](https://arxiv.org/html/2609.19288#S6.F3)shows 5 representative runs along with the DMFT prediction after the feedback and learning are implemented\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/ztrajRaw.png)Figure 3:Network output as a function of time for five representative realizations with different onset times\. The dashed curve shows the DMFT prediction\. The histogram in \(f\) shows the distribution of onset times over 33 realizations\.We note that at finite N, for each run there is a variable time lag before macroscopic learning commences\. We also show the histogram for the onset time of learning over 33 runs in Fig\.[3](https://arxiv.org/html/2609.19288#S6.F3)\. The DMFT gives a deterministic trajectory and so does not capture run to run variability in the onset time of learning\. We can clearly identify the critical timetcrt\_\{\\text\{cr\}\}from the output trajectory as the time when the fluctuations completely disappear and the output rises smoothly\. This happens aty^=0\.2\\hat\{y\}=0\.2for g=1\.3 as predicted by the DMFT analysis\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/CttRaw.png)Figure 4:Diagonal correlation C\(t,t\) for five representative realizations\. The dashed curve shows the DMFT prediction\.Fig\.[4](https://arxiv.org/html/2609.19288#S6.F4)shows the diagonal entries of the correlation matrixC⁡\(t,t\)C\(t,t\)for the same representative runs\. We see the the same time lag after which the correlation matrix starts rising above the SCS value of 0\.23\. We also overlay the DMFT prediction for comparison\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/shifted.png)Figure 5:\(a\) Network output as a function of time for representative runs shifted by their onset times\. Also shown is the critical timetcrt\_\{\\text\{cr\}\}as predicted by the DMFT \(b\) The Diagonal Correlation matrixC⁡\(t,t\)C\(t,t\)as a function of time for the same representative runs shifted by their onset times\.We now shift each run by its corresponding onset time\. We determine the onset times of the various runs by demanding that they reach the outputy^=0\.8\\hat\{y\}=0\.8at the same time\. In Fig\.[5](https://arxiv.org/html/2609.19288#S6.F5), we show the overlap of the 5 representative runs shifted by their onset times, along with the DMFT prediction\. We also show the critical timetcrt\_\{\\text\{cr\}\}and the correspoding critical value of the feedbacky^c\\hat\{y\}\_\{c\}as predicted by DMFT, which does agree quite well with the data\. Finally we plot the mean and the standard deviation for 33 runs in Fig\.[6](https://arxiv.org/html/2609.19288#S6.F6)for the network output and the diagonal correlation matrixC⁡\(t,t\)C\(t,t\)\. We find excellent agreement with the DMFT prediction which validates the theory\.

![Refer to caption](https://arxiv.org/html/2609.19288v1/MeanPlot.png)Figure 6:Mean and standard deviation over 33 realizations for the network output \(left\) and diagonal correlation C\(t,t\) \(right\), after alignment by the realization\-dependent onset time\. The dashed curves show the DMFT predictions\.![Refer to caption](https://arxiv.org/html/2609.19288v1/OffDc.png)Figure 7:Off\-diagonal elements of the correlation matrixC⁡\(t,s\)C\(t,s\)along with the standard deviation for representative times before \(\(a\) and \(b\)\) and after \(\(c\) and \(d\)\) critical time\. Dashed curves are the DMFT predictions\.Finally, we plot the off diagonal correlation functionC⁡\(t,s\)C\(t,s\)for representative choices of t and overt−s∈\{0,t\}t\-s\\in\\\{0,t\\\}in Fig\.[7](https://arxiv.org/html/2609.19288#S6.F7)\. We look at two slices of the off\-diagonal elements fort<tcrt<t\_\{\\text\{cr\}\}\(Fig\.[7](https://arxiv.org/html/2609.19288#S6.F7)\(a\) and \(b\)\)\. As predicted by the DMFT, in this chaotic regime, as we move away from the diagonal, we clearly see a fast decay to a plateau over a time scaletdecay∼20​τct\_\{\\text\{decay\}\}\\sim 20\\tau\_\{c\}followed by slow evolution of the plateau to 0 at s=0\. The off\-diagonal elements are noisier especially in the chaotic regime, but the mean agrees quite well with our approximation\. Fort\>tcrt\>t\_\{\\text\{cr\}\}, the fast decay is absent, signifying frozen fast dynamics and we only see a slow learning induced evolution of the plateau\.

## 7Conclusion

In this work, we have developed a non\-equilibrium dynamical mean\-field theory for recurrent neural networks undergoing slow, feedback\-driven learning\. Going beyond characterizing the final trained network, we followed the evolution of the network dynamics throughout learning\. This allowed us to describe how plasticity progressively reorganizes the statistical state of the recurrent network and drives a transition between dynamical regimes\. A central result of our analysis is that learning generates an evolving effective feedback strength that acts as a dynamical control parameter\. As the learned feedback develops, the effective dynamical landscape is progressively deformed, driving the network from an initially chaotic regime towards a dynamically stable state\. The theory identifies a well\-defined critical feedback strength and a corresponding learning time marking the transition between these regimes\. By retaining the full two\-time structure of the correlation function, rather than assuming stationarity, the theory captures both the rapid relaxation of fluctuations and the slow evolution of the correlation plateau associated with the learning process\.

The transition identified here is a non\-equilibrium dynamical transition rather than an equilibrium thermodynamic transition\. The effective potential provides a convenient representation of the dynamical mean\-field equations and should not be interpreted as a thermodynamic free energy\. The two dynamical phases are distinguished by the presence or absence of persistent fast fluctuations and are separated by a bifurcation of the self\-consistent DMFT solutions\.

We further demonstrated that the resulting DMFT predictions quantitatively describe the learning dynamics observed in finite\-size simulations\. In particular, after accounting for realization\-dependent fluctuations in the onset time, the ensemble\-averaged trajectories and two\-time correlations agree with the theoretical prediction\. These finite\-size variations are consistent with fluctuations around the macroscopic learning trajectory described by the DMFT, but instead lead to fluctuations in the time at which individual realizations commence macroscopic learning\.

More broadly, this work suggests that dynamical mean\-field theory can provide a useful framework for connecting microscopic plasticity rules to macroscopic changes in the dynamical state of recurrent networks\. The two\-time description developed here makes it possible to characterize learning trajectories beyond the stationary states typically considered in analyses of trained networks, and provides a route toward studying how different forms of plasticity shape the temporal organization of recurrent computation\.

An important direction for future work is to extend this framework to learning rules that modify the recurrent connectivity itself, rather than only the readout and feedback pathways considered here\. Such plasticity would allow the network to develop persistent internal structure and long\-term memory, raising the question of how previously acquired dynamical states influence subsequent learning dynamics and the ability of the network to acquire new tasks\. Extending the present non\-equilibrium DMFT framework to this setting could provide a way to study the interaction between memory formation, dynamical phase transitions, and continual learning in recurrent networks\.

## Acknowledgements

V\.V\. would like to thank Dr\. Rodrigue Rizk for his valuable comments on the manuscript\. V\.V\. is supported by startup funds from the University of South Dakota and by the U\.S\. Department of Energy, EPSCoR program under contract No\. DE\-SC0025545\.

## Appendix ADerivation of DMFT equation

The generating functional which enforces the equation of motion of the neurons can be written as a path integral over the neuron degrees of freedom:

Z⁡\(J,Wf​b\)=∫∏i𝒟​xi​𝒟​x~i​exp⁡\{i​∑i∫d​t​x~i​\(t\)​\[τc​x˙i​\(t\)\+xi​\(t\)−∑j\(Ji​j\+Wif​b​wjo​u​t\)​ϕ​\(xj​\(t\)\)\]\}Z\(J,W^\{fb\}\)=\\int\\prod\_\{i\}\\mathcal\{D\}x\_\{i\}\\mathcal\{D\}\\tilde\{x\}\_\{i\}\\;\\exp\\Bigg\\\{i\\sum\_\{i\}\\int dt\\,\\tilde\{x\}\_\{i\}\(t\)\\Big\[\\tau\_\{c\}\\dot\{x\}\_\{i\}\(t\)\+x\_\{i\}\(t\)\-\\sum\_\{j\}\\left\(J\_\{ij\}\+W\_\{i\}^\{fb\}w\_\{j\}^\{out\}\\right\)\\phi\(x\_\{j\}\(t\)\)\\Big\]\\Bigg\\\}\(35\)
we can perform an ensemble average over the network parametersJi​j,Wif​bJ\_\{ij\},W\_\{i\}^\{fb\}, using the the distributions

P⁡\(Ji​j\)=g2​π​N​exp​\{−N2​g2​Ji​j2\},P⁡\(Wif​b\)=σfb2​π​N​exp​\{−N2​σfb2​\(Wif​b\)2\}\\displaystyle P\(J\_\{ij\}\)=\\frac\{g\}\{\\sqrt\{2\\pi N\}\}\\exp\{\\\{\-\\frac\{N\}\{2g^\{2\}\}J\_\{ij\}^\{2\}\\\}\},\\ \\ \\ \\ P\(W\_\{i\}^\{fb\}\)=\\frac\{\\sigma\_\{\\text\{fb\}\}\}\{\\sqrt\{2\\pi N\}\}\\exp\{\\Big\\\{\-\\frac\{N\}\{2\\sigma\_\{\\text\{fb\}\}^\{2\}\}\(W^\{fb\}\_\{i\}\)^\{2\}\\Big\\\}\}\(36\)
Z¯\\displaystyle\\bar\{Z\}=\\displaystyle=∫D​J​∫D​Wf​b​P​\(J\)​P​\(Wf​b\)​Z​\(J,Wf​b\)\\displaystyle\\int DJ\\int DW^\{fb\}P\(J\)P\(W^\{fb\}\)Z\(J,W^\{fb\}\)\(37\)=\\displaystyle=∫∏i𝒟xi𝒟x~iexp\{i∑i∫dtx~i\(t\)\[τcx˙i\(t\)\+xi\(t\)\]\\displaystyle\\int\\prod\_\{i\}\\mathcal\{D\}x\_\{i\}\\mathcal\{D\}\\tilde\{x\}\_\{i\}\\;\\exp\\Bigg\\\{i\\sum\_\{i\}\\int dt\\,\\tilde\{x\}\_\{i\}\(t\)\\Big\[\\tau\_\{c\}\\dot\{x\}\_\{i\}\(t\)\+x\_\{i\}\(t\)\\Big\]−\\displaystyle\-∑i∫dt∫dt′x~i\(t\)x~i\(t′\)\[∑jg22​Nϕ\(xj\(t\)\)ϕ\(xj\(t′\)\)\+σf​b22​Ny\(t\)y\(t′\)\]\}\\displaystyle\\sum\_\{i\}\\int dt\\int dt^\{\\prime\}\\tilde\{x\}\_\{i\}\(t\)\\tilde\{x\}\_\{i\}\(t^\{\\prime\}\)\\Big\[\\sum\_\{j\}\\frac\{g^\{2\}\}\{2N\}\\phi\(x\_\{j\}\(t\)\)\\phi\(x\_\{j\}\(t^\{\\prime\}\)\)\+\\frac\{\\sigma^\{2\}\_\{fb\}\}\{2N\}y\(t\)y\(t^\{\\prime\}\)\\Big\]\\Bigg\\\}The path integral is done using a saddle point approximation, where the saddle point is given by

x~i∗​\(t\)\\displaystyle\\tilde\{x\}\_\{i\}^\{\*\}\(t\)=\\displaystyle=i​∫d​s​C¯−1​\(t,s\)​\(x˙​\(s\)​τc\+x⁡\(s\)\)whereC¯−1​\(t,s\)solves\\displaystyle i\\int ds\\bar\{C\}^\{\-1\}\(t,s\)\(\\dot\{x\}\(s\)\\tau\_\{c\}\+x\(s\)\)\\ \\ \\text\{where\}\\ \\ \\bar\{C\}^\{\-1\}\(t,s\)\\ \\ \\text\{solves\}∫d​s​C¯​\(t,s\)​C¯−1​\(s,u\)=δ⁡\(t−u\)and\\displaystyle\\int ds\\bar\{C\}\(t,s\)\\bar\{C\}^\{\-1\}\(s,u\)=\\delta\(t\-u\)\\ \\ \\ \\text\{and\}C¯​\(t,s\)\\displaystyle\\bar\{C\}\(t,s\)=\\displaystyle=g2N​∑jϕ⁡\(xj​\(t\)\)​ϕ​\(xj​\(s\)\)\+σfb2N​y​\(t\)​y​\(s\)\\displaystyle\\frac\{g^\{2\}\}\{N\}\\sum\_\{j\}\\phi\(x\_\{j\}\(t\)\)\\phi\(x\_\{j\}\(s\)\)\+\\frac\{\\sigma^\{2\}\_\{\\text\{fb\}\}\}\{N\}y\(t\)y\(s\)\(38\)The saddle point is on the imaginary axis but the integral overx~i​\(t\)\\tilde\{x\}\_\{i\}\(t\)is along the real line\. We therefore choose a contour that passes through the saddle point and orients along the path of steepest descent which in this case points along direction of the real line\. So we simply shift the contour to be parallel to the real line passing throughx~i​\(t\)∗\\tilde\{x\}\_\{i\}\(t\)^\{\*\}noting that we do not cross any poles and the contributions from the segments at infinity vanish\. This yields the result

Z¯=∏i∫D​xi​exp⁡\{i​∑i∫d​t​𝑑s​\[τc​x˙i​\(t\)\+xi​\(t\)\]​C¯−1​\(t,s\)​\[τc​x˙i​\(s\)\+xi​\(s\)\]\}\\displaystyle\\bar\{Z\}=\\prod\_\{i\}\\int Dx\_\{i\}\\exp\\Bigg\\\{i\\sum\_\{i\}\\int dtds\\Big\[\\tau\_\{c\}\\dot\{x\}\_\{i\}\(t\)\+x\_\{i\}\(t\)\\Big\]\\bar\{C\}^\{\-1\}\(t,s\)\\Big\[\\tau\_\{c\}\\dot\{x\}\_\{i\}\(s\)\+x\_\{i\}\(s\)\\Big\]\\Bigg\\\}\(39\)Next we introduce conjugate fieldC^​\(t,s\)\\hat\{C\}\(t,s\)which enforces the definition for the collective fieldC¯​\(t,s\)\\bar\{C\}\(t,s\),

Z¯\\displaystyle\\bar\{Z\}=\\displaystyle=∫D​C^​D​C¯​exp⁡\{i​∑∫d​t​𝑑s​C¯​\(t,s\)​C^​\(t,s\)\}\\displaystyle\\int D\\hat\{C\}D\\bar\{C\}\\exp\\Bigg\\\{i\\sum\\int dtds\\bar\{C\}\(t,s\)\\hat\{C\}\(t,s\)\\Bigg\\\}\(40\)∏iDηiexp\{−12∑i∫dtdsηi\(t\)C¯−1\(t,s\)ηi\(s\)−i∫dtdsC^\(t,s\)\[g2N∑jϕ\(xj\(t\)\)ϕ\(xj\(s\)\)\+σfb2Ny\(t\)y\(s\)\]\}\\displaystyle\\prod\_\{i\}D\\eta\_\{i\}\\exp\\Bigg\\\{\-\\frac\{1\}\{2\}\\sum\_\{i\}\\int dtds\\eta\_\{i\}\(t\)\\bar\{C\}^\{\-1\}\(t,s\)\\eta\_\{i\}\(s\)\-i\\int dtds\\hat\{C\}\(t,s\)\\Big\[\\frac\{g^\{2\}\}\{N\}\\sum\_\{j\}\\phi\(x\_\{j\}\(t\)\)\\phi\(x\_\{j\}\(s\)\)\+\\frac\{\\sigma^\{2\}\_\{\\text\{fb\}\}\}\{N\}y\(t\)y\(s\)\\Big\]\\Bigg\\\}≡\\displaystyle\\equiv∫D​C^​D​C¯​exp⁡\{i​∑∫d​t​𝑑s​C¯​\(t,s\)​C^​\(t,s\)\}​exp​\{ln⁡Z^​\(C¯,C^\)\}\\displaystyle\\int D\\hat\{C\}D\\bar\{C\}\\exp\\Bigg\\\{i\\sum\\int dtds\\bar\{C\}\(t,s\)\\hat\{C\}\(t,s\)\\Bigg\\\}\\exp\\Bigg\\\{\\ln\\hat\{Z\}\(\\bar\{C\},\\hat\{C\}\)\\Bigg\\\}Next we find the saddle point solution by minimizing the action overC^\\hat\{C\}which yields

C¯​\(t,s\)=g2N​∑i⟨ϕ⁡\(xi​\(t\)\)​ϕ​\(xi​\(s\)\)⟩Z^​\(C¯,C^\)\+σfb2N​⟨y⁡\(t\)​y​\(s\)⟩Z^​\(C¯,C^\)\\displaystyle\\bar\{C\}\(t,s\)=\\frac\{g^\{2\}\}\{N\}\\sum\_\{i\}\\langle\\phi\(x\_\{i\}\(t\)\)\\phi\(x\_\{i\}\(s\)\)\\rangle\_\{\\hat\{Z\}\(\\bar\{C\},\\hat\{C\}\)\}\+\\frac\{\\sigma^\{2\}\_\{\\text\{fb\}\}\}\{N\}\\langle y\(t\)y\(s\)\\rangle\_\{\\hat\{Z\}\(\\bar\{C\},\\hat\{C\}\)\}\(41\)i\.e\., the matrixC¯\\bar\{C\}is set to its average value computed self consistently\. When we substitute this back into our path integral, using the central limit theorem, we can ignore the variance inC¯\\bar\{C\}in the large N limit, which then effectively reduces our generating functional upto a normalizing factor to

Z^=∏iDηiexp\{−12∑i∫dtdsηi\(t\)C¯−1\(t,s\)ηi\(s\)\}\\displaystyle\\hat\{Z\}=\\prod\_\{i\}D\\eta\_\{i\}\\exp\\Bigg\\\{\-\\frac\{1\}\{2\}\\sum\_\{i\}\\int dtds\\eta\_\{i\}\(t\)\\bar\{C\}^\{\-1\}\(t,s\)\\eta\_\{i\}\(s\)\\Bigg\\\}\(42\)whereC¯​\(t,s\)\\bar\{C\}\(t,s\)is average value in Eq\.[41](https://arxiv.org/html/2609.19288#A1.E41)that is determined self consistently using the Gaussian actionZ^\\hat\{Z\}\.

## References

- \(1\)H\. Sompolinsky, A\. Crisanti, and H\. Sommers, “Chaos in random neural networks,”Physical Review Letters61no\. 3, \(1988\) 259–262\.
- \(2\)E\. Daucé, M\. Quoy, B\. Cessac, B\. Doyon, and M\. Samuelides, “Self\-organization and dynamics reduction in recurrent networks: stimulus presentation and learning,”[Neural Networks11no\. 3, \(1998\) 521–533](http://dx.doi.org/10.1016/S0893-6080(97)00131-7)\.
- \(3\)R\. Laje and D\. V\. Buonomano, “Robust timing and motor patterns by taming chaos in recurrent neural networks,”[Nature Neuroscience16\(2013\) 925–933](http://dx.doi.org/10.1038/nn.3405)\.
- \(4\)D\. Sussillo and L\. F\. Abbott, “Generating coherent patterns of activity from chaotic neural networks,”[Neuron63no\. 4, \(2009\) 544–557](http://dx.doi.org/10.1016/j.neuron.2009.07.018)\.
- \(5\)F\. Mastrogiuseppe and S\. Ostojic, “Linking connectivity, dynamics, and computations in low\-rank recurrent neural networks,”[Neuron99no\. 3, \(2018\) 609–623\.e29](http://dx.doi.org/10.1016/j.neuron.2018.07.003)\.
- \(6\)A\. Rivkind and O\. Barak, “Local dynamics in trained recurrent neural networks,”[Physical Review Letters118no\. 25, \(2017\) 258101](http://dx.doi.org/10.1103/PhysRevLett.118.258101)\.
- \(7\)D\. G\. Clark, B\. Bordelon, J\. A\. Zavatone\-Veth, and C\. Pehlevan, “Structure, disorder, and dynamics in task\-trained recurrent neural circuits,”[bioRxiv\(2026\)](http://dx.doi.org/10.64898/2026.03.02.708943)\.
- \(8\)H\. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks,” Tech\. Rep\. 148, German National Research Center for Information Technology, 2001\.
- \(9\)D\. Koryakin, J\. Lohmann, and M\. V\. Butz, “Balanced echo state networks,”[Neural Networks36\(2012\) 35–45](http://dx.doi.org/10.1016/j.neunet.2012.08.008)\.
- \(10\)J\. M\. Murray, “Local online learning in recurrent networks with random feedback,”[eLife8\(2019\) e43299](http://dx.doi.org/10.7554/eLife.43299)\.
- \(11\)T\. Miconi, “Biologically plausible learning in recurrent neural networks reproduces neural dynamics observed during cognitive tasks,”[eLife6\(2017\) e20899](http://dx.doi.org/10.7554/eLife.20899)\.
- \(12\)T\. Asabuki and C\. Clopath, “Taming the chaos gently: a predictive alignment learning rule in recurrent neural networks,”[Nature Communications16\(2025\) 6784](http://dx.doi.org/10.1038/s41467-025-61309-9)\.
- \(13\)P\. C\. Martin, E\. D\. Siggia, and H\. A\. Rose, “Statistical dynamics of classical systems,”[Physical Review A8\(1973\) 423–437](http://dx.doi.org/10.1103/PhysRevA.8.423)\.
- \(14\)H\. K\. Janssen, “On a lagrangean for classical field dynamics and renormalization group calculations of dynamical critical properties,”[Zeitschrift für Physik B23\(1976\) 377–380](http://dx.doi.org/10.1007/BF01316547)\.
- \(15\)C\. De Dominicis, “Technics of field renormalization and dynamics of critical phenomena,”Journal de Physique Colloques37\(1976\) C1–247–C1–253\.
- \(16\)A\. Crisanti and H\. Sompolinsky, “Path integral approach to random neural networks,”[Physical Review E98\(2018\) 062120](http://dx.doi.org/10.1103/PhysRevE.98.062120)\.