Multireference Ab Initio Studies of Magnetic Properties of Terbium-Based Single-Molecule Magnets

TL;DR

This study uses relativistic multireference ab initio calculations to analyze how chemical environment and molecular geometry influence the magnetic properties and quantum tunneling behaviors of terbium-based single-molecule magnets.

Contribution

It provides detailed computational insights into the electronic structure and magnetic anisotropy of Tb-based SMMs, highlighting the effects of ligand type, geometry, and oxidation state.

Findings

Ligand type and geometry minimally affect zero-field splitting.

Oxidation number weakly influences zero-field splitting.

Molecular symmetry and ligand distortion impact tunnel splitting and quantum tunneling.

Abstract

We investigate how different chemical environment influences magnetic properties of terbium(III) (Tb)-based single-molecule magnets (SMMs), using first-principles relativistic multireference methods. Recent experiments showed that Tb-based SMMs can have exceptionally large magnetic anisotropy and that they can be used for experimental realization of quantum information applications, with a judicious choice of chemical environment. Here, we perform complete active space self-consistent field (CASSCF) calculations including relativistic spin-orbit interaction (SOI) for representative Tb-based SMMs such as TbPc$_2$ and TbPcNc in three charge states. We calculate low-energy electronic structure from which we compute the Tb crystal-field parameters and construct an effective pseudospin Hamiltonian. Our calculations show that ligand type and fine points of molecular geometry do not affect the zero-field splitting, while the latter varies weakly with oxidation number. On the other hand, higher-energy levels have a strong dependence on all these characteristics. For neutral TbPc$_2$ and TbPcNc molecules, the Tb magnetic moment and the ligand spin are parallel to each other and the coupling strength between them does not depend much on ligand type and details of atomic structure. However, ligand distortion and molecular symmetry play a crucial role in transverse crystal-field parameters which lead to tunnel splitting. The tunnel splitting induces quantum tunneling of magnetization by itself or by combining with other processes. Our results provide insight into mechanisms of magnetization relaxation in the representative Tb-based SMMs.

Tables

Table 1. Table 1: Structural Parameters a for the Considered Experimental Geometries. The Errors Correspond to the Standard Deviation from the Average. Numbers in Parentheses Denote the Range of Observed Values where the Tb ion is at the origin. The Distances Are in Å and Angles in Degree

	[TbPc2]0 M1	[TbPc2]0 M2	[TbPcNc]0 M3	[TbPc2]- M4	[TbPc2]- M5	[TbPcNc]+ M6
Counter ionb	n/a	n/a	n/a	[TBA]+$=$[N(C4H9)4]+	[TBA]+	(PF6)-
Solvent or	CH2Cl2	n/a	CHCl3	CH3OH $\cdot$ H2O	3([TBA]+Br-)	2 CH2Cl2
diultionb					$\cdot$ 3$H_2$O
Crystal packingc	$P_{nma}$ (62)	$P$212121 (19)	$P$4212 (90)	$P_{na}$21 (33)	$P$21$c$ (56)	$P$4212 (90)
Symmetryd	$C_s$	$C_1$	$C_4$	$C_1$	$C_1$	$C_4$
$D_{{\text{TbN}}_{\text{nn}}}$e	2.42$\pm$0.01	2.41$\pm$0.01	2.42$\pm$0.02	2.43$\pm$0.01	2.44$\pm$0.10	2.41$\pm$0.01
$D_{{\text{TbC}}_{\text{far}}}$f	6.79$\pm$0.05	6.82$\pm$0.04	6.85$\pm$0.03	6.88$\pm$0.07	6.91$\pm$0.41	6.85$\pm$0.03
${\theta }_{{\text{N}}_{\text{nn}}{\text{TbN}}_{\text{nn}}}$g	(44.4,45.6)	(40.6,49.5)	(43.7,46.4)	(43.3,46.8)	(34.2,49.8)	(43.4,46.6)
${\theta }_{{\text{C}}_{\text{far}}{\text{TbC}}_{\text{far}}}$h	(44.6,45.4)	(40.7,49.9)	(43.5,46.5)	(42.8,46.8)	(34.5,51.8)	(43.4,46.6)
$Z_{{\text{N}}_{\text{nn}}^{\text{up}}}$j	1.41$\pm$0.00	1.39$\pm$0.01	1.39$\pm$0.00	1.41$\pm$0.02	1.41$\pm$0.03	1.38$\pm$0.00
$Z_{{\text{N}}_{\text{nn}}^{\text{down}}}$i	-1.41$\pm$0.01	-1.40$\pm$0.02	-1.41$\pm$0.00	-1.41$\pm$0.01	-1.42$\pm$0.04	-1.40$\pm$0.00
$Z_{{\text{C}}^{\text{up}}}$k	(1.49,1.66)	(1.51,2.00)	(1.51,1.78)	(1.53,2.68)	(1.46,2.45)	(1.50,1.81)
$Z_{{\text{C}}^{\text{down}}}$l	(-2.38,-1.47)	(-2.26,-1.49)	(-2.12,-1.57)	(-2.83,-1.52)	(-2.99,-1.36)	(-2.19,-1.55)
$Z_{{\text{C}}_{\text{down}}^{\text{ave}}}$;$Z_{{\text{C}}_{\text{up}}^{\text{ave}}}$m	-1.70;1.59	-1.68;1.69	-1.54;1.66	-1.90;1.92	-1.84;1.91	-1.57;1.67
Reference	32	31	33	27	29	33

Table 2. Table 2: Calculated Magnetic Energy Scales for Different TbPc 2 -Type SMMs in Units of cm -1

	M1	M2	M3	M4	M5	M6
$E_{\text{TS}}^{1,2}$a	n/a	n/a	n/a	0.000	0.007	0.000
$E_{\text{TS}}^{3,4}$b	n/a	n/a	n/a	0.015	0.090	0.000
$E_{\text{TS}}^{5,6}$c	n/a	n/a	n/a	0.530	7.969	0.499
$\Delta E_{\text{ex}}$d	8.2	6.6	5.1	n/a	n/a	n/a
$J_{\text{ex}}$e	0.8	0.6	0.6	n/a	n/a	n/a
$E_{\text{ZFS}}$f	308	305	308	289	292	298
$E_{\text{MAB}}$g	658	639	622	582	734	592
$\Delta E_J$h	2083	2068	2068	2036	2045	2049

