Fokker-Planck formalism approach to Kibble-Zurek scaling laws and non-equilibrium dynamics

TL;DR

This paper employs the Fokker-Planck formalism to analyze non-equilibrium dynamics and Kibble-Zurek scaling laws in phase transitions, validating theoretical predictions and extending applicability to microscopic models and ion Coulomb crystals.

Contribution

It introduces a Fokker-Planck based framework for deriving and validating Kibble-Zurek scaling laws in various dynamical regimes of phase transitions.

Findings

Fokker-Planck method confirms Kibble-Zurek predictions in overdamped and underdamped limits.

Framework predicts finite-size scaling functions.

Applied to ion Coulomb crystals, demonstrating broader applicability.

Abstract

We study the non-equilibrium dynamics of second-order phase transitions in a simplified Ginzburg-Landau model using the Fokker-Planck formalism. In particular, we focus on deriving the Kibble-Zurek scaling laws that dictate the dependence of spatial correlations on the quench rate. In the limiting cases of overdamped and underdamped dynamics, the Fokker-Planck method confirms the theoretical predictions of the Kibble-Zurek scaling theory. The developed framework is computationally efficient, enables the prediction of finite-size scaling functions and is applicable to microscopic models as well as their hydrodynamic approximations. We demonstrate this extended range of applicability by analyzing the non-equilibrium linear to zigzag structural phase transition in ion Coulomb crystals confined in a trap with periodic boundary conditions.

Equations

$$\mathcal{F}=\frac{1}{2}\int dx\left[h^2({\partial }_x\varphi {)}^2+V(\varphi )\right],$$

$$V(\varphi )=\frac{\epsilon }{2}{\varphi }^2+\frac{g}{4}{\varphi }^4.$$

$$\epsilon (t)={\epsilon }_0+\frac{t}{{\tau }_Q}\left({\epsilon }_1-{\epsilon }_0\right),0\le t\le {\tau }_Q.$$

$$\left(\frac{{\partial }^2}{\partial t^2}+\eta \frac{\partial }{\partial t}\right)\varphi (x,t)=h^2\frac{{\partial }^2}{\partial x^2}\varphi (x,t)-\frac{\delta V(\varphi )}{\delta \varphi }+\zeta (x,t)$$

$$⟨\zeta (x,t)⟩$$

$$⟨\zeta (x,t)\zeta (x^{\prime },t^{\prime })⟩$$

$$\left(\frac{{\partial }^2}{\partial t^2}+\eta \frac{\partial }{\partial t}\right)\varphi (x,t)=\left(h^2\frac{{\partial }^2}{\partial x^2}-\epsilon (t)\right)\varphi (x,t)+\zeta (x,t).$$

$$\varphi (x,t)=\sum _n{\phi }_n(t)e^{i{k}_{n}x},$$

$$\left(\frac{{\partial }^2}{\partial t^2}+\eta \frac{\partial }{\partial t}\right){\phi }_n(t)=\left(-k_n^2{h}^2-\epsilon (t)\right){\phi }_n(t)+{\zeta }_n(t),$$

$$⟨\mathcal{A}⟩=\prod {\int }_{-\infty }^{+\infty }d{\phi }_n\prod {\int }_{-\infty }^{+\infty }d{\dot{\phi }}_n\mathcal{A}\prod P_n(t,{\phi }_n,{\dot{\phi }}_n),$$

$${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }d{\phi }_{n}d{\dot{\phi }}_n{P}_n(t,{\phi }_n,{\dot{\phi }}_n)=1.$$

$$\frac{\partial P_n(t,{\phi }_n,\dot{{\phi }_n})}{\partial t}=[-\frac{\partial }{\partial {\phi }_n}{\dot{\phi }}_n+\frac{\partial }{\partial {\dot{\phi }}_n}(\eta {\dot{\phi }}_n+$$

$$+\left(h^2{k}_n^2+\epsilon (t)\right){\phi }_n)+\frac{\eta }{\beta }\frac{{\partial }^2}{\partial {\dot{\phi }}_n^2}]P_n(t,{\phi }_n,\dot{{\phi }_n}),$$

$$G(x_1,x_2,t)=⟨\varphi (x_1,t)\varphi (x_2,t)⟩-⟨\varphi (x_1,t)⟩⟨\varphi (x_2,t)⟩.$$

$${\xi }_L(t)=\frac{\sqrt{{\int }_0^{L/2}dx{x}^{2}G(x,t)}}{\sqrt{2{\int }_0^{L/2}dxG(x,t)}}.$$

$$g_L(t)=L\frac{{\int }_0^{L}dx⟨{\left({\partial }_x\varphi (x,t)\right)}^2⟩}{{\int }_0^{L}dx⟨{\left(\varphi (x,t)\right)}^2⟩}.$$

$$\eta \frac{\partial }{\partial t}\varphi (x,t)=\left(h^2\frac{{\partial }^2}{\partial x^2}-\epsilon (t)\right)\varphi (x,t)+\zeta (x,t).$$

$$\eta \frac{\partial }{\partial t}{\phi }_n(t)=\left(-h^2{k}_n^2-\epsilon (t)\right){\phi }_n(t)+{\zeta }_n(t),$$

$$\varphi (x,t)=\frac{{\phi }_0^R(t)}{\sqrt{N_c}}$$

$$+\sqrt{\frac{2}{N_c}}\sum _{n=1}^{(N_c-1)/2}{\phi }_n^R(t)\cos (k_{n}x)+{\phi }_n^I(t)\sin (k_{n}x),$$

$$\eta \frac{\partial P_{o,n}(t,{\phi }_n^{R,I})}{\partial t}=\frac{\partial }{\partial {\phi }_n^{R,I}}[\frac{1}{\beta }\frac{\partial }{\partial {\phi }_n^{R,I}}+$$

$$+\left(h^2{k}_n^2+\epsilon (t)\right){\phi }_n^{R,I}]P_{o,n}(t,{\phi }_n^{R,I}).$$

