Influence of the Contact Geometry and Counterions on the Current Flow and Charge Transfer in Polyoxometalate Molecular Junctions: A Density Functional Theory Study.

Polyoxometalates (POMs) are promising candidates for molecular electronic applications because (1) they are inorganic molecules, which have better CMOS compatibility compared to organic molecules; (2) they are easily synthesized in a one-pot reaction from metal oxides (MO x ) (where the metal M can be, e.g., W, V, or Mo, and x is an integer between 4 and 7); (3) POMs can self-assemble to form various shapes and configurations, and thus the chemical synthesis can be tailored for specific device performance; and (4) they are redox-active with multiple states that have a very low voltage switching between polarized states. However, a deep understanding is required if we are to make commercial molecular devices a reality. Simulation and modeling are the most time efficient and cost-effective methods to evaluate a potential device performance. Here, we use density functional theory in combination with nonequilibrium Green's function to study the transport properties of [W18O54(SO3)2]4-, a POM cluster, in a variety of molecular junction configurations. Our calculations reveal that the transport profile not only is linked to the electronic structure of the molecule but also is influenced by contact geometry and presence of ions. More specifically, the contact geometry and the number of bonds between the POM and the electrodes determine the current flow. Hence, strong and reproducible contact between the leads and the molecule is mandatory to establish a reliable fabrication process. Moreover, although often ignored, our simulations show that the charge balancing counterions activate the conductance channels intrinsic to the molecule, leading to a dramatic increase in the computed current at low bias. Therefore, the role of these counterions cannot be ignored when molecular based devices are fabricated. In summary, this work shows that the current transport in POM junctions is determined by not only the contact geometry between the molecule and the electrode but also the presence of ions around the molecule. This significantly impacts the transport properties in such nanoscale molecular electronic devices.

Introduction

The semiconductor industry is driven by the ambition to extend Moore’s law. However, further miniaturization of CMOS technology is reaching its inherent limitations.[1] There is a significant research effort to overcome these limitations through a new approach. One such approach is from “the bottom up”, creating electronic circuits from single molecules. This is the goal of the field of molecular electronics.[2] While a lot of progress has been made since the discovery of molecular rectifiers in 1974,[3] a commercial molecular based electronic device is still out of reach. Moreover, the understanding of electronic transport at the molecular level is very challenging from experimental point of view.[4] For this reason various computational modeling techniques such as density functional theory (DFT)[5] are an efficient way to probe molecular electronic systems and gain a more fundamental understanding of electron transport and charge distribution.

A promising candidate for molecular electronic applications is polyoxometalates (POMs). POMs are inorganic anionic clusters typically made from high valent group V and VI transition metals (predominantly Mn and W) linked together by oxo ligands to make 3D networks.[6] The structural variety of POMs in addition to their rich redox chemistry, high conductivity, and diverse reactivity has found them wide-ranging applications in photochemistry,[7] medicinal chemistry,[8] sensors,[9] sustainable energy,[9−11] catalysis,[12] nanotechnology,[13] and electronics.[14,15]

The idea of a POM-based memory device was presented in two ground-breaking concept papers published in 2006 and 2013, which described their unique characteristics and their great potential as elements of molecular memory.[13,15] This concept was not only very well received, but it also made the community think in a new way about how nanosystems might look in the future.[16,17] As devices shrink and the molecular/atomic level approaches, vast problems with variability and lifetime are becoming more of a reality. In 2014, a new approach to the design of flash memory was proposed, incorporating polyoxometalate clusters within the floating gate to enhance the performance parameters.[18] Further device simulations have been explored since based on this design.[19−21] However, a fundamental understanding of electron transport through POM molecules is lacking but pivotal if we are to reach our goal of designing working commercial molecular based memories.

In this work, the transport properties of a specific tungsten-based POM molecule, [W18O54(SO3)2]4–, are explored theoretically by computational modeling of POM molecular junctions. Single molecule junctions are the simplest devices for exploring the properties of single molecule electronics both experimentally and theoretically, with many examples in the literature.[22−27]

We employ DFT in combination with the nonequilibrium Green’s function (NEGF) method to probe the transport properties of [W18O54(SO3)2]4– molecular junctions. DFT has been used successfully to study the structure of POMs previously.[28] NEGF has proven to be a successful method for understanding the electron transport at the molecular level for a wide range of molecular systems.[26,27]

