text:
# Decoherence effects on the tunneling rate of paramagnetic and superparamagnetic particles

## Abstract

We analyze the effects of the environment on the spin tunneling process of paramagnetic and superparamagnetic particles and conclude that the “bare” macroscopic tunneling rate is hardly affected in such case, but others more effective processes come up to change the magnetization state of the particle. We measure the degree of coherence loss by the linear entropy which will account for all correlations present in the process. We conclude that for both, paramagnetic and superparamagnetic particles, the decoherence time scale is extremely short ($`10^{8..16}s`$), indicating that coherent tunneling should be strongly suppressed in favor of incoherent tunneling, i.e., the population of higher levels with subsequent decay. Interestingly enough, the ground state tunneling rate is hardly affected by such dissipative mechanisms. However other processes immediately give rise to new (incoherent) possibilities of crossing the barrier.

, , ,

Recently discovered magnetic molecules are of interest because of the gigantic relaxation times of their magnetization. Recent experiments on magnetic relaxation of molecular crystals of Mn12 and Fe8 find strong evidence for tunneling mediated relaxation at low temperatures. Several theoretical studies have been performed, leading to tunneling rates of the order of a few months. The situation with superparamagnetic particles is less clear, however. A change in magnetization direction at low temperatures is indeed observed, the question remaining as to whether this is due to a unitary tunneling effect. Contradictory statements regarding the observability of the tunneling effect can be found in the literature.

In the present letter, we set up a collective model for the tunneling of magnetic molecules, derive the corresponding tunneling rates and include dissipation via coupling with a Caldeira-Leggett type of environment. As expected, the decoherence time scale is found to be much shorter than the characteristic tunneling time. However, surprisingly enough, we find that the tunneling rate is hardly affected by the environment induced dissipation. We therefore conclude, within the scope of the present model, that the spin-phonon interaction does not significantly interfere with the coherences which are essential for the tunneling process. It however very quickly opens up several alternative routes for the system to change its magnetization, rendering the unitary tunneling mechanism but a possibility among many others. Our result strongly supports the picture of incoherent tunneling, since coherence is destroyed just about immediately.

We start by considering a Hamiltonian which describes the magnetization properties of such molecules,

$$H_S=AS_z^2+B(S_+^2+S_{}^2).$$
(1)
The anisotropy constants satisfy $`AB>0`$. $`S_z`$ and $`S_\pm =S_x\pm iS_y`$ are spin operators. The energy splitting of the ground state due to the second term on the r.h.s. of Eq. (1) can be calculated by means of path integral techniques and the instanton method. Details of this calculation can be found in Ref. . We get

$$\mathrm{\Delta }E_{\text{inst}}=\frac{8As^{3/2}}{\sqrt{\pi }}\left(\frac{B}{A}\right)^s\mathrm{cos}(\pi s).$$
(2)
This result can be compared with that obtained by Hartmann-Boutron by means of others methods. For magnetic molecules, like Mn12, $`s=10`$ and for typical superparamagnetic particles, $`s3000\mathrm{}`$. For both we have $`BA`$, and therefore $`\mathrm{\Delta }E_{\text{inst}}0`$, suggesting then that th observation of this phenomenon in the laboratory is impossible.

Obviously thermally activated over-barrier turnover of the magnetization direction at sufficiently high temperatures (normal superparamagnetism) is possible. In the following, however, we analyze the possibility of spin-tunneling, i.e., under-barrier changes of the magnetization direction.

The model of a single collective variable is apparently too naive. We therefore next take into account the fact that the magnetic particles are not isolated and estimate the characteristic time of environmental effects. We consider phonon degrees of freedom, represented as a set of harmonic oscillators described by the Hamiltonian

$$H_R=\underset{\alpha }{}\mathrm{}\omega _\alpha a_\alpha ^{}a_\alpha $$
(3)
and coupled to $`H_S`$ (this model has been proposed in )

$$H_{\text{int}}=B_{\text{int}}(S_+^2+S_{}^2)\underset{\alpha }{}\sqrt{\frac{\mathrm{}}{2M\omega _\alpha }}k_\alpha (a_\alpha ^{}+a_\alpha ),$$
(4)
where $`B_{\text{int}}`$ is an overall interaction constant, $`M`$ the particle mass, $`k_\alpha `$ the wave number associated to the phonon with frequency $`\omega _\alpha `$ and $`a_\alpha ^{}(a_\alpha )`$ creation (annihilation) operators of the environmental degrees of freedom.

Assuming this interaction to be weak enough with respect to the collective dynamics, Eq. (1), we use perturbation theory to evaluate the reduced density of the spin degree of freedom,

$`\stackrel{~}{\rho }(t)`$ $`=`$ $`\stackrel{~}{\rho }(0){\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^t}𝑑t^{}[\stackrel{~}{H}_{\text{int}}(t^{}),\stackrel{~}{\rho }(0)]`$ (5)
$`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t^{}}𝑑t^{\prime \prime }[\stackrel{~}{H}_{\text{int}}(t^{}),[\stackrel{~}{H}_{\text{int}}(t^{\prime \prime }),\stackrel{~}{\rho }(0)]]+\mathrm{}`$
$`=`$ $`\stackrel{~}{\rho }^{(0)}+\stackrel{~}{\rho }^{(1)}+\stackrel{~}{\rho }^{(2)}+\mathrm{}`$
where $`\stackrel{~}{H}_{\text{int}}=\mathrm{}^{iH_0t}H_{\text{int}}\mathrm{}^{iH_0t}`$ and $`H_0H_S+H_R`$. The reduced density is then

$$\stackrel{~}{\rho }_S=\mathrm{Tr}_R[\stackrel{~}{\rho }].$$
(6)
$`\mathrm{Tr}_R`$ represents the trace over the reservoir (phonons) variables. From the above expressions it is not difficult to obtain the two quantities which lead us to our final conclusions:

1. The tunneling probability, i.e., the probability as a function of time to find our system in state $`|s`$, having been initially prepared in $`|+s`$

