1. Introduction

A heterodyne receiver, consisting of a mixer and a local oscillator (LO), is of crucial importance for astronomical observation and atmospheric remote sensing in the submillimeter wavelength range, due to its high spectral resolution (ν/Δν>10⁶) and quantum limited sensitivity. Nowadays, multi-beam heterodyne receivers play a vital role in single dish astronomical observations, due to their improved mapping speed with an enhanced number of pixels (e.g. SuperCam, upGREAT, GUSTO [1–3]). One of the major challenges in large format heterodyne array receivers lies in the difficulty of achieving a high efficiency LO power distribution system. Approaching the submillimeter wavelength range, due to the difficulty in micro-fabrication of metallic waveguide structures, phase grating was developed as the free space coupling geometry to split the incoming beam into equally distributed multiple beams for pumping multiple mixers. Over the past few years, pseudo-2D Fourier phase gratings have undergone immense progress. Based on the continuous phase modulation, Fourier phase gratings generate a smooth surface profile, and have been demonstrated from 345 GHz up to 4.7 THz [4–8].

For next generation large format THz multi-beam heterodyne array receivers, two-dimensional (2D) mixer arrangement is necessary. Typical 2D phase gratings are designed by the superposition of two orthogonal one-dimensional gratings, which are regarded as pseudo two-dimensional gratings. Although the pseudo-2D grating is easy to design, it also has several major drawbacks. Firstly, some atypical beam arrangements (such as the hexagonal configuration used in the upGREAT spectrometer onboard of SOFIA [2]) cannot be generated in this pseudo way, since it can only form standard rectangular beam arrangements (such as 2×2 or 4×4 patterns). Secondly, the superposition of two orthogonal one-dimensional gratings creates a certain correlation between the phase of the diffracted beams, which is unnecessary for many applications. More importantly, such dispensable constraint in phase will restrain the maximum diffraction efficiency of the grating been obtained.

In contrast, true two-dimensional gratings rely on a complete two-dimensional phase retrieval design approach, which in principle can generate arbitrary beam profiles and release unnecessary phase constraints between diffracted beams. Therefore, as demonstrated in [9], such true two-dimensional grating brings significantly higher efficiency than conventional pseudo two-dimensional gratings. The design of a phase grating lies in solving the phase retrieval problem. The most commonly used method is the Gerchberg-Saxton (GS) algorithm, which holds the advantage of achieving a predetermined complex distribution with high diffraction efficiency [10]. However, one major limitation of using the GS algorithm for phase grating design is to achieve high intensity accuracy of the beam distribution, especially in the presence of a large number of diffracted beams.

In this study, a 2×2 beams true two-dimensional grating at 0.85 THz with 20 degree incident angle of the LO beam, has been realized. The far-field beam pattern has been measured, where a measured diffraction power efficiency of 81.9% confirmed good agreement with the simulation results. Furthermore, we report the design of a 100-pixel two-dimensional phase grating at terahertz wavelengths, based on the Mixed-Region-Amplitude-Freedom (MRAF) algorithm [11]. Rigorous full wave simulation proves the efficiency and accuracy of the design, where an improvement by a factor of four in accuracy has been obtained compared with the standard GS scheme.

2. Gerchberg-Saxton algorithm design procedure

The Gerchberg-Saxton algorithm is an efficient solution for retrieving the phase distributions of light at the hologram plane and the image plane. In our case, it is applied to find the phase distribution of the reflection grating, which transforms the given input intensity distribution at the hologram plane (the grating surface) into the desired intensity distribution at the image plane (the reflected far-field plane), as indicated in Fig. 1. In our case, the field at the image plane can be treated as the Fourier transformation of the field at the hologram plane.

 

Fig. 1. Optical scheme of the generation of the desired intensity distributions at the image plane, obtained by Fourier transform of the optical field at the hologram plane.

Download Full Size | PDF

The first step of the design is to treat the unit cell of the grating surface (hologram plane) as an N×N segments grid (N=4/8/16/32). Therefore, the image plane contains the diffraction orders between -N/2 to N/2 in each dimension. Later, at the image plane, the amplitude of the desired diffraction orders are set to A₀ and all other unwanted diffraction beams are set to zero. In the optimization process of the GS algorithm, this diffraction field distribution is set as the target of the image plane. The algorithm uses the defined amplitude at the hologram plane and the image plane to alternate between the two-dimensional Fourier transform and the inverse Fourier transform at each iteration to improve the phase estimation. For the initial iteration step, a random phase is applied to all diffracted beams at the image plane, and no symmetry phase constraint at the hologram plane is applied in the algorithm.

