Visualization of internal structure and internal stress in visibly opaque objects using full-field phase-shifting terahertz digital holography

: We construct a full-ﬁeld phase-shifting terahertz digital holography (PS-THz-DH) system by use of a THz quantum cascade laser and an uncooled, 2D micro-bolometer array. The PS-THz-DH enables us to separate the necessary diﬀraction-order image from unnecessary diﬀraction-order images without the need for spatial Fourier ﬁltering, leading to suppress the decrease of spatial resolution. 3D shape of a visibly opaque object is visualized with a sub-millimeter lateral resolution and a sub-µm axial resolution. Also, the digital focusing of amplitude image enables the visualization of internal structure with the millimeter-order axial selectivity. Furthermore, the internal stress distribution of an externally compressed object is visualized from the phase image. The demonstrated results imply a possibility for non-destructive inspection of visibly opaque non-metal materials.


Introduction
Terahertz (THz) imaging [1][2][3] has attracted attention for three-dimensional (3D) visualization of visibly opaque non-metal materials such as plastics, ceramics, rubbers, food, and so on, because THz radiation benefits from good immunity to optical scattering, good penetration to visibly opaque objects, and low invasiveness to human body. One promising method for 3D THz imaging is THz reflection tomography with pulsed THz radiation [4,5], in which the depth information is obtained by a time-of-flight measurement of the THz pulse echo. However, the requirement for mechanical time-delay scanning and mechanical sample-position scanning often hampers its practical use due to slow image acquisition rate. Although the line-field THz reflection tomography boosted the image acquisition rate up to the real-time cross-sectional imaging of the moving object [6], use of a bulky, complicated, expensive, femtosecond regenerative amplifier limits its wide use. THz computed tomography (THz-CT) [7][8][9] is another interesting approach for 3D THz imaging, in which a series of transmitted images are measured at different projection angles and the internal structure of the object is reconstructed by analyzing these images with specific reconstruction algorithms. However, since it is based on the point-to-point measurement with a single THz detector, it is hampered by the low imaging acquisition rate resulting from multi-dimensional mechanical scanning of a sample. Although the line-field THz-CT [10,11] was proposed similar to the line-field THz refection tomography, use of femtosecond regenerative amplifier is still limiting factors.
Recently, THz digital holography (THz-DH) has emerged as a new mode for 3D THz imaging because the recent progress of continuous-wave (CW) THz laser sources and THz digital imagers allows full-field DH to be easily performed in the THz region [12][13][14][15][16][17]. In usual DH, an interference image is generated by interference between a diffracted beam in an object and a reference one, and then is recorded as a hologram by a digital imaging device. By performing a diffraction calculation of the acquired hologram with a computer, the amplitude and phase images of the object can be reconstructed at an arbitrary plane. In particular, the phase image can be used for 3D imaging of geometrical shapes or optical thicknesses in an object with a sub-wavelength axial resolution. Combination of DH with THz radiation, namely THz-DH, enables us to expand the application of DH to visually opaque objects or objects covered with visibly rough surface. As an early work of THz-DH, Gabor-type THz-DH has been demonstrated for phase imaging of objects at a transmission configuration [12][13][14]; while benefitting from a simple optical configuration without the need for a two-arm interferometer, the necessary first-diffraction-order image is spatially overlapped on unnecessary diffraction-order images, such as zero-diffraction-order and conjugate first-diffraction-order, in spatial-frequency domain, leading to the degraded quality of the reconstructed image. To overcome this problem, off-axis THz-DH (OA-THz-DH) [15][16][17] was proposed, enabling us to spatially separate the necessary first-diffraction-order image from unnecessary diffraction-order images in spatial-frequency domain by Fourier filtering. However, such spatial Fourier filtering decreases the spatial resolution and the number of pixels in the reconstructed image.
One possible approach to extract the necessary first-diffraction-order image without the need for spatial Fourier filtering is use of phase-shifting method in DH, namely phase-shifting DH or PS-DH, which is widely used in the visible region [18][19][20]. In PS-DH, a few holograms with known different phase shifts are used to reconstruct the complex amplitude of the object beam of an object without the decrease of pixels. Effectiveness of PS-DH has been demonstrated in THz region as PS-THz-DH [21]; however, this frontier work of PS-THz-DH was based on sample-scanning acquisition of interferogram with a THz point-detector, hampering the long image acquisition time. More recently, full-field PS-THz-DH was achieved by a CO 2 -laserpumped CW-THz laser; however, versality of full-field PS-THz-DH was hampered by a bulky and complicated laser [22]. There are no attempts to perform full-field PS-THz-DH with a combination of a compact CW-THz source and a THz imaging device.
In this article, we construct full-field PS-THz-DH equipped with a THz quantum cascade laser (THz-QCL) [23,24] and an uncooled, 2D micro-bolometer array (THz imager) [25]. The proposed full-field PS-THz-DH is applied for visualization of internal structure and internal stress of visibly opaque objects. Furthermore, we discuss limitation factors of spatial resolution in PS-THz-DH, OA-THz-DH, and visible PS-DH.

