Resonance-like enhancement of forced nonlinear diffusion as a nonequilibrium phase transition

We describe the phenomenon of a resonance-like, giant enhancement of diffusion in a basic model of nonlinear diffusion featured by a nonlinear in velocity friction and the corresponding multiplicative thermal noise. The model is consistent with thermal equilibrium in the absence of driving. Different from previous studies of this phenomenon, where the crucial nonlinearity originates from a periodic external potential while friction is linear, we focus on the case of a constant force driving, whereas the crucial nonlinearity stems from the friction. The basic model of such friction considered interpolates between linear viscous Stokes friction at small velocities and dry Coulomb-like friction at large velocities corresponding to a stress plateau in some nonlinear viscoelastic materials. Recently, a nonequilibrium phase transition to super-diffusion and super-transport was discovered within this basic model. We show that adding a tiny viscous friction part to major nonlinear friction regularizes in part this behavior. Diffusion becomes asymptotically normal. However, the phase transition translates into a giant enhancement of normal diffusion and mobility of particles at the transition point over the intuitively expected large force limit, where the linearization of friction occurs. Such a giant enhancement of diffusion is closely related to the largely enhanced kinetic temperature of the particles at and beyond the critical point. We provide analytical results obtained within an effective mass approximation which nicely agree with stochastic numerics.


Introduction
Giant enhancement of diffusion is a phenomenon whereby a driven beyond thermal equilibrium diffusion can be hugely enhanced over the limit of force-free thermal diffusion in a resonance-like manner at some critical force [1][2][3][4][5][6][7][8][9][10][11]. It was first theoretically discovered for overdamped diffusion in tilted washboard potentials [1,2] and later confirmed experimentally [12,13]. The corresponding critical force is associated with a critical potential tilt separating localized and sliding inertialess solutions at zero temperature [1,2]. The akin phenomenon was also discovered for underdamped diffusion in the same nonlinear potentials, where its physical origin is entirely different [4][5][6]8]. For very small friction coefficients, it occurs at much smaller values of forces which correspond to equal probabilities of localized and running inertial solutions [5,6,8] within a bistable velocity pseudo-potential picture [14], which can be traced back to pioneering contributions by Risken [15] with colleagues [16]. The diffusion enhancement over the inertialess limit becomes much more prominent due to inertial effects. Moreover, taking the hydrodynamic memory effects into account makes this enhancement dramatically even larger [9,10] and prolongs profoundly transient regimes of anomalously fast diffusion and transport [6,9,11]. In the inertial case, it is generally associated with a significant enhancement of the kinetic temperature of particles [5,9,17,18] defined by the width of the velocity distribution [17][18][19][20]. For an intermediate friction strength, it can become multi-peaked, as discovered in references [9][10][11]. Moreover, hydrodynamic memory effects can make the inertial effects in the giant enhancement of diffusion profound even for seemingly (critically) overdamped diffusion [10]. All this previous work was done for a linear in velocity model of friction.
Nonlinear models of friction lead to the so-called nonlinear Brownian motion [15,21,22], which is much less investigated until now than the Brownian motion featured by linear velocity friction in nonlinear spatial potentials. The latter one is also nonlinear but in a profoundly different way. The former emerges in quite different research areas from diffusion in nonlinear condensed media, like various polymer and colloidal solutions [23][24][25][26][27][28][29], relativistic Brownian motion [22,30], which are at thermal equilibrium unless being driven out of it by external forces, to cold atoms diffusion in optical potentials [31][32][33][34][35] and active Brownian particles [36], where it is profoundly out of thermal equilibrium always. For such active Brownian motion featured by bistable velocity fluctuations, a giant enhancement of diffusion was also found at the condition when the potential minima of a velocity pseudo-potential are nearly equal [36], similar to the inertial case in tilted washboard potential [5,6,8].
Our focus in this paper is on nonlinear Brownian motion consistent with thermal equilibrium in the absence of driving force, where at a small constant driving force, linear mobility and zero-force thermal diffusion coefficient are related by the thermal fluctuation-dissipation theorem (FDT) [15,[37][38][39]. The nonlinearity of friction greatly influences the diffusion coefficient, particularly its temperature dependence, but not the linear-response FDT. For a significant force, however, the system becomes strongly driven beyond thermal equilibrium, and the particles become kinetically heated, with their stationary velocity distribution (when exists) broadening beyond the Maxwell distribution at the environment's temperature [39]. This broadening can be characterized by nonequilibrium kinetic temperature [17,19,20]. We consider the model of friction similar to the one in reference [39] which bridges between linear viscous friction at small velocities and dry Coulomb-like friction at large velocities. The later feature corresponds to the stress plateau at large shear velocities in some viscoelastic nonlinear media [23][24][25][26][27][28]. The shear rate and particle velocity, as well as the nonlinear medium's viscosity and the friction acting on the particle, are proportional. In viscoelastic media, the memory effects are generally substantial. They will be, however, neglected in our model study. We consider a general theoretical model, which is of profound interest per se, as a basic model of nonlinear Brownian motion and statistical physics [40]. Recently, we have shown [39] that there exists a nonequilibrium phase transition [40,41] in this model at a critical driving force, where the stationary kinetic temperature diverges. Both diffusion and transport become asymptotically anomalously fast and nonergodic [35] at the edge of phase transition and beyond it [39]. The model serves as a vivid illustration of the supposition [40] that anomalous diffusion [31][32][33][34]42] can emerge by a phase transition. Earlier, a similar anomalous diffusion regime was found (for a different model of friction) in a model of nonlinear Brownian motion applied to ultracold atoms in shallow optical potentials in the absence of bias, where the unbiased diffusion is profoundly out of thermal equilibrium per se [31][32][33][34][35].
For systems featured by a nonlinear shear-thinning viscosity, it is generally expected that diffusion will be enhanced by forcing over the equilibrium force-free limit. It is because the corresponding velocity-dependent friction coefficient becomes smaller upon forcing, and the diffusion coefficient is inversely proportional to friction by the Einstein relation [15,38]. This statement is evident if to think in terms of an effectively linearized problem with a driving-renormalized diffusion coefficient. However, can the diffusion enhancement become even more significant, giant, beyond this expectation, and resonance-like? This question we address and answer with this work.
Below we consider the same basic model of nonlinear friction as in reference [39] augmented, however, by a small viscous friction force added to the nonlinear friction, together with the corresponding thermal noise part to fulfill the FDT. Such an addition will be shown to regularize the asymptotic behavior of the transport and diffusion at and beyond the critical force. They become asymptotically normal, with anomalous transport and diffusion present as a transient phenomenon. However, a vast enhancement of mobility and diffusion coefficient occurs at the critical point. The smaller is an add-on, the larger are diffusion at the peak, and the mobility of particles after the phase transition took place, and the longer lasts anomalous diffusion and transport transient. Simple approximate, however, accurate analytical expressions for the mobility and diffusion coefficient are found beyond the transition point. Moreover, we show that the diffusion enhancement at and beyond the critical force is physically explained and related to the behavior of the kinetic temperature, which for vast forces reduces to the temperature of the environment, like in the opposite limit of small forces. Conceptually, thinking in terms of a dynamical or nonequilibrium phase transition remains useful and predictive also in the case considered.

