Rapid and Precise Determination of Zero-Field Splittings by Terahertz Time-Domain Electron Paramagnetic Resonance Spectroscopy

Zero-field splitting (ZFS) parameters are fundamentally tied to the geometries of metal ion complexes. Despite their critical importance for understanding the magnetism and spectroscopy of metal complexes, they are not routinely available through general laboratory-based techniques, and are often inferred from magnetism data. Here we demonstrate a simple tabletop experimental approach that enables direct and reliable determination of ZFS parameters in the terahertz (THz) regime. We report time-domain measurements of electron paramagnetic resonance (EPR) signals associated with THz-frequency ZFSs in molecular complexes containing high-spin transition-metal ions. We measure the temporal profiles of the free-induction decays of spin resonances in the complexes at zero and nonzero external magnetic fields, and we derive the EPR spectra via numerical Fourier transformation of the time-domain signals. In most cases, absolute values of the ZFS parameters are extracted from the measured zero-field EPR frequencies, and the signs can be determined by zero-field measurements at two different temperatures. Field-dependent EPR measurements further allow refined determination of the ZFS parameters and access to the g-factor. The results show good agreement with those obtained by other methods. The simplicity of the method portends wide applicability in chemistry, biology and material science.


Introduction
Transition-metal or rare-earth molecular complexes and biological molecules assume well-defined symmetries, and may undergo structural distortion to lower symmetries due to interactions of the central metal ions with each ligand. 1 The molecular orbitals consequently undergo crystal field splitting, and valence electrons rearrange to form an energetically stable ground state. Zero-field splitting (ZFS) refers to the magnetic sublevel fine structure of unpaired electrons in such molecular orbitals in the absence of an external magnetic field. ZFS originates from spin-spin interactions mediated by the ligand field and from spin-orbit coupling. 1 The behavior of ZFS can be adequately captured by two parameters, D and E, which are the axial and transverse components of the magnetic anisotropy, respectively. 2 Accurate determination of the ZFS parameters for atomic centers with unpaired spins is critical in many research fields. For instance, in nanoscale thermometry based on nitrogen vacancy centers in diamond, the temperature dependence of D, the axial magnetic anisotropy, enables precise measurements of local temperature 3 . In metalloproteins, ZFS parameters are essential for interpreting EPR signals and understanding the electronic structure of metallo-cofactors 4 . In molecular magnets, D and E are related to both magnetic susceptibility and magnetization 5 ; measuring these parameters independently anchors the structure-function relationship that enables the design of molecules with higher effective barriers for spin inversion. Indeed, in such molecules, a negative D parameter indicates Ising-type magnetic anisotropy, and the potential energy surface of the spin sublevels has a double-well shape, with a spin inversion barrier that is proportional to D. 6 Due to their long relaxation and coherence times at low temperatures, single-molecule and single-ion magnets have been proposed to serve as the smallest units in high-density memory, quantum information processing and spintronics, 7 as such systems have two stable states that can be switched reversibly by magnetic fields, analogous to nitrogen vacancy centers in diamond 8 .
Conventional pulsed EPR spectroscopy using microwave technology has been widely applied to measure ZFS parameters between 1 and 100 GHz in the time domain. However, ZFS parameters of many metalloproteins and molecular complexes, in particular molecules with high magnetic anisotropy, lie in the THz frequency range between 0.1 and 10 THz (1 THz = 33.3 cm -1 ). Access to these higher frequencies is not available currently with routine magnetic resonance techniques. In high-frequency, high-field EPR (HFEPR), an external magnetic field is continuously swept to shift the spin resonances to the frequencies of the narrowband THz radiation sources used. 9,10 Leading work in this area is conducted with magnetic fields that typically can be varied continuously from 0 to 25 T. 11 Compounds with particularly large values of the ZFS parameters would require even higher magnetic fields which scale up with increasing ZFS parameters 12 . Because these measurements are conducted with applied magnetic fields, they also do not allow for direct measurements of the ZFS parameters, which are instead inferred from fitting the field-swept resonance data at discrete frequencies. Alternatively, frequency-domain Fourier transform EPR using THz sources from coherent synchrotron radiation 13 or blackbody radiation 14,15 has been used to measure ZFS parameters in single-molecule magnets 13,16 and biological systems 14,17 at zero and nonzero external magnetic field. Although this allows direct and broadband measurements of ZFS parameters and enhanced precision when the fieldswept broadband spectra are fitted by theoretical modeling 18 , the synchrotron THz sources employed thus far place severe limitations on wide applicability in the community. Inelastic neutron scattering has been used to directly measure ZFSs, but a typical experiment requires large amounts of sample (often more than 1 g), and is only available at specialized facilities with neutron sources. Most commonly, temperature-dependent magnetometry measurements have been used to determine ZFS parameters. Although more widely available, this technique is also arguably the most imprecise, because the measurements inevitably convolve ZFS parameters with other magnetic responses (e.g., exchange coupling, which is often of similar frequency to ZFS) and average over the temperature variation of ZFS values. Indeed, ZFS parameters determined by magnetometry often deviate from those measured by other methods. Clearly, the development of a simple technique employing benchtop equipment to directly and reliably measure absolute values of ZFS parameters could have a transformative effect in physical inorganic and organic chemistry, leading to more facile development of new magnetic materials and greatly facilitating the understanding of structure and function in metalloenzymes.
In this work, we present direct characterization of THz-frequency ZFSs in transition-metal complexes by measuring their EPR spectra using THz time-domain spectroscopy both at zero and nonzero external magnetic fields 19 . The broadband spectral coverage of our THz generation and detection methods readily allows for the measurement of EPR signals in transition-metal complexes with THz-frequency ZFSs at zero field. In cases where there is more than one spin transition, the ZFS parameter absolute values can be derived from a single EPR spectrum, and the spectral linewidths are typically narrower than those at nonzero magnetic fields, yielding optimal resolution. 20 The sign of the D parameter can be determined by zero-field EPR measurements at two different temperatures. The measurements can be supplemented in important ways with the addition of an applied magnetic field. Field-dependent EPR measurements allow the determination of ZFS parameters in spin systems where only one spin transition is present at zero field. Refined determination of the ZFS parameter values as well as access to the g-factor are also provided by fielddependent measurements. At the qualitative level, field-dependent measurements can unambiguously distinguish spin transitions from other low-frequency resonances including molecular or lattice vibrational transitions.

