1. Introduction

Two-beam interferometers are used to measure the flow speed of plasma via the Doppler-shifted wavelength of an emitted spectral line [1–4]. According to the principles of Fourier transform spectroscopy, flow of the light-emitting plasma along the interferometer line-of-sight causes a proportionate shift in the phase of the measured interference fringes. Birefringent interferometers offer an alternative to traditional double-path interferometers, with the advantage that they are less sensitive to vibration and easier to assemble [5]. In the last decade, imaging birefringent interferometers have found routine use in nuclear fusion energy research, on magnetically-confined plasma experiments such as the tokamak [4,6–11]. In this field, the technique is better known as “coherence imaging”. These instruments typically target the spectral lines in the visible wavelength range that are emitted by ions in the edge plasma region. The 2-D flow data, captured at $<1$ kHz rates, supports the study of particle and energy transport in the reactor. The instruments are compact and affordable to research institutions, and can provide better spatial coverage than is typically practical for grating spectrometers. The same technique has been demonstrated for measuring wind speeds in the upper atmosphere [12,13]. However the lower speeds ($<200$ m/s) compared to the edge of fusion experiments ($<40$ km/s) make this application more challenging. A large interferometer delay, on the order of $10^3$ or $10^4$ waves, is required to achieve sensitivity to Doppler shifts of this size.

Calibrating the zero-point of these flow measurements requires a fringe image measured at the rest-frame wavelength of the targeted emission line. For some lines it is practical to take this measurement directly using a gas-discharge lamp. For other lines, the measurement requires a tuneable laser and precision wavemeter [14–16]. If an appropriate light source is not available, then another approach is required. In the earth’s atmosphere, this can mean measuring along anti-parallel lines of sight and taking an average [3], or else by assuming that the wind speed in the vertical direction is negligible [2]. In previous work on fusion experiments, the plasma has been assumed to be stationary in the reactor’s initial breakdown phase [4] or else assumed stationary along a particular line of sight [7]. These assumptions can be difficult to verify. They can also constrain the instrument view, require periodic changes in the view to calibrate, or can require that an additional view be incorporated into the instrument design.

This paper presents a new method for absolute wavelength calibration of birefringent interferometers, useful when a tunable monochromatic light source at the appropriate wavelength is not available. The method is analogous to the standard calibration procedure for a grating spectrometer: measure spectral lines at known wavelengths and constrain a model of the instrument response across the wavelength range. A key difference with this work is the need to measure each spectral line separately, each measurement being a phase image that must be reproduced by the model. Previous work on model-fitting to birefringent interferometer fringe images has been limited to single wavelengths or to a narrow wavelength range by the inherent “$2\pi$ ambiguity” of the phase data [14,17]. We handle this ambiguity using circular statistics, allowing measurements to be made over a wide wavelength range relative to the unambiguous measurement range of the interferometer. Our work is motivated by application to the UK’s MAST Upgrade tokamak [18], where we seek to measure flows to $\pm 1$km/s accuracy. While we restrict our discussion and demonstration to interferometers producing fringes in the spatial domain (across a camera’s sensor), the principles of the calibration method apply similarly to fringes produced in the time domain.

This paper proceeds as follows. Section 2 outlines the flow measurement principle for imaging birefringent interferometers. Section 3 describes the interferometer hardware used in this work. Section 4 then describes a framework for model-fitting to wrapped phase data. Section 5 presents an experimental demonstration of the calibration method. Section 6 is a critical discussion of the method compared to the existing alternative and finally Section 7 summarises the paper.

2. Measurement principle

Figure 1(a) defines the relevant geometry of a plane-parallel uniaxial birefringent crystal plate: thickness $L$, cut angle $\gamma$ and orientation $\rho$. Figure 1(b) then defines the relevant geometry of a light ray incident on such a plate: incidence angle $\alpha$ and azimuthal angle $\beta$. A simple imaging birefringent interferometer consists of a uniaxial crystal sandwiched between two co-aligned (or crossed) polarizers [5,19] with orientations as shown in Fig. 1(c). After the first polarizer, the crystal resolves the light into orthogonally polarized ordinary (o) and extraordinary (e) rays which, experiencing refractive indices $n_{\mathrm {o}}$ and $n_{\mathrm {e}}$, propagate at different speeds and along different paths through the crystal. The two rays are resolved into a common polarization state by the second polarizer and are brought to a focus on a camera sensor by an imaging lens that is focused to infinity. When the rays interfere at the sensor, the phase delay $\phi$ between them is a function [20] of $n_{\mathrm {e}}$, $n_{\mathrm {o}}$, $\gamma$, $L$, $\alpha$, $(\beta -\rho )$ and wavelength $\lambda$.

 

