Unconventional Current Scaling and Edge Effects for Charge Transport through Molecular Clusters

Metal–molecule–metal junctions are the key components of molecular electronics circuits. Gaining a microscopic understanding of their conducting properties is central to advancing the field. In the present contribution, we highlight the fundamental differences between single-molecule and ensemble junctions focusing on the fundamentals of transport through molecular clusters. In this way, we elucidate the collective behavior of parallel molecular wires, bridging the gap between single molecule and large-area monolayer electronics, where even in the latter case transport is usually dominated by finite-size islands. On the basis of first-principles charge-transport simulations, we explain why the scaling of the conductivity of a junction has to be distinctly nonlinear in the number of molecules it contains. Moreover, transport through molecular clusters is found to be highly inhomogeneous with pronounced edge effects determined by molecules in locally different electrostatic environments. These effects are most pronounced for comparably small clusters, but electrostatic considerations show that they prevail also for more extended systems.

M olecular electronics aims at realizing electronic devices by contacting nanoscale assemblies of molecules between metallic electrodes. 1,2 A key goal is to meet the increasing technical demands of miniaturization. Beyond that, new device types are sought after, exploiting the enormous variety of conceivable systems arising from the flexibility of chemical design. Conceptually, the field of molecular electronics can be divided into two branches, namely single molecule electronics, where junctions ideally consist of an individual molecule 3,4 and molecular ensemble electronics 5,6 comprising junctions with large numbers of molecules or even extended monolayers with a quasi-infinite number of molecules in parallel. While there has been tremendous progress in understanding charge-transport in each of these fields, the transition between single molecule and monolayer junctions is rarely addressed. 7−10 The main conceptual difference between idealized molecular and ensemble (respectively, monolayer) junctions is that in the latter case the interaction of individual molecules becomes important. It determines the "scaling" of charge-transport properties with the number of molecules the junction contains. This issue has been quite controversially discussed in literature: While in some experimental studies, 11−13 the conductance per molecule has been observed to scale directly with the number of molecules in the junction, in other cases the current per molecule in single molecule junctions has been found to be orders of magnitude larger than for the monolayer. 14 Also, several theoretical studies 7−9,15−20 reported that the interaction between molecules can "help or hamper" 8 charge-transport through molecular junctions.
The situation is further complicated by the fact that in singlemolecule junctions the molecules might not exist as isolated entities and several molecules bond to the electrodes simultaneously, even if they might not fully bridge the gap. Conversely, even when studying transport through monolayers, often only a few hundred molecules are contacted simultaneously (e.g., in conducting-probe atomic force microscopy experiments 21 ) or transport occurs only through relatively small regions of the monolayer rendering only a tiny fraction electrically active (like in EGaIn junctions 22,23 ). In some experiments, contact is also intentionally made to comparably small clusters of organic molecules, for example, when growing them on metal nanoparticles. 24 When discussing the collective properties of a molecular ensemble, one generally has to distinguish between effects originating primarily from quantum-mechanical interactions 7−9,15,17−19 and effects caused by mere electrostatics due to polar elements within the junction, so-called collective electrostatic effects. 10,25,26 The latter generally arise from the combined electric fields of periodically assembled neighboring dipoles, which significantly affect the electrostatic potentialenergy landscape. 27−30 In molecular junctions they emerge "naturally" from polar docking groups and from interface dipoles arising due to the bonding of the molecules to the leads; 26 alternatively, they can be intentionally triggered by incorporating polar elements into the molecular backbones. 10 They crucially affect the alignment between the molecular transport channels and the Fermi level of the metal, which massively changes the current per molecule, potentially even switching the charge-transport polarity. 10,25 Collective electrostatic effects have also been found to shift measured transition voltages in monolayer junctions and through the spatial localization of charge-transporting states they can even cause rectification. 31−33 In this Letter, we address the "transition-region" between single molecule and ensemble junctions focusing on charge transport through molecular clusters of varying size. We show that a number of potentially unexpected effects occur even for idealized systems. In these, to isolate electrostatically triggered effects from other factors impacting ballistic transport, we assume flat electrodes and clusters consisting of equivalently oriented molecules in identical configurations and at equivalent docking sites. These junctions mimic a defect-free low-temperature situation, which is crucial for proving fundamental insights. Thus, fluctuations in the above-mentioned geometrical parameters are not considered here despite previous observations that they can have a profound impact on the ballistic transport properties. 34−40 In this way, on the basis of first-principles quantum transport simulations for junctions containing up to 16 molecules and electrostatic considerations for larger clusters we find that (i) the scaling of charge transport properties with the number of molecules in the junction is far from trivial, highly nonlinear, and massively affected by the above-mentioned collective electrostatic effects; (ii) collective electrostatic effects crucially impact transport already for comparably small cluster sizes; and (iii) transport in clusters is highly inhomogeneous with pronounced edge effects, which are linked to molecules in different electrostatic environments due to different numbers of neighbors.
The molecular junctions we study are based on derivatives of the prototypical "Tour-wire" 41 molecule (1,2,-bis(2phenylethynyl)benzene), which are bonded to Au electrodes via the commonly used anchoring groups pyridine, thiolate (−S), and isocyanide (−NC), as shown in Figure 1a.
The type of anchoring group is known to have a profound impact on the degree of coupling between molecules and leads. They also determine the level alignment, that is, the relative energetic position of the transmission channels with respect to the Fermi level, 42−47 due to interfacial potential steps arising from the superposition of the electric fields caused by the dipoles of the anchoring groups and the interfacial charge rearrangements. 26 For modeling densely packed monolayers employing periodic boundary conditions, we considered one molecule in a p(2 × 2) Au(111) surface unit cell (see solid box in Figure Figure S1 in the Supporting Information.

Nano Letters
Letter 1c), while the isolated (i.e., single) molecule is approximated by one molecule in a (8 × 8) Au(111) surface unit cell (i.e., at 1/ 16 coverage; see dashed box in Figure 1c). To set up the differently sized clusters of molecules, we start from the single molecule unit cell and consecutively add molecules while enlarging the electrodes and, concomitantly, also the unit cell. In this way, the distance between the outermost molecules of the clusters remains independent of the cluster size (see Figure  S1 in the Supporting Information). The clusters explicitly considered in the quantum-mechanical simulations contain 1, 2, 3, 4, 9, and 16 molecules. The structure of the latter is shown in Figure 1d and contains ∼2900 atoms per unit cell, which turns out to be a practical upper limit for doing the transport calculations (see Supporting Information). Ballistic charge transport calculations were performed using a recently improved version of the TranSIESTA/TBTrans suite 48,49 in conjunction with geometries optimized using VASP. 50 Further details can be found in Methods and in the Supporting Information.
To discuss the main effects we choose the pyridine-linked system shown in Figure 1a. This choice is motivated by the observation that the comparably weak coupling between the pyridines and the electrodes 26,51,52 allows a clear distinction between well-resolved transport channels. Note that the used geometry corresponds to the "low-conductance" mode of pyridine-linked junctions. 51,52 Conceptually similar trends would be expected for the "high-conductance" mode containing adatoms directly below the pyridine group. The main difference there is, however, that Fermi-level pinning (discussed below for the 16 molecule cluster and the monolayer) would set in at already significantly smaller cluster sizes due to the smaller barrier between the lowest unoccupied transport channel and the Fermi energy. 26 The current−voltage (I−V) characteristics calculated on the basis of the zero-bias transmission functions are shown in Figure 2a for the differently sized pyridine-linked clusters. For comparison, also the data for the single molecule and monolayer junctions are included. Overall, one observes a sharp onset of the current that shifts to smaller voltages with increasing cluster size. For the monolayer and the two largest clusters, an immediate steep increase of the current is obtained. These data clearly show that the current per molecule at a given voltage changes dramatically with cluster size.
To understand that evolution, it is useful to analyze the (zero-bias) transmission functions, which directly reflect the energetic alignment of the transmission peaks corresponding to the molecular states relative to the Fermi level of the electrodes. These are shown in Figure 2b,c in the energy range of the lowest unoccupied molecular orbital (LUMO), as there one finds the dominant conductance channels for the pyridinelinked junction. Transmission functions over a wider range can be found in Figure S5 in the Supporting Information. For the single molecule junction (black line), the transmission feature associated with the LUMO is represented by a narrow peak at around 0.8 eV. For two molecules in the junction, it splits and the lower-energy feature is shifted by as much as 0.22 eV. Upon further increasing the number of molecules in the cluster this trend continues. For 16 molecules in the junction, the transmission features spread out over a wide energy range (about 0.3 eV) resembling the situation for the monolayer (Figure 2c) and the net shift between the lowest-energy transmission features in the single molecule limit and the 16 molecule cluster amounts to 0.8 eV.
The energetic shifts of the transmission peaks are a consequence of changes of the electrostatic energy in the region of the molecules caused by the fields arising from the dipoles associated with the polar docking groups and bonding-

Nano Letters
Letter induced charge rearrangements. These are particularly pronounced when several dipoles are arranged in an ordered fashion, as it is the case for the monolayer 27−29 and to a lesser extent also in the molecular clusters. This is a mere consequence of the Helmholtz solution to Poisson's equation, which shows that electrons have to overcome a step in the electrostatic energy, when passing through a regular assembly of dipoles. At this point it is important to realize that this effect is not a peculiarity of pyridine-docked systems but is universal, as docking groups containing heteroatoms will always generate local dipoles and bonding to a substrate is inevitably associated with interfacial charge rearrangements. The situation for molecular clusters is, in fact, reminiscent of what has been observed when reducing the coverage of a homogeneous monolayer. 26 There, the shift has been found to be directly proportional to the dipole density. Upon closer inspection, however, one notices a fundamental difference between clusters and homogeneous low-coverage layers. While in the latter case each molecule feels an identical field generated by the neighboring molecules, in the cluster case the number and position of neighbors and, concomitantly, the net field varies depending on a molecule's position. 28 The variation of the electrostatic energy landscape each molecule is exposed to results in variations of the energetic positions of the transmission features associated with specific molecules within the cluster.
These electrostatic "edge-effects" together with the quantummechanical coupling between the molecules are responsible for the emergence of multiple peaks in the transmission functions shown in Figure 2b,c. In the following, the situation is discussed in more detail for the pyridine-linked 9-molecule cluster (see Figure 3a). There, the transmission peaks derived from the lowest unoccupied molecular states group into three main features centered at 0.07, 0.17, and 0.35 eV above the Fermi level.
To analyze the origin of that peak splitting, we will first discuss the impact of electrostatic effects. Subsequently, we will analyze the role played by the quantum-mechanical coupling between the molecules. The local densities of states (LDOS) associated with the transmission peaks are shown in the top region of Figure 3a. They clearly show that the lowest-energy feature is mostly associated with transport through the central molecule. This molecule is surrounded by eight neighboring molecules with their associated dipoles and, hence, the electrostatically induced shift is largest. The next feature can be associated primarily with molecules at the edges of the cluster, while the molecules at the corners have the smallest number of neighboring molecules and, consequently, experience also the smallest fields resulting in the least shifted transmission feature. Interestingly, exactly the same succession of transmission peaks is repeated around 1 eV for the transmission features derived from the next unoccupied molecular state (see Figure 3a).
The pronounced edge-effects are also clearly visible in Figure  3b, where we plot the cluster-size dependent shifts of the lowest and highest LUMO-derived transmission features relative to the shift between single molecule and the monolayer (for details see figure caption). One sees that the energy of the transmission peaks associated with the central molecule(s) (i.e., the squares in Figure 3b) approaches the monolayer limit faster. This results in a continuous increase of the energy range in which transmission features exist (see shaded yellow area in Figure 3b). ES is the electrostatic energy difference between the pair of dipoles and two infinitely extended dipole sheets. The energies are determined in the plane between the dipole sheets; squares refer to positions between the centers of the dipole clusters, circles to positions between corners (see Figure S7). These data represent the electrostatic analogue to the quantum-mechanical results depicted in panel (b) albeit extended to much larger cluster sizes and neglecting the extent of the electronic states perpendicular to the planes of the dipoles. The red vertical line indicates the 16 dipole cluster, which corresponds to the largest cluster calculated also quantum-mechanically. As only relative quantities are reported, this plot is generally valid independent of the magnitude of the dipoles or the distance