EPR spectra in zero and nonzero magnetic fields
The Hamiltonian 4,21 for a single electron spin includes the ZFS and electron Zeeman interaction (EZI) terms given by, The ZFS Hamiltonian is commonly written as, where ; ( = , , ) are spin matrices, is the total spin quantum number, and D and E are the axial and transverse ZFS parameters ( ≤ 3). 21 Diagonalization of the ZFS Hamiltonian yields the eigenenergies and eigenvectors of the magnetic sublevels of the spin system at zero static magnetic field. With a nonzero E value, the eigenstates are linear combinations of B states. Though B is a "good" quantum number only when E is zero, the magnetic sublevels are usually denoted by B . For integer spin systems, the degeneracy among the magnetic sublevels is completely removed by nonzero D and E. For half-integer spin systems, degeneracy is removed among different B states. Kramers degeneracy between ± B doublets remains, but can be removed by applying an external static magnetic field.
The static EZI term accounts for the spin system under a static magnetic field " . The EZI term is usually written as $%& = F 7 7 "7 + 8 8 "8 + 2 2 "2 , (1c) where F is the Bohr magneton, "; ( = , , ) is the static magnetic field, and ; ( = , , ) is the gfactor along each molecular axis. The application of a field induces splittings and shifts of the magnetic sublevels derived from the ZFS Hamiltonian. The new eigenenergies and eigenvectors of the spin system under a static external magnetic field can be sufficiently described by %() and $%& .
The magnetic field of electromagnetic radiation I interacts with the spin system also via the Zeeman interaction. It induces resonant magnetic dipole-allowed transitions between the magnetic sublevels. By measuring the transition frequencies, the ZFS parameters in the static spin Hamiltonian can be accessed at zero field. The Zeeman splittings and frequency shifts of the magnetic dipole-allowed transitions under a nonzero static magnetic field can also be measured, allowing determination of the g-factor.

Experimental
The experimental setup is a free-space THz time-domain spectroscopy system in transmission geometry as shown schematically in Fig. 1. (Details of the experimental setup are presented in the Supplementary  Information Fig. S1.) The THz emitter, which was a 1-mm thick (110)-cut zinc telluride (ZnTe) crystal, was illuminated by 800-nm pulses with 100 fs duration from a commercial Ti:Sapphire amplifier at a 5-kHz repetition rate. Single-cycle THz pulses with usable bandwidths spanning from 0.1 to 2.5 THz were generated by optical rectification of the laser pulses in the ZnTe crystal. 22 The broadband, linearly polarized THz pulses were focused onto the sample where the THz magnetic fields induced magnetic dipole-allowed transitions between the magnetic sublevels. 23 The resulting spin coherences radiated electromagnetic signals, known as free-induction decays (FIDs), at the resonance frequencies. We determined the electric field profiles of the FID signals with sub-ps resolution in the time domain by measuring the THz field-induced depolarization of a variably delayed 800-nm probe pulse in another ZnTe crystal. 22,24 The sample was placed in a helium cryostat with a split superconducting magnet that could be used to apply a static magnetic field B 0 in the 0-5.5 T range. The orientation of B 0 was perpendicular to the polarization of the THz magnetic field B 1 .

Fig. 1.
Schematic illustration of experimental geometry. Single-cycle THz pulses are generated by illuminating the THz emitter with fs laser pulses. The THz pulses are incident onto the sample to excite the spin transitions. The transmitted THz pulses and FID signals are directed into the THz detector. Fs laser pulses are overlapped with the THz fields at the THz detector for phase-resolved detection of the THz signals. The sample is placed in a cryostat with an external magnetic field B 0 perpendicular to the THz magnetic field B 1 (e.g. the Voigt geometry shown here). Fourier transformation of the THz signals yields the EPR spectra. In this work, the THz emitter and detector are both ZnTe crystals as described in the Supplementary Information.
We derive the THz Fourier-transform (FT) amplitude spectra KLM of the sample both at zero field and at discrete B 0 levels through numerical Fourier transformation of the FID signals. Absorbance spectra of the samples are obtained by comparing the sample spectra with the reference spectra PQR which were measured without the presence of the sample. The absorbance in units of optical density (OD) is given by (2)

