%This program produces the ionisation threshold distrbution, weighted by %the probability of ionisation at different electric field strengths in the %trap. This weighting depends on the laser profiles. %In this program only a rf electric field, applied to the radial %electrodes of a typical linear Paul trap, is considered. The distribution %of molecules is assumed to be constant across the spatial extent of the %ionisation lasers. %By changing the parameters listed below, the inhomogeneous broadening of %the ionisation threshold for different trapping and laser parameters can %be explored. %Authors: Laura Blackburn (laura.blackburn@sussex.ac.uk), Matthias Keller %A more detailed explanation of the methods used can be found here: %www.nature.com/articles/10.1038/s41598-020-74759-6 %Record start time tic %Parameters %Basic parameters alpha = 6.1*100; %constant from adiabaticity of process, V^{-1/2} m^{-1/2} Omega = 2 * pi * 20e6; %frequency of rf trap drive field, 2\pi Hz q = 0.1; %trap q parameter m = 2 * 2.3258671e-26; %mass of single molecule of ^{14}N_2, kg e = 1.6e-19; %charge of produced ion, C %Parameters for ionisation probability equation P_0 = 1e300; %arbitrary probability multiplication factor - this factor is very large as a brute-force solution to numerical integration artefacts N = 2; %number of photons for excitation step in REMPI process M = 1; %number of photons for ionisation step in REMPI process w_N = 50e-6; %spatial width of REMPI laser pulse for excitation step, m w_M = 50e-6; %spatial width of REMPI laser pulse for ionisation step, m tau_N = 5e-9; %temporal width of REMPI laser pulse for excitation step, m tau_M = 5e-9; %temporal width of REMPI laser pulse for ionisation step, m x_d = 0e-6; %x displacement of lasers from trap centre, m t_d = 1.25e-8; %delay between rf start and laser pulse arrival time, s %Voltage on radial electrodes (dc = 0) V_rad_rf = (q .* m .* Omega.^2) ./ (2 * e); %Find E D = 500; %number of wavenumber shifts to simulate for Delta = linspace(-150000, 0, D); %array of wavenumber shifts from the zero-field ionisation threshold for the lowest-lying state of the ion, m^{-1} E = Delta.^2 / alpha.^2; %electric field needed for each wavenumber shift, V m^{-1} %Make empty array to hold ionisation threshold distribution I = zeros(D); %The coefficient a from parameterisation of ellipsoid surface %a = ((E(i)) ./ (V_rad_rf.* cos(Omega .* t))); for i = 1:D %Function to integrate f = @(u, t) P_0.*sqrt(((E(i)) ./ (V_rad_rf.* cos(Omega .* t))).^2.*sin(u).^2 + ((E(i)) ./ (V_rad_rf.* cos(Omega .* t))).^2.*cos(u).^2).*exp(-(t - t_d).^2.*(M./tau_M.^2 + N./tau_N.^2)).*exp((M./w_M.^2 + N./w_N.^2).*(-(((E(i)) ./ (V_rad_rf.* cos(Omega .* t))).*sin(u)).^2 - (((E(i)) ./ (V_rad_rf.* cos(Omega .* t))).*cos(u) - x_d).^2)); %Numerical integration du dt to find ionisation probability. I(i) = integral2(f, 0, 2*pi, t_d - 10 * (tau_M + tau_N), t_d + 10 * (tau_M + tau_N), 'AbsTol', 1e-15); end figure() plot(Delta/100, I/max(I)) xlabel('Wavenumber shift from zero-field ionisation threshold/ cm^{-1}') ylabel('Normalised ionisation probability') %Print time taken to run program toc