Two channel model for optical conductivity of high mobility organic crystals

TL;DR

This paper models the optical conductivity of high mobility organic crystals like Pentacene and Rubrene, revealing that their carrier dynamics can be described by a superposition of diffusive and localized charge behaviors.

Contribution

It introduces an approximation-free Monte Carlo approach to analyze carrier scattering by impurities and phonons in organic crystals, providing new insights into their conductivity mechanisms.

Findings

Carrier dynamics differ from conventional metals.

Optical conductivity can be modeled as a superposition of Drude and Lorentz terms.

The model accurately describes temperature-dependent conductivity.

Abstract

We show that the temperature dependence of conductivity of high mobility organic crystals Pentacene and Rubrene can be quantitatively described in the framework of the model where carriers are scattered by quenched local impurities and interact with phonons by Su-Schrieffer-Hegger (SSH) coupling. Within this model, we present approximation free results for mobility and optical conductivity obtained by world line Monte Carlo, which we generalize to the case of coupling both to phonons and impurities. We find fingerprints of carrier dynamics in these compounds which differ from conventional metals and show that the dynamics of carriers can be described as a superposition of a Drude term representing diffusive mobile particles and a Lorentz term associated with dynamics of localized charges.

Equations

$$H_1=\hslash {\omega }_0\sum _k{a}_k^{†}{a}_k+\frac{1}{\sqrt{N}}\sum _{q,k}{\delta }_k{c}_{q+k}^{†}{c}_q$$

$$H_2=\sum _q\epsilon (q)c_q^{†}{c}_q+\frac{1}{\sqrt{N}}\sum _{q,k}M(k,q)\left(a_{-k}^{†}+a_k\right)c_{q+k}^{†}{c}_q,$$

$$Z=\sum _{|{\varphi }_1⟩}\cdots \sum _{|{\varphi }_L⟩}⟨{\varphi }_1|e^{-(\beta /L)H}|{\varphi }_2⟩⟨{\varphi }_2|e^{-(\beta /L)H}|{\varphi }_3⟩\cdots$$

$$\frac{{\beta }^2}{N{L}^2}M(k_{2i-1},q_{2i-1})M^{\ast }(k_{2i},q_{2i})\times$$

$$\left[(n_k+1)e^{-\hslash {\omega }_0(\beta n_k+|{\tau }_{2i}-{\tau }_{2i-1}|)}+n_k{e}^{-\hslash {\omega }_0(\beta n_k-|{\tau }_{2i}-{\tau }_{2i-1}|)}\right],$$

$$\frac{{\beta }^2}{N{L}^2}M(k_{2i-1},q_{2i-1})M^{\ast }(k_{2i},q_{2i})\frac{\cosh \left[\hslash {\omega }_0\left(\frac{\beta }{2}-|{\tau }_{2i}-{\tau }_{2i-1}|\right)\right]}{\sinh (\beta \hslash {\omega }_0/2)}$$

$$\tilde{W}[q(\tau )]=\frac{n!!m!!}{N^{(n+m)/2}{L}^{n+m}}W[q(\tau )],$$

$$W[q(\tau )]=\exp \left[2t_h\underset{0}{\overset{\beta }{\int }}\cos q(\tau )d\tau \right]\frac{{\beta }^m}{m!!}\frac{{\beta }^n}{n!!}\times$$

$$\prod _{j=1}^{m/2}{\delta }_{k_{2j-1}}{\delta }_{k_{2j}}\prod _{i=1}^{n/2}M(k_{2i-1},q_{2i-1})M^{\ast }(k_{2i},q_{2i})\times$$

$$\frac{\cosh \left[\hslash {\omega }_0\left(\frac{\beta }{2}-|{\tau }_{2i}-{\tau }_{2i-1}|\right)\right]}{\sinh (\beta \hslash {\omega }_0/2)}.$$

$$A(x\to y)R(x\to y)\tilde{W}[x]=A(y\to x)R(y\to x)\tilde{W}[y],$$

$$R(x\to y)W[x]=R(y\to x)W[y].$$

$${\sigma }_D(\omega )=\frac{{\sigma }_D}{1+({\tau }_D\omega {)}^2},$$

$${\sigma }_L(\omega )=\frac{{\sigma }_L(\omega {\tau }_L{)}^2}{(1-\frac{{\omega }^2}{{\omega }_L^2}{)}^2+(\omega {\tau }_L{)}^2}.$$

