#### Abstract

A simulation method for acquiring spectrometer’s Spectral Response Function (SRF) based on Huygens Point Spread Function (PSF) is suggested. Taking into account the effects of optical aberrations and diffraction, the method can obtain the fine SRF curve and corresponding spectral bandwidth at any nominal wavelength as early as in the design phase. A prism monochromator is proposed for illustrating the simulation procedure. For comparison, a geometrical ray-tracing method is also provided, with bandwidth deviations varying from 5% at 250 nm to 25% at 2400 nm. Further comparison with reported experiments shows that the areas of the SRF profiles agree to about 1%. However, the weak scattered background light on the level of 10^{−4} to 10^{−5} observed by experiment could not be covered by this simulation. This simulation method is a useful tool for forecasting the performance of an underdesigned spectrometer.

#### 1. Introduction

The spectral calibration [1, 2] of a spectrometer mainly consists of the accurate measurement of the Spectral Response Function [3, 4] (SRF, also known as Instrument Function [5]), from which the peak response wavelength and the spectral bandwidth of each detection channel can also be derived. Traditionally, experiments are usually performed to conduct the spectral calibration. Generally, a tunable monochromatic source such as a tunable laser [6] or a monochromator with higher resolution [7] is used to stimulate the spectrometer from wavelength to wavelength. For monochromators, the inverse process is also carried out, in the way of rotating the dispersion elements to scan a wavelength-fixed source. In addition, white-light interferometry [8], observing the spectral lines from a tungsten/discharge lamp [9, 10], or comparing reflectance of two different doped diffuse panels [11, 12] is also used by (imaging) spectrometers. The tunable wavelength scanning method theoretically can get all the SRF curves of each detection channel, but the operation is always so time consuming that only some typical positions or wavelengths can be measured in practice. Spectral line lamp and doped diffuse panel are only used for calibrating the peak response wavelength and are more effective for linear dispersion instruments. The white-light interferometry technology is not yet fully mature. Besides, experimental methods cannot be implemented until the instrument is assembled and adjusted. A simulation of SRF based on Line Spread Function (LSF) is reported by Mouroulis et al. [13]. The contour of the SRF is evaluated but it does not establish the relationship between the simulated SRF and the wavelength. Some parameters such as central wavelength and spectral bandwidth are still needed. Therefore, based on Huygens PSF, which is more accurate than Fast Fourier Transform (FFT) PSF, a modified simulation method for SRF and spectral bandwidth is suggested. Even in the optical design phase, the simulation procedure can also be carried out immediately at need to forecast the designer’s concerned specifications. This simulation method takes account of the effects of optical aberrations and diffraction, which are the two most important disturbances to SRF. The alignment and fabrication errors can be added as well for a postsimulation. In this paper, a spectrograph using array detectors is considered as a monochromator with multiple consecutive exit slits, and the detector characteristics are not considered here as they are relatively independent from the instrument’s intrinsic dispersion properties.

#### 2. Basic Principles in Spectral Response Function Simulation

The spectral integrated energy received by a spectrometer’s exit slit can be regarded as the slit’s response to the superposition of all the shifted monochromatic images of entrance slit [3, 13]. Considered in the wavelength domain, the received energy can be transformed to spectral density distribution versus wavelength, which is, namely, the SRF function of the exit slit. The SRF in fact is a spectral filer function of the spectrometer. Slit function can be simply assumed to be a rectangular one; image of entrance slit is the superposition integral of the entrance slit with the optical system Point Spread Function (PSF); if the optical system satisfies spatially invariant condition, the superposition integral can be replaced by convolution [14]. SRF of traditional spectrometer only focuses on the dispersion direction (assumed to be tangential direction), so the slit image will be obtained by using tangential Line Spread Function (). Typically, with the entrance slit having the same width as the exit slit, say, , the SRF of linear space invariant optical system can be expressed as where denotes cross-correlation and denotes convolution; is the tangential coordinate on image plane, which can be mapped onto wavelength domain coordinate by chief ray iterative tracings. Cross-correlation has an equivalent result on condition that the images of wavelengths near central are of the same profiles except being shifted. If is symmetric, the cross-correlation can be replaced by convolution. If is narrow compared with the slit width or ideally close to a Dirac delta function, then the SRF approximates to a triangle function with a bandwidth (FWHM, Full Width at Half Maximum) equivalent to in image plane coordinate. Otherwise, if the size is comparable with the exit slit, then it approaches a Gaussian function with a bandwidth broader than . If SRF can be approximated by a triangle or Gaussian function, all the spectral response characteristics can be identified by peak response wavelength and FWHM accurately.