Model and theory
Let us consider a nonlinear Brownian motion (NLBM) in an environment at thermal equilibrium characterized by temperature T, governed by the Langevin equation [15,21,22,43] for the velocity variable v =ẋ of a particle with one-dimensional coordinate x and inertial mass m. The particle is subjected to a constant force f, nonlinear frictional force f fr (v) = η(v)v and a thermal noise f th (v, t) = √ 2k B Tη(v)ξ(t) related to a formally introduced white Gaussian noise ξ(t) of unit intensity. It is worth noting that thermal noise is not Gaussian for nonlinear friction, η(v) = const. Nevertheless, the Gaussian model of ξ(t) within a multiplicative noise setting remains useful, despite another line of thinking including non-Gaussian thermal environment is also essential [44]. It is, however, beyond the NLBM scope of this work which provides a generic enough theoretical framework.
The well-known problem with equations like equation (1) [45][46][47][48] arises at evaluating stochastic integral t+Δt t f th (v(t ), t )dt involving singular delta-correlated noise ξ(t), which is featured by an infinite root-mean-square (rms) amplitude at any integration time step Δt. Without some pertinent interpretation such equations make no perfect sense [45][46][47][48]. Generally, the proper interpretation depends on the physics of the problem considered as greatly advocated by van Kampen and others [45][46][47][48][49]. For consistency with thermal equilibrium at zero-forcing, f = 0, this Langevin equation must be interpreted in the Hänggi-Klimontovich [21,22,46,47,49], or kinetic sense, named also anti-Ito [50,51] or post-point interpretation to distinguish it [52,53] from the Ito (pre-point) and Stratonovich (middle-point) interpretations [45,48]. Otherwise, at thermal equilibrium, the energy dissipated by friction would not be compensated on average by the energy gained from the thermal noise of the environment, with mean kinetic energy being equal to k B T/2 on average (1D model is considered) and satisfying the thermal equipartition theorem [37,38,54]. Moreover, a directed motion can emerge at f = 0 in violation of the second law of thermodynamics. It becomes obvious in the case of an asymmetric friction, η(−v) = η(v), which we show in appendix A. The second law of thermodynamics forbidding a 'Maxwell demon' [54][55][56] current at f = 0 fixes the physically correct mathematical interpretation of this stochastic equation. Its ensemble, master equation counter-part, the Fokker-Planck equation (FPE) [15] can be written in a kinetic form [20,21,39] where K(v) = k B Tη(v)/m 2 , upon introducing a velocity pseudo-potential 1 The kinetic form (2) of FPE and our way to introduce the velocity pseudo-potential make obvious that the stationary solution is always Divergence of statistical integral means that no stationary velocity distribution exists 2 . It signals a phase transition to superdiffusion, which may happen at some critical force f c [39]. However, for any model of η(v), the equilibrium velocity distribution at f = 0 is Maxwellian, Hence, at f = 0, in the non-driven case, the particle current is always absent, in agreement with thermodynamics and kinetic theory. Both Ito and Stratonovich interpretations of equation (1) can erroneously entail a particle current 3 .
We discuss and demonstrate it numerically in the appendix A for an asymmetric friction model, which is consistent with thermal equilibrium within the correct kinetic interpretation at use.
Using p st (f, v), any moment of stationary velocity distribution can be found as, For f = 0, we have immediately, v(f = 0) st = 0, at any T, consistently with thermal equilibrium, independently of a particular model of η(v). Henceforth, other interpretations of equation (1) are expected to generally imply a current at f = 0. Surely, for a symmetric η(−v) = η(v) it may be forbidden by this symmetry (if it is not spontaneously broken for some particular symmetric models of η(v)). Nevertheless, if this symmetry is violated, a nonphysical current does emerge immediately within physically incorrect interpretations of the Langevin equation considered, see appendix A. The second moment of distribution p(v, t) can be used to define the kinetic temperature T k (t) by is fluctuation of the velocity from its mean value. It is a common notion in classical kinetic theory [19,20] and has the same kinetic meaning as equilibrium temperature. Beyond equilibrium, the kinetic temperature of driven Brownian particles can largely exceed T [9,17]. The stationary limit of the kinetic temperature st ]/k B , can be found from equation (4) if this limit exists.
Also the diffusion coefficient D = (1/2) lim t→∞ δx 2 (t) /t characterizing the asymptotic behavior of the particles position variance, δx 2 2 can be found, when it exists, in terms of quadratures as [39] 4 where where η eff (T) = η −1 (v) −1 st is an effective friction coefficient obtained by averaging the inverse friction coefficient over equilibrium velocity fluctuations. It is especially useful in the linear response (LR) regime, where the inverse of it defines the LR-mobility μ(0, T) = lim f→0 ∂ v(f ) st /∂f = 1/η eff (T). This is also a very general result, independently of the model of nonlinear friction η(v).

