Single Qubit Photon Sensor with Quantum Pulses¶

  1. Quantum non-demolition detection of an itinerant microwave photon
  2. Input-Output Theory with Quantum Pulses
  3. Quantum interactions with pulses of radiation
  4. Interactions of quantum systems with pulses of quantized radiation: From a cascaded master equation to a traveling mode perspective
In [1]:
import matplotlib.pyplot as plt
import numpy as np

from mpl_toolkits.mplot3d import Axes3D
from qutip import *
from qutip_qip.operations import ry, hadamard_transform
from scipy.integrate import quad, trapezoid, cumulative_trapezoid
from scipy.interpolate import CubicSpline
from typing import Callable, Union
In [2]:
# I choose you Bulbasaur!
plt.style.use('../../styles/bulbasaur.mplstyle')

Definition of pulses¶

The input pulse is a Gaussian pulse $$ \begin{aligned} u(t) = \frac{1}{\sqrt[4]{\pi}\sqrt{\sigma_t}}e^{-\frac{(t-t_0)^2}{2\sigma_t^2}}e^{-i\omega_0(t-t_0)} \qquad , \qquad u(\omega) = \frac{1}{\sqrt[4]{\pi}\sqrt{\sigma_\omega}}e^{-\frac{(\omega-\omega_0)^2}{2\sigma_\omega^2}}e^{i\omega t_0} \qquad , \qquad \sigma_\omega \sigma_t = 1 \end{aligned} $$ while the output pulse is obtained from the reflection relation $$ \begin{aligned} r(\omega) = \frac{i(\omega - \omega_c) + \frac\gamma2}{i(\omega - \omega_c) - \frac\gamma2} \end{aligned} $$ The effective couplings are defined as $$ \begin{aligned} g_u(t) = \frac{u^*(t)}{\sqrt{1 - \int_0^t\mathrm{dt}' |u(t')|^2}} \qquad , \qquad g_v(t) = -\frac{v^*(t)}{\sqrt{\int_0^t\mathrm{dt}' |v(t')|^2}} \end{aligned} $$

In [3]:
def u_omega(
    time_0: float, omega_0: float, sigma_omega: float
) -> Callable[[float], complex]:
    prefactor = 1 / (np.pi**0.25 * np.sqrt(sigma_omega))

    def func(omega: float) -> complex:
        return (
            prefactor
            * np.exp(-((omega - omega_0) ** 2) / (2 * sigma_omega**2))
            * np.exp(1j * omega * time_0)
        )

    return func


def u_time(
    time_0: float, omega_0: float, sigma_omega: float
) -> Callable[[float], complex]:
    sigma_time = 1 / sigma_omega
    prefactor = 1 / (np.pi**0.25 * np.sqrt(sigma_time))

    def func(time: float) -> complex:
        return (
            prefactor
            * np.exp(-((time - time_0) ** 2) / (2 * sigma_time**2))
            * np.exp(-1j * omega_0 * (time - time_0))
        )

    return func


def u_integral(
    time_0: float, omega_0: float, sigma_omega: float, times: np.ndarray
) -> CubicSpline:
    u_values = u_time(time_0, omega_0, sigma_omega)(times)
    u_abs_2 = np.abs(u_values) ** 2
    u_abs_2_int = cumulative_trapezoid(u_abs_2, times, initial=0)
    return CubicSpline(times, u_abs_2_int)


def r_omega(omega_c: float, gamma: float) -> Callable[[float], complex]:
    def func(omega: float) -> complex:
        return (1j * (omega - omega_c) + gamma / 2) / (
            1j * (omega - omega_c) - gamma / 2
        )

    return func

def v_omega(
    time_0: float, omega_0: float, sigma_omega: float, omega_c: float, gamma: float
) -> Callable[[float], complex]:
    u_w = u_omega(time_0, omega_0, sigma_omega)
    r_w = r_omega(omega_c, gamma)
    return lambda omega: r_w(omega) * u_w(omega)


def get_v_time(v_w: Callable[[float], complex], times: np.ndarray) -> CubicSpline:
    v_t = np.array(
        [inverse_fourier_transform(v_w, t) for t in times], dtype=np.complex128
    )
    v_t /= np.sqrt(trapezoid(np.abs(v_t) ** 2, times))
    return CubicSpline(times, v_t)


def v_time(
    time_0: float,
    omega_0: float,
    sigma_omega: float,
    omega_c: float,
    gamma: float,
    times: np.ndarray,
) -> CubicSpline:
    v_w = v_omega(time_0, omega_0, sigma_omega, omega_c, gamma)
    return get_v_time(v_w, times)


def inverse_fourier_transform(f: Callable[[float], complex], t: float) -> complex:
    def f_with_factor(w: float) -> complex:
        return f(w) * np.exp(-1j * w * t) / np.sqrt(2 * np.pi)

    return (
        quad(lambda w: np.real(f_with_factor(w)), -np.inf, np.inf)[0]
        + 1j * quad(lambda w: np.imag(f_with_factor(w)), -np.inf, np.inf)[0]
    )


def v_integral(
    time_0: float,
    omega_0: float,
    sigma_omega: float,
    omega_c: float,
    gamma: float,
    times: np.ndarray,
) -> CubicSpline:
    v_values = v_time(time_0, omega_0, sigma_omega, omega_c, gamma, times)(times)
    v_abs_2 = np.abs(v_values) ** 2
    v_abs_2_int = cumulative_trapezoid(v_abs_2, times, initial=0)
    return CubicSpline(times, v_abs_2_int)


class Pulse:
    shape: str
    in_going: bool
    _u: Callable[[float], complex]
    _g: Callable[[float], complex]

    def __init__(self, shape: str, in_going: bool, *args: Union[float, np.ndarray]):
        self.shape = shape
        self.in_going = in_going
        if shape == "u":
            self._u = u_time(*args[:3])
            self._g = u_integral(*args[:4])
        elif shape == "v":
            self._u = v_time(*args[:6])
            self._g = v_integral(*args[:6])
        else:
            raise ValueError("Unsupported shape")

    def u(self, t: float) -> complex:
        return self._u(t)

    def g(self, t: float) -> complex:
        temp = np.conj(self.u(t))
        real = np.imag(temp) == 0
        if abs(temp) >= EPSILON:
            factor = np.sqrt(
                EPSILON + (1 - self._g(t)) if self.in_going else EPSILON + self._g(t)
            )
            temp = temp / factor if self.in_going else -temp / factor
        if not self.in_going and abs(temp) >= 1 / EPSILON:
            temp = 0
        return np.real(temp) if real else temp