Results and Discussion
To demonstrate the utility of the method, we chose prototypical samples that cover all of the common transition metal high-spin (HS) states and that are of interest in a diverse spectrum of research fields in which the knowledge of ZFS parameters is crucial. The following proof-of-principle samples were chosen: hemin, a compound that is related to heme-based enzymes, containing HS Fe(III) in a square-pyramidal environment 14,17 ; CoX 2 (PPh 3 ) 2 (X = Cl or Br, PPh 3 = triphenylphosphine), a known series of compounds exhibiting single-ion single-molecule magnetic behavior stemming from HS Co(II) in a pseudo-tetrahedral coordination environment 25 ; [Fe(H 2 O) 6 ](BF 4 ) 2 , a well-known integer-spin compound with HS Fe (II) in an octahedral environment 26 ; and NiCl 2 (PPh 3 ) 2 , an integer-spin compound with HS Ni(II) in pseudo-tetrahedral rather than the usual square planar geometry 18 . Microcrystalline powders of each compound were pressed into pellets which were used in the measurements (for details see Supplementary Information). Hemin has been under extensive study in EPR, as it is related to heme, which is the functional group in heme-based metalloproteins such as hemoglobin and myoglobin. The square-pyramidal structure of hemin is  Fig. 2(d). By measuring the frequency shifts of the transitions as a function of applied field strength, the values of the g-factor components can be obtained. The crystallites in the pellet sample are oriented randomly with respect to the applied magnetic field, so all three components can be determined. The raw time-domain FID signal and the FT amplitude spectrum measured from hemin at zero field and 20 K are shown in Fig. 3. The time-domain waveform of the THz pulse transmitted through the sample shows attenuation and slight broadening due to the THz absorption and dispersion in the sample. The FID signal is identified as the complex waveform profile following the transmitted THz pulse shown in Fig. 3(a). Numerical Fourier transformation of the THz time-domain signals yielded the complex FT spectra of the reference and sample. The amplitude spectra are shown in Fig. 3(b), where two dips are assigned as the magnetic dipole-allowed transitions indicated in Fig. 2(b). The absorption spectrum of hemin at 20 K is plotted in Fig. 3(c). In the raw spectrum, these two peaks sit on the wing of broad higher-lying absorptions that may be due to low-frequency vibrations. The background was subtracted manually to yield Fig. 3(d). The lineshapes were fitted to two Gaussian functions, yielding the peak frequencies 0.404 ± 0.001 THz and 0.809 ± 0.001 THz. The = 5 2 spin Hamiltonian was used to calculate the frequencies with variable ZFS parameters D and E to determine the absolute values = 6.74 ± 0.01 cm bI and = 0.048 ± 0.048 cm bI which show good agreement with previous measurements using frequency-domain FT THz EPR. 18 Fig. 4. Temperature-dependent zero-field absorbance spectra of hemin after background subtraction. The data are color-coded according to the temperatures. The higher-frequency peak disappears and the lowerfrequency peak shifts slightly at reduced temperatures.

High-spin Fe(III): Spin-5/2 system
The absorbance spectra of hemin at several temperatures are shown in Fig. 4. The time-domain data and FT spectra from which the absorbance spectra were determined are shown in Fig. S2 of the Supplementary Information. At temperatures above 20 K, both peaks resulting from the two magnetic dipole-allowed transitions are present. As the temperature decreases to 3 K, the lower-frequency peak becomes more pronounced while the higher-frequency peak disappears. This shows that at low temperature the spin populations are concentrated in the B = ± 1 2 states which must therefore be lowest in energy, indicating that the sign of D is positive. If D were negative, the B = ± 5 2 states would be lowest in energy and the higher-frequency transition would predominate at low temperature. Details about determination of the sign of D can be found in the Supplementary Information. The spectral peak at ~0.40 THz at 20 K shifts to ~0.41 THz at 3 K. This shift indicates a change in the ZFS parameter that may arise from subtle changes in the geometric or electronic structure of hemin at different temperatures.
Field-dependent measurements were conducted at 3 K and 20 K. B 0 was oriented along the THz propagation direction (i.e., Faraday geometry). The experimental absorbance spectra at zero field and six discrete B 0 levels are shown in Figs. 5(a) and 5(b) at 3 K and 20 K, respectively. The background was manually subtracted from all spectra. The time-domain data and FT amplitude spectra are shown in Fig. S3 of the Supplementary Information. The Kramers degeneracy is lifted and the ± B doublets are split at nonzero B 0 . Due to the anisotropy in the g-factor, the spectral peaks are shifted by the Zeeman interaction along all three molecular axes in the powder sample used, corresponding to the field-dependent frequency shifts shown in Fig. 2(d). A new magnetic dipole-allowed transition between the B = ± 1 2 states emerges at nonzero B 0 and is shifted from zero frequency to the spectral window shown as B 0 increases. The spectra at 20 K become somewhat crowded due to the presence of two peaks at zero field which both split to produce four peaks, plus the new peak starting from zero frequency, with some peaks merging to produce complicated lineshapes as the field strength is increased. The spectra at 3 K are simpler since only transitions originating from the B = ± 1 2 states appear. To quantitatively analyze the field-dependent spectra of our pellet samples with randomly oriented crystallites, we employed Easyspin, an EPR simulation software package 27 . A frequency-domain EPR simulation program was used to calculate the spectra for an = 5 2 spin system at the experimental temperatures and magnetic field levels. The input parameters include the total spin quantum number (set to 5/2), the ZFS parameters D and E, the g-factor elements 7 , 8 , and 2 , and the spectral lineshape and linewidth. The simulated spectra are shown in Fig. 5(c) and 5(d) for comparison with the experimental data. The frequency shifts and spectral lineshapes show good agreement between the experimental and simulated EPR spectra.
To eliminate possible errors introduced by the background subtraction, we plot in Figs. 6(a) and 6(b) the difference absorbance spectra between the spectra (without background subtraction) at successive B 0 levels.
The difference absorbance spectra can also provide enhanced sensitivity to spectral changes induced by B 0 . The Easyspin simulation results at both temperatures are plotted in thick dashed lines and are overlaid with the experimental difference spectra plotted in thin solid lines in Figs. 6(a) and 6(b). The simulation results at 20 K show good agreement with the experimental data in terms of the frequency shifts and the lineshapes. The simulation results at 3 K capture the transition peak between the B = ± 1 2 doublets due to 7 and 8 . However, the blue-shifted spectral peak originating from the transitions between B = ± 1 2 and B = ± 3 2 states due to 7 and 8 shows slightly smaller frequency shifts than the simulated ones. This discrepancy warrants further study. Based on the comparison between the simulated and experimental fielddependent EPR spectra, we can refine the quantitative determination of the parameters in the spin Hamiltonian of hemin. The parameters of the spin Hamiltonian determined by the simulations of the difference spectra are summarized in Table 1.