Table 3. Table 3: Calculated Principal Values a of the g 𝑔 g -Tensor for Different TbPc 2 -Type SMMs

	M1b	M2b	M3b	M4	M5	M6
$g_x$	1.505	1.505	1.504	1.505	1.510	1.504
$g_y$	1.504	1.503	1.504	1.502	1.498	1.504
$g_z$	1.476	1.477	1.477	1.480	1.477	1.478

Table 4. Table 4: CF Parameters a , b in cm -1 for Tb Ion for M4, M5 and M6 SMMs .

$B_k^q$	M4	M5	M6
$B_2^{-2}$	-0.12802616	-0.02332695	0.00000240
$B_2^{-1}$	-0.03939183	-0.00727754	0.00001437
$B_2^0$	-5.05068413	-5.52638125	-5.33499223
$B_2^1$	-0.02663218	-0.00634158	-0.00001624
$B_2^2$	0.80463902	3.45788034	0.00003865
$B_4^{-4}$	-0.00332496	-0.00213682	0.00610215
$B_4^{-3}$	0.00793656	-0.00424813	-0.00000026
$B_4^{-2}$	0.00410350	0.00096348	-0.00000007
$B_4^{-1}$	-0.00375441	0.00049609	0.00000289
$B_4^0$	-0.01406960	-0.01300444	-0.01412822
$B_4^1$	-0.00514061	0.00707655	-0.00000337
$B_4^2$	-0.00151449	-0.02369506	-0.00000004
$B_4^3$	-0.00030393	0.00188537	0.00000259
$B_4^4$	-0.00079989	0.01117566	0.00540901

Table 5. Table 5: Measured effective magnetic anisotropy barrier (in cm -1 ) for TbPc 2 -like single-molecule magnets

Molecule	TbPc2a	TbPc2b	TbPc2b	TbL2c,d	TbL2c,d	TbPcNce	TbPcNce
Charge	-1	-1	-1	-1	+1	+1	0
Dilution	Yes	No	Yes	No	No	No	No
MAB	260	584	641	594	550	584	342
Relevance	n/a	M5	M5	n/a	n/a	M6	M3

Table 6. Table 6: CF Parameters a , b in cm -1 for Tb Ion for M1, M2 and M3 SMMs

$B_k^q$	M1c	M2c	M3c
$B_2^{-2}$	0.00000001	-0.04824674	-0.00000478
$B_2^{-1}$	-0.00000003	0.03000187	0.00000119
$B_2^0$	-5.93432163	-5.68334414	-5.65157592
$B_2^1$	-0.00444594	-0.04186354	-0.00001973
$B_2^2$	0.32444131	0.44453996	0.00019387
$B_4^{-4}$	0.00000000	-0.01188166	-0.00302990
$B_4^{-3}$	0.00000000	-0.00089983	0.00000000
$B_4^{-2}$	0.00000000	0.00162901	0.00000014
$B_4^{-1}$	0.00000000	-0.00034396	0.00000135
$B_4^0$	-0.01336535	-0.01369266	-0.01380654
$B_4^1$	-0.00043523	-0.00806695	-0.00000275
$B_4^2$	-0.00143680	-0.00160371	-0.00000081
$B_4^3$	0.00160576	-0.00288249	-0.00000487
$B_4^4$	-0.00027466	-0.01198852	-0.00720075

Equations

$$\hat{H}=\sum _{k=2,\mathrm{even}}\sum _{q=-k}^k{B}_k^q{\hat{O}}_k^q(\hat{\mathbf{J}})-{\mu }_B\mathbf{B}\cdot \mathbf{g}\cdot \hat{\mathbf{J}},$$

$$\hat{H}=\sum _{kq}{B}_k^q{\hat{O}}_k^q(\hat{\mathbf{J}})+J_{\text{ex}}\hat{\mathbf{s}}\cdot \hat{\mathbf{J}}-{\mu }_B\mathbf{B}\cdot \left(\mathbf{g}\cdot \hat{\mathbf{J}}+g_e\hat{\mathbf{s}}\right),$$

Full text

\SectionNumbersOn

1 Introduction

Single-molecule magnets1, 2, 3 (SMMs) are magnetic molecules typically composed of one or several transition metal or lanthanide ions surrounded by ligands. The key feature of SMMs is existence of inherent magnetic anisotropy due to interplay between the ligand crystal field (CF) and spin-orbit interaction (SOI). This feature allows one to explore interesting phenomena and potential applications of SMMs for magnetic information storage,4, 5 spintronics6, 7, 8, 9, 10 and quantum information processing.11, 12, 13, 14, 15, 16, 17 Recently, it was reported that monometallic lanthanide-based SMMs can have an effective energy barrier of over 1000 cm*-1* with magnetic hysteresis above liquid nitrogen temperature.5, 18 For reviews of lanthanide-based SMMs, see Refs.19, 20, 21, 22, 23. Furthermore, since the first proposal of the implementation of Grover’s algorithm into the prototype SMM Mn12,13 quantum bits (qubits) and quantum gates based on SMMs14, 15, 16, 17 have been experimentally realized. One promising SMM candidate is a double-decker TbPc2 family (Pc = phthalocyanine),24 where Rabi oscillations and Grover’s algorithm were experimentally implemented by tuning the nuclear spin states of the Tb ion with AC electric fields within single-molecule transistor set-ups.14, 17

SMM TbPc224, 25 consists of a Tb3+ ion sandwiched between two Pc ligands, as shown in Fig. 1a. Three different charge states were experimentally realized for TbPc2 molecules such as anionic [TbPc2]-,26, 27, 28, 29, 30 (Fig. 1(b)) neutral [TbPc2]0,31, 32 (Fig. 1(a)) and cationic [TbPc2]+.28 In the latter case, however, X-ray crystallography data have not been reported yet. Recently, [TbPcNc]0,+ molecules, where one of the Pc ligand rings was replaced by a larger Nc (naphthalocyaninato) ligand (Fig. 1c), were experimentally studied.33 The synthesized TbPc2 (TbPcNc) molecules have only approximate $D_{4d}$ ($C_{4v}$) symmetry. The degree of symmetry deviation varies with crystal packing, diamagnetic dilution molecules, or solvent molecules used in synthesis processes.