The Huygens wavelets method and FFT are the common method to compute the diffraction PSF. Among them, the Huygens wavelets method imagines each point on a wavefront as a perfect point source with an amplitude and phase. Each of these point sources radiates a spherical “wavelet.” The diffraction intensity at any point on the image surface is the complex sum of all these wavelets, squared. The PSF is computed this way for every point on the image grid. So the Huygens PSF accounts for any local tilt in the image surface caused by the image surface slope, the chief ray incidence angle, or both. The optical system of a prism or grating spectrometer is not always axial symmetric as tilt surfaces often used, so using direct integration of Huygens wavelets method in place of FFT of the complex amplitude of wavefront at exit pupil has a higher accuracy in computing system PSF or . Generally, PSF is a nonelementary function with respect to either wavelength or spatial coordinates, so only the numerical expression form can be obtained through numerical analysis methods. A practical way to get the results is to extract them from optical design software (such as Zemax), where the calculation algorithm has been tested and verified in practice.

The numerical calculation procedure of SRF is based on (1). In order to establish the relationship between the contour of SRF and the corresponding wavelength, before the convolution, the coordinate system in wavelength space should be mapped onto a position coordinate system on image plane. It can be achieved by iteratively tracing chief rays on axis FOV (Field of View) at a serial of wavelengths. The wavelength corresponding to a chief ray intersected with the image plane at a certain position is referred to as the position’s* nominal wavelength*, while the position is called as the wavelength’s* nominal position*. This operation is performed only along the dispersion direction. For the wavelength , the in (1) can be calculated by integrating the discrete direction FOV Huygens . If the linear space invariant condition is satisfied, then only the zero FOV is needed. The monochromatic images of can be obtained by convoluting the with the rectangular object function. Then the calculation of SRF could be finished through a further convolution with the exit slit function. The SRF function corresponding to this exit slit could be obtained through inverse mapping the results from position coordinate system back onto the wavelength coordinate system.

The mutual mapping between wavelength and position coordinates is introduced to solve the nonlinear dispersion problem, as cross-correlation of exit slit with monochromatic images is carried out on a uniformly spaced position grid, while the SRF is a function of wavelength described in the wavelength domain. Calculation speed is the main drawback of Huygens PSF method. As the of specific wavelength changes continuously and smoothly with position coordinates on image plane, an algorithm of spline interpolating extracted of larger intervals will accelerate calculations effectively.

#### 3. Demonstration of Simulation Instance

A Féry prism monochromator [15] is proposed for illustrating the simulation algorithm. The monochromator adopts only one refractive-reflective prism with curved surfaces for light dispersion, while the wavelength scanning from 250 nm to 2500 nm is achieved via rotating the prism. The only exit slit is coplanar with the entrance slit and of the same size. Optical layout and construction parameters of the monochromator are showed in Figure 1 and Table 1, respectively. At each angle within the ±2.5° prism rotation range, the transverse aberration along dispersion direction is smaller than 8 micrometers. This optical system is used for solar spectral irradiance measurement [16].

##### 3.1. Analysis of Monochromatic Imaging Characteristics

Geometrical ray traced images of entrance slit at four different wavelengths are showed in Figure 2. Ideally, at each arbitrary wavelength, entrance slit should be imaged onto the spectrum plane without any scaling or distortion. In fact, aberrations more or less exist and will directly affect the imaging process through ray tracing. As can be seen from the simulation results, along the dispersion direction, or, namely, the slit width direction, edges of the monochromatic images are clearly sharp which indicates good wavelength separating capability, whereas the edges along the cross dispersion direction are of a certain degree of diffusion, which however will not directly affect the wavelength resolution property. It also can be seen that the width decreases slightly as wavelength increases, and the image appears to distort in some sort at the same time.

**(a)**

**(b)**

**(c)**

**(d)**

