Silicon – single molecule – silicon circuits

Abstract

In 2020, silicon – molecule – silicon junctions were fabricated and shown to be on average one third as conductive as traditional junctions made using gold electrodes, but in some instances to be even more conductive, and significantly 3 times more extendable and 5 times more mechanically stable. Herein, calculations are performed of single-molecule junction structure and conductivity pertaining to blinking and scanning-tunnelling-microscopy (STM) break junction (STMBJ) experiments performed using chemisorbed 1,6-hexanedithiol linkers. Some strikingly different characteristics are found compared to analogous junctions formed using the metals which, to date, have dominated the field of molecular electronics. In the STMBJ experiment, following retraction of the STM tip after collision with the substrate, unterminated silicon surface dangling bonds are predicted to remain after reaction of the fresh tips with the dithiol solute. These dangling bonds occupy the silicon band gap and are predicted to facilitate extraordinary single-molecule conductivity. Enhanced junction extendibility is attributed to junction flexibility and the translation of adsorbed molecules between silicon dangling bonds. The calculations investigate a range of junction atomic-structural models using density-functional-theory (DFT) calculations of structure, often explored at 300 K using molecular dynamics (MD) simulations. These are aided by DFT calculations of barriers for passivation reactions of the dangling bonds. Thermally averaged conductivities are then evaluated using non-equilibrium Green's function (NEGF) methods. Countless applications through electronics, nanotechnology, photonics, and sensing are envisaged for this technology.

---

Introduction

Recently, synthetic techniques compatible with technologies used in silicon fabrication plants were developed that can assemble molecules with either thiol¹ or disulfide² terminating groups between two silicon-electrode contacts. This can, in principle, pave the way for the inclusion of single-molecules, or else nanoscopic self-assembled monolayers (SAMs), to be imbedded into silicon diodes and transistors. Recently reported conductances for silicon – S(CH₂)₆S – silicon junctions were observed to be, on average, one third of analogous values measured on traditional gold – S(CH₂)₆S – gold junctions, a useful result. More significantly, the silicon devices displaced mechanical stability lasting 5 times longer, with junction extendibility up to 3 times greater,¹ and some single-molecule junctions were found to conduct even better than gold junctions. Unlike with traditional gold electrodes, the alignment of the Fermi energy with the molecular orbitals, and consequently the junction's electric properties, can be dopant controlled. Hence these junctions appear attractive compared to standard approaches used in Molecular Electronics.

A key feature of the synthesis conditions recently developed^1,2 is that they require no heating, pressure, irradiation, or external catalysis. This differs from the wide variety of approaches previously used for grafting molecules onto silicon³ using, e.g.,^4–6 Lewis acids,⁷ Grignard reagents,^8,9 electrografting,¹⁰ and microwave¹¹ or UV-Visible irradiation,^12–18 involving perhaps ultra-high vacuum technologies,^19,20 high-temperature solution chemistry,^21,22 or high-temperature high-pressure processes in supercritical CO₂.²³ Techniques that form oxide-free surfaces have been of central importance,^8,9,23–30 and the newly developed techniques not only avoid oxide but also avoid fabrication-induced SAM damage.

Utilising control over both the properties of the silicon contacts and the bridging molecule,^8,28,30–33 applications can be envisaged to field-effect transistors^25,30,34,35 perhaps for biomedical applications,²⁹ electrochemical applications^36,37 including sensing,³⁸ polymer engineering,³⁹ hydrophobicity,⁴⁰ quantum-dot photonics,⁴¹ photoluminescence,⁴² light harvesting and usage,^43,44 bioimaging, biosensing, and cancer treatment,^45,46 as well as molecular-electronics applications.^{31,36,47–51} Silicon–molecule–metal junctions can also be envisaged and would have useful properties, by analogy to results found for GaAs–molecule–Au junctions.^52–54

Thiol SAMs on gold have dominated previous applications in molecular electronics as they are also easy to prepare and have properties of widespread interest,^32,55–62 but suffer from drawbacks as the Au–S bonding is weak⁶³ and dispersion controlled.^64–66 To make robust devices of atomic dimensions, structural regularity and stability is important, and hence covalent bonding of molecules to silicon offers new technology directions; related covalent-bonding applications involving, e.g., graphene point contacts are also of modern interest.^67–69

To date, two types of silicon – molecule – silicon junctions have been prepared using scanning-tunnelling microscopy (STM) technology: “blinking” junctions, formed by holding a silicon STM tip fixed above a SAM of the molecule pre-prepared on a Si(111)–H substrate, and “break-junction” (STMBJ) junctions formed by crashing a silicon tip into a silicon substrate and then withdrawing the tip. Fig. 1 illustrates the two approaches. The blinking method is adapted from the current–time approach of Nichols and co-workers,⁷⁰ whilst the STMBJ method is adapted from the approach of Tao and co-workers.⁷¹ Such approaches drive many modern applications in Molecular Electronics.⁷²

Fig. 1 Sketches of (a) the blinking experiment, and (b) the STMBJ experiment. Blinking experiments are performed on pre-formed regular SAMs (at ca. 75% coverage), bringing up an STM tip that is assumed to be flat on the atomic scale, detecting single-molecule junctions forming and breaking. STMBJ experiments involving crashing an STM tip into the substrate, and withdrawing it until first single-atom-wide silicon–silicon junctions form and then further until single molecules from solution bridge the formed gap.

The blinking junctions (Fig. 1a) are so-called as the connection between the bridging molecules and one of the electrodes can break and reform on the seconds timescale. The STM tips used in such experiments are regarded as being “sharp”, yet present curvatures of the order of 100 nm and hence can be considered as being atomically flat on the length scale of the experiment. Note that mostly one molecule at a time bridges the two silicon electrodes, despite the large surface area of interaction. The blinking experiment is therefore one in which the atomic structures are controlled: the SAM is ordered, and the two silicon surfaces are ordered. Experimental conditions can be varied to control factors such as SAM coverage and structure.

On the other hand, understanding the atomic structure of STMBJs involving silicon (Fig. 1b) demands solution of a significant number of issues. When the tip is crashed into the substrate, damage to both structures over multiple atomic layers is expected. The original Si(111)–H surface termination will be destroyed by the impact. Covalent bonds will form between the atoms originally in the tip and those originally in the substrate, and extraction of the tip will result in the fracturing of enough covalent bonds to allow separate structures to reform. This will expose silicon atoms on both tip and substrate with “dangling bonds” associated with the collision and separation processes. As these experiments are performed in solution, dangling bonds will then react with re-entering solvent and/or solute, plus any contaminant molecules that it may bring. Relevant to this work, STMBJ conditions involve 1,3,5-trichlorobenzene as the solvent, with the 1,6-hexanedithiol reactant present at 4 μM concentration.¹ The timescale of individual STMBJ measurements is of the order of ms, and so to modify outcome, such reactions must be completed by then. Reported experimentally is conductance information averaged over thousands of such collisions.

