%% Physical constants
PlanckConstant = 6.62607015E-34;
PlanckConstantReduced = 6.62607015E-34/(2*pi);
FineStructureConstant = 7.2973525698E-3;
ElectronMass = 9.10938291E-31;
GravitationalConstant = 6.67384E-11;
ProtonMass = 1.672621777E-27;
AtomicMassUnit = 1.660539066E-27;
BohrRadius = 5.2917721067E-11;
BohrMagneton = 9.274009994E-24;
BoltzmannConstant = 1.38064852E-23;
StandardGravityAcceleration = 9.80665;
SpeedOfLight = 299792458;
StefanBoltzmannConstant = 5.670373E-8;
ElectronCharge = 1.602176634E-19;
VacuumPermeability = 1.25663706212E-6;
DielectricConstant = 8.8541878128E-12;
ElectronGyromagneticFactor = -2.00231930436153;
AvogadroConstant = 6.02214076E23;
ZeroKelvin = 273.15;
GravitationalAcceleration = 9.80553;
VacuumPermittivity = 1 / (SpeedOfLight^2 * VacuumPermeability);
HartreeEnergy = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3
% Dy specific constants
Dy164Mass = 163.929174751*AtomicMassUnit;
Dy164IsotopicAbundance = 0.2826;
DyMagneticMoment = 9.93*BohrMagneton;
%% Lattice spacing
Wavelength = 532e-9;
theta = linspace(1.5, 18.0, 100);
LatticeSpacing = Wavelength ./ (2.*sin((theta*pi/180)/2));
figure(1);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
plot(theta, LatticeSpacing * 1E6, LineWidth=2.0)
xlim([0 19]);
ylim([0.5 21]);
xlabel('Angle (deg)', FontSize=16)
ylabel('Lattice spacing (µm)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(LatticeSpacing * 1E6),1)) ' µm ; \bf Lower bound = ' num2str(round(min(LatticeSpacing * 1E6),1)) ' µm'], FontSize=16)
grid on
%% Scaling of vertical trap frequency with lattice spacing
Wavelength = 532e-9;
a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
Power = 5;
waist_y = 250E-6;
waist_z = 50E-6;
thetas = linspace(1.5, 18.0, 100);
LatticeSpacings = zeros(1, length(thetas));
Omega_z = zeros(1, length(thetas));
for idx = 1:length(thetas)
theta = 0.5 * thetas(idx) .* pi/180;
LatticeSpacings(idx) = Wavelength ./ (2.*sin(theta));
Omega_z(idx) = sqrt(((16 * a * Power) / (pi * Dy164Mass * waist_y * waist_z)) * ...
((2 * (cos(theta)/waist_z)^2) + ((Wavelength * sin(theta)/pi)^2 * ...
((1/waist_y^4) + (1/waist_z^4))) + (pi / LatticeSpacings(idx))^2));
end
nu_z = Omega_z ./ (2*pi);
figure(2);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
plot(LatticeSpacings * 1E6, nu_z * 1E-3, LineWidth=2.0)
xlim([0.5 21]);
xlabel('Lattice spacing (µm)', FontSize=16)
ylabel('Trap frequency (kHz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(nu_z * 1E-3),2)) ' kHz ; \bf Lower bound = ' num2str(round(min(nu_z * 1E-3),2)) ' kHz'], FontSize=16)
grid on
%% Scaling of trap depth with power
a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y = 250E-6;
waist_z = 50E-6;
Powers = linspace(0.1, 5, 100);
TrapDepths = ((8 * a .* Powers) ./ (pi * waist_y * waist_z));
TrapDepthsInHz = TrapDepths ./ PlanckConstant; % in Hz
TrapDepthsInmicroK = TrapDepths ./ BoltzmannConstant; % in µK
TwoPhotonRecoilEnergy = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
TrapDepthsInUnitsOfRecoilEnergy = TrapDepths ./ TwoPhotonRecoilEnergy;
% TrapDepthsToPlot = TrapDepthsInHz * 1E-3;
% units = ' kHz';
% TrapDepthsToPlot = TrapDepthsInmicroK * 1E6;
% units = ' µK';
TrapDepthsToPlot = TrapDepthsInUnitsOfRecoilEnergy;
units = ' E_r';
figure(3);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
plot(Powers, TrapDepthsToPlot, LineWidth=2.0)
xlim([0.0 5.25]);
xlabel('Powers (W)', FontSize=16)
ylabel(['Trap depth (' units ' )'], FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(TrapDepthsToPlot),2)) units '; \bf Lower bound = ' num2str(round(min(TrapDepthsToPlot),2)) units], FontSize=16)
grid on
%% Scaling of the lattice recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy
LatticeSpacing = linspace(2E-6, 20E-6, 100);
RecoilEnergy = PlanckConstant^2 ./ (8 .* Dy164Mass .* LatticeSpacing.^2);
figure(4);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
semilogy(LatticeSpacing * 1E6, RecoilEnergy/PlanckConstant, LineWidth=2.0)
xlim([0.5 21]);
xlabel('Lattice spacing (µm)', FontSize=16)
ylabel('Recoil Energy (Hz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(RecoilEnergy / PlanckConstant),1)) ' Hz; Lower bound = ' num2str(round(min(RecoilEnergy / PlanckConstant),1)) ' Hz'], FontSize=16)
grid on
%% Interference pattern spacing in ToF - de Broglie wavelength associated with the relative motion of atoms
ExpansionTime = linspace(1E-3, 20.0E-3, 100);
figure(5);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
labels = [];
for ls = [2E-6:2E-6:5E-6 6E-6:6E-6:20E-6]
InteferencePatternSpacing = (PlanckConstant .* ExpansionTime) ./ (Dy164Mass * ls);
plot(ExpansionTime*1E3, InteferencePatternSpacing* 1E6, LineWidth=2.0, DisplayName=['\bf Lattice spacing = ' num2str(round(max(ls * 1E6),1)) ' µm'])
hold on
end
xlim([0 22]);
xlabel('Free expansion time (milliseconds)', FontSize=16)
ylabel('Interference pattern period (µm)', FontSize=16)
title('\bf Interference of condensates - Fringe period is the de Broglie wavelength associated with the relative motion of atoms', FontSize=12)
legend(labels, 'Location','NorthWest', FontSize=12);
grid on
legend show
%% Scaling of frequency of oscillation in the first order in Kapitza-Dirac scattering
a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y = 250E-6;
waist_z = 50E-6;
Powers = linspace(0.001, 0.4, 100);
TrapDepths = ((8 * a .* Powers) ./ (pi * waist_y * waist_z));
TwoPhotonRecoilEnergy = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
RabiOscillationFrequency = (1/PlanckConstantReduced) .* (sqrt(TrapDepths.^2/2 + TwoPhotonRecoilEnergy^2));
TrapDepthsInHz = TrapDepths./ PlanckConstant;
TwoPhotonRecoilEnergyInHz = TwoPhotonRecoilEnergy / PlanckConstant;
TrapDepthsInUnitsOfRecoilEnergy = TrapDepthsInHz ./ TwoPhotonRecoilEnergyInHz;
figure(6);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
plot(TrapDepthsInUnitsOfRecoilEnergy, RabiOscillationFrequency .* 1E-3, LineWidth=2.0)
xlim([0 4]);
xlabel('Trap depths (E_r)', FontSize=16)
ylabel('Rabi oscillation frequency (kHz)', FontSize=16)
title(['\bf Upper bound = ' num2str(round(max(RabiOscillationFrequency .* 1E-3),1)) ' kHz; Lower bound = ' num2str(round(min(RabiOscillationFrequency .* 1E-3),1)) ' kHz'], FontSize=16)
grid on
%% Rabi oscillations in the first order in Kapitza-Dirac scattering
a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
waist_y = 250E-6;
waist_z = 50E-6;
Power = 0.2;
TrapDepth = ((8 * a .* Power) ./ (pi * waist_y * waist_z));
TwoPhotonRecoilEnergy = (2*PlanckConstantReduced*2*pi/Wavelength)^2 / (2 * Dy164Mass);
TrapDepthsInUnitsOfRecoilEnergy = TrapDepth ./ TwoPhotonRecoilEnergy;
RabiOscillationFrequency = (1/PlanckConstantReduced) .* (sqrt(TrapDepth.^2/2 + TwoPhotonRecoilEnergy^2));
alpha = TwoPhotonRecoilEnergy / PlanckConstantReduced;
beta = TrapDepth / PlanckConstantReduced;
C = beta^2 / ((2*beta^2) + (4*alpha^2));
PulseDurations = linspace(1E-6, 150E-6, 1000);
PopulationInFirstOrders = C .* (sin(0.5 .* PulseDurations .* (RabiOscillationFrequency)));
figure(7);
set(gcf,'Position',[100 100 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
plot(PulseDurations .* 1E6, PopulationInFirstOrders, LineWidth=2.0, DisplayName=['\bf Power = ' num2str(Power) ' W / Trap depth = ' num2str(round(TrapDepthsInUnitsOfRecoilEnergy, 1)) ' E_r'])
xlabel('Pulse duration (µs)', FontSize=16)
ylabel('Fraction of atoms in first order', FontSize=16)
title('\bf Expected Rabi oscillation', FontSize=16)
grid on
legend(FontSize=16)
%% Scaling of vertical trap frequency with power for fixed lattice spacing
Wavelength = 532e-9;
a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
LatticeSpacing = 13.89e-6;
waist_y = 250E-6;
waist_z = 75E-6;
thetaindeg = 1.0973;
Powers = linspace(0.5, 5.0, 100);
Omega_z = zeros(1, length(Powers));
for idx = 1:length(thetas)
theta = 0.5 * thetaindeg .* pi/180;
Omega_z(idx) = sqrt(((16 * a * Powers(idx)) / (pi * Dy164Mass * waist_y * waist_z)) * ...
((2 * (cos(theta)/waist_z)^2) + ((Wavelength * sin(theta)/pi)^2 * ...
((1/waist_y^4) + (1/waist_z^4))) + (pi / LatticeSpacing)^2));
end
nu_z = Omega_z ./ (2*pi);
figure(8);
set(gcf,'Position',[100 100 950 750])
plot(Powers, nu_z * 1E-3, LineWidth=2.0)
% xlim([0.5 21]);
xlabel('Powers (W)', FontSize=16)
ylabel('Trap frequency (kHz)', FontSize=16)
title(['\bf Maximum = ' num2str(round(max(nu_z * 1E-3),2)) ' kHz ; \bf Minimum = ' num2str(round(min(nu_z * 1E-3),2)) ' kHz'], FontSize=16)
grid on