Both TbPc2 and TbPcNc molecules in different charge states were experimentally shown to exhibit SMM behavior.24, 26, 28, 31, 30, 32, 27, 29, 33 For the TbPcNc molecule, the measured effective energy barrier is in the range of 340-580 cm*-1* depending on oxidation number,33 while for the TbPc2 molecule, the measured barrier is in the range of 230-640 cm*-1* depending on oxidation number.24, 26, 28, 30, 29 There are no theoretical studies of the origin of this wide range of the energy barrier for TbPcNc or TbPc2 molecules. Observed magnetization relaxation in SMMs may arise from combined contributions of different relaxation mechanisms such as quantum tunneling of magnetization (QTM), Raman process, Orbach process, hyperfine coupling, and/or intermolecular interaction.34 In extracting the experimental barrier, it is often assumed that there are no intermolecular interactions and hyperfine interactions. The experimental barrier may also depend on magnetization relaxation mechanisms considered in the fitting of experimental data. Considering the complexity of the relaxation mechanisms and the assumption and ambiguity in the fitting process, it may be difficult to unambiguously determine the magnetization relaxation mechanisms solely from measurements. Furthermore, depending on experimental set-ups, additional environmental factors qualitatively affect the magnetic properties of the SMMs.

In TbPc2,24, 25 the Tb3+ ion has the 4f8 electronic configuration. According to Hund’s rules, its ground multiplet corresponds to spin, orbital and total angular momentum quantum numbers of $S=3$, $L=3$ and $J=6$, respectively. The unquenched orbital angular momentum gives rise to strong SOI. CF of Pc ligands in conjunction with the SOI splits the Tb $J$-multiplets (Fig. 2), leading to large magnetic anisotropy. CF parameters are often referred to as magnetic anisotropy parameters. As shown in Fig. 2, for low-energy $J$-multiplets, different $J$-multiplets are well separated from one another. In this case, assuming uniaxial magnetic anisotropy, total angular momentum projected onto the magnetic easy axis ($M_{J}$) remains a good quantum number. For a given $J$-multiplet, states with the same magnitude of $M_{J}$ are degenerate, whereas states with different $|M_{J}|$ values are split. Here we define magnetic anisotropy barrier (MAB) as the energy difference between the lowest and highest magnetic levels within the ground $J$-multiplet. Note that this barrier can differ from experimental effective energy barrier. The energy difference between the ground-state and the first-excited doublets for a given $J$-multiplet is referred to as zero-field splitting (ZFS).

For cationic or anionic TbPc2 or TbPcNc SMMs, the Tb $J=6$ multiplet structure constitutes the entire low-energy spectrum. In this case, the transverse CF may split the two-fold degeneracy of nonzero $|M_{J}|$ levels by mixing states with different $M_{J}$. This phenomenon is referred to as tunnel splitting. For neutral TbPc2 SMMs, one unpaired electron (with the spin $s=1/2$) is delocalized (or shared) within the two Pc ligands. This ligand spin interacts with the Tb magnetic moment by exchange coupling ($J_{\text{ex}}$) and doubles the number of low-energy levels. In addition, this extra electron makes neutral TbPc2 a Kramers system, and irrespective of symmetry of the ligand field, it ensures at least two-fold degeneracy of all electronic levels enforced by time-reversal symmetry. Low-energy levels within the ground multiplet, ZFS, tunnel splitting, $J_{\text{ex}}$, and separation between the ground and first-excited multiplets are important energy scales that control the magnetic properties of TbPc2 and TbPcNc SMMs. They play an important role in elucidation of magnetization relaxation mechanisms.

Despite the great interests and the experimental efforts and ambiguity, there are quite few computational studies of TbPc2 molecule or its derivatives. In Ref. 35 a neutral TbPc2 molecule was investigated using density-functional theory (DFT) calculations with and without an on-site Coulomb repulsion $U$ term in the absence of SOI. However, SOI on Tb ion is much stronger than the ligand CF (Fig. 2) and, therefore, it is imperative to include SOI for an even qualitative description of TbPc2. Equally importantly, lanthanides atoms are known to have nearly degenerate electronic configurations demanding multireference treatments. The addition of the $U$ term alone does not suffice to describe the electronic structure and magnetic properties of the TbPc2 molecule even qualitatively. Some multireference calculations of the TbPc2 molecule using complete active space self-consistent field (CASSCF) method including SOI within restricted active space state-interaction (RASSI) have also been reported.36, 37 In particular, in Ref. 36 a neutral TbPc2 molecule was studied using a DFT-optimized structure when the molecule is adsorbed on a Ni substrate. In Ref. 37, ZFS and MAB of anionic TbPc2 were calculated with a particular experimental geometry.29 There are no multireference calculations of the magnetic properties of the TbPcNc molecule. Therefore, there is still a lack of multireference ab initio studies of magnetic energy scales of TbPc2-type SMMs as a function of charge state, type of ligand, and details of molecular geometry.

Here we investigate the effects of ligand and oxidation on the electronic and magnetic properties of the TbPc2 and TbPcNc molecules in a gas phase, by using the CASSCF multireference method including SOI within RASSI. Using experimental geometries, we study the electronic levels characteristics and analyze the dependence of important magnetic energy scales on oxidation, ligand type and details of molecular structure. Furthermore, we construct an effective pseudospin Hamiltonian that includes Tb CF parameters and the Zeeman interaction (as well as exchange coupling between the Tb magnetic moment and the ligand spin for neutral molecules). The paper is structured as follows. Geometries of interest and methods used in this study are described in Sec. 2 and Sec. 3, respectively. Results of the TbPc2 molecule in charged forms are followed by those for neutral forms in Sec. 4. We make conclusions in Sec. 5.

2 Geometries of Study

We perform calculations for TbPc2-type molecules with different charge states and ligand types for which experimental atomic structures are available. The use of experimental molecular geometries is preferred over theoretically optimized geometries since the latter, due to prohibitive computational costs, cannot include counter ions and solvent or dilution molecules which can be quite sizable. Indeed, the type of solvent molecules or diamagnetic dilution molecules as well as crystal packing, can significantly affect the ligand geometry of TbPc2-type molecules. For example, two [TbPc2]- molecules in Refs. 27 and 29 (referred to as M4 and M5 later) have significantly different geometries despite the same oxidation and ligand type (Table 1).