$${\partial }_t{P}_{o,n}^{\mathrm{th}}=0.$$

$$P_{o,n}^{\text{th}}({\phi }_n^{R,I})$$

$$f_n^{\mathrm{th}}$$

$$P_{o,n}(t)=\frac{f_n(t)}{\sqrt{\pi }}{e}^{-f_n^2(t){\phi }_n^2}.$$

$$\frac{\partial }{\partial t}{f}_n(t)$$

$$G(x_1,x_2,t)=\frac{1}{2N_c{f}_0^2(t)}+\sum _{n=1}^{N_c-1}\frac{\cos (k_n(x_1-x_2))}{N_c{f}_n^2(t)}.$$

$${\xi }_L(t)=\frac{L}{2\sqrt{6}}\sqrt{1+12f_0^2(t)\sum _{n=1}^{(N_c-1)/2}\frac{(-1{)}^n}{f_n^2(t)n^2{\pi }^2}}.$$

$$g_L(t)=L\frac{{\sum }_{n=1}^{(N_c-1)/2}\frac{k_n^2}{f_n^2(t)}}{\frac{1}{2f_0^2(t)}+{\sum }_{n=1}^{(N_c-1)/2}\frac{1}{f_n^2(t)}}.$$

$${\xi }_L^{\text{th}}(\epsilon >0)=\sqrt{\frac{h^2}{\epsilon }-\frac{hL}{2{\epsilon }^{1/2}\sinh ({\epsilon }^{1/2}L/2h)}},$$

$${\xi }_L^{\text{th}}({\epsilon }_c=0)$$

$$g_{L\to \infty }^{\mathrm{th}}(\epsilon >0)/L$$

Full text

I Introduction

Non-equilibrium dynamics involving critical phenomena, such as phase transitions, is an important area of statistical physics Landau80 . The physical phenomena that arise when traversing a symmetry breaking second-order phase transitions at finite rate are of particular interest. Specifically, symmetry breaking at finite rate promotes the formation of non-equilibrium excitations that can stabilize forming topological defects in a process known as the Kibble-Zurek (KZ) mechanism Kibble:76 . When the quench is performed at finite rates the symmetry is broken locally, and spatially separated regions can select different symmetry-broken states within the ground state manifold, which results in defects forming at spatial locations where phases of different symmetries meet. A major achievement of KZ theory is the prediction that the average number of defects exhibits a power-law dependence on the quench rate, whose scaling exponents are determined by the equilibrium critical exponents of the phase transition.

KZ mechanism has been studied in a number of experiments (see delCampo:14 for a recent review). While the standard KZ argument applies in spatially homogeneous systems, in some experimental systems, such as Bose-Einstein condensates Sadler:06 ; Weiler:08 ; Lamporesi:13 and ion Coulomb crystal Pyka:13 ; Ulm:13 inhomogeneities, as well as finite-size effects need to be accounted for. Thus measured scaling may not agree with the prediction of KZ scaling exponents in the thermodynamic limit. In such cases numerical simulations are particularly valuable tools for gaining insights into the non-equilibrium dynamics. Simulations of KZ experiments typically involve the numerical evaluation of many stochastic trajectories to allow for the calculation of an accurate estimate of any statistical quantity, including the expected density of defects. The tracking of stochastic trajectories of individual quench realizations, followed by averaging over the obtained ensemble is known as the Langevin approach of stochastic thermodynamics. In stochastic thermodynamics, there exists the different but equivalent approach of studying expectations of observables known as the Fokker-Planck approach Risken:84 . Fokker-Planck equations are deterministic partial differential equations specifying the time-evolution of the probability distribution of the configuration of the system that interacts with a Markovian heat bath. Thus the Fokker-Planck approach aims to solve the non-equilibrium dynamics problem at the ensemble rather than at the individual realization level, as is the case in the Langevin approach. The aim of the current paper is to apply the Fokker-Planck formalism to the KZ problem. We develop the Fokker-Planck method for the KZ problem and show that it can reproduce the known non-equilibrium scaling laws. The advantages of this approach include a computationally fast evaluation of the scaling laws and access to numerically exact probability distributions.

The paper is structured as follows. In Section II, we formulate the problem by introducing the Ginzburg-Landau model of phase transition, the equations of motion within the Langevin and Fokker-Planck formulation and the observables relevant in the context of KZ scenario. In Section III and IV, we solve the Fokker-Planck equations in, respectively, the overdamped regime and underdamped regime. In Section V, we apply the method to a non-equilibrium structural phase transition between linear and zigzag configurations in Coulomb crystals.

II Non-equilibrium dynamics in Ginzburg-Landau theory

Ginzburg-Landau (GL) theory provides a good model of second-order phase transitions. We consider a scalar one-dimensional order parameter $\phi(x,t)$. In the GL theory of second-order phase transitions, the free energy of the systems is given by

[TABLE]

where the Ginzburg-Landau potential $V(\phi)$ reads

[TABLE]

The constants $g$ and $h$ are parameters of the model that depend on the microscopic structure of the system. The parameter $\varepsilon$ quantifies the distance to the critical point of the phase transition, located at $\varepsilon=\varepsilon_{c}=0$. The model specified by Eqs. (1) and (2) is ubiquitous in physics as it describes a symmetry-breaking phase transition: the position of the minimum of $V(\phi)$ changes from being found at $\phi=0$ for $\varepsilon>0$ to two energetically equivalent choices at $\phi=\pm\sqrt{-\varepsilon/g}$ for $\varepsilon<0$.

The main purpose of the present article resides in analyzing the non-equilibrium dynamics resulting from the finite-rate symmetry breaking induced by an externally controlled time-dependent parameter $\varepsilon(t)$. We consider linear quenches in $\varepsilon(t)$ with functional dependence given by

[TABLE]

where $\varepsilon(0)=\varepsilon_{0}>0$ and $\varepsilon(\tau_{Q})=\varepsilon_{1}<0$, so that the systems is in the symmetric phase at the start of the quench protocol and in the symmetry broken phase at the end of the quench protocol. The rate at which the critical point is traversed is $d{\varepsilon}(t)/dt=(\varepsilon_{1}-\varepsilon_{0})/\tau_{Q}$, and thus, it is determined by the quench time $\tau_{Q}$ once $\varepsilon_{0}$ and $\varepsilon_{1}$ are fixed.

Langevin approach.— The dynamics of the one-dimensional order parameter $\phi(x,t)$ is described by the following general stochastic equation of motion Hohenberg:77 ; Chiara:10 ; Nigmatullin:16

[TABLE]

where $\eta$ is the friction parameter and $\zeta(x,t)$ the stochastic force, which fulfills

[TABLE]

where $\left\langle\textstyle{\ldots}\right\rangle$ denotes the ensemble average, $\beta\equiv(k_{B}T)^{-1}$, $T$ is the temperature and $k_{B}$ is the Boltzmann constant. For $\varepsilon\geq 0$, the field $\phi(x)$, close to the ground state, has small amplitude such that $|g\phi^{3}|\ll|\varepsilon\phi|$. Therefore, to a good approximation, the higher order terms of $\phi$ in $V(\phi)$ can be neglected resulting in $V(\phi)\approx\varepsilon\phi^{2}/2$, as depicted in Fig. 1. Despite of this apparently naive simplification, the linearized stochastic equations of motion still reproduce the dynamics of realistic models Moro:99 ; Chiara:10 ; Nigmatullin:16 . The linearized version of Eq. (4) reads

[TABLE]

It is convenient to express the field $\phi(x,t)$ in the Fourier space

[TABLE]

where $k_{n}=2\pi n/L$. For simplicity, we consider periodic boundary conditions, $\phi(0,t)=\phi(L,t)$, and a real field $\phi(x,t)$, which implies $\varphi_{n}(t)=\varphi_{-n}^{*}(t)$. Substituting Eq. (8) in Eq. (7) results in decoupled equations of motion for each mode. The equation for the $n$th mode reads

[TABLE]

with $\langle\zeta_{n}(t)\rangle=0$ and $\langle\zeta_{n}(t)\zeta_{m}(t^{\prime})\rangle=2\eta\beta^{-1}\delta_{nm}\delta(t-t^{\prime})$.

As usual within the framework of statistical mechanics, we are interested in the ensemble averages of observable macroscopic quantities (e.g. correlation length) that are functions of the microstates of the system. The ensemble averaged value $\left<\mathcal{A}[\phi(t),\dot{\phi}(t);t]\right>$ of a physical quantity $\mathcal{A}$ of interest at time $t$ may be obtained by averaging over a number of stochastic trajectories generated by the Langevin dynamics. Typically, the smaller the system, the more stochastic trajectories are necessary to obtain a reliable and meaningful average. Formally, the ensemble average can be obtained by integrating over the field distribution expressed in the Fourier space as

[TABLE]

where $P_{n}(t,\varphi_{n},\dot{\varphi}_{n})$ is the time-dependent probability distributions for $n$th mode in the phase space, i.e. it defines the probability of obtaining a value $\varphi_{n}$, and its velocity, $\dot{\varphi}_{n}$, at time $t$ for a mode with momentum $k_{n}=2\pi n/L$. The probability distribution is normalized according to

[TABLE]

In the Langevin approach, ensemble averaging involves evaluating approximations to these probabilities from the repeated solutions of the stochastic dynamical equations. The Fokker-Planck approach provides an analytic expression for $P_{n}(t,\varphi_{m},\dot{\varphi}_{n})$ as a solution to fully deterministic partial differential equations, as we explain in the following.

Fokker-Planck approach.— The Fokker-Planck formalism is a well-known approach to handle stochastic dynamics, which, in contrast to the Langevin approach focuses from the start on the probability distributions of the stochastic variables. The dynamical equations for the probability distributions are deterministic partial differential equations Risken:84 . The Fokker-Planck counterpart of Eq. (9) that specifies the dynamics of the $n$th mode, reads

[TABLE]

which is known as the Kramers equation Risken:84 . Hence the full probabilistic dynamics is acquired solving Eq. (II).

Quantities of interest.— In the spirit of KZ mechanism, we characterize the dynamics by means of the correlations induced in the system as it traverses the second-order phase transition at $\varepsilon_{c}=0$ at different rates $d\varepsilon(t)/dt\propto\tau_{Q}^{-1}$. To quantify such correlations, we introduce the usual two-point correlation function

[TABLE]

As a consequence of the periodic boundary conditions $\phi(0,t)=\phi(L,t)$, the two-point correlation function depends only in the distance $x_{1}-x_{2}$, i.e. $G(x_{1},x_{2},t)\equiv G(x_{1}-x_{2},t)$. The correlation length of the field $\phi(x,t)$ can be defined as

[TABLE]

We also attempt to quantify the defect density formed during the evolution, that is, the number of domains or regions per unit length within $\phi(x,t)$ with an equivalent choice of the broken symmetry. In the defect region the field interpolates rapidly but smoothly between the chosen configurations and thus in those regions the field has large spatial variations. For that reason the density of defects may be quantified by the gradient of the field Laguna:97 . Hence, to quantify such spatial variations, we introduce the density $g_{L}$ as

[TABLE]

The quantities $\xi_{L}(t)$ and $g_{L}(t)$ contain important non-equilibrium dynamical information, allowing us to test the emergence of universal behavior such as KZ scaling.