The effective mass approximation
Some further general results can be obtained within an effective mass approximation. For this we notice that for any model of friction and a sufficiently small f, the minimum of U(v, f) at v = 0, will be shifted to some v m (f), and we can approximate in the vicinity of the minimum, by introducing an effective force-dependent mass m eff (f). Below we assume that it is the only minimum of U(v, f). Then, it is easy to show that within our model for any model of η(v) leading to monostable U(v, f). It is an effective mass approximation thus introduced. It is expected to work at sufficiently low temperatures. Then, any velocity moment can be found as U m (f) does not matter, as any additive constant in U(v, f). In particular, v(f) st = v m (f) within this approximation, which is clearly a T → 0 result. Moreover, from which, Hence, kinetic heating of Brownian particles with mass m can alternatively be interpreted as a thermal equilibrium behavior of particles with the effective mass m eff (f). 5 However, the first interpretation is more physical because driven Brownian particles are not at thermal equilibrium. Nevertheless, the concept of quasi-particles with some effective masses is of great importance in general in the condensed matter physics [61]. Hence, this interpretation is also helpful. Furthermore, for the diffusion coefficient within this approximation, we obtain where Unlike the averaged velocity, it is temperature-dependent within the same approximation. Moreover, in the limit of very low temperatures we can approximate η m (v) eff by η(v m ) to obtain a nice compact result whose importance will be clarified below. We shall name it the second effective mass approximation. One more general result can be obtained within this latter approximation for the autocorrelation time st . Indeed by the Kubo-Green relation [38], Hence, from this and (14) we obtain