The aim of this work is as follows: first, to establish a link between the underlying electronic structure of the POM molecule and its predicted transport properties; second, to explore the effects of different electrode–molecule contact geometries on the transport; and third, to predict how the POM will behave under an applied bias.

Finally, because of the anionic nature of POMs, countercations are required for charge balance and to isolate the molecule in the pure phase. The roles of countercations in POM chemistry and the chemistry of other metal–oxo clusters have been overlooked as pointed out in the review by Nyman et al.[30] However, several studies show that countercations play a larger role than simple charge balance. In some applications, such as heterogeneous catalysis, their role is critical.[31] We believe that the POM counterions are not insignificant and could influence the properties of the molecule in the context of single molecule electronic transport of POMs. This is explored theoretically for the first time here in this work.

Computational Methodology

The calculations of the POM molecule and POM junctions were performed by using the QuantumATK-2018.06 software. For the counterion-containing systems the newer QuantumATK-2019.12 software was used.[32] The electronic properties of the molecule were calculated by using the generalized gradient approximation (GGA) and the BP86 functional with spin polarization. This was chosen as it has successfully described POMs in previous publications.[15,28]

The SG15 pseudopotential[33] and a medium basis set (comparable to double-zeta polarized) were employed and satisfactorily reproduced the geometry and electronic properties of [W18O54(SO3)2]4–. Figure shows the structure of the [W18O54(SO3)2]4– POM cluster. For the gold electrodes, single-zeta polarized basis set was used to reduce the computational cost. The geometry of the molecule was optimized in the gas phase, and its energy levels were calculated.

Figure 1

POM molecule in this study, with atom labels for bond data. W = blue, O = red, and S = yellow. O1, O2, and O3 are all oxygen atoms which are labeled according to the oxo band in the molecule. O1 = capping, O2 = central, and O3 = upper oxo bands.

The molecule was then sandwiched between two Au (111) surfaces in several possible configurations. The main molecule–electrode contact was the binding of an oxygen atom to a face-centered cubic (fcc) site in the gold surface; because of the geometry of the POM molecule when this contact is made, neighboring oxygen atoms are close enough to also interact with the surface, resulting in multiple molecule–electrode contacts. The Au (111) surfaces act as the source and drain electrodes of the molecular junction.

For the transport calculations the NEGF method was then used to calculate the transmission spectra and the current flow through the molecule under applied drain bias. For the transport simulations a large transverse k-point sampling is required for accuracy. Thus, a Monkhorst–Pack grid of 3 × 3 × 200 k-point sampling was implemented.

According to the NEGF formalism, the transmission coefficients (T(ε)) at a particular energy are calculated via the equation[34]where ΓL/R is the broadening function due to the coupling between the left/right electrode and the central scattering region and GR is the retarded Green’s function. The current through the junction is computed over a range of applied bias and is calculated by using the Landauer–Buttiker formula:[35]where fL/R is the Fermi function for the electrodes, e is the charge of the electron, and h is Planck’s constant.

For the counterion-containing junctions, the same procedure was repeated. An increased electrode size was required to encapsulate the whole system. The positioning of the counterions was taken from X-ray diffraction data.

Results and Discussion

Justification of DFT Settings

All DFT studies rely on the choice of the functional and basis set and its ability to accurately describe the system of interest. In this case, it is vital that the HOMO–LUMO gap is accurately calculated and the molecular geometry is successfully reproduced. Several functionals and basis sets were tested on the molecule to reproduce the theoretical HOMO–LUMO gap and the experimental bond lengths and angles reported by Vilà-Nadal et al.[15]

As is evident from Table , the BP86 functional, with SG15 pseudopotential and medium basis set, satisfactorily reproduces the bond lengths and bond angles of the theoretical and experimental benchmarks. There are only slight variations, most notably with the bond angle. Furthermore, DFT at this level reproduces the HOMO–LUMO gap with a difference of only 0.02 eV. On the basis of these results, we are satisfied that this level of DFT theory can describe the POM system for transport calculations.

Table 1