Principle of operation
In PS-DH, after acquiring three or more holograms with known phase difference and then calculating the complex amplitude of object beam from them, we obtain the amplitude and phase images of an object at an arbitrary place by the wavefront propagation calculation (diffraction integral) of the complex amplitude. In this article, the four-step phase-shifting method was selected from the viewpoint of balance between easy implementation and measurement precision [26]. In this method, four holograms were acquired while shifting the phase in a step-by-step manner (0 rad, π/2 rad, π rad, and 3π/2 rad). We here define (x, y) coordinate as a lateral plane of a sample and z coordinate as an axial direction of the sample corresponding to a propagation direction of THz radiation. The complex amplitude distribution of object beam and reference beam, u o (x, y) and u r (x, y), is given by where a o (x, y) and a r (x, y) are real amplitude distribution of object beam and reference beam and φ o (x, y) and φ r (x, y) are phase distribution of them. When the phase of the reference beam is shifted by δ, the corresponding hologram I(x, y;δ) is given by where u o *(x, y) and u r *(x, y) are conjugates of u o (x, y) and u r (x, y). For example, when δ is set to 0 rad, π/2 rad, π rad, and 3π/2 rad, the corresponding four holograms are given by From Eqs. (4), (5), (6), and (7), u o (x, y) at the detector plane can be obtained by where the reference beam is assumed to be a plane wave. Assuming that the reference beam is a plane wave, Eq. (8) can be rewritten as To obtain u o (x, y) at an arbitrary plane, we used the angular spectrum method [27] together with the following equation where U o (f x ,f y ;0) is Fourier spectrum of u o (x,y;0). Figure 1 shows a schematic diagram of a transmission PS-THz-DH system. A continuouswave (CW) THz-QCL (EasyQCL, LongWave Photonics LLC, center frequency = 3 THz, center wavelength = 100 µm, average power = 1.78 mW, operation temperature = 48 K) was used for CW-THz radiation. The emitted CW-THz radiation was collimated with a diameter of 6 mm by an off-axis parabolic mirror (OA-PM, off-axis angle = 90°, diameter = 25.4 mm, focal length = 50.8 mm) and passes through an optical chopper (OC, 3501, New Focus, modulation frequency = 7.5 Hz). Then, the CW-THz radiation was split into an object beam and a reference beam of a Mach-Zehnder interferometer by a silicon beam splitter (BS, diameter = 101.6 mm, thickness = 500 µm, reflection = 54%, transmittance = 46%). A sample was placed at an optical path of the object beam. An optical path of the refence beam includes a combination of a retro reflector (RR) and a stepping-motor-driven translation stage (KS261-85, SURUGA SEIKI CO., LTD., travel range = 85 mm, incremental motion = 1 µm, accuracy = ±30 µm for travel range of 85 mm, uni-directional repeatability = ± 3 µm) for stage translation (0 µm, 25 µm, 50 µm, and 75 µm), corresponding to four-step phase shift (0 rad, π/2 rad, π rad, and 3π/2 rad). The object beam and the reference beam were collinearly overlapped again by another BS to form the hologram. The formed hologram was acquired by an uncooled, 2D micro-bolometer array (THz imager, IRV-T0831, NEC Inc., sensor area = 7.52 mm × 5.64 mm, pixel number = 320 pixels × 240 pixels, pixel size = 23.5 µm × 23.5 µm, frame rate = 7.5 frames per second, exposure time = 66.7 ms, digital output resolution = 14 bit). To improve the image SNR, we adopted a dynamic subtraction technique [17], in which the chopped CW-THz radiation by OC was synchronized with the frame timing of the THz imager.