This work concerns, firstly, the elucidation of the atomic structure of Si – S(CH₂)₆S – Si junctions formed under either blinking or STMBJ conditions. Some basic interfacial structural models are postulated and optimised using density-functional theory (DFT), mostly utilising molecular dynamics (MD) simulations at 300 K.⁷³ Structures considered include: simple SAMs on regular silicon-surface structures with or without hydrogen termination; fractured silicon structures after STMBJ tip retraction; transition-state structures for reaction of dangling bonds with solvent, solute, and possible contaminants; and SAMs reformed to fractured silicon structures after tip retraction. These structures are obtained as a function of gradual retraction of the STM tip. They are then tested using DFT-based non-equilibrium Green's function (NEGF) calculations⁷⁴ of the electric properties of the electronic device, partitioned as semi-infinite doped silicon – junction region – semi-infinite doped silicon, leading to the evaluation of the electrical conductance. How this conductance varies as a function of tip retraction is then compared with observed conductances,¹ allowing the principles that lead to the observed exceptional device characteristics to be determined. Such knowledge of atomic structure and dynamics, electronic structure, and conductivity is essential for the design and interpretation of almost all experiments performed on applications systems and devices. NEGF calculations are performed for the situations in which the silicon tip and substrate are P-doped (1.15 × 10²⁰ atoms cm⁻³) and N-doped (7 × 10¹⁹ atoms cm⁻³) – these doping levels correspond to the same silicon bulk resistivity of 0.001 Ω cm.⁷⁵

Results

(a) Structure and conductivity of the regular SAMs examined in blinking experiments

The atomic model used for the blinking experiments involves a single S(CH₂)₆S group bridging regular Si(111)–H surfaces representing the substrate (lower surface) and tip (upper surface). Owing to the large curvature of the STM tip (Fig. 1), the upper electrode is assumed flat on the atomic length scale, and the same surface model is therefore used for both surfaces. The thiol hydrogen atoms are regarded as being lost by the chemisorption process¹ that binds the solvent molecule to the two surfaces. Structures are obtained by taking an already-prepared¹ low-coverage (1 : 9) SAM of S(CH₂)₆SH on (3 × 3) Si(111)–H, removing the tail thiol H atom, and optimising the coordinates obtained when a second Si(111)–H surface is brought up that contains one silicon dangling bond. The optimised structure is shown in Fig. 2a; full structural and computational details for this, and indeed all structures reported in figures, are provided in the ESI.†

Fig. 2 Modelling the blinking experiment. (a) 2D (3 × 3) model of two flat Si(111)–H surfaces spanned by S(CH₂)₆S at 1 : 9 coverage; bulk silicon is emulated by H-terminating the outer layer of each silicon slab, with this and the next inner layer frozen at bulk coordinates. The vertical height inside the junction is taken to be the distance z from the uppermost frozen atom in the lower electrode. The junction length is taken to be the vertical extent, d = 24.6 Å, of the region containing optimised atoms. Si-brown, S-yellow, C-cyan, H-white. (b) PDOS (in (eV Å)⁻¹) from NEGF calculations of the conductivity, showing the relative density of electronic states (colour coded) as a function of the vertical height z inside the junction and the electron energy E. See also ESI† Fig. S1.

NEGF calculations of the conductance of this blinking-experiment model yields transmission curves (see Fig. S1 (ESI†)) and projected densities of states (PDOS, Fig. 2b) for the situations in which the silicon tip and substrate are P-doped and N-doped. The PDOS tells the electronic state density at height z above the start of the geometrically optimised region of the junction and at electronic energy E away from the Fermi energy. Fig. 2b shows the PDOS for P-doped silicon, for which the Fermi energy is at the top of the valence band. Within the molecular region, features in the PDOS are seen at z = 8.3 Å and 16.3 Å arising from the sulfur atoms. Between them are three bands each attributable to two carbons. This electronic structure can interact with the valence band of the silicon and portrays a significant conduction pathway that is close to the Fermi energy, supporting hole conductivity in the P-doped junction. As Fig. 2b implies and ESI† Fig. S1 details, no such analogous pathway for electron conduction is apparent, implying that the junction conductance will be much greater for P-type Si than for N-type. Indeed, our NEGF calculations predict a ratio of 10 : 1 for the conductance in P-type silicon compared to N-type.

Fig. 2b depicts the electronic properties of the electrodes and bridging molecule that are critical to the ability to perform realistic conductance predictions. However, as one sees from the figure (and also ESI† Fig. S1), the native silicon bandgap is predicted to be 0.6 eV, significantly removed from the actual value⁷⁶ of 1.17 eV. The calculations could by enhanced by using either improved density functionals,⁷⁷ or GW-based NEGF theory^78,79 to open up the silicon bandgap; however, these calculations are hardly feasible for the size and number of geometries considered herein. Subtle but sometimes significant effects⁸⁰ associated with the stabilisation of (partial) molecular charges at the interface are also not included in these calculations. The difficulty in making quantitatively reliable predictions of junction conductance are likely to be enhanced when irregular and unterminated silicon contacts are considered in subsequent results subsections.

To date, conductance for this junction has only been measured¹ for P-type Si, with a summary given in Table 1. Those experiments were performed using boron P-doped silicon with a bulk resistivity of 0.001 Ω cm, which corresponds to the dopant concentration used in the modelling. The observed conductance at zero voltage is on average G = 70 μG₀, where G₀ = 2e²/ℏ = 77.5 μS is the quantum of conductance. At half maximum, the conductance range HM₋–HM₊ is 65–85 μG₀, a narrow range with ratio η = HM₊/HM₋ = 1.3. In comparison, the calculated conductance is 4 μG₀; an order of magnitude difference between observed and calculated conductances is indicative of reasonable qualitative agreement.

Table 1 Comparison of observed^a and calculated conductances (μG₀) at zero voltage for P-type silicon and for gold junctions bridged by S(CH₂)₆S

Junction	Calc.	Peak	HM−	HM+	η
a Peak, values at half the maximum on the high (HM+) and low (HM−) sides, and their ratio η = HM+/HM−, extracted from the probability distribution^1 of log(G/G0) so that the peak corresponds to the geometric mean of the observed conductances. b Geometric mean from Fig. 6e at 300 K. c Long tail, perhaps extending beyond the detection limit of 1000 μG0. d At T = 0 K.
Si blinking	4	70	65	85	1.3
Si STMBJ	28^b	60	18	250^c	14.3
Au blinking	230 (ref. 1)	180	150	210	1.5
Au STMBJ	230^d	200	90	350	3.9

---

A gauge of the usefulness of the calculated result is its comparison with analogous data pertaining to HS(CH₂)₆SH tethered between gold electrodes. The calculated conductance at zero voltage is 230 μG₀, whilst the observed¹ conductance is similar, displaying a peak at 180 μG₀ and a range of 150–210 μG₀ (η = 1.5). Of note, the narrower observed conductance range for silicon junctions (η = 1.3) implies that the blinking silicon junctions are more reproducible and more stable than the analogous gold junctions. Concerning the calculations, increased errors are anticipated for the modelling of silicon junctions compared to metal junctions as results will depend significantly on how doping, band-gap, and other effects are handled within the DFT and NEGF calculations.

(b) A crude atomic model of an STMBJ: small H-terminated tips on regular Si(111)–H surfaces

Little is known about the atomic structures prepared after a silicon STM tip is crashed into a silicon substrate and then withdrawn, as occurs in an STMBJ experiment. Nevertheless, the unique properties of silicon STM tips after they are crashed into a surface and retracted have long been exploited to make high-resolution STM images,⁸¹ with retracted tips shown to retain function on the day's timescale. It is likely that the differences between standard silicon STM tips and such “super tips” are captured in the caricatures of the blinking and STMBJ experiments presented in Fig. 1. For silicon tips produced by surface collisions, allowed to equilibrate and then used in subsequent experiments, the assumption of small apex-shaped H-terminated tips has led to useful modelling of high-resolution STM images,⁸² and hence this model is taken as a starting point.