Comparison of Selected Bond Lengths (Å) and Bond Angles (deg) of DFT Theory Used in This Work and the Vilà et al.[15] Paper Used as a Benchmark

	W1–O1	W1–O3	W2–O2	W2–O2–W2	ΔEHOMO–LUMO (eV)
experimental	1.89	1.93	1.92	145.58	N/A
theoretical	1.92	1.93	1.92	147.55	1.87
BP86, SG15 medium basis set	1.92	1.93	1.92	148.75	1.89

Single POM Molecular Junction Studies

One of the major challenges in studying single molecule junctions is a lack of structural detail at the molecular level, that is, how the molecules are actually attached to the electrodes. The first goal in this work is to assess the effect of molecule orientation and contact geometry of different junctions with the same molecule on the current and link this to the underlying electronic structure. In experimental conditions, because of the thermal vibrations of the molecule and the stochastic nature of the binding process of the molecule to the contacts, it is practically impossible to accurately determine the exact position of the molecule in relation to the contacts. In reality, there are multiple orientations which all have similar probabilities of forming. From the simulation point of view, it is computationally expensive and time-consuming to study all possible of configurations.

Therefore, we report three configurations that we believe cover the different “extremes” of effects in the junctions. It is known that the nature of the molecule to electrode contact has a significant impact on the current flow in molecular electronic systems.[36,37] Hence, to explore this, junction setups with 1, 2, and 3 direct contacts of oxygen to gold (Au) on either side have been modeled. The configurations discussed here are shown in Figure . For all configurations we have calculated the binding energies which can be found in the Supporting Information. Our results show that for all systems the binding energy is positive, and has much higher value than the thermal energy which suggests a very strong absorption of the molecule to the Au surface.

Figure 2

Three Au–POM–Au configurations studied.

The first step to understanding the molecular junctions is by exploring its electronic structure. To that end, the energy levels, specifically the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), were calculated and visualized. It must be stressed that the analysis is qualitative. However, GGA functionals have been shown to be reasonably accurate in calculating HOMO–LUMO gap.[38] As the method is consistent for all configurations, their energetic properties can be compared.

Figure shows the relative alignment of the HOMO, HOMO–1, HOMO–2, and the LUMO for all three configurations with respect to the Fermi energy of the electrodes. It is clear from the figure that the molecule to electrode configuration influences the energy levels of the molecule. First, the HOMO–LUMO gap in all cases is very close to the value of the isolated molecule, 1.89 eV (see Table ). Second, the alignment of the energy levels to the Fermi energy (Ef) is different for each configuration. For configurations 1 and 2, the Fermi energy lies closer to the HOMO than the LUMO, whereas for configuration 3 it lies approximately halfway between the HOMO and LUMO. This suggests the increased contact strength to the electrodes stabilizes the HOMO level and pulls the LUMO closer to the Fermi energy. This observation is directly linked to the number of bonds between the POM cluster and the electrode. For instance, configuration 1 has three contacts (bonds) between the POM–Au electrode, and this would result in a stronger contact.

Figure 3

Comparison of HOMO–LUMO relative energy positioning and their respective frontier Kohn–Sham molecular orbitals. The blue dotted line marks the positioning Fermi energy of the electrodes.

Third, an important observation is that the symmetry of the molecular orbitals of the HOMO and LUMO levels for configuration 2 and 3 is identical with the energy levels in the isolated molecule.[15] The HOMO is delocalized around the [SO3]2– moieties. The LUMO level is delocalized around the tungsten cage and show d-like symmetry and character. For configuration 1 the HOMO–2 resembles the HOMO of the isolated molecule, while the LUMO level is identical. However, the HOMO–2, HOMO–1, and HOMO all lie within 0.02 eV of each other and thus are degenerate states. Similar degeneracy is seen for configuration 3, suggesting that the “horizontal” orientation brings these energy levels closer. This is consistent with the fact that for these configurations the energy states show that the electron density lies on part of the central oxo bands. Whereas for configuration 2 the densities for these levels lie on the upper and lower oxo bands (see Figure ), and thus the levels are more discrete and as a result different in energy.

