Current saturation and Coulomb interactions in organic single-crystal transistors

Electronic transport through rubrene single-crystal field effect transistors (FETs) is investigated experimentally in the high carrier density regime (n ~ 0.1 carrier/molecule). In this regime, we find that the current does not increase linearly with the density of charge carriers, and tends to saturate. At the same time, the activation energy for transport unexpectedly increases with increasing n. We perform a theoretical analysis in terms of a well-defined microscopic model for interacting Frohlich polarons, that quantitatively accounts for our experimental observations. This work is particularly significant for our understanding of electronic transport through organic FETs.


Introduction
The use of single-crystalline material for the fabrication of organic field-effect transistors (FETs) has given, over the past few years, the experimental control needed for the investigation of the intrinsic transport properties of dielectric/organic interfaces [1]. This has resulted in the observation of anisotropic transport [2], a metallic-like temperature dependence of the mobility [3], the Hall effect [4] and quasiparticle response in the infrared conductivity [5]. Following this progress, the successful quantitative analysis of experiments in terms of a simple microscopic model has recently been possible in single-crystal FETs with highly polarizable gate dielectrics [6]. In these devices, charge carriers accumulated electrostatically at the surface of the organic semiconductor were shown to interact strongly with the polarization in the dielectric, leading to the formation of Fröhlich polarons [7]. The quantitative agreement between the behavior predicted by the Fröhlich Hamiltonian and the experimental data demonstrates that the microscopic processes responsible for transport in organic FETs can be qualitatively different from that in bulk organic semiconductors.
Here, we investigate the transport properties of rubrene (C 42 H 28 ) single-crystal FETs with highly polarizable Ta 2 O 5 gate dielectrics, in the as yet unexplored high-carrier-density regime. We find that in this regime the electrical characteristics of the devices exhibit pronounced deviations from those of conventional FETs. Specifically, in all devices the source-drain current I sd stops increasing linearly with the gate voltage V g , and shows a clear saturation. Concomitantly, the activation energy for temperature-dependent transport increases when V g increases, a trend opposite to that usually observed in organic transistors. We build on the Fröhlich model, which was used in our previous study to describe the dressing of the carriers by the polarizability of the gate dielectric [6], and find that the observed behavior can be quantitatively explained by considering the effects of the Coulomb interactions between holes. Our results confirm that organic single-crystal transistors are suitable for the experimental investigation of the intrinsic electronic properties of dielectric/organic interfaces, and extend our fundamental understanding of transport in organic transistors to the high-carrier-density regime.

Experimental results
The rubrene transistors used in our experiments are fabricated by laminating thin (<1 µm thick) rubrene single crystals onto a substrate with pre-fabricated FET circuitry, as described in [8,9]. The gate dielectric is a ≃500 nm-thick layer of Ta 2 O 5 (dielectric constant ǫ s = 25; breakdown voltage ≃6.5 MV cm −1 ), enabling the accumulation of up to 10 14 holes cm −2 in the FET channel. During device operation, however, we never exceed gate voltage values corresponding to half the breakdown field, to minimize the chance of device failure [10] (we reach hole densities of 5 × 10 13 cm −2 , i.e. one order of magnitude higher than in [6]). The transistor electrical characteristics were measured in the vacuum (10 −6 mbar) chamber of a flow cryostat, between 300 and 125-210 K depending on the specific device (at low temperature the devices can easily break due to the difference in thermal expansion between crystal and substrate). Figure 1 shows the source-drain current I sd measured on two different devices as a function of gate voltage V g , at different temperatures. Since the carrier density n in the channel is linearly proportional to the gate voltage V g (n = C(V g − V th )/e, where C is the capacitance per unit area and V th is the threshold voltage), it is normally expected that I sd also increases linearly with V g . The data, however, show a pronounced deviation from the conventional linear behavior, since at high V g the source-drain current tends to saturate. The effect is reversible and reproducible: the devices can be measured many times without noticeable difference in the data, which implies that the saturation of I sd is not due to device degradation [10]. We have measured tens of similar devices (room-temperature mobility values between 1.0 and 1.5 cm 2 V s −1 ) and found a similarly pronounced saturation of the source-drain current at high gate voltage in virtually all cases, albeit with some differences in the details of the I sd -V g characteristics (illustrated by the data in figures 1(a) and (b)).
The temperature dependence of the current I sd is also unusual. At any fixed value of carrier density n, the current I sd decreases in a thermally activated way [6]. The inset of figure 1(a) shows that the activation energy E a depends on n and that, for densities larger than ≃0.02 holes molecule −1 , E a increases with increasing n. This observation is surprising, because normally (e.g. in organic thin-film FETs) E a exhibits the opposite trend, which is attributed to the effect of disorder (filling of traps) in the organic material.
A behavior of the I sd versus V g curves such as the one just described has not been reported earlier in FETs realized on other conventional gate dielectrics (e.g. SiO 2 , Al 2 O 3 and Si 3 N 4 ). Only recently, a strongly nonlinear I sd -V g relation was observed in organic FETs with gate 4 electrolytes [11,12]. Gate electrolytes enable the accumulation of a large density of carriers, 5 × 10 13 carriers cm −2 or more, similar to our Ta 2 O 5 dielectrics. At these densities, at which the charges accumulate in the uppermost molecular layer of the crystal [6], the average distance between the carriers is only a few molecules and the resulting (bare) Coulomb interaction is a few hundreds of millielectronvolts, much larger than the thermal energy at room temperature. This, together with the poor screening associated with the thermally activated motion of charge carriers, suggests that interactions may be responsible for the observed anomalous FET characteristics.

