Speckle patterns formed by broadband terahertz radiation and their applications for ghost imaging

Speckle patterns can be very promising for many applications due to their unique properties. This paper presents the possibility of numerically and experimentally formation of speckle patterns using broadband THz radiation. Strong dependence of the statistical parameters of speckles, such as size and sharpness on the parameters of the diffuser are demonstrated: the correlation length and the mean square deviation of the phase surface inhomogeneity. As the surface correlation length is increasing, the speckle size also increases and its sharpness goes down. Alternatively, the magnification of the standard deviation of the surface height leads to the speckle size diminishing and growth of the speckle sharpness. The dimensions of the experimentally formed speckles correspond to the results of numerical simulation. The possibility of utilizing formed speckle patterns for the implementation of the ghost imaging technique has been demonstrated by methods of numerical modeling.

Speckle patterns can be very promising for many applications due to their unique properties. This paper presents the possibility of numerically and experimentally formation of speckle patterns using broadband THz radiation. Strong dependence of the statistical parameters of speckles, such as size and sharpness on the parameters of the diffuser are demonstrated: the correlation length and the mean square deviation of the phase surface inhomogeneity. As the surface correlation length is increasing, the speckle size also increases and its sharpness goes down. Alternatively, the magnification of the standard deviation of the surface height leads to the speckle size diminishing and growth of the speckle sharpness. The dimensions of the experimentally formed speckles correspond to the results of numerical simulation. The possibility of utilizing formed speckle patterns for the implementation of the ghost imaging technique has been demonstrated by methods of numerical modeling.
Speckle pattern formation is a phenomenon that occurs in situations where a coherent light reflects or transmits through a rough surface or a turbid medium 1,2 . Despite the fact that the occurrence of speckle-noise can affect negatively for some imaging techniques (e. g. optical coherence tomography and LIDAR technologies), there are fields in which usage of speckle patterns can be very promising: it applied for image retrieving in speckle metrology 3 , single pixel imaging 4 , for ghost imaging [5][6][7][8] , cryptography technique 9 etc. In turn, terahertz (THz) radiation, due to the optical properties of various materials in this frequency range, as well as the features of individual properties of the radiation itself, is becoming more popular in quality control systems [10][11][12] , safety [13][14][15][16] and biomedicine application [17][18][19] , communication [20][21][22][23] .
The formation of speckle patterns for THz radiation begins to be studied in detail only at the present time. The first speckle structures in the THz region were first obtained using monochromatic continuous THz radiation and discussed in papers 24,25 . All these results were experimentally tested using a free electron laser (FEL). In [25][26][27] formation of speckle patterns by reflection FEL THz radiation from a random rough surface was also considered. The use of other sources of monochromatic radiation is presented in 28 . In this paper, for the formation of speckles, metal disk with a set of various random amplitude images around the disk circumference was used. In works [25][26][27][28] , the possibility of application such speckle structures for the ghost imaging technique was demonstrated. High energies, various amplitude and phase objects transmission are the main advantages of monochromatic continuous THz radiation.
Usually, pulsed THz radiation has low energies 29,30 , but there are already widespread methods of generating high-energy radiation, for example, such as in organic crystals 31,32 or by filamentation in gases or liquids 33,34 . The generation of speckle structures by pulsed THz radiation has not been studied yet. However, it should be noted that mask formation methods using nonlinear transformations in the crystal for generating the pulsed THz radiation exist and could be found in 35 . In this case, spatial light modulator, which modulate pump radiation is used. The modified beam is directed to the crystal for generating the THz field. This allows the ghost imaging technique in the THz range realization. A variant of amplitude modulation of probe radiation was proposed as a separate method for implementing THz imaging with sub-wavelength accuracy 36 . Producing of speckle patters

