Quantum computation of the lowest-energy Kramers states and magnetic g -factors of rare earth ions in crystals

We present the results of the quantum calculation of the ground state energies and magnetic g -factors of two rare earth (RE) ions: Yb 3+ in Y 2 Ti 2 O 7 crystal and Er 3+ in YPO 4 crystal. The Variational Quantum Eigensolver (VQE) algorithm has been performed on 5-qubit IBM super-conducting quantum computers via IBM Quantum Experience cloud access. The Hamiltonian of the lowest spectroscopic multiplet of each RE ion, containing crystal field and Zeeman interaction, has been projected onto the collective states of three (Yb 3+ ) and four (Er 3+ ) coupled transmon qubits. The lowest-energy states of RE ions have been found by minimizing the mean energy in ∼ 250 – 350 iterations of the algorithm: the first part was performed on a quantum simulator, and the last 25 iterations were conducted on the real quantum computing hardware. All the calculated ground-state energies and magnetic g -factors agree well with their exact values, while the estimated error of 2 ÷ 15% is mostly attributed to the decoherence associated with the two-qubit operations.


Introduction
During the past decade, the development of working physical realizations of multiqubit quantum computers has been quite rapid.In 2011, the D-Wave One system was announced: a 128-qubit quantum annealing computer, which, however, was restricted to certain tasks and lacked the possibility to implement arbitrary quantum gates necessary to realize the "iconic" quantum algorithms [1].In 2015, D-Wave Systems announced the general availability of the D-Wave 2X, a 1000+ qubit quantum annealing computer.In 2016, IBM introduced their first superconducting quantum computer accessible by the vast scientific community through their web cloud [2].In 2017, another company, Rigetti Computing, announced the public beta availability of a quantum cloud computing platform [3].In 2019, Google claimed that their quantum 54-qubit processor achieved quantum supremacy [4].
However, quantum technology is still far from fault-tolerant computations, and most of algorithms we know are extremely sensitive to noise.Quantum computers with > 50 qubits may be able to perform tasks which surpass the capabilities of today's classical digital computers, but noise in quantum gates will limit the size of quantum circuits that can be executed reliably.In this sense, modern quantum devices can be related to Noisy Intermediate-Scale Quantum (NISQ) technology [5].
Modern quantum computing platforms enable the simulation of various coupled electronic quantum systems, including Ising and central spin models [6,7].However, possible implementations are obviously not restricted to exchange-coupled spin aggregates.Another interesting problem tackled by the modern quantum computing devices is simulation of the lowest-energy state of simple atomic clusters like the hydrogen molecule [8].A possible generalization would be a study of a multi-electronic ion (a transition metal or rare earth ion) in a crystal.Currently, we are unaware of any research on quantum computer simulations of spectroscopic properties of rare earth ions in crystals conducted so far.Each rare earth impurity ion has approximately 70 electrons interacting with each other, their nucleus and the nearby ions of the host crystal, which makes almost impossible to find an exact solution using any classical algorithm.
In this work, we deal with a more simple problem of finding the lowest-energy states of two Kramers ions (Yb 3+ and Er 3+ ) taking into account only their lower-energy multiplet states.The electrostatic interaction between the electrons of the valence the 4f shell, together with their spin-orbit interaction, produce a number of well-separated spectroscopic multiplets.The total angular momentum of 4f shell and its projection are considered as good quantum numbers.The electrostatic interaction produces some splitting between these levels within the multiplet, which can be modelled by introducing the so-called crystal field Hamiltonian [9].The lowestenergy multiplets of Yb 3+ ( 2 F 7/2 ) and Er 3+ ( 4 I 15/2 ) contain 8 and 16 states respectively, and can be simulated using 3 and 4 coupled qubits.