The microscopic model
As we proceed to show, the inclusion of Coulomb interactions between holes does indeed successfully account for our experimental findings. To introduce the theoretical framework needed for the quantitative description of the device characteristics, we recall the conclusions of [6], namely that holes in rubrene FETs with highly polarizable gate dielectrics form small Fröhlich polarons due to their strong interaction with the ionic polarization at the interface. This conclusion was established by analyzing transport measurements in the low-density regime, performed on FETs with a range of different gate dielectrics. Specifically, it was shown that, using gate dielectrics of increasing polarizability, the 'metallic-like' transport previously observed in rubrene [3,13] is progressively suppressed, turning into an activated behavior characteristic of self-localized particles. The experimentally observed behavior was consistently explained in terms of the Fröhlich model-which describes the polar interaction with interface phonon modes-by including explicitly the narrow-band nature of the organic crystal [6]. In FETs with Ta 2 O 5 devices, the most polarizable dielectric used in our previous work, the interaction with the interface modes completely dominates the carrier motion, reducing the mobility at 200 K by two orders of magnitude with respect to isolated rubrene.
Here, we extend the analysis of [6] to the high-density regime that had not been analyzed previously, by adding a term which accounts for the mutual Coulomb interactions between the carriers [14]. The resulting Hamiltonian: is the simplest model that catches the main microscopic mechanisms that determine the transport in our devices: finite electronic bandwidth, polar coupling of the carriers with the gate dielectric and mutual interaction between carriers. The first three terms describe free holes in a tightbinding scheme (c + i and c i are the corresponding creation and annihilation operators on site i), interacting with lattice deformations X ℓ through a non-local coupling g iℓ . Physically, g iℓ represents the polar interaction between the holes in the two-dimensional (2D) conducting channel and the phonons at the organic/dielectric interface [6,15]. When this interaction is sufficiently strong, as is the case for the rubrene/Ta 2 O 5 interface, it leads to the formation of small polarons, i.e. self-trapped states whose characteristic radius is comparable with the lattice spacing. This is demonstrated theoretically in the appendix, and is confirmed by the experimental data, which exhibit a thermally activated mobility characteristic of small polarons (see [6] and figure 2 below). As in [6], we explicitly neglect the kinetic energy of the phonons, which is appropriate in the adiabatic regime, for which the lattice dynamics is slow compared to the band electrons. . Red squares are obtained by fitting the data with the theory for interacting polarons described in section 3, while black circles correspond to linear fits of the I sd -V g curves restricted to the low-density regime, as done in [6]. The lines are best fits to the polaronic thermally activated behavior, yielding respectively = 53 and 55 meV (see text). The inset is the polaron mobility for the sample of figure 1(b) (units are the same as in the main panel; the temperature range is too small to extract reliable values for ).
The last term-not considered previously-represents the mutual interaction between holes through the unscreened Coulomb repulsion V i j , which is long-ranged (this contrasts with the Hubbard model, where a purely on-site interaction is considered, leading to the appearance of interaction effects only at densities close to 1 carrier molecule −1 ). In the following, we evaluate quantitatively the effects of Coulomb interactions on the transport characteristics. We shall only present the main formulae needed for the comparison with the experiments, leaving the full details of the calculations to a further publication 4 .
In the temperature range of the experiment the motion of small Fröhlich polarons at the rubrene/Ta 2 O 5 interface takes place through successive uncorrelated hops from a given molecule to a nearest neighbor. In this case, the mobility can be rigorously determined by analyzing the hopping rate of a given carrier between two neighboring molecules (say i = 1, 2) [17]. Integrating out the electronic degrees of freedom, it is found that the hole dynamics is coupled to the long-range lattice polarization solely through the collective coordinate Q = 1 g ℓ (g 1ℓ − g 2ℓ )X ℓ . Taking advantage of the adiabatic assumption, we end up with the following 4 It should be noted that the electron-lattice interaction in equation (1) can mediate an effective attraction between holes, leading under specific conditions to the stabilization of bipolarons, i.e. two holes bound together by a common lattice deformation. This requires in particular that the phonon-mediated attraction is sufficiently strong to overcome the instantaneous Coulomb repulsion between holes. This condition can be cast in terms of the effective dielectric constants of the interface at high and low frequencies, and corresponds to η = (ǫ ∞ + κ)/(ǫ s + κ) being less than a given critical value η c . An accurate estimate η c = 0.131 was given in [16] for a pure 2D Fröhlich interaction in the continuum limit. Such a critical value should be further reduced in the present case, because the finite distance z of the conducting channel to the interface suppresses the electron-phonon interaction at large momenta, M q ∼ e −qz / √ q (see the appendix). For the present rubrene/Ta 2 O 5 interface, we have η = 0.26, which lies safely above the critical value η c , so that bipolaron formation can in principle be excluded. 6 double-well potential: where we have defined g 2 = ℓ g 1ℓ (g 1ℓ − g 2ℓ ). The physical meaning of the above equation is that in the adiabatic regime, the hole follows instantaneously the dynamics of the slow coordinate Q in the effective potential E ad . The two minima at Q ≃ ∓g/2K are then associated with the hole being at sites 1 and 2, respectively. The corresponding hopping rate can be evaluated by standard techniques [18], by studying the escape rate of the coordinate Q over the barrier. The variable ξ that appears in equation (2) (defined below) accounts for the instantaneous repulsion of all the other carriers in the conducting layer, and vanishes in the lowdensity limit. In that case the double-well potential equation (2) is symmetric, with a barrier given by = g 2 /4K − t (the polaronic gap, valid in the strong coupling regime g 2 /4K ≫ t), all the carriers diffuse with the same hopping rate w 0 = (ω s /2π ) e − /k B T , and we recover the mobility of independent polarons obtained in [6], namely µ P (T ) = (ea 2 /hT )w 0 (here ω s is the frequency of the interface optical phonons and a is the intermolecular distance).
In the high-density regime accessible with Ta 2 O 5 gate dielectrics, it becomes essential to include the Coulomb interactions that were not analyzed previously, which requires considering the case of a finite ξ . From the model equation (1) it can be shown that: whereV i j = 2e 2 /(ǫ s + κ)/R i j is the Coulomb potential, appropriately screened by the ionic polarization at the interface (κ = 3 is the dielectric constant of rubrene). A finite ξ causes an energy unbalance between the initial and final hopping sites (the double-well potential E ad (Q) becomes asymmetric) so that the hopping rate changes to w(ξ ) = w 0 e −ξ/2k B T . Now each carrier diffuses with a different hopping rate, determined by its own environment. Correspondingly, the mobility of the sample is obtained as the statistical average w(ξ ) over all the possible values of the local electronic potential ξ . Since the polaronic barrier sets the dominant energy scale, due to the high polarizability of the gate dielectric Ta 2 O 5 used in our devices, we always have ξ ≪ in the explored density range. As a result, correlations between subsequent hops can be neglected to a first approximation, which corresponds to replacing w(ξ ) → w( ξ ). This yields the leading term of the mobility in the presence of electron-electron interactions as The density-dependent quantity ξ represents the average extra energy cost for hopping from site to site, induced by the long-range Coulomb interactions between the carriers. As the comparison between theory and experiments will show (see below), it is this dependence of µ on the density (through ξ ) that causes the saturation of I sd observed in figure 1.
In the linear response regime, the calculation of ξ follows, through equation (3), from the statistical distribution of the occupation numbers {n j } constrained to n 1 = 1 (site 1 is initially occupied by the hole, which then hops to site 2). Such a constraint reflects the fact that the relaxation of the remaining carriers occurs on a timescale ∼w −1 0 , much longer than that of the hopping event under consideration, ∼ω −1 s , and can be neglected. In this case, the occupation numbers {n j } are given, by definition, by the pair distribution function g (2) (R 1 j ) of the electronic system at equilibrium [19]. Furthermore, due to the diffusive nature of their hopping motion (transport is thermally activated), the small polarons behave as classical interacting particles. The two-body correlations of such a classical fluid of charged particles are determined solely by the Coulomb interactions, and are therefore equivalent to those of a classical 2D plasma. Observing that in the explored density range the mean interparticle distance is always larger than the lattice periodicity, we can use for g (2) the known pair distribution function of the 2D plasma in the continuous limit. Correspondingly, we can replace the discrete sum in equation (3) by an integral, which finally yields In the above equation, the universal function F(Ŵ) is an intrinsic property of the 2D electron plasma. It depends on a single dimensionless parameter Ŵ that measures the ratio of the Coulomb interactions to the thermal energy. F(Ŵ) can be evaluated directly by using the data of g (2) (r ) available from extensive Monte Carlo simulations [20]. In the range 1 < Ŵ < 20, we find F(Ŵ) = 1 + 0.85Ŵ to within 1% accuracy. In the present FET geometry, the plasma interaction parameter is given by Ŵ = 2 √ π ne 2 /(k B T )/(ǫ s + κ), and can attain the value ∼9 at the lowest temperatures/highest V g , placing our devices in the range of weak to moderate interactions (electron crystallization is expected at Ŵ = 125).