$$\begin{array}{cc}\hfill P_s(t)& =\mathrm{Tr}[|ss|\stackrel{~}{\rho }_S(t)]\hfill \\ & =\mathrm{sin}^2\left(\frac{\mathrm{\Delta }\omega _{\text{inst}}t}{2}\right)\left[1\frac{18B_{\text{int}}^2s_2}{\mathrm{}Mc^2\omega _D^3}\left(\frac{(\omega _{s,s2})^3}{\mathrm{}^{\mathrm{}\omega _{s,s2}/kT}1}\right)t\right]\hfill \end{array}$$
(7)
2. The coherence loss as measured by the linear entropy (or idempotency defect )

$$\begin{array}{cc}\hfill \delta _S(t)& =1\mathrm{Tr}_S[(\stackrel{~}{\rho }_S)^2]\hfill \\ & =\frac{3B_{\text{int}}^2}{4\pi ^2\mathrm{}\rho c^5}\left[\frac{s_{+2}^m(\omega _{m,m+2})^3}{1\mathrm{}^{\mathrm{}\omega _{m,m+2}/k_BT}}+\frac{s_2^m(\omega _{m,m2})^3}{1\mathrm{}^{\mathrm{}\omega _{m,m2}/k_BT}}\right]t\hfill \end{array}$$
(8)
From this formula we also define the “decoherence time” ($`\tau _{\text{dec}}`$) as being the inverse of the factor multiplying $`t`$, or in others words, the time at which $`\delta _S(t)=1`$.

In the above expressions, we have used a spectral function for the phonons given by

$$J(\omega )=\frac{\pi }{2}B_{\text{int}}^2\underset{\alpha }{}\frac{k_\alpha ^2}{M\omega _\alpha }\delta (\omega \omega _\alpha )$$
(9)
and a linear dispersion relation for the phonons $`\omega _\alpha =ck_\alpha `$, where $`c`$ is the sound velocity in the medium. Moreover, Debye’s model is used to transform the sum into an integral, i.e.

$$\underset{\alpha }{}Q(\omega _\alpha )_0^{\omega _D}𝑑\omega g(\omega )Q(\omega )$$

with $`g(\omega )=9\omega ^2/\omega _D^3`$ and $`\omega _D^3=3V/(18\pi ^2c^3)`$, $`V`$ being the particle’s volume. Nevertheless it is only necessary to know the value of the density $`\rho `$ since $`V`$ and $`M`$ always appear in formulas in the form $`M/V`$. With these ingredients, Eq. (9) becomes

$$J(\omega )=\frac{9\pi B_{\text{int}}^2\omega ^3}{2Mc^2\omega _D^3}$$
(10)
corresponding to superohmic dissipation, according to Caldeira and Leggett’s model. In Eq. (8), $`s_{\pm 2}^m(sm)(s\pm m+1)(sm1)(s\pm m+2)`$ and $`\omega _{m,m\pm 2}=(E_mE_{m\pm 2})/\mathrm{}`$. In Eq. (7), $`s_2`$ is equal $`s_2^m`$, for $`m=s`$.

Let us begin by analyzing the tunneling rate. The environmental effects are contained in the brackets of Eq. (7) and are, for typical values of the constants (see Table 1), much smaller than 1. We can therefore conceive that the terms in brackets correspond to the first terms in the expansion of the exponential $`\mathrm{}^{\gamma (\beta )t}`$, with

$$\gamma (\beta )=\frac{18B_{\text{int}}^2s_2}{\mathrm{}Mc^2\omega _D^3}\frac{(\omega _{s,s2})^3}{\mathrm{}^{\mathrm{}\omega _{s,s2}/k_BT}1}.$$

This result is thus very analogous to the tunneling rate obtained in the context of the spin-boson model, i.e., the “free” tunneling frequency is damped. Our model hypothesis, which we believe valid for magnetic particles, indicate that the tunneling rate remains essentially unaffected by the presence of the other degrees of freedom. This means, in particular, that the unitary quantum mechanical tunneling process will be hardly affected by this interaction. However, a careful look at the decoherence time (see Fig. 1) shows that very quickly ($`\tau _{\text{dec}}10^{16}`$s ($`10^8`$s) for superparamagnetic (Mn12) particles), many new channels will be open, the density matrix will be “contaminated” by other available states and lose purity. Physically this means that, although the particular coherences necessary for the macroscopic quantum coherent tunneling of the ground state remains essentially unaffected by the dissipation, states, candidate to tunnel the barrier, will not do it in a coherent way since others mechanisms for magnetization change will come into play very quickly, thus giving support to the idea that the experimentally observed magnetization change must be phonon assisted.

This is work was supported by Brazilian research agencies CNPq, FAPESP and FAPEMIG.


format: text

text:
# Quantum Simulations of Physics Problems

## 1 introduction

Quantum simulation of physical systems on a QC has acquired importance during the last years since it is believed that QCs can simulate quantum physics problems more efficiently than their classical analogues : The number of operations needed for deterministically solving a quantum many-body problem on a classical computer (CC) increases exponentially with the number of degrees of freedom of the system.

In quantum mechanics, each physical system has associated a language of operators and an algebra realizing this language, and can be considered as a possible model of quantum computation . As we discussed in a previous paper , the existence of one-to-one mappings between different languages (e.g., the Jordan-Wigner transformation that maps fermionic operators onto spin-1/2 operators) and between quantum states of different Hilbert spaces, allows the quantum simulation of one physical system by any other one. For example, a liquid nuclear magnetic resonance QC (NMR) can simulate a system of 4He atoms (hard-core bosons) because an isomorphic mapping between both algebras of observables exists.

The existence of mappings between operators allows us to construct quantum network models from sets of elementary gates, to which we map the operators of our physical system. An important remark is that these mappings can be performed efficiently: we need a number of steps that scales polynomially with the system size. However, this fact alone is not sufficient to establish that any quantum problem can be solved efficiently. One needs to show that all steps involved in the simulation (i.e., preparation of the initial state, evolution, measurement, and measurement control) can be performed with polynomial complexity. For example, the number of different eigenvalues in the two-dimensional Hubbard model scales exponentially with the system size, so QC algorithms for obtaining its energy spectrum will also require a number of operations that scales exponentially with the system size .