What is also highlighted by Figure is the influence of the configuration of the device on the HOMO–LUMO gap. From the engineering point of view the HOMO level can be associated with the valence band edge (EV) and the LUMO level with the conduction band edge (EC). For configuration 2, “the vertical” orientation results in a smaller HOMO–LUMO gap (1.71 eV) compared to the 1.89 eV of the isolated molecule. However, for the “horizontal” orientation, the HOMO–LUMO gap increases (2.17 and 2.12 eV). This suggests the orientation of the molecule to the electrode influences whether the gap increases or decreases. The orientation of the molecule seems to affect mostly the HOMO and near-HOMO levels, resulting in them being pulled closer together and stabilized; hence, the gap to the LUMO widens. This can be expected because these energy levels are full of valence electrons which contribute to the forming of chemical bonds between the molecule and the electrode.

A useful way of analyzing the electronic structure of the molecular junction is by computing the local density of states (LDOS). It shows where in the molecular device there are available energy states and the relative occupation. The LDOS for all configurations projected onto the device coordinates is shown in Figure . These were calculated at zero drain bias, and hence the Fermi levels for the source (εL) and the drain (εR) are the same. Also, the source Fermi level is positioned at 0 eV to facilitate comparisons between the structures.

Figure 4

Projected local density of states calculated at zero bias for (a) configuration 1, (b) configuration 2, and (c) configuration 3. The region spanning 10–23 Å is the central region of the molecule. The pink regions are where there is high electron density. The green dotted line is the Fermi energy, which is set at 0 eV.

In general, the pink regions are regions of high electron density and the black regions are without any states. In Figure the Z-axis is the transport direction, and the continuous regions at 0–10 and 22–32 Å are the LDOS for the electrodes, while the region in-between (10–22 Å) is the molecule.

The PLDOS echoes the results for the energy level alignment in Figure and shows significant differences across all three configurations. This highlights the effect of molecule–electrode contact and orientation on the electronic structure of the device. Increased contact to the electrodes shifts the energy levels and subsequently the available density of states downward in energy. This is evident as for configurations 1 and 2 the LDOS lie closer to the Fermi energy. The gaps between the states on either side of the Fermi energy are consistent with the calculated HOMO–LUMO gaps.

What is also clear by analyzing the PLDOS is the available DOS are highly localized around the molecule. Thus, the positioning of the localization is highly dependent on the configuration. Specifically, in the states around the Fermi level, which are critical to transport, there are clear differences between the configurations, with split resonances across the molecule. It shows that the contact to the electrodes and molecular orientation dictates both the energy and atomic localization of the DOS of the molecule.

Figure shows the theoretical transmission spectra calculated at zero bias. It shows the different energies at which electrons can scatter through the device and the relative probability (T(E)). The NEFG formalism treats the electrons as waves, so the peaks in the transmission spectra can be associated with constructive interference between the electrons in the electrodes and the molecule when they are interacting. The transmission spectra are linked to the LDOS data in Figure .

Figure 5

(a) Transmission spectra for all three configurations. (b) Transmission spectra with the dominant transmission eigenstates visualized for selected peaks. All energies are relative to EF, which is at zero.

In Figure , all energies are calculated relative to the Fermi energy (which is set as zero) at zero drain bias. Despite differences in T(E) and energy positions of the peaks, the spectra are similar for all three configurations. This demonstrates a clear link between the electronic structure of the POM and its transport. However, because of the different electrode–molecule configurations, differences arise in the transmission spectra, consistent with the differences in energy levels discussed previously. Plotting the transmission spectra as E – EF (eV) vs log(T(E)) provides a clear visualization of resonance broadening near the Fermi energy. Interested readers can see the log scale transmission spectra in the Supporting Information.

There is an evident trend in Figure a regarding the effect of molecule–electrode contact strength. Focusing on the first peaks above the Fermi energy (EF = 0.0 V) for each configuration, the closest peak is for configuration 1, followed by configuration 2 and then configuration 3. A similar trend is seen for the other peaks, and below EF, the peak for configuration 1 is lowest in energy. This reflects the HOMO–LUMO positioning in Figure . Configuration 1 has three oxygen atoms bonded to the gold electrode on either side, while only two for configuration 2 and just one contact either side for configuration 3 (see Figure ). More points of contact to the electrodes lowers the energy of the transmission modes intrinsic to the molecule. Or this can be explained by creating more than one conducting channel which allow the electrons to move from the source to the molecule. As a result, the transmission spectra and transport are strongly dependent on the molecule–electrode geometry.