We consider the following six molecules (Table 1): (M1) neutral TbPc2 with experimental geometry from Ref. 32, (M2) neutral TbPc2 with experimental geometry from Ref. 31, (M3) neutral TbPcNc with experimental geometry from Ref. 33, (M4) anionic TbPc2 with experimental geometry from Ref. 27, (M5) anionic TbPc2 with experimental geometry from Ref. 29, (M6) cationic TbPcNc with experimental geometry from Ref. 33. As is a common practice, for all experimental geometries we correct the carbon-hydrogen bond length to 1.09 Å since this distance cannot be reliably extracted from X-ray measurements. Solvent molecules, diamagnetic dilution molecules, or counter ions are not included in our calculations.

Figure 1 shows the atomic geometry of M1, M5 and M3 molecules. The structure can be viewed in terms of the Tb ion sandwiched between two approximately flat Pc or Nc ligand planes. The ligands are rotated with respect to each other by roughly 45∘ angle. Each ligand has four roughly identical branches that form $\approx$90∘ angle with each other. Each branch starts with a nitrogen atom that is the nearest neighbor of the Tb ion and is denoted as N${}_{\text{nn}}$ (Fig. 3). Away from the Tb ion, N${}_{\text{nn}}$ is connected to two intermediate carbon atoms which, in turn, are connected to a benzene-like carbon ring. Two carbon atoms from the ring which are the farthest from the Tb ion are denoted as C${}_{\text{far}}$ (Fig. 3). For Nc, C${}_{\text{far}}$ atoms are further connected to an additional carbon ring. The neighboring branches are linked by bridging nitrogen atoms which form bonds with the intermediate carbon atoms.

In order to characterize the experimental geometries we introduce several structural parameters (see Table 1 and Fig. 3). The analysis of these parameters allows us to quantify deviations of the molecular structure from the ideal $D_{4d}$ (or $C_{4v}$ for TbPcNc) symmetry induced by counter cations/anions, solvent molecules, diamagnetic dilution molecules, or crystal packing. For M1, M3, and M6 molecules, both $\theta_{\text{N}_{\text{nn}}\text{TbN}_{\text{nn}}}$ and $\theta_{\text{C}_{\text{far}}\text{TbC}_{\text{far}}}$ angles are close to 45∘ which indicates that the four-fold symmetry is approximately preserved. To a somewhat lesser degree, this is also true for M2 and M4 molecules. This is consistent with a small standard deviation of the $D_{\text{TbN}_{\text{nn}}}$ and $D_{\text{TbC}_{\text{far}}}$ bond lengths for these molecules. On the other hand, for M5 molecule, we have strong deviations from the four-fold symmetry as reflected in the fact that both the angles and the $\theta_{\text{N}_{\text{nn}}\text{TbN}_{\text{nn}}}$ and $\theta_{\text{C}_{\text{far}}\text{TbC}_{\text{far}}}$ as well as the $D_{\text{TbN}_{\text{nn}}}$ and $D_{\text{TbC}_{\text{far}}}$ bond lengths vary significantly. For all considered geometries, we observe curving of the carbon parts of the ligand planes away from each other (this effect can be also seen in Fig. 1). The curving can be different for different ligand branches and it is most pronounced for the M5 molecule. The strong deviations of the M5 geometry from the ideal TbPc2 structure is likely a result of bulky diamagnetic dilution molecules used in the synthesis process.29

3 Methods

The multireference calculations are performed using the Molcas quantum chemistry code (version 8.2).38 Scalar relativistic effects are included based on the Douglas-Kroll-Hess Hamiltonian39, 40 using relativistically contracted atomic natural orbital (ANO-RCC) basis sets.41, 42 In particular, polarized valence triple-$\zeta$ quality (ANO-RCC-VTZP) is used for the Tb ion, polarized valence double-$\zeta$ quality (ANO-RCC-VDZP) is used for the nitrogen and carbon atoms, and valence double-$\zeta$ quality (ANO-RCC-VDZ) is used for the hydrogen atoms. Such a choice of the basis set is made to maintain a high accuracy and to not exceed computational capabilities. More details on the basis set dependence are discussed in Tables S4 and S5 in Supporting Information.

First, in the absence of SOI, for a given spin multiplicity, the spin-free eigenstates are obtained using state-averaged CASSCF method.43, 44 The valence electronic configurations of Tb*+3* ion consists of eight electrons at 4$f$ orbitals which must be included in the active space. It would be desirable to include all ligand $\pi$/$\pi^{\star}$-type orbitals as well. This is, however, computationally prohibitive and, therefore, only several near-in-energy ligand orbitals are included in the active space. In order to identify such ligand orbitals, we consider molecules in the cationic state and perform CASSCF calculations with eight electrons and seven Tb 4$f$-type orbitals in the active space. We find that HOMO (L1) is always energetically separated ($\sim$0.1 a.u.) from other occupied ligand orbitals. For unoccupied orbitals we find that LUMO-1 (L2) as well as nearly degenerate LUMO-2 (L3) and LUMO-3 (L4) are separated from higher unoccupied states. Therefore, we include L1, L2, L3, and L4 ligand orbitals in the active space. Altogether, as illustrated in Fig. 4, we consider eleven active orbitals with ten, eleven, and twelve active electrons for cationic, neutral, and anionic molecules, respectively. The effect of the choice of the active space on the final results is discussed in Supporting Information (Tables S7 and S8).

In the case of charged molecules we consider only the $S=3$ configuration since other spin configurations such as $S=2$, $S=1$, and $S=0$ are much higher in energy and their inclusion changes the energy levels of the ground $J$-multiplet only by about a few cm*-1*. See Table S6 in Supporting Information for details. For neutral molecules, depending on whether the ligand spin ($s=1/2$) is parallel or antiparallel to the Tb spin $S=3$, we have two possible values of the total spin of the molecule: $S_{\text{tot}}=7/2$ or $S_{\text{tot}}=5/2$ (see Fig. 4). We consider both $S_{\text{tot}}$ values since they lie close in energy. For a given spin configuration, we evaluate seven lowest spin-free states (roots) that correspond to different configurations of eight electrons in Tb 4$f$-type orbitals with Tb spin of $S=3$. These seven spin-free states are used in the state-averaged procedure.

In the next step, we include SOI, within the atomic mean-field approximation,45 in the space of aforementioned spin configurations and spin-free eigenstates, using the RASSI method.46 With SOI, all possible $J$-multiplets from the addition of ${\mathbf{L}}$ and ${\mathbf{S}}$ are generated as illustrated in Fig. 2. For the calculation of the CF parameters, we use the methodology implemented in the SINGLE_ANISO47 module of the Molcas code.