Typically, the degrees of freedom of the physical system over which we have quantum control constitute the model of computation. In this paper, we consider the simulation of any physical system by the standard model of quantum computation (spin-1/2 system), since this might be the language needed for the practical implementation of the quantum algorithms (e.g., NMR). Therefore, the complexity of the quantum algorithms is analyzed from the point of view of the number of resources (elementary gates) needed for their implementation in the language of the standard model. Had another model of computation being used, one should follow the same qualitative steps although the mappings and network structure would be different.

The main purpose of this work is to show how to simulate any physical process and system using the least possible number of resources. We organized the paper in the following way: In section 2 we describe the standard model of quantum computation (spin-1/2 system). Section 3 shows the mappings between physical systems governed by a generalized Pauli’s exclusion principle (fermions, etc.) and the standard model, giving examples of algorithms for the first two steps (preparation of the initial state and evolution) of the quantum simulation. In section 4 we develop similar steps for the simulation of quantum systems whose language has an infinite-dimensional representation, thus, there is no exclusion principle (e.g., canonical bosons). In section 5 we explain the measurement process used to extract information of some relevant and generic physical properties, such as correlation functions and energy spectra. We conclude with a discussion about efficiency and quantum errors (section 6), and a summary about the general statements (section 7).

## 2 standard model of quantum computation

In the standard model of quantum computation, the fundamental unit is the qubit, represented by a two level quantum system $`|𝖺=a|0+b|1`$. For a spin-1/2 particle, for example, the two “levels” are the two different orientations of the spin, $`|=|0`$ and $`|=|1`$. In this model, the algebra assigned to a system of $`N`$-qubits is built upon the Pauli spin-1/2 operators $`\sigma _x^j`$, $`\sigma _y^j`$ and $`\sigma _z^j`$ acting on the $`j`$-th qubit (individual qubit). The commutation relations for these operators satisfy an $`\underset{i=1}{\overset{N}{}}su(2)_i`$ algebra defined by ($`\mu ,\nu ,\lambda =x,y,z`$)

$$[\sigma _\mu ^j,\sigma _\nu ^k]=2i\delta _{jk}ϵ_{\mu \nu \lambda }\sigma _\lambda ^j,$$
(1)
where $`ϵ_{\mu \nu \lambda }`$ is the totally anti-symmetric Levi-Civita symbol. Sometimes it is useful to write the commutation relations in terms of the raising and lowering spin-1/2 operators

$$\sigma _\pm ^j=\frac{\sigma _x^j\pm i\sigma _y^j}{2}.$$
(2)
Any operation on a QC is represented by a unitary operator $`U`$ that evolves some initial state (boot-up state) in a way that satisfies the time-dependent Schrödinger equation for some Hamiltonian $`H`$. Any unitary operation (evolution) $`U`$ applied to a system of $`N`$ qubits can be decomposed into either single qubit rotations $`R_\mu (\vartheta )=e^{i\frac{\vartheta }{2}\sigma _\mu }`$ by an angle $`\vartheta `$ about the $`\mu `$ axis or two qubits Ising interactions $`R_{z^j,z^k}=e^{i\omega \sigma _z^i\sigma _z^j}`$. This is an important result of quantum information, since with these operations one can perform universal quantum computation. It is important to mention that we could also perform universal quantum computation with single qubit rotations and C-NOT gates  or even with different control Hamiltonians. The crucial point is that we need to have quantum control over those elementary operations in the real physical system.

In the following, we will write down our algorithms in terms of single qubit rotations and two qubits Ising interactions, since this is the language needed for the implementation of the algorithms, for example, in a liquid NMR QC. Again, had we used a different set of elementary gates our main results still hold but with modified quantum networks.

As an example of such decompositions, we consider the unitary operator $`U(t)=e^{iHt}`$, where $`H=\alpha \sigma _x^1\sigma _z^2\sigma _x^3`$ represents a time-independent Hamiltonian. After some simple calculations  we decompose $`U`$ into elementary gates (one qubit rotations and two qubits interactions) in the following way

$$U(t)=e^{i\alpha \sigma _x^1\sigma _z^2\sigma _x^3t}=e^{i\frac{\pi }{4}\sigma _y^3}e^{i\frac{\pi }{4}\sigma _z^1\sigma _z^3}e^{i\frac{\pi }{4}\sigma _x^1}e^{i\alpha \sigma _z^1\sigma _z^2t}e^{i\frac{\pi }{4}\sigma _x^1}e^{i\frac{\pi }{4}\sigma _z^1\sigma _z^3}e^{i\frac{\pi }{4}\sigma _y^3}.$$
(3)
This decomposition is shown in Fig. 1, where the quantum network representation is displayed. In the same way, we could also decompose an operator $`U^{}(t)=e^{i\alpha \sigma _y^1\sigma _z^2\sigma _y^3t}`$ using similar steps, by replacing $`\sigma _x^i\sigma _y^i`$ in the right hand side of Eq. 3.

## 3 simulation of fermionic systems

As discussed in the Introduction, quantum simulations require simulations of systems with diverse degrees of freedom and particle statistics. Fermionic systems are governed by Pauli’s exclusion principle, which implies that no more than one fermion can occupy the same quantum state at the same time. In this way, the Hilbert space of quantum states that represent a system of fermions in a solid is finite-dimensional ($`2^N`$ for spinless fermions, where $`N`$ is the number of sites or modes in the solid), and one could think in the existence of one-to-one mappings between the fermionic and Pauli’s spin-1/2 algebras. Similarly, any language which involves operators with a finite-dimensional representation (e.g., hard-core bosons, higher irreps of $`su(2)`$, etc.) can be mapped onto the standard model language .

