Is there a finite mobility for the one vibrational mode Holstein model? Implications from real time simulations

Charge carrier mobility is an important property of organic semiconducting materials and is critical in many of their applications.^1–3 It is now well known that charge carrier transport in these materials is significantly affected by the electron–phonon interaction.^4–7 Depending on the strengths of the intermolecular electronic coupling and the electron–phonon coupling, the charge carrier transport mechanism can range from band-like behavior^2,8 to hopping transitions⁹ and cases in between.^10–12 As a unified theory is not yet available, charge carrier transport mechanisms in organic materials continue to be an active area of theoretical research.

The Holstein model^4,13 has been widely applied to study charge carrier transport in organic molecular crystals and has been indispensable in understanding such processes.^14,15 Analytically, solutions of the Holstein model are only available in the strong and weak coupling limits.^16–18 In the case of arbitrary coupling, approximate^19–21 and numerically exact^22–26 methods can be employed. For example, quantum Monte Carlo (QMC)^27–29 simulation can be used to obtain imaginary time properties, while charge carrier transport mobility can be obtained from an inverse Laplace transformation of the imaginary time results. The Green function methods^{16,23,30–34} can also be used to numerically simulate the Holstein model and obtain charge carrier mobility.

Charge carrier mobility can also be obtained by calculating the real time current autocorrelation functions using methods such as the density matrix renormalization group (DMRG)^25,35,36 approach. Since long time simulations are still very challenging, the correlation functions are usually truncated at a finite time. For this treatment to be valid, it is required that the long time behavior of the autocorrelation function either is not important or can be extrapolated using simple approximations.

A particularly interesting problem in the literature is the one vibrational mode Holstein model, where each site couples to a local vibrational mode. It is still not very clear whether there is a well-defined mobility for this problem when there is no static disorder or dissipation. In the literature, charge carrier mobilities have been obtained using both imaginary time and real time methods.^{10,15,37–39} However, in a recent study of the friction-free Holstein model by Kloss et al.,⁴⁰ it was found that the mean square displacement (MSD) does not show a linear dependence of time t for a wide range of parameters. This result indicates that a diffusion constant may not exist and, thus, the charge carrier transport mobility cannot be defined for the standard one vibrational mode Holstein model.

In previous studies, our group has employed the hierarchical equation of motion (HEOM) method to simulate the Holstein and Holstein–Peierls models.^39,41,42 The HEOM approach is a “numerically exact” method to simulate quantum dynamics in condensed phases^43,44 and is widely applied in different fields, including excitation energy transfer (EET) dynamics and related spectroscopic phenomena,^45,46 as well as non-equilibrium charge carrier transport dynamics.⁴⁷ In all our previous studies of the Holstein model,^39,41,42 dissipative modes were employed for the vibrational degrees of freedom (DOFs), and well-defined diffusion constants and mobilities can be obtained.

In this work, we investigate whether a well-defined mobility exists for the 1D one vibrational mode Holstein model by starting from a 1D Holstein model with dissipation and systematically decreasing the dissipation strength toward the friction-free limit. More specifically, we apply the methods previously developed in Refs. 48 and 49 to a Holstein model with a Brownian oscillator spectral density, where the strength of dissipation is controlled by a single parameter. Charge carrier mobilities for the dissipative model are then calculated from both the MSD and the real time current autocorrelation functions and extrapolated to the friction-free limit.

In Eq. (12), the reorganization energy is determined by the vibrational frequency Ω and the coupling constant C in Eqs. (3) and (4) as λ = C²/Ω². So, the parameter γ in Eq. (12) actually controls the strength of the dissipative effects. To study the effects of dissipation on the charge carrier mobility, we chose to fix the vibrational frequency Ω and reorganization energy λ and investigate how γ affects the charge carrier transport. Figure S1 in the supplementary material shows the BO spectral density J(ω) for different values of γ, with λ = 1.0 and Ω = 1.0. It can be seen that the shape of J(ω) changes significantly as γ varies. When γ = 0, J(ω) is a delta function, and η = 0 in Eq. (7), the coupling strength between the vibrational mode and the ohmic bath now vanishes, and the standard friction-free 1D Holstein model is recovered. As the coupling between the vibrational mode and the ohmic bath increases, γ also increases. Consequently, J(ω) is broadened, and the position of the peak shifts to lower frequencies. When 0 < γ < Ω, the system is underdamped, and when γ > Ω, the system is overdamped.

