Quantum operators for the computation of exponential weighted integrals of expectation values

TL;DR

This paper introduces an analytical 'integral operator' method for efficiently computing exponential weighted expectation values of quantum operators, enhancing speed and accuracy in NMR and Muon Spectroscopy simulations.

Contribution

The paper presents a novel analytical approach that improves the efficiency and precision of calculating exponential weighted expectation values in quantum systems.

Findings

Faster computation of expectation values compared to numerical methods

Higher precision in simulation results for NMR and Muon Spectroscopy

Applicable to a range of quantum expectation value calculations

Abstract

An analytically derived 'integral operator' approach is introduced to estimate the expectation value of a quantum operator for an evolving state weighted with an exponential function. This allows to compute quantities useful in Nuclear Magnetic Resonance or Muon Spectroscopy experiment simulations with greater speed and precision than the standard numerical integration approach.

Equations

$$I_{ALC}(B)\propto {\int }_0^{\infty }⟨S_z⟩(B,t)e^{-at}dt$$

$${\int }_0^t⟨\psi |e^{\frac{iH{t}^{\prime }}{\hslash }}{S}_z{e}^{-a{t}^{\prime }}{e}^{-\frac{iH{t}^{\prime }}{\hslash }}|\psi ⟩d{t}^{\prime }=⟨\psi |e^{\frac{iHt}{\hslash }}O(t)e^{\frac{-iHt}{\hslash }}|\psi ⟩-⟨\psi |O(0)|\psi ⟩$$

$$S_z{e}^{-at}=\frac{i}{\hslash }[H,O]+\frac{\partial O}{\partial t}$$

$$s_{ij}{e}^{-at}=\frac{i}{\hslash }{o}_{ij}({\lambda }_i-{\lambda }_j)+{\dot{o}}_{ij}$$

$$\frac{i}{\hslash }{o}_{ij}({\lambda }_i-{\lambda }_j)+\dot{o_{ij}}=0⟹o_{ij}=o_{ij}^0{e}^{-i\frac{{\lambda }_i-{\lambda }_j}{\hslash }t}$$

$$s_{ij}{e}^{-at}=\frac{i}{\hslash }{o}_{ij}({\lambda }_i-{\lambda }_j)+\dot{o_{ij}}⟹o_{ij}=\frac{s_{ij}}{\frac{i}{\hslash }({\lambda }_i-{\lambda }_j)-a}{e}^{-at}$$

$${\int }_0^t⟨\psi (t)|S_z{e}^{-a{t}^{\prime }}|\psi (t)⟩d{t}^{\prime }=⟨\psi (0)|P(t)|\psi (0)⟩$$

$$p_{ij}(t)=\frac{s_{ij}}{\frac{i}{\hslash }({\lambda }_i-{\lambda }_j)-a}\left(e^{-at+i\frac{{\lambda }_i-{\lambda }_j}{\hslash }t}-1\right)$$

$${\overline{S}}_N=\frac{{\int }_0^t⟨S_N⟩(t^{\prime })e^{-\frac{t^{\prime }}{\tau }}d{t}^{\prime }}{{\int }_0^t{e}^{-\frac{t^{\prime }}{\tau }}d{t}^{\prime }}$$

Full text

1 Introduction

A common technique used in muon spectroscopy is avoided level crossing resonance [1, 2]. In this kind of experiment, a beam of polarised muons is injected into a sample under a strong external magnetic field, and their subsequent decay causes emission of positrons whose momentum is directly related to the expectation value of $S_{z}$. The final spectrum is determined by the integral of this quantity over time, weighted with the decaying density of muons, for a given external magnetic field $B$:

[TABLE]

with $a=\frac{1}{\tau_{\mu}}$ here is the natural decay time of the muons. To compute the quantity seen in equation 1 in a simulation the most straightforward method involves evolving the state of the system with the time-dependent Schrödinger equation, tracing out the eigenvalue and numerically integrating. In this note a different approach to compute this kind of integral is proposed. This involves finding an ‘integral operator’ such that its expectation value will match exactly the result of equation 1. It should also be noted that many similar problems can be described in the same general form; for example, if $a=i\omega$, we get the Fourier transform at a given frequency $\omega$. This method could thus apply to exploring specific frequencies of the spectrum of an expectation value without having to compute its full time evolution.

2 Derivation

Let us consider the case of a quantum system evolving under a general Hamiltonian. In the Schrödinger picture, the wavefunction of such a system would be time dependent and the operator whose expectation value we wish to average (in this case, $S_{z}$) would be time independent. However we will be looking for a Heisenberg picture operator $O$ such that

