Charge and spin transport in single and packed ruthenium-terpyridine molecular devices: Insight from first-principles calculations

Using first-principles calculations, we study the electronic and transport properties of rutheniumterpyridine molecules sandwiched between two Au(111) electrodes. We analyse both single and packed molecular devices, more amenable to scaling and realistic integration approaches. The devices display all together robust negative differential resistance features at low bias voltages. Remarkably, the electrical control of the spin transport in the studied systems implies a subtle distribution of the magnetisation density within the biased devices and highlights the key role of the Au(111) electrical contacts.

the first-principles calculations we used the SIESTA and SMEAGOL codes. Electronic exchange and correlation were treated within the generalized gradient approximation. The two-dimensional Brillouin zone sampling was limited to a uniform 2 x 2 grid of k-points and the convergence threshold for the residual forces during atomic-structure relaxation was set to 0.01 eV/Å. The relaxed structures can be characterized by the bonding of the molecular units to the Au(111) surfaces and the molecular geometry deformation. 1 The atomic displacements at the electrode surface are presented below.
To quantify the distortion of the electrode surface in the studied systems, we computed the displacement of the Au atoms along the z-coordinate (i.e. the transport direction) from their averaged positions in the relaxed system: where N is the number of the atoms in the top layer of the surface, z i is the current atomic position and z is the average value of the z-coordinate of the atoms in the top layer. For all models, the surface distortion is found to be negligible, with δ z ≈ 0.002Å, which is much smaller than the calculated distancesd S−Au , as summarized in Table S1.

Computed electronic properties of the devices
For the interpretation of the transmission spectra T(E) in the considered systems, related to the zero-bias conductance, we are referring to the projected density of states (PDOS) of the Ru, N, C and S atoms. The latter is computed using SIESTA on the whole system and applying a broadening of 0.1 eV. We represent the PDOS and the T(E) results for Models 1 -4 in Figure S1 and for Models 5, 6, 2B, and 2P in Figure S2.

Single and packed ruthenium-terpyridine molecules between two
Au(111) electrodes: Negative differential resistance The nonequilibrium transport method used for the SMEAGOL calculations is described in detail in Ref. 2. Basically, the current I flowing through the molecular junction is given by: Model 1

Model 2
Model 3 Model 4 Figure S1: Transmission T(E) and atom-projected density of states (PDOS) for Models 1 -4. We use red curves for spin-up and blue curves for spin-down data. The Fermi level was set to zero. For each atom-projected density of states, the units correspond to 3 states/eV, while the transmission results are displayed between 0 and 1.

Model 6
Model 2B Model 2P Figure S2: Transmission T(E) and atom-projected density of states (PDOS) for Models 5, 6, 2B and 2P. We use red curves for spin-up and blue curves for spin-down data. The Fermi level was set to zero. For each atom-projected density of states, the units correspond to 3 states/eV, while the transmission results are displayed between 0 and 1.
where T(E, V) is the total (spin-up plus spin-down) transmission coefficient represented in Figure S3, f (E) is the Fermi function at a given temperature, and µ L and µ R are the chemical potentials of the left and right leads, respectively. Due to the difference between the two Fermi functions, the integration range is limited to the energy window between µ L and µ R (with a broadening related to the temperature). When a bias V is applied to the molecular junction, µ L (resp. µ R ) is shifted up (resp. down) by eV/2: as illustrated in Fig The I(V) characteristics in Fig. 2 of the main text can be easily understood by inspect-ing Figure S3. For Model 2 (representative for Models 1-6), the peak P 1 in T(E, V) around (0.1 eV, -0.2 V) is responsible for the peaks at V = ±0.2 V in the I(V) curve. The peaks P 2 around (0.7 eV, 1.0 V) and P 3 around (0.9 eV, 1.4 V) which lie slightly outside the integration range, translate into a plateau in the current. So, the NDR is related to the valley between the peaks P 1 and P 2 . For the Model 2P, the peak P 1 around (0.3 eV, 0.3 V) is quite extended and gives rise to the peaks at V = −0.2 and 0.4 V (due to its asymmetric shape).
On the positive bias side, the current peak is followed by a dip at ∼0.6 V, which induces a NDR, due to the valley between P 1 and P 2 around (0.5 eV, 0.9 V). In contrast with Model 2, the separation between the peaks P 2 and P 3 around (1.0 eV, 1.9 V) is pronounced and induces a new dip in the current, and hence a NDR behaviour.

