A single-walker approach for studying quasi-nonergodic systems

The jump-walking Monte-Carlo algorithm is revisited and updated to study the equilibrium properties of systems exhibiting quasi-nonergodicity. It is designed for a single processing thread as opposed to currently predominant algorithms for large parallel processing systems. The updated algorithm is tested on the Ising model and applied to the lattice-gas model for sorption in aerogel at low temperatures, when dynamics of the system is critically slowed down. It is demonstrated that the updated jump-walking simulations are able to produce equilibrium isotherms which are typically hidden by the hysteresis effect characteristic of the standard single-flip simulations.

study low-temperature dynamics of lattice models of relatively large size (~10 5 particles) using just a single-core simulations as contrasted to massively parallel (> ∼ 10 3 cores 36, 37 ) simulations by REM or similar methods.
The structure of the paper is the following. The description of the original JW algorithm and its development are presented in Section Method. The performance of the improved JW algorithm is discussed for the Ising and lattice-gas models in Section Results and Discussion. Some technical details are shown in Supplementary Information (SI).

Method
The original JW algorithm. Consider a system described by a Hamiltonian H = H(Q) which is in equilibrium with the thermal bath at temperature T (measured in energy units) where Q denotes a state from the set {Q} of all possible states the system can be in (canonical ensemble). Assume that the target temperature T is low enough so that the system is quasi-nonergodic at T. The JW algorithm 18,19 explores the state space by means of the standard MC sampling at T augmented by jumps to the states which have been sampled at higher temperature T′ > T. The temperature T′ is chosen high enough to ensure that the system is ergodic at T′. The states, to which the jumps are attempted to, are sampled according to the equilibrium Boltzmann distribution (thus constructing the j-distribution), where Z(β′) stands for the partition function at the inverse temperature β′ = (T′) −1 . These additional jumps from states Q i (conventionally sampled at T) to states Q j (sampled at T′) allow for the ergodic sampling of the state space at temperature T if they preserve the detailed balance, j being the transition probability between states Q i and Q j . The transition probability,  → Q Q ( ) i j , can be split into sampling and acceptance probabilities, where P s (Q j , β′) is proportional to the Boltzman probability at T′. The detailed balance condition given by Eq. (2) imposes a constrain on the acceptance probability P acc , The standard Metropolis expression for P acc is recovered from Eq. (5) as T′ goes to infinity (β′ → 0), since the j-distribution is then uniformly random. On the other hand, if the state space is sampled at the same temperature as that of the j-distribution (β = β′), the detailed balance is maintained with the acceptance probability independent of the state parameters, since the j-distribution itself would already encompass the required transition probability.
Development of the algorithm. Although the original JW algorithm successfully tackled the quasi-nonergodicity problem, its implementation faced technical constraints at time of its creation. As mentioned in the Introduction, they were related to (i) inability to maintain the detailed balance, (ii) inefficiency of jumps between low-and high-temperature states with increasing temperature difference and (iii) low access rate of states in the j-distribution stored on hard drive. In this section, we discuss solutions to these problems and suggest an updated protocol for the JW algorithm, which addresses the mentioned issues. Namely, we show that (i) the detailed balance can be maintained exactly within the new protocol presented below, (ii) the standard temperature-sequence approach can be employed for the JW to increase the jump-acceptance probability and (iii) RAM can be used for accessing the j-distribution. We start with reviewing of the MC algorithms using "jumpy" protocols for dealing with quasi-nonergodicity which can be considered as developments of the original JW algorithm aiming to improve its performance and avoid its limitations and constrains. They can be broadly classified into two groups one based on multiple walkers (multiple replicas) and another one using one or two walkers. Multiple-replica methods with replica exchange are used to simulate a broad set of systems ranging from spin glasses to bio-systems 3, 17, 38, 39 while single-(two-)walker approaches are less popular and mainly used for modelling small atomic clusters [20][21][22][27][28][29] . Multiple-walkers algorithms such as multiple-stage JW 40 , REM 30,31 (parallel tempering 23,24,41 which originated from simulated tempering 42,43 ) and annealed swapping 26 run in parallel several walkers at temperatures covering the gap between high (ergodic) and low (target) temperature. This allows to increase the jump acceptance probability and preserve the detailed balance in REM and annealed swapping although demanding massively parallel computations 36,37 .
Smart-walking (SW) method 28 , in contrast to REM and similar to JW, uses only two walkers running at high and low temperatures. In contrast to JW, SW before making a jump to a high-temperature distribution performs energy minimisation for high-temperature state and then makes a jump. SW significantly improves the acceptance probability for jumps and permits effective exploration of the energy landscapes especially for protein systems although it does not maintain the detailed balance. This drawback of SW (lack of detailed balance) is addressed by smart-darting (SD) 29 algorithm which makes jumps only to neighbouring energy minima, although finding all relevant minima can be problematic for SD in complex landscapes. Cool walking 27 suggests a different method for achieving detailed balance in two-walker approach by means of statistical quench performed on configurations sampled by high-temperature walker.
Technically, the detailed-balance problem for the JW algorithm was tackled by two approaches. The first one 40 requires multiple parallel simulations at high temperature T′, as well as one run at the target temperature, T < T′. The j-distribution is sampled at runtime by randomly choosing one of the systems simulated at T′ and attempting a transition to its current state. Use of several systems for sampling at high temperature prevents possible correlations, which may appear if sampling from the same simulation is done before it loses memory of the previously sampled state. While this method does not use large amount of memory since the states are sampled at runtime, it is based, similarly to REM, on multiple simulations running at the same time and it was replaced by well established REM maintaining the detailed balance 25 .
The second approach 18 is to perform the initial sampling of the j-distribution at T′ prior to the main simulation and save the recorded states. During the main simulation at T, a random state is periodically chosen from the saved j-distribution. The state could be accepted or rejected according to the Metropolis rule (Eq. (5)). Due to insufficient amount of RAM, the set of states, {Q rec (T′)}, recorded at temperature T′, was kept on hard-drive storage, thus severely limiting access speed and manipulation capabilities. In order to maintain the detailed balance, one had to ensure that the same state from {Q rec (T′)} is not chosen twice 25,44 . Within this strategy, the probability of picking the same element in {Q rec (T′)} more than once, decreases as the ratio of the initially recorded and eventually required states increases. However, even if significantly more states were recorded than the simulation at T ultimately needed, this probability would remain finite and the detailed balance would be only approximately obeyed.
Below we suggest a new algorithm which resolves the detailed-balance problem for a single-walker JW method. The key point of our approach is in using RAM for j-distribution which permits easy manipulations with recorded states. In particular, the recorded states stored in RAM can be easily removed once chosen for a possible jump-transition, whether it occurred or not. Moreover, only necessary number of states from j-distribution for its unbiased sampling can be stored thus avoiding any storage overhead.
A further challenge faced by JW as well as REM and related methods arises due to the transition acceptance probability vanishing as the difference between temperatures T and T′ increases 25 . This hinders access to the states sampled at T′ and thus prevents efficient equilibration of the system. The standard solution to this problem for the JW algorithm 45 is to run the simulation sequentially at several temperatures, . The initial temperature T 1 , is set high enough so that the system is fully ergodic at T 1 (see SI 1). Any two consecutive temperatures T α and T α+1 are chosen in such a way that the Boltzmann distributions at these temperatures overlap, i.e. there exist states which the system can visit with non-vanishing probabilities at both temperatures. At each temperature T α , the target states Q j for the jump transitions are sampled from {Q rec (T α−1 )}, and simultaneously the new j-distribution is constructed by recording the states to {Q rec (T α )}. This way the system remains ergodic at each successive temperature down to T n . In case of REM, a typical solution to the equivalent problem is to employ multiple different temperature replicas running concurrently 46,47 .
There has been a lot of work done to optimise the set of selected temperatures for REM 48,49 as well as the JW 44,45 and other methods 50 . Typically, the temperatures are selected to follow geometric sequence 51 . If the free-energy landscapes are sufficiently similar across the temperature range (which is not necessary the case e.g. for lattice polymers 52 and Lennard-Jones clusters 12,25 ), the resultant overlaps between the subsets of the state space explored at different consecutive temperatures are comparable. This, in turn, leads to approximately the same acceptance probabilities for replica exchanges or JW jumps. However, in general, the free-energy landscape can be significantly different for different problems (e.g. for proteins 39 and structural glasses 53,54 ). Therefore, the temperature-sequence optimisation strategy can depend on a particular system and preliminary simulations are often necessary to examine the subsets of available states at each temperature and optimise {T α } 24 . The question of such optimisation is outside the scope of this work. Instead, we provide a general-purpose multi-stage protocol for the JW algorithm given an optimised set of temperatures {T α }. Using this protocol the JW MC simulation run on a single thread can ergodically explore the state space of the system exhibiting quasi-nonergodic behaviour below a certain characteristic temperature T c .
The protocol is as follows: (i) Run the standard single-flip MC simulation at temperature T 1 > T c . Record the current state of the system at equal intervals between the regular MC steps, so that R states are recorded into array {Q rec (T 1 )} until the total required number, S, of MC steps is executed. (ii) Decrease the temperature to T α , where T α is the next element in the optimised set of temperatures {T α }, and run the standard single-flip MC simulation again. Choose R step indexes at random from 1 to S and once each of them is reached in the simulation, attempt a transition to a state randomly picked from {Q rec (T α−1 )} with the acceptance probability P acc defined by Eq. (5). Whether the transition is accepted or not, remove the suggested state from {Q rec (T α−1 )}, so that it cannot be chosen again. In the meantime, continue to record the states as in (i) into the new array {Q rec (T α )}, which now contains an ergodic sample of the state space at T α . Proceed until S MC steps are executed. (iii) Set {Q rec (T α )} as the new {Q rec (T α−1 )} from which to sample the j-distribution. Repeat (ii) and (iii) until the final temperature T n in the list {T α } is reached.
The algorithm provides us with two free parameters: R and S. The constraint on precision naturally sets the required number of MC steps, S, to be completed. The interval between the attempted jumps, S/R, is determined by the desired frequency of the jumps, (R/S) × P acc , and limited by how many states R it is feasible to store in memory. For REM, it has been shown that given there were no computational constraints, equilibration accelerates as the replica-exchange frequency increases 55,56 . However, since the JW algorithm simulates only one copy of the system at any time, the jump transition, unlike the replica-exchange process in REM, requires to modify larger amounts of data associated with the state of the system. Computational cost of the transition is thus system dependent and could impose an upper bound on the optimal jump rate. Since the attempted states are removed from the storage independently of their acceptance, it is efficient to keep P acc as high as possible, minimising the number of stored states. Therefore, we are only free to adjust R in order to set the required jump rate.
For the purpose of clarity, the above analysis assumes that the average acceptance probability P acc is the same for each {T α }. However, in general, P acc depends on temperature and the constant jump-rate can be achieved by adjusting the number of stored states {R α } for each temperature T α . All simulation parameters used for testing the algorithm are given in SI 1.
Finally, it is straightforward to generalise the JW method for parameter-dependent Hamiltonians in a similar way to that used in the case of multidimensional REM 57 . For a system described by a Hamiltonian H = H(Q, {λ m }), where {λ m } is a set of parameters (e.g. chemical potential in the case of grand-canonical ensemble), an equivalent formalism leads to a generalised version of the expression for P acc , This allows us to sample the non-local transitions from a distribution constructed not only at a different temperature, but also based on Hamiltonian characterised by different values of parameters λ ′ m . An example of such analysis is presented below for the lattice-gas model.
To summarise, the updated JW algorithm described above maintains the detailed balance by storing the j-distribution in RAM and, given an optimised set of temperatures, is able to explore ergodic behaviour of the system for relatively low temperatures at least for some lattice models as demonstrated in Sec.