Fig. 1. (a) A plane-parallel, uniaxial birefringent plate of thickness $L$. A principal section contains the crystal optic axis (COA) and the surface normal $\boldsymbol {\hat {n}}$. Orientation angle $\rho$ is between the principal section and the $x$-axis. Cut angle $\gamma$ is between the COA and the front face. (b) The plane of incidence contains the incident ray and $\boldsymbol {\hat {n}}$. Incidence angle is $\alpha$ and azimuthal angle is $\beta$. (c) The layout of a simple imaging birefringent interferometer. Orientation markings indicate polarizer transmission axes and the crystal’s principal section. (d) The corresponding interferometer delay contours across the camera sensor plane modeled for $\gamma =0$° (a waveplate). (e) For $\gamma =45$°.

Download Full Size | PDF

To produce an interferogram in the spatial domain, $\phi$ is modulated across the sensor plane. If the imaging lens has focal length $f_l$ and is considered thin, then point $\boldsymbol {p}=(x,y)$ on the sensor maps to $\alpha$ as

Here, $\mathrm {atan}2(\dots )$ is the two-argument arctangent function and the $\pi$ term appears since the lens produces an inverted image. Terms $\Delta _x$ and $\Delta _y$ account for small rotational misalignment of the crystal: $\Delta _x\equiv {}f_l\psi _y$ and $\Delta _y\equiv {}f_l\psi _x$. Crystal tilt angle about the $x$-axis is $\psi _x$, in the clockwise direction when the $x$-axis points towards the observer. Crystal tilt angle about the $y$-axis is $\psi _y$, in the anti-clockwise direction when the $y$-axis points towards the observer. For a waveplate ($\gamma =0$°) in air the expression for $\phi$ is [20]:

At normal incidence ($\alpha =0$°), this reduces to the familiar

Combining Eqs. (1), (2) and (3), we can model $\phi (\boldsymbol {p},\lambda )$ across the sensor plane. Figure 1(d) shows the hyperbolic $\phi$ contours characteristic of a waveplate. Figure 1(e) then shows the roughly linear $\phi$ contours of a crystal with $\gamma =45$°. Linear phase ramps like this are commonly used as a spatial carrier wave in fusion applications [7–9].

For the setup shown in Fig. 1(c) observing an unpolarized source, it can be shown using Mueller calculus (or Jones calculus) that the spectral flux reaching a pixel centered at point $\boldsymbol {p}$ on the sensor plane is

Here, $I_\lambda (\boldsymbol {p},\lambda )$ is the spectral flux reaching the first polarizer, relevant to the pixel in question. Delay $\phi$ is assumed constant across each pixel’s area. If the observed source has velocity $\boldsymbol {v}$ and emits a spectral line with rest-frame wavelength $\lambda _0$, then the flux reaching the pixel $S(\boldsymbol {p})=\int _0^\infty {}I'_\lambda (\boldsymbol {p},\lambda )d\lambda {}$ can be written:

Here, we have defined the shorthand $\phi _0(\boldsymbol {p})\equiv \phi (\boldsymbol {p},\lambda _0)$ and $I(\boldsymbol {p})\equiv \int _0^{\infty }I_\lambda (\boldsymbol {p},\lambda )d\lambda {}$. Unit vector $\boldsymbol {\hat {l}}(\boldsymbol {p})$ is the line-of-sight vector at the source and $c$ is the speed of light. The group delay

The phase $\Phi$ recovered from the interferogram, corresponding to the pixel centered at $\boldsymbol {p}$, is inherently wrapped into the interval $(-\pi, \pi ]$ rad. Rearranging Eq. (6), flow is recovered from $\Phi$ using:

Here, $\mathcal {W[\dots ]}$ represents the operation of wrapping within the $(-\pi, \pi ]$ rad interval. To avoid ambiguity, $\phi _0$ is chosen such that the phase shift does not exceed $\pi$ rad for the maximum expected $|\boldsymbol {v}\cdot {}\boldsymbol {\hat {l}}(\boldsymbol {p})|$. Choice of $\phi _0$ is also informed by fringe visibility $\zeta$ and by data analysis tractability for measurements of an inhomogeneous source volume [21].

From Eq. (8), the flow measurement requires two calibration images: phase delay $\phi _0(\boldsymbol {p})$ and group delay $\hat {\phi }_0(\boldsymbol {p})$. Since $\phi _0(\boldsymbol {p})$ only needs to be known modulo $2\pi$ rad, it can be measured directly with an extended monochromatic source emitting at $\lambda _0$. To measure $\hat {\phi }_0(\boldsymbol {p})$, that source can be tuned about $\lambda _0$. These procedures can be realized with a tunable laser and precision wavemeter [7,14–16]. Alternatively, both parameters can, in principle, be measured via the interferometer’s white light interference pattern using a high-resolution grating spectrometer [21]. However, this is a challenging measurement to make for an imaging system. Since $|\boldsymbol {v}|\ll {}c$, the flow measurement is far more sensitive to absolute error in $\phi _0(\boldsymbol {p})$ than in $\hat {\phi }_0(\boldsymbol {p})$. It follows that regular calibration of $\phi _0(\boldsymbol {p})$ is required to track changes in ambient temperature and alignment, while a single measurement of $\hat {\phi }_0(\boldsymbol {p})$ for each assembled interferometer is typically sufficient [7]. Our calibration method proposes that both parameters be simulated using an optimized instrument model rather than measured directly, removing the need for a tunable coherent source.

The next section describes the experimental setup and instrument model used in this work.

3. Instrument model

The interferometer used to demonstrate our calibration method has a pixelated phase mask (PPM) design [22–24]. This design generates a pixelated interferogram that maximizes the spatial resolution of the recovered phase image [23]. Figure 2(a) is a labelled photograph of the instrument, with the schematic in Fig. 2(b) showing the interferometer components and their orientations. The polarized camera sensor used here (Sony IMX250MZR CMOS sensor [25]) has an array of pixel-sized wire-grid polarizers bonded directly to the chip, oriented in the repeating $2\times {}2$ pixel pattern shown in Fig. 2(c). We use $m$ to index the camera pixels according to their polarizer orientation and refer to each $2\times {}2$ pixel block as a superpixel. If Eq. (6) is re-derived for this configuration, interferogram phase now includes an additional pixel-dependent term: $m\pi /2$ rad. Figure 2(d) shows an $8\times {}8$ pixel crop of the resulting interferogram image measured using an extended, monochromatic source. Since the source is uniform in its brightness and wavelength, the observed signal variation is due to delay modulation only. Interferogram phase $\Phi$ for a superpixel can be recovered, wrapped within $(-\pi, \pi ]$ rad, using the ‘$4$-bucket algorithm’ [22]:

Here, $S_m$ is the signal of the $m$-th pixel within the given superpixel. Phase images are recovered at superpixel resolution using this algorithm (for a detector-limited system) and example images will be shown in Section 5. It should be noted that our calibration method applies similarly to interferometer designs that instead modulate delay using arbitrarily-cut crystal plates [8,9,14,26], Savart plates [6–8,11] or Wollaston prisms [5,13,27].

 

Fig. 2. (a) The imaging birefringent interferometer used to demonstrate the wavelength calibration technique in this work. (b) Schematic showing the transmission axis orientation of the polarizer and the fast axis orientations of the waveplates. (c) The repeating $2\times {}2$ pattern of pixelated polarizers on the camera sensor, with orientation indicated by index $m$. (d) An $8\times {}8$ pixel crop from a fringe image corresponding to a uniform, monochromatic scene created using an integrating sphere and a Cd spectral lamp.

Download Full Size | PDF

