Ab Initio Design of Molecular Qubits with Electric Field Control

Current scalable quantum computers require large footprints and complex interconnections due to the design of superconducting qubits. While this architecture is competitive, molecular qubits offer a promising alternative due to their atomic scale and tuneable properties through chemical design. The use of electric fields to precisely, selectively and coherently manipulate molecular spins with resonant pulses has the potential to solve the experimental limitations of current molecular spin manipulation techniques such as electron paramagnetic resonance (EPR) spectroscopy. EPR can only address a macroscopic ensemble of molecules, defeating the inherent benefits of molecule-based quantum information. Hence, numerous experiments have been performed using EPR in combination with electric fields to demonstrate coherent spin manipulation. In this work, we explore the underlying theory of spin-electric coupling in lanthanide molecules, and outline ab initio methods to design molecules with enhanced electric field responses. We show how structural distortions arising from electric fields generate coupling elements in the crystal field Hamiltonian within a Kramers doublet ground state and demonstrate the impact of molecular geometry on this phenomenon. We use perturbation theory to rationalize the magnetic and electric field orientation dependence of the spin-electric coupling. We use pseudo-symmetry point groups to decompose molecular distortions to understand the role that symmetry has on spin-electric coupling. Finally, we present an analytical electric field model of structural perturbations that provides large savings in computational expense and allows for the investigation of experimentally accessible electric field magnitudes which cannot be accessed using common ab initio methods.


S1 Analytic Model
Due to limitations in both computational capability and cost the need for an analytical model of electric field distortions is required to aid in high throughput calculations and in calculations for electric field strengths smaller than the smallest available electric field ab initio (1 × 10 8 V m −1 ) to replicate experimental results.
We first define the potential energy of the system to be given as: Where ⃗ r is the atomic coordinate vector containing 3N entries, H is the Hessian matrix obtained from a frequency calculation at the equilibrium geometry and ⃗ g is the potential energy gradient vector which can be assumed to be zero at the equilibrium geometry.
The force associated with the displacement of atoms is given as the negative gradient of the potential energy and is given by: To compute the gradient of the potential energy with respect to atomic coordinate displacements, equation S1 is written in terms of its components.
Where terms the same indices are taken outside of the summation.Taking the derivative the potential with respect to a general element k we get.

∂V ∂r
Adding the k th element of the derivative back into the summations we get.

∂V ∂r
Writing equation (S6) back into matrix form gives the following equations.
Where l is the number of modes.The Hessian is a symmetric matrix (H T = H) such that the gradient of the potential energy can be written as H • r as such the force is given by.
To calculate the distortion due to an applied electric field, we calcualate the potential energy due to the applied electric field and the associated force.The potential energy due to an applied electric field can be written as a power series expansion of the molecular potential energy as a function of electric field.
Where V 0 is the molecular potential energy at zero field, µ is the dipole moment, ⃗ α is the polarisability, ⃗ β is the hyperpolarisability, and ⃗ E is the electric field vector.Here we truncate the Taylor series expansion of the potential energy to fist order as the term that depends on the electric field squared has negligible contribution to structural distortions.
Hence, we take the potential energy due to an electric field to be: The force exerted of the atomic coordinated due to an electric field, F E , is therefor given as the negative of gradient of the electric dipole moment with respect the the displacement of atomic coordinates.
Alternatively this can be written in terms of its components such that the α th component of the force exerted on atom i due to an applied electric field is given by: The atomic electric dipole derivatives are obtained from directly from the initial geometry optimisation maintaining the the single shot approach of this method.The force due to electric field and the force due to atomic displacements are set equal to each other giving the following total equation for the distortion due to an electric field.
This forms a set of exact linear equations for ⃗ r.However, the the Hessian still contains the energy invariant rigid body translations and rotations (3 and 6 for the solid state and gas phase respectively) which should not be included when finding a solution.We Where the solution to the set of linear equations is given by U T r, where U T ⃗ r is the projection of ⃗ r in to the column space of U .Such that our full rank solution is projected back into the atomic displacement vector space.
Where U U T = 1 as U is an orthonormal basis.This yields the distortion vector r.
The units of quantities extracted from Gaussian 16 can be somewhat ambiguous and therefore we provide a table of extracted units and the conversions we perform form the analytical model.All unit conversions are performed using the physical constants library Table S1: All unit conversions from the extracted Gaussian16 quantities to the units used in the analytical model.

Quantity Gaussian16 unit Model unit
as part of the scipy package for python.
The calculated displacements are then converted from m to Å to be added to the equilibrium coordinates extracted from the formatted checkpoint file.
Our electric field model can be used to obtain the derivatives of the crystal field Hamiltonian with respect to an applied electric field ∂H Using these analytical derivatives we can evaluate the electric field Hamiltonian for any given electric field without the requirement to evaluate the geometric distortion due to an applied electric field first.