Mathematical model
In this work we study formation of speckle patterns using broadband pulsed THz radiation. There are a lot of approaches for obtaining of random light fields [37][38][39] . One of the most widespread methods propose using random refractive inhomogeneity for light modulation 40 . In our paper we expand this approach for broadband THz radiation and we form speckle patterns by means of transmission of THz pulses through the transparent plate with random rough surface on its back side. This random phase screen provides a phase modulation of THz pulses, which leads to the appearance of speckle patterns at some distance behind the screen, as shown in the Fig. 1a. In this section, a numerical model of speckle formation is presented and the dependence of the speckle parameters on the characteristics of the transparent plate with random rough surface is investigated.
Speckle pattern formation in the monochromatic and pulsed THz radiation. In the first part of the mathematical modeling, the formation of speckle structures by THz radiation, both monochromatic and pulsed, is considered. To describe speckle pattern formation for the monochromatic THz radiation we use model of Gaussian beam with the field distribution E in which can be described by: where x, y, z are the spatial coordinates, a is the spatial parameter of Gaussian beam related to the beam waist, E 0 is the initial amplitude of beam, k is the wave number. . Speckle patterns produced by (d) monochromatic THz radiation with frequency ν = 0.7 THz, produced by pulsed THz radiation with symmetric spectrum (e) and spectrum with shifted central frequency (f). Spectrum (g) with a blue curve has a symmetric Gaussian shape, spectrum with red curve has gaussian shape with shifted central frequency. As a material of plate with random rough surface on its back side with refraction index n = 1.55 for THz spectral range was used, parameters of the surface are h RMS = 0.3 mm and Cl = 0.9 mm. The distance between random phase screen and speckle pattern observation plane is = 5 cm www.nature.com/scientificreports/ Suppose, that a Gaussian beam that is incident to the random phase screen at z = 0 . Using the phase shift introduced by random phase screen, the field distribution of beam right after the random phase screen is evaluated in the following way 41 : where is a wavelength, n is the refraction index of material of plate with random rough surface on its back side, h(x, y) is the function of height profile of random rough surface. Gaussian random rough surface model is used for formation of rough structure 42 . This surface has two parameters which describe the height profile. The first parameter is the root mean square deviation of the surface height h RMS , and the second parameter is the horizontal correlation length Cl of the surface, which shows the frequency of spatial inhomogeneities of the surface. Formation of Gaussian random rough surface h(x, y) is described by the following equation: where L is the length of the surface side, N is a size of the matrix h(x i , y j ) ( i, j = 1, 2, . . . , N ) obtained due to discretization, H 1 is N × N matrix of random numbers with Gaussian distribution ϕ and standard deviation h RMS : H 1 = h RMS · ϕ , g denotes a Gaussian filter g(x, y) = exp(−2(x 2 + y 2 )/Cl 2 ) with correlation length Cl; F and F −1 are direct and inverse Fourier transforms, respectively. While propagating along the z axis over the distance , the intensity profile of the modulated beam becomes random, which results in a speckle pattern formation at the observation plane z = . Propagation of the modulated beam along z axis was calculated by means of the well-known two methods: Fresnel propagation formula (see "Methods" section) and the angular spectrum approach (plane wave decomposition, see "Methods" section). It is worth remarking that both methods provided the same numerical results for speckle pattern formation.
Consider formation of speckle patterns for pulsed THz radiation. Taking into account all spectral components of THz pulse with spectrum S(ν) , we have: Thus, we obtain time-dependent speckle pattern distribution E(x, y, �, t) . The overall speckle pattern intensity I r (x, y) that is evaluated by integrating E r (x, y, �, t) over the pulse duration T is given by: Depending on the shape of the time-domain spectrum the resulting speckle image can change considerably. Below we present a comparison of speckle images obtained using monochromatic and pulsed THz radiation. We confine our modelling by using two types of pulses. First type has Gaussian spectrum shape with central frequency ν 0 = 0.7 THz and can be described as (the blue curve in Fig. 1g): where 4τ p log 2 is the duration of the pulse, ν 0 is the central frequency of the pulse. The second type of pulse has asymmetric spectrum shape (an asymmetric spectrum is usually observed in experiments) with central frequency ν 0 = 0.5 THz (the red curve in Fig. 1g): Some typical speckle patterns for monochromatic and pulsed THz radiation are shown in Fig. 1d-f. The frequency of the monochromatic THz radiation was taken as ν = 0.7 THz, the duration of THz Gaussian pulse was 2 ps. For the plate with random rough surface on its back side with refraction index n = 1.55 for THz spectral range was used, and parameters of the rough surface were chosen as h RMS = 0.3 mm and Cl = 0.9 mm. The distance between random phase screen and speckle pattern observation plane was = 5 cm.
As it can be seen from Fig. 1, speckle patterns significantly vary depending on type of radiation [monochromatic ( Fig. 1d) and pulsed (Fig. 1e,f)] and shape of pulse spectrum (Fig. 1g). Comparing the speckle images in Fig. 1d-f it is noticeable that the number of speckles in Fig. 1e,f which is produced by THz pulses with spectrum S 1 , has been increased dramatically. It can be explained by the fact that different spectral components of the pulse produce different speckle patterns which overlap each other. This increases the overall number of speckles in the image which is the advantage of the pulsed radiation over the monochromatic one. However, the shift of the central frequency of the pulse from ν 0 = 0.7 THz in S 1 to ν 0 = 0.5 THz in S 2 leads to the decrease of contrast of speckle image, forming a lighter background comparing to speckle images for monochromatic THz radiation and THz pulses with spectrum S 1 . Shift of the central frequency of the pulse from ν 0 = 0.7 THz in S 1 to ν 0 = 0.5 THz in S 2 does not make significant changes because of using the same total spectrum range (Fig. 1g). Broadband radiation forms a lighter background comparing to speckle images for monochromatic THz radiation and THz pulses with spectrum S 1 . www.nature.com/scientificreports/ Statistical characteristics of speckle patterns. Due to the extensive use of speckle patterns in the optical range 1-8 , such possibilities of application of speckle patterns in the THz frequency range is needed to be considered. For the qualitative applications of speckle structures in the THz range, it is necessary to analyze their characteristics. In this section we describe some characteristics of speckles formed by the propagation of broadband THz radiation through the random phase screen. For statistical characteristics of speckle patterns generated by THz pulsed radiation, we considered the two parameters. The first parameter is an equivalent circular diameter (ECD) of a single speckle. This parameter is usually responsible for the spatial resolution of the reconstructed picture by various methods using speckles, for example 28 . The second parameter is described by the following equation: where I peak is a peak intensity of a speckle, �·� denotes an ensemble average over m single speckles. Parameter P can be considered as an average degree of sharpness of the speckle pattern. This parameter is more responsible for the contrast of the speckle patterns. The quantity P should increase in order to obtain single speckles that are more distinct in intensity (amplitude). The dependence of these parameters on root mean square deviation of the surface height h RMS and correlation length Cl is presented in Fig. 2.
As it can be seen from Fig. 2b, the increase of correlation length leads to the magnification of the speckle size 41 , wherein the saturation occurs at the certain value of the correlation length, after which the speckle size remains approximately constant. It can be explained by the fact that by considering broadband pulsed radiation, the use of a field with certain wavelengths affects the speckle size. In this case, the speckle size becomes comparable with the maximum wavelength containing in the spectrum of THz pulses. Moreover, according to Fig. 2c, increasing the root mean square deviation of the surface height leads to a decrease of the speckle size. An increase in the surface roughness after a certain value leads to a decrease of the speckle size, when the phase difference of the field at the minimum and maximum heights of surface described by Eq. (3) becomes greater than 2 π . This was observed in the work 41 . For speckles formed by broadband radiation, the contribution of adjacent frequencies takes place. After passing the element of inhomogeneity the steepness of the front increases so much that a strong change in the phase difference arises at neighboring points, which leads to the decrease of the speckle size. Thus, a strong increase in surface roughness (phase greater than 2 π ) is equivalent to the decrease in the correlation length. As one could expect, the sharpness of the speckles is decreasing with the growth of the correlation length as it is shown in Fig. 2d. At the same time the value of speckle sharpness grows with magnification of the root mean square deviation of the surface height due to the strong phase shift for the whole spectrum, see Fig. 2e. The study of these parameters makes an important contribution to the possibility of using formed speckle patterns in various applications. It is essential that the speckle diameter (ECD) is responsible for the spatial resolution of the images, and P is responsible for the image contrast. Thus, it is possible to control changes in the formed speckles by changing various characteristics in the phase object.