Discussion
We are now in a position to compare the experimental data with our calculations. By inserting equation (5) into (4), and using the definition of Ŵ, we obtain the explicit functional dependence of the mobility µ on density n = C(V g -V th )/e and temperature. The I sd -V g curves can then be calculated from I sd = A|V g -V th | e − ξ /2k B T , with A = µ P C V sd W/L having the dimensions of conductance (here W and L are the width and length of the channel, V sd the applied source-drain voltage and e the electron charge). To analyze the data, we fix the temperature and fit the I sd -V g curve using µ P and V th as fitting parameters. The polaron microscopic parameters are obtained by analyzing the temperature dependence of µ P extracted in this way.
The values of µ P extracted for the different devices from the theory of interacting polarons are reported in figure 2 (red squares). They are close to the ones that one would obtain from linear fits of the I sd -V g curves in the low-density regime (the conventional definition of the mobility; data shown as black circles in figure 2) and they exhibit the same trend. Indeed, by comparing the temperature dependence of the fitted µ P to the theoretical relation µ P (T ) = (ea 2 ω s /2π k B T )exp(− /k B T ) we obtain, for the device of figure 1(a), the following values of the polaron parameters: ω s /(2π ) = 390 cm −1 and = 53 meV. These values are comparable to those obtained by assuming a linear I sd versus V g dependence (ω s /(2π ) = 315 cm −1 and = 55 meV). This shows that the results of the interpretation based on the theory with interactions is compatible with the analysis performed by only looking at the low-density part of the I sd -V g curves, and that the effect of the interactions is determined by the carrier density only (i.e. without the need to include any additional parameter). For the second device (the inset of figure 2), a precise quantitative analysis of the temperature dependence of the mobility is prevented by the restricted temperature range of the available data, and it is not possible to extract reliable values for the parameters ω s /(2π ) and . Still, the fitted values of µ P and V th are well defined for all the curves of figure 1(b), thus allowing for an accurate analysis of the density dependence of the mobility at each given temperature.  The values of V th obtained from the interacting theory also follow the same trend as those that one obtains by extrapolating the linear part of the I sd -V g curves to zero current, i.e. the usual definition of V th . The two estimates differ by a systematic offset of ∼7 V, which is approximately the same for all samples and at all temperatures. Such a deviation does not have a particular physical significance, because the exact absolute value of V th defined in the usual way does not have a precise physical meaning.
The continuous lines in figures 1(a) and (b) are the results of the theoretical fits of the I sd -V g curves, using the theory for interacting polarons. The theory reproduces the saturation of the source-drain current at high carrier density and the quantitative agreement with the data is remarkable, at all temperatures and in both the devices analyzed. From equation (4), we see that theory predicts an activation energy which increases with increasing density (via the term ξ /2). This is illustrated in the inset of figure 1(a), where the circles are the experimental values extracted from the data and the continuous line is the theoretical curve computed using equation (5), again without any new free parameters. Also here the agreement is very satisfactory, and at small n we recover the value given by the analysis of the linear regime in terms of non-interacting polarons [6].
Notably, the theoretical derivation presented above implies that the density-induced increase of the hopping activation barrier (cf equation (4)) is a scaling function of the parameter Ŵ of the interacting plasma. This scaling behavior is checked in figure 3, where we use the experimental data to directly plot the quantity log(I sd /A|V g − V th |)/|V g − V th |, which is proportional to the universal function F(Ŵ) (see equations (4) and (5)) versus |V g − V th |/T ∝ Ŵ. In terms of these variables and at sufficiently high gate voltages, the data do indeed tend to collapse on a single linear curve corresponding to F(Ŵ) = 1 + 0.85Ŵ (full line in the figure). Considering the high sensitivity to details of such a scaling plot, the agreement between theory and experiment is satisfactory. It confirms the validity of the assumptions underlying our derivation, giving a strong indication that the anomalous behavior of the electrical characteristics observed at finite carrier concentrations is indeed due to the long-range Coulomb interactions between the carriers.