$$\sigma (\omega )\simeq {\sigma }_D+({\sigma }_L{\tau }_L^2-{\sigma }_D{\tau }_D^2){\omega }^2.$$

$$D(t)=\frac{1}{2}\frac{d\Delta x^2}{dt}=\frac{\hslash }{\pi e^2}\underset{0}{\overset{\infty }{\int }}\frac{\sigma (\omega )\sin (\omega t)}{\tanh (\beta \hslash \omega /2)}d\omega ,$$

$$\sigma (z)=\frac{i}{z}\left(\Pi (z)-\Gamma \right),$$

$$\Gamma =-\frac{1}{\hslash }{\int }_0^{\beta \hslash }ds⟨j(s)j(0)⟩,$$

$$\Pi (z)=-\frac{i}{\hslash }{\int }_0^{\infty }d\tau e^{iz\tau }⟨[j(\tau ),j(0)]⟩.$$

$$\Pi (s)={\int }_{-\infty }^{\infty }d\omega \frac{1}{\pi }\frac{\omega e^{-\omega s}}{1-e^{\beta \hslash \omega }}\mathrm{Re}\sigma (\omega ).$$

$$\underset{0}{\overset{\infty }{\int }}\sigma (\omega )d\omega =-\frac{\pi }{2}{\left(\frac{ea}{\hslash }\right)}^2⟨H_2⟩,$$

$$H_K$$

$$H_{EPI}$$

Full text

I Introduction

Numerous technological applications Hasegawa and Takeya (2009); Gershenson et al. (2006); Brédas et al. (2009); Mei et al. (2017); MEN (2017) of high mobility organic crystals Pentacene and Rubrene attracted considerable interest to its transport properties. However, in spite of intense efforts, there are only few theoretical approaches capable of providing reasonable descriptions of the physical properties of these compounds Troisi and Orlandi (2006); Troisi (2011); Fratini and Ciuchi (2009); Ciuchi and Fratini (2012); Cataudella et al. (2011); Ciuchi et al. (2011); Fratini et al. (2016, 2014, 2017); De Filippis et al. (2015). Moreover, a complete description of experiments is achieved only in the high-temperature range $150\,\mbox{K}<T<350\,\mbox{K}$. In this paper we present a novel model which is able to provide a unified description of the transport properties of these compounds in the whole temperature range. The model is based on the traditionally used Su-Schrieffer-Hegger model, which is generalized in our study by adding a coupling of carriers to quenched local impurities. We solve this model by a novel approximation free modification of the world line Monte Carlo (WLMC) approach, demonstrate agreement with experimental mobility at all temperatures, and give physical interpretation of the optical conductivity (OC) and diffusivity of these compounds.

To review the current theoretical understanding, we start with the pioneering approach by Troisi-Orlandi Troisi and Orlandi (2006); Troisi (2011), using the one dimensional (1D) Su-Schrieffer-Hegger (SSH) model, where optical phonons of very low energies ($\hbar\omega_{0}\sim\,\mbox{few meV}$) are coupled to the electrons characterized by a hopping of the order of $t_{h}\simeq 100\,\mbox{meV}$, so that the relation $t_{h}>k_{B}T\gg\hbar\omega_{0}$ holds between the model parameters.

Indeed, the mobility in these materials is characterized by a pronounced anisotropy (Ref. Gershenson et al. (2006) and references therein), so that a one dimensional model, where carrier hopping in the other two directions is completely neglected, is a good approximation. Note that this is the main approximatoin of our approach, as we extract the mobility using a method which is approximation free within the model Hamiltonian chosen.

Actually a more detailed description of these semiconductors would require a contribution à la Holstein, that is a local electron-phonon coupling, describing the interaction between the charge carriers and the intramolecular deformations. On the other hand, in literature, the characteristic frequency of these vibrations has been estimated to be of the order of 1-2 $t_{h}/\hbar$ Brédas et al. (2009); Coropceanu et al. (2007), much greater than the frequency associated to the intermolecular displacements, modelled by means of the SSH Hamiltonian, that is about 0.05 $t_{h}/\hbar$. The low energy physics, at $k_{B}T$ much less than $t_{h}$, is determined mainly by the coupling between the electrons and the lattice excitations within the SSH model (nonlocal electron-phonon interaction). Note that $t_{h}/k_{B}$ is about 1200 kelvin.