A similar dependence of the molecule–electrode geometry on the properties of tunnel junctions has been shown in a series of studies exploring predominantly organic and organometallic molecules including junctions with a series of different electrodes. Our results show that the same is true for POM-based junctions.[39−41]

Another key disagreement between the transmission spectra of the configurations is the differences in T(E), i.e., the intensity of the peaks. Focusing on the region of 0–1.5 eV in Figure a, the peaks for configuration 1 are much larger than that of 2 but similar to the size of the peaks for configuration 3. This indicates that at these energies the “horizontal” orientation provides more favorable alignment of the transport modes than the “vertical” configuration. The peaks closest to EF will be the transmission modes which will contribute to the transport under applied bias; as transmission is directly related to conductance, it can be expected that the horizontal orientation results in a higher conductance.

Below the Fermi energy, the peaks correspond to the filled energy levels of the molecule and gold. The differences in T(E) are not so prominent. This can be explained by the fact that these states are already filled, and so the probabilities of the electrons scattering at these energies are all similar. Given that they are low-lying, a high applied bias would be required for these peaks to contribute to the current flow.

In Figure b, the dominant transmission eigenstates responsible for the transmission of the first peaks above and below the Fermi energy have been visualized for all configurations. The transmission eigenstates provide insight into where in the molecule the transmission is occurring and also the underlying electronic structure (the molecular orbitals responsible)—in effect, helping understand the link between electronic structure with the transport at the molecular level. The first peaks below the EF for all configurations (in the range −1 to −2 eV) unsurprisingly resemble the HOMO energy level of the molecule for all configurations.

The peaks above EF for all three configurations provide the most interesting results with regards to transport. All peaks resemble the LUMO of the molecule, indicating that the main transport pathways in the POM-based molecular junctions will be LUMO dominated at low drain bias. This is consistent with other studies of other POM-based molecular electronic transport.[17,42] In addition, the eigenstates are identical for all configurations, providing further evidence that the transport of a molecular system is dominated by the molecular electronic structure of the POMs. However, other features that influence the alignment of the molecular energy levels will influence the transport characteristics. Specifically, differences in molecule to electrode strength and in molecule–electrode contact orientation.

Figure shows the theoretical I–V curves for the three junction configurations and reaffirms the predictions from the zero bias transmission spectra. The computed curves predict low current flow between the bias of −1.5 to 1.5 V, with magnitudes ranging from 19.5 to 425 nA. The low current flow is consistent with the transmission spectra of the devices. This can be explained first by the size of the HOMO–LUMO gap, where there are no transmission modes. In addition, the transmission peaks are all sharp and narrow, as only a few energy levels/states contribute to the transmission in the junction.

Figure 6

Computed I–V characteristics for all configurations. The current is electronic current rather than conventional; i.e., positive current describes electron movement from the source (left electrode) to the drain (right electrode).

Unsurprisingly, configuration 1 shows the highest predicted current flow between the bias of −1.5 to 1.5 V, with it reaching 427 nA at 1.5 V and 387 nA at −1 V. Configuration 3 shows the lowest current flow, reaching 95 nA at 1.5 V and −67.5 nA at 1.5 V bias. Configuration 2 shows a higher current flow than 3 but significantly less than 1, with the predicted current reaching 184 nA. These results are highly suggestive that the contact strength between electrode and molecule most strongly influences the current flow of the molecular device, with more contacts resulting in a higher current.

The computed current is directly related to the number and magnitude of transmission modes within an energy range due to the applied bias (see the Landauer formula above). Hence, analyzing how the transmission spectrum changes under bias helps rationalize the I–V characteristics. Figure shows how the transmission changes under applied bias with peaks highlighted at ±0.5 and ±1 V. In general, with increased magnitude of the applied bias the intensity of the peaks (value of T(E)) decreases. This results in a lower than expected current. By comparison across the three configurations, in configuration 1 the energy window from the applied bias (yellow line) reaches transmission peaks under an applied bias of 1 and −1 V. Under negative applied bias the transmission peaks are not shifted as much as under positive bias; as a result, there are two small peaks included in the window. This explains why the current for configuration 1 at −1 V (387 nA) is predicted to be more than double the current at +1 V (179 nA).

