Contents

Scalar Difference Frequency Generation

Simulates DFG for monochromatic, CW, plane waves.

Input: TolAbs Absolute tolerence used in RungeKutta45 TolRel Relative tolerence used in RungeKutta45 Pp1 Signal seed power Pp2 Idler seed power Pp3 Pump power zsteps Number of steps that shall be returned from the RungeKutta45 zlength Length of crystal QPM Order of QPM structure QPMtype Type of QPM simulation; 'eff' (effeective) or 'domain'

Output: output0 z axis output1 Signal power; output2 Idler power; output3 Pump power output4 Total power output5 Signal photon number output6 Idler photon number output7 Pump photon number output8 Total photon number

function [output0, output1, output2, output3, output4, output5, output6,...
                    output7, output8] = ScalarDFG(TolAbs, TolRel,Pp1,...
                    Pp2, Pp3, zsteps,zlength, QPM, QPMtype)

Independent Parameters

Idler.Lambda  = 1.570e-6;           %wavelength of idler in [m]
Pump.Lambda   = 1.064e-6;           %wavelength of pump in [m]
Param.Length  = zlength;            %Length of nonlinear crystal in [m]
Param.Area    = pi*(40e-6)^2;       %Area of the nonlinear crystal in [m^2]
Param.Temp    = 88;                 %Temperature in degree
Param.QPM     = QPM;                %QPM: 0=NON, 1=1st order, 2=2th order,...
Param.QPMtype = QPMtype;            %Type of QPM: domain og eff.
Param.zsteps  = zsteps;             %Number of z-steps

Error using ScalarDFG (line 34)
Not enough input arguments.

Boundary conditions

Signal.PeakPower    = Pp1;          %power in watt - signal
Idler.PeakPower     = Pp2;          %power in watt - idler
Pump.PeakPower      = Pp3;          %power in watt - pump

constants

const.c           = 299792458;            %Speed of light in [m/s]
const.hbar        = 1.05457173e-34; %const.hbar in [m2*kg/s]
const.epsilon0    = 8.854187817e-12;%in [F/m]
const.mu0         = pi*4e-7;        %in [V*s/A*m]
const.chi         = 22e-12;       %eff. nonlinear suceptibility in [m/V]

Dependent Parameters

Idler.Omega     = 2*pi*const.c / Idler.Lambda;
Pump.Omega      = 2*pi*const.c / Pump.Lambda;
Signal.Omega    = Pump.Omega - Idler.Omega;
Signal.Lambda   = 2*pi*const.c / Signal.Omega;

Signal.Omega / Pump.Omega;
Pump.Lambda / Signal.Lambda;

Refractive index from Sellmeier equations

[n , kVector] = sellmeier([Signal.Lambda Idler.Lambda Pump.Lambda],Param.Temp);
% k vector
Signal.kVector  = kVector(1);
Idler.kVector   = kVector(2);
Pump.kVector    = kVector(3);
% refractive index
Signal.n        = n(1);
Idler.n         = n(2);
Pump.n          = n(3);

deltak  = Pump.kVector - Signal.kVector - Idler.kVector;
lcoh    = 2 * pi / deltak;
if( strcmp(Param.QPMtype , 'eff') )
   const.chi  = 2/(pi*Param.QPM)*const.chi;
   lcoh       = 1e10;
   deltak     = 0;
end

Nonlinear parameters

eta             = sqrt( const.mu0 / const.epsilon0 ) / nthroot(Signal.n*Idler.n*Pump.n,3);
g               = const.epsilon0 * const.chi * sqrt(1/2 * const.hbar * eta^3 ...
            * Signal.Omega * Idler.Omega * Pump.Omega);

%calculate amplitudes
Signal.A    = sqrt(Signal.PeakPower / (const.hbar*Signal.Omega)/Param.Area);
Idler.A     = sqrt(Idler.PeakPower  / (const.hbar*Idler.Omega)/Param.Area);
Pump.A      = 1i*sqrt(Pump.PeakPower   / (const.hbar*Pump.Omega)/Param.Area);

Solve set of ODEs

Sets the global and local maximum fault tolerance.

options = odeset('AbsTol',TolAbs,'RelTol',TolRel);
    % Solve ODE: ode45(solve for , interval , init.values , options).

[z, A]  = ode45(@(z,A) ode45equation(z,A,g,deltak,Param.QPM,lcoh), ...
    linspace(0,Param.Length,Param.zsteps), [Signal.A, Idler.A, Pump.A], options);

Find number of photons instead of amplitudes

Signal.N    = A(:,1).*conj(A(:,1)) ;
Idler.N     = A(:,2).*conj(A(:,2)) ;
Pump.N      = A(:,3).*conj(A(:,3)) ;
Total.N     = Signal.N + Idler.N + Pump.N;

Find electric fields

Signal.E    = A(:,1) * sqrt( 2 * eta * const.hbar * Signal.Omega) .* exp(1i .* Signal.kVector .* z);
Idler.E     = A(:,1) * sqrt( 2 * eta * const.hbar * Idler.Omega ) .* exp(1i .* Idler.kVector .* z);
Pump.E     = A(:,1) * sqrt( 2 * eta * const.hbar * Pump.Omega ) .* exp(1i .* Pump.kVector .* z);

Find powers instead of amplitudes

Signal.P    = Signal.N .* Param.Area .* const.hbar .* Signal.Omega;
Idler.P     = Idler.N .* Param.Area .* const.hbar .* Idler.Omega;
Pump.P      = Pump.N .* Param.Area .* const.hbar .* Pump.Omega;
Total.P     = Signal.P + Idler.P + Pump.P;

Z axis

Axis.zmili  = z*1e3;

output0 = Axis.zmili;
output1 = Signal.P;
output2 = Idler.P;
output3 = Pump.P;
output4 = Total.P;
output5 = Signal.N;
output6 = Idler.N;
output7 = Pump.N;
output8 = Total.N;

end

Published with MATLAB® R2013a