The instrument shown in Fig. 2(a) is a modified version of the coherence imaging instrument deployed on the MAST and MAST Upgrade tokamaks [7,24]. Lens 1 is the objective, forming an intermediate image at its back focal plane. Lens 2 is a relay, sending collimated light through the interferometer, and Lens 3, focused at infinity, forms an image on the sensor. The lenses have focal lengths $f_1=40$ mm, $f_2=105$ mm and $f_3=150$ mm. A motorized, programmable filter wheel houses the bandpass interference filters required to isolate single spectral emission lines for measurement. All components are mounted within $50$ mm cage mounts, except the camera which is mounted directly to Lens 3. The camera sensor has format $2448\times {}2048$ and a pixel dimension of $3.45$ $\mu$m [25]. The birefringent material used is barium borate BaB$_2$O$_4$ in its alpha crystalline phase ($\alpha$-BBO). Two waveplates are used in this work, with measured thicknesses of $4.48\pm 0.05$ mm and $9.79\pm 0.05$ mm.

Dispersion in birefringent materials is appreciable — $\alpha$-BBO birefringence $n_{\mathrm {e}}-n_{\mathrm {o}}\approx {}-0.125$ varies almost $10\%$ over the visible wavelength range. This variation is modeled using the Sellmeier equation

Combining Eqs. (1), (2), (3) and (10) with the positions of the camera superpixel centers allows phase images captured by the PPM interferometer to be modeled for any wavelength. However, nominal values for the model parameters are not accurate enough to generate usable calibration images and so must be optimized using measurements. The choice of which parameters to optimize in the fit is determined by the specifics of the instrument setup and by the accuracy required. For example, for an ideal setup the crystal tilt angle model values would not need to be optimized. For the Sellmeier coefficients, multiple sets are available for $\alpha$-BBO and $\beta$-BBO, from the academic literature [29–33] and from manufacturers [34,35]. However, no information is available as to the accuracy of these coefficients. If the spread in available coefficient values is taken as being representative of their uncertainty, then modeling suggests that the coefficients are not sufficiently accurate to model phase shifts over a $>5$ nm wavelength range to the $\pm 1$ km/s accuracy target. This means that one or more of the coefficient values must be optimized to match data measured over a wider wavelength range. Correlations between model parameters also informs which parameters should be optimized. For example, crystal thickness $L$ and Sellmeier coefficient values inferred from a fit to phase will be highly correlated with one another. In other words, the data gives relatively little information about their true values. This shouldn’t worry us though as our concern here is the predictive power of the model, not the the universal applicability of our parameter values. For the instrument and targeted accuracy of this work, a fit to 8 parameters was sufficient to reproduce the measured data: $\rho$, $\psi _x$, $\psi _y$, $f_3$, $A_\mathrm {e}$, $B_\mathrm {e}$, $C_\mathrm {e}$ and $D_\mathrm {e}$. Initial values and assumptions made about these parameters will be given in Section 5.

Next, we describe the method used to fit this instrument model to measured phase data.

4. Fitting to wrapped phase data

Measurements of interferogram phase are wrapped into the interval $(-\pi, \pi ]$ rad, introducing ambiguity and discontinuities. Unwrapping to recover an unambiguous, continuous signal requires a sufficiently low noise level and high sampling rate [19,36]. For an optimized interferometer, the unambiguous measurement range for wavelength shifts will approximate the maximum anticipated Doppler shift. This is $\sim {}0.1$ nm for the visible light emitted at the edge plasma of fusion experiments—smaller than the typical wavelength separation of the bright spectral lines used in spectrometer calibration. We therefore need to handle this ambiguity carefully to interpret measurements made over a wide ($\gg 0.1$ nm) wavelength range. A Bayesian framework is well-suited for this, since prior knowledge of the model parameters can be used to remove the ambiguity in a self-consistent way.

For a set of data points $\boldsymbol {\mathcal {D}}=\{\mathcal {D}_1,\mathcal {D}_2\dots {}\}$ that is well described by a model $f(\boldsymbol {\theta })$ with parameters $\boldsymbol {\theta }=\{\theta _1,\theta _2\dots {}\}$, the likelihood function is written $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$. In the absence of prior knowledge, the most plausible parameter values are those that maximize $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$. Prior knowledge is accounted for using a prior probability density function (PDF) $P(\boldsymbol {\theta }){}$. The maximum in the posterior PDF $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}})\propto \mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}P(\boldsymbol {\theta }){}$ then gives the most plausible parameter values, given all of the available information [37]:

This is the maximum a posteriori (MAP) estimator.

Consider a single, wrapped phase data point $\mathcal {D}_i$ with noise variance $\sigma _i^2$. The ambiguity can be accounted for in the likelihood function for $\mathcal {D}_i$ using the von Mises distribution [38]:

Here, $\kappa _i \approx 1 / \sigma _i^2$ is the concentration parameter, $I_0(\dots )$ is the modified Bessel function of order 0 and model $f_i(\boldsymbol {\theta })$ corresponds to the $i$th data point. The von Mises distribution, plotted in Fig. 3, is commonly used in circular statistics as a close approximation to the wrapped normal distribution, which is less numerically tractable. We assume that the data points are statistically independent, so the total likelihood is $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}=\prod _{i}\mathcal {L}(\boldsymbol {\theta }\mid {}\mathcal {D}_i)$. By combining $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ with a suitable prior and the model described in Sections 2 and 3, we can calculate the relative probability of any possible $\boldsymbol {\theta }$ values, given a set of measurements. For wrapped data, $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ is infinitely multi-modal in model parameter space. The distance between neighboring (global) modes along a particular dimension of parameter space can be considered the unambiguous measurement range (UMR) of that parameter. Increasing the UMR reduces the amount of prior information needed to leave a single, pronounced global mode in $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ and therefore a single $\hat {\boldsymbol {\theta }}_{\mathrm {MAP}}{}$ estimate.

 

Fig. 3. The von Mises distribution is a close approximation to the wrapped normal distribution. We use it as the likelihood function for a single wrapped data point. Here, it is plotted as a function of the residual between the data point and the model for two values of concentration parameter $\kappa _i$.

Download Full Size | PDF

For models with more than a few parameters, it is not practical to evaluate $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ on a regular grid to find $\hat {\boldsymbol {\theta }}_{\mathrm {MAP}}{}$. In this work we use Gibbs sampling instead—a Markov chain Monte-Carlo (MCMC) technique where each step through parameter space is either accepted or rejected using a criterion that favors travel towards regions of high probability [37]. The sample density then converges to $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$. To perform the sampling we use the Python package “inference-tools” [39]. For wrapped data, $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}})$ possesses many local maxima. To avoid mistaking a local maximum for the global maximum, we use parallel tempering: running multiple MCMC “chains” with different settings in parallel, and occasionally swapping their positions to encourage transitions between local and global maxima [37]. In our problem, a useful way to reduce the density of local maxima in parameter space is to model the phase difference between wavelengths rather than the absolute phase. This can be seen by considering the simplified case of normal ray incidence given in Eq. (4). Consider the problem of determining optical path difference (OPD) $L(n_{\mathrm {e}}{}-n_{\mathrm {o}}{})$ from phase measurements made at two wavelengths $\lambda _1$ and $\lambda _2$. Substituting Eq. (4) into Eq. (12) and evaluating $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}})$, it can be shown that the distance between neighboring local maxima in OPD parameter space is $\sim \lambda _1$ ($\sim \lambda _2$) for a fit to absolute phase. If, instead, the model used is the phase difference

An important question for our calibration method is this: at how many wavelengths does the phase need to be measured to unambiguously constrain the model parameter values? Firstly, for the problem to be well-posed, the number of wavelengths must match or exceed the number of dispersion model parameters (e.g., Sellmeier coefficients) to be constrained. Following previous work in multi-wavelength interferometry [40,41], an estimate for the UMR of the OPD is given by the largest (finite) beat wavelength (“synthetic wavelength”) across the set of wavelengths measured:

This can be seen by deriving an expression for $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}})$ using the model in Eq. (4). Let $\delta$ be a width characteristic of the prior uncertainty in the OPD. If $\delta \ll {}u$ then the fitting problem is likely well-posed. For the set of wavelengths measured in Section 5, $u=10^{-3}$ m and, for the $\alpha$-BBO waveplates used in this work, $\delta <10^{-4}$ m. While this is a useful heuristic for gauging ambiguity, a numerical simulation of $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ as a function of OPD (using the nominal $\boldsymbol {\theta }$ values to generate a synthetic data set) is more conclusive. Clearly then, the wavelengths chosen must contain a sufficiently close pair to remove the ambiguity, while also spanning the wavelength range of interest for the calibration. The choice of wavelengths is also determined by the availability of bright spectral sources in the range of interest.

We now turn to the full problem of fitting to phase images. The spatial variation in the images generally contains important information about the model parameters that would be lost if an image measured at one wavelength were to be subtracted from all of the the other images prior to fitting. This information can be retained by instead subtracting a scalar phase value corresponding to normal ray incidence at a single wavelength. The model used for fitting can then be written

Here, $j$ indexes superpixel center position $\boldsymbol {p}_j$, $k$ indexes wavelength $\lambda _k$, $\boldsymbol {p}_{\alpha =0\text {\textdegree {}}}$ is the position on the sensor plane corresponding to normal incidence through the waveplate and $\lambda _\mathrm {ref}$ is the chosen reference wavelength. The dependence of delay $\phi$ due to a waveplate (Eq. (3)) on $\boldsymbol {\theta }$ has here been made explicit. The corresponding data used for fitting can be written