Experimental verification
THz imaging setup was used for experimental verification of broadband THz field speckle patterns formation. Schematic diagram of detection system is shown in Fig. 3a (A detailed description of the experimental setup is presented in the work 43 ). The terahertz radiation was generated due to optical rectification of femtosecond www.nature.com/scientificreports/ pulses in lithium niobate crystal 44 . The femtosecond pulses from a regenerative Ti-sapphire amplifier (pulse duration 35 fs, energy 1.1 mJ, 1 kHz repetition rate, 800 nm central wavelength) were used as pump radiation. The pulse energy was divided into two parts using a beam splitter (1:49 for probe and pump beam respectively). The probe pulse passed through the delay line and was detected by an electro-optical system (EOS). In this case, a 18 × 18 × 1 mm ZnTe crystal was used. The generated THz field with the parameters given below (see "Methods" section) had a Gaussian profile [ Fig. 3c(1), a(I)] and passed through the phase diffused plates which were printed from the plastic with the following parameters: a refractive index in the THz range n THz = 1, 55 (photos of phase objects are shown in Fig. 3b), h RMS = 0.8 mm and Cl = 1.5 mm (shown at the top of Fig. 3b), h RMS = 0.6 mm and Cl = 0.9 mm (shown at the bottom of Fig. 3b). One of the key requirements for the material from which the diffusers are made of by 3D printing to form speckle patterns with THz radiation is its transparency and a sufficiently high refractive index. In principle, other polymeric materials can also be used for these purposes, for example Tsurupica 45 , polylactide (PLA) 46 , Teflon 47 . The parameters of the phase object were selected in relation to the studies carried out in the previous section (Statistical characteristics of speckle patterns) for later comparison of the size of the speckles with the numerically calculated ones. Then, the modulated THz field (shown on the II Fig. 3a) was detected by ZnTe crystal, where the real valued spatial distribution of THz pulse in temporal domain 48 was measured. After that EOS allowed to completely detect the entire picture of THz radiation by using a probe femtosecond field (shown on the III Fig. 3a) and register it on the CMOS camera (CCD). Processing of the experimental results included: averaging in each pixel of the frame by 10 scans to increase the signal to noise ratio. Then Fourier transform was done for each pixel in order to obtain spectra at each point of the frame. Finally the radiation energy in each pixel was calculated by multiplying the amplitude at the corresponding frequency.
As it can be seen from Fig. 3c(1), a THz field with a slightly shifted Gaussian profile was suppressed at the input. However, after passing through the plates, a changes in the spatial distribution of the THz radiation intensity can be noticed . Local maxima are clearly observed, which can be attributed to individual parts of the speckle structure. Moreover, the sizes of these peaks are within the segment of 1-2 mm (Fig. 3c2,c3) which corresponds to the estimates given in Fig. 2b,c. For example, typical numerical calculation result for the parameters selected in Fig. 3c3 is shown in (Fig. 3c4). It can be seen from the figure that typical speckle sizes are also in the range of 1-2 mm. Both results (numerical calculation and experimental one) correspond to the 1-2 mm segment of speckle size calculated on the basis of the equations given in works 49,50 . The experimental data plots are similar and have a correlation coefficient of 0.9. Due to the inhomogeneous Gaussian profile of the probe optical beam that affects the speckle patterns in the THz range during registration, sizes of individual speckles (1-2 mm) and recorded field (5 mm) the high correlation coefficient is obtained. Thus, the possibility of speckle formation by broadband THz radiation has been experimentally confirmed in the work.  (1) and resulting THz speckle patterns obtained from the radiation intensity recorded on an imaging setup for the case in which THz radiation propagates through the phase plates (2 is b (top) and 3 is b (bottom)) located at a distance of 5 cm from the ZnTe crystal, 4-typical result of numerical calculation of speckle pattern formation for the same conditions as in c3