For practical applications, oblique angle of the incident beam is considered. The GS algorithm optimization only considers plane wave incident situation, while at certain oblique angle, the phase front of the beam will create additional optical path difference at different positions of the grating. Therefore, it is necessary to perform phase correction depending on the angle θ of the incident beam, to compensate the optical path difference. Assuming the incident beam lives in the plane xz, the dimension of the phase grating has to be modified along the x and z directions of the incident beam. The modification is performed with [8]:

3. 2×2 beams phase grating

In our work, a true two-dimensional grating operating at 0.85 THz with 2×2 beams distribution is considered, where the target far field beam pattern is composed of four diffraction beams with identical intensity. The GS algorithm is employed to obtain the phase distribution of the grating, which is defined as a grid of N×N segments. In principle, larger number of segments for each phase grating unit leads to higher diffraction efficiency, but it also brings more difficulties in simulation and manufacture. In our case, a group of scenarios with 4×4, 8×8, 16×16 and 32×32 segments for one phase grating unit are considered. As shown in Fig. 1, the phase distribution of one grating unit with different segments is obtained from the GS algorithm. Subsequently, based on the calculated phase profile matrix, a model in FEKO consisting of nine grating units (3×3) is developed. For the most complex case (32×32 segments), the entire design and simulation process takes about 5 hours on a 3.5 GHz quad-core CPU, 40 GB RAM, and 16 GB GPU desktop computer.

Figure 2 summarizes the simulated far field beam patterns of different segments, where the color of the beam pattern is on a logarithmic scale. In all cases, four main diffraction beams are obtained. Obviously, with 4×4 segments, additional higher order diffraction beams are presented. As the number of segments increases, the unwanted diffraction orders are less pronounced. The beam intensity is calculated by integrating the power in each divergent beam with a noise threshold level setting to 1% of the maximum intensity, while the total power is set as the power integrated over the full hemispherical plane. The total power efficiency is calculated as the sum of the power of the four beams. For the 4×4, 8×8, 16×16 and 32×32 segments design in FEKO, they are 73.2%, 82.3%, 83.9% and 84.6%, respectively. For comparison, the total power efficiency of a typical pseudo-2D Fourier grating with 2×2 beams distribution is 65% [6], which reflects 10-20% improvement for the true two-dimensional phase gratings at least. Obviously, larger number of segments with finer resolved phase hologram will lead to even higher diffraction efficiency, which is limited by our computing power.

 

Fig. 2. Phase distribution (upper panel) and the simulated far field beam patterns (lower panel) of two-dimensional grating with 4×4 (a,e), 8×8 (b,f), 16×16 (c,g) and 32×32 (d,h) segments. The color of the phase distribution is on a linear scale from –π to π, while the intensity of the beam pattern is on a logarithmic scale.

Download Full Size | PDF

To verify the design and simulations, we fabricated a prototype true two-dimensional phase grating. Due to the compromise between the diffraction efficiency and the difficulty in the physical realization of the grating, an 8×8 segments grid is selected. As shown in Fig. 3, the grating is composed of 9 unit cells (3×3), and the unit cell length is 4 mm. In this case, 20 degree incident angle is considered. As explained, the phase profile matrix is obtained from the GS algorithm and tuned according to the FEKO simulation results and the angle compensation formulas. Overall, the size of the manufactured grating is 12.8 mm × 12 mm. The entire grating structure takes about 130 hours to be processed, by a computer controlled micro-milling machine, where the short flute end mill with a diameter of 100 μm is employed. A laser confocal microscope (OLS5000 from Olympus Corporation) is used to measure the surface profile of the manufactured phase grating, which is presented in Fig. 3. The height of the periodic segments grid of the grating is displayed in linear color scale. The designed and manufactured surface profile of the grating in two orthogonal directions are also drawn. The maximum discrepancy between the two curves is about 5 μm, which results in a phase difference of less than 6°. The average roughness of the structured surface is less than 2 μm.

 

Fig. 3. (a) The far-field beam pattern experimental setup for characterization of the two-dimensional phase grating. (b) The picture of the manufactured grating composing 3×3 unit cells, shown with the lengths in two orthogonal directions of the grating. (c) The surface profile of the manufactured phase grating measured with a laser confocal microscope, and the height in μm is shown on a linear scale. The left and upper panels show the comparison between the designed profile (grey curve) and manufactured profile (blue/red curve) in two orthogonal directions. Also shown are the discrepancies between the designed and manufactured profiles (black curves) in two directions which are enlarged by a factor of 10.

Download Full Size | PDF

For the measurement, as indicated in Fig. 3, a solid-state multiplier chain is employed to generate a linear polarized incident beam at 0.85 THz. The output power of the LO chain is about 50 μW. The beam is coupled by two off-axis parabolic mirrors, and a Gaussian beam with a beam-width of 6 mm is obtained at the position of the grating surface. The reflected beam is measured by a pyro-electrical detector scanned with two PC-controlled rotation stages in spherical coordinates. A pinhole of 2 mm is located in front of the detector to define the angular resolution of the system. The distance between the detector and the grating is 250 mm. The scanning area is set to 26°×26° to cover four diffraction beams.

