Data Structure

Replication data for said article is provided in the form of a HDF5 file with the NetCDF4 file format. This choice allows us to embed the metadata and hierarchical structure right within the dataset. The data is best read with the use of the `xarray` Python library, with minimal examples below showcasing usage and code segments for replicating article plots.

Dependencies

• Python
• numpy
• matplotlib
• xarray
• lmfit

Usage

As mentioned, the data structure has been carefully prepared for top transparency and ease of exploration in mind. Its loading can be done as follows:

import xarray as xr

dt = xr.open_datatree("ARTICLE_DATA.h5")

print(dt)

and gives us an immediate overview of the entire data structure:

Group: /

├── Group: /HWP

│    Dimensions:       (θ: 8, z: 3420, uncertainty: 2, Δ1269: 80)

│    Coordinates:

│     * θ          (θ) int64 64B -24 -12 0 12 24 36 48 60

│     * z          (z) float64 27kB -19.3 -19.29 -19.28 ... 18.79 18.8

│     * uncertainty     (uncertainty) int64 16B 0 1

│     * Δ1269        (Δ1269) float64 640B -58.5 -57.0 -55.5 ... 58.5 60.0

│    Data variables:

│      Rabi Fits      (θ, z, uncertainty) float64 438kB ...

│      Points to Evaluate (z) float64 27kB ...

│      Fluorescence    (θ, Δ1269, z) uint16 4MB ...

└── Group: /REFLECTOR

    Dimensions:    (mirror_d: 14, z: 3420, uncertainty: 2, Δ1269: 80)

    Coordinates:

     * mirror_d   (mirror_d) float64 112B -1.2 -0.9 -0.45 -0.3 ... 0.9 1.05 1.2

     * z       (z) float64 27kB -19.3 -19.29 -19.28 ... 18.77 18.79 18.8

     * uncertainty  (uncertainty) int64 16B 0 1

     * Δ1269     (Δ1269) float64 640B -58.5 -57.0 -55.5 ... 57.0 58.5 60.0

    Data variables:

      Rabi Fits   (mirror_d, z, uncertainty) float64 766kB ...

      Fluorescence (mirror_d, Δ1269, z) float32 15MB ...

In interactive environments, for example a Jupyter notebook, the data structure features rich HTML previews with interacting controls, making this a recommended environment for processing the attached dataset.

The HWP group features a dataset with fluorescence data for half waveplate attenuated measurements. This includes both the fluorescence data from the camera, after dark count subtraction and vertical averaging, as well as master equation Rabi fits results performed with the methodology from the article. The fits have and added uncertainty flag dimension, as we store both the nominal value as well as its uncertainty obtained from the fit.

Similarily, the REFLECTOR group contains analogus data, but with the HIPS mirror position being varied instead of the half waveplate angle.

All of the parameters are preserved as Xarray coordinates, with long names and units for plotting.

Imports

import matplotlib.pyplot as plt

import matplotlib as mpl

import numpy as np

import xarray as xr

from lmfit import Model, Parameters, minimize

dt = xr.open_datatree("ARTICLE_DATA.h5")

Figure 3 plot

points_to_evaluate = dt["HWP/Points to Evaluate"].dropna("z")

dt["HWP/Rabi Fits"].sel(uncertainty=0, z=slice(-15, 15)).plot.line(hue="θ")

for z in points_to_evaluate:

plt.axvline(z.item())

plt.show()

Figure 4 plot

This example showcases also the multifit of a composite model. As described in the article, the model components share all parameters apart from the scaling so a custom objective function is required.

points_to_evaluate = dt["HWP/Points to Evaluate"].dropna("z")

da_at_points = dt["HWP/Rabi Fits"].sel(z=points_to_evaluate).isel(uncertainty=0)

err_at_points = dt["HWP/Rabi Fits"].sel(z=points_to_evaluate).isel(uncertainty=1)

def objective(params: Parameters, da: xr.DataArray, err_da: xr.DataArray, model: Model):

"""Calculate total residual for fits of a model sum to several data sets."""

resid = xr.zeros_like(da)

# make residual per data set

for i in range(da.z.size):

components = model.components[i].eval(params=params, angle=da.θ.values)

resid.isel(z=i).values[:] = (da.isel(z=i).values - components) / err_da.isel(

z=i

).values

# now flatten this to a 1D array, as minimize() needs

return resid.values.flatten()

def hwp_sine(angle, amplitude, shift, c):