S2 Perturbation Theory Derivation
To rationalise the electric and magnetic field orientation dependence observed we used perturbation theory to approximate the dependencies of the coupling elements in the electric field Hamiltonian, H ⃗ E .The crystal field Hamiltonian is comprised of degenerate Kramers doublet states which are even under time-reversal symmetry.Hence, structural perturbations alone will not allow for intra-Kramers doublet mixing.Therefore, the degeneracy of the doublet states must be lifted by a small perturbation in the form of a Zeeman term to allow intra-kramers doublet mixing.The ab initio crystal field parameters were projected into a model Hamiltonian and diagonalised into its eigenstate basis, H 0 .Where the eigenvalues are doubly degenerate such that Degenerate perturbation theory was used to express the ground doublet |1⟩, | 1⟩ as the states |1 ± ⟩ such that there is an angle, θ, between the basis vectors, they remain orthonormal is this basis, and they remain a time-reversal conjugate doublet.The transformed basis vectors are expressed in equations (S20) and (S21).
such that energy of states |1 ± ⟩ can be expressed as where ∆E is the splitting generate from the small Zeeman perturbation to the basis states as given in equation S22.
Non-degenerate perturbation theory was subsequently used to derive the corrections the states and energy from a small Zeeman term added to the total Hamiltonian.Such that in this basis the Schrodinger equation is expressed by equation (S23).Where state |n⟩ is an eigenstate of the Hamiltonian.
To obtain the first order correction to the energy, equation ( S23) is acted on by the bra, ⟨n 0 |, which is the zeroth order correction the state |n⟩.Acting with bra, ⟨n 0 |, separating the terms of the total Hamiltonian, and expanding |n⟩ yields.
Given that we are only going to compute the first order corrections for this analysis the terms of the equations of higher order corrections are neglected.Evaluating the expression gives.
As expected the correction to the energy of a state is approximately the expectation value of the Zeeman Hamiltonian in that state.To compute the first order correction to the states, equation ( S23) is acted on by the bra, ⟨m 0 |.Where ⟨m 0 | is the zeroth order correction of an excited state such that m ̸ = 1, 1. Again, the terms of the Hamiltionian are separated and expanded giving.
Given that H 0 is in its eigenstate basis off diagonal matrix elements are zero.Therefore, the expression above can be simplified and then rearranged to give an expression for the correction to the state.
Using the the fact: 1 = m̸ =1, 1 |m 0 ⟩ ⟨m 0 |, The first order correction to the state can be expressed by equation ( S30) Using the first order correction to the state and energy, The states |1 ± ⟩ can be expressed in a new basis |1 ′ ± ⟩ such that the states and energies are corrected for the perturbation as given in equation S32.
The correction to the state is dependent on the mixing of the ground doublet states, |1 ± ⟩ with excited states divided by the difference in energy of those states.To compute the relative orientation dependence of the magnetic and electric fields, we project the total perturbed Hamiltonian, H ′ onto the corrected ground doublet derived using perturbation theory.For simpler notation, we define operator, Q, that has the form Em−En .This allows us to write equation (S32) as equation ( S34) The projection operator, P , projects the total perturbed Hamiltonian, H dist onto the states |1 ′ ± ⟩ and is given by equation (S35).
Projecting and evaluating the known terms gives equation (S36).
The off-diagonal spin-electric coupling terms are only found in the projection of the electric field Hamiltonian onto the ground doublet.Therefore, the decomposition of the spin-electric coupling can be found by expressing the off-diagonal matrix elements in the unperturbed basis using the equation (S34).
Evaluating inside the braket gives (noting that on the left hand side the basis states are primed and they are non-primed on the right hand side) Where A can be expanded using equation (S34) to reform the full expression for the off-diagonal matrix elements.O(H 2 Zee ) is omitted as its second order in the Zeeman term.
Combining the summation gives the perturbative expression seen in the main text.

