University of Groningen Spectroscopy of Trapped ¹³⁸Ba⁺ Ions for Atomic Parity Violation and Optical Clocks Dijck, Elwin

(1)

Spectroscopy of Trapped ¹³⁸Ba⁺ Ions for Atomic Parity Violation and Optical Clocks

Dijck, Elwin

DOI:

10.33612/diss.108023683

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2020

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Dijck, E. (2020). Spectroscopy of Trapped ¹³⁸Ba⁺ Ions for Atomic Parity Violation and Optical Clocks.

University of Groningen. https://doi.org/10.33612/diss.108023683

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the

author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately

and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the

number of authors shown on this cover page is limited to 10 maximum.

(2)

endices

Additional Appendix

Optical Bloch Equations Code

This appendix contains the matlab code to calculate fluorescence spectra of heavy

alkaline earth metal ions based on the optical Bloch equations detailed in Section 1.3.3.

Our code is based on the work by Oberst

1

, extending it to include several new features:

• Easy and transparent switching between different ions, laser setups and

calcula-tion models (i.e. three-level or eight-level);

• Support for arbitrary propagation directions and polarization states of both

lasers (instead of assuming linear polarization under a fixed angle to the B-field);

• Allowing also contra-propagating lasers instead of assuming co-propagation

(relevant for the effect of ion motion on spectra);

• Using fast matrix exponentiation in transient calculations if the Liouville operator

is time-independent;

• Possibility to sum over detuning classes or oscillator modes to model the effect

of a thermal (Maxwell–Boltzmann) distribution;

• Infrastructure for fitting experimental data using nlinfit (non-linear fitter).

In addition, two bugs in the code were fixed:

• Time coordinate in transient calculation corrected (factor 2π);

• Micromotion velocity amplitude in steady-state calculation corrected (factor 2).

We include several examples reproducing figures from this and other works to show

how the various functions and parameters in the code are used.

This appendix is part of: E. A. Dijck, Spectroscopy of Trapped138Ba+ Ions for Atomic Parity Violation and Optical Clocks, PhD Thesis, University of Groningen (2020). The code files are

available for download with the online version of this thesis, doi:10.33612/diss.108023683.

1_{H. Oberst, Resonance Fluorescence of Single Barium Ions, Diplomarbeit, University of Innsbruck}

(3)

Core code files

The core of the calculation is formed by model_lambda3 and model_lambda8 (adapted

from Oberst’s b311 and b811). These functions return the Liouville operator as a

matrix so that the time derivative of the density matrix (written as a column vector)

is given by matrix multiplication:

function_{L = model_lambda3(S1, S2, l1, l2, D1, D2, ~, ~, ~, lshift , atomdataf, ~)} % Calculate Liouville operator for a three- level system in Lambda configuration

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 based on code by H. Oberst %

% S1, S2 ... saturation parameters defined as I/Isat % l1, l2 ... laser half-linewidths [MHz]

% D1, D2 ... laser detunings [MHz]

% lshift ... optional: extra level shifts (e.g. light shift ) [MHz] % atomdataf ... function providing atomic properties (decay constants) % Level labeling:

% |1>: S % |2>: P % |3>: D

[~, g1, ~, g2] = atomdataf(); % atomic decay constants [MHz] O1 = g1∗sqrt(S1/2); % Rabi frequencies [MHz]

O2 = g2∗sqrt(S2/2);

E = eye(3); ii = (1:3)'∗ones(1,3);

