Computational Hyperspectral Terahertz Holography

Kokoteksti

(1)MAKSIM KULYA. Computational Hyperspectral Terahertz Holography. Tampereen yliopiston väitöskirjat 387.

(2)

(3) Tampere University Dissertations 387. MAKSIM KULYA. Computational Hyperspectral Terahertz Holography. ACADEMIC DISSERTATION To be presented, with the permission of the Faculty of Information Technology and Communication Sciences of Tampere University, for online public discussion, on 5 March 2021, at 12 o’clock..

(4) ACADEMIC DISSERTATION Tampere University, Faculty of Information Technology and Communication Sciences Finland Responsible supervisor and Custos. Prof. Karen Eguiazarian Tampere University. Finland. Supervisors. Prof. Vladimir Katkovnik Tampere University Finland. Dr. Nikolay Petrov ITMO University Russia. Pre-examiners. Dr. Lorenzo Valzania École Normale Supérieure France. Dr. Lu Rong Beijing University of Technology China. Opponents. Dr. Konstantinos Falaggis University of North California at Sharlotte USA. Dr. Lu Rong Beijing University of Technology China. The originality of this thesis has been checked using the Turnitin OriginalityCheck service.. Copyright ©2021 author. Cover design: Roihu Inc.. ISBN 978-952-03-1883-3 (print) ISBN 978-952-03-1884-0 (pdf) ISSN 2489-9860 (print) ISSN 2490-0028 (pdf) http://urn.fi/URN:ISBN:978-952-03-1884-0. PunaMusta Oy – Yliopistopaino Joensuu 2021.

(5) To Anna of the Netherlands "Knowledge is real knowledge only when it is acquired by the efforts of your intellect, not by memory". Leo Tolstoy. iii.

(6) iv.

(7) PREFACE/ACKNOWLEDGEMENTS. I would like to express my gratitude to my advisor Prof. Karen Eguiazarian for the continuous support of my Ph.D study and related research, for his patience, motivation, and immense knowledge. My sincere thanks also goes to my co-advisor Prof. Vladimir Katkovnik who directed me in the various scientific problems. I would like to thank also my chef Nikolay Petrov for multi-year leadership of my scientific activities.. v.

(8) vi.

(9) ABSTRACT. The progress in terahertz (0.1-10 THz) technologies has opened access to electromagnetic waves that are extremely short in the number of oscillations. Moreover, the result of THz technologies is not the envelope or radiation power but the spatial distribution of the electric field itself. This allows to access the amplitude-phase information for all spectral components. It makes possible the development of the holographic approach which involves measurements in a wide collimated beam. Moreover, in the practical implementation of measurements, it is possible to decrease the number of optical elements for THz range, since the field is recorded by a scanning diaphragm or a wide-aperture crystal. In the thesis the method of digital hyperspectral THz holography was developed and improved, including the description of the physical, mathematical and computational aspects. This method allows to study the spatio-temporal and spatio-spectral dynamics of electric wavefields in linear isotropic homogeneous non-magnetic medium. The developed THz holography allows also to analyze the arbitrary vector beams including the behavior of transverse and longitudinal components of the field and polarization state. The practical application of THz holography was demonstrated on example of data transmission using the interference of pairs of THz vortex beams. The low signal-to-noise ration and high level of noise problems were solved by using the block-matching algorithms adapted for three dimensional complex THz fields, thus increasing the resolution and image quality of THz hyperspectral holograms. The developed approach was also expanded to the near infrared range corresponding to femtosecond laser pulses. Here the features of correlative registration of broadband laser fields of femtosecond duration was realized using the hyperspectral holography.. vii.

(10) viii.

(11) CONTENTS. 1. 2. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 37. 1.1. Thesis contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 39. 1.2. Thesis and published results . . . . . . . . . . . . . . . . . . . . . . . . . . .. 39. Terahertz hyperspectral holography . . . . . . . . . . . . . . . . . . . . . . . . .. 41. 2.1. Theoretical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 41. 2.1.1. General case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 41. 2.1.2. Specific cases of Fourier transforms . . . . . . . . . . . . . . . . .. 44. Numerical calculation peculiarities . . . . . . . . . . . . . . . . . . . . . .. 45. 2.2.1. Angular spectrum method . . . . . . . . . . . . . . . . . . . . . . .. 45. 2.2.1.1. Discrete Fourier transform peculiarities . . . . . . .. 48. 2.2.2. Convolution method . . . . . . . . . . . . . . . . . . . . . . . . . .. 49. 2.2.3. The choice of the calculation method . . . . . . . . . . . . . . . .. 51. Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 52. 2.3.1. Pulsed THz wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 52. 2.3.2. THz wavefront . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 54. 2.3.3. Results of THz vectorial wavefront propagation . . . . . . . .. 55. 2.4. Layouts for THz hyperspectral holography . . . . . . . . . . . . . . . . .. 60. 2.5. Data processing and visualization . . . . . . . . . . . . . . . . . . . . . . .. 61. 2.6. Peculiarities of the dynamics of THz wave field propagation . . . . . .. 63. 2.6.1. Diffraction effects . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 64. 2.6.2. Dispersion effects . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 66. Application of THz hyperspectral holography for beams analysis . . . . . .. 75. 2.2. 2.3. 3. ix.

(12) 3.1. 3.2. 4. Gauss-Bessel beams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 75. 3.1.1. Bessel-Gaussian beam’s paraxial diffraction . . . . . . . . . . . .. 76. 3.1.2. Peculiarities of Bessel-Gaussian beam’s non-paraxial diffraction 81. Vortex beams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 83. 3.2.1. Vortex beam’s diffraction . . . . . . . . . . . . . . . . . . . . . . .. 84. 3.2.2. Application of THz vortex beams for wireless data transfer .. 86. 3.2.3. Coding of information based on the interference of pairs of THz beams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 91. Hyperspectral data denoising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 97. 4.1. Photo conductive antenna system . . . . . . . . . . . . . . . . . . . . . . .. 97. 4.1.1. Mathematical model of denoising . . . . . . . . . . . . . . . . . .. 98. 4.1.2. Numerical simulation scheme . . . . . . . . . . . . . . . . . . . . 100. 4.1.3. Experimental and simulation results . . . . . . . . . . . . . . . . 102. 4.2. 4.1.3.1. Experimental setup . . . . . . . . . . . . . . . . . . . . . 102. 4.1.3.2. Spectral domain denoising . . . . . . . . . . . . . . . . 104. 4.1.3.3. Relief reconstruction . . . . . . . . . . . . . . . . . . . . 107. 4.1.3.4. Time-domain denoising . . . . . . . . . . . . . . . . . . 108. Balance detection system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.2.1. Mathematical model of denoising . . . . . . . . . . . . . . . . . . 111 4.2.1.1. 5. Complex domain cube filter . . . . . . . . . . . . . . . 114. 4.2.2. Numerical simulation scheme . . . . . . . . . . . . . . . . . . . . 115. 4.2.3. Experimental and simulation results . . . . . . . . . . . . . . . . 117 4.2.3.1. Experimental setup . . . . . . . . . . . . . . . . . . . . . 117. 4.2.3.2. Time-domain denoising . . . . . . . . . . . . . . . . . . 118. 4.2.3.3. Spectral-domain denoising . . . . . . . . . . . . . . . . 119. Peculiarities of correlation measurements of pulsed hyperspectral optical fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127. 5.2. Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 x.

(13) 5.3 6. Mathematical model and results . . . . . . . . . . . . . . . . . . . . . . . . 130. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.1. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138. References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141. List of Programs and Algorithms matlab/main.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii matlab/E_initial.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv matlab/read_masks.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi matlab/GaussBeam.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi matlab/phasefromrelief.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii matlab/PROP.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii matlab/PROPz.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii matlab/PROPz_polar.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix matlab/AS.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. xx. matlab/ASz.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi matlab/ASz_polar.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii matlab/C.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii matlab/Cz.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxv matlab/Cz_polar.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvi matlab/Gshift.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii matlab/back_E.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii matlab/visual_Axnu_y.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxix matlab/visual_PHxnu_y.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxx matlab//visual_Ext_y.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxxi Here the listing of working codes with necessary comments inside is presented. Hope that these realization will help all researches to have a progress in the field of xi.

(14) wavefront propagation, holography, hyperspectral imaging and other areas of wave optics. % Version "v .5.2.1 Vanilla " includes THz pulse initialisation , wavefront modulation, its propagation for transverse and longitudinal components and visualization for temporal and spectral patterns . %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % main.m file which execute wavefront propagation opengl (’ save ’, ’ software ’) path = ’C:\SCIENCE\2020\testt’;. % create main path. if ~exist(path, ’ dir ’) mkdir(path); end pathAM = strcat(path,’\\white128.bmp’); % indicate the path for amplitude mask. Use white square 8 bit pathPH = strcat(path ,’\\ white128.bmp’); % indicate the path for phase mask. Use white square 8 bit c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/ps %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Initialization of the pulse N=1024; % number of points in time window T=100; % time window in ps tau=0.4; % pulse duration [E,G,t,nu,ind] = E_initial (N,T,tau,path) ; % initialize. initial E(t ) from. parameters %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Initialization of G(x,y,nu) in object plane [Mask_Am, Mask_Ph] = read_masks(pathAM,pathPH); % read and normalize masks. xii.