In Fig. 3 is shown the atomic model used to depict this scenario. It consists of regular two-layer tips (9 and 4 Si atoms in each layer) that are fully H-terminated and sit on top of Si(111)–H surfaces. The horizontal alignment of the top and bottom electrodes is adjusted to minimise the energy. This is done by first optimising the coordinates of a constrained atomic-cluster model,⁸³ then mapping out a potential-energy surface as a function of the separation between the tips, see ESI† Fig. S2. Overall, this separation is varied over 13 Å, manifesting junction structures that present at least five distinct conformations of the alkane chain. Eight representative structures from this cluster model were then inserted into 2D periodic models of the junction and re-optimised, freezing the outside two layers of each silicon electrode at 3^rd-layer (from the surface) to 3^rd-layer distance d, see Fig. 3 and ESI† Fig. S3.

Fig. 3 Modelling the STMBJ experiment assuming regular H-terminated silicon tips form on each electrode. (a) 2D (3 × 3) model of two flat Si(111)–H surfaces with two added H-terminated Si layers in the shape of a tip, spanned by S(CH₂)₆S, at d = 36.6 Å; Si-brown, S-yellow, C-cyan, H-white. (b) Conductivity at zero voltage evaluated along a 800 fs MD trajectory for P-type silicon. (c) Average energy of each MD simulation. (d) Conductivity at zero voltage for P-type and N-type silicon for static geometries (T = 0 K) and averaged over T = 300 K MD trajectories; the shaded region shows the observed conductance found throughout extensions of 3–10 Å for P-type Si at 300 K. See ESI† Fig. S2 and S3 for more information.

A primary result apparent in Fig. 3 is that, through alkane conformational changes and bond-extension possibilities, these STMBJs can form over 13 Å length extension. This is consistent with key experimental results¹ depicting junctions stable over typically 3–10 Å. Contrary to experiment, however, the calculated conductances are shown in Fig. 3c and depict four orders of magnitude variation in the conductance. Observed¹ conductance histograms (of log(G/G₀)) depict circuits formed with G = 10–1000 μG₀, with a conductance maximum (approximate geometric mean) at 60 μG₀ and conductances at half-maximum of 18 and 250 μG₀ (Table 1) giving η = 14. To compare with the logarithmic scale reported in Fig. 3 and subsequent figures, the observed full conductance range is −5.0 < log₁₀(G/G₀) < −3.0, and representative shaded regions are drawn on the figures.

In summary, the calculated conductances for the regular H-terminated tip model are too low and too variable to explain the observed conductance plateaus. It also appears unlikely that STMBJ structures of this type on regular structures could account for the very wide range of conductances observed in the experiment (η = 14, Table 1) and its sharp contrast to that observed in blinking experiments (η = 1.3). In traditional STMBJ experiments using gold electrodes, it is customary to expect greater conductance variations in STMBJ experiments, e.g. for HS(CH₂)₆SH, η = 1.5 in STMBJ experiments compared to 1.3 in blinking ones (Table 1),¹ but the effect using silicon is dramatically magnified, increasing from 4 to 14.

Nevertheless, insight gained by consideration of the calculated conductances for the regular H-terminated junction provides important information needed to understand more complicated junction models. ESI† Fig. S3b shows transmission as a function of extension for thus model. The conductivity variation shown in Fig. 3d throughout 27 Å ≤ d ≤ 35 Å arises from direct through-space tunnelling between the silicon electrodes (see also ESI† Fig. S5); it takes a similar form if the bridging molecule is removed and the conductance calculations simply repeated. Within the range 35 Å ≤ d ≤ 40 Å, the conductance increases owing to alkane-chain conformational straightening. The conductance is predicted to increase significantly as gauche-type conformations are eliminated at larger junction extensions.⁷³ To interpret the small observed variations in conductance over extensions of 3–10 Å, it is apparent that only straight-chain conformations, or conformations with at most one gauche linkage, can be involved. Hence structural variations other than chain conformational changes must be responsible for a substantial part of the observed effect.

The relatively high conductances indicated in Fig. 3d at d = 40 Å arise as the molecule has become highly stretched. This reduces the molecular band gap and hence increases the conductivity. In practice, such effects may not be observable as stretched structures have a high energy cost of production (Fig. 3c), and hence are likely to embody unsustainable force magnitudes. Indeed, the average force between the last two data points in Fig. 3c is 0.7 eV Å⁻¹, larger than the value of ∼0.5 eV Å⁻¹ observed⁸⁴ for physisorbed sulfur–gold junctions. Hence, of the structures considered, the ones that are most indicative of through-molecule conductivity through linear (or near linear) chains are those in Fig. 3d in the range of 35–37.4 Å. These are 100 times smaller than those calculated for the blinking experimental configuration. Added fully H-terminated tips therefore act as insulators that block conductivity from the silicon through the molecule.

In general,⁸⁵ molecules can show considerable enhancements of conductance between metal electrodes using STMBJ experiments compared to blinking experiments, although much small effects are common,^86–89 and specifically the effect for HS(CH₂)₆SH between gold electrodes is quite small (Table 1). Complementary calculations for junctions with Au electrodes were performed, with the conductance calculated for HS(CH₂)₆SH between flat Au surfaces being G = 227 μG₀, compared to 229 μG₀ when 4-atom tips are inserted on each side of the junction (Table 1). The relative insensitivity to interface structure can be attributable to the ability of gold SAMs to adapt to their circumstances,^90–92 although for some properties such as inelastic electron-tunnelling spectroscopy (IETS) structural details can be critical.⁹³ That the calculations predict 100 times reduction in conductance for regular H-terminated silicon STMBJ tips therefore presents strong evidence suggesting that such a chemical model is inappropriate, a conclusion supported by the observed large conductance range (e.g., η = 14, Table 1).

One effect of note is that the experiments are performed at room temperature, whereas the junction structural models so far reported pertain to 0 K and do not include thermal fluctuations of the interface. The eight representative structures shown Fig. S3a† (and Fig. 3a) are actually those following 1 ps of MD at 300 K. Shown in Fig. 3b is the conductance calculated for the (perhaps most important) trajectory at d = 37.4 Å. The predicted variation in conductance attributable to thermal fluctuations is small compared to other effects that have been noted. As a function of junction extension, the calculated conductances at 0 K and 300 K are very similar, see Fig. 3d. These results are consistent with the identification of the structure of “super tips” used in high-resolution STM imaging as having this basic structural form.^81,82

(c) A crude atomic model of an STMBJ: small bare tips on regular Si(111) bare surfaces

A variant crude structural model was then used in which the bare silicon tips produced after breaking of the contacts in an STMBJ experiment (Fig. 1) are assumed not to react extensively with solvent or solute on the ms timescale of the experiment, leaving all dangling bonds bare. The simplest model for this is to take the H-terminated surface and tips from Fig. 3, remove the terminating hydrogens, and re-optimise the structures. This was followed by MD simulations for 1 ps at 300 K, an operation that induced significant structural changes to the interface, see Fig. 4a for d = 37.4 Å; results at other tip extensions are shown in ESI† Fig. S4. As shown in Fig. 4b, conductances at most junction separations are significantly enhanced compared to the analogous fully H-terminated model in Fig. 3, often now encroaching on the region in which conductance is observed to be sustained over 3–10 Å extensions. Nevertheless, the conductance is predicted to change erratically over 5 orders of magnitude during tip extension, differing from experiment¹ not only in the range of the variation but also in that individual conductance traces are usually observed with at most factors of 3 variation.