Basic performance
We first performed amplitude and phase imaging of a visibly opaque object, in which two letters, "P" and "S", were impressed in a convex manner (height = 228 µm) on a polystyrene plate (thickness = 0.8 mm, refractive index = 1.67), as shown in Fig. 2(a). The measured region is indicated as a red box in Fig. 2(a). The acquired holograms I 0 (x, y), I π/2 (x, y), I π (x, y), and I 3π/2 (x, y) are shown in Figs. 2(b), 2(c), 2(d), and 2(e), respectively (image size = 7.52 mm × 5.64 mm, pixel number = 320 pixels × 240 pixels, number of image integrations = 128). From these holograms, we reconstructed the amplitude and phase images of the sample at a reconstruction distance z of 31.6 mm as shown in Figs. 2(f) and 2(g) (image size = 9.49 mm × 9.49 mm). The reason why the image size is larger than the sensor area size in the THz imager is due to the extended calculation region at the object plane by use of null data padding in Fourier transform in the reconstructed process. The letter 'P' was clearly visualized in the amplitude contrast and the phase contrast by four-step PS-THz-DH. The phase image contrast in the PS-THz-DH [see Fig. 2(g)] was lower than that in the previous OA-THz-DH [17] because the image SNR in the former equipped with two BSs was lower than that in the latter equipped with one BS. (a) Optical photograph of a visibly opaque, plastic sample with two impressed letters, "P" and "S". THz digital holograms with a phase-shift of (b) 0 rad, (c) π/2 rad, (d) π rad, and (e) 3π/2 rad. Reconstructed images of (f) amplitude and (g) phase.
We next evaluated lateral resolution of amplitude imaging in PS-THz-DH. To this end, we prepared a sufficiently fine, plastic line (diameter = 330 µm) for sample and placed it as intersecting in a cross shape. After reconstructing the amplitude image of this sample, we obtained the intensity image by squaring the amplitude image, as shown in Fig. 3(a) (image size = 9.49 mm × 9.49 mm, number of image integrations = 128). Then, we extracted the intensity profile across horizontal and vertical lines of the cross shape [see blue lines in Fig. 3

(a)]. Figures 3(b) and 3(c)
show intensity profiles across vertical and horizontal lines for the intensity image, respectively. Then, we calculated the half width at 1-1/e 2 of difference of the maximum and minimum intensity values in them by curve fitting analysis with a Gaussian function: 657 µm and 842 µm for horizontal and vertical dimensions in amplitude image. These values were in good agreement with theoretical value for lateral resolution, LR x and LR y , given by [28] where λ is the wavelength of the CW-THz radiation ( = 100 µm), z is the reconstruction distance, M and N are the number of pixels along the horizontal and vertical directions of the THz imager, ∆ x and ∆ y are the size of pixel along the horizontal and vertical directions of the THz imager. Such the diffraction-limited spatial resolution is similar to that in the previous research in OA-THz-DH [17]. For more detailed evaluation of spatial resolution, 'Siemens star' target will be a powerful tool [29]. When the phase image is used for 3D shape measurement, the phase resolution, defined as the fluctuation of the phase value at each pixel between different phase images, determines the axial resolution in 3D shape measurement. To evaluate the phase resolution of PS-THz-DH, we repeated 10 acquisitions of phase image without a sample (number of image integrations = 128). Then, we calculated the standard deviation of the phase values at each pixel in 10 phase images to provide the phase resolution. Figure 4 shows the spatial distribution of the phase resolution (image size = 7.52 mm × 5.64 mm, pixel number = 320 pixels × 240 pixels). The spatial dependence of the phase resolution was relatively flat. The mean of the phase resolution in a region of interest (ROI) shown by the yellow square region of Fig. 4 (size = 2.2 mm × 2.2 mm) was 0.041 rad, corresponding to a phase resolving power of λ/151. From these values, the axial resolution in the 3D shape measurement was estimated to be 0.66 µm for the ROI.