The result of coordinate mapping from 294 nm to 306 nm is shown in Figure 3. It is achieved by iteratively tracing chief rays on axis FOV from wavelengths of 294 nm to 306 nm. Figure 4 shows four monochromatic image intensity distributions along the dispersion direction. They are obtained through the convolution in (1). As is computed based on Huygens wavelet integration which has already considered both the optical path difference and wave character of light, thus the aberration and diffraction effects are included in the image acquisition. The image profile is near rectangular at short wavelength and then reshapes to trapezoidal and at last to triangular gradually as wavelength increases. Also, it should be noted that, in the wings of the longer wavelength image, there are two more and more obvious bulges which will not be possible to appear by geometrical ray-tracing method. The bulges are mainly attributed to the diffraction effect. The relationship between the image profile’s FWHM width and wavelength is presented in Figure 5. Clearly the widening trend with wavelength is distinct from the slightly decreasing condition of geometrically ray traced spots as showed in Figure 2. Diffraction effect is believed to directly contribute to the widening condition and be evaluated quantitatively here. These monochromatic images can be directly recorded by array detectors or films when the separate photosensitive element is infinitely smaller than the image size. This is easily inferred by setting the exit slit width to Dirac delta function according to (1).

**(a)**

**(b)**

**(c)**

**(d)**

##### 3.2. SRF Curves at Nominal Peak Wavelengths

According to (1), the SRF is considered as a cross-correlation of the exit slit function and the image function (the convolution result of terms in the square brackets). Figure 6 shows two SRF curves at nominal peak wavelengths of 299.0 nm and 1912.2 nm, respectively. Obviously the SRF profile of 299.0 nm is close to a triangular shape, though the peak and bottom corners are smoothed to some degree. The maximum response is about 0.95 other than the ray traced result of almost 1, which means about 5% of input energy at this peak wavelength is not caught by the exit slit but redistributed to the neighborhood area. An inverse case is when the vicinal wavelengths will also contribute extra light to this slit. Further data processing is conducted as follows: several data points near the nominal peak are parabolic fitted to get the finer resolving peak wavelength and the corresponding maximum response; then a superposition integral of the SRF over the horizontal wavelength coordinate is performed, while the wavelength coordinate may be not linear as it is nonlinearly inverse mapped from the computing grid; finally, the effective FWHM is obtained by dividing the throughput of the SRF by the maximum response. The parabolic fitting model is set as follows:where is the fitted peak wavelength and and are another two fitting parameters. The term “effective FWHM” is used to distinct it from the one that is originally defined as the wavelength interval between the geometrically half maximums of the SRF profile. The difference between the two FWHMs may be slight at most time, but when in spectroradiometry applications, the effective FWHM acts more in deriving the exact spectral radiance intensity.

**(a)**

**(b)**

Through the method described above, the fitted peak wavelength corresponding to nominal 299.0 nm SRF is 299.0 nm and the effective FWHM is 1.4 nm; the deviation of the fitted peak wavelength from the nominal 299.0 nm is less than 0.1 nm. The RMS error of 1.5 × 10^{−4} in the parabolic fitting is preferred. The fitted curves are plotted in solid lines. Then the similar fitting process is applied to the longer nominal wavelength of 1911.4 nm with a RMS fitting error lower than 10^{−6}. The fitted peak wavelength and effective FWHM are 1912.2 nm and 37.5 nm, respectively. The remarkable 0.8 nm deviation of fitted peak wavelength from the nominal 1911.4 nm is believed to be a result of the nonrotationally symmetric optical structure and the exit lights not exactly perpendicularly passing through the exit slit. Since the whole SRF profile at this wavelength is so deeply smoothed, a Gaussian model can be tried to fit the main part of the SRF. In the actual operation, the middle 5/9 of the entire 401 data points is selected. The fitted peak wavelength is 1912.2 nm which is the same as the parabolic fitting result; the RMS fitting error is 0.4%. The two same fitted peak wavelengths are just coincidence but the Gaussian RMS error is so larger compared with parabolic fitting. This is the best fitting result when using Gaussian model all over the wavelength range. The geometrical FWHM value is 35.0 nm. However, the triangular and Gaussian shapes have both emerged in the simulation which is consistent with the former analysis of (1). Since the SRF always has a rounded peak and the effective FWHM usually acts more in practice, the method of parabolic fitting and throughput-equivalence is more general in describing SRF than simple triangular or Gaussian model fitting.

##### 3.3. Spectral Resolution over the Full Response Wavelength Range

