首先将影响泛函写成另一种容易处理的表达式.定义如下两类厄米算符

{X_μ=1/2(f_μ+f^_μ),

P_μ=1/(2i)(f_μ-f^_μ).(7)

将其带入式(6)中并将Φ_μ(τ,τ')分成两部分

Φ_μ(τ,τ')=Φ⁰_μ(τ,τ')+Φ^cross_μ(τ,τ'),(8)

其中:

Φ⁰_μ(τ,τ')=[X⁽⁺⁾_μ(τ)-X^(-)_μ(τ)]·

[α_μ(τ-τ')X⁽⁺⁾_μ(τ')-α^*_μ(τ-τ')X^(-)_μ(τ')]+

[P⁽⁺⁾_μ(τ)-P^(-)_μ(τ)]·

[α_μ(τ-τ')P⁽⁺⁾_μ(τ')-α^*_μ(τ-τ')P^(-)_μ(τ')],(9)

包含了X_μ、P_μ与自身的耦合; 而另一项

Φ^cross_μ(τ,τ')=[X⁽⁺⁾_μ(τ)-X^(-)_μ(τ)]·

[α^～_μ(τ-τ')P⁽⁺⁾_μ(τ')-α^～^*_μ(τ-τ')P^(-)_μ(τ')]-

[P⁽⁺⁾_μ(τ)-P^(-)_μ(τ)]·

[α^～_μ(τ-τ')X⁽⁺⁾_μ(τ')-α^～^*_μ(τ-τ')X^(-)_μ(τ')].(10)

则包含了两者之间的耦合.其中,

{α_μ(t)=(2n_μ+1)cos(ω_μt)-i·sin(ω_μt),

α^～_μ(t)=(2n_μ+1)sin(ω_μt)+i·cos(ω_μt),(11)

其中,α^*_μ与α^～^*_μ分别是α_μ与α^～_μ的复共轭.在很多情况下(如当f_μ为厄米算符或当总哈密顿量具有平移对称性时),∑_μΦ^cross_μ的值为零,对动力学没有贡献.因此,为了方便,我们将不再考虑Φ^cross_μ.对于Φ^cross_μ不为零的情况,下文中对Φ⁰_μ采用的所有数学技巧均可直接用于处理Φ^cross_μ.

对影响泛函进行完全随机去耦合可能导致动力学方程收敛性很差,只进行部分随机去耦合则可大大减轻这个问题.为此,将关联函数α_μ分成主体部分α^pri_μ和剩余部分α^res_μ:

α_μ(t)=α^pri_μ(t)+α^res_μ(t).(12)

其中,

{α^pri_μ(t)=γ_μcos(ω_μt),

α^res_μ(t)=(2n_μ+1-γ_μ)cos(ω_μt)-i·sin(ω_μt).(13)

这里,0≤γ_μ≤2n_μ+1是一个待定参数,它决定了α_μ的划分方式.由此,我们将Φ⁰_μ划分为两个部分:

Φ⁰_μ(τ,τ')=Φ^pri_μ(τ,τ')+Φ^res_μ(τ,τ').(14)

其中,主体部分为

Φ^pri_μ(τ,τ')=(X⁽⁺⁾_μ(τ)X⁽⁺⁾_μ(τ')+P⁽⁺⁾_μ(τ)·

P⁽⁺⁾_μ(τ')+X^(-)_μ(τ)X^(-)_μ(τ')+

P^(-)_μ(τ)P^(-)_μ(τ'))·α^pri_μ(τ-τ')-

(X^(-)_μ(τ)X⁽⁺⁾_μ(τ')+P^(-)_μ(τ)P⁽⁺⁾_μ(τ'))·

α_μ(τ-τ')-(X⁽⁺⁾_μ(τ)X^(-)_μ(τ')+

P⁽⁺⁾_μ(τ)P^(-)_μ(τ'))·α^*_μ(τ-τ'),(15)

剩余部分为

Φ^res_μ(τ,τ')=(X⁽⁺⁾_μ(τ)X⁽⁺⁾_μ(τ')+P⁽⁺⁾_μ(τ)·

P⁽⁺⁾_μ(τ'))·α^res_μ(τ-τ')+

(X^(-)_μ(τ)X^(-)_μ(τ')+P^(-)_μ(τ)P^(-)_μ(τ'))·

(α^res_μ(τ-τ'))^*.(16)

可以看到,向前、向后两条路径相耦合的项全都在主体部分中,因此若对主体部分进行随机去耦合,不同路径之间就不再有交叉项.具体地,引入一系列含时随机场{ξ^±,L_μ(t)}(+或-分别对应向前或向后路径,L代表X或P),它们的平均值为零,〈ξ^±,L_μ(t)〉=0,并且服从如下统计性质

{〈ξ^+,L_μ(t)ξ^+,L'_μ'(t')〉=δ_LL'δ_μμ'α^pri_μ(t-t'),

〈ξ ^-,L_μ(t)ξ ^-,L'_μ'(t')〉^*=δ_LL'δ_μμ'α^pri_μ(t-t'),

〈ξ ^+,L_μ(t)(ξ ^-,L'_μ'(t'))^*〉=δ_LL'δ_μμ'α^*_μ(t-t').(17)

利用Hubbard-Stratonovich变换^[74-75],我们可以将影响泛函部分随机去耦合并得到

F(t)=M[F_ξ(t)],(18)

其中,M[·]表示对随机场进行平均,F_ξ(t)与随机场有关,

F_ξ(t)=exp[-i∫^t₀∑_μ[ξ ^+,X_μ(τ)X⁽⁺⁾_μ(τ)+

ξ ^+,P_μ(τ)P⁽⁺⁾_μ(τ)]dτ+i∫^t₀∑_μ[(ξ ^-,X_μ(τ))^*·

X^(-)_μ(τ)+(ξ ^-,P_μ(τ))^*P^(-)_μ(τ)]dτ-

∫^t₀∫^τ₀∑_μΦ^res_μ(τ,τ')dτ'dτ].(19)

可以看到,主体部分Φ^pri_μ已经被时间上局域的随机项取代.将式(18)和(19)代入式(4)中,并假设ρ_S(0)=|ψ₀〉〈ψ₀|,可以得到

ρ_S(t)=M[|ψ⁺_ξ(t)〉〈ψ^-_ξ(t)|],(20)

其中,

|ψ^±_ξ(t)〉=U₀(t)U^±_ξ(t)|ψ₀〉(21)

是向前或向后随机波函数,而

U^±_ξ(t)=T₊exp[-i∫^t₀∑_μ(ξ^±,X_μ(τ)X_μ(τ)+

ξ^±,P_μ(τ)P_μ(τ))dτ-∫^t₀∫^τ₀∑_μ(X_μ(τ)·

X_μ(τ')+P_μ(τ)P_μ(τ'))α^res_μ(τ-τ')dτ'dτ].(22)

将式(21)对时间进行求导,便得到随机波函数的演化方程

其中,R^±_ξ的形式表达式为

R^±_ξ(t)=T₊[∫^t₀∑_μ(X_μ(t)X_μ(τ)+

P_μ(t)P_μ(τ))α^res_μ(t-τ)dτ·

exp[-i∫^t₀∑_μ(ξ^±,X_μ(τ)X_μ(τ)+ξ^±,P_μ(τ)

P_μ(τ))dτ-∫^t₀∫^τ₀∑_μ(X_μ(τ)X_μ(τ')+

P_μ(τ)P_μ(τ'))α^res_μ(τ-τ')dτ'dτ]].(24)

R^±_ξ和(U^±_ξ)^-1的复杂形式使得我们无法直接对式(23)进行数值求解.注意到随机去耦合所引入的随机场严格处理了影响泛函的主体部分,可以预见剩余部分对动力学的影响较小,因此可对式(23)中的R^±_ξ(t)(U^±_ξ)^-1(t)项采用微扰近似,将其以H_I进行展开并截断到二阶项,得到

R^±_ξ(t)(U^±_ξ)^-1(t)≈∫^t₀∑_μ(X_μ(t)X_μ(τ)+

P_μ(t)P_μ(τ))α^res_μ(t-τ)dτ.(25)

将式(25)代入式(23)中便得到一个非常利于数值求解的动力学方程

其中,H^±_ξ(t)等于式(26)中方括号内所有项的和.不难发现NMSSE与一般的含时薛定谔方程形式完全相同,唯一的区别在于前者的哈密顿量具有随机性,系统-热库相互作用对系统动力学的影响经由随机场和非马尔可夫项(式(26)中的积分项)包含于H^±_ξ(t)中.因此在式(26)数值求解的过程中,我们可以将时间离散化,先计算t时刻的H^±_ξ,然后用常规的算法如四阶龙格库塔法等将t时刻的向前、向后波函数分别演化至下一时刻,接着再重复这个过程.

式(26)实际上代表了无穷多种微扰形式的NMSSE,因为我们并未指定γ_μ的值.对于目前文献中已被提出的几种NMSSE,γ_μ=0对应于NMQSD^[51];

对应于含时波包扩散方法^[58]; γ_μ=2n_μ+1对应于最近被提出的一种新的NMSSE^[59].理论上,γ_μ的选取对NMSSE的数值表现有非常大的影响:若γ_μ取较小的值,则NMSSE中的非马尔可夫项有可能变得很大,使得微扰近似失效; 若γ_μ取较大的值,则随机场的收敛性有可能变差.经过大量数值测试,在权衡计算精度和计算资源消耗的平衡后,我们发现一个比较令人满意的选择是

其中,a¹_μ,L和a²_μ,L是相互独立的标准正态分布随机变量,〈aⁱ_μ,La^j_μ',L'〉=δ_ijδ_LL'δ_μμ'.容易验证,这种方式生成的随机场确实服从式(17).此外,由于向前和向后随机场完全相等,每一时刻的向前和向后波函数也必然相等,|ψ⁺_ξ(t)〉=|ψ^-_ξ(t)〉,因此在实际计算中我们只需演化其中一条路径即可.在下文中,我们将默认按照上述方式产生随机场,并略去方程中的±号.

在实际计算过程中,需要根据式(29)对随机场进行大量采样,每次采样都对NMSSE进行一次求解得到动力学结果,采样结束后根据式(20)对结果进行平均.由于NMSSE采用了部分微扰近似,最后得到的约化密度矩阵一般不是归一化的,因此通常需要在得到平均结果后进行人为归一化,以保证数值结果的稳定性.

数值求解式(26)需要我们选择一个具体的表象.在很多情况下,简正模和系统-热库相互作用在空间上局域,因此实空间表象是一个很自然的选择.例如,在传统Frenkel激子模型中,电子态哈密顿量在实空间(位点)表象下可写成

H_S=∑_n≠mV_nm|n〉〈m|,(30)

其中,|n〉代表处在第n个位点的分子被激发而其他分子处在基态,V_nm是激子耦合强度.分子振动及激子-声子相互作用采用独立热库模型,

这里,式(3)中的指标{μ}被指标{n,j}代替,对应于第n个位点上第j个简正模,g_nj代表它与电子态|n〉之间的耦合强度,常常通过谱密度函数

J_n(ω)=π∑_jg²_njδ(ω-ω_nj)(33)

进行描述.将式(32)与式(3)相比较可知,此时

f_nj=f^_nj=X_nj=g_nj|n〉〈n|,(34)

且P_nj=0.因此式(8)中的交叉项Φ^cross_nj确实为零,而且H_ξ化简为

H_ξ(t)=H_S+∑_nξ^X_n(t)|n〉〈n|-i∑_n|n〉〈n|·

∫^t₀U₀(τ)|n〉〈n|U₀(-τ)α^res_n(τ)dτ,(35)

其中,ξ^X_n(t)≡∑_jg_njξ^X_nj(t),α^res_n(t)≡∑_jg²_njα^res_nj(t).

假设系统希尔伯特空间的维度(在Frenkel激子模型中等于分子位点的个数)为N,容易看出,由于H_B和H_I的高度空间局域性,式(35)中随机项的计算量与N成正比,而非马尔可夫项的计算量与N³成正比.考虑到H_ξ作用在波函数上的计算量与N²成正比,可以预见,当N比较大时非马尔可夫项的计算最耗时.但是注意到这一项实际上与随机场无关,因此当N很大时,一个可行的策略是事先计算好每个时间点的非马尔可夫项并存储起来,这样便无需在每次随机取样时都重复计算一遍,其计算成本将远小于方程的其他部分.这使得我们可以处理包含几百个电子态的系统.

上述策略虽然大大减小了计算量,但对内存却有很高的要求,所以无法用于处理包含上千个电子态的系统.此外,当H_B和H_I在空间上离域时(例如,Holstein-Peierls模型),随机项的计算量将与N²成正比,而非马尔可夫项的计算量将与N⁴成正比,此时即使只对后者进行单次计算也非常耗时.这种困难产生的根本原因在于计算非马尔可夫项时涉及到多个稠密矩阵的乘法运算.例如在上一小节的例子中,即使f_nj在实空间表象下非常稀疏(只有一个矩阵元非零),但U₀(t)却是稠密矩阵,其计算过程涉及到H_S的本征表象与实空间表象之间的变换,计算量与N³成正比.这启发我们寻找一个特殊的表象,使得每一个矩阵都是稀疏矩阵,从而大大减少非马尔可夫项的计算量.

当哈密顿量具有平移对称性时,k空间(倒易空间)表象便是电子态哈密顿量的本征表象,且由电声耦合引起的电子态的跃迁必须满足晶格动量守恒,因此计算非马尔可夫项所涉及的矩阵均为稀疏矩阵.在这种情况下,相较于实空间而言,k空间表象是一个更好的选择.例如,在k空间表象下,研究载流子输运过程可以考虑如下电子态哈密顿

H_S=∑_k_ka^_ka_k,(36)

其中,a^_k和a_k分别是晶格动量为k的轨道上载流子的产生和湮灭算符,_k是对应的轨道能量,k=(2πm)/N,m=-N/2,-N/2+1,…,N/2-1(这里N代表第一布里渊区内可取k值的个数).晶格振动和载流子-声子相互作用的一般形式可分别写成

此外,式(3)中的指标{μ}被指标{ν,q}代替,对应于第ν个声子谱分支上晶格动量为q的振动,g^ν_kq表征了电声耦合强度.与式(3)相比较可知,此时

f_{ν q}=N^-1/2∑_kg^ν_{k q}a^_ka_k+q.(39)

可以看到,由于电声耦合满足晶格动量守恒,f_{ν q}在k空间表象下是稀疏矩阵.此外,由于布里渊区的对称性,k与k±2π等价、q与q±2π等价,故由H_I的厄米性可得g^ν_{k q}=(g^ν_k+q,q)^*.根据这些性质,容易推得f_{ν q}=f_ν,q+2π=f^_ν,-q.因此,虽然式(8)中单个交叉项Φ^cross_{ν q}不为零,但是Φ^cross_{ν q}+Φ^cross_ν,-q=0,所以它们的总和也为零.

考虑单载流子的输运过程,定义电子态|k〉≡a^_k|vac〉,其中|vac〉代表真空态.将上述哈密顿量代入NMSSE中可得,H_ξ中的随机项为

容易看出,其计算量也与N²成正比.因此,当总哈密顿量具有平移对称性时,在k空间表象下数值求解NMSSE无需占用额外内存便可大大减少计算量,因此可以实现包含上千个电子态的系统的动力学模拟.

此外,在式(37)和(38)中,不同位置局域振动之间的相互作用可以用ω_νq的色散关系完全表示,局域、非局域的电声耦合也在g^ν_kq的函数关系中被统一处理,因此k空间表象不受独立热库近似限制,可以考虑任意形式的线性电声耦合.作为一个特例,独立热库模型对应于无色散的爱因斯坦声子模型,ω_νq=ω_ν,且对于局域的电声耦合,g^ν_kq=g^ν_loc与晶格动量无关,局域电声耦合的谱密度函数为

J_loc(ω)=π∑_ν(g^ν_loc)²δ(ω-ω_ν).(42)

若在独立热库模型中考虑非局域的电声耦合,该电声耦合引起相邻位点电子耦合强度的涨落,则g^ν_{k q}=g^ν_non(e^-ik+e^i(k+q)),其对应的谱密度函数为

J_non(ω)=π∑_ν(g^ν_non)²δ(ω-ω_ν).(43)

k空间下电声耦合的统一性使得我们可以结合NMSSE与密度泛函理论、密度泛函微扰理论等方法,将在k空间下从头算(ab initio)得到的电子能带结构、声子谱以及电声耦合矩阵直接作为NMSSE的输入参数,研究真实无机或有机材料中载流子的超快动力学过程.

通常将能量超过k_BT的激子态称为热激子.在有机光伏材料中,热激子能量弛豫过程对后续发生的电荷分离、多重激子产生等过程具有极大影响.实验上发现热激子能量弛豫是一个先快后慢的“两步”过程^[77-84],且有证据表明其中的快速步可能和激子的相干运动有关^[85],然而人们对其机理的理解仍然有限.最近一项针对热激子弛豫过程的研究对此进行了明确的理论描述,并阐明了多重时间尺度的来源^[86].这里,我们以一维Frenkel激子模型为例,运用NMSSE研究热激子弛豫与动力学相干性的关联,揭示多重时间尺度背后的物理本质.

具体而言,我们采用式(30)～(33)的哈密顿量模拟由200个单体组成的有机分子链,选择分子链中心100个位点参与形成的相干态作为初始激子波包(具体细节见

参考文献[86]).设每个位点的能量相同,且只考虑相邻分子间的激子耦合,将其作为能量参照值,V_n,n±1=1.不同位点的分子振动特性由同一欧姆谱密度函数表征,即J_n(ω)=J(ω)=1/2πξωe^-ω/ω_c,其中ξ和ω_c分别是Kondo参数和截断频率,二者的乘积是重组能λ.这里,我们设ω_c=1,温度为β=10.

图1(a)是热激子平均能量E _-(t)≡tr{H_Sρ_S(t)}的时间演化曲线.可以看到,激子平均能量呈现多指数衰减,且衰减速度与电声耦合强度正相关.通过对实空间下激子在不同位点上的布

图1 200个单体构成的有机分子链的E -(t)(a)和平均相干长度L-(t)(b)的时间演化,退相干时间td以及快速弛豫时间t1随λ的变化曲线(c)[86]
Fig.1 Time evolution of E-(t)(a)and the averaged coherent length L-(t)(b), variation of the dephasing time(td)and the fast energy relaxation time(t1)with respect to λ(c)in a molecular chain consisting of 200 monomers[86]

E^-(t)=E_e+Ae^-k₁t+Be^-k₂t,(44)

其中,A=(((E_d-E_h)k₁+(E_h-E_e)k₂)C)/(k₂-k₁),B=((E_e-E_d)k₁C)/(k₂-k₁),C为激子的初始布居数,k₁和k₂分别为快速退相干和缓慢跳跃过程的速率.容易看出,B始终大于零,而当(k₁)/(k₂)>(E_e-E_h)/(E_d-E_h)时,A值为负.

为了验证上述机理的合理性,我们需要进一步研究激子的退相干动力学,为此,定义相干长度序列L_l(t):

L_l(t)=∑^N+1-l_n=1|〈n|ρ_S(t)|n+l-1〉|,(45)

其中,1≤l≤N,N为位点个数.L_l(t)的大小反映了相隔l-1个分子的两个位点间的相干程度.图1(b)展示了不同重组能下平均相干长度L_-(t)=∑^N_l=1L_l(t)/(N-l+1)随时间的变化,可以看到,λ越大,退相干速度越快,这与相关实验结果一致^[92].此外,图1(c)展示了利用多个指数函数对不同λ下的E_-(t)和L_-(t)进行拟合所得到的快速能量弛豫时间t₁≡1/k₁和退相干时间t_d.可以看到,t_d非常接近t₁,这证实了快速能量弛豫对应于激子退相干过程的假设.

为了证明上述结论在真实有机材料中的有效性,我们进一步选取有机太阳能电池电子给体聚苯并[1,2-b:4,5-b']二噻吩-噻吩并[3,4-c]吡咯-4,6-二酮(PBDTTPD)作为研究对象.PBDTTPD共聚物及其单体的结构如图2所示.计算动力学所需的重组能以及激子耦合强度通过第一性原理计算得到,其中重组能采用振动模贡献法λ=∑_lλ_l=∑_l1/2ω²_lΔQ²_l计算^[93],ω_l是单体第l个简正模的频率,ΔQ_l是对应的单体基态和第一激发态之间的位移.PBDTTPD重组能的振动模分析数值结果见图2.

插图为PBDTTPD的分子及其单体结构示意图.

图2 PBDTTPD振动分辨重组能(因峰值过低,在3 000～3 300 cm-1频率范围内的结果没有画出)[86]
Fig.2 The mode-resolved reorganization energy for PBDTTPD (frequencies between 3 000 and 3 300 cm-1 are omitted because their peaks are too weak to exhibit)[86]

图3 室温(298 K)下PBDTTPD激子布居数演化图(a)和相干长度序列{Ll(t)}演化图(b)[86]
Fig.3 Evolution of the exciton population of PBDTTPD at room temperature(298 K)(a)and the coherence length sequence {Ll(t)}(b)

实验上发现在部分有机材料中载流子呈现很强的局域性^[96-100],这暗示着材料中存在较强的电声耦合,载流子输运服从跳跃机制.然而,光导率检测结果^[101-102]却表明其迁移率随温度升高而下降,载流子的有效质量较小,属于能带型传输机制.针对这个矛盾,许多理论研究提出了可能的解释^[103-110].一部分研究认为载流子传输遵循相干的能带型传输,其局域特征由非局域的电声耦合导致^[107]; 一部分研究认为其传输本质是低频的分子间振动导致的瞬态定域波包的扩散^[110]; 还有一部分研究发现在跳跃模型的基础上考虑了核隧穿效应后也能得到迁移率随温度负相关变化的传输特征^[103].由此可见,关于局域和非局域电声耦合对载流子迁移率的影响,人们还未有一个统一的结论.这里,我们利用k空间下的NMSSE方法在较大跨度参数范围内计算一维周期性有机半导体的迁移率,研究这两种电声耦合作用对载流子迁移率的影响^[70].

我们采用式(36)～(38)的哈密顿量进行研究,其中电子能带结构由_k=2Vcos k表征,V=100 cm^-1作为能量参照值.声子部分采用无色散的爱因斯坦声子模型,局域和非局域电声耦合强度分别由式(42)和(43)描述,并采用如下两套独立的超欧姆谱密度:

J_loc(ω)=(πλ)/2(ω/(ω_c,loc))³exp[-ω/ω_c,loc],(46)

J_non(ω)=(πη)/2(ω/(ω_c,non))³exp[-ω/ω_c,non],(47)

其中,λ和η分别代表两种电声耦合作用的强度,ω_c,loc和ω_c,non是其对应的振动特征频率.利用NMSSE模拟局域载流子的长时间扩散过程便可得到其扩散系数及迁移率.在k空间下,局域的载流子初始态为|ψ₀〉=N^-1/2∑_k|k〉,均方位移通过下式计算

〈q²(t)〉=N^-1l²∑_kk'∑_nn²exp[i(k-k')n]

〈k|ρ_S(t)|k'〉,(48)

其中,实空间下相邻晶格的间距l设为0.38 nm.扩散系数D和迁移率μ可通过公式D=lim_{t→SymboleB@}〈q²(t)〉/2t以及μ=eD/k_BT计算得到.

ω_c,loc=288 cm^-1,ω_c,non=14 cm^-1,T=300 K时,迁移率随非局域电声耦合强度η的变化曲线如图4所示.在室温下,随着η从零开始不断增加,迁移率先是急剧减小,再缓慢回弹,对应的转折流子传输机制的转变.理论上普遍认为,在弱电声耦合区域,载流子遵循能带型相干传输^[106,111],此时非零的电声耦合会引起载流子的非弹性散射,破坏其相干,进而导致迁移率下降.另一方面,许多研究表明非局域电声耦合对载流子有很强的定域作用^{[107-108,112]}.随着非局域电声耦合强度的增加,载流子在空间上越来越定域,波包相干性不断弱化,载流子逐渐转变为跳跃型传输.在这种传输机制下,非局域电声耦合在一定程度上增加了有效电子耦合强度,从而起到协助载流子迁移的作用^[113-116],这解释了为什么在转折点后迁移率逐渐增加.因此,我们把迁移率随η变化的反转现象归因于载流子从能带型到跳跃型传输的转变.

图4 μ随η的变化曲线[70]
Fig.4 Cuves of μ vs η [70]

虚线表示μ∝T^-α的线性拟合结果,相应的α值在虚线旁标出

图5 载流子迁移率随温度变化曲线[70]
Fig.5 Cuves of μ vs T[70]

图6 瞬时扩散率随时间的演化(只考虑非局域电声耦合)[70]
Fig.6 Evolution of the instantaneous diffusivity(only the nonlocal electron-phonon coupling is included)[70]

图7 瞬时扩散率随时间的演化(只考虑局域电声耦合)[70]
Fig.7 Evolution of the instantaneous diffusivity(only the local electron-phonon coupling is included)[70]