I am currently participating in research on transient absorption spectroscopy and four wave mixing. In the experimental design, an extreme ultraviolet (XUV) laser pulse and infrared (IR) laser pulse are sent through a finite gas sample with a certain time delay between them. My current task is to develop a function that, given data on the system and input pulses, provides the transmission spectrum of the XUV pulse electric field after exiting the sample.
The Procedure
The gist of the procedure is to take two matrices describing the sample susceptibility, combine them into a single matrix, diagonalize said matrix, normalize the left and right eigenvector matrices according to each other, apply an exponential function to the eigenvalue matrix, and finally compute the exiting spectrum given an input spectrum.
The Code
function E_out = FiniteSampleTransmissionXUV(w, X_0, X_nd, E_in, a, N, tau)
% Function takes as input:
% > w - A vector that describes the frequency spectrum range
% > X_0 - A vector that describes the susceptibility
% without the IR pulse
% > X_nd - A matrix that describes the susceptibility with the
% IR pulse
% > E_in - A vector that describes the frequency spectrum of the
% incoming XUV pulse electric field
% > a - A constant that describes the intensity of the IR pulse
% > N - A constant that describes the optical depth of the
% sample
% > tau - A constant that describes the time delay between the
% IR pulse and the XUV pulse
%
% Function provides output:
% > A vector that describes the frequency spectrum of the
% XUV pulse electric field at the end of the sample
% determines number of frequency steps from input frequency range
n = size(w); % constant
% determines frequency step size
delta_w = 1 / (n(2) - 1); % constant
% create matrix sqrt(w_i*w_j)
sqrt_w = w.^(1/2).' * w.^(1/2); % nxn matrix
% combine X_0 and X_nd into total suscptibility matrix
X_ij = (a * delta_w * sqrt_w .* X_nd) + ...
diag(diag(sqrt_w).' .* (X_0)); % nxn matrix
% diagonalize X_ij where sum(U_R_i^2) = 1
[U_R, LAMBDA, U_L] = eig(X_ij); % nxn matrices
% attain the function that scales U_L'*U_R
F = sum(U_L' * U_R, 1); % row vector
sqrt_F = F.^2;
% scale U_R and U_L so that U_L'*U_R = 1
U_R_bar = U_R ./ sqrt_F; % nxn matrix
U_L_bar = U_L ./ sqrt_F; % nxn matrix
% apply exponential transform to eigenvalues % diagonal nxn matrix
exp_LAMBDA = diag(exp(1i * (2*pi*N / (3*10^8)) * diag(LAMBDA)));
% create phase shift due to the time delay
tau_phase = exp(1i * w * tau); % row vector
% recombine transformed susceptibility matrix
ULAMBDAU = U_R_bar * exp_LAMBDA * U_L_bar'; % nxn matrix
% apply effects of finite propagation and pulse
% time delay to input electric field spectrum
E_out = ULAMBDAU * (tau_phase .* E_in).'; % vector
end
Testing
The following is a script to test and demonstrate the function.
n = 100; % number of frequency steps
w = linspace(0,1,n); % linearly spaced 1D array
X_0 = (1 + 1i) * (sin(5*pi*w)).^2;
X_nd = (1 + 1i) * ((sin(1*pi*w')) * sin(1*pi*w)).^2;
E_in = GaussianSpectrum(w, 1, 0.5, 0.1);
a = 1;
N = 10^8;
tau = 0;
E_out = FiniteSampleTransmissionXUV(w, X_0, X_nd, E_in, a, N, tau);
figure(1)
plot(w, E_in,'b')
hold on
plot(w, abs(E_out),'r')
hold off
function E_w_x0 = GaussianSpectrum(w, E_0, w_0, sigma)
% creates a quasi-gaussian spectrum defined by the input
% frequency range w, amplitude E_0, center frequency w_0,
% and spectral width sigma
E_w_x0 = E_0 * sin(pi*w) .* ...
(exp(-((w - w_0) / (2^(1/2) * sigma)).^2));
end
In running the test, we see the frequency spectrum of the XUV pulse before entering the sample in blue with the exiting spectrum in red. One sees that frequencies of the XUV pulse have been absorbed, while new frequencies have been generated.
Please let me know what violates best practices and how the function could be refactored to be more robust and clean. Thank you!