High-spin Co(II): Spin-3/2 systems
The pseudo-tetrahedral structure of CoX 2 (PPh 3 ) 2 (X = Cl or Br) is shown in Fig. 7(a). 25 The valence electrons (d 7 ) of the HS Co(II) indicate a total spin quantum number S = 3/2 due to three unpaired electrons.
The magnetic sublevels are derived by diagonalization of the ZFS Hamiltonian with S = 3/2 and are shown in Fig. 7(b) where a negative D value is assumed. The eigenstates are denoted by B . Magnetic dipoleallowed transitions are denoted by the double-sided arrows, and the frequency is given by 2 3 + 3 3 . Due to Kramers degeneracy between the ± B states, only one transition is observable at zero field, which is insufficient for separate determination of D and E or determination of the sign of D. The application of B 0 lifts the Kramers degeneracy and shifts the energies of the B states as shown in Fig. 7(c). The fielddependent transition frequencies are calculated and plotted in Fig. 7(d). Several magnetic dipole-allowed transitions emerge under a nonzero B 0 and allow separate determination of D and E. Measuring the transition between B = ± 1 2 states at different temperatures allows determination of the sign of D. The zero-field absorbance spectra of the two compounds at low temperatures are shown in Figs. 8(a) and 8(b). The time-domain signals and FT amplitude spectra are shown in Fig. S4 of the Supplementary Information. Several strong spectral peaks are all due to vibrations. One spectral peak due to the spin transition is indicated by the arrow in each figure. The magnetic origin of the transitions is confirmed in the field-dependent measurements discussed later. Each spectral peak is fitted to a Lorenztian function, with central frequencies 0.901 ± 0.001 THz for CoCl 2 (PPh 3 ) 2 Fig. 8. Zero-field absorbance spectra of CoCl 2 (PPh 3 ) 2 at 6 K (a) and CoBr 2 (PPh 3 ) 2 at 2 K. The spin resonance peak in each figure is indicated by the arrow. Other strong absorption peaks are due to vibrations.
Field-dependent measurements from 0 to 5.5 T were conducted on each sample. B 0 was oriented perpendicular to the THz propagation direction (i.e., Voigt geometry). The time-domain FID data and spectra are shown in Fig. S4 of the Supplementary Information. The field-dependent absorbance spectra are shown in Fig. 9 for the two compounds. The modulations are artifacts from the Fourier transformation as discussed in the Supplementary Information. At 1 T, the spin resonance peak in each compound exhibits a decrease in amplitude and broadening due to the anisotropic Zeeman interaction. At higher B 0 levels, the peaks split into several peaks, following the trends shown in Fig. 7(d). The transition between the B = ± 1 2 doublets emerges at nonzero B 0 and is expected to be around 0.3 THz at 5.5 T. The absence of these absorption peaks at 5.5 T at such low temperatures implies that the populations in the B = ± 1 2 states are depleted and the populations are concentrated in the lower-energy B = ± 3 2 states. The sign of the D parameter is therefore determined to be negative in each compound.
Quantitative analysis of the field-dependent spectra is conducted by comparing the experimental difference absorbance spectra with the simulated ones. The results are shown in Fig. 10. The simulation assumes negative D values for both compounds, which further confirms the determination of the sign of the D parameter. Separate determination of the D and E parameters is possible by analyzing the field-dependent spectra for = 3 2 systems. The relevant parameters of the spin Hamiltonian determined by the simulation are listed in Table 1. The measurements yield 3 + 3 3 values of 15.02 cm bI for CoCl 2 (PPh 3 ) 2 at 6 K and 14.00 cm bI for CoBr 2 (PPh 3 ) 2 at 2 K. In these two compounds, the spin transition peaks lie between two stronger absorptions due to vibrations. The peak vibrational absorptions are so strong that even small imperfections in their subtraction result in artifacts in the difference absorbance spectra.