In literature, the 1D model has been studied by using different approximate approaches, reproducing the experimental power law $\mu\sim T^{-\gamma}$ temperature dependence, where the exponent $\gamma$ depends on materials and disorder Gershenson et al. (2006). These approaches are based on: (i) mixed quantum-classical simulations based on the Ehrenfest coupled equations Troisi and Orlandi (2006); Ciuchi et al. (2011); (ii) neglecting vertex corrections in the OC calculation Fratini and Ciuchi (2009); Ciuchi and Fratini (2012); (iii) introducing an ad hoc energy broadening of the system energy levels Cataudella et al. (2011); (iv) using the relaxation time approximation Ciuchi et al. (2011); Fratini et al. (2016, 2014); and (v) methods based on realistic chemical models Landi and Troisi (2018); Landi et al. (2018).

The scenario emerging in these approaches, i.e. the “transient localization scenario”, has been largely adopted to describe these materials Fratini et al. (2017).

It assumes that, on a time scale less than $\tau\approx 1/\omega_{0}$, the electrons tend to localize due to Anderson theorem Anderson (1958) (the lattice fluctuations are considered ”frozen” on this time scale), while for times larger than $\tau>1/\omega_{0}$ the electrons acquire a diffusion.

Recent approximations free study De Filippis et al. (2015) of the 1D SSH Troisi-Orlandi model confirmed power law of mobility temperature dependence, and established ultimate understanding of the high-temperature behavior of mobility.

However, the approaches based solely on the SSH model Troisi and Orlandi (2006); Troisi (2011); Fratini and Ciuchi (2009); De Filippis et al. (2015); Ciuchi et al. (2011) are unable to describe low temperature ($T<180K$) transport of Pentacene and Rubrene because this model does not account for static extrinsic disorder of organic crystals Gershenson et al. (2006); Liu et al. (2017); Matsui et al. (2010); Mishchenko et al. (2012); Di Sante et al. (2017). To resolve this drawback of disorder-free scenario we add the effect of quenched static local impurities into 1D SSH Hamiltonian, to end up with the model giving comprehensive description of the $\omega$-dependent OC $\sigma(\omega,T)$ and mobility of the crystals under study. Rewriting long-established WLMC Kornilovitch (1998) into momentum space, and generalizing it to the interaction of polaron with impurities, we use analytic continuation to real frequencies based on Maximum Entropy Jarrell and Gubernatis (1996) and Stochastic Optimization Mishchenko et al. (2000) methods to obtain OC from the current-current correlation function. Our analysis reveals two contributions to mobility and OC: one is associate with the diffusive motion of a mobile quasiparticles, and the second arises from the excitations of localized carriers.

II Model and methods

The model Hamiltonian is given by $H=H_{1}+H_{2}$, in momentum space

[TABLE]

and

[TABLE]

where operators $c_{k}^{\dagger}$ ($a_{k}^{\dagger}$) represent the electron (phonon) creation operators, Einstein optical phonons have frequency $\omega_{0}$, and ${\delta_{k}}$ are random variables representing the static on site disorder, distributed according to a Gaussian with $\langle\delta_{k}\rangle=0$, and $\langle\delta_{k}\delta_{k^{\prime}}\rangle=\Delta^{2}$ if $k^{\prime}=-k$, zero otherwise. The electron band is $\varepsilon(q)=-2t_{h}\cos(q)$, with $t_{h}$ the hopping parameter, while $M(k,q)=2it_{h}\lambda\left[\sin q-\sin(q+k)\right]$ is the EPI vertex, with $\lambda$ the dimensionless coupling constant.

We study the model with a WLMC approach, where the trace over the state vectors (we work in the thermodynamic limit $N\to\infty$) becomes a sum over all the paths in imaginary times (world lines). A world line is defined by a function $q(\tau)$ of the electron momentum, piecewise constant and periodic in the interval $[0,\beta]$. The function has a finite number of discontinuities (hoppings) due to the interaction with phonons and with disorder. Each hopping is paired with a corresponding one having opposite momentum.

It is convenient to work in the Suzuki-Trotter approximation, in which the operator $e^{-\beta H}$ is written as the product of a finite number $L$ of “time slices”, $\left[e^{-(\beta/L)H}\right]^{L}$, and take the limit $L\to\infty$ at the end. The partition function, after insertion of the appropriate number of identities, can be written as

[TABLE]