The exit slit of the monochromator locates at 35 mm above the entrance slit, and the wavelength scanning of exit slit is achieved through prism rotation. A further simulation of the spectral resolution or, namely, the effective FWHMs over the full wavelength range is carried out based on the Huygens SRF simulation results. Additionally, a geometrical ray-tracing method is also conducted for a comparison. Both of the two resolution results versus rotation angle/nominal wavelength are plotted in Figure 7. At a given rotation angle, the geometrical method traces 401 monochromatic images at a series of wavelengths corresponding to the position grid in Huygens method; for each wavelength 5000 rays are firstly traced with FOVs and pupil positions randomly generated; then, the rays trapped by the exit slit are counted. The distribution of the ray numbers with the wavelength is the SRF at this rotation angle. Optical aberrations are inherently included in the ray-tracing process but diffraction limit is ignored. The relative resolution deviation defined as the ratio of the absolute deviation to the average of the two is also plotted in Figure 7. The relative deviation keeps increasing with wavelength. At 250 nm, the relative deviation is 5%; then it reaches up to 10% at 600 nm and 25% at 2400 nm. The enlarging deviation trend indicates the enhanced diffraction effect in the longer wavelength region, as well as the importance and necessity of diffraction analysis especially for infrared optical system.

##### 3.4. Comparisons with Reported Experiments

Ideally, a tunable monochromatic source covering all the wavelengths is used to obtain the SRF function directly. References [6, 17] have been validated by experiments that the inverse process is a preferable approximation, in which dispersion element is rotated for scan but source wavelength is fixed. This process can be considered as the slit moving relative to the image and scanning it. The response can be expressed aswhere denotes the source wavelength, *λ* is the nominal wavelength of the “shifted” slit center, denotes the slit’s relative displacement, and and , respectively, represent image function and slit function. The nominal position of is selected as the displacement’s origin, and the longer wavelength’s nominal position is assumed to move along the positive direction. When illuminated by a chromatic source, the monochromatic images near the exit slit can be approximated as the translations of the image of the slit center’s nominal wavelength. Then the image of a vicinal wavelength *λ* can be expressed asSo the general response function through direct scan method isClearly, (3) and (5) gain the same results. Then a comparison simulation is conducted according to the real working process rather than the approximate expressions in (3) and (5), and the results are shown in Figure 8. At wavelengths of 628.5 nm and 1912.2 nm, the response residuals of different scan methods are between ±0.01, and the area differences are lower than 0.6% so also the effective SRF deviations. The residual level is coincident with reported result in [17].

**(a)**

**(b)**

A new design was created by Monte Carlo analysis when tolerancing the original design. This new design has included randomly generated alignment and fabrication errors which have resulted in the RMS spot radiuses being doubled. Figure 9 shows four Huygens SRF curves by direct scan with logarithmic scale for the -axis; the ray traced ones are also plotted for comparison. The area differences of the Huygens PSF and the ray trace modeled profiles are between 0.9% and 1.2%. When response is below 10^{−2}, two shoulders grow in the Huygens SRF profile. Vertical grid lines are added as mirror-symmetric in the vicinity of each curve’s maximum. So it is easy to see the asymmetry in the shoulders of each Huygens SRF profile. But below 10^{−4} or 10^{−5}, the Huygens SRFs quickly drop to zero without any background level. This case disagrees with experiments in [17]. Comparisons with stray light-free measurements by Kristensson et al. [18] show that the stray light also plays an important role when the response is below 10^{−2} or even 10^{−5}. However, the simulation based on Huygens SRF fails to include the stray light. Between 10^{−2} and 10^{−4}, the Huygens SRF characteristics agree well with that experiment, as the diffraction effect is dominant in this response interval.

#### 4. Conclusions

The SRF and spectral resolution simulation method based on Huygens PSF is proposed to quantificationally characterize the spectroscopic properties of a spectrometer. A wide-spectral prism monochromator is taken to illustrate the simulation procedure. The results are compared with that from geometrical ray-tracing method, and the resolution deviations from 5% to 25% indicate the necessity and importance of diffraction analysis. Finally, comparison with reported experiment shows that the simulation method agrees well with that except in the far wing regions where scattered background light on the level of 10^{−4} or 10^{−5} is dominant. The simulation algorithm is a useful complement to current commercial optical design software in forecasting the performance of underdesigned spectrometer.

#### Competing Interests

The authors declare that there are no competing interests related to this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant no. 41474161.