In the second quantization representation, the (spinless) fermionic operators $`c_i^{}`$ ($`c_i^{}`$) are defined as the creation (annihilation) operators of a fermion in the $`i`$-th mode ($`i=1,\mathrm{},N`$). Due to the Pauli’s exclusion principle and the antisymmetric nature of the fermionic wave function under the permutation of two fermions, the fermionic algebra is given by the following commutation relations

$$\{c_i,c_j\}=0,\text{ }\{c_i^{},c_j\}=\delta _{ij}$$
(4)
where $`\{,\}`$ denotes the anticommutator.

The Jordan-Wigner transformation  is the isomorphic mapping that allows the description of a fermionic system by the standard model

$`c_j\left({\displaystyle \underset{l=1}{\overset{j1}{}}}\sigma _z^l\right)\sigma _{}^j`$ (5)
$`c_j^{}\left({\displaystyle \underset{l=1}{\overset{j1}{}}}\sigma _z^l\right)\sigma _+^j,`$ (6)
where $`\sigma _\mu ^i`$ are the Pauli operators defined in section 2. One can easily verify that if the operators $`\sigma _\mu ^i`$ satisfy the $`su(2)`$ commutation relations (Eq. 1), the operators $`c_i^{}`$ and $`c_i^{}`$ obey Eqs. 4.

We now need to show how to simulate a fermionic system by a QC. Just as for a simulation on a CC, the quantum simulation has three basic steps: the preparation of an initial state, the evolution of this state, and the measurement of a relevant physical property of the evolved state. We will now explain the first two steps, postponing the third until section 5.

### 3.1 Preparation of the initial state

In the most general case, any quantum state $`|\psi `$ of $`N_e`$ fermions can be written as a linear combination of Slater determinants $`|\varphi _\alpha `$

$$|\psi =\underset{\alpha =1}{\overset{L}{}}g_\alpha |\varphi _\alpha ,$$
(7)
where

$$|\varphi _\alpha =\underset{j=1}{\overset{N_e}{}}c_j^{}|\mathrm{𝗏𝖺𝖼}$$
(8)
with the vacuum state $`|\mathrm{𝗏𝖺𝖼}`$ defined as the state with no fermions. In the spin language, $`|\mathrm{𝗏𝖺𝖼}=|\mathrm{}`$.

We can easily prepare the states $`|\varphi _\alpha `$ by noticing that the quantum gate, represented by the unitary operator

$$U_m=e^{i\frac{\pi }{2}(c_m^{}+c_m^{})}$$
(9)
when acting on the vacuum state, produces $`c_m^{}|0`$ up to a phase factor. Making use of the Jordan-Wigner transformation (Eqs. 5, 6), we can write the operators $`U_m`$ in the spin language

$$U_m=e^{i\frac{\pi }{2}\sigma _x^m\underset{j=1}{\overset{m1}{}}\sigma _z^j}.$$
(10)
The successive application of $`N_e`$ similar unitary operators will generate the state $`|\varphi _\alpha `$ up to an irrelevant global phase.

A detailed preparation of the fermionic state $`|\psi =\underset{\alpha =1}{\overset{L}{}}g_\alpha |\varphi _\alpha `$ can be found in a previous work . The basic idea is to use $`L`$ extra (ancilla) qubits, then perform unitary evolutions controlled in the state of the ancillas, and finally perform a measurement of the $`z`$-component of the spin of the ancillas. In this way, the probability of successful preparation of $`|\psi `$ is $`1/L`$. (We need of the order of $`L`$ trials before a successful preparation.)

Another important case is the preparation of a Slater determinant in a different basis than the one given before

$$|\varphi _\beta =\underset{i=1}{\overset{N_e}{}}d_i^{}|\mathrm{𝗏𝖺𝖼},$$
(11)
where the fermionic operators $`d_i^{}`$’s are related to the operators $`c_j^{}`$ through the following canonical transformation

$$\stackrel{}{d}^{}=e^{iM}\stackrel{}{c}^{}$$
(12)
with $`\stackrel{}{d}^{}=(d_1^{},d_2^{},\mathrm{},d_N^{})`$, $`\stackrel{}{c}^{}=(c_1^{},c_2^{},\mathrm{},c_N^{})`$, and $`M`$ is an $`N\times N`$ Hermitian matrix. Making use of Thouless’s theorem , we observe that one Slater determinant evolves into the other, $`|\varphi _\beta =U|\varphi _\alpha `$, where the unitary operator $`U=e^{i\stackrel{}{c}^{}M\stackrel{}{c}}`$ can be written in spin operators using the Jordan-Wigner transformation and can be decomposed into elementary gates , as described in section 2. Since the number of gates scales polynomially with the system size, the state $`|\varphi _\beta `$ can be efficiently prepared from the state $`|\varphi _\alpha `$.

### 3.2 Evolution of the initial state

The second step in the quantum simulation is the evolution of the initial state. The unitary evolution operator of a time-independent Hamiltonian $`H`$ is $`U(t)=e^{iHt}`$. In general, $`H=K+V`$ with $`K`$ representing the kinetic energy and $`V`$ the potential energy. Since we usually have $`[K,V]0`$, the decomposition of $`U(t)`$, written in the spin language through the Jordan-Wigner transformation (Eqs. 5,6), in terms of elementary gates (one qubit rotations and two qubits interactions), becomes complicated. To avoid this problem, we instead use a Trotter decomposition, so the evolution during a short period of time ($`\mathrm{\Delta }t=t/𝒩`$ with $`\mathrm{\Delta }t0`$) is approximated. To order $`𝖮(\mathrm{\Delta }t)`$ (first order Trotter breakup)