(15) x_size=50; % screen size in mm y_size=50; % screen size in mm rhogauss=5; % Gaussian parameter rho in mm [ Gauss ] = GaussBeam( Mask_Am,x_size,y_size,rhogauss ); % create Gaussian beam Mask_Am=Mask_Am.Gauss; % multiply initial Mask_Am to the gaussian beam %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % wavefront modulation n=1.46; % refractive index of the object height=0; % height of the object in mm. Use 0 for "no phase modulation" [G_object] = phasefromrelief (Mask_Am,Mask_Ph,n,height,nu,G); % zombilizer %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % remove extra frequencies from the G(x,y,nu) with zero amplitude nucrop=200; % crop limit G_object=G_object (:,:, size (G_object,3)/2+1:size (G_object,3)/2+nucrop); % crop G(x,y,nu) nu=nu((size(nu,2) /2)+1:size (nu,2)/2+nucrop); % crop nu array %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % wavefront propagation z=10; % propagation distance nu0=cz/x_size /( x_size/ size (Mask_Am,1)); % set boundary frequency nu0 which determines AS or C method [ GG ] = PROP( G_object,nu,nu0,x_size,y_size,z ) ; % propagation of transverse component [ GGz ] = PROPz( G_object,nu,nu0,x_size,y_size,z ) ;GGz (:,:,1) =0; % propagation of longitudinal component [GG_shift] = Gshift(GG,nu,z,Mask_Am); [EE] = back_E(GG_shift,N); % transition back to E(t ) form G(nu) ( for transverse component). xiii.

(16) [GGz_shift] = Gshift(GGz,nu,z,Mask_Am); % phase shift of propagated waveront to visualize it at the same temporal scale [EEz] = back_E(GGz_shift,N); % transition back to E(t ) form G(nu) ( for longitudinal component) %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % visualizers % for transverse component y_coord=64;nu_min=1;nu_max=200;prfx=0;type=’trans’;caxis1=0;caxis2=15; xlim1=0;xlim2=2;ylim1=−20;ylim2=20; visual_Axnu_y(GG,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max,prfx,path, type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0); visual_PHxnu_y(GG,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max,prfx,path, type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0); y_coord=64;t_min=1;t_max=N;prfx=0;type=’trans’;caxis1=−1;caxis2=1;xlim1 =−5;xlim2=5;ylim1=−20;ylim2=20; visual_Ext_y (EE,t,z,Mask_Am,x_size,y_coord,t_min,t_max,prfx,path,type, caxis1, caxis2 ,xlim1,xlim2,ylim1,ylim2) ; % for longitudinal component y_coord=64;nu_min=1;nu_max=200;prfx=0;type=’long’;caxis1=0;caxis2=15; xlim1=0;xlim2=2;ylim1=−20;ylim2=20; visual_Axnu_y(GGz,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max,prfx,path, type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0); visual_PHxnu_y(GGz,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max,prfx,path, type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0); y_coord=64;t_min=1;t_max=N;prfx=0;type=’long’;caxis1=−1;caxis2=1;xlim1 =−5;xlim2=5;ylim1=−20;ylim2=20; visual_Ext_y (EEz,t,z,Mask_Am,x_size,y_coord,t_min,t_max,prfx,path,type, caxis1, caxis2 ,xlim1,xlim2,ylim1,ylim2) ; % E_initial .m: initialize. initial E(t ) and G(nu) from parameters. function [ E,G,t,nu,ind ] = E_initial (N,T,tau,path). xiv.

(17) if ~exist( strcat (path ,’\\ E_initial ’) ,’ dir ’) % create path E_initial to save figures mkdir( strcat (path ,’\\ E_initial ’) ) ; end t=linspace(−T/2,T/2−T/N,N); % create time grid nu=linspace(−N/2/T,N/2/T−1/T,N); % create frequency grid E=exp(−(t/tau).^2) . t/tau ; % initialize E(t ) E=E/(max(E)); % normalization [−1;1] f1=figure (’ visible ’,’ off ’) ; % save figure E(t ) plot ( t ,E ,’−’,’ LineWidth’,2) ; xlabel (’ t , ps ’) ; ylabel (’ E/Eo’); set (gca ,’ YLim’,[−1.1 1.1]) ; set (gca ,’ XLim’,[−5 5]); set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f1 , strcat (path ,’\\ E_initial ’,’\\ Et.png’) ) ; close ( gcf ) ; G=fftshift ( fft ( fftshift (E)) ) ; % FFT solution of E(t ) to get G(nu) [G0,ind]=max(abs(G)); % find maximum of G(nu) f2=figure (’ visible ’,’ off ’) ; % save figure amplitude of G(nu) plot (nu(N/2+1:N/2+400),abs(G(N/2+1:N/2+400)/G0),’−’,’LineWidth’,2); xlabel(’\nu, THz’);ylabel(’|G/Go|’); set (gca ,’ XLim’,[0 3.5]) ;%set (gca ,’ YLim’,[0 1.1]) ; set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f2 , strcat (path ,’\\ E_initial ’,’\\ abs(G(nu)).png’) ) ;. xv.

(18) close ( gcf ) ; f3=figure (’ visible ’,’ off ’) ; % save figure phase of G(nu) plot (nu(N/2+2:N/2+400),angle(G(N/2+2:N/2+400)),’−’,’LineWidth’,2);xlabel (’\nu, THz’);ylabel(’phi,rad .’) ; set (gca ,’ XLim’,[0 3.5]) ; set (gca ,’ YLim’,[−2 −1]); set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f3 , strcat (path ,’\\ E_initial ’,’\\ angle (G(nu)).png’) ) ; close ( gcf ) ; end % read_masks.m: read and normalize masks function [ Mask_Am, Mask_Ph ] = read_masks( Am_file, Ph_file ) Am=imread(Am_file); Mask_Am=im2double(Am); Mask_Am=Mask_Am/max(max(Mask_Am));max(Mask_Am(:)) Ph=imread(Ph_file); Mask_Ph = im2double(Ph); Mask_Ph=Mask_Ph/max(max(Mask_Ph)); end % GaussBeam.m: create Gaussian beam function [ Gauss ] = GaussBeam( Mask_Am,x_size,y_size,rho ) x=linspace(−x_size/2, x_size/2−x_size/size (Mask_Am,1),size(Mask_Am,1)); y=linspace(−y_size/2, y_size/2−y_size/size(Mask_Am,2),size(Mask_Am,2)); for i=1:size (Mask_Am,1) for j=1:size (Mask_Am,2) xvi.

(19) Gauss(i , j )=exp(−(x(i).^2+y(j ) .^2) /rho^2); end end end % phasefromrelief .m: set amplitude and phase modulation according Mask_Am and Mask_Ph function [ G_object,am_object, phi_object ] = phasefromrelief (Mask_Am, Mask_Ph,n,height,nu,G) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in m/s phi_object=zeros( size (Mask_Ph,1),size (Mask_Ph,2),size (nu,2) ) ; % preallocation am_object=zeros(size(Mask_Am,1),size(Mask_Am,2),size(nu,2)); % preallocation G_object=zeros(size(Mask_Am,1),size(Mask_Am,2),size(nu,2)); % preallocation for k=1:size(nu,2) am_object (:,:, k)= Mask_Am.abs(G(k)); % amplitude modulation phi_object (:,:, k)= Mask_Ph .(2 pi ( n−1) height nu(k)/c)+angle(G(k)); % phase modulation G_object (:,:, k)=am_object (:,:, k) . exp(1 j phi_object (:,:, k) ) ; % set initial wavefront as 3D array G(x,y,nu) end end % PROP.m: wavefront propagation of transverse component function [ GG ] = PROP( G_object,nu,nu0,x_size,y_size,z ) if z==0 % for z=0 use initial waveront GG=G_object; else xvii.

(20) aa=find(nu<nu0);nu0_ind=max(aa)+1; % find nu0 in array nu and devide it by 2 subarrays if mod(nu0_ind,2) == 1 nu0_ind_odd=nu0_ind+1; else nu0_ind_odd=nu0_ind; end if nu0_ind_odd >= size(G_object,3) [GG] = C(G_object,nu,x_size,y_size,z) ; else G_object_C=G_object(:,:,1:nu0_ind_odd); % set subaaray for C method [GG_C] = C(G_object_C,nu(1:nu0_ind_odd),x_size,y_size,z); % C method G_object_AS=G_object(:,:,nu0_ind_odd+1:size(G_object,3)); % set subaaray for AS method [GG_AS] = AS(G_object_AS,nu(nu0_ind_odd+1:size(G_object,3)),x_size, y_size,z); % AS method GG=cat(3,GG_C,GG_AS); % concatenate 2 subarrays end end end % PROPz.m: wavefront propagation of longitudianal component function [ GGz ] = PROPz( G_object,nu,nu0,x_size,y_size,z ) if z==0 % for z=0 use only AS method for longitudinal component GGz = ASz(G_object,nu,x_size,y_size,z) ; else aa=find(nu<nu0);nu0_ind=max(aa)+1; % find nu0 in array nu and devide it by 2 subarrays if mod(nu0_ind,2) == 1 nu0_ind_odd=nu0_ind+1; xviii.

