Cotunnelling and polaronic effect in granular systems

TL;DR

This paper presents a theoretical study of conductivity in granular metallic systems, highlighting how polaronic effects due to stray charges influence the transition from variable-range hopping to Arrhenius-like behavior at low temperatures.

Contribution

It introduces a unified theoretical framework that accounts for both weak and strong polaron effects and describes the crossover from elastic to inelastic cotunneling in granular metals.

Findings

Polaronic effects induce Arrhenius-like conductivity at low temperatures.

Effective activation energy depends logarithmically on temperature.

The theory covers both weak and strong polaron regimes and elastic-inelastic cotunneling crossover.

Abstract

We theoretically study the conductivity in arrays of metallic grains due to the variable-range multiple cotunneling of electrons with short-range (screened) Coulomb interaction. The system is supposed to be coupled to random stray charges in the dielectric matrix that are only loosely bounded to their spatial positions by elastic forces. The flexibility of the stray charges gives rise to a polaronic effect, which leads to the onset of Arrhenius-like conductivity behavior at low temperatures, replacing conventional Mott variable-range hopping. The effective activation energy logarithmically depends on temperature due to fluctuations of the polaron barrier heights. We present the unified theory that covers both weak and strong polaron effect regimes of hopping in granular metals and describes the crossover from elastic to inelastic cotunneling.

Figures

Equations

$$\sigma \propto \exp \{-(T_0/T{)}^{\alpha }\}$$

$$T_0\sim \mathcal{L}(T)E_c$$