$`U(t)`$ $`=`$ $`{\displaystyle \underset{g=1}{\overset{𝒩}{}}}U(\mathrm{\Delta }t),`$ (13)
$`U(\mathrm{\Delta }t)`$ $`=`$ $`e^{iH\mathrm{\Delta }t}=e^{i(K+V)\mathrm{\Delta }t}e^{iK\mathrm{\Delta }t}e^{iV\mathrm{\Delta }t}.`$ (14)
The potential energy $`V`$ is usually a sum of commuting diagonal terms, and the decomposition of $`e^{iV\mathrm{\Delta }t}`$ into elementary gates is straightforward. However, the kinetic energy $`K`$ is usually a sum of noncommuting terms of the form $`c_i^{}c_j^{}+c_j^{}c_i^{}`$ (bilinear fermionic operators), so we need again to perform a Trotter approximation of the operator $`e^{iK\mathrm{\Delta }t}`$. As an example of such a decomposition, we consider a typical term $`e^{i(c_i^{}c_j^{}+c_j^{}c_i^{})\mathrm{\Delta }t}`$ ($`i<j`$), when mapped onto the spin language gives

$$e^{\frac{i}{2}(\sigma _x^i\sigma _x^j+\sigma _y^i\sigma _y^j)\underset{k=i+1}{\overset{j1}{}}(\sigma _z^k)}=e^{\frac{i}{2}\sigma _x^i\sigma _x^j\underset{k=i+1}{\overset{j1}{}}(\sigma _z^k)}e^{\frac{i}{2}\sigma _y^i\sigma _y^j\underset{k=i+1}{\overset{j1}{}}(\sigma _z^k)}.$$
(15)
The decomposition of each term on the right hand side of Eq. 15 into elementary gates was already described in previous work . In section 2 and Fig. 1, we also showed an example of such a decomposition for $`i=1`$ and $`j=3`$. It is important to mention that the required number of elementary gates scales polynomially with the length $`|ji|`$. Notice that this step is not necessary for bosonic systems since no string of $`\sigma _z^k`$ operators is involved (see section 4).

The accuracy of this method increases as $`\mathrm{\Delta }t`$ decreases, so we might require a large number of gates to perform the evolution with small errors. To overcome this problem, one could use Trotter approximations of higher order in $`\mathrm{\Delta }t`$ .

### 3.3 Generalization: simulation of anyonic systems

The concepts described in sections 3.1 and 3.2 can be easily generalized to other more general particle statistics, namely hard-core anyons. By “hard-core”, we mean that only zero or one particle can occupy a single mode (Pauli’s exclusion principle).

The commutation relations between the anyonic creation and annihilation operators $`a_i^{}`$ and $`a_i`$, are given by

$`[a_i^{},a_j^{}]_\theta `$ $`=`$ $`[a_i^{},a_j^{}]_\theta =0,`$
$`[a_i^{},a_j^{}]_\theta `$ $`=`$ $`\delta _{ij}(1(e^{i\theta }+1)n_j),`$ (16)
$`[n_i,a_j^{}]`$ $`=`$ $`\delta _{ij}a_j^{},`$
($`ij`$) where $`n_j=a_j^{}a_j^{}`$, $`[\widehat{A},\widehat{B}]_\theta =\widehat{A}\widehat{B}e^{i\theta }\widehat{B}\widehat{A}`$, with $`0\theta <2\pi `$ defining the statistical angle. In particular, $`\theta =\pi `$ mod($`2\pi `$) corresponds to canonical spinless fermions, while $`\theta =0`$ mod($`2\pi `$) represents hard-core bosons.

In order to simulate this problem with a QC made of qubits, we need to apply the following isomorphic and efficient mapping between algebras

