# Local Plastic Response and Slow Heterogeneous Dynamics of Supercooled
Liquids

Yan-Wei Li<EMAIL_ADDRESS>Key Laboratory of Advanced Optoelectronic
Quantum Architecture and Measurement (MOE), School of Physics, Beijing
Institute of Technology, Beijing, 100081, China  Yugui Yao<EMAIL_ADDRESS>Key Laboratory of Advanced Optoelectronic Quantum Architecture and Measurement
(MOE), School of Physics, Beijing Institute of Technology, Beijing, 100081,
China  Massimo Pica Ciamarra<EMAIL_ADDRESS>Division of Physics and
Applied Physics, School of Physical and Mathematical Sciences, Nanyang
Technological University, Singapore 637371, Singapore  CNR–SPIN, Dipartimento
di Scienze Fisiche, Università di Napoli Federico II, I-80126, Napoli, Italy
CNRS@CREATE LTD, 1 Create Way, #08-01 CREATE Tower, Singapore 138602

###### Abstract

We demonstrate, via numerical simulations, that the relaxation dynamics of
supercooled liquids correlates well with a plastic length scale measuring a
particle’s response to impulsive localized perturbations and weakly to
measures of local elasticity. We find that the particle averaged plastic
length scale vanishes linearly in temperature and controls the super-Arrhenius
temperature dependence of the relaxation time. Furthermore, we show that the
plastic length scale of individual particles correlates with their typical
displacement at the relaxation time. In contrast, the local elastic response
only correlates with the dynamics on the vibrational time scale.

When supercooled, liquids develop a slow relaxation dynamics [1, 2, 3] that
proceeds through a sequence of activated events [4]. These relaxation events
occur in preferential locations fixed by the particle configuration [5],
rather than randomly in space. Furthermore, these events cluster in space and
propagate through a facilitation mechanism [6, 7, 8], making the dynamics
spatially heterogeneous on the relaxation time scale [1, 2, 3, 9, 10, 11, 12].
It is currently unclear what features of a particle configuration identify the
regions where relaxation events occur with higher probability.

Previous works suggested that the local structure surrounding each particle
determines its possible involvement in irreversible relaxation events and
analysed this structure introducing physically motivated structural parameter
[13, 14, 15, 16, 17, 13, 18] or via machine learning approaches [19, 20].
These structural parameters are system dependent as affected, e.g., by the
particles’ shape. Alternatively, local elastic rather than geometric
properties may have predictive ability, if the established correlation [21,
22, 23] between macroscopic elastic properties and relaxation dynamics also
holds at the microscopic scale. Considered measures of local elasticity
include the Debye-Waller factor (DWF), its harmonic approximation, and
parameters probing particles’ vibrational motion [24, 9, 25, 26, 27]. However,
the correlation of a particle’s local elastic properties with its involvement
in local relaxation processes is debated. Indeed, a particle’s vibrational
amplitude correlates highly with its short-time displacement [9] but weakly
with its displacement evaluated at the relaxation time scale [15]. Besides,
correlations between local elastic vibrations and plastic relaxation events
entail non-obvious assumptions [28] on the shape of the free energy basin the
system is transiently confined.

The response of individual particles to externally applied forces is an
alternative approach to characterize local mechanical properties. The response
to small forces, which probes the linear elastic regime, revealed the spatial
heterogeneity of the local elastic moduli in metallic glasses [29] and dense
colloidal suspensions [30]. Larger forces induce plastic flow and allow for
microrheology investigations [31, 32, 33, 34]. Recently, the response to
transient applied forces has been used to probe the emergence of caging in
colloidal suspensions [35] and energy absorption in soft colloids [36, 37].

In this Letter, we demonstrate that the plastic response induced by transient
localized perturbations defines a length, $\xi$, that correlates with the
relaxation dynamics at the macroscopic and at the individual particle level.
The average plastic length vanishes linearly with the temperature and
regulates the super-Arrhenius divergence of the relaxation time,
$\tau=\tau_{0}\exp(\xi_{0}/\xi)$. A particle’s local plastic response
correlates with its typical displacement at the relaxation time. On the
contrary, the local elastic response only correlates with a particle’s
displacement at shorter times, on the vibrational time scale. Our findings
demonstrate that structural relaxation proceeds via local rearrangements that
are weakly related to the local elastic response but well connected to the
local plastic response induced by transient perturbations.

We consider the Kob-Andersen [38, 39], 65(A):35(B) binary Lennard-Jones
mixture of $N=2000$ particles in two spatial dimensions. The interaction
potential is
$U_{\alpha\beta}(r)=4\epsilon_{\alpha\beta}[(\frac{\sigma_{\alpha\beta}}{r})^{12}-(\frac{\sigma_{\alpha\beta}}{r})^{6}+C_{\alpha\beta}]$
for $r\leq r_{\alpha\beta}^{c}=2.5\sigma_{\alpha\beta}$ and
$U_{\alpha\beta}(r)=0$ otherwise, $\alpha,\beta\in$ $\\{\rm{A,B}\\}$. We set
$\sigma_{AB}/\sigma_{AA}=0.8$, $\sigma_{BB}/\sigma_{AA}=0.88$,
$\epsilon_{AB}/\epsilon_{AA}=1.5$, $\epsilon_{BB}/\epsilon_{AA}=0.5$, and fix
$C_{\alpha\beta}$ so that the potential vanishes continuously at the cutoff.
Length, energy and time are recorded in units of $\sigma_{AA}$,
$\epsilon_{AA}$ and $\sqrt{m\sigma_{AA}^{2}/\epsilon_{AA}}$, respectively. We
perform Langevin dynamics simulations using lammps [40], numerically
integrating the equations of motion,
$m\mathbf{\ddot{r}}_{i}=-\mathbf{\nabla}_{i}\sum_{j(\neq
i)}U_{\alpha\beta}(r_{ij})-\gamma\mathbf{\dot{r}}_{i}+\mathbf{\eta}_{i}(t)$,
where $\mathbf{r}_{i}$ is the position of the $i$th particle, $r_{ij}$ is the
inter-particle distance, and $\gamma=1$ is the friction coefficient.
$\mathbf{\eta}_{i}$ is a random noise, satisfying
$\langle\mathbf{\eta}_{i}(t)\mathbf{\eta}_{j}(t^{\prime})\rangle=2k_{B}T\gamma\delta_{ij}\delta_{t-t^{\prime}}\mathbf{1}$,
with $\mathbf{1}$ the unit tensor. This model has been shown to be a good
glass former [39], the relaxation time developing a super-Arrhenius
temperature dependence below the onset temperature $T_{\rm onset}=0.7$ [41,
42], for number density $\rho=1.2$. We recap the main features of the
relaxation dynamics of this model in the Supplementary Material (SM) [42].

We probe the local mechanical response by exerting to a particle $i$ of an
equilibrium configuration a force of magnitude $f$ and random orientation in
the time interval $0<t<t_{\rm p}$, mimicking recent experimental studies [35,
30, 37, 36]. The response of the system depends on $f$, $t_{p}$ and the
damping parameter $\gamma$. In the following, we concentrate on the $f$
dependence of the response at fixed $\gamma=1$ and $t_{p}=0.1$ corresponding
to $\approx 1/10$ of the vibrational timescale. In the SM [42], we show that
results are robust with respect to changes in $\gamma$ and $t_{p}$ as long as
the perturbation triggers an irreversible response.

We illustrate in Fig. 1(a) the time dependence of the displacement $\Delta
x(t)=\langle\Delta x_{i}(t_{\rm obs})\rangle_{\theta,N_{\rm r},c}$ induced by
a local force of magnitude $f$. The displacement is averaged over $160$ forces
differing in their random orientation $\theta$ for each particle, and then
further averaged over $N_{\rm r}=200$ randomly selected particles and $c=8$
independent configurations. The displacement attains a constant value $L$ for
a transient which approximately ends when the mean square displacement becomes
of order $L^{2}$. Below the onset temperature, this transient begins at
$t\simeq 2\times 10^{2}$, corresponding to roughly $t_{\rm obs}\simeq 200$
vibrational times.

The time dependence of the displacement informs on the elastic properties of
the systems. Indeed, in a viscous medium, the displacement induced by a
transient localized force increases monotonically and saturates to a plateau
value, $L\propto f$. In Fig. 1(a), where dashed lines mark the value of $L$,
we observe this behaviour at the largest of the considered temperatures. In
contrast, in a viscoplastic medium, the displacement reaches a maximum and
then saturates to a smaller $L$ value, as we observe at the other considered
temperatures. In the presence of a purely elastic response, the final
displacement $L$ would vanish.