i1 = reshape(ii',1,9); % First index (row) i2 = reshape(ii,1,9); % Second index (column) % Coherent Hamiltonian (rotating frame) [MHz]

H = [D1 -O1/2 0; -O1/2 0 -O2/2; 0 -O2/2 D2];

% Extra shifts of the energy levels (e.g. light shift ) if size_{( lshift ) == 3, H = H + diag(lshift); end} % Relaxation terms

C1 = sqrt(g1) ∗ E(:,1)∗E(2,:); % Radiative decay C2 = sqrt(g2) ∗ E(:,3)∗E(2,:);

Cl1 = sqrt(2∗l1) ∗ E(:,1)∗E(1,:) ; % Laser linewidth Cl2 = sqrt(2∗l2) ∗ E(:,3)∗E(3,:) ;

CC = C1'∗C1 + C2'∗C2 + Cl1'∗Cl1 + Cl2'∗Cl2;

% Effective Hamiltonian and conjugate

H = H -1i/2∗CC; Hc = H'.';

% Liouville [MHz]

L = -1i∗(H(i1,i1).∗E(i2, i2) - Hc(i2,i2).∗E(i1, i1));

% Feeding terms L = L + C1(i1,i1).∗C1(i2,i2); L = L + C2(i1,i1).∗C2(i2,i2); L = L + Cl1(i1,i1).∗Cl1(i2, i2); L = L + Cl2(i1,i1).∗Cl2(i2, i2); end

(4)

endices

Core code files

functionL = model_lambda8(S1, S2, l1, l2, D1, D2, pol1, pol2, u, lshift , atomdataf, setupdataf) % Calculate Liouville operator for an eight-level system in Lambda configuration

% (e.g. an alkaline earth metal ion with Zeeman splitting, no nuclear spin)

% S1, S2 ... saturation parameters defined as I/Isat % l1, l2 ... laser half-linewidths [MHz]

% D1, D2 ... laser detunings [MHz] % pol1, pol2 ... laser polarization states

% u ... magnetic field vector [MHz] u = B ∗ (muB/h) = B ∗ 1.4 MHz/G % lshift ... optional: extra level shifts (e.g. light shift ) [MHz]

% atomdataf ... function providing atomic properties (decay constants) % setupdataf ... function providing setup properties (laser propagation vector) % Level labeling:

% |1>: S1/2 m-1/2 |5>: D3/2 m-3/2 % |2>: S1/2 m+1/2 |6>: D3/2 m-1/2 % |3>: P1/2 m-1/2 |7>: D3/2 m+1/2 % |4>: P1/2 m+1/2 |8>: D3/2 m+3/2

[~, g1, ~, g2] = atomdataf(); % Atomic decay constants [MHz] O1 = g1∗sqrt(S1/2); % Rabi frequencies [MHz]

O2 = g2∗sqrt(S2/2);

gS = 2; gP = 2/3; gD = 4/5; % Lande factors gJ E = eye(8); ii = (1:8)'∗ones(1,8);

i1 = reshape(ii',1,64); % First index i2 = reshape(ii,1,64); % Second index % Atomic Hamiltonian (rotating frame) [MHz]

H = diag([D1 D1 0 0 D2 D2 D2 D2] + norm(u) ...

∗ [( -1/2)∗gS (+1/2)∗gS (-1/2)∗gP (+1/2)∗gP (-3/2)∗gD (-1/2)∗gD (+1/2)∗gD (+3/2)∗gD]);

% Extra shifts of the energy levels (e.g. light shift ) if length(lshift) == 8, H = H + diag(lshift); end

[ ltheta , lphi , lcoprop] = setupdataf(); % Define laser propagation direction

% Determine normalized transverse polarization (Jones) vectors (Etheta, Ephi) % from polarization parametrization (angle, phase). E.g. for a laser propagating % in -y direction, this corresponds to [-Ez Ex]

er1 = [cos(pol1(1)) - 1i∗sin(pol1(1))∗sin(pol1(2)) ...

sin(pol1(1)) + 1i∗cos(pol1(1))∗sin(pol1(2))] / sqrt(1 + sin(pol1(2))^2);

er2 = [cos(pol2(1)) - 1i∗sin(pol2(1))∗sin(pol2(2)) ...

sin(pol2(1)) + 1i∗cos(pol2(1))∗sin(pol2(2))] / sqrt(1 + sin(pol2(2))^2); % Normalized electric field vectors in lab frame

e1 = [er1(1)∗cos(ltheta)∗cos(lphi)-er1(2)∗sin(lphi) ...

er1(1)∗cos(ltheta)∗sin(lphi)+er1(2)∗cos(lphi) -er1(1)∗sin(ltheta) ]; e2 = lcoprop ∗ [er2(1)∗cos(ltheta)∗cos(lphi)-er2(2)∗sin(lphi) ...

er2(1)∗cos(ltheta)∗sin(lphi)+er2(2)∗cos(lphi) -er2(1)∗sin(ltheta) ];

% The magnetic field defines the z-axis, so rotate electric fields to match

ax = cross(u, [0 0 1]) ; tht = asin(norm(ax) / norm(u)); ax = ax / norm(ax);

if _{tht > 0}

e1 = cos(tht)∗e1 + sin(tht)∗(cross(ax,e1)) + (1-cos(tht))∗dot(ax,e1)∗ax; e2 = cos(tht)∗e2 + sin(tht)∗(cross(ax,e2)) + (1-cos(tht))∗dot(ax,e2)∗ax;

(5)

% Interaction Hamiltonian [MHz]

q3 = sqrt(3);

HSP = O1/2 ∗ [dot(e1, [0 0 1/q3]) dot(e1, [ -1/q3 -1i/q3 0]) ; ... dot(e1, [ -1/q3 1i/q3 0]) dot(e1, [0 0 -1/q3]) ];

HDP = O2/2 ∗ [dot(e2, [1/2 1i/2 0]) 0 ; ...

dot_{(e2, [0 0 1/q3])} dot_{(e2, [1/2/q3 1i/2/q3 0]); ...} dot_{(e2, [ -1/2/q3 1i/2/q3 0]) dot(e2, [0 0 1/q3])} _{; ...}

0 dot_{(e2, [ -1/2 1i/2 0])} _];

% Add interaction Hamiltonian

H(1:2,3:4) = HSP; H(3:4,1:2) = HSP'; H(5:8,3:4) = HDP; H(3:4,5:8) = HDP';

% Relaxation terms (Lindblad terms)

C1 = (2/3∗g1)^0.5∗E(:,1)∗E(4,:); C2 = (2/3∗g1)^0.5∗E(:,2)∗E(3,:); C3 = (1/3∗g1)^0.5∗(E(:,1)∗E(3,:) - E(:,2)∗E(4,:)); C4 = (1/2∗g2)^0.5∗E(:,5)∗E(3,:) + (1/6∗g2)^0.5∗E(:,6)∗E(4,:); C5 = (1/6∗g2)^0.5∗E(:,7)∗E(3,:) + (1/2∗g2)^0.5∗E(:,8)∗E(4,:); C6 = (1/3∗g2)^0.5∗(E(:,6)∗E(3,:) + E(:,7)∗E(4,:)); C7 = (2∗l1)^0.5∗(E(:,1)∗E(1,:) + E(:,2)∗E(2,:));

C8 = (2∗l2)^0.5∗(E(:,5)∗E(5,:) + E(:,6)∗E(6,:) + E(:,7)∗E(7,:) + E(:,8)∗E(8,:)); CC = C1'∗C1 + C2'∗C2 + C3'∗C3 + C4'∗C4 + C5'∗C5 + C6'∗C6 + C7'∗C7 + C8'∗C8;

% Effective Hamiltonian (including non-feeding damping terms)

H = H -1i/2 ∗ CC; Hc = H'.';

% Liouville [MHz]

L = -1i∗(H(i1,i1).∗E(i2, i2) - Hc(i2,i2).∗E(i1, i1));

% Feeding terms L = L + C1(i1,i1).∗C1(i2,i2); L = L + C2(i1,i1).∗C2(i2,i2); L = L + C3(i1,i1).∗C3(i2,i2); L = L + C4(i1,i1).∗C4(i2,i2); L = L + C5(i1,i1).∗C5(i2,i2); L = L + C6(i1,i1).∗C6(i2,i2); L = L + C7(i1,i1).∗C7(i2,i2); L = L + C8(i1,i1).∗C8(i2,i2); end

Both preceding functions accept the same set of arguments, including atomdataf

which should be a function providing the decay constants of the ion of interest, as well

as the transition wavelengths (used in calculating the effect of ion motion):

function[lam1, g1, lam2, g2] = atom_caion % Define calcium ion (Ca+) atomic properties

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

lam1 = 396.8; % cooling laser wavelength [nm] lam2 = 866.2; % repump laser wavelength [nm]

g1 = 22.5; % P->S decay constant [MHz] source: Appl. Phys. B 81, 5 (2005)

g2 = 1.35; % P->D decay constant [MHz] end

(6)

endices

Core code files

function[lam1, g1, lam2, g2] = atom_srion % Define strontium ion (Sr+) atomic properties

g1 = 20.4; % P->S decay const [MHz] source: New J. Phys. 18, 123021 (2016)

g2 = 1.18; % P->D decay const [MHz] end

function[lam1, g1, lam2, g2] = atom_baion % Define barium ion (Ba+) atomic properties

g1 = 14.8; % P->S decay const [MHz] source: Phys. Rev. A 82, 012506 (2010)

g2 = 5.38; % P->D decay const [MHz] source: J. Opt. Soc. Am. B 10, 1330 (1993) end

function_{[lam1, g1, lam2, g2] = atom_raion} % Define radium ion (Ra+) atomic properties

g1 = 16.9; % P->S decay const [MHz] source: Phys. Rev. A 79, 052512 (2009)

g2 = 2.64; % P->D decay const [MHz] end

The argument setupdataf of model_lambda8 should be a function defining the

prop-agation direction of both laser beams, for example under a right angle to the z-axis:

function[ltheta, lphi , lcoprop] = setup_rightangle % Define co-propagating laser direction as -y

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Physics spherical coordinates in the frame of the magnetic field

ltheta = 90 ∗ pi/180; lphi = -90 ∗ pi/180; lcoprop = +1;

end

function_{[ltheta, lphi , lcoprop] = setup_rightangle_contra} % Define contra-propagating laser direction as -y

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Physics spherical coordinates in the frame of the magnetic field

ltheta = 90 ∗ pi/180; lphi = -90 ∗ pi/180; lcoprop = -1;

end

With the Liouville operator, transient or steady-state density matrices can be

calculated. The former is implemented in the function transient (adapted from

(7)

Oberst’s dndynfa), which uses either numerical integration or matrix exponentiation,

depending on whether the Liouville operator is time-dependent (due to ion motion):

function_{rt = transient(par, t, r0)}

% Calculate density matrices as function of time

% par ... parameter set

% t ... vector of time samples [us]

% r0 ... initial density matrix (as column vector)

S1 = par{1}; % saturation parameter pump laser

S2 = par{2}; % saturation parameter repump laser

D1 = par{3}; % detuning pump laser [MHz]

D2 = par{4}; % detuning repump laser [MHz]

l1 = par{5}; % half-linewidth pump laser [MHz]

l2 = par{6}; % half-linewidth repump laser [MHz]

pol1 = par{7}; % polarization state pump laser (angle, phase)

pol2 = par{8}; % polarization state repump laser (angle, phase)

u = par{9}; % magnetic field vector [MHz]

dv = par{10}; % velocity amplitude [m/s]

tf = par{11}; % trap frequency [MHz]

modelf = par{12}; % model function: 3-level or 8-level

atomdataf = par{13}; % atom data function

setupdataf = par{14}; % experimental setup data function

lshift = par{15}; % optional: extra energy shifts of levels [MHz]

[lam1, ~, lam2, ~] = atomdataf(); [~, ~, lcoprop] = setupdataf();

% Liouville matrix components

L0 = modelf(S1, S2, l1, l2 , 0, 0, pol1, pol2, u, lshift , atomdataf, setupdataf); dL1 = modelf(S1, S2, l1, l2 , 1, 0, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0; dL2 = modelf(S1, S2, l1, l2 , 0, 1, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0; L = L0 + D1∗dL1 + D2∗dL2;

if dv > 0

% Numerically integrate the Liouville equation ...

DL = 1000∗dv∗(1/lam1∗dL1 + lcoprop ∗ 1/lam2∗dL2); th = [0 t(end)];

[th, rhoh] = ode45('transient_dvec', th, r0', [], L, DL, tf); rt = interp1(th, rhoh, t, ' spline ') ';

else

% ... or use matrix exponentiation (faster when operator is time-independent)

rt = zeros(length(L0), length(t));

for_{kk = 1:length(t)}

rt (:, kk) = expm(L ∗ 2∗pi∗t(kk)) ∗ r0; % Note factor 2pi

end end end

function_{rprime = transient_dvec(t, r, ~, L, DL, tf)} % Calculate the derivative of the density matrix vector

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 based on code by H. Oberst

rprime = 2∗pi∗(L + DL∗cos(2∗pi∗t∗tf))∗r; % Note overall factor 2pi

end

Calculating steady-state density matrices as function of detuning (fluorescence spectra)

is implemented in steadystate_spec (adapted from Oberst’s anspecfa2) with helper

functions for only changing the pump (‘blue’) or repump (‘red’) laser detuning:

(8)

endices

Core code files

functionrd = steadystate_spec(par, D1, D2)

% Calculate steady-state density matrices for varying laser detuning(s) % The effect of the oscillation of the ion in the trap is considered using % the matrix continued fraction formalism.

% D1 ... vector of pump laser detuning settings [MHz] % D2 ... vector of repump laser detuning settings [MHz]

D1 = par{3} + D1; % detuning pump laser [MHz]

D2 = par{4} + D2; % detuning repump laser [MHz]

dv = par{10}; % velocity amplitude [m/s]

tf = par{11}; % trap frequency [MHz]

L0 = modelf(S1, S2, l1, l2 , 0, 0, pol1, pol2, u, lshift , atomdataf, setupdataf); dL1 = modelf(S1, S2, l1, l2 , 1, 0, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0; dL2 = modelf(S1, S2, l1, l2 , 0, 1, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0;

% Frequency modulation amplitude [MHz] and max modulation index to consider

DL = 1000∗dv ∗ (1/lam1 ∗ dL1 + lcoprop ∗ 1/lam2 ∗ dL2) / 2;

maxmodidx = 0; if abs(dv)>0, maxmodidx = ceil(2∗1000∗abs(dv)/(lam1∗tf)); end z = length(L0); zh = sqrt(z); E = eye(z);

rd = zeros(z, length(D1));

formm = 1:length(D1)

L = L0 + D1(mm)∗dL1 + D2(mm)∗dL2; % Liouville matrix [MHz] Ah = L;

if abs(dv) > 0 % Calculate S+ and S- for n=maxmodidx

Sp = -(L - 1i∗(maxmodidx+1)∗tf∗E)\DL; Sm = -(L + 1i∗(maxmodidx+1)∗tf∗E)\DL;

fornn = maxmodidx:-1:1 % Iteratively calculate S+ and S- for n=0

Sp = -(L - 1i∗nn∗tf∗E + DL∗Sp)\DL; Sm = -(L + 1i∗nn∗tf∗E + DL∗Sm)\DL;

end

Ah = Ah + DL∗(Sp+Sm);

end

% Add normalization condition and solve for steady-state

v = zeros(1,z); v (1,1:( zh+1):z) = ones(1,zh); A = [v; Ah]; c = [1 zeros(1,z) ]';

rd (:, mm) = A\c;

end end

(9)

functionrd = steadystate_specblue(par, D1)

% Calculate steady-state density matrices for varying pump laser detuning % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

rd = steadystate_spec(par, D1, zeros(size(D1)));

end

functionrd = steadystate_specred(par, D2)

% Calculate steady-state density matrices for varying repump laser detuning % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

rd = steadystate_spec(par, zeros(size(D2)), D2);

end

The function steadystate_spec includes the possible effect of a single coherent

oscillation on the density matrix (e.g. excess micromotion of a trapped ion). The

following second version allows specifying any number of motion components with

differing velocity amplitudes and frequencies (e.g. for including macromotion and

intrinsic micromotion):

functionrd = steadystate_spec_multimotion(par, D1, D2, dv, tf)

% Calculate steady-state density matrices including multiple motion components % Note: the usual micromotion velocity in the parameter set is ignored!

% D1 ... vector of pump laser detuning settings [MHz] % D2 ... vector of repump laser detuning settings [MHz] % dv ... vector of motion velocities [m/s]

% tf ... vector of motion frequencies [MHz]

D1 = par{3} + D1; % detuning pump laser [MHz]

D2 = par{4} + D2; % detuning repump laser [MHz]

L0 = modelf(S1, S2, l1, l2 , 0, 0, pol1, pol2, u, lshift , atomdataf, setupdataf); dL1 = modelf(S1, S2, l1, l2 , 1, 0, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0; dL2 = modelf(S1, S2, l1, l2 , 0, 1, pol1, pol2, u, lshift , atomdataf, setupdataf) - L0;

% Frequency modulation amplitudes [/um] and max modulation indices

DL = 1000∗(1/lam1 ∗ dL1 + lcoprop∗1/lam2 ∗ dL2) / 2; maxmodidx = ceil(2∗1000∗abs(dv)./(lam1∗tf)); z = length(L0); zh = sqrt(z); E = eye(z);

(10)

endices

Core code files

rd = zeros(z, length(D1));

formm = 1:length(D1)

L = L0 + D1(mm)∗dL1 + D2(mm)∗dL2; % Liouville matrix [MHz] Ah = L;

forvv = 1:length(dv) % Calculate S+ and S- for n=maxmodidx

Sp = -(L - 1i∗(maxmodidx(vv)+1)∗tf(vv)∗E)\(dv(vv)∗DL); Sm = -(L + 1i∗(maxmodidx(vv)+1)∗tf(vv)∗E)\(dv(vv)∗DL);

for_{nn = maxmodidx(vv):-1:1 % Iteratively calculate S+ and S- for n=0}

Sp = -(L - 1i∗nn∗tf(vv)∗E + dv(vv)∗DL∗Sp)\(dv(vv)∗DL); Sm = -(L + 1i∗nn∗tf(vv)∗E + dv(vv)∗DL∗Sm)\(dv(vv)∗DL);

end

Ah = Ah + dv(vv)∗DL∗(Sp+Sm);

end

% Add normalization condition and solve for steady-state

v = zeros(1,z); v (1,1:( zh+1):z) = ones(1,zh); A = [v; Ah]; c = [1 zeros(1,z) ]';

rd (:, mm) = A\c;

end end

A few convenience functions are provided to extract values of interest (populations

or coherences) from an array of calculated density matrices:

functionfl = get_s_population(rr)

% Return row matrix of summed S-level populations % from an array of density matrices (3-level or 8-level). % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Location of S-level populations in the density matrix if size_{(rr ,1) == 3∗3, slevels = [1];}

elseif size_{(rr ,1) == 8∗8, slevels = [1 10];} end

fl = real(sum(rr(slevels,:) , 1));

end

functionfl = get_p_population(rr)

% Return row matrix of summed P-level populations % from an array of density matrices (3-level or 8-level). % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Location of P-level populations in the density matrix if size_{(rr ,1) == 3∗3, plevels = [5];}

elseif size_{(rr ,1) == 8∗8, plevels = [19 28];} end

fl = real(sum(rr(plevels,:), 1));

end

functionfl = get_d_population(rr)

% Return row matrix of summed D-level populations % from an array of density matrices (3-level or 8-level). % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Location of D-level populations in the density matrix if size_{(rr ,1) == 3∗3, dlevels = [9];}

elseif size_{(rr ,1) == 8∗8, dlevels = [37 46 55 64];} end

fl = real(sum(rr(dlevels,:), 1));

(11)

functionfl = get_sd_coherences_mag(rr) % Return array of S-D coherence magnitudes

% from an array of density matrices (3-level or 8-level). % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

% Take norm of complex values and take average of above and below diagonal % (should be complex conjugates except for numerical uncertainty)

if size_{(rr ,1) == 3∗3} fl = [(sqrt(real(rr (3,:) ).^2 + imag(rr(3,:)).^2) ... + sqrt(real(rr(7,:) ).^2 + imag(rr(7,:)).^2)) / 2]; elseif size_{(rr ,1) == 8∗8} fl = [(sqrt(real(rr (5,:) ).^2 + imag(rr(5,:)).^2) ... + sqrt(real(rr(33,:)).^2 + imag(rr(33,:)).^2)) / 2; ... % |1><5| (sqrt(real(rr (6,:) ).^2 + imag(rr(6,:)).^2) ... + sqrt(real(rr(41,:)).^2 + imag(rr(41,:)).^2)) / 2; ... % |1><6| (sqrt(real(rr (7,:) ).^2 + imag(rr(7,:)).^2) ... + sqrt(real(rr(49,:)).^2 + imag(rr(49,:)).^2)) / 2; ... % |1><7| (sqrt(real(rr (8,:) ).^2 + imag(rr(8,:)).^2) ... + sqrt(real(rr(57,:)).^2 + imag(rr(57,:)).^2)) / 2; ... % |1><8| (sqrt(real(rr (13,:) ).^2 + imag(rr(13,:)).^2) ... + sqrt(real(rr(34,:)).^2 + imag(rr(34,:)).^2)) / 2; ... % |2><5| (sqrt(real(rr (14,:) ).^2 + imag(rr(14,:)).^2) ... + sqrt(real(rr(42,:)).^2 + imag(rr(42,:)).^2)) / 2; ... % |2><6| (sqrt(real(rr (15,:) ).^2 + imag(rr(15,:)).^2) ... + sqrt(real(rr(50,:)).^2 + imag(rr(50,:)).^2)) / 2; ... % |2><7| (sqrt(real(rr (16,:) ).^2 + imag(rr(16,:)).^2) ... + sqrt(real(rr(58,:)).^2 + imag(rr(58,:)).^2)) / 2]; % |2><8| end end

The effect of thermal motion on fluorescence spectra is calculated by summing

over detuning classes (implemented in steadystate_spec_maxwellboltzmann) or

oscillator modes (implemented in steadystate_spec_maxwellboltzmann_osc):

function_{rd = steadystate_spec_maxwellboltzmann(par, D1, D2, vscale, nbins)} % Apply Maxwell-Boltzmann (thermal) velocity distribution to steady-state spectrum % by averaging over Doppler detuning classes

% PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 %

% par ... parameter set

% D1 ... vector of pump laser detuning settings [MHz] % D2 ... vector of repump laser detuning settings [MHz] % vscale ... sqrt(kT/m): scale factor of the M-B distribution [m/s] % nbins ... number of velocity classes to include

[lam1, ~, lam2, ~] = par{13}(); % laser wavelengths [nm] [~, ~, lcoprop] = par{14}(); % co- or contra-propagation?

% The 1D (projected) Maxwell-Boltzmann velocity distribution is a normal % distribution with std. dev. = sqrt(kT/m). Get velocity classes of equal % probability by using the inverse CDF.

cdfs = conv(linspace(0, 1, nbins+1), [0.5 0.5], ' valid ');

dv = vscale∗sqrt(2) ∗ erfinv(2∗cdfs - 1); % velocity classes [m/s]

% Initialize with first velocity class

D1v = 1000 ∗ dv(1) / lam1; % Doppler detuning [MHz]

D2v = lcoprop ∗ 1000 ∗ dv(1) / lam2; % Doppler detuning [MHz] rd = (1/nbins) ∗ steadystate_spec(par, D1+D1v, D2+D2v);

(12)

endices

Consistency check code files

for ii = 2:nbins

D1v = 1000 ∗ dv(ii) / lam1; % Doppler detuning [MHz]

D2v = lcoprop ∗ 1000 ∗ dv(ii) / lam2; % Doppler detuning [MHz] rd = rd + (1/nbins) ∗ steadystate_spec(par, D1+D1v, D2+D2v);

end end

function_{rd = steadystate_spec_maxwellboltzmann_osc(par, D1, D2, T, m, fsec, nbins)} % Apply Maxwell-Boltzmann (thermal) velocity distribution to steady-state spectrum % by averaging over harmonic oscillator modes

% D1 ... vector of pump laser detuning settings [MHz] % D2 ... vector of repump laser detuning settings [MHz] % T ... temperature [mK]

% m ... ion/atom mass [kg] % fsec ... secular frequency [MHz]

% nbins ... number of velocity classes to include

% The effect of a thermal distribution can also be calculated by summing % over 1D harmonic oscillator modes (the secular frequency should be small % compared to the linewidth)

kB = 1.380648e-23 / 1000; % Boltzmann constant [J/mK] h = 6.626070e-34 ∗ 1e6; % Planck constant [J/MHz]

cdfs = conv(linspace(0, 1, nbins+1), [0.5 0.5], ' valid ');

n = kB∗T ∗ log(1 ./ (1 - cdfs)) / (h∗fsec) - 1; % oscillator eigenstates dv = sqrt(h∗fsec ./ m ∗ (2∗n+1)); % velocity classes [m/s] % Initialize with first velocity class

rd = (1/nbins)∗steadystate_spec_multimotion(par, D1, D2, dv(1), fsec);

% Sum over rest of distribution for _{ii = 2:nbins}

rd = rd + (1/nbins)∗steadystate_spec_multimotion(par, D1, D2, dv(ii), fsec);

end end

Consistency check code files

To verify the internal consistency and normalization of calculation results, we include

a few tests. Part of these rely on running the calculation on a simplified two-level

system, which is done by defining an ‘atom’ with one decay constant set to zero:

function[lam1, g1, lam2, g2] = atom_twolevel_10meg

% Set P->D decay to zero to form a two-level atom with 10 MHz decay width % Note: this produces warnings about rank deficient matrices, but it does work % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

lam1 = 0; % not used lam2 = 0;

g1 = 10; % P->S decay constant [MHz]

g2 = 0; % P->D decay constant [MHz] end

(13)

The first test verifies the normalization of the time coordinate in transient and the

definition of the Rabi frequencies in model_lambda3 and model_lambda8 by plotting

transient calculation results for spontaneous decay and Rabi oscillations:

functiontest_transient_equations

% Verify the definition of time coordinate and Rabi frequency % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Select either model_lambda3 or model_lambda8:

modelf = @model_lambda3;

%modelf = @model_lambda8;

% Take two-level atom so we know what the result should be

atomf = @atom_twolevel_10meg;

par = {0, 0, 0, 0, 0, 0, [0; 0], [0; 0], [0 0 0], 0, 0, modelf, atomf, @setup_rightangle, []}; [~, g, ~, ~] = atomf();

% Initial density matrix

if isequal_{(modelf, @model_lambda3)}

sstates = [1]; pstates = [5]; r0 = zeros(9, 1); else % model_lambda8 sstates = [1 10]; pstates = [19 28]; r0 = zeros(64, 1); end

% Check time normalization with decay rate from P level

r0p = r0; r0p(pstates) = 1 / length(pstates); t = linspace(0, 0.1, 25); % Time in microseconds

par{1} = 0; % Saturation parameter: zero laser intensity

rt = transient(par, t, r0p);

figure_{; hold on;}

plot_{(t, get_s_population(rt)); plot(t, get_p_population(rt));}

plot_{(t, exp(-t∗(2∗pi∗g)), 'o'); % Plot expected decay curve (g = 10MHz)} xlabel_{('Time␣(\mus)'); ylabel('Population');}

legend_{('S␣population', 'P␣population', 'Expected␣decay␣curve');} % Check Rabi frequency definition vs. calculated Rabi flops

r0s = r0; r0s( sstates ) = 1 / length(sstates); t = linspace(0, 0.02, 50); % Time in microseconds

par{1} = 300; % Saturation parameter: large enough for some oscillations

rt = transient(par, t, r0s);

figure; hold on;

plot(t, get_s_population(rt)); plot(t, get_p_population(rt)); if isequal(modelf, @model_lambda3)

% Rabi flopping matching the model_lambda3 Rabi frequency definition

Oflop = g ∗ sqrt(par{1}/2);

else

% Rabi flopping matching the model_lambda8 Rabi frequency definition % with a factor sqrt(3) from coupling strength

Oflop = g ∗ sqrt(par{1}/2) / sqrt(3);

end

plot_{(t, 0.5+0.5∗cos(2∗pi∗t∗Oflop), 'o');} xlabel_{('Time␣(\mus)'); ylabel('Population');}

legend_{('S␣population', 'P␣population', 'Expected␣Rabi␣oscillations');} end

(14)

endices

Consistency check code files

The following test shows the definition of the laser saturation parameters in

model_lambda3 and model_lambda8 by comparing steady-state calculation results for

a two-level system with the expected saturation curves of excited state population

and power broadening:

functiontest_saturation_equations

% Check the definition of laser saturation parameters

% Using atom_twolevel_10meg gives a two-level system where we know the % expected saturation curves, allowing us to check the definition of the % saturation parameters.

% Select either model_lambda3 or model_lambda8:

modelf = @model_lambda3;

%modelf = @model_lambda8;

% Varying laser power (only laser 1 for two-level system)

S1 = 0.1:0.1:5;

% Define fitted range of spectrum sufficiently wide

x = linspace(-100, +100, 200); % [MHz]

% Let's ignore these warnings due to forcing a two-level system

wrdm = warning('query', 'MATLAB:rankDeficientMatrix');

warning('off', 'MATLAB:rankDeficientMatrix');

par = @(s1) {s1, 0, 0, 0, 0, 0, [0; 0], [0; pi/2], 1.3996∗[0 0 0.1], 0, ... [], modelf, @atom_twolevel_10meg, @setup_rightangle, []};

for _{ii = 1:length(S1)} y{ii} = get_p_population(steadystate_specblue(par(S1(ii)), x)); % Do a simple FWHM determination halfmax = max(y{ii}) / 2; w = []; for _{jj = 2:length(y{ii})}

if sign(y{ii}( jj ) - halfmax) ~= sign(y{ii}(jj-1) - halfmax)

w(end+1) = abs((halfmax-y{ii}(jj-1))/(y{ii}(jj)-y{ii}( jj -1))) ... ∗ x( jj ) + abs((halfmax-y{ii}(jj))/(y{ii}( jj )-y{ii}( jj -1)))∗x(jj -1); end end fwhm(ii) = w(2) - w(1); chm(ii) = (w(1) + w(2))/2; fm(ii) = 2∗halfmax; end figure;

subplot(1,3,1); hold on;

for ii = 1:floor((length(S1)-1)/2):length(S1) % show just 3 spectra plot_{(x, get_p_population(steadystate_specblue(par(S1(ii)), x)), ...}

'LineWidth', 1, 'Color', [0 .7∗ ii /length(S1) 0]);

end

xlabel_{('Laser␣detuning␣(MHz)'); ylabel('P-level␣population');} title_{('Calculated␣spectra');}

subplot_{(1,3,2); hold on;} plot_{(S1, fm, 'k. ');} xlim([0 ceil(max(S1))]);

(15)

% The fits require the Curve Fitting Toolbox

v = ver;

if any(strcmp({v.Name}, 'Curve␣Fitting␣Toolbox'))

% a is the saturated excited state population, x/b is the usual saturation parameter

fo = fitoptions('Method', 'NonlinearLeastSquares', 'Lower', [0 0], ... 'Upper', [inf inf ], 'StartPoint' , [max(fm) 1]);

ft = fittype('a∗(1-1/(1+x/b))', 'Options', fo); f = fit(S1', fm', ft );

cv = coeffvalues(f);

plot_{(f, 'b-- ');}

plot_{([cv(2) cv(2) ], ylim, 'b-- ');}

legend_{('Simple␣max.', 'Fit␣saturation␣curve', 'Location', 'Southeast');} else

legend('Simple␣max.', 'Location', 'Southeast'); end

xlabel('Laser␣saturation␣param'); ylabel('Peak␣height'); title('Peak␣height');

subplot(1,3,3); hold on; plot(S1, fwhm, 'k.');

xlim([0 ceil(max(S1))]); ylim(min(ylim, [0 inf])) if any(strcmp({v.Name}, 'Curve␣Fitting␣Toolbox'))

% a is the natural linewidth, x/b is the usual saturation parameter

fo = fitoptions('Method', 'NonlinearLeastSquares', 'Lower', [0 0], ... 'Upper', [inf inf ], 'StartPoint' , [15 1]) ;

ft = fittype('a∗sqrt(1+x/b)', 'Options', fo); f = fit(S1', fwhm', ft);

cv = coeffvalues(f);

plot_{(f, 'b-- ');}

plot_{([cv(2) cv(2) ], ylim, 'b-- ');}

legend_{('Simple␣FWHM', 'Fit␣saturation␣curve', 'Location', 'Southeast');} else

legend_{('Simple␣FWHM', 'Location', 'Southeast');} end

xlabel('Laser␣saturation␣param'); ylabel('FWHM␣width␣(MHz)'); title('Peak␣width');

warning(wrdm.state, 'MATLAB:rankDeficientMatrix');

% The saturation parameters are defined as S = I/Isat, this gives the % expected behavior with model_lambda3 and model_lambda8, where in the % latter case the electric field interacts only with 1/3 of the full dipole % moment, so the effective saturation intensity is higher by a factor 3. end

The next function checks whether the definition of the velocity amplitude parameter

is consistent between the transient (transient) and steady-state (steadystate_spec)

calculations, which use different methods of modeling the effect of ion motion:

function_{test_micromotion_equations}

% Verify the consistency of the micromotion velocity amplitude parameter % between transient and steady-state calculations

% Calculation parameters, varying only laser 2 detuning and velocity amplitude

(16)

endices

Consistency check code files

par = @(D2off, dv) {2∗(10/g1)^2, 2∗(5/g2)^2, -20, D2off, 0, 0, [], [], ... [], dv, tf , @model_lambda3, @atom_baion, @setup_rightangle, []};

figure('Position' , [100 100 1600 300]); subplot(1,3,2); hold on;

D2 = -50:0.1:50;

plot_{(D2, get_p_population(steadystate_specred(par(0, 0), D2)), 'b');} plot_{(D2, get_p_population(steadystate_specred(par(0, 5), D2)), 'g');} plot_{(D2, get_p_population(steadystate_specred(par(0, 10), D2)), 'r');} xlim_{([-50 50]) ; ylim([0 0.1]) ;}

xlabel_{('\Delta_2␣(MHz)'); ylabel('Steady-state␣\rho_{22}');} % Now calculate the transient solution at two points in the spectrum

x1 = -0.8; x2 = 3.5;

line([x1 x1], ylim, 'Color', 'k'); annotation('Arrow', [.51 .47], [.8 .8]) ; line([x2 x2], ylim, 'Color', 'k'); annotation('Arrow', [.53 .57], [.8 .8]) ; % Integrate to 1 microsecond starting with population in S state

t = linspace(0, 1, 1000); r0 = [1; zeros(8, 1) ];

% Define running average function for plotting the time-averaged solution

windowsz = ceil(1/tf / (t(2)-t(1))); b = (1/windowsz)∗ones(1,windowsz); a = 1; pt1 = get_p_population(transient(par(x1, 0), t, r0)); pt2 = get_p_population(transient(par(x1, 5), t, r0)); pt3 = get_p_population(transient(par(x1, 10), t, r0));

subplot_{(1,3,1); hold on;} plot_{(t, pt1, 'b');}

plot_{(t, pt2, 'g: '); plot(t, filter (b, a, pt2), 'g');} plot_{(t, pt3, 'r: '); plot(t, filter (b, a, pt3), 'r');} ylim_{([0 0.1]) ;}

xlabel('t␣(\mus)'); ylabel('\rho_{22}');

pt1 = get_p_population(transient(par(x2, 0), t, r0)); pt2 = get_p_population(transient(par(x2, 5), t, r0)); pt3 = get_p_population(transient(par(x2, 10), t, r0));

subplot(1,3,3); hold on; plot(t, pt1, 'b');

plot(t, pt2, 'g: '); plot(t, filter (b, a, pt2), 'g'); plot(t, pt3, 'r: '); plot(t, filter (b, a, pt3), 'r'); ylim([0 0.1]) ;

xlabel('t␣(\mus)'); ylabel('\rho_{22}');

% Both calculations should result in the same time-averaged steady-state % solution for the same parameter set (i.e. equal micromotion amplitude). end

The final test verifies that the calculated effect of a Maxwell–Boltzmann distribution

for a two-level system produces the expected Voigt lineshape:

function_{test_maxwellboltzmann_equations}

% Verify the lineshapes produced by the steady-state Maxwell-Boltzmann calculations % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

(17)

% expected Lorentz and Voigt profiles % Set up our two-level atom for testing

g = 10; % 10 MHz decay width

lam = 500; % 500 nm transition wavelength atomf = @() deal(lam, g, inf, 0);

% Set a low laser intensity to prevent saturation broadening

par = {g/1000, 0, 0, 0, 0, 0, [], [], [], 0, [], @model_lambda3, atomf, @setup_rightangle, []};

% Settings for calculating Doppler broadening (sum over velocity classes)

T = 100; % temperature [mK]

m = 50 ∗ 1.660539e-27; % atomic mass [kg]

kB = 1.380648e-23 / 1000; % Boltzmann constant [J/mK] vbins = 100; % number of velocity bins to use

% Doppler broadening can also be calculated by summing over harmonic % oscillator modes (secular frequency should be small compared to linewidth)

fsec = 0.1; % secular frequency [MHz]

fbins = 20; % number of oscillator modes to use

D = -25:1:25; % detuning [MHz]

% Let's ignore these warnings due to forcing a two-level system

warning('off', 'MATLAB:rankDeficientMatrix');

yp1 = get_p_population(steadystate_spec(par, D, zeros(size(D)))); yp2 = get_p_population(steadystate_spec_maxwellboltzmann(par, D, ...

zeros_{(size(D)), sqrt(kB∗T/m), vbins));}

yp3 = get_p_population(steadystate_spec_maxwellboltzmann_osc(par, D, ...

zeros_{(size(D)), T, m, fsec , fbins));}

warning_{(wrdm.state, 'MATLAB:rankDeficientMatrix');} figure_{; hold on;}

plot(D, yp1, 'bo'); plot(D, yp2, 'ro'); plot(D, yp3, 'r. ');

% Check results with the expected Lorentz and Voigt profiles

lorentz = @(x, mu, gam) 1/(pi∗gam) ∗ gam^2./((x-mu).^2+gam^2); gaus = @(x, mu, sig) 1/sqrt(2∗pi∗sig^2) ∗ exp(-(x-mu).^2./(2∗sig^2));

% Use a somewhat larger range in the convolution to avoid edge effects

Dconv = min(D)-10:1:max(D)+10; voigt = conv(lorentz(Dconv, 0, g/2), ...

gaus(Dconv, 0, 1000∗sqrt(kB∗T/m)/lam), 'same'); voigt = voigt(1+10:end-10);

plot_{(D, max(yp1)∗pi∗g/2 ∗ lorentz(D, 0, g/2), 'b');} plot_{(D, max(yp2)/max(voigt) ∗ voigt, 'r');} xlabel_{('Frequency␣detuning␣\Delta\nu␣(MHz)');} ylabel_{('P-level␣population');}

legend_{('Spectrum␣(0␣mK)', 'Spectrum␣(100␣mK,␣vel.␣classes)', ...}

'Spectrum␣(100␣mK,␣osc.␣modes)', 'Analytic␣Lorentz␣profile', ... 'Analytic␣Voigt␣profile ' , 'Location', 'Northwest');

(18)

endices

Example code files

Here we provide several more examples of using the code. The first example illustrates

how the figures in this thesis are produced:

% Generate example figures optical Bloch equations from thesis % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

par = {5, 25, - 5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], 1.3996∗[0 0 10], 0, ... [], @model_lambda8, @atom_baion, @setup_rightangle, []};

% - Figure 1.6 (page 42)

---t = linspace(0, 0.5, 201); % ---time (us) sstates = [1 10];

r0 = zeros(64, 1);

r0( sstates ) = 1 / length(sstates); rt = transient(par, t, r0); rd = steadystate_specred(par, [0]);

plot(t, get_s_population(rt), 'g');

plot(xlim, [get_s_population(rd) get_s_population(rd)], 'g--'); plot(t, get_p_population(rt), 'b');

plot(xlim, [get_p_population(rd) get_p_population(rd)], 'b--'); plot(t, get_d_population(rt), 'r');

plot_{(xlim, [get_d_population(rd) get_d_population(rd)], 'r--');} ylim_{([0 1]) ;}

title_{('Thesis␣Figure␣1.6'); xlabel('t␣(\mus)'); ylabel('Populations');} % - Figure 1.7 (page 42) ---D2 = linspace(-50, 50, 401); rd = steadystate_specred(par, D2); figure; subplot(2,1,1); plot(D2, get_p_population(rd));

line([par{4}; par{4}], ylim, 'LineStyle' , ' - . '); ylabel('\rho_{33}␣+␣\rho_{44}');

ysdr = [real(rd (5,:) ); real(rd (6,:) ); real(rd (7,:) ); real(rd (8,:) ); ...

real(rd (13,:) ); real(rd (14,:) ); real(rd (15,:) ); real(rd (16,:) ) ];

ysdi = [imag(rd(5,:)); imag(rd(6,:)); imag(rd(7,:)); imag(rd(8,:)); ...

imag(rd(13,:)); imag(rd(14,:)); imag(rd(15,:)); imag(rd(16,:))]; plot(D2, ysdr); gca.ColorOrderIndex = 1; plot(D2, ysdi, '--');

legend({'S_{-1/2}D_{-3/2}', 'S_{-1/2}D_{-1/2}', 'S_{-1/2}D_{+1/2}', ...

'S_{-1/2}D_{+3/2}', 'S_{+1/2}D_{-3/2}', 'S_{+1/2}D_{-1/2}', ... 'S_{+1/2}D_{+1/2}', 'S_{+1/2}D_{+3/2}'}, 'Location', 'East');

ylim_{([-0.4 0.4]) ;}

xlabel_{('\Delta\nu_{2}␣(MHz)'); ylabel('Coherences');} if exist_{(' sgtitle '), sgtitle('Thesis␣Figure␣1.7'); end}

% - Figure 1.8 (page 44) ---figure;

(19)

rdm = steadystate_specred({0.5∗5, 25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ... 1.3996∗[0 0 10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2); rdp = steadystate_specred({2.0∗5, 25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ...

1.3996∗[0 0 10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2);

plot_{(D2, get_p_population(rd), 'b:');} plot_{(D2, get_p_population(rdm), 'c');} plot_{(D2, get_p_population(rdp), 'b');} ylim_{([0 0.12]) ; gca.YTick = [0 0.04 0.08 0.12];} ylabel_{('\rho_{33}␣+␣\rho_{44}')}

subplot_{(2,2,2); hold on;}

rdm = steadystate_specred({5, 0.5∗25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ... 1.3996∗[0 0 10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2); rdp = steadystate_specred({5, 2.0∗25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ...

plot(D2, get_p_population(rd), 'b:'); plot(D2, get_p_population(rdm), 'c'); plot(D2, get_p_population(rdp), 'b'); ylim([0 0.12]) ; gca.YTick = [0 0.04 0.08 0.12]; subplot(2,2,3); hold on;

rdm = steadystate_specred({5, 25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ...

1.3996∗[0 0 0.75∗10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2); rdp = steadystate_specred({5, 25, -5, 0, 0.1/2, 0.1/2, [0; 0], [0; pi/2], ...

1.3996∗[0 0 1.50∗10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2);

plot(D2, get_p_population(rd), 'b:'); plot_{(D2, get_p_population(rdm), 'b');} plot_{(D2, get_p_population(rdp), 'c');} ylim_{([0 0.12]) ; gca.YTick = [0 0.04 0.08 0.12];}

xlabel_{('\Delta\nu_2␣(MHz)'); ylabel('\rho_{33}␣+␣\rho_{44}')} subplot_{(2,2,4); hold on;}

rdm = steadystate_specred({5, 25, -5, 0, 0.0/2, 0.0/2, [0; 0], [0; pi/2], ...

1.3996∗[0 0 10], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2); rdp = steadystate_specred({5, 25, -5, 0, 0.5/2, 0.5/2, [0; 0], [0; pi/2], ...

plot(D2, get_p_population(rd), 'b:'); plot(D2, get_p_population(rdm), 'c'); plot(D2, get_p_population(rdp), 'b'); ylim([0 0.12]) ; gca.YTick = [0 0.04 0.08 0.12]; xlabel('\Delta\nu_2␣(MHz)');

if exist(' sgtitle '), sgtitle('Thesis␣Figure␣1.8'); end

% - Figure 1.9 (page 44) ---figure;

% Ignore warnings when all population is pumped into dark states

warning_{('off', 'MATLAB:rankDeficientMatrix');}

D2 = linspace(-50, 50, 401);

pol = {[0; 0]∗pi/180, [90; 0]∗pi/180, [0; 90]∗pi/180, [0; -90]∗pi/180};

for _{ii = 1:4} for _{jj = 1:4}

nn = 4∗(ii-1) + jj ;

subplot(4,4,nn);

(20)

endices

Example code files

pol{ ii }, pol{jj }, 1.3996∗[0 0 5], 0, [], @model_lambda8, ... @atom_baion, @setup_rightangle, []}, D2));

plot(D2, pp(nn,:)); ylim([0 0.1]) ; end

end

warning_{(wrdm.state, 'MATLAB:rankDeficientMatrix');} if exist_{(' sgtitle '), sgtitle('Thesis␣Figure␣1.9'); end}

% - Figure 5.4 (page 108) ---% Light shift template, red-detuned laser 649.82nm/6mW/100um/Bz/theta0/phi0

lshifttemplate = [-0.000 -0.000 -0.495 -0.495 -0.000 +0.495 +0.495 -0.000]; D2 = linspace(-50, 50, 201);

rd = steadystate_specred({4, 8, -10, 0, 0, 0, [0; 0], [0; pi/2], 1.3996∗[0 0 8], 0, ... [], @model_lambda8, @atom_baion, @setup_rightangle, []}, D2);

yp = get_p_population(rd); ysd = get_sd_coherences_mag(rd);

% Set scalar P level light shift to +5 MHz

rdshift = steadystate_specred({4, 8, -10, 0, 0, 0, [0; 0], [0; pi/2], ... 1.3996∗[0 0 8], 0, [], @model_lambda8, @atom_baion, @setup_rightangle, ... -5.0∗ lshifttemplate /0.495}, D2);

ypshift = get_p_population(rdshift); ysdshift = get_sd_coherences_mag(rdshift);

figure_;

subplot_{(2,1,1); hold on;}

plot_{(D2, yp, 'r: '); plot(D2, ypshift, 'r');} ylim_{([0 0.06]) ;}

ylabel_{('\rho_{33}␣+␣\rho_{44}');} subplot(2,1,2); hold on;

plot(D2, ysd, ' : '); gca.ColorOrderIndex = 1; plot(D2, ysdshift);

legend({'S_{-1/2}D_{-3/2}', 'S_{-1/2}D_{-1/2}', 'S_{-1/2}D_{+1/2}', ...

'S_{-1/2}D_{+3/2}', 'S_{+1/2}D_{-3/2}', 'S_{+1/2}D_{-1/2}', ... 'S_{+1/2}D_{+1/2}', 'S_{+1/2}D_{+3/2}'}, 'Location', 'East');

ylim([0 0.35]) ;

xlabel('\Delta\nu_2␣(MHz)'); ylabel('Coherences␣(magnitudes)'); if exist(' sgtitle '), sgtitle('Thesis␣Figure␣5.4'); end

%

---The following example reproduces several figures from Oberst’s Diplomarbeit,

show-ing how the parameter values correspond to the definitions in the present code. In

particular, we define the laser saturation parameters differently:

% Reproduce figures from [Oberst, Diplomarbeit, Innsbruck (1999)] % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

% Note that the definition of the saturation parameters S_k has been changed

% to the conventional S_k = 2∗(Omega_k/Gamma_k)^2, where Omega_k and Gamma_k % are the corresponding Rabi frequency and decay constant.

(21)

[~, g1, ~, g2] = atom_baion(); % [MHz] decay constants

% - Figure 2.2 (page 16) ---figure; hold on;

t = 0:0.005:1;

rt = transient(params_simple(0, 0, 2∗(10/g1)^2, 0, 0, 0, [], [], 0, [], ... @model_lambda3), t, [1 0 0 0 0 0 0 0 0]') ;

plot_{(t, get_s_population(rt), 'k'); text(0.1, 0.8, '\rho_{11}');} plot_{(t, get_p_population(rt), 'k'); text(0.2, 0.1, '\rho_{22}');} plot_{(t, get_d_population(rt), 'k'); text(0.4, 0.8, '\rho_{33}');} ylim_{([0 1]) ;}

title_{('Oberst␣Figure␣2.2'); xlabel('t␣(\mus)');}

% - Figure 2.3 (page 16) ---figure; D2 = -100:0.1:100; plot(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ... 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 0, [], @model_lambda3), D2)), 'k'); ylim([0 5]) ;

title('Oberst␣Figure␣2.3'); xlabel('\Delta_r␣(MHz)'); ylabel('\rho_{22}∗100'); % - Figure 2.4 (page 18) ---figure; hold on;

t = 0:0.005:1;

rt = transient(params_simple(0, 0, 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 0, [], ... @model_lambda3), t, [1 0 0 0 0 0 0 0 0]') ;

plot_{(t, get_s_population(rt), 'k'); text(0.8, 0.25, '\rho_{11}');} plot_{(t, get_p_population(rt), 'k'); text(0.2, 0.1, '\rho_{22}');} plot_{(t, get_d_population(rt), 'k'); text(0.8, 0.85, '\rho_{33}');} plot_{(t, get_sd_coherences_mag(rt), 'k'); text(0.8, 0.45, ' |\rho_{13}|');} ylim_{([0 1]) ;}

title_{('Oberst␣Figure␣2.4'); xlabel('t␣(\mus)');}

% - Figure 2.5 (page 19) ---figure; hold on;

D2 = -40:0.1:0; plot(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ... 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 0, [], @model_lambda3), D2)), 'k'); plot(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ... 2∗(10/g1)^2, 2∗(5/g2)^2, 50, 50, [], [], 0, [], @model_lambda3), D2)), 'k'); plot(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ... 2∗(10/g1)^2, 2∗(5/g2)^2, 500, 500, [], [], 0, [], @model_lambda3), D2)), 'k'); ylim([0 5]) ;

title('Oberst␣Figure␣2.5'); xlabel('\Delta_r␣(MHz)'); ylabel('\rho_{22}∗100'); % - Figure 3.2 (page 30) ---figure_; D2 = -80:0.1:80; plot_{(D2, 100∗get_p_population(steadystate_specred(params_simple(-10, 0, ...} 8∗(0.5)^2, 8∗(1.0)^2, 0, 0, 5, 90, 0, [], @model_lambda8), D2)), 'k'); ylim_{([0 5]) ;} title_{('Oberst␣Figure␣3.2');}

xlabel_{('\Delta_r␣(MHz)'); ylabel('(\rho_{33}+\rho_{44})∗100');}

(22)

---endices

Example code files

figure; D2 = -60:0.1:40; plot(D2, 100∗get_p_population(steadystate_specred(params_simple(-10, 0, ... 8∗(0.25)^2, 8∗(1.0)^2, 0, 0, 5, 30, 0, [], @model_lambda8), D2)), 'k'); ylim([0 2.5]) ; title('Oberst␣Figure␣3.3');

xlabel_{('\Delta_r␣(MHz)'); ylabel('(\rho_{33}+\rho_{44})∗100');}

% - Figure 3.4 (page 32) ---% Note: the phase convention has been changed.

figure_;

D2 = -40:0.1:20;

rd = steadystate_specred(params_simple(-10, 0, 8∗(0.5)^2, 8∗(1.0)^2, 0, 0, 5, ... 90, 0, [], @model_lambda8), D2);

subplot(5,1,1); plot(D2, 100∗real(sum(rd([19 28],:), 1)), 'k'); ylim([0 5]) ; title('Oberst␣Figure␣3.4');

plot(D2, real(rd(8,:)), 'k'); plot(D2, imag(rd(8,:)), 'k'); subplot(5,1,3); hold on;

plot(D2, real(rd(13,:)), 'k'); plot(D2, imag(rd(13,:)), 'k'); xlabel('\Delta_r␣(MHz)');

% - Figure 4.2 (page 37) ---% Note: labeling of the time scale and trap frequency both miss a factor 2pi. figure_{; hold on;}

t = (0:0.002:3)/(2∗pi); plot_{(t, 100∗get_p_population(transient(params_simple(-20, -10, 2∗(10/g1)^2, ...} 2∗(5/g2)^2, 0, 0, [], [], 0, 2∗pi∗20, @model_lambda3), t, ... [1 0 0 0 0 0 0 0 0]') ), 'b'); plot(t, 100∗get_p_population(transient(params_simple(-20, -10, 2∗(10/g1)^2, ... 2∗(5/g2)^2, 0, 0, [], [], 10, 2∗pi∗20, @model_lambda3), t, ... [1 0 0 0 0 0 0 0 0]') ), 'g'); ylim([0 7]) ;

title('Oberst␣Figure␣4.2'); xlabel('t␣(\mus)'); ylabel('\rho_{22}∗100'); % - Figure 4.3 (page 38) ---% Note: the code as given in Oberst's Diplomarbeit has a bug where in steady-% state calculations the velocity amplitude is effectively doubled.

D2 = -100:0.1:100; plot_{(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ...} 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 0, 20, @model_lambda3), D2)), 'b'); plot_{(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ...} 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 5, 20, @model_lambda3), D2)), 'g'); plot_{(D2, 100∗get_p_population(steadystate_specred(params_simple(-20, 0, ...} 2∗(10/g1)^2, 2∗(5/g2)^2, 0, 0, [], [], 10, 20, @model_lambda3), D2)), 'r'); ylim_{([0 6]) ;}

title_{('Oberst␣Figure␣4.3'); xlabel('\Delta_r␣(MHz)'); ylabel('\rho_{22}∗100');} %

(23)

---The preceding makes use of the function params_simple which allows to conveniently

work with linear laser polarization under a specified angle to the B-field:

function_{par = params_simple(D1off, D2off, S1, S2, l1, l2, u, alpha, dv, tf , varargin)} % Convert parameter set for simple parametrization

% Specify saturation parameters directly and set the polarization of both % lasers by specifying the angle between the magnetic field and the % (linear) polarization direction .

par = {abs(S1), ... % pump laser saturation parameter abs(S2), ... % repump laser saturation parameter

D1off, ... % pump laser detuning (zero) [MHz]

D2off, ... % repump laser detuning (zero) [MHz] abs(l1) / 1000, ... % pump laser half-linewidth [MHz] abs(l2) / 1000, ... % repump laser half-linewidth [MHz]

[alpha; 0] ∗ pi/180, ... % pump laser polarization state

[alpha; 0] ∗ pi/180, ... % repump laser polarization state

[0 0 u ], ... % magnetic field vector [MHz] abs_{(dv), ...} _{% velocity amplitude [m/s]} abs_{(tf), ...} _{% trap frequency [MHz]}

@model_lambda8, ... % model function

@atom_baion, ... % atomic data function

@setup_rightangle, ... % experimental setup function

[]}; % light shifts [MHz]

if length_{(varargin) >= 1}

par{12} = varargin{1}; % optionally specify model function

end end

The following produces figures similar to those in the work of Lisowski et al.

2

,

showing the use of steadystate_spec_multimotion to calculate the combined effect

of macromotion and micromotion, and the difference that co- and contra-propagating

laser configurations make:

% Reproduce figures from [Lisowski et al., Appl. Phys. B 81, 5 (2005)] % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683 % Calibration parameters

l1 = 1.0; % [MHz] laser 1 full linewidth

l2 = 0.2; % [MHz] laser 2 full linewidth

B = [1.3996 0 0]; % [MHz/G] magnetic field vector calibration

tf = 11.6; % [MHz] trap frequency

fsec = 1.0; % [MHz] secular frequency

[~, g1, ~, g2] = atom_caion(); % [MHz] decay constants

% Note on saturation parameters: it appears Lisowski et al. use a different % normalization of the Rabi frequencies. Additionally, setting a magnetic % field of 0.75 G instead of 1.0 G gives better agreement with the figures. % - Figure 2 ---% The secular motion of a trapped ion produces motion components at the % secular frequency and at the trap frequency +/- the secular frequency

2_{C. Lisowski et al., Dark resonances as a probe for the motional state of a single ion, Appl.}

(24)

endices

Example code files

% (intrinsic micromotion). Additional (excess) micromotion at the trap % frequency may be caused by ion displacement from the trap center. Using % steadystate_spec_multimotion these four motion components can be modeled, % specifying their amplitudes in [m/s].

D2 = linspace(-200, 200, 401);

yp1 = 100∗get_p_population(steadystate_spec_multimotion( ...

{2∗(1.25∗120/g1)^2, 2∗(1.07∗30/g2)^2, -50, 0, l1/2, l2/2, [0; 0], [0; 0], 0.75∗B, ... [], [], @model_lambda8, @atom_caion, @setup_rightangle_contra, []}, ...

zeros_{(size(D2)), D2, [1.0 0.77 0.65 0.0], [ fsec tf+fsec tf - fsec tf ]) );}

zeros(size(D2)), D2, [1.0 0.77 0.65 4.4], [ fsec tf+fsec tf - fsec tf ]) );

zeros(size(D2)), D2, [1.0 0.77 0.65 8.9], [ fsec tf+fsec tf - fsec tf ]) ); figure; hold on;