Nano Letters
Letter To further analyze the situation and to extend the discussion to much larger clusters, we devised a simple electrostatic point dipole model. As shown in Figure S7 in the Supporting Information, we describe the electrostatic situation by two opposite square 2D point dipole arrays, where each array mimics the dipoles due to the docking groups and the bond dipoles. Then, we calculated the shift in the electrostatic energy an electron would experience in the middle between the two arrays at the position of the central molecule as well as at the position of a corner molecule. Figure 3c shows those shifts relative to the shift between the single molecule and the monolayer equivalents. The obtained data confirm several of the trends discussed already above for the small clusters calculated quantum-mechanically. The energy at the position of the central molecule gradually shifts toward the monolayer limit. The shift at the corner of the cluster is much smaller highlighting again the boundary effects expected for transport through molecular clusters. In this context it should, however, be noted that in extended clusters the effect is at least partially offset by a decreasing ratio of molecules (dipoles) at the border of the cluster relative to molecules inside the cluster (see Figure  S8 in the Supporting Information).
While the qualitative conclusions from the electrostatic model match those of the quantum-mechanical calculations, the actual situation is more complex. In the quantum-mechanical simulations the monolayer limit is nearly reached for the 16molecule cluster. In the electrostatic model, much larger structures would be needed to achieve that. One reason are depolarization effects, 28,30 which are less pronounced for border molecules and which are not considered in the electrostatic model. The main reason, however, lies in an electronic peculiarity of the pyridine docking group. For pyridine-linked junctions beyond a certain cluster size, the LUMO aligns with E F , 26,53 an effect known as Fermi level pinning. 54,55 In the present case that means that for clusters containing more than 16 molecules on purely electrostatic grounds the shift of transmission features to lower energies would continue. In the actual junction this is, however, prevented by interfacial charge rearrangement avoiding that unoccupied states of the cluster are pushed below the Fermi level.
A further complication is that electrostatic effects cannot be the only reason for peak splitting and the energetic broadening of the features. This can, for example, be inferred from the observations in Figure 2b,c that already in the two-molecule cluster consisting of "electrostatically equivalent" molecules a pronounced peak-splitting occurs and that in the limit of a continuous monolayer consisting of identical molecules the low-energy transmission feature has a width of ∼0.3 eV (similar to the peak splitting in the cluster; see Figure 2c). In those cases, quantum-mechanical coupling between the molecules has to be responsible for the broadening. To distinguish between the impact of that coupling and the electrostatic edge effects discussed above for the 9-molecule cluster, we pursued a dual   Figure S5. Equivalent plots for the pyridine-docked system can be found in the Supporting Information.

