Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation

The simulation of open quantum dynamics is a critical tool for understanding how the non-classical properties of matter might be functionalised in future devices. However, unlocking the enormous potential of molecular quantum processes is highly challenging due to the very strong and non-Markovian coupling of ‘environmental’ molecular vibrations to the electronic ‘system’ degrees of freedom. Here, we present an advanced but general computational strategy that allows tensor network methods to effectively compute the non-perturbative, real-time dynamics of exponentially large vibronic wave functions of real molecules. We demonstrate how ab initio modelling, machine learning and entanglement analysis can enable simulations which provide real-time insight and direct visualisation of dissipative photophysics, and illustrate this with an example based on the ultrafast process known as singlet fission.


Supplementary
(a) DP-Mes structure with symmetry elements of the D 2d point group, and (b) frontier molecular orbitals used in calculations and this supplemental note. Note that extremal hydrogens of mesitylene side-groups in (a) weakly break the point group symmetry. However, this is irrelevant for the low-energy excitations since the frontier orbitals shown in (b) firmly reside on the pentacene units.

Supplementary Note 1 -Electronic structure calculations
Off-diagonal couplings and symmetries. To first order in nuclear displacements, the matrix elements of the diabatic potential matrix between electronic states i and j are given by 1 where the Q n represent normal-mode coordinates of the ground state. The diagonal elements are simply the adiabatic forces (expressed in terms of normal coordinates), which we obtain directly from DFT. For the purposes of the symmetry analysis and evaluation of the off-diagonal coupling vectors we employ an approximate multi-Supplementary Reference description of the electronic states, following Supplementary Reference 2. This uses a small active space comprised of the four orbitals h A , h B , A , B , namely the monomer-localised HO-MOs/LUMOs corresponding to the two pentacene units that make up DP-Mes (see Supplementary Figure 1). As can be seen from the figure this set of localised orbitals is already orthogonal due to symmetry, meaning that a separate Löwdin orthogonalisation step is not necessary. In this set the ground state can be represented as the single determinant |h A αh A βh B αh B β|, where α and β label the electron spins. Furthermore, we have the locally excited states of each pentacene unit: Supplementary as well as the corresponding CT states where the excited electron is transferred to the other pentacene unit, respectively: To obtain states that transform as irreducible representations of the point group we need to perform an (anti)symmetrisation step: Finally, the correlated triplet pair TT can be expressed as: Using this picture one can now classify the electronic states according to irreducible representations of the point group. To evaluate the diabatic potential matrix W to first order we can now go ahead and substitute the frontier orbital representations of the electronic states. Application of the Slater-Condon rules and some algebra then yields the following off-diagonal coupling contributions: The matrix elements W (ab) between single orbitals a, b are now simply the Coulombic forces on the nuclei due to the overlap density a(r)b(r). It should be pointed out that the first-order coupling vectors cannot have contributions from two-electron matrix elements since the electron-electron Coulomb interaction is independent of the nuclear coordinates and vanishes when taking the derivative in Supplementary equation (1). Any vibrational mode can only give a non-zero contribution to the coupling between two electronic states if the product of the respective irreducible representations yields the totally symmetric representation A 1 . For example modes contributing to the coupling between LE + (B 2 ) and CT + (A 2 ) necessarily transform as B 1 (see Supplementary Figure 2 and Supplementary Table II).
Limitations. There are a number of limitations of the methods and approximations described. These include treating the normal modes in the harmonic regime, and only accounting for the coupling between vibrations and electronic states to first order in the nuclear displacements. Furthermore, the evaluation of the off-diagonal elements of the potential matrix is based on a description of the electronic states in terms of a small active space of frontier-orbitals only. This is necessary to correctly capture the spin-structure of the states involved in fission, at the expense of the representational power of the basis used. In addition, these wave functions are constructed from Kohn-Sham orbitals which are eigenfunctions of the non-interacting Kohn-Sham system rather than the full electronic Hamiltonian. Since these orbitals are expressed in terms of Gaussian basis functions centered on the nuclei we also do not capture Pulay terms when evaluating the off-diagonal elements of W according to Supplementary equation (1). Finally, it should be noted that CT energies from TDDFT, while not severely underestimated thanks to the use of range-separation, instead come out quite significantly higher than LE. It can be suspected that the lack of diffuse functions in the basis set and the linear-response nature of the TDDFT calculations play a role here. Both hamper the ability to capture the intra-molecular relaxation of a CT configuration.
Given all these limitations, our parametrisation of the vibronic coupling Hamiltonian should be considered semiquantitative.