plot(D2, yp1, 'k: '); plot(D2, yp2, 'k-- '); plot(D2, yp3, 'k-'); xlabel('\Delta_R␣[MHz]');

ylabel('Probability␣of␣population␣of␣4P_{1/2}-state␣[%]'); ylim([0 16]); title('Lisowski2005␣Figure␣2');

% - Figure 3 ---% Varying Rabi frequencies (saturation parameters).

D2 = linspace(-200, 200, 401);

{2∗(1.25∗120/g1)^2, 2∗(1.07∗30/g2)^2, -50, 0, 1/2, 2/2, [0; 0], [0; 0], 0.75∗B, ... [], [], @model_lambda8, @atom_caion, @setup_rightangle_contra, []}, ...

zeros_{(size(D2)), D2, [1.0 0.77 0.65 4.4], [ fsec tf+fsec tf - fsec tf ]) );}

zeros(size(D2)), D2, [1.0 0.77 0.65 4.4], [ fsec tf+fsec tf - fsec tf ]) );

plot(D2, yp1, 'k: '); plot(D2, yp2, 'k-- '); plot(D2, yp3, 'k-'); xlabel('\Delta_R␣[MHz]');

ylabel('Probability␣of␣population␣of␣4P_{1/2}-state␣[%]'); ylim([0 14]); title('Lisowski2005␣Figure␣3');