Nano Letters
Letter approach that is described in detail in the Supporting Information: On the one hand, we devised a tight-binding model parametrized on the basis of the quantum-mechanical data. There one can test the impact of the two effects by selectively switching off either the electrostatically triggered asymmetries (by choosing identical on-site energies) or the quantum-mechanical coupling (by setting the transfer integral to 0 eV). On the other hand, we calculated orbitals in molecular clusters isolated from the electrodes in which we manipulated the electrostatic effects by changing the type of polar substituents. Both approaches yield the same conclusions, namely that already the quantum-mechanical coupling gives rise to a peak splitting, respectively, broadening. The electrostatic edge effects are, however, crucial for the final situation, as they increase the magnitude of the splitting and, most importantly, determine the localization and energetic order of the states. The latter occurs fully in line with the above discussion of variations of the energetic shifts caused by different numbers of neighboring molecule-related dipoles. The only case in which, naturally, edge effects are not relevant is the monolayer, where the entire broadening arises from quantum-mechanical coupling, which has the strongest impact on the bandwidth for the infinitely extended system (see Supporting Information).
As a last step, the trends obtained for the pyridine-linked system shall be compared to those of thiolate-and isocyanidelinked junctions. The corresponding transmission functions for increasing cluster sizes are shown in Figure 4a,b in the range of the states closest to E F , as those determine conduction at small biases. For the thiolate-linked system, these are the highest occupied states (resulting in p-type transport) and for the isocyanide-linked system the lowest unoccupied states (yielding n-type transport). Compared to the pyridine-docked systems, one observes two differences: the shifts between single molecule and monolayer are smaller, which is a consequence of the smaller interfacial dipoles; 26 additionally, due to the stronger quantum-mechanical coupling between molecular and metal states, the transmission peaks are significantly broadened, which makes a reliable determination of the positions of individual transmission features difficult, if not impossible. It also complicates the identification of edge effects. Thus, to obtain quantitative trends also for these systems, we calculate the (zero-bias) conductance, G(E F ), per molecule, which is the product of the transmission at E F and the quantum of conductance (2e 2 /h). The results are plotted as a function of cluster size Figure 4c,d (an equivalent plot for the pyridinelinked system can be found in Figure S11 in the Supporting Information).
For all systems, the (zero-bias) conductance per molecule displays a significant dependence on the cluster size. Interestingly, G(E F ) per molecule increases with the number of molecules in the junction for the isocyanide-linked junctions, while it decreases for the thiolate-linked ones. The reason for that is that for all docking groups studied here, collective electrostatic effects shift transmission features to lower energies. Consequently, whenever transport occurs through unoccupied states, this results in an increase of the conductance with cluster size; conversely, a decrease is observed, when transport ensues through occupied states. In passing we note that a possible strategy for breaking this pattern would be the incorporation of polar elements with opposite dipole orientations into the molecular backbone. This has, for example, been shown for pyrimidine-containing, tour-wire based junctions. 10 In conclusion, we performed first-principles charge-transport calculations through molecular clusters of increasing size in order to bridge the gap between single molecule and large-area molecular junctions. At a given voltage, we find hugely differing currents per molecule as a function of cluster size, where it depends on the docking group whether larger clusters are more or less conductive. Moreover, pronounced edge effects are observed, which results in highly inhomogeneous transport through the clusters. These observations are identified as a consequence of so-called collective electrostatic effects, which arise from the superposition of the fields caused by dipoles associated with docking groups and interfacial charge rearrangements. They shift the molecule-derived electronic states within the junction in energy, which directly translates into a shift of the associated transmission features. The edge effects are a consequence of the dependence of these shifts on the location of a specific molecule within the cluster in combination with the quantum-mechanical coupling between the molecules. These considerations show that even for idealized junctions the electrostatic environment of the conducting molecules is a very important factor for understanding transport in any junction, especially in those that contain more than a single molecule.
Methods. To determine the electronic properties of the systems and optimize their geometries, we performed periodic calculations within the framework of density functional theory (DFT) using the VASP 50 code in conjunction with the PBE 56 functional and a plane-wave basis set (cutoff: 274 eV). We optimized the structure of the monolayer junction (considering three Au layers on each side of the junction) according to the procedure described in ref 10. The geometries of the clusters were not optimized since no significant changes in molecular conformation are expected considering that the molecules are suspended between two electrodes. To obtain current−voltage characteristics and zero-bias transmission functions we used the recently improved TranSIESTA/TBTrans suite. 48,49 We employed a double-ζ polarized orbital basis set (DZP) accompanied by a PAO energy shift of 0.001 Ry, which we found to be crucial to correctly reproduce the level alignment obtained from highly converged VASP calculations (see the Supporting Information of ref 10). Due to the system size (up to 2900 atoms) and the related computational effort (see Figure S4 in the Supporting Information), we only present transmission and current calculations using the Kohn−Sham Hamiltonians as calculated by SIESTA as input to the transport calculations in TBtrans. We have performed equilibrium Green function calculations with TranSIESTA to assert that the physics are unchanged due to sufficient screening towards the bulk gold electrodes, see Figure S3 in the Supporting Information. The junction consists of the molecule and three Au layers on each side. In the terminology of transport calculations this is called the "central region". For electronic transport calculations, on each side three more Au layers are added, which represent the semi-infinite leads in the Green function calculation scheme, also referred to as "electrodes" (this setup is depicted in Figure S2 in the Supporting Information). The zero-bias conductance, G(E F ) = T(E F )·G 0 , was evaluated from the zero-bias transmission function T at the Fermi level E F ; G 0 = (2e 2 /h) is the quantum of conductance. Xcrysden 57 and VMD 58 were used for graphical visualization. For full details on the applied computational methodology and numerical parameters used in our calculations, see the Supporting Information.