Model of friction
Below we consider the model of friction with a = 1/2 similar to one in reference [39], which is augmented, however, by a constant friction term η 0 η a . The limit v v c , where v c is a critical velocity, corresponds to viscous Stokes friction η a + η 0 ; whereas in the limit v v c , it corresponds to a much smaller Stokes friction η 0 . Our model friction presents a particular case of friction corresponding to the Yasuda-Carreau model of nonlinear viscosity [23,25]. The striking feature of this model at η 0 = 0 is that it corresponds to a plateau in the elastic stress dependence on shear rate in some nonlinear viscoelastic materials [23][24][25][26][27][28], as the particle velocity is proportional to the shear rate, and the friction to viscosity. This plateau can be very extended in some materials [26]. However, the elastic stress starts to grow again for further increasing shear rate, which is reflected by a small η 0 in equation (17), cf figure B1. The ratio η a /η 0 can easily reach several thousand [26].
2.2.1. The case η 0 = 0 As a basic physical model of nonlinear friction, the considered model interpolates, at η 0 = 0, between viscous Stokes friction f fr ≈ η a v at small v and a dry-like Coulomb friction, f fr ≈ f c v/|v|, with f c ≈ η a v c at large velocities. Then, the corresponding velocity potential has a single minimum at [39] v 2 (18) 5 The kinetic heating of Brownian particles by a constant force in nonlinear media should not be confused with another interesting phenomenon of hot Brownian motion due to the heating of the Brownian particles and their microenvironment in laser fields [60].
for a subcritical tilt f < f c , which is also the mean transport velocity v(f) st in the approximation of effective mass. Let us improve now this approximation by replacing η a with the above-introduced temperature-dependent η eff (T), to obtain an approximation better suitable at finite T. Making this improvement, while replacing η(0) by η eff (T) in the results obtained within the second effective mass approximation will be named the improved effective mass approximation (IEMA).
In the case considered, the effective mass is exactly. At the critical tilt f c , it vanishes, and a nonequilibrium phase transition occurs to anomalous fast transport and diffusion, which also persist beyond the critical point. The transport can be interpreted so that, for f > f c , the frictional force becomes exactly f c v/|v|, i.e., a critical transition to the Coulomb-like friction occurs. Below critical point, both diffusion and transport are normal. Within the IEMA, we obtain the following analytical results for nonlinear mobility Notice that the mobility enhancement is larger that one would expect from the (12), (14) within the IEMA reads Notice that the diffusion enhancement is much more significant than one would expect from an effectively linearized theory relationship D LT (f, T) = k B Tμ(f, T) ad hoc applied to our profoundly nonlinear problem. Such a naive expectation is simply wrong. These insightful results nicely agree with numerics [39], where a power law exponent γ with a numerical fitting value close to γ = 5/2 was used instead of the exact value 5/2 derived in the present work. Moreover, we have within the same approximation, It also excellently agrees with numerics in [39] (denoted as T max therein), where unfortunately a misprint is present: power two at f/f c was missed. It is corrected here. Furthermore, from (16) we obtain within the IEMA. At the critical point, f = f c , transport becomes anomalously fast asymptotically with δx(t) ∼ t α , where α = 4/3, and diffusion is strongly superballistic with δx 2 (t) ∼ t 2α [39]. Beyond the critical point, f > f c , transport becomes asymptotically ballistic, δx(t) = at 2 /2, under acceleration a = (f − f c )/m. Next, diffusion is asymptotically weakly super-ballistic δx 2 (t) ∼ k B T r t 2 ln(t/t 1 )/m, with a temperature parameter, T r ≈ 2T/(f/f c − 1), which diverges at f = f c [39]. This weakly superballistic dependence is related to a logarithmic law of the kinetic temperature growth, T k (t) ∼ T r ln(t/t 0 ), with t 0 being some numerical scaling time parameter and t 1 = t 0 exp(3/2) [39].