In [17]:
# To avoid numeric precision
EPSILON = 1e-10

# Cavity dissipation rate
gamma = (2 * np.pi) * 0.00332 # GHz

# Qubit-cavity coupling
chi = gamma / 2

# Cavity
omega_c = 0 # GHz

# Pulse
time_0 = 0
omega_0 = 0
sigma_time = 1500 / (2 * np.sqrt(2 * np.log(2))) # ns
sigma_omega = 1/sigma_time

omegas = np.linspace(omega_0 - 5 * sigma_omega, omega_0 + 5 * sigma_omega, 1000)
times = np.linspace(time_0 - 5 * sigma_time, time_0 + 5 * sigma_time, 1000)
In [18]:
u_pulses = Pulse("u", True, time_0, omega_0, sigma_omega, times)
v_pulses = Pulse("v", False, time_0, omega_0, sigma_omega, omega_c, gamma, times)

u_w_values = u_omega(time_0, omega_0, sigma_omega)(omegas)
v_w_values = v_omega(time_0, omega_0, sigma_omega, omega_c, gamma)(omegas)

u_data = [u_pulses.u(t) for t in times]
v_data = [v_pulses.u(t) for t in times]
gu_data = [u_pulses.g(t) for t in times]
gv_data = [v_pulses.g(t) for t in times]
In [19]:
colors = plt.cm.managua(np.linspace(0, 1, 10))

fig, ax = plt.subplots(2, 1, figsize=(14, 8))

ax[0].fill(omegas, np.real(u_w_values), label=r"$\mathrm{Re}\,u(\omega)$", color=colors[1], alpha=0.7)
ax[0].fill(omegas, np.imag(u_w_values), label=r"$\mathrm{Im}\,u(\omega)$", color=colors[2], alpha=0.7)
ax[0].fill(omegas, np.real(v_w_values), label=r"$\mathrm{Re}\,v(\omega)$", color=colors[-2], alpha=0.7)
ax[0].fill(omegas, np.imag(v_w_values), label=r"$\mathrm{Im}\,v(\omega)$", color=colors[-3], alpha=0.7)
ax[0].set_xlim(ax[0].get_xlim()[0], ax[0].get_xlim()[1]*1.2)
ax[0].set_xlabel("Angular frequency")
ax[0].set_ylabel(r"Mode")
ax[0].grid()
ax[0].legend()

ax[1].fill(times, np.real(u_data), label=r"$\mathrm{Re}\,u(t)$", color=colors[1], alpha=0.7)
ax[1].fill(times, np.imag(u_data), label=r"$\mathrm{Im}\,u(t)$", color=colors[2], alpha=0.7)
ax[1].fill(times, np.real(v_data), label=r"$\mathrm{Re}\,v(t)$", color=colors[-2], alpha=0.7)
ax[1].fill(times, np.imag(v_data), label=r"$\mathrm{Im}\,v(t)$", color=colors[-3], alpha=0.7)
ax[1].set_xlim(ax[1].get_xlim()[0], ax[1].get_xlim()[1]*1.2)
ax[1].set_xlabel("Time")
ax[1].set_ylabel(r"Mode")
ax[1].grid()

ax2 = ax[1].twinx()
ax2.plot(times, np.abs(gu_data)**2, label=r"$|g_u(t)|^2$", color=colors[2],  linewidth=3)
ax2.plot(times, np.abs(gv_data)**2, label=r"$|g_v(t)|^2$", color=colors[-3],  linewidth=3, linestyle="dashed")
ax2.set_ylabel(r"Coupling", labelpad=2)

handles1, labels1 = ax[1].get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()

handles = handles1 + handles2
labels = labels1 + labels2

new_order = [0, 1, 4, 2, 3, 5]

handles = [handles[i] for i in new_order]
labels = [labels[i] for i in new_order]

ax2.legend(handles, labels, loc="upper right")

plt.tight_layout()
#plt.savefig("plots/pulse.pdf")
plt.show()
No description has been provided for this image

Replicating Kono's paper¶

In [20]:
N = 7
M = 7
C = 7
Q = 2

au = destroy(N)
av = destroy(M)
c = destroy(C)

Iu = qeye(N)
Iv = qeye(M)
Ic = qeye(C)
Iq = qeye(Q)

au_tot = tensor(au, Iv, Ic, Iq)
av_tot = tensor(Iu, av, Ic, Iq)
c_tot = tensor(Iu, Iv, c, Iq)
sigmaz_tot = tensor(Iu, Iv, Ic, sigmaz())

The density matrix of the total system evolves according to a master equation on the Lindblad form $$ \begin{aligned} \frac{\mathrm{d}\hat \rho}{\mathrm{d}t} = -i[\hat H, \hat \rho] + \sum_{i=0} \left(\hat L_i \hat \rho \hat L_i^\dagger - \frac12 \{\hat L_i^\dagger \hat L_i, \hat \rho\}\right) \end{aligned} $$ the Hamiltonian contains terms that represent coherent exchange of energy between the three components $$ \begin{aligned} \hat H(t) = \chi \hat c^\dagger \hat c \hat \sigma_z + \frac i2 \left(\sqrt{\gamma} g_u^*(t)\hat a_u^\dagger \hat c + \sqrt{\gamma^*} g_v(t) \hat c ^\dagger \hat a_v + g_u^*(t)g_v(t)\hat a_u^\dagger \hat a_v - \mathrm{h.c.}\right) \end{aligned} $$ The damping terms include a single Lindblad operator of the form representing the output loss from the last cavity $$ \begin{aligned} \hat L_0(t) = \sqrt{\gamma}\hat c + g_u(t) \hat a_u + g_v(t) \hat a_v \end{aligned} $$