(21) else nu0_ind_odd=nu0_ind; end if nu0_ind_odd >= size(G_object,3) [GGz] = Cz(G_object,nu,x_size,y_size,z) ; else G_object_C=G_object(:,:,1:nu0_ind_odd); % set subarray for Cz method [GG_Cz] = Cz(G_object_C,nu(1:nu0_ind_odd),x_size,y_size,z); % Cz method G_object_AS=G_object(:,:,nu0_ind_odd+1:size(G_object,3)); % set subarray for ASz method [GG_ASz] = ASz(G_object_AS,nu(nu0_ind_odd+1:size(G_object,3)),x_size, y_size,z); % ASz method GGz=cat(3,GG_Cz,GG_ASz); % concatenate 2 subarrays end end end % PROPz_polar.m: wavefront propagation of longitudinal component Gx when 2 transverse components Gx % and Gy exist function [ GGz ] = PROPz_polar( G_object_x,G_object_y,nu,nu0,x_size,y_size ,z ) if z==0 % for z=0 use only ASz_polar method for longitudinal component GGz = ASz_polar(G_object_x,G_object_y,nu,x_size,y_size,z); else aa=find(nu<nu0);nu0_ind=max(aa)+1; % find nu0 in array nu and devide it by 2 subarrays if mod(nu0_ind,2) == 1 nu0_ind_odd=nu0_ind+1; else xix.

(22) nu0_ind_odd=nu0_ind; end if nu0_ind_odd >= size(G_object_x,3) [GGz] = Cz(G_object_x,nu,x_size,y_size,z);%GG_C(:,:,1)=0; else G_object_x_C=G_object_x(:,:,1:nu0_ind_odd); % set subarray for Cz_polar method for Gx component G_object_y_C=G_object_y(:,:,1:nu0_ind_odd); % set subarray for Cz_polar method for Gy component [GG_Cz] = Cz_polar(G_object_x_C,G_object_y_C,nu(1:nu0_ind_odd), x_size,y_size,z); % Cz_polar method G_object_x_AS=G_object_x(:,:,nu0_ind_odd+1:size(G_object_x,3)); % set subarray for ASz_polar method for Gx component G_object_y_AS=G_object_y(:,:,nu0_ind_odd+1:size(G_object_y,3)); % set subarray for ASz_polar method for Gy component [GG_ASz] = ASz_polar(G_object_x_AS,G_object_y_AS,nu(nu0_ind_odd +1:size(G_object_x,3)),x_size,y_size,z); % ASz_polar method GGz=cat(3,GG_Cz,GG_ASz); % concatenate 2 subarrays end end end % AS.m: execute AS method for transverse component function [ GG ] = AS( G_object,nu,x_size, y_size ,z ) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/fs G_plwave=fftshift( fft2 ( fftshift (G_object)) ) ; % plane wave decomposition H=zeros(size(G_object,1) , size (G_object,2) , size (G_object,3) ) ; % preallocation xx.

(23) U=zeros(size(G_object,1) , size (G_object,2) , size (G_object,3) ) ; % preallocation Nx=size(G_object,1); Ny=size(G_object,2); fx=linspace(−Nx/2/x_size,Nx/2/x_size−1/x_size,Nx); % set fx grid fy=linspace(−Ny/2/y_size,Ny/2/y_size−1/y_size,Ny); % set fy grid for k=1:size(G_object,3) for i=1:size (G_object,1) for j=1:size (G_object,2) U(i, j ,k)=1−((c fx ( i )/nu(k)) ^2)−((c fy ( j )/nu(k)) ^2) ; H(i, j ,k)=exp(−1jz (2 pi nu(k)/c) sqrt (U(i, j ,k) ) ) ; % transfer function end end end H(U<0)=0; % condition for non negative radicant U G_plane_H=G_plwave.H; % multiplication of plane waves to transfer function GG=ifftshift( ifft2 ( ifftshift (G_plane_H))); % plane wave composition end % ASz.m: execute AS method for longitudinal component function [ GGz ] = ASz( G_object,nu,x_size, y_size ,z ) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in m/s G_plwave=fftshift( fft2 ( fftshift (G_object)) ) ; % plane wave decomposition Hz=zeros(size(G_object,1) , size (G_object,2) , size (G_object,3) ) ; % preallocation U=zeros(size(G_object,1) , size (G_object,2) , size (G_object,3) ) ; % preallocation xxi.

(24) Nx=size(G_object,1); Ny=size(G_object,2); fx=linspace(−Nx/2/x_size,Nx/2/x_size−1/x_size,Nx); % set fx grid fy=linspace(−Ny/2/y_size,Ny/2/y_size−1/y_size,Ny); % set fy grid for k=1:size(G_object,3) for i=1:size (G_object,1) for j=1:size (G_object,2) U(i, j ,k)=1−((c fx ( i )/nu(k)) ^2)−((c fy ( j )/nu(k)) ^2) ; Hz(i, j ,k)= exp(−1j z (2 pi nu(k)/c) sqrt (U(i, j ,k) ) ) fx ( i ) c/(nu(k) sqrt (U(i , j ,k) ) ) ; % transfer function end end end Hz(U<0)=0; % condition for non negative radicant U G_plwave_Hz=G_plwave.Hz; % multiplication of plane waves to transfer function GGz=ifftshift( ifft2 ( ifftshift (G_plwave_Hz))); % plane wave composition end % ASz_polar.m: execute ASz_polar method for longitudinal component Gx when 2 transverse components Gx % and Gy exist function [ GGz ] = ASz_polar( G_object_x,G_object_y,nu,x_size,y_size ,z ) c=(1E−9)physconst(’LightSpeed ’) ;%speed of light in m/s G_plwave_x=fftshift( fft2 ( fftshift (G_object_x))) ; % plane wave decomposition of Gx component G_plwave_y=fftshift( fft2 ( fftshift (G_object_y))) ; % plane wave decomposition of Gy component. xxii.

(25) Hxz=zeros(size(G_object_x,1), size (G_object_x,2) , size (G_object_x,3)) ; % preallocation Hyz=zeros(size(G_object_x,1), size (G_object_x,2) , size (G_object_x,3)) ; % preallocation U=zeros(size(G_object_x,1) , size (G_object_x,2) , size (G_object_x,3)) ; % preallocation Nx=size(G_object_x,1); Ny=size(G_object_x,2); fx=linspace(−Nx/2/x_size,Nx/2/x_size−1/x_size,Nx); % set fx grid fy=linspace(−Ny/2/y_size,Ny/2/y_size−1/y_size,Ny); % set fy grid for k=1:size(G_object_x,3) for i=1:size (G_object_x,1) for j=1:size (G_object_x,2) U(i, j ,k)=1−((c fx ( i )/nu(k)) ^2)−((c fy ( j )/nu(k)) ^2) ; Hxz(i, j ,k)= exp(−1j z (2 pi nu(k)/c) sqrt (U(i, j ,k) ) ) fx ( i ) c/(nu(k) sqrt (U( i, j ,k) ) ) ; % transfer function Hxz Hyz(i, j ,k)= exp(−1j z (2 pi nu(k)/c) sqrt (U(i, j ,k) ) ) fy ( j ) c/(nu(k) sqrt (U( i, j ,k) ) ) ; % transfer function Hyz end end end Hxz(U<0)=0; % condition for non negative radicant U Hyz(U<0)=0; % condition for non negative radicant U G_plwave_Hz=G_plwave_x.Hxz+G_plwave_y.Hyz; % multiplication of plane waves to transfer function GGz=ifftshift( ifft2 ( ifftshift (G_plwave_Hz))); % plane wave composition end % C.m: execute C method for transverse component function [ GG ] = C( G_object,nu,x_size, y_size ,z ). xxiii.

(26) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/fs G_object_zp=padarray(G_object,[size(G_object,1)/2 size (G_object,2)/2 0],0+0 j ,’ both’) ; % twice zeroppading x_size=x_size 2; y_size=y_size 2; % twice increasing of grid size G_plwave=fftshift( fft2 ( fftshift (G_object_zp))) ; % plane wave decomposition Nx=size(G_object_zp,1); Ny=size(G_object_zp,2); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid y=linspace(−y_size/2, y_size/2−y_size/Ny,Ny); % set y grid h=zeros(size (G_object_zp,1), size (G_object_zp,2), size (G_object_zp,3)) ; % preallocation for k=1:size(G_object_zp,3) for i=1:size (G_object_zp,1) for j=1:size (G_object_zp,2) h( i , j ,k)=((x_size/ size (G_object_zp,1)) ^2) exp(−1j sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2) 2 pi nu(k)/c) z/(−1j c/nu(k) sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2).^2) ; end end end HH=fftshift(fft2 ( fftshift (h) ) ) ; % transfer function G_plwave_HH=G_plwave.HH; % multiplication of plane waves to transfer function GG=ifftshift( ifft2 ( ifftshift (G_plwave_HH))); % plane wave composition GG=GG(size(G_object,1)/2+1:size(G_object,1)/2+1+size(G_object,1)−1,size( G_object,2)/2+1:size(G_object,2)/2+1+size(G_object,2)−1,:) ; % crop propagated array to the initial size end. xxiv.