Conclusion
We have studied electronic transport through rubrene single-crystal FETs in the high-density regime and found new unexpected phenomena, such as a saturation of the current I sd versus V g and an increase of the activation energy with increasing density. The experimental data are accurately reproduced by a theoretical analysis based on interacting small Fröhlich polarons, which naturally extends previous studies in the low-density regime to include the effect of the long-range Coulomb repulsion between charge carriers. Remarkably, the quantitative description of the high-density regime does not require the introduction of any new parameter, as the effect of the interactions is fully determined by the (known) density of charge carriers. Our results demonstrate that, by devoting sufficient effort to control the experimental systems, single-crystal organic FETs do allow a detailed quantitative study of the intrinsic transport properties of organic FETs.

Acknowledgment
We thank R W I de Boer and N Iosad for experimental help and the NWO and Nanoned for financial support.

Appendix. Calculation of the polaron radius
The theoretical analysis developed in the text treats the charge carriers as small Fröhlich polarons, having dimensions comparable to the lattice spacing. To demonstrate that this is the case, in this appendix we evaluate the radius of the polarons that form at a rubrene/Ta 2 O 5 interface. To this aim we solve the adiabatic model equation (1) variationally for a single hole in the HOMO band, and evaluate the polaron radius R P as a function of the electron-phonon coupling strength. The calculations indeed show that R P is of the order of the lattice spacing.
The electron-phonon interaction g i j is obtained from the Fourier transform of the electron-phonon matrix element M q = M 0 e −qz / √ q [6,15]. Here, M 2 0 = 2πhω s e 2 β/S, where S is the total surface of the device, and z is the distance of the conducting layer to the polar dielectric. The parameter β = (ǫ s − ǫ ∞ )/(ǫ s + κ)/(ǫ ∞ + κ) is a combination of the known dielectric constants of the gate dielectric and of the organic semiconductor, which determines the strength of the electron-phonon coupling. In principle, the interaction g i j defined above must be cut off at short distances to account for the discreteness of the rubrene lattice. A precise prescription to carry out this procedure will be reported elsewhere. For the present purposes it suffices to say that the overall behavior of the polaron radius does not depend on the particular choice of the short-distance cutoff.
As is customary for lattice models, we introduce the dimensionless electron-phonon coupling defined as the energy of a fully localized polaron in units of the half bandwidth D ∝ t. The energy E P can be written in general as: where the sum spans all the lattice sites R i . The model is solved by taking the following Gaussian trial wavefunction: where N is a normalization factor. The parameter α is obtained after minimization of the ground-state energy. The polaron radius is then given by which tends to 1/α in the large polaron limit (small α, R P ≫ a). The results for R P at different values of the electron-phonon coupling are summarized in table A.1. As expected, the polaron radius becomes comparable to the lattice spacing around λ ∼ 1, i.e. when the polaron energy is of the order of the half bandwidth. Performing the sum in equation (A.2) with the electron-phonon coupling parameters appropriate to our devices (β = 0.1 and z = 2 Å), we obtain E P = 170 meV. The width of the HOMO band in the rubrene crystal has been evaluated through a semi-empirical ab initio calculation in [21] as 2D = 341 meV. Such a value should be considered as an upper bound to the actual bandwidth, which was shown in [22] to be sizably reduced by the effects of molecular polarization. A further effective reduction of the bandwidth can arise from thermal fluctuations, as considered in [13]. From the above arguments and from the definition equation (A.1), we conclude that λ > 1, which, as can be seen from table A.1, corresponds to a small polaron.