Crystal Field and Spectroscopic Data
Over the last two decades, the pyrochlore crystals with various impurity RE ions, being common physical realizations of a geometrically frustrated system, have attracted the attention of many researchers [10].Here, we consider Y 2 Ti 2 O 7 crystal doped with Yb 3+ ions containing 13 electrons in the valence 4f shell.Their electrostatic interaction with the surrounding Y 3+ , Ti 4+ and O 2− ions of the host crystal projected onto the lowest-energy multiplet 2 F 7/2 (orbital momentum L = 3, spin S = 1/2, total angular momentum J = 7/2) can be described by the effective crystal field operator, which, in the case of trigonal symmetry of the Yb 3+ site, is defined as: Above, O k p are Stevens operators (linear combinations of spherical tensor operators [9,11]), the parameters α p define the reduced matrix elements [9].The crystal field parameters B k p determined previously [11] are presented in Table 1.
The second system that we consider in this work is the yttrium orthophosphate crystal YPO 4 activated with Er 3+ ions.This system is a promising telecom-wavelength material for applications in quantum electronics and quantum information processing [12].Er 3+ ions substitute for Y 3+ ions on the sites of D 2d symmetry.The crystal-field interaction projected onto the subspace of 4 I 15/2 multiplet (L = 6, S = 3/2, momentum J = 15/2) is expressed as: The parameters B k p determined previously [12] are presented in Table 1.The interaction with the static magnetic field, when projected onto the ground multiplet 2S+1 L J , results in the following Zeeman interaction term: where µ B is the Bohr magneton, B is the magnetic field vector, and g J is the Landé g-factor (below, g e = 2.0023 is electronic g-factor) [9] g J = (g e − 1) • J(J + 1) + S(S + 1) − L(L + 1) 2J(J + 1) The matrices of the operators (1)-(3) in the basis of |J, M J ⟩ states have been computed on a classical computer, and then expanded into series of either three-or four-fold tensor products of Pauli matrices, see Section 4.

Computational Hardware and Software
We conduct our experiments using cloud access to IBM devices via the Qiskit framework [13].
IBM processors utilize transmon qubit technology [14,15].The transmon qubit is a form of superconducting charge qubit with additional capacity resistant to charge noise.The device comprises of fixed-frequency Josephson-junction-based transmon qubits, with individual superconducting coplanar waveguide (CPW) resonators for qubit control and readout, and another pair of CPW resonators providing the qubit connectivity.This fixed-frequency architecture is favorable for obtaining long coherence times, and the qubit control and readout is achieved using only microwave pulses.Each quantum chip is calibrated on a daily basis, thus minimizing its single-qubit and CNOT error rates.Since our problems do not involve more than 4 qubits, we mainly use a 5-qubit devices depicted in Figure 1.A quantum algorithm for the ground-state problem, as proposed in works [16,17], relies on preparation of the initial state that has a large overlap with the ground state and solves the problem using the quantum phase estimation algorithm (PEA).This approach requires sufficient circuit depths and long coherence times, and such requirements cannot be fully handled by contemporary hardware.In our work, we use the VQE algorithm introduced in 2014 by Peruzzo and others [18,19], which has less strict hardware requirements and is suitable for eigenvalue calculations using NISQ devices.The low-depth circuits utilized by VQE can help to avoid decoherence errors during computation.The algorithm is based on Ritz variational principle.If we parametrize an eigenfunction |ψ(θ)⟩ of the Hamiltonian H by θ, and if we minimize left side of inequality ( 5), we will be able to approach the ground state energy E 0 of H: The VQE is a semi-classical algorithm omprising of two parts: (i) the quantum-computer part that evaluates the expectation value (5) of the target Hamiltonian, and (ii) the classical-computer part that seeks the minimum of the expectation value with respect to θ using a classical optimization algorithm.
In order to parametrize the eigenstate, we have to adopt one of the strategies used to build so-called ansatz circuits to prepare trial state.The first one relies on the knowledge of properties and symmetries of system's Hamiltonian.Such ansatzes are widely used in quantum chemistry problems, e.g. the Unitary Coupled Cluster Ansatz (UCCA) in recent studies [20][21][22].Second strategy is called Hardware-Efficient Ansatz: this one utilizes quantum gates that are naturally available on the quantum hardware [23].We have chosen the second strategy for our work, because it requires less circuit depth, as we should exclude as many error sources as possible due to the limited coherence times of contemporary IBM quantum devices.
Our ansatz acts on this state as unitary transformation U (θ): where θ = {θ 1 , θ 2 , . . ., θ m } is a set of m parameters.A generalized unitary transformation performed on a single qubit modeling the states of a spin-1/2 particle could be written as a combination of three spin rotations around z and y axes in the form [24]: where q identifies the qubit and d identifies the algorithm depth layer of the circuit.For N qubits, we can write: where U ENT represents an entangling layer of two-qubit gates.Since the qubits are all initialized in their ground state |0⟩, the first R Z -set of rotations U q,0 has no effect on it and is not implemented, this gives the algorithm N (3d + 2) parameters to be optimized, where d is the depth of the algorithm.After each trial state is prepared, we estimate the associated mean energy (5) by measuring the expectation values of the individual Pauli tensor products where each A, B, . . ., C is either one of the three Pauli operators X q , Y q , Z q of the qubit q, or the identity operator I q .The Hamiltonian H is presented as a sum of such products, see section 4. Since we work in the basis of the eigenstates of Z q operators, Z q values can be obtained by directly measuring the final state of a qubit q, and then averaged over a large number of runs (8192 in our case) to produce ⟨. . .Z q . ..⟩.In order to make measurements for another two Pauli operators X q , Y q , we apply additional π/2 rotations prior to measurement of Z q .
After the preparation of the trial state is finished and all the measurements are done, the algorithm continues with its classical part.This part involves the summation of all average values ⟨P ⟩ comprising the expectation value of H, and a search for a better set of parameters that would lower ⟨H⟩.After that, the whole cycle repeats again until the convergence conditions are fulfilled.