The technique from Ref. 47 that we use for the charged molecules cannot be directly applied to the neutral molecules. Indeed, for lanthanides this method finds CF parameters and $g$-tensor elements for a given $J$-multiplet using lowest $(2J+1)$ ab initio eigenvalues and the corresponding eigenfunctions. For neutral molecules, however, the low-energy levels do not only originate solely from the ground $J$-multiplet but also they involve the unpaired ligand electron spin-flip states. Both types of excitations are entangled and there is no obvious way how to extract multiplet levels and their wave functions. In fact, the entanglement is essential for formation of Kramers doublet. In order to circumvent this problem, we take the experimental geometries of the neutral molecules and consider them in the cationic state for the calculation of the Tb CF parameters. In this way we remove the unpaired ligand electron so that the calculated low-energy spectrum can be put into correspondence with the Tb $J=6$ multiplet and the CF parameters as well as the $g$-tensor elements can be calculated. Assuming that the unpaired ligand electron has a small contribution to the Tb CF and $g$-tensor, these parameters as well as the exchange coupling constant can be used in effective pseudospin Hamiltonian for the neutral molecules.

4 Results and Discussion

In each subsection, we present the calculated magnetic energy levels obtained from CASSCF-RASSI-SOI method and construct the effective pseudospin Hamiltonian with the calculated CF parameters and $g$ tensor. We then compare with relevant theoretical results and experimental data when they are available.

4.1 Charged molecules

For charged TbPc2-type molecules (M4, M5 and M6), we find that among eleven active molecular orbitals, seven are Tb 4$f$-type. These orbitals are similar for all considered molecules and are shown in Fig. S1 in Supporting Information. The occupation of each 4$f$-type orbital is approximately 1.14 which corresponds to about eight electrons occupying these orbitals (note that this is a result of state-averaged calculations). The remaining four active orbitals (L1-L4) are ligand orbitals and are shown in Figs. S5-S7 in Supporting Information. For cationic M6 molecule, L1 is almost doubly occupied while L2-L4 are almost empty. All ligand orbitals arise mostly from inner C atoms and they have in-plane symmetry due to $C_{4}$ symmetry of the molecule (Table 1). There is, however, asymmetry between orbitals in the top Pc and the bottom Nc ligands. For anionic M4 molecule, the occupation of nominally filled L1 and L2 orbitals is somewhat lower than two while the nominally empty L3 and L4 orbitals have a significant occupation (Fig. S5). The effect is even stronger for the M5 molecule (Fig. S6) and it suggests sizable correlations. The orbitals have significant in-plane asymmetry as well as asymmetry between the two Pc planes. This is a consequence of deviations from ideal $D_{4d}$ symmetry as illustrated by large spread and uneven distribution of $Z_{\text{C}^{\text{up}}}$ and $Z_{\text{C}^{\text{down}}}$ structural parameters (Table 1).

For the charged molecules, the low-energy spectrum is largely determined by Tb atomic levels which are split in the ligand CF. As listed in Table 2, the separation between the ground and the first-excited $J$-multiplet, $\Delta E_{J}$, is almost a factor of three larger than the MAB of the ground $J$-multiplet, $E_{\text{MAB}}$. Thus, we can project the lowest $J$ multiplet onto a multiplet of effective pseudospin $\tilde{J}=6$. Figure 5 shows the thirteen lowest calculated energy levels corresponding to the $J=6$ ground multiplet for all considered charged molecules (numeric data together with spin-free energies are provided in Tables S1 and S2 in Supporting Information). In order to elucidate the properties of these levels we construct an effective pseudospin Hamiltonian

[TABLE]

where $\hat{O}_{k}^{q}(\hat{\mathbf{J}})$ are $k$-th rank extended Stevens operators48 and $B_{k}^{q}$ are corresponding CF parameters ($q=-k,...,0,...,k$). Here $\hat{\mathbf{J}}$ is the Tb pseudospin operator, and $\hat{O}_{k}^{-q}(\hat{\mathbf{J}})=[\hat{O}_{k}^{q}(\hat{\mathbf{J}})]^{\star}$, where $\star$ means complex conjugate. For example, second-rank Stevens operators are $\hat{O}_{2}^{0}=(3J_{z}^{2}-J(J+1)I)$, $\hat{O}_{2}^{1}=(J_{z}(J_{+}+J_{-})+(J_{+}+J_{-})J_{z})/4$ and $\hat{O}_{2}^{2}=(J_{+}^{2}+J_{-}^{2})/2$, where $J_{\pm}$ are raising and lowering operators and $I$ is a 2$\times$2 identity matrix. Higher-rank Stevens operators are listed in Table S9 in Supporting Information. Here $\hat{O}_{k}^{q=0}$ terms represent uniaxial or diagonal contributions, while $\hat{O}_{k}^{q\neq 0}$ terms are transverse or off-diagonal contributions. Time-reversal symmetry enforces only even integer $k$ in the summation. Molecular symmetry dictates allowed nonzero $B_{k}^{q}$ values. When the second-order uniaxial term, $\hat{O}_{2}^{0}$, is dominant, ZFS and MAB are approximated to be $3|B_{2}^{0}|(2J-1)$ and $3|B_{2}^{0}|J^{2}$, respectively, for integer $J$. The second term in Eq. (1) is the Zeeman interaction with the $g$-tensor, $\mathbf{g}$.

We evaluate the elements of the $g$-tensor and $B_{k}^{q}$ values using the thirteen ab initio energy levels and the corresponding eigenfunctions.47, 49 The results are shown in Tables 3 and 4 (higher order CF parameters are shown in Table S10 in Supporting Information). We use the coordinate system along the principal axes of the $g$-tensor. In this coordinate system, the $z$ axis is roughly perpendicular to the ligand planes. Note that the principal values of the $g$-tensor are close to the ideal Lande $g$-factor value of $3/2$. The calculated $|B_{2}^{0}|$ value is the largest for M5 and the smallest for M4. Due to the absence of any symmetry, however, M4 and M5 molecules have significant all-order off-diagonal CF parameters (Table 4). In particular, for M5, the $B_{2}^{2}$ value is about 63% of the $|B_{2}^{0}|$ value, while the $|B_{4}^{4}|$ value is similar to the $|B_{4}^{0}|$ value. On the other hand, for M6, only off-diagonal terms with $q={\pm 4}$ are significant, which is dictated by the molecular $C_{4}$ symmetry. Including the full set of CF parameters up to sixth order ($k=6$) and the diagonal eighth order term ($B_{8}^{0}$), diagonalization of Eq. (1) reproduces the ab initio energy levels up to 0.5 cm*-1*. This indicates that Eq. (1) is a proper Hamiltonian for the low-energy spectrum of charged TbPc$2$-type SMMs.