$`a_j^{}`$ $`=`$ $`{\displaystyle \underset{i<j}{}}[{\displaystyle \frac{e^{i\theta }+1}{2}}+{\displaystyle \frac{e^{i\theta }1}{2}}\sigma _z^i]\sigma _+^j,`$
$`a_j`$ $`=`$ $`{\displaystyle \underset{i<j}{}}[{\displaystyle \frac{e^{i\theta }+1}{2}}+{\displaystyle \frac{e^{i\theta }1}{2}}\sigma _z^i]\sigma _{}^j,`$ (17)
$`n_j`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\sigma _z^j),`$
where the Pauli operators $`\sigma _\mu ^j`$ where defined in section 2, and since they satisfy Eq. 1, the corresponding commutation relations for the anyonic operators (Eqs. 3.3) are satisfied, too.

We can now proceed in the same way as in the fermionic case, writing our anyonic evolution operator in terms of single qubit rotations and two qubits interactions in the spin-1/2 language. As we already mentioned, anyon statistics have fermion and hard-core boson statistics as limiting cases. We now relax the hard-core condition on the bosons.

## 4 Simulation of Bosonic systems

Quantum computation is based on the manipulation of quantum systems that possess finite number of degrees of freedom (e.g., qubits). From this point of view, the simulation of bosonic systems appears to be impossible, since the non existence of an exclusion principle implies that the Hilbert space used to represent bosonic quantum sates is infinite-dimensional; that is, there is no limit to the number of bosons that can occupy a given mode. However, sometimes we might be interested in simulating and studying properties such that the use of the whole Hilbert space is unnecessary, and only a finite sub-basis of states is sufficient. This is the case for physical systems with interactions given by the Hamiltonian

$$H=\underset{i,j=1}{\overset{N}{}}\alpha _{ij}b_i^{}b_j^{}+\beta _{ij}n_in_j,$$
(18)
where the operators $`b_i^{}`$ ($`b_i^{}`$) create (destroy) a boson at site $`i`$, and $`n_i=b_i^{}b_i^{}`$ is the number operator. The space dimension of the lattice is encoded in the parameters $`\alpha _{ij}`$ and $`\beta _{ij}`$. Obviously, the total number of bosons $`N_P`$ in the system is conserved, and we restrict ourselves to work with a finite sub-basis of states, where the dimension depends on the value of $`N_P`$.

The respective bosonic commutation relations (in an infinite-dimensional Hilbert space) are

$$[b_i,b_j]=0,[b_i,b_j^{}]=\delta _{ij}.$$
(19)
However, in a finite basis of states represented by $`\{|n_1,n_2,\mathrm{},n_N`$ with $`n_i=0,\mathrm{},N_P\}`$, where $`N_P`$ is the maximum number of bosons per site, the operators $`b_i^{}`$ can have the following matrix representation

$$\overline{b}_i^{}=1\mathrm{l}\mathrm{}1\mathrm{l}\underset{i^{th}\text{ factor}}{\underset{}{\widehat{b}^{}}}1\mathrm{l}\mathrm{}1\mathrm{l}$$
(20)
where $``$ indicates the usual tensorial product between matrices, and the $`(N_P+1)\times (N_P+1)`$ dimensional matrices $`1\mathrm{l}`$ and $`\widehat{b}^{}`$ are

$$1\mathrm{l}=\left(\begin{array}{ccccc}1& 0& 0& \mathrm{}& 0\\ 0& 1& 0& \mathrm{}& 0\\ 0& 0& 1& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& 1\end{array}\right)\text{ , }\widehat{b}^{}=\left(\begin{array}{cccccc}0& 0& 0& \mathrm{}& 0& 0\\ 1& 0& 0& \mathrm{}& 0& 0\\ 0& \sqrt{2}& 0& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& \sqrt{N_P}& 0\end{array}\right).$$
(21)
It is important to note that in this finite basis, the commutation relations of the bosons $`\overline{b}_i^{}`$ differ from the standard bosonic ones (Eq. 19)

$$[\overline{b}_i^{},\overline{b}_j^{}]=0,\text{ }[\overline{b}_i^{},\overline{b}_j^{}]=\delta _{ij}\left[1\frac{N_P+1}{N_P!}(\overline{b}_i^{})^{N_P}(\overline{b}_i^{})^{N_P}\right],$$
(22)
and clearly $`(\overline{b}_i^{})^{N_P+1}=0`$.

As we mentioned in the Introduction, our idea is to simulate any physical system in a QC made of qubits. For this purpose, we need to map the bosonic algebra into the spin-1/2 language. However, since Eqs. 22 imply that the linear span of the operators $`\overline{b}_i^{}`$ and $`\overline{b}_i^{}`$ is not closed under the bracket (commutator), a direct mapping between the bosonic algebra and the spin-1/2 algebra (such as the case of the Jordan-Wigner transformation between the fermionic and spin-1/2 algebra) is not possible. Therefore, we could think in a one-to-one mapping between the bosonic and spin-1/2 quantum states, instead of an isomorphic mapping between algebras. Let us show a possible mapping of quantum states.

We start by considering only the $`i`$-th site in the chain. Since this site can be occupied with at most $`N_P`$ bosons, it is possible to associate an $`N_P+1`$ qubits quantum state to each particle number state, in the following way

$`|0_i`$ $``$ $`|_0_1_2\mathrm{}_{N_P}_i`$
$`|1_i`$ $``$ $`|_0_1_2\mathrm{}_{N_P}_i`$
$`|2_i`$ $``$ $`|_0_1_2\mathrm{}_{N_P}_i`$ (23)
$`\mathrm{}`$ $`\mathrm{}`$
$`|N_P_i`$ $``$ $`|_0_1_2\mathrm{}_{N_P}_i`$
where $`|n_i`$ denotes a quantum state with $`n`$ bosons in site $`i`$. Therefore, we need $`N(N_P+1)`$ qubits for the simulation (where $`N`$ is the number of sites). In Fig. 2 we show an example of this mapping for a quantum state with 7 bosons in a chain of 5 sites.

By definition (see Eqs. 20, 21) $`\overline{b}_i^{}|n_i=\sqrt{n+1}|n+1_i`$, so the operator

$$\overline{b}_i^{}=\underset{n=0}{\overset{N_P1}{}}\sqrt{n+1}\sigma _{}^{n,i}\sigma _+^{n+1,i},$$
(24)
where the pair $`(n,i)`$ indicates the qubit $`n`$ that represents the $`i`$-th site, acts in the $`N_P+1`$ qubits states of Eqs. 4 as $`\overline{b}_i^{}|_0\mathrm{}_{n1}_n_{n+1}\mathrm{}_{N_P}_i=\sqrt{n+1}|_0\mathrm{}_n_{n+1}_{n+2}\mathrm{}_{N_P}_i`$. Then, its matrix representation in this basis is the same matrix representation of $`b_i^{}`$ in the basis of bosonic states. Similarly, the number operator can be written

$$\overline{n}_i=\underset{n=0}{\overset{N_P}{}}n\frac{\sigma _z^{n,i}+1}{2},$$
(25)
and act as $`\overline{n}_i|_0\mathrm{}_{n1}_n_{n+1}\mathrm{}_{N_P}_i=n|_0\mathrm{}_n_{n+1}_{n+2}\mathrm{}_{N_P}_i`$. Notice that $`[\overline{b}_i^{},_{n=0}^{N_P}\sigma _z^{n,i}]=0`$, which means that these operators conserve the total $`z`$-component of the spin and, thus, always keep states within the same subspace.

We can now write down the Hamiltonian in Eq. 18 in the spin-1/2 algebra as

$$H=\underset{i,j=1}{\overset{N}{}}\alpha _{ij}\overline{b}_i^{}\overline{b}_j^{}+\beta _{ij}\overline{n}_i\overline{n}_j,$$
(26)
where the operators $`\overline{b}_i^{}`$ ($`\overline{b}_i^{}`$) are given by Eq. 24 and $`\overline{n}_i`$ bt Eq. 25. In this way, we are able to obtain physical properties of the bosonic system (such as the mean value of an observable, the mean value of the evolution operator, etc.) in a QC made of qubits. It is important to note that the type of Hamiltonian given by Eq. 18 is not the only one that can be simulatable using the described method. The only constraint is a fixed maximum number of bosons per site (or mode).

### 4.1 Preparation of the initial state

As in the fermionic case, the most general bosonic state of an $`N`$ sites quantum system with $`N_P`$ bosons can be written as a linear combination of product states like

$$|\varphi _\alpha =𝖪(b_1^{})^{n_1}(b_2^{})^{n_2}\mathrm{}(b_N^{})^{n_N}|\mathrm{𝗏𝖺𝖼},$$
(27)
where $`𝖪`$ is a normalization factor, $`n_i`$ is the number of bosons at site $`i`$ ($`\underset{i=1}{\overset{N}{}}n_i=N_P`$), and $`|\mathrm{𝗏𝖺𝖼}`$ is the boson vacuum state (no particle state). Using the mapping described in Eq. 4, we can write the vacuum state in the spin language as $`|\mathrm{𝗏𝖺𝖼}=|_0_1\mathrm{}_{N_P}_1\mathrm{}|_0_1\mathrm{}_{N_P}_N`$ and $`|\varphi _\alpha =|_0\mathrm{}_{n_1}\mathrm{}_{N_P}_1\mathrm{}|_0\mathrm{}_{n_N}\mathrm{}_{N_P}_N`$ (see Fig. 2 for an example). Therefore, the preparation of $`|\varphi _\alpha `$ in a QC made of qubits is an easy process: only $`N`$ spins are flipped from the fully polarized state, where all spins are pointing down.

The preparation of a bosonic initial state of the form $`|\psi =\underset{\alpha =1}{\overset{L}{}}g_\alpha |\varphi _\alpha `$ is realized as in the fermionic case. Again, we need to add $`L`$ ancillas (extra qubits), perform controlled evolutions on their states, and finally perform a measurement of an spin component .

### 4.2 Evolution of the initial state

The basic idea is to use the first order Trotter approximation (see the fermionic case) to separate those terms of the Hamiltonian that belong to the kinetic energy $`K`$, from the ones that belong to the potential energy $`V`$ ($`H=K+V`$, $`[K,V]0`$), i.e.,

$$e^{iH\mathrm{\Delta }t}e^{iK\mathrm{\Delta }t}e^{iV\mathrm{\Delta }t}.$$
(28)
In general, $`K`$ is a sum of non commuting terms of the form $`b_k^{}b_l^{}+b_l^{}b_k^{}`$, and we need to perform another first order Trotter approximation to decompose it into elementary gates (in the spin language). Then, a typical term $`e^{i(b_i^{}b_j^{}+b_j^{}b_i^{})t}`$ when mapped onto the spin language (Eq. 24) gives

$`\mathrm{exp}[{\displaystyle \frac{it}{8}}{\displaystyle \underset{n,n^{}=0}{\overset{N_P1}{}}}\sqrt{(n+1)(n^{}+1)}[(\sigma _x^{n,i}\sigma _x^{n+1,i}+\sigma _y^{n,i}\sigma _y^{n+1,i})(\sigma _x^{n^{},j}\sigma _x^{n^{}+1,j}+\sigma _y^{n^{},j}\sigma _y^{n^{}+1,j})`$
$`+(\sigma _x^{n,i}\sigma _y^{n+1,i}\sigma _y^{n,i}\sigma _x^{n+1,i})(\sigma _x^{n^{},j}\sigma _y^{n^{}+1,j}\sigma _y^{n^{},j}\sigma _x^{n^{}+1,j})]],`$ (29)
where $`N_P`$ is the number of bosons. The terms in the exponent of Eq. 4.2 commute with each other, so the decomposition into elementary gates becomes straightforward. As an example (see Fig.3), we consider a system of two sites with one boson. We need then $`2(1+1)=4`$ qubits for the simulation, and Eq. 24 implies that $`\overline{b}_1^{}=\sigma _{}^{0,1}\sigma _+^{1,1}`$ and $`\overline{b}_2^{}=\sigma _{}^{0,2}\sigma _+^{1,2}`$. Then, $`e^{i(b_i^{}b_j^{}+b_j^{}b_i^{})t}`$ becomes

$`\mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _x^{0,1}\sigma _x^{1,1}\sigma _x^{0,2}\sigma _x^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _x^{0,1}\sigma _x^{1,1}\sigma _y^{0,2}\sigma _y^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _y^{0,1}\sigma _y^{1,1}\sigma _x^{0,2}\sigma _x^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _y^{0,1}\sigma _y^{1,1}\sigma _y^{0,2}\sigma _y^{1,2})`$ (30)
$`\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _y^{0,1}\sigma _x^{1,1}\sigma _y^{0,2}\sigma _x^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _y^{0,1}\sigma _x^{1,1}\sigma _x^{0,2}\sigma _y^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _x^{0,1}\sigma _y^{1,1}\sigma _y^{0,2}\sigma _x^{1,2})\times \mathrm{exp}({\displaystyle \frac{it}{8}}\sigma _x^{0,1}\sigma _y^{1,1}\sigma _x^{0,2}\sigma _y^{1,2}),`$
where the decomposition of each of the terms in Eq. 30 in elementary gates can be done using the methods described in previous works . In particular, in Fig. 3 we show the decomposition of the term $`\mathrm{exp}\left(\frac{i}{8}\sigma _x^{0,1}\sigma _y^{1,1}\sigma _y^{0,2}\sigma _x^{1,2}t\right)`$, where the qubits were relabeled as $`(n,i)n+2i1`$ (e.g., $`(0,1)1`$).

On the other hand, it is important to mention that the number of operations involved in the decomposition is not related to the distance between the sites $`i`$ and $`j`$, as in the fermionic case.

## 5 Measurement: Correlation functions and Energy spectra

In previous work  we introduced an efficient algorithm for the measurement of correlation functions in quantum systems. The idea is to make an indirect measurement, that is, we prepare an ancilla qubit (extra qubit) in a given initial state, then interact with the system whose properties one wants to measure, and finally we measure some observable of the ancilla to obtain information about the system. Particularly, we could be interested in the measurement of dynamical correlation functions of the form

$$G(t)=\psi |T^{}A_i^{}TB_j\psi $$
(31)
where $`A_i`$ and $`B_j`$ are unitary operators (any operator can be decomposed in a unitary operator basis as $`A=\underset{i}{}\alpha _iA_i`$, $`B=\underset{j}{}\beta _jB_j`$), $`T=e^{iHt}`$ is the time evolution operator of a time-independent Hamiltonian $`H`$, and $`|\psi `$ is the state of the system whose correlations one wants to determine. If we were interested in the evaluation of spatial correlation functions, we would replace the evolution operator $`T`$ by the space translation operator. In Fig. 4 we show the quantum algorithm (quantum network) for the evaluation of $`G(t)`$. As explained before , the initial state (ancilla plus system) has to be prepared in the quantum state $`|+_𝖺|\psi `$ (where $`𝖺`$ denotes the ancilla qubit and $`|+=\frac{|0+|1}{\sqrt{2}}`$). Additionally, we have to perform an evolution (unitary operation) in the following three steps: i) a controlled evolution in the state $`|1`$ of the ancilla $`\text{C-B}=|00|I+|11|B_j`$, ii) a time evolution $`T`$, and iii) a controlled evolution in the state $`|0`$ of the ancilla $`\text{C-A}=|00|A_i+|11|I`$. Finally we measure the observable $`2\sigma _+^𝖺=\sigma _x^𝖺+i\sigma _y^𝖺=G(t)`$.

On the other hand, sometimes we are interested in obtaining the spectrum (eigenvalues) of a given observable $`\widehat{Q}`$ (i.e., an hermitian operator). A quantum algorithm (network) for this purpose was also given in previous work . Again, the basic idea is to perform an indirect measurement using an extra qubit (see Fig. 5). Basically, we prepare the initial state (ancilla plus system) $`|+_𝖺|\varphi `$, then apply the evolution $`e^{i\widehat{Q}\sigma _z^𝖺\frac{t}{2}}`$, and finally measure the observable $`2\sigma _+^𝖺(t)=\varphi |e^{i\widehat{Q}t}\varphi `$. Since the initial state of the system can be written as a linear combination of eigenstates of $`\widehat{Q}`$, $`|\varphi =\underset{n=0}{\overset{L}{}}\gamma _n|\psi _n`$, where $`\gamma _n`$ are complex coefficients and $`|\psi _n`$ are eigenstates of $`\widehat{Q}`$ with eigenvalue $`\lambda _n`$, the classical Fourier transform applied to the function of time $`2\sigma _+^𝖺(t)`$ gives us $`\lambda _n`$

$$\widehat{F}(\lambda )=\underset{n=0}{\overset{L}{}}2\pi |\gamma _n|^2\delta (\lambda \lambda _n).$$
(32)
Without loss of generality, we can choose $`\widehat{Q}=H`$, with $`H`$ some particular Hamiltonian.

It is important to note that in order to obtain the different eigenvalues of $`\widehat{Q}`$, the overlap between the initial state and the eigenstates of $`\widehat{Q}`$ must be different from zero. One can use different mean-field solutions of $`\widehat{Q}`$ as initial states $`|\varphi `$ depending on the part of the spectrum one wants to determine with higher accuracy.

## 6 Algorithm efficiency and errors

An algorithm is considered efficient if the number of operations involved scales polynomially with the system size, and if the effort required to make the error $`ϵ`$ in the measurement of a relevant property smaller, scales polynomially with $`1/ϵ`$.

While the evolution step involves a number of unitary operations that scales polynomially with the system size (such is the case for the Trotter approximation) whenever the Hamiltonian $`H`$ is physical (e.g., is a sum of a number of terms that scales polynomially with the system size), the preparation of the initial state could be inefficient. Such inefficiency would arise, for example, if the state $`|\psi `$ defined in Eq. 8 or Eq. 27 is a linear combination of an exponential number of states ($`Lx^N`$, with $`N`$ the number of sites in the system and $`x`$ a positive number). However, if we assume that $`|\psi `$ is a finite combination of states ($`L`$ scales polynomially with $`N`$), its preparation can be done efficiently. (Any (Perelomov-Gilmore) generalized coherent state can be prepared in a number of steps that scales polynomially with the number of generators of the respective algebra.) On the other hand, the measurement process described in section 5 is always an efficient step, since it only involves the measurement of the spin of one qubit, despite the number of qubits or sites $`N`$ of the quantum system.

Errors $`ϵ`$ come from gate imperfections, the use of the Trotter approximation in the evolution operator, and the statistics in measuring the spin of the ancilla qubit (sections 3.2, 4.2, and 5). A precise description and study of the error sources can be found in previous work . The result is that the algorithms described here, for the simulation of physical systems and processes, are efficient if the preparation of the initial state is efficient, too.

## 7 conclusions

We studied the implementation of quantum algorithms for the simulation of an arbitrary quantum physical system on a QC made of qubits, making a distinction between systems that are governed by Pauli’s exclusion principle (fermions, hard-core bosons, anyons, spins, etc.), and systems that are not (e.g, canonical bosons). For the first class of quantum systems, we showed that a mapping between the corresponding algebra of operators and the spin-1/2 algebra exists, since both have a finite-dimensional representation. On the other hand, the operator representation of quantum systems that are not governed by an exclusion principle is infinite-dimensional, and an isomorphic mapping to the spin-1/2 algebra is not possible. However, one can work with a finite set of quantum states, setting a constraint, such as fixing the number of bosons in the system. Then, the representation of bosonic operators becomes finite-dimensional, and we showed that we can write down bosonic operators in the spin-1/2 language (Eq. 24), mapping bosonic states to spin-1/2 states (Eq. 4) .

We also showed how to perform quantum simulations in a QC made of qubits (quantum networks), giving algorithms for the preparation of the initial state, the evolution, and the measurement of a relevant physical property, where in the most general case the unitary operations have to be approximated (sections 3.2,4.2).

The mappings explained are efficient in the sense that we can perform them in a number of operations that scales polynomially with the system size. This implies that the evaluation of some correlation functions in quantums states that can be prepared efficiently is also efficient, showing an exponential speed-up of these algorithms with respect to their classical simulation. However, these mappings are insufficient to establish that quantum networks can simulate any physical problem efficiently. As we mentioned in the introduction, this is the case for the determination of the spectrum of the Hamiltonian in the two-dimensional Hubbard model , where the signal-to-noise ratio decays exponentially with the system size.

Finally, in Fig. 6 a table displays the advantages of simulating some known algorithms with a QC than with a CC, concluding that QCs behave as efficient devices for some quantum simulations.