Application in ghost imaging algorithm
In ghost imaging technique classical laser beams typically are used to reconstruct ghost images. In this work we employ speckle patterns obtained numerically for both monochromatic and broadband pulsed THz radiation ("Mathematical model" section). The scheme of ghost imaging utilized in our model is a classical pseudo thermal light ghost imaging setup shown in Fig. 6 (see "Methods" section). In this scheme THz radiation is incident onto the rotating plate with random rough surface on its back side, thus representing a random phase screen, and it forms speckle-picture.
In this paper the ghost image reconstruction is considered for both monochromatic and broadband pulsed THz radiation. For the monochromatic case we use frequency-dependent speckle patterns I r (x, y, ν) = E r (x 2 , y 2 , z = �) 2 from the Eq. (11), and for the case of the pulsed THz radiation we use speckle patterns described by Eq. (5) by taking into account all spectral components of the pulse. For ghost image reconstruction we calculate cross-correlation function 51 between the speckle pattern I r and the intensity transmitted through the object B r , which is measured by a point detector and evaluated by means of the following equation: where T(x, y) is the transmission function of the object. Using a large number of realizations of speckle structures, one can calculate the image estimation function G(x, y) in the following form: Here �·� ≡ 1/M r deotes an ensemble average over M realizations of ghost imaging algorithm.
As an object to be imaged it was used a rectangular ring, shown in Fig. 4a, defined by the inequalities a x < |x| < b x , and a y < |y| < b y , in which the geometric parameters are determined by the following values a x = 4 mm, b x = 6 mm, a y = 3 mm, b y = 5 mm, respectively. Numerical results of ghost image reconstruction for the monochromatic THz radiation are obtained for frequency ν = 0.7 THz and shown in Fig. 4b,d for different number of realizations of the ghost imaging algorithm, M = 10,000 and M = 20,000 , respectively. Parameters of the random phase screen were chosen as follows: h RMS = 1.0 mm, Cl = 1.0 mm, n = 1.55 . The choice of parameters was carried out on the basis of the estimates made in previous section (Statistical characteristics of speckle patterns). It was necessary to obtain such a ratio of parameters at which the speckle diameter (ECD) was close to the minimum and the value of parameter P was high. The distance between the random phase screen and the registration plane is = 5 cm. The signal-to-noise ratio is evaluated in the same manner as it was suggested in 51 . The ghost image reconstruction results for the pulsed THz radiation for the same parameters of the random phase screen were obtained for the Gaussian pulse with duration of 2 ps and spectral profile with central frequency ν 0 = 0.7 THz as shown with blue curve in Fig. 1g. Reconstructed images for different number of realizations of ghost imaging algorithm, M = 10,000 and M = 20,000 , respectively, are presented in Fig. 4c,e. Figure 4c,e demonstrates the fact that speckle patterns formed by broadband THz radiation can be used for ghost imaging technique. Now we can see clearly that the increase of the number of realizations of ghost imaging algorithm in Fig. 4, as one could expect, leads to a growth of the quality of reconstructed image that is also in accordance with rise of SNR 8 . In case of image reconstruction with pulsed THz radiation the values of SNR (9) B r = dxdyI r (x, y, z = �)T(x, y), www.nature.com/scientificreports/ of obtained images are about 10% lower comparing with the corresponding frequency domain results. When comparing monochromatic and pulsed radiation, the same contrast ratio is observed for other central radiation frequencies. It can be explained by the fact that speckle patterns from different spectral components overlap each other. This leads to a decrease of the sharpness of the reconstructed images.
The results presented in this paragraph were made for a speckle pattern correlation coefficient below 0.01. In the experiment we have demonstrated that the correlation coefficient between the speckle patterns is equal to 0.9. For this purpose, in the methods, we showed a comparison of using speckle patterns with correlation coefficient below 0.01 and equal to 0.9. As it was shown by the simulation results, the image reconstruction by the ghost imaging method is also possible even if the correlation coefficient is 0.9. However, in this case, as the correlation coefficient increases, the image quality is getting lower.