Figure 4 represents the experimental results of the far-field beam distribution pattern, at an incident beam angle of 20°. The color of the beam pattern is on a logarithmic scale and the step size of the measurement is 0.1 degree. The beam efficiency is calculated by integrating the power in each divergent beam, and compared to the total power measured by scanning the incident beam at the position of the grating. According to the FEKO simulation, the power distributions of four divergent beams are 19.8%, 20.0%, 20.6% and 20.9%, and the measurement results show that the power efficiency of each beam varies less than 2.2% (17.6%, 19.6%, 21.9% and 22.8%) at this frequency. The total power efficiency is calculated as the sum of the power of the four beams, which is 81.3% in the simulation and 81.9% (±2%) in the measurement. The slightly higher power efficiency in the measurement could be due to the uncertainty caused by the power variations of the LO source during the scanning measurement. Overall, good agreement between measurements and simulations validates the design and the manufacture process. Simulation results demonstrate a bandwidth up to over 350 GHz (750-1000 GHz), with a total power efficiency over 50%. As the operating frequency deviates from the center frequency, not only the total power efficiency drops, but the uniformity of the power distribution among all the four beams also deteriorates. Compared with the standard 2D Fourier phase grating, the measured efficiency is improved by 17%. New manufacture technologies such as femtosecond laser machines and 3D printing will be able to create phase grating with larger number of segments, where even higher diffraction efficiency is feasible.

 

Fig. 4. Measured far-field beam pattern at 0.85 THz with an incident angle of 20°, and the color of the beam pattern is on a logarithmic scale. The beam is scanned in spherical coordinates at a distance of 250 mm from the grating.

Download Full Size | PDF

4. 100 pixels phase grating by the mixed-region-amplitude-freedom algorithm

The hundred-pixel level heterodyne array receiver is essential for the next generation submillimeter radio telescopes to improve the observation efficiency. As an important step, Fourier phase grating generating 81 beams at 3.86 THz has been demonstrated in [12]. By superimposing two orthogonal one-dimensional Fourier gratings with nine diffraction beams, a total power efficiency of 93% is achieved. In general, pseudo-2D Fourier phase grating maintains lower diffraction efficiency with even diffraction beam numbers, due to the additional half-wave phase step required to suppress the 0^th order diffraction beam. For instance, the optimal total power efficiency of a 10×10 beams pseudo-2D grating is around 75%.

For true two-dimensional gratings, the removal of unnecessary phase constraints will further improve the total power efficiency. Taking a 10×10 beams true phase grating as an example, the total power efficiency can reach 96% from simulation results in FEKO. However, as illustrated in Fig. 5 (b), the main challenge in designing a large-pixel 2D phase grating by the GS algorithm is the intensity accuracy of the diffraction beams, which is not an issue for the 2×2 beams distribution case. In order to characterize the accuracy of the intensity distribution, the accuracy metric (σ) is defined as:

 

Fig. 5. 10×10 beams true two-dimensional phase design with the standard GS algorithm (upper panel) and the MRAF algorithm (lower panel), including the phase distribution (a,c), and simulated intensity distribution (b,d). The signal region (SR) is indicated by the white square zone, and the noise region (NR) is the rest of the image plane.

Download Full Size | PDF

To further improve the accuracy metric, the standard GS algorithm has been modified by various methods, including the non-zero padding GS algorithm [13], the Mixed-Region-Amplitude-Freedom (MRAF) algorithm [11], and the offset-MRAF (OMRAF) algorithm [14]. For all these schemes, the image plane can be defined in two subsets, namely the signal region (SR) and the noise region (NR). The idea of the non-zero padding GS algorithm is to treat the noise region as the padding zone, and the padding zone is filled with fixed non-zero values. In contrast, in the MRAF algorithm, the field at the image plane is defined as the combination of the target field (in the SR) with the propagated field from the previous iteration (in the NR). In this configuration, the phase freedom is permitted in the entire region at the image plane, while the amplitude freedom is allowed only in the NR. The core of the MRAF algorithm is to further improve the convergence over the SR by giving up control of the NR, which in turn can provide better performance in terms of accuracy of the intensity distribution. The OMRAF algorithm is built on the MRAF algorithm. By including non-zero background light intensity in the SR, the OMRAF algorithm presents a great advantage for generating arbitrary two-dimensional holographic traps in cold atom experiments [14]. For large-pixel two-dimensional THz phase grating applications, zero background intensity distribution is required. Therefore, in this work, we will focus on the MRAF algorithm.

