Computing fidelity of quantum state

I write my program, however, i cannot plot the correct figure. When i use the ket psi_0 as initial state, i can plot the correct figure same Matlab. However, when i use the density  psi_0⊗dagger(psi_0), i cannot plot the same figure as Matlab. Further i cannot use the command ‘fidelity(rho, sigma)’ to plot the figure when i use ket as initial state.
This is my program

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms

# baisis
b = NLevelBasis(dim) # Basis for a single atom
s00 = nlevelstate(b,1) # |0> state
sD = nlevelstate(b,2) # |1> state
srr = nlevelstate(b,3) # |r> state
# Identity4 = identityoperator(b) # Identity operator  one or identityoperator

# parameters
Omega = 1;
Ud = 0.5;
delta = 0.2;
tf = 20*pi/Omega;

Hami = sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD) + dagger(sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD)) + Ud*srr ⊗ dagger(srr) - delta*sD ⊗ dagger(sD);

# initial and final state
psi_0 = s00;
psi_f = s00;

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

# evolution
tlist_pump = range(0, stop=tf, length=1000);

exp_psi0 = [];

tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami);
# exp_psi0 = real(expect(rho_final, rho_XZ))
exp_psi0 = [fidelity(rhoXZ, rho_final) for rhoXZ in rho_XZ]

plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(exp_psi0)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()


When i use ket psi_0 as initial state (i believe it is correct, same as Matlab), i get

while i use density dm0 = psi_0⊗dagger(psi_0) as initial state (i believe it is wrong, not same as Matlab), i get

I can’t figure out why？

Furthermore, i also consider fidelity, i use exp_psi0 = real(expect(rho_final, rho_XZ)). First, it cannot use rho_final of type ket. Second, i think its figure is also wrong currently:

Could you elaborate on the solver+plotting commands you use for each of the three plots?

If I understand correctly:

1. The first plot is timeevolution.schroedinger with a Ket initial state. But then I do not understand what is being plotted, because fidelity does not have methods that work with Ket types.

2. The second plot is timeevolution.schroedinger with an Operator initial state which is simply the outer product of the initial Ket from case 1. It is actually the code you have shared and it plots fidelity(...).

3. The third plot is timeevolution.schroedinger with an Operator initial state which is simply the outer product of the initial Ket from case 1. It plots real(expect(...)) instead of fidelity(...).

I guess it makes sense that 2 and 3 are different. In 2 you plot F(ρ, σ) = Tr\left(\sqrt{\sqrt{ρ}σ\sqrt{ρ}}\right) while in 3 you plot E(ρ, σ) = Tr(ρ, σ^\dagger) which are not the same unless you have some commutativity property fulfilled (and even then, they will be not equal, rather one of them will be the square of the other).

I am not sure why case 1 does not match, but mostly because I am not sure what case 1 is.

Sorry, my description may not be very clear to you. Here, I will explain my question in detail. (I also slightly changed the parameter values to make the different in the figure more obvious.)

My first code is

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms

# baisis
b = NLevelBasis(dim) # Basis for a single atom
s00 = nlevelstate(b,1) # |0> state
sD = nlevelstate(b,2) # |1> state
srr = nlevelstate(b,3) # |r> state
# Identity4 = identityoperator(b) # Identity operator  one or identityoperator

# parameters
Omega = 1;
Ud = 10;
delta = 0;
tf = pi/Omega;

Hami = sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD) + dagger(sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD)) + Ud*srr ⊗ dagger(srr) - delta*sD ⊗ dagger(sD);

# initial and final state
psi_0 = s00;
psi_f = s00;

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

# evolution
tlist_pump = range(0, stop=tf, length=1000);

ffidelity = [];

tout, psi_XZ = timeevolution.schroedinger(tlist_pump, psi_0, Hami);
ffidelity = real(expect(rho_final, psi_XZ))

# plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(ffidelity)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()


Here I write tout, rho_XZ = timeevolution.schroedinger(tlist_pump, psi_0, Hami) and use psi_0 (ket) as intial state. The fidelity is defined as exp_psi0 = real(expect(rho_final, rho_XZ)). I can plot

Here because I cannot find the source code of expect , i guess it can be used as |\langle \psi_{XZ}|\rho_{final}|\psi_{XZ}\rangle| (here \rho_{final}=|\psi_{final}\rangle\langle \psi_{final}|). And the image is same as that in Matlab (which does use the definition abs(psi_final'*rho_xz*psi_final)) of fidelity. I want to know is the definition of expect |\langle \psi_{XZ}|\rho_{final}|\psi_{XZ}\rangle|?

My second code is

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms

# baisis
b = NLevelBasis(dim) # Basis for a single atom
s00 = nlevelstate(b,1) # |0> state
sD = nlevelstate(b,2) # |1> state
srr = nlevelstate(b,3) # |r> state
# Identity4 = identityoperator(b) # Identity operator  one or identityoperator

# parameters
Omega = 1;
Ud = 10;
delta = 0;
tf = pi/Omega;

Hami = sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD) + dagger(sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD)) + Ud*srr ⊗ dagger(srr) - delta*sD ⊗ dagger(sD);

# initial and final state
psi_0 = s00;
psi_f = s00;

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

# evolution
tlist_pump = range(0, stop=tf, length=1000);

ffidelity = [];

tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami);
ffidelity = real(expect(rho_final, rho_XZ))

# plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(ffidelity)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()


Here tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami) I use initial state dm0 (density operator) and write real(expect(rho_final, rho_XZ)) to compute fidelity. And I plot

I hope here expect acts as trace(\rho_{final}\rho_{XZ}) = \langle \psi_{final}|\rho_{XZ}|\psi_{final}\rangle . However, it must not. Because it is different as that from my first code.