The magnitude of the average force-induced displacement $L(T,f)$ (see Fig. S4
[42] for the distribution of the individual particle displacements) depends on
the applied force and temperature. Below the onset temperature the asymptotic
displacement increases linearly with the temperature, $L(T,f)\propto(T-T_{0})$
and with the applied force, $L\propto(f-f_{0})$, as illustrated in Fig. 1(b)
and (c). Here, $T_{0}\simeq 0.12$ and $f_{0}$ is a temperature independent
threshold value. This scaling does not hold at higher temperatures (see Fig.
S3 [42]), and indeed, a small deviation already occurs at the onset
temperature in Fig. 1(c). The need of applying a minimum force to induce
irreversible particle motion is consistent with previous results [34].

Figure 1:  (a) Time dependence of the average displacement of a particle in
response to a force acting on the time window identified by the shaded region.
The dashed lines mark the asymptotic displacement values, $L(T,f)$. The filled
red stars mark the $\Delta x(\tau_{\rm CR})$ for $T\geq 0.5$. (b) The
asymptotic displacement grows linearly in temperature below the onset
temperature marked by the vertical dashed line, $L\propto(T-T_{0})$. (c) A
plot of $A(f)=L(T,f)/(T-T_{0})$ reveals that, below the onset temperature, the
asymptotic displacement grows linearly in $f$ for large enough forces.

To rationalize the force dependence of $L$, we investigate the response to the
application of localized forces of systems brought in energy minima, in the SM
[42]. This investigation suggests that the temperature independence of $f_{0}$
reflects the weak dependence of the energy landscape on the temperature of the
parent liquid, in the supercooled regime, and that $L$ has a linear dependence
on $f$ as it results from the accumulation of $n(f)\propto(f-f_{0})$ plastic
displacements of typical size $\langle l\rangle$ [42].

Summarizing, Fig. 1 shows that in the regime of low-temperature and forces
large enough to trigger plastic rearrangements, the force-induced displacement
varies as $L(T,f)=\xi(T)(1-f/f_{0})$, with $\xi(T)\propto(T-T_{0})$ a
temperature dependent length scale. Next, we show that this plastic length
regulates the relaxation dynamics.

We probe the relaxation dynamics using cage-relative (CR) measures to filter
out the influence of long-wavelength fluctuations [43, 44, 45, 46, 41] that
lead to an underestimation of the structural relaxation time in two spatial
dimensions. While standard relaxation measures focus on the actual
displacement of each particle, CR measures focus on the displacement of a
particle with respect to the average displacement of its neighbours, which we
identify via the Voronoi construction. This definition makes CR measures
unaffected by the collective translations that characterise long-wavelength
fluctuations.

We estimate the CR relaxation time, $\tau_{\rm CR}(T)$, which is proportional
to the shear viscosity [41, 47], from the decay of the CR self-scattering
function evaluated at the wave-vector corresponding to the first peak of the
structure factor. Figure 2 illustrates that, below the onset temperature, the
CR relaxation time grows exponentially with the inverse plastic length scale
$\xi(T)$,

$\tau_{\rm CR}(T)=\tau_{0}\exp\left(\frac{\xi_{0}}{\xi(T)}\right)$ (1)

with $\tau_{0}$ a microscopic vibrational timescale and $\xi_{0}$ a constant.
Given that $\xi\propto(T-T_{0})$, Eq. 1 corresponds to a Vogel-Fulcher-Tammann
(VFT) dependence of the relaxation time on temperature, which we illustrate in
Fig. S2 of the SM [42].

We stress that other functional forms may describe the super-Arrhenius
temperature dependence of the relaxation time. As such, we do not want to give
fundamental significance to the fitting parameters. Yet, since
$\tau=\tau_{0}\exp(\Delta F/T)$, where $\Delta F(T)$ is the free energy
barrier of the elementary relaxation events, our finding demonstrates a
connection between the plastic length scale $\xi(T)$ and $\Delta F(T)$, i.e.,
$\Delta F\propto 1/\xi$. This is our first main result.

Figure 2:  The cage-relative relaxation time grows exponentially with the
reciprocal plastic length $1/\xi(T)=(1-f/f_{0})(1-T/T_{0})/L(T,f)$. Data refer
to $\rho=1.2$.

A VFT-like divergence of the relaxation time is predicted by entropic theories
of the glass transition, such as the early Adam-Gibbs [48] scenario and the
subsequent Random First Order Transition theory [49]. According to these
theories, the super-Arrhenius dynamics follows from the vanishing of the
configurational entropy, or equivalently from the divergence of a static
length which measures the size of loosely defined cooperative rearranging
regions [50, 51]. Previous works supported the existence of a growing static
length on cooling [52, 13, 18, 53, 11, 54, 55]. These theories are consistent
with our results only if the dynamical slowdown involves both a decreasing and
an increasing length scale. This scenario may occur if relaxation events have
a core-corona structure, with the core and corona size respectively decreasing
and increasing on cooling. In this respect, we remark that particles
undergoing large displacements in a relaxation event become increasingly
localized [56, 57] on cooling.

Elastic models of the glass transitions [58, 21] relates $\Delta F$ to a
vanishing length scale, the DWF, e.g., $\Delta F\propto 1/\langle
u^{2}\rangle^{\alpha/2}$ with $\alpha$ a phenomenological parameter describing
cage anisotropy [59, 23] or related functional forms [24, 22]. While both
$\xi$ and $\langle u^{2}\rangle$ decrease on cooling, these lengths are
conceptually different. $\xi$ relates to the plastic response and it appears
to vanish at a finite temperature. In contrast, $\langle u^{2}\rangle$ probes
the elastic, possibly anharmonic, vibrational motion and only vanishes at zero
temperature. Indeed, to recover a super-Arrhenius behaviour, earlier works
[24] suggested the relaxation time relates to the difference between $\langle
u^{2}\rangle$ and its expected crystalline value.

We now turn our attention to the spatial and temporal heterogeneity of the
relaxation dynamics [1, 2, 3], which we quantify by resorting to the CR mean
square displacement in the iso-configurational ensemble [5]. We evaluate
$\langle\Delta r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ by averaging the mean
square displacement of each particle over $128$ equilibrium simulations
sharing the same initial particle configuration and differing in the
particles’ momenta, we randomly sample from the relevant equilibrium
distribution. Particles with a large $\langle\Delta r^{2}_{i,\rm
CR}(t)\rangle_{\rm iso}$ at the relaxation time identify the locations more
prone to structural rearrangements [5, 9].

We investigate the correlations between $\langle\Delta r^{2}_{i,\rm
CR}(t)\rangle_{\rm iso}$ and a local plastic length scale
$\xi_{i}=\langle\Delta x_{i}(t_{\rm obs})\rangle_{\theta}$, we define as the
displacement induced by the application of a randomly oriented force of
magnitude $f=140$ on particle $i$, averaged over $200$ realizations. To
compare with previous works [9, 15, 6], we also investigate how $\langle\Delta
r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ correlates with the CR-DWF $\langle
u_{i}^{2}\rangle=\langle\Delta r^{2}_{i,\rm CR}(\tau_{\beta})\rangle_{\rm
iso}$, the mode participation $p_{i}$ [9], and the harmonic mean square
displacement $\psi_{i}$ [15]. Here $\tau_{\beta}$ is the time at which the
logarithmic derivative of the mean square displacement acquires its minimum,
implying that the system is maximally subdiffusive [22]. To evaluate $p$ and
$\psi$, we diagonalize the Hessian matrix of the disordered solid obtained by
minimizing the energy of the $t=0$ configuration via the Tinker package [60],
obtaining $2N-2$ eigenmodes $\mathbf{e}(\omega_{a})$ with non-zero
eigenfrequencies $\omega_{a}$. We define the participation fraction of
particle $i$ as
$p_{i}=\langle|\mathbf{e}_{i}(\omega_{j})|^{2}\rangle_{N_{m}}$, where $j$ runs
over the $N_{m}=50$ modes with lowest non-zero frequency. The mean square
displacement in the harmonic approximation is
$\psi_{i}=\sum_{a}\omega_{a}^{-2}\left|\mathbf{e}_{i}(\omega_{a})\right|^{2}$.
At each time $t$, we resort to Spearman’s rank correlation coefficient $S(t)$
[61, 62, 11, 63] to evaluate how $\langle\Delta r^{2}_{\rm CR}(t)\rangle_{\rm
iso}$ correlates to each of the different considered static quantities. $S=1$
($-1$) for monotonically increasing (decreasing) dependence of two quantities,
while $S=0$ in the absence of correlations.