% - Figure 4 ---% Showing the effect of co-/contra-propagating lasers.

D2 = linspace(-200, 200, 401);

(25)

{2∗(1.25∗60/g1)^2, 2∗(1.07∗15/g2)^2, -50, 20, 1/2, 2/2, [0; 0], [0; 0], 0.75∗B, ... [], [], @model_lambda8, @atom_caion, @setup_rightangle, []}, ...

plot_{(D2, yp1, 'k-'); plot(D2, yp2, 'k-. ');} xlabel_{('\Delta_R␣[MHz]');}

ylabel_{('4S_{1/2}␣-␣4P_{1/2}␣fluorescence␣signal␣[a.u.]');} title_{('Lisowski2005␣Figure␣4');}

%

---The next example produces figures similar to those in the work of Sikorsky et al.

3

.

This illustrates the use of steadystate_spec to calculate the effect of varying both

laser detunings simultaneously (as a Doppler shift) and some micromotion effects:

% Reproduce figures from [Sikorsky et al., Phys. Rev. A 96, 012519 (2017)] % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

% Calibration parameters

[lam1, g1, lam2, g2] = atom_srion(); % wavelengths and decay constants

Bz = [0 0 1.3996]; % [MHz/G] magnetic field vector normalization

tf = 26.51; % [MHz] trap frequency