In [21]:
h_0 = chi * c_tot.dag() * c_tot * sigmaz_tot
h_int = QobjEvo(
    [
        [
            0.5j * np.sqrt(gamma) * (au_tot.dag() * c_tot),
            lambda t, args: np.conjugate(u_pulses.g(t)),
        ],
        [
            0.5j * np.sqrt(np.conjugate(gamma)) * (c_tot.dag() * av_tot),
            lambda t, args: v_pulses.g(t),
        ],
        [
            0.5j * (au_tot.dag() * av_tot),
            lambda t, args: np.conjugate(u_pulses.g(t)) * v_pulses.g(t),
        ],
        [
            - 0.5j * np.sqrt(np.conjugate(gamma)) * (c_tot.dag() * au_tot),
            lambda t, args: u_pulses.g(t),
        ],
        [
            - 0.5j * np.sqrt(gamma) * (av_tot.dag() * c_tot),
            lambda t, args: np.conjugate(v_pulses.g(t)),
        ],
        [
            - 0.5j * (av_tot.dag() * au_tot),
            lambda t, args: u_pulses.g(t) * np.conjugate(v_pulses.g(t)),
        ]
    ]
)

h_tot = h_0 + h_int

L_0 = [
    np.sqrt(gamma) * c_tot
    + QobjEvo([[au_tot, lambda t, args: u_pulses.g(t)]])
    + QobjEvo([[av_tot, lambda t, args: v_pulses.g(t)]])
]

Fock state¶

In [22]:
# Input mode: Fock state
options = {"progress_bar": True}
psi_0 = tensor(basis(N, 1), basis(M, 0), basis(C, 0), (basis(2, 0) + basis(2, 1)).unit())
result = mesolve(h_tot, psi_0, tlist=times, c_ops=L_0, options=options)

10.0%. Run time:   3.79s. Est. time left: 00:00:00:34
20.0%. Run time:   9.15s. Est. time left: 00:00:00:36
30.0%. Run time:  14.61s. Est. time left: 00:00:00:34
40.0%. Run time:  19.95s. Est. time left: 00:00:00:29
50.1%. Run time:  25.81s. Est. time left: 00:00:00:25
60.1%. Run time:  32.63s. Est. time left: 00:00:00:21
70.1%. Run time:  39.36s. Est. time left: 00:00:00:16
80.1%. Run time:  45.53s. Est. time left: 00:00:00:11
90.1%. Run time:  51.56s. Est. time left: 00:00:00:05
100.0%. Run time:  58.27s. Est. time left: 00:00:00:00
Total run time:  58.27s
In [23]:
ry_tot = tensor(Iu, Iv, Ic, ry(-np.pi / 2))
results_rotated = [ry_tot * state * ry_tot.dag() for state in result.states]
results_rotated_protocol = [state for state in result.states]
results_rotated_protocol[-1] = ry_tot * result.states[-1] * ry_tot.dag()
In [24]:
au_population = expect(au_tot.dag() * au_tot, results_rotated_protocol)
av_population = expect(av_tot.dag() * av_tot, results_rotated_protocol)
c_population = expect(c_tot.dag() * c_tot, results_rotated_protocol)
x_population = expect(tensor(Iu, Iv, Ic, sigmax()), results_rotated_protocol)
y_population = expect(tensor(Iu, Iv, Ic, sigmay()), results_rotated_protocol)
z_population = expect(tensor(Iu, Iv, Ic, sigmaz()), results_rotated_protocol)

I_out = []

for i in range(len(times)):
    L_0_t = L_0[0](times[i])
    expectation_value = expect(L_0_t.dag() * L_0_t, results_rotated_protocol[i])
    I_out.append(expectation_value)

I_out_int = trapezoid(I_out, times)
In [25]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, au_population, label=r"$\langle \hat a_u^\dagger \hat a_u\rangle$", color=colors[1], linewidth=2)
ax.plot(times, av_population, label=r"$\langle \hat a_v^\dagger \hat a_v\rangle$", color=colors[-2], linewidth=2)
ax.plot(times, c_population, label=r"$\langle \hat c^\dagger \hat c\rangle$", color="black", linewidth=2, linestyle="dashed")
ax.set_xlabel("Time")
ax.set_ylabel("Population")
ax.legend()
ax.grid()
#plt.savefig("plots/fock/population.pdf")
plt.show()
No description has been provided for this image
In [26]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, x_population, label=r"$\langle \hat \sigma_x \rangle$", color=colors[1], linewidth=2, linestyle="dashed")
ax.plot(times, y_population, label=r"$\langle \hat \sigma_y \rangle$", color=colors[4], linewidth=2, linestyle="dashed")
ax.plot(times, z_population, label=r"$\langle \hat \sigma_z \rangle$", color=colors[-2], linewidth=2, linestyle="dashed")
ax.set_xlabel("Time")
ax.set_ylabel("Qubit components")
ax.legend()
ax.grid()
#plt.savefig("plots/fock/population.pdf")
plt.show()
No description has been provided for this image
In [27]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, I_out, label=r"$I_\mathrm{out}(t) = \langle \hat L_0^\dagger(t) \hat L_0(t) \rangle$", color=colors[5], linewidth=2)
ax.scatter(0.00002, 0, label=rf"$\int_{{t_1}}^{{t_2}}\,\mathrm{{d}}t\, I_\mathrm{{out}}(t) =  {round(I_out_int, 3)}$", alpha=0)
ax.set_ylabel("Output flux")
ax.set_xlabel("Time")
ax.grid()
ax.legend()
#plt.savefig("plots/fock/output_flux.pdf")
plt.show()
No description has been provided for this image
In [28]:
indices = list(range(0, len(times), 15))

if 0 not in indices:
    indices.insert(0, 0)
if (len(times) - 1) not in indices:
    indices.append(len(times) - 1)
indices = sorted(set(indices))

colors_bs = plt.cm.managua(np.linspace(0, 1, len(times)))
selected_colors = [colors_bs[i] for i in indices]

b = Bloch()
b.point_color = selected_colors
b.sphere_color = "white"
b.frame_color = 'black'
b.frame_width = 1
b.point_size[0] = 30
b.add_points([
    [x_population[i] for i in indices],
    [y_population[i] for i in indices],
    [z_population[i] for i in indices]
])

fig = plt.figure(figsize=(6, 7))
ax = fig.add_subplot(111, projection='3d')
b.fig = fig
b.axes = ax
b.render()
ax.view_init(elev=30, azim=-65)

cbar_ax = fig.add_axes([0.2, 0.15, 0.6, 0.03])
norm = plt.Normalize(vmin=times[0], vmax=times[-1])
sm = plt.cm.ScalarMappable(cmap='managua', norm=norm)
sm.set_array([])
cbar = fig.colorbar(sm, cax=cbar_ax, orientation='horizontal')
cbar.set_label('Time')