To obtain the current autocorrelation function $⟨ĵ(t)ĵ(0)⟩$⁠, we first use the imaginary time HEOM to calculate the correlated initial state, i.e., the equilibrium $ρ̃n$⁠. The current operator $ĵ$ is then multiplied by the equilibrium $ρ̃n$ and the new quantity $ĵρ̃n$ is used as the initial state to propagate the real time HEOM. Finally, we multiply the current operator by the reduced density operator at time t and trace over the electron DOFs to obtain the correlation function $⟨ĵ(t)ĵ(0)⟩$⁠. The exact form and details of the real and imaginary time HEOM can be found in the supplementary material.

By applying the real time HEOM, we first calculate the MSD for the dissipative Holstein model. The periodic boundary condition is employed, and two different sets of parameters are used in the simulation. In the first set of parameters, λ = 1.0, Ω = 1.0, β = 1.0, and J = 0.1. Since the reorganization energy λ is significantly larger than the electronic coupling constant J, it is supposed to be in the strong electron–phonon coupling regime. The relatively large λ parameter is used as a relatively small number of sites are needed to obtain convergent results. Five different values for the friction parameter γ = 0.05, 0.2, 0.5, 0.8, 1.5 are used in the simulation, where the first four are in the underdamped regime (γ < 1.0) and the last one is in the overdamped regime (γ > 1.0).

The initial state of the real time HEOM simulation is $|0〉〈0|⊗e−βĤvib$⁠. The total number of sites used in the simulation is up to 21, depending on different values of γ, ensuring that the MSD has reached the linear growth region and the calculated diffusion constant is not affected by boundary effects. The HEOM simulation is based on the matrix product state (MPS) method with a bond dimension of 80. The simulated MSD is shown in Fig. 1(a). It can be seen that with an increase in the friction constant γ, the slope of the MSD curve becomes smaller. In addition, the time to reach the linear growth region of the MSD also increases with the decrease in γ.

In another set of parameters, we use λ = 1.0, Ω = 1.0, β = 1.0, and J = 0.5. Since the reorganization energy λ and the electronic coupling constant J are now comparable in magnitude, this is a case in the so-called intermediate coupling regime. As the coupling constant J is now much larger than that in the first set of parameters, the real time HEOM propagation becomes more challenging. To obtain converged results, the number of sites needed is 51 for γ = 0.8 and 61 for γ = 0.2. The bond dimension for the MPS also increases to 120. The time dependent MSD obtained for the second set of parameters is presented in Fig. 1(b). As in the J = 0.1 case, for smaller γ values, it takes a longer time for the MSD to reach the linear growth region. The diffusion constants are also larger for small γ.

We then apply the imaginary and real time HEOM to simulate the current autocorrelation function for the aforementioned two sets of parameters. The MPS method is applied to propagate an imaginary time HEOM, with a bond dimension of 90 for all different γs. The cosine fitting scheme with seven cosine terms is used to fit the imaginary time bath correlation function.⁴⁹ To obtain converged results, the numbers of sites used in the simulations are from 8 (γ = 0.8) to 10 (γ = 0.05) for J = 0.1, and 10 (γ = 0.8) to 26 (γ = 0.05) for J = 0.5. Figure S3 in the supplementary material shows the convergence of the current autocorrelation function with respect to different numbers of sites when γ = 0.05 and J = 0.5.

In the J = 0.1 case, the current autocorrelation functions for different values of γ are shown in Fig. 2(a), and the comparison with the second order perturbation is shown in Fig. S4 in the supplementary material. The current autocorrelation functions for J = 0.5 are shown in Figs. 2(b) and S5. In both cases, a faster decay of the correlation function is observed when increasing γ. It is also shown that, for a larger coupling constant J, the current autocorrelation function becomes less oscillatory.