$$\mathcal{L}(T)\sim \{\begin{array}{rl}\ln \left(E_c/g\delta \right),& T\ll T_{c0}\\ \ln \left(E_c^2/g{T}^2{\mathcal{L}}^2\right),& T\gg T_{c0}\end{array}$$

$$T\sim T_{c0}=\sqrt{E_C\delta }/\mathcal{L}$$

$$\sigma \sim \exp \{-E_H/T\}$$

$$E_C(\vec{N},\vec{q})=\frac{1}{2}\left\{(\vec{N}-\vec{Q})\hat{U}(\vec{N}-\vec{Q})-\vec{Q}\hat{U}\vec{Q}\right\}$$

$${\hat{E}}_{\mathrm{dev}}=\frac{1}{2}(\vec{Q}-{\vec{Q}}^{(0)})\hat{K}(\vec{Q}-{\vec{Q}}^{(0)})$$

$${\hat{\mathcal{H}}}_{\mathrm{kin}}=-\frac{{\hslash }^2}{2}{\partial }_{\vec{Q}}{\hat{M}}^{-1}{\partial }_{\vec{Q}}$$

$$K_{i{i}^{\prime }}(\vec{N})=K_{i{i}^{\prime }}^{(0)}+\sum _j{K}_{i{i}^{\prime }j}^{(1)}{N}_j+\sum _{j{j}^{\prime }}{K}_{i{i}^{\prime }j{j}^{\prime }}^{(2)}{N}_j{N}_{j^{\prime }}$$

$$K_{i{i}^{\prime }}\equiv K_{i{i}^{\prime }}^{(0)}+\sum _j{K}_{i{i}^{\prime }j}^{(1)}{N}_j^{(\mathrm{eq})}+\sum _{j{j}^{\prime }}{K}_{i{i}^{\prime }j{j}^{\prime }}^{(2)}{N}_j^{(\mathrm{eq})}{N}_{j^{\prime }}^{(\mathrm{eq})}$$

$$E_{\mathrm{ch}}(\vec{N},\vec{Q})=E_C(\vec{N},\vec{Q})+E_{\mathrm{dev}}(\vec{Q})=$$

$$=\frac{1}{2}\left\{(\vec{N}-\vec{Q})\hat{U}(\vec{N}-\vec{Q})-\vec{Q}\hat{U}\vec{Q}\right\}+$$

$$+\frac{1}{2}(\vec{Q}-{\vec{Q}}^{(0)})\hat{K}(\vec{Q}-{\vec{Q}}^{(0)})$$

$${\vec{Q}}_{\min }(\vec{N})={\vec{Q}}^{(0)}+{\hat{K}}^{-1}\hat{U}\vec{N},$$

$$\underset{\vec{Q}}{\min }{E}_{\mathrm{ch}}(\vec{Q},\vec{N})=\frac{1}{2}\vec{N}\hat{\mathcal{U}}\vec{N}-\vec{N}\hat{U}{\vec{Q}}^{(0)}$$

$$\hat{\mathcal{U}}=\hat{U}-\hat{U}{\hat{K}}^{-1}\hat{U}$$

$${\vec{N}}^{(\mathrm{eq})}={\hat{\mathcal{U}}}^{-1}\hat{U}{\vec{Q}}^{(0)}+\vec{\gamma },$$

$${\vec{Q}}^{(\mathrm{eq})}={\vec{Q}}^{(0)}+{\hat{K}}^{-1}\hat{U}{\vec{N}}^{(\mathrm{eq})}$$

$$-1/2\le {\gamma }_k\le 1/2$$

$$\vec{N}={\vec{N}}^{(\mathrm{eq})}+\vec{n},\vec{Q}={\vec{Q}}^{(\mathrm{eq})}+\vec{q}$$

$$E_{\mathrm{ch}}=\frac{1}{2}\left\{(\vec{n}-\vec{q})\hat{U}(\vec{n}-\vec{q})+\vec{q}(\hat{K}-\hat{U})\vec{q}\right\}+\vec{n}\hat{\mathcal{U}}\vec{\gamma }$$

$$\tilde{E}(\vec{n})\equiv \underset{\vec{q}}{\min }{E}_{\mathrm{ch}}(\vec{n},\vec{q})=\frac{1}{2}\vec{n}\hat{\mathcal{U}}(\vec{n}+2\vec{\gamma })$$

$${\tilde{E}}_k^{(\pm )}=\frac{1}{2}{\vec{z}}_k\hat{\mathcal{U}}({\vec{z}}_k\pm 2\vec{\gamma })$$

$$\displaystyle\varepsilon_{k}=\left\{\begin{aligned} \tilde{E}_{k}^{(+)},\qquad&\mbox{for $\vec{z}_{k}\hat{\cal U}\vec{\gamma}<0$}\\ -\tilde{E}_{k}^{(-)},\qquad&\mbox{for $\vec{z}_{k}\hat{\cal U}\vec{\gamma}>0$}\end{aligned}\right.$$

$${\epsilon }_{lr}\equiv \frac{1}{2}\left\{|{\epsilon }_l-{\epsilon }_r|+|{\epsilon }_l|+|{\epsilon }_r|\right\}$$

$$E_k^{(\pm )}(\vec{q})\equiv E_{\mathrm{ch}}(\pm {\vec{z}}_k,\vec{q})-E_{\mathrm{ch}}(0,\vec{q})=$$

$$={\tilde{E}}_k^{(\pm )}+\frac{1}{2}{\vec{z}}_k\hat{U}{\hat{K}}^{-1}\hat{U}{\vec{z}}_k\mp {\vec{z}}_k\hat{U}\vec{q}$$

$${\tilde{E}}_{lr}^{(\pm \pm )}={\tilde{E}}_l^{(\pm )}+{\tilde{E}}_r^{(\pm )}$$

$${\tilde{E}}_{lr}^{(\pm \mp )}={\tilde{E}}_l^{(\pm )}+{\tilde{E}}_r^{(\mp )}-{\vec{z}}_l\hat{\mathcal{U}}{\vec{z}}_r$$

$$E_{lr}^{(\pm \pm )}(\vec{q})\equiv E_{\mathrm{ch}}(\pm {\vec{z}}_l\pm {\vec{z}}_r,\vec{q})-E_{\mathrm{ch}}(0,\vec{q})=$$

$$=E_l^{(\pm )}+E_r^{(\pm )}$$

$$E_{lr}^{(\pm \mp )}(\vec{q})\equiv E_{\mathrm{ch}}(\pm {\vec{z}}_l\mp {\vec{z}}_r,\vec{q})-E_{\mathrm{ch}}(0,\vec{q})=$$

$$=E_l^{(\pm )}+E_r^{(\mp )}-{\vec{z}}_l\hat{U}{\hat{K}}^{-1}\hat{U}{\vec{z}}_r$$

$$W=\frac{1}{8}({\vec{n}}_f-{\vec{n}}_i)\hat{U}{\hat{K}}^{-1}\hat{U}({\vec{n}}_f-{\vec{n}}_i)$$

$${\mathcal{E}}_{\mathrm{act}}=\frac{{\tilde{E}}_i+{\tilde{E}}_f}{2}+W+\frac{({\tilde{E}}_f-{\tilde{E}}_i{)}^2}{16W}$$

Full text

I Introduction

In this paper we are discussing a multiparticle cotunneling mechanism of conductivity in a granular metal with flexible charges, randomly placed in the insulating matrix. The flexibility of disorder gives rise to a sort of “random polaronic effect”, which dramatically affects the temperature dependence of conductivity. A study of this subject requires a combination of different physical concepts and methods, so we start from a brief review of the necessary ingredients.

I.1 Variable range hopping

The stretched-exponential temperature dependence of conductivity

[TABLE]

is characteristic of variable range hopping (VRH) in homogeneously disordered materials (such as amorphous solids mott or doped semiconductors ES_book ) at low temperatures. In the case when the long-range interactions do not play any significant role, the exponent $\alpha=1/(d+1)$ (where $d$ is the space dimensionality), $T_{0}\equiv T_{\rm M}\sim[\nu_{F}\xi^{d}]^{-1}$ (where $\nu_{F}$ is the density of states at the Fermi level and $\xi$ is the inverse decrement of electronic wave-functions), and the corresponding dependence is known as the Mott law mott0 . In the opposite case, when the long range Coulomb interaction is crucial and gives rise to the soft Coulomb gap in the electronic density of states ES_gap , the exponent $\alpha=1/2$, $T_{0}\equiv T_{\rm ES}\sim e^{2}/\kappa\xi$ (where $\kappa$ is the dielectric constant) and this dependence is known as the Efros-Shklovskii law.

I.2 Variable range cotunneling in granular systems

A similar behaviour of conductivity was also observed in arrays of metallic and semiconducting quantum dots in the temperature range 1 – 200K (see reviews abeles ; beloborodov and some old VRHold and more recent VRH experimental papers). An explanation of the stretched-exponential $T$-dependence of conductivity in granular materials attracted interest of theorists. Many of early theories earlytheory1 were based on special assumptions about the distribution of the random parameters of grains (sizes etc) and were criticised – see, e.g., PollakAdkins92 , because of their ad hoc character and the lack of universality. Some other theories ZhangShklovskii ; earlytheory2 correctly indicated the important role of sequential tunneling of electrons through a chain of intermediate grains, but did not give a correct multi-particle description for this tunneling and a valid recipe for evaluation of the tunneling amplitude. Such a prescription was worked out in iosel2005 , where the simple idea, introduced in ZhangShklovskii , was generalised to take into account both multi-particle character (i.e., the cotunneling, see Nazarov1990 ; NazarovBlanter ) of the process and Coulomb effects.

In contrast with the standard single-particle-like VRH scenario, where a particle would travel from one end of the chain to the other, consequently hopping through all intermediate grains, the multiple cotunneling scenario, developed in iosel2005 involves all possible sequences of hops between neighbouring grains in the chain. In general, these hops are executed by different electrons; the intermediate states of the process involve many grains with altered charges. However, upon the completion of the process there is only one net electron that is transferred between the terminal grains of the chain, while the charges of all intermediate grains return to their initial values. It does not mean that all the processes with different sequences of hops lead to the same final state of the system: the final states of some grains may be identical to the initial ones (elastic cotunneling), while for other grains the initial and final states may differ by an electron-hole pair (inelastic cotunneling). As a result (see iosel2005 for details), the law (1) is reproduced with

[TABLE]

where

[TABLE]

and the Coulomb charging energy is $E_{c}\sim e^{2}/2C$, $C$ being typical capacitance in the system of grains. A typical level spacing in a grain $\delta\sim(\nu_{F}a^{3})^{-1}$, $a$ being the grains size, and $g\ll 1$ is the typical dimensionless conductance between adjacent grains. Note that the logarithmic factor ${\cal L}$ is large. The crossover between the elastic and inelastic cotunneling regimes takes place at

[TABLE]

Thus, the stretched exponential law in the case of granular materials is slightly modified in the intermediate temperature range ($T>T_{c0}$) due to additional logarithmic $T$-dependence of $T_{0}$. This deviation from the Mott-Efros-Shklovskii law is, however, not easy to detect experimentally.

I.3 Hard gap and polaronic effect in homogeneously disordered systems

It is well known that at relatively high temperatures the VRH is not operative, it is changed to the nearest neighbour hopping (NNH), so that the stretched-exponential law (1) is replaced by the Arrhenius law for conductivity $\sigma\propto\exp\{-\varepsilon_{3}/T\}$ (see ES_book ). What is much less trivial, in some cases hard-gap the reentrance of the Arrhenius law

[TABLE]

is observed also at low temperatures! This reentrance is usually attributed to polaronic effect: a “hard gap” $E_{H}$ is supposed to be related to the energy, necessary for the creation in advance of a polaronic cloud, that then will accommodate a hopping electron at a new position.

In principle, the hopping electron can take along the necessary energy while hopping from the initial position to the new one, therefore at still lower temperatures the activational mechanism of polaron hopping is substituted by the tunneling one (see mott ) and the stretched exponential law (1) is again restored at $T\lesssim\omega$, where $\omega$ is the characteristic frequency of phonons (or some other species – magnons, localised electronic excitations etc), that constitute the polaronic cloud. The case of magnetic polarons is special: due to the local conservation of the magnetisation the process of the tunneling transfer of the polaronic cloud is strongly suppressed MHG , and the classical-quantum crossover is shifted from $T\sim\omega$ to much lower temperatures. It explains why the hard gap phenomenon was experimentally observed predominantly (but not exclusively!) in magnetic systems. In general, the theory of VRH with account for polaronic effect was developed in MHG ; polaron-theory ; iossel86 ; foygel for different types of polarons.

It is important that, if the strength of the polaronic effect randomly varies from place to place, then it doesn’t necessarily have to lead to the activation Arrhenius law. If the distribution of barrier heights has a power-law tail at zero, one can expect the dependence (1) with $\alpha\neq 1$. In particular, for homogeneously disordered solids it was shown in foygel that, if the barrier distribution is constant in the vicinity of zero, the Mott conductivity should have the exponent $\alpha=2/(d+2)$.

I.4 Polaronic effect in granular systems

The onset of Arrhenius-like behaviour of conductivity at low temperatures in granular materials is reported far less often. For example, in experiments polaron exp the Arrhenius behaviour was observed below $\sim 5K$ in two-dimensional arrays of semiconducting Ge/Si quantum dots. Above that temperature the conductivity followed Efros-Shklovskii law. In paper polaron exp 2 the observation of “almost” Arrhenius conductivity was reported in the temperature range $20K<T<30K$ in the array of metallic Co nanoparticles. This dependence was well described by (1) with $\alpha=1.1$. At higher temperatures $45K<T<80K$ Efros-Shklovskii law was found. There are also experiments polaron exp 3 in which $\alpha=2/3$ was observed in the broad temperature range $7K<T<200K$ in the arrays of ZnO nanocrystals.

It is tempting to attribute these findings to some kind of random polaronic effect.

But what does the polaron effect mean in the case, where we deal not with single charge carriers, but with grains, containing many electrons? And what are the effective degrees of freedom that constitute here a polaronic cloud? As to our knowledge, no models of the polaronic effect in granular systems have been discussed so far. In this paper we introduce such a model for metallic grains and study its implications for the transport properties within the framework of multiple cotunnelling concept.

I.5 Structure of the article

The article is organized in the following way. We discuss the basic concepts of our model in Sections II and III. The detailed calculation of electron transition rate between the pair of distant resonant grains in the case of short-range interaction is carried out in Sections IV-VII. We deal with different parameter ranges and also provide physical interpretation of results. Based on these findings we analyse the temperature behaviour of conductivity in Sections VIII and IX. Section X contains the summary of our results and in Section XI we discuss the limitations and possible future directions of research.

II Dynamical fluctuations of the offset charges: General Model

The main source of disorder in granular systems is the “stray charges” – hardly removable charged impurities and defects, trapped in the insulating part of the system, see Fig. 1. They produce random Coulomb fields acting on the grains, so that the Coulomb energy of the system is

[TABLE]

where we have introduced vectorial notation $\vec{N}\hat{U}\vec{N}\equiv\sum_{ii^{\prime}}U_{ii^{\prime}}N_{i}N_{i^{\prime}}$ Here integer $N_{i}$ denotes a number of excess electrons on $i$-th grain, $\hat{U}\equiv e^{2}\hat{C}^{-1}$ is the inverse matrix of capacitances, and the components $Q_{i}$ of the vector $\vec{Q}$ are the so called “offset charges” (not necessarily integers!). It should be noted that each offset charge $Q_{i}$ can not be identified with certain unique impurity: all impurities that effectively interact with a given grain $i$ contribute to $Q_{i}$. Vice versa, each impurity may effectively contribute to many different variables $Q_{i}$.

In the context of VRH the offset charges are usually treated as static random variables, but in this paper we are going to take into account their dynamics. Indeed, $Q_{i}$ depend on the positions of the charged impurities, that are not absolutely rigid, but can deviate from their equilibrium places. In the harmonic approximation these deviations are governed by the hamiltonian $\hat{H}_{\rm dev}=\hat{E}_{\rm dev}+\hat{\cal H}_{\rm kin}$,

[TABLE]

where the positively defined matrix of effective masses $\hat{M}$ is related to the masses of impurities. The vector $\vec{Q}^{(0)}$ describes the set of equilibrium values of offset charges for the case of neutral grains: $\vec{N}\equiv 0$. The matrix $\hat{K}$ in (7) is related to the stiffness of the system with respect to displacements of the charged impurities. It contains both the “mechanical” part (due to deformation of surrounding medium) and the “electrostatic” part (the variation of the electrostatic energy of grains (6) due to the displacement of impurities). The second part depends on the set $\vec{N}$; namely, since $K_{ii^{\prime}}(\vec{N})$ is linearly related to the second derivative of the energy with respect to the coordinates of impurities, $\hat{K}$ should be a quadratic polynomial in $\vec{N}$:

[TABLE]

The mechanical part contributes only to the first, $N$-independent, term in (9). Therefore, in the most natural case, when the mechanical stiffness dominates over the electrostatic one, the $N$-dependence of $\hat{K}$ can be neglected. Anyway, even if the electrostatic part of $\hat{K}$ is considerable, at low temperatures it seems reasonable to ignore the “live” $\vec{N}$-dependence of $\hat{K}$, and replace the function $K_{ii^{\prime}}(\vec{N})$ by its equilibrium value

[TABLE]

where $\vec{N}^{\rm(eq)}$ are the charges, that the grains acquire at the equilibrium.

II.1 The classical ground state

Since the effective masses are large, in the leading approximation the kinetic energy term (8) in the hamiltonian can be neglected, and the ground state of the system corresponds to the minimum of the total potential energy of the system of charges

[TABLE]

Minimising (11) with respect to $\vec{Q}$ at fixed $\vec{N}$, we get

[TABLE]

where the symmetric matrix

[TABLE]

has a meaning of the inverse capacitance matrix, renormalised due to the effects of finite elasticity of the system.

For the stability of the system’s ground state the matrix $\hat{\cal U}$ must be positively defined. A violation of this requirement would mean that we have incorrectly chosen the ground state set $\vec{N}^{\rm(eq)}$, which turned out to be unstable; the system eventually will move to a different state, where the stability will be restored due to nonlinearity of the problem, expressed in the $\hat{K}(\vec{N})$ dependence (9).

Further minimisation of (13) with respect to $\vec{N}$ gives the equilibrium values of charges as integers, closest (in a sense, see below) to $\hat{\cal U}^{-1}\hat{U}\vec{Q}^{(0)}$. As a result

[TABLE]

where $\vec{\gamma}$ is the set of “effective residual offset charges”. Note, that in the case of rigid impurities (when $\hat{K}\to\infty$) $\hat{\cal U}\to\hat{U}$, and the effective offset charges are reduced to the standard ones: $\vec{\gamma}\to$ noninteger part of $\vec{Q}^{(0)}$

The common restriction, usually imposed on the residual offset charges, reads

[TABLE]

Strictly speaking, this is not correct in general case: for symmetric matrix $\hat{\cal U}$ the all-integer-$\vec{N}$ minimum of the energy functional may, in principle, lay quite far from the unrestricted one, so that some components of $\vec{\gamma}$ may be quite large. Examples are easy to produce (say, a highly anisotropic potential profile with valleys, looking in low-symmetry directions) but all these examples seem to be exotic, if not pathological. At least we were not able to construct any physically relevant matrix $\hat{\cal U}$ for that the restriction (17) would be violated. Anyway, it is definitely valid for the case which we are going to study in detail below: the screened Coulomb interaction with diagonal matrix $\hat{\cal U}$. Therefore, in what follows we will consider the restriction (17) granted.

II.2 Low-energy hamiltonian

At low temperatures both the grains’ charges $\vec{N}$ and the offset charges $\vec{Q}$ only slightly deviate from the equilibrium values, so one can write

[TABLE]

where $n_{k}$ may take values $(-1,0,1)$. Now we are prepared to rewrite the hamiltonian in terms of deviations $\vec{n}$ and $\vec{q}$. Substituting (18) into (11) and omitting the terms, that do not contain deviations, we get

[TABLE]

II.3 Thermodynamic excitation energies

Besides the ground state, the low-lying excited states are of great importance for the transport properties of the system. For the state $\vec{n}$ the thermodynamic (i.e., minimized with respect to $\vec{q}$) excitation energy is

[TABLE]

Different branches of spectrum (19) can be visualized as multi-dimensional paraboloids in the $\vec{q}$-space. The paraboloids are indexed by vector $\vec{n}$. Excitation energy $\tilde{E}(\vec{n})$ is nothing else, but the energetic distance between the bottom of corresponding paraboloid and the global ground state energy, see Fig.2. For brevity we have denoted $\tilde{E}_{i}\equiv\tilde{E}(\vec{n}_{i})$ and $\tilde{E}_{f}\equiv\tilde{E}(\vec{n}_{f})$.

II.3.1 Single-particle excitations

Below we will obtain some crucial characteristics of single-particle excited states, in which only one entry $n_{k}$ in $\vec{n}$ is nonzero, while all other entries are zeroes (namely, $\vec{n}=\pm\vec{z}_{k}$ where $[\vec{z}_{k}]_{m}=\delta_{km}$).

The thermodynamic excitation energies

[TABLE]

correspond to variation of the “relaxed” energy due to creation (annihilation) of one electron at grain $k$. By definition, the inequality $\tilde{E}_{k}^{(\pm)}\geq 0$ should hold for all $k$, which is ensured by (17). In granular systems it is convenient to introduce

[TABLE]

that has the meaning of the energy of “charged ground state”, counted from the global ground state of the grain. Density of such states is sometimes called the density of ground states (DOGS) in the literature – see ZhangShklovskii ; nanocrystals for more details.

As we will see soon, another useful combination that enters the activation exponent of the conductance between two distant grains $l$ (by agreement “left”) and $r$ (by agreement “right”) is

[TABLE]

Note that it has a conventional analogue in the standard hopping conductivity theory ES_book .

In the following we will assume that electronic transitions proceed very fast, so that the slow variables $\vec{q}$ do not have time to change during the process – the Franck-Condon principle (see more discussion on this topic in Section IV.1). Thus we introduce additional Franck-Condon excitation energies

[TABLE]

II.3.2 Two-particle excitations

There are four classes of possible two-particle processes, corresponding to an act of the charge $-e$ transfer from grain $l$ to grain $r$, see Fig. 3.

1. a)

($++$ process): Transfer of an electron-type single-particle excitation from $l$ to $r$: $(\vec{n}=\vec{z}_{l})\longrightarrow(\vec{n}=\vec{z}_{r})$; 2. b)