Figure 3 illustrates (right axis) the time evolution of the average CR mean
square displacement at $\rho=1.2$ and $T=0.4$, and (left axis) the time
evolution of Spearman’s rank correlation coefficients. By definition,
$\langle\Delta r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ and $\langle
u_{i}^{2}\rangle$ are perfectly correlated at $t=\tau_{\beta}$
($\tau_{\beta}\simeq 3\times 10^{-4}\tau_{\rm CR}$ at $\rho=1.2$ and $T=0.4$).
Figure 3 clarifies that the correlation sharply drops at longer times [15],
equaling $\simeq 0.1$ at the relaxation time scale. Similarly, the correlation
between $\langle\Delta r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ and both
participation fraction and mean square displacement in the harmonic
approximation reaches a maximum at short times and sensibly drops on the
relaxation time scale.

These results agree with previous investigations [9] on the relation between
iso-configurational mean square displacement and participation fraction, that
found correlations at $t\simeq 0.3\tau$, with $\tau$ the relaxation time
evaluated without resorting to CR measures [64]. Indeed, for the system
considered in Fig. 3, $t=0.3\tau$ corresponds to $t\simeq 0.1\tau_{\rm CR}$, a
short time where the dynamics is still dominated by the vibrational motion.
The critical conclusion we reach from this investigation is that participation
fraction and DWF correlate well with the vibrational dynamics but weakly with
the dynamics on the relaxation time scale. Henceforth, these quantities are
not a good proxy for the location of the relaxation events.

Figure 3: Time dependence of Spearman’s rank correlation coefficient $S$
between the CR mean square displacement of the particles at time $t$ and 1)
their CR mean square displacement at $t=\tau_{\beta}$; 2) their participation
fraction $p$; 3) their mean-square displacement evaluated in the harmonic
approximation and 4) their plastic length scale $\xi$. The blue dashed line
represents the average CR mean square displacement.

In Figure 3, we also illustrate Spearman’s rank correlation coefficient
between $\langle\Delta r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ and the plastic
length scale $\xi_{i}$. This correlation coefficient is small at short times
and reaches a maximum at the relaxation time. Hence, it critically differs
from the other correlation coefficients we have investigated that peak on the
vibrational time scale. These findings demonstrate that the plastic length
scale is a robust indicator of the relaxation dynamics at the particle level.
This result is our second significant result.

To further test the validity of our results, we have replaced $\langle\Delta
r^{2}_{i,\rm CR}(t)\rangle_{\rm iso}$ with alternative measures of structural
relaxation, a particle-based self-scattering function and the fraction of
neighbors the particle has lost at time $t$. Fig S11 [42] shows that our
finding are robust with respect to the adopted definition of local structural
relaxation. In addition, in the SM we show that the average plastic length
scale correlates with the relaxation time, and the locally defined plastic
length scale identifies the particles more prone to rearrange, also when the
dynamical slowdown is induced by an increase of the density at a constant
temperature. These investigations support the robustness of our findings and
clarify that our predictions could be experimentally tested in suspensions of
colloidal particles [35, 30, 37, 36], where density is the main control
parameter.

Our results show that the response of individual particles to transient,
almost impulsive perturbations defines a plastic correlation length that
regulates supercooled liquids’ relaxation dynamics. The average plastic
correlation length $\xi$ is inversely proportional to the free energy barrier
$\Delta F(T)$ governing the dynamic slowdown. The plastic response of
individual particles informs on their participation in localized relaxation
events, and hence on the spatial heterogeneity of the relaxation process at
the relaxation time scale. Our findings provide novel insights into the origin
of slow and heterogeneous dynamics, showing that the local elastic response
mostly correlates with the vibrational dynamics rather than with the dynamics
at the relaxation time. Our results may inspire novel glass transition
theories and stimulate experimental research on the plastic response to
external forces induced, e.g., by laser pulses.

We thank M. Wyart for fruitful discussions. Y.-W.L. acknowledges support of
the start-up funding of Beijing Institute of Technology and support from the
National Natural Science Foundation (NSF) of China (Grants No. 12105012). Y.Y.
is supported by the NSF of China (Grants Nos. 11734003, 12061131002). MPC
acknowledges support from the Singapore Ministry of Education through the
Singapore Academic Research Fund (MOE2019-T1-001-03 & MOE2019-T1-001-22), and
the National Supercomputing Centre Singapore (NSCC) for the computational
resources.