Fig. 4 Modelling the STMBJ experiment, using regular unterminated silicon tips on each unterminated electrode. (a) 2D (3 × 3) model for tips spanned by S(CH₂)₆S at d = 37.37 Å; Si-brown, S-yellow, C-cyan, H-white. (b) Conductivity at zero voltage for P-type and N-type silicon at 0 K; the shaded region shows the observed conductance found throughout extensions of 3–10 Å for P-type Si. See ESI† Fig. S4 for more information. (c) Comparison of the transmission for P-type Si, as a function of electron energy from the Fermi energy, at d = 32.17 Å, from this unterminated series to that for the analogous terminated structure (Fig. 3, ESI† Fig. S3).

Fig. 4c demonstrates the cause of this erratic behaviour, the major qualitative difference found between calculations using regular tips with and without H-termination. Shown is the electronic transmission probability as a function of energy away from the Fermi energy for P-type silicon at d = 32.2 Å, evaluated for both H-terminated tips and surfaces and for their unterminated analogues. For the unterminated tip, silicon dangling-bond states appear within the silicon bandgap at seemingly random intervals. If perchance one of these states appears close to the Fermi energy and the molecular energies, then it will dramatically enhance conductivity, giving rise to the profound differences in conductivity perceived between Fig. 3d and 4b, as well as the erratic tip-extension dependence shown in Fig. 4b.

The variation in the electronic structure associated with dangling bonds will be very sensitive to the details of both the silicon geometry and the silicon–sulfur interface. Hence the conductance is expected to vary significantly during MD simulations. Realistic calculations of junction conductance therefore require understanding of what actual tip shapes could be, how these can react with solvent, solute, etc. on the experimental timescale, and understanding of the thermal effects that are likely to modulate them. Needless to say, quantitatively accurate calculations will also require a much better description of the silicon band gap than the current, PBE-based, calculations provide.

Inspiration concerning the basic nature of the problems encountered in understanding the nuclear and electronic structure of unterminated Si(111) can be found in the literature of the bare Si(111) surface itself. For this, following decades of research, controversy remains concerning the basic surface nuclear and electronic structure, as well as the determination of appropriate computational methods for its simulation.^94,95 The central issue pertains to how silicon dangling bonds interact with their neighbours, concerning both the alleviation of the intrinsic chemical instability that arises and the resulting electronic structure and its influence on conductivity. Unless full passivation of the STMBJ tips occur after the solution re-enters the broken junction, such issues will also arise. As full passivation acts as an electrical insulator (Fig. 3), it is unlikely to have been observed in silicon STMBJ experiments.

(d) Structures made using MD to fracture a fused substrate-tip structure

In search of a more realistic silicon atomic topology, MD simulations were performed, at 300 K, in which an initial structure representing the STM tip crashed into the substrate was taken and pulled to breaking point. The initial structure, somewhat unrealistically, was taken as an 11-layer Si(111) slab with random atoms removed from the central layer. An odd number of total atoms was chosen to assure asymmetry in the products obtained. Steps during the fracturing process, as well as in its compression to the bulk silicon density, are illustrated in ESI †Fig. S5. These MD simulations do not take into account any large-scale tip deformations that could occur on the experimental timescale⁸⁵ as they explore only the chemical processes occurring at the junction interface. The conductance as a function of the fracturing is also shown in ESI† Fig. S5.

The two irregular tips so produced, one with a clear apex atom, the other with 3 atoms in its uppermost plane, were not H-terminated and a straight-chained single bridging unit S(CH₂)₆S was added to connect to the top and bottom electrodes. MD calculations were then run for 1 ps, and the structures were slowly pulled apart and the dynamics repeated for each step (Fig. 5b). At the shortest separation investigated, a kink resulted in the alkane chain, but otherwise a linear chain persisted over an 8 Å retraction of the top electrode. To facilitate this, a binding location change occurred, with the S atom binding to an Si atom below the top of the tip at short distances, and to the top at larger ones.

Fig. 5 Modelling the STMBJ experiment, following MD at 300 K on regular unterminated silicon tips made by simulating a tip-surface crash and withdrawal (see ESI† Fig. S5). (a) Conductivity at zero voltage for P-type and N-type silicon at 300 K; the shaded region shows the observed conductance found throughout extensions of 3–10 Å for P-type Si. (b) 2D (3 × 3) model for tips spanned by S(CH₂)₆S; Si-brown, S-yellow, C-cyan, H-white. See ESI† Fig. S6 for more information.

Fig. 5a shows that this model predicts conductances that strongly encroach on the observed range, sustained for extensions of up to 8 Å as also observed. These were averaged over 1 ps of MD and do not show the same erratic behaviour found in Fig. 4 as the accidental resonances between dangling-bond levels and molecular levels that drive conduction are averaged over. From a quantitative perspective, the thermal fluctuations in the conductance are so large that the performed MD averaging is insufficient, and hence order-of-magnitude error bars may still be anticipated for the thermally averaged results.

(e) Likely chemical reactions on the fractured structures

As dangling bonds dominate conduction properties, it would seem essential that the best-possible descriptions of the nuclear and electronic structures of STMBJ junctions be obtained. This demands understanding how dangling bonds become passivated on the ms STMBJ timescale.

In the experiment being modelled, exposed silicon dangling bonds would react rapidly with any dissolved O₂ in the solution to form surface oxides. The facile synthesis used to make the junctions starting from thiol reactants relies on the presence of ambient levels of dissolved O₂, which effectively provides a catalyst for the process.¹ In this way, all dissolved O₂ is believed to be consumed, and indeed no trace of surface oxide can be found after SAMs were left standing. Hence the available experimental evidence indicates that oxide formation by this mechanism is not a likely process.

An alternative route to surface oxides is through reactions of dangling bonds with contaminate water in the 1,3,5-trichlorobenzene solvent. The solvent had been extensively dried,¹ but the concentration of the dithiol solute is only 4 μM and so water may still present possible reactive pathways. According to standard bond-energy tables, the reaction

2Si˙ + H₂O → SiOH + SiH (1)

should be exothermic by 3.1 eV. A MD simulation was run in which the gap between the two unterminated electrodes was widened and liquid water inserted, see ESI† Fig. S7. Water reacted with the dangling bonds within 1 ps. Hence the water concentration is an important parameter pertaining to the chemical stabilisation of the bare tips.

The solute molecules HS(CH₂)₆SH clearly do react with the silicon tips on the STMBJ timescale as the observed conductance signal stems from the resultant bridged-electrode structure. Calculations¹ predict that thiols react barrierlessly with silicon surface radicals, hence facilitating sub-ps reaction times akin to that calculated above for the reaction with water. At 4 μM concentration, the reaction time would be expected to be in the ns–μs range. Diffusion of solvent and solute into the nano-cavity formed by tip retraction is therefore likely to be a critical aspect in determining the actual reaction kinetics. The relative importance of water reacting instead of the dithiol would scale with the relative concentrations of the water contaminant to the dithiol solute. If the reaction with water was dominant, then break junctions would be difficult to form. Hence the experiment suggest that reactions with water are not significant, and we neglect this possibility henceforth.

A reaction of surface dangling bonds was noted in some MD runs when the starting geometry was at high energy. This involved the hydrogen abstraction reaction

2Si˙ + HS(CH₂)₆SH → 2SiH + HS(CH=CH)(CH₂)₄SH, (2)

but standard bond energies suggest this reaction is endothermic and therefore not of great concern.