plt.show()
No description has been provided for this image
In [29]:
# Projectors
P0 = tensor(qeye(N), qeye(C), qeye(M), basis(2, 0) * basis(2, 0).dag())
P1 = tensor(qeye(N), qeye(C), qeye(M), basis(2, 1) * basis(2, 1).dag())

# Density matrices
rhou = ptrace(results_rotated[0], 0)
rhov = ptrace(results_rotated[-1], 1)
rhov_0 = ptrace((P0 * results_rotated[-1] * P0) / (P0 * results_rotated[-1] * P0).tr(), 1)
rhov_1 = ptrace((P1 * results_rotated[-1] * P1) / (P1 * results_rotated[-1] * P1).tr(), 1)

values = np.linspace(-3, 3, int(1e2))

w_rhou = wigner(rhou, values, values)
w_rhov = wigner(rhov, values, values)
w_rhov_0 = wigner(rhov_0, values, values)
w_rhov_1 = wigner(rhov_1, values, values)

fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)

c1 = ax[0,0].pcolor(values, values, w_rhou, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[0,0].set_title("Input mode", pad=8)
ax[0,0].set_ylabel(r"$Q$")
ax[0,0].set_xlim(-3, 3)
ax[0,0].set_ylim(-3, 3)
ax[0,0].set_aspect('equal')
ax[0,0].grid(True)

c2 = ax[0,1].pcolor(values, values, w_rhov, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[0,1].text(2.2, 2.3, r"$\boldsymbol{|?\rangle}$", color="w", fontsize=12)
ax[0,1].set_title("Output mode", pad=8)
ax[0,1].set_xlim(-3, 3)
ax[0,1].set_ylim(-3, 3)
ax[0,1].set_aspect('equal')
ax[0,1].grid(True)

c3 = ax[1,0].pcolor(values, values, w_rhov_0, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[1,0].text(2.2, 2.3, r"$\boldsymbol{|g\rangle}$", color="w", fontsize=12)
ax[1,0].text(-2.6, 2.25, rf"$\boldsymbol{{\mathcal{{F}}={round(fidelity(rhov_0, basis(M, 0)),2)}}}$", color="w", fontsize=12)
ax[1,0].set_title("Output mode", pad=8)
ax[1,0].set_xlabel(r"$I$")
ax[1,0].set_ylabel(r"$Q$")
ax[1,0].set_xlim(-3, 3)
ax[1,0].set_ylim(-3, 3)
ax[1,0].set_aspect('equal')
ax[1,0].grid(True)

c4 = ax[1,1].pcolor(values, values, w_rhov_1, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[1,1].text(2.2, 2.3, r"$\boldsymbol{|e\rangle}$", color="w", fontsize=12)
ax[1,1].text(-2.6, 2.25, rf"$\boldsymbol{{\mathcal{{F}}={round(fidelity(rhov_1, basis(M, 1)),2)}}}$", color="w", fontsize=12)
ax[1,1].set_title("Output mode", pad=8)
ax[1,1].set_xlabel(r"$I$")
ax[1,1].set_xlim(-3, 3)
ax[1,1].set_ylim(-3, 3)
ax[1,1].set_aspect('equal')
ax[1,1].grid(True)

cbar = fig.colorbar(c1, ax=ax, orientation='vertical')
cbar.ax.set_position([0.67, 0.12, 0.15, 0.752])
cbar.set_ticks(ticks=[-1/np.pi, 0, 1/np.pi], labels=[r'$-1/\pi$', r'$0$', r'$1/\pi$'])
cbar.ax.set_aspect('auto')
plt.subplots_adjust(left=0.1, right=0.65)
#plt.savefig("plots/fock/wigner.pdf")
plt.show()
No description has been provided for this image
In [31]:
fidelities = {
    "0": [
        fidelity(ptrace((P0 * state * P0) / (P0 * state * P0).tr(), 1), basis(M, 0))
        if (P0 * state * P0).tr() != 0 else 0
        for state in results_rotated
    ],
    "1": [
        fidelity(ptrace((P1 * state * P1) / (P1 * state * P1).tr(), 1), basis(M, 1))
        if (P1 * state * P1).tr() != 0 else 0
        for state in results_rotated
    ]
}
In [37]:
fig, ax = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 1], 'hspace': 0}, sharex=True)

ax[0].plot(times, fidelities["0"], label=r"$\mathcal{F}(|0\rangle\,|\,|g\rangle)$", color=colors[1], linewidth=2, zorder=1)
ax[0].axhline(np.min(fidelities["0"]), label=rf"$\mathcal{{F}}_\mathrm{{min}} = {round(np.min(fidelities['0']),2)}$", color=colors[2], linewidth=2, linestyle="dotted", zorder=0)
ax[0].set_ylabel("Fidelity")
ax[0].legend(loc=7)

ax[1].plot(times, fidelities["1"], label=r"$\mathcal{F}(|1\rangle\,|\,|e\rangle)$", color=colors[-2], linewidth=2, zorder=1)
ax[1].axhline(np.max(fidelities["1"][-500:]), label=rf"$\mathcal{{F}}_\mathrm{{max}} = {round(np.max(fidelities['1'][-500:]),2)}$", color=colors[-3], linewidth=2, linestyle="dotted", zorder=0)
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Fidelity")
ax[1].legend(loc=7)

plt.tight_layout()
#plt.savefig("plots/fock/fidelity.pdf")
plt.show()
No description has been provided for this image
In [33]:
negativities = {
    "oc": [negativity(ptrace(state, [1,2]),0) for state in results_rotated],
    "oq": [negativity(ptrace(state, [1,3]),0) for state in results_rotated]}
In [34]:
fig, ax = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 1], 'hspace': 0}, sharex=True)

ax[0].plot(times, negativities["oc"], label=r"$\mathcal{N}^{\mathrm{\,output\,\,field}}_{\mathrm{cavity}}$", color=colors[1], linewidth=2, zorder=1)
ax[0].axhline(np.max(negativities["oc"]), label=rf"$\mathcal{{N}}_\mathrm{{max}} = {round(np.max(negativities['oc']),2)}$", color=colors[2], linewidth=2, linestyle="dotted", zorder=0)
ax[0].set_ylabel("Negativity")
ax[0].legend(loc=6)