## References

  * Ediger [2000] M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000).
  * Debenedetti and Stillinger [2001] P. Debenedetti and F. Stillinger, Nature 410, 259 (2001).
  * Berthier and Biroli [2011] L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
  * Pica Ciamarra  _et al._ [2015] M. Pica Ciamarra, R. Pastore, and A. Coniglio, Soft Matter 12, 358 (2015).
  * Widmer-Cooper  _et al._ [2004] A. Widmer-Cooper, P. Harrowell, and H. Fynewever, Phys. Rev. Lett. 93, 135701 (2004).
  * Candelier  _et al._ [2010] R. Candelier, A. Widmer-Cooper, J. K. Kummerfeld, O. Dauchot, G. Biroli, P. Harrowell, and D. R. Reichman, Phys. Rev. Lett. 105, 135702 (2010).
  * Keys  _et al._ [2011a] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Phys. Rev. X 1, 021013 (2011a).
  * Pastore  _et al._ [2015] R. Pastore, A. Coniglio, and M. Pica Ciamarra, Soft Matter 11, 7214 (2015).
  * Widmer-Cooper  _et al._ [2008] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. Reichman, Nat. Phys 4, 711 (2008).
  * Keys  _et al._ [2011b] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Phys. Rev. X 1, 021013 (2011b).
  * Tong and Tanaka [2018] H. Tong and H. Tanaka, Phys. Rev. X 8, 011041 (2018).
  * Hu  _et al._ [2018] Y.-C. Hu, Y.-W. Li, Y. Yang, P.-F. Guan, H.-Y. Bai, and W.-H. Wang, Proc. Natl. Acad. Sci. U. S. A. (2018), 10.1073/pnas.1802300115.
  * Kawasaki  _et al._ [2007] T. Kawasaki, T. Araki, and H. Tanaka, Phys. Rev. Lett. 99, 215701 (2007).
  * Patrick Royall  _et al._ [2008] C. Patrick Royall, S. R. Williams, T. Ohtsuka, and H. Tanaka, Nat. Mater. 7, 556 (2008).
  * Jack  _et al._ [2014] R. L. Jack, A. J. Dunleavy, and C. P. Royall, Phys. Rev. Lett. 113, 095703 (2014).
  * Furukawa [2018] A. Furukawa, Phys. Rev. E 97, 022615 (2018).
  * Li  _et al._ [2016] Y.-W. Li, Y.-L. Zhu, and Z.-Y. Sun, Phys. Rev. E 94, 062601 (2016).
  * Watanabe  _et al._ [2011] K. Watanabe, T. Kawasaki, and H. Tanaka, Nat Mater 10, 512 (2011).
  * Cubuk  _et al._ [2015] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, Phys. Rev. Lett. 114, 108001 (2015).
  * Schoenholz  _et al._ [2016] S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, Nat. Phys. 12, 469 (2016).
  * Dyre [2006] J. C. Dyre, Rev. Mod. Phys. 78, 953 (2006).
  * Larini  _et al._ [2008] L. Larini, A. Ottochian, C. De Michele, and D. Leporini, Nat. Phys. 4, 42 (2008).
  * Pazmiño Betancourt  _et al._ [2015] B. A. Pazmiño Betancourt, P. Z. Hanakata, F. W. Starr, and J. F. Douglas, Proc. Natl. Acad. Sci. U. S. A. 112, 2966 (2015).
  * Buchenau and Zorn [1992] U. Buchenau and R. Zorn, Europhys. Lett. 18, 523 (1992).
  * Smessaert and Rottler [2014] A. Smessaert and J. Rottler, Soft Matter 10, 8533 (2014).
  * Mosayebi  _et al._ [2014] M. Mosayebi, P. Ilg, A. Widmer-Cooper, and E. Del Gado, Phys. Rev. Lett. 112, 105503 (2014).
  * Widmer-Cooper and Harrowell [2006a] A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett. 96, 185701 (2006a).
  * Kumar  _et al._ [2021] A. Kumar, I. Procaccia, and M. Singh, Europhys. Lett. 135, 66001 (2021).
  * Wagner  _et al._ [2011] H. Wagner, D. Bedorf, S. Küchemann, M. Schwabe, B. Zhang, W. Arnold, and K. Samwer, Nat. Mater. 10, 439 (2011).
  * Anderson  _et al._ [2013] D. Anderson, D. Schaar, H. G. E. Hentschel, J. Hay, P. Habdas, and E. R. Weeks, J. Chem. Phys. 138, 12A520 (2013).
  * Habdas  _et al._ [2004] P. Habdas, D. Schaar, A. C. Levitt, and E. R. Weeks, Europhys. Lett. 67, 477 (2004).
  * Puertas and Voigtmann [2014] A. M. Puertas and T. Voigtmann, Journal of Physics: Condensed Matter 26, 243101 (2014).
  * Furst and Squires [2018] E. M. Furst and T. M. Squires,  _Microrheology_, Vol. 1 (Oxford University Press, 2018).
  * Şenbil  _et al._ [2019] N. Şenbil, M. Gruber, C. Zhang, M. Fuchs, and F. Scheffold, Phys. Rev. Lett. 122, 108002 (2019).
  * Li  _et al._ [2020] B. Li, K. Lou, W. Kob, and S. Granick, Nature 587, 225 (2020).
  * van Doorn  _et al._ [2020] J. M. van Doorn, R. Higler, R. Wegh, R. Fokkink, A. Zaccone, J. Sprakel, and J. van der Gucht, Sci. Adv. 6, eaba6601 (2020).
  * Buttinoni  _et al._ [2017] I. Buttinoni, J. Cha, W.-H. Lin, S. Job, C. Daraio, and L. Isa, Proc. Natl. Acad. Sci. U. S. A. 114, 12150 (2017).
  * Kob and Andersen [1994] W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 (1994).
  * Brüning  _et al._ [2009] R. Brüning, D. A. St-Onge, S. Patterson, and W. Kob, J. Phys.: Condens Matter 21, 035117 (2009).
  * Plimpton [1995] S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Li  _et al._ [2019] Y.-W. Li, C. K. Mishra, Z.-Y. Sun, K. Zhao, T. G. Mason, R. Ganapathy, and M. Pica Ciamarra, Proc. Natl. Acad. Sci. U. S. A. 116, 22977 (2019).
  * [42] See Supplemental Material at http://… for additional information on the slow dynamics, plastic length distribution, role of damping parameter, force duration and force magnitude, plastic length at $T=0$, alternative measure of structural relaxation and density dependence of the plastic length.
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
  * Shiba  _et al._ [2016] H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, Phys. Rev. Lett. 117, 245701 (2016).
  * Vivek  _et al._ [2017] S. Vivek, C. P. Kelleher, P. M. Chaikin, and E. R. Weeks, Proc. Natl. Acad. Sci. U. S. A. 114, 1850 (2017).
  * Illing  _et al._ [2017] B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, Proc. Natl. Acad. Sci. U. S. A. 114, 1856 (2017).
  * Shiba  _et al._ [2019] H. Shiba, T. Kawasaki, and K. Kim, Phys. Rev. Lett. 123, 265501 (2019).
  * Adam and Gibbs [1965] G. Adam and J. H. Gibbs, J. Chem. Phys. 43, 139 (1965).
  * Kirkpatrick  _et al._ [1989] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev. A 40, 1045 (1989).
  * Stevenson  _et al._ [2006] J. D. Stevenson, J. Schmalian, and P. G. Wolynes, Nat. Phys. 2, 268 (2006).
  * Donati  _et al._ [1998] C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 80, 2338 (1998).
  * Biroli  _et al._ [2008] G. Biroli, J.-P. Bouchaud, A. Cavagna, T. S. Grigera, and P. Verrocchio, Nat. Phys. 4, 771 (2008).
  * Karmakar  _et al._ [2014] S. Karmakar, C. Dasgupta, and S. Sastry, Annu. Rev. Condens. Matter Phys. 5, 255 (2014).
  * Hima Nagamanasa  _et al._ [2015] K. Hima Nagamanasa, S. Gokhale, A. K. Sood, and R. Ganapathy, Nat. Phys. 11, 403 (2015).
  * Zheng  _et al._ [2021] Z. Zheng, R. Ni, Y. Wang, and Y. Han, Science Advances 7, eabd1958 (2021).
  * Coslovich  _et al._ [2019] D. Coslovich, A. Ninarello, and L. Berthier, SciPost Physics 7, 077 (2019).
  * Ji  _et al._ [2020] W. Ji, T. W. J. De Geus, M. Popović, E. Agoritsas, and M. Wyart, Phys. Rev. E 102, 062110 (2020).
  * Hall and Wolynes [1998] R. W. Hall and P. G. Wolynes, J. Chem. Phys. 86, 2943 (1998).
  * Simmons  _et al._ [2012] D. S. Simmons, M. T. Cicerone, Q. Zhong, M. Tyagi, and J. F. Douglas, Soft Matter 8, 11455 (2012).
  * Ponder and Richards [1987] J. W. Ponder and F. M. Richards, J. Comput. Chem. 8, 1016 (1987).
  * Sammut and Webb [2011] C. Sammut and G. I. Webb,  _Encyclopedia of Machine Learning_ (Springer-Verlag, Berlin, 2011).
  * Widmer-Cooper and Harrowell [2006b] A. Widmer-Cooper and P. Harrowell, J. Non-Cryst. Solids 352, 5098 (2006b).
  * Hocky  _et al._ [2014] G. M. Hocky, D. Coslovich, A. Ikeda, and D. R. Reichman, Phys. Rev. Lett. 113, 157801 (2014).
  * Perera and Harrowell [1999] D. N. Perera and P. Harrowell, J. Chem. Phys. 111, 5441 (1999).


# On the scaling properties of (2+1) directed polymers in the high temperature
limit

Victor Dotsenko${}^{\,a,b}$ and Boris Klumov${}^{\,c}$ aLPTMC, Sorbonne
Université, Paris, France  bLandau Institute for Theoretical Physics, Moscow,
Russia  cJoint Institute for High Temperatures, Russian Academy of Sciences,
Moscow, Russia

###### Abstract

In this paper in terms of the replica method we consider the high temperature
limit of (2+1) directed polymers in a random potential and propose an approach
which allows to compute the scaling exponent $\theta$ of the free energy
fluctuations as well as the left tail of its probability distribution
function. It is argued that $\theta=1/2$ which is different from the zero-
temperature numerical value which is close to $0.241$. This result implies
that unlike the $(1+1)$ system in the two-dimensional case the free energy
scaling exponent is non-universal being temperature dependent.

###### pacs:

05.20.-y 75.10.Nr 74.25.Qt 61.41.+e

##  I Introduction

At present the statistical properties of one-dimensional directed polymers as
well as the other systems belonging to the so called KPZ universality class
KPZ  are sufficiently well studied (for the reviews see e.g. Corwin ; Borodin
; Rev ). In contrast to that, not so much is achieved in the studies of the so
called (2+1) model of directed polymers which describes the fluctuations of an
elastic string directed along the time axes which passes through a random
medium in the three-dimensional space. Due to extensive numerical simulations
it is rather convincingly established that at the zero-temperature and in the
limit of large times $t$ the free energy fluctuations of such directed
polymers scale as $t^{\theta}$ with the scaling exponent $\theta\simeq 0.241$
numerics1 ; numerics2 ; numerics3 . The third digit here is still disputed but
there is unanimous agreement that $\theta$ is not equal to $1/4$ (as every
theoretician would hope). On the other hand, on the theoretical side, there is
a rigorous proof that in the high temperature limit the disorder potential
remains relevant Lacoin ; Berger-Lacoin ; Comets . In other words, there is
effective localization of the polymers trajectories around disorder defined
favorite ”corridors” so that (like in the (1+1) system) even at high
temperatures the polymers statistics do not reduce to simple diffusion (as it
takes place e.g. in the high temperature phase of the (3+1) modelBolthausen ).

In this paper in terms of the replica technique we propose an approximate
(mean field) method which in the high temperature limit allows to estimate the
replica partition function in the limit of large number of replicas. This in
turn, makes possible to derive the left tail asymptotics of the free energy
fluctuations distribution function. Assuming that this (unknown) distribution
function is defined by only one energy scale one eventually finds that the
typical value of the disorder defined free energy fluctuations scale as
$\sim\sqrt{u}\,t^{1/2}$ (where $u$ is the parameter which describes the
strength of the disorder) which means that the time scaling exponent
$\theta=1/2$. This value is remarkably different from the zero temperature
numerical results. If correct, this statement implies that unlike the $(1+1)$
system (where it is rigorously proved that $\theta=1/3$ both at $T\to 0$ and
at $T\to\infty$) in the two-dimensional case the free energy scaling exponent
is non-universal being temperature dependent. Note also that with $\theta=1/2$
for the mean square fluctuations of the polymer trajectory
$\boldsymbol{\phi}(\tau)\;(0\leq\tau\leq t)$ simple energy balance arguments
yields the scaling $\overline{\phi^{2}}\sim t^{\theta+1}=t^{3/2}$.

It is also worthing to note that comparing the features of the high
temperature limits in (1+1) and (2+1) systems one finds one more interesting
point. In (1+1) dimension there is the high temperature regime in which the
effective temperature parameter is rescaled with time as
$\beta\to\beta\,t^{-1/4}$ and intermediate exponents are established Albert-
Khanin-Quastel-1 ; Albert-Khanin-Quastel-2 . In particular, in the limit
$t\to\infty$ this scaling interpolates between the weak disorder ($\beta=0$)
and strong disorder regimes ($\beta>0$). In a sense, this result shares some
features of the problem considered in the present paper.

In Section II we define the model and describe the general ideas of the
present approach. In Section III we present the systematic mean-field method
(hopefully valid in the limit of large number of replicas) which eventually
can be reduced to the solution of two-dimensional one-particle non-local
(integral) non-linear  differential equation (23). In the high temperature
limit this equation can be (numerically) solved providing the value of the
exponent $\theta=1/2$ (Section IV). Section V is devoted to the brief
discussion of the obtained results.

##  II The model and the replicas approach

We consider the model of directed polymers defined in terms of an elastic
string described by the two-dimensional vector
$\boldsymbol{\phi}(\tau)\equiv\bigl{(}\phi_{x}(\tau),\,\phi_{y}(\tau)\bigr{)}$
directed along the $\tau$-axes within an interval $[0,t]$ which passes through
a random medium described by a random potential $V(\boldsymbol{\phi},\tau)$.
The energy of a given polymer’s trajectory $\boldsymbol{\phi}(\tau)$ is

$H[\boldsymbol{\phi}(\tau);V]=\int_{0}^{t}d\tau\biggl{\\{}\frac{1}{2}\bigl{[}\partial_{\tau}\boldsymbol{\phi}(\tau)\bigr{]}^{2}+V[\boldsymbol{\phi}(\tau),\tau]\biggr{\\}};$
(1)

Here the disorder potential $V[\boldsymbol{\phi},\tau]$ is supposed to be
Gaussian distributed with a zero mean $\overline{V(\boldsymbol{\phi},\tau)}=0$
and the correlation function

$\overline{V(\boldsymbol{\phi},\tau)V(\boldsymbol{\phi}^{\prime},\tau^{\prime})}=u\,\delta(\tau-\tau^{\prime})U(\boldsymbol{\phi}-\boldsymbol{\phi}^{\prime})$
(2)

The parameter $u$ is the strength of the disorder and $U(\boldsymbol{\phi})$
is a smooth function characterized by the correlation length $\epsilon$. For
simplicity we take

$U(\boldsymbol{\phi})\;=\;\frac{1}{2\pi\,\epsilon^{2}}\;\exp\Bigl{\\{}-\frac{\boldsymbol{\phi}^{2}}{2\epsilon^{2}}\Bigr{\\}}$
(3)

One-dimensional, or the so called (1+1) version of this problem (when instead
of the vector we have a scalar field $\phi(\tau)$) with the
$\delta$-correlated random potential has been the focus of intense studies
during past three decades hhf_85 ; numer1 ; numer2 ; kardar_87 ; bouchaud-
orland ; hh_zhang_95 ; Johansson ; Prahofer-Spohn ; Ferrari-Spohn1 ; KPZ-TW1a
; KPZ-TW1b ; KPZ-TW1c ; KPZ-TW2 ; BA-TW2 ; BA-TW3 ; LeDoussal1 ; LeDoussal2 ;
goe ; end-point ; LeDoussal3 . At present it is well established that the
fluctuations of the free energy of this system are described by the Tracy-
Widom (TW) distribution TW-GUE  and their typical value scale with time as
$t^{1/3}$.

The general formulation of the considered (2+1) problem in terms of the
replica approach looks quite similar to the (1+1) one. For a given realization
of the random potential $V(\boldsymbol{\phi},\tau)$ the partition function of
the considered system (with fixed boundary conditions) is

$Z({\bf r},t)=\int_{\boldsymbol{\phi}(0)=\bf{0}}^{\boldsymbol{\phi}(t)={\bf
r}}{\cal D}\boldsymbol{\phi}(\tau)\exp\bigl{\\{}-\beta
H[\boldsymbol{\phi}(\tau),V]\bigr{\\}}\;=\;\exp\bigl{\\{}-\beta F({\bf
r},t)\bigr{\\}}$ (4)

where $\beta$ is the inverse temperature, $F({\bf r},t)$ is the free energy
which is a random quantity and the integration is taken over all trajectories
$\boldsymbol{\phi}(\tau)$ starting at $\bf{0}$ (at $\tau=0$) and ending at a
point ${\bf r}$ (at $\tau=t$). Note that this problem is equivalent to the KPZ
equation KPZ

$\partial_{t}F({\bf r},t)\;=\;\frac{1}{2\beta}\nabla^{2}F({\bf
r},t)\;-\;\frac{1}{2}\Bigl{(}\boldsymbol{\nabla}\,F({\bf
r},t)\Bigr{)}^{2}\;+\;V({\bf r},t)$ (5)

which describe the time evolution of the two-dimensional manifold $F({\bf
r},t)$ in a random potential $V({\bf r},t)$.

For simplicity, in what follows we are going to consider the problem with the
zero boundary conditions:
$\boldsymbol{\phi}(\tau=0)=\boldsymbol{\phi}(\tau=t)={\bf 0}$. The free energy
probability distribution function $P(F)$ of this system can be studied in
terms of the integer moments of the above partition function, eq.(4):

$\overline{Z^{N}}\equiv
Z(N,t)\;=\;\prod_{a=1}^{N}\int_{\boldsymbol{\phi}_{a}(0)=0}^{\boldsymbol{\phi}_{a}(t)=0}{\cal
D}\boldsymbol{\phi}_{a}(\tau)\;\overline{\Biggl{(}\exp\Bigl{\\{}-\beta\sum_{a=1}^{N}H[\boldsymbol{\phi}_{a}(\tau),V]\Bigr{\\}}\Biggr{)}}\;=\;\int_{-\infty}^{+\infty}dF\,P(F)\,\exp\bigl{\\{}-\beta
NF\bigr{\\}}$ (6)

where $\overline{(...)}$ denotes the averaging over the random potentials
$V[\boldsymbol{\phi},\tau]$ Performing this simple Gaussian averaging we get

$Z(N,t)\;=\;\prod_{a=1}^{N}\int_{\boldsymbol{\phi}_{a}(0)=0}^{\boldsymbol{\phi}_{a}(t)=0}{\cal
D}\boldsymbol{\phi}_{a}(\tau)\;\exp\Bigl{\\{}-\beta
H_{N}[\boldsymbol{\phi}_{1}(\tau),\,\boldsymbol{\phi}_{2}(\tau),\,...\,,\boldsymbol{\phi}_{N}(\tau)]\Bigr{\\}}$
(7)

where

$\beta
H_{N}[\boldsymbol{\phi}_{1}(\tau),\,\boldsymbol{\phi}_{2}(\tau),\,...\,,\boldsymbol{\phi}_{N}(\tau)]=\int_{0}^{t}d\tau\Biggl{[}\frac{1}{2}\beta\sum_{a=1}^{N}\Bigl{(}\partial_{\tau}\boldsymbol{\phi}_{a}(\tau)\Bigr{)}^{2}\;-\;\frac{1}{2}\beta^{2}\,u\,\sum_{a,b=1}^{N}U\bigl{(}\boldsymbol{\phi}_{a}(\tau)-\boldsymbol{\phi}_{b}(\tau)\bigr{)}\Biggr{]};$
(8)

is the replica Hamiltonian which describes $N$ elastic strings
$\bigl{\\{}\boldsymbol{\phi}_{1}(\tau),\,\boldsymbol{\phi}_{2}(\tau),\,...\,,\boldsymbol{\phi}_{N}(\tau)\bigr{\\}}$
with the attractive interactions
$U\bigl{(}\boldsymbol{\phi}_{a}-\boldsymbol{\phi}_{b}\bigr{)}$, eq.(3). To
compute the replica partition function $Z(N,t)$, eq.(7), one introduces the
function:

$\Psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N};\;t)\;=\;\prod_{a=1}^{N}\int_{\boldsymbol{\phi}_{a}(0)=\bf{0}}^{\boldsymbol{\phi}_{a}(t)={\bf
r}_{a}}{\cal D}\boldsymbol{\phi}_{a}(\tau)\;\exp\Bigl{\\{}-\beta
H_{N}[\boldsymbol{\phi}_{1}(\tau),\,\boldsymbol{\phi}_{2}(\tau),\,...\,,\boldsymbol{\phi}_{N}(\tau)]\Bigr{\\}}$
(9)

such that

$Z(N,t)\,=\,\Psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N};\,t)\big{|}_{{\bf r}_{a}=0}$ (10)

Here the spatial arguments of this function are $N$ two-dimensional vectors
$\\{{\bf r}_{a}\\}$. One can easily show that $\Psi({\bf r}_{1},\,{\bf
r}_{2},\,...\,{\bf r}_{N};\,t)$ is the wave function of N quantum bosons
defined by the imaginary time Schrödinger equation

$\beta\frac{\partial}{\partial t}\Psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N};\;t)\;=\;\frac{1}{2}\sum_{a=1}^{N}\,\Delta_{a}\Psi({\bf r}_{1},\,{\bf
r}_{2},\,...\,{\bf
r}_{N};\;t)\;+\;\frac{1}{2}\,\beta^{3}u\,\sum_{a,b=1}^{N}U({\bf r}_{a}-{\bf
r}_{b})\,\Psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N};\;t)$ (11)

where $\Delta_{a}$ is the two-dimensional Laplacian with respect to the
coordinate ${\bf r}_{a}$. The corresponding eigenvalue equation for the
eigenfunctions $\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})$, defined
by the relation

$\Psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N};\;t)\;=\;\psi({\bf
r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})\,\exp\bigl{\\{}-t\,E_{N}\bigr{\\}}$
(12)

reads:

$-2\beta\,E_{N}\,\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N})\;=\;\sum_{a=1}^{N}\,\Delta_{a}\psi({\bf r}_{1},\,{\bf
r}_{2},\,...\,{\bf r}_{N})\;+\;\beta^{3}u\,\sum_{a,b=1}^{N}U({\bf r}_{a}-{\bf
r}_{b})\,\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})$ (13)

It is at this stage that we are facing the crucial difference of the
considered problem with the corresponding (1+1) one. The general solution of
the one-dimensional counterpart of eq.(13) is given by the Bethe ansatz wave
function which is valid only for $U(x)\,=\,\delta(x)$ and which is based on
the exact two-particle wave functions ($N=2$) solution exhibiting finite value
energy $E_{N=2}$. It is this fundamental property of (1+1) problem which
eventually allows to derive the Tracy-Widom distribution for the free energy
fluctuation.

The situation in the (2+1) case, eq.(13), is much more complicated. First of
all, in two dimensions there exists no finite two-particle solution for
$U({\bf r})\,=\,\delta({\bf r})$. One can easily construct an approximate
ground state solution of the two-particle problem for the finite-size function
$U({\bf r})$, eq.(3), but then one finds that in the limit $\epsilon\to 0$
(when $U({\bf r})$ turns into the $\delta$-function) the ground state energy
of this solution $E_{N=2}\,\to\,-\infty$. In other words, in two dimensions
(unlike one-dimensional case) one can not consider the problem with
$\delta$-correlated random potential. We have to study the system with finite
size function $U({\bf r})$ and the value of its spatial size $\epsilon$ must
explicitly enter into the final results.

Second, one can easily demonstrate that in the two-dimensional case the
construction of the $N$-particle wave function à la Bethe ansatz structure
based on the approximate two-particle solution (for finite $\epsilon$) doesn’t
work. So that, unlike one-dimensional case, here even the ground state energy
$E_{N}$ as well as $N$-particle ground state wave function $\psi({\bf
r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})$ of eq.(13) are not known. All that
makes the perspective to find the exact solution of the (2+1) problem rather
doubtful.

Here we would like to propose somewhat different strategy which, at least in
the high-temperature limit, makes possible to estimate the replica partition
function $Z(N,t)$ at $N\gg 1$ which in turn allows to derive the scaling
exponent of the free energy fluctuations.

##  III Mean Field Approach

First, to simplify notations, let us eliminate the parameter $\epsilon$ of the
correlation function $U\bigl{(}{\bf r}\bigr{)}$. Redefining

${\bf r}\;\to\;\epsilon\,{\bf r}$ (14)

and

$E_{N}\;=\;\frac{1}{\epsilon^{2}}\,{\tilde{E}}_{N}\,,$ (15)

instead of eq.(13) we get

$-2\beta\,\tilde{E}_{N}\,\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N})\;=\;\sum_{a=1}^{N}\,\Delta_{a}\psi({\bf r}_{1},\,{\bf
r}_{2},\,...\,{\bf r}_{N})\;+\;\beta^{3}u\,\sum_{a,b=1}^{N}U_{0}({\bf
r}_{a}-{\bf r}_{b})\,\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})$ (16)

where

$U_{0}({\bf r})\;=\;\frac{1}{2\pi}\;\exp\biggl{\\{}-\frac{1}{2}\,{\bf
r}^{2}\biggr{\\}}$ (17)

It is evident that in a general case eq.(16) can not be solved. However in the
limit of large number of particles, $N\gg 1$, one hopefully can use the
standard trick of the mean field approximation, in which the $N$-particle wave
function $\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf r}_{N})$ factorizes into
the product of $N$ one-particle functions, namely

$\psi({\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N})\;\simeq\;\prod_{a=1}^{N}\psi({\bf r}_{a})$ (18)

Substituting this into eq.(16), we obtain

$-2\beta\,\tilde{E}_{N}\,\prod_{a=1}^{N}\psi({\bf
r}_{a})\;=\;\sum_{a=1}^{N}\,\Delta_{a}\psi({\bf
r}_{a})\;\prod_{b\not=a}^{N}\psi({\bf
r}_{b})\;+\;\beta^{3}u\,\sum_{a\not=b}^{N}U_{0}({\bf r}_{a}-{\bf
r}_{b})\,\psi({\bf r}_{a})\psi({\bf r}_{b})\;\prod_{c\not=a,b}^{N}\psi({\bf
r}_{c})\,+\,\frac{1}{2\pi}\beta^{3}uN\prod_{a=1}^{N}\psi({\bf r}_{a})$ (19)

Introducing notations

$\displaystyle\tilde{E}_{N}$ $\displaystyle=$
$\displaystyle-\frac{N}{4\beta}\,\lambda$ (20) $\displaystyle\beta^{3}u\,N$
$\displaystyle=$ $\displaystyle\frac{1}{2}\kappa$ (21)

and integrating eq.(19) over all ${\bf r}_{1},\,{\bf r}_{2},\,...\,{\bf
r}_{N}$ (taking into account that $\int d^{2}r\Delta\psi({\bf r})=0$) we get

$\kappa(N-1)C_{0}^{-2}\int d^{2}r_{1}\int d^{2}r_{2}U_{0}({\bf r}_{1}-{\bf
r}_{2})\,\psi({\bf r}_{1})\psi({\bf
r}_{2})\;=\;\lambda\,N\;-\;\frac{1}{2\pi}\kappa$ (22)

where $C_{0}\;=\;\int d^{2}r\,\psi({\bf r})$. Now, integrating eq.(19) over
${\bf r}_{2},\,...\,,{\bf r}_{N}$, using the above relation, eq.(22),
redefining $\psi({\bf r})\to C_{0}\,\psi({\bf r})$ and neglecting terms of
order $N^{-1}$ we get the following non-linear mean-field equation for the
one-particle function $\psi({\bf r})$:

$\Delta\psi({\bf r})\;-\;\lambda\,\psi({\bf r})\;+\;\kappa\,\psi({\bf
r})\,\int d^{2}r^{\prime}\,U_{0}({\bf r}-{\bf r}^{\prime})\,\psi({\bf
r}^{\prime})\;=\;0$ (23)

where

$\int d^{2}r\,\psi({\bf r})\;=\;1$ (24)

and the function $U_{0}({\bf r})$ is given in eq.(17).

Further strategy is the following. For given values of the parameters
$\lambda$ and $\kappa$, eq.(20), we have to find smooth non-negative solution
of eq.(23) such that $\psi({\bf r}\to\infty)\to 0$. Next, substituting this
solution into the constraint (24) we can find $\lambda$ as a function of
$\kappa$, which eventually gives us the dependence of the ground state energy,
eqs.(15) and (20), on the replica parameter $N$. First, let us demonstrate how
this strategy works in the well studied one-dimensional case.

###  III.1 The example of $(1+1)$ system

The one-dimensional version of eqs.(23)-(24) reads