The lower part of the spectrum (six lowest levels for M4 and M6 and four lowest levels for M5) is composed of approximate doublets, as shown in Fig. 5. The former three doublets correspond to states $M_{J}=\pm 6$, $\pm$5, and $\pm$4, while the latter two doublets are states $M_{J}=\pm 6$ and $\pm$5, respectively. For further level characteristics, see Table S14-S16 in Supporting Information. The transverse CF, however, mixes states with different $|M_{J}|$ values and breaks the degeneracies of the doublets leading to tunnel splitting. It is important to understand tunneling splitting of the levels within the ground multiplet since magnetization can be relaxed via phonon-assisted tunneling. For low-energy levels with large $|M_{J}|$, the transverse CF has a small effect since large powers of its matrix elements are needed to connect states with opposite $M_{J}$. For high-energy levels with small $|M_{J}|$, the transverse CF becomes more important and tunnel splitting is more pronounced. In fact, for higher part of the spectrum, $|M_{J}|$ ceases to be a good quantum number and the doublet structure disappears. The tunnel splitting and $|M_{J}|$ mixing for M5 are significantly larger than those for M4 and M6. This is because the substantial distortions and curving of the ligand planes for M5 lead to significant transverse CF (Table 4). The geometrical distortions are shown in the larger spread of $D_{\text{TbN}_{\text{nn}}}$ and $D_{\text{TbC}_{\text{far}}}$ values and of $\theta_{\text{N}_{\text{nn}}\text{TbN}_{\text{nn}}}$ and $\theta_{\text{C}_{\text{far}}\text{TbC}_{\text{far}}}$ values and in the larger range of the $Z_{\text{C}^{\text{up}}}$ and $Z_{\text{C}^{\text{down}}}$ values in Table 1, compared to the other molecules of interest. For M5, the tunnel splitting values are of the order of $10^{-2}$, $10^{-1}$, and 10 cm*-1* for the ground state and the first- and second-excited states, respectively (Table 2). For M6, the $C_{4}$ symmetry allows mixing of $M_{J}$ levels with \Delta M_{J}$$=$$\pm 4 only. This explains the significant tunnel splitting only for a pair of states 5 and 6 and a pair of states 9 and 13 (Table 2 and Table S2 in Supporting Information).

The calculated ZFS and MAB values are shown in Table 2. We find that ZFS lies in the range of 289-292 cm*-1* for the anionic molecules (M4 and M5), while it is somewhat larger for the cationic molecule (M6). As seen from Fig. 5, ZFS for M6 is similar to ZFS for other cationic TbPc2 or TbPcNc molecules. Therefore, we conclude that ZFS does not depend on ligand type and geometry details, and it only shows a weak dependence on oxidation number. This reflects the fact that ZFS characterizes the lower part of the energy spectrum where the transverse CF has a small effect. However, higher-energy part of the spectrum where the transverse CF plays a significant role reveals a much stronger dependence on the ligand type and geometry details, as shown in Fig. 5.

Our calculated results for M5 molecule can be directly compared with similar calculations from Ref. 37 where the same experimental geometry was used. ZFS and MAB reported in Ref. 37 are 308 cm*-1* and 809 cm*-1*, respectively. While ZFS is reasonably close to our value (292 cm*-1*), the MAB differs more significantly from our value (770 cm*-1*). We check that both the choice of the active space and the effect of higher spin-free states are not responsible for the significant difference (see Supporting Information). The most likely reason is the basis set difference between our study and Ref. 37.

On the experimental side, the Tb CF parameters for [TbPc2]-TBA+ 24 estimated by Ishikawa et al.50 are widely used in the community. These CF parameters were obtained by fitting both the experimental nuclear magnetic resonance (NMR) shift and magnetic susceptibility data to a ligand-field model, assuming perfect C4 symmetry. The estimated CF parameters are listed in Thiele et al.14. In this estimate, three uniaxial terms $B_{2}^{0}$, $B_{4}^{0}$, and $B_{6}^{0}$ as well as only one transverse term $B_{4}^{4}$ were considered. The estimated values are B_{2}^{0}$$=$$-4.18193, B_{4}^{0}$$=$$-0.02790, $B_{4}^{4}=$0.00122, and B_{6}^{0}$$=$$-0.00004 cm-1. Since X-ray crystallography data was not reported and perfect C4 symmetry was assumed in Ishikawa et al.,50 our calculated CF parameters cannot be directly compared to their fitted values. Nonetheless, it is worthwhile providing some remarks. Note that M4 molecule has the same counter cation without bulky dilution molecules as in Ishikawa et al.50. Thus, we discuss major differences between our calculated CF parameters for M4 molecule with those by Ishikawa et al.50. For M4 molecule, the $|B_{2}^{0}|$ value is somewhat greater than the fitted value, while the $|B_{4}^{0}|$, $|B_{4}^{4}|$, and $|B_{6}^{0}|$ values are comparable to the fitted values. See Table 6 and Table S10 in Supporting Information. The major difference is that we find large low-order transverse CF parameters such as $|B_{2}^{q\neq 0,2}|$, $|B_{4}^{q\neq 0,4}|$ and $|B_{6}^{q\neq 0}|$, whereas the literature considered only one transverse CF parameter $B_{4}^{4}$. As a result, our result predicts a much larger tunnel splitting and qualitatively different level characteristics above the first-excited states.

Table 5 shows the experimental effective barrier for several TbPc2-like molecules. The measured barrier for M6 was reported, considering Raman and Orbach processes, while the barrier for a diluted crystal of neutral TbPcNc (M3) molecules was obtained, considering Orbach and quantum tunneling processes.33 It is interesting to compare the measured barrier for M5 with and without bulky diamagnetic dilution molecules.29 These experimental values are much larger than the earlier reported values from Refs. 24 and 25, for anionic TbPc2 molecular crystals. Comparison with the experimental data indicates that for M6 molecule, the experimental barrier is close to states 11 and 12, and for M5 molecule, the range of the experimental barrier falls on states 10 and 11. See Table S2 in Supporting Information.