(27) % Cz.m: execute C method for longitudinal component function [ GGz ] = Cz( G_object,nu,x_size, y_size ,z ) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/fs G_object_zp=padarray(G_object,[size(G_object,1)/2 size (G_object,2)/2 0],0,’ both’) ; % twice zeroppading x_size=x_size 2; y_size=y_size 2; % twice increasing of grid size G_plwave=fftshift( fft2 ( fftshift (G_object_zp))) ; % plane wave decomposition Nx=size(G_object_zp,1); Ny=size(G_object_zp,2); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid y=linspace(−y_size/2, y_size/2−y_size/Ny,Ny); % set y grid fx=linspace(−Nx/2/x_size,Nx/2/x_size−1/x_size,Nx); % set fx grid fy=linspace(−Ny/2/y_size,Ny/2/y_size−1/y_size,Ny); % set fx grid U=zeros(size(G_object_zp,1), size (G_object_zp,2), size (G_object_zp,3)) ; % preallocation h=zeros(size (G_object_zp,1), size (G_object_zp,2), size (G_object_zp,3)) ; % preallocation kxkz=zeros(size(G_object_zp,1), size (G_object_zp,2), size (G_object_zp,3)) ; % preallocation for k=1:size(G_object_zp,3) for i=1:size (G_object_zp,1) for j=1:size (G_object_zp,2) U(i, j ,k)=1−((c fx ( i )/nu(k)) ^2)−((c fy ( j )/nu(k)) ^2) ; h( i , j ,k)=((x_size/ size (G_object_zp,1)) ^2) exp(−1j sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2) 2 pi nu(k)/c) z/(−1j c/nu(k) sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2).^2) ; kxkz(i , j ,k)= fx( i ) c/(nu(k) sqrt (U(i, j ,k) ) ) ; end end. xxv.

(28) end kxkz(U<0)=0; % condition for non negative radicant U HHz=fftshift(fft2 ( fftshift (h) ) ) . kxkz; % transfer function G_plwave_HHz=G_plwave.HHz; % multiplication of plane waves to transfer function GGz=ifftshift( ifft2 ( ifftshift (G_plwave_HHz))); % plane wave composition GGz=GGz(size(G_object,1)/2+1:size(G_object,1)/2+1+size(G_object,1)−1,size( G_object,2)/2+1:size(G_object,2)/2+1+size(G_object,2)−1,:); % crop propagated array to the initial size end % Cz_polar.m: execute C method for longitudinal component Gx when 2 transverse components Gx % and Gy exist function [ GGz ] = Cz_polar( G_object_x,G_object_y,nu,x_size,y_size ,z ) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/fs G_object_x_zp=padarray(G_object_x,[size(G_object_x,1)/2 size (G_object_x,2) /2 0],0,’ both’) ; % twice zeroppading G_object_y_zp=padarray(G_object_y,[size(G_object_y,1)/2 size (G_object_y,2) /2 0],0,’ both’) ; % twice zeroppading x_size=x_size 2; y_size=y_size 2; % twice increasing of grid size G_plwave_x=fftshift( fft2 ( fftshift (G_object_x_zp))) ; % plane wave decomposition for Gx component G_plwave_y=fftshift( fft2 ( fftshift (G_object_y_zp))); % plane wave decomposition for Gy component Nx=size(G_object_x_zp,1); Ny=size(G_object_x_zp,2); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid y=linspace(−y_size/2, y_size/2−y_size/Ny,Ny); % set y grid xxvi.

(29) fx=linspace(−Nx/2/x_size,Nx/2/x_size−1/x_size,Nx); % set fx grid fy=linspace(−Ny/2/y_size,Ny/2/y_size−1/y_size,Ny); % set fy grid U=zeros(size(G_object_x_zp,1), size (G_object_x_zp,2), size (G_object_x_zp,3)) ; % preallocation h=zeros(size (G_object_x_zp,1), size (G_object_x_zp,2), size (G_object_x_zp,3)) ; % preallocation kxkz=zeros(size(G_object_x_zp,1), size (G_object_x_zp,2), size (G_object_x_zp ,3)) ; % preallocation kykz=zeros(size(G_object_x_zp,1), size (G_object_x_zp,2), size (G_object_x_zp ,3)) ; % preallocation for k=1:size(G_object_x_zp,3) for i=1:size (G_object_x_zp,1) for j=1:size (G_object_x_zp,2) U(i, j ,k)=1−((c fx ( i )/nu(k)) ^2)−((c fy ( j )/nu(k)) ^2) ; h( i , j ,k)=((x_size/ size (G_object_x_zp,1))^2) exp(−1j sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2) 2 pi nu(k)/c) z/(−1j c/nu(k) sqrt (( x( i ) .^2) +(y(j ) .^2) +z^2) .^2) ; kxkz(i , j ,k)= fx( i ) c/(nu(k) sqrt (U(i, j ,k) ) ) ; kykz(i , j ,k)= fy( j ) c/(nu(k) sqrt (U(i, j ,k) ) ) ; end end end kxkz(U<0)=0; % condition for non negative radicant U kykz(U<0)=0; % condition for non negative radicant U HHxz=fftshift(fft2 ( fftshift (h) ) ) . kxkz; % transfer function Hxz HHyz=fftshift(fft2( fftshift (h) ) ) . kykz; % transfer function Hyz G_plwave_HHz=G_plwave_x.HHxz+G_plwave_y.HHyz; % multiplication of plane waves to transfer functions GGz=ifftshift( ifft2 ( ifftshift (G_plwave_HHz))); % plane wave composition GGz=GGz(size(G_object_x,1)/2+1:size(G_object_x,1)/2+1+size(G_object_x,1). xxvii.

(30) −1,size(G_object_x,2)/2+1:size(G_object_x,2)/2+1+size(G_object_x,2)−1,:); % crop propagated array to the initial size end % Gshift.m: makes phase shift of propagated waveront to visualize it at the same temporal scale function [GGout] = Gshift(GGin,nu,z,Mask_Am) c=(1E−9)physconst(’LightSpeed ’) ; %speed of light in mm/fs GGout=zeros(size(GGin,1),size(GGin,2),size (GGin,3)); % preallocation for i=1:size (Mask_Am,1) for j=1:size (Mask_Am,2) for k=1:size(nu,2) GGout(i,j,k)=GGin(i,j,k) . exp(1 j 2 pi nu(k) z/c) ; end end end end % back_E.m: transition back to E(t ) form G(nu) function [EE] = back_E(GG,N) GG_ansig=padarray(GG,[0 0 N/2−size(GG,3)],0,’post’); % zeropadding according analytical signal theory GG_ansig2=padarray(GG_ansig,[0 0 N/2],0,’pre’); % zeropadding according analytical signal theory EE=zeros(size(GG_ansig2,1), size (GG_ansig2,2), size (GG_ansig2,3)); % preallocation for i=1:size (GG,1) for j=1:size (GG,2) EE(i, j ,:) =2 real ( fftshift ( ifft ( fftshift (GG_ansig2(i,j ,:) ) ) ) ) ; % get E(t ) xxviii.

