Here, we derive the exact analytical expression for Ramsey fringes.
from sympy import Function, symbols, sqrt, Eq, Abs, exp, cos, sin, Matrix, dsolve, solve, S, lambdify
from sympy import I as 𝕚, pi as π
g = Function('g')
e = Function('e')
t, τ = symbols('t, tau', positive=True)
Δ = symbols('Delta', real=True)
η = symbols('eta', positive=True)
Ω = symbols('Omega', real=True)
g0, e0 = symbols('g_0, e_0')
Ĥ = Matrix([
[-Δ/2, η ],
[ η, Δ/2]
])
Ĥ
TDSE_system = [
g(t).diff(t) + 𝕚 * (Ĥ[0,0] * g(t) + Ĥ[0,1] * e(t)),
e(t).diff(t) + 𝕚 * (Ĥ[1,0]* g(t) + Ĥ[1,1] * e(t)),
]
sols_gen = dsolve(TDSE_system, [g(t), e(t)])
effective_rabi_freq = {
Ω: sqrt(Δ**2 + (2*η)**2)
}
Eq(Ω, effective_rabi_freq[Ω])
find_Ω = {
sqrt(-Δ**2 - 4 * η**2): 𝕚 * Ω,
sqrt(Δ**2 + 4 * η**2): Ω
}
sols_gen[0].subs(find_Ω)
sols_gen[1].subs(find_Ω)
We specialize the above general solution to an arbitrary initial state with complex amplidues $g_0$, $e_0$.
boundary_conditions = {
t: 0,
g(t): g0,
e(t) : e0,
}
find_4ηsq = {
Ω**2 - Δ**2: 4 * η**2
}
C1, C2 = symbols("C1, C2")
_integration_constants = solve(
[sol.subs(find_Ω).subs(boundary_conditions) for sol in sols_gen],
[C1, C2]
)
integration_constants = {
k: v.subs(find_4ηsq)
for (k, v) in _integration_constants.items()
}
Eq(C1, integration_constants[C1])
Eq(C2, integration_constants[C2])
def simplify_sol_g(sol_g):
rhs = (
sol_g
.rhs
.rewrite(sin)
.expand()
.collect(𝕚)
.subs(effective_rabi_freq)
.simplify()
.subs(find_Ω)
)
return Eq(sol_g.lhs, rhs)
sol_g = simplify_sol_g(sols_gen[0].subs(find_Ω).subs(integration_constants))
sol_g
sol_e = (
sols_gen[1]
.subs(find_Ω)
.subs(integration_constants)
.expand()
.rewrite(sin)
.expand()
)
sol_e
For example, when starting from the ground state:
sol_g.subs({g0: 1, e0: 0})
sol_e.subs({g0: 1, e0: 0})
For the population dynamics, we find:
def abs_sq(eq):
lhs = eq.lhs
rhs = eq.rhs
return Eq(Abs(lhs)**2, (rhs * rhs.conjugate()).expand())
abs_sq(sol_g.subs({g0: 1, e0: 0}))
abs_sq(sol_e.subs({g0: 1, e0: 0}))
T = π / (4*η)
For the ground state amplitude, we find:
sol_g_πhalf = sol_g.subs({g0: 1, e0: 0}).subs({t:T})
sol_g_πhalf
sol_g_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: 0, η: 15})
sol_g_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: 0, η: 15}).evalf()
sol_g_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: S(1)/5, η: 15}).evalf()
For the excited state amplitude, we find:
sol_e_πhalf = sol_e.subs({g0: 1, e0: 0}).subs({t:T})
sol_e_πhalf
sol_e_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: 0, η: 15})
sol_e_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: 0, η: 15}).evalf()
sol_e_πhalf.subs({g0: 1, e0: 0}).subs(effective_rabi_freq).subs({Δ: S(1)/5, η: 15}).evalf()
timeshift_free = {g(τ): g(τ + T), e(τ): e(τ + T)}
sol_g_free = (
sol_g
.expand()
.subs(effective_rabi_freq)
.subs({η: 0})
.subs({Δ: symbols("Delta", positive=True)})
.subs({t:τ})
.rewrite(exp)
.expand()
.subs({symbols("Delta", positive=True): Δ})
.subs(timeshift_free)
)
sol_g_free
sol_e_free = (
sol_e
.expand()
.subs(effective_rabi_freq)
.subs({η: 0})
.subs({Δ: symbols("Delta", positive=True)})
.subs({t:τ})
.rewrite(exp)
.expand()
.subs({symbols("Delta", positive=True): Δ})
.subs(timeshift_free)
)
sol_e_free
Ψ_πhalf_ideal = {g0: sqrt(2)/2, e0: -𝕚*sqrt(2)/2}
sol_g_free_ideal = sol_g_free.subs(Ψ_πhalf_ideal)
sol_g_free_ideal
abs_sq(sol_g_free_ideal)
sol_e_free_ideal = sol_e_free.subs(Ψ_πhalf_ideal)
sol_e_free_ideal
abs_sq(sol_e_free_ideal)
Ψ_πhalf_actual = {g0: sol_g_πhalf.rhs, e0: sol_e_πhalf.rhs}
sol_g_free_actual = sol_g_free.subs(Ψ_πhalf_actual).expand()
sol_g_free_actual
sol_e_free_actual = sol_e_free.subs(Ψ_πhalf_actual).expand()
sol_e_free_actual
Ψ_free_actual = {g0: sol_g_free_actual.rhs, e0: sol_e_free_actual.rhs}
timeshift_final = {g(T): g(2*T + τ), e(T): e(2*T + τ)}
sol_g_final = sol_g.subs(Ψ_free_actual).subs({t:T}).expand().subs(timeshift_final)
sol_g_final
pop_g_final = abs_sq(sol_g_final)
pop_g_final
sol_e_final = sol_e.subs(Ψ_free_actual).subs({t:T}).expand().subs(timeshift_final)
sol_e_final
pop_e_final = abs_sq(sol_e_final)
pop_e_final
We can find an "ideal" expression in the limit $\eta \gg \Delta$. That is, $\Omega \approx 2\eta$ and $\Delta/\eta \rightarrow 0$:
pop_g_final_ideal = Eq(
symbols('p_g'),
pop_g_final
.subs({Ω: 2*η})
.rewrite(sin)
.expand()
.subs({Δ/η: 0})
.rhs
)
pop_g_final_ideal
pop_e_final_ideal = Eq(
symbols('p_e'),
pop_e_final
.subs({Ω: 2*η})
.rewrite(sin)
.expand()
.subs({Δ/η: 0})
.rhs
)
pop_e_final_ideal
Note that his is simply $\cos(\Delta\cdot\tau/2)^2$:
assert (pop_e_final_ideal.rhs - cos(Δ*τ/2)**2).simplify() == S(0)
The "Ramsey fringes" result from plotting the excitated state popultion for varying detunging $\Delta$ and fixed pulse amplitude $\eta$ and time-of-flight $\tau$.
from matplotlib import pylab as plt
import numpy as np
pop_e_func = lambdify([Δ, η, τ], pop_e_final.subs(effective_rabi_freq).rhs)
pop_e_ideal_func = lambdify([Δ, η, τ], pop_e_final_ideal.subs(effective_rabi_freq).rhs)
def plot_ramsey_fringes(
Δ_min=-1,
Δ_max=1,
η=1,
τ=10.0,
N=101,
ax=None,
label=None,
figsize=(10, 6),
func=pop_e_func,
legend=None,
):
if ax is None:
fig, ax = plt.subplots(figsize=figsize)
Δ_vals = np.linspace(Δ_min, Δ_max, N)
p_vals = func(Δ_vals, eta=η, tau=τ).real
ax.plot(Δ_vals, p_vals, label=label)
ax.set_xlabel("Δ")
ax.set_ylabel("pₑ")
if legend is None:
legend = label is not None
if legend:
ax.legend()
return ax
ax = plot_ramsey_fringes(Δ_min=-20, Δ_max=20, N=1001, label="Ramsey")
#plot_ramsey_fringes(Δ_min=-20, Δ_max=20, N=1001, label="ideal", ax=ax, func=pop_e_ideal_func)
ax = plot_ramsey_fringes(Δ_min=-5, Δ_max=5, N=1001, η=1, label="Ramsey");
ax = plot_ramsey_fringes(Δ_min=-5, Δ_max=5, N=1001, η=1, label="ideal", ax=ax, func=pop_e_ideal_func);
ax.figure.suptitle("Actual vs ideal fringes for small η");
ax.figure.tight_layout()
ax = plot_ramsey_fringes(Δ_min=-5, Δ_max=5, N=1001, η=10, label="Ramsey");
ax = plot_ramsey_fringes(Δ_min=-5, Δ_max=5, N=1001, η=10, label="ideal", ax=ax, func=pop_e_ideal_func);
ax.figure.suptitle("Actual vs ideal fringes for large η");
ax.figure.tight_layout()
ax = plot_ramsey_fringes(Δ_min=0.9, Δ_max=1.1, N=1001, η=1, label="η=1")
plot_ramsey_fringes(Δ_min=0.9, Δ_max=1.1, N=1001, η=10, label="η=10", ax=ax)
plot_ramsey_fringes(Δ_min=0.9, Δ_max=1.1, N=1001, η=100, label="η=100", ax=ax)
plot_ramsey_fringes(Δ_min=0.9, Δ_max=1.1, N=1001, label="ideal", ax=ax, func=pop_e_ideal_func);
ax.figure.suptitle("Probing around Δ=1");
ax.figure.tight_layout()
ax = plot_ramsey_fringes(Δ_min=-1, Δ_max=1, N=1001, τ=10, label="τ=10");
plot_ramsey_fringes(Δ_min=-1, Δ_max=1, N=1001, τ=1, label="τ=1", ax=ax);
plot_ramsey_fringes(Δ_min=-1, Δ_max=1, N=1001, τ=20, label="τ=20", ax=ax);
ax.figure.suptitle("Varying the time of flight");
ax.figure.tight_layout()
Conclusions: