# Computing fidelity of quantum state
**URL:** https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264
**Category:** Quantum
**Tags:** question, package, plotting
**Created:** [June 6, 2024, 2:29am UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264 "2024-06-06T02:29:45Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@fqguo](https://discourse.julialang.org/u/fqguo)
#### Post date: [June 6, 2024, 2:29am UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/1 "2024-06-06T02:29:45Z")
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
```julia
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:
---
### Author:  [@Krastanov](https://discourse.julialang.org/u/Krastanov)
#### Post date: [June 6, 2024, 3:16am UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/2 "2024-06-06T03:16:20Z")
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.
---
### Author:  [@fqguo](https://discourse.julialang.org/u/fqguo)
#### Post date: [June 6, 2024, 6:45am UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/3 "2024-06-06T06:45:12Z")
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
```julia-auto
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
```julia-auto
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
```julia-auto
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?
---
### Author:  [@david-pl](https://discourse.julialang.org/u/david-pl)
#### Post date: [June 6, 2024, 1:18pm UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/4 "2024-06-06T13:18:48Z")
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
```julia
tout, rho_XZ = timeevolution.schroedinger(tlist_pump, dm0, Hami);
```
to
```julia
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.
---
### Author:  [@Krastanov](https://discourse.julialang.org/u/Krastanov)
#### Post date: [June 6, 2024, 2:01pm UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/5 "2024-06-06T14:01:10Z")
You can find documentation about these functions in a couple of ways:
- Press `?` in the Julia REPL and type out the function name
```julia
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](https://docs.qojulia.org/api/#QuantumOpticsBase.fidelity)
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](https://github.com/qojulia/QuantumOpticsBase.jl/blob/v0.4.20/src/metrics.jl#L176-L188)
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.**
---
### Author:  [@fqguo](https://discourse.julialang.org/u/fqguo)
#### Post date: [June 6, 2024, 2:17pm UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/6 "2024-06-06T14:17:35Z")
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!
---
### Author:  [@fqguo](https://discourse.julialang.org/u/fqguo)
#### Post date: [June 6, 2024, 2:24pm UTC](https://discourse.julialang.org/t/computing-fidelity-of-quantum-state/115264/7 "2024-06-06T14:24:15Z")
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!