Reactions with the 1,3,5-trichlorobenzene solvent provide another possible mechanism to stabilise surface silicon dangling bonds. Some products obtained computationally by bringing a solvent molecule up to either the top or bottom tip fragments from ESI† Fig. S5 are shown in ESI† Fig. S8. A reaction to break a C=C π bond in the benzene and add two Si–C bonds in its place was calculated to be exothermic, with ΔE = −0.39 eV, whilst the insertion of a Si biradical inside a C–Cl bond was calculated to be more so, with ΔE = −1.70 eV. A subsequent H-shift onto a neighbouring Si radical lowered the energy further to −2.88 eV. Analogous insertions of a silicon biradical into a solvent CH bond were predicted to be endothermic. All reactions considered were found to require activation, with the lowest-energy transition-state obtained being that for the Si biradical insertion reaction at ΔE^† = 0.49 eV. To be competitive against a barrierless reaction with the solute at 4 μM concentration, the reaction barrier would need to be less than 0.3 eV.

Of all the reaction mechanisms considered, the most likely one under the experimental conditions is therefore reaction of silicon radicals with the dithiol solute molecules. This reaction is expected to be fast on the ms timescale of the experiments, and is known to occur to some extent as only the products of this reaction can give rise to the observed conductivity. Whilst the calculations do not rule out the possibility of reactions of silicon dangling bonds with the solvent, this possibility is neglected in the subsequent modelling.

(f) Molecular dynamics to model reformed SAMs

The calculations suggest that the bare silicon tips formed after fracturing of a crashed tip and substrate junction will spontaneously react with the dithiol solute molecules to reform SAMs on both the top and bottom electrodes. Three possible SAMs attached to the fractured tips from Fig. S5† were manually constructed and named “l”, “m”, and “n”. This construction was done at a short electrode spacing such that the SAMs attached to the top and bottom electrodes overlapped, delivering a combined SAM coverage of 6 : 9, slightly less than the value of 75% ± 8% observed¹ for dense regular S(CH₂)₆S SAMs on flat Si(111)–H surfaces. Results following optimisation at 0 K are shown in see Fig. 6 and ESI† Fig. S9a–c, including calculated conductances in both N-type and P-type silicon. The structures were then stretched and the geometry optimisations repeated. The total extension applied was 6.5 Å, encompassing some alkane-chain conformational changes at short distances and junction rapture at the largest separation used, d = 41.2 Å. The calculated conductances shown in Fig. 6b and c show regions of regular conductance extending over 4 Å of stretching, with magnitudes similar to those observed. Combined, the three series motivate how seemingly similar chemical structures can sustain conductance over the lower-region of the range observed, with conductance values spreading over the wide observed range.

Fig. 6 Modelling the STMBJ experiment, following MD at 300 K on SAMs made by simulating a tip-surface crash and withdrawal and then chemisorption of HS(CH₂)₆SH solute molecules, with one molecule bridging top and bottom. (a) 2D (3 × 3) model for tips spanned by S(CH₂)₆S; Si-brown, S-yellow, C-cyan, H-white. (b) Conductivity variations at zero voltage for P-type and N-type silicon at 300 K along a MD trajectory. (c) Conductivity at zero voltage for P-type silicon at 0 K. (d) Conductivity at zero voltage for N-type silicon at optimised geometries (T = 0 K). (e) Conductivity at zero voltage for P-type and N-type silicon averaged over T = 300 K MD trajectories. The shaded region shows the observed conductance found throughout extensions of 3–10 Å for P-type Si. See ESI† Fig. S5 and S9–S12 for more information.

The “n” series was then chosen and MD simulations run for 1 ps. Sample structures obtained are highlighted in ESI† Fig. S10; at large distances, fracturing of the junction, at least for small times, are predicted, sometimes associated with Si–S bond breakage and sometimes with Si–Si cleavage. Sometimes, dynamics resulted in chemical reactions between the breaking fragments and the surrounding SAM, forming S–S or Si–S bonds, indicating a new mechanism for chemical reactions occurring inside molecular junctions.^89,96,97 Averaged conductances are reported in Fig. 6e as a function of extension and overall in Table 1 (900 μG₀), many averaged junction properties are reported in ESI† Fig. S12, and individual conductance time histories are reported in ESI† Fig. S11 and Fig. 6b. The variations in conductance with distance averaged over the MD trajectories appears enhanced compared to those reported in Fig. 6c and d for T = 0, though the sampling performed may not be sufficient to justify this conclusion. Indeed, the calculated time dependence of the conductivity is large (Fig. 6b, ESI† Fig. S10), as found for previous MD simulations. The number of passivated dangling bonds per cell in these junctions is 14, whereas the number of remnant dangling bonds is 18, so only 43% are passivated, and the unpassivated bonds continue to exert a dominant influence over the conductivity.

Note that none of the investigated structures could facilitate slippage of the Si–S contact from an inner Si row to an outer row, as was found in Fig. 5 to significantly enhance the extendibility junctions beyond the ca. 4 Å extent perceived in Fig. 6. Another limitation of the current calculations is that they do not allow for extended surface passivation by SAM growth as the tips are pulled apart and the overlap between the SAMs bound to the top and bottom electrodes diminishes. Such SAM growth would act to reduce the coverage of dangling bonds, but the coverage would presumably remain sufficient to maintain the predicted high conductances.

Conclusions

The calculations presented herein confirm the observed relatively high average conductivity of the silicon junctions, explain how it can be that some silicon junctions can be more conductive than typical gold ones, and reveal the mechanisms through which the increased extendibility arises: a combination of Si–S bond-slippage effects, as well as Si-tip reconstruction and elongation, rather than alkyl-chain conformational flexibility. Overall, the calculations motivate how it is that STMBJ junctions can display conductances varying over an extremely wide range, and that conductance can be sustained over unprecedented mechanical extensions.

Whereas specific quantitative predictions arising from the calculations performed are likely to be wayward owing to the many simplifications used, the overall qualitative scenarios concerning junction performance and properties are expected to be robust. Calculation improvements could involve using a more accurate method than PBE-D3(BJ) to model structure and reactivity, significantly increased MD simulation times, and use of improved NEGF techniques. In particular, the PBE method used in the NEGF calculations underestimates bandgaps and poorly treats image charges;⁸⁰ the use of empirical modifications or more advanced and expensive methods^77,80,98 could be used to reduce this effect, if warranted.

A significant difference found between the properties of metal – molecule – metal junctions and the silicon – S(CH₂)₆S – silicon ones considered herein is that the conductivity mechanism appropriate to STMBJ experiments differs from that appropriate to blinking experiments. The different experiments realise different junction atomic structures, a feature that produces an often small effect on conductance, whereas the effects on silicon junctions are profound.

The silicon STMBJ junction model most similar to the analogous blinking-experiment model involves fully H-terminated tips, but this is predicted to be insulating rather than conductive. The conductance of STMBJ junctions could only reasonably be modelled when chemical models included silicon junctions with surface dangling bonds. Their conductance is predicted to be very sensitive to the location of the dangling-bonds within the silicon band gap, making them also very sensitive to nuclear structure and associated thermal motion. All dangling-bond chemical scenarios considered herein as models for STMBJs led to features consistent with the primary experimental observations, but the details varied considerably. The minimum level of calculation needed to produce qualitatively reasonable results would appear to be simple bare tips made following the fracture of STM-tip and substrate, as considered in Fig. 5, with sophisticated models like the reformed SAMs considered in Fig. 6 being the next level of improvement.