Supplementary Note 2 -Transformation of vibrational environments
For the efficient application of our TTNS method we want the Hamiltonian to assume a star topology, with the exciton system at its centre and each arm representing one independent vibrational environment mapped onto a 1D-chain. In order to achieve this we need to further analyse all coupling matrices W n . A collection of modes ω n belong to the same independent environment i and can be mapped onto a 1D-chain if all modes couple with the same system operator W i . This means their coupling matrix elements W n,(jl) can be factorised as where W (jl) is the vectorisation of the system operator and λ n is the relative coupling strength for each mode n. This definition is equivalent to requiring rank(W n,(jl) ) = 1, or that it should have only one singular value/principal component. In particular, to ensure that resulting W are symmetric and stay zero where required by symmetry, in the following we restrict (jl) strictly to the upper right triangular part including the diagonal and excluding entries which are zero across all W n,(jl) . This yields the normalisation W A1 = 1 for the diagonally coupled A 1 and W = √ 2 for the off-diagonally coupled environments.
Clustering. For a chain mapping of arbitrary W n we need to find the correct partition of the modes into subsets, such that each subset fulfills above condition independently. We can find this through clustering methods, e.g. kmeans clustering 3 , where the cluster centroids represent the approximate (mean) optimal system coupling operator W i . We modified the algorithm to adjust the centroids weighted by λ n to reduce approximation errors for strongly coupled modes. In particular we need to find an optimal number of i clusters, such that the approximation error is minimised. An objective t-SNE projection 4 of the couplings (Supplementary Figure 3) shows that four clusters are Supplementary the minimum, which also agree with the irreps A 1 , A 2 , B 1 , and B 2 as obtained from group theory. However, since the clusters are fairly broad, our cluster analysis yields that at least seven clusters (coloured as in Figure 3) are necessary for an accurate representation, thus subdividing irreps A 1 , B 1 , and B 2 .
After the clustering we can map the initial linear vibronic coupling Hamiltonian onto the star Hamiltonian by numerically applying orthogonal polynomials or the Lanczos method 5,6 to each subset i independently. This transformation finally leads to the star Hamiltonian with the coupling vector λ i = (λ i1 , . . . , λ in ), containing all couplings of assigned modes, and with the optimal coupling matrices (columns ordered as TT, LE + , LE − , CT + , CT − ): Figure 3. 2D-Projection of the normalised vectorised coupling matrices Wn using t-SNE already reveals environments of different symmetry. The crosses represent the cluster centroids W i and colours the assigned modes. Especially the ring structure makes it necessary to split B2 into at least two baths.
After removing the strongly anharmonic low-frequency modes and irrelevant high frequency C-H stretches (time scales much faster than observed dynamics), 252 out of 318 modes are left. In Supplementary Table III we give all modes out of the 252 which contain at least 90% of the total coupling in each environment. In Supplementary Table IV we give the frequencies and coupling of the first sites, the reaction coordinates (RC) of each environment after chain transformation. This star Hamiltonian model minimises the number of nearest neighbours to two for all vibrational modes and 7 for the electronic system. This is an essential optimisation needed to allow efficient TTNS simulations. Furthermore, Supplementary the RCs contain valuable information about the collective action and effect of their respective environment. The RC coupling and frequency indicate total amount of coupling and weighted average frequency of all represented modes. This allows estimation of time scales and coupling regime. Additionally, expectation values of the RC's occupation and displacement represent weighted mean values for all constituent modes and can be used as a good proxy for the environment's collective state. This greatly reduces the amount of observables that have to be considered and supports the correlation with observables of the electronic system.