KZ scaling laws in the thermodynamic limit.— The equilibrium correlation length in the thermodynamic limit diverges as $\xi\propto\left|\varepsilon-\varepsilon_{c}\right|^{-\nu}$, where $\nu$ is the corresponding critical exponent. Additionally, a second-order phase transition is also characterized by a diverging relaxation time, $\tau\propto\left|\varepsilon-\varepsilon_{c}\right|^{-\nu z}$, where $z$ is the dynamical critical exponent Hohenberg:77 . If the fourth and higher order terms of $\phi$ in Eq. (2) are negligible then $\nu=1/2$ (mean field exponent), while depending on the dynamical regime, $z=2$ or $z=1$, for overdamped or underdamped dynamics Laguna:98 ; Chiara:10 . The former regime is found when $\left|\eta\partial_{t}\phi\right|\gg\left|\partial^{2}_{t^{2}}\phi\right|$, while the latter takes place in the opposite limit. To derive scaling laws, we resort to the KZ argument, which states that due to the diverging relaxation time near the critical point, there will be a freeze-out instance, $\hat{t}$, at which the system is not able any longer to adjust its correlation length to its equilibrium value Kibble:76 ; delCampo:14 during the quench. Accordingly, traversing the critical point at a finite rate $\tau_{Q}^{-1}$ provokes the formation of defects or excitations, whose typical size scales as $\hat{\xi}\sim\tau_{Q}^{\nu/(1+z\nu)}$, where $\hat{\xi}$ corresponds to the correlation length at the freeze-out instant. Therefore, the density of excitations will scale as $L^{d}/\hat{\xi}^{d}\sim\tau_{Q}^{-d\nu/(1+z\nu)}$, where $d$ is the dimension of the system. The KZ scaling laws can also be derived using rescaling transformations of the equations of motion Gor:16 , which does not rely on the physical arguments of transition between adiabatic and impulsive dynamics Kibble:76 ; delCampo:14 . In this paper, we show how these rescaling transformations can be applied to our system, elucidating a set of non-equilibrium scaling functions, where KZ scaling laws appear as a special case.

III Overdamped regime: Smoluchowski equation

The general equation of motion which governs the dynamics is given in Eqs. (7) and (II) for Langevin and Fokker-Planck formalism, respectively. However, when the term $\eta\partial_{t}\phi$ dominates, $|\eta\partial_{t}\phi|\gg|\partial^{2}_{t^{2}}\phi|$, the dynamics of $\phi(x,t)$ is overdamped or pure relaxational Hohenberg:77 . In the overdamped regime, the Langevin equation of motion reads

[TABLE]

The Fourier decomposition Eq. (8) results in decoupled equations for each of the normal modes. The dynamical equation for the $n$th mode is

[TABLE]

with $\langle\zeta_{n}(t)\rangle=0$ and $\langle\zeta_{n}(t)\zeta_{m}(t^{\prime})\rangle=2\eta\beta^{-1}\delta_{nm}\delta(t-t^{\prime})$. The Fokker-Planck equation in the overdamped regime is known as the Smoluchowski equation. Let $P_{o,n}(t,\varphi_{n})$ denote the probability distribution for the $n$th mode with $k_{n}=2\pi n/L$; the subscript $o$ emphasizes the overdamped nature of the dynamics. Since $\varphi_{n}(t)$ is in general complex, it is more convenient to express the field as

[TABLE]

where we have introduced a momentum cut-off $k_{c}$ which sets a maximum number of modes $(N_{c}-1)/2$, and $\varphi_{n}^{R}(t)\equiv\textrm{Re}(\varphi_{n}(t))$ and $\varphi_{n}^{I}(t)\equiv\textrm{Im}(\varphi_{n}(t))$ due to the condition $\varphi_{-n}(t)=\varphi_{n}^{*}(t)$. Note that, without loss of generality, $N_{c}$ is chosen to be odd. Then, the Smoluchowski equation for $P_{o,n}(t,\varphi_{n}^{R,I})$ reads Risken:84

[TABLE]

We will assume that the system is initially in thermal equilibrium at $\varepsilon_{0}$. In thermal equilibrium, a probability distribution, $P^{\rm{th}}_{o,n}$, must fulfill

[TABLE]

Substituting Eq. (20) in (III) and solving the resulting differential equation gives

[TABLE]

As expected, these probabilities correspond to the Boltzmann distribution at given temperature $1/\beta$ of the bath. Note that they exist only for $\varepsilon_{0}>\varepsilon_{c}=0$ as a consequence of the harmonic approximation $V(\phi)\approx\varepsilon\phi^{2}/2$. Having determined the initial state, we now solve Eq. (III) to find the time-dependent probability distribution $P_{o}(t)$. For that, we make use of a Gaussian Ansatz

[TABLE]

Substituting Eq. (22) into Eq. (III) gives

[TABLE]

with the initial condition determined by the thermal equilibrium, $f_{n}(0)=f_{n}^{\rm th}$. Thus the full knowledge of the probabilistic dynamics is captured in a set of uncoupled differential equations for the variance of each probability distribution.

The knowledge of the functional form of the probability distributions, allows us to explicitly calculate the quantities of interests, namely $\xi_{L}(t)$ and $g_{L}(t)$. Substituting Eq. (22) into the expression for the two-point correlation function given by Eq. (13), and simplifying the resulting expression gives

[TABLE]

The correlation length is then obtained using Eq. (14),

[TABLE]

Similarly, the expression for $g_{L}(t)$ evaluates to

[TABLE]

Thus, we have obtained analytic expressions for the correlation length and the number of defects as a function of time for non-equilibrium Kibble-Zurek dynamical scenario. In the analytical expressions, we can now look for physical meaning such as the presence of non-equilibrium scaling laws. We break down the discussion into three parts, each of which corresponds to a different quench rate regime: infinitely slow or isothermal quench (the case of thermal equilibrium), sudden quench and finite-rate quench.