Using metal junctions like gold, the conductance measured using STMBJ experiments shows greater variability than that measured using blinking experiments, but the results are qualitatively similar. That analogous experiments give very different results using silicon, demonstrating STMBJ conductance ranges that are ten times broader than those obtained from blinking experiments, is explained by the calculations in terms of the necessary presence of remnant silicon surface dangling bonds after cleavage of fused silicon units and subsequent reaction with the solute. It is by manipulating the silicon dangling bonds that some junctions can conduct better than do analogous gold junctions. This dangling-bond effect therefore presents a new way by which molecular-electronic devices could be manipulated and controlled. Added to this is the result apparent from the calculations that circuits can be engineered by varying the silicon doping. These opportunities present many ways by which future devices could be designed.

Overall, the properties of silicon – single molecule – silicon circuits are found to be very different to the traditional metal – single molecule – metal circuits. To expand on this, the potential properties of silicon P–N junctions bridged by molecules are of great interest, and there is an immediate need for their construction and characterisation.

Methods

(a) Geometries of model clusters

Initial structures were obtained by taking two H-terminated tip-like clusters Si₁₀H₂₁, each with one silicon dangling bond site, and bridging S(CH₂)₆S between them. The outside rows of Si atoms were frozen parallel at distance d′ (see ESI† Fig. S2) ranging from 6.8 Å to 28.0 Å apart. These structures were then optimized by Gaussian-16 (ref. 99) using the B3LYP density function,¹⁰⁰ the 6-31G* basis set,¹⁰¹ and the D3(BJ) dispersion correction.¹⁰²

(b) Geometries at 0 K of 2D-periodic Si – S(CH₂)₆S – Si junctions

2D-periodic models for Si – molecule – Si junctions embody the molecule attached to two slabs representing bulk silicon. A total of 4 layers of Si(111) atoms were used in each slab, with the outside edges H-terminated to mimic the bulk material. The outer two rows on each slab were frozen in all calculations, with the distance d (Fig. 3–6) between the second-row silicon atoms in the top and bottom slabs reported as a measure of the inter-electrode vertical separation. A (3 × 3) model of the surfaces was used containing 18 silicon atoms per row, with a lattice parameter of 3 × 3.826862 Å (bulk Si–Si bond length 2.34346 Å). At least 15 Å of vacuum region was placed between each vertical image of the system. Treatment of the inner surfaces varied with sample: sometimes atoms were left bare, sometimes H-terminated, sometimes regular Si tips were added to the innermost layer, and sometimes irregular-shaped tips formed by breaking apart a single piece of silicon were used. Also, sometimes only one slab was used and chemical reaction and activation energies calculated for the attack of various solvent or solute molecules on silicon dangling bonds.

2D-geometry optimisations were performed by VASP^103,104 5.4.1 using the PBE density functional¹⁰⁵ with the D3(BJ) dispersion correction¹⁰² applied in the in-layer directions. Control parameters used include: normal precision, an energy convergence criterion of 10⁻⁵ eV, and a gradient convergence of 0.01 eV Å⁻¹. Thermal broadening of the occupied orbital energies was performed at an electronic temperature corresponding to 0.2 eV. Dipole corrections were not employed.

(c) Molecular dynamics studies of 2D-periodic Si – S(CH₂)₆S – Si junctions

Molecular dynamics simulations at 300 K were performed on various junctions by VASP^103,104 5.4.1 using the PBE density functional¹⁰⁵ with the D3(BJ) dispersion correction¹⁰² applied in the in-layer directions. Sample and unit cell details were as per the 2D calculations performed at 0 K (see above). Control parameters used include: low precision, the “very fast” algorithm, an energy convergence criterion of 10⁻⁴ eV, and a time step of 1 fs (results using 2 fs were insignificantly different, but 3 fs proved unacceptable), and a trajectory length of 1 ps. All of the data was used in calculating averaged properties, except that configurations were sampled after every 16 fs for conductivity analysis.

(d) NEGF calculations

Electronic-transport properties were calculated with the NEGF technique,⁷⁴ as implemented in Nanodcal package (version 2016).^74,106 The molecular electronic junction was divided into three regions, bottom electrode, contact region, and top electrode. The contact region includes parts of the physical electrodes (4 layers of silicon or 3 layers of gold), possible additional electrode atoms above the surface, and the bridging molecule. All NEGF calculations of molecular junction electron transport properties were performed self-consistently at T = 300 K electronic temperature. The contact region Green's functions are used to compute electronic transmission spectra and, consequently, molecular junction conductivities. The PBE density functional was used, along with double-zeta with polarization basis set, k-space grids of 3 × 3 × 1, and 50 a.u. energy cutoff for the real-space grid. Pseudopotentials for silicon atoms were fitted separately to reproduce the required experimental phosphorous (N-type) and boron (P-type) doping, utilizing DFT bulk silicon calculations. These are set to mimic doping levels of 1.15 × 10²⁰ atoms cm⁻³ for P-doped silicon, and 7 × 10¹⁹ atoms cm⁻³ for N-doped.

Data availability

The most significant data concerns the atomic coordinates and calculation properties pertaining to the 80 structural images shown in figures in the main text and ESI, for which full details are provides in ESI. Additional information is available from the authors.

Author contributions

J. R. R. and J. Y. performed the DFT calculations, whilst D. S. K. performed the NEGF calculations. J. R. R., N. D., and D. S. K. designed the research.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank National Computational Infrastructure (NCI) Australia for the provision of computer resources, as well as the Shanghai University ICQMS computational facility.

References

