Wrong matrix element of a string of creation and annihilation operators

[udevd]

Since the recent refactor, the function coefficient() of manybody.jl apparently does not calculate the right matrix elements. The issue seems to be that in the following portion of the code the vm/vn can only decrement by 1. If the string of annihilation and creation operators contains repeated entries, this is not the right behavior:

    for i in 1:length(occ_m)
        vm = occ_m[i]
        vn = occ_n[i]
        i in at_indices && (vm -= 1) # <- vm should decrease by the number of times i appears in at_indices, not 1
        i in a_indices && (vn -= 1) # <- similar issue
        vm != vn && return zero(C)
    end

The following code should yield a coefficient of 2 (it corresponds to ⟨2,0|a₁⁺a₁⁺ a₂a₂|0,2⟩), but 0 is returned:

init_fock=[0,2]
end_fock=[2,0]
annihilation_indices=[2,2]
creation_indices=[1,1]
coefficient(end_fock,init_fock,creation_indices,annihilation_indices)

The bug propagates further to manybodyoperator(), leading to entirely wrong results. It appears in version 1.0.14 of QuantumOptics.jl.

[Krastanov]

@udevd thank you for tracking this down, it is very much appreciated!

0.4.9 works fine but the problem appears in 0.4.10

@aryavorskiy it seems this happened due to the optimizations in #121 (we did not have a test case to detect this issue). Would you be able to submit a fix and a unit test in the near future? If it would take you some time (which is perfectly reasonable), I would prefer to reverse that commit, back to the much slower version that works correctly for the edge case.

[Krastanov]

@udevd , there is a fix in #128 -- it will be merged and released soon (together with another significant performance improvement thanks to Alexander)

[udevd]

Thanks! I'm not using the function directly, was just trying out the many body functionality or QuantumOptics.jl. The results of manybodyoperator() were inconsistent with my intuitions so I went digging through the code looking for the cause.

[Krastanov]

0.4.18 will be registered in the next few minutes. Could you give it a try to confirm this is indeed fixed?

JuliaRegistries/General#92348

[udevd]

Yes, it does work properly now. Checked the output of manybodyoperator() against the matrix elements (calculated explicitly in the full Hilbert space) ⟨Dⁿₖ|O|Dⁿₗ⟩ for all symmetric two-body operators O and Dicke states |Dⁿₖ⟩ for up to 10 qubits and they are the same, barring some numerical noise.

[Krastanov]

Wonderful! If you feel like it, do not hesitate to paste the code you used for tests here and we can add it to the test suite. Or feel free to send a pull request if that is easier.

[udevd]

Sure! I'll try to find the time to tidy up the code, was really busy for the last two weeks.

[udevd]

So, below is the code – I'm not sure how to perform the pull request in a manner consistent with your local style in this repo, so it's easier this way. It's far from optimal (and since it has to do the full n-qubit space calculations, it possibly does not matter much), but it does work. Feel free to extract from it what you need.

using QuantumOptics
using Combinatorics
using Printf
"""
Calculate the tensor product of operators described by opsdict in a multipartite space. 
For indices which are not the keys of opsdict, the default operator is used (e.g. identity).

For example, operators_at(Dict(1=>x,2=>y), 3, id) is equal to tensor(x,y,id).
"""
function operators_at(opsdict,nsubsystems,default)
    ops=[get(opsdict,i,default) for i in 1:nsubsystems]
    return tensor(ops...)
end

"""
Calculate the symmetrized multiqubit operator described by ops. 
Roughly speaking, if ops==[op1,op2,...], the output is a sum of operators_at(Dict(idx1=>op1,idx2=>op2,...),...), summed over all pairwise different [idx1,idx2,...]. 
Note: it does not check whether the ops==[op1,op2,...] are pairwise different. 
Thus, symmetric_nqubit_operator(2,[x,x])=2*tensor(x,x),
 since both permutations [idx1,idx2]=[1,2] and [2,1] are considered.
This behavior seems to be consisent with the results of manybodyoperator(...).
"""
function symmetric_nqubit_operator(nqubits::Integer,ops)
    b = SpinBasis(1//2)
    identity=identityoperator(b)
    subsystemsize=length(ops)
    indices=multiset_permutations(1:nqubits,subsystemsize)
    basis=[]
    result=sum(indices) do ids
    	opsdict=Dict(zip(ids,ops))
        operators_at(opsdict,nqubits,identity)
    end
    return result
end
"""
Determines the matrix elements of an operator given a basis (useful for a subspace basis).
"""
function matrix_elements(operator,basis)
    @assert operator.basis_l==operator.basis_r
    result = DenseOperator(basis)
    N=length(basis)
    for m=1:N, n=1:N
        result.data[m,n]+=dagger(basis.basisstates[m])*operator*basis.basisstates[n]
    end
    return result
end
"""
Calculates the Dicke state D(nqubits,num) of a multiqubit system -- the symmetric state built as a (normalized) sum of kets, each with (num) number of spin down/0 states and (nqubits-num) spin up/1 states. 
For example, the Dicke state D(3,1) is equal to (|100⟩+|010⟩+|001⟩)/√3.
"""
function dicke_state(nqubits,num)
    b=SpinBasis(1//2)
    basis=vcat(fill(spinup(b),nqubits-num),fill(spindown(b),num))
    v=sum(map(p->tensor(p...),multiset_permutations(basis,nqubits)))
    v/norm(v)
end

"""
Check whether the matrix elements of an symmetric operator in the full n-qubit state space 
between the Dicke states (which span the fully symmetric space) 
and the operator built using manybodyoperator(...) do match.
"""
function check_2qubit_ops(nqubits,eps=10^-10)
    b=SpinBasis(1//2)
    x=sigmax(b)
    y=sigmay(b)
    z=sigmaz(b)
    pauli_dict=Dict("x"=>x,"y"=>y,"z"=>z)
    dickevecs=[dicke_state(nqubits,m) for m in 0:nqubits];
    dickebasis=SubspaceBasis(dickevecs)
    bosonbasis=ManyBodyBasis(b,bosonstates(b,nqubits))
    for s1 in ["x","y","z"], s2 in ["x","y","z"] 
    	op1=pauli_dict[s1]
    	op2=pauli_dict[s2]
        matrix_explicit=matrix_elements(symmetric_nqubit_operator(nqubits,(op1,op2)),dickebasis)
        matrix_manybody=manybodyoperator(bosonbasis,tensor(op1,op2))|>dense
        matdiff=DenseOperator(bosonbasis)
        matdiff.data=matrix_explicit.data-matrix_manybody.data
        if !all(element->abs(element)<eps,matdiff.data)
             @printf("Wrong matrix elements for %d-qubit symmetric operator built from [%s,%s].\n",nqubits,s1,s2)
        end
    end
end

check_2qubit_ops(10)