The vectors $|\phi_{k}\rangle$ represent the state of the system (electron + phonons) at imaginary time $\tau_{k}=(\beta/L)k$, with $k=1\ldots,L$. The electron kinetic term of the Hamiltonian gives simply $\exp[2t_{h}\int\limits_{0}^{\beta}\!\cos q(\tau)\,d\tau]$. Of course, before taking the limit $L\to\infty$, the integral over $\tau$ has to be considered as a sum over the discrete times $\tau_{i}$. If the world line has $n$ hoppings due to phonons, let us call $\tau_{2i-1}<\tau_{2i}$ the imaginary times at which the hoppings occur, with $i=1,\ldots,n/2$ and $n$ even, $q_{2i-1}$ and $q_{2i}$ the momentum of the electron just before the hopping, and $k_{2i-1}$ and $k_{2i}$ the momentum of the phonon emitted or absorbed, with $k_{2i}=-k_{2i-1}$. While diagonal terms in the Hamiltonian contribute to the product over all the time slices where the electron does not hop, non-diagonal terms like electron-phonon interaction $H_{\text{epi}}$ contribute only when the electron hops. Therefore, for $L\to\infty$, the term $\exp[-(\beta/L)H_{\text{epi}}]$ can be replaced by $1-(\beta/L)H_{\text{epi}}$. The matrix elements of this operator, and of the diagonal operators $\exp[-(\beta/L)\hbar\omega_{0}a_{k}^{\dagger}a_{k}]$ relative to the phonon being emitted or absorbed over the whole interval $[0,\beta]$, are given by

[TABLE]

where $n_{k}$ is the number of phonons before time $\tau_{2i-1}$, and the two terms correnspond respectively to emission or absorption. Summing over the number of phonons, apart from a factor $(1-e^{-\beta\hbar\omega_{0}})^{-1}$ that corresponds to the partition function of free phonons and can be dropped, we have

[TABLE]

Note that the weight of world lines in which more than one phonon with the same momentum is emitted or absorbed vanishes as $N^{-2}$ or faster, and as we will see later is negligible in the thermodynamic limit. Calling $k_{2j-1}$ and $k_{2j}$, with $j=1,\ldots,m/2$ and $m$ even, the hoppings of the electron due to the disorder term of the hamiltonian, the matrix elements over a pair of hoppings can be computed similarly, and give $\frac{\beta^{2}}{NL^{2}}\delta_{k_{2j-1}}\delta_{k_{2j}}$. Putting all together, we obtain that the weight of the world line is given by

[TABLE]

where

[TABLE]

Monte Carlo updates change the values of the hopping momenta, times, and the number of hoppings. When doing the Monte Carlo simulation, we have to impose the detailed balance condition

[TABLE]

where $x$ and $y$ are two different world lines, $A(x\to y)$ is the probability to propose the Monte Carlo move $x\to y$, and $R(x\to y)$ is the probability to accept it. If for example $x$ is a world line with $n$ hoppings due to phonons, and $y$ a world line with $n+2$ hoppings, then $A(x\to y)=N^{-1}L^{-2}$, because we have to choose a value of the momentum transferred and two values of the imaginary times of emission and absorption, and $A(y\to x)=\frac{2}{n+2}$, because we have to choose which pair of hoppings to remove. Therefore, it is not difficult to see that we can choose $R(x\to y)$ such that

[TABLE]

In this way, all the factors $N$ and $L$ have disappeared, and we can take the limit $L\to\infty$ and $N\to\infty$, extracting momenta from the whole interval $[-\pi,\pi]$, and times from $[0,\beta]$. Note that, if a world line contains more than a pair of hoppings with the same value of the absolute momentum, then a factor $N^{-1}$ will still appear in $W[y]$, and the weight will vanish in the thermodynamic limit. This corresponds to the fact that it is not possible to extract two times the same real number from a continuous interval. Therefore, the simulation of the $N\to\infty$ system is actually simpler than that of the finite system, because we have not to deal with different possible pairings of the hoppings. Note also that, despite we are doing a quenched average over the disorder, the factor $\delta_{k_{2j-1}}\delta_{k_{2j}}$ can be replaced by its average $\Delta^{2}$. Indeed, in the thermodynamic limit, one can average the weight of a world line over an arbitrarily small interval around $k_{2j-1}$ and $k_{2j}$, that contains an infinite number of possible momenta. This means in practice that quenched and annealed averages coincide.