($--$ process): Transfer of a hole-type single-particle excitation from $r$ to $l$: $(\vec{n}=-\vec{z}_{r})\longrightarrow(\vec{n}=-\vec{z}_{l})$; 3. c)

($+-$ process): Annihilation a two-particle excitation, consisting of an electron-type excitation in grain $l$ and a hole-type one in grain $r$: $(\vec{n}=\vec{z}_{l}-\vec{z}_{r})\longrightarrow(\vec{n}=0)$; 4. d)

($-+$ process): Creation of a two-particle excitation, consisting of an electron-type excitation in grain $r$ and a hole-type one in grain $l$: $(\vec{n}=0)\longrightarrow(\vec{n}=\vec{z}_{r}-\vec{z}_{l})$;

While in the first two processes only single-particle excitations are involved, in the third and fourth ones the two-particle complexes (intergrain electron-hole pairs) are created or annihilated. Because of generally long-range character of interaction matrix $\hat{\cal U}$ the components of two-particle excitations interact with each other, so that their energies generally are not additive. Thermodynamic excitation energies for the two-particle excitations can be obtained from (20):

[TABLE]

and the Franck-Condon energies are

[TABLE]

We will see in Section III that the interaction parts in (26) and (28) vanish in the case of short-range interaction.

II.3.3 Potential barrier between resonant grains and activation energy

A very important role in low temperature physics is played by the resonant grains, for which either $\tilde{E}_{k}^{(+)}$ or $\tilde{E}_{k}^{(-)}$ is anomalously close to zero. The transition between such resonant states $\vec{n}_{i}=\vec{z}_{i}$ and $\vec{n}_{f}=\vec{z}_{f}$ requires, however, a considerable change of the surrounding (i.e., the vector $\vec{q}$), which can only be done continuously. In the course of this change the potential energy of the system also changes – first increases, then decreases, so that the system has to overcome the potential barrier. This can be accomplished either by means of activation over the barrier, or by tunneling. For both processes the height of the barrier $W$ is crucial. To find it we should minimise the energy $E_{\rm ch}(\vec{n}_{i},\vec{q})$ over $\vec{q}$ with additional condition $E_{\rm ch}(\vec{n}_{i},\vec{q})=E_{\rm ch}(\vec{n}_{f},\vec{q})$, implying that the states are resonant. The result is

[TABLE]

Note that $W$ depends only on the difference $\vec{n}_{f}-\vec{n}_{i}$ between the final and initial states. In particular, it is the same for all kinds of processes described in Section II.3.2.

We can also find the activation energy ${\cal E}_{\rm act}$ for a transition between arbitrary (not necessarily resonant) states $\vec{n}_{i}$ and $\vec{n}_{f}$. It can be defined as the lowest point of interception of two paraboloids (see Fig. 2):

[TABLE]

where $W$ is defined by (29). If the states are resonant, we have ${\cal E}_{\rm act}=W$.

III Short range interaction model

Those granular system, where the long range part of interaction is screened (say, because of the presence of metallic gate, or due to residual conductivity of the insulating matrix) can be roughly described by the simplest model of “short range Coulomb interaction” (see, e.g., NazarovBlanter ). In this model we assume the matrices $\hat{U}$ and $\hat{K}$ to be diagonal

[TABLE]

but their diagonal entries in general are not identical, since different grains have different capacitances etc. The electrostatic energy (19) can be rewritten in a simple way $E_{\rm ch}=~{}\sum_{i}E_{i}(n_{i},q_{i})$, where

[TABLE]

and each grain is characterised by three constants:

• •

The standard charging energy $E^{c}_{i}\equiv U_{i}/2=e^{2}/2C_{i}$ with $C_{i}$ being the capacitance of the grain $i$,

• •

The “polaronic coupling constant” $\alpha_{i}\equiv U_{i}/K_{i}$. The inequality $\alpha_{i}<1$ is the stability condition.

• •

The random “effective offset charge” $\gamma_{i}$, distributed in the interval $[-1/2,1/2]$.

The thermodynamic and Franck-Condon excitation energies for this model are

[TABLE]

IV Transition rate

As we know from iosel2005 , the main features of variable-range hopping in granular systems can be revealed already in the simplest model of short range Coulomb interaction, described in Section III, so in this paper we restrict our consideration to this model.

IV.1 Franck-Condon principle

We will also ignore the quantum aspects of the offset charges dynamics (i.e., put $M\to\infty$). In particular, the last assumption means that we can use the Franck-Condon principle in the calculation of transition rates. According to this principle, the set $\vec{q}$ remains unchanged on the cotunneling time scale $\sim E_{c}^{-1}$. We should calculate the transition rates $w_{ii^{\prime}}(\vec{q}\;)$ between relevant grains at fixed configuration of $\vec{q}$ and only afterwards perform the thermodynamic averaging of the result with respect to $\vec{q}$. The cotunneling time $\tau_{\rm cotun}$ can be roughly estimated as the inverse scale of energetic denominators appearing in the perturbation theory

[TABLE]

On the other hand, characteristic time scale $\tau_{q}$ for the dynamics of $\vec{q}$ is nothing but the inverse frequency of the impurity oscillations $\omega_{\rm imp}$, which apparently should be of the order of phonon frequencies $\omega_{\rm ph}$. As a result, the applicability criterion for the Franck-Condon principle $\tau_{\rm cotun}\ll\tau_{q}$ becomes

[TABLE]

Note that in typical granular systems the rough estimate $\sim$ 100 K holds for $E_{c}$ as well as for $\omega_{\rm imp}$. Thus, the condition (36) can not be considered automarically fulfilled, the opposite case is well probable. However, as we already noted, we are going to neglect the kinetic term for $\vec{q}$ in the hamiltonian, which means that $\vec{q}$ can overcome potential barriers only by activation, and not by tunneling. This implies even more stringent condition

[TABLE]

that we will consider granted throughout this paper.

IV.2 Hamiltonian and some qualitative considerations

Under the above approximations, the hamiltonian of the system may be written in the form

[TABLE]

where $E_{i}(n_{i},q_{i})$ is given by expression (31), and

[TABLE]

is the hamiltonian of electrons within the $i$-th grain. The index $\lambda_{i}$ denotes electronic eigenstates with eigenenergies $\epsilon_{\lambda_{i}}$, which are supposed to be $\sigma$-independent ($\sigma$ is a spin projection). The level spacing $\delta_{i}$ for electrons at the Fermi level in grain $i$ is small: $\delta_{i}\ll E_{i}^{c}$.

The tunneling Hamiltonian

[TABLE]

describes the hops of electrons between neighbouring grains $\langle ii^{\prime}\rangle$. We are interested in the case, when the typical dimensionless intergrain conductances

[TABLE]

are small, and the tunneling hamiltonian may be treated perturbatively; the necessary order of the perturbation theory, however, appears to be high – the lower the temperature, the higher the order!

The principal idea of any VRH-type calculation is the famous observation of Mott mott0 : at low temperatures hopping electrons prefer to visit only resonant sites (in our case for sites stand the grains), where their energies are confined to a narrow strip of width $\varepsilon$ near the Fermi energy. Decreasing $\varepsilon$ loosens the factor $\exp\{-\varepsilon/T\}$ that suppresses the conductivity due to small number of available excitations. On the other hand the resonant sites are rare (the typical distance $r(\varepsilon)$ between them grows with decreasing $\varepsilon$), so that the overlap of corresponding wave-functions $I$ is small: $I\propto\exp\{-2r(\varepsilon)/\xi\}$ and this small factor becomes still smaller with decreasing $\varepsilon$. Thus, one has to find a compromise between the two exponentially small factors, that results in certain optimal $\varepsilon_{\rm opt}(T)$ and the conductivity $\sigma\propto\exp\{-\varepsilon_{\rm opt}(T)/T\}$. In the presence of polaronic effect the above calculation scenario is somewhat modified, but the main idea remains the same.

While in case of single electron tunneling evaluation of the overlap exponential factor $I$ is straightforward, and $\xi$ is simply related to the decrement of electronic wave-function, in the case of metallic grains the origin of the exponential dependence and explicit form of $\xi$ is more sophisticated. In this section we will explore this problem, incorporating the additional physics that arises from the effects of the stray-charges flexibility.

IV.3 Amplitudes of multiple cotunneling: perturbation theory

Thus, we are interested in the amplitude of electronic transition between two distant resonant grains $l$ (left) and $r$ (right). Unless there is a tunnelling term in the Hamiltonian, the occupation numbers $\vec{n}$ are preserved. The tunneling term $\hat{\cal H}_{\rm tun}$ allows hops of electrons between neighbouring grains; to accomplish a hop between distant grains, $\hat{\cal H}_{\rm tun}$ has to be applied in some $({\cal N}+1)$-st order of perturbation theory, where ${\cal N}$ is a number of intermediate grains, constituting a continuous chain between $l$ and $r$. A multi-particle process, described by this perturbational approach, is generally known as cotunneling. It was introduced in Nazarov1990 and applied to the theory of transport in quantum dots NazarovBlanter ; Glazman and granular metals iosel2005 ; bel-hopping2 .

Let’s consider the transfer of an electron via virtual states on ${\cal N}$ intermediate grains, $l\equiv 0,1,\ldots,{\cal N},{\cal N}+1\equiv r$ being the path of adjacent grains, starting at the initial grain $l$ and terminating at the final grain $r$, so that $k$ is a neighbour of $k-1$ and $k+1$ (see Fig.4). In principle, one should sum over all possible paths connecting $l$ and $r$, but in the case of small tunnelling amplitudes $t$ we can expect the sum to be dominated by shortest paths – those, with minimal possible number ${\cal N}$ of intermediate grains. Moreover, in most situations only one particular path will be important. The amplitude of such multiple cotunneling event is given by a proper matrix element of

[TABLE]

calculated at given static set of deviations $\vec{q}$.

The amplitude $A_{\{h_{k},e_{k}\}}$ describes the process, at the end of which a hole with the set of quantum numbers $h_{0}$ is created in grain $l$ and an electron with the set $e_{{\cal N}+1}$ is created in grain $r$; generally speaking, each of ${\cal N}$ intermediate grains $k=1,\ldots,{\cal N}$ acquires an electron-hole pair with quantum numbers $\{e_{k},h_{k}\}$ (ineleastic cotunneling). However, it is possible to have $e_{k}=h_{k}$ for certain $k$’s, then no electron-hole pairs are created in the corresponding grains (elastic cotunneling). Let us denote the set of such $k$’s as $m_{1},m_{2},\ldots,m_{\cal M}$. If this set is empty (${\cal M}=~{}0$), one speaks about purely inelastic multiple cotunneling; if it includes all the intermediate grains (${\cal M}={\cal N}$), we deal with purely elastic multiple cotunneling.

Consider the set of indices $\left\{\alpha_{1}\beta_{1};...;\alpha_{{\cal N}+1}\beta_{{\cal N}+1}\right\}$, which describes certain sequence of individual tunnelings between pairs of neighbouring grains in the chain (in the $k$-th entry a tunnelling event occurs from the state $\alpha_{k}$ to the state $\beta_{k}$). Such a set is some permutation of “the natural sequence” $\left\{h_{0}e_{1};h_{1}e_{2};...;h_{\cal N}e_{{\cal N}+1}\right\}$, corresponding to progressive motion of electron from the left end of the chain to its right end. All the permuted sets contribute to the amplitude $A_{\{h_{k},e_{k}\}}$ alongside with the natural one. From (42) we get

[TABLE]

where we have used the shortcut notation $|\rangle$ for initial and $\langle|$ for final state. The tunneling matrix element $t_{e_{k}h_{k-1}}$ describes the transition of electron from the state $h$ in the $k-1$-th grain to the state $e$ in the $k$-th grain.

Within the short-range interaction model we only need to time-order the $\hat{\psi}$-operators acting in each single-grain subspace:

[TABLE]

As a result, we obtain (45) for inelastic grains and (46) for elastic ones:

[TABLE]

Here $n_{s}$ is the fermionic occupation number for state $s$ (not to be mixed with the occupation numbers of grains!). The two terms in (45), (46) correspond to different sequences of tunnelings: in one case the hole is created first and the electron afterwards, and in another case the order is reversed. Let us now collect all the factors containing $\tau_{k}$ for each $k\in 1...{\cal N}$ in (43) and integrate them out. Then we eventually arrive at

[TABLE]

The factors $B_{k}$ for inelastic and elastic grains respectively are

[TABLE]

where we have introduced

[TABLE]

Note that the explicit summation over all tunnelling sequences, which we have just performed, is only possible in the short-range interaction model.

According to Fermi golden rule the transition rate is

[TABLE]

where the Franck-Condon energy difference

[TABLE]

is the difference between electrostatic energies of final and initial states. We should also define the thermodynamic energy difference