My third code is

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms

# baisis
b = NLevelBasis(dim) # Basis for a single atom
s00 = nlevelstate(b,1) # |0> state
sD = nlevelstate(b,2) # |1> state
srr = nlevelstate(b,3) # |r> state
# Identity4 = identityoperator(b) # Identity operator  one or identityoperator

# parameters
Omega = 1;
Ud = 10;
delta = 0;
tf = pi/Omega;

Hami = sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD) + dagger(sqrt(2)*Omega*sD ⊗ dagger(s00) + sqrt(2)*Omega*srr ⊗ dagger(sD)) + Ud*srr ⊗ dagger(srr) - delta*sD ⊗ dagger(sD);

# initial and final state
psi_0 = s00;
psi_f = s00;

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

# evolution
tlist_pump = range(0, stop=tf, length=1000);

ffidelity = [];

tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami);
ffidelity = [fidelity(rhoXZ, rho_final) for rhoXZ in rho_XZ]

# plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(ffidelity)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()


Here I use tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami) and  ffidelity = [fidelity(rhoXZ, rho_final) for rhoXZ in rho_XZ] with initial state dm0 (density operator). And I get

I believe it is the definition of fidelity, however, I do not want to use trace(\sqrt{\rho*\sigma}). Instead, I want to use trace(\rho*\sigma). I cannot find a direct way to define fidelity like this.

conclusion

1. I do not know the knowledge and source code of code expect.
My ideal guess is when I use initial state with ket, it is \langle \psi_{target}|\rho(t)|\psi_{target}\rangle, and when I use initial state with density operator, it is trace(\rho(t)*\rho_{target}).

2. Compared to fidelity(rho,sigma) with definition trace(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}), I want to use trace(\rho(t)*\rho_{target}) and \langle \psi_{target}|\rho(t)|\psi_{target}\rangle. It is precisely because there is no such two definitions and I cannot directly write it myself, so I used expect. But I find when I use initial state with density operator, i cannot get the correct image (using expect) same as Matlab.

3. Do you have any way to write directly the fidelity in the way trace(\rho(t)*\rho_{target}) and \langle \psi_{target}|\rho(t)|\psi_{target}\rangle?

The issue here is this: timeevolution.schroedinger doesn’t work they way you think with operators. It always solves the equation of the form \frac{d}{dt}X = -i H X. Whether you now supply |\psi\rangle \equiv X or \rho \equiv X, doesn’t matter, the solver does the same thing. However, in the latter case, you end up with a non-Hermitian matrix (not a state!). Rather, this is used to calculate the propagate, i.e. you can obtain the unitary U(t) that evolves a pure state, |\psi(t)\rangle = U(t)|\psi(0)\rangle .

What you actually want is to solve the von Neumann equation, which is equivalent to a Master equation without any jump operators. Taking your second coding example and just changing the line

tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami);


to

tout, rho_XZ = timeevolution.master(tlist_pump, dm0, Hami, []);


produces the same output as you get in the first case.

I didn’t check for the third case, but it should be the same. That and the square root factor in the fidelity as @Krastanov explained should make up the difference.

1 Like

You can find documentation about these functions in a couple of ways:

• Press ? in the Julia REPL and type out the function name
julia> using QuantumOptics

help?> fidelity
search: fidelity avg_gate_fidelity

fidelity(rho, sigma)

Fidelity of two density operators.

The fidelity of two density operators ρ and σ is defined by

F(ρ, σ) = Tr\left(\sqrt{\sqrt{ρ}σ\sqrt{ρ}}\right),

where \sqrt{ρ}=\sum_n\sqrt{λ_n}|ψ⟩⟨ψ|.

help?> expect
search: expect onebodyexpect ExponentialBackOff

expect(index, op, state)

If an index is given, it assumes that op is defined in the subsystem
specified by this number.

────────────────────────────────────────────────────────────────────────────

expect(op, state)

Expectation value of the given operator op for the specified state.

state can either be a (density) operator or a ket.


• The various documentation panels built in tools like Pluto.jl and VSCode (with the julia plugin), depending on what you are using

• The website for the package API · QuantumOptics.jl

From these we can see that expect and fidelity as defined in this package are very different, so that should answer your questions about differences between point 2 and 3.

If you want to look at the code, it can be found on the public repository for the package, but it can also be seen by clicking the “source” button in the online documentation I shared (bottom right corner of the panel with the function description) or in the Julia REPL with the edit function and macro (and you can even directly modify the code yourself in that manner). Here is the link to see it: QuantumOpticsBase.jl/src/metrics.jl at v0.4.20 · qojulia/QuantumOpticsBase.jl · GitHub

But that does not explain why you have differences in your case 1 and case 2.

There is a much bigger semantic misconception that leads to this issue. You already know that schroedinger(tspan, psi::Ket, H::Operator) gives you the Schroedinger evolution of the ket psi under the Hamiltonian H.

However, as you can see in the documentation, schroedinger(tspan, psi::Operator, H::Operator) does NOT interpret the matrix psi as a density matrix, rather it interprets it as a “each column is a separate state (ket) that we want to evolve”. This is a useful choice if one wants to do optimal control or to compute time evolution operators. This function is for getting a unitary operator, not a density operator.

1 Like

Yeah, I get what you mean. In fact, I misunderstood timeevolution.schroedinger. As you said, I understand it as a density operator. According to you suggestion, I have get the correct result. Thank you very much!

Thank you very much! I have get the correct result. In fact, I can find the source code, however, it only has Languages description about expect, so I think my question stems from an unknown expression for the expect (unlike fidelity which has mathematical formula that can be intuitively seen), now I know I misunderstand the timeevolution.schroedinger. Thank for your help again!