2.2.2.
The case η 0 > 0 In this case and for η a > η 0 , the velocity pseudo-potential reads, This bulky expression precludes finding μ(f, T) and D(f, T) in an exact analytical form. However, good-quality approximate analytical results are still possible. Moreover, equation (24) is crucial for  1(a). These approximations do not work only in some small neighborhood of f c , which shrinks with ever-smaller η 0 . These features indicate that the phase transition from viscous to dry friction in the nonlinear friction part takes place also at a small η 0 , with f c in equation (25) corresponding to the dry friction force in the supercritical regime. However, a finite η 0 makes the stationary velocity for f > f c now finite, which explains the structure of equation (25). Asymptotically anomalous diffusion and transport disappear. However, one expects a peak of the diffusion coefficient and a sharp increase of the particles' mobility at some critical force, which is close to f c and depends in addition on η 0 and temperature. Hence, a sort of phase transition takes place also for a small but finite η 0 . The IEMA allows also in this case arrive at the insightful results. Indeed, we obtain exactly. It is a very cumbersome expression because of a bulky v m (f) entering it. Nevertheless, it fortunately allows for very simple approximations, in the parameter regime η a η 0 and |f/f c − 1| η 0 /η a . The smaller η 0 /η a , the better are approximations. In the subcritical regime, the approximation is given by equation (19), as expected. It works pretty well, see figure 1(b). In this case, the above approximate results for nonlinear mobility μ(f, T) and diffusion coefficient D(f, T), obtained for η 0 = 0, remain valid with the corresponding η eff (T) wherein accounting now for a small η 0 contribution as well.
In the supercritical regime, we obtain approximately, which immediately implies Notice that instead of m eff (f c ) = 0 at f c in both approximations, the true m eff (f) > 0 arrives at a minimum at some f m > f c near f c , see in figure 1(b), and then starts to grow in accordance with (27) and reaches asymptotically m eff (∞) = m. This growth is, however, algebraically slow. This approximation works also very well, cf figure 1(b). The nonlinear mobility in supercritical regime is accordingly It gradually approaches μ(∞) = 1/η 0 , i.e., a linearization occurs in the limit of vast driving forces. Furthermore, for the diffusion coefficient in the supercritical regime, we obtain with T st (f) in equation (28). It implies that (i) a resonance-like enhancement of the diffusion coefficient near to f c in supercritical regime over the asymptotic value, D(∞, T) = k B T/η 0 , is entirely on the count of increased kinetic temperature (one, major interpretation), or decreased effective mass (alternative interpretation), and (ii) VACF time does not depend on f being equal τ cor (∞) = m/η 0 in the supercritical regime. It is in a sharp contrast with the behavior of τ cor (f) in the subcritical regime, which diverges with the exponent −3/2, cf reference [39] and equation (23) at η 0 = 0, but is getting the final value m/η 0 at a finite η 0 . It is also worth mentioning that in the supercritical regime independently of f within this model, as deduced from equations (28)- (30). It is very different from the LR-regime, where both the mobility and the diffusion coefficient do not depend on the force applied and μ(0, T) = D(0, T)/(k B T). This linear proportionality within LR-theory rules an intuition in the classical transport theory [61]: the larger the diffusion coefficient, the larger is the mobility of particles. They are proportional. With equation (31), we see, however, that in the case of nonlinear friction, this intuition can be quite wrong in a strongly nonlinear response regime: nonlinear mobility can be inversely proportional to the diffusion coefficient.