[TABLE]

which is the energy distance between two local minima corresponding to initial $\vec{n}_{i}$ and final $\vec{n}_{f}$ states.

Note that it is impossible to distinguish the final states of the systems with the same sets of inelastic quantum numbers and different sets of elastic quantum numbers. For this reason we should add the amplitudes of such processed rather than probabilities. However, because of violent sign fluctuations of tunneling elements we will neglect such interference terms, which is implied in (50).

V Averaging of the transition rate

The transition rate (50) depends on both dynamic random variables $\vec{q}$ and static ones – $\vec{\gamma},\vec{\alpha}$ etc. There is an important difference between these two groups of variables: we should perform thermodynamic averaging of (50) over $\vec{q}$, while the frozen static variables do not imply averaging, so that, in principle, they remain “live” and characterise a specific surrounding of a particular chain. However, because of typically large number ${\cal N}$ of grains in the chain the static disorder is partly self-averaged and washed out. As we will see, only the dependence on few static random variables (mainly those, characterising the terminal grains $l$ and $r$) remains live.

V.1 Thermodynamic averaging

Gibbs thermodynamic averaging over $\vec{q}$ has the form

[TABLE]

where $Z_{c}$ is the partition function:

[TABLE]

In averaging over electronic states we will assume an equilibrium noncorrelated distribution, so that $\langle n_{s_{k}}\rangle=f_{F}(\epsilon_{s_{k}})$ and $\langle n_{s_{k}}n_{s_{l}}\rangle=\langle n_{s_{k}}\rangle\langle n_{s_{l}}\rangle$, with $f_{F}$ being the Fermi function.

Further simplification can be made if we note that the characteristic value of inelastic energies $\epsilon_{h_{l}}$, $\epsilon_{e_{l}}$ is controlled by the combination of Fermi-functions and delta-function in (50). We can see that $\epsilon_{\rm inel}\sim\Delta/L$, where $L={\cal N}-{\cal M}$ is the number of inelastic grains in the string. This means that typically $\epsilon_{\rm inel}\ll E_{c}$ . The same is also true for $\Delta_{k}$. Therefore, these quantities in the denominators in (48) can be neglected compared to $E_{k}^{(\pm)}$. However, this is not true for the elastic energies $\epsilon_{m}$, since they don’t enter the delta-function. These energies are of order $E_{c}$ and should be kept in the denominators. For the same reason (as $\epsilon_{h_{m}}\sim E_{c}\gg T$) we can substitute $1-\theta(\epsilon_{h_{m}})$ instead of $f_{F}(\epsilon_{h_{m}})$.

The summation over spin projections should be performed only for inelastic grains. It is straightforward and yields the factor $2^{L+1}$.

Also, we can substitute $t_{e_{k}h_{k-1}}$ by its “coarse grained” value at the Fermi level $t_{k,k-1}$ (i.e., averaged over an interval of energies, large compared to the level spacing $\delta$, but small compared to any other relevant scale). It allows for replacement of the summation over electronic states by integration $\sum_{\epsilon_{s}}\rightarrow\int d\epsilon_{s}/\delta_{k}$, which immediately allows to integrate out the elastic energies, and we arrive at:

[TABLE]

The variables $q_{k}$ for intermediate grains enter $\overline{w}_{if}$ only through the denominators in square brackets of both types (elastic and inelastic), and the result of integration over $dq_{k}$ in thermal averaging formally diverges at resonances, where either $E^{(-)}_{k}$ or $E^{(+)}_{k}$ goes to zero. This divergency, however, should be cut off at $E^{\pm}_{k}\sim|t|$, where the perturbation theory ceases to be valid, thus the integration in the vicinity of the resonance gives finite result. So, each integration over $dq_{k}$ gives three contributions: one from the vicinity of thermal equilibrium $q_{k}\approx 0$, and the other two from the vicinities of the resonances $q_{k}^{(\pm)}$, defined by conditions $E^{(\pm)}_{k}(q_{k}^{(\pm)})\approx 0$. The relative magnitudes of the resonant contributions are exponentially suppressed $\sim E_{c}/|t|\exp\{-E_{c}/T\}$ and thus are negligible in the low-temperature VRH regime. As a result, thermodynamic fluctuations of $q$ at intermediate grains are not relevant, and one can simply neglect them, putting $q_{k}\approx 0$, so that $E^{(\pm)}_{k}(q_{k})\to E^{(\pm)}_{k}(0)$ in (55).

V.2 Self-averaging of intermediate grains

To treat the products of large number of random factors, occurring in (55), we will apply the CLT self-averaging rule in the form

[TABLE]

and hence

[TABLE]

[TABLE]

where

[TABLE]

The factors $A_{1}$, $A_{2}$ depend on the specific distribution functions of random parameters $\alpha$ and $\gamma$ and their correlation. They are model-dependent numbers of order unity and we will never refer to their exact values in this paper.

V.3 Back to thermodynamic averaging

Finally, we have to perform the remaining summation over the energies $\epsilon_{e_{l}},\epsilon_{h_{l}}$ of the components of electron-hole pairs, created in the inelastic grains, and over the vector $\vec{q}$.

Let us introduce the new temperature scale

[TABLE]

and more convenient notations $\epsilon_{l}=\epsilon_{e_{l}}$, $\epsilon_{{\cal N}+l}=-\epsilon_{h_{l}}$ for all ”inelastic” intermediate grains, and $\epsilon_{2{\cal N}+1}=\epsilon_{e_{{\cal N}+1}}$, $\epsilon_{2{\cal N}+2}=-\epsilon_{h_{0}}$ for terminal grains. Then we can rewrite (55) as

[TABLE]

where

[TABLE]

The binomial coefficient $C_{\cal N}^{L}\equiv{\cal N}![({\cal N}-L)!L!]^{-1}$ appears in the formula (61) as a number of possible partitions of the string into elastic and inelastic subsets. Note that in (62) we have used the relation (41) and expressed the result in the terms of the average dimensionless conductance between adjacent grains $g\equiv(|t|/\delta)^{2}\ll 1$.

As it was shown in iosel2005 , in the absence of the polaronic effect (for $W\to 0$) the characteristic scale $\epsilon_{\rm inel}$ for the energies of electron-hole pairs, created in the acts of inelastic cotunneling, is much larger than temperature. It allows for an evaluation of the multiple integral in (62), leading to the result:

[TABLE]

A physical meaning of this result is clear: if the number of inelastic grains is $L$ and their total energy is $\Delta$, then the characteristic energy of one electron (or one hole), created in the process is $\epsilon_{\rm inel}\sim\Delta/L$. The phase volume, corresponding to the processes with $2L+2$ particles with energies $\epsilon\sim\epsilon_{\rm inel}$ and finite density of states is then proportional to $\epsilon_{\rm inel}^{2L+2}$, which (with the help of Stirling formula) explains formula (63).

We will see soon, that in the presence of the polaronic effect (namely, for $\Delta\lesssim W$) $\epsilon_{\rm inel}$ becomes comparable to temperature. Finding $I_{L}(\Delta)$ for $\epsilon_{\rm inel}\lesssim T$ is a much more difficult problem, which, however, is possible to resolve, using a trick, proposed in melnikov .

We should Fourier-transform the $\delta$-function in (62) to decouple the integrals over $\epsilon_{l}$:

[TABLE]

The trick is to multiply it by $1=\exp(-\frac{\Delta}{2T}-\frac{1}{2T}\sum\epsilon_{l})$, which will make the integrals convergent, and they can be easily calculated via residue theory:

[TABLE]

[TABLE]

As a result,

[TABLE]

The integral (66) can be evaluated exactly, but it is more convenient first to perform summation over $L$, which in this representation turns out to be trivial:

[TABLE]

We are left with the thermodynamic averaging $\langle\exp\{-\Delta(1/2T+it)\}\rangle_{\vec{q}}$, which is reduced to the gaussian integration that can be easily performed. As a result

[TABLE]

where

[TABLE]

and we have omitted all preexponential factors.

VI Case studies: Transition rates at different strength of polaronic effect

The general formula (68), in principle, contains the answers for all possible questions concerning different modes of the charge transfer between two particular grains. However, for understanding of the physical origin of each particular mode and the corresponding $T$-dependences, it is necessary to consider the limiting cases. In this Section we undertake such a case study.

To evaluate the integral over $t$ one can use the steepest descent method. The saddle point is located below the lowest pole on the imaginary axis of the complex $t$ plane at $t=i(1-\xi)/2T$, where $\xi$ satisfies the equation

[TABLE]

and

[TABLE]

determines the typical number $\overline{L}$ of inelastic grains in the string. In terms of $\xi$ the transition rate can be written as

[TABLE]

We have omitted in (72) the preexponential factor $\sin^{-2}\pi\xi$, because it becomes essential only at extremely low temperatures in the case of purely elastic cotunneling. This factor may be important for systems where the length of strings ${\cal N}$ is restricted from above (as in small arrays of quantum dots, or single-electron transistors), which is not the case as long the VRH conductivity of a large sample is concerned.

We will see that the entire range of the parameter $D$ ($-4W<D<\infty$) may be split into three intervals with different types of approximations applicable:

Weak polaronic effect ($D>0$, $D\gg\Delta D$) 2. 2.

Strong polaronic effect ($D<0$, $|D|\gg\Delta D$) 3. 3.

