Quantum equilibration of the double-proton transfer in a model system porphine †

There is a renewed interest in the derivation of statistical mechanics from the dynamics of closed quantum systems. A central part of this program is to understand how closed quantum systems, i.e. , in the absence of a thermal bath, initialized far-from-equilibrium can share a dynamics that is typical to the relaxation towards thermal equilibrium. Equilibration dynamics has been traditionally studied with a focus on the so-called quenches of large-scale many-body systems. We consider here the equilibration of a two-dimensional molecular model system describing the double proton transfer reaction in porphine. Using numerical simulations, we show that equilibration indeed takes place very rapidly ( B 200 fs) for initial states induced by pump–dump laser pulse control with energies well above the synchronous barrier. The resulting equilibration state is characterized by a strong delocalization of the probability density of the protons that can be explained, mechanistically, as the result of (i) an initial state consisting of a large superposition of vibrational states, and (ii) the presence of a very effective dephasing mechanism.


Introduction
There is currently a renewed interest in the derivation of statistical mechanics from the dynamics of a closed quantum system. 1 In this approach, instead of assuming a priori that the system is in some mixed state, such as e.g. a micro-canonical ensemble, one describes it at all times using a pure state. One then seeks to show that, under reasonable conditions, the system behaves as if it was described by a statistical ensemble.
In this way the use of statistical mechanics can be justified without introducing additional external degrees of freedom, such as e.g. thermal ''baths''.
A central part of this programme has been to understand the process of equilibration, i.e., how a constantly-evolving closed quantum system can share aspects typical to the relaxation towards a static equilibrium. The main insight relies on the fact 2-4 that, if measurements are limited to small subsystems or restricted sets of observables, then ''typical'' pure states of large quantum systems are essentially indistinguishable from thermal states. It can then be shown 5,6 that under very general conditions on the Hamiltonian and for nearly all initial states, the system will eventually equilibrate, in the sense that an (again, restricted) set of relevant physical quantities will remain most of the time very close to fixed, ''equilibrium'' values.
Quantum equilibration has been recently put within reach of experimental verification. Fueled by enormous improvements in experimental techniques it is now feasible to control quantum systems with many degrees of freedom. This is particularly true for the development of techniques to cool and trap ultracold atoms with optical lattices (giving raise to effective lattice systems) [7][8][9][10][11] or suitable confinement potentials (that lead to effective low-dimensional continuous systems). [12][13][14][15] Similarly, systems of trapped ions 16,17 allow us to precisely study the physics of interacting systems in the laboratory. [18][19][20][21][22][23] In such highly controlled settings, equilibration and thermalisation dynamics has been studied both experimentally [24][25][26][27][28] and numerically, often with a focus on the so-called quenches, i.e. rapid changes of the Hamiltonian. [29][30][31][32][33][34][35][36][37] While the most common notion of equilibration applies to thermodynamically large systems, we are here interested on the more general notion of equilibration on average 1 (viz., a timedependent property equilibrates on average if its value is for most times during the evolution close to some equilibrium value). This notion of equilibration applies also to small systems and we use it here to characterize the dynamics of an isolated molecular system. We consider a model system describing the double proton transfer reaction in the electronic ground state of porphine, a paradigmatic system in which the making and breaking of H-bonds occurs in an anharmonic Born-Oppenheimer potential-energy surface. Despite the approximations involved, the two dimensional analytical potential energy surface describing the coupled linear motions of the two protons serves well for the proof of principle discussed here. As it will be shown, for initial states with energies above the synchronous double-proton transfer energy barrier involving a large coherent sum of vibrational states, the isomerization dynamics can reach a long-lived quasi-stationary regime with associated physical quantities, such as the mean position or momentum, that remain close to fixed, ''equilibrium'', values.