% The Rabi frequency definitions match Oberst's; the laser linewidths are % described as on the order of hundreds of kHz.

% - Figure 2 ---% Laser half-widths of 500 kHz seem to give decent agreement. Here only the % blue laser frequency is scanned.

D1 = linspace(-80, 80, 201);

yp = get_p_population(steadystate_specblue({8∗(9.3/g1)^2, 8∗(7.0/g2)^2, 0, ... - 9.0, 0.5, 0.5, [0; 0], [pi/4; 0], 3.0∗Bz, 0, tf , @model_lambda8, @atom_srion, ... @setup_rightangle, []}, D1));

figure;

plot(D1, yp, 'b');

xlabel('\delta_{422}␣[MHz]'); ylabel('\rho␣(P_{1/2})'); ylim([0 0.15]); title('Sikorsky2017␣Figure␣2');

% - Figure 3 ---% Now both laser frequencies are scanned to model a Doppler shift.

D = linspace(-80, 80, 201);

yp = get_p_population(steadystate_spec({8∗(9.3/g1)^2, 8∗(7.0/g2)^2, 0, -9.0, ... 0.5, 0.5, [0; 0], [pi/4; 0], 3.0∗Bz, 0, tf , @model_lambda8, @atom_srion, ... @setup_rightangle, []}, D, (lam1/lam2)∗D));