$\psi^{\prime\prime}(x)\;-\;\lambda\,\psi(x)\;+\;\kappa\,\psi(x)\,\int
dx^{\prime}\,U_{1}(x-x^{\prime})\,\psi(x^{\prime})\;=\;0$ (25)

where

$\int_{-\infty}^{+\infty}dx\,\psi(x)\;=\;1$ (26)

and

$U_{1}(x)\;=\;\frac{1}{\sqrt{2\pi}}\;\exp\biggl{\\{}-\frac{1}{2}\,x^{2}\biggr{\\}}$
(27)

Redefining

$\psi(x)\;=\;\frac{\lambda}{\kappa}\,\phi\bigl{(}\sqrt{\lambda}\,x)$ (28)

and denoting $\sqrt{\lambda}\,x\;=\;z$, instead of eqs.(25)-(26) we get

$\phi^{\prime\prime}(z)\;-\;\phi(z)\;+\;\phi(z)\,\int
dz^{\prime}\,U_{\lambda}(z-z^{\prime})\,\phi(z^{\prime})\;=\;0$ (29)

where the function $\phi(z)$ satisfy the constraint

$\frac{\sqrt{\lambda}}{\kappa}\,\int_{-\infty}^{+\infty}dz\,\phi(z)\;=\;1$
(30)

and

$U_{\lambda}(z)\;=\;\frac{1}{\sqrt{2\pi\lambda}}\;\exp\biggl{\\{}-\frac{1}{2\lambda}\,x^{2}\biggr{\\}}$
(31)

According to eq.(30),

$\lambda(\kappa)\;=\;I^{-2}\,\kappa^{2}$ (32)

where

$I\;=\;\int_{-\infty}^{+\infty}dz\,\phi(z)$ (33)

According to eq.(32) in the high temperature limit both
$\kappa\propto\beta^{3}uN\to 0$ and $\lambda\to 0$. Thus, in this limit,
according to eq.(31),

$\lim_{\beta\to 0}\;U_{\lambda}(z)\;\to\;\delta(z)$ (34)

so that eq.(29) reduces to

$\phi^{\prime\prime}(z)\;-\;\phi(z)\;+\;\phi^{2}(z)\;=\;0$ (35)

One can easily check (numerically) that this equation has an instanton-like
solution with $\phi(0)\simeq 1.50\;,\phi^{\prime}(0)=0$ and
$\phi(z\to\infty)\to 0$ (see Fig.1)

Figure 1: Instanton solution of eq.(35)

Substituting this solution into eq.(33) we find

$I\;=\;\int_{-\infty}^{+\infty}dz\,\phi(z)\;\simeq\;6.00$ (36)

Thus, according to eqs.(20), (21), (32) and (36) we find

$E_{N}\;\simeq\;-\frac{1}{36}\beta^{5}\,u^{2}\,N^{3}\;\propto\;-N^{3}$ (37)

This result, except for the numerical prefactor, perfectly fits with the exact
value of ground state energy $-\frac{1}{24}\beta^{5}\,u^{2}\,N^{3}$ of the
one-dimensional $N$-particle boson system (see e.g. Rev ) and correspondingly
provide the well known value of the free energy scaling exponent $\theta=1/3$.

###  III.2 $(2+1)$ directed polymers

The situation in the $(2+1)$ case is more complicated. As there are no reasons
to expect that the ground state solution of the original Schrödinger equation
(13) is anisotropic, in what follows it will be assumed that the function
$\psi({\bf r})$ is radially symmetric: $\psi({\bf r})=\psi(|{\bf
r}|)\equiv\psi(r)$. In this case equations (23)-(24) take the form

$\psi^{\prime\prime}(r)\;+\;\frac{1}{r}\psi^{\prime}(r)\;-\;\lambda\psi(r)\;+\;\kappa\,\psi(r)\,\int
d^{2}r^{\prime}\,U_{0}(|{\bf r}-{\bf r}^{\prime}|)\,\psi(r^{\prime})\;=\;0$
(38) $2\pi\int_{0}^{\infty}dr\,r\,\psi(r)\;=\;1$ (39)

Redefining

$\psi(r)\;=\;\frac{\lambda}{\kappa}\,\phi\bigl{(}\sqrt{\lambda}\,r)$ (40)

and denoting $\sqrt{\lambda}\,{\bf r}\;=\;{\bf z}$, instead of eqs.(38)-(39)
we get

$\phi^{\prime\prime}(z)\;+\;\frac{1}{z}\phi^{\prime}(z)\;-\;\phi(z)\;+\;\phi(z)\,\int
d^{2}z^{\prime}\,U_{\lambda}(|{\bf z}-{\bf
z}^{\prime}|)\,\phi(z^{\prime})\;=\;0$ (41)
$2\pi\int_{0}^{+\infty}dz\,z\,\phi(z)\;=\;\kappa$ (42)

where

$U_{\lambda}(|{\bf
z}|)\;=\;\frac{1}{2\pi\lambda}\;\exp\biggl{\\{}-\frac{1}{2\lambda}\,|{\bf
z}|^{2}\biggr{\\}}$ (43)

The main difference with the one-dimensional case is that now in the limit
$\lambda\to 0$ the value of $\kappa$, eq.(42) remains finite. Indeed,
according to eq.(43),

$\lim_{\lambda\to 0}\,U_{\lambda}(|{\bf z}|)\;=\;\delta({\bf z})$ (44)

In this case eq.(41) reduces to

$\phi^{\prime\prime}(z)\;+\;\frac{1}{z}\phi^{\prime}(z)\;-\;\phi(z)\;+\;\phi^{2}(z)\;=\;0$
(45)

Figure 2: Instanton solution of eq.(45)

This equation has an instanton-like solution with $\phi(0)\simeq
2.39\;,\phi^{\prime}(0)=0$ and $\phi(z\to\infty)\to 0$ (see Fig.2).
Substituting this solution into eq.(42) we find that at $\lambda=0$

$\kappa(\lambda=0)\;\equiv\;\kappa_{0}\;\simeq\;31.00$ (46)

At non-zero $\lambda\ll 1$, for $\kappa>\kappa_{0}$ numerical solution of
eqs.(41)-(42) demonstrate perfect linear dependence (see Fig.3)

$\lambda(\kappa)\;=\;\gamma\,(\kappa-\kappa_{0})$ (47)

with

$\gamma\;\simeq\;0.050$ (48)

Figure 3: Dependence of $\lambda$ on $(\kappa-\kappa_{0})$, eq.(47)

On the other hand, in the region $\kappa<\kappa_{0}$ by numerical methods we
find no non-negative solutions of equations (41)-(42) (such that $\psi({\bf
r}\to\infty)\to 0$) for any value of $\lambda$.

##  IV Free energy scaling

Substituting eqs (15) and (20) into eq.(47) for the mean-field ground state
energy of considered two-dimensional $N$-particle boson system we find

$E_{N}\;\simeq\;-\frac{\gamma}{4\beta\epsilon^{2}}\,N\,\bigl{(}2\beta^{3}uN\;-\;\kappa_{0}\bigr{)}$
(49)

where $\epsilon$ is the correlation size of the random potential, eq.(3), and
the numerical (approximate) values of the factors $\gamma$ and $\kappa_{0}$
are given in eqs.(46) and (48). Note that the above result is valid only for

$N\;>\;N_{*}\,\equiv\,\frac{\kappa_{0}}{2\beta^{3}u}\;\gg\;1$ (50)

Thus, according to eqs.(49), (10) and(12), in the limit $t\to\infty$ for the
replica partition function we find the following estimate

$Z(N,t)\;\sim\;\exp\biggl{\\{}\frac{\gamma}{4\beta\epsilon^{2}}\,N\,\bigl{(}2\beta^{3}uN\;-\;\kappa_{0}\bigr{)}\,t\biggr{\\}}$
(51)

Correspondingly, we see that in the high-temperature limit and at large $N$
eq.(6) reads

$\int_{-\infty}^{+\infty}dF\,P(F)\,\exp\bigl{\\{}-\beta
NF\bigr{\\}}\;\sim\;\exp\biggl{\\{}\frac{\gamma u}{2\epsilon^{2}}\,(\beta
N)^{2}\,t\;-\;\frac{\gamma\kappa_{0}}{4\beta^{2}\epsilon^{2}}\,\beta
N\,t\biggr{\\}}$ (52)

As in the exponential in r.h.s of the above relation there is linear in $N$
term the total free energy $F$ can be redefined such that it splits into two
independent parts: $F\,=\,\overline{F}\,+\,\tilde{F}$, where
$\overline{F}=\gamma\kappa_{0}t/(4\beta^{2}\epsilon^{2})$ is an extensive non-
random (selfaveraging) part while $\tilde{F}$ is the fluctuating contribution
described by a distribution function $\tilde{P}\bigl{(}\tilde{F}\bigr{)}$
which according to eq.(52) is defined by the relation