ax[1].plot(times, negativities["oq"], label=r"$\mathcal{N}^{\mathrm{\,output\,\,field}}_{\mathrm{qubit}}$", color=colors[-2], linewidth=2, zorder=1)
ax[1].axhline(np.max(negativities["oq"]), label=rf"$\mathcal{{N}}_\mathrm{{max}} = {round(np.max(negativities['oq']),2)}$", color=colors[-3], linewidth=2, linestyle="dotted", zorder=0)
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Negativity")
ax[1].legend(loc=6)

plt.tight_layout()
#plt.savefig("plots/fock/negativity.pdf")
plt.show()
No description has been provided for this image
In [38]:
qubit_state = [ptrace(state, 3) for state in results_rotated]

prob_g = expect(qubit_state, basis(2,0))
prob_e = expect(qubit_state, basis(2,1))
In [39]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, prob_g, label=r"$P(|g\rangle)$", color=colors[-2], linewidth=2, zorder=1)
ax.plot(times, prob_e, label=r"$P(|e\rangle)$", color=colors[1], linewidth=2, zorder=2)
ax.set_xlabel("Time")
ax.set_ylabel("Probability")
ax.legend()
ax.grid()
#plt.savefig("plots/fock/probability.pdf")
plt.show()
No description has been provided for this image

Coherent state¶

In [40]:
# Input mode: Coherent state
options = {"progress_bar": True}
psi_0 = tensor(coherent(N, np.sqrt(0.165)), basis(M, 0), basis(C, 0), (basis(2, 0) + basis(2, 1)).unit())
result = mesolve(h_tot, psi_0, tlist=times, c_ops=L_0, options=options)

10.0%. Run time:   3.96s. Est. time left: 00:00:00:35
20.0%. Run time:  11.26s. Est. time left: 00:00:00:45
30.0%. Run time:  35.54s. Est. time left: 00:00:01:22
40.0%. Run time:  61.85s. Est. time left: 00:00:01:32
50.1%. Run time:  90.01s. Est. time left: 00:00:01:29
60.1%. Run time: 116.87s. Est. time left: 00:00:01:17
70.1%. Run time: 144.07s. Est. time left: 00:00:01:01
80.1%. Run time: 170.58s. Est. time left: 00:00:00:42
90.1%. Run time: 196.54s. Est. time left: 00:00:00:21
100.0%. Run time: 221.75s. Est. time left: 00:00:00:00
Total run time: 221.76s
In [41]:
ry_tot = tensor(Iu, Iv, Ic, ry(-np.pi / 2))
results_rotated = [ry_tot * state * ry_tot.dag() for state in result.states]
results_rotated_protocol = [state for state in result.states]
results_rotated_protocol[-1] = ry_tot * result.states[-1] * ry_tot.dag()
In [42]:
au_population = expect(au_tot.dag() * au_tot, results_rotated_protocol)
av_population = expect(av_tot.dag() * av_tot, results_rotated_protocol)
c_population = expect(c_tot.dag() * c_tot, results_rotated_protocol)
x_population = expect(tensor(Iu, Iv, Ic, sigmax()), results_rotated_protocol)
y_population = expect(tensor(Iu, Iv, Ic, sigmay()), results_rotated_protocol)
z_population = expect(tensor(Iu, Iv, Ic, sigmaz()), results_rotated_protocol)

I_out = []

for i in range(len(times)):
    L_0_t = L_0[0](times[i])
    expectation_value = expect(L_0_t.dag() * L_0_t, results_rotated_protocol[i])
    I_out.append(expectation_value)

I_out_int = trapezoid(I_out, times)
In [43]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, au_population, label=r"$\langle \hat a_u^\dagger \hat a_u\rangle$", color=colors[1], linewidth=2)
ax.plot(times, av_population, label=r"$\langle \hat a_v^\dagger \hat a_v\rangle$", color=colors[-2], linewidth=2)
ax.plot(times, c_population, label=r"$\langle \hat c^\dagger \hat c\rangle$", color="black", linewidth=2, linestyle="dashed")
ax.set_xlabel("Time")
ax.set_ylabel("Population")
ax.legend()
ax.grid()
#plt.savefig("plots/coherent/population_c.pdf")
plt.show()
No description has been provided for this image
In [44]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, I_out, label=r"$I_\mathrm{out}(t) = \langle \hat L_0^\dagger(t) \hat L_0(t) \rangle$", color=colors[5], linewidth=2)
ax.scatter(0.00002, 0, label=rf"$\int_{{t_1}}^{{t_2}}\,\mathrm{{d}}t\, I_\mathrm{{out}}(t) =  {round(I_out_int, 3)}$", alpha=0)
ax.set_ylabel("Output flux")
ax.set_xlabel("Time")
ax.grid()
ax.legend()
#plt.savefig("plots/coherent/output_flux.pdf")
plt.show()
No description has been provided for this image
In [45]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, x_population, label=r"$\langle \hat \sigma_x \rangle$", color=colors[1], linewidth=2, linestyle="dashed")
ax.plot(times, y_population, label=r"$\langle \hat \sigma_y \rangle$", color=colors[4], linewidth=2, linestyle="dashed")
ax.plot(times, z_population, label=r"$\langle \hat \sigma_z \rangle$", color=colors[-2], linewidth=2, linestyle="dashed")
ax.set_xlabel("Time")
ax.set_ylabel("Qubit components")
ax.legend()
ax.grid()
#plt.savefig("plots/coherent/population_q.pdf")
plt.show()
No description has been provided for this image
In [46]:
indices = list(range(0, len(times), 15))

if 0 not in indices:
    indices.insert(0, 0)
if (len(times) - 1) not in indices:
    indices.append(len(times) - 1)
indices = sorted(set(indices))

colors_bs = plt.cm.managua(np.linspace(0, 1, len(times)))
selected_colors = [colors_bs[i] for i in indices]

b = Bloch()
b.point_color = selected_colors
b.sphere_color = "white"
b.frame_color = 'black'
b.frame_width = 1
b.point_size[0] = 30
b.add_points([
    [x_population[i] for i in indices],
    [y_population[i] for i in indices],
    [z_population[i] for i in indices]
])

fig = plt.figure(figsize=(6, 7))
ax = fig.add_subplot(111, projection='3d')
b.fig = fig
b.axes = ax
b.render()
ax.view_init(elev=30, azim=-65)