figure_;

plot_{(D, yp, 'Color', [0 0.5 0], 'Linewidth', 3);}

xlabel_{('\delta_D␣[MHz]'); ylabel('\rho␣(P_{1/2})'); ylim([0 0.15]);} title_{('Sikorsky2017␣Figure␣3');}

3_{T. Sikorsky et al., Doppler cooling thermometry of a multilevel ion in the presence of micromotion,}

(26)

endices

Example code files

% - Figure 6 (inset) ---% The dominant effect is (inherent) micromotion in y-direction at the trap % frequency +/- the corresponding secular frequency, each component with a % projected velocity amplitude of about 4.5 m/s. Laser half-linewidth appears % to be about 200 kHz.

fsecy = 0.73; % [MHz] secular frequency D = linspace(-80, 80, 201);

yp1 = get_p_population(steadystate_spec({8∗(12/g1)^2, 8∗(6/g2)^2, 0, -7.5, ... 0.2, 0.2, [0; 0], [pi/4; 0], 3.0∗Bz, 0, tf , @model_lambda8, @atom_srion, ... @setup_rightangle, []}, D, (lam1/lam2)∗D));

yp2 = get_p_population(steadystate_spec_multimotion({8∗(12/g1)^2, ...

8∗(6/g2)^2, 0, - 7.5, 0.2, 0.2, [0; 0], [pi/4; 0], 3.0∗Bz, [], [], @model_lambda8, ... @atom_srion, @setup_rightangle, []}, D, (lam1/lam2)∗D, [4.5 4.5], ...

[ tf+fsecy tf- fsecy ]) );

