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

PCCP View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. PAPER View Journal | View Issue Quantum equilibration of the double-proton transfer in a model system porphine† Cite this: Phys. Chem. Chem. Phys., 2020, 22, 22332 Guillermo Albareda, Josep Maria Bofill, Ivano Tavernelli f *ab Arnau Riera,c Miguel González, bd Iberio de P. R. Moreira, bd Rosendo Valero be bd and 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 Received 3rd June 2020, Accepted 7th September 2020 (B200 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 DOI: 10.1039/d0cp02991b 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 rsc.li/pccp dephasing mechanism. 1 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. a Max Planck Institute for the Structure and Dynamics of Matter and Center for Free-Electron Laser Science, Luruper Chaussee 149, 22761 Hamburg, Germany. E-mail: guillermo.albareda@mpsd.mpg.de, guillealpi@gmail.com b Institute of Theoretical and Computational Chemistry, Universitat de Barcelona, Martı́ i Franquès 1, 08028 Barcelona, Spain c ICFO-Institut de Ciències Fotòniques, The Barcelona Institute of Science and Technology, Castelldefels (Barcelona), 08860, Spain d Departament de Ciència de Materials i Quı́mica Fı́sica, Universitat de Barcelona, Martı́ i Franquès 1, 08028 Barcelona, Spain e Departament de Quı́mica Inorgànica i Orgànica, Universitat de Barcelona, Martı́ i Franquès 1, 08028 Barcelona, Spain f IBM Quantum, IBM Research – Zurich, 8803 Rüschlikon, Switzerland † Electronic supplementary information (ESI) available: (S1) Dynamics of the twoproton density for the initial state C0(R1,R2;2.2a0). (S2) Dynamics of the twoproton density for the initial state C0(R1,R2;1.9a0). (S3) Dynamics of the two-proton density for the initial state C0(R1,R2;1.75a0). (S4) Dynamics of the two-proton density for the initial state CG(R1 + 2.2a0,R2  2.2a0). (S5) Dynamics of the two-proton density for the initial state CG(R1 + 2.2a0,R2). (S6) Dynamics of the two-proton density for the initial state CG(R1 + 2.2a0,R2  1.0a0). See DOI: 10.1039/ d0cp02991b 22332 | Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 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 fact2–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 shown5,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–11 or suitable confinement potentials (that lead to effective low-dimensional continuous systems).12–15 Similarly, systems of trapped ions16,17 allow us to precisely study the physics of interacting systems in the laboratory.18–23 In such highly controlled settings, equilibration and thermalisation dynamics has been studied both experimentally24–28 and This journal is © the Owner Societies 2020 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Paper PCCP numerically, often with a focus on the so-called quenches, i.e. rapid changes of the Hamiltonian.29–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 average1 (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. 2 Definitions ^ ^ hAðtÞi ¼ hcðtÞjAjcðtÞi ¼ ðT 0 ^ hcðtÞjAjcðtÞidt: (1) 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 dT and whose Hamiltonian has a dE P ^¼ Ek jEk ihEk j, where Ek are its spectral representation H k¼1 X k;l   ^ l : ck cl eiðEk El Þt Ek jAjE (2) The long-time average in (1) is then: A ¼ X k   ^ k ; jck j2 Ek jAjE (3) 1 Ð T iðEk El Þt e dt. Therefore, T 0 if the system relaxes at all, it must be to the following density matrix: where we have used that dk;l ¼ limT !1 o¼ X k jck j2 jEk ihEk j: (4) As noted in ref. 29, the result in (4) can be thought as stating the prediction of a ‘‘diagonal ensemble’’, |ck|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: deff 1 : ¼ 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 1 A ¼ lim T!1 T of  can be written as: X k   jck j4 ¼ Tr o2 (5) 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 dT energy eigenstates to different energies is dT. Provided that a given observable equilibrates, it was recently shown40 that the equilibration time can be related with a dephasing time td as: Teq B td = p/sG. (6) where sG is the standard deviation of the energy gaps, Ga = Ej  Ei, when weighted by their respective relevances qa, i.e.: sG2 : ¼ X qa Ga2 ; (7) a where a A G = {(i, j):i, j A {1,. . .dE}, i a j}, and qa : ¼ jva j2 P  2 vb  b energies and |Eki are eigenstates of Ĥ.‡ Choosing units such that h = 1, the state at time t is then given by jcðtÞi ¼ P iE t ck e k jEk i with ck  hEk|c(0)i, and the expectation value with va = cj*Ajici/sA, and Aij := hEi|A|Eji are the matrix elements of A in the energy eigenbasis. The denominator sA = amax  amin is the range of possible outcomes, being amax(min) the largest (smallest) occupied eigenvalue of Â. ‡ Note that the sum runs over dE r dT 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 |Eki 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. k This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 | 22333 View Article Online PCCP Paper Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. 3 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 R1 and R2, respectively, from the domains of the reactant (R) to the product (P) (see Fig. 1). The PES model is41 i 2  2 U0 h V ðR1 ; R2 Þ ¼ 4 R21  D20 þ R22  D20 4GD20 R1 R2 D0 (8) þ 2Gð2 þ GÞU0 ; where G is the coupling parameter which formally represents the interaction between the two hydrogen bonds, 2D0 is the transfer distance of each proton, i.e., the distance between its two equilibrium positions, and U0 is the barrier height for single-proton transfer. The parameter U0 = 0.473 eV has been fitted in ref. 41 in order to account for the experimental results of nuclear magnetic resonance and laser-induced fluorescence measurements of ref. 44–46. The other two parameters, D0 = 1.251a0 and G = 0.063 are based on density functional theory calculations of Smedarchina et al.47 at the B3LYP/6-31G* level. The resulting 2D model PES is illustrated in the background of Fig. 1. The barriers are labeled TS (‘‘transition states’’) for two alternative reaction paths. 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 Fig. 1 Double proton tranfer of the model porphine. The protons move along coordinates R1 and R2. The four snapshots represent the transfer of the two protons from reactant (R) to product (P), sequentially along intermediate states (I) involving four transition states (TS), or simultaneously through a second order saddle point (SP2). In the background: Potential energy surface for the model porphine, eqn (8), adopted from ref. 41–43. The equidistant values of the contours range from 0 eV, for the potential minima of the reactant (R) and product (P) configurations, to 5 eV. The corresponding energies of the local minima for the intermediates (I), of the four barriers labeled TS, and of the second order saddle point (SP2) are 0.238, 0.600, and 1.069 eV, respectively. 22334 | Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 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 R1 = R2. It accommodates nearly degenerate doublets of eigenstates Cv+(R1,R2) and Cv(R1,R2), with energies below the barriers TS, plus higher excited states. We then chose our initial state to be of the general form: C0(R1,R2;DR) = C0,R(R1 + DR,R2 + DR), (9) where DR represents a shift along the diagonal R1 = R2 and 1 C0;R ðR1 ; R2 Þ ¼ pffiffiffiðC0þ þ C0 Þ is a superposition state that 2 represents the localized ground state wavefunction of the reactant, with C0+(R1,R2) and C0(R1,R2) representing the lowest doublet (n = 0), i.e., HCn(R1,R2) = EnCn(R1,R2), where H is the full Hamiltonian of the system, i.e.: H(R1,R2) = TR1 + TR2 + V(R1,R2), (10) 2 @ 2 h the kinetic energies of the two protons, 2M @Rk2 M = mp the proton mass, and V(R1,R2) the scalar potential given with TRk ¼  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 C4N 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 pump– dump laser pulse control,48 as designed by Tannor and Rice.49–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 Fig. 2 Due to the symmetry of the model potential used in this work, i.e., eqn (8), the states in eqn (9) are representative of the two equivalent hydrogen configurations in (a) and (b). This journal is © the Owner Societies 2020 View Article Online Paper PCCP Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. 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 4 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: R̂ = |R1ihR1| # I2 + I1 # |R2ihR2|. (11)    Ð1 In eqn (11), I1ð2Þ ¼ 1 dR1ð2Þ R1ð2Þ R1ð2Þ  are identity operators in the position representation acting, respectively, in the subspaces defined by R1 and R2. 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 deff in eqn (5) and the equilibration time Teq in (6) for different values of DR A [2.7,2.7] (see Fig. 3). For that, we first computed va = cj*Rjici/sR, where Rij := hEi|R̂|Eji and sR = Rmax  Rmin has been chosen to pffiffiffiffiffi be sR = 3sgs, and sgs ¼ 1 25 is the dispersion of the localized ground state wavefunction of the reactant. In addition to deff, which represents an equilibration criterion, and Teq, 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): T~eq ¼ Teq : deff  1 (12) Note that, by definition, even if Teq is small, T̃eq can still be close to infinite for initial states being eigenstates of the observable of interest (i.e., deff B 1). The resulting values of deff, Teq 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 C0 is just the sum of two energy eigenstates C0+ and C0. 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 Teq and deff. Note, for example, the relative minimum shown by Teq 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 Teq This journal is © the Owner Societies 2020 Fig. 3 Effective dimension deff (top panel), equilibration time Teq (mid B panel), and effective equilibration time T eq (bottom panel), all in dashed red line. The results are for the initial states defined in eqn (9) and different values of DR for the type of measurement defined in eqn (11) using a square grid of size 151a0 with a grid spacing of 0.08a0. 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. 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 deff, Teq 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. 5 Equilibration dynamics We now aim to study the dynamics that leads to equilibration for initial states induced by pump–dump laser pulse control of Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 | 22335 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. PCCP the form in eqn (9). Specifically, we choose DR = 2.2a0 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 C0(R1,R2;2.2a0) (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 R1 = R2 (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̂ = |P1ihP1| # I2 + I1 # |P2ihP2|. 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 Fig. 4 Dynamics of the two-proton density |C(R1,R2,t)|2 at four different time snapshots for the initial state C0(R1,R2;2.2a0). 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. 22336 | Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 Paper Fig. 5 Equilibration of the position (solid black line) and momentum (dashed blue line). In the inset: Ensemble average of the phase space as a function of time. Upper panel: The initial state, C0(R1,R2;2.2a0), has an energy of 4.885 eV, i.e., 3.7 eV above the SP2. Middle panel: The initial state, C0(R1,R2;1.9a0), has an energy of 1.905 eV, i.e., 0.8 eV above the SP2. Bottom panel: The initial state, C0(R1,R2;1.75a0), has an energy of 1.185 eV, i.e., 0.1 eV above the SP2. The values of deff, Teff, and B T eff for the three different initial conditions DR(a0) = {2.2,1.9,1.75}, are respectively deff = {29.89,33.39,48.2}, Teq = {35.12 fs,59.27 fs,72.75 fs}, and B T eq = {1.17 fs,1.77 fs,1.51 fs}. differences R%  hR̂(t)i and P%  hP̂(t)i as a function of time for the initial state C0(R1,R2;2.2a0). 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., C0(R1,R2;1.9a0) and C0(R1,R2;1.75a0), 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 This journal is © the Owner Societies 2020 View Article Online Paper PCCP Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. 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 C0(R1,R2;2.2a0), C0(R1,R2;1.9a0), and C0(R1,R2;1.75a0) can be found in S1, S2, and S3 (ESI†) respectively. 6 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–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–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 {Ra(t)} = a {R1(t),Ra2(t)} initially sampled from the probability distribution P(R1,R2) = |C0(R1,R2,2.2a0)|2 and whose time evolution is defined as ðt (13) Ra ðtÞ ¼ Ra ðt0 Þ þ vðRa ðt 0 Þ; t 0 Þdt 0 ; t0 This journal is © the Owner Societies 2020 Fig. 6 Time evolution of a randomly chosen set of five hundred Bohmian trajectories Ra1(t) sampled from |C0(R1,R2,2.2a0)|2. where vðRa ðtÞ; tÞ ¼   h 1 @C 1 @C  ; : Im M C @R1 C @R2 Ra ðtÞ (14) In Fig. 6 we show the time evolution of an ensemble of five hundred Bohmian trajectories Ra1(t) sampled from |C0(R1,R2,2.2a0)|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. Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 | 22337 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. PCCP 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 CG(R1 + DR1,R2 + DR2), where CG(R1,R2) represents a Gaussian approximation to the localized ground state wavefunction pffiffiffiffiffi of the reactant, with a dispersion sgs ¼ 1 25 and centered at the origin of coordinates. In particular, for DR1 = DR2, the wavepacket CG(R1 + DR1,R2 + DR2) is centered along the diagonal R1 = R2 that leads from one Intermediate to the other (i.e., perpendicular to the diagonal R1 = R2 that leads from the reactant to the product). Due to the symmetry of the model potential in eqn (8), the state CG(R1 + DR1,R2  DR1) around the Intermediate minima, may correspond to four equivalent hydrogen atoms distributions (see Fig. 7). Starting from the initial state CG(R1 + 2.2a0,R2  2.2a0), 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 R1 = R2 instead of R1 = R2 and there is a notable increase of the probability density along the sequential paths. The associated behaviour of hR̂i and hP̂i 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 Fig. 7 Due to the symmetry of the model potential defined in eqn (8), the initial states CG(R1 + DR1,R2  DR1) localized around the intermediate minima may correspond to the four equivalent proton distributions in (a–d). 22338 | Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 Paper Fig. 8 Equilibration of the position (solid black line) and momentum (dashed blue line) for the initial state CG(R1 + 2.2a0,R2  2.2a0). In the inset: Ensemble average of the phase space as a function of time. The initial state, CG(R1 + 2.2a0,R2  2.2a0), has an energy of 5.10 eV, i.e., 4.03 eV above the SP2. to the potential minima. The full dynamics for the initial state CG(R1 + 2.2a0,R2  2.2a0) 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 hR̂i and hP̂i (and the corresponding expectation value dynamics of hR̂1,2i and hP̂1,2i) that result from the initial states CG(R1 + 2.2a0,R2) and CG(R1 + 2.2a0,R2  1.0a0) respectively. The full dynamics for the initial states CG(R1 + 2.2a0,R2) and CG(R1 + 2.2a0,R2  1.0a0) can be found respectively in S5 and S6 (ESI†). Fig. 9 Top panel: Equilibration of the position (solid black line) and momentum (dashed blue line) for the initial state CG(R1 + 2.2a0,R2). In the inset: Ensemble average of the phase space as a function of time. The initial state, CG(R1 + 2.2a0,R2), has an energy of 3.08 eV, i.e., 2.01 eV above the SP2. Bottom panel: Dynamics of hR̂1,2i and hP̂1,2i. This journal is © the Owner Societies 2020 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Paper Fig. 10 Top panel: Equilibration of the position (solid black line) and momentum (dashed blue line) for the initial state CG(R1 + 2.2a0,R2  1.0a0). In the inset: Ensemble average of the phase space as a function of time. The initial state, CG(R1 + 2.2a0,R2  1.0a0), has an energy of 2.85 eV, i.e., 1.78 eV above the SP2. Bottom panel: Dynamics of hR̂1,2i and hP̂1,2i. 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.: X X    ^ k ¼ 1 ^ a ; A ¼ jck j2 Ek jAjE Ea jAjE (15) N E0 ;DE a2DE k where E0 is the mean energy of the initial state, DE is the half-width of an appropriately chosen energy window centred at E0, and the normalization is the number of energy eigenstates with energies in the window [E0  DE,E0 + 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 ck, 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, E0 and DE do depend on the height of the SP2 barrier. 7 Conclusions To summarize, we have shown that quantum equilibration, mostly studied for many-body systems made of a large number This journal is © the Owner Societies 2020 PCCP 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 values64,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 spectroscopy73,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–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 Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 | 22339 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. PCCP 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–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 anchoring81,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 B13+.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. Acknowledgements The authors acknowledge financial support from the Spanish Ministerio de Economı́a y Competitividad (Projects No. CTQ201676423-P and No. PID2019-109518GB-I00), the Generalitat de Catalunya (Project No. 2017 SGR 348), and the Spanish Structures of Excellence Marı́a de Maeztu program (through grant MDM-2017-0767). G. A. acknowledges also financial support from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 752822. Open Access funding provided by the Max Planck Society. Notes and references 1 C. Gogolin and J. Eisert, Rep. Prog. Phys., 2016, 79, 056001. 2 P. Reimann, Phys. Rev. Lett., 2007, 99, 160404. 3 S. Goldstein, J. L. Lebowitz, R. Tumulka and N. Zangh, Phys. Rev. Lett., 2006, 96, 050403. 4 S. Popescu, A. J. Short and A. Winter, Nat. Phys., 2006, 2, 754. 5 P. Reimann, Phys. Rev. Lett., 2008, 101, 190403. 6 N. Linden, S. Popescu, A. J. Short and A. Winter, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 061103. 7 M. Greiner, O. Mandel, T. W. Hänsch and I. Bloch, Nature, 2002, 419, 51. 22340 | Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 Paper 8 M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch and I. Bloch, Nature, 2002, 415, 39. 9 A. Tuchman, C. Orzel, A. Polkovnikov and M. Kasevich, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 74, 051601. 10 I. Bloch, Nat. Phys., 2005, 1, 23. 11 T. Langen, R. Geiger and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys., 2015, 6, 201–217. 12 L. Sadler, J. Higbie, S. Leslie, M. Vengalattore and D. StamperKurn, Nature, 2006, 443, 312. 13 T. Kinoshita, T. Wenger and D. S. Weiss, Nature, 2006, 440, 900. 14 S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm and J. Schmiedmayer, Nature, 2007, 449, 324. 15 A. Weller, J. Ronzheimer, C. Gross, J. Esteve, M. Oberthaler, D. Frantzeskakis, G. Theocharis and P. Kevrekidis, Phys. Rev. Lett., 2008, 101, 130401. 16 D. Porras and J. I. Cirac, Phys. Rev. Lett., 2004, 92, 207901. 17 A. Friedenauer, H. Schmitz, J. T. Glueckert, D. Porras and T. Schätz, Nat. Phys., 2008, 4, 757. 18 H. Häffner, W. Hänsel, C. Roos, J. Benhelm, M. Chwalla, T. Körber, U. Rapol, M. Riebe, P. Schmidt and C. Becher, Nature, 2005, 438, 643. 19 P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt and C. F. Roos, Nature, 2014, 511, 202. 20 B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach and G. Kirchmair, Science, 2011, 334, 57–61. 21 P. Schindler, M. Müller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller and R. Blatt, Nat. Phys., 2013, 9, 361. 22 R. Blatt and C. F. Roos, Nat. Phys., 2012, 8, 277. 23 J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk and J. J. Bollinger, Nature, 2012, 484, 489. 24 S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert and I. Bloch, Nat. Phys., 2012, 8, 325. 25 M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath and S. Kuhr, Nature, 2012, 481, 484. 26 T. Langen, R. Geiger, M. Kuhnert, B. Rauer and J. Schmiedmayer, Nat. Phys., 2013, 9, 640. 27 M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler and J. Schmiedmayer, Science, 2012, 337, 1318–1322. 28 J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch and U. Schneider, Phys. Rev. Lett., 2013, 110, 205301. 29 M. Rigol, V. Dunjko and M. Olshanii, Nature, 2008, 452, 854. 30 M. Rigol, Phys. Rev. Lett., 2009, 103, 100403. 31 M. Moeckel and S. Kehrein, Phys. Rev. Lett., 2008, 100, 175702. 32 C. Kollath, A. M. Läuchli and E. Altman, Phys. Rev. Lett., 2007, 98, 180601. 33 S. Deng, G. Ortiz and L. Viola, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 094304. 34 L. C. Venuti and P. Zanardi, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 81, 032113. This journal is © the Owner Societies 2020 View Article Online Open Access Article. Published on 14 September 2020. Downloaded on 12/15/2021 3:54:39 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Paper 35 F. Goth and F. F. Assaad, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 085129. 36 M. Rigol and M. Fitzpatrick, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 033640. 37 E. Torres-Herrera and L. F. Santos, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 043620. 38 A. J. Short, New J. Phys., 2011, 13, 053009. 39 A. J. Short and T. C. Farrelly, New J. Phys., 2012, 14, 013063. 40 T. R. de Oliveira, C. Charalambous, D. Jonathan, M. Lewenstein and A. Riera, New J. Phys., 2018, 20, 033032. 41 Z. Smedarchina, W. Siebrand and A. Fernández-Ramos, J. Chem. Phys., 2007, 127, 174513. 42 A. Accardi, I. Barth, O. Kühn and J. Manz, J. Phys. Chem. A, 2010, 114, 11252–11262. 43 G. Albareda, J. M. Bofill, I. Tavernelli, F. Huarte-Larrañaga, F. Illas and A. Rubio, J. Phys. Chem. Lett., 2015, 6, 1529–1535. 44 T. J. Butenhoff and C. B. Moore, J. Am. Chem. Soc., 1988, 110, 8336–8341. 45 J. Braun, M. Koecher, M. Schlabach, B. Wehrle, H.-H. Limbach and E. Vogel, J. Am. Chem. Soc., 1994, 116, 6593–6604. 46 J. Braun, H.-H. Limbach, P. G. Williams, H. Morimoto and D. E. Wemmer, J. Am. Chem. Soc., 1996, 118, 7231–7232. 47 Z. Smedarchina, M. Z. Zgierski, W. Siebrand and P. M. Kozlowski, J. Chem. Phys., 1998, 109, 1014–1024. 48 M. K. Abdel-Latif and O. Kühn, Theor. Chem. Acc., 2011, 128, 307–316. 49 D. J. Tannor and S. A. Rice, J. Chem. Phys., 1985, 83, 5013–5018. 50 D. J. Tannor, R. Kosloff and S. A. Rice, J. Chem. Phys., 1986, 85, 5805–5820. 51 M. Shapiro and P. Brumer, Quantum Control of Molecular Processes, John Wiley & Sons, 2012. 52 W. Li, X. Zhou, R. Lock, S. Patchkovskii, A. Stolow, H. C. Kapteyn and M. M. Murnane, Science, 2008, 322, 1207–1211. 53 A. S. Sanz and S. Miret-Artés, J. Chem. Phys., 2007, 126, 234106. 54 A. S. Malabarba, N. Linden and A. J. Short, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 062128. 55 D. Jia, J. Manz and Y. Yang, AIP Adv., 2018, 8, 045222. 56 A. Benseny, G. Albareda, Ã. S. Sanz, J. Mompart and X. Oriols, Eur. Phys. J. D, 2014, 68, 1–42. 57 X. Oriols and J. Mompart, Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology, Pan Stanford, Florida, USA, 2012. 58 G. Albareda, H. Appel, I. Franco, A. Abedi and A. Rubio, Phys. Rev. Lett., 2014, 113, 083003. 59 G. Albareda, A. Abedi, I. Tavernelli and A. Rubio, Phys. Rev. A, 2016, 94, 062511. 60 L. González and R. Lindh, Quantum Chemistry and Dynamics of Excited States: Methods and Applications, Wiley, 2020. This journal is © the Owner Societies 2020 PCCP 61 G. Albareda, A. Kelly and A. Rubio, Phys. Rev. Mater., 2019, 3, 023803. 62 S. Kocsis, B. Braverman, S. Ravets, M. J. Stevens, R. P. Mirin, L. K. Shalm and A. M. Steinberg, Science, 2011, 332, 1170–1173. 63 D. H. Mahler, L. Rozema, K. Fisher, L. Vermeyden, K. J. Resch, H. M. Wiseman and A. Steinberg, Sci. Adv., 2016, 2, e1501466. 64 F. Traversa, G. Albareda, M. Di Ventra and X. Oriols, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 87, 052124. 65 W. H. Miller, Faraday Discuss. Chem. Soc., 1977, 62, 40–46. 66 G. Albareda, X. Saura, X. Oriols and J. Suné, J. Appl. Phys., 2010, 108, 043706. 67 Y. Goldfarb, J. Schiff and D. J. Tannor, J. Chem. Phys., 2008, 128, 164114. 68 B. Doyon, 2019, arXiv preprint arXiv:1912.08496. 69 R. P. Feynman, Feynman’s Thesis – A New Approach to Quantum Theory, World Scientific, 2005, pp. 71–109. 70 D. Pandey, R. Sampaio, T. Ala-Nissila, G. Albareda and X. Oriols, 2018, arXiv preprint arXiv:1812.10257. 71 H. Wiseman, New J. Phys., 2007, 9, 165. 72 D. Pandey, R. Sampaio, T. Ala-Nissila, G. Albareda Piquer and X. Oriols, 2019, arXiv preprint arXiv:1812.10257. 73 I. Hertel and W. Radloff, Rep. Prog. Phys., 2006, 69, 1897. 74 D. J. Tannor, Introduction to Quantum Mechanics: a TimeDependent Perspective, University Science Books, 2007. 75 H.-G. Duan, V. I. Prokhorenko, R. J. Cogdell, K. Ashraf, A. L. Stevens, M. Thorwart and R. D. Miller, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 8493–8498. 76 L. Bruder, U. Bangert, M. Binz, D. Uhl, R. Vexiau, N. Bouloufa-Maafa, O. Dulieu and F. Stienkemeier, Nat. Commun., 2018, 9, 4823. 77 F. Ma, E. Romero, M. R. Jones, V. I. Novoderezhkin and R. van Grondelle, Nat. Commun., 2019, 10, 933. 78 Y. Litman, J. O. Richardson, T. Kumagai and M. Rossi, J. Am. Chem. Soc., 2019, 141, 2526–2534. 79 P. Liljeroth, J. Repp and G. Meyer, Science, 2007, 317, 1203–1206. 80 M. Jurow, A. E. Schuckman, J. D. Batteas and C. M. Drain, Coord. Chem. Rev., 2010, 254, 2297–2310. 81 W. Auwärter, K. Seufert, F. Bischoff, D. Ecija, S. Vijayaraghavan, S. Joshi, F. Klappenberger, N. Samudrala and J. V. Barth, Nat. Nanotechnol., 2012, 7, 41. 82 G. Bussetti, M. Campione, M. Riva, A. Picone, L. Raimondo, L. Ferraro, C. Hogan, M. Palummo, A. Brambilla and M. Finazzi, Adv. Funct. Mater., 2014, 24, 958–963. 83 Q. Zhang, X. Zheng, G. Kuang, W. Wang, L. Zhu, R. Pang, X. Shi, X. Shang, X. Huang and P. N. Liu, J. Phys. Chem. Lett., 2017, 8, 1241–1247. 84 E. Y. Li and N. Marzari, J. Phys. Chem. Lett., 2013, 4, 3039–3044. Phys. Chem. Chem. Phys., 2020, 22, 22332--22341 | 22341