Narrow transition range ($|D|\lesssim\Delta D$)

The transition range width

[TABLE]

The latter inequality is valid because, as we will see in the Section VIII, $\max\{W,|\tilde{\Delta}|\}/T\gg{\cal N}$ for typical strings of grains that contribute to the conductivity. Below we discuss these three intervals separately.

VI.1 Weak polaronic effect: “Electron Hopping”

For $D>0$ the parameter $\xi\ll 1$, so that all trigonometric functions in (70) and (71) can be expanded. Besides that, as long as $D\gg\Delta D$ one can neglect the second term on the left hand side of (70), so that the latter is reduced to the cubic equation

[TABLE]

and one can easily express $\xi$ via $y$:

[TABLE]

The transition rate then can be rewritten as

[TABLE]

Substituting (75) to (74), we arrive at the equation

[TABLE]

which implicitly determines the function $y(z)$. Finally, for the transition rate we obtain

[TABLE]

the function $\Theta_{1}(z)$ being defined as

[TABLE]

The asymptotics of $\Theta_{1}(z)$ are

[TABLE]

Note, that $\Theta_{1}(z)$ coincides with $\varphi(z)$ which was described in iosel2005 , so that for $D>0$ their result coincides with the result of the present paper up to a slight modification $|\tilde{\Delta}|\to D$ in the definition of the argument $z$.

One can easily check, that the parameter $\xi$ is indeed small under conditions $D>0$, $\max\{W,|\tilde{\Delta}|\}\gg\overline{L}T$, no matter whether $T<T_{*}$ or $T>T_{*}$.

VI.2 Strong polaronic effect: “Polaron Hopping”

For negative and not very small $D$ we can neglect the right hand side of equation (70) (it can definitely be done if $|D|\gg\delta D$), and then

[TABLE]

where $\xi$ is not necessarily small. Substituting (81) into (72), we get

[TABLE]

and also

[TABLE]

In order to write the transition rate in a compact form, like (78), we have introduced the new function $\Theta_{2}(z)$

[TABLE]

VI.3 Narrow transition range

Looking at the results of the two preceding subsections, we conclude that $\xi\ll 1$ for all positive $D$ and also for small negative $D$, (such that $|D|\ll 4W$). However, for very small $|D|$ (both positive and negative) the first term on the left hand side of (70) can not be neglected. The condition of “very small” $|D|$ reads $|D|\sim\xi W$, where $\xi$ is given by the solution of (74). It yields $|D|\sim\Delta D$, where $\Delta D$ is given by (73). To find $\xi$ in the narrow transition range $|D|\lesssim\Delta D$ one would have to solve a quartic equation

[TABLE]

To justify the numerical coefficients in (73), we note that the saddle point equation (70) for small $\xi$ can be written in the form

[TABLE]

where $\overline{L}\equiv y(\xi){\cal N}$ itself depends on $\xi$. A formal “solution” of this equation is

[TABLE]

From the result (87) immediately follows the estimate for the width of transition range $\Delta D\sim 8\sqrt{W\overline{L}T}$. Since the requirement $\overline{L}\leq\cal N$ is always satisfied, we arrive at the inequality in (73).

We will not discuss the transition region in detail, since the corresponding range of random parameters $D,W$ is narrow and does not give any considerable contribution to physical observables.

VII Physical interpretation

Activation exponential factors in formulas (78) and (82) coincide with the corresponding factors from standard polaron hopping theory. Additional modifications, specific for our problem, arise only in the power-law factors (effective overlap integrals) due to the many-particle nature of the cotunnelling process.

VII.1 Main exponential dependence

Activation exponential factors for conventional polarons were obtained previously by other researchers. Still, we would like to present some physical arguments that qualitatively explain the origin of these factors.

Let’s firstly consider the exothermic transition with $\Delta<0$, i.e. $E_{\rm ch}(\vec{n}_{f},\vec{q})<E_{\rm ch}({\vec{n}_{i},\vec{q}})$, (see Fig.5 a). Since the total electrostatic energy decreases, the electron-hole pairs will be created in the inelastic intermediate grains to ensure energy conservation. Thus, in this situation with exponential accuracy the probability of transition is just the probability to find the system in state $\vec{n}_{i}$. Its maximum value is

[TABLE]

Now let’s consider the endothermic transition with $\Delta>0$, i.e. $E_{\rm ch}(\vec{n}_{f},\vec{q})>E_{\rm ch}({\vec{n}_{i},\vec{q}})$, (see Fig.5 b). Since the electrostatic energy increases, the electron-hole pairs will be annihilated in some intermediate grains to make up the shortfall. Such pairs are difficult to find at low temperature, which is accounted for by the additional exponential factor $\exp(-\Delta/T)$. The transition rate also contains $\exp\left\{-E_{\rm ch}(\vec{n}_{i},\vec{q})/T\right\}$ – the probability to find the system in state $\vec{n}_{i}$. Together these two contributions yield $\overline{w}\propto\exp\left\{-E_{\rm ch}(\vec{n}_{f},\vec{q})/T\right\}$, and we should take its maximum value:

[TABLE]

But what if $\vec{q}$, that delivers the minimum of $E_{\rm ch}(\vec{n}_{f},\vec{q})$, violates the requirement $\Delta>0$? The result (89) will be incorrect in this case, since there will be no exponential factor $e^{|\Delta|/T}$ for $\Delta<0$! In this case the true minimum of the activation energy should lay on the boundary of the regions with $\Delta>0$ and $\Delta<0$, i.e. at $\Delta=0$, (see Fig.5 c).

[TABLE]

Finding the minimum, we obtain

[TABLE]

where ${\cal E}_{\rm act}$ is given by (30).

To determine which of the above described solutions gives the largest contribution to the probability of transition, we should simply compare the corresponding exponents. It is easily seen that the contribution of the boundary minimum (91) dominates when $|\tilde{E}_{i}-\tilde{E}_{f}|<4W$ and explains the main exponential dependence in Polaron Hopping regime (82). In the opposite case $|\tilde{E}_{i}-\tilde{E}_{f}|>4W$ we have to choose the maximum of (88) and (89), which explains the result $\overline{w}\propto\exp\left\{-\varepsilon_{if}/T\right\}$ from (78).

VII.2 Energies of electron-hole pairs

In the Electron Hopping regime the Franck-Condon transition energy $\Delta$ (that coincides with the aggregate energy of electron-hole pairs involved) doesn’t differ much from its “relaxed” value $\tilde{\Delta}$. Hence

[TABLE]

where $\cal L$ is the large logarithmic factor (see later).

On the contrary, according to the above described physical picture, in Polaron Hopping regime the transition occurs at $\Delta=0$, see (90). However, presented in Section VI more careful calculation (which takes into account the $\Delta$-dependence of the power-law factor) shows that $\Delta$ is indeed small, but finite, in contrast with conventional polaron transitions (where it is exactly zero). In fact, $\Delta\sim\overline{L}T$, which means that in the Polaron Hopping regime

[TABLE]

Another way to see this is to note that in the presence of polaronic effect the integral (68) converges at $t\lesssim 1/T$, which means that the integrals over the inelastic energies (64) converge solely due to the denominators, which means at $\varepsilon_{\rm inel}\sim T$.

Thus, we conclude that during the crossover from electron to polaron hopping the characteristic energies of electron-hole pairs decrease from $\varepsilon_{\rm inel}\sim{\cal L}T$ to $\varepsilon_{\rm inel}\sim T\ll{\cal L}T$.

VIII Conductivity: general consideration

We now finally turn to the calculation of conductivity of a sample of granular metal. According to the general philosophy of the hopping, we introduce the Miller-Abrahams network of conductances $g_{ij}\propto\overline{w}_{ij}$, connecting each pair of grains. As usual, at low temperatures the conductivity is dominated by distant pairs of resonant grains, so that we can use the results of the preceding Section, and represent the conductances in (almost) standard form:

[TABLE]

where ${\bf r}_{i}$, ${\bf r}_{j}$ are the positions of the centres of the grains and the “effective decay length” $a_{ij}$ is given by

[TABLE]

where $a_{0}$ is the average distance between neighbouring grains, so that $|{\bf r}_{i}-{\bf r}_{j}|\equiv a_{0}{\cal N}_{ij}$. The functions $\Theta_{1}$ and $\Theta_{2}$ are given by (80) and (84).

The dependence of $a_{ij}$ on both $r_{ij}$ and characteristics of grains seems to be unusual, however it is a distinct feature of hopping in granular materials. This dependence is only logarithmic and therefore it can be taken into account perturbatively (see later).

The activation energy is standard for a polaronic problem:

[TABLE]

Thus, we have come to a system of grains, each of them being characterised by the position ${\bf r}_{i}$ of it’s centre, by the energy $\varepsilon_{i}$, and by the barrier $W_{i}$. Obviously, at low $T$ hopping electrons prefer to “make stops” only at resonant grains – those with small $\varepsilon$ and $W$, while the nonresonant grains with typically large $\varepsilon$ and $W$ may serve only as intermediate places, where electrons occur only virtually, staying there for very short times, governed by the quantum uncertainty, see Fig. 6. The random variables $\varepsilon_{i}$ and $W_{i}$ have the distribution function $\tilde{\nu}(\varepsilon,W)$ and may be correlated or uncorrelated, depending on the underlying physics. In this paper we will consider only the noncorrelated case, in which for small enough $\varepsilon$