cbar_ax = fig.add_axes([0.2, 0.15, 0.6, 0.03])
norm = plt.Normalize(vmin=times[0], vmax=times[-1])
sm = plt.cm.ScalarMappable(cmap='managua', norm=norm)
sm.set_array([])
cbar = fig.colorbar(sm, cax=cbar_ax, orientation='horizontal')
cbar.set_label('Time')

plt.show()
No description has been provided for this image
In [47]:
# Projectors
P0 = tensor(qeye(N), qeye(M), qeye(C), basis(2, 0) * basis(2, 0).dag())
P1 = tensor(qeye(N), qeye(M), qeye(C), basis(2, 1) * basis(2, 1).dag())

# Density matrices
rhou = ptrace(results_rotated[0], 0)
rhov = ptrace(results_rotated[-1], 1)
rhov_0 = ptrace((P0 * results_rotated[-1] * P0) / (P0 * results_rotated[-1] * P0).tr(), 1)
rhov_1 = ptrace((P1 * results_rotated[-1] * P1) / (P1 * results_rotated[-1] * P1).tr(), 1)

values = np.linspace(-3, 3, int(1e2))

w_rhou = wigner(rhou, values, values)
w_rhov = wigner(rhov, values, values)
w_rhov_0 = wigner(rhov_0, values, values)
w_rhov_1 = wigner(rhov_1, values, values)

fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)

c1 = ax[0,0].pcolor(values, values, w_rhou, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[0,0].set_title("Input mode", pad=8)
ax[0,0].set_ylabel(r"$Q$")
ax[0,0].set_xlim(-3, 3)
ax[0,0].set_ylim(-3, 3)
ax[0,0].set_aspect('equal')
ax[0,0].grid(True)

c2 = ax[0,1].pcolor(values, values, w_rhov, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[0,1].text(2.2, 2.3, r"$\boldsymbol{|?\rangle}$", color="w", fontsize=12)
ax[0,1].set_title("Output mode", pad=8)
ax[0,1].set_xlim(-3, 3)
ax[0,1].set_ylim(-3, 3)
ax[0,1].set_aspect('equal')
ax[0,1].grid(True)

c3 = ax[1,0].pcolor(values, values, w_rhov_0, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[1,0].text(2.2, 2.3, r"$\boldsymbol{|g\rangle}$", color="w", fontsize=12)
ax[1,0].text(-2.6, 2.25, rf"$\boldsymbol{{\mathcal{{F}}={round(fidelity(rhov_0, basis(M, 0)),2)}}}$", color="w", fontsize=12)
ax[1,0].set_title("Output mode", pad=8)
ax[1,0].set_xlabel(r"$I$")
ax[1,0].set_ylabel(r"$Q$")
ax[1,0].set_xlim(-3, 3)
ax[1,0].set_ylim(-3, 3)
ax[1,0].set_aspect('equal')
ax[1,0].grid(True)

c4 = ax[1,1].pcolor(values, values, w_rhov_1, cmap="managua", vmin=-1/np.pi, vmax=1/np.pi)
ax[1,1].text(2.2, 2.3, r"$\boldsymbol{|e\rangle}$", color="w", fontsize=12)
ax[1,1].text(-2.6, 2.25, rf"$\boldsymbol{{\mathcal{{F}}={round(fidelity(rhov_1, basis(M, 1)),2)}}}$", color="w", fontsize=12)
ax[1,1].set_title("Output mode", pad=8)
ax[1,1].set_xlabel(r"$I$")
ax[1,1].set_xlim(-3, 3)
ax[1,1].set_ylim(-3, 3)
ax[1,1].set_aspect('equal')
ax[1,1].grid(True)

cbar = fig.colorbar(c1, ax=ax, orientation='vertical')
cbar.ax.set_position([0.67, 0.12, 0.15, 0.752])
cbar.set_ticks(ticks=[-1/np.pi, 0, 1/np.pi], labels=[r'$-1/\pi$', r'$0$', r'$1/\pi$'])
cbar.ax.set_aspect('auto')
plt.subplots_adjust(left=0.1, right=0.65)
#plt.savefig("plots/coherent/wigner.pdf")
plt.show()
No description has been provided for this image
In [48]:
fidelities = {
    "0":[fidelity(ptrace((P0 * state * P0) / (P0 * state * P0).tr(), 1), basis(M, 0)) for state in results_rotated],
    "1":[fidelity(ptrace((P1 * state * P1) / (P1 * state * P1).tr(), 1), basis(M, 1)) for state in results_rotated]
}
In [49]:
fig, ax = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 1], 'hspace': 0}, sharex=True)

ax[0].plot(times, fidelities["0"], label=r"$\mathcal{F}(|0\rangle\,|\,|g\rangle)$", color=colors[1], linewidth=2, zorder=1)
ax[0].axhline(np.min(fidelities["0"]), label=rf"$\mathcal{{F}}_\mathrm{{min}} = {round(np.min(fidelities['0']),2)}$", color=colors[2], linewidth=2, linestyle="dotted", zorder=0)
ax[0].set_ylabel("Fidelity")
ax[0].legend(loc=7)

ax[1].plot(times, fidelities["1"], label=r"$\mathcal{F}(|1\rangle\,|\,|e\rangle)$", color=colors[-2], linewidth=2, zorder=1)
ax[1].axhline(np.max(fidelities["1"]), label=rf"$\mathcal{{F}}_\mathrm{{max}} = {round(np.max(fidelities['1']),2)}$", color=colors[-3], linewidth=2, linestyle="dotted", zorder=0)
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Fidelity")
ax[1].legend(loc=7)

plt.tight_layout()
#plt.savefig("plots/coherent/fidelity.pdf")
plt.show()
No description has been provided for this image
In [50]:
negativities = {
    "oc": [negativity(ptrace(state, [1,2]),0) for state in results_rotated],
    "oq": [negativity(ptrace(state, [1,3]),0) for state in results_rotated]}
In [51]:
fig, ax = plt.subplots(2, 1, gridspec_kw={'height_ratios': [1, 1], 'hspace': 0}, sharex=True)

ax[0].plot(times, negativities["oc"], label=r"$\mathcal{N}^{\mathrm{\,output\,\,field}}_{\mathrm{cavity}}$", color=colors[1], linewidth=2, zorder=1)
ax[0].axhline(np.max(negativities["oc"]), label=rf"$\mathcal{{N}}_\mathrm{{max}} = {round(np.max(negativities['oc'])*1e5,2)} \times 10^{{-5}}$", color=colors[2], linewidth=2, linestyle="dotted", zorder=0)
ax[0].set_ylabel("Negativity")
ax[0].legend()

