Subcritical jump probability and anomalous order parameter autocorrelations

Critical phenomena are crucial for the description of physical world's structure at different scales. Their natural realization is usually related to the spontaneous breaking of symmetries of the Hamiltonian describing the dynamics of the considered physical system when a control parameter is varied. In systems at thermal equilibrium the spontaneous symmetry breaking (SSB) occurs at macroscopic scales, at which the relevant dynamical behaviour is described effectively through the associated free energy. The later depends on macroscopic quantities which are either extensive (order parameters), obtained through spatial sum of the microscopic degrees of freedom, or intensive (control parameters) determined by the system's environment. For thermal systems one of the control parameters is the temperature T. Symmetries of the microscopic Hamiltonian are transferred also to the free energy, expressed through its dependence on the order parameter(s). Then, the appearance of SSB is encoded in the modification of the free energy minima with varying T. In the unbroken phase (usually high temperature regime) the free energy possesses a single minimum at a particular value of the order parameter(s) and there is a critical temperature value T_c , below which a degenerate set of minima emerges inducing the spontaneous breaking of the free-energy symmetries. An important condition for the emergence of SSB is the lack of communication between the free-energy degenerate minima through the spontaneous dynamical evolution of the system. For a system of infinite size the free-energy minima are disconnected for T < T_c and the SSB is exact, since the equilibrated system has to choose a specific free-energy minimum [1]. The situation is more complicated for finite systems. In fact, the dynamics in a finite system allow the communication between the free-energy minima even for arbitrary small temperatures in the range 0 < T < T_c . In this sense the SSB in a finite system, despite the emergence of degenerate minima in the associated free energy, is never completed [2]. The system is not forced to choose a free energy minimum at its equilibrium. Nevertheless, as we will show in the present manuscript, there is a way to define rigorously SSB also for finite systems. To illustrate our proposal we will consider the magnetization properties of the Ising classical spin model in two $(D=2)$ and three $(D=3)$ spatial dimensions. The microscopic degrees of freedom are the binary spin variables $\sigma_{\mathbf{i}}=\pm 1$ where the index i defines their position in the respective lattice. The free energy $\mathcal{F}(\Omega,M,T)$ of the extended system with volume $\Omega=N^D$ depends on the total magnetization $M={\sum_{\mathbf{i}}} \sigma_{\mathbf{i}}$, which is the macroscopic extensive order parameter, and the temperature T which is the intensive control parameter. The free energy is symmetric under the transformation $M \to -M$ reflecting the $\mathbb{Z}_2{\text{-symmetry}}$ of the microscopic Hamiltonian with respect to $\sigma_{\mathbf{i}} \to -\sigma_{\mathbf{i}}$.

In the high temperature regime (symmetric phase) the free energy possesses a single minimum at M = 0 while for low temperatures (broken phase) two degenerate minima at $\pm M_0(T) \neq 0$ occur. To distinguish the case of a finite system from that of an infinite one, we will use the pseudo-critical temperature T_pc whenever we refer to a finite system, in contrast to the critical temperature T_c used exclusively for the infinite system. According to this notation the pseudocritical temperature T_pc is defined as the temperature below which the two degenerate free-energy minima first occur in the finite system. Certainly T_pc depends on the system's size N. Our first task is to quantify the "communication" between these two minima through the dynamical evolution of the finite system in the temperature range $T \lesssim T_{pc}(N)$. To achieve this we will employ the Metropolis algorithm for the description of the system's dynamics. For each temperature T, in the range mentioned previously, we generate an ensemble of magnetization time-series. For each of these time-series we construct symbolic sequences consisting of two symbols L and R using the map $M \to L$ when M < 0 and $M \to R$ when M > 0. Based on these symbolic sequences we can define a jump probability P_LR from the subspace with M < 0 to the subspace M > 0, as well as the reverse jump probability P_RL , for the considered ensemble. Clearly, for symmetry reasons it holds $P_{LR}=P_{RL}$. We obtain a general law describing the dependence of P_LR on T which is valid both for the 2D as well as the 3D case. These jump probabilities can be used to define a mean number of jumps $\langle n_{LR} \rangle=\langle n_{RL} \rangle$ in the ensemble which in turn provides a measure quantifying the communication between the two free-energy minima in the broken symmetry phase. An accurate definition of these quantities will be given in the next section.

In fact, when the time interval of the simulation $\mathcal{N}_{sim}$ (measured in lattice sweeps) is given, we can strictly define a temperature T_SSB (N) below which the mean number of jumps between the subspaces with different magnetization sign in the ensemble becomes less than ${1 \over 2}$, indicating practically the disconnection of the two free-energy minima with respect to the system's evolution up to the time $\mathcal{N}_{sim}$. We use the term SSB completion temperature for $T_{SSB}(N,\mathcal{N}_{sim})$. Searching for an additional indicator distinguishing the temperature $T_{SSB}(N,\mathcal{N}_{sim})$ we also calculate the autocorrelation functions of the magnetization time-series in the temperature range $T_{SSB}(N,\mathcal{N}_{sim}) \leq T \leq T_{pc}(N)$ and we observe an anomalous enhancement of correlations as $T \to T_{SSB}(N,\mathcal{N}_{sim})$. To explain this behaviour we use the Method of Critical Fluctuations (MCF), developed in [3], which has revealed the presence of type I intermittent dynamics [4,5] in the magnetization time-series of the 3D Ising model at $T=T_{pc}(N)$. With the MCF tool we reveal the emergence of on-off intermittency [6] as we approach $T_{SSB}(N,\mathcal{N}_{sim})$. This is certainly an interesting scenario which needs deeper understanding (work in progress).

Our findings are eventually relevant when searching for critical points in experimental data, particularly in mesoscopic systems [7,8]. In such a case, depending on N, the distance $T_{pc}(N)-T_{SSB}(N,\mathcal{N}_{sim})$ in temperature space can be significantly large and therefore even the term "critical point" could be misleading. Nevertheless, phenomenological characteristics discussed in the present work enable an identification of traces of criticality going beyond thermodynamic response.

The remaining part of the paper is organized as follows. In the next section we present the definition of the jump probability P_LR and its calculation as a function of temperature ($T \lesssim T_{pc}(N)$) for the 2D and 3D Ising models with nearest neighbour interactions. In the third section we calculate the corresponding autocorrelation functions of the magnetization time-series for temperatures $T \lesssim T_{pc}(N)$ focusing on their anomalous enhancement as T approaches the SSB completion temperature $T_{SSB}(N,\mathcal{N}_{sim})$. We also use the MCF to analyse and explain the observed anomalous behaviour. Finally in the last section we present our concluding remarks.

To illustrate our claims we will consider the magnetization dynamics of the ferromagnetic Ising model with nearest neighbour interactions in 2 and 3 dimensions. The lattice geometry is assumed to be the conventional one (square and cubic respectively). The Hamiltonian of the system is given as

$$\begin{equation} H=-J \displaystyle{\sum_{\langle \mathbf{i},\mathbf{j} \rangle}} \sigma_{\mathbf{i}} \sigma_{\mathbf{j}}, \quad\quad J > 0, \end{equation} \tag{ 1 }$$

where $\langle \mathbf{i},\mathbf{j} \rangle$ denotes summation over nearest neighbours. As usual, the coupling J can be absorbed into the inverse temperature $\beta=T^{-1}$. The system obeys periodic boundary conditions. The Metropolis algorithm is used in all simulations, taking into account autocorrelations in order to produce statistically independent spin configurations. Our calculations are performed using an ensemble of 100 magnetization time series each having the time length of $\mathcal{N}_{sim}=10^5$ sweeps/spin. We have checked the convergence of our results for specific values of the size N and the temperature T, in both the 2D as well as the 3D model, using several $\mathcal{N}_{sim}$ values ranging from $5 \cdot 10^4$ to $5 \cdot 10^5$ sweeps/spin.

We focus on the properties of the average magnetization density $m={\sum_{i} \sigma_i \over \Omega}$ at equilibrium. With $\Omega=N^D\,(D=2,~3)$ we denote the volume of the considered lattice in lattice constant units. Firstly, we investigate the magnetization distribution $\rho(m)$ for various temperature values. It is calculated through the ensemble of the 100 magnetization density time-series. To achieve this we divide in 100 bins the interval $[-1,1]$ and we determine the number of times each magnetization density time-series visits a bin. Then, we average over the ensemble of 100 time-series. At the end we normalize the obtained distribution to have integral equal to one. In fig. 1(a) we show $\rho(m)$ at $\beta^{-1}=2.3$ for a square lattice with N = 128. This temperature value is slightly higher than the pseudo-critical one for this system's size, notated as T_pc (128). In fact, T_pc (128) is above but close to the critical temperature $T_c =\frac{2}{\ln(1 +\sqrt{2})} \approx 2.269$ [9] of the infinite system. In our approach T_pc (N) is determined as the temperature value for which a plateau region around the single maximum at m = 0 of $\rho(m)$ is formed. It turns out that it is in general slightly higher than the corresponding pseudocritical temperature obtained in the literature through specific heat expansion [10]. In fig. 1(b) we show $\rho(m)$ for $\beta^{-1} = 2.27$ where the formation of two maxima in $\rho(m)$ is clearly seen, indicating the emergence of two free-energy minima at this temperature value for the N = 128 system. Notice that the minima of the free energy appear as maxima in the distribution $\rho(m)$ since $\rho(m) \propto e^{-\beta \mathcal{F}(N,m,T)}$. Due to this relation $\rho(m)$ displays the same $\mathbb{Z}_2$ symmetry as the free energy, namely $\rho(m)=\rho(-m)$. In the inset of this plot we provide an example of a part of one of the magnetization density time-series after equilibrium is reached. Notice the jumps from the region of negative m to the region of positive m and vice versa.

Zoom In Zoom Out Reset image size

As a next step, for fixed lattice size N, we vary the temperature in our simulations within a region just below T_pc (N). For every temperature value we form an ensemble of 100 magnetization density time-series with length of $\mathcal{N}_{sim}=10^5$ sweeps/spin each. For any such time-series we construct a symbolic sequence consisting of two symbols L and R by applying the mapping f given as

$$\begin{equation} f(m)=\begin{cases} L, \quad m < 0, \\ R, \quad m > 0. \end{cases} \end{equation} \tag{ 2 }$$

We end up with 100 symbol sequences of the type (LRLLLLLR...RRL) per temperature value each having the length $\mathcal{N}_{sim}$. For each of these symbol sequences we calculate $n_{LR}\,(n_{RL})$ which is the total number of the strings $LR\,(RL)$ in this sequence. Then, we define the jump probabilities P_LR and P_RL as

$$\begin{equation} P_{LR}=\displaystyle{\frac{\langle n_{LR} \rangle}{\mathcal{N}_{sim}}}, \quad P_{RL}=\displaystyle{\frac{\langle n_{RL} \rangle}{\mathcal{N}_{sim}}}, \end{equation} \tag{ 3 }$$