[TABLE]

where $P(W)$ is the barrier distribution function and $\nu_{0}$ is the density of states (DOGS) near $\varepsilon=0$

[TABLE]

where $\tilde{P}(\gamma=\pm 1/2)$ is the probability density to have $\gamma=1/2$ (the same as $\gamma=-1/2$), and $n_{g}$ is the concentration of grains.

VIII.1 Statistics of barriers

Within the paradigm of short range interaction it is natural to assume that the interaction between grains and impurities is also a short range one. It means that only the impurities, situated immediately at the surface of particular grain $i$ contribute to both $\varepsilon_{i}$ (through the offset charge $\gamma_{i}$) and $W_{i}$ (through finite elasticity $\alpha_{i}$).

In general, the number of such impurities $p_{i}$ fluctuates from grain to grain. Different impurities give to $\varepsilon_{i}$ contributions of different signs, depending on their charges and positions. Therefore $\varepsilon_{i}$ would fluctuate from grain to grain already in the hypothetical case, where all $p_{i}$ are the same.

Because of the random signs of different contributions to $\varepsilon_{i}$ the density of states $\nu_{0}(\varepsilon)$ can easily be nonzero at $\varepsilon=0$. However, it is not the case for the barriers $W_{i}$. Each barrier is additive with respect to different impurities

[TABLE]

Since the individual contributions are positively defined, the random charges of different impurities do not matter, and the distribution function $P(W)$ should vanish at $W=0$. At the same time, for any positive barrier height $P(W)$ should be nonzero: there is no physical reason to expect a hard gap in $P(W)$.

VIII.1.1 Gaussian distribution of barriers

The contributions of individual impurities $W_{k_{i}}$ are assumed to be noncorrelated random variables with identical relatively narrow distributions. Thus, for the most interesting case, in which the average number of impurities coupled to each grain is large $\overline{p}_{i}\gg 1$, the central part of the distribution $P(W)$ should be a Gaussian one:

[TABLE]

where $\overline{W}$ is the average barrier, $\delta W$ is its standard deviation, $\delta W^{2}\equiv\overline{(W-\overline{W})^{2}}$. We will stick to the case $\overline{p}\gg 1$, when $\overline{W}\propto\overline{p}$ and $\delta W\propto\overline{p}^{1/2}$, so that the distribution is narrow: $\delta W\ll\overline{W}$.

The condition $\overline{p}_{i}\gg 1$ also allows to neglect the correlations between $\varepsilon_{i}$ and $W_{i}$. Note that for $p_{i}\sim 1$ the correlations are strong: for instance, the absence of impurities, associated with given grain (i.e., $p_{i}=0$) leads to $W_{i}=0$ and, simultaneously, $\gamma_{i}=0$.

Thus, in the present paper we are going to analyze the case where the joint density of states (97) has the polaronic factor $P(W)$ defined in (100). We have physically justified such choice, although other distribution functions can also be considered.

VIII.1.2 Rectangular distribution of barriers

For example, in the paper foygel a factorized rectangular distribution

[TABLE]

was studied by the method similar to the original percolational approach to the Mott VRH, proposed in halperin . The most spectacular result was obtained for the case $W_{\min}\ll\overline{\varepsilon}(T)\ll W_{\max}$, when for relevant $\varepsilon\sim\overline{\varepsilon}(T)$ within the Mott strip of width $\overline{\varepsilon}(T)$ the“generalized density of states” $\tilde{\nu}(\varepsilon,W)=\tilde{\nu}_{0}$ may be treated as $\varepsilon$- and $W$-independent. Under this condition the temperature dependence of conductivity is given by

[TABLE]

For low temperatures, such that $W_{\min}>\overline{\varepsilon}(T)$ the low-energy gap in the distribution (101) becomes essential and the conductivity acquires a hard gap as well:

[TABLE]

As we have already noted, in reality the gap in the distribution $P(W)$ is not hard, it is most likely to have an exponential tail at zero. In the rest of this paper we will adopt the gaussian distribution of barriers (100).

VIII.2 Percolation problem

The conductivity of Miller-Abrahams network with exponentialy wide distribution of conductances should be found by means of percolation theory halperin ; ES_book .

We have come to a kind of coloured percolation problem on random sites, homogeneously distributed in $({\bf r},\epsilon)$ space with the density (DOS) $\nu_{0}$. Each site $i$, besides its position $({\bf r}_{i},\epsilon_{i})$, is characterised by a positive variable $W_{i}$, distributed according to (100). The variables $\epsilon_{i}$ and $W_{i}$ are interpreted as two components of a composite colour of a site $i$, and $n(\epsilon,W)\equiv\nu(\epsilon)P(W)$ is the density of sites with given colour.

By definition, a pair of sites $\langle ij\rangle$ is “$\xi$-connected” if

[TABLE]

where

[TABLE]

As it is usual for the problems of VRH type, the dc charge transfer processes at low temperatures are confined to certain critical subnetwork of resonant grains with small energies $\varepsilon$ within certain narrow strip of width $\overline{\varepsilon}(T)$. This width is $T$-dependent and has to be defined self-consistently. It is more convenient to proceed with calculations in different temperature ranges separately.

IX Conductivity in different temperature ranges

There are three principal energy scales in our problem: (i) The average barrier $\overline{W}$, (ii) The barrier dispersion $\delta W\ll\overline{W}$, and (iii) the width of the Mott strip $\overline{\varepsilon}(T)$. Correspondingly, there are three ranges of temperature with different dominating physics.

IX.1 Electron Hopping – standard Mott VRH

In this case $\overline{W}\ll\overline{\varepsilon}(T)$, so that the polaron effect is negligible, and the standard Mott law is valid

[TABLE]

The width of effective energy strip is

[TABLE]

The condition $\overline{W}\ll\overline{\varepsilon}(T)$ is equivalent to

[TABLE]

being the temperature of crossover between the electron and the polaron regimes.

IX.2 Polaron Hopping, grains with typical barriers

Here

[TABLE]

so that the polaron effect is dominant, but the fluctuations of the barriers are still negligible. In this range the results of the model with identical barriers for all grains are applicable:

[TABLE]

Note that the second term in the exponent is small compared to the first one. The condition $\delta W\ll\overline{\varepsilon}(T)\ll\overline{W}$ is equivalent to $T_{c2}\ll T\ll T_{c1}$, where

[TABLE]

In principle, one can also write a general formula that describes the behaviour of the conductivity in the entire range $T\gg T_{c2}$:

[TABLE]

$F_{1}(x)$ being a universal function, with a shape, depending only on the space dimensionality $d$. Its asymptotics are

[TABLE]

where $c_{1},c_{2}$ are some universal constants. An interpolation formula for the function $F_{1}(x)$ was proposed in foygel .

IX.3 Polaron Hopping, grains with the barriers in the Gaussian tail of distribution

Here

[TABLE]

and the grains with anomalously low barriers $W_{i}$ from the tail of the distribution (100) dominate the critical network.

[TABLE]

Namely, the critical subnetwork consists of grains with $W_{i}$ from a narrow strip

[TABLE]

The derivation of results (118) and (119) is given in the Appendix A. It is based on the possibility to reduce the initial multiparametrical percolation problem to certain universal one, governed by a single parameter.

Note, that the result (118) is only valid, if the second term in the exponent is small compared to the first one. In the terms of temperature it means

[TABLE]

In principle, one can also combine the results in the entire range $T_{\min}\ll T\ll T_{c1}$ in one formula:

[TABLE]

where $F_{2}(x)$ is, again, a universal function, depending only on the space dimensionality $d$. Its asymptotics are

[TABLE]

where $c_{3}$ is a constant. Note that the high-temperature asymptotics of (121) coincides with the low-temperature asymptotics of (115).

The inequality $T\gg T_{\min}$ was imposed to secure the condition of relatively small relevant fluctuations with $W-\overline{W}\ll\overline{W}$. This condition is needed to justify their universal gaussian distribution. At $T\sim T_{\min}$ the relevant fluctuations are large: $W-\overline{W}\sim\overline{W}$ and their distribution is strongly model-dependent.

The range of exponentially low temperatures $T\lesssim T_{\min}$ may only be of academic interest, because at such low $T$ the the barrier would rather be penetrated due to quantum tunneling mechanism, than due to the activational one, which we are discussing here.

IX.4 Effective radius of the wavefunctions

Now we have to specify the “effective radius of the wave function” $a$, entering the criterion (105). In the Electron Hopping regime, when $|\varepsilon|\gg W$, the actual “radius” $a=a_{ij}$ depends on the characteristics of particular grains $i$ and $j$. This logarithmic dependence can be easily taken into account. This was done in iosel2005 with the help of the perturbational method in the percolation theory. As a result, the effective radius $a$ appeared to be inversely proportional to the large logarithm ${\cal L}(T)$:

[TABLE]

where we define $a_{0}$ as the mean distance between the grains: $r_{ij}=a_{0}{\cal N}_{ij}$. As a consequence, the parameter $T_{M}$ in (109) acquires a logarithmic temperature dependence

[TABLE]

Fortunately, for $|\varepsilon|\ll W$, in all the Polaron Hopping regimes, including both weak fluctuations and strong fluctuations cases, $a$ doesn’t depend on $i$ and $j$, but only on temperature. This can be seen from the expression for transition rate, which contains the factor