ax[1].plot(times, negativities["oq"], label=r"$\mathcal{N}^{\mathrm{\,output\,\,field}}_{\mathrm{qubit}}$", color=colors[-2], linewidth=2, zorder=1)
ax[1].axhline(np.max(negativities["oq"]), label=rf"$\mathcal{{N}}_\mathrm{{max}} = {round(np.max(negativities['oq']),2)}$", color=colors[-3], linewidth=2, linestyle="dotted", zorder=0)
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Negativity")
ax[1].legend(loc=6)

plt.tight_layout()
#plt.savefig("plots/coherent/negativity.pdf")
plt.show()
No description has been provided for this image
In [52]:
qubit_state = [ptrace(state, 3) for state in results_rotated]

alpha = np.sqrt(0.165)
prob_1 = alpha**2 * np.exp(-(alpha**2))
prob_e = expect(qubit_state, basis(2,1)) / prob_1
prob_g = 1 - prob_e
In [53]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times, prob_g, label=r"$P(|g\rangle\,|\,|1\rangle)$", color=colors[-2], linewidth=2, zorder=1)
ax.plot(times, prob_e, label=r"$P(|e\rangle\,|\,|1\rangle)$", color=colors[1], linewidth=2, zorder=2)
ax.set_xlabel("Time")
ax.set_ylabel("Conditional probability")
plt.legend()
plt.grid()
#plt.savefig("plots/coherent/conditional_prob.pdf")
plt.show()
No description has been provided for this image
In [56]:
rho_oq = ptrace(results_rotated[-1], [1, 3])
u_gate_tot = tensor(Iv, -ry(np.pi / 2))
rho_pm = u_gate_tot * rho_oq * u_gate_tot.dag()

rho_real = np.real(rho_pm.full())[:6,:6]
rho_imag = np.imag(rho_pm.full())[:6,:6]


labels = [r"$|0\,+\rangle$", r"$|0\,-\rangle$", r"$|1\,+\rangle$", r"$|1\,-\rangle$", r"$|2\,+\rangle$", r"$|2\,-\rangle$",]

fig = plt.figure(figsize=(16, 8))
ax = fig.add_subplot(111, projection='3d')
x, y = np.meshgrid(range(rho_real.shape[0]), range(rho_real.shape[1]))
ax.bar3d(x.flatten(), y.flatten(), np.zeros_like(rho_real.flatten()), 0.6, 0.6, rho_real.flatten(), color=colors[-2], alpha=0.7, edgecolor="black")
ax.set_xticks(range(len(labels)))
ax.set_yticks(range(len(labels)))

ax.set_xticklabels(labels, fontsize=14)
ax.set_yticklabels(labels, fontsize=14)
ax.tick_params(axis='z', labelsize=14)

ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r"$\mathrm{Re}\,\rho_{oq}$", rotation=0, labelpad=-40)
ax.set_zlim(0, 1)
ax.view_init(elev=30, azim=-45)

#plt.savefig("plots/coherent/qst_real.pdf")
plt.show()


fig = plt.figure(figsize=(16, 8))
ax = fig.add_subplot(111, projection='3d')
x, y = np.meshgrid(range(rho_imag.shape[0]), range(rho_imag.shape[1]))
ax.bar3d(x.flatten(), y.flatten(), np.zeros_like(rho_imag.flatten()), 0.6, 0.6, rho_imag.flatten(), color=colors[0], alpha=0.7, edgecolor="black")
ax.set_xticks(range(len(labels)))
ax.set_yticks(range(len(labels)))

ax.set_xticklabels(labels, fontsize=14)
ax.set_yticklabels(labels, fontsize=14)
ax.tick_params(axis='z', labelsize=14)

ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r"$\mathrm{Im}\,\rho_{oq}$", rotation=0, labelpad=-40)
ax.set_zlim(0, 1)
ax.view_init(elev=30, azim=-45)

#plt.savefig("plots/coherent/qst_imag.pdf")
plt.show()
No description has been provided for this image
No description has been provided for this image
In [57]:
import matplotlib

# Ottieni la colormap correttamente
cmap = matplotlib.colormaps.get_cmap('managua')

rhou_bar = np.real(np.diag(rhou.full()))
rhov_bar = np.real(np.diag(rhov.full()))
rhov_0_bar = np.real(np.diag(rhov_0.full()))
rhov_1_bar = np.real(np.diag(rhov_1.full()))
N = np.arange(len(rhou_bar))

fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)

