Fourier optics
HW 3

define functions for normalized FFTs in the t and ω domains

FFT[fLst_,npts_,dt_]:=(Sqrt[npts (dt^2)])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
InverseFFT[FLst_,npts_,dt_]:=(Sqrt[(1/(npts (dt^2)))])RotateLeft[InverseFourier[RotateRight[FLst,npts/2]],npts/2];

Problem 1

Using the FFT and convolve.nb program as a guide, create a pulse function, put it on a grid, and calculate and plot the final pulse shape for pure 2nd, 3rd and 4th order spectral phase. Do this for a Gaussian and a rectangular pulses of FWHM of 10fs and 100fs and a center wavelength of 800nm. 
Symbolically: start with (Subscript[E, in])(t), transform to (Subscript[E, in])(ω), multiply by Exp[ i ϕ(ω) ], then inverse FT to (Subscript[E, out])(t), then plot ((Abs[ (Subscript[E, out])(t) ])^2)

The Taylor expansion of spectral phase is defined as: 
ϕ(ω)=(Subscript[ϕ, 0])+(ω-(Subscript[ω, 0]))(Subscript[ϕ, 1])+(1/2)(((ω-(Subscript[ω, 0])))^2)(Subscript[ϕ, 2])+(1/(3!))(((ω-(Subscript[ω, 0])))^3)(Subscript[ϕ, 3])+(1/(4!))(((ω-(Subscript[ω, 0])))^4)(Subscript[ϕ, 4])+...,
so specifically, you’ll be plotting the output pulse shape for different values of (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4]), with (Subscript[ϕ, 0])=(Subscript[ϕ, 1])=0.
It is conventional to represent the (Subscript[ϕ, n])'s in units of (fs^n). 

So you should make 6 plots for each of the two pulse shapes: 
τ = 10fs: (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4])
τ = 100fs: (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4])
When plotting the temporal pulse shapes, determine the amount of each type of spectral phase that approximately triples the duration of the pulse. For each plot, make it clear what value of phase you are using, and plot the x-axis in units of fs. 
To do this, set up a grid for the time points of your array, and do a 
ListPlot[ Transpose[ {timeLst / fs, data} }, Joined  -> True], 
potentially using PlotRange to zoom in to get a better display.  
Keep an eye on whether the pulse broadens so much that the edges get close to the boundaries of the grid. If it does, you can increase the number of points or increase the dt spacing.

problem 2:

A signal f(t)=(sinc^2)((Subscript[ω, 1]) t), where (Subscript[ω, 1])=100, is sampled by a comb function, to produce the list (Subscript[f, s])(t)=f(t)comb((Subscript[ω, s])t), where (Subscript[ω, s]) is the sampling frequency (in rad/s). The sampled function is then passed through an LSI system with transfer function H(ω)=rect(ω/(Subscript[ω, s])) to produce an output g(t). 

For the following questions, I want you to explore the effect of different sampling rates.

part a:

Find the analytic transform, F(ω)=FT{f(t)} and make plots of f(t) and F(ω). Use the internally defined function Sinc[ ].

part b:

Find the minimum sampling rate (Subscript[ω, s])=(Subscript[ω, Ny]) (Nyquist rate) and the corresponding sampling interval (Subscript[t, Ny])=2π/(Subscript[ω, Ny]) that will permit an exact recovery of f(t) from (Subscript[f, s])(t); i.e. such that g(t)=f(t) within a multiplicative constant.

part c:

Generate a list (Subscript[f, s])(t) (using the Table[ ] function in Mathematica) with a sampling period (Subscript[t, s])=2π/(Subscript[ω, Ny]) for 32 points.  

Use the format of ListPlot[ ] shown below to plot (Subscript[f, s])(t). These options for ListPlot[ ] allow you to see dots for the points as well as connecting the lines.

ListPlot[fInLst,Mesh->All,Joined->True,MeshStyle->PointSize[Medium],PlotRange->All]

Take the FFT (using the FFT[ ] function defined in fft and convolve.nb), and plot the spectrum. 