Numerical results and discussion
Next, we performed two sorts of numerics to check the validity of our approximate analytical theory against the exact analytical results evaluated numerically and accurate stochastic simulations. First, we did a numerical evaluation of the theoretical results using equations (4)- (7), (17) and (24). Second, we also performed stochastic simulations of the Langevin equation (1)   , also the stationary kinetic temperature is plotted. The numerically-evaluated for η 0 = 1 exact theory results are plotted in (a)-(c) with green asterisk symbols. They nicely agree with the corresponding stochastic numerics. In (a), we also plotted numerical results for both values of η 0 together with the subcritical and supercritical approximations. They work very well except for the vicinity of the transition points. Far from it, the mean v st depends on f linearly, however, with two very different slopes and the offset f c above the critical point, as the inset in (a) demonstrates. In (b), the numerical results for the diffusion coefficient D are compared with the numerical values of the scaled stationary kinetic temperature T st /η 0 and the corresponding two analytical approximations. Both nicely work. The inset in (b) resolves a sharp diffusion and kinetic temperature peak for η 0 = 0.1. In the supercritical regime, up to the peak, the diffusion coefficient reflects completely the change of the kinetic temperature. It is not so in the subcritical regime, where D is nicely described by equation (21) far from the vicinity of the critical point. In (c), we resolve the D peak structure and compare the results with the full effective mass approximation. Astoundingly enough, it works very well at η 0 = 0.1 even for T = 1 in our numerics, which is not small. The simple analytical approximations also work well, apart from the vicinity of the peak in this case. At η 0 = 1, the agreement worsens, and, nevertheless, the approximation still works fairly well. Surprisingly, simple analytical result in equation (30) well estimates the maximal value of D at peak upon using the numerically found f max in it. done on professional graphics processing units (GPUs) NVIDIA V100 propagating 10 5 trajectories in parallel with double precision, where it takes about 92 min to integrate until t max = 1000.
An instance of the ensemble-averaged mean displacement of particles and their variance is shown in figure 2 for the tilting force f = 250, which is the critical force at η 0 = 0. The results for four different values of η 0 are depicted. For very small t 0.01, we have universally δx(t) = ft 2 /(2m) and δx 2 (t) = (v T t) 2 , as the friction plays initially no role.
For η 0 = 0, transport and diffusion quickly arrive at anomalously fast transport, t 4/3 , and diffusion, t 8/3 , asymptotic regimes [39]. Finite η 0 makes these regimes of finite duration. Asymptotically, both transport and diffusion are normal. However, the intermediate anomalous regime can last long, depending on η 0 . The smaller η 0 , the longer this regime is. For the smallest η 0 = 0.01 in figure 1 the asymptotic normal regime is still not achieved. However, for η 0 = 0.1 and especially η 0 = 1, it is clearly seen. By fitting the corresponding numerical regime with δx(t) = v(f) st t + c 1 and δx 2 (t) = 2D(f)t + c 2 , we extract the corresponding numerical values of v(f) st and D(f). Here, c 1,2 are just some fitting constants required to improve the fitting accuracy at finite t. They can be omitted for very large t.
Symbols depict the corresponding numerical values of mean velocity and diffusion coefficient for two values η 0 = 0.1 and η 0 = 1 in figure 3. For η 0 = 1, we also depict the results obtained from the numerical evaluation of the theory results. Both numerics remarkably agree. In this respect, one should mention that the numerical evaluation of the first two moments of velocity distribution is not problematic in general. However, the numerical evaluation of D(f) involves a double-quadrature and is numerically troublesome. We did it using a Gauss-Legendre quadrature in mpmath library [63], Python, while (i) shifting v by v st and U(v, f) by U( v st , f) in the corresponding integrals, (ii) using finite integration limits in (5) instead of infinite, and (iii) getting convergent results upon gradual extension of these finite limits. It was not successful for all parameters attempted and for many parameter values this procedure took more time than finding D(f) from the corresponding stochastic numerics. In such cases, one can perceive the numerical solution of equation (1), as a sort of indirect Monte Carlo integration of the double-quadrature in equations (5) and (6) with (24). The agreement of both numerics for η 0 = 1 in figure 3 is convincing.
From figure 3, one can deduce that the approximate analytical results obtained within the second effective mass approximation (improved in the subcritical regime) work very well except for some vicinity of the critical point. Of course, they cannot describe the value f max at which the diffusion coefficient arrives at its maximum. Here, the result with the exact m eff in equations (26) and (12) (30), when used with f max found numerically, also allows to estimate D(f max ) with a reasonable accuracy, see in figure 3(c). Indeed, in accordance with this approximation, D(f max ) ≈ 17.66 for η 0 = 1 and D(f max ) ≈ 724.28 for η 0 = 0.1. f max is always larger f c and gradually approaches f c with lowering η 0 . D(f max ) increases with ever smaller η 0 and tends to infinity, when η 0 → 0, following equation (30).
Next, in the supercritical regime, the enhancement of diffusion is indeed entirely on the count of the enhanced kinetic temperature, as the theory predicts. Here, we see a remarkable universality: independently of η 0 , the enhancement of the diffusion coefficient over its f → ∞ limit is always given by the factor Furthermore, in figure 4, we plotted dynamics of T k (t) at η 0 = 1 and several values of f shown in the plot. Several remarkable features are observed. First, for a subcritical force, e.g., f = 250 in this plot, T k (t) grows monotonously and arrives finally at the corresponding T st which value is used in figure 3(b). T st grows with f and arrives at it maximal value near to f max . In figure 4, it happens at f = 263 somewhat smaller than f max = 265. With further growing f, T st becomes smaller and T k (t) becomes nonmonotonous in time having a maximum at some intermediate time. In other words, T st is not the maximal kinetic temperature in supercritical regime. Both the value of this maximum and T st algebraically slow decay to the temperature of the bath with ever increasing f.
The stationary distribution of velocity is also insightful. It is shown in figure 5 for η 0 = 1 and four different values of f near-to or not-far-away-from the critical point. The exact analytical results are plotted by dashed-dotted red lines using the velocity potential (24), whereas the numerical results obtained from stochastic simulations are plotted with empty circle symbols. The non-Gaussian form of the velocity distribution is obvious, except from the supercritical case of f = 300, part (d), where it is nearly Gaussian. It indicates that p st (v, f) becomes quickly Gaussian behind the critical point in the supercritical regime.
We fitted non-Gaussian dependencies with an asymmetric stretched-to-compacted Gaussian or a generalized exponential distribution where i = 1 corresponds to the backward part and i = 2 to the forward part of this distribution correspondingly, with two half-width parameters v i and power-exponents b i . C is a normalization constant. It turns out that the forward part is always nearly Gaussian, b 2 = 2, i.e., it is a half-Gaussian distribution, in fact. It seems to be the reason why  (24) and also their fitting discussed in the text, with the parameters given in the plots. the effective mass approximation works unexpectedly well even for a non-small T = 1 in our numerics. However, the backward part can strongly deviate from the Gaussian b 1 = 2. For example, b 1 = 2.67 at f = 200, i.e. the Gaussian becomes compressed, see in part (a). At f = f c = 250, it becomes strongly compressed with b 1 = 4.68 (part (b)), and at f = f max = 265, the left central part becomes stretched, b 1 = 1.64. Its tail, however, decays very fast to zero, see in part (c). We reiterate that the velocity distribution rapidly becomes nearly Gaussian after the critical force f max is passed. Then, it is looking like a standard Maxwell distribution. However, the kinetic temperature entering it can strongly exceed the environment's temperature. For f = 300 in (d), T st ≈ 6, i.e., if T = 1 would correspond to a room temperature 300 K, then T st would be 1800 K, correspondingly! The particles are kinetically very hot, which is reflected by their largely enhanced diffusion, whereas the decay time constant of VACF remains the same m/η 0 . It is not affected by force at all.