Supplementary Note 3 -Tree Tensor Network States
The presented simulations were performed with the time-dependent variational principle (TDVP) for matrix product states (MPS) with an optimised boson basis (OBB) as described in detail in Refs. 6-8. To accommodate the extended structure of the DP-Mes Hamiltonian H star , with seven chains coupling to the electronic system, we generalised our previously described method to Tree Tensor Network States (TTNS), the higher dimensional generalisation of an MPS, which is structurally similar to the hierarchical tucker format 9 , commonly used in the multi-layer multiconfigurational time-dependent Hartree-Fock (ML-MCTDH) algorithm 10 . We will briefly outline the features of the algorithm, for more information we refer to Refs. 6-8, 11-13. An MPS can approximate a many-body wave function where n k indexes the local Hilbert space of site k with local dimension d k , by decomposing the single rank-L tensor Ψ n1,...,n L , containing all probability amplitudes, into L rank-3 tensors by means of a low-rank tensor approximation, for the case of a 1D chain Hamiltonian of length L. The bond dimensions D k directly relate to the maximum amount of entanglement between the site k and its neighbours which can be encoded by the MPS. In order to optimally capture the entanglement arising in the time-evolution of our model, with the smallest possible D k , we advanced our MPS tensor network to a TTNS (see Supplementary Figure 4). Entanglement Renormalisation Tensors. In particular, we shape the TTNS such that each tensor represents one physical (local) site of H star . The tensor connections follow the coupling structure in H star to minimise inter-site entanglement. Since the central tensor (blue shaded in Supplementary Figure 4) connecting the system (red) with all seven environments (green) grows rapidly with increasing bond dimension (D 7 ), we decompose it into a tree-network of six entanglement renormalisation tensors (ER-nodes) (blue). The optimal network structure is determined by analysing the entanglement within the original core tensor in a preceding, expensive calculation. The key insight is that although this calculation of the star MPS network with D ∈ [5,8] is highly inaccurate, the core tensor still exhibits the major entanglement structure between the environments. By calculating the von Neumann entanglement entropy of all possible decompositions of this tensor, it is possible to find an optimal tree-network which minimises the occurring entropies (see Supplementary Figure 4).
The ER-nodes help to capture inter-environment entanglement and pass on only the most relevant effective multienvironment states. This way we can reduce the number of free parameters from D 7 down to 6D 3 and turn the initial exponential scaling into an almost linear scaling in the number of environments. This technique is key to enable efficient simulation for > 3 environment models with high accuracy. For comparison, while in the star MPS without ER-nodes we can reasonably simulate with up to D ∈ [5, 10] between the system and each environment, in the TTNS we can use D N ode ∈ [50, 80] while still gaining a speed-up and higher accuracy (see Supplementary Figure 5).
Time-evolution with TDVP. For the time-evolution of the TTNS we decided to use the TDVP, despite previous work reporting the time-evolving block decimation (TEBD) method for TTNS 13 . TEBD involves a substantial amount of excess operations for index swaps in order to apply nearest neighbour Trotter gates which would be less feasible for the ER-tree network. Additionally TDVP offers a better complexity scaling than TEBD due to its accurate 1-site integrator, allowing larger bond dimensions. Importantly, the TTNS does not contain loops, which are typical problems for PEPS and MERA 14 , it can be properly normalised and TDVP can be applied. Since we perform TDVP in a similar way as already described in 6,8,15 , we will only give a brief outline here.
Although the TTNS seems to consist of functionally different components (ER-tensors, MPS chains, OBB matrices), we can identify them as nodes having one parent and multiple child tensors. Since the quasi 1D structure of the tree  Figure 4. The TTNS tree used to simulate the model for DP-Mes. The root node (red) represents the five-level exciton system comprising TT, LE ± , and CT ± . The seven vibrational chain environments (green) are connected to the system through a tree network of entanglement renormalising tensors (blue). The long-time entanglement entropies of each bond are indicated next to the edges, as well as two example entanglement compression rates. This tree structure optimises the inter-chain entanglement such that only a small number of relevant environmental states are coupled to the exciton system, leading to a reduced bond dimension without reducing accuracy.
ensures correct TTNS normalisation we can use the insights from Supplementary Reference 6 and generalise it to nodes with more than one child. We recursively define the time-evolution of a single node as: 1. Evolve node t → t + ∆t 2 2. for i = {1 to N child } (a) Create C to child node i and evolve t + ∆t 2 → t (b) Contract C with child i and evolve t → t + ∆t (c) Create C from child i and evolve t + ∆t → t + ∆t 2 (d) Contract C from child with node 3. Evolve node t + ∆t 2 → t + ∆t Where we denote C the center matrix obtained by SVD. Starting this recursive scheme from the root node with a pre-defined ordering of child nodes constitutes a sweep. To obtain a symmetric second-order integrator with error of order O(∆t 3 ), we sweep through the tree such that each tensor is updated twice by a time step of ∆t/2. The effective Hamiltonians needed for the time-evolution of a node are obtained in a very similar way as for the propagation of a chain tensor, by contracting the Hamiltonian with all MPS tensors which are not being propagated in the current update. This algorithm generalises and applies very naturally to ER-tensors since they can also be seen as time-dependent tensors for which we can calculate an effective Hamiltonian. To apply the matrix exponential of the effective Hamiltonians to a tensor during the update, we use a customised Krylov subspace integrator based on EXPOKIT 16 .
Finally, this scheme can also be applied recursively to any tree tensor network with more layers. This allows in principle the simulation of arbitrary Hamiltonians which are representable as a tree without loops.