return amplitude * (np.sin(4 * np.deg2rad(angle) + shift) + c)

model_arr = [Model(hwp_sine, prefix=f"z{i}_") for i in range(da_at_points.z.size)]

model = sum(model_arr[1:], start=model_arr[0])

fit_params = model.make_params()

for i in range(da_at_points.z.size):

fit_params.add(

f"z{i}_amplitude",

value=(

da_at_points.isel(z=i).max().item() - da_at_points.isel(z=i).min().item()

)

/ 2,

)

if i == 0:

fit_params.add(f"z{i}_shift", value=1.93)

fit_params.add(

f"z{i}_c",

value=3.5,

# vary=False,

)

else:

fit_params.add(f"z{i}_shift", expr="z0_shift")

fit_params.add(f"z{i}_c", expr="z0_c")

out = minimize(

objective,

fit_params,

args=(

da_at_points,

err_at_points,

model,

),

nan_policy="omit",

)

hwp_space = np.linspace(da_at_points.θ.min().item(), da_at_points.θ.max().item(), 200)

components = model.eval_components(params=out.params, angle=hwp_space)

for i in range(da_at_points.z.size):

color = mpl.rcParams["axes.prop_cycle"].by_key()["color"][i % len(mpl.rcParams["axes.prop_cycle"])]

plt.plot(

hwp_space,

components[f"z{i}_"],

color=color,

)

plt.errorbar(

da_at_points.θ,

da_at_points.isel(z=i),

yerr=err_at_points.isel(z=i),

fmt=".",

label=f"z={da_at_points.z.isel(z=i).item():.2f} mm",

markersize=5,

ecolor="grey",

elinewidth=1,

capsize=2,

color=color,

)

plt.legend(title=r"$\mathcal{R}(θ, z)$ [MHz]")

plt.xlabel(r"HWP Angle $\Theta$ [deg]")

plt.ylabel(r"$\hat\mathcal{R}(θ, z)$ [MHz]")

plt.show()

Figure 5

This code cell needs to be executed after figure 4 as it reuses the fit data.

for i in range(da_at_points.z.size):

plt.plot(

sorted(

model.components[i].eval(params=out.params, angle=da_at_points.θ.values)

),

sorted(da_at_points.isel(z=i).values),

label=f"z={da_at_points.z.isel(z=i).item():.2f} mm",

)

plt.legend()

plt.show()

Figure 6

fig, axs = plt.subplots(

nrows=2,

ncols=2,

sharey="row",

sharex=True,

height_ratios=[1, 1],

layout="constrained",

)

# Pick mirror positions to plot

ds_min = dt["REFLECTOR"].sel(mirror_d=1.05)

ds_max = dt["REFLECTOR"].sel(mirror_d=-0.3)

for i, ds in enumerate([ds_min, ds_max]):

fax, eax = axs.T[i]

# Fluorescence map

pcm = ds["Fluorescence"].plot.imshow(

ax=fax, interpolation="nearest", add_colorbar=False, cmap=cm.rainbow4

)

fax.set_xlabel("")

# AT fits

eax.errorbar(

ds.z,

ds["Rabi Fits"].sel(uncertainty=0),

yerr=ds["Rabi Fits"].sel(uncertainty=1),

label="Rabi fit",

)

eax.set_xlabel("Distance from camera focus point [mm]")

axs.T[i][0].set_title(f"HIPS mirror position: {ds.mirror_d.item():.2f} mm")

fig.colorbar(pcm, ax=axs[0, :], label="Fluorescence [a.u.]")

axs[1, 0].set_ylabel(r"$\Omega_{\text{mm}}$ [MHz]")

axs[0, 1].set_ylabel("")

axs[1, 0].legend()

fig.show()

Figure 7

mirror_slice = dt["REFLECTOR/Rabi Fits"].sel(mirror_d=[-0.3, 0.9, 1.2], uncertainty=0)

mirror_slice.plot.line(x="z", hue="mirror_d")

plt.show()

Funding

Funded by the European Union (Grant Agreement: 101120422). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or REA. Neither the European Union nor the granting authority can be held responsible for them. The "Quantum Optical Technologies" (FENG.02.01-IP.05-0017/23) project is carried out within the Measure 2.1 International Research Agendas programme of the Foundation for Polish Science, co-financed by the European Union under the European Funds for Smart Economy 2021-2027 (FENG). This research was funded in whole or in part by the National Science Centre, Poland grant No. 2021/43/D/ST2/03114.

(2026-04-07)