Magnetic and Electric Field Orientation Dependence of the Spin-Electric Coupling
To understand the magnetic field orientation dependence of the spin-electric coupling the magnetic anisotropy (Figure S2) of excited Kramers doublets need to be considered as perturbation theory shows that the spin-electric coupling between the ground Kramers doublet is dependent on the mixing of the ground Kramers doublets with excited Kramers doublets.The magnetic anisotropy of excited Kramers doublets becomes increasingly easy-axis along x for as the doublet energy increases (Table S2).The Zeeman Hamiltonian's for each magnetic field orientation (Figure S7b) shows that when the magnetic field is orientated along the main molecular axis x, there are intra-Kramers doublet mixing which break the degeneracy of the Kramers doublets and near-zero inter-Kramers doublet mixing.For excited Kramers doublets that have easy-axis anisotropy we can chose that the states |m⟩ and | m⟩ each correspond to |J x ; m J ⟩ states such that they are eigenstates of J x which have the freedom to rotate within Kramers doublet forming linear combinations and hence there is only off diagonal matrix elements within each Kramers doublets.On the other hand, when the magnetic field is orientated in the easy-plane (y, z) we observe only inter-Kramers doublet mixing as the excited states are not eigenstates of J y or J z and we see off diagonal matrix elements between the the ground Kramers doublet and Kramers doublets with |J x ; m J ⟩ representations of |J x ; m J ± 1⟩.This coincides with the magnetic field orientation dependence seen in ab initio, as the perturbation theory shows that the spin-electric coupling is dependent on the inter-Kramers doublet mixing.As there is only intra-Kramers doublet mixing when the magnetic field is orientated along x this explains the smallest spin-electric coupling.
Perturbation theory can also be used to rationalise the electric field orientation dependence.The electric field Hamiltonian in the eigenbasis of the equilibrium geometry crystal field Hamiltonian (Figure S7a) for each electric field orientation all show zero intra-Kramers doublet mixing as the doublets in this basis are degenerate and are related by time-reversal conjugate symmetry.In the case where the electric field is orientated along (pseudo)-high-symmetry axis (y) where the distortion is symmetry preserving the electric field Hamiltonian is near-diagonal such that this distortion purely shifts the energies of the crystal field eigenstates.However, when the electric field is orientated perpendicular to the (pseudo)-high-symmetry axis where the distortion are always symmetry breaking there are inter-Kramers doublet mixing elements which support the ab initio findings the electric field orientated perpendicular to the high-symmetry axis generate larger spinelectric coupling.
The comparison of spin-electric couplings derived from perturbation theory with the spinelectric couplings from ab initio and the electric field model are shown in figure S6.Perturbation theory correctly predicts the magnetic field orientation dependence for E x and E z when comparing to ab initio with B y and B z having very good agreement.However perturbation theory under-predicts the spin-electric coupling for B x by approximately a)  The proposed experiment would manipulate the molecular spin using a resonant frequency electric field pulse.Where the resonant frequency (ω) is calculated from the splitting of the ground doublet due to the weak static magnetic field.

H Ez
where ω is the resonant frequency that drives spin population between states |1 ′ − ⟩ and |1 ′ + ⟩.To model these dynamics, we use evolve a spin population density matrix using the Liouville-Von Neumann equation which is defined as the commutator between the time dependent Hamiltonian and the spin density matrix.
Where H(t) is the total time dependent Hamiltonian of the system at time t, which includes the equilibrium Hamiltonian (H eq ) and the oscillating electric field Hamiltonian (H ⃗ E ) as shown in equation (S45).
We write equation (S43) such that the density matrix ρ(t) is a column vector noting that Hence the dynamics of the system can be solved using an ordinary differential equation (ODE) solver rather than using explicit numerical integration at each time step.The modified Liouville-Von Neuamn equation is given by: The transition probability between states |i⟩ and |j⟩ (|c ij (t)| 2 ) shown in figure S8c as a function of pulse duration and driving frequency detuning was calculated using equation Where Ω ij is the Rabi frequency of the transition between |i⟩ and |j⟩, and δ is the detuning of the drive frequency.
Where the Rabi frequency Ω ij of the transition between states |i⟩ and |j⟩ is calculated using equation (S49).We note our dynamics do not fully reproduce experimental conditions for a number of reasons.Firstly, in our π-pulse simulation our dynamics are unitary and neglect spin decoherence due to spin-phonon coupling.Secondly, we make a number of assumptions in our oscillating electric field Hamiltonian term.Optimisation of the molecular geometry or using our analytical model produce the equilibrium geometry in the applied electric field.However, these optimisations do not account for the time scale of molecular distortion.Hence, when the electric field Hamiltonian oscillates in our simulation we make the assumption that in the period of ω our geometry can distort from the equilibrium zerofield geometry to the maximum distorted geometry at that specific electric field strength.
Finally, our simulation makes the assumption of a perfectly uniform pulse, i.e. we do not include any shaping to our pulse.