Figure 7

Theoretical transmission spectra under applied bias (−1.0 to 1.0 V) for (a) configuration 1, (b) configuration 2, and (c) configuration 3. The yellow line corresponds to the energy range of the applied bias window.

In summary, in this section the transport properties of three different junction configurations for a POM molecule [W18O54(SO3)4]4– have been studied computationally. The DFT calculations demonstrate that the number of contacts between the molecule and the Au electrode determines the bond strength. A higher molecule to electrode contact strength (number of chemical bonds) lowers the energy levels, bringing the transmission peaks closer to the EF and, as a result, increasing current flow in the junction under bias. The molecular orientation to electrode has been shown to affect the magnitude of T(E) at energies close to the EF, with the calculations suggesting the “horizontal” configuration is the most favorable for transport.

Exploring the Effects of Counterions

Commonly in theoretical explorations of POM systems, the counterbalancing ions are ignored and assumed to be spectator ions important only for isolation of the system. However, the reality is that in a real device system the POMs will be surrounded by relatively large organic cations (see Figure ). We anticipate that from an electronic device point of view the influence of these counterions is not negligible. It has been shown that counterion dynamics are vitally important in molecular tunnel junctions. Studies by Han et al. and Goswami et al. report the stabilizing effect of counterions, which in turn have a crucial influence on the charge transport of organometallic based tunnel junctions.[43,44]

Figure 8

Illustration of [W18O54(SO3)2]4– POM molecules surrounded by tetrapropylammonium (TPA) counterions.

Because of this, it is vitally important to understand how counterions will affect the transport properties of [W18O54(SO3)2]4–and by extension other POM molecules. This is explored here theoretically for the first time.

[W18O54(SO3)2]4– is commonly isolated as a salt with four tetrapropylammonium cations (TPA).[45] In a proposed device structure these cations or similar ones would be present. To reduce the computational burden, the size of the alkyl chains was reduced to methyl groups, and thus the influence of tetramethylammonium (TMA) on the POM junction is explored. Previous device calculations modeled the presence of counterions as point charges.[18] To explore something similar, we have also investigated the influence of Cs+ counterions. Cs+ was chosen for two main reasons: Its ionic radius is similar to that of ammonium (which has similar properties to TMA).[46] Also, there are several examples of Cs+ POM salts in the literature.[30,31,47,48]

The main challenge with modeling counterions is realistic positioning of the counterions with respect to the molecule. In this study, we have taken the counterion positioning directly from X-ray crystallography coordinates of four TPA cations with respect to one molecule of [W18O54(SO3)2]4– within a unit cell. We accept that this is one arrangement of many possible and that because of the different sizes of TPA, TMA, and Cs+ there will be differences in position. However, for the purposes of this study it is a reasonable approximation to gauge the influence of counterions on POM transport properties.

A junction geometry like configuration 1 was chosen to study the counterion effects as it showed the best transport properties. To accommodate the presence of the counterions, the electrodes for the junction were increased (see Figure ). Benchmarking calculations show that increasing the size of the electrodes has a negligible effect on the transport properties of the system. Hence, any changes in predicted transport properties can be credited to the presence of counterions.

Figure 9

Junction configuration for systems including the counterions TMA (left) and Cs+ (right).

The effect of counterions on the energy levels of the system was assessed. A comparison between the two counterion systems and configuration 1 is shown in Figure . The Fermi energy has increased by 0.04 eV to −2.36 eV due to the increased size of electrodes. All comparisons have been adjusted for this. The spatial arrangement of the HOMO for the Cs+ and TMA system is identical with the HOMO for configurations 2 and 3 and to the HOMO–2 of configuration 1.

Figure 10

Comparison of HOMO–LUMO energies and the respective frontier Kohn–Sham molecular orbitals. The blue dotted line marks the Fermi energy dominated by the electrodes.

Interestingly, the results for the TMA+ and Cs+ are almost identical with respect to energy level positioning and symmetry of the orbitals. The results show that the presence of counterions influence the positioning of the energy levels, resulting in a stabilization, bringing down the energy much closer to the EF in comparison to the structures without counterions. The LUMO level for these systems lies much closer to the Fermi energy, while the HOMO level lies deeper. The results suggest a slight lowering of the HOMO–LUMO gap; however, the difference is not significant.

