From coherent quasi-irreversible quantum dynamics towards the second law of thermodynamics: The model boron rotor B13+

The planar boron cluster B13+ provides a model to investigate the microscopic origin of the second law of thermodynamics in a small system. It is a molecular rotor with an inner wheel that rotates in an outer bearing. The cyclic reaction path of B13+ passes along thirty equivalent global minimum structures (GMi, i = 1, 2, ..., 30). The GMs are embedded in a cyclic thirty-well potential. They are separated by thirty equivalent transition states with potential barrier Vb. If the boron rotor B13+ is prepared initially in one of the thirty GMs, with energy below Vb, then it tunnels sequentially to its nearest, next-nearest etc. neighbors (520 fs per step) such that all the other GMs get populated. As a consequence, the entropy of occupying the GMs takes about 6 ps to increases from zero to a value close to the maximum value for equi-distribution. Perfect recurrences are practically not observable.

The importance of molecular rotors has been recognized by the 2016 Nobel Prize in Chemistry awarded to B. L. Feringa, J.-P. Sauvage and J. F. Stoddart. [1][2][3] Fundamental developments have discovered fascinating phenomena and enabled a wealth of applications in the gas phase, [4][5][6][7][8] in solution, 4,[9][10][11] in crystals, 4,12 on surfaces, [4][5][6][13][14][15] and in biosystems. 16 Here we use a molecular model rotor to demonstrate the way from molecular quantum dynamics to explicitly time dependent statistical thermodynamics. Specifically, we discover the relation between the tunneling quantum dynamics of a system with cyclic multi-well potential and the time evolution of the entropy of occupying the wells, from zero to close to the maximum value for equidistribution. Our model system is the planar molecular boron rotor B 13 +17-20 which has a trilateral inner wheel and a ring-shaped outer bearing with N w = 3 and N b = 10 boron atoms, respectively. One cycle of the wheel rotating with respect to the bearing corresponds to motion along a cyclic potential energy surface with N w ×N b = 30 equivalent wells. The potential wells support thirty equivalent global minimum structures that are arranged in cyclic order, GM 1 , GM 2 , ..., GM 30 . These are separated from each other by thirty equivalent transition states, again in cyclic order, TS 1,2 , TS 2,3 , ..., TS 30,1 , and all with the same low potential barrier (V b ≈ 106.4 h c cm -1 ).
The choice of B 13 + among similar candidates, B 11 -, 21,22 B 15 + , 23 B 19 -24,25 (for systematic investigations and reviews of boron clusters and rotors, see , is motivated by the fact that it is the first planar boron rotor whose structure and fluxionality have been probed by means of infrared (IR) photodissociation spectroscopy. 19 Experiment and theory suggest that the rotor is prepared in one of the global minimum structures. [17][18][19][20] Here we consider this type of preparation of B 13 + , with energy E below the barrier V b , and with zero angular momentum. Preparation of a single GM means that the entropy for occupying the GMs is equal to zero. Since this initial state is not an eigenstate, B 13 + will tunnel out of the potential well that supports the initial GM. This work aims to discover the mechanism of tunneling in the cyclic multi-well potential and the consequence for the statistical thermodynamics. There are many questions such as: Will the tunneling proceed sequentially from the initial GM to its nearest neighbors, then to the next nearest ones, etc.? Will the tunneling in the cyclic multi-well potential exhibit new phenomena that are fundamentally different from the familiar tunneling in cyclic double-well potentials observed in examples such as B 2 Cl 2 F 2 29 and T 2 captured in a metal-organic framework with active Cu(I) sites. 30 Will the quantum dynamics induce the increase of the entropy for occupying the GMs? On what time scale? These are of course rather general questions; the model system B 13 + just serves as an example which has been extensively investigated experimentally 19 and theoretically. [17][18][19][20] The results will be obtained by means of quantum molecular dynamics simulations, using the model of Ref. 20. The model parameters (e.g. V b ) are obtained by quantum chemistry calculations at the PBE0/def2-TZVP level, 31,32 using the Gaussian 09 suite of programs. 33 Before starting, let us recall some important properties of the model B 13 + : All the equivalent GMs have the same properties such as the same symmetry (C 2v ) and the same normal mode frequencies ν 1 , ν 2 , ..., ν 33 . Likewise, all the equivalent TSs also have the same properties such as the same barrier height V b (see above), the same symmetry (again C 2v ) and the same normal mode frequencies i ν 1 , ν 2 , ..., ν 33 . Motions along the normal mode with the lowest frequency ν 1 of the GMs or with the imaginary frequency i ν 1 of the TSs correlate with motions along the reaction path, corresponding to the rotation of the wheel with respect to its bearing. The rotational angle ϕ represents the corresponding reaction coordinate. Actually when the nuclei of the wheel rotate with respect to the bearing, all the nuclei of the bearing move concertedly along equivalent small periodic orbits such that the shape of the bearing adapts to the instantaneous orientation of the wheel, as illustrated in Fig. 1. The shape of the bearing thus rotates synchronously with the rotating wheel, without any direct rotation of the bearing. This is called "pseudo-rotation", analogous to rotating molecules in pseudo-rotating cages. 34 The low barrier of the TSs together with the pseudo-rotation of the bearing support nearly free rotations of the wheel along the angle ϕ.
The values of ν 1 =135.2 c cm 1 and ν 1 = 134.8 c cm −1 (or ω 1 =2πν 1 and ω 1 = 2πν 1 ) are almost equal to each other. Likewise, all other normal mode frequencies ν 2 , ν 3 ,...,ν 33 of the GMs correlate closely with the ν 2 , ν 3 , ..., ν 33 at the TSs. This suggests a simple one-dimensional model, i.e. the rotation/pseudo-rotation along ϕ is essentially decoupled from the other vibrations. Marginal effects of small couplings may be considered in a refined model. The corresponding model Hamiltonian is Tunneling dynamics of the model boron rotor B 13 + with inner triangular wheel rotating in outer pseudo-rotating bearing (schematic). Clockwise and anti-clockwise tunneling from the global minimum GM 16 to the neighbors GM 15 and GM 17 , respectively (left), and then sequentially to GM 14 , GM 13 , GM 12 , GM 11 , etc, and simultaneously also to GM 18  of inertia (see below). The cyclic thirty-well potential is modeled as ..,30 that range from one transition state TS i-1,i to the next one TS i,i+1 . The potential well in domain d i supports GM i . The force constants at the potential minima and maxima have the same absolute values, k = V b ×30 2 = I eff ·ω 2 . This is in accord with the near equality of the absolute values of the normal mode frequencies ω 1 and ω 1 (see above). This suggests the estimate of the effective moment of inertia as the mean value at GM and at TS, The result can be rationalized as the (approximate) sum of two contributions for the moments of inertia for the rotation of the wheel in its bearing and for the pseudo-rotation of the bearing (see supplementary material).
The rotational/pseudo-rotational eigenstates of the model B 13 + are obtained as solutions of the Schrödinger equation H m ψ m (ϕ) = E m ψ m (ϕ). The wavefunctions ψ m (ϕ) and their levels are illustrated in the supplementary material. The lowest 30 levels form a rather narrow band, from E 0 = 45.24 hc cm -1 to E 29 = 67.52 hc cm -1 , all below the barrier V b . The wave function that represents the initial preparation of a single GM, for example GM 16 , can be written as superposition of the eigenfunctions ψ m (ϕ) with expansion coefficients c m . This wavepacket evolves as The resulting density is ρ(ϕ, t) = |ψ(ϕ, t)| 2 . It is carried by the rotational/pseudo-rotational flux 35 Preparation of the initial GM 16 without any angular momentum implies the boundary condition j(ϕ b , t) = 0 at the bottom (ϕ b = 0.5 (ϕ 15,16 +ϕ 16,17 )) of its potential well. 29,35,36 This causes equivalent clockwise and anti-clockwise branches of the flux along -ϕ and +ϕ, respectively. When the rotor tunnels out of the potential well that supports the initial GM 16 , it will populate other GMs. This will increase the entropy for occupying the GMs. In general, the probability of observing GM i is obtained by integrating the density in the domain d i , The corresponding entropy of occupying the GMs (in units of the Boltzmann constant k B ) is For the ideal preparation of a single GM, the initial entropy is equal to zero. The maximum value of the entropy of occupation of the GMs is S max (t)=k B ln N m =k B ln 30. This applies to the hypothetical case of equal populations of all GM i , P i (t) = 1/30. In general, the value of the entropy, 0 ≤ S(t) ≤ k B ln N m may be considered as a measure of the degree of delocalization of the rotor in various GMs. Smaller or larger values of S(t) indicate better localization in few GMs, or more equal populations of all GMs, respectively. According to the second law of thermodynamics, the entropy should increase from zero to the maximum value for occupying all GMs equally. But this second law does not tell anything about the time scale. The present model system allows to discover this missing link between quantum dynamics and statistical mechanics, in the present case for the entropy of occupying the GMs, eq. (4). Alternative definitions are considered in the supplementary material.
The key result of this investigation is condensed in Fig. 2. Accordingly, the boron rotor B 13 + is prepared initially in one of its global minimum structures (e.g. GM 16 ), with energy E below the barrier V b , and with zero angular momentum. Note that the nuclear spins prohibit the preparation of this initial state by laser pulses. 37 The tunneling dynamics in its thirty-well potential proceeds such that the entropy S(t) of occupation of the GMs increases rapidly. Within approximately 6 ps, S(t) increases from its very small initial value to a rather large value close to the maximum. . a, b, Coherent quasi-irreversible tunneling of the boron rotor B 13 + in its cyclic thirty-well potential, starting from a single global minimum structure, e.g. GM 16 , with energy below the barrier, and zero angular momentum. c, d, Coherent periodic tunneling of a system with cyclic double-well potential. 30 For reference, the horizontal lines mark the values of the entropies ln2, ln3,..., ln30 for equal populations of 2,3,...,30 GMs, respectively.
The ideal initial value S(t=0) is equal to 0, slightly larger values are also possible if the initial localization in one of the GMs is not 100% perfect -this case arises because the 30 eigenfunctions ψ m (ϕ) of the tunneling band in the expansion (1) constitute an excellent but not perfectly complete set of basis functions. The rather large value of S(t) which is reached after ca 6 ps is quite robust, except for rather small oscillations. Systematic search during 50 ns yields the mean value S/k B = 3.017 = ln 20.4, and a very sharp distribution of the entropy centered at this value, with variance ∆ S/k B = 0.101. The mean value S/k B is slightly below the maximum value ln 30 = 3.401. This points to incomplete equilibration; analogous incomplete explorations of molecular phase space have been discovered by Schlag and Levine. 38 Exceptional deviations from the mean value S/k B are possible, but they are extremely rare. For example, the smallest value during 50 ns after quasi-equilibration is Min(S/k B ) = 1.303 = ln 3.7 at t = 41241.2 ps, but the probability of this rare event is only 1.9×10 -5 , see supplementary material. No perfect recurrence of S(t) back to the very small initial value occurs within any practical limits of observing the system. This is coherent quasi-irreversible tunneling in the cyclic multi-well potential. The results shown in Figs. 2a and 2b document the quantum dynamical results for a small system that reflect the trends of the second law of thermodynamics.
The apparent absence of recurrences of B 13 + from the initial global minimum structures (for example GM 16 ) back to the same GM during tunneling in the multi-well potential can be explained as follows: Essentially, the first hypothetical time T of such recurrences would correspond to identical phase factors e −i E m T/ in the expansion (1). The energies E m are incommensurable, i.e. this condition will never be satisfied. Even if one allows sufficiently small deviations from perfectly identical phases modulo 2π it would call for times T beyond all practical limits.
Coherent quasi-irreversible tunneling in cyclic thirty-well potential of B 13 + is in marked contrast with the familiar coherent tunneling in a cyclic periodic double-well potential such as what is used for the model T 2 captured in a metal-organic framework with active Cu(I) sites. 30  (as discussed above) and the bottom panels of Fig. 2, respectively. In the case of the cyclic double-well potential, the ideal time dependent probabilities of occupying the initial (GM 1 ) and the other (GM 2 ) global minimum structures are known analytically, 35 P 1 (t) = cos 2 (ω t), P 2 (t) = sin 2 (ω t) where the tunneling frequency ω is defined by the tunneling splitting, ∆E = ω. Their periodicity implies that the entropy S(t)/k B = P 1 (t) ln P 1 (t) P 2 (t) ln P 2 (t) oscillates periodically between its minimum and maximum (k B ln 2). Recurrences to occupations of a single GM are also observed periodically, with period T = π/ω (Figs. 2c-2d).
The entirely different coherent quasi-irreversible tunneling of the model systems B 13 + in the thirty-well potential versus coherent periodic tunneling of T 2 in double-well potential can also be rationalized by the different time evolutions of the densities ρ(ϕ, t) and the fluxes j(ϕ, t), starting from the initial occupation of a single global minimum structure. This is illustrated in the left and right panels of Fig. 3 for the examples of the initial preparations of GM 16 (B 13 + ) and GM 1 (T 2 ), 30 respectively. Apparently, the boron rotor B 13 + (left) tunnels from the initial GM 16 in equivalent clockwise (along -ϕ) and anti-clockwise (along +ϕ) manners sequentially to its "left" (GM 15 ) and "right" (GM 17 ) neighbors, then to the next nearest neighbors (GM 14 and GM 18 ), etc. see Fig. 1. Each step takes about 520 femtoseconds. The density of B 13 + which was localized initially in a single GM thus penetrates into the domains of all GMs. The two equivalent clockwise and anti-clockwise delocalizations of the density (Fig. 3a) are carried by the corresponding fluxes (Fig. 3c). Sequential fluxes into new domains are accompanied with fluxes back to previous domains. This effect of bi-directional tunneling supports continuous delocalization of the wave function. Delocalization is also supported by the effect of interfering partial waves, or lobes of the wavefunction. A prominent example is observed at ca 8 ps when the frontier lobes of the two branches of the wavefunction that tunnel in clockwise and anti-clockwise directions "collide" in approximately the domains of GM 30 , GM 1 , GM 2 , opposite to the initial domain of GM 16 . Subsequently, this collision causes rather complex interference patterns of the wave function that tend to populate all GMs: this contributes to the quasi-irreversible delocalization of the density in all GMs and, therefore to the rather large mean value of the entropy S(t) of occupying the GMs, which is recalled in the time domain of ca 6 ps, as documented in Fig. 2. Close inspection of the initial rise of S(t) from (almost) zero to the large mean value reveals an overall smooth behavior which is superimposed by tiny steps, or shoulders: These correlate with the step-wise penetration of the density from the initial GM to its nearest, next-nearest etc. neighbors (ca 520 fs per step). For comparison, the density and flux for coherent periodic tunneling in the cyclic double-well potential 30   In conclusion, the present model B 13 + offers a rather simple demonstration of the challenging phenomenon of coherent quasi-irreversible quantum molecular dynamics, analogous to the second law of thermodynamics, and even beyond: The law says that the entropy should increase, but does not tell us on which time scale. In the present example, quantum molecular dynamics yields the missing information, i.e. the entropy of occupying the GMs increases from zero to the mean value S/k B within ca 6 ps. This value can be rationalized by the underlying tunneling mechanism i.e. step-bystep penetration into neighboring, next-neighboring etc GMs, with ca 520 fs per step. It is amazing that the emergence of the second law with respect to the entropy of occupying equivalent GMs is observed already for a rather small isolated system such as the boron rotor B 13 + , quite different from small open quantum systems which are coupled to large heat baths or environments; for a review, see Ref. 39. This should stimulate extended investigations, e.g. search for analogous coherent quasiirreversible tunneling in different molecules or molecular clusters with multi-well potentials that support many equivalent global minimum structures. For the present system, B 13 + , and likewise for analogous planar boron rotors, we suggest to monitor the time evolution for the mean value of the dipole moment, as a signature of the coherent quasi-irreversible tunneling. Specifically, the initial dipole should be "washed-out" due to the nearly perfect delocalization to all other GMs, on the same time scale of ca 6 ps for the increase of the entropy of population of the GM, cf.