Thermal equilibrium.— In the limit $\tau_{Q}\rightarrow\infty$, thermal equilibrium is achieved at any $\varepsilon(t)$. Although this scenario is just a limiting case of a finite-rate quench, it allows us to gather some interesting equilibrium properties which will be helpful later on. We stress that the derived expressions given in Eqs. (25) and (26) are valid for either equilibrium or non-equilibrium. Equilibrium properties are recovered simply by considering thermal probability distributions $P_{o,n}^{\rm th}$ at any $\varepsilon$, that is, $f_{n}^{\rm th}$, given in Eq. (III). In the limit $N_{c}\rightarrow\infty$, the correlation length at thermal equilibrium for finite $L$ reads

[TABLE]

while in the thermodynamic limit, $L,N_{c}\rightarrow\infty$, keeping the cut-off $k_{c}$ finite, is just $\xi^{\rm th}_{L\rightarrow\infty}(\varepsilon>0)=h\sqrt{1/\varepsilon}$. As we expected for Ginzburg-Landau theory, we obtain a critical exponent $\nu=1/2$ for the diverging correlation length at the critical point, i.e., $\xi\propto|\varepsilon-\varepsilon_{c}|^{-\nu}$. At $\varepsilon=\varepsilon_{c}=0$, the resulting expression is particularly simple for finite $L$,

[TABLE]

As one expects for a finite system, the correlation length can not exceed the system size, reaching its maximum at the critical point $\xi_{L}\lesssim L$. The previous result gives precisely its saturation value in the harmonic approximation of the Ginzburg-Landau model, as well as the scaling $\xi_{L}(\varepsilon_{c}=0)\propto L$ in agreement with finite-size scaling theory Fisher:72 ; Hohenberg:77 .

In a similar way, we can calculate the $g_{L}^{\textrm{th}}(\varepsilon)$ in the thermodynamic limit as

[TABLE]

and hence, as $\varepsilon\rightarrow\varepsilon_{c}=0$, it vanishes as

[TABLE]

revealing its critical exponent, which turns out to be $1/2$.

Sudden quench limit.— We briefly comment on the limit of sudden quenches, that is, when $\tau_{Q}\rightarrow 0$. In this case, as the system has no time to react to external perturbations, the corresponding properties of the system remain unchanged from its initial thermal state. Therefore, the results of sudden quenches are simply given by the thermal initial state at $\varepsilon_{0}$, i.e., $\xi_{L}^{\rm th}(\varepsilon_{0})$ and $g_{L}^{\rm th}(\varepsilon_{0})$, which for the former the expression is explicitly given in Eq. (27).

Finite-rate quenches.— We consider now the non-equilibrium dynamics for finite $\tau_{Q}$ in the overdamped regime, where KZ theory predicts scaling laws as a function of the quench rate. That is, we quench linearly the parameter $\varepsilon(t)$ in a time $\tau_{Q}$ according to Eq. (3), and solve the equations of motion (23), by numerical integration. Numerical solutions can be easily done by means of standard Runge-Kutta techniques. We emphasize that solving Eq. (23) immediately allows us to calculate precise average quantities, while the Langevin approach requires evaluating many realizations, which is far more costly from a computational point of view.

We set a momentum cut-off of $k_{N_{c}}=5\pi$, which leads to a maximum number of modes $N_{c}$ for a given length $L$. Initially, the system is prepared in thermal equilibrium at $\varepsilon_{0}=100$ and quenched in a time $\tau_{Q}$ towards $\varepsilon_{1}=-10$. The other parameters are set to $h=5$, $\eta=10$ and $\beta=1$. The results for three different system sizes at $\varepsilon(t)=0$, $L=10$, $20$ and $40$, are presented in Fig. 2, where $g_{L}$ and $1/\xi_{L}$ exhibit KZ scaling at intermediate quench rates. As the system size increases, the region of universal power-law scaling gets broader. For very fast quenches, $\tau_{Q}\rightarrow 0$, $1/\xi_{L}$ and $g_{L}$ saturate to their initial value, while for $\tau_{Q}\rightarrow\infty$, they tend to its value at thermal equilibrium, $g_{L}\rightarrow 0$ and $1/\xi_{L}\rightarrow 2\sqrt{6}/L$, as explained previously. For $L=40$ we perform a fit to obtain the power-law exponent which agrees well with the KZ prediction, $\tau_{Q}^{-\nu/(1+z\nu)}=\tau_{Q}^{-1/4}$ for the overdamped regime. Furthermore, we illustrate the finite-size scaling at intermediate quench rates, which predicts that

[TABLE]

where $y_{\xi}(x)$ and $y_{g}(x)$ are non-equilibrium scaling functions which fulfill $y(x\ll 1)\sim$ constant (where KZ scaling law emerges) and $y(x\gg 1)\sim x^{\nu/(1+z\nu)}$, for nearly adiabatic quenches Nigmatullin:16 . For that, we plot $L/\xi_{L}$ which is expected to follow a functional form $L\tau_{Q}^{-\nu/(1+z\nu)}y_{\xi}(\tau_{Q}L^{-1/\nu-z})$ or simply $x^{-\nu/(1+z\nu)}y_{\xi}(x)$ being $x=\tau_{Q}L^{-1/\nu-z}$. Thus, $L/\xi_{L}$ and $g_{L}L$ depend only on the scaling variable $\tau_{Q}L^{-1/\nu-z}$. The collapse of the data onto a single curve shown in Fig. 3 confirms this non-equilibrium scaling hypothesis.