Another interesting effect of counterions is the breaking of the degeneracy of HOMO, HOMO–1, and HOMO–2. The HOMO is clearly higher in energy than HOMO–1 and HOMO–2 for the counterion containing systems. The symmetry of the orbitals is similar, localized on several of the oxygens of the cage, but they are not identical with the HOMO and HOMO–1 of configuration 1.

The LUMO for the systems are identical with configurations 1–3 and are delocalized over the central tungsten cage. The counterions do, however, stabilize the energy with respect to the Fermi energy. As a result, it lies ∼0.08 eV above EF. Given that, the transport of the POM has been shown to be LUMO dominated. This is significant for the transport properties of the system.

The PLDOS for counterion-containing systems at zero bias is shown in Figure . The calculated PLDOS is effectively identical for Cs+ and TMA with only slight differences in the first resonances below the Fermi energy. Compared to the single molecule systems, there are more available states above the Fermi energy and are all closer to the Fermi energy. Below EF, the density of states is found a lot lower in energy. The presence of counterions influence the local density of states of the junction, with available states being lower in energy, analogous to the effect on the energy levels.

Figure 11

PLDOS for the counterion-containing systems: (a) TMA-containing system; (b) Cs-containing system.

The theoretical transmission spectra for the counterion-containing junctions are shown in Figure a. Like the previous results, the spectra for TMA and Cs are almost identical. All peaks are in the same position, with only slight differences in T(E), possibly due to the respective sizes of the counterions. These results suggest that the two cations influence the transport properties of the POM in the same way or the nature of the counterion does not matter and the transport is influenced by the counterion positioning.

Figure 12

(a) Zero bias transmission spectra for configuration 1, with Cs+ and TMA counterions. (b) Visualized transmission eigenstates for the first four peaks above the Fermi energy. Comparison of the energy states for the TMA-containing system (left) with the Cs+-containing system (right).

The calculations show that the presence of the counterion clearly affects the transport properties of the POM junctions, as evidenced by the difference in the transmission spectra compared to configuration 1 (“no counterions”). The counterions lower the energy of the transmission peaks, bringing them closer to EF. The peaks remain narrow, suggesting that there are still only a few energy states involved in the conduction of electrons. There are significantly more peaks between 1.5 and 2 eV, likely caused due to interaction of counterions with the electrodes. However, given the high bias required to activate these channels, this influence is largely unimportant.

To further understand how counterions influence the transport, the transmission peaks closest to the Fermi energy (0.00 eV) in the range 0–0.51 were analyzed by calculating and visualizing the dominant transmission eigenstates. These are shown in Figure b. First, the symmetry and location of the eigenstates for both systems are identical. There are slight differences in the size of the electron clouds when visualized at the same isovalue—most notably, at 0.135 eV, where the orbitals are much smaller for the Cs+ system. This is consistent with the lower value of T(E). Because the transmission eigenstates for the peaks are the same for both species of counterion, it can be concluded that they influence the transport of POM junctions the same way.

The calculations also predict that all the dominant transmission eigenstates for the analyzed peaks are localized exclusively on the POM molecule. They are all LUMO-like in character, with delocalization over parts of the tungsten atoms in the cage. This implies that the presence of counterions does not create new conductance channels but makes existing ones more energetically stable, which brings them close to the Fermi level and as a result contribute significantly to transport. This reduction in energy makes the transmission modes more likely to be in the bias window and hence increasing the conductance of the system.

Along with the results of configurations 1–3 it highlights the fact the main feasible conductance pathways are focused on the tungsten d-orbitals, clearly like the LUMO level which is delocalized over the whole tungsten cage. These conductance pathways are influenced by the contact to the electrodes, which shifts the density toward other tungsten areas. In addition, the presence of counterions influences the energy levels of the molecule with respect to the Fermi energy of the electrodes. As a result, the transmission channels are affected and moved closer to the Fermi level. The delocalization of the transmission eigenstates may offer a way to tune the transmission pathways by chemically influencing these tungsten atoms or by changing counterions.