For THz phase grating applications, the finite length of the grating unit cell defines the number of diffraction orders that can be reflected. By placing the NR beyond the maximum reflection order, the target diffraction pattern can be well preserved. In this work, the segments size (64×64) defines the diffraction orders between -32 and 32 calculated in both the horizontal and the vertical directions at the image plane. The divergent angle between each diffraction order can be defined as $\theta = si{n^{ - 1}}\left( {\frac{\lambda }{D}} \right)$, where λ is the wavelength and D is the length of the unit cell. In our case, a 5° divergent angle is obtained with a unit cell length of 4 mm and a center frequency of 0.85 THz. Consequently, a diffracted beam with a maximum diffraction order of 19^th can be obtained at the entire reflection plane. Therefore, the diffraction field beyond this order is defined as the NR, and the diffraction field within this diffraction order is defined as the SR, as indicated in Fig. 5 (d). The GS algorithm uses the amplitude defined at the hologram plane and the image plane to alternate between the two-dimensional Fourier transform and the inverse Fourier transform to achieve the phase estimation. As described in [11], the MRAF algorithm modifies the field at the image plane as the combination with the target field in the SR and the propagated field from the previous iteration in the NR, according to:

To quantitatively characterize the performance of the MRAF algorithm and the GS algorithm, we evaluate the accuracy metric and the diffraction efficiency. The accuracy metric defines the accuracy of the intensity distribution at the image plane, and the total diffraction efficiency η is defined as:

We examine the performance of the MRAF algorithm with a variation of mixing parameter m. For the case with 10×10 beams diffraction pattern as the target field, the signal region is outlined in Fig. 5, and the rest of the pattern is defined as the noise region. The MRAF algorithm would improve the accuracy metric of the diffraction beams in the SR, by the release of the amplitude freedom in the NR. Thus, in general as plotted in Fig. 6, lower m value defines higher power ratio in the NR, which in turn achieves better accuracy of the intensity distribution but lower total diffraction efficiency. For the best scenario, the MRAF algorithm on average outperforms the standard GS algorithm by at least a factor of four in accuracy (σ). For the case with a mixing parameter of 0.5, the MRAF algorithm yields a total diffraction efficiency of 83%, and an accuracy metric below 1E-5. For comparison, the other two conditions with 8×11, and 6×6 beams diffraction patterns are also calculated (described in Supplement 1).

 

Fig. 6. Calculated accuracy metric and diffraction efficiency as a function of the mixing parameter m for the 10×10 beams diffraction pattern. The inset shows the comparison of the accuracy metric between the standard GS algorithm and the MRAF algorithm with different mixing parameter m for 10×10 beams diffraction pattern.

Download Full Size | PDF

A model in FEKO consisting of nine grating units (3×3) is developed, with 64×64 segments for each unit. As shown in Fig. 7, far-field beam pattern for a 10×10 beams phase grating at 0.85 THz is obtain. Full wave electromagnetic simulation results confirms the performance of the MRAF design. We notice that for large-pixel LO beam multiplexer applications, small angles between diffraction orders are preferred. Thus, the current definition of the NR as the diffraction orders beyond the reflection zone may be impractical. Therefore, the NR can be defined as the high diffraction order area based on the segment size, the number of pixels, and the operating wavelength.

 

Fig. 7. Simulated far-field beam pattern for a 10×10 beams phase grating at 0.85 THz designed by the MRAF algorithm, and the color of the beam pattern is on a logarithmic scale.

Download Full Size | PDF

With regard to the MRAF algorithm, similar principle has been applied to digital holography [15] and ptychography [16] to create superresolved images. Other principles has also been shown with the potential to generate multiple beam spots at terahertz frequencies, including the photonic nanojets method [17,18] and the array of solid immersion lenses method [19]. In comparison, the true phase grating method holds the advantage of achieving high diffraction efficiency by avoiding dielectric loss.

5. Conclusion

A true two-dimensional phase grating generating 2×2 diffraction beams at 0.85 THz has been realized. A 2D phase retrieval method is developed with oblique angle incident illumination. The measured total power efficiency is 81.9%, which is in good agreement with the simulation results from FEKO. It is at least 17% improvements over a pseudo-2D Fourier phase grating. This true two-dimensional phase grating not only holds higher efficiency, but in principle can generate arbitrary far-field beam patterns. Moreover, for high accuracy large-pixel terahertz multi-beam heterodyne receivers, the MRAF algorithm is analyzed. By analyzing the algorithm with different mixing parameters, the performance of the phase grating in the total diffraction efficiency and the intensity distribution accuracy is evaluated. Confirmed with rigorous full wave electromagnetic simulation results, a 100-pixel true two-dimensional phase grating is designed with an improvement in accuracy by at least four orders of magnitude compared with standard GS algorithm, and the total power efficiency is 83%, higher than the optimal value from a pseudo-2D grating. The results outline a path for the development of high accuracy large-pixel terahertz multi-beam heterodyne receivers.