$\int_{-\infty}^{+\infty}d\tilde{F}\,\tilde{P}\bigl{(}\tilde{F}\bigr{)}\,\exp\bigl{\\{}-\beta
N\tilde{F}\bigr{\\}}\;\sim\;\exp\biggl{\\{}\frac{\gamma
u}{2\epsilon^{2}}\,t\,(\beta N)^{2}\biggr{\\}}$ (53)

In the case the above relation would be valid for any $N$ we would find that
$\tilde{P}\bigl{(}\tilde{F}\bigr{)}$ is just simple Gaussian distribution
function. In fact, as eq.(53) is valid only for $N\,>\,N_{*}\,\gg\,1$,
eq.(50), it gives us only the left tail of this distribution:

$\tilde{P}\bigl{(}\tilde{F}\to-\infty\bigr{)}\;\sim\;\exp\biggl{\\{}-\frac{\epsilon^{2}}{2\gamma\,u\,t}\,\tilde{F}^{2}\biggr{\\}}$
(54)

Note that this asymptotics sets in at $|\tilde{F}|\gg(\beta N_{*})\,\gamma
ut/\epsilon^{2}\,\sim\,t/(\beta^{2}\epsilon^{2})$.

It would be natural to suppose that the entire (unknown) probability
distribution $\tilde{P}(\tilde{F})$ is a universal function in a sense that it
is defined by the only energy scale so that it can be represented in the form
$\tilde{P}(\tilde{F})=G\bigl{(}\tilde{F}/F_{0}(t))$, where $G(x)$ is some
unknown universal function and $F_{0}(t)$ is the typical energy scale of the
free energy fluctuations. It is supposed that this form holds for all scales
of $\tilde{F}$ including far tails where $|\tilde{F}|\gg F_{0}(t)$. If the
above hypothesis is correct then using the explicit expression for far left
tail of the distribution function $\tilde{P}\bigl{(}\tilde{F}\bigr{)}$,
eq.(54), one immediately finds that in the considered high-temperature limit
the typical value of the free energy fluctuations scale as

$\tilde{F}\;\sim\;F_{0}(t)\;=\;\frac{\sqrt{\gamma\,u}}{\epsilon}\,t^{1/2}$
(55)

where $\gamma\simeq 0.05$, eq.(48), and $u$ and $\epsilon$ are the strength
and the correlation length of the random potential, eqs.(2)-(3). The above
eqs.(54)-(55) constitute the main results of the present research.

##  V Conclusions

Regardless of their somewhat ”trivial” form, the results presented in this
paper, eqs. (54) and (55), imply several rather non-trivial conclusions.

First of all, it should be noted that the prefactor in the time scaling of the
free energy fluctuations, eq.(55), is defined by the parameters of the
disorder potential. Thus the estimate, eq.(55), is in a full agreement with
the rigorous proof Lacoin ; Berger-Lacoin ; Comets  which states that even in
the high temperature limit the statistical properties of (2+1) directed
polymers are defined by the random potential.

Second, the fact that in the high temperature limit the free energy time
scaling exponent $\theta=1/2$ is different from the one at the zero
temperature (which is close to $0.241$) means that in two dimensions this
scaling exponent is temperature dependent. Maybe this takes place because the
dimension $d=2$ is critical for ($d$ \+ 1) directed polymers.

Third, it is also interesting to note that in the high-temperature limit the
prefactor in time scaling of the free energy fluctuations, eq.(55), turns out
to be temperature independent (unlike the one-dimensional case, where it is
proportional to $\beta^{2/3}$).

It should be stressed however, that all the results presented in this paper
are based on two crucial assumptions. The first one is pure heuristic mean-
field ansatz, eq.(18). It looks quite reasonable in the limit of large number
of replicas (which according to eq.(50) corresponds to the high-temperature
limit). Moreover, it works very well in the one-dimensional case (see Section
III.A). The second one is the hypothesis that in the considered system the
entire probability distribution function of the free energy fluctuations
$\tilde{P}(\tilde{F})$ reduces to a universal function:
$\tilde{P}(\tilde{F})=G\bigl{(}\tilde{F}/F_{0}(t))$ and this is valid for all
energy scales including far tails where $|\tilde{F}|\gg F_{0}(t)$. For the
moment the only support for this hypothesis is that usually in physical
systems exhibiting scaling phenomena it is correct. In particular, it is
certainly correct for (1+1) polymers where the exact solution provide us with
the universal function $G(x)$ which is the Tracy-Widom distribution function.
On the other hand, although both assumptions looks quite reasonable this does
not guarantee that they are correct… In view of that, further analytic studies
as well as numerical simulations of the considered system at finite (or high)
temperatures would be extremely helpful.

###### Acknowledgements.

VD is thankful to Quentin Berger for fruitful discussion of mathematical and
probabilistic aspects of (d+1) directed polymers problem. VD is also grateful
to Pierre Le Doussal, Jeremy Quastel, Lev Ioffe, Sergei Nechaev and Maxim
Dolgushev for numerous illuminating discussions.

## References

  * (1) Kardar M, Parisi G, Zhang Y-C, Phys. Rev. Lett. 56, 889 (1986)
  * (2) Corwin I, Random Matrices: Theory Appl. 1, 1130001 (2012)
  * (3) Borodin A, Corwin I and Ferrari P, Comm. Pure Appl. Math. 67, 1129–1214 (2014)
  * (4) V.Dotsenko Statistical properties of one-dimensional directed polymers in a random potential, arXiv:1703.04305 (2017)
  * (5) Halpin-Healy T, Phys. Rev. Lett. 109, 170602 (2012)
  * (6) Halpin-Healy T, Phys. Rev. E 109, 042118 (2013)
  * (7) Halpin-Healy T and Palasantzas G, EPL, 105, 50001 (2013)
  * (8) Lacoin H, Commun. Math. Phys. 294 471 (2020)
  * (9) Berger Q and Lacoin H, Ann. Inst. Henri Poincaré: Probab. Stat. 53 430 (2017)
  * (10) Comets F, Lecture notes from the 46th Probability Summer School, Saint-Flour, 2016, Springer International Publishing (2017)
  * (11) Bolthausen E, Commun. Math. Phys. 123 529 (1989)
  * (12) Alberts T, Khanin K and Quastel J, Phys. Rev. Lett. 105 090603 (2010)
  * (13) Alberts T, Khanin K and Quastel J, Ann. Probab. 43 1212 (2014)
  * (14) Huse D A , Henley C L and Fisher D S, Phys. Rev. Lett. 55, 2924 (1985)
  * (15) Huse D A and Henley C L, Phys. Rev. Lett. 54, 2708 (1985)
  * (16) Kardar M and Zhang Y-C, Phys. Rev. Lett. 58, 2087 (1987)
  * (17) Kardar M, Nucl. Phys. B 290, 582 (1987)
  * (18) Bouchaud J P and Orland H, J. Stat. Phys. 61, 877 (1990)
  * (19) Halpin-Healy T and Zhang Y-C, Phys. Rep. 254, 215 (1995)
  * (20) Johansson K, Comm. Math. Phys. 209, 437 (2000)
  * (21) Prahofer M and Spohn H, J. Stat. Phys. 108, 1071 (2002)
  * (22) Ferrari P L and Spohn H, Comm. Math. Phys. 265, 1 (2006)
  * (23) Sasamoto T and Spohn H, Phys. Rev. Lett. 104, 230602 (2010)
  * (24) Sasamoto T and Spohn H, Nucl. Phys. B834, 523 (2010)
  * (25) Sasamoto T and Spohn H, J. Stat. Phys. 140, 209 (2010)
  * (26) Amir G , Corwin I and Quastel J, Comm. Pure Appl. Math. 64, 466 (2011)
  * (27) Calabrese P, Le Doussal P and Rosso A, EPL, 90, 20002 (2010)
  * (28) Dotsenko V, EPL, 90, 20003 (2010)
  * (29) Dotsenko V, J.Stat.Mech. P07010 (2010)
  * (30) Calabrese P and Le Doussal P, Phys. Rev. Lett. 106, 250603 (2011)
  * (31) Dotsenko V, J. Stat. Mech. P11014 (2012)
  * (32) Dotsenko V, J. Stat. Mech., P02012 (2013)
  * (33) Gueudré T and Le Doussal P, EPL, 100, 26006 (2012)
  * (34) Tracy C A and Widom H, Commun.Math Phys., 159, 151 (1994)