High-spin Fe(II): Spin-2 system
The structure of Fe(H 2 O) 6 2+ is shown in Fig. 11(a). 29 The HS Fe(II) ion is in octahedral coordination with six water ligands. Four unpaired electrons of the valence electrons (d 6 ) indicate a total spin quantum number S = 2. The magnetic sublevel energy diagram assuming a positive D and a nonzero E is shown in Fig. 11(b). The new eigenstates are denoted by Φ ; with their eigenenergies labeled in Fig. 11(b). Magnetic dipoleallowed transitions are denoted by the double-sided arrows. As the Φ ; states are superpositions of the B states, six magnetic dipole-allowed transitions exist at zero field. As ≫ is typical, the splitting between Φ e and Φ f is small. The transitions to Φ e and Φ f states are often merged, which results in four distinct transitions. The values of the ZFS parameters D and E are adequately determined by a zero-field measurement of the frequencies I3 and Ig . The sign of D can be determined by temperature-dependent measurements at zero field. The application of B 0 further shifts the magnetic sublevels as shown in Fig. 11(c). The frequency shifts of the magnetic dipole-allowed transitions as a function of B 0 are shown in Fig. 11(d).  The zero-field absorbance spectra of Fe(H 2 O) 6 (BF 4 ) 2 at 1.8 K and 20 K are shown in Fig. 12. The timedomain signals and FT amplitude spectra are shown in Fig. S5 of the Supplementary Information. At 1.8 K, only the lowest state Φ I is populated. Two absorption peaks located through Lorentzian fits at 0.243 ± 0.001 THz and 0.423 ± 0.001 THz are assigned as the spin transitions at frequencies I3 and Ig between Φ I and Φ 3 and between Φ I and Φ g , respectively, as shown in Fig. 11(b). The assignments were confirmed by field-dependent measurements discussed below. Based on the values of I3 and Ig , the ZFS parameters are calculated to be = 10.82 ± 0.03 cm bI and = 1.00 ± 0.01 cm bI . At 20 K, four transitions are observed as the intermediate states Φ 3 and Φ g are also populated. The two absorption peaks at frequencies I3 and Ig become weaker, and two additional spectral peaks at ~0.90 THz and ~1.07 THz emerge due to the transitions between the higher-lying states. Hence the sign of the D parameter is determined to be positive. Due to the relative simplicity of the zero-field EPR spectra at 1.8 K where only two transitions appear, field-dependent measurements from 0 to 5.5 T were conducted at 1.8 K in the Voigt geometry. The experimental absorbance spectra are shown in Fig. 13(a). The time-domain signals and FT amplitude spectra are shown in Fig. S6 of the Supplementary Information. The spectral peaks assigned as spin resonances show splittings and shifts as a function of B 0 , which confirm their magnetic origins. To eliminate the background due to the wing of higher-lying absorptions by vibrations, difference absorbance spectra as a function of B 0 are shown in Fig. 13(b). The spectra are overlapped with simulated ones for comparison, which show excellent agreement. The spin Hamiltonian parameters determined by the simulations are listed in Table 1.

High-spin Ni(II): Spin-1 system
The pseudo-tetrahedral structure of NiCl 2 (PPh 3 ) 2 , which is similar to the structure of CoCl 2 (PPh 3 ) 2 , is shown in Fig. 14(a). Two unpaired electrons of the valence electrons (d 8 ) of the HS Ni(II) indicate a total spin quantum number S = 1. Similar to the case of the S = 2 spin system, with a zero E value, the degeneracy between the ± B states remains and the eigenstates are the B states. With a nonzero E value, the degeneracy between the B = ±1 states is lifted and the new eigenstates are superpositions of the B states. The magnetic sublevel energy diagram with positive D and nonzero E values is shown in Fig. 14(b) with new eigenstates denoted by Φ ; whose eigenenergies are indicated. Magnetic dipole-allowed transitions are denoted by the double-sided arrows. At zero field, two transitions are observable with frequencies denoted as I3 and Ig . The values of the ZFS parameters D and E are adequately determined by measuring I3 and Ig . The sign of D can be determined by temperature-dependent measurements at zero field. The application of B 0 shifts the magnetic sublevels as shown in Fig. 14(c). The frequency shifts of the magnetic dipoleallowed transitions as a function of B 0 are shown in Fig. 14(d) which allows the determination of the g-factor. The zero-field absorbance spectra of NiCl 2 (PPh 3 ) 2 at 2 K and 10 K are shown in Fig. 15. The time-domain signals and FT amplitude spectra are shown in Fig. S7 of the Supplementary Information. The dominant peak is likely due to a lattice vibrational mode. It is at a frequency similar to that of a peak in the spectra of CoCl 2 (PPh 3 ) 2 due to their similar molecular structures and atomic masses. Two weaker absorption peaks are assigned as the spin resonances. The ZFS parameters can be calculated based on the frequencies of these two peaks, which correspond to the transitions at I3 and Ig in Fig. 14(b). The frequencies were obtained by fitting the spectral lineshapes to two Gaussian functions, yielding I3 = 0.337 ± 0.001 THz and Ig = 0.459 ± 0.001 THz. The ZFS parameters were calculated to be = 13.27 ± 0.03 cm bI and = 2.03 ± 0.03 cm bI at 2 K, which show good agreement with previous HFEPR results 28 yielding = 13.196 cm bI and = 1.848 cm bI . The appearance of two peaks at both 2 K and 10 K implies that the lowest energy state is Φ I , where the population is concentrated. If Φ g were the lowest energy state, the population at Φ 3 would be mostly depleted and only one transition peak would be expected at 2 K. Hence the sign of the D parameter is determined to be positive. Fig. 15. Zero-field absorbance spectra of NiCl 2 (PPh 3 ) 2 at 2 K (red) and 10 K (blue). The two spin resonance peaks are shown by the arrows and are present at both temperatures. The stronger peak is due to a vibration. The field-dependent absorbance spectra of NiCl 2 (PPh 3 ) 2 at 2 K (Voigt geometry) are shown in Fig. 16(a). The time-domain signals and FT amplitude spectra are shown in Fig. S8 of the Supplementary Information. The two spin resonance peaks show field-dependent frequency shifts. The dependence of the spectral peaks on the magnetic field is analyzed by comparing the difference absorbance spectra with those obtained from simulations. The results are shown in Fig. 16(b), which shows excellent agreement between the experimental and simulated spectra. The spin Hamiltonian parameters determined by the simulations are listed in Table 1.