Supplementary Note 4 -Convergence of Simulation
We assured sufficient convergence of presented simulations by sweeping the parameters limiting the maximum bond dimensions D Node,max and D Chain,max and comparing the results. As measure of accuracy we compare simulations with pairs (D Node,max , D Chain,max ) to an expensive benchmark simulation by taking the relative mean absolute error ε of the TT population kinetics. To find the performance optimum as balance between accuracy and computation time we present the error ε against the CPU time needed to propagate the TTNS to the next timestep (see Supplementary Figure 5). We chose as benchmark (120, 30) since results converged sufficiently and t CP U /Sweep = 10.2 min was the upper limit to finish a simulation within 2 weeks. We found that the accuracy is mainly limited by D Node,max since any D Chain,max ∈ [20, 30] yielded very similar errors. With ε < 3%, the optimal truncation bounds were set to (80, 30) which gave a total CPU time of 4 days on a single core of an Intel Xeon E5-2670 for the simulation of 1 ps.

Supplementary Note 5 -Dynamic energy surfaces
One of the central operators relevant for the time-evolution of the TTNS and for finding the ground state in DMRG is the effective Hamiltonian of a site. This operator, constructed by contractingĤ with all tensors except for the tensor of the chosen site, contains all information relevant for the time-evolution of the site tensor. It is this remarkable property which allows us to gain deeper insights of the dynamics and its pathways. Here we will outline how we can derive the potential energy surface (PES) and the total energy surface (TES) from the effective Hamiltonian for the exciton system, as well as the corresponding adiabatic and non-adiabatic basis states.
To introduce the concept we derive the energy surfaces of a simplified example (quantum double well). In a single mode spin-boson model (SBM)Ĥ with spin Hamiltonian H S = ∆ 2 σ x + 2 σ z , tunnelling rate ∆, coupling λ, and mode frequency ω 0 , which is similar to the quantum Rabi model in case of vanishing energy level bias , we can obtain the operator for the potential energŷ V by subtracting the kinetic energyT = ω0 2 (P 2 − 1 2 ) where we used the mass-frequency normalised coordinatesQ =â † +â √ 2 =x √ mω,P = iâ † −â √ 2 =p √ mω and divided the zero point energy equally onT andV . In the next step we replace bath operators by their scalar valuesQ → Q,P → P At this point, it would not be useful to replace the bath operators by scalar phase space coordinates (Q i,k ,P i,k ) → (Q i,k , P i,k ) since this would yield a 500-dimensional coordinate space which is impossible to analyse. Even the restriction to the RC only would still result in a 7D PES or 14D TES. Instead, we propose to only evaluate the operators along the trajectory taken by the time-evolved wave function which unravels the multi-dimensional energy surfaces into a single dimension, with time as the generalised reaction coordinate. Here, A l 1k is the TTNS tensor for the electronic system and we know from the Schmidt decomposition that there are as many bath states |φ k (t) as there are system states |l for maximum system-bath entanglement.
Each |φ k (t) can be interpreted as a particular bath configuration (displacement) or superposition of multiple bath configurations. We evaluateĤ star andV star by contracting them with the time-dependent bath-states where is the effective Hamiltonian of the system tensor A l 1k . To obtain surfaces similar to the ones described in the single-mode SBM, we would need to diagonalise H ef f ii (t) and V ef f ii (t), yielding the TES and PES created by |φ i (t) only. However, to understand the observed dynamics with strong system-environment entanglement, where the Born-Oppenheimer approximation breaks down, this picture is insufficient since there are also non-vanishing terms for j = lcoupling different bath states. In order to account for them, we diagonalise the full yielding the non-adiabatic basis U H (t) and corresponding TES in E H (t), and the adiabatic basis With access to the basis maps U H (t) and U V (t) we can also obtain the TES in the adiabatic basis and PES in the non-adiabatic basis as the diagonals ofẼ which may also contain couplings since they are not purely diagonal any more. In addition to the PES and TES landscape we want to obtain information about the population of each surface in order to track the pathway of SF. Each system wave function associated with a bath state j can be interpreted as an individual wave packet on the respective potential. By rotating the vectorised system tensor wave function A ∈ C D·ds into the basis defining the PES/TESÃ we can get the probability amplitudeÃ k on each surface k. The populationP k is then given by the element-wise productP Furthermore we can calculate the contribution of each diabatic system state i on each surface k by using Since this expression does not include coherences between surface (adiabatic) states, it yields incorrect diabatic populations when summing over k and can only be taken as an estimate. Vice-versa, we can calculate the contribution of each surface k to each diabatic population P i via which will yield incorrect surface populations when summing over i since the coherences between diabatic states are not included.
To give an example, we present in Supplementary Figure 6 the TES (a,c) and PES (b,d) plotted with respect to two different basis states. In particular, we distinguish between the non-adiabatic energy eigenstates U H of H ef f (a,b), the adiabatic states U V (c,d), and the diabatic basis of the electronic system TT, LE ± , CT ± (colour of fills). The filled areas indicate by how much the surface is populated, while the colours give an estimate of diabatic mixing.