Within the approach described in Ref. Mishchenko et al. (2015), we determine mobility $\mu=\lim_{\omega\mapsto 0}\sigma(\omega)/e$ from the OC $\sigma(\omega)$ which is obtained from current-current correlation function (see Supplemental Material sm that includes Refs. ).

III The mobility

Our results reproduce the experimental data for rubrene Gershenson et al. (2006) in the entire range of temperatures $100\,\mbox{K}<T<300\,\mbox{K}$ (Fig. 1), supporting the assumption that the interaction with low energy phonons and static on site disorder are the key ingredients needed to describe the experiments, and showing that a finite coupling with quantum phonons is able to overcome the Anderson localization and provide a non-vanishing mobility Bonča et al. (2018).

At low $T$ the disorder plays the major role, driving the system from metallic to insulating behavior due to localization on impurities. The effect of disorder is very important even in metallic regime at $T>180K$, where a power law $\mu\propto T^{-1}$ is observed. Note that in the present context, where the carrier density is very small, with “metallic regime” we intend to refer to the region where the mobility exhibits the typical behavior of the metals, i.e. power-law temperature dependence: the mobility decreases by increasing the temperature with an exponent depending on the material and the disorder. On the other hand with insulating regime we refer to the temperature window where mobility increases by increasing temperature. Depending on the disorder strength, the apparent exponent of the power law ($\mu\propto T^{-\gamma}$) can vary from $\gamma=1$ (see Fig. 1) to $\gamma=2$ (negligible disorder), which is in agreement with experimental situations in different samples Gershenson et al. (2006).

To identify the transport regime in high mobility organic crystals Pentacene and Rubrene, we estimated the diffusion time $\tau_{D}$ and the diffusion lenght $L_{d}$, using the exact relation $\mu=-ea^{2}\hbar^{-2}\tau_{D}\left<H_{2}\right>$ provided by Mori formalism Mori (1965), and the Einstein equation $\mu=eD/(k_{B}T)$, with $D=L_{d}^{2}/(2\tau_{D})$. Here $e$ is the electron charge, $a$ the lattice constant, $k_{B}$ the Boltzmann constant, and $\left<H_{2}\right>$ is straightforwardly obtained from the WLMC calculations. Figure 2 shows that, in the entire temperature range, the diffusion time $\tau_{D}$ is of the order of $\tau_{0}=\hbar/t_{h}$, which is many times smaller than that in typical metals, and its temperature dependence is very similar to that of mobility, pointing to diffusive behavior. The diffusion length is of the order of the lattice parameter $a$, and even smaller than $a$ at low temperatures (Fig. 2), violating the Ioffe-Regel condition Ioffe and Regel (1960) and showing that the mobility is dominated by a strong incoherent regime. Such small $L_{d}$ corresponds to diffusivity $D=L_{d}^{2}/(2\tau_{D})$ which is two order of magnitude smaller than in typical metals.

IV Optical conductivity: the two channel model

The profound difference of the low- and high-temperature nature of mobility can be highlighted by the analysis of the dynamic OC in the frequency range $0<\omega<{\tilde{\omega}}/2$, where ${\tilde{\omega}}=t_{h}/\hbar=1/\tau_{0}$ and OC can be represented as a linear combination of two contributions (Fig. 3). The first one is a Drude term describing the OC of mobile particles at low frequency

[TABLE]

and the second one is the Lorentz term describing the dynamical conductivity associated with localized charges

[TABLE]

Furthermore, for comparison, we show in Fig. 3 also the fully adiabatic conductivity, whose limit for $\omega\mapsto 0$ is zero, as expected by the Anderson theorem. It is worth noticing that the peak in the exact and fully adiabatic approaches are very close. Actually for low temperature the agreement is perfect, whereas at room temperature the exact calculation peak is located at a slightly higher frequency, and exhibits a smaller intensity. This signals the presence of quantum effects even at this high temperature.

The ratio of the two coefficients, $\sigma_{D}/\sigma_{L}$, quickly increases with the temperature in the insulating phase, and saturates in the metallic phase at a value around $8.0$ (Fig. 4). Indeed, temperature increasing releases carriers trapped at low temperatures and transfers the spectral weight at low frequencies from the Lorentz to the Drude component, which dominates at high temperatures. In the inset of Fig. 4 we plot the spectral weight of the two components in the interval $0\leq\omega\leq{\tilde{\omega}}/2$, together with their sum which remains almost equal as changing the temperature, showing that the exchange of the weights between Drude and Lorentz component is the leading source of the strong temperature dependence of mobility.