Conclusions
We demonstrate using four representative molecules a fast, facile, and reliable technique to measure ZFS parameters directly, in the absence of a magnetic field, and to refine the determination through fielddependent measurements. We use THz time-domain FID measurements to obtain EPR signals associated with THz-frequency ZFSs in molecular complexes. We fully characterized the values and signs of the ZFS parameters for several compounds belonging to = 1, S = 3 2, = 2 and = 5 2 spin systems based on the zero-field and/or field-dependent EPR spectra. Values of the g-factor are also obtained from the fielddependent measurements. This technique permits unambiguous assignment of THz-frequency ZFS parameters at different temperatures, which is difficult to accomplish by magnetometry measurements. More specifically, integer THz-frequency spin systems were termed "EPR silent" as the magnetic dipole-allowed transitions are not accessible in traditional EPR measurements due to the large ZFS parameters. HFEPR can access complexes with moderate ZFSs that are shifted into the excitation bandwidth via the Zeeman interaction, but the magnetic fields required for such measurements with monochromatic excitation sources can be quite large for complexes with high ZFS. The THz time-domain EPR measurement provides a direct way to measure large ZFSs in the THz-frequency region which is characteristic for molecular magnets.
In our current experiments, the spin number density was on the order of 10 21 cm -3 . Considering the ~5-mm THz beam diameter throughout the 2-mm thick pellet samples, the measured signals emerged from ~10 20 spins (see Supplementary Information for details). The estimated sensitivity of our current measurement configuration is 10 19 spins which is close to that required from a reasonable amount of large biomolecules 15,17 . In this work, the THz source and detector utilized allow us to measure resonances between about 6 and 80 cm -1 . The ranges of D and E that are accessible in this work depend on the specific spin system under study. Although we used homebuilt THz systems for our measurements, existing THz technologies 30 including some commercially available tabletop instruments can provide broader THz spectral coverage throughout the farinfrared range 31 , higher resolution 32 , higher sensitivity [33][34][35][36] , and faster data acquisition times 32,37 .These will allow determination of ZFS parameters and other magnetic fine structure revealed through THz-frequency spin transitions.
Because the technique is independent of the identity of the spin center and of the spin value, it can be applied broadly to systems with non-zero ZFS. We expect that the general approach in using THz timedomain spectroscopy to characterize transitions in the spin manifold of open-shell systems can be further elaborated to include other magnetic interactions, including magnetic exchange, for instance. Apart from the demonstrated advantages of THz time-domain spectroscopy over its frequency-domain counterpart, 38 the time-domain oscillation periods of ~1 ps will allow measurement of the dynamic evolution of unpaired electron spins as revealed by time-dependent changes in the ZFS spectra following, for example, pulsed optical excitation of a molecular state from which charge transfer or spin crossover occurs 39,40 . Twodimensional (2D) THz magnetic resonance spectroscopy of collective spin waves (magnons) has recently been demonstrated 41 , and 2D THz EPR measurements of HS compounds may also prove possible. THz EPR echoes (as observed in the 2D THz measurements of magnons) could prove useful for separation of spin transitions in biomolecules from low-frequency vibrations which may undergo very rapid dephasing, after which echo signals may be dominated by otherwise obscured spin coherences. The methodology presented herein and its future developments will find wide-ranging applications in characterizing magnetic properties of molecules, biological systems, and condensed matter. Our experimental setup for the tabletop THz time-domain EPR spectroscopy system is shown in basic form in Fig. 1 of the main paper and in more detail schematically in Fig.  S1. The laser was a commercial Ti:sapphire amplifier (Spitfire Pro, Spectra Physics) delivering 800 nm pulses with duration of 100 femtoseconds at a repetition rate of 5 kHz. The total output power of 400 mW was split into two optical paths by a 92/8 pellicle beamsplitter. The stronger pulses were modulated at 2.5 kHz by an optical chopper and were incident onto a 1-mm (110)-cut ZnTe crystal to generate THz pulses via optical rectification 1-3 . Single-cycle THz pulses were generated from the ZnTe crystal with a useable bandwidth spanning from 0.1 to 2.5 THz. The residual laser light transmitting through the ZnTe crystal was blocked by a black Teflon sheet. The THz pulses were collimated by a 45-degree off-axis parabolic mirror (PM) and focused onto the sample by a 90-degree off-axis PM (PM1 and PM2, respectively, in Fig. S1). The THz pulses transmitted through the sample and the FID signals that followed them were collimated and focused into a 2-mm ZnTe detection crystal by a pair of 90-degree off-axis PMs (PM3 and PM4 in Fig. S1). The weaker laser pulses from the beamsplitter were time delayed by a delay line (a mechanical translation stage) and attenuated by a half wave-plate (λ/2) and a polarizer (P). They were subsequently focused and overlapped with the THz beam in the ZnTe detection crystal to measure the phase-resolved THz signals via electro-optic sampling 4 . In this measurement, THz electric fields induced a modulation of the refractive indices of the ZnTe crystal along two orthogonal directions. The laser pulses experienced a transient birefringence due to this modulation. The THz-induced birefringence was measured as intensity modulations of the two optical polarization components which were separated by a quarter wave-plate (λ/4) and a Wollaston prism (WP) and detected by a pair of photodiodes (PD1 and PD2). The difference between the measured intensities was detected by a data acquisition card (DAQ) triggered by the chopper. The temporal profiles of the THz pulses and the FID signals were measured with sub-picosecond time resolution by scanning the delay line. The THz beam path was kept under dry air purge, which suppressed THz absorption due to water vapor in the atmosphere. The dynamic range of spectral amplitude of the system was in excess of 10 3 (corresponding to a dynamic range of spectral intensity in excess of 10 6 ). Further discussions of the underlying principles of THz generation and detection by optical rectification and electro-optic sampling can be found in References 1-3. The samples were placed in a helium cryostat with a split superconducting magnet (SuperOptiMag, Janis) which could provide static magnetic fields B 0 ranging from 0 to 5.5 T. The orientation of B 0 was perpendicular to the polarization of the THz magnetic field B 1 . B 0 can be either parallel or perpendicular to the propagation direction of the THz pulse. The former geometry is called Faraday geometry and the latter is called Voigt geometry 5 . These two geometries do not result in any differences in the EPR measurements. All the elements used in the setup are commercially available.