Supplementary Note 6 -Supporting energy surface results
The SF dynamics mainly occupy two relevant energy surfaces, one with major LE + and the other with TT character, which we will focus our discussions on.  Figure 7c) we can identify the actual crossing of the two coupled surfaces. In the following we detail main observations leading to SF and compare them to environmental dynamics and RC properties to investigate the mechanism. At 20 fs the unoccupied TT TES is pushed above the LE + (Supplementary Figure 7a) through interaction with the B 22 and A 12 environments, as confirmed by comparison with their TT-projected RC displacements ( Supplementary  Figure 7b). While a larger x RC in B 22 raises the TT energy, the effect of A 12 is opposite as evident from the bottom three TT surfaces which move anti-parallel to the displacement. We can find similar timescales in the eigenfrequencies of the respective RCs of ω RC,A12 ≈ 1377 cm −1 ∼ 24 fs and ω RC,B2 ≈ 1100 cm −1 ∼ 30 fs, confirming our observations. The AC at 20 fs with a gap of 6 meV is traversed rapidly by the TT, leading to almost no population exchange, while the AC at 80 fs with a gap of 12 meV seems to be the main driver of SF. The relatively weak coupling of 3 meV (half AC gap) is a result of slow B 11 reorganisation which happens on a longer ω B1 ≈ 400 cm −1 ∼ 83 fs timescale leading to a larger effective vibronic coupling of 6 meV at 80 fs. Also the lowering of the TT TES is correlated with the excitation of the B 11 RC, which leads to traversing the crossing at the same time as the effective vibronic coupling is maximised. The crossing velocity of the two surfaces is low enough to keep the majority of initial LE + population on   very similar dynamics (Supplementary Figure 8). This observation suggests that there is only little qualitative difference between the 252 modes and 75 modes model which is also confirmed by the parameters of the RCs in Supplementary Table V. Although this might be a reasonable argument against simulating full model dynamics, it has to be pointed out that there is only a negligible simulation speed-up with the 75 modes model. The reason is that these dynamics are causing the same amount of inter-environment entanglement that has to be compressed in the ER-network, which is the major scaling bottleneck, such that decrease in the size of each environment can be neglected.
We studied the effect of an external dissipative environment on the SF dynamics by broadening the vibrational spectra. We broaden each vibrational mode by 10 cm −1 with a Lorentzian lineshape while conserving the reorganisation energy. The resulting spectral densities are mapped onto chain environments using orthogonal polynomials as described in Supplementary Reference 6. We truncate each semi-infinite chain to a length of L = 70, giving a total of 490 modes, such that the exciton dynamics are unaffected by reflections from the end of the chains. The resulting dynamics show a lower TT yield of 80% along with a three times larger LE + population of 15%. The absence of further kinetic energy seems to limit SF to this upper value. However, as soon as the reflection from the end of the chains impact on the system past 1 ps, the TT increases again by a few percent. This suggests, that even after an equilibrium is reached, excitation of the environment, such as heat or a push pulse, could increase the yield of fission.
The super-exchange character of LE − B1 − − → CT − B1 − − → TT is confirmed by the appearance of Im(ρ TT,LE − ) and the large static coherence Re(ρ TT,LE − ) (Supplementary Figure 9). The latter manifests as a coherent superposition of TT, LE − , and CT − in the long-time state (> 300 fs). The mixture of this strong vibronic (dressed) state is modulated by B 1 and A 1 (Supplementary Figure 8c). All other coherences are small compared to ρ TT,LE − , as expected for linear vibronic coupling.
The von Neumann entropy S of the RDM increases until 90 fs which shows the strong non-adiabatic character of the SF dynamics (see Supplementary Figure 9). At this point the electronic state is maximally distributed among LE ± , CT ± , and TT. As the TT population increases subsequently, S decreases again due to purification of the state and demixing with the environment.