Additionally, we demonstrate the universality of the phase transition dynamics by transforming physical quantities in such a way that any $\tau_{Q}$-dependence is removed from the equations of motion, following the theory developed in Gor:16 . This is achieved by performing transformations $x\rightarrow x\tau_{Q}^{-1/4}$ and $t\rightarrow(t-t_{c})\tau_{Q}^{-1/2}$, where $t_{c}=\tau_{Q}\frac{\varepsilon_{0}}{\varepsilon_{0}-\varepsilon_{1}}$ is the instance at which the critical point is crossed. In this rescaled frame the dynamics is universal and hence the functional dependence of $\xi_{L}\tau_{Q}^{-1/4}$ and $g_{L}\tau_{Q}^{1/4}$ on $(t-t_{c})\tau_{Q}^{-1/2}$ is expected to be the same irrespective of the value of $\tau_{Q}$. This universality during the whole evolution is demonstrated in Fig. 4, which shows the collapse of the results for three different values of $\tau_{Q}$ onto two curves, one for $\xi_{L}$ and one for $g_{L}$.

IV General and underdamped regime: Kramers equation

Let us recall that the Kramers equation Eq. (II) describes the general dynamical regime that includes both dissipative and inertial terms. In our particular case, the Fokker-Planck equation which describes the dynamics reads

[TABLE]

where $P_{n}(t,\varphi_{n}(t),\dot{\varphi}_{n}(t))$ is now a two-dimensional probability distribution at time $t$. The analysis of Eq. (IV) is more intricate, but nevertheless the procedure is similar to the one presented in Sec. III for the overdamped dynamics.

Thermal equilibrium states are obtained from $\partial_{t}P_{n}^{\textrm{th}}=0$, whose solution in terms of $\varphi_{n}^{R}\equiv{\rm Re}(\varphi_{n})$ and $\varphi_{n}^{I}\equiv{\rm Im}(\varphi_{n})$ reads

[TABLE]

where

[TABLE]

The time evolution of the probability distributions $P_{n}(t)$ is given by time-dependent coefficients $A_{n}(t)$, $B_{n}(t)$ and $C_{n}(t)$. In this way, three coupled differential equations per mode under the protocol $\varepsilon(t)$ determine the dynamics,

[TABLE]

The average quantities are obtained in the same way as in the overdamped regime, but with the time-dependent probability distributions also dependent on $\dot{\varphi}_{n}$. Indeed, we can define the probability distribution $Q_{n}(t,\varphi_{n}^{R,I})$ once the velocity dependence is integrated out,

[TABLE]

where $F_{n}(t)$ depends on the coefficients $A_{n}(t)$, $B_{n}(t)$ and $C_{n}(t)$ as

[TABLE]

This allows us to directly apply the same expressions as those derived for the overdamped regime. Eqs. (25) and (26) for correlation length $\xi_{L}(t)$ and density $g_{L}(t)$ can be directly applied by just replacing $f_{n}(t)$ with $F_{n}(t)$.

Thermal equilibrium and sudden quenches.— Clearly, thermal equilibrium does not depend on the considered dynamical regime. Therefore, the same thermal equilibrium probability distributions are retrieved from Eq. (IV) and we refer to Sec. III for the discussion on equilibrium features, as well as the opposite limit, $\tau_{Q}\rightarrow 0$, of sudden quenches.

Finite-rate quenches.— As in the case of overdamped dynamics, the KZ scaling laws are observed at finite quench rates. We solve numerically the equations of motion Eq. (38) to obtain the time-dependent probability distributions for different $\tau_{Q}$. Then, from $F_{n}(t)$ we calculate $\xi_{L}(t)$ and $g_{L}(t)$ using Eqs. (25) and (26), respectively. As in the overdamped regime, we set a momentum cut-off of $k_{c}=5\pi$. The initial thermal state at $\varepsilon_{0}=100$ is quenched in a time $\tau_{Q}$ towards $\varepsilon_{1}=-10$. Additionally, we set $h=5$ and $\beta=1$. To illustrate the KZ scaling in the underdamped regime, we select a small friction coefficient $\eta=0.1$. The results are presented in Fig. 5 for three different system sizes, $L=10$, $20$ and $40$, where the latter already exhibits a power-law scaling $\sim\tau_{Q}^{\alpha}$ for wide range of quench times. The performed fit gives an exponent $\alpha=-0.34(1)$ for $1/\xi_{L}$, in agreement with the predicted KZ scaling $\tau_{Q}^{-\nu/(1+z\nu)}=\tau_{Q}^{-1/3}$ since $\nu=1/2$ and $z=1$. A more pronounced deviation is found for $g_{L}$, with $\alpha=-0.36(1)$, which might be caused by finite-size effects. Moreover, in Fig. 6 the finite-size scaling at intermediate quench rates is verified. The data collapse onto a single curve corroborates the relations $1/\xi_{L}\sim\tau_{Q}^{-\nu/(1+z\nu)}y_{\xi}(\tau_{Q}L^{-1/\nu-z})$ and $g_{L}\sim\tau_{Q}^{-\nu/(1+z\nu)}y_{g}(\tau_{Q}L^{-1/\nu-z})$. Recall that $L/\xi_{L}$ and $Lg_{L}$ are expected to follow $x^{-\nu/(1+z\nu)}y_{\xi}(x)$ and $x^{-\nu/(1+z\nu)}y_{g}(x)$ where $x=\tau_{Q}L^{-1/\nu-z}$ is a scaling variable and $y(x)$ are non-equilibrium scaling functions.

Finally, we exemplify the universality of the dynamics in the underdamped regime. If $\left|\eta\partial_{t}\phi\right|\ll\left|\partial^{2}_{t^{2}}\phi\right|$, one removes the $\tau_{Q}$ dependence by performing the transformation $x\rightarrow x\tau_{Q}^{-1/3}$ and $t\rightarrow(t-t_{c})\tau_{Q}^{-1/3}$ Gor:16 . Note that the transformation is different to the overdamped case. The quantities $\xi_{L}\tau_{Q}^{-1/3}$ and $g_{L}\tau_{Q}^{1/3}$ are expected to collapse for different quench times $\tau_{Q}$ when plotted against the rescaled time $(t-t_{c})\tau_{Q}^{-1/3}$. This collapse is shown in Fig. 7, where we plot the results of the calculations of $\xi_{L}(t)$ and $g_{L}(t)$ for three different values of $\tau_{Q}$ in the rescaled coordinates.

