Accessing the Ion position after time evolution

[pkalaczynski]

Hi,

it is not clear to me how can I access the position of the ion after performing time evolution of the system, as described in https://examples.ionsim.org/

I guess it should also be done via expect(), but how should I construct (and use) the right ion position operator if the solution has a composite basis like:

sol[1].basis

[⁴⁰Ca ⊗ ⁴⁰Ca ⊗ VibrationalMode(ν=3.0e6, axis=x̂, N=10) ⊗ VibrationalMode(ν=3.0e6, axis=x̂, N=10) ⊗ VibrationalMode(ν=3.0e6, axis=x̂, N=10)]

Sorry if it's a noob question, but this is who I am in this context 🙈

[neil-glikin]

You can construct a position operator from the creation and annihilation operators of some chosen vibrational mode: $\hat{x} = \sqrt{\hbar/2m\omega}(\hat{a}+\hat{a}^\dagger)$. In IonSim you can use the functions create and destroy to define the ladder operators: a = destroy(v::VibrationalMode), adag = create(v::VibrationalMode), and the function frequency to get the frequency: ω = 2π*frequency(v::VibrationalMode).

So if m is the ion mass, something like x = √(ℏ/(2mω))*(a + adag) (given the above definitions) should give the position operator you're looking for.

[pkalaczynski]

Thanks! What about the ionposition() function? Can it be used as well or is it meant for something else?

By the way the result I get:

seems to be independent of the number of ions I put in the trap (or at least it seems so to me), i.e. I get the same with a single ion and two or three. I imagine this should not be like that, even if the result is for the center of mass of the whole set of ions, plus I'd really like to have individual positions, or even better probability densities, similar to the QuantumOptics.jl example in https://docs.qojulia.org/examples/particle-in-harmonic-trap/
Is this also doable?

What I did was basically:

# define the x position operator in the composite basis of sol:
m = mass(ca)
ω = 2π*ν
ion1=ions(chamber)[1]
ion2=ions(chamber)[2]
ion3=ions(chamber)[3]
mode_x=modes(chamber)[1]
mode_y=modes(chamber)[2]
mode_z=modes(chamber)[3]

x̂_op = sqrt(ħ/(2. *m*ω))*(IonSim.destroy(mode_x) + IonSim.create(mode_x))
x̂₍composite₎ = one(ion1) ⊗ one(ion2) ⊗ one(ion3)⊗  x̂_op ⊗ one(mode_y) ⊗ one(mode_z); # we need to adapt to the basis of sol

ŷ_op = sqrt(ħ/(2. *m*ω))*(IonSim.destroy(mode_y) + IonSim.create(mode_y))
ŷ₍composite₎ = one(ion1) ⊗ one(ion2) ⊗ one(ion3)⊗  one(mode_x) ⊗ ŷ_op ⊗ one(mode_z); # we need to adapt to the basis of sol

ẑ_op = sqrt(ħ/(2. *m*ω))*(IonSim.destroy(mode_z) + IonSim.create(mode_z))
ẑ₍composite₎ = one(ion1) ⊗ one(ion2) ⊗ one(ion3)⊗  one(mode_x) ⊗ one(mode_y) ⊗ ẑ_op; # we need to adapt to the basis of sol

x_vals = real.(expect(x̂₍composite₎, sol))
y_vals = real.(expect(ŷ₍composite₎, sol))
z_vals = real.(expect(ẑ₍composite₎, sol));

p=plot(T,x_vals*1e12,label=L"x")
plot!(T,y_vals*1e12,label=L"y")
plot!(T,z_vals*1e12,label=L"z")
xlabel!("time [μs]")
ylabel!("position [fm]")
savefig("plot.png")
@show p

[pkalaczynski]

Hmm, okay, I think I understood that ionposition() will return me the initial position for each ion (after I've put it into the trap). But still I an unsure if there's a proper way to use it after time evolving the system? Or do I just add the vibrational mode results to it to obtain the position at a certain time?

The 2nd thing which I find confusing: what dimension are these positions along? Do we just always assume a ### trap and then it's along the long axis? Looking at https://github.com/HaeffnerLab/IonSim.jl/blob/master/src/iontraps.jl#L112, the distances are basically just a function of the number of ions ...

[neil-glikin]

You are correct; there is not built-in position operator, so you would need to construct it yourself for each ion as the appropriate sum of contributions from all motional modes. ionposition is simply the equlibrium position of the ion.

To answer your second question, yes, in IonSim the ions are always along a single axis, the z-axis. The position is always along that axis.

[pkalaczynski]

Alright, thanks, this already brought some clarity :)