Here, $\Phi _{jk}$ is the phase measured at the $j$th superpixel and $k$th wavelength. $\Phi _\text {ref}$ is an estimate of phase at normal incidence ($\alpha =0$°) extracted from the phase image measured at $\lambda _\mathrm {ref}$ by averaging over a $10\times 10$ superpixel window centred approximately at $\boldsymbol {p}_{\alpha =0\text {\textdegree {}}}$. Substituting Eqs. (15) and (16) into Eq. (12) gives us the point likelihood $\mathcal {L}(\boldsymbol {\theta }\mid {}\mathcal {D}_{jk})$ and total likelihood $\mathcal {L}(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ is then calculated as a product over indices $j$ and $k$.

With the fitting methodology described, we next present an experimental demonstration.

5. Experimental demonstration

We will test the wavelength calibration method using spectral lines emitted by standard cadmium and zinc hollow cathode gas-discharge lamps. Figure 4(a) plots the relative intensities of eight lines emitted by atomic Cd and Zn in the $465$–$515$ nm wavelength range, as measured using a Czerny–Turner spectrometer. The six lines targeted for measurement with the interferometer have their precise wavelengths shown. For our purposes, these lines may be treated as thin singlet lines with known wavelengths [42], though uncertainty in the wavelength could be factored into the analysis. To provide an extended light source, the lamps illuminated an integrating sphere with a diameter of $6$ inches, into whose entrance port the interferometer viewed. Four interference filters with appropriate bandpass transmission profiles were chosen, as plotted in Fig. 4(a). Nominal out-of-band transmission for these filters is $10^{-4}$. To isolate a line for measurement, each filter must sufficiently suppress the neighboring lines such that their influence on the measured phase is negligible. Although filters $1$ and $4$ were not specifically procured for this work, they do meet the criterion for isolation. The $2$-inch diameter filters were mounted in the motorized filter wheel shown in Fig. 2(a). An interference fringe image was measured at each of the six wavelengths in sequence for both the $4.48$ mm and $9.79$ mm $\alpha$-BBO waveplates.

 

Fig. 4. (a) Measured transmission profiles of the bandpass filters and measured intensities (normalized) of the spectral lines used in the wavelength calibration. The precise wavelength values (in air) shown are from NIST [42] (b) and (c) plot the phase data and fit as a function of wavelength, for the centre point of each measured image, for two $\alpha$-BBO waveplates.

Download Full Size | PDF

Ambient temperature drift causes a proportionate fringe phase drift via the waveplate’s thermal expansion and the thermo-optic effect. Without thermal stabilization, a short measurement window was important to minimize systematic error. Camera exposures and filter selection were automated, allowing data collection within a $3$ minute window. Automated monitoring of the phase at a single wavelength in the $20$ minutes before and after the acquisition of the calibration data suggests a linear drift of $\approx 0.01$ rad during the acquisition. This corresponds to a temperature change of $\sim 0.02$°C [30]. A first order correction to the calibration data was applied by interpolating the monitoring data in time.

The six raw phase images for the $4.48$ mm waveplate are plotted in the left column of Fig. 5. The standard deviation of the noise on a single phase data point was $\sigma \approx 0.01$ rad, estimated using the spread of values in the flat central region of the phase images. The phase at the centre of the images is plotted for a reduced wavelength range for both waveplates in Fig. 4(b) and 4(c).

 

Fig. 5. Left column: measured fringe phase images at six wavelengths. Black contours show the data and white contours show the fit. Right column: the corresponding residuals between data and fit, with histograms on the colorbars. Data corresponds to the $4.48$ mm $\alpha$-BBO waveplate.

Download Full Size | PDF

To test the predictive power of the optimized model, data from 5 wavelengths was used to constrain the fit while the data measured at $472.2$ nm was withheld. Table 1 lists the starting values and prior PDF widths used for each parameter. The starting values used for the extraordinary $\alpha$-BBO Sellmeier coefficients, and the fixed values for the ordinary coefficients, come from [34]. The prior PDF width for the extraordinary coefficients was set at $40\%$ of the starting value, which comfortably reproduces the spread in available coefficient values. The phase images were downsampled in each dimension by a factor $50$ prior to fitting to improve performance. Each fit was run for 200,000 steps on a desktop PC and typically converged after $\sim$50,000 steps, taking around $30$ hours total on a standard desktop PC. Each fit was repeated ten times and found to be reproducible for a range of starting parameters selected from the range of available Sellmeier coefficient sets.

Table 1. Model parameter values before and after fitting to the measured phase images for both waveplates tested.^a

View Table

The $\hat {\boldsymbol {\theta }}_{\mathrm {MAP}}{}$ parameter values produced by the fits are listed in Table 1. Figure 5 plots the contours of the optimized model atop the measured phase in the left column and plots the residuals between measured and modeled phase in the right column. Figure 4 plots the modeled phase at the image centre as a function of wavelength, clearly showing the large number of phase wraps between some of the wavelengths measured. The residuals have standard deviation $\sigma \approx 0.02$ rad corresponding to a flow speed of $\approx {}0.7$ km/s for the wavelengths considered. The residual mean is smaller than this for all wavelengths except $467.8$ nm, where it is comparable. Importantly, the residual is consistent for the $472.2$ nm data that was not used to constrain the fit, indicating the model’s predictive accuracy. Results for the $9.79$ mm waveplate have comparable flow-equivalent residuals to the $4.48$ mm waveplate. When the process was repeated with different wavelengths withheld from during the fit, the results were comparable. Table 1 lists the $\hat {\boldsymbol {\theta }}_{\mathrm {MAP}}$ parameter values for fits to both the waveplate thicknesses tested. Figure 6 then plots the kernel density estimate of the posterior PDF $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ for the $4.48$ mm waveplate, from which the parameter (co)variance values can be inferred. The expected correlation between Sellmeier coefficient model parameter is observed.

 

Fig. 6. A kernel density estimate of the posterior PDF $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ for the fit to the $4.48$ mm $\alpha$-BBO waveplate data shown in Figs. 4 and 5, generated using the Markov chain Monte-Carlo samples from the distribution. The $8$-dimensional $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ is shown as $28$ marginalized joint distributions, one across each parameter pair. Parameter symbols and units are defined in Table 1. The red dot indicates the maximum a posteriori estimate $\hat {\boldsymbol {\theta }}_{\mathrm {MAP}}{}$.

Download Full Size | PDF

The optimized model parameters also allow us to calculate group delay $\hat {\phi }_0$ for the instrument. Measurements of $\hat {\phi }_0$ for the two waveplates tested have been published in previous work [7], where they were measured using a tuneable laser at $460.9$ nm. Values for $\hat {\phi }_0$ at this wavelength generated using the posterior PDF $P(\boldsymbol {\theta }\mid {}\boldsymbol {\mathcal {D}}){}$ are $1412\pm 1$ waves and $3091\pm 2$ waves for the $4.48$ mm and $9.79$ mm waveplates respectively. The corresponding values from the previous work are $1401\pm 10$ waves and $3142\pm 16$ waves respectively. Roughly $10$ waves of discrepancy can be expected due to instrument misalignment evident in the data from the previous work. However, this cannot fully explain the discrepancy in the two results for the $9.79$ mm waveplate. The uncertainty values in these $\hat {\phi }_0$ estimates were found by modeling $\hat {\phi }_0$ for $1000$ samples of $\boldsymbol {\theta }$ from the fits.

6. Discussion

The method described here must be compared to the standard for absolute wavelength calibration of Doppler flow measurements: using a tuneable laser to directly illuminate the instrument at the appropriate wavelength. The systematics in this work limited the calibration accuracy to $\pm 0.7$ km/s flow equivalent. This is comparable to the accuracy achieved by laser-calibrated interferometers used in fusion experiments [10,16]. The advantage of using lamps (and the motivator for this work) is the substantially lower hardware cost. As well as being a cheaper source than a laser, lamps do not require precision wavelength measurement. That said, this method is less direct and has limitations that will now be discussed.

Accuracy is limited by the agreement between modeled and measured phase. An etalon effect in the filter is the likely cause of the periodic residual pattern clearly visible at $508.6$ nm and somewhat visible at $467.8$ nm and $468.0$ nm (Fig. 5). Other contributors across all wavelengths could include crystal inhomogeneity, non-flat surfaces, and image distortion. An investigation into non-ideal crystal behaviour could feasibly be carried out and incorporated into the model. Alternatively (or additionally), higher tolerances for homogeneity and flatness could be specified during procurement. Similarly, distortion could be incorporated using a standard distortion model, or by replacing the thin-lens approximation with a ray-traced model. If the systematics could be confidently considered to be unchanging wavelength, as was observed between $472.2$ nm and $481.1$ nm in Fig. 5, then they could be subtracted from from the modeled phase to improve accuracy.

It was found in this work that the available BBO Sellmeier coefficients could not accurately reproduce the measured phase shift data over the wavelength range tested. As mentioned in Section 3, this is to be expected based on the spread in values across the available coefficient sets. For other birefringent materials used in interferometry (e.g., lithium niobate or yttrium orthovanadate) the available Sellmeier coefficients may be sufficiently accurate.

Although the MCMC sampling required to perform the fit took $\sim 10$ hours, this is for a single desktop PC running non-optimized Python code. It is also likely that the optimization algorithm used here —the parallel-tempered Gibbs sampling discussed in Section 4— is sub-optimal at sampling a complex PDF with many local maxima like this one. Figure 6 appears to show some evidence of local maximum sampling for the $D_\mathrm {e}$ parameter. There is therefore large scope for improving the computational and algorithmic efficiency of the method, which could reduce analysis time significantly.

While the accuracy demonstrated here meets our $\pm 1$km/s target, it does not approach the accuracy required for atmospheric windspeed measurements.

7. Summary

We have presented a new method for absolute wavelength calibration of imaging birefringent interferometers for the spectroscopic measurement of ionised gas flows. These instruments are currently deployed on nuclear fusion experiments around the world and are being investigated for application to measurements of the upper atmosphere. Instead of measuring the “zero-flow” reference image directly with a source at the rest-frame wavelength of the targeted spectral line, we propose using an instrument model to simulate the reference image. Due to uncertainties in the instrument parameters, this requires fitting the model to measurements made at a small number of nearby wavelengths. The key development here is a framework for model-fitting to fringe phase measurements made over a wavelength range much wider than the unambiguous measurement range of the instrument. This allows standard lamp sources to be used for the calibration. The fitting framework is based on established statistical techniques and requires making modest assumptions about about the instrument model parameters. The method was validated by showing that measurements of 5 lines across a $40$ nm wavelength range were enough to accurately predict the phase at a wavelength within that range. The accuracy achieved corresponded to $\pm 0.7$ km/s flow equivalent. This is comparable to the accuracy reported in fusion experiments that use a tuneable laser source to calibrate the measurement directly. The method represents a significant reduction in calibration hardware costs, but requires a sophisticated data analysis routine. Difficulties are foreseen reaching the higher accuracy required to measure the slower flows of the upper atmosphere.