Definitions
Equilibration is achieved when a measurable quantity, after having been initialised at a non-equilibrium value, evolves towards some value and then stays close to it for an extended amount of time. More precisely, given some observable A and a system of finite but arbitrarily large size, if its expectation value hÂ(t)i equilibrates, then it must do so around the infinite time average, 1 A ¼ lim If the infinite-time average fluctuation of hÂ(t)i around Ā is small, then we say that the observable A equilibrates. This notion of equilibration is compatible with the recurrent and time reversal invariant nature of unitary quantum dynamics in finite dimensional systems. To see that, let us assume a closed system whose state is described by a vector in a Hilbert space of dimension d T and whose Hamiltonian has a spectral representationĤ ¼ where E k are its energies and |E k i are eigenstates of Ĥ . ‡ Choosing units such that h = 1, the state at time t is then given by jcðtÞi ¼ P k c k e ÀiE k t E k j i with c k hE k |c(0)i, and the expectation value of Â can be written as: The long-time average in (1) is then: where we have used that d k; if the system relaxes at all, it must be to the following density matrix: As noted in ref. 29, the result in (4) can be thought as stating the prediction of a ''diagonal ensemble'', |c k | 2 being the corresponding weights. For a non-integrable system, this ensemble is generally expected to reduce to a thermal ensemble. Sufficient conditions for equilibration in the sense defined above can be then defined in terms of the so-called effective dimension. Roughly speaking, the effective dimension is a measure of the number of significantly occupied energy eigenstates, and is defined as: In ref. 5, 6 and 38 it is shown that for any Hamiltonian with non-degenerate gaps, § a large effective dimension is sufficient to guarantee that equilibration will be attained. Note that if the initial state is taken to be an energy eigenstate, the resulting effective dimension is one, while the one resulting from a uniform coherent superposition of d T energy eigenstates to different energies is d T .
Provided that a given observable equilibrates, it was recently shown 40 that the equilibration time can be related with a dephasing time t d as: where s G is the standard deviation of the energy gaps, G a = E j À E i , when weighted by their respective relevances q a , i.e.: where with v a = c j *A ji c i /s A , and A ij := hE i |A|E j i are the matrix elements of A in the energy eigenbasis. The denominator s A = a max À a min is the range of possible outcomes, being a max(min) the largest (smallest) occupied eigenvalue of Â. ‡ Note that the sum runs over d E r d T terms, since some eigenspaces can be degenerate. If the Hamiltonian has degenerate energies, we choose an eigenbasis of Ĥ such that the initial state, |c(0)i, has non-zero overlap with only one eigenstate |E k i for each distinct energy. § Concerning the assumption of the Hamiltonian having non-degenerate gaps, it is shown in ref. 39 that as long as there are not exponentially many degeneracies the argument stays the same.