where $\langle n_{LR} \rangle\,(\langle n_{RL} \rangle)$ is the average number of the strings $LR\,(RL)$ in the ensemble. As discussed previously, in the considered temperature range, due to $\mathbb{Z}_2$ symmetry, it holds $P_{LR}=P_{RL}$. Thus, in the following we restrict our study on the jump probability P_RL without loss of generality. Our goal is to calculate P_RL at different subcritical temperatures close to the pseudocritical T_pc (N) and for different lattice sizes N to obtain P_RL (N, T). The results of our calculations are presented in figs. 2(a), (b). In fig. 2(a) we present the data collapse plot for $Q(T,N)=-{1 \over b(N)}\ln{\left(\frac{P_{RL}(T,N)}{P_{RL}(T_{pc}(N),N)}\right)}$ for the 2D Ising model and various sizes of the square lattice N = 64 (black open circles), N = 80 (red crosses), N = 96 (green up triangles), N = 112 (blue stars) and N = 128 (violet triangles down). The scaling factor b(N) is used to reveal the universal behaviour obeyed by the calculated jump probabilities. On the x-axis we set the temperature difference $T_{pc}(N)-T$ with T_pc (N) the pseudocritical temperature corresponding to each N-value, calculated with separate simulations and cross-checked with the literature values [10,11]. The dashed black line is the function $(T_{pc}(N)-T)^{3 \over 2}$ indicating the master curve for this data collapse plot. As shown in fig. 2(b), similar results for Q(T, N) are also obtained for the 3D Ising case using cubic lattices with N = 16 (black open circles), N = 20 (red crosses), N = 24 (green up triangles), N = 28 (blue stars) and N = 32 (violet triangles down).

Zoom In Zoom Out Reset image size

The universal behaviour in figs. 2(a) and (b) is captured by a modified exponential decay law given by the function

$$\begin{equation} P_{RL}(T,N)=P_{RL}(T_{pc}(N),N) e^{-b(N) (T_{pc}(N) - T)^{3 \over 2}}. \end{equation} \tag{ 4 }$$

From eq. (4) it is possible to obtain an estimate for the SSB completion temperature $T_{SSB}(N,\mathcal{N}_{sim})$. For given $\mathcal{N}_{sim}$ the condition $\langle n_{RL} \rangle < \frac{1}{2}$ is used to define the temperature at which the m < 0 and the m > 0 regions (or the associated free-energy minima) can be considered as disconnected. This leads to the determination of $T_{SSB}(N,\mathcal{N}_{sim})$ through the implicit relation

$$\begin{equation} P_{RL}(T_{SSB}(N,\mathcal{N}_{sim}),N) = \frac{1}{2 \mathcal{N}_{sim}} \end{equation} \tag{ 5 }$$

leading to

$$\begin{equation} T_{SSB}(N,\mathcal{N}_{sim})=T_{pc}(N)-\left[\frac{\ln(s(N,\mathcal{N}_{sim}))}{b(N)}\right]^{\frac{2}{3}} \end{equation} \tag{ 6 }$$

with $s(N,\mathcal{N}_{sim})=2 P_{RL}(T_{pc}(N),N) \cdot \mathcal{N}_{sim}$ which describes very well the results obtained through the numerical simulations (up to a $2\%$ deviation). When considering experimental data it is more useful to calculate the transition rate w_RL instead of P_RL . In the context of lattice simulations it is defined by the expression $w_{RL}=\frac{P_{RL}}{\mathcal{N}_{sim}}$. Then, the condition for SSB completion becomes $w_{RL}\cdot \mathcal{N}_{sim}^2 < \frac{1}{2}$. Of course in this case $\mathcal{N}_{sim}$ should be replaced by the experimental observation time. In fact, we can write $b(N)=\frac{\sigma(N)}{T_{pc}(N)^{3 \over 2}}$ with $\sigma(N)$ dimensionless. As shown in fig. 3 the coefficient $\sigma(N)$ has a power-law dependence on N in the form $\sigma(N) \sim N^{3 \over 2}$ for the 2D Ising and $\sigma(N) \sim N^{\frac{7}{3}}$ in the 3D Ising case.

Zoom In Zoom Out Reset image size

Combining eq. (6) with the results shown in fig. 3 we obtain that the size $\Delta T(N) = T_{pc}(N)-T_{SSB}(N,\mathcal{N}_{sim})$ of the SSB completion region in temperature space, for fixed $\mathcal{N}_{sim}$, depends on the lattice size N, shrinking with increasing N according to the power laws