V Coulomb crystals: linear to zigzag phase transition

The analysis in the previous sections has been done for phase transitions described by a one-dimensional Ginzburg-Landau field theory. However, the Fokker-Planck approach is valid beyond the Ginzburg-Landau theory; the knowledge of the quench function and the dispersion relation of the system is sufficient to predict the expected density of defects or any other statistical observable. We will illustrate this by applying our method to the problem of dynamic structural phase transition in Coulomb crystals Fishman:08 ; Chiara:10 .

Coulomb crystals are ordered structures that form when charged particles in a global confining potential are cooled below a critical temperature. An example of the physical realization of Coulomb crystals are ion crystals in Paul traps. Structural transitions in Coulomb crystal can be induced by varying the global confining potential Retzker:08 ; Partner:15 . The KZ mechanism of defect formation was studied numerically and experimentally using linear to zigzag non-equilibrium phase transition in ion traps Pyka:13 ; Ulm:13 . The analysis in references delCampo:10 ; Chiara:10 ; Nigmatullin:16 relied on the mapping of the linear to zigzag transition to a Ginzburg-Landau field theory model. In this section, we show how to use the methods developed in this paper to analyze the dynamic linear to zigzag transition without resorting to Ginzburg-Landau theory.

We consider $N$ charged particles moving in a periodic cell of size $L$. The periodic boundary conditions simplify the analysis since they result in a homogeneous Coulomb crystal. Moreover, periodic boundary condition can be realized with the existing technology of ring ion traps Wang:15 ; Li:16 . The potential energy of the $N$ particles reads

[TABLE]

where $\textbf{r}_{i}=(x_{j},z_{j})$ is the coordinate of the $j$th ion, $m$ is the mass of the ions, $\omega_{t}$ is the transverse trapping secular frequency and $Q^{2}\equiv e^{2}/4\pi\epsilon_{0}$. There exists a critical frequency value $\omega_{t}^{c}=\sqrt{7\zeta_{R}(3)/2}\omega_{0}$, where $\omega_{0}=\sqrt{Q^{2}/ma^{3}}$ and $\zeta_{R}(x)$ is the Riemann zeta function. For $\omega_{t}>\omega_{t}^{c}$ the lowest energy configuration is a linear chain and for $\omega_{t}<\omega_{t}^{c}$ the lowest energy configuration is a two-row zigzag chain.

Initially, the $N$ ions are in thermal equilibrium in a linear chain configuration. The transverse frequency $\omega_{t}$ is then quenched linearly in time through the critical point $\omega_{t}^{c}$, thereby inducing a transition from a linear to zigzag configuration at a rate proportional to $1/\tau_{Q}$

[TABLE]

where $\omega_{i}>\omega_{t}^{c}$ and $\omega_{f}<\omega_{t}^{c}$. Since the quench is performed at finite rate, the system is driven out of equilibrium and there is a non-zero probability of formation of a number of structural defects.

Langevin approach.— The expectation of any observable $\left\langle\mathcal{A}(t)\right\rangle$ (including the number of defects) can be evaluated by repeatedly solving the stochastic equations of motion that describe the dynamics of the systems, and then estimating the expectation from the obtained sample of trajectories. The dynamics of the system is determined by the following Langevin equations of motion

[TABLE]

There, $\eta$ is the friction coefficient and $\zeta_{j}^{x,z}$ is the stochastic force that satisfies the following statistical relations

[TABLE]

where $\alpha,\gamma\in\{x,z\}$ and $j,k\in\{1,2,...,N\}$.

We solve equations of motion, Eqs. (45)-(46), using the Langevin Impulse integration method with a timestep of $10$ ns Skeel:02 . We consider a system of $N=21$ ions, with inter-ion spacing in the linear configuration $a=10\ \mu$m, mass $m=172$ amu, which corresponds to $\textrm{Yb}^{+}$ ions, temperature $T=5$ mK and friction coefficient $\eta=1.5\times 10^{-21}$ kg s*-1* Pyka:13 . The initial and final transverse frequency is set to $\omega_{i}=2\pi\times 477.5$ kHz and $\omega_{f}=2\pi\times 159$ kHz. Note that $\omega_{0}=2\pi\times 143$ kHz and hence, the critical frequency is $\omega_{t}^{c}=2\pi\times 293.4$ kHz. To ensure that the system is initially in thermal equilibrium, the system is evolved under fixed trap parameters for $100\ \mu$s before starting the quench protocol. The quench time $\tau_{Q}$ is varied from $10\ \mu$s to $400\ \mu$s. For each value of $\tau_{Q}$, we perform $2000$ simulations in order to obtain accurate estimations of statistical observables such as two-point correlation function $G(x,t)$, correlation length $\xi(t)$, number of defects at the end of the quench $\left<N_{d}\right>$, and probability distributions of the transverse displacement $z(t)$, which due to translational symmetry does not depend on the ion position.