Interesting indication of the profound role of disorder even in the high-temperature metallic region follows from the low-frequency expansion of $\sigma(\omega)$ in powers of $\omega^{2}$:

[TABLE]

The sign of the constant before $\omega^{2}$ is given by the balance between the Drude metallic contribution (always negative) and the Lorentz insulating contribution (always positive). For the high mobility organic materials studied in this work it is found that the positive contribution always overcomes that negative (Fig. 4) at all the studied temperatures. Positive curvature around $\omega=0$ is very large in the insulating phase due to trapping of carriers, but it remains positive even at high temperatures, which distinguishes it form typical ordinary metals.

We note that the Lorentz contribution remains nonzero and saturates even at high temperatures where the majority of trapped carriers are released. This persistence of the Lorentz contribution follows from the Feynman polaron scenario, where a carrier interacting with a phonon field is described in an effective way as a charge linked, through a spring, to a massive particle (simulating the phonon field). As a result the system behavior can be decomposed into the free motion of the center of mass of the effective quasiparticle (polaron) and the excitations associated to the relative motion of the pair (charge - massive particle) which are bound oscillations in the harmonic potential due to the spring. In this case one is able to decompose the low frequency behavior of the system into two contributions: the Drude one, associated with the center of mass motion, and the Lorentz contribution describing the bound oscillations of the internal degrees of freedom Shinozuka and Toyozawa (1979). Here we note a substantial change introduced by the disorder into the balance between Drude and Lorentz contributions. In the crystals without disorder,the Lorentz term is negligible at low temperatures and gains a finite weight only at high temperatures De Filippis et al. (2015). To the contrary, the Lorentz term dominates at low temperatures in system with notable disorder, because bound oscillations around real impurities contribute to the Lorents term much more effectively that those around the fictitious Feynman charge.

V Time dependent diffusivity

The role of the Drude and Lorentz component in transport can be highlighted studying their contribution to the time dependent diffusivity, that is related to the $\sigma(\omega)$ by the following expression Kubo (1957); Ciuchi et al. (2011); Fratini et al. (2016)

[TABLE]

where $\beta=1/k_{B}T$, and $\Delta x^{2}=\langle[x(t)-x(0)]^{2}\rangle$ is the mean square displacement of the position operator.

Indeed, the long time behavior is determined solely by the fitted Drude contribution (3) [dashed green curve in Fig. 5], whereas the diffusivity due to the fitted Lorentz-term (4) [dotted magenta line] is seen only at small times because localized particles can not contribute into long-time large distance diffusion. Moreover, the Drude term dominates even at small times, showing a cross-over from a quadratic in time balistic contribution at very short times to a dissipative behavior at long times going through a maximum at a time around $\tau_{0}$. The differences between the full calculations (blue points) and the sum of the fitted Drude and Lorentz contributions (solid red curve) are due to the high frequencies oscillations obtained in the full calculations, that are not reproduced by the sum of Drude and Lorentz contributions. We note that the diffusivity at long times (mobility) is characterized by a scattering time $\tau_{D}\simeq\tau_{0}\simeq 10^{-15}\,\mbox{s}$, which is few times smaller than that typical of good metals. In this regime the Drude term gives rise to a peak in the time depending diffusivity at a time around $\tau_{0}$ (see dashed green curve in Fig. 5) signaling the quantum nature of the transport. Indeed, in the classical case, given by Eq. (6) in the limit $\hbar\to 0$, the time dependent diffusivity would be a monotonically increasing function of time, as shown by the solid black curve in Fig. 5.

It is worth to note that correlations and memory effects have been extensively studied in other contexts, in particular in statistical physics and studies of dieletric polarization relaxation Jonscher (1999).

VI Conclusions

We suggested a model, where transport properties originate from the interplay of two different interactions. Charge carriers are coupled to phonons by SSH-type interaction, and scattered by quenched local impurities. Within the model proposed by Troisi and Orlandi Troisi and Orlandi (2006), our approximation free results provide quantitative theoretical description of the conductivity of high mobility organic crystals Pentacene and Rubrene from the insulating low- to the metallic high-temperature regime. Essential ingredient of our model is the coupling to impurities which was missing in previous approaches. We are not aware of any other results giving theoretical description of transport data for such wide temperature range.