Quantum computation of the Lowest Energy State of a RE ion
To conform with VQE approach, we map the problem's Hamiltonian H to an N -qubit Hamiltonian containing Pauli tensor products.We use the Hilbert-Schmidt inner product decomposition: Above, i, j, . . ., k refer to {X, Y, Z, I} operators.In the case of Yb 3+ ion's 2 F 7/2 state, we map its Hamiltonian (1) to 8 states of three qubits, and obtain (omitting the upper indices for simplicity) Next, we utilize the fact that the crystal field Hamiltonian contains only real terms.In this case, one can simplify the one-qubit rotation (7) to the form that contains only one parameter [25].With this simplification, our algorithm has only N (d + 1) parameters to be optimized.We choose d = 1 for both crystals.
For the classical optimization part, we choose Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm, since it is robust to stochastic fluctuations and requires only two cost-function evaluations [23], irrespective of the dimensionality of the parameter space.Because we did not have continuous access to IBM quantum devices, to avoid time consumption, the first part of the optimization procedure was performed on the simulator of a quantum computer.The energy was estimated on the actual hardware during the last 25 optimization iterations, and then averaged over several iterations to improve results.
It is possible to reduce the number of terms in the Hamiltonian using simple qubit-wise commutativity (QWC) [26] for Pauli operators.Two Pauli strings QWC-commute if, at each index, the corresponding two Pauli operators commute.For instance, {XX, XI, IX, II} is a QWC partition, so all these Pauli strings can be measured simultaneously, and the results then re-calculated straightforwardly to obtain the corresponding expectation values.Since each additional measurement is likely to increase the error, the use of QWC is advantageous.To this extent, we can rewrite the Hamiltonian (10) in the form: For Er 3+ ion, the Hamiltonian (2) maps to the Hilbert space of four qubits: Again, using QWC approach, we reduce the number of independently measured terms to 4. The Zeeman interaction (3) has been decomposed into the Pauli strings in the same manner.
The magnetic g-factors of the lowest Kramers doublet of each multiplet have been calculated by minimizing the energy of the lowest state of the ion subjected to the external magnetic field B directed along (g ∥ ) or perpendicular to (g ⊥ ) the crystal's symmetry axis.In either case, where θ * represents the parameter set minimizing the expectation value of total Hamiltonian For the case of Er 3+ , we used the different way to calculate g ⊥ because the expansion of Hamiltonian H (Er) CF + H ⊥ Z consisted of more than 25 terms, which could not be reduced by QWC method.This led to the accumulation of significant statistical error and making it impossible to obtain g ⊥ with an acceptable accuracy.
We took advantage of the fact that Er 3+ is a Kramers ion with odd number of electrons on valence shell, and its multiplet states are doubly degenerate and form so-called Kramers doublet {|ψ 1 ⟩, |ψ 2 ⟩.We replace these states by the states of an effective spin 1/2 {|ψ 1 ⟩ ≡ | ↓⟩, |ψ 2 ⟩ ≡ | ↑⟩} and write an effective Hamiltonian: Up to a global phase factor we write: The relation between these degenerate Kramers states can be expressed with time-reversal operator T : The state |ψ 1 ⟩ can be obtained after the minimization of , where H ∥ Z does not require the introduction of additional Pauli strings.The expansion of J x T in (16) into the Pauli tensor products, in combination with the QWC, consists only of 8 terms.We have also calculated the standard error where σ is a standard deviation of the result and n is the number of measurement repetitions.The standard deviation in the measurement of one of the Pauli strings in the Hamiltonian is given by The standard error in the energy measurement is then upper bounded by [23] α This quantity is shown in the inset of Figure 3 and Figure 2.

Results and Discussion
We performed a quantum computations of ground state energies and g-factors of two rare-earth ion impurities via cloud access to IBM quantum devices.We used Hilbert-Schmidt decomposition along with QWC-method to obtain a weighted sum of Pauli operator tensor products.Using a hardware-efficient ansatz to prepare trial states we found ground state energies of impurity ions.By adding an interaction with an external magnetic field to our initial Hamiltonians we determined the projections of ground state g-factors.
Some of the results obtained for the YPO 4 :Er 3+ and Y 2 Ti 2 O 7 :Yb 3+ crystals are shown in Figures 2 and 3 respectively.The SPSA minimization was performed with qiskit qasm simulator [13], and the ground-state energy ⟨ψ(θ * )|H|ψ(θ * )⟩ has been calculated on the real hardware during the last 25 steps of optimization with 5-qubit quantum chip, and then averaged (dashed green line).The variance for 25 estimated energy values is Σ 2 = 1.32 cm −2 for Er 3+ and Σ 2 = 184.98cm −2 for Yb 3+ .The estimated error is attributed mainly to the decoherence associated with the two-qubit operations.The exact eigenvalues were calculated with standard diagonalization procedure (dashed black line).As shown in Figure 2, the difference between the final energy estimate and the "exact" ground state energy is ∼ 10%.g-factors were calculated according to Eq. 13.One should bear in mind that, due to nonlin-  earity of the average state energy with respect to the field value B, the calculated value (13) would depend on B. We chose B=2000 G and 10000 G for Er 3+ and Yb 3+ , respectively, to assure that these nonlinear contributions played insignificant role, and, at the same time, to allow the reliable calculation of the ground doublet splitting.The calculated g-factor error values, if compared with the exact classical diagonalization, are 2 ÷ 15 % (see Table 2).
These results show that IBM's near-term quantum computers are capable to simulate simple physical systems like an impurity ion in a crystal.However it's too early to say that these devices can operate efficiently.The amount of errors is still sufficient, the biggest problem here is short coherence time of qubits that limits the length of a quantum circuit.In turn, this limits the number of variational parameters in the ansatz circuit, As a result, we can't consider quantum systems with large state spaces (the largest dimension of the state space in our article is 16).This fact does not allow us to scale this approach to more complex physical systems.
Scaling our task to address problems involving interactions with more than four particles will require developing a new approach for constructing ansatz state, one should consider symmetries that limit the size of Hilbert space for more precise search.Straightforward use of Hardware-efficient ansatz will lead to accumulating errors due to increasing ansatz depth, and the minimization procedure will probably stuck in local minima.However, in our small-scale problem, the SPSA optimization algorithm demonstrated good performance, achieving convergence within 100-150 iterations in all cases.We used the measurement error mitigation method to reduce the errors that occure during the process of qubits state measurement (such as errors associated with qubit signal recognition).Most errors that result in this difference between exact or experimental values and the values obtained with quantum computer should be attributed to non-ideal one-qubit and two-qubit operations.Besides, results can be influenced by the current calibration quality since two-qubit, one-qubit and readout errors vary with each daily calibration.
Furthermore, scaling will lead to large number of Pauli strings in expansion of target Hamiltonian, which suggests the need to significantly decrease measurement errors with new methods of quantum data post-processing.Neural networks are supposed to be a good way to address this problem, according to recent studies [27].For NISQ devices, such semi-classical algorithms like VQE are expected to be widely used in the coming years.We expect that quantum computers may give us a computing speed advantage since in the case of variational algorithms and n-particle Hamiltonian, all quantum resources (quantum memory, number of qubits, number of measurements) grow polynomially with the number of particles.Classical procedure of the Hamiltonian diagonalization requires an exponential number of operations with respect to the number of particles (∼ 2 n ), and so do classical resources.

Figure 1 .
Figure 1.Layout of a 5-qubit IBM quantum computer.The arrows indicate the coupling arranged between the adjacent qubits.

Figure 2 .
Figure 2. Results of the minimization procedure in the case of Er 3+ crystal field Hamiltonian.The first part (red and blue dots) correspond to the quantum simulator results obtained for the two close sets of parameters determining the next-iteration set in SPSA.The last part of the curve (green dots) corresponds to the quantum hardware part of the calculation.Vertical bars represent the upper error bound (19).

Figure 3 .
Figure 3. Results of the minimization procedure in the case of Yb 3+ crystal field Hamiltonian.The first part (red and blue dots) correspond to the quantum simulator results obtained for the two close sets of parameters determining the next-iteration set in SPSA.The last part of the curve (green dots) corresponds to the quantum hardware part of the calculation.Vertical bars represent the upper error bound (19).

Table 1 .
Crystal field parameters B k p (cm −1 ) and the reduced matrix elements for Yb 3+ in Y 2 Ti 2 O 7 crystal and Er 3+ in YPO 4 crystal.