[TABLE]

Thus, no perturbational method is needed in this case and

[TABLE]

so that the parameter $T_{M}^{\prime}$ in (113) also becomes logarithmically $T$-dependent:

[TABLE]

[TABLE]

X Summary

We have considered a granular metal with charged impurities (stray charges) that are not rigidly fixed at certain positions, but can be slightly displaced from their equilibrium positions due to electrostatic forces and thermal fluctuations. The flexibility of the system of impurities leads to a “polaronic effect”, that can be either weak (for $T\gg T_{c1}$) or strong (for $T\ll T_{c1}$), where the crossover temperature

[TABLE]

and $T_{M}\approx T_{M}(T_{c1})$ is given by (135) below.

In this paper we have studied a model with short range interaction (e.g., screened Coulomb) and taken into account only the thermoactivational mechanism of the penetration of the polaronic barrier by the configurational degrees of freedom, not the tunneling one. The latter may become relevant at very low temperatures (the larger the effective mass $M$ of the impurities, the lower the corresponding crossover temperature).

Experimentalists often introduce the differential activation energy according to

[TABLE]

Below we discuss the temperature dependence of $E_{\rm act}(T)$ in different temperature ranges.

X.1 Weak polaronic effect

In the range $T\gg T_{c1}$ the hopping can be with minor modifications described by the electronic multiple cotunneling iosel2005

[TABLE]

X.2 Strong polaronic effect

The range $T\ll T_{c1}$, where the cotunneling is dominated by polaronic effect, is split into two subranges with the crossover temperature

[TABLE]

• •

For $T_{c2}\ll T\ll T_{c1}$ the spatial fluctuations of polaronic barriers $W_{i}$ can be neglected and the choice of resonant grains that constitute the effectively conducting network is dictated exclusively by the values of $\varepsilon_{i}$: they should lie in the Mott-like strip $|\varepsilon_{i}|\lesssim\overline{\varepsilon}(T)$, while the particular values of $W_{i}$ are irrelevant. Here the approach, similar to the one, proposed in foygel , is applicable, and

[TABLE]

where the second term is relatively small, compared to the first one.

• •

For $T\ll T_{c2}$ both $\varepsilon_{i}$ and $W_{i}$ are important, the effective network is formed by the grains for which both $\varepsilon_{i}$ and $W_{i}$ are anomalously small: $|\varepsilon_{i}|\lesssim\overline{\varepsilon}(T)$ and simultaneously $|W_{i}-\frac{1}{2}W_{\rm eff}(T)|\lesssim\overline{\varepsilon}(T)$. Here

[TABLE]

In this temperature range

[TABLE]

Thus, due to the tail in the distribution of the barrier fluctuations, the effective activation energy continues to decrease even at lowest temperatures, though this decrease becomes very slow.

X.3 Elastic vs inelastic cotunneling

The above mentioned crossovers discriminate different modes of hopping with respect to the strength of polaronic effect (crossover at $T\sim T_{c1}$) and to the importance of the fluctuations of the latter (crossover at $T\sim T_{c2}$). There is, however, one additional crossover at $T\sim T_{c0}$ that discriminates different modes with respect to the character of cotunneling: elasic or inelastic. For $T\gg T_{c0}$ the “constants” $T_{M}(T)$ and $T^{\prime}_{M}(T)$ logarithmically increase with the lowering of temperature, while at $T\lesssim T_{c0}$ they saturate and become $T$-independent:

[TABLE]

[TABLE]

The explicit form of $T_{c0}$ depends on the relation between $T_{c0}$ and $T_{c1}$:

[TABLE]

For $T_{c0}\gg T_{c1}$ the crossover between elastic and inelastic cotunneling takes place within the weak polaron effect domain, while for $T_{c0}\ll T_{c1}$ it happens within the strong polaron effect domain.

XI Conclusion

In this paper we have introduced a concept of polaronic effect in granular systems, related to the flexibility of the random charges, trapped in the insulating matrix. We have explained, how this effect is manifested in the conductivity of the system, the latter being controlled by multiple cotunneling of electrons through long “strings” of adjacent grains. The basic line of our reasoning was similar to that of every VRH-like calculation, and could be split into two basic steps:

(i) We calculate the carrier transition rate between distant resonant sites, taking into account all the necessary physics (the cotunneling and coupling to flexible impurities). The transition rate is given by the expression (68) and its simplifications (78) and (82).

(ii) Afterwards, we find the conductivity of Miller-Abrahams network of conductances with the help of percolation theory. We distinguish two important temperature ranges, namely $T>T_{c1}$ (Electron Hopping) and $T<T_{c1}$ (Polaron Hopping), which differ in both physical characteristics of transport and $T$-dependence of conductivity.

Besides the crossover between Electron and Polaron hopping we have also studied another crossover – between elastic and inelastic cotunneling regimes. It takes place at the temperature $T=T_{c0}$, given by (137). In the presence of strong polaronic effect this temperature turns out to be the same as for single quantum dot, in contrast with the Electron Hopping regime, where it is much lower.

A few important questions remain open and are subject to future research.

(i) In the present work we have treated the configurational degrees of freedom (i.e., coordinates of the charged impurities) as classic ones. This can be justified only if the temperature $T$ is higher than characteristic frequencies $\omega_{i}$ of the impurities vibrations; at lower $T$ the polaronic barriers would be penetrated by means of quantum tunneling. In principle, it should lead to the reentrance of Mott law at $T<\overline{\omega}$. However, since $\omega_{i}$ may be different for different impurities, the crossover to tunneling may occur not simultaneously at all grains, and that may give rise to some interesting new physics.

(ii) We have studied here only the case of short range interactions. Some qualitative ideas about the effects of the long range Coulomb interactions for the Electron Hopping were presented in iosel2005 .

(iii) The lack of experimental data on polaron effect in granular systems does not allow to reliably choose the distribution function of the polaronic barriers $W$. It is not clear whether it should be wide or narrow, have a power-law tail at $W\to 0$ or exponential one.

(iv) It is also unclear if the present model, in which the occupation numbers $n$ are coupled to oscillators, is adequate for real world applications. Another possible option would be the coupling to the tunneling two-level systems (TTLS, see TTLS1 ; TTLS2 ). Physically the two-level system may be represented by impurity atoms that can tunnel between two adjacent potential wells. Such systems play important role in physics of glasses, they also contribute to dephasing in qubits TLS .

We believe that the formalism developed in this paper will allow us to answer at least the first question from the list above. As for the other questions, some novel approaches may be required.

We are indebted to M.V.Feigel’man and L.B.Ioffe for helpful discussions. This work was supported by 5top100 grant of Russian Ministry of Education and Science.

Appendix A Percolation with gaussian distribution of barriers. Derivation of the formula (118)

We are interested in the low-$T$ behaviour of the system, when the critical subnetwork, responsible for the conductivity, is comprised of the rare grains with anomalously small $\epsilon_{i}$ and anomalously low barriers $W_{i}$ from the tail of the distribution (100). As we will see later, these parameters are confined in narrow strips

[TABLE]

Now we will adopt the inequalities (138) as a conjecture, the real values of $\overline{\varepsilon}_{\epsilon},\overline{\varepsilon}_{W}$ will be determined later, and these values will justify (138).

Based on the assumption (138), the expression (108) for $\Lambda_{ij}$ can be expanded and we get

[TABLE]

The result (139) means that indeed each site may be characterised by one “composite colour” $\tilde{W}_{i}\equiv W_{i}+|\epsilon_{i}|/2>0$ not by two independent colours $\epsilon_{i}$ and $W_{i}$. The density of this composite colour is

[TABLE]

(the factor 2 arises due to symmetry $\varepsilon\to-\varepsilon$), and

[TABLE]

Now we introduce dimensionless variables

[TABLE]

and arrive at the dimensionless percolation problem

[TABLE]

with the density in the $({\bf x},\delta)$-space

[TABLE]

where

[TABLE]

In particular, we will need the asymptotics

[TABLE]

It is convenient to write

[TABLE]

and note that $w\equiv 2\overline{W}-W_{\rm eff}(T)\ll\overline{W}$ and $\xi T\delta\ll\overline{W}$, so that

[TABLE]

As we will see soon, for typical $\delta$ and $w$ a hierarchy $\overline{W}\delta\ll\delta W\ll w$ holds. Therefore only the first two terms in the exponent should be kept, while last term is much less than unity and can be neglected. As a result, we can rewrite (148) in a form

[TABLE]

After the renormalisation of variables

[TABLE]

we arrive at the universal percolation problem, where the sites of colour $u$ are randomly distributed with the density

[TABLE]

while the percolation criterion reads

[TABLE]

Note that this problem is characterised by single constant

[TABLE]

and therefore the percolation should be established at

[TABLE]

where $A_{d}\sim 1$ is some universal constant, depending only on the space dimensionality $d$.

So, the dependence of the effective barrier $W_{\rm eff}(T)\equiv 2\overline{W}-w$ on the parameters may be found from the equation for $w$:

[TABLE]

In the leading logarithmic approximation the solution of (155) reads

[TABLE]

Now we can estimate $\overline{\varepsilon}_{\varepsilon}$ and $\overline{\varepsilon}_{W}$. As it follows from (150)

[TABLE]

which justifies our conjecture (138).