Conclusion
In this study, the possibility of experimental formation of speckle patterns using pulsed broadband THz radiation was demonstrated. Numerical modeling was used to describe the conditions for the formation of a speckle structure with the help of a transparent phase object. It was shown that changes in the statistical characteristics of speckles (sizes and contrast of speckles) are significantly determined by the spectrum of the pulsed THz field and the characteristics of the phase object. The characteristics of the phase object include the horizontal correlation length of surface and the root-mean-square deviation of the surface height. We guess that the minimum speckle size is still limited by the maximum frequency in the THz field spectrum. An increase in the correlation length leads to an increase of the speckle size, while saturation occurs at a certain value of the correlation length, after which the speckle size remains approximately constant. Moreover, an increase in the root-mean-square deviation of the surface height leads to the decrease of the speckle size. At the same time, the value of speckle sharpness increases with an increase of the root-mean-square deviation of the surface height due to a strong phase shift for the entire spectrum. Using theoretical estimates, an experimental verification of the formation of speckle patterns for a pulsed THz radiation was carried out. According to the statistical characteristics, the experimentally obtained speckle patterns correspond to the numerically calculated results. Experimental verification demonstrates, in principle, the possibility of forming speckle patterns using broadband THz radiation. Moreover, the work has numerically demonstrated the opportunity of using such speckles in the image restoration system using the ghost imaging method. It must be said that the use of broadband radiation reduced the contrast of the restored images by only 10%. But at the same time, a broadband radiation will allow in the future to restore not only the image of an object itself, it can also control its spectral properties. This feature is very widely used in various applications of science and technology.