Experimental setup
The time-domain signals measured experimentally typically had 390 time steps of 67 fs. At each time point, the signal was typically averaged for 1000 laser shots. Under these conditions, the data acquisition time for a single absorbance spectrum was roughly 4 minutes, including measuring the time-domain signals for both the reference and the sample. The absorbance spectra reported below were averaged data from 10 measurements, which took roughly 40 minutes to collect. The time window of 26 ps was limited by the THz double reflection in the 1-mm ZnTe crystal used for THz generation. The resulting instrument-limited frequency resolution was approximately 39 GHz. Each spectrum reported here was interpolated by zero padding of the time-domain signals to 4096 data points in the data processing.
In some cases, the oscillatory FID signals emerging from the samples last longer than the instrument-limited time window of 26 ps. The spectra show ringing artifacts due to Fourier transformation of the raw time-domain signals with square windows. A Hamming apodization function was applied to the absorbance spectra, which reduced the ringing effects and had minimal effects on the frequencies of the spin resonance signals of interest.

Spin Hamiltonian
The spin Hamiltonian 6,7 discussed in the main text consists of the ZFS and EZI terms. The general form of the ZFS parameter is a second-rank tensor, which is set traceless (the sum of the diagonal components is zero) and symmetric ( ;i = i; ). The general form of the spin Hamiltonian describing the ZFS for a single spin system is written as where is the Kronecker delta. The ZFS energy levels of the spin systems studied in the main text can be derived from these equations.
The general form of the EZI term is written in the tensor form as where F is Bohr magneton and " = "7 "y "z is the applied static magnetic field vector, and is the g-factor, which is a tensor. The g-factor is usually symmetric and can be transformed into diagonal form, with the diagonal elements ; ( = , , ) determined from experimental measurements.

Determination of D and E parameters from zero-field and field-dependent EPR measurements
As shown in the main paper, zero-field EPR measurements yield the absolute values |D| and |E| of the ZFS parameters for integer spin systems, in which the spin sublevels are nondegenerate. For S = 3/2 systems, degeneracies among the sublevels limit the zero-field EPR measurements to determination of combinations of the ZFS parameters, as D and E cannot be separately determined from one doubly degenerate transition at 3 + 3 3 measured in zero-field EPR. Application of an external magnetic field shifts the levels, enabling separate determination of the absolute values |D| and |E| for S = 3/2 systems. For other half-integer spin systems, zero-field EPR can measure more than one (doubly degenerate) magnetic dipole-allowed transition derived from the magnetic sublevels, so |D| and |E| can be determined without application of an external field.
Variation of the temperature in the presence of a magnetic field allows determination of the sign of D for S = 3/2 systems and for S = 1 systems with zero E parameter values. In all other spin systems, variation of the temperature at zero magnetic field is sufficient for determination of the sign of D. The details with examples of S = 3/2, S = 5/2, S = 1 and S = 2 systems are elaborated briefly below. Though the magnetic sublevels are labeled with B , note that B is a "good" quantum number only for = 0. In each case, applying an external magnetic field along the molecular z-axis, the Kramers doublets are split which allows determination of the sign of D.
In order to determine the sign of the D parameter in S = 3/2 systems, EPR spectra with varying external magnetic field and temperature are necessary. As an example, we show in Fig. S9 the magnetic sublevel energy diagrams for S = 3/2 systems with positive or negative D. Without splitting the Kramers doublets, the EPR transition frequencies are both 3 + 3 3 . We assume an external magnetic field B z along the molecular z-axis, which can be achieved with a single-crystal sample. (For powder samples, the molecular orientation is random with respect to the magnetic field. Zeeman interactions with an anisotropic g tensor result in somewhat complex lineshapes and need to be analyzed numerically. ZFS and g tensor parameters can be obtained with high precision from field-swept EPR measurements combined with spin Hamiltonian simulations 5,8 .) The applied field splits the Kramers doublets, separating the two transitions that are overlapped in the zero-field spectrum and also introducing a new magnetic dipole-allowed transition at frequency I3 between the Φ I and Φ 3 levels (derived from ) = ± 1 2 states that were degenerate under zero field). The spectral amplitude of the new peak exhibits different trends for different signs of D as the temperature is varied, due to the Boltzmann factor. For positive D, states Φ I and Φ 3 have lower energies than states Φ g and Φ e and the I3 transition has an increasing spectral amplitude as temperature is reduced and population is increased in Φ I . For negative D, states Φ I and Φ 3 have the higher energies and the I3 transition has a decreasing spectral amplitude as temperature decreases. The external magnetic field also allows separate determination of D and E. The frequencies of the magnetic dipole-allowed transitions shown in Fig. S9