x = np.arange(0, np.shape(rhou[:3,:3])[0])
y = np.arange(0, np.shape(rhou[:3,:3])[0])
z_real = np.abs(rhou[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

ax[0,0].bar(N, rhou_bar, color=cmap(np.diag(z_real)), edgecolor="black")
ax[0,0].set_title("Input mode", pad=8)
ax[0,0].set_ylabel("Occupation probability")
ax[0,0].grid(True)

x = np.arange(0, np.shape(rhov[:3,:3])[0])
y = np.arange(0, np.shape(rhov[:3,:3])[0])
z_real = np.abs(rhov[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

ax[0,1].bar(N, rhov_bar, color=cmap(np.diag(z_real)), edgecolor="black")
ax[0,1].set_title("Output mode", pad=8)
ax[0,1].text(5.8, 0.18, r"$\boldsymbol{|?\rangle}$", color="black", fontsize=12)
ax[0,1].grid(True)

x = np.arange(0, np.shape(rhov_0[:3,:3])[0])
y = np.arange(0, np.shape(rhov_0[:3,:3])[0])
z_real = np.abs(rhov_0[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

ax[1,0].bar(N, rhov_0_bar, color=cmap(np.diag(z_real)), edgecolor="black")
ax[1,0].text(5.8, 0.18, r"$\boldsymbol{|g\rangle}$", color="black", fontsize=12)
ax[1,0].set_title("Output mode", pad=8)
ax[1,0].set_xlabel("Fock number")
ax[1,0].set_ylabel("Occupation probability")
ax[1,0].set_xticks(N)
ax[1,0].grid(True)

x = np.arange(0, np.shape(rhov_1[:3,:3])[0])
y = np.arange(0, np.shape(rhov_1[:3,:3])[0])
z_real = np.abs(rhov_1[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

ax[1,1].bar(N, rhov_1_bar, color=cmap(np.diag(z_real)), edgecolor="black")
ax[1,1].text(5.8, 0.18, r"$\boldsymbol{|e\rangle}$", color="black", fontsize=12)
ax[1,1].set_title("Output mode", pad=8)
ax[1,1].set_xlabel("Fock number")
ax[1,1].set_xticks(N)
ax[1,1].grid(True)

from mpl_toolkits.axes_grid1.inset_locator import inset_axes

x = np.arange(0, np.shape(rhov_0[:3,:3])[0])
y = np.arange(0, np.shape(rhov_0[:3,:3])[0])
z_real = np.abs(rhov_0[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

axins = inset_axes(ax[1, 0], width="100%", height="100%", bbox_to_anchor=(253, 170, 100, 100))

c1 = axins.pcolor(x, y, z_real, cmap="managua", vmin=0, vmax=1)
axins.set_xticks(N[:3])
axins.set_yticks(N[:3])
axins.tick_params(labelsize=8)
axins.set_aspect('equal')
cbaxes = inset_axes(axins, width="5%", height="100%", loc='right',
                    bbox_to_anchor=(0.15, 0., 1, 1),
                    bbox_transform=axins.transAxes,
                    borderpad=0)

cbar = plt.colorbar(c1, cax=cbaxes)
cbar.ax.tick_params(labelsize=6)


x = np.arange(0, np.shape(rhov_1[:3,:3])[0])
y = np.arange(0, np.shape(rhov_1[:3,:3])[0])
z_real = np.abs(rhov_1[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

axins2 = inset_axes(ax[1, 1], width="100%", height="100%", bbox_to_anchor=(590, 170, 100, 100))

c2 = axins2.pcolor(x, y, z_real, cmap="managua", vmin=0, vmax=1)
axins2.set_xticks(N[:3])
axins2.set_yticks(N[:3])
axins2.tick_params(labelsize=8)
axins2.set_aspect('equal')
cbaxes2 = inset_axes(axins2, width="5%", height="100%", loc='right',
                    bbox_to_anchor=(3.53, 0., 1, 1),
                    bbox_transform=axins.transAxes,
                    borderpad=0)

cbar2 = plt.colorbar(c2, cax=cbaxes2)
cbar2.ax.tick_params(labelsize=6)

x = np.arange(0, np.shape(rhov[:3,:3])[0])
y = np.arange(0, np.shape(rhov[:3,:3])[0])
z_real = np.abs(rhov[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

axins3 = inset_axes(ax[0, 1], width="100%", height="100%", bbox_to_anchor=(590, 422, 100, 100))

c3 = axins3.pcolor(x, y, z_real, cmap="managua", vmin=0, vmax=1)
axins3.set_xticks(N[:3])
axins3.set_yticks(N[:3])
axins3.tick_params(labelsize=8)
axins3.set_aspect('equal')
cbaxes3 = inset_axes(axins3, width="5%", height="100%", loc='right',
                    bbox_to_anchor=(3.53, 2.52, 1, 1),
                    bbox_transform=axins.transAxes,
                    borderpad=0)

cbar3 = plt.colorbar(c3, cax=cbaxes3)
cbar3.ax.tick_params(labelsize=6)

x = np.arange(0, np.shape(rhou[:3,:3])[0])
y = np.arange(0, np.shape(rhou[:3,:3])[0])
z_real = np.abs(rhou[:3,:3].copy().real)
z_real[z_real < 1e-4] = 0

axins4 = inset_axes(ax[1, 1], width="100%", height="100%", bbox_to_anchor=(253, 422, 100, 100))

c4 = axins4.pcolor(x, y, z_real, cmap="managua", vmin=0, vmax=1)
axins4.set_xticks(N[:3])
axins4.set_yticks(N[:3])
axins4.tick_params(labelsize=8)
axins4.set_aspect('equal')
cbaxes4 = inset_axes(axins4, width="5%", height="100%", loc='right',
                    bbox_to_anchor=(0.15, 2.52, 1, 1),
                    bbox_transform=axins.transAxes,
                    borderpad=0)

cbar4 = plt.colorbar(c4, cax=cbaxes4)
cbar4.ax.tick_params(labelsize=6)
No description has been provided for this image
In [1]:
!pip install pandoc

Defaulting to user installation because normal site-packages is not writeable
Collecting pandoc
  Downloading pandoc-2.4.tar.gz (34 kB)
  Installing build dependencies: started
  Installing build dependencies: finished with status 'done'
  Getting requirements to build wheel: started
  Getting requirements to build wheel: finished with status 'done'
  Preparing metadata (pyproject.toml): started
  Preparing metadata (pyproject.toml): finished with status 'done'
Collecting plumbum (from pandoc)
  Downloading plumbum-1.9.0-py3-none-any.whl.metadata (10 kB)
Collecting ply (from pandoc)
  Downloading ply-3.11-py2.py3-none-any.whl.metadata (844 bytes)
Requirement already satisfied: pywin32 in c:\users\giachero\appdata\local\packages\pythonsoftwarefoundation.python.3.12_qbz5n2kfra8p0\localcache\local-packages\python312\site-packages (from plumbum->pandoc) (308)
Downloading plumbum-1.9.0-py3-none-any.whl (127 kB)
Downloading ply-3.11-py2.py3-none-any.whl (49 kB)
Building wheels for collected packages: pandoc
  Building wheel for pandoc (pyproject.toml): started
  Building wheel for pandoc (pyproject.toml): finished with status 'done'
  Created wheel for pandoc: filename=pandoc-2.4-py3-none-any.whl size=34897 sha256=74b9aee51037cb9c4443f4cb3b07780f097c225d852ca1ba0889f2b40ea9990a
  Stored in directory: c:\users\giachero\appdata\local\packages\pythonsoftwarefoundation.python.3.12_qbz5n2kfra8p0\localcache\local\pip\cache\wheels\9c\2f\9f\b1aac8c3e74b4ee327dc8c6eac5128996f9eadf586e2c0ba67
Successfully built pandoc
Installing collected packages: ply, plumbum, pandoc
Successfully installed pandoc-2.4 plumbum-1.9.0 ply-3.11

[notice] A new release of pip is available: 25.0.1 -> 25.1.1
[notice] To update, run: C:\Users\giachero\AppData\Local\Microsoft\WindowsApps\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\python.exe -m pip install --upgrade pip
In [ ]: