Tuning the effective spin-orbit coupling in molecular semiconductors

The control of spins and spin to charge conversion in organics requires understanding the molecular spin-orbit coupling (SOC), and a means to tune its strength. However, quantifying SOC strengths indirectly through spin relaxation effects has proven difficult due to competing relaxation mechanisms. Here we present a systematic study of the g-tensor shift in molecular semiconductors and link it directly to the SOC strength in a series of high-mobility molecular semiconductors with strong potential for future devices. The results demonstrate a rich variability of the molecular g-shifts with the effective SOC, depending on subtle aspects of molecular composition and structure. We correlate the above g-shifts to spin-lattice relaxation times over four orders of magnitude, from 200 to 0.15 μs, for isolated molecules in solution and relate our findings for isolated molecules in solution to the spin relaxation mechanisms that are likely to be relevant in solid state systems.


Supplementary Note 1: Theoretical DFT calculations
This section provides additional detail on the theoretical calculations performed. As stated in the main text, the theoretical work is intended to support the key message that g-shifts in the studied class of molecules may be understood in terms of the effective spin-orbit coupling, or equivalently, the degree of overlap between the orbital-and spin angular momentum distributions. We argue that, given simple geometries and / or weaker SOC, we may predict the dominant ∆g OZ/SOC term of the g-tensor shifts from a linear model of atomically localized SOC interactions. The latter amounts to a linear fit of atomically localized spin populations to the ∆g OZ/SOC terms, on the form of Eq. 3. In other words, if a linear model with the distribution of atomic spin population as only independent variable is predictive, we consider our hypothesis confirmed. We validate the model by studying the correlation between the fitted ∆g OZ/SOC values, and those calculated from first-principles.
The theoretical methodology follows three steps, i) DFT calculations of molecular geometries, cation spin densities, and g-tensors, ii) calculations of atomic spin populations using an atom-inmolecule (AIM) decomposition method, and iii) a multivariate linear regression of the calculated dominant ∆g OZ/SOC g-shift term versus atomic spin populations. These three steps are described in reproducible detail in the following paragraphs.

DFT Calculations
We use DFT to predict g-tensor shifts from first-principles theory. These shifts depend on spin-orbit coupling (SOC), ranging from weak and nearly insignificant in the pure hydrocarbons studied, to significant in the Se-substituted molecules. We therefore require an accurate and transferable level of theory, able to describe a wide range of SOC strengths on a consistent footing.
Consequentially, we opted for an all-electron DFT method, with nuclear relativistic effects described by the zeroth-order regular approximation (ZORA 1 ) using the standard point-charge approximation for the atomic nuclei. g-tensors were calculated using the method 2 developed by Neese et al. and related techniques 3 as implemented in the ORCA software package 4 , version 3.0.3. Tests showed hybrid exchange-correlation functionals in general, and the PBE0 5 parameterization in particular, to perform very well for the prediction of the experimentally measured g-tensor shifts, in agreement with previous findings 2 . Exhaustive tests of the SARC basis set family 6 revealed that no qualitative features improved past the triple-zeta level (a.k.a. the 'ma-zora-def2-tzvp' basis set). The benchmark level quadruple-zeta 'ma-zora-def2-qzvp' basis set gives slightly larger (on the order of hundreds of [ppm]) g-shifts. Therefore, while the 'ma-zora-def2-tzvp' basis set, with linear dependencies eliminated by removing the diffuse (a.k.a. 'minimal augmentation') functions on the carbon atoms represents the optimal balance between accuracy of predictions and computational cost, the accuracy of predicted g-shifts is partially due to error cancellation between the basis set and other, over-estimating factors.
All molecules were fully geometry optimized in a charge neutral state and the above described level of theory including all-electron spin-orbit coupling, using the NWChem quantum chemistry software package 7 , version 6.5. Compared to a non-relativistic, spin-orbit free level of theory, the influence of the ZORA on the geometry is significant, but SOC effects negligible. Geometry differences between molecules optimized in a neutral and charged state were small, allowing for greater comparability by obtaining all geometries in the neutral state. The optimized pentacene and rubrene geometries have D 2H point-group symmetry. All other optimized molecules belong to the C 2H point-group. While experimental measurements were performed on molecules of varying alkyl chain lengths, e. g. C10, C12 etc., all molecules in the theoretical calculations were given 8-membered alkyl chains. Tests showed the alkyl chain influence to be well converged at this length, meaning that the theoretical results presented here are comparable to experimental measurements at chains equal to or longer than 8 methylene units.
For each optimized geometry, a single-point calculation with a total charge of +1 and multiplicity 2 was performed, with a subsequent calculation of the g-tensor and dumping of the DFT spin density onto a dense regular grid. The resulting calculated g-tensor shifts are shown in Supplementary Table  1 in parts per million (ppm). Evidently, the total g-tensor shift ∆g is strongly dominated by the orbital Zeeman / SOC-dependent term ∆g OZ/SOC in the studied class of molecules. That is, as mentioned in the main text, the gauge correction and relativistic mass correction terms add up to a very small, mostly positive term (67 ppm on average in Supplementary Table 1).

Atomic Spin Calculations
In the main text we show, that g-tensor shifts in the studied class of molecules can be understood in terms of a linear relationship between the dominant ∆g OZ/SOC term, and the spin density distribution among the nuclei in the molecule. More specifically, ∆g OZ/SOC is approximated as a sum over fitted proportionality constants representing the effective orbital angular momentum coupling for a given element, multiplied with sums of atomic spin populations at that element. This requires calculating atomic spin populations from the molecular spin density, in other words, finding an appropriate solution to the atoms-in-molecules (AIM) problem.
The AIM problem is commonly solved using methods partitioning either a) wave-functions based on projections of basis functions, or b) charge-or spin-densities based on their spatial distribution. The first option, exemplified by the Mulliken and Löwdin 8 population analysis methods, are conceptually simpler and readily implemented, but suffer from significant basis set effects 9,10 . Tests using the arguably more rigorous Löwdin method confirmed significant variations in calculated atomic spin populations with basis set type and size. The basis set dependence also manifested in a sensitive dependence on relative nuclear coordinates, appearing as an unexplained scatter in correlation plots of data fitted according to the procedure below. The second option, exemplified by the Bader method 11 , does not depend on the basis set used in the calculation of the molecular spin density. However, for the current set of molecules, the Bader method did not prove significantly more stable with respect to small variations in molecular structure than the above mentioned Löwdin method, nor did using the alternative weighting method proposed by Yu and Trinkle 12 alleviate the problem.
The Voronoi Deformation Density (VDD) method 10 proved significantly more stable and consistent, however, and is the method used in all fits of ∆g OZ/SOC presented in this manuscript. All density-based partitioning methods were applied using the 'bader' program developed by Henkelman et al. 13 , version 0.95a.

Fit of ∆g OZ/SOC vs Atomic Spins
For completeness, the fitted c e are presented in Supplementary Table 2. Since each molecule is excluded from the fitting set when fitted for, the c e are in principle unique for each molecule.
Consequentially, Supplementary Table 2 shows the average c e for each element in the complete set of 32 molecules studied, and the corresponding standard deviation.
As expected of a parameter representing the net effect of orbital angular momentum, c e increases dramatically going from carbon via sulphur to selenium. The reasons for the precise magnitude of the c e are convoluted, and we will here refrain from a detailed interpretation. It is however noteworthy that the fit is relatively consistent for the heavier elements, with standard deviations below a percent of the average value. Because of it's very small contributions to the ∆g OZ/SOC combined with the limited precision of input data, the carbon fits are considerably less numerically consistent, resulting in a relative standard deviation of approximately 70 percent.

Supplementary Note 2: Deviation of g-shifts for BSBS
All g-factor factors show an excellent agreement between experiment and theory apart from the values determined for BSBS. A systematic error in the DFT calculations is unlikely since only BSBS shows a large discrepancy. We consequently tried to eliminate all error sources for the ESR measurements. To rule out aggregation of the molecules in solution, we reduced the concentration of BSBS close to the detection limit of the instrument (0.05 × 10 −3 mol L −1 or 4 × 10 14 spins in the sample tube) and heated up the solution to 38 • C (boiling point of DCM: 39.6 • C) during the measurement. Neither resulted in any changes of the g-factor. In a control experiment, we investigated the solutions at different concentrations using diffusion-ordered nuclear magnetic resonance spectroscopy (DOSY NMR). These spectra show that the diffusion constant of BSBS increases marginally from 1.8 × 10 −9 m 2 s −1 to 2.1 × 10 −9 m 2 s −1 when reducing its concentration from 28 × 10 −3 mol L −1 to 0.28 × 10 −3 mol L −1 . This indicates the absence of aggregation which should otherwise result in a strong decrease of diffusion constants with increasing concentrations. Consistent with the difference in g-factors, we also observe different HFI coupling constants in theory and experiment: BSBS shows a root mean square (rms) deviation of ∼ 3 between theoretical and experimental HFI couplings compared to rms deviations < 0.8 for all other molecules. Aggregation seems unlikely as the cause for this difference, especially since the spread of the spin density over more than one molecule would increase the number of protons participating in the HFI. At present, we are not able to pinpoint the source of the discrepancy, apart from possible defects in the molecular structure. The hydrogen NMR spectrum ( Supplementary Fig. 3) reveals small impurity peaks in the aliphatic region though it remains unclear at present if they could cause such a strong shift in spin density.