Outlook and conclusion
The theory developed in this work is based essentially on the assumption of a single minimum of the velocity potential U(v, f) (corresponding to the unique solution of the force balance equation η(v)v = f) and the effective mass approximation which yields a series of insightful analytical results confirmed by numerics. The occurrence of a nonequilibrium phase transition in driven Brownian particles transport in conjunction with a greatly enhanced diffusion is related to the occurrence of a plateau in the frictional force vs the particle velocity, cf figure B1. It, in turn, is related to a corresponding stress plateau vs the shear rate in the viscoelastic media [26,27], where the phenomenon of a huge diffusion enhancement is expected to be observable. However, it should be stressed that in the corresponding real viscoelastic media, the memory effects are essential. For simplicity, they were omitted in our model study. We shall address them in a subsequent work, which requires a generalized Langevin equation approach [37,38,64,65] generalized toward a nonlinear in velocity friction, which is not a trivial task. Nevertheless, the experimental work is most welcome even without such memory effects in our current theory. It is expected to reveal the phenomenon of a resonance-like enhanced driven diffusion in such media predicted in this work. Its basic mechanism differs much from the original phenomenon of hugely enhanced diffusion in overdamped physical systems and its subsequent inertial generalizations.
Next, in some nonlinear viscoelastic media, the phenomenon of stress overshooting is also observed [28], i.e., the stress first arrives at a maximum vs the shear rate, then diminishes and rises again. It can be accounted for by a > 1/2 in equation (17). Then, however, the force balance equation can have three solutions for a range of f variation, see figure B1 for two illustrative cases, a = 2/3 and a = 1. In those cases, the corresponding velocity potential U(v, f) becomes bistable, which invalidates the single effective mass approximation and requires a further generalization of the theory. Nevertheless, even without a detailed theory to be developed in the future, it is clear, given the results in references [6,8,9,36] for two different models, that the phenomenon of a hugely enhanced diffusion will also persist in this case. It will occur at the force value when both minima of U(v, f) will be equal [6,8,9,36]. In such a case, a huge diffusion enhancement phenomenon will be reminiscent in physical mechanism to one of the underdamped diffusion enhancement in a tilted washboard potential and other setups leading to bistable or multistable velocity fluctuations. Combined with the memory effects, which should be further considered in nonlinear viscoelastic media, it opens intriguing perspectives for future research.
In conclusion, this paper introduced a physically well-grounded model predicting a giant enhancement of diffusion in forced Brownian motion in nonlinear physical media, seen as a nonequilibrium or dynamic phase transition. Further generalizations and ramifications of this basic model are expected to be crucial to describe stochastic transport phenomena in nonlinear viscoelastic media. The predicted phenomenon of a resonance-like enhanced forced diffusion in such media calls for experimental verification.
corresponding to the noise entering the discussed equation, switched in series, make together a closed electrical circuit, with the voltage noise source presenting voltage fluctuations in it. By this analogy, equation (1) describes the stochastic current I(t) dynamics in this electrical circuit. In the absence of emf, U 0 = 0, the averaged current must be absent at thermal equilibrium. The resistance can, in principle, be asymmetric, R(−I) = R(I) with diodes providing a pertinent example. No matter how big this asymmetry is, no mean current can emerge from a diode rectifying thermal voltage fluctuations across its ends. Its emergence would mean a Perpetuum mobile of the second kind violating the second law of thermodynamics. This physical demand fixes the physically correct kinetic interpretation of the Langevin equation. Namely, it leads to FPE (2) in the kinetic form, which allows making these statements precise.
Indeed, let us consider some asymmetric friction model (or a diode model in the electric circuit setup), e.g., which is (17) at η 0 = 0 and shifted by a constant v 0 . In this case we obtain, exactly, after some straightforward calculation. We recall that it corresponds to v(f) st within the effective mass approximation, which can be improved by replacing η a with the corresponding η eff (T). This dependence is very asymmetric for a sufficiently large v 0 , e.g., v 0 = v c . However, at any v 0 = 0, v(f = 0) st = 0: no rectification of thermal fluctuations is possible. Striking asymmetry of the considered setup implies immediately that in the Ito and Stratonovich interpretations of (1) a Maxwell demon current [54][55][56] will emerge at f = 0, which is forbidden by the second law of thermodynamics.