4.2 Neutral molecules

Neutral TbPc2-type molecules (M1, M2 and M3) have much richer electronic structure than their charged counterparts due to presence of an extra unpaired electron that is expected to reside in the ligands. As in the case of the charged molecules, seven of the eleven active orbitals are Tb 4$f$-like (Fig. S1 in Supporting Information). These orbitals are similar for all considered molecules and occupation of each of these orbitals is roughly 1.14 consistent with eight electrons occupying the 4$f$ shell. The remaining four active orbitals (L1-L4) are ligand orbitals and are shown in Figs. S2-S4 in Supporting Information. For all neutral molecules L1 is almost doubly occupied, L2 is singly occupied (SOMO) while L3 and L4 have small occupation numbers. For M3 molecule, the SOMO (L2) has top-down plane asymmetry due to different ligand type but shows in-plane symmetry. Interestingly, for M1 and M2 molecules, despite the same ligand type, the SOMO is highly asymmetric with large weight on the bottom Pc plane only (Fig. S2 and S3). This is due to the asymmetric curvature within the two Pc planes, as listed in Table 1. Compare the $Z_{\text{C}^{\text{up}}}$ and $Z_{\text{C}^{\text{down}}}$ values in Table 1. In order to check on this, we enforce mirror symmetry about the $xy$ plane ($\sigma_{h}$) in M1 molecule and compute the ligand orbitals in the active space. For this flattened M1 molecule, we find that top-bottom plane symmetry is more or less restored in the SOMO. Compare Fig. 6(a) with (b).

The spin of the ligand electron can be either parallel ($S_{\text{tot}}=7/2$) or antiparallel ($S_{\text{tot}}=5/2$) to the Tb spin. The calculated spin density for both cases is shown in Fig. 7 for M1 and M3 molecules. The Tb spin density is localized in the vicinity of Tb ionic core. For M1, the ligand spin density is mostly shared by inner carbon atoms in one side of the bottom Pc plane. This is consistent with the SOMO (Fig. 6(a)) that is primarily delocalized on these carbon atoms. The same behavior is observed for M2 molecule (not shown). For M3, consistently with SOMO, the majority of the ligand spin density is symmetrically shared by the inner carbon atoms from the Nc ligand. Table S1 in Supporting Information shows the calculated seven lowest spin-free states for both values of $S_{\text{tot}}$ for different neutral molecules. For all systems, the parallel configuration of the Tb and ligand spin ($S_{\text{tot}}=7/2$) has a lower energy than the antiparallel configuration ($S_{\text{tot}}=5/2$).

Since neutral TbPc2 and TbPcNc molecules have an odd number of electrons, the Kramers theorem dictates that all electronic levels are at least doubly degenerate. The low symmetry of the considered molecules prevents from appearance of higher degeneracy. The low-energy part of the calculated spectra for M1, M2, and M3 molecules consists of a group of thirteen Kramers doublets (Fig. 8 and Table S3 in Supporting Information).

The thirteen lowest Kramers doublets result from Tb $J=6$ multiplet that is coupled with the unpaired ligand electron spin. The next electronic level belongs to the first-excited multiplet and it lies more than 1425 cm*-1* higher than the highest level in the first multiplet (Table 2). Therefore, for analysis of the magnetic properties we can focus only on the thirteen lowest doublets. Such an energy spectrum can be described by the following generalization of Eq. (1) in order to include the ligand spin:

[TABLE]

where we introduce Heisenberg exchange coupling $J_{\text{ex}}$ between the ligand spin ($\hat{\mathbf{s}}$) and the Tb total angular momentum operators. In addition, the Zeeman term is generalized to describe interaction of the ligand spin with magnetic field ($g_{e}\approx 2$ is a free-electron $g$-factor). Note that in the presence of strong SOI, anisotropic and antisymmetric exchange couplings may play an important role and rigorous treatment of exchange interaction requires more general formalism.52, 53, 54 We expect, however, that the simple isotropic Heisenberg form in Eq. (2) is a good approximation for our molecules since the exchange interaction is very small. The anisotropic and antisymetric exchanges (as well as higher order spin interactions) are, in general, smaller than the isotropic exchange and, therefore should not play a significant role for TbPc2-type molecules. As discussed below, we find that Hamiltonian (2) provides a good representation of the low-energy spectrum of the neutral molecules.

First-principles evaluations of the Tb $B_{k}^{q}$ coefficients and the $g$-tensor elements are not straightforward for the neutral molecules. As discussed in Sec. 3, we obtain the Tb $B_{k}^{q}$ coefficients and the $g$-tensor by using the calculated low-energy spectra (Fig. 5(a)-(c) and Table S2 in Supporting Information) of the neutral molecular geometries with one electron removed. Here we assume that the contribution of the unpaired ligand electron to the Tb CF and $g$-tensor is negligible. Later we show that this is indeed a valid assumption. The calculated elements of the $g$-tensor and the Tb CF parameters are shown in Tables  3 and 6 (see Table S10 in Supporting Information for a full set of $B_{k}^{q}$). Note that as in the case of M4, M5 and M6 molecules, we use a coordinate system along the principal axes of the $g$-tensor for which the $z$ axis (easy axis) points approximately in the direction perpendicular to the ligand planes. We find all the principal values of the $g$-tensor being close to the ideal Lande $g$-factor value of $3/2$ (Table 3). The calculated $|B_{2}^{0}|$ value for M1 is a bit larger than that for M2 and M3, which is reflected in the largest ZFS among the three neutral molecules. The calculated $|B_{2}^{0}|$ value for M1, M2, and M3 is consistently larger than that for the charged molecules. Compared to the charged TbPc2 molecules, the degree of the structural distortion is much less for the neutral TbPc2 molecules. See $D_{\text{TbC}_{\text{far}}}$, $Z_{\text{C}^{\text{up}}}$, and $Z_{\text{C}^{\text{down}}}$ values in Table 1. This explains overall smaller transverse CF parameters in the neutral TbPc2 molecules than in the charged TbPc2 molecules.

