Mean first passage times reconstruct the slowest relaxations in potential energy landscapes of nanoclusters
Teruaki Okushima, Tomoaki Niiyama, Kensuke S. Ikeda, Yasushi Shimizu
aa r X i v : . [ phy s i c s . a t m - c l u s ] S e p Mean first passage times reconstruct the slowest relaxations in potential energylandscapes of nanoclusters
Teruaki Okushima, ∗ Tomoaki Niiyama, † Kensuke S. Ikeda, ‡ and Yasushi Shimizu § College of Engineering, Chubu University, Matsumoto-cho, Kasugai, Aichi 487-8501, Japan Graduate School of Natural Science and Technology,Kanazawa University, Kakuma-cho, Kanazawa, Ishikawa 920-1192, Japan College of Science and Engineering, Ritsumeikan University,Noji-higashi 1-1-1, Kusatsu, shiga 525-8577, Japan Department of Physics, Ritsumeikan University,Noji-higashi 1-1-1, Kusatsu, shiga 525-8577, Japan (Dated: September 25, 2019)Relaxation modes are the collective modes in which all probability deviations from equilibriumstates decay with the same relaxation rates. In contrast, a first passage time is the required timefor arriving for the first time from one state to another. In this paper, we discuss how and whythe slowest relaxation rates of relaxation modes are reconstructed from the first passage times.As an illustrative model, we use a continuous-time Markov state model of vacancy diffusion inKCl nanoclusters. Using this model, we reveal that all characteristics of the relaxations in KClnanoclusters come from the fact that they are hybrids of two kinetically different regions of the fastsurface and slow bulk diffusions. The origin of the different diffusivities turns out to come fromthe heterogeneity of the activation energies on the potential energy landscapes. We also developa stationary population method to compute the mean first passage times as mean times requiredfor pair annihilations of particle-hole pairs, which enables us to obtain the symmetric results ofrelaxation rates under the exchange of the sinks and the sources. With this symmetric method,we finally show why the slowest relaxation times can be reconstructed from the mean first passagetimes.
I. INTRODUCTION
Recently, the dynamics of complex systems, such as therelaxation of glass-forming materials [1–12], the kineticsof biomolecules [13–23], and diffusion in nanoclusters [24–30], were studied in a unified way for Markov state models[31–35]. The slowest relaxation modes of these systemsdescribe the bottleneck processes, and hence they are themost crucial, e.g., for understanding glass transitions andrapid formations of mixed crystals [36–43].The relaxation rates and modes are the eigenvalues andeigenvectors, respectively, of the transition rate matrixof a Markov state model. In general, physical quanti-ties are expressed in terms of the eigenvalues and eigen-vectors. The resulting expressions, called spectral repre-sentations, give useful formulas that enable us to eval-uate the physical quantities with use of the eigenvaluesand eigenvectors [44–52]. We can compute the relax-ation rates and modes of realistic, complicated Markovstate models using numerical diagonalization algorithms[53]. However, it is hard to understand why the eigenvec-tors are formed in the shapes of the numerical diagonal-ization results because the eigenvectors are quite high-dimensional and complicated. To extract the essence ofthe relaxation properties of Markov state models, there ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] have been many studies, such as lumping or renormal-izing Markov state models [35–39, 42, 43], and applica-tions of network algorithms, such as Dijkstra’s shortestpath algorithm [40]. Although there are many pioneer-ing works concerning this problem [44–52], to the best ofour knowledge, this problem has not yet been completelyclarified.As a more specific indicator of diffusive transport thanthe slowest relaxation rates, the first passage time iswidely studied mainly for analyzing the kinetics in com-plex networks [15, 16, 54–58]. The first passage timesfrom one state to another target state in a kinetic networkare the required times of stochastic realizations for travel-ing from the former to the latter state for the first time.The corresponding mean first passage time is given bythe mean value of the first passage times of the stochas-tic realizations.Intuitively, we may interpret the slowest relaxation of asystem as the process that transports the excess probabil-ity to the maximum probability states of the equilibriumdistribution along the unavoidable and slowest transportroutes in order to achieve the equilibrium distribution.Therefore, it may be possible to understand the forma-tion of the slowest relaxation mode by searching for thestates, where the first passage times to the maximumprobability states are the largest, and then by findingout the main routes connecting the former to the lat-ter states. To the best of our knowledge, however, therehave been no such studies that search for the slowest re-laxations with this idea. Instead, all pioneering works,e.g., Refs. [44–51], concern mainly how the mean firstpassage times are expressed with the relaxation modesvia renewal theorems. It should be noted that in this pa-per we study the inverse problem, i.e., how and why theslowest relaxation modes are reconstructed by the meanfirst passage times.As a realistic problem, we analyze a KCl nanoclustermodel having one vacancy [29, 30, 59]. The vacancy dif-fuses in the cluster and introduces the mixing of atoms inthe cluster. As for the pioneering works on vacancy diffu-sion, the equilibrium vacancy concentration [60] and therelaxation process using a stochastic process simplifiedby a uniform diffusion equation [61] have been studied.Nevertheless, there are no studies in which the surfaceeffects of nanoclusters on the relaxations of the vacancydiffusion are taken into consideration. The most sub-stantial reason that makes such approaches difficult isthat it is hard to estimate the transition rates betweenall adjacent states on the high-dimensional potential en-ergy surface from interatomic interactions. Fortunately,in Ref. [30] we have successfully enumerated all statesand all transition rates between adjacent states in nan-oclusters of various sizes, and we elucidated the specificproperties, such as migration energies of vacancies, aris-ing from the surfaces of nanoclusters. In this paper, weuse these transition data to construct the Markov statemodel of KCl nanoclusters and investigate the relation-ship between the slowest relaxation modes and the firstpassage times in the Markov state model equipped withthe cluster surfaces of characteristic transition regions.The purpose of this paper is twofold. One is to un-derstand the formations of the slowest relaxation modesin terms of the first passage times of the vacancy diffu-sion in KCl nanoclusters. The other is to elucidate thetheoretical basis for why such a mean first passage timeanalysis applies to the slowest relaxations.In Sec. II, we introduce a Markov state model, its re-laxation modes, and its mean first passage times in ageneral setting. For the mean first passage times, we de-velop a stationary population method that enables us tocompute the first passage times from the stationary pop-ulations of the Markov state models that connect sinkswith sources. In Sec. III, we introduce the interatomicinteraction of the KCl nanoclusters and then constructthe Markov state models of the vacancy diffusion model.In Sec. IV, we compute the slowest relaxation mode andthe mean first passage times of a KCl nanocluster. Wefind there that the shape of the energy landscape [32, 33]tells us why its relaxation makes effective use of the short-est routes of the vacancy from the center to the surface.In Sec. V, we study the second slowest relaxation modeand the corresponding mean first passage times of theKCl nanocluster. In Sec. VI, we first confirm that, un-der exchanging sinks and sources, the mean first passagetime approximation for the relaxation times in Sec. IVis asymmetric, while that in Sec. V is suitably symmet-ric. We then develop a symmetric stationary populationmethod for the mean first passage times, where they areinterpreted as the mean first encounter times of particle- hole pairs. The iterative use of the symmetric methodturns out to be equivalent to an inverse power methodfor diagonalization of matrices. We show that the meanfirst passage time approximations of the relaxation timesare good approximations that converge to the exact re-laxation times with the iterative use of the symmetricmethod. II. MARKOV STATE MODEL, RELAXATIONRATES, AND MEAN FIRST PASSAGE TIMES
In this section, we introduce a Markov state model, andwe describe a popular method of calculating first passagetimes for this model according to Refs. [50, 52]. We alsoshow that the mean first passage times obey stationarypopulation equations, which will be used to develop asymmetrized version of the population method later inSec. VI.
A. Continuous-time Markov state model
We start with a continuous-time Markov state modeldescribed by a transition rate matrix K with finite di-mension, n , of the state space. The kinetic equation isgiven by d P dt = K P , (1)where P is the probability distribution P =( p , . . . , p n ) T , with p i denoting the probability of the i thstate and the superscript T denoting the transpose. Weassume that K is time-independent and satisfies K ij > i = j ) and the probability conservation condition of P ni =1 ( K ) ij = 0 ( j = 1 , , . . . , n ). Further, we assumethat the equilibrium, lim t →∞ P ( t ), is a unique vector P eq satisfying the detailed balance conditions ( K ) i,j ( P eq ) j =( K ) j,i ( P eq ) i [50, 62]. Then, the eigenvalues of K satisfy0 = λ > λ > · · · > λ n − . (2)The equilibrium P eq coincides with the zeroth eigenvec-tor of K , and the first, second, . . . eigenvectors P i of K represent the slowest relaxation modes with the relax-ation times of | λ | − > | λ | − > . . . , respectively. B. Mean first passage times
The mean first passage times, t i,j , from a state j to i are evaluated by connecting perfect absorbers to all thefinal destinations of i . The resulting equation is given by d P dt = K P − S − , (3)where S − represents the perfect absorbers that alwayskeep ( P ) i = 0 for the sink states of i . Without the lossof generality, the sink states are assumed to be the statesof i = 1 , . . . , m , and the other states, which are free fromthe absorbers, are the remainders of i = m +1 , . . . , n . Theperfect absorber conditions are represented as follows: S − = ( s , . . . , s m , , . . . , T ≡ (cid:18) s − (cid:19) , (4) P = (0 , . . . , , p m +1 , . . . , p n ) T ≡ (cid:18) p (cid:19) . (5)By substituting Eqs. (4) and (5) for Eq. (3), we havethe following solution with the initial condition of P = (cid:18) p (cid:19) : p ( t ) = e tK FF p , (6) s − ( t ) = K SF e tK FF p , (7)where p satisfies k p k ≡ n − m X i =1 | ( p ) i | = n − m X i =1 ( p ) i = 1 . (8) K F F is the submatrix with dimension ( n − m ) × ( n − m )formed by selecting the rows of K from m + 1 to n andthe columns from m + 1 to n , and K SF is the submatrixwith dimension m × ( n − m ) formed by selecting the rowsfrom 1 to m and the columns from m + 1 to n . Theprobability conservation property of the rate matrix of K can be represented by m X i =1 ( K SF ) ij + n − m X i =1 ( K F F ) ij = 0 (9)for j = 1 , , . . . , n − m . Multiplying both sides of Eq. (9)by ( K − F F ) jk and adding the resultant equations from j =1 to n − m , we have the following equations m X i =1 ( − K SF K − F F ) ik = 1 ( k = 1 , , . . . , n − m ) . (10)The i th element of s − ( t ) describes the first passagetime distribution of being absorbed in the i th sink attime t . Therefore, by integrating Eq. (7) from t = 0 to ∞ , the probability of being absorbed in the i th sink for0 t < ∞ is given by the i th element of¯ s − = Z ∞ s − ( t ) dt = − K SF K − F F p , (11)where we use Eq. (7) and lim t →∞ e tK FF = 0, which holdsbecause all eigenvalues of K F F are negative values. FromEq. (11) and ( s − ( t )) i >
0, we see that (¯ s − ) i >
0. More- over, with the use of Eqs. (8) and (10), we have k ¯ s − k ≡ m X i =1 | (¯ s − ) i | = X i (¯ s − ) i = X i (cid:0) − K SF K − F F p (cid:1) i = X i,k (cid:0) − K SF K − F F (cid:1) ik ( p ) k = X k ( p ) k = k p k = 1 , whence k ¯ s − k = k p k = 1 . (12)The conditional probability distribution ρ i ( t ) of beingabsorbed at time t when the system is known to be ab-sorbed in the state of i is given by ρ i ( t ) = [ s − ( t )] i (¯ s − ) i = ( K SF e tK FF p ) i ( − K SF K − F F p ) i . (13)Therefore, the mean first passage time, t i,j , from thestate j to the sink state i is given by t i,j = Z ∞ tρ i ( t ) dt = ( K SF K − F F p ) i ( − K SF K − F F p ) i (14)with ( p ) k = δ k,j − m ( k = 1 , , . . . , n − m ). Moreover,the mean first passage time, t j , from the state j to anyabsorbing states is given by t j = X i (¯ s − ) i t i,j = X i (cid:0) K SF K − F F p (cid:1) i = X i,k (cid:0) − K SF K − F F (cid:1) ik (cid:0) − K − F F p (cid:1) k = X k (cid:0) − K − F F p (cid:1) k , that is t j = k − K − F F p k , (15)where Eqs. (10), (11), and (14) are used.Equations (14) and (15) are the basic formulas for cal-culating the mean first passage times.Next, we show that the mean first passage times canbe evaluated from a stationary population distribution.Let us consider the following mean residence time distri-bution, ¯ p = Z ∞ p ( t ) dt = − K − F F p , (16)where (¯ p ) i is the mean residence time in the ( i + m )thstate for i = 1 , . . . , n − m . Hence, the mean residencetime in the whole system is given by the sum, k ¯ p k , of(¯ p ) i from i = 1 to n − m , which is, of course, equivalentto the mean first passage time t j of Eq. (15).Equation (16) enables us to confirm that ¯ p satisfies thefollowing non-equilibrium stationary state equation: d ¯ p dt = K F F ¯ p + ¯ s + = (17)with ¯ s + = p . Hence, we can interpret ¯ s + as the sourceterm that adds one particle with distribution p per unittime, ¯ p as the stationary population of Eq. (17), k ¯ p k as the total population contained in ¯ p , and k ¯ p k − as theprobabilistic flow carried by one particle. Namely, we canalso compute the mean first passage times as the totalnumbers, k ¯ P k , of particles in the stationary population¯ P obeying the following stationary equation: ddt ¯ P = K ¯ P + S + − S − = , (18)where¯ P = (cid:18) ¯ p (cid:19) , S + = (cid:18) p (cid:19) , S − = (cid:18) ¯ s − (cid:19) . (19)Note that the stationary population equation (18) willbe used in Sec. VI. III. KCL NANOCLUSTER VACANCYDIFFUSION MODEL
In this section, according to Ref. [30], we first presentthe vacancy diffusion model of KCl nanoclusters as anexample of a practical problem, and then we introducethe corresponding Markov state model of the vacancydiffusion.
A. Local minima and saddle points on thepotential energy surface of a KCl nanocluster
Let us assume that one chlorine ion is extracted from acube of ionic crystal with equal N L -atom edges and fur-ther that N L is an odd number 2 n L +1, and the resultantcluster with N ≡ N L − v ( r ij ) = Q i Q j πǫ r ij + A ij exp (cid:18) R i + R j − r ij ρ (cid:19) , (20)where Q i , Q j are the charges of the i th and j th atoms, ǫ is the vacuum permittivity, and r ij is the distancebetween the i th and j th atoms. We use the values of thethree parameters A ij , R i , and ρ that were introduced byTosi and Fumi in Ref. [63]: A ij = 0 . . . R i = 1 .
463 and 1 .
585 ˚A for K and Cl, respectively; and ρ = 0 .
337 ˚A. Then, the total potential energy of thecluster is given by V ( r , . . . , r N ) = N − X i =1 N X j = i +1 v ( r ij ) . (21)In the course of the time evolution, the vacancy movesaround the cluster, which introduces atomic mixing tothe cluster. Note that the cubic form of the clusteris kept with the time evolution when the temperatureis sufficiently low [29]. At such low temperatures, theposition of the vacancy is specified by the cubic latticepoint n = ( n x , n y , n z ) with − n L n x , n y , n z n L .Moreover, we are able to find the atomic structure corre-sponding to the vacancy lattice point n as follows: First,atoms are arranged at d ( m x , m y , m z ) with lattice con-stant d = 3 .
147 ˚A for KCl, where ( m x , m y , m z ) = n and − n L m x , m y , m z n L . Then, the configuration of theatoms is relaxed to the local minimum (LM) configura-tion, r = ( r , . . . , r N ), of the potential energy surface,e.g., by the conjugate gradient method [53]. In this way, n is assigned to the LM atomic structure as r n = r .We compute the LM configurations r n and the energies V ( r n ) for all n . The LM datasets of V ( r n ) and r n arestored in a file in nondecreasing order of energy V ( r n ).For the sake of notational simplicity, the i th lowest en-ergy is denoted as E i , and the corresponding LM, atomicconfiguration, and vacancy lattice point are denoted as i , r i , and n i , respectively. Then, we proceed to find out allof the saddle points (SPs) connecting the adjacent LMs,e.g., by the nudged elastic band method [32]. The corre-sponding saddle point connecting the i th and j th LMs,atomic configuration, and potential energy are denotedas ij , r ij , and E ij , respectively. For the computationaldetails of enumerating all of the LMs and SPs, we referthe reader to Ref. [30]. B. Markov state model of KCl vacancy diffusion
Let f ( r ) denote the probability density function at aconfiguration r . We suppose that the intra-LM relax-ations are so fast that f ( r ) is represented as f ( r ) = p f ( r ) + p f ( r ) + · · · + p n f n ( r ) , (22)where n is the number of the LMs, f i ( r ) is the localequilibrium in the i th-LM basin, and p i is the probabilitythat r is in the i th-LM basin. We identify p i with theprobability in the i th state of the Markov state model.Then, the probability vector P of the Markov state modelis given by P = ( p , p , . . . , p n ) T .Next, we evaluate the transition rate k i,j from the j thto the adjacent i th state, when the potential barrier en-ergies are sufficiently larger than the average kinetic en-ergy of k B T / T , where k B denotes the Boltzmann constant. In thiscase, the transition rate k i,j from the j th to the i th stateis given by k i,j = ν i,j exp (cid:18) − E ij − E j k B T (cid:19) . (23)Here, the prefactor ν i,j , called a frequency factor, is givenby ν i,j = Q ′ k ( ν i ) k Q ′ k ( ν ij ) k , (24)where ν i and ν ij are vibrational frequency vectors thatare calculated from the Hessians at r i and r ij , respec-tively. The product Q ′ k ( ν ∗ ) k denotes the partial prod-uct of the positive frequency modes ( ν ∗ ) k >
0, wherethe imaginary frequency modes and the zero frequencymodes are left out from the products.Finally, the transition rate matrix K is given by( K ) i,j = k i,j ( i = j ) and ( K ) j,j = − X i = j k i,j , (25)where the probability conservation equations P i ( K ) i,j =0 for j = 1 , , . . . , n and the detailed balance conditions( K ) i,j ( P eq ) j = ( K ) j,i ( P eq ) i for i , j = 1 , , . . . , n aresatisfied since ( P eq ) i ∝ exp( − E i /k B T ).In the following, we examine the slowest relaxationmodes of the KCl vacancy diffusion model of N L = 13.To this end, we searched for the LMs and the SPs of thecluster, thereby finding 1099 (= 1 + ( N L − /
2) LMs and5472 SPs. Then, with the use of Eqs. (23), (24), and (25),we formed its rate matrix of K at k B T = 0 .
03 eV, whosematrix dimension n is 1099 and the number of nonzerooffdiagonal elements is 10948 (= 5472 × K , we obtained the eigenvalues λ i and the corresponding relaxation modes P i for i =0 , , , . . . , n −
1, where P = P eq . In Sec. IV andSec. V, we study the properties of the slowest relaxation P and the second slowest relaxations P , P , and P ( λ = λ = λ ), respectively. IV. THE SLOWEST RELAXATION MODE
In this section, we show that the slowest relaxationmode of the KCl nanocluster makes effective use of thefast surface diffusion of the cluster, in terms of the meanfirst passage times, the free energy landscapes, and theatomic interactions.
A. Dominant pathways
Figure 1(a) shows the slowest relaxation mode of P ,from which we see that P has the probability excessesat around the origin (0 , ,
0) and the probability short-ages at around the eight vertices of ( ± n L , ± n L , ± n L ). In ● ● ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ - - - - - - FIG. 1. The slowest relaxation mode, P , at k B T = 0 .
03 eV:(a) Red (light gray) and blue (dark gray) balls are depictedat vacancy lattice points of n i with radii ∝ | ( P ) i | / for( P ) i > P ) i <
0, respectively. Note that P hascubic symmetry around the x -, y -, and z -axes. (b) ( P ) i alongpathways ℓ of Eq. (26) and ℓ of Eq. (27) are plotted as afunction of the number of steps from the origin (0 , ,
0) withblue circles connected by the lower line and orange squaresconnected by the upper line, respectively.
Fig. 1(b), we also plot the values of ( P ) i along two path-ways from the origin to the vertex ( n L , n L , n L ), whichclearly shows that they have the maximum values at theorigin and positive values up to three steps from the ori-gin and negative values at the vertices. P decays with the rate of λ = − . × s − overthe course of time. Hence, the probabilistic flow from thecenter to the vertices is expected in the relaxation pro-cess of P . To confirm this, we compute all of the prob-abilistic flows f i,j from j to i , generated by P , where f i,j is given by f i,j = k i,j p j − k j,i p i with p i = ( P ) i for1 i, j n . In Fig. 2, f i,j are represented by the arrowsfrom n j to n i when f i,j >
0. We see that the probabilis-tic flows from the center to the vertices are generated.More precisely, the probabilistic flows are not uniformbut mostly along the two dominant pathways of ℓ and ℓ as depicted in Fig. 2. FIG. 2. Probabilistic flows f i,j of P at k B T = 0 .
03 eV arerepresented by arrows from n j to n i with cylinder radii ∝ p | f i,j | for f i,j >
0. The flows have the same cubic symmetryas in Fig. 1. Hence, only the flows in a reduced zone n y > n x > n z > ℓ and ℓ in blue (dark gray) carry the dominant flows: ℓ iscomposed of the straight move from the origin (0 , ,
0) to theedge center (6 , ,
0) and the succeeding zigzag move to thevertex (6 , , ℓ is composed of the zigzag move from theorigin to the face center (0 , ,
0) and the succeeding straightmove to the vertex. [See Eqs. (26) and (27).]
The dominant pathways of ℓ and ℓ are defined bythe following algorithm that searches for the maximumflow pathways flowing into the terminals. First, we startwith the terminal of the vertex ( n L , n L , n L ) = (6 , , , ,
5) flows into the ter-minal (6 , , , ,
5) and(4 , ,
4) flow into (5 , , , ,
5) and (4 , ,
4) are denoted as ℓ and ℓ , respectively. We then search for the source flow of ℓ as follows. The probabilistic flow from (5 , ,
4) is themaximum flow flowing into (6 , , , ,
3) isthe maximum flow flowing into (5 , , , ,
0) appears andgives the dominant pathway as ℓ =(0 , , → (1 , , → · · · → (6 , , → (5 , , → (6 , , → (5 , ,
3) (26) → (6 , , → (5 , , → (6 , , . The dominant paths of this kind are composed of sixstraight steps from the origin to the 12 centers of theedges ( ± , ± , ± , , ± , ± , ± , ,
4) and obtain ℓ as ℓ =(0 , , → (1 , , → (0 , , → (1 , , → (0 , , → (1 , , → (0 , , → (27)(1 , , → · · · → (5 , , → (6 , , . The dominant paths of the second kind are composedof six zigzag steps from the center to the six centers ofthe faces, ( ± , , , ± , , , ± ℓ -type paths, which climb along thethree edges connected to the vertex, and the three ℓ -typepaths, which move across the three faces containing thevertex. Note that these observations are consistent withour previous results from Ref. [42]. There, all states aredivided into groups, called metabasins, that are locatedaround the vertices, the edges, the faces, and the centerpart, and then the relaxation processes are described ac-curately by the renormalized transitions between thesemetabasins. That is to say, we have reconfirmed herethat the essential pathways connecting the vertices, theedges, the faces, and the center part are indispensable fordescribing the slowest transport of probabilities. B. Mean first passage times
Here, we consider why the dominant paths carryinglarge probabilistic flows are not almost straight, nine-step shortest paths from the origin to the vertices, suchas (0 , , → (0 , , → (1 , , → (2 , , → (2 , , → (2 , , → (3 , , → (4 , , → (5 , , → (6 , , ℓ and ℓ in Fig. 2. To this end,we examine the mean first passage times of t j from var-ious initial states of j to the sink states of the vertices( ± n L , ± n L , ± n L ).Using Eq. (15), we compute t j for various initial statesof j . The resulting t j are plotted in Fig. 3(a). We see thatthe states in the central part have large t j values, whilethe states on the surface have quite small values. Thatis, the KCl nanocluster is a hybrid system that combinesentirely different microscopic diffusive regions: The cen-tral part is the region that is hard to move stochastically,whereas the surface part is the region that is quite easyto move. To see this more closely, we plot t j along thepaths of ℓ and ℓ in Fig. 3(b), where both of t j havethe maximum value at the origin and they are negligiblysmall compared to the maximum value in two layers fromthe surface.Now, we see the reason why the detoured pathwaysare selected to be the dominant pathways, as depicted inFig. 2. Namely, it is because all the dominant pathwaysprefer to pass the slow diffusion region of the central partas soon as possible, with the fewest steps of n L = 6, inorder to make the most effective use of the fast diffusionin the surface region.Next, we show that λ can be evaluated approximatelyfrom t j . The longest mean first passage time to the ● ● ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ FIG. 3. Mean first passage times of t j at k B T = 0 .
03 eV,with perfect sinks connected to the vertices ( ± , ± , ± t j are represented by balls of radii ∝ | t j | / located at n j .The mean first passage times have cubic symmetry aroundthe x -, y -, and z -axes. (b) t j along ℓ and ℓ are plottedas a function of the number of steps from the origin withblue circles connected by the lower line and orange squaresconnected by the upper line, respectively. Both of t j have themaximum value of 8 . × − s at the origin of (0 , , t j in the center part, where the number of steps isfrom 0 to 4, have the same order of magnitude, while thosein the two layers from the surface, where the number of stepsis from 5 to 12, are quite small values. vertices is t (0 , , = 8 . × − s. The probability ofbeing at the vertices in equilibrium is P eq (vertices) = P i ∈ vertices ( P eq ) i = 0 . P eq (vertices).Then, the equilibration time is approximately given by P eq (vertices) × t (0 , , = 7 . × − s . (28)The corresponding equilibration rate is given by the in-verse of the equilibration time, 1 . × s − . The FIG. 4. Free energy sequences in units of eV are plottedat k B T = 0 .
03 eV, as functions of steps counted from theorigin (0 , ,
0) to a vertex along the geometric shortest path( N , green), along the dominant path, ℓ , of Eq. (26) via anedge center ( • , blue) and the dominant path, ℓ , of Eq. (27)via a face center ( (cid:4) , orange). The integer steps indicate thefree energies of the local minima, and the half-integer stepsindicate the free energies of the saddle points that connectthe basins of adjacent local minima. estimate agrees qualitatively with the values of | λ | =1 . × s − , although it is a smaller value than | λ | .This discrepancy arises because, although the actualexcess probabilities in P are distributed in the centralpart as depicted in Fig. 1, the excess probability is ap-proximated to the distribution concentrated on the ori-gin, for the mean first passage time approximation of λ .We will revisit this point in Sec. VI. C. Free energy sequences
Here, we consider the physical reason why the bottle-neck of diffusion in the Markov state model of the KClnanocluster is located at the central part. To this end,we examine the following free energies for the LMs of j and the SPs of ij , respectively: F j = E j − k B T ln Y k ′ ( ν j ) k , (29) F ij = E ij − k B T ln Y k ′ ( ν ij ) k , (30)where k i,j = exp[ − β ( F ij − F j )] holds. Then, the freeenergy sequence of local minima and saddle pointsalong a pathway i → i → · · · → i s is given by F i , F i i , F i , F i i , F i , . . . F i s .In Fig. 4, we plot the free energy sequences of thedominant pathways of ℓ [Eq. (26)] and ℓ [Eq. (27)].Along these dominant pathways, the activation energiesfor the inner transitions i k → i k +1 are about ∆ F i k +1 ,i k = F i k i k +1 − F i k ≈ . k i k +1 ,i k ≈ × s − at k B T = 0 .
03 eV. In contrast, those for the surface transi-tions are ∆ F i k +1 ,i k ≈ . k i k +1 ,i k ≈ × s − , which are about 10 timeshigher than the inner rates, at the same temperature.For comparison, we also plot the free energy sequencealong the nine-step geometric shortest path in Fig. 4.We see that the first seven steps are in the slow diffusionregion and the last two steps are in the faster surfacediffusion region. Therefore, the geometric shortest pathcannot be dominant, because the extra steps in the slowerdiffusion region reduce its diffusive flow drastically.Note that the activation free energies of ∆ F i k ,i k +1 & . k B T = 0 .
03 eV andhence the harmonic approximation (23) used in thisstudy is accurate.
D. Activation energies on potential energylandscapes
Next, we show that it can be understood in terms ofthe interatomic interaction energies why the surface ac-tivation free energies are so small compared to the innerones.First, from Eqs. (24), (29), and (30), ∆ F i,j = E ij − E j + k B T ln ν i,j holds. At the low temperature of k B T =0 .
03 eV, k B T ln ν i,j ≈ .
03 eV is negligible compared to∆ F i,j ≈ . F i,j ≈ E ij − E j holds. Hence, inthe following we consider the activation energies, ∆ E i,j = E ij − E j , of various transitions.We examine the activation energy of ∆ E i,j when thevacancy lattice points n i and n j are present in the innerpart of the cluster. In this case, E i and E j are almostthe same, as shown in Fig. 4, and hence ∆ E i,j is deter-mined by the energy increase from E i ≃ E j , due to thedeformation of the crystal structure near the lattice de-fect. To quantify the deformation energies, we introducethe individual potential energy of the k th atom as V k = 12 X l = k v ( r kl ) . (31)Then, the total potential energy of Eq. (21) is representedas V ( r ) = X k V k . (32)Note that V k is half of the required energy to remove the k th atom from the cluster, since V ( r , . . . , r k , . . . r N ) − V ( r , . . . , r k − , r k +1 . . . , r N ) = 2 V k holds. Figures 5(a)and 5(c), respectively, show the values of V k for theLMs of the vacancy lattice points of n = (0 , ,
0) and n = (1 , , V k are con-centrated in the vicinities of the vacancies at around r = d n and r = d n with lattice constant d , andthat V k of Cl and K atoms decrease and increase, re-spectively, when approaching the vacancy positions. Fig-ure 5(b) shows the values of V k for the SP connectingthe LMs of n and n . The SP has the high energyCl atom as the lattice defect at around the midpoint, FIG. 5. Individual energies V k of k th atoms contained inthe z = 0 planes are shown for (a) the LM of the vacancylattice point of n = (0 , , n and n = (1 , , n . The Cl and K atoms of individual energies V k are represented by blue (dark gray) and yellow (light gray)balls with radii ∝ p V k − min k { V k } . The defect neighborsare defined by the regions inside the green (gray) frames of | x − y | < . d , − . d < x + y < . d , and | z | < .
2. (See thetext.) r , = d (1 / , / , V k of Cl and K atoms decrease and increase,respectively, when approaching the defect of the high-energy Cl atom.Next, we show that ∆ E i,j can be estimated withuse of the local V k values around the defects. To thisend, we obtain the local energies E loc.1 = − .
34 eV, E loc.1 , = − .
804 eV, and E loc.2 = − .
334 eV, whichare the sums of V k inside the local regions surrounded bythe green (gray) rectangle frames in Figs. 5(a), 5(b), and5(c), respectively. Hence, the activation energies eval-uated from these local energies are given by ∆ E loc.2 , = E loc.1 , − E loc.1 = 0 .
530 eV and ∆ E loc.1 , = E loc.1 , − E loc.2 =0 .
536 eV, which agree qualitatively with the exact acti-vation energies of ∆ E , = 0 .
58 eV and ∆ E , = 0 . E face ← face ≈ ∆ E inner ← inner / E edge ← face ≈ ∆ E inner ← inner /
4, which are impliedin Table I. To understand these relations, we assumefor simplicity that the deformation energy is uniformlydistributed inside the ball of radius a ≈ d located atthe defect point. Assuming further that the deforma-tion energy per unit volume is given by ǫ , then theactivation energies for inner vacancies are estimated as∆ E inner ← inner = 4 πa ǫ/ .
58 eV).Next, we consider the activation energies for vacan-cies in a face. In this case, the energies of adjacent localminima are also supposed to be the same for simplic-ity. Since the deformed regions are half of the innercase, the deformation energies of the vacancies in thefaces are estimated to be ∆ E face ← face = 4 πa ǫ/ / E inner ← inner / .
29 eV).Moreover, when a vacancy inside a face moves toan adjacent edge, the deformation energy reduces tohalf of ∆ E face ← face , since the deformed region is halvedfrom that of the transition in a face. Hence, we have∆ E edge ← face ≈ ∆ E inner ← inner / .
15 eV). We seethat the estimated values agree quantitatively with theexact values. Lastly, when the vacancy in a face moves toan adjacent vertex, the deformation energy is evaluatedto be halved to ∆ E edge ← face . Thus, we have the follow-ing approximation: ∆ E vertex ← face ≈ ∆ E inner ← inner / .
07 eV), which agrees qualitatively with the exact valueof 0 .
027 eV.Here, we have revealed that the activation energiesfor the system of n L = 6 are determined by thelocal deformation energies of ∆ E loc .i,j around the de-fects. Accordingly, the activation energies are supposedto be almost independent of the system size of n L .In fact, we have ∆ E inner ← inner = ∆ E (0 , , ← (0 , , =0 .
58 eV, ∆ E face ← face = ∆ E (4 , , ← (4 , , = 0 .
33 eV,∆ E edge ← face = ∆ E (4 , , ← (4 , , = 0 .
16 eV, and∆ E vertex ← face = ∆ E (4 , , ← (4 , , = 0 .
03 eV for n L = 4.These results show that all types of activation energiesare indeed almost independent of the system size when n L > n L = 2, we have∆ E inner ← inner = ∆ E (0 , , ← (0 , , = 0 .
47 eV,∆ E face ← face = ∆ E (2 , , ← (2 , , = 0 .
39 eV,∆ E edge ← face = ∆ E (2 , , ← (2 , , = 0 . E vertex ← face = ∆ E (2 , , ← (2 , , = 0 .
05 eV, whichshows that the uniform local deformation assumption forthe activation energies employed above does not hold for n L = 2. In other words, the cluster of n L = 2 is too smallto separate the deformations of the surface from thoseof the central portion, and thus some non-negligiblecouplings are generated between the inner and surfacedeformations. As a result of the couplings, the relativelyhigh activation energies between inner transitions aredecreased, while the other relatively low activationenergies between surface transitions are increased for n L = 2.In addition, the saddle connectivity graphs of n L = 4,6, and 8 depicted in Ref. [30] also show visually thatthe activation energies of ∆ E inner ← inner , ∆ E face ← face ,∆ E edge ← face , and ∆ E vertex ← face are almost independentof the sizes n L of the clusters. V. THE SECOND SLOWEST RELAXATIONS
In this section, we examine the second slowest relax-ations of λ , λ , λ = − . × s − .In Fig. 6, we plot P in the same way as in Fig. 1.The probability deviations of P are polarized in the x -direction. Here, the probability excess is in the regionof x >
0, the probability shortage is in x <
0, and theprobabilities are zero in x = 0. From this observation, therelaxation process is expected as follows: the probabilityexcess moves in the opposite x -direction, the probabilityshortage moves in the x -direction, and these pairs meet FIG. 6. We plot the second slowest relaxation mode of P in the same manner as in Fig. 1. P has the four-fold rota-tional symmetry around the x -axis. The relaxation mode ispolarized in the x -direction, where the excess and shortage ofprobability are distributed, respectively, in x > x < x = 0. Similarly to P , P and P are the relaxation modes, which are polarizedin the y - and z -directions, respectively. TABLE I. Activation energies ∆ E i,j and local activation energies ∆ E loc .i,j of vacancy transitions in the KCl cluster of n L = 6are enumerated in units of eV. ∆ E , = E , − E and ∆ E loc.1 , = E loc.1 , − E loc.2 , where E loc.1 , and E loc.2 are sums of individualatomic energies around the defect points. (See the text.)Activation type n ← n ∆ E , ∆ E loc . , Inner ← Inner (1 , , ← (0 , ,
0) 0 .
58 0 . ← Face (6 , , ← (6 , ,
0) 0 .
32 0 . ← Face (6 , , ← (6 , ,
1) 0 .
15 0 . ← Face (6 , , ← (6 , ,
5) 0 .
027 0 . with each other in the region of x = 0, to be annihilated.To confirm this expectation, we evaluated the meanfirst passage times with sinks connected to (0 , n y , n z )for − n L n y , n z n L . The longest passage time is t (2 , , = 4 . × − s. The resulting rate of this pro-cess is 2 . × s − . Also here, the estimated values ofthe rate agree qualitatively with | λ | , but they are some-what smaller than the exact rate of | λ | = 3 . × s − ,because this approximate rate is evaluated only from thelongest mean first passage time, as discussed in Sec. IV B.Finally, we show that the approximate relation of λ ≈ λ holds. Here, the value of λ evaluated from themean first passage time from (2 , ,
0) to x = 0 is approx-imated by that from (2 , ,
0) to (0 , , λ is approximately evaluated from the meanfirst passage time from (0 , ,
0) to (4 , , , ,
0) to (2 , , / | λ | : 1 / | λ | = 4 : 2 holds, and thus λ = 2 λ holds. Similarly, we can derive the approxi-mate relations of λ = 2 λ and λ = 2 λ .In this section, we have confirmed that the second slow-est relaxations of P , P , and P , respectively, smoothout the nonequilibrium distribution deviations in the x -, y -, and z -directions in the course of time, and the bottle-neck processes for the second slowest relaxations are alsothe slow diffusions inside the cluster. This fact allows usto derive the approximate relation of λ , λ , λ ≈ λ . VI. SYMMETRIC EVALUATION OFRELAXATION RATESA. Symmetric population method
In the previous sections, we have successfully evalu-ated the values of the slowest and second slowest relax-ation rates with the use of the mean first passage times.Recall that the usages of the mean first passage timesfor the slowest and second slowest relaxations were dif-ferent. That is, for the slowest relaxation, the mean firstpassage times concerning particles were used, whereas,for the second slowest relaxations, the mean annihilationtimes of particle-hole pairs were evaluated with the usesof the mean first passage times. The difference manifests itself in the symmetry of re-laxation times under exchanging the sinks and sources.Namely, the approaches to evaluating the slowest andsecond slowest relaxation rates, respectively, give the re-laxation times that are asymmetric and symmetric for ex-changing sinks and sources. In fact, as shown in Eq. (28),the slowest relaxation time is 7 . × − s with the sourceand the sinks being connected to the origin and the ver-tices, respectively, whereas, with the sources and the sinkbeing connected to the vertices and the origin, the meanfirst passage time is given by τ ′ = 3 . × − s and thusthe relaxation time is ( P eq ) (0 , , × τ ′ = 8 . × − s.The symmetry corresponds to the property that P i and − P i have the same value of λ i , and hence it is requiredfor a consistent treatment.Here, we develop an alternative population method forestimating mean first passage times that is symmetricunder the exchange of the sinks and sources.To this end, we consider the following stationary pop-ulation equation: K Q + S = , (33)with k Q + k = k Q − k , k S + k = k S − k = 1 . (34)Here, Q = Q + − Q − and S = S + − S − , where the pos-itive population Q + and the negative population − Q − are given by Q + = Q + | Q | , Q − = − Q − | Q | , (35)with | Q | ≡ ( | q | , | q | , . . . , | q n | ). The source part S + andthe sink part − S − are, respectively, given by S + = S + | S | , S − = − S − | S | , (36)with | S | ≡ ( | s | , | s | , . . . , | s n | ). With use of the stationarysolution ¯ P , satisfying Eqs. (18) and (19), the stationarysolution Q of Eqs. (33) and (34) is given by Q = ¯ P − k ¯ P k P eq , (37)In fact, Q given in Eq. (37) satisfies the constraint k Q + k = k Q − k , because k Q + k − k Q − k = k Q k = k ( ¯ P − k ¯ P k P eq ) k = k ¯ P k − k ¯ P k = 0 holds.As illustrated in Fig. 7(b), the negative population − Q − can be interpreted as the hole population of Q − .1 source hole source hole source sourcesinksinksinksource FIG. 7. Schematic illustration of the population methods. (a) The stationary particle population of ( ¯ P ) i satisfying Eqs. (18)and (19) is shown as a function of states i . The red (light gray) arrow indicates how the particles injected at the source movediffusively to the sink. We also show k ¯ P k P eq with a dashed line, to illustrate the latter term in Eq. (37). (b) The symmetricstationary population Q with a source and a sink connected to the left and right end states, respectively, is given by Eq. (37).The negative population − Q − is interpreted as the stationary population Q − of holes. The blue (gray) arrow indicates howthe holes injected at the hole source move diffusively. Particles from the source and the holes meet at the center and areannihilated by pair annihilation. (c) The stationary population under the exchange of the sink and source is shown. In thiscase, particles of population Q − move to the left and holes of population Q + move to the right. The particles and the holesmeet at the center to be annihilated. Both of the mean annihilation times in (b) and (c) agree with the mean first passagetime, k Q + k = k Q − k , of particles and holes. (See the text.) Hence, the mean annihilation times of particles and holesare, respectively, given by the first passage times of k Q + k and k Q − k as discussed in Sec. V. Similarly, S − is inter-preted as a hole source part that adds one hole per unittime. Hence, the constraint of k S + k = k S − k = 1 meansthat S + adds one particle per unit time, and S − addsone hole per unit time. Note that, since P eq satisfies thedetailed balance condition, the particle flows of ¯ P and Q are the same, and so are their dominant pathways, asillustrated in Figs. 7(a) and 7(b).The constraint k Q + k = k Q − k means that the meanfirst passage times are symmetric under exchanging thesinks and sources. In fact, by exchanging the sinks andsources, S and Q are, respectively, converted to − S and − Q , and hence Q + and Q − are, respectively, convertedto Q − and Q + , as illustrated in Fig. 7(c). Hence, themean first passage times of Eqs. (33) and (34) satisfy k Q − k = k Q + k with sinks and sources exchanged, whichis the same value as the value before the exchange.Here, with this symmetric population method, we eval-uate the slowest relaxation rate of λ for the vacancydiffusion model of the KCl nanocluster. Setting( S + ) (0 , , = 1 , ( S − ) ( ± n L , ± n L , ± n L ) = 1 / , (38)and otherwise ( S ± ) i,j,k = 0, we evaluated the symmet-ric stationary solution of Eq. (41), thereby obtaining theslowest relaxation time of k Q + k = k Q − k = 7 . × − s.Namely, with this symmetric method, we obtained anapproximation of the exact slowest relaxation time of − /λ = 5 . × − s, which is symmetric under ex-changing the sinks and sources and as accurate as theasymmetric result of 7 . × − s given in Sec. IV. For the second slowest relaxation, we set( S + ) (2 , , = 1 , ( S − ) ( − , , = − , (39)and otherwise ( S ± ) i,j,k = 0, thereby obtaining the sym-metric result of k Q + k = k Q − k = 4 . × − s, which ofcourse agrees with the result given in Sec. V.In this subsection, we have developed the symmetricpopulation method for mean first passage times, whichenables us to approximately evaluate the slowest relax-ation times symmetrically by exchanging the sinks andsources. B. Symmetric population method as an inversepower method
Here, we show that the iterative use of the symmetricpopulation method enables us to compute the slowestrelaxation times accurately.First, the n × n matrix P = ( P eq , P , P , . . . , P n − )is invertible, where P i is the eigenvector correspondingto the i th relaxation mode. We expand S and Q as S = s ′ P eq + s ′ P + s ′ P + . . . , (40) Q = q ′ P eq + q ′ P + q ′ P + . . . , (41)where the coefficients s ′ i and q ′ i are defined as follows: s ′ i = ( P − S ) i , q ′ i = ( P − P ) i . (42)With the use of s ′ i and q ′ i , Eq. (33) is represented as q ′ = s ′ = 0 , (43) q ′ i = s ′ i − λ i ( i = 1 , , . . . ) . (44)2Substituting Eqs. (43) and (44) into Eq. (41), we have Q = s ′ − λ P + s ′ − λ P + . . . . (45)Equations (42) and (45) define the procedure to obtain Q from S , which is denoted as Q = Q ( S ).Then, the mean first passage time approximation τ ofthe relaxation time is written as follows: τ = k Q + kk S + k = k Q − kk S − k = k Q kk S k = k Q ( S ) kk S k , (46)as discussed in the previous subsection.Now, we consider the effect of the iterative use of thisprocedure. We apply the m -times function compositionof Q to S , and the resultant vector Q ( m ) ( S ) is given by Q ( m ) ( S ) = s ′ ( − λ ) m P + s ′ ( − λ ) m P + . . . , (47)which shows that, as m → ∞ , if s ′ = 0, then Q ∝ P ,else if s ′ = · · · = s ′ i − = 0 and s ′ i = 0, then Q ∝ P i .From this, we understand that the symmetric populationmethod can be interpreted as an inverse power methodfor the eigenvalue problem [53].In Fig. 8(a), we plot the mean first passage times τ = k Q ( m ) ( S ) kk Q ( m − ( S ) k = k Q ( Q ( m − ( S )) kk Q ( m − ( S ) k for the slowest and second slowest relaxation modes withthe settings of Eqs. (38) and (39), respectively. At m = 1,both are larger than the corresponding exact relaxationtimes, because S of Eqs. (38) and (39) are selected so asto maximize the mean first passage times of k Q ( S ) k / k S k .At m = 2 and 3, the mean first passage time approxi-mations almost converge to the corresponding relaxationtimes. These findings show that the mean first passagetime approximations of the slowest and the second slow-est relaxation times satisfy s ′ = 0 and s ′ = 0 , s ′ = 0, re-spectively. Moreover, both of S are sufficiently accurate,so as to converge to the slowest and second slowest relax-ation modes, respectively, with a few iterations. Namely, S of Eqs. (38) and (39) closely approximate the exacteigenvectors of P i ( i = 1 , S are the drasticsimplifications of P i with very few sinks and sources.To see the convergence of Q ( m ) ( S ) to the eigenvectorswith the symmetric population method, we plot Q ( S )for the slowest and second slowest relaxation modes inFigs. 8(b) and 8(c), respectively. Comparing these graphswith Figs. 1(a) and 6, respectively, we see that Q ( S )with Eqs. (38) and (39) are almost the same with P and P . That is, we can obtain accurate approximateeigenvectors by applying this procedure just once to thequite simplified sinks and sources of S .We remark finally that when it is difficult to set S tobe in the convergence region of P , we can obtain thefirst and second slowest relaxation modes simultaneouslyby iteratively applying Q to two vectors that span a two-dimensional subspace and orthogonalizing the vectors, asin the general diagonalization algorithms [53]. ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ FIG. 8. The mean first passage times of the symmetric pop-ulation method converge to the slowest relaxation times. (a)The mean first passage times, τ = k Q ( m ) ( S ) k / k Q ( m − ( S ) k ,are plotted as a function of iteration number m ( m =1 , , . . . ,
10) with circles for the slowest relaxation and withsquares for the second slowest relaxation. The dashed anddotted lines show the values of the exact relaxation times of1 / | λ | and 1 / | λ | , respectively. Q ( S ) are plotted in the samemanner as in Fig. 1(a), for the slowest relaxation [panel (b)]and for the second slowest relaxation [panel(c)]. VII. SUMMARY
We studied the slowest and second slowest relaxationsof vacancy diffusion in a KCl nanocluster.In Sec. IV, we found that the slowest relaxation modeof P has cubic symmetry around the origin (0 , , ± n L , ± n L , ± n L ) over the course oftime evolution. We also found that the dominant path-ways that carry large diffusive flows are classified into twotypes of pathways from the origin to the vertices. Oneis through the face centers, and the other is through theedge centers.To understand why these pathways are selected asdominant pathways, we estimated the mean first pas-sage times from various states to the vertex sinks. As aresult, the surface diffusion turned out to be about 10 times faster than the surface diffusion at room temper-ature of k B T = 0 .
03 eV. Hence, the dominant pathwaysturned out to be the shortest pathways to the surfaces.There, we also gave an approximation of the slowest re-laxation rate λ with the use of the mean first passagetimes.Next, the reasons for the slow inner and fast surfacediffusions were studied in terms of the free energy land-scape. The sequences of free energies at minima andsaddle points along the two types of pathways were ex-amined. We found that the activation free energies inthe inner region are about twice as large as those in thesurface region, which explains the drastic slow inner dif-fusion. We also gave an intuitive explanation for the ratioof the activation energies that are leading terms of theactivation free energies, with the use of the individualatomic energies of V k .In Sec. V, we considered the second slowest relaxationmodes. With use of the three-dimensional plot of Fig. 6,the second relaxation modes of P , P , and P turnedout to correspond to relaxation of the excesses and thedeficiencies of probability in the x - y -, and z -directions,respectively. The second slowest relaxation rate of λ is also success-fully estimated by use of the mean first passage timeswith sinks connected to the region of x = 0. There,the intuitive explanation for the approximate relation λ , λ , λ ≈ λ was given in terms of the free energylandscapes along the dominant pathways.In Sec. VI, we have developed a symmetric popula-tion method, which computes the approximate relax-ation rate as the mean first passage times of particlesand holes. The symmetric population method has a rea-sonable property in that both P i and − P i are the eigen-vectors of the same eigenvalue. We have also shown thatiterative use of the symmetric population method enablesus to obtain the accurate slowest relaxation times, simi-larly to the inverse power method of matrix diagonaliza-tion.In summary, we have shown that the properties of theslowest relaxation modes are reconstructed by mean firstpassage times in Markov state models suitably connectedwith sinks and sources. The mean first passage timesare useful to extract the bottleneck processes buried inMarkov state models. We have also shown that the for-mation of the bottlenecks can be understood from thephysical basis of potential energy landscapes that sup-port the networks of the Markov state models. ACKNOWLEDGMENTS
The authors are very grateful to Shoji Tsuji andKankikai for the use of their facilities at Kawaraya duringthe early stage of this study. Y.S and T.O. are supportedby a Grant-in-Aid for Challenging Exploratory Research(Grant No. JP 15K13539) from the Japan Society forthe Promotion of Science. [1] M. Goldstein, J. Chem. Phys. , 3728 (1969).[2] F. H. Stillinger and T. A. Weber, Science , 983(1984).[3] F. H. Stillinger, Science , 1935 (1995).[4] A. Heuer, Phys. Rev. Lett. , 4051 (1997).[5] L. Angelani, G. Parisi, G. Ruocco, and G. Viliani, Phys.Rev. Lett. , 4648 (1998).[6] P. G. Debenedetti and F. H. Stillinger, Nature (London) , 259 (2001).[7] S. Sastry, Nature (London) , 164 (2001).[8] R. A. Denny, D. R. Reichman, and J.-P. Bouchaud, Phys.Rev. Lett. , 025503 (2003).[9] B. Doliwa and A. Heuer, Phys. Rev. Lett. , 235501(2003).[10] G. A. Appignanesi, J. A. Rodr´ıguez Fris, R. A. Montani,and W. Kob, Phys. Rev. Lett. , 057801 (2006).[11] S. De, B. Schaefer, A. Sadeghi, M. Sicher, D. G. Kanhere,and S. Goedecker, Phys. Rev. Lett. , 083401 (2014).[12] Y. Yang and B. Chakraborty, Phys. Rev. E , 011501(2009).[13] O. M. Becker and M. Karplus, J. Chem. Phys. , 1495 (1997).[14] D. Shukla, C.X.Hern´andez, J.K. Weber, and V. S. Pande,Acc. Chem. Res. , 414 (2015).[15] A.B. Kolomeisky, E. B. Stukalin, and A. A. Popov, Phys.Rev. E , 031902 (2005).[16] K. R. Ghusingaa, J. J. Dennehyb, and A.Singh,Proc. Natl. Acad. Sci. (U.S.A.) , 693 (2017).[17] N.-V. Buchete and G. Hummer, J. Phys. Chem. B ,6057 (2008).[18] G. Hummer and A. Szabo, J. Phys. Chem. B ,9029(2015).[19] S. S. Cho, Y. Levy, and P. G. Wolynes, Proc. Natl. Acad.Sci. (U.S.A.) , 586 (2006).[20] G. R. Bowman and V. S. Pande, Proc. Natl. Acad. Sci.(U.S.A.) ,10890 (2010).[21] J. Wang, R.J. Oliveira, X. Chu, P. C. Whitford, J.Chahine, W. Han, E. Wang, J. N. Onuchic, and V.B.P.Leite, Proc. Natl. Acad. Sci. (U.S.A.) , 15763 (2012).[22] F. Pontiggia, D.V. Pachov, M.W. Clarkson, J. Villali,M.F. Hagan, V.S. Pande, and D. Kern, Nat. Commun. , 7284 (2015). [23] B. Zhang, W. Zheng, G.A. Papoian, and P.G. Wolynes,J. Am. Chem. Soc. , 8126 (2016).[24] G. A. Breaux, R. C. Benirschke, T. Sugai, B. S. Kinnear,and M. F. Jarrold, Phys. Rev. Lett. , 215508 (2003).[25] H. Haberland, T. Hippler, J. Donges, O. Kostko, M.Schmidt, and B. von Issendorff Phys. Rev. Lett. ,035701 (2005).[26] K. Joshi, S. Krishnamurty, and D. G. Kanhere, Phys.Rev. Lett. , 135703 (2006).[27] C. Hock, S. Straßburg, H. Haberland, B. v. Issendorff, A.Aguado, and M. Schmidt, Phys. Rev. Lett. , 023401(2008).[28] C. Hock, C. Bartels, S. Straßburg, M. Schmidt, H. Haber-land, B. von Issendorff, and A. Aguado Phys. Rev. Lett. , 043401 (2009).[29] T. Niiyama, S.-I. Sawada, K. S. Ikeda, and Y. Shimizu,Eur. Phys. J. D ,1 (2014).[30] T. Niiyama, T. Okushima, K. S. Ikeda, and Y. Shimizu,Chem. Phys. Lett. , 52 (2016).[31] C. L. Brooks III, J.N. Onuchic, and D. J. Wales, Science , 612(2001).[32] D. J. Wales, Energy Landscapes: Applications to Clus-ters, Biomolecules and Glasses (Cambridge UniversityPress, Cambridge, UK, 2003).[33] F. H. Stillinger,
Energy Landscapes, Inherent Structures,and Condensed-Matter Phenomena (Princeton Univer-sity Press, Princeton, New Jersey, 2016).[34]
An Introduction to Markov State Models and Their Ap-plication to Long Timescale Molecular Simulation , editedby G. R. Bowman, V. S. Pande, and F. No´e (Springer,New York, 2013).[35] J.G. Kemeny and J.L. Snell,
Finite Markov Chains ,(Springer-Verlag, New York, 1976).[36] J.P. Tian and D. Kannan, Stoch. Anal. Appl. , 685(2006).[37] G.R. Bowman, J. Chem. Phys. , 134111 (2012).[38] D. J. Wales, Mol. Phys. , 3285 (2002).[39] K. Klemm, C. Flamm, and P. F. Stadler, Eur. Phys. J.B 63, 387 (2008).[40] T. Okushima, T. Niiyama, K. S. Ikeda, and Y. Shimizu,Phys. Rev. E , 036112 (2009).[41] T. Okushima, T. Niiyama, K. S. Ikeda, and Y. Shimizu,Phys. Rev. E , 036109 (2007).[42] T. Okushima, T. Niiyama, K. S. Ikeda, and Y. Shimizu, Phys. Rev. E , 021301(R) (2018).[43] T. Okushima, T. Niiyama, K. S. Ikeda, and Y. Shimizu,Phys. Rev.E , 032304 (2018).[44] A. J. F. Siegert, Phys. Rev. , 617 (1951).[45] J. Keilson, J. Appl. Prob. , 247 (1964).[46] J. Keilson, J. Appl. Prob. , 405 (1965).[47] D. A. Darling and A. J. F. Siegert, Ann. Math. Stat. ,624 (1953).[48] Z. Schuss and B. J. Matkowsky, SIAM J. Appl. Math. , 604 (1979).[49] B. J. Matkowsky and Z. Schuss, SIAM J. Appl. Math. , 242 (1981).[50] N.G. van Kampen, Stochastic Processes in Physics andChemistry, 3rd ed. , (Elsevier, Amsterdam, 2007).[51] D. Hartich and A. Godec, J. Stat. Mech. , 024002(2019).[52] D. Aldous,
Probability Approximations via the PoissonClumping Heuristic (Springer-Verlag, New York, 1989).[53] G.H. Golub and C.F. Van Loan,
Matrix Computations ,4th ed. (Johns Hopkins University Press, Baltimore, MD,2012).[54] E.W. Montroll, K.E. Shuler, Adv. Chem. Phys. , 12787 (2013).[56] W. Dai, A. M. Sengupta, and R. M. Levy, Phys. Rev.Lett. , 048101 (2015).[57] Z. Q. Lee, W. Hsu, M. Lin, PLoS ONE , e93348 (2014).[58] V. Tejedor, O. B´enichou, and R. Voituriez, Phys. Rev. E , 065104(R) (2009).[59] Y. Kimura, Y. Saito,T. Nakada, and C. Kaito, PhysicaE , 11 (2002).[60] M. M¨uller and K. Albe, Acta Mater. , 3237 (2007).[61] T. Shibata, B. A. Bunker, Z. Zhang,D. Meisel, C. F.Vardeman II, and J. D. Gezelter, J. Am. Chem. Soc. , 11989 (2002).[62] H. Haken, Synergetics, An Introduction: NonequilibriumPhase Transitions and Self-Organization in Physics,Chemistry, and Biology , 3rd rev. enl. ed. (Springer,Berlin, 1983).[63] M. P. Tosi, F. G. Fumi, J. Phys. Chem. Solids25