Results and Discussion
In order to demonstrate the updated JW MC method and further clarify its performance, the algorithm has been implemented for two systems: a standard 3d Ising model and lattice-gas model for condensation of fluid in porous media. Both models, their implementation details and simulation results are described in the following subsections.
Ising model. The Ising model is one of the simplest models exhibiting phase transitions, making it well suited to test the performance of the JW algorithm. Since the memory required to record the states of such system grows linearly with the system size N, the total memory needed to store the multiple recorded states has been prohibitively high and the original JW method 18,19 has never been implemented for the Ising or other spin models. Here, we demonstrate that the updated JW algorithm is capable of handling relatively large systems (N ~ 10 5 ) using only the memory of a regular contemporary personal computer and a single CPU core.
As a test example we study a 3d Ising model on a body centred cubic (bcc) lattice of N = 2 × L 3 ferromagnetically interacting spins (with periodic boundary conditions), where L is the number of unit cells along the edge of a cubic sample. The dimensionality and lattice type are convenient for comparison with the lattice-gas models considered below.
The Hamiltonian I  of the Ising model with the interaction strength between two neighbouring spins J > 0 (J = 1) in external field H (measured in units of J) is given by, where the first sum is taken over all the nearest-neighbour spin pairs 〈ij〉 and τ ∈ + − { 1, 1} i is the spin variable for a spin on site i.
In order to visualise the multi-dimensional state space of the spin system, it is convenient to study its projection onto two dimensions representing magnetisation . As MC sampling is running, the number of MC steps (or total time for kinetic MC 58 ) the system spends in a macrostate  Q M ( , ) I is recorded. In the ergodic regime, this quantity is proportional to the probability distribution for system to be in a certain macrostate and below we refer to it as distribution of visits. The probability to visit some macrostates can be very small and domains of the state space containing such states (far tails of the probability distribution) cannot be explored by the JW algorithm for finite value of parameter S 12 (see Section Development of the algorithm). The results for the distribution of visits are displayed on the M-I  plane using the following colour-scheme. On the white background of unexplored state space, the most to the least visited macrostates are coloured in yellow to dark blue, respectively. Figure 1a illustrates the results of the updated JW MC simulations on a 18 × 18 × 18 bcc lattice (N = 11664) for zero external field H = 0 and several temperatures. At sufficiently high temperatures, all transitions between microstates are approximately equally likely, and therefore mainly the degeneracy of the macrostates determines which part of the state space is visited within the simulation time (see the single-peak distribution of visits around zero magnetisation for β = 0.13 in the upper panel of Fig. 1a). As temperature decreases and occupation of the low-energy states becomes exponentially more probable (according to Eq. (1)), the subspace of explored states shifts to the lower energies (cf. distributions at β = 0.13 and β = 0.16). Moreover, with decreasing temperature, the distribution of visits becomes bi-modal (see e.g. the distribution for β = 0.16) reflecting the symmetry of the spin system with respect to the change of the spin orientation. This is due to the fact that the degeneracy of macrostates with approximately zero magnetisation decreases (with decreasing energy) more rapidly than of those with finite magnetisation. The two peaks in the distribution of visits become more pronounced with further decrease in temperature (cf. distributions at β = 0.17 and β = 0.19) and at T = 0 only two lowest energy states retain non-zero probability, each corresponding to all spins aligned either up or down. This picture can be described by using the language of the free-energy landscapes, referring to the dependence of the relative free energy For temperatures above some critical value, T c , the relative free energy f(M) has a single minimum at zero magnetisation (see yellow line marked by solid squares in the bottom panel of Fig. 1a). At ≈ T T c , two symmetric minima at finite values of M separated by a barrier at  M 0 appear (see red line marked by open squares). With decreasing temperature, the barrier between minima increases, the minima become deeper (see lines marked by the circles) and their positions tend to M = ±1 at zero temperature.
The picture described above cannot be seen by means of MC simulations governed only by the standard single-flip mechanism, because the transitions through the free-energy barrier become exponentially rare with decreasing temperature. Eventually, within the standard single-flip MC simulations, the system settles down in one of the free-energy minima shown in Fig. 1a, thus violating the equilibrium and breaking ergodicity. This does not happen for the JW MC simulations for which the multiple-flip jump transitions to the states recorded at higher temperatures are possible. These transitions enable exploration of all the states proportionally to the canonical distribution at the current temperature thus maintaining an ergodic sampling of the state space and all free-energy minima present at sufficiently low temperatures become achievable within a single JW simulation.  Fig. 1b). The negative external field breaks the symmetry and makes the minimum at M − deeper. The JW algorithm samples the states at both free-energy minima according to the canonical probability distribution. This means that only the energy difference at these minima is significant for their relative occupation and the barrier between the minima is irrelevant for the JW sampling. As the macrostates around M + become ever less probable with decreasing temperature, the jump-transitions from this minimum to the minimum at M − are favoured over the reverse ones (see the asymmetric bi-modal distributions of visits at β = 0.16 and β = 0.17 in the upper panel of Fig. 1b). Therefore, the amount of time that the system spends exploring the M + minimum gradually decreases with decreasing temperature (because the energy difference between minima increases) until it eventually vanishes completely for the precision set up within the JW algorithm (see the single-peaked distribution at β = 0.19 in Fig. 1b).
An example of equilibrium behaviour of the Ising model obtained by means of the JW MC simulations contrasted with the results of the single-flip MC is presented in Fig. 2. The standard way to investigate M(H) is to run a conventional single-flip MC simulations at fixed temperature gradually incrementing the external field from sufficiently negative (so that −  M 1) to high positive values (  M 1) and then back again. Under this protocol, a hysteresis loop is observed 59 (the areas between solid and dashed lines in Fig. 2 for β = 0.17, 0.19, 0.21). In this hysteresis loop (in some systems found experimentally as well), the magnetisation of the sample starts to align with the external field and then reverse in response to the change in H. Such a hysteresis loop is still observed even if the change of the field takes place adiabatically slowly. This implies that the magnetic material is not in equilibrium and is stuck in a state which it cannot leave on the single-flip MC (or experimental) time-scales with the thermal energy available, i.e. the system gets stuck in a metastable state.
The JW MC algorithm provides a possibility to avoid trapping in metastable states, i.e. to remove the hysteresis loop, by sampling the state space according to the canonical equilibrium distribution. Indeed, Fig. 2 (see the solid symbols corresponding to the JW data) shows a clear and expected phase transition at H = 0 below the critical temperature T c (for bcc lattice, . −  T 0 157371(1) c 1 60 ). Therefore, the JW MC simulations produce equilibrium isotherms, i.e. M(H), which are obfuscated by hysteresis phenomenon in the single-flip MC simulations or experiment. The inset in Fig. 2 magnifies the transition region. In this region for fixed temperature (below critical) and H close to zero, the JW MC simulations detect two minima of the free-energy, energetically favourable (solid symbols) and unfavourable (open symbols), which are also shown in lower panel of Fig. 1b. The probability to visit these minima follows the canonical equilibrium distribution and thus the accessibility of the energetically unfavourable minimum becomes exponentially small with increasing energy difference between the minima. Eventually, the energetically unfavourable minimum cannot be detected by the JW sampling with fixed precision. This occurs for relatively small deviations of the external field from zero, e.g. |H| ~ 10 −3 for β = 0.17. For the single-flip MC simulations, the height of the barrier between the free-energy minima matters and energetically unfavourable minima, in this case, act as traps, i.e. they become metastable states from which the system can escape only due to action of relatively large external fields, e.g. H ~ 10 −1 for β = 0.17. The behaviour of the Ising model described above is well established and served as a test for the JW MC simulations. As we can see from the results presented in Figs 1 and 2, the JW MC simulations reproduce all expected phenomena and, in addition, allow the magnetisation equilibrium isotherms, i.e. M(H), to be obtained for sufficiently low temperatures. Lattice-gas model. In order to test the applicability of the updated JW MC algorithm to relatively large and complex systems, we have implemented the method for a lattice-gas model to study fluid sorption in porous media. Condensation of gas in porous materials occurs differently from that in free space, e.g. it takes place at a lower pressure (or chemical potential) compared to the bulk saturation value 61,62 . While studies of fluids confined in single pores of simple geometry have clarified the mechanism for such shifted transitions [62][63][64][65][66] , the situation in real materials such as mesoporous glasses and silica gels, that consist of an interconnected network of pores of various shape and size, is still under active investigation [67][68][69][70][71][72][73][74] .
In this work, porous media is represented by two models: (i) a small lattice toy model for which exact numerical solution is available (see SI 2.2) and (ii) structural model of silica aerogel (see SI 2.3 for structural details of the model).
In the lattice-gas model 4, 75-77 , each out of N tot lattice site can be occupied by a fluid or a matrix particle, as described by the occupancy variables τ i and 1 − η i , respectively, which are equal to unity (zero) for occupied (unoccupied) sites. The matrix sites do not change their state, i.e. cannot be occupied by the fluid, and their concentration is quantified by porosity, tot . In the simulation setup, the porous material is assumed to be an open system of N = φN tot pore sites connected to a reservoir of fluid particles. The number of the fluid particles η τ = ∑ N i N i i f tot in the system and thus the energy  (see SI 2.1 for lattice-gas Hamiltonian) of the system can vary, but the chemical potential μ, temperature T and volume V (or equivalently N) of the system are fixed. For a given set of μ, T and N, the system can be in a set of states which form a grand canonical ensemble described by the following grand partition function,  6)) for such a system is given by the following expression, where the non-local multiple-flip transitions from a current microstate of the system simulated at β and μ to a state in the j-distribution previously sampled at β′ and μ′ are allowed. However, to assure clarity and simplicity of the illustrative examples, the results presented below refer to a value of μ fixed within a single JW simulation. The quantity of interest in our analysis is the mean equilibrium fluid density (order parameter), 〈ρ〉 = 〈N f 〉/N as a function of chemical potential and temperature. The angular brackets mean thermodynamic averaging over the states {τ} of the grand canonical ensemble, i.e.
The lattice-gas model is a discrete model and thus the density ρ takes N + 1 discrete values, The probability, ρ ∼ P Q ( ( , ))  , introduced in Eq. (13), to find the system in a certain macrostate ρ ∼ Q( , )  is given by  stands for the degeneracy of macrostate  ρ ∼ Q( , ), and thus the probability of the system to be in a state with a particular ρ is  ). Similarly to F(M) in the case of the Ising model, the grand-potential landscape can be defined as the ρ-dependence of the grand potential ρ β ρ , where Z μ = Z μ P(ρ). The shape of ρ Ω( ) can describe possible phases of the system. For example, if ρ Ω( ) has one minimum then the system is characterised by a single phase. However, if two minima appear in ρ Ω( ) this might be an indication of coexistence of two phases with different densities. In order to visualize the grand-potential landscape, similarly to f(M) for the Ising model, we plot ω ρ β ρ ρ ≡ Ω − Ω = − P ( ) ( ( ) ) ln( ( )) 78 (calling this quantity as grand-potential landscape), where Ω = −β −1 ln Z μ is the total grand-potential of the system (see the bottom panels in Fig. 3).
The state space of the lattice-gas model can be sampled by means of the JW MC algorithm in a similar way as that of the Ising model. The distribution of visits of different macrostates is displayed in the ρ- ∼ plane (see the upper panels in Fig. 3), where the same colour-scheme is used as in Fig. 1. Approximate high-and low-temperature Scientific RepoRts | 7: 2242 | DOI:10.1038/s41598-017-01704-5 boundaries of the state space are also computed (see SI 3 for details) and displayed for reference in the upper panels of Fig. 3 (the red (dashed) and green (solid) lines, respectively). Figure 3 presents the results of the JW MC simulations for sorption of fluid in a model sample of aerogel for two different values of μ and several temperatures. At low temperatures for both values of μ shown in Fig. 3, the grand-potential landscape has two minima, as can be inferred from the shape of the low-temperature boundary,  ρ ∼ ( ) min (the solid green line in the upper panels). The low-density minimum (e.g. at ρ .  0 15 for μ = −4.18) represents the fluid distribution in the system for which the sites occupied by fluid are concentrated only around and nearby the quenched matrix sites (see SI 2.4). The minimum at the higher values of ρ corresponds to the aerogel sample filled with fluid, and only the added surface layers left unoccupied (see SI 2.4). As temperature is decreased, the explored part of the state space shifts down until the system eventually settles in a minimum (see the sequences of distribution of visits corresponding to gradually decreasing temperature in the top panels of Fig. 3). For sufficiently low or high values of μ the system straightforwardly descents to the low-or high-density minimum, respectively (the descent to the low-density minimum is shown in Fig. 3a). However, at the intermediate values of μ (see Fig. 3b) the state space sampled by the JW MC simulation splits into several regions and thus the distribution of visits has several peaks (see e.g. the distribution at β = 0.91 in the top panel of Fig. 3b). In this case, the right peak in the distribution (see colour-map for β = 0.80 in Fig. 3b) becomes dominant with decreasing temperature and eventually the system settles down in the high-density minimum of the grand-potential. The single-flip MC simulations would not be able to follow such equilibration of the system for μ = −4.145. This is due to development of the large entropic barrier between two minima leading to significant slowing down of single-flip dynamics. Without the jump-transitions provided by the JW MC simulations, the region of the state space corresponding to the high-density peak could not be sufficiently visited and system would remain in the low-density state, i.e. in the metastable state. Figure 3b also shows the existence of the secondary split of the high-density peak in the distribution of visits into two peaks (see the two neighbouring bright spots on the right in the colour map for β = 0.91) corresponding to two minima in the grand-potential landscape at ρ . The JW MC simulations performed for fixed temperature and different values of chemical potential provide equilibrium isotherms, i.e. ρ(μ), for sorption in aerogel (see the data points labelled by circles, squares and diamonds in Fig. 4a). The solid symbols refer to the densities corresponding to the deepest minimum of ω(ρ) while the open symbols (see the inset in Fig. 4a) are related to other minima of ω(ρ) available for given values of μ and β. A typical evolution of the grand-potential landscape with chemical potential for relatively low temperature (β = 0.89) is shown in Fig. 4b. In two bottom panels (μ = −4.147 and −4.146), the low-density minimum (marked by solid square) is the deepest one and system mainly has the density corresponding to this minimum. One (for μ = −4.147) and two (for μ = −4.146) higher-density minima marked by open squares are weakly occupied in equilibrium. They are separated by a large barrier (not fully shown) from the low-density minimum and cannot be detected by the single-flip MC simulations due to large height of this barrier as compared to temperature. For μ = −4.145, the low-and intermediate-density minima become of approximately the same depth (they both marked by solid squares) and the system, in equilibrium, can be in one of them with approximately equal probability but, again, this cannot be detected by the single-flip MC simulations due to high barrier between the minima. In the relatively narrow range of chemical potential, −4.145 ⩽ μ ⩽ −4.144, the most likely fluid density is related to the middle minimum of ω(ρ) (see the inset in Fig. 4a) which becomes approximately of the same depth as the high-density minimum at µ − .  4 144 (both minima are marked by solid squares in Fig. 4b). For greater values of chemical potential, μ ⩾ −4.144, the system in thermal equilibrium is in the high-density state (see the top panel in Fig. 4b where the high-density minimum is marked by solid square and two other states with lower density by open squares).
In contrast to the JW MC simulations, both experimental studies 73,79 and single-flip MC simulations 80,81 are not able to access the equilibrium for sufficiently low temperatures due to presence of the hysteresis effect. Indeed, the results of the single-flip MC simulations for adsorption (solid lines in Fig. 4a) and desorption (dashed lines) for the same aerogel model exhibit hysteresis (area between solid and dashed lines) for β = 0.89 and β = 1.39, although both adsorption and desorption curves coincide producing the equilibrium isotherm (solid line marked by circles representing the JW MC data) for β = 0.71. Therefore, obtaining the equilibrium isotherms for sorption in aerogel and other porous materials is often problematic, and requires techniques such as mean-field scanning curves 68 . However, despite those efforts conclusive answers have not yet been obtained to even the key questions such as about existence of the discontinuous phase transition in aerogel samples of different porosity 73 . The JW algorithm is designed to avoid hysteresis and thus can be an invaluable tool in detecting phase transitions and resolving such problems. A detailed analysis of equilibrium sorption isotherms in aerogel samples of different porosity is outside the scope of this paper and will be presented elsewhere.