1. C. R. Peiris, S. Ciampi, E. M. Dief, J. Zhang, P. J. Canfield, A. P. Le Brun, D. S. Kosov, J. R. Reimers and N. Darwish, Chem. Sci., 2020, 11, 5246.
2. E. Dief, Y. Vogel, C. Peiris, A. Le Brun, V. Gonçales, S. Ciampi, J. Reimers and N. Darwish, Langmuir, 2020, 36, 14999.
3. H. Sugimura, in Nanocrystalline Materials: Their Synthesis-Structure-Property Relationships and Applications, ed. S.-C. Tjong, Elsevier, 2 edn, 2013, p. 161.
4. M. R. Linford and C. E. D. Chidsey, J. Am. Chem. Soc., 1993, 115, 12631.
5. A. B. Sieval, A. L. Demirel, J. W. M. Nissink, M. R. Linford, J. H. Van Der Maas, W. H. De Jeu, H. Zuilhof and E. J. R. Sudhölter, Langmuir, 1998, 14, 1759.
6. J. M. Buriak and M. D. H. Sikder, J. Am. Chem. Soc., 2015, 137, 9730.
7. J. M. Buriak, M. P. Stewart, T. W. Geders, M. J. Allen, H. C. Choi, J. Smith, D. Raftery and L. T. Canham, J. Am. Chem. Soc., 1999, 121, 11491.
8. S. R. Puniredd, O. Assad and H. Haick, J. Am. Chem. Soc., 2008, 130, 13727.
9. E. J. Nemanick, P. T. Hurley, L. J. Webb, D. W. Knapp, D. J. Michalak, B. S. Brunschwig and N. S. Lewis, J. Phys. Chem. B, 2006, 110, 14770.
10. E. G. Robins, M. P. Stewart and J. M. Buriak, Chem. Commun., 1999, 2479.
11. R. Boukherroub, A. Petit, A. Loupy, J. N. Chazalviel and F. Ozanam, J. Phys. Chem. B, 2003, 107, 13459.
12. A. I. A. Soliman, T. Utsunomiya, T. Ichii and H. Sugimura, Langmuir, 2018, 34, 13162.
13. Q. Y. Sun, L. C. P. M. de Smet, B. van Lagen, M. Giesbers, P. C. Thüne, J. van Engelenburg, F. A. de Wolf, H. Zuilhof and E. J. R. Sudhölter, J. Am. Chem. Soc., 2005, 127, 2514.
14. Q. Y. Sun, L. C. P. M. de Smet, B. van Lagen, A. Wright, H. Zuilhof and E. J. R. Sudhölter, Angew. Chem., Int. Ed., 2004, 43, 1352.
15. Y. B. Vogel, L. Zhang, N. Darwish, V. R. Gonçales, A. Le Brun, J. J. Gooding, A. Molina, G. G. Wallace, M. L. Coote and J. Gonzalez, Nat. Commun., 2017, 8, 2066.
16. N. Gergel-Hackett, C. D. Zangmeister, C. A. Hacker, L. J. Richter and C. A. Richter, J. Am. Chem. Soc., 2008, 130, 4259.
17. M. Hu, F. Liu and J. M. Buriak, ACS Appl. Mater. Interfaces, 2016, 8, 11091.
18. Y. S. Huang, C. H. Chen, C. H. Chen and W. H. Hung, ACS Appl. Mater. Interfaces, 2013, 5, 5771.
19. S. K. Coulter, M. P. Schwartz and R. J. Hamers, J. Phys. Chem. B, 2001, 105, 3079.
20. Y. H. Lai, C. T. Yeh, C. C. Yeh and W. H. Hung, J. Phys. Chem. B, 2003, 107, 9351.
21. H. Sano, K. Ohno, T. Ichii, K. Murase and H. Sugimura, Jpn. J. Appl. Phys., 2010, 49, 01.
22. C. A. Hacker, Solid-State Electron., 2010, 54, 1657.
23. B. Bhartia, S. R. Puniredd, S. Jayaraman, C. Gandhimathi, M. Sharma, Y. C. Kuo, C. H. Chen, V. J. Reddy, C. Troadec and M. P. Srinivasan, ACS Appl. Mater. Interfaces, 2016, 8, 24933.
24. H. G. Steinrück, A. Schiener, T. Schindler, J. Will, A. Magerl, O. Konovalov, G. Li Destri, O. H. Seeck, M. Mezger, J. Haddad, M. Deutsch, A. Checco and B. M. Ocko, ACS Nano, 2014, 8, 12676.
25. L. J. Webb, D. J. Michalak, J. S. Biteen, B. S. Brunschwig, A. S. Y. Chan, D. W. Knapp, H. M. Meyer III, E. J. Nemanick, M. C. Traub and N. S. Lewis, J. Phys. Chem. B, 2006, 110, 23450.
26. J. W. McPherson and H. C. Mogul, J. Appl. Phys., 1998, 84, 1513.
27. A. B. Sieval, R. Linke, H. Zuilhof and E. J. R. Sudhölter, Adv. Mater., 2000, 12, 1457.
28. A. Vilan and D. Cahen, Chem. Rev., 2017, 117, 4624.
29. O. Assad and H. Haick, IEEE Int. Symp. Ind. Electron., 2008.
30. G. Ashkenasy, D. Cahen, R. Cohen, A. Shanzer and A. Vilan, Acc. Chem. Res., 2002, 35, 121.
31. A. Vilan, O. Yaffe, A. Biller, A. Salomon, A. Kahn and D. Cahen, Adv. Mater., 2010, 22, 140.
32. A. Vilan, D. Aswal and D. Cahen, Chem. Rev., 2017, 117, 4248.
33. J. M. Buriak, Chem. Rev., 2002, 102, 1271.
34. J. W. McPherson, Microelectron. Reliab., 2012, 52, 1753.
35. K. T. Wong and N. S. Lewis, Acc. Chem. Res., 2014, 47, 3037.
36. O. S. Kang, J. P. Bruce, D. E. Herbert and M. S. Freund, ACS Appl. Mater. Interfaces, 2015, 7, 26959.
37. A. A. Yasseri, D. Syomin, R. S. Loewe, J. S. Lindsey, F. Zaera and D. F. Bocian, J. Am. Chem. Soc., 2004, 126, 15603.
38. J. Wang, Y. Zhou, M. Watkinson, J. Gautrot and S. Krause, Sens. Actuators, B, 2015, 209, 230.
39. M. Dasog, J. Kehrle, B. Rieger and J. G. C. Veinot, Angew. Chem., Int. Ed., 2016, 55, 2322.
40. E. S. Brown, S. Hlynchuk and S. Maldonado, Surf. Sci., 2016, 645, 49.
41. M. Dasog, K. Bader and J. G. C. Veinot, Chem. Mater., 2015, 27, 1153.
42. W. J. I. Debenedetti, S. K. Chiu, C. M. Radlinger, R. J. Ellison, B. A. Manhat, J. Z. Zhang, J. Shi and A. M. Goforth, J. Phys. Chem. C, 2015, 119, 9595.
43. A. Fermi, M. Locritani, G. Di Carlo, M. Pizzotti, S. Caramori, Y. Yu, B. A. Korgel, G. Bergamini and P. Ceroni, Faraday Discuss., 2015, 185, 481.
44. W. Peng, S. M. Rupich, N. Shafiq, Y. N. Gartstein, A. V. Malko and Y. J. Chabal, Chem. Rev., 2015, 115, 12764.
45. F. Peng, Y. Su, Y. Zhong, C. Fan, S. T. Lee and Y. He, Acc. Chem. Res., 2014, 47, 612.
46. X. Cheng, S. B. Lowe, P. J. Reece and J. J. Gooding, Chem. Soc. Rev., 2014, 43, 2680.
47. L. H. Yu, N. Gergel-Hackett, C. D. Zangmeister, C. A. Hacker, C. A. Richter and J. G. Kushmerick, J. Phys.: Condens. Matter, 2008, 20, 374114.
48. S. Ciampi, M. H. Choudhury, S. A. B. A. Ahmad, N. Darwish, A. L. Brun and J. J. Gooding, Electrochim. Acta, 2015, 186, 216.
49. Y. B. Vogel, L. Zhang, N. Darwish, V. R. Gonçales, A. Le Brun, J. J. Gooding, A. Molina, G. G. Wallace, M. L. Coote, J. Gonzalez and S. Ciampi, Nat. Commun., 2017, 8, 2066.
50. A. C. Aragonès, N. Darwish, S. Ciampi, F. Sanz, J. J. Gooding and I. Diéz-Pérez, Nat. Commun., 2017, 8, 15056.
51. A. Allain, J. Kang, K. Banerjee and A. Kis, Nat. Mater., 2015, 14, 1195.
52. A. Vezzoli, R. J. Brooke, N. Ferri, S. J. Higgins, W. Schwarzacher and R. J. Nichols, Nano Lett., 2017, 17, 1109.
53. A. Vezzoli, R. J. Brooke, S. J. Higgins, W. Schwarzacher and R. J. Nichols, Nano Lett., 2017, 17, 6702.
54. A. Vezzoli, R. J. Brooke, N. Ferri, C. Brooke, S. J. Higgins, W. Schwarzacher and R. J. Nichols, Faraday Discuss., 2018, 210, 397.
55. H. Rieley and G. K. Kendall, Langmuir, 1999, 15, 8867.
56. J. C. Love, L. A. Estroff, J. K. Kriebel, R. G. Nuzzo and G. M. Whitesides, Chem. Rev., 2005, 105, 1103.
57. G. Heimel, L. Romaner, E. Zojer and J. L. Brédas, Nano Lett., 2007, 7, 932.
58. J. G. Van Alsten, Langmuir, 1999, 15, 7605.
59. H. Haick and Y. Paz, J. Phys. Chem. B, 2003, 107, 2319.
60. C. Vericat, M. E. Vela, G. Benitez, P. Carro and R. C. Salvarezza, Chem. Soc. Rev., 2010, 39, 1805.
61. N. Darwish, P. K. Eggers, S. Ciampi, Y. Tong, S. Ye, M. N. Paddon-Row and J. J. Gooding, J. Am. Chem. Soc., 2012, 134, 18401.
62. A. H. Pakiari and Z. Jamshidi, J. Phys. Chem. A, 2010, 114, 9212.
63. N. Darwish, P. K. Eggers, S. Ciampi, Y. Zhang, Y. Tong, S. Ye, M. N. Paddon-Row and J. J. Gooding, Electrochem. Commun., 2011, 13, 387.
64. Y. Wang, Q. J. Chi, J. D. Zhang, N. S. Hush, J. R. Reimers and J. Ulstrup, J. Am. Chem. Soc., 2011, 133, 14856.
65. J. Yan, R. Ouyang, P. S. Jensen, E. Ascic, D. A. Tanner, B. Mao, J. Zhang, C. Tang, N. S. Hush, J. Ulstrup and J. R. Reimers, J. Am. Chem. Soc., 2014, 136, 17087.
66. R. Ouyang, J. Yan, P. S. Jensen, E. Ascic, S. Gan, D. Tanner, B. Mao, L. Niu, J. Zhang and C. Tang, ChemPhysChem, 2015, 16, 928.
67. Y. Cao, S. Dong, S. Liu, L. He, L. Gan, X. Yu, M. L. Steigerwald, X. Wu, Z. Liu and X. Guo, Angew. Chem., Int. Ed., 2012, 51, 12228.
68. C. Yang, Z. Liu, Y. Li, S. Zhou, C. Lu, Y. Guo, M. Ramirez, Q. Zhang, Y. Li, Z. Liu, K. N. Houk, D. Zhang and X. Guo, Sci. Adv., 2021, 7, eabf0689.
69. Z. Liu, X. Li, H. Masai, X. Huang, S. Tsuda, J. Terao, J. Yang and X. Guo, Sci. Adv., 2021, 7, eabe4365.
70. W. Haiss, C. Wang, I. Grace, A. S. Batsanov, D. J. Schiffrin, S. J. Higgins, M. R. Bryce, C. J. Lambert and R. J. Nichols, Nat. Mater., 2006, 5, 995.
71. B. Xu and N. J. Tao, Science, 2003, 301, 1221.
72. P. Gehring, J. M. Thijssen and H. S. J. van der Zant, Nat. Rev. Phys., 2019, 1, 381.
73. M. Paulsson, C. Krag, T. Frederiksen and M. Brandbyge, Nano Lett., 2009, 9, 117.
74. J. Taylor, H. Guo and J. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 245407.
75. W. R. Thurber, R. L. Mattis and Y. M. Liu, The relationship between resistivity and dopant density for phosphorous- and boron-doped silicon, US Government Printing Office, Washington, 1981.
76. O. Madelung, Semiconductors: Data Handbook, Springer-Verlag, Berlin, 2004.
77. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
78. K. S. Thygesen and A. Rubio, Phys. Rev. B, 2008, 77, 115333.
79. M. Strange, C. Rostgaard, H. Häkkinen and K. S. Thygesen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 115108.
80. M. Thoss and F. Evers, J. Chem. Phys., 2018, 148, 030901.
81. J. E. Demuth, U. Koehler and R. J. Hamers, J. Microsc., 1988, 152, 299.
82. Ó. Paz, I. Brihuega, J. M. Gómez-Rodríguez and J. M. Soler, Phys. Rev. Lett., 2005, 94, 056103.
83. J. Yang, D. S. Kosov and J. R. Reimers, Shenzhen Daxue Xuebao, Ligongban, 2021 , in press, https://journal.szu.edu.cn/en/default.aspx.
84. M. Frei, S. V. Aradhya, M. S. Hybertsen and L. Venkataraman, J. Am. Chem. Soc., 2012, 134, 4003.
85. F. Evers, R. Korytár, S. Tewari and J. M. van Ruitenbeek, Rev. Mod. Phys., 2020, 92, 035001.
86. W. Haiss, R. J. Nichols, H. van Zalinge, S. J. Higgins, D. Bethell and D. J. Schiffrin, Phys. Chem. Chem. Phys., 2004, 6, 4330.
87. Z. Li and D. S. Kosov, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 035415.
88. I. Díez-Pérez, J. Hihath, Y. Lee, L. Yu, L. Adamska, M. A. Kozhushner, I. I. Oleynik and N. Tao, Nat. Chem., 2009, 1, 635.
89. A. C. Aragonès, N. L. Haworth, N. Darwish, S. Ciampi, N. J. Bloomfield, G. G. Wallace, I. Diez-Perez and M. L. Coote, Nature, 2016, 531, 88.
90. Y. Wang, Q.-J. Chi, J.-D. Zhang, N. S. Hush, J. R. Reimers and J. Ulstrup, J. Am. Chem. Soc., 2011, 133, 14856.
91. R. Ouyang, J. Yan, P. S. Jensen, E. Ascic, S. Gan, D. A. Tanner, B. Mao, L. Ni, J. Zhang, C. Tang, N. S. Hush, J. R. Reimers and J. Ulstrup, ChemPhysChem, 2015, 16, 928.
92. J. R. Reimers, M. J. Ford, S. M. Marcuccio, J. Ulstrup and N. S. Hush, Nat. Rev. Chem., 2017, 1, 0017.
93. G. C. Solomon, A. Gagliardi, A. Pecchia, T. Frauenheim, A. Di Carlo, J. R. Reimers and N. S. Hush, J. Chem. Phys., 2006, 124, 094704.
94. J. E. Demuth, Phys. Status Solidi B, 2020, 257, 2000229.
95. J. E. Demuth, J. Phys. Chem. C, 2020, 124, 22435.
96. N. S. Hush, A. T. Wong, G. B. Bacskay and J. R. Reimers, J. Am. Chem. Soc., 1990, 112, 4192.
97. C. Yang, L. Zhang, C. Lu, S. Zhou, X. Li, Y. Li, Y. Yang, Y. Li, Z. Liu, J. Yang, K. N. Houk, F. Mo and X. Guo, Nat. Nanotechnol., 2021 DOI:10.1038/s41565-021-00959-4.
98. M. S. Hybertsen and S. G. Louie, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 5390.
99. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford, CT, 2016.
100. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
101. W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 56, 2257.
102. L. Goerigk and S. Grimme, Phys. Chem. Chem. Phys., 2011, 13, 6670.
103. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
104. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15.
105. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
106. D. J. Oliver, J. Maassen, M. El Ouali, W. Paul, T. Hagedorn, Y. Miyahara, Y. Qi, H. Guo and P. Grutter, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 19097.

---

Footnote

† Electronic supplementary information (ESI) available: Additional figures and tables, as well as listings of optimised coordinates. See DOI: 10.1039/d1sc04943g

---