3D thickness measurement of visibly opaque object
We performed 3D thickness imaging of a visibly opaque, semiconductor sample with known dimension. Figure 5(a) shows an optical photograph of this sample, in which a checker pattern with 1-mm square was carved within a size of 20 mm by 20 mm on a silicon substrate (n-type, thickness = 300 µm, electric resistivity > 1 Ω•cm) by a photolithography method. Black square and white square in an inset of Fig. 5(a) respectively show convex surface and concave surface, corresponding to top surface and bottom surface. Difference between convex surface and concave surface was beforehand determined to be 21.24 ± 0.02 µm using 3D optical profilometer (Profilm 3D, Filmetrics, Inc., reflection configuration, axial resolution = 0.05 µm). Figures 5(b) and 5(c) show the amplitude and phase images of this sample reconstructed at z of 15 mm (image size = 9.49 mm × 9.49 mm, number of image integrations = 128). The checker pattern was visualized with different image contrast in these images. Since the contrast of the amplitude image is given by attenuation of CW-THz radiation, scattering at edges of the checker pattern results in high contrast. A little difference of image contrast was confirmed between convex surface and concave surface, which is mainly due to difference of surface roughness between optically-polished convex surface and caved concave surface rather than difference of absorption between thicker convex surface and thinner concave surface. On the other hand, the phase image gives high contrast between convex surface and concave surface because its image contrast depends on the optical phase and hence the optical thickness between them. This phase image enables us to obtain 2D distribution of the relative thickness t(x, y) for this sample under the assumption of a known refractive index. t(x, y) is calculated from the phase image φ(x, y) by where λ is the wavelength of the CW-THz radiation ( = 100 µm), n Si is the refractive index of silicon ( = 3.1), and n air is the refractive index of air ( = 1.0). Figure 5(d) shows the calculated t(x, y) for ROI as shown by a red box region in Fig. 5(c) (size = 4 mm × 4 mm). Surface profile of this sample was clearly visualized as the caved checker pattern. From this surface profile, we determined the difference between convex surface and concave surface to be 21 ± 1 µm, which is in good agreement with that ( = 21.24 ± 0.02 µm) by using 3D optical profilometer. This comparison indicates high potential of phase imaging in PS-THz-DH for 3D thickness measurement. We next evaluated the robustness of amplitude and phase imaging to the low SNR. To this end, we obtained a series of the amplitude and phase images of the same checker-pattern sample with different THz power incident onto the object as shown in Figs. 6(a) and 6(b). The amplitude image decreased its contrast with respect to the reduced power of the object beam, and almost disappeared at 26% power. On the other hand, the phase images maintained the same contrast value regardless of the reduced power. This comparison clearly indicated that the phase imaging is more robust to the low SNR than the amplitude imaging. To investigate difference of low-SNR robustness, we measured the mean value for the selected ROI (size = 1.17 mm × 1.17 mm) of amplitude image under no samples with respect to the normalized THz power, as shown in Fig.  6(c); the mean amplitude is proportional to the THz power. We next measured the spatial phase noise for the same ROI with respect to with respect to the normalized THz power, as shown in Fig. 6(d); the spatial phase noise maintains the same level over a certain threshold of THz power. From these results, we consider that difference of the low-SNR robustness between amplitude imaging and phase one is due to different contributions of signal and noise to the reduced THz power in amplitude and phase.

Visualization of internal structure of visibly opaque object
Visualization of internal structure of a visibly opaque object is another interesting imaging modality in PS-THz-DH. Although reconstruction of a target object through visibly opaque, nonconducting material was demonstrated based on the phase image [16,30], we here demonstrated visualization of internal structure with the axial selectivity by digital focusing of amplitude image. The digital focusing enables 3D imaging with an axial resolution over a wide axial range without changing the optical system, in which a series of amplitude images are acquired while changing the reconstruction distance z. Figure 7(a) shows an optical photograph of a sample. We pierced the foamed polystyrene block with two metal nails (diameter = 1.25 mm) orthogonal to each other with a separation of 13 mm. Then, we reconstructed the amplitude images at a 1-mm interval of z. Visualization 1 shows a series of amplitude images with respect to z. Image contrast of the vertical and horizontal nails changes in different way with respect to z due to different axial positions of them. Figures 7(b) and 7(c) show amplitude images of the sample with best focusing at different z positions: 36 mm and 49 mm. The horizontal nail was visualized in focus at 36 mm and out of focus at 49 mm; on the other hand, the vertical nails show the opposite behavior in these images. To evaluate the out of focus quantitatively, we extracted the line profile across the vertical and horizontal nails in these images and determined the diameter of the line profile.