In order to evaluate the exchange coupling $J_{\text{ex}}$, we fit the ab-initio energies of the lowest thirteen doublets to eigenvalues of Eq. (2) with a fitting parameter $J_{\text{ex}}$. For all molecules, we obtain a high-quality fit. This indicates that the method we use for evaluation of CF parameters is reliable and that Eq. (2) provides a reasonable description of the low-energy spectrum of neutral TbPc2-type molecules. We find a small ferromagnetic exchange coupling with $J_{\text{ex}}\approx 0.6-0.8$ cm*-1* (Table 2). An increase of the active space or the atomic basis sets does not affect this number significantly (see Tables S4 and S7 in Supporting Information). The insensitivity of $J_{\text{ex}}$ to ligand type and synthesis process indicates that the exchange coupling does not depend on details of molecular geometry much as long as structural changes are moderate. This conclusion is only applied to molecules in crystals or films on substrates rather than within single-molecule transistors where the molecules experience much stronger structural changes.

Using the pseudospin Hamiltonian, Eq. (2), we can determine the characteristics of energy levels shown in Fig. 8. For all considered molecules, four lowest Kramers doublets have a well-defined $|M_{J}|$ value. For each pair, the lower (higher) energy level has the ligand spin parallel (antiparallel) to the Tb angular momentum. Within each Kramers doublet, it is convenient to choose such a basis set in which each Kramers partner state is characterized by $M_{J}$ and $m_{s}$ quantum numbers. Since Kramers partner states are related by time reversal symmetry, they have opposite values of $M_{J}$ and $m_{s}$. If we focus on the Kramers partner state with positive $M_{J}$, we find ($M_{J}=6$, $m_{s}=1/2$) for the ground state level (level 1), ($M_{J}=6$, $m_{s}=-1/2$) for the first-excited level (level 2), ($M_{J}=5$, $m_{s}=1/2$) for the second-excited level (level 3), ($M_{J}=5$, $m_{s}=-1/2$) for the third-excited level (level 4) (Fig. 8). For the most symmetric M3 molecule, level 5 consists of a majority of ($M_{J}=4$, $m_{s}=1/2$) and small contributions from ($M_{J}=-4$, $m_{s}=1/2$) and ($M_{J}=0$, $m_{s}=1/2$), while level 6 consists of a majority of ($M_{J}=4$, $m_{s}=-1/2$) with small contributions from ($M_{J}=-4$, $m_{s}=-1/2$), ($M_{J}=-4$, $m_{s}=1/2$), ($M_{J}=0$, $m_{s}=-1/2$) and ($M_{J}=3$, $m_{s}=1/2$).

For the neutral molecules, we define ZFS as an energy difference between $|M_{J}|=6$ and $|M_{J}|=5$ levels with the same direction of ligand spin with respect to the Tb total angular momentum (values for parallel and antiparallel configurations are very close). As seen in Table 2 and Fig. 8, ZFS is similar for M1, M2 and M3 molecules. This indicates that ZFS is not sensitive to the ligand type and geometry details, similarly to $J_{\text{ex}}$. However, the ZFS for the neutral molecules is somewhat larger than the values for the anionic and cationic molecules. Higher-energy part of the spectrum shows some dependence on ligand type and geometry details but we have less sensitivity than in the case of charged molecules.

We now compare our calculated CASSCF-RASSI-SO results to other theoretical calculations and experimental data. Regarding $J_{\text{ex}}$, a similar value to our value has been obtained from CASSCF-RASSI-SO calculations for a neutral TbPc2 molecule adsorbed on a Ni substrate using a DFT-optimized atomic structure.36 Compared to the experimental $J_{\text{ex}}$ from electron paramagnetic resonance spectra for a crystal of M1 molecules,32 the sign of our calculated $J_{\text{ex}}$ agrees with experiment but the magnitude is a bit larger than the experimental value. Experiments on TbPc2-based single-molecule transistors have shown both antiferromagnetic 14 and ferromagnetic 55 couplings between the ligand spin and the Tb spin. These seemingly conflicting experimental results are not surprising, considering the small energy scale of $J_{\text{ex}}$ and possible large configurational changes of the ligand planes of the TbPc2 molecule bridged between gold electrodes in single-molecule transistor set-ups. Experimental data for the effective energy barrier is rare for neutral TbPc2 molecules. There exists an experimental report on a neutral TbPcNc molecular crystal diluted with YPcNc.33 In this case, the experimental barrier is 342 cm*-1* (Table 5). The measured value seems to be close to our calculated ZFS for M3 molecule. See Tables 2 and 6.

5 Conclusions

We investigate electronic structure and magnetic properties of six different TbPc2 and TbPcNc molecules in different charge states using relativistic multireference methods and pseudospin Hamiltonian technique. For the charged molecules, we evaluate CF parameters and $g$-tensor elements by projecting the low-energy spectrum onto effective $J=6$ pseudospin. For the neutral molecules, we consider exchange coupling of the Tb magnetic moment with the ligand spin and extract Tb CF parameters and $g$-tensor elements by separating the Tb ground multiplet from the unpaired ligand electron spin-flip states using artificial oxidation. The key findings are as follows:

• •

For the neutral molecules, the exchange coupling constant between the Tb magnetic moment and the ligand spin does not depend much on ligand type and geometry details. This result is valid as long as molecular structures are more or less controlled such as in crystals or layers on substrates.

• •

Geometry details and ligand type do not affect ZFS.

• •

ZFS weakly depends on oxidation number. The neutral molecules have somewhat higher ZFS than the charged molecules.

• •

The higher-energy levels and associated tunnel splittings strongly depend on ligand type, oxidation number, and geometry details.

• •

Comparison to experimental effective barrier suggests that in some cases higher-energy levels rather than just ZFS may contribute to the magnetization relaxation through phonon-assisted tunneling.

These results provide insights in separating the effects intrinsic to individual molecules from extrinsic effects on magnetization relaxation and in interpretation of reported experimental data and stimulating new experiments on TbPc2-type molecules.

{acknowledgement}

This work was funded by the Department of Energy (DOE) Basic Energy Sciences (BES) grant No DE-SC0018326. Computational support by Virginia Tech ARC and San Diego Supercomputer Center (SDSC) under DMR060009N. We also thank Dr. Benjamin Pritchard for helpful discussion and insight.

{suppinfo}

The Supporting Information is available free of charge: Molecular orbitals, spin-free energies, numeric data for energy levels, basis set and active space dependence, definitions of extended Stevens operators, full set of crystal field parameters, and character of electronic levels.