Repeat the above for 64 points with the same sampling interval (Subscript[t, s])=2π/(Subscript[ω, Ny]). Comment on what the effect is in both domains when you increase npts while maintaining dt.

part d:

As we showed in the class notes, the analytic calculation of the FT of the sampled function results in the convolution of the spectrum with another comb function. In terms of equations: FT{ (Subscript[f, s])(t) } = FT{ f(t) comb[t/(Subscript[t, s])] } = (Subscript[F, s])(ω)=(1/(2π))F(ω)⊗[(Subscript[ω, s])comb((ω/(Subscript[ω, s])))]. 
Aside: to obtain the same amplitude scaling for this result, so that the peak value of (Subscript[F, s])(ω) equals the peak value of F(ω), I find that we need to divide by dω=((2π)/(N dt))=((Subscript[ω, s])/N), so that
(Subscript[F, s])(ω)=(1/(2π dω))F(ω)⊗[(Subscript[ω, s])comb((ω/(Subscript[ω, s])))]=(N/(2π))F(ω)⊗[comb((ω/(Subscript[ω, s])))]

Make a plot of this analytic (Subscript[F, s])[ω] over a range from -2(Subscript[ω, s]) to 2(Subscript[ω, s]). To represent the comb[ω/(Subscript[ω, s])] function, you can use the following definition,

comb[ω_,ωs_]:=Sum[DiracDelta[ω-n ωs],{n,-3,3}];

then you can use the Mathematica function Convolve[ ] to do the analytic convolution. 

How does the result of the FFT numerical transform relate to the analytic result?

part e:

Repeat the plots made in parts c and d above, but for under-sampling (Subscript[t, s])=(3/2)((2π)/(Subscript[ω, Ny])),  and over-sampling, (Subscript[t, s])=(1/2)((2π)/(Subscript[ω, Ny])).

Self-phase modulation and spectral broadening

Review class notes on the dispersion length and the nonlinear refractive index and self-phase modulation.

constants and unit definitions: internally, all numbers in SI

ps=(10^(-12));fs=(10^(-15));
mm=(10^(-3));cm=(10^(-2));cm3=(10^6);um=(10^(-6));nm=(10^(-9));angs=(10^(-10));pm=(10^(-12));
mJ=(10^(-3));uJ=(10^(-6));
c=2.9979 10^8;                          (*speed of light,in m/s*)
eps0=8.8542 (10^(-12));            (*permittivity of free space*)
e=1.6 (10^(-19));(*Coul electron charge*)

refractive index for BK7: nbk7

nbk7[lambda_]:=√(2.2718929-(1.0108077 (lambda/um)^2)/10^2+1.0592509/(10^2 (lambda/um)^2)+2.0816965/(10^4 (lambda/um)^4)-7.6472538/(10^6 (lambda/um)^6)+4.9240991/(10^7 (lambda/um)^8));

refractive index for fused silica: nSiO2

as1=.6961663;
as2=.4079426;
as3=.8974794;
bs1=.0684043^2;
bs2=.1162414^2;
bs3=9.896161^2;
nSiO2[lambda_]:=√(1+(as1 (lambda/um)^2)/((lambda/um)^2-bs1)+(as2 (lambda/um)^2)/((lambda/um)^2-bs2)+(as3 (lambda/um)^2)/((lambda/um)^2-bs3))

Problem 3. Experiment design: dispersion length

First we need to do some back-of-the-envelope calculations to determine what experimental conditions we would need to observe spectral broadening from SPM.

a. For fused silica and BK7 glass, plot (Subscript[k, 2])(λ) over the range 400nm to 1800nm, both curves on the same plot. Identify the zero dispersion wavelength for each material.

b. Next, assuming that the laser pulse duration is (Subscript[τ, p])=50fs (1/(e^2) half-width), plot (Subscript[L, D])(λ) over the same wavelength range for both glasses on the same plot. Near the zero dispersion point, the dispersion length goes to infinity, and third-order dispersion becomes the next-most important order of dispersion.