$$\begin{equation} \Delta T_{2D} \sim N^{-1}, \quad \Delta T_{3D} \sim N^{-14/9}. \end{equation} \tag{ 7 }$$

Notice that the exponents appearing in eq. (7) are compatible with the value ${1 \over \nu}$, where ν is the correlation length critical exponent for the 2D and the 3D systems, respectively. In fact, the exponent ${1 \over \nu}$ is the finite-size scaling exponent determining the extension $T_{pc}(N)-T_c \sim N^{-{1 \over \nu}}$ (with $T_c=T_{pc}(\infty)$) of the critical region in finite systems [12–17]. Notice that $\nu=1$ for the 2D [9] and $\nu \approx {2 \over 3}$ for the 3D Ising model, respectively [18].

Our next goal is to explore in more detail the magnetization dynamics within the SSB completion region. Fur this purpose, in this part of the present work we will consider a single magnetization density time-series (instead of the ensemble of 100 time-series used in the previous section) of length $\mathcal{N}_{sim}=2 \cdot 10^6$ sweeps/spin. In a previous work we have shown, in the context of the 3D Ising model, that the magnetization dynamics at the pseudocritical point possesses characteristics of type I intermittency [3]. The term type I intermittency is used to characterize dynamics describing the transition —through tangent bifurcation— of a stable fixed point to an unstable one, tuning a control parameter. In 3D Ising model at the pseudocritical point the fixed point undergoing the transition is at m = 0 and the intermittent dynamics refers to the case when it becomes marginally unstable. Then, the trajectory of magnetization density is expected to spend long time in the immediate neighbourhood of m = 0. In dynamical systems terminology the time intervals, for which the trajectory remains close to the marginally unstable fixed point (here m = 0), are called laminar lengths and they provide a benchmark of intermittent dynamics. However, in lattice simulations there is always stochastic noise which in general prohibits a direct observation of this behaviour. Nevertheless, there is a typical signature of intermittent dynamics which is quite robust against noise. It is the distribution of the associated laminar lengths [4]. As shown in [3], for critical thermal systems this distribution attains the form

$$\begin{equation} S(\lambda) \sim \lambda^{-\frac{\delta + 1}{\delta}}, \end{equation} \tag{ 8 }$$

where δ is the isothermal critical exponent taking the values $\delta \approx 15 (5)$ for the 2D (3D) Ising model, respectively [9,19]. In the simulations the laminar lengths λ are found by counting the number of successive sweeps for which $\vert m(t) \vert < \epsilon$ with $\epsilon \ll 1$. The actual value of is not strictly defined. There is an entire interval of values which lead to a power-law form for the laminar length distribution originating from the magnetization density timeseries. Usually, within this interval the estimated exponent possess very small fluctuations (less than $5\%$). Nevertheless, for sufficient statistics an optimal can be determined through minimization of the residua with respect to power-law form.

Here, we argue that, for given size N and simulation time $\mathcal{N}_{sim}$, varying the temperature in the SSB completion region from T_pc (N) to $T_{SSB}(N,\mathcal{N}_{sim})$, a gradual transition from type I intermittent dynamics to on-off intermittency [6] takes place. In on-off intermittency the dynamics evolve in a higher dimensional chaotic phase space containing an invariant subspace, usually protected by some symmetry, which may act as an attractor. Then the system remains in this attractor for long time intervals which are interrupted by high amplitude abrupt bursts following a path along the unstable phase space manifolds induced by a stochastic force [6]. In the case of the Ising model this scenario is realized as follows. Below $T_{SSB}(N,\mathcal{N}_{sim})$ there are two symmetric free-energy minima for the magnetization density values $\pm \frac{M_0}{\Omega}$. These two minima provide the invariant subspace of the available phase space and they form a set which is invariant under $\mathbb{Z}_2$. In contrast to the case at $T=T_{pc}(N)$ when the effective dynamics can be described exclusively in the one dimensional magnetization density subspace [20], for $T=T_{SSB}(N,\mathcal{N}_{sim})$ the dynamics is not one-dimensional any more. Our proposal is to consider the energy density as an additional variable influencing the dynamics of the system. Thus, the effective dynamics become two-dimensional at $T=T_{SSB}$. An argument supporting our proposal is shown in figs. 4(a), (b) where we present the time-series of the magnetization (fig. 4(a)) and the energy (fig. 4(b)) densities for the 3D Ising model (lattice size N = 32) at the temperatures $T=4.515 \approx T_{pc}(32)$ (red lines) and $T=4.45 \approx T_{SSB}(32,2 \cdot 10^6)$ (blue lines). In fig. 4(b) we clearly observe that the energy density is almost constant for $T \approx T_{pc}(32)$ (red line) and therefore it is not expected to influence the magnetization density dynamics. On the contrary, the energy density for $T \approx T_{SSB}(32,2 \cdot 10^6)$ (blue line) possesses a characteristic peak, exactly in the time interval in which the magnetization density at the same temperature (blue line in fig. 4(a) jumps from the one part of the invariant set to the other. Thus, in this case one could relate the energy density fluctuations to the stochastic force appearing in the on-off intermittency process.

Zoom In Zoom Out Reset image size

After having presented the general idea it is worth to try to connect the effective dynamical behaviour in the SSB-zone with the microscopic spin dynamics. At $T=T_{pc}(N)$ the equilibrium state is dominated by microstates containing different clusters of aligned spins having practically all possible sizes. The total volume of clusters with positive orientation is exactly equal to that of negatively oriented clusters. As long as the temperature is within the SSB-zone there is a communication between the two minima of the free energy. As the temperature decreases, approaching $T=T_{SSB}(N,\mathcal{N}_{sim})$, the equally oriented clusters start to coalesce building larger clusters. It is achieved through the m = 0 channel which acts as the unstable manifold in the effective magnetization density dynamics. However, as the size of the clusters of aligned spins increases (for temperatures T < T_pc ) this channel becomes more and more narrow (statistically suppressed). As discussed previously, for the system to follow this channel a systematic enhancement of the energy density must take place (as shown in the peak of fig. 4(b), blue line). As a consequence the on-off intermittency becomes established. A clear signature of this transition is the modification of the laminar length distribution $S(\lambda)$ as $T_{SSB}(N,\mathcal{N}_{sim})$ is approached. As shown in [21,22] for noisy on-off intermittency the laminar length distribution becomes a power law with exponent $\gtrsim \frac{3}{2}$. We have calculated the laminar length distribution around m = 0 for the 3D Ising system at the temperatures $T_{pc}(N=32)=4.515$ and $T=4.45$ which lies very close to $T_{SSB}(N=32,\mathcal{N}_{sim}=2 \cdot 10^6)$ according to the relation in eq. (4). We pick out in fig. 5 the bulk of the distribution where the power-law behaviour is clearly seen and the saturation region close to $\lambda=1$ (notice that the time is measured in sweeps) as well as the exponential tail of the distribution, induced by the finite size N of the considered system, are avoided. The increase of the slope of the power-law decay from 1.2 to 1.5 approximately, as T_SSB is approached, can be clearly seen.

Zoom In Zoom Out Reset image size

Related to this transition from type I to on-off intermittency in the magnetization density dynamics —as the temperature is varied from T_pc (N) to $T_{SSB}(N,\mathcal{N}_{sim})$— is an anomalous behaviour of the magnetization autocorrelation function displayed in figs. 6(a), (b). We use the standard definition of the autocorrelation function,

$$\begin{equation} C(\tau)=\frac{1}{\Delta m^2}\displaystyle{\sum_{t=1}^{\mathcal{N}_{sim}-\tau}} \bar{m}(t + \tau) \bar{m}(t) \end{equation} \tag{ 9 }$$

with $\bar{m}(t)=m(t)- \langle m \rangle$, $\langle m \rangle$ the mean value of the magnetization density in the time interval $1 \leq t \leq \mathcal{N}_{sim}$ (in sweeps) and $\Delta m^2 = \langle m^2 \rangle - \langle m \rangle^2$ the associated variance.

Zoom In Zoom Out Reset image size

We clearly observe in both cases an increase of the autocorrelations as T approaches $T_{SSB}(N,\mathcal{N}_{sim})$. This anomalous enhancement of autocorrelations in the on-off intermittency is due to the fact that, although the mean magnetization in this case remains close to zero, the corresponding time-series contains mainly small oscillations around the minima of the free energy (maxima of the distribution $\rho(m)$) which in turn may lead to large waiting times (laminar lengths) in this region. Thus, in the on-off intermittency the laminar phase is displaced in the neighbourhood of the free energy minima while the chaotic bursts correspond to the jump processes between them. This is in contrast to the type I intermittent dynamics where the laminar part of the magnetization density time series is located around the marginal free energy maximum at m = 0. The rare jump processes are responsible for the maximization of the autocorrelations at the border of the SSB zone. Certainly, if the jumps were very frequent the autocorrelations would decrease since the corresponding time-series would resemble an oscillatory behaviour in the presence of strong noise.

The analysis presented in this work aims to illuminate aspects of criticality as it appears in mesoscopic, thermally equilibrated, physical systems possessing a finite size. All results have been obtained through simulations of the conventional 2D and 3D Ising model using the Metropolis algorithm. In particular, we have shown the existence of a temperature region $[T_{SSB}(N,\mathcal{N}_{sim}),T_{pc}(N)]$, depending on system's size N and the simulation time $\mathcal{N}_{sim}$, within which the spontaneous symmetry breaking in finite systems takes place. The upper limiting temperature T_pc (N) (pseudocritical) is related to the emergence of degenerate minima in the system's free energy while the lower limiting temperature $T_{SSB}(N,\mathcal{N}_{sim})$ is related to the suitably defined disconnection of these two minima with respect to the underlying order parameter dynamics. We called this temperature range as SSB completion region. A characteristic property of this region is the increase of its size $T_{pc}(N) - T_{SSB}(N,\mathcal{N}_{sim})$ with decreasing system's size N. As a consequence, in a mesoscopic system this distance can be large. For example, such a case could appear in the Lattice QCD simulations at zero chemical potential. There, the critical point is estimated through the formation of a non-zero $\langle \bar{q} q \rangle$ condensate (sigma condensate), which in fact takes place at T_SSB (N). However, since the calculations are performed in quite small lattices (systems) the corresponding pseudocritical temperature could be significantly higher. Such a scenario could explain inconsistencies occurring in the attempt to apply Lattice QCD results to the phenomenology of ultra-relativistic ion collisions, related to the experimental search for the QCD critical endpoint [7,8]. Even more, assuming that the sigma condensate is formed in the fireball created in the ion collision experiments, it certainly possesses a finite life time and therefore the arguments developed in the present work should be taken into account when analysing the recorded data.

Furthermore, we have demonstrated that the order parameter autocorrelations in time possess an anomalous enhancement as T approaches $T_{SSB}(N,\mathcal{N}_{sim})$. We have argued that this anomalous behaviour is related to a transition from type I intermittent dynamics ($T \approx T_{pc}(N)$) to on-off intermittency ($T \approx T_{SSB}(N,\mathcal{N}_{sim})$).

We expect that the qualitative aspects of our analysis are independent of the algorithm (Metropolis here) used in the simulations. However, it is not clear if the values of the exponents appearing in eqs. (4), (7) are universal. To answer this question a very demanding computational effort going beyond the aim of the present letter is required.

Data availability statement: All data that support the findings of this study are included within the article (and any supplementary files).