A.1. Stochastic simulations and the spurious drift
Maybe even more instructive is to show this by direct stochastic simulations of equation (1). Here, depending on the numerical method at use, one must add one or another spurious drift term to be consistent with this equation's kinetic interpretation. Only in the case of an implicit stochastic Euler method, there is no need for such spurious drift because such a method does correspond to the kinetic interpretation of the noise term integral at each integration time-step. The working horse of stochastic simulations and the simplest numerical method-the stochastic Euler method-does correspond to the Ito interpretation of (1) [62]. Hence, one must add a spurious drift term to (1) while using the stochastic Euler method to solve it, to have a physically meaningful solution of (1). Within the realm of overdamped Brownian dynamics, this fact was probably first realized by Ermak and McCammon [67]. Since then, it has been commonly used in Brownian molecular dynamics simulations. In our simulations, we use the second-order stochastic Heun method, which naturally corresponds to the Stratonovich interpretation of (1) [62]. Hence, one must add the spurious drift term f ghost (v) = (m/2)K (v) to f [39], to have (1) in the kinetic interpretation, while solving it by using the Heun method or subtract this term, to arrive at its Ito interpretation in effect. The results of the corresponding simulations at f = 0 and v 0 = v c = 1 and T = 1 are shown in figure A1. A nonphysical current of particles is absent only within the kinetic interpretation, in agreement with the theory.  Furthermore, as noticed in [39], for a > 1/2, the velocity potential may have a maximum for certain values of f and become then bistable for a finite η 0 . It is reflected by a nonmonotonous behavior of the frictional force η(v)v in figure B1 for two illustrative values of α = 2/3 and α = 1. At some f, there are three crossing points: two correspond to the minima and one to the maximum of the velocity potential, which becomes bistable. These interesting cases will be studied somewhere else.