For all values of γ and J, the current autocorrelation function decays to zero, which indicates that a finite mobility can be obtained via the Green–Kubo relation in Eq. (15). The second order perturbation results in Figs. S4 and S5 also show that, although the accuracy of the second order autocorrelation functions depends on specific model parameters, they also decay to zero and lead to a finite approximate mobility.

The mobilities calculated by using the MSD method and the Green–Kubo relation are consistent with each other. Figure 3(a) shows the calculated mobilities via the MSD and the Green–Kubo relation, compared with the second order perturbation theory for J = 0.1. It can be seen that mobility decreases when the dissipation becomes stronger. The second order perturbation actually works very well for large frictions but overestimates the mobility at low friction. The same behavior is observed in the results for J = 0.5 shown in Fig. 3(b), although the overestimation is rather pronounced for small friction (γ = 0.05). It is thus shown that the friction constant γ plays an important role in charge transport. Increasing γ leads to stronger decoherence, which makes the current autocorrelation function decay faster to zero and reduces mobility.

However, when we try to extrapolate the results to the friction-free limit by taking γ → 0, the γ = 0 limit cannot be defined in both the J = 0.1 and 0.5 cases because of the sharp decrease in the mobility for small γs. In particular, when we investigate the asymptotic behavior of the mobility when γ → 0, all curves show power law scaling μ ∝ γ^−α with different α values. We thus conclude that mobility tends to diverge when γ → 0. The divergence of the second order perturbation result is actually easy to understand. As shown in Fig. S6, in the case γ = 0, the second order current autocorrelation function calculated via Eq. (17) a periodic function with a period of 2π/Ω, so the integration in the Green–Kubo relation goes to infinity.

It would then be interesting to investigate how the current autocorrelation functions in the friction-free Holstein model behave beyond the second order perturbation. We now apply the combined imaginary and real time HEOM toward this problem. For the friction-free Holstein model, the imaginary time bath (vibrational) correlation function is expanded analytically by hyperbolic functions.⁴⁹ 

The current autocorrelation functions for the friction-free model with different numbers of sites are shown in Fig. 4. The parameters used in the simulation are λ = 1.0, Ω = 1.0, β = 1.0, γ = 0, and J = 0.5. A similar set of results is shown for the J = 0.1 case in Fig. S7 in the supplementary material. It can be seen from Fig. 4 that the HEOM current autocorrelation function does not converge even for a large number of sites (N = 50). However, compared with the second order result shown in Fig. S6, the HEOM correlation function also does not show significant recurrences. We thus calculate the integrated current autocorrelation functions in Fig. 4 and show the results in Fig. S8, where the long time behavior for N = 40 is shown in the inset. It can be seen that, in the simulation time up to t = 200, the autocorrelation function has not reached zero, and the integrated results still show small oscillations.

Since the oscillations are not significant, we can define an approximate mobility by truncating the integration in Eq. (15) at some large t. For example, the straight line in the inset of Fig. S8 shows the result for the integration truncated at t = 60. The approximate mobility as a function of the number of sites is shown in Fig. S9. Again, we cannot see a trend of convergence with the increase in the number of sites. So, for the parameter regimes explored in this work, based on the characteristics of the friction-free current autocorrelation function and the failure to extrapolate the mobilities in the dissipative model to the zero friction limit, we conclude that a finite mobility may not be obtained for the friction-free one vibrational mode Holstein model.

In summary, we have applied the imaginary and real time HEOM to simulate the MSD and current autocorrelation functions for the 1D single mode Holstein model in the strong electron–phonon coupling regime. The main observation is that, with dissipative vibrational modes, the current autocorrelation function decays to zero such that a finite mobility can be obtained. On the other hand, in the friction-free Holstein model, the current autocorrelation function does not decay to zero, and a finite mobility cannot be defined either from the Green–Kubo relation or via extrapolation from finite friction results. This is consistent with the recent finding that no diffusive behavior can be observed from the MSD of the friction-free Holstein model.⁴⁰ Of course, in realistic systems, dissipation and static disorder⁵⁷ are unavoidable, and a finite mobility should always be obtained.