c. Using our square waveguide approximation for optical fiber dispersion, estimate the waveguide dimension that will push the zero dispersion wavelength to 800nm.

Problem 4. Experiment design: nonlinear length

Now we need to estimate the spot size and pulse energy to get sufficient intensity to see a nonlinear effect. We want the beam to be approximately collimated throughout the length of the sample.

a. Calculate the minimum focal spot size ((Subscript[w, 0]),intensity 1/(e^2) radius), so that the confocal parameter is equal to the window thickness. (The confocal parameter is 2(Subscript[z, R]), where (Subscript[z, R])=((π (Subscript[n, 0])((Subscript[w, 0])^2))/λ) is the Rayleigh range).

b. Given this minimum spot radius, we can estimate the intensity as I=(E/(τ π ((Subscript[w, 0])^2))). For future use, write a function that calculates corresponding intensity as a function of the pulse duration and energy.

c. Now calculate the approximate pulse energy such that the B-integral is equal to one.  (B = ((Subscript[ω, 0])/c)(Subscript[n, 2])I L). At B = 1, the nonlinear effect on the spectrum begins to be noticeable. In fused silica, the nonlinear coefficient has been measured to be (Subscript[n, 2])=2.5(x10^(-16))(cm^2)/W. Also then calculate the intensity for the focused beam under these conditions in W/(cm^2).

Problem 5. instantaneous frequency shift

a. The nonlinear phase shift is (Subscript[ϕ, NL])(t)=((Subscript[ω, 0])/c)(Subscript[n, 2])I(t)L, where I(t)=(Subscript[I, 0])exp[-2(t^2)/((Subscript[τ, 0])^2)]. 

Write a function for the nonlinear phase shift as a function of time, pulse duration and pulse energy. (We will hold the spot size and sample length constant. ) Plot ϕNL[t] for the energy calculated above. If your pulse energy is correct, the peak value should be 1.

b. Now calculate the instantaneous frequency shift: 
(Subscript[Δω, inst])(t)=-(∂_t) (Subscript[ϕ, NL])(t)
Calculate an expression for (Subscript[Δω, inst])(t), then plot the relative frequency shift: (((Subscript[Δω, inst])(t))/(Δω)), where Δω=2/(Subscript[τ, 0]) is the bandwidth of the incident pulse. Again, if your pulse energy is correct, the peak value should be of the order 1.

Problem 6. Numerical calculation of SPM spectrum

To see what the self-phase modulated spectrum looks like, we need to calculate it numerically. Set up time and frequency grids, sample the self-phase modulated pulse in the time domain, then calculate the pulse spectral intensity. Feel free to use the FFT functions and code for setting up the grids from my examples.  

To calculate a list for the nonlinear phase at a given pulse energy, you can use the function you defined above with the list of time points ϕNL[tLst, τ, E ]. The peak B-integral is just the maximum value of this list. 
 
Make plots of the spectrum vs frequency for peak B-integrals of 1, 5, and 20.  On each plot, also show the original spectrum for comparison. Normalize the plots by dividing the array you are plotting by its maximum value.

Problem 7. simple pulse compression

a. The self-phase modulated pulse is not any different in duration than the input pulse, but the if the right spectral phase is applied, the broader bandwidth can be compressed to get a shorter pulse. 

The best compression would be a pulse with spectral phase that is constant.  To see what the Fourier transform limit is for the broadened spectrum, calculate the inverse FFT of the absolute value of the spectrum for the case where the B-integral = 20.

b. Next calculate a list Exp[i(1/2)(((ωLst-ω0))^2)ϕ2Tst] and use this to try to compress the pulse. Calculate the inverse FFT to get the pulse in the time domain. Plot the intensity in the time domain for the compressed pulse along with the original pulse shape.  Vary the value of the (Subscript[ϕ, 2]) used for compression so that you get a shorter pulse than the initial value. 

Normalize the plots to better compare the original with the compressed pulse.