Visualization of internal stress distribution of externally compressed, visibly opaque object
Another potential application of the PS-THz-DH-based phase imaging is in visualization of internal or residual stress distribution of visibly opaque objects because it is often mentioned in failure analysis as a contributing factor in fracture. Strain gauge is widely used for stress measurement due to its cost effectiveness, however, it does not meet the requirement for non-contact measurement; optical polarimetry [31] enables us to visualize internal stress in a non-contact manner, however, its application is limited to transparent objects in the visible region. We here applied the PS-THz-DH-based phase imaging for visualization of internal stress distribution of a visibly opaque object caused by external force. A sample is a cube-block piece of high-density white polyethylene (size = 36-mm wide × 20-mm height × 10-mm depth, corresponding to x-y-z coordinate in Fig. 1, refractive index = 1.54, Young's modulus = 1.3 GPa, compressive yield stress = 22 MPa, Poisson's ratio = 0.3). External compression force was applied to the sample along its height direction parallel to y axis in Fig. 1 within the range of elastic deformation by a small hydraulic press machine (GS03940, Specac, Inc., maximum load = 2 ton). MPa was applied to the sample, respectively. The differential phase images ∆φ(x, y) were obtained by calculating difference between phase images with and without the compression stress. Then, we calculated the mean of the phase difference for ROI shown by the black square region of Figs. 8(a) to 8(d) (size = 5.9 mm × 5.9 mm), as shown in Fig.  8(e). This graph clearly indicates the dependence of phase difference on the compression stress. The slope coefficient between them was 0.0616 rad/MPa although such non-linear behavior requires further investigation for experimental methodology and material property. We next calculated the internal stress of this compressed sample assuming that its geometrical dimension (∆x, ∆y, and ∆z) was changed without the change of its refractive index. In the elastic deformation of isotropic material, the compression stress along the y axis (σ y ) induces its orthogonal strain along x-axis and z-axis (ε x and ε z ) in addition to its parallel strain along y-axis (ε y ). The differential phase image obtained by PS-THz-DH is sensitive to ∆z due to the change of the geometrical thickness of the sample. The strain distribution ε z (x, y) is given by where z is the thickness of the sample, n wp is the refractive index of white polyethylene, and ∆φ(x, y) is the phase difference distribution. Using the Poisson's ratio, the strain distribution ε y (x, y) is given by where ν is the Poisson's ratio of the sample material. From Hooke's Law, the internal stress distribution σ y (x, y) is given by where E is the Young's modulus of the sample material. Figure 8(f) shows the internal stress distribution σ y (x, y) under the externally applied compression stress of 7.35 MPa, which was  Fig. 4, the measurement resolution of internal stress was estimated to be 0.667 MPa. Such visualization of internal stress in visibly opaque object will be useful for production process of a bulk material and/or machining process of a product in the field of soft materials.