Supplementary Note 3: ESR power saturation measurements
When the decay of spin polarization is an exponential process, it is possible to define spin lifetimes T 1 and T 2 and describe the time evolution of the spins by Bloch equations. The absorption rate of microwaves then follows a Lorentzian profile 14 where n 0 is the population difference between Zeeman split levels, 2π is the Planck constant, γ e the gyromagnetic ratio of an electron, ω the microwave frequency, B 1 the microwave magnetic field and ω 0 = γ e B 0 the Larmor frequency in an external magnetic field B 0 . P (ω) scales with B 2 1 in the absence of power saturation. The ESR signal is detected by a Schottky diode crystal and the detector current varies with the square root of the microwave power, i.e., linearly with B 1 : Here, χ 0 is the paramagnetic susceptibility of the sample. All measurements are recorded by a lock-in amplifier with a modulated external field B 0 and are thus recorded as a derivative of Eq. 2 with respect to B 0 . It is possible to determine the product T 1 T 2 of the two lifetimes by saturating the ESR signal with increasing microwave powers. The integral of Eq. 2 over all fields B 0 , or the double integral over the recorded derivative spectrum, scales with the microwave magnetic field B 1 as 15,16 The measured power saturation curves together with the fits to Eq. 3 are shown Fig. 3b. Before the onset of power saturation and line width broadening, the coherence time is given by T 2 = (γ e ∆B 1/2 ) −1 where ∆B 1/2 is the half width at half maximum of the Lorentzian resonance 15 . In practice, we use the line widths of the individual Lorentzian resonances that make up the spectrum, determined by the least-squares fitting of the hyperfine structure. Dividing the product of spin lifetimes by T 2 then gives the spin lattice relaxation time T 1 .
Microwave magnetic field in the cavity: Determining the product T 1 T 2 relies on the knowledge of B 1 over the sample volume while only the applied power P mw is recorded during an experiment. The scaling of B 1 with √ P mw depends on the resonator architecture and losses inside the cavity.
In the absence of power saturation, the integrated absorption signal of a paramagnetic sample scales only with its magnetic susceptibility and the field amplitude B 1 at the sample position. This allows us to determine B m from the integrated absorption signal of a reference sample with a known susceptibility. For the center of the cavity, we find as a function of the cavity Q-value and the incident microwave power P m : The vertical microwave magnetic field distribution f (y) is given by Bruker as a 9th order polynomial and stored with each data acquisition file ( Supplementary Fig. 2a). We additionally confirmed the values by monitoring the ESR signal intensity for a point like marker sample for different positions in the cavity (Fig. 3a). Since our samples are confined in a capillary tube with a radius of 0.5 mm around the cavity center, we only correct for the vertical field distribution by calculating the average microwave field across the sample length L and neglected any radial variations: A small error is introduced when using the average over the sample volume to fit the power saturation curves with Eq. 3 because we do not take into account that parts of the sample at positions with a lower field will contribute less to the total power saturation curve. In fact, the precise expression for the observed power saturation behavior would be as opposed to our approximation: However, the resulting error for our sample length of L = 20 mm remains well below 1% because the volume fraction at positions y = 0 is double that of the volume fraction at y = 0. A comparison of estimated and exact power saturation curves for a 25 mm long sample is shown in Supplementary  Fig. 2b, any visible differences disappear completely for sample lengths ∼ 20 mm. We also neglect the (small) distortion of the microwave field due to the sample's dielectric properties and therefore calculate with an error of ±8% in T 1 T 2 from the uncertainty in B 1 in addition to the standard error of from fitting.

Supplementary Note 4: DOSY NMR: diffusion constants and correlation times
DOSY NMR spectra were recorded on a Bruker AV-500 spectrometer at 298 K in deuterated dichloromethane. The diffusion experiments were performed using the Bruker ledbpgp2s pulse sequence. The number of data points per experiment was 32000 (td2) and the number of scans was 192. 24 experiments (td1) were collected using a diffusion delay (∆) of 28 ms (d20) and δ/2 pulse of 1.3 ms (p30).
The translational diffusion constants D tr from DOSY NMR allow us to calculate the hydrodynamic volume and Stokes radius r S of the respective molecule by the Stokes-Einstein equation with the Boltzmann constant k B , the sample temperature T = 295 K and the viscosity of DCM η = 4.13 × 10 −4 kg m −1 s −1 at room temperature. Knowledge of the hydrodynamic volume then allows us to estimate the rotational diffusion constant and the corresponding rotational correlation time, i.e., the time scale for which the rotational motion of the molecule is completely randomized, which is given by 14,15 :