Porphine model
The above introduced concepts are now used to characterize the non-integrable model Porphine designed by Smedarchina et al. 41 to describe the switch from synchronous (or concerted) to sequential (or stepwise) double-proton transfer. 42,43 This model accounts for the motions of two protons (labeled 1 and 2) along coordinates R 1 and R 2 , respectively, from the domains of the reactant (R) to the product (P) (see Fig. 1). The PES model is 41 where G is the coupling parameter which formally represents the interaction between the two hydrogen bonds, 2D 0 is the transfer distance of each proton, i.e., the distance between its two equilibrium positions, and U 0 is the barrier height for single-proton transfer. The parameter U 0 = 0.473 eV has been fitted in ref. The reaction can lead from the reactant R via alternative transitions states TS to the intermediates (I), and subsequently via the other two TS to the product P. In addition, Fig. 1 shows a central maximum (that corresponds to a saddle point of second order) labeled SP2. The competing synchronous reaction mechanism leads from the reactant R via SP2 to the product P. The model potential in eqn (8) is symmetric with respect to the diagonals R 1 = AER 2 . It accommodates nearly degenerate doublets of eigenstates C v+ (R 1 ,R 2 ) and C vÀ (R 1 ,R 2 ), with energies below the barriers TS, plus higher excited states. We then chose our initial state to be of the general form: where DR represents a shift along the diagonal R 1 = R 2 and ð Þis a superposition state that represents the localized ground state wavefunction of the reactant, with C 0+ (R 1 ,R 2 ) and C 0À (R 1 ,R 2 ) representing the lowest doublet (n = 0), i.e., HC nAE (R 1 , where H is the full Hamiltonian of the system, i.e.: with @ 2 @R k 2 the kinetic energies of the two protons, M = m p the proton mass, and V(R 1 ,R 2 ) the scalar potential given in eqn (8).
Due to the symmetry of the model potential in eqn (8), the localized ground state wavefunction of the reactant is equivalent to the localized ground state wavefunction of the product. Therefore, the initial states in eqn (9) are, in fact, representative of two equivalent configurations where the two protons are placed along the reactant-product diagonal (see Fig. 2).
The localization, eqn (9), of the reactant wavefunction may be generated by symmetry breaking, for example, by selective deuteration of the C 4 N rings for the reactant; this would not change the PES, but it would induce localization of the reactant wavefunction due to the slightly heavier mass that is associated with the reactant compared to the product. 42 Also, shifts of initial wavefunctions from equilibrium to non-equilibrium positions may be induced, for example, by means of pumpdump laser pulse control, 48 as designed by Tannor and Rice. [49][50][51] Essentially, the ultrashort pump pulse transfers the molecule from the electronic ground state to an excited state. Here, the system evolves from the original configuration until  it is shifted to the target position. Finally, the dump pulse brings the wavepacket back to the electronic ground state, thus preparing the initial state for the subsequent reaction. Analogous shifts of the original wavefunctions to nonequilibrium positions have been demonstrated by means of laser pulse control, by Kapteyn, Murnane, and co-workers. 52

Equilibration time
To gain some insight into the behaviour of the effective dimension and the equilibration time in porphine, we consider the case of initial states defined in eqn (9) and the type of projective position measurements: In eqn (11), are identity operators in the position representation acting, respectively, in the subspaces defined by R 1 and R 2 . Whereas the position operator in eqn (11) involves an infinite Hilbert space, the relevant energy subspace is finite in practice (for numerical purposes) and hence it is possible to apply the machinery introduced above. We evaluate the effective dimension d eff in eqn (5) and the equilibration time T eq in (6) for different values of DR A [À2.7,2.7] (see Fig. 3). For that, we first computed v a = c j *R ji c i /s R , where R ij := hE i |R|E j i and s R = R max À R min has been chosen to be s R = 3s gs , and s gs ¼ 1 ffiffiffiffiffi 25 p is the dispersion of the localized ground state wavefunction of the reactant.
In addition to d eff , which represents an equilibration criterion, and T eq , which provides an estimate of the velocity of the equilibration process, we found convenient to show also these two concepts in a compact way by introducing an effective equilibration time (or effective equilibration speed): Note that, by definition, even if T eq is small, T eq can still be close to infinite for initial states being eigenstates of the observable of interest (i.e., d eff B 1).
The resulting values of d eff , T eq and T eq are plotted in Fig. 3. We identified three different regions: (i) in green, regions where equilibration is expected to occur fast, (ii) in blue, regions where equilibration is expected to occur slowly, and (iii) in red, regions where equilibration is not expected to occur. Interestingly, region (iii) corresponds to the regime where our initial state can be written as the superposition of a very small number of energy eigenstates. Note that for DR E 0 then C 0 is just the sum of two energy eigenstates C 0+ and C 0À . The number of energy eigenstates required to define the initial wavefunction increases as one departs from the global minima. Regions of slow equilibration (in blue) are the result of a nontrivial interplay between T eq and d eff . Note, for example, the relative minimum shown by T eq at DR = 0 that corresponds to the situation where the initial wavefunction would fall equally to the left and right wells. In regions where equilibration is expected to occur fast (in green), the equilibration time T eq decreases almost monotonically. This is due to the fact that the barrier leads to a very effective dephasing mechanism at energies close or above the global maximum.
It is important to emphasize that the quantities d eff , T eq and T eq are to be understood only as qualitative descriptors of equilibration, as they are defined using rather heuristic arguments. Therefore, as it will be shown in the next section, these quantities do not provide an accurate quantitative estimation of the equilibration time, but do provide us with a good guess.

Equilibration dynamics
We now aim to study the dynamics that leads to equilibration for initial states induced by pump-dump laser pulse control of  (9) and different values of DR for the type of measurement defined in eqn (11) using a square grid of size 151a 0 with a grid spacing of 0.08a 0 . The Born-Oppenheimer potential-energy surface is also plot in solid black line for reference. Three different regions are identified in the bottom panel: (i) in green, regions where equilibration is expected to occur rapidly, (ii) in blue, regions where equilibration is expected to be attained slowly, and (iii) in red, regions where equilibration is not expected to occur. the form in eqn (9). Specifically, we choose DR = À2.2a 0 such that the initial state lies in the rapid equilibration regime (see Fig. 3). The mean energy of this initial state is rather high (4.885 eV) and well above the TS and SP2 barriers (respectively 0.600 eV and 1.069 eV).
Starting with C 0 (R 1 ,R 2 ;À2.2a 0 ) (see Fig. 4a), the initial synchronous mechanism of the first forward reaction is characterized by a rapid dispersion of the wavepacket and relief reflections of the broadened wavepacket from wide regions of the steep repulsive wall of the PES close to the minimum of the product that leads to the switch from the sequential to the synchronous mechanism (see Fig. 4b). This situation is reminiscent of the near field interference effect arising when periodic diffracting structures are illuminated by highly coherent light or particle beams. 53 The time dilation supported by continuous wavepacket dispersion leads to a strong proton delocalization. An apparently ''chaotic'' flux is however coherently directing the recovery of the concerted double proton transfer at t \ 100 fs. Due to the dephasing between partial waves, the grid structure of the probability density associated to the sequential double proton transfer progressively dilutes into what is reminiscent of a stationary state, showing a series of minima disposed perpendicular to the diagonal R 1 = R 2 (see Fig. 4c). The long-lived quasi-stationary state, formed already at B150 fs, lasts beyond 500 fs (see Fig. 4d) and it gives rise to the equilibration of the position operator in eqn (11) as well as of the momentum operator P = |P 1 ihP 1 | # I 2 + I 1 # |P 2 ihP 2 |. The full evolution dynamics can be found in S1 (ESI †).
The above dynamics result into an effective equilibration of the phase-space variables. In Fig. 5 (top panel) we show the differences % R À hR(t)i and % P À hP(t)i as a function of time for the initial state C 0 (R 1 ,R 2 ;À2.2a 0 ). Equilibration is achieved once these differences become small ({1) and stay small for a lapse of time that is comparable to the relevant time-scale of the dynamics. 54 The equilibration of the position and momentum yields a phase-space portrait (see the inset in Fig. 5 top panel) that is reminiscent of the phase-space portrait of a damped harmonic oscillator. Interestingly, this phase-space dynamics is attained here without any loss (nor gain) of energy.
In Fig. 5 (middle and bottom panels) we show the dynamics of % R À hR(t)i and % P À hP(t)i for another two different initial states, viz., C 0 (R 1 ,R 2 ;À1.9a 0 ) and C 0 (R 1 ,R 2 ;À1.75a 0 ), with energies 1.905 eV and 1.185 eV (i.e., B0.8 eV and B0.1 eV above the SP2 barrier) respectively. As we decrease DR and thus approach the equilibrium position DR = 0, the energy of the initial state gets closer to the top of the synchronous barrier. This results into a dynamics that departs more and more from an equilibration process. Note, for instance, that while the phase-space portrait in Fig. 5 (mid panel) is still reminiscent of some Fig. 4 Dynamics of the two-proton density |C(R 1 ,R 2 ,t)| 2 at four different time snapshots for the initial state C 0 (R 1 ,R 2 ;À2.2a 0 ). The initially welllocalized nuclear density gives way later to a strong proton delocalization along the synchronous pathway that is maintained for at least 1 ps. damped dynamics, the phase-space diagram in Fig. 5 (bottom panel) is clearly not, as it does not converge to equilibrium. For even smaller DR values, the initial states fall down into the minimum of the well, where equilibration is no longer expected to occur. The full dynamics for the initial states C 0 (R 1 ,R 2 ;À2.2a 0 ), C 0 (R 1 ,R 2 ;À1.9a 0 ), and C 0 (R 1 ,R 2 ;À1.75a 0 ) can be found in S1, S2, and S3 (ESI †) respectively.

Discussion
One might think that provided there are at least two degrees of freedom with non-negligible coupling between each other, one can expect one to act as a bath for the other and vice versa, both leading to the damping of the mutual dynamics. As a result, one may encounter that observables reach some equilibrium value. While this hypothesis might sound reasonable, it is incorrect. First of all, two particles with identical masses (e.g., two protons) can hardly act one as a bath for the other. Second, even a system with a single degree of freedom can equilibrate. See for example the works in ref. 54 and 55. In order a quantum system to equilibrate there is no need for a bath. That is the main point that differentiates equilibration from thermalization. Also, working with initial conditions far-from-equilibrium within anharmonic regions of the potentialenergy surface does not guarantee equilibration. Otherwise equilibration should be observed for almost every molecule out-of-equilibrium. This is definitely not the case.
To gain more insight into the physical mechanism that yields equilibration, we here take advantage of the interpretation value of Bohmian mechanics. [56][57][58][59][60][61] It is possible to reinterpret the quantum formalism as describing particles following definite trajectories, each with a precisely defined position at each instant in time. In this interpretation, called Bohmian mechanics, the trajectories of the particles are guided by the wavefunction. This allows for phenomena such as double-slit interference, as has been investigated experimentally for single photons. [62][63][64] Note that this is different from semi-classical trajectories, 65-67 other generalized hydrodynamic approaches, 68 and also from the Feynman path formalism of quantum mechanics. 69 For example, in contrast to the Feynman formalism, Bohmian mechanics says that each quantum particle in a given experiment follows a trajectory in a deterministic manner. Thus, much of the intuition of classical mechanics is regained. 70 In any case, let us insist that we use an ensemble of exact Bohmian trajectories (evaluated from the wavefunction solution of the timedependent Schrödinger equation) with the only intention to provide a different viewpoint on the process of equilibration that is more ''digestible'' for an adherent of classical mechanics and that it is at the same time compatible with the expectation values associated to the exact quantum dynamics.
We consider an ensemble of Bohmian trajectories {R a (t)} = {R a 1 (t),R a 2 (t)} initially sampled from the probability distribution P(R 1 ,R 2 ) = |C 0 (R 1 ,R 2 ,À2.2a 0 )| 2 and whose time evolution is defined as where In Fig. 6 we show the time evolution of an ensemble of five hundred Bohmian trajectories R a 1 (t) sampled from |C 0 (R 1 ,R 2 ,À2.2a 0 )| 2 . An initially orchestrated oscillatory dynamics with an amplitude that extends over the two wells is rapidly (at around B50 fs) substituted by apparently incoherent oscillations with a very small amplitude. Not much later, at around 150 fs, the trajectories fall into a dynamics that is highly localized and homogeneously distributed. Of particular interest is the fine structure of the Bohmian trajectories, characterized by high frequency components and rapid shifts in their direction. This quantum mechanical motion arises from strong timedependent quantum force fields (originating from quantum interferences) which essentially dominate the dynamics and are responsible for shifting the direction of trajectories away from that dictated by the classical force field.
As already anticipated in Section (4), the above dynamics can be seen as the result of a very efficient dephasing mechanism that leads to a diagonal ensemble. The barrier in between the two wells induces a continuous time dilatation supported by wavepacket dispersion that yields strong interference effects between different parts of the wavepacket. Eventually, these interferences are responsible for freezing the transfer of the two protons and lead to the localization of the corresponding trajectories within the concerted pathway (see Fig. 6). This equilibration dynamics corresponds to the picture where the two protons become localized at certain configurations along the concerted transfer reaction path, as if they were associated to a stationary state.
In view of the above results, we may conclude that there are two effects that must coexist to find equilibration: (i) Initial states consisting of a large superposition of energy eigenstates.
(ii) The presence of a very effective dephasing mechanism. Note that these two conditions are not easily met in general, but can be found for the porphine molecule and probably also for other related molecules, under certain initial conditions. More specifically, for the porphine model considered in this work, condition (i) is met for initial states with energies well above the SP2 barrier, and the dephasing mechanism in condition (ii) is originated by the SP2 barrier. Concerning condition (i), there are other types of initial states, different from those defined in eqn (9), that can be considered as potentially leading to the equilibration of the canonical variables R and P. Very general initial conditions can be defined, for example, by introducing the wavefunction C G (R 1 + DR 1 ,R 2 + DR 2 ), where C G (R 1 ,R 2 ) represents a Gaussian approximation to the localized ground state wavefunction of the reactant, with a dispersion s gs ¼ 1 ffiffiffiffiffi 25 p and centered at the origin of coordinates. In particular, for DR 1 = ÀDR 2 , the wavepacket C G (R 1 + DR 1 ,R 2 + DR 2 ) is centered along the diagonal R 1 = ÀR 2 that leads from one Intermediate to the other (i.e., perpendicular to the diagonal R 1 = R 2 that leads from the reactant to the product). Due to the symmetry of the model potential in eqn (8), the state C G (R 1 + DR 1 ,R 2 À DR 1 ) around the Intermediate minima, may correspond to four equivalent hydrogen atoms distributions (see Fig. 7). Starting from the initial state C G (R 1 + 2.2a 0 ,R 2 À 2.2a 0 ), the resulting dynamics is analogous to that shown in Fig. 4 but the probability density that corresponds to Fig. 4c and d extends now along the diagonal R 1 = ÀR 2 instead of R 1 = R 2 and there is a notable increase of the probability density along the sequential paths. The associated behaviour of hRi and hPi is not significantly different from the results in Fig. 5. Specifically, the results in Fig. 8 illustrate that the present definition of ''quantum equilibration'' in eqn (1) allows relaxation to states that are far from the traditional equilibrium scenarios, which are centered close to the potential minima. The full dynamics for the initial state C G (R 1 + 2.2a 0 ,R 2 À 2.2a 0 ) can be found in S4 (ESI †).
Other initial states, with resulting dynamics occurring outside the two main diagonals (concerted mechanism of proton transfer), can be also considered. Should the wavepacket be displaced along only one degree of freedom, one would likely observe the sequential mechanism of the proton transfer. In Fig. 9 and 10, we show the dynamics of hRi and hPi (and the corresponding expectation value dynamics of hR 1,2 i and hP 1,2 i) that result from the initial states C G (R 1 + 2.2a 0 ,R 2 ) and C G (R 1 + 2.2a 0 ,R 2 À 1.0a 0 ) respectively. The full dynamics for the initial states C G (R 1 + 2.2a 0 ,R 2 ) and C G (R 1 + 2.2a 0 ,R 2 À 1.0a 0 ) can be found respectively in S5 and S6 (ESI †).   As it is apparent from all the above dynamics, the symmetry of the initial state does not seem to play a crucial role. Instead, it seems to be the energy of the initial state in comparison with the height of the SP2 barrier what determines whether or not equilibration will take place. As the energy of the initial wavepacket is decreased and gets closer to the height of the SP2 barrier, it takes longer to equilibrate. Overall, this is in consonance with the theoretical predictions in Fig. 3.
In general, our results are in accordance with the idea that the diagonal and microcanonical ensembles give the same predictions for the relaxed value, 29 i.e.: where E 0 is the mean energy of the initial state, DE is the half-width of an appropriately chosen energy window centred at E 0 , and the normalization is the number of energy eigenstates with energies in the window [E 0 À DE,E 0 + DE]. Thermodynamical universality is evident in this equality: although the left-hand side depends on the details of the initial conditions through the set of coefficients c k , the right-hand side depends only on the total energy, which is the same for many different initial conditions. Explanations for the above universality have been proposed based on at least three possible scenarios. 29 Here we have proposed a possible explanation based on (i) and (ii) above. Specifically, in our example, E 0 and DE do depend on the height of the SP2 barrier.

Conclusions
To summarize, we have shown that quantum equilibration, mostly studied for many-body systems made of a large number of particles, can also play a role in isolated molecular systems that involve a few number of degrees of freedom. Specifically, we have seen that a two-dimensional model system describing the double proton transfer in porphine can present equilibration dynamics for a wide range of far-from-equilibrium initial states that are compatible with a pump-dump preparation scheme. We have seen that the resulting equilibration state is reached in B200 fs and that it is associated with zero mean position and momentum of the two protons. This equilibration state is characterized by a strong delocalization of the probability density of the protons, which can appear either along the concerted or sequential transfer paths depending on the symmetry of the initial conditions. By picking up a number of very general initial conditions, we have shown that the symmetry of the initial state does not play a crucial role. Instead, it is the energy of the initial state in comparison with the height of the SP2 barrier what determines whether or not equilibration will take place and how long it will take to develop. We have thus concluded that there are two main facts that must coexist to find equilibration: 40 (i) initial states consisting of a large superposition of energy eigenstates, and (ii) the presence of a very effective dephasing mechanism. For the porphine model considered in this work, condition (i) is met for initial states with energies well above the SP2 barrier, and it is the SP2 barrier itself that facilitates the fulfillment of (ii). This conclusion is compatible with the idea that the diagonal and microcanonical ensembles give the same predictions for the relaxed value for a wide range of initial conditions. 29 Relying on the de Broglie-Bohm interpretation of quantum mechanics, we have shown that the equilibration process in porphine can be associated with the picture where, due to strong interferences, the two protons become localized at certain configurations along the concerted and sequential transfer paths as if they where associated to a stationary state. This result may help, for example, in the interpretation of experimental data. That is, the close connection between Bohmian mechanics and weak values 64,71,72 provides an operational interpretation of the above described equilibration mechanism based on trajectories that could be put within reach of experimental verification. 62,63 Let us insist that Bohmian trajectories have been introduced in this work with the intention to gain insight into the physical mechanism that yields equilibration. However, the conclusion that equilibration is actually taking place is based solely on the analysis of the time-evolution of the position and momentum expectation values (which are not dependent on the particular interpretation choice).
Witnessing equilibration dynamics in the laboratory for an arbitrary type of molecules could be feasible thanks to the development of femtosecond spectroscopy 73,74 and to the most recent technological advances that allow characterizing wavepacket dynamics at an experimental level in a wide variety of molecular systems. [75][76][77] However, whether quantum equilibration dynamics in an isolated porphine molecule might be actually observed will strongly depend on the validity of the approximations considered in the employed model. The model porphine considered here does not include the effect of the heavy atom vibrations of the molecular scaffold, 78 and therefore, while it served well for a proof of principle, the excitation of specific heavy atom vibrations might have a sizable influence on the proton transfer dynamics. Yet, avoiding possible effects coming from the scaffold dynamics could be possible. Because of the bistability in the position of the two protons, which can occupy different, energetically equivalent positions (tautomerization), porphines have been devised as molecular switches. [79][80][81][82] In this respect, there exist an extensive literature (see e.g. ref. 80) on the control of the scaffold structural freedom by means of surface anchoring 81,83 or molecular assembling. 82,84 Finally, let us mention a recent work where the third law of thermodynamics is investigated for a 1D model boron rotor B 13 + . 55 In ref. 55 the entropy production reaches a stationary (or quasi-stationary) value (3.017). This value, however, does not correspond to the equilibrium value (3.401), which is far beyond the variance of the fluctuations (0.101). In the authors words: ''this points to incomplete equilibration''. Furthermore, this quasi-equilibration process occurs in approximately 10 picoseconds. At the time scale of a few picoseconds it is difficult to imagine that the coupling of the rotation/pseudo-rotation of the boron rotor to the rest of degrees of freedom will not modify the dynamics.

Conflicts of interest
There are no conflicts to declare.