added calculation of the adjasent sideband sucsectibilty in for wave mixing

4 files changed, 38 insertions, 13 deletions

diff --git a/field_description.m b/field_description.m
index ffe6bcc..5d8aa77 100644
--- a/field_description.m
+++ b/field_description.m
@@ -1,13 +1,16 @@
 1;
 
 %EM field definition
-Ep=0.01;
+Ep=0.01; %probe
 Epc=conj(Ep);
-Ed=1;
+Ed=.1; %drive
 Edc=conj(Ed);
+Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
+Emc=conj(Em);
 wd=w13;
 wp=w12;
-modulation_freq=[0,  wp, wd, -wp, -wd, wp-wd, wd-wp];
-E_field        =[0,  Ep, Ed, Epc, Edc, 0,     0    ];
+wm=wd-(wp-wd);
+modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
+E_field        =[0,  Ep, Ed, Em,   Epc, Edc, Emc,  0,     0    ];
 %modulation_freq=[0,  wp, -wp];
 %E_field        =[0,  Ep,  Epc];
diff --git a/liouville.m b/liouville.m
index 4812775..6fe83be 100644
--- a/liouville.m
+++ b/liouville.m
@@ -1,3 +1,6 @@
+1;
+t0 = clock (); % we will use this latter to calculate elapsed time
+
 
 % load useful functions;
 useful_functions;
@@ -19,17 +22,19 @@ Nfreq=length(modulation_freq);
 %tune probe frequency
 detuning_p=0;
 N_detun_steps=100;
-detuning_p_min=-15;
+detuning_p_min=-1;
 detuning_p_max=-detuning_p_min;
 detuning_freq=zeros(1,N_detun_steps+1);
 kappa_p      =zeros(1,N_detun_steps+1);
+kappa_m      =zeros(1,N_detun_steps+1);
 detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps;
 for detuning_p_cntr=1:N_detun_steps+1;
 
 wp0=w12;
 detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1);
 wp=wp0+detuning_p;
-modulation_freq=[0,  wp, wd, -wp, -wd, wp-wd, wd-wp];
+wm=wd-(wp-wd);
+modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
 
 % now we create Liouville indexes list
 [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
@@ -49,22 +54,27 @@ L=Liouville_operator_matrix(
 		N, L,
 		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
 
+
 %solving for density matrix vector
 rhoLiouville=L\rhoLiouville_dot;
 
 
-rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq);
-rho_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq);
-rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq);
-%rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq)
+%rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); % 0 frequency, 
+%rho_p=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); % probe frequency
+%rho_d=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); % drive frequency
+%rho_m=rhoOfFreq(rhoLiouville, 4, Nlevels, Nfreq); % opposite sideband frequency
 
-kappa_p(detuning_p_cntr)=sum(sum(rho_1));
+kappa_p(detuning_p_cntr)=sucseptibility(2, rhoLiouville, dipole_elements, Nlevels, Nfreq);
+kappa_m(detuning_p_cntr)=sucseptibility(4, rhoLiouville, dipole_elements, Nlevels, Nfreq);
 detuning_freq(detuning_p_cntr)=detuning_p;
 
 %kappa_p_re=real(kappa_p);
 %kappa_p_im=imag(kappa_p);
 
 endfor
-figure(1); plot(detuning_freq, real(kappa_p));
-figure(2); plot(detuning_freq, imag(kappa_p));
+figure(1); plot(detuning_freq, real(kappa_p)); title("probe dispersion");
+figure(2); plot(detuning_freq, imag(kappa_p)); title("probe absorption");
+figure(3); plot(detuning_freq, real(kappa_m)); title("off resonant sideband dispersion");
+figure(4); plot(detuning_freq, imag(kappa_m)); title("off resonant absorption");
 
+elapsed_time = etime (clock (), t0)
diff --git a/useful_constants.m b/useful_constants.m
index 115efb4..82d1d1a 100644
--- a/useful_constants.m
+++ b/useful_constants.m
@@ -3,3 +3,4 @@
 % some physical constants
 hbar=1;
 im_one=0+1i;
+el_charge=1;
diff --git a/useful_functions.m b/useful_functions.m
index 5522afb..2ae4ced 100644
--- a/useful_functions.m
+++ b/useful_functions.m
@@ -141,3 +141,14 @@ function [rhoLiouville_dot, L]=constran_rho_and_match_L(
 	rhoLiouville_dot(1)=1;
 endfunction
 
+% calculate sucseptibility for the field at given frequency index
+function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
+
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq);
+	kappa=0;
+	for i=1:Nlevels
+		for j=1:Nlevels
+			kappa+=dipole_elements(j,i)*rho(i,j);
+		endfor
+	endfor
+endfunction