The scenario emerging from our calculations is based on the existence of two additive contributions to the optical conductivity, characterized by two different times. One ($\tau_{D}$) controls the Drude contribution, and is of the order of $\hbar/t_{h}$. It is the time when the diffusive motion sets in, signaling the onset of a fully incoherent motion, whereas for shorter times the motion is of course ballistic. This contribution is attributed to the motion of the center of mass of the relevant quasi-particle. The other contribution is characterized by a longer time scale ($\tau_{L}$), and controls the maximum of the optical conductivity ($\omega_{L}\simeq 1/\tau_{L}$). This time is related to the relative motion of the internal degrees of freedom of the quasi-particle.

We revealed several fingerprints which discern transport in the high mobility organic crystals investigated in this paper, from that in ordinary metals: extremely small scattering times, diffusion length comparable to lattice spacing, and substantial contribution of localized charges to the high temperature optical conductivity.

Acknowledgements.

This work was supported by the ImPACT Program of the Council for Science, Technology and Innovation (Cabinet Office, Government of Japan).

I Optical conductivity

Within the linear response theory, the optical conductivity (OC) Berg et al. (2008), $\sigma(z)$, can be obtained calculating the current-current correlation function:

[TABLE]

where $z$ lies in the complex upper half-plane, $z=\omega+i0^{+}$. The quantity $\Gamma$ is:

[TABLE]

and $\Pi(z)$ represents the current-current correlation function

[TABLE]

In Eq. (9) (Eq. (8)) $j(\tau)$ ($j(s)$) is the real-time (imaginary-time) Heisenberg representation of the current operator, $[$,$]$ denotes the commutator, $\left\langle\right\rangle$ indicates the thermal average, and $\beta=1/k_{B}T$. The real part of OC is related to the imaginary-time current-current correlation function:

[TABLE]

The function $\Pi(s)$ has been calculated by using World Line Monte Carlo Filippis et al. ; De Filippis et al. (2015) methods. The dynamical spectra, then, is extracted from the integral equation, Eq. (10), through two different approaches i) the maximum entropy method Jarrell and Gubernatis (1996) and ii) the stochastic optimization method Mishchenko et al. (2000). In the former case, in particolar, we followed the Bryan’s method choosing, as default model, the OC obtained through exact diagonalization on a lattice of $20$ sites with periodic boundary conditions (at the investigated temperatures the mean free path (MFP) is less than $6a$, so that such a small lattice provides a very good starting point).

II Sum rule for optical conductivity

The optical conductivity obeys the sum rule

[TABLE]

where $\langle H_{2}\rangle$ is the thermal average of $H_{2}=H_{K}+H_{EPI}$, and

[TABLE]

are the kinetic energy of the electron, and the electron-phonon interaction (see main text for the definition of $\varepsilon(q)$ and $M(k,q)$). In Fig. 6 we show the thermal average of $H_{K}$ (red circles), $H_{EPI}$ (blue squares), and the sum $H_{2}$ (magenta triangles). The average $\langle H_{2}\rangle$ is almost constant in the considered temperature range, while there is an exchange between the two components $\langle H_{K}\rangle$ and $\langle H_{EPI}\rangle$.

III Adiabatic approach

In Fig. 7 we compare time dependent diffusivity given by the Monte Carlo full calculation (blue line), with the so-called adiabatic limit (green line), at temperature $T=270\,\mbox{K}$. In the adiabatic approach the quantum phonon dynamics is completely neglected and they are approximated as static lattice deformations Cataudella et al. (2011). As it can be noted the adiabatic result closely follows the full result at short times ($t<7\,\tau_{0}$, where $\tau_{0}$ is $\hbar$ divided by the hopping parameter) but then it deviates tending to zero as it is expected by the Anderson theorem.

We emphasize that, within our approach, the numerical exact solution is characterized by two times, $\tau_{D}$ controlling the Drude contribution, which is always close to $\tau_{0}$, and $\tau_{L}$, controlling the localized Lorentz contribution which assumes larger values ranging from $16.2\,\tau_{0}$ at $T=124\,\mbox{K}$ to $4\,\tau_{0}$ at $T=270\,\mbox{K}$. It is worth noting that the time emerging from Fig. 7, the time when the full calculation deviates from the fully adiabatic results, does not coincide neither with $\tau_{D}$, nor with $\tau_{L}$.