Figure shows the computed I–V characteristics for the counterion-containing systems. Because of the increased size of the systems, the computational cost for calculating the transport under an applied bias is much higher. As a result, the calculated bias range is smaller than for the non-counterion-containing system. Nonetheless, the results clearly reveal a strong influence of the counterions on the transport.

Figure 13

Computed I–V curves for the counterion-containing system.

The first thing to note from Figure is that the current for counterion-containing systems rises much more quickly than for configuration 1. Additionally, the magnitude of the calculated current is much higher at lower bias than for configuration 1. Without counterions, at a bias of −1 V the current reaches ∼400 nA. However, for the TMA system, the current reaches this value at a bias of −0.3 V and is more than double by −0.5 V. This is a clear indication of an increase in conductance and therefore current flow for the TMA-containing system. Similarly, for the Cs+ system, the current rises more quickly than in the absence of counterions. By a bias of 0.26 V it is already at around 100 nA, whereas a bias of nearly 1 V is required to reach the same current with no counterions.

The increased current at lower bias is consistent with the effects of the counterion on the transmission peaks relative to the Fermi energy. As these conductance channels have been pushed down in energy, a smaller applied bias is required for them to be within the energy window caused by an applied bias.

Given the almost identical transmission spectra for Cs and TMA system shown in Figure a, the slight differences in the IV curve for the Cs system compared to the TMA system are surprising. Clearly, in positive bias the current rises more slowly in the presence of Cs compared to TMA. It is possible that this is due to the differences in the value of the T(E) in the peaks between 0 and 0.5 eV. From Figure a the TMA spectra show peaks with higher transmission in this region.

The results suggest that the organic cation (TMA) provides better charge transport than the inorganic one (Cs). Although surprising given the identical transmission curves, when considering the relative size/polarization of the two cations studied, it is not surprising that TMA would provide better charge transfer. Further experimental work would be key in confirming this prediction.

Overall, the DFT and NEGF simulations illustrate that the effect of the countercations on the transport properties of POMs is not negligible. In fact, the presence of counterions pushes the unoccupied energy levels of the POM closer to the Fermi energy. As a result, the transmission peaks are also found closer to the Fermi energy. Therefore, a lower bias is required for significant current flow in the POM molecular junction.

Conclusions

We have thoroughly explored the fundamental transport properties of [W18O54(SO3)2]4–, a POM molecule, using DFT and NEGF methods. Our study has provided insight into what factors will influence the transport of this molecule and by extension other POM molecules.

The calculations demonstrate that the electronic structure of the molecule dictates the transport profile of the molecular system. However, the transport is influenced by some key aspects of the device setup. We showed that the contact strength of the molecule to the electrodes greatly influences the current flow through the molecule, with stronger contacts resulting in higher predicted currents. It was also shown that the orientation of the molecule to the electrodes had a key influence on the device’s transmission spectra. The horizontal geometry results in the most favorable transmission profile for this molecule.

Lastly, we illustrated the importance of the counterions on the transport. By exploring two potential counterions, TMA and Cs+, the simulations suggest that the presence of these ions pushes the unoccupied energy levels closer to the Fermi level of the device. Thus, the conductance channels are more readily available at lower biases, and hence the predicted current is much greater than for systems where the counterions are ignored. The simulations suggest that TMA and Cs+ influence the device in much the same way, with no difference in transmission spectra and PLDOS. However, TMA provides better charge transport as illustrated by the computed I–V curves.

This is a purely theoretical and computational work. However, from an experimental point of view we believe that is possible to measure POM-based junction’s I–V curves using the standard experimental techniques in the field. Indeed, this is one of the current efforts for the authors in this work. In addition, future work will focus on further theoretical exploration of the POM molecule on different surface geometries (nonideal structures) and on SiO2 and HfO2 to explore flash memory architectures and further understand the transport of POMs from a theoretical perspective.

We believe that our theoretical exploration provides key insights into the transport behavior of this molecule which will inevitably aid device design for molecular electronic device applications. Namely, the influence of counterions on potential device operation cannot be ignored. Here, it has briefly been studied with one case, but further work could explore the details on how different counterion arrangements, number of ions, and identity of counterion influence the system’s transport and whether they can be taken advantage off from a device engineering point of view.