Fig. 8 shows an example of the results of the calculations for a selected quench time $\tau_{Q}=41.7\ \mu$s, namely the probability distribution of the transverse displacement and the two-point correlation function at the critical point. In Fig. 9, we plot the scaling of several quantities as a function of the quench time $\tau_{Q}$. These quantities are the correlation length $\xi$ and the averaged square displacement $\left<z^{2}\right>$ at $\omega_{t}(t_{c})=\omega_{t}^{c}$, and the number of defects $\left<N_{d}\right>$ at the end of the quench. We find that the scaling exponent is $\sim 1/3$, which is in agreement with the existing results in Refs. Chiara:10 ; Nigmatullin:16 , which predict this scaling by mapping the problem to GL theory and using the KZ relation $\tau_{Q}^{-\nu/(1+z\nu)}$ with $\nu=1/2$ and $z=1$. Furthermore, following the theory developed in Gor:16 , we can transform physical quantities to remove their dependence on the quench time $\tau_{Q}$, as explained and demonstrated in previous sections for GL theory. This entails the collapse of the correlation length $\xi$ into a single curve for different $\tau_{Q}$ values when $\xi\tau_{Q}^{-1/3}$ is plotted against the rescaled time $(t-t_{c})\tau_{Q}^{-1/3}$ with $\omega_{t}(t_{c})=\omega_{t}^{c}$, as shown in Fig. 10. Note that since $\nu=1/2$ and $z=1$, the used transformation to obtain the collapse is equivalent as GL theory in underdamped regime (see Sec. IV). In addition to the results of the Langevin dynamics simulations, Figs. 8, 9 and 10 include the results of the Fokker-Planck approach. The results show a good agreement with the Fokker-Planck description of the problem even when non-linear terms are neglected, as we explain in the following.

Fokker-Planck approach.— We now apply the Fokker-Planck approach to the problem of non-equilibrium quenches from the linear to zigzag configuration. At the start of the quench the system is in the symmetric linear phase. The equilibrium configuration of the ions is given by $\textbf{r}_{j}^{(0)}=(x_{j}^{(0)},0)$, where for convenience we take $x_{i}>x_{j}$ for $i>j$. Due to the periodic boundary conditions, the equilibrium inter-particle distance is constant i.e. $a=x_{j+1}^{(0)}-x_{j}^{(0)}$. The linearized equations of motions for small displacements around the equilibrium configurations $q_{j}=x_{j}-x_{j}^{(0)}$ and $z_{j}$ are obtained by Taylor expanding the potential. In the second-order Taylor expansion the axial and transverse motion decouple Chiara:10 . The equations of motion for the transverse displacements in this limit are

[TABLE]

where $\mathcal{K}_{i,j}\equiv\left.-\partial^{2}V/\partial x_{j}\partial x_{i}\right|_{x_{j}^{(0)}}$ is given by

[TABLE]

Eq. (49) describes the motion of coupled oscillators, which can be decoupled by rewriting it in terms of the normal modes. The relation between the transverse coordinate vector $\vec{z}$ and the normal mode vector $\vec{\Psi}$ can be written as

[TABLE]

where $k_{n}=2\pi n/Na$ and the sign $+$ ($-$) indicates the parity under $k_{n}\rightarrow-k_{n}$. Substituting Eq. (51) in (49) gives

[TABLE]

where $\omega_{n}(t)$ defines the frequency of the normal modes,

[TABLE]

and $\zeta_{n}^{\pm}(t)$ represents the stochastic force in the normal mode space, which again fulfills $\left<\zeta_{n}^{p}(t)\zeta_{m}^{q}(t^{\prime})\right>=2\eta/\beta\delta_{pq}\delta_{nm}\delta(t-t^{\prime})$. The Fokker-Planck equations corresponding to the Langevin Eq. (52) are

[TABLE]

Therefore, we have reduced the problem to the solution of $N$ deterministic Fokker-Planck equations that determine the mode population probabilities at a chosen time $t$. Using Eq. (54) and the expressions for $P_{n}$, following the same procedure as in Sec. IV, allows the determination of any statistical observable, as well for example, the probability distributions for $z_{j}$ at time $t$. Note that, from Eq. (51) and since $\Psi_{n}^{\pm}$ are statistically independent and Gaussian distributed, $P(t,z_{j})$ adopts also a Gaussian form and independent of $j$,

[TABLE]

with a time-dependent variance

[TABLE]

where $F_{n}(t)$ is obtained from Eq. (54) in the same way as explained in Sec. IV. In particular, we calculate the non-equilibrium correlation length $\xi$ and the mean square transverse displacement $\left<z^{2}\right>$ for the same set of parameters as was used previously in the Langevin approach. A comparison between the two approaches is shown in Fig. 8 and 9. We emphasize that the Fokker-Planck results which have been obtained under a simplified description of the realistic model, where non-linear terms and fluctuations in the longitudinal coordinates have been neglected, still reproduce essential features of the considered non-equilibrium scenario in a quantitative way.

VI Conclusions

We have studied the emergence of universal scaling laws in non-equilibrium second-order phase transitions using Fokker-Planck formalism. We verify that the developed approach reproduces Kibble-Zurek scaling laws in one dimensional Ginzburg-Landau model in overdamped and underdamped dynamical regimes. Additionally, we use this approach to obtain the universal finite-size scaling functions and demonstrate the universality of the dynamics.

There are several advantages of the developed method. It allows us to determine universal scaling laws in a efficient way in comparison to Langevin approach, where ensemble averages must be computed numerically. It provides analytic results that are easily amenable to further analysis. Moreover, it has an extended range of applicability - it can be used to derive insights into the non-equilibrium symmetry breaking phase transitions in overdamped, underdamped and intermediate dynamical regimes; finite as well as infinite systems and systems that are not described directly by the Ginzburg-Landau model. We have illustrated the power of the developed framework by analyzing the non-equilibrium linear to zigzag structural phase transition of an ion chain with periodic boundary conditions. We find an excellent agreement between the results obtained using the Fokker-Planck approach and the non-linear Langevin dynamics simulations.

One challenge for future theoretical work is to include the coupling between the normal modes of the system during quench protocol and find a way to predict the resulting corrections to the scaling laws. Another interesting direction of research would be apply this framework to investigation of scaling laws of other important quantities in stochastic thermodynamics such as entropy production and work done.

Acknowledgements.

This work is supported by an Alexander von Humboldt Professorship, by DFG through grant ME 3648/1-1 and the EU STREP project EQUAM. This work was performed on the computational resource bwUniCluster funded by the Ministry of Science, Research and the Arts Baden-Württemberg and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC. R. P. thanks P. Fernández-Acebal and A. Smirne for useful discussions.