plot(D, yp1, 'r' , 'Linewidth', 2); plot(D, yp2, 'b', 'Linewidth', 2); xlabel('\delta_D␣[MHz]'); ylabel('\rho␣(P_{1/2})'); ylim([0 0.15]); title('Sikorsky2017␣Figure␣6␣(inset)');

% - Figure 10

---D1 = linspace(-60, 60, 201);

% The excess micromotion only has a single motion component at the trap % frequency. Laser half-linewidth appears to be about 350 kHz.

yp1 = 700∗get_p_population(steadystate_specblue({8∗(15.6/g1)^2, ...

8∗(10.7/g2)^2, - 18.5, - 0.66, 0.35, 0.35, [6; 0]∗pi/180, [35; 0]∗pi/180, 3.0∗Bz, ... 0.0, tf , @model_lambda8, @atom_srion, @setup_rightangle, []}, D1));

plot(D1, yp1, 'b'); plot(D1, yp2, 'g'); plot(D1, yp3, 'r');

xlabel('\delta_{422}␣[MHz]'); ylabel('PMT␣counts␣in␣1␣ms'); ylim([0 120]); title('Sikorsky2017␣Figure␣10');

%

---The last example produces a figure similar to that in the work of Tugayé et al.

4

,

combining a custom Liouville operator model_tugaye2019 with calculating the effect

of a thermal velocity distribution using steadystate_spec_maxwellboltzmann:

% Reproduce figures from [Tugaye et al., Phys. Rev. A 99, 023412 (2019)] % PhD Thesis E. A. Dijck (2020), doi:10.33612/diss.108023683

[~, g1, ~, g2] = atom_srion();

m = 88 ∗ 1.660539e-27; % atomic mass [kg]

kB = 1.380648e-23 / 1000; % boltzmann constant [J/mK]

% - Figure 3