Conclusions
To conclude, we have revisited the JW algorithm 18 which has been originally designed to tackle quasi-nonergodicity problem, a known problem in dynamics of interacting particles. However, due to demanding (at the time of its creation) memory requirements the algorithm lost in competition with other approaches such as REM. Bearing in mind significant improvements in RAM available in contemporary computers we have updated the JW algorithm, tested it for the Ising model and demonstrated its performance for sorption in lattice-gas models (a toy model and aerogel model).
The efficiency of the JW algorithm stems from the fact that only a single non-local (multiple-flip) jump is needed to cross the free-energy barriers as opposed to multiple replica exchanges in the case of REM. Even though each particular replica within REM is allowed to cross energy barriers by configuration exchange with a replica at a slightly different temperature, REM overall simply swaps the two replicas, and computer resources continue to be used for exploring the same free-energy minima. Therefore, if the minimum is deep enough the algorithm based on REM can spend vast amount of time randomly diffusing in the temperature dimension while exploring the same local free-energy minimum. In contrast, as temperature decreases the JW algorithm naturally drives the system to escape from the energetically unfavourable local minima and therefore, the system is maintained in the equilibrium throughout the entire duration of the simulation.
The principal differences between the two methods, the JW and REM, suggest the types of problems that can be tackled by each approach. Since the REM replicas tend to explore the local minima for extended periods of time throughout the simulation, the method is well suited for tasks such as mapping of the entire free-energy landscape, weight factor estimation for multicanonical algorithms and, more generally, for investigating the properties of the whole state space of the system. The JW algorithm is more appropriate for studying equilibrium behaviour of the system. As temperature decreases, the system simulated by the JW algorithm naturally approaches the most favourable explored minimum, thus this method is well suited to search for the global free-energy minimum and investigate system's behaviour at the lowest temperatures. Since significantly less CPU resources is required by the JW algorithm as discussed above, it can either enable to study systems that have previously been prohibited by the lack of computing power, or to allow a more efficient investigation of the parameter space, provided that the analysis is mainly focused on the equilibrium behaviour of the system. In particular, the JW MC algorithm can be used on a single computing thread to investigate equilibrium properties of a system, especially near discontinuous phase transitions where the standard single-flip MC simulations are inefficient due to critical slowing down of dynamics and REM is too expensive in terms of computational resources.
Successful application of the updated JW algorithm to lattice-spin and lattice-gas models might promise its relevance to other models facing quasi-nonergodic behaviour, such as spin-facilitated models for glassy dynamics 54 and lattice models for proteins 52 . Another possible area for application of the updated JW method could be in studying the systems exhibiting weak ergodicity breaking 14 with non-Boltzmann limiting distributions 82 .
The description of the improved JW MC method was intentionally given as simple and general as possible. Therefore, despite certain advantages in its current form, the algorithm leaves much room for further optimisation and improvements. Moreover, due to key similarities with REM, numerous developments that studies of REM provided throughout the last two decades are readily available to implement for the updated JW algorithm (e.g. the temperature set {T α } and jump frequency optimisation as well as combinations with Wang-Landau, multicanonical or simulated tempering algorithms).