Supplementary Note 8 -Rabi Dynamics
The < 20 fs dynamics are dominated by very fast vacuum Rabi-like oscillations. However, correlating system and environment observables, such as n RC (see Figure 3) to visualise these is indicative, since expectation values are averages over microscopic states of competing dynamics, such as environmental reorganisation. By taking n only with respect to one of the five effective (microscopic) bath states φ j |n|φ j we can single out one particular type of process. The system tensor wave function A ∈ C D×ds tells which pairs of system and bath states are dominant.
The Rabi dynamics LE + B1 ← → CT + are happening between |LE + , φ 1 and |CT + , φ 3 , as evident from the antiparallel dynamics of ρ(CT + ) and φ 3 |n B1,RC |φ 3 (Supplementary Figure 10b). We can see that B 1,2 only participates to 2% while B 1,1 has a contribution of 98%. Since the driving is very off-resonant and comparably weak, with a RC frequency ω = 450 cm −1 , an LE + -CT + energy gap ∆ ≈ 5500 cm −1 − 6500 cm −1 , and coupling λ ≈ 1350 cm −1 , the Rabi frequency Ω R is very large (T ≈ 5.1 fs → Ω R ≈ 6500 cm −1 ). It should be noted that in this regime the counter-rotating terms cannot be neglected such that the usual rotating-wave approximation (RWA) is invalid. Thus we obtain the Rabi frequency as 18  Figure 11. A2 RC occupation projected onto TT. Most of the initial < 10 fs TT was generated via pathway I as indicated by the large nRC,A 2 > 0.7. At the equilibrium state nRC,A 2 ≈ 0.07 shows that only 7% of the total TT came from path I, while the rest has taken path II.

Supplementary Note 9 -Exclusion of A 2 pathway
We can confidently prove that pathway I coupled via the A 2 environment is not contributing substantially to the observed SF dynamics by looking at the TT projected occupation n RC,A2,TT = Ψ| P TTnRC,A2 |Ψ , where P TT = |TT TT| is the projector onto TT. Since the vibrationally mediated transport CT + → TT can only happen through the excitation of A 2 as |CT + |0 A2 → |TT |1 A2 , while pathway II via LE − → CT − → TT happens at |0 A2 , the magnitude of this TT projected occupation gives direct evidence of how much TT has come via pathway I. Initially, n RC,A2,TT ≈ 1 indicates that almost all TT was generated through path I. Subsequently this occupation drops rapidly, indicating a larger population of |TT |0 A2 which can only happen via path II. The final value n RC,A2 ≈ 0.07 shows that 93% of TT has been generated through path II.
Supplementary Note 10 -LE + → TT as Landau-Zener transition We can interpret the passage through the avoided crossing in terms of a Landau-Zener passage which allows the adiabatic transition LE + → TT if the decrease in energy gap is slow enough. Based on the Landau-Zener formula for the probability of diabatic transition with a = 6 meV being the effective coupling between the surfaces and α = 0.727 meV/fs the rate of change of the energy gap, we can calculate a 62% transition probability which is not properly matching the observed kinetics. The formula would predict a coupling of 10 meV to yield the expected 25% transition probability. Given that the Landau-Zener formula is based on multiple assumptions which are all violated in these real dynamics, such as a time-independent diabatic coupling, linearly varying energy gap, and the perturbation parameter of the Hamiltonian being a linear function of time, it is ensuring that the found and expected effective coupling are at a comparable order of magnitude.