%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 and molecular density %in the trap. %In this program only a dc electric field, applied to the axial endcap %electrodes of a typical linear Paul trap, is considered. The molecules are %ionised from a homogeneous background gas. %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 %Constants P_0 = 1; %probability multiplier alpha = 6.12 * 100; %diabaticity constant, m^{-1/2} V^{-1/2} N = 2; %number of photons for excitation step of REMPI process M = 1; %number of photons for ionisation step of REMPI process w_N = 50e-6; %width of excitation laser at centre of trap, m w_M = 50e-6; %width of ionisation laser at centre of trap, m x_d = 0e-6; %displacement of laser focusses in x direction, m y_d = 0e-6; %displacement of laser focusses in y direction, m z_d = 0e-6; %displacement of laser focusses in z direction, m %Radial confinement omega_z = 2 * pi * 200e3; %axial secular frequency, Hz m = 2 * 2.3258671e-26; %mass of nitrogen molecule, kg Q = 1.6e-19; %charge of nitrogen ion created, C %Rayleigh lengths lambda_N = 255e-9; %wavelength of excitation laser, m lambda_M = 212e-9; %wavelength of ionisation laser, m z_RN = pi * w_N^2 / (2 .* lambda_N); %Rayleigh length of excitation laser, m z_RM = pi * w_M^2 / (2 .*lambda_M); %Rayleigh length of ionisation laser, m %Make array of Delta values to determine ionisation probability for K = 1000; %number of wavenumber shifts to find ionisation probability for Delta = linspace(-200000, 0, K); %array of wavenumber shifts, m^{-1} %Make empty array to hold ionisation threshold distributions I = zeros(K); for i=1:K %The a and b parameters from the circular parameterisation of the electric field a = (Delta(i).^2 ./ alpha.^2) .* ((2 .* Q) ./ (m .* omega_z.^2)); b = (Delta(i).^2 ./ alpha.^2) .* ((Q) ./ (m .* omega_z.^2)); %Make function for integration f = @(u, v) P_0 .* sqrt(a.^2.*b.^2.*sin(u).^2.*sin(v).^4 + a.^2.*b.^2.*sin(v).^4.*cos(u).^2 + (-a.^2.*sin(u).^2.*sin(v).*cos(v) - a.^2.*sin(v).*cos(u).^2.*cos(v)).^2) .* (1./(1 + (b.*cos(v) - z_d).^2./z_RM.^2)).^M.*(1./(1 + (b.*cos(v) - z_d).^2./z_RN.^2)).^N.*exp((M./w_M.^2 + N./w_N.^2).*(-(a.*sin(u).*sin(v) - y_d).^2 - (a.*sin(v).*cos(u) - x_d).^2)); %Integrate function between u = 0 and 2\pi and v = 0, pi I(i) = integral2(f, 0, 2*pi, 0, pi, 'AbsTol', 1e-17, 'Method', 'iterated'); end %Plot results figure() plot(Delta/100, I/max(I)) xlabel('Wavenumber shift from zero-field ionisation threshold/ cm^{-1}') ylabel('Normalised ionisation probability') %Print runtime toc