[TABLE]

where we took care to make the time dependence of the wave function explicit. By recalling the expression for the derivative of an expectation value in the Heisenberg picture we can write

[TABLE]

It should be noted that since we moved from expectation values to operators this is not a necessary condition, but it is sufficient to guarantee that equation 2 holds.

For systems with a finite number of states we can describe the operators as matrices without loss of generality. For convenience, let us consider a Hilbert space in which the Hamiltonian is diagonal and has eigenvalues $\lambda_{i}$. We will also rewrite $S_{z}$ as a general operator $S$, since we don’t know anything about its form in this basis, and element-wise we can impose:

[TABLE]

where lowercase letters refer to elements of the correspondingly named uppercase operators. This differential equation can be solved easily. First we consider a solution for the associated homogeneous equation:

[TABLE]

where we notice the presence of a matrix of integration constants $O^{0}$. Then we find a particular solution:

[TABLE]

The integration constants may give us some thought, but in fact do not matter. Since this is a Heisenberg operator, in order to account for the time evolution it will need to be multiplied by a factor of $e^{i\frac{\lambda_{i}-\lambda_{j}}{\hbar}t}$. But for equation 5 that means the entire term becomes a constant, and since equation 2 shows that the final value of interest only depends on the difference between the operator at time $t$ and [math], such constants always cancel each other. Therefore equation 6 gives in fact the full solution we need. Including the time dependence factor we can define a new operator $P$ and write

[TABLE]

where

[TABLE]

3 Performance

We now proceed to test the performance of the proposed algorithm by generating random Hermitian Hamiltonians of variable size. The simulations were performed using Python and the numerical libraries Numpy and Scipy [3, 4] and run on a desktop PC. The script used to produce the data is given for reference as supplementary material and to provide an implementation of the algorithm. The generated Hamiltonians are all created as uniform random numbers between -0.5 and 0.5 for both real and imaginary part. Units are arbitrary. The first test focuses on evaluating the time evolution of the number operator $S_{N}$ for a density matrix initialized in the highest state and computing the weighted time average

[TABLE]

for $\tau=10.0$. Figure 1a shows the example of this integral computed at different times with the integral operator as well as with a numerical integration method for different amount of steps. The numerical integration was carried out with Numpy’s default trapezoidal integration function; the number of steps is fixed for all points (meaning that the resolution is better at short times). The Hamiltonian diagonalization is not counted in the time as it is a necessary step for both methods. It is possible to compute the time evolution numerically without diagonalizing the Hamiltonian, but that introduces a further source of errors and requires smaller time steps, therefore is not desirable for small systems and has not been tried. It can be seen how increasing the number of steps makes the numerical solution progressively closer to the analytical result. In figure 1b we can however take a look at the average performance of these approaches, computed over 200 averages for different system sizes and step numbers. It becomes obvious that even for the lowest integration step counts the numerical approach is far inferior to the integral operator, taking almost an order of magnitude longer for worse results.

Figure 2 shows instead the example of an application to Fourier transform of the time evolution generated by a random Hamiltonian of size 4. While it’s not strictly necessary, in this example a damping time $\tau=50$ has also been used to broaden the peaks and represent the most general case. After performing the FFT of the numerically computed time evolution a peak is identified and a small window around it is explored with the integral operator. The results are in excellent agreement. Use for Fourier transform is not time efficient: in this example, the FFT on 1000 points has taken approximately $0.01\,s$, while each point computed with the integral operator has taken almost $45\,\mu s$. At this rate, if we had to compute the whole spectrum that way it would take almost 5 times longer. However, this method allows one to explore single values or small frequency windows offset from the origin, and in the case that is all that is needed it could result convenient.

4 Conclusions

A simple method that to our knowledge was not previously documented to compute exponentially weighted time integrals of an expectation value in quantum mechanics simulation has this proposed. This approach is shown to be significantly faster and more precise than traditional ones. Numerical simulations in the field of muon spectroscopy and nuclear magnetic resonance could potentially benefit from using this method in select cases of interest.

5 Acknowledgments

Thanks to Stewart Clark, Dominik Jochym, Barbara Montanari, Leonardo Bernasconi and Leandro Liborio for useful conversations. This work was supported by the Collaborative Computational Project for NMR Crystallography, funded by the EPSRC (UK) Grants EP/J010510/1 and EP/M022501/1.