Methods of THz radiation beam propagation. The famous Fresnel propagation equation is described by:
where x 1 , y 1 are the source plane coordinates, x 2 , y 2 are the observation plane coordinates, is the distance between source and observation plane (distance of the beam propagation). Here E r (x 1 , y 1 ) and E r (x 2 , y 2 ) is the field distribution in the source and observation plane, respectively. The integral [main text Eq. (3)] was evaluated numerically by means of FFT algorithm.
In the second method we utilized angular spectrum approach (plane wave decomposition 1 ).
where ν 1 = (ν x 1 , ν y 1 ) is a spatial frequency of source plane, F and F −1 are direct and inverse two dimensional Fourier transforms, respectively, H is the transfer function of free-space propagation Experimental parameters of THz radiation. The timeform of THz pulse and its spectrum are shown in Fig. 5. The duration of a THz pulse is 1 ps, and the spectrum range is from 0.05 to 2 THz, pulse energy 300 nJ, the spatial dimension of the Gaussian profile at half-height (FWHM) is 5 mm (Shown Fig. 3c1).
Model of ghost imaging for pulsed THz radiation. The scheme of ghost imaging has a classical pseudo thermal light ghost imaging setup. Here THz pulses are incident onto the plate with changing random rough surface on its back side, representing a random phase screen. The random phase screen introduces a phase shift which subsequently leads to the appearance of random intensity distributions (speckle patterns). After that the beam is divided into two parts by beam splitter. The first beam interacts with the object and propagates to the single pixel detector. The second beam is recorded by multi pixel detector without interacting with the object. Then the cross correlation function between speckle patterns registered by multi pixel detector and overall intensity of speckle patterns transmitted through the object is calculated 51 (Fig. 6). Formation of a large number of speckle patterns can be organized using one or two rotating discs, as in work 28 . At the same time, in contrast to this article, which used monochromatic THz radiation and a metal disk creating amplitude modulation of the transmitted radiation, our study demonstrates the possibility of forming speckle patterns using broadband THz radiation by propagating through transparent stochastic phase screens. Another www.nature.com/scientificreports/ alternative option for arranging spatially inhomogeneous modulation for the formation of speckle structures is the possibility of using the optical range spatial light modulator (SLM) to pump beam modulation which generate THz radiation in crystal 35 as well as modulation of the probe beam 36 , or fabricating SLM for monochromatic THz radiation based on metamaterials 52,53 . Also, up today, for THz radiation, the use of liquid crystal cells as spatially homogeneous phase shifters is being considered 54-56 . Comparing ghost imaging for different correlation coefficients between speckle patterns. It is worth noting that the correlation coefficient between experimentally obtained speckle patterns in Fig. 3c2,c3 is approximately equal to 0.9 as they look relatively similar to each other. However, using the numerical model described above, we demonstrate that ghost image reconstruction for the case with the correlation coefficient between two speckle patterns being equal to 0.9 still provides reasonable results. Thus, in Fig. 7 for comparison two reconstructed ghost images are shown with correlation coefficients equal to 0.01 and 0.9, respectively, for . Model of ghost imaging for pulsed THz radiation. THz pulses are incident onto the plate with changing random rough surface (PP) on its back side, representing a random phase screen. The random phase screen introduces a phase shift which subsequently leads to the appearance of random speckle patterns. After that beam is divided into two parts by beam splitter (BS), the first beam interacts with the object (O) and propagates to the single pixel detector (SPD). The second beam is recorded by CCD camera without interacting with the object Figure 7. Results of a numerical experiment for ghost imaging. Reference object (a), reconstructed images from data with a wide spectrum for 10 thousand realizations using a speckle pattern with correlation coefficients equal to 0.01 (b) and 0.9 (c)