Discussion
We first discuss the spatial resolution in PS-DH between the visible region and THz region. Although the spatial resolution in DH is given in detail elsewhere [32], it is important to note that a size relation between wavelength and camera pixel size is different between THz and visible regions: wavelength < pixel size in visible region and wavelength > pixel size in THz region. Figures 9(a) and 9(b) show the schematic drawing of spatial resolution for visible PS-DH and PS-THz-DH in spatial frequency domain, respectively. Blue square region indicates the available region of spatial frequency determined by the camera pixel size ∆ x and ∆ y . On the other hand, red circle region indicates the spatial-frequency spectrum of the diffraction-order image, whose radius depends on the numerical aperture (NA) of optical system and wavelength λ. In the visible PS-DH, since the blue square region is smaller than the red circle region, the spatial resolution is determined by the camera pixel size [see Fig. 9(a)]. In contrast, in the PS-THz-DH, since the red circle region is smaller than the blue square region, the spatial resolution is determined by the NA/λ [see Fig. 9(b)]. In this way, the spatial resolution of PS-DH was limited by different factors in the visible region and THz region due to difference of the size relation between wavelength and camera pixel size. We next discuss the difference of spatial resolution between PS-THz-DH and OA-THz-DH. As shown in Fig. 9(b), the spatial resolution of PS-THz-DH is limited by NA/λ rather than the camera pixel size. In the case of OA-THz-DH, one has to further consider spatial frequency of interference fringe, size of the spatial-frequency spectra of diffraction and non-diffraction order images, and their spectral overlapping. Figure 10 shows a schematic drawing of the spatial-frequency spectrum when (a) 2NA ≥ sinθ, (b) 2NA < sinθ < 3NA, and (c) 3NA ≤ sinθ, respectively, where θ is the off-axis angle of two beams in the OA-THz-DH. The green circle indicates spatial-frequency spectrum of the zero-diffraction-order image, whose radius is given by 2NA/λ. The red and purple circles indicate spatial-frequency spectra of the first-diffraction-order image and the conjugate first-diffraction-order image, respectively; their radius is given by NA/λ whereas distance of their center from the spatial-frequency origin is given by sinθ/λ, corresponding to the spatial frequency of interference fringe in OA-THz-DH. When 2NA = sinθ, the spectral center of the first-diffraction-order image (red circle) is located on the spectral periphery of the zero-diffractionorder image (green circle). On the other hand, when 3NA = sinθ, the spectral periphery of the first-diffraction-order image circumscribes that of the zero-diffraction-order image. In Fig. 10(a), the center of the red circle is inside of the green circle; it is impossible to extract only the spectrum of the first-diffraction-order image by the spatial Fourier filtering because the spatial Fourier filtering is performed around the spectral center of the first-diffraction-order image and always includes the spectrum of the zero-diffraction-order image. In Fig. 10(b), the center of the red circle is outside of the green circle; the spatial Fourier filtering enables us to separate the spectrum of the first-diffraction-order image from that of the zero-diffraction-order image although the partially spectral overlapping of them reduces the size of the separated spectrum for the first-diffraction-order image. For example, if NA of the OA-THz-DH is equal to that of the PS-THz-DH, the spatial resolution of the PS-THz-DH is improved by NA/(sinθ -2NA) from that of the OA-THz-DH. In Fig. 10(c), the whole red circle is outside of the green circle; the whole spectrum of the first-diffraction-order image can be extracted without the spectral interference of the zero-diffraction-order image. If NA of the OA-THz-DH is equal to that of the PS-THz-DH, the spatial resolution of the former is also equal to that of the latter. However, the reduced spectral area of the first-diffraction-order image to avoid the spectral overlapping leads to the lower spatial resolution. We simulated the amplitude images of vertically stripe lines spaced with 235 µm interval [see Fig. 10(d)] for the PS-THz-DH with NA = 0.4 and for the OA-THz-DH with θ = 60°and NA = 0.4, corresponding to the case of Fig. 10(b). The resulting amplitude images in the PS-THz-DH and the OA-THz-DH are shown in Figs. 10(e) and 10(f). In this case, the spatial resolution of the PS-THz-DH is 6-times better than that of the OA-THz-DH. In this way, PS-THz-DH has the advantage of spatial resolution over the OA-THz-DH.

Conclusions
We constructed full-field THz-DH with the four-step phase-shifting method. Use of a 3-THz CW-THz-QCL and an uncooled THz imager enables us to simplify its optical setup. Using the proposed setup, 3D shape of visibly opaque objects was visualized as amplitude and phase images; the lateral resolution of amplitude image was 657 µm and 842 µm for the horizontal and vertical dimensions whereas the phase resolving power of phase image was λ/151 corresponding to the axial resolution of 0.66 µm. Furthermore, we demonstrated two potential applications of the full-field phase-shifting THz-DH for visibly opaque objects: visualization of internal structure with the millimeter-order axial selectivity and mapping of internal stress distribution with MPa order. The full-field phase-shifting THz-DH method will be a powerful tool for non-destructive inspection of visibly opaque materials.

Funding
Exploratory Research for Advanced Technology (JPMJER1304); Japan Society for the Promotion of Science (18K04981); Cabinet Office, Government of Japan (Subsidy for Reg. Univ. and Reg. Ind. Creation).