(31) according analytical signal theory end end end % visual_Axnu_y.m: visualization of spatio −spectral amplitude of G(x,y= y_coord,nu) function visual_Axnu_y(GG,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max, prfx,path,type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0) pt=strcat (path ,’\\’, type ,’\\ A(x,nu)_y’) ; % create path A(x,nu)_y to save figures if ~exist(pt ,’ dir ’) mkdir(pt) ; end Nx=size(GG,1); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid % save figure f=figure (’ visible ’,’ off ’) ; imagesc(nu(nu_min:nu_max),x,squeeze(abs(GG(:,y_coord,nu_min:nu_max)))) axis square title ( strcat (’ z=’, string (z) ,’ mm ’,’//’,’ nu0=’,string (nu0) ,’ THz’)) xlabel (’\ nu, THz’) xlim([xlim1 xlim2]) ; % set limits for x scale ylabel (’ x, mm’) ylim([ylim1 ylim2]) ; % set limits for y scale %caxis([ caxis1 caxis2 ]) % set color scale . comment it to use autoscale c=colorbar; c .Label. String =’|G|,a.u .’; cpos=get(c ,’ Position ’) ; cpos (1) =0.76;cpos (2) =0.205;cpos(3)=0.03;cpos (4) =0.685; set (c ,’ Position ’, cpos) ; xxix.

(32) %colormap(cmocean(’amp’)); % set special colormap set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f , strcat (pt ,’\\’, num2str(prfx) ,’ _ ’, num2str(z) ,’. png’) ) ; close ( gcf ) ; end % visual_PHxnu_y.m: visualization of spatio −spectral phase of G(x,y= y_coord,nu) function visual_PHxnu_y(GG,nu,z,Mask_Am,x_size,y_coord,nu_min,nu_max, prfx,path,type,caxis1,caxis2,xlim1,xlim2,ylim1,ylim2,nu0) pt=strcat (path ,’\\’, type ,’\\ PH(x,nu)_y’); % create path PH(x,nu)_y to save figures if ~exist(pt ,’ dir ’) mkdir(pt) ; end Nx=size(GG,1); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid % save figure f=figure (’ visible ’,’ off ’) ; imagesc(nu(nu_min:nu_max),x,squeeze(angle(GG(:,y_coord,nu_min:nu_max)))) axis square title ( strcat (’ z=’, string (z) ,’ mm ’,’//’,’ nu0=’,string (nu0) ,’ THz’)) xlabel (’\ nu, THz’) xlim([xlim1 xlim2]) ; % set limits for x scale ylabel (’ x, mm’) ylim([ylim1 ylim2]) ; % set limits for y scale xxx.

(33) %caxis([ caxis1 caxis2 ]) % set color scale . comment it to use autoscale c=colorbar; c .Label. String =’Phase,rad ’; cpos=get(c ,’ Position ’) ; cpos (1) =0.76;cpos (2) =0.205;cpos(3)=0.03;cpos (4) =0.685; set (c ,’ Position ’, cpos) ; %colormap(cmocean(’balance’)); % set special colormap set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f , strcat (pt ,’\\’, num2str(prfx) ,’ _ ’, num2str(z) ,’. png’) ) ; close ( gcf ) ; end % visual_Ext_y.m: visualization of spatio −temporal amplitude of E(x,y= y_coord,t) function visual_Ext_y (EE,t,z,Mask_Am,x_size,y_coord,t_min,t_max,prfx,path, type,caxis1, caxis2 ,xlim1,xlim2,ylim1,ylim2) pt=strcat (path ,’\\’, type ,’\\ E(x, t )_y’) ; % create path E(x, t )_y to save figures if ~exist(pt ,’ dir ’) mkdir(pt) ; end Nx=size(EE,1); x=linspace(−x_size/2, x_size/2−x_size/Nx,Nx); % set x grid % save figure f=figure (’ visible ’,’ off ’) ; imagesc( t (t_min:t_max),x,squeeze(EE(:,y_coord,t_min:t_max))) xxxi.

(34) axis square title ( strcat (’ z=’, string (z) ,’ mm’)) xlabel (’ t , ps ’) xlim([xlim1 xlim2]) ; % set limits for x scale ylabel (’ x, mm’) ylim([ylim1 ylim2]) ; % set limits for y scale %caxis([ caxis1 caxis2 ]) % set color scale . comment it to use autoscale c=colorbar; c .Label. String =’E/E0’; cpos=get(c ,’ Position ’) ; cpos (1) =0.76;cpos (2) =0.205;cpos(3)=0.03;cpos (4) =0.685; set (c ,’ Position ’, cpos) ; %colormap(cmocean(’balance’)) % set special colormap set (gca ,... ’FontSize ’,10,’ FontWeight’,’Bold ’,... ’FontName’,’Times ’,... ’Box ’,’ on’) set ( gcf ,’ PaperUnits ’,’ inches ’,’ PaperPosition ’,[0 0 3 2]) saveas ( f , strcat (pt ,’\\’, num2str(prfx) ,’ _ ’, num2str(z) ,’. png’) ) ; close ( gcf ) ; end. xxxii.

(35) ABBREVIATIONS. B M 3D. operator of BM3D denoising algorithm (block matching 3D). C. cross-correlation function. CCF. operator of CCF denoising algorithm (complex-domain cube filter). C DB M 3D. operator of CDBM3D denoising algorithm (complex-domain block matching 3D). Cx. boundary condition for x component. Cy. boundary condition for y component. F. lens focus. G. temporal spectrum. H. radiation transfer function. He. Heaviside function. I1. intensity in the first arm of balance detector. I2. intensity in the second arm of balance detector. M. opaque mask. M. number of points in spatial grid. N. number of points in time grid. O(x, y, ν). spatio-spectral function of the object wave. Pµ 1. distribution of Poissonian variable for I1. Pµ 2. distribution of Poissonian variable for I2. QT (x, y, t , z r ). spatio-temporal array as a collection of real-valued noisy measurements xxxiii.

(36) R(x, y, ν). spatio-spectral function of the reference wave. T. amplitude transmittance of the object. V B M 3D. operator of VBM3D denoising algorithm (video block matching 3D). Vph. phase velocity. ∆. Laplace operator. ∆I. difference between two Poissonian variables described by the Skellam distribution. ∆ν. sample intervals in frequency grid. ∆t. sample intervals in time grid. ∆x. sample intervals in space grid. ∆k x. sample intervals in spatial frequency grid. δ(q). Dirac delta function. ϵ̂. random value in accordance to the normalized Gaussian distribution N (0, 1). µ. permeability. ∇. Del operator. ν. oscillation linear frequency. ν0. critical frequency defining AS and C methods. ω. oscillation angular frequency. ρ. transverse size of the beam. ρ. extraneous electric charge density. σ. noise standard deviations. τ. pulse duration. t̃ 0. time window size. ϵ. permittivity. ϕ. phase of the pulse. B⃗. vector of magnetic induction xxxiv.

(37) ⃗ D. vector of dielectric flux density. ⃗ E. vector of electric field. ⃗ H. vector of magnetic field. ⃗j. vector of electric current density. b (x, y, ν). wavefront distortion. c. speed of light in vacuum. d. object height. g. spatio-temporal spectrum. h. system pulse response in C method. k. cyclic wave number. kx. cyclic spatial frequency in x component. ky. cyclic spatial frequency in y component. kz. cyclic spatial frequency in z component. n(ν). refractive index dispersion. o. Object transmittance function. q. frequency in cross-correlation function. r. distance between the object plane and registration plane. t hs. threshold factor in CDBM3D. t ht. threshold factor in VBM3D. QV (x, y, ν, z). spatio-spectral array as a collection of complex-valued noisy measurements. n o var k̂. variance. F−1 2D. inverse two-dimensional Fourier transform. F12D. direct two-dimensional Fourier transform. AS. Angular spectrum method. BD. balance detector. BUTCH beam. broadband uniformly topologically charged (BUTCH) beams. xxxv.

(38) C. Convolution method. FFT. Fast Fourier Transform. fs. femtosecond. HI. hyperspectral interferometry. NRMSE. Normalized Root Mean Square Error. OAM. Orbital angular momentum. PCA. photoconductive antenna. RMSE. Root Mean Square Error. RRMSE. Ralative Root Mean Square Error. SPP. spiral phase plate. THz. terahertz. ToF. Time-of-flight data. xxxvi.

(39) 1. INTRODUCTION. Terahertz radiation, T-rays, T-waves, T-light, T-lux or THz consists of electromagnetic waves in frequency range 0.1 - 10 THz. However, the upper range is arbitrary and could be expanded up to 30 THz in some papers [18]. One terahertz is 1012 Hz or 1000 GHz. Wavelength of radiation in THz range corresponds from 30 µm to 3 mm (1 THz = 300 µm). Sometimes THz radiation is known as submillimeter waves (especially in astronomy). This type of electromagnetic waves could be also regarded as microwave or far infrared waves. Progress in femtosecond (fs) laser technology, microelectronics and digital signal processing allows to obtain and detect electromagnetic terahertz (THz) waves which are extremely short in the number of oscillations. Pulsed THz waves are usually recorded using the techniques of coherent detection with use of fs laser pulses. The result of such measurements is not the envelope or radiation power but the temporal form of the electric field itself. Typically, THz pulse represents a wave of several electric field oscillations with a corresponding spectrum in interval 0.110 THz. Technologies based on the use of broadband THz radiation, generated by different methods [10, 68, 77], have become an important tools for many practical applications nowadays, such as non-invasive diagnostic and control methods in industry [9, 66], medicine [20, 64] and biology [65, 67], safety and security [33]. In the past decade, another large-scale application of the THz range has been actively explored – THz wireless communications [24, 52]. In addition, unlike continuous THz radiation, pulsed has an ultra-wide spectrum. Its features, such as extremely low absorption by dielectrics, and the presence of characteristic frequencies of various organic molecules make it possible to carry out spectroscopic diagnostics of objects and materials. Due to the development of techniques for obtaining complete information about the structure and dynamics of THz fields, the studying of the evolution of structured beams of broadband THz radiation has become urgent [72, 74]. The existing 37.

(40) coherent detection technique allows the direct registration of the temporal form of electric field amplitude, whence the phase information also becomes available. The valuable result of the progress in THz radiation was the development of many new scientific directions, an important place among which is given to techniques for obtaining complete information about the wavefield in space and time. It can be named THz holography in accordance with its literal meaning - from the Greek words "holos" - whole and "graphe" - writing. One of the most characteristic features of visualization techniques today is the transition from analog to digital registration, which is implemented in most cases using matrix photodetectors. With regard to holography, this has led to the development of digital holography methods, as well as the development of new approaches implemented primarily through numerical processing methods. These prerequisites made it possible to develop the method of digital pulsed hyperspectral holography with time resolution and its applications. Historically, this method was proposed in [23], where authors present THz pulse holography for imaging of binary amplitude objects. They present referenceless THz pulse holography based on direct measurements of temporal form of THz pulse in each points of wide-aperture collimated beam passed through the investigated object. This become a huge advantage of THz holography against THz time-domain spectroscopy (TDS). For different tasks previously solved by TDS, in particular for determining the dielectric properties of materials, it is important to take into account effects associated with dynamics of THz wavefront propagation. However, TDS does not allow this accounting, but THz holography does. Also holographic approach allows to manage another valuable aspect: usually, the dependence of a pulsed electric field on spatial coordinates during its propagation is often considered separately, under the assumption that the temporal properties of the pulse remain unchanged along the beam. This is similar to the representation of an electric field through the product of time and space factors independent of each other. However, such a representation is often incorrect. For the pulses containing only a several oscillations of a wavefield and having a transverse size comparable with the central wavelength, unusual properties are manifested associated with diffraction and dispersion of the propagation media. Dispersion effects and angular spatio-temporal rearrangement of the spectral components during its propagation in dispersive elements lead to the interconnection of. 38.

(41) temporal and spatial structures. Such a phenomenon is called “spatio-temporal couplings” (STCs) [4]. Thus, it is important to consider simultaneously the spatial and temporal reshapes of the wavefield of such extremely short pulses (in this case, THz pulse is the limiting case of a pulse consisting of only one oscillation of the electric field).. 1.1 Thesis contributions The method of digital hyperspectral THz holography is developed, which allows to study the spatio-temporal and spatio-spectral dynamics of vector wavefields, including transverse and longitudinal components of the THz field and its polarization state. Digital processing of wavefields is implemented, which provides information on the spatial structure of the complex electric field taking into account the effects such as wave diffraction and dispersion of refractive index of the object and propagation medium. The possibility of digital hyperspectral THz holography to study the dynamics of complex-structured wavefields is demonstrated. The application of data transmission by dynamical interference of pairs of vortex beams is demonstrated. The methods to increase the resolution and image quality reconstruction in THz hyperspectral holography are investigated with use block-matching algorithms adapted for three dimensional complex-valued wavefields. The features of correlative measurements of broadband laser fields of femtosecond pulsed is realized using hyperspectral holography approach.. 1.2 Thesis and published results The main results on the topic of the thesis are presented in 16 papers. Results of Chapter 2 about method of terahertz hyperspectral holography are presented in [7, 43, 47, 57]. Results of Chapter 3 about applications of THz hyperspectral holography for beams analysis are presented in [40, 41, 44, 45, 58, 59, 60]. Results of Chapter 4 about hyperspectral data denoising are presented in [37, 38, 39, 42]. 39.

(42) Results of Chapter 5 about peculiarities of correlation measurements of pulsed hyperspectral optical fields are presented in [46].. 40.

(43) 2. TERAHERTZ HYPERSPECTRAL HOLOGRAPHY. This chapter is devoted to the method of digital hyperspectral THz holography, its theoretical foundations, computational features, and physical aspects related to wavefield propagation including diffraction and dispersion of the object’s refractive index and propagation medium.. 2.1 Theoretical model 2.1.1 General case Now consider this model in general case. In wave optics theory the propagation of electromagnetic wave is described by Maxwell’s equations. In differential form these equations could be formalized as follows (in CGS system): ⎧ ⃗ = − 1 ∂ B⃗ ⎪ ∇×E ⎪ c ∂t ⎪ ⎪ ⎪ ⎪ ⎨ ∇×H ⃗ = 1 ∂ D⃗ + 4π ⃗j c ∂t c ⎪ ⃗ ⎪ ∇D = 4πρ ⎪ ⎪ ⎪ ⎪ ⎩ ⃗ ∇B = 0.. (2.1). ⃗ is dielectric flux density, ⃗ and H ⃗ are vectors of electric and magnetic field, D Here E B⃗ is magnetic induction, ⃗j is electric current density, ρ is extraneous electric charge density, ∇ is Del operator, t is time; c speed of light in vacuum. Further we will consider only dielectric homogeneous isotropic medium where external currents and ⃗ = 0 corresponds charges are absent: ⃗j = 0 and ρ = 0. Note, that the case when ∇D to the solenoidal vector field. 41.

(44) To solve Maxwell’s equations for specific medium they should be supplemented with material equations. Consider the linear case when: ⎧ ⎨ D ⃗ = ϵE ⃗, ⎩ B⃗ = µH ⃗.. (2.2). Here ϵ is permittivity and µ is permeability of the propagation medium. Then first two equations in (2.1) could be rewritten: ⎧ ⎨ ∇×E ⃗ = − µ ∂ H⃗ , c ∂t ⎩ ∇×H ⃗ =. ⃗ ϵ∂E c ∂t.. Now consider the harmonic wave: ⎧ ⎨ E ⃗ =E ⃗ e −iωt , 0 ⎩ H ⃗ =H ⃗ e −iωt . 0. (2.3). (2.4). ⃗ and H ⃗ are electric and magnetic amplitude, ω is oscillation angular freHere E 0 0 quency. Thus, eq.(2.3) can be rewritten according harmonic wave equations: ⎧ ⎨ ∇×E ⃗= ⎩ ∇×H ⃗. iµω ⃗ c H, i ϵω ⃗ =− c E .. (2.5). ⃗ from the Then applying operator (∇×) to the first equation in (2.5) and using ∇×H ⃗ = 0 from the third equation in second equation and solenoidal vector field with ∇E system(2.1) we obtain: ⃗ = 0. ⃗ + ω ϵµ E ∆E c2 2. Here ∆ is Laplace operator and. ω 2 ϵµ c2. (2.6). = k 2 , k is a wave number. Note, that eq.(2.6) is. known as a spectral analogue of wave equation. The eq.(2.6) is linear and k is scalar, thus we can solve it for each Cartesian component of the field independently. Now consider that wave propagates along the direction z and x, y are transverse. 42.

(45) axis. It can be formalized as follows: ⎧ ⎪ E →0 ⎪ ⎪ ⎨ x,y,z. (2.7). x → ±∞ ⎪ ⎪ ⎪ ⎩ y → ±∞. In general case pulsed THz or femtosecond wavefront can have both a wide temporal and spatial spectrum. In this way, the spatio-temporal field distribution E(x, y, t , z) could be represented as its corresponding spectrum g k x , ky , ω, z or vice versa by direct and inverse Fourier transform: ⎧ ∞ ∫︁ ∞ ∫︁ ∞ ∫︁ ⎪ ⎪ g k , k , ω, z = E x,y,z (x, y, t , z) · ⎪ x,y,z x y ⎪ ⎪ ⎪ −∞ −∞ −∞ ⎪ ⎪ ⎪ ⎨ ·e −i (ωt +kx x+ky y ) d x d y d t , ∞ ∫︁ ∞ ∫︁ ∞ ∫︁ ⎪ 1 ⎪ , k , ω, z · E y, t , z) = g k (x, ⎪ x,y,z x y x,y,z ⎪ (2π)3 ⎪ ⎪ −∞ −∞ −∞ ⎪ ⎪ ⎪ ⎩ ·e i (ωt +kx x+ky y ) d k d k d ω. x. (2.8). y. Firstly consider the solutions for transverse g x,y components. Applying this expression of E x,y (x, y, t , z) to eq. (2.6) one can obtain the equation describing the dynamics of spatio-temporal spectrum propagation: ∂ 2 g x,y ∂ z2. + g x,y k 2 − k x2 − ky2 = 0.. (2.9). The solution of this equation gives the following: g x,y = C0 e −i. ⎷. k 2 −k x2 −ky2 z. + C1 e i. ⎷. k 2 −k x2 −ky2 z. .. (2.10). Here the first term describes the propagation of direct wave and the second term is reverse wave. Further only the propagation dynamics of the direct wave will be considered. In eq. (2.10) C0 is the integration constant and further we will assume C x as a boundary condition for x component and Cy for y components. After the solutions of scalar equations for the transverse components it is possible to find solution ⃗ = 0 it is derived as for longitudinal component. From solenoidal field equation ∇E. 43.

(46) follows:. ∂ E x ∂ Ey ∂ Ez =− − . ∂z ∂z ∂z. (2.11). Thus, applying the expression of E x,y,z (x, y, t , z) from the second eq. in system (2.8) to eq. (2.11) one can obtain: ∂ gz = −i k x g x − i ky gy . ∂z. (2.12). The solution of this equation gives: k x C x + ky Cy −i ⎷k 2 −k 2 −k 2 z x y . gz = Æ e 2 2 2 k − k x − ky. (2.13). Thus, in this section we considered the wavefronts whose spatial and temporal spectrum could be broad. Such waves have transverse size comparable with the central wavelength, and the temporal wave duration is comparable with the central oscillation period. Resuming the equations for the general case of nonparaxial propagation of an individual monochromatic component of the spectrum for the transverse and longitudinal components over a distance z could be described as follows: ⎧ ⎪ ⎪ ⎪ ⎨. g x k x , ky , ω, z = C x k x , ky , ω e −i kz z , gy k x , ky , ω, z = Cy k x , ky , ω e −i kz z , ⎪ ⎪ ⎪ ⎩ g k , k , ω, z = kx ·Cx (kx ,ky ,ω)+ky ·Cy (kx ,ky ,ω) · e −i kz z . z x y k. (2.14). z. Here we denote k z =. q. k 2 − k x 2 − ky 2 . Note, that solutions when radicand is < 0. correspond to evanescent waves.. 2.1.2 Specific cases of Fourier transforms For practical reasons, it may be useful to take a one-dimensional Fourier transform to return from spectral G(x, y, ν, z) to the temporal representation: E (x, y, t , z) =. ∫︂∞. G (x, y, ν, z) e i2πν t d ν.. −∞. 44. (2.15).

(47) One of the necessary conditions for the correct operation of the engine is the conservation of energy during the transition from the time representation of the signal to the spectral one and back. Energy conservation is checked using the Parseval ratio, which in integral form is written as: ∫︂∞. ∫︂∞. 2. |E (t )| d t =. −∞. |G (ν)|2 d ν. (2.16). −∞. There are also more specific cases of transition between dimensions in 2.8. For example, the analysis in time-angular coordinates (k x , t ), in angles-spectral (k x , ν) [4] or in Wiener function space (k x , x) are available.. 2.2 Numerical calculation peculiarities 2.2.1 Angular spectrum method In this section we consider the numerical peculiarities of wavefront propagation. This method (further ‘Angular spectrum’ - AS) is formalized in equation system (2.14). Note, that in the numerical calculation a linear frequency ν is more convenient than cyclic ω. In addition, in THz research it is also common to apply a linear frequency: ω = 2πν; k x = 2π f x ; ky = 2π fy .. (2.17). In system (2.14) C x and Cy are boundary conditions corresponding to the initial wavefront distribution in plane z = 0. However, for general objective of wavefront propagation we can consider the wavefront g k x , ky , ω, z0 in any arbitrary plane z0 as initial condition. In numerical calculations system (2.14) is realized as follows:. 45.

(48) ⎡ ⎧ ⎪ ⎪ ⎪ −1 ⎣ ⎪ Gx (x, y, ν, z) = F2D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎡ ⎪ ⎪ ⎪ ⎪ ⎪ −1 ⎣ ⎪ ⎪ Gy (x, y, ν, z) = F2D ⎪ ⎪ ⎨ ⎡ ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ −1 ⎢ ⎪ ⎪ G y, ν, z) = F (x, ⎢ ⎪ z 2D ⎢ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎣ ⎪ ⎪ ⎩. 1 F2D [Gx (x, y, ν, z0 )] ·. ⎤. ⎦, ·H x (x, y, ν, z) ⎤ 1 F2D Gy (x, y, ν, z0 ) · ⎦, ·Hy (x, y, ν, z) 1 F2D. [Gx (x, y, ν, z0 )] ·. ⎤. (2.18). ⎥ ⎥ ⎥ ·H x z (x, y, ν, z)+ ⎥ ⎥. 1 +F2D Gy (x, y, ν, z0 ) · ⎥ ⎥ ⎦ ·Hy z (x, y, ν, z). where F12D and F−1 is the direct and inverse two-dimensional Fourier transform. 2D Physically this operation corresponds to the decomposition of wavefront into plane waves. Then the operator H in eq.(2.18) describes the radiation transfer function: ⎧ ⎪ H x (x, y, ν, z) = e −i kz z , ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ H (x, y, ν, z) = e −i kz z , y k ⎪ ⎪ H x z (x, y, ν, z) = kx e −i kz z , ⎪ ⎪ z ⎪ ⎪ k ⎩ Hy z (x, y, ν, z) = ky e −i kz z .. (2.19). z. It is useful to analyze the applicability limits of the above methods, due to the discretization, which have to face in numerical calculations. In practice, the calculation of the above equations is carried out with a discrete step along a rectangular grid. Thus, following the convention of Fast Fourier Transform we produce the discretization of eq.(2.15) according to the relations: N N N t = ∆ t p; p ∈ − 2 , − 2 + 1, ..., 2 − 1 , N N N ν = ∆ν u; u ∈ − 2 , − 2 + 1, ..., 2 − 1 .. (2.20). where p and u are integers. Besides, sample intervals ∆ t and ∆ν should be such that: ∆ t ∆ν =. 1 . N. (2.21). where N is a number of points in time grid. The similar way is to be in spatial grid. 46.

(49) for spatial Fourier transform: M M M x = ∆ x m x ; m x ∈ − 2 , − 2 + 1, ..., 2 − 1 , M M M y = ∆y my ; my ∈ − 2 , − 2 + 1, ..., 2 − 1 , M M M k x = ∆kx n x ; n x ∈ − 2 , − 2 + 1, ..., 2 − 1 , M M M ky = ∆ky ny ; ny ∈ − 2 , − 2 + 1, ..., 2 − 1 .. (2.22). where M is a number of points in spatial grid. Here the sample intervals ∆ x and ∆kx should also be such that: ∆ x ∆k x =. 1 . M. (2.23). The equation for wavefront propagation calculation in discrete grid looks as: ⎡ ⎧ ⎪ ⎪ ⎪ −1 ⎣ ⎪ Gx m x , my , u, z = F2D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎡ ⎪ ⎪ ⎪ ⎪ ⎪ −1 ⎣ ⎪ ⎪ Gy m x , my , u, z = F2D ⎪ ⎪ ⎨ ⎡ ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ −1 ⎢ ⎪ ⎪ Gz m x , my , u, z = F2D ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎣ ⎪ ⎪ ⎩. ⎤ Gx m x , my , u, z0 · ⎦, ·H x n x , ny , u, z ⎤ 1 F2D Gy m x , my , u, z0 · ⎦, ·Hy n x , ny , u, z ⎤ 1 F2D Gx m x , my , u, z0 · ⎥ ⎥ ⎥ ·H x z (n x , ny , u, z)+ ⎥ ⎥ . 1 +F2D Gy m x , my , u, z0 · ⎥ ⎥ ⎦ ·Hy z (n x , ny , u, z) 1 F2D. . . (2.24). Here H is transfer function defined in eq. (2.18) as: ⎧ ⎪ H n , n , u, z = e −i nz ·z , ⎪ x x y ⎪ ⎪ ⎪ ⎪ ⎨ H n , n , u, z = e −i nz ·z , y x y n ⎪ ⎪ H n , n , u, z = nx e −i nz ·z , ⎪ x z x y ⎪ ⎪ nz ⎪ ⎩ Hy z n x , ny , u, z = ny e −i nz ·z . z. Here n z =. r. 2πu 2 − n x2 c. − ny2 .. 47. (2.25).

(50) 2.2.1.1. Discrete Fourier transform peculiarities. Consider the discrete analogue of the eq.(2.15) in direct and inverse form is: ⎧ ⎪ ⎪ ⎪ ⎨ E ( p) = ⎪ ⎪ ⎪ ⎩ G (u) =. 1 N. N ∑︁ /2−1. u=−N /2 N ∑︁ /2−1. p=−N /2. 2π. G (u) e i N. E ( p) e. −i. up. 2π N up. , (2.26). .. Parseval relation 2.16 in discrete form is written as: N −1 ∑︂ p=0. |E ( p)|2 =. −1 1 N∑︂ |G (u)|2 N u=0. (2.27). This equation controls the conservation of energy and serves as one of the criteria for the correctness of the software implementation. This transition from the temporary representation of the data to spectral and vice versa is presented in Fig. 2.1. The energy conservation according to the eq. (2.27) is proved.. Figure 2.1 The transition from a temporary representation of the data to spectral and vice versa with an estimate of energy conservation according to the eq.(2.27).. Since only positive frequencies have the physical meaning it is redundant to carry out a numerical calculation for all frequencies. Thus, the information contained in the positive frequencies is sufficient to reconstruct the temporal shape of the signal due to the Hermitian symmetry of the complex spectrum: G (−u) = G ∗ (u). To restore the initial distribution of the temporal shape of the pulse without losses of temporal resolution, it is possible to use the mathematical representation of the real48.

(51) valued signal E( p) in the form of a complex-valued analytical time function A[E( p)] (see for example [63]). This analytical signal is represented by the following expression: A(E ( p)) = i F F T (H (u) · G (u)). (2.28). where H is Heaviside piecewise constant function, equal to zero for negative frequencies and one for positive ones. The doubled value of the real part of the resulting complex-valued analytical signal will correspond to the signal amplitude in a temporal representation: E ( p) = 2Re [A(E ( p))]. (2.29). In numerical modeling, it is enough to determine the complex signal spectrum by zero elements in the negative frequency branch, then use the inverse FFT and consider the double value of real part of the obtained complex number.. 2.2.2 Convolution method As an alternative, the diffraction propagation of a broadband field into the z plane can be represented as a convolution of the field in the plane z0 with the pulse response of the system: G (x, y, ν, z) = G (x, y, ν, z0 ) ∗ h(x, y, ν, z).. (2.30). Here h is the system’s pulse response, defined as: ν e −i2πr ν c z , −i c r 2 −1. h(x, y, ν, z) = where r =. p. (2.31). z 2 + x 2 + y 2 is the distance between the object plane and registration. plane. According to the property of Fourier transforms, we can rewrite the convo-. 49.

(52) lution (eq. 2.30) into the form similar to eq. (2.18): ⎡ ⎧ ⎪ ⎪ ⎪ −1 ⎣ ⎪ Gx (x, y, ν, z) = F2D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎡ ⎪ ⎪ ⎪ ⎪ ⎪ −1 ⎣ ⎪ ⎪ Gy (x, y, ν, z) = F2D ⎪ ⎪ ⎨ ⎡ ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ −1 ⎪ ⎪ G y, ν, z) = F (x, ⎢ ⎪ z 2D ⎢ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎣ ⎪ ⎪ ⎩. 1 F2D [Gx (x, y, ν, z0 )] ·. ⎤. ⎦, 1 ·F2D [h(x, y, ν, z)] ⎤ 1 F2D Gy (x, y, ν, z0 ) · ⎦, 1 ·F2D [h(x, y, ν, z)] 1 F2D. [Gx (x, y, ν, z0 )] ·. (2.32). ⎤. ⎥ ⎥ ⎥ [h(x, y, ν, z)] ⎥ ⎥. 1 +F2D Gy (x, y, ν, z0 ) · ⎥ ⎥ ⎦ ky 1 ·F2D [h(x, y, ν, z)] k kx + kz. 1 ·F2D. z. where the transfer function is the Fourier transform of the impulse response of the system h. Accordingly, all x, y, z field components can also be calculated by the convolution method, similarly to the system of equations (2.18). In a discrete form, eq. (2.32) could be rewritten as: ⎡ ⎧ ⎪ ⎪ ⎪ −1 ⎣ ⎪ Gx m x , my , u, z = F2D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎡ ⎪ ⎪ ⎪ ⎪ ⎪ −1 ⎣ ⎪ ⎪ Gy m x , my , u, z = F2D ⎪ ⎪ ⎨ ⎡ ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ ⎢ ⎪ ⎪ −1 ⎢ ⎪ ⎪ G m , m , u, z = F ⎢ ⎪ x y 2D ⎢ ⎪ z ⎪ ⎪ ⎢ ⎪ ⎪ ⎣ ⎪ ⎪ ⎩. ⎤ 1 F2D Gx m x , my , u, z0 · ⎦, 1 ·F2D h m x , my , u, z ⎤ 1 F2D Gy m x , my , u, z0 · ⎦ , 1 ·F2D h m x , my , u, z 1 F2D Gx m x , my , u, z0 · n 1 ·F2D h m x , my , u, z nx + z 1 +F2D Gy m x , my , u, z0 · n 1 ·F2D h m x , my , u, z ny. ⎤. (2.33). ⎥ ⎥ ⎥ ⎥ ⎥. ⎥ ⎥ ⎦. z. and pulse response is: . . h m x , my , u, z =. u e −i 2π −i c. ⎷. z 2 +m x2 +my2 u c −1. z 2 + m x2. + my2. z .. (2.34). Note, that calculation of z component by convolution method supposes to use k x /k z and ky /k z from general approach. Thus the transfer functions from both numerical 50.

(53) methods are related to each other: ⎧ 1 ⎪ H x = F2D [h] , ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ H = F 1 [h] , y 2D k 1 ⎪ ⎪ H x z = F2D [h] kx , ⎪ ⎪ z ⎪ ⎪ k ⎩ 1 Hy z = F2D [h] ky .. (2.35). z. Further the numerical realization of this wavefront propagation methods will be demonstrated. We will use the notations in system of Eq. (2.18) as angular spectrum (AS) method and system (2.32) as convolution (C) method.. 2.2.3 The choice of the calculation method Pulse THz radiation has a wide spectrum, and therefore, often the calculation of wavefront propagation performed for all spectral components using only one approach leads to numerical errors. In the case when the wavefront is not very curved, the sampling theorem allows the usage of two different methods depending on the ratio between frequencies and distances. The choice of calculation method for an individual spectral component is determined depending on the critical frequency ν0 : νc r =. c·z Dx · ∆x. ,. (2.36). where D x is the linear size of the rectangular calculation grid, ∆ x - pixel size. Calculation through the angular spectrum can be done if ν ⩾ ν0 , while the convolution representation of the field with an impulse response is appropriate if ν ⩽ ν0 . Fig. 2.2 illustrates this principle of choosing a calculation method.. 51.

(54) 1.00. (a). 0.75 |G/G0|. (b). Initial wavefront i<. cr. 0.50. i >=. i. AS. C. 0.25 AS. C. 0.00 0. 1. cr. cr. 2 / 0. 3. 4. Propagated wavefront. Figure 2.2 Explanation of the use of calculation methods in depending on the critical frequency νc r in the spectrum: (a) -visual concept, (b) - block diagram.. 2.3 Numerical Results 2.3.1 Pulsed THz wave Now consider the wavefront propagation on example of THz single-period wave defined as in [35]:. 2 t t E(t ) = E0 exp − . τ τ2. (2.37) t̃. t̃. Here E0 is the electric field amplitude, τ is pulse duration, t ∈ [− 20 ; 20 ], and t̃ 0 is time window size. The eq. (2.37) corresponds to the approximation of THz pulses generated by optical rectification of fs pulses in non-linear crystals [76] or by a photoconductive emitters [10] or in air plasma [5] irradiated by pulsed fs lasers. Thus, eq. (2.37) describes the THz pulse temporal form E(t ). This ideal single-period THz waveform is close to the pulses obtained in different experiments (see for example pulses in Fig. 2.3). The signal described by eq. 2.37 is real-valued and has no temporal phase φ(t ) in explicit form. But the entire temporal phase of any signal could be calculated as the complex phase of its analytical signal (see for example [63]). This approach also allows to estimate the instantaneous frequency of the signal and its temporal envelope. However, for a single-cycle THz pulse the temporal phase and envelope are not of significant interest and remain beyond the scope of this consideration, 52.

(55) Figure 2.3 Single-period THz wave given by eq. (2.37) and pulses obtained in several experiments.. since the concept of the envelope and temporal phase for a pulse with one oscillation loses its physical content [36]. More important in assessing the physical properties of the time signal given by Eq. 2.37 is its spectral characteristics. The solution of the Fourier integral for the field E(t ) allows to obtain an expression for its corresponding complex spectrum G(ν): G (ν) =. ∫︂∞. 3 E (t ) exp (−i2πν t ) d t = −iπ 2 E0 τ 2 ν exp −π2 τ 2 ν 2 .. (2.38). −∞. In this case, the obtained spectral density is complex, which allows to consider separately its amplitude |G(ν)| and phase ϕ(ν). Note, that phase components are synchronized for all frequencies in the considered THz pulse. Spectral amplitude and spectral phase of considered THz pulse are presented in Fig. 2.4. 53.

(56) Figure 2.4 Spectral amplitude (left) and spectral phase (right) of THz pulse.. 2.3.2 THz wavefront Consider the model of THz wavefront. Firstly we assume that initial THz wavefield emitted from the source is ideally collimated. Then, in general case, the spatial structure of THz field could be formalized by Gaussian distribution with transverse size ρ:. 2 x + y2 G(x, y, ν) = G (ν) · exp − ρ2. (2.39). In the limiting case when ρ = ∞ the eq. (2.39) defines a plane transversely uniform wave. Then, this Gaussian beam is exposed to amplitude and phase modulation (Fig. 4.13). These modulation components could be represented by the arbitrary amplitudephase object or by any system of elements for the objective of beam shaping. For example, these components could be phase axicons [44, 56], spiral phase plates [59] or specially composed systems like a THz broadband uniformly topologically charged vortex beam generator [40] or even more complicated devices. In plane z0 the wavefield gets a determinate structure according to the modulation components. Then the wavefront propagates through the medium and at the arbitrary distance z the electric field is registered. Now consider the frequency dependent amplitude-phase modulation (caused by arbitrary object). Object transmittance function will be considered as o (x, y, ν) = 54.

(57) Figure 2.5 Schematic representation of pulsed THz vectorial wavefield propagation.. T (x, y, ν) exp(i · ϕo b j (x, y, ν)). Here T represents the amplitude transmittance and ϕo b j is the phase shift induced by the object. This phase delay is defined by the object height d (x, y) and refractive index no b j (x, y, ν): ϕo b j (x, y, ν) =. 2πν no b j (x, y, ν) − 1 d (x, y) c. (2.40). Thus, the impact of the object to the resultant spatio-spectral wavefield is describe as: G (x, y, ν, z0 ) = G (x, y, ν) · o (x, y, ν) = 2πν = G (x, y, ν) · T (x, y, ν) exp i c no b j (x, y, ν) − 1 d (x, y) ,. (2.41). Thus we obtain the initial wavefront G (x, y, ν, z0 ) to be propagated.. 2.3.3 Results of THz vectorial wavefront propagation Further consider the numerical realization of AS and C methods on the example of a circular aperture diffraction. However, the peculiarities are associated with THz wavefront with broad temporal and spatial spectrum (see Fig. 2.6). The circular aperture can be formalized as an opaque mask M (x, y) with a radius 55.

(58) Figure 2.6 The setup of THz wavefront diffraction on the circle. (a, d) show the transverse slices of THz wavefront; (b, e) are time-domain form of THz pulse in the beam center; (d, f) are corresponding spectra.. R: M (x, y) =. ⎧ ⎨ 1, x 2 + y 2 ≤ R2 ; ⎩ 0, x 2 + y 2 > R2 .. (2.42). Mask is located in the initial plane z=0 and it modulates two-dimensional complex spectrum as G (x, y, ν, z0 ) = G (x, y, ν) · M (x, y). The following is the results of calculating of wavefield propagation for the angular spectrum method, convolution method, and their combination depending on the limiting frequency ν0 . The calculation is shown for the spatial-frequency distribution of the modulus of the spectral density of radiation (Fig. 2.7). The AS method gives an accurate result in the range above ν0 , while below ν0 it leads to calculation errors. Similarly, the C method gives a correct result within less than ν0 and gives numerical errors in the region higher than critical frequency. Combining of these two methods along the boundary ν0 gives the correct calculation result for all frequencies. This is also clearly confirmed on one-dimensional graphs displaying a point in the center of the beam (Fig. 2.8). It is useful to show how the corresponding spatial-temporal structure behaves. 56.

(59) Figure 2.7 Spatial-frequency distribution of the spectral density modulus of diffracted radiation (transverse field component) calculated by the angular spectrum method (upper row), convolution method (middle row) and by combining two methods (lower row). Propagation distances correspond to the columns 0.25z0 , 0.5z0 , 1z0 , 2z0 . z0 = 180m m . Critical frequency ν0 is marked by dot line.. Fig. 2.9 shows the wavefields for the angular spectrum method, convolution method, and their combination. It can be seen how the combination of the two methods compensates the errors caused by each method used outside its applicability. Further, the result of calculating the wavefield propagation for the angular spectrum method, convolution method, and their combination in the region ν0 (in this case, for the longitudinal field component) is presented. The calculation is shown for the spatial-frequency distribution of the modulus of the spectral density of radiation (Fig. 2.10). The angular spectrum method gives an accurate result in the range above ν0 , while below ν0 it leads to calculation errors. Similarly, the convolution method gives a correct result within less than ν0 and gives numerical errors in the region of a higher critical frequency. Combining the two methods along the boundary ν0 gives the correct calculation result for all frequencies. 57.