Au(111) electrodes: Polarization
Insight into the behaviour of the molecular junctions can be obtained by re-plotting the I(V) curves as polarization traces P(V), as shown in Figure S4. We observe that for Models 1 -3 the geometry of the devices is not strongly correlated with the polarization, which does not change if the devices are probed at negative or positive biases. The calculated P(V) for the packed system could be understood as arising from the lowest unoccupied molecular orbital (LUMO) given by the superposition of the two LUMOs of the constituent molecules (see below).

Au(111) electrodes: Magnetization
The calculated magnetization curves m(V) are displayed in Figure S5 for all investigated cases. The shape of the magnetization peak around zero bias is indicative for a peculiar distribution of the magnetization density across the molecular constituents of the devices.   . We use red curves for spin-up, blue curves for spin-down data, and green curves for their difference. The Fermi level was set to zero. For each atomprojected density of states, the units correspond to 3 states/eV, while the transmission results are displayed between 0 and 1.

The space-resolved spin density
The space-resolved spin density for the Model 2P is represented at various bias voltages (V = 0.0, 0.2, 1.0 and 2.0 V) in Figure S7 . The overall behaviour of the packed rutheniumterpyridine molecules linked between two Au(111) electrodes is imposed by the Molecule 2, which is tightly linked to the bottom electrode. At large bias configuration, the magnetisation density is essentially localised on the terpyridines linked to the bottom electrode.  Figure 1 (main text). The isosurfaces in green and violet correspond to a magnetization density of ±2 × 10 −3 µ B /Bohr 3 . Figure   HOMO LUMO Figure S8: Energy dependence of the density of states for the free radical. The vacuum energy was set to zero. Contour plots for the HOMO and LUMO orbitals are drawn in the inset. The isosurfaces in green and violet correspond to ±2 × 10 −3 e/Bohr 3 . Note that, due to the symmetry of the molecule, the HOMO-1 (resp. LUMO+1) is degenerate with the HOMO (resp. LUMO). The corresponding contour plots (not reported here for sake of clarity) are very similar to the HOMO and LUMO shown in the figure, but spatially localized on the other terpyridine unit.

Computed electronic properties of the isolated molecule
Since the conventional DFT exchange-correlation functionals (such as PBE) tend to underestimate the difference between the HOMO and LUMO energies (E g ), we have also performed the calculations using hybrid functionals (PBE0 and B3LYP) and computed many-body corrections (using both the one-shot G 0 W 0 method and the self-consistent evGW approach with self-consistency on the eigenvalues) for the isolated molecule. The calculations with the hybrid functionals where performed using the NWCHEM code, while the many-body corrections were obtained with the FIESTA code. The many-body corrections were obtained starting from the PBE0 eigenvalues and wavefunctions. For both DFT and many-body calculations, changing the basis set from cc-pVTZ to aug-pVTZ affected the value of E g by less than 0.1 eV. The results are presented in Table S2.
Table S 2: Difference between the HOMO and LUMO energies (E g ) as calculated using DFT with the PBE, PBE0, and B3LYP exchange-correlation functionals and including many-body corrections (using both the one-shot G 0 W 0 method and the self-consistent evGW approach with self-consistency on the eigenvalues).
Method PBE PBE0 B3LYP G 0 W 0 evGW E g 0.5 0.9 0.7 2.3 2.5 These values can be used to correct the eigenvalues for the devices using the the DFT+Σ method. [3][4][5] This requires to also take into account an "image charge" term, accounting for the polarization energy associated with static nonlocal correlations between the electrons on the molecule and in the metal electrodes that closes the gap of the molecule upon absorption. Typically, this will reduce the difference between the HOMO and LUMO energies of the isolated molecule by 35-45%. 6 As a result, we expect a separation of ∼1.5 eV for the system considered here.