S4 Pseudo C 2 Irrep Symmetry Decomposition
Distortions can not only be written as a linear combination of normal mode displacements but can be written as a linear combination of any arbitrary orthonormal basis.To provide further insight into the type of distortion that has the largest influence on spinelectric coupling.The distortion due to an applied electric field can be decomposed into transformed into a symmetry adapted coordinate basis.Where a set of displacement vectors are constructed such that (pseudo)-equivalent atoms related by (pseudo)-symmetry are simultaneously displaced in in-phase and out-of-phase displacements, i.e. each displacement is categorised by its definite pseudo-symmetry irreducible representation of the point group.These explicit displacements are orthonormalised using the Gram-Schmidt procedure, which takes the set of linear independent displacement vectors that span the 3N nuclear degrees of freedom and generates a set of orthonormal basis vectors of the same 3N dimensions.This forms the new symmetry adapted coordinate basis, Z.Where the basis vectors are categorised into their corresponding pseudo-symmetry irrep, that consists of symmetry breaking (B symmetry) and symmetry preserving (A symmetry) orthogonal coordinate basis vectors.Furthermore, energy invariant rigid-body rotations and translations are removed as these displacements do not contribute to geometric distortions.The distorted structure due to an electric field, r E , can be written as a liner combination of both A and B symmetry distortions.
The symmetry adapted coordinate basis Z which has the matrix form.
Using the new orthonormal Z basis the magnitude of each displacement vector, ⃗ z, can be projected. ( As the matrix Z is an orthogonal basis ZZ T = I such that we recover a vector of displacement magnitudes z which are redistributed into ⃗ a, ⃗ b using equation (S51).A set of 3N − 6 structures are then produced with a single distortion displaced.
The crystal field parameters were extracted using the Linear Vibronic Coupling model and projected onto a crystal field Hamiltonian.Rotating into the equilibrium geometries eigenbasis allows for the computation of the contribution each distortion has to the total spin-electric coupling and total contribution from symmetry preserving A distortions and symmetry breaking B distortions.Table S2: Crystal field states of the equilibrium geometry of Tm(N † † ) 2 .Where the principal g-values are given in the molecular reference to aid in the discussion of the magnetic field orientation dependence of the spin-electric coupling in the main text.

Molcas
Crystal Field Principal G values Energy (cm −1 ) Energy (cm −1 ) Gx Gy Gz Crystal field wavefunction 0.00 0.00  S3: Spin electric coupling values calculated from ab initio for each electric and magnetic file orientation, for each electric field strength, and for a magnetic field strength of 320 mT construct a new orthonormal basis of internal coordinate displacement vectors U such they have dimensions [3N × 3N − k] where the energy invariant rigid body displacements are removed.We can then transform the set of linear equations into the new reduced rank U basis to form a new set of linear equations with the energy invariant rigid body rotations and translations projected out.

E
using the derivatives obtained from the LVC model.The electric field model gives structural distortions as a function of electric field strength.We use the autodiff library which is a part of the JAX library to obtain the analytical derivative of the structural distortion as a function of electric field strength ∂⃗ r ∂ ⃗ E .Using the derivatives of the crystal field Hamiltonian with respect to atomic coordi-

Figure S1 :
Figure S1: The atomic coordinate Root Mean Square Distance (RMSD) of the ab initio and electric field model geometries in comparison to the equilibrium geometry as a function of electric field strength.

Figure S2 :Figure S3 :
Figure S2: The decomposition of the ab initio and analytical model predicted structures into the symmetry adapted coordinated basis with the contribution each symmetry has to the total spin-electric coupling magnitude.For a magnetic field of B y = 320 mT and an electric field of E y = 3 × 10 9 V m −1 .

Figure S4 :Figure S5 :
Figure S4: The real and imaginary derivatives of the distorted crystal field Hamiltonian in the SO basis for each electric field orientation (x, y, z).8 Figure S7: Heat-map representations of the electric field (a) and Zeeman (b) Hamiltonian's in the eigenstate basis of the equilibrium geometry crystal field Hamiltonian.Such that the matrix elements ⟨m|H Zee |1 ∓ ⟩ and ⟨1|H ⃗ E |m⟩ from perturbation theory are the off diagonal matrix elements of the electric field and Zeeman Hamiltonian's shown.
Figure S8: a) The pulse duration required to perform a π-pulse as a function of field strength for E z and all magnetic field orientations using the spin-electric couplings calculated ab initio.b) The Rabi oscillation of state |1 ′ − ⟩.The simulations used an applied electric field magnitude of 3.08 × 10 9 V m −1 and a magnetic field magnitude of 320 mT both orientated along the molecular z axis and a driving frequency of 1 × 10 11 Hz.c) The transition probability of the transition from state |1 ′ − ⟩ to state |1 ′ + ⟩ as a function of time and detuning of the resonant drive frequency.

Figure S9 :Figure S10 :
FigureS9: The contribution to the total spin-electric coupling magnitude due to an applied electric field along each Cartesian orientation (3×10 9 V m −1 ) between the ground Kramers doublet split by the Zeeman effect due to a magnetic field aligned along Cartesian x (320 mT) for each 'A' symmetry mode and each 'B' symmetry mode in the symmetry adapted coordinate basis for the pseudo-C 2 point group.
states are well represented by as projections of the total angular momentum, ⟨J x ⟩, and show a bistable m j = ±1/2 ground state.