. (S5c)
By analyzing the field-dependences of these three transition frequencies from field-swept EPR spectra, we can obtain the g z factor and the D and E parameters. Though the magnetic sublevels are labeled with ) , note that ) is a "good" quantum number only for = 0. The transition frequencies are calculated through second-order perturbation theory. D and E can both be determined from the two transition frequencies measured at zero field. The sign of D can be determined through temperaturedependent measurements at zero field.
In half-integer spin systems other than S = 3/2, the degeneracies of magnetic dipoleallowed transitions are partially removed by nonzero D, which allows separate determination of |D| and |E| and determination of the sign of D at zero field. The magnetic sublevel energy diagrams for S = 5/2 systems are shown in Fig. S10. Zero-field EPR measurements at a fixed temperature measure the magnetic dipole-allowed transitions between B = ± 1 2 and B = ± 3 2 and between B = ± 3 2 and B = ± 5 2 , whose frequencies are calculated by second-order perturbation theory 9 to be 2 + 11.2 3 and 4 + 11.2 3 in both cases. Hence for spin-5/2 systems, D and E can be separately determined by measuring these two frequencies at zero external magnetic field. The sign of D can be deduced from zero-field EPR spectra with varying temperature. If D is positive, B = ± 1 2 states are lowest in energy. As temperature decreases, the lower-frequency transition amplitudes between B = ± 1 2 and B = ± 3 2 states increase monotonically while the higher-frequency transition amplitude between B = ± 3 2 and B = ± 5 2 states increases first and then decreases as the thermal population in states B = ± 3 2 decreases. If D is negative, B = ± 5 2 states are lowest in energy. As temperature decreases, the higher-frequency transition amplitude between B = ± 5 2 and B = ± 3 2 states increases monotonically while the lower-frequency transition amplitude between B = ± 3 2 and B = ± 1 2 states increases first and then decreases as the thermal population in states B = ± 3 2 decreases. Applying an external magnetic field further splits the degenerate doublets, allowing determination of the g-factor through field-dependent EPR measurements. For S = 1 systems with a zero (or nonzero) E parameter, zero-field EPR measurement yields in both cases. An external magnetic field is required to split the doublet with ) = ±1 as shown in Fig. S11. The sign of D can then be determined by the temperature-dependent changes in spectral amplitudes of the magnetic dipole-allowed transitions. As temperature decreases, the spectral amplitudes at I3 and Ig both increase monotonically if D is positive. If D is negative, the higher-frequency transition amplitude at Ig increases monotonically while the lower-frequency transition amplitude at I3 increases first and then decreases as the thermal population in state Φ 3 decreases. Nonzero E values also can be determined through zero-field EPR measurements as illustrated in the main paper.
The magnetic sublevel energy diagrams for S = 2 systems are shown in Fig. S12. With E = 0, zero-field EPR measurements at a fixed temperature measure the magnetic dipole-allowed transitions between B = 0 and B = ±1 and between B = ±1 and B = ±2 , whose frequencies are and 3 respectively regardless of the sign of D. As temperature decreases, the lower-frequency transition amplitudes between B = 0 and B = ±1 states increase monotonically while the higher-frequency transition amplitude between B = ±1 and B = ±2 states increases first and then decreases as the thermal population in states B = ±1 decreases if D is positive. If D is negative, as temperature decreases, the higher-frequency transition amplitude between B = ±1 and B = ±2 states increases monotonically while the lower-frequency transition amplitude between B = 0 and B = ±1 states increases first and then decreases as the thermal population in states B = ±1 decreases. With nonzero E, all the degeneracy is removed as shown in the main paper. Once again, based on different trends of the spectral amplitudes as functions of temperature one can determine the sign of D. Figure S12. Magnetic sublevel energy diagrams for S = 2 systems with (a) positive or (b) negative D. With either zero E or nonzero E, the sign of D can be determined from the temperature-dependent spectral amplitudes.
Approximately 200 mg of hemin and CoX 2 (PPh 3 ) 2 were pressed into pellets of 13-mm diameter and ~2 mm thickness. Approximately 200 mg of Fe(H 2 O) 6 (BF 4 ) 2 and NiCl 2 (PPh 3 ) 2 were mixed with 100 mg high-density polyethylene and the solid solutions were pressed into pellets of 13-mm diameter and ~2 mm thickness. The measurements were conducted on these pellets and the results are discussed as follows. Pellets of hemin (nominally pure powders), CoX 2 (PPh 3 ) 2 (nominally pure powders) and NiCl 2 (PPh 3 ) 2 (nominally pure powders mixed with HDPE) were pressed to ~2 mm thickness in air with a manual hydraulic press (MTI) using a 13 mm stainless steel pellet die (Specac) with 3 tons of applied mass. Pellets of Fe(H 2 O) 6 (BF 4 ) 2 (nominally pure powders mixed with HDPE) were pressed to ~2 mm thickness in a N 2 -filled glovebox (Innovative Technology) with a manual hydraulic press (Specac) using a 13 mm stainless steel pellet die (Specac) with 3 tons of applied mass.
Though the diameters of the pellets were 13 mm, the focused THz beam with a spot size of approximately ~5 mm diameter throughout the 2-mm samples defined the ~40 mm 3 sample volumes that were measured. The spin number density in the samples was on the order of 10 23 cm -3 . The measured signals emerged from approximately 10 20 spins within the spot of the focused THz beam.

Powder X-ray diffraction
Powder X-ray diffraction (PXRD) patterns were recorded with a Bruker D8 Advance diffractometer equipped with a Göbel mirror, rotating sample stage, LynxEye detector and Cu Kα (λ = 1.5405 Å) X-ray source in a θ/2θ Bragg-Brentano geometry. An anti-scattering incident source slit (2 mm) and an exchangeable steckblende detector slit (8 mm) were used. The tube voltage and current were 40 kV and 40 mA, respectively. Knife-edge attachments were used to remove scattering at low angles. Samples for PXRD were prepared by placing a thin layer of the designated materials on a zero-background silicon (510) crystal plate.