Bayesian framework for THz-TDS plasma diagnostics

Terahertz time-domain spectroscopy (THz-TDS) is an optical diagnostic used to noninvasively measure plasma electron density and collision frequency. Conventional methods for analyzing THz-TDS plasma diagnostic data often do not account for measurement artifacts and do not quantify parameter uncertainties. We introduce a novel Bayesian framework that overcomes these deficiencies. The framework enables computation of both the density and collision frequency, compensates for artifacts produced by refraction and delay line errors, and quantifies parameter uncertainties caused by noise and imprecise knowledge of unmeasured plasma properties. We demonstrate the framework with sample measurements of a radio frequency inductively-coupled plasma discharge. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Optical diagnostics play a crucial role in experimental plasma physics. In contrast to electrostatic probes, optical diagnostics do not perturb plasma conditions and are not susceptible to damage by plasma sputtering [1,2]. The application of optical plasma diagnostics, such as optical emission spectroscopy, laser Thomson scattering, and microwave interferometry, has led to advances in fusion energy [3,4], materials processing [5], and spacecraft propulsion [6,7].
Terahertz time-domain spectroscopy (THz-TDS) is a relatively novel optical plasma diagnostic that utilizes pulsed THz radiation to measure plasma electron density and collision frequency [8][9][10][11][12][13][14][15][16]. The picosecond-scale duration and broadband nature of the THz pulses enables THz-TDS to make time-resolved measurements across a broad range of plasma conditions [17]. Further, because THz radiation transmits through many common plasma insulators (e.g., boron nitride, alumina, and silica), THz-TDS is a candidate for probing plasma discharges that are inaccessible to other optical techniques [18].
However, despite the many advantages of (transmission-mode) THz-TDS for (gaseous) plasma diagnostics, the capability of the technique is limited by conventional data analysis methods. The governing equations are often simplified to preclude determination of the collision frequency, data are rarely corrected for measurement artifacts, and parameter uncertainties are rarely reported [8][9][10][11][12][13][14][15][16]. These simplifications are often motivated by difficulties associated with the extraction of plasma parameters from THz-TDS measurements. Time-domain data are Fourier-transformed to produce a large set of nonlinear frequency-domain equations that relate complex refractive indices to the plasma parameters. In conventional methods, the electron density and collision frequency are recovered through numerical inversion of these nonlinear equations. The analysis is complicated by the inclusion of the collision frequency, implementation of data corrections, quantification of parameter uncertainties, and appropriate incorporation of data measured across the frequency spectrum.
An analysis method that is especially well-suited for parameter estimation and uncertainty quantification (UQ) in this context is Bayesian inference [19]. We present a novel Bayesian framework that compensates for measurement artifacts and appropriately weights measured data to enable robust computation of both the electron density and collision frequency, as well as the quantification of their respective uncertainties. We begin by presenting our experimental configuration for THz-TDS measurements of a radio frequency (RF) inductively-coupled plasma (ICP) discharge. We then review the THz-TDS plasma equations and introduce a novel correction factor that offsets the effects of refraction. Next, we detail the Bayesian framework for analyzing THz-TDS plasma measurements. Finally, we apply our approach to experimental measurements of the RF ICP discharge and demonstrate the benefits relative to conventional signal analysis.

THz-TDS system
Shown schematically in Fig. 1, the THz-TDS system generates THz radiation via the optical pumping of a photoconductive antenna (PCA) [20] and detects THz radiation electro-optically by utilizing the Pockels effect in a <110> ZnTe crystal [21]. The pump and probe beams consist of ultrashort laser pulses emitted from a Ti:sapphire laser (Coherent Vitara-T HP) at a repetition rate of 100 MHz. The pulses have a measured temporal FWHM of approximately 50 fs and center wavelength of 800 nm. Group velocity dispersion (GVD) is mitigated through the use of commercially available low-GVD optics and pre-compensated for by a pulse compressor (Coherent CPC-II). The PCA is a low-temperature-grown gallium arsenide parallel-line antenna with a dipole gap of 5 µm (BATOP PCA-60-05-10-800-h). A DC bias of 10 V is applied across the dipole gap with an external power supply. The pump beam is focused onto the dipole gap with an aspheric lens and limited to an average power of 5 mW to avoid damaging the PCA. Photoelectrons produced in the PCA substrate by each pump pulse respond to the PCA voltage bias by generating a time-varying current that, in conjunction with the short recombination times of the charge carriers, results in the emission of a THz pulse.
The THz pulses emitted by the PCA are coupled to air with a hyper-hemispherical float-zone silicon lens. The THz beam is collimated, steered through the plasma discharge, and focused onto the ZnTe crystal by a pair of gold parabolic mirrors. The ZnTe crystal (Eskma Optics ZnTe-1000H) has a thickness of 1 mm and is oriented such that the pump and THz beams propagate along the [110] axis with polarizations 45 • from the [001] axis [22]. Each THz pulse induces birefringence in the ZnTe crystal, thereby causing a differential phase retardation in each probing laser pulse that is linearly proportional to the electric field strength of the THz pulse [23]. The differential phase retardation in the probe pulse is measured using a quarter-wave plate, Wollaston prism, and balanced photodetector connected to a lock-in amplifier (LIA) and computer running a custom LabVIEW Virtual Instrument. The 1 kHz reference signal for the LIA (Zurich Instruments MFLI) is provided by an optical chopper in the pump path. The femtosecond-scale probe pulses are scanned across the picosecond-scale THz pulses with a mechanical delay line in the pump path.

RF ICP discharge
The RF ICP discharge produces a stable, uniform plasma with an approximate length of 12 cm. The discharge consists of a custom quartz tube cross (50 mm outer diameter, 3 mm wall thickness) connected to steel KF vacuum flanges by quick-connect couplings. A rotary-vane mechanical pump (Pfeiffer Adixen PASCAL 2021SD) evacuates the discharge chamber to a base pressure of 1 mTorr, as measured by a convection pressure gauge (Kurt J. Lesker KJL275807LL). Ultra-high purity (99.999%) argon is fed into the discharge by a regulator and precision flow meter. Pressure values outputted by the pressure gauge are corrected for argon according to fits made from argon calibration data provided by the manufacturer.
RF power is coupled to the plasma via a three-turn hollow copper antenna wrapped around the quartz tube. The unbalanced 13.56 MHz RF signal is generated by an RF power supply (Materials Science, Inc. RF-3-XIII), tuned by an RF radio antenna tuner (Palstar HF-AUTO), and converted to a balanced signal by a custom 1000 Ω current balun connected to each end of the antenna. Power is carried between devices by 50 Ω RG400 coaxial patch cables with UHF connectors (Pasternack PE3743 series). The plasma discharge is always operated with a standing wave ratio of 1.05 or less, and forward power drift was observed to remain within ±1 W throughout the duration of all THz-TDS measurements.
An electrically grounded Faraday cage surrounding the antenna prevents stray electromagnetic radiation or electrical arcs from interfering with THz-TDS measurement equipment. The discharge viewports (Torr Scientific BKVPZ50NQZ) feature Z-cut crystalline quartz windows that, compared to conventional amorphous quartz windows, exhibit relatively low absorption in the THz regime [24]. The distance between the visible edge of the plasma and the nearest quartz window was always greater than 6 cm.

Standard governing equations
THz-TDS measures the electric fields of sample and reference THz pulses to determine the complex plasma refractive index [17]. The sample field is recorded with the plasma in the pulse path, and the reference field is recorded without the plasma in the pulse path. Assuming the plasma is bounded by vacuum, reflections at the plasma-vacuum boundary may be ignored, and pulse scattering is negligible, the sample (Ê s ) and reference (Ê r ) electric field spectra are related to the complex plasma refractive index (ñ) bŷ where ω is the angular frequency of the respective harmonic component comprising the THz pulse, c is the vacuum speed of light, L is the plasma length, and z is the position along the pulse propagation axis [25]. As noted by Jamison et al. [8], because the real component of the plasma refractive index is typically near unity, the impact of reflections at the plasma-vacuum boundaries are generally below the detection limit of conventional THz-TDS systems and can therefore be ignored. Other reflections, such as those at the viewport-vacuum interfaces, are present in both the sample and reference signals and are therefore cancelled in the transfer function [26]. If the plasma is uniform along the z-direction, Eq. (1) simplifies tô In many cases, Eq. (2) is also a valid expression for the average plasma refractive index in non-uniform plasmas. The error incurred for ignoring the distribution is generally negligible, as long as the THz pulse frequency components at which the complex refractive index is measured are not near the plasma cut-off frequency [18].
Assuming electron temperature and applied magnetic field effects are negligible, the complex refractive index is related to the angular plasma frequency (ω p ) and angular electron collision frequency (ν) byñ The plasma frequency is a function of the electron density (n e ), according to where e is the elementary charge, ε 0 is the permittivity of free space, and m e is the mass of an electron [27,28]. It is convenient to decompose the measured spectral ratio of Eq. (2) into a magnitude (A) and phase (Φ):Ê According to the governing relations presented in Eqs. (2) and (3), the model magnitude (A m ) and phase (Φ m ) are explicitly related to the plasma parameters by and It is important to emphasize that, in this work, we utilize the following convention: all frequencies presented in the equations are angular frequencies (rad/s) and all frequencies plotted in the figures are linear frequencies (Hz or THz).

Modified model equations
The model equations presented above do not account for experimental artifacts that can skew estimates of the electron density and collision frequency. One such artifact, recently investigated by Meier et al. [13], is inaccuracy in the temporal axis caused by delay line registration error. Measurements recorded in the time domain are Fourier-transformed to obtain the THz pulse electric field spectrum. As a result, delay line registration error in the time axis (δt reg ) produces a systematic error in the frequency domain according to where F is the Fourier operator, E is the time-domain electric field, andÊ is the true spectrum. The modified phase model (Φ mm ) includes this effect and is therefore given by We determine the delay line registration error using the Bayesian framework introduced in the next section. By comparison, Meier et al. [13] use a regression to eliminate the influence of delay line errors within the conventional framework for THz-TDS analysis.
Another experimental artifact is refraction of the THz pulse. Even if the mismatch between the real components of the vacuum and plasma refractive indices is small, non-perpendicular propagation of the THz beam across plasma-vacuum interfaces (caused by curvature in the plasma boundary, frequency-dependent THz beam divergence, errors in THz beam alignment, etc.) can cause the plasma to refract a non-negligible portion of the THz beam power away from the detector. This refraction artificially decreases the measured spectrum magnitude and results in the calculation of erroneous electron collision frequencies. As a result, the computed electron collision frequency increases with THz radiation frequency. This effect can be seen in the uncorrected data shown in Fig. 2. The impact of refraction has already been well-documented by researchers working with microwave interferometry, an optical plasma diagnostic that also relies on refractive index measurements to infer plasma properties [27]. However, because microwave interferometry typically uses continuous-wave radiation with a single resolvable frequency component, correcting for signal losses due to refraction often requires precise ray tracing.
THz-TDS, on the other hand, utilizes broadband THz pulses comprising a broad bandwidth of multiple resolvable frequencies. Data recorded across the spectrum bandwidth can therefore be collectively leveraged to correct for refraction effects. The real portion of the complex plasma refractive index, which governs beam steering, does not vary significantly across the measured spectrum. Refraction thus attenuates the various frequency components of the THz beam by approximately the same fraction. We therefore introduce a new, frequency-independent correction factor (r) equal to unity minus the fraction by which the magnitude is reduced via refraction. Our resulting modified magnitude model (A mm ) is We estimate r with our Bayesian framework, but the correction factor can also be computed in the conventional analysis framework by determining the value that minimizes the variance in the computed electron collision frequency. It should be noted that computation within the conventional framework requires careful weighting of the data provided at various THz frequency components -a weighting that is provided automatically in the Bayesian framework. Figure 2 compares the electron collision frequency computed as a function of the radiation frequency with and without correction for refraction. Data used to construct the plot were taken with the RF ICP discharge operating at 200 W and 1 Torr. Correction of the magnitude results in an order of magnitude decrease in the span of the calculated electron collision frequency. The span of computed values among THz pulse frequency components with the largest signal-to-noise ratio (SNR), ranging from approximately 0.6 to 1.2 THz (with a few exceptions caused by water vapor absorption [29]), is most significantly reduced. The apparent trends and scatter in the corrected data outside this range are caused by low SNR. The Bayesian framework enables consideration and appropriate weighting of the data at each frequency component for estimation of the correction factors and plasma parameters, as well as their respective uncertainties.

Bayesian inference and UQ
The parameter estimation problem outlined in Sect. 3 features complex inter-dependencies between THz-TDS data and numerous quantities of interest (QoI), i.e., electron density and collision frequency. Estimates of the QoI are highly sensitive to noise and variables that are difficult to precisely characterize, such as the plasma length and necessary correction factors. Bayesian inference is a statistical framework for parameter estimation that combines prior and measured information from multiple sources to generate accurate estimates of the QoI and make a comprehensive account of the uncertainty about those estimates [19]. We employ Bayesian inference to determine the plasma electron density and collision frequency using THz-TDS data, and to quantify their respective uncertainties.

Overview
Traditional methods for UQ are based on the propagation of standard errors using a Taylor series expansion of the measurement equations [30]. This approach fails when errors are non-Gaussian, the measurement model is highly nonlinear, or regularization is required to obtain accurate estimates of the QoI [31]. The Bayesian framework for UQ accounts for these effects. In Bayesian inference, all variables (the measured data, QoI, and other physical parameters) are conceived of as random variables, characterized by a probability density function (PDF) that reflects one's knowledge of said variables. Narrow distributions represent a high degree of confidence about the value of a given parameter, whereas a diffuse distribution represents ignorance thereof. The variables in question include vectors of the measured data (b), QoI (x), and nuisance parameters (θ). Nuisance parameters are uncertain quantities that affect the inference of x but are not of primary interest, e.g., the plasma length. (See Sect. 4.2 for an overview of b, x, and θ.) The goal of Bayesian inference is to compute the PDF of x subject to a measurement, called the posterior (p(x|b)). The posterior PDF is considered a comprehensive solution to the parameter estimation problem because p(x|b) carries all prior and measured information about the QoI. Bayes' equation provides an expression for the posterior in terms of the likelihood (p(b|x, θ)) and prior (p pr (x, θ)) PDFs: The likelihood describes the chance of observing b for a hypothetical set of QoI and nuisance parameters, accounting for noise and errors in the measurement equations, and the prior encodes one's antecedent knowledge of the QoI and nuisance parameters, i.e., independent of the current measurement. Bayes' equation is normalized by the evidence (p(b)), which ensures that the posterior PDF integrates to unity. The evidence is obtained by marginalizing the numerator of Bayes' equation: Similarly, while the posterior on the left side of Eq. (11) contains the nuisance parameters, we seek the distribution of x that accounts for the distribution of θ, which also is obtained by marginalization: Despite the fact that the posterior PDF is a full solution to the parameter estimation problem, it is often useful to report a singular estimate of x along with a credible interval (CI) that reflects the width of the posterior density. We summarize the posterior of individual QoI using the conditional mean (CM), which is computed for each QoI to construct x CM . Uncertainties about the CM reported in this paper are based on an equal-tailed Pth-percentile CI, [x CM − γ, x CM + γ], where P is a percentage and γ satisfies Using the subjective Bayesian interpretation of this interval, there is a P% chance that x is in the Pth-percentile CI, subject to the data, the distribution of θ, and one's prior knowledge of x. The remainder of this section describes the variables considered in Bayesian THz-TDS, construction of the likelihood and prior PDFs, and computation of the posterior.

Variable definitions
Quantities of interest are the line-averaged electron density and collision frequency, inferred parameters of secondary importance (i.e., nuisance parameters) are the plasma length, delay line registration error, and refraction correction factor, θ = [L, δt reg , r] T ; (17) and the data vector contains the measured frequency-dependent magnitudes and phases of the THz-TDS transfer function, where N is the number of resolved angular frequencies. The phases in b are unwrapped, but the phase and magnitude data are otherwise uncorrected. Finally, the discrete Fourier transform convention used in this paper iŝ where t m are the discrete time values for all m ∈ 1, 2, . . . , M.

Likelihood
The likelihood describes the chance of observing a measurement for a hypothetical set of plasma conditions. This calculation is based on a measurement model (b m ) that relates x and θ to b; in THz-TDS, b m consists of the expressions for A mm and Φ mm from Eqs. (10) and (9) applied to each resolved angular frequency: If the measurement equations perfectly described the plasma and THz-TDS apparatus, then p(b|x, θ) would simply be a Dirac delta distribution centered at b m . However, there are always numerous sources of noise that corrupt b, as well as errors and inadequacies in the measurement model. The relationship between b and b m is expressed in terms of additive errors (e), where the error vector accounts for all discrepancies between the measured and modeled data. Consequently, the likelihood PDF describes the distribution of e, which is considered to be a random variable along with b, x, and θ. This section presents one approach to characterizing errors in THz-TDS plasma diagnostic data. The relations can be modified and extended as necessary for application to other THz-TDS plasma diagnostic scenarios.
We model e using a centered multivariate Gaussian distribution, based on the analyses of Withayachumnankul et al. [32] and Meier et al. [13]. It is important to note that e is centered because corrections to known biases in the measurement equations are incorporated into b m . The resulting likelihood is where Γ is the error covariance matrix. In general, the structure of Γ accounts for the magnitude of errors and correlations between individual errors. However, error covariances in THz-TDS are generally negligible in both the time and frequency domains, so Γ is diagonal with elements σ 2 A,tot for amplitude measurements and σ 2 Φ,tot for phase measurements. These quantities are given by and where σ 2 A,Ê (ω k ) and σ 2 Φ,Ê (ω k ) denote uncertainties in the magnitude and phase of the spectral transfer function, respectively, and σ 2 Φ,t (ω k ) is the phase uncertainty caused by uncertainty in the temporal axis of the time-domain measurement.
The temporal axis uncertainty term is equal to the uncertainty in the time axis multiplied by the resolved frequency. Assuming the probe laser pulse has a Gaussian profile in the time domain with known full width at half maximum (δt FWHM ), the frequency-domain phase error variance (σ 2 Φ,t ) can be approximated as (2) . (25) Uncertainties in the spectrum magnitude and phase caused by signal noise (σ 2 noise ) and drift (σ 2 drift ) are inherently measured as time-domain sample (σ 2 E,s ) and reference (σ 2 E,r ) electric field uncertainties. The necessary frequency-domain uncertainties are computed from the time-domain uncertainties with the expressions developed by Withayachumnankul et al. [32]: and The time-domain errors, in turn, are given by and Drift errors are associated with the uncertainty about E r during the measurement of E s , so the latter is not subject to drift errors. Drift affects one's estimate of the reference signal when E s is recorded. For measurements of a solid target, the delay line can be repeatably scanned to minimize the effect of drift errors [32]. However, repeated scanning is not preferred for plasma targets due to challenges associated with maintaining a perfectly steady plasma. Instead, the delay line is typically scanned just once for each measurement of E, with multiple acquisitions at every delay line position; a method is therefore required to quantify the drift.
To do so, reference measurements are made before and after E s , at times t r,1 and t r,2 , such that the time of the sample measurement (t s ) is in the interval [t r,1 , t r,2 ]. The reference field used to evaluate Eq. (2) is interpolated between these measured references, and the uncertainty due to drift is conservatively estimated to be The uncertainty due to noise (σ 2 noise ) is estimated as the square of the standard deviation of the time-domain electric field at each delay line position. In our experiments, we average together 1000 samples at each delay line position and typically recover a noise standard deviation value on the order of 1% of the peak electric field value.

Prior
The prior incorporates information known about the plasma that is independent of the current measurement. This information can be derived from previous measurements, results reported in the literature, simulations, or fundamental constraints on the QoI and nuisance parameters. In principle, the prior should encode all of one's general knowledge about x and θ. We employ minimal assumptions about these variables to build p pr (x, θ) in this paper, which reflects our limited knowledge of the RF ICP parameters prior to our measurements.
To start, we assume that prior information about the QoI and nuisance parameters is uncorrelated, p pr (x, θ) = p pr (n e ) · p pr (ν) · p pr (L) · p pr (δt reg ) · p pr (r).
It is important to note that the use of independent priors for each QoI and nuisance parameter does not imply that these quantities are physically independent -we merely do not impose correlations between the variables. However, as such relationships are observed and confirmed through repeated testing, they can be incorporated into the prior to reduce the uncertainty of subsequent measurements. We put scale-invariant priors on n e and ν, i.e., and where [n e,min , n e,max ] and [ν min , ν max ] are the assumed ranges of viable electron densities and collision frequencies. These ranges are conservatively estimated as [10 16 m −3 , 10 21 m −3 ] and 2π × [10 6 Hz, 10 13 Hz], respectively. Scale-invariant PDFs are akin to uniform PDFs for variables whose range spans several orders of magnitude, with equal weight assigned to each decibel of the range. The plasma length is modeled using a Gaussian distribution, where µ L is our best estimate of the plasma length (12 cm) and σ L is the corresponding standard deviation of that estimate (1 cm), determined in accordance with the measurement technique. Finally, we use uniform priors for the correction factors, and where δt reg,min and δt reg,max are -1 and 1 ps, respectively, and r min and r max are 0.5 and 1, respectively. Our priors and corresponding parameter limits are conservative and may be further narrowed to better describe specific experiments. Except for the case of the plasma length, we found that prior widths had little impact on estimated uncertainties.

Posterior
The posterior PDF produced by our likelihood and prior is not analytically tractable, so we use a Markov chain Monte Carlo algorithm to sample from the joint density, i.e., p(b, x, θ) = p(b|x, θ)p pr (x, θ). Since b is fixed, p(b, x, θ) is proportional to p(x, θ|b) (see Eq. (11)), and sampling from the joint density is equivalent to sampling from the posterior. We use MATLAB's implementation of Neal's [33] slice sampling algorithm to recover the posterior PDF. To improve numerical stability, we directly sample log-scale PDFs and employ the Cholesky factorization of the error covariance matrix. Marginal PDFs are obtained by simply binning the samples along the dimensions of interest, and our equal-tailed CIs are computed via numerical integration.

Results
This section provides parameter estimations determined via application of the Bayesian framework to measurements of the RF ICP discharge. Figure 3 shows the measured electron density and collision frequency as a function of discharge power for discharge pressures of 0.5 and 1 Torr. At each condition, the CM is indicated by the center line, the 50% CI is bounded by the box, and the 90% CI is contained within the error bars. Faded colors in the density plot correspond to conditions for which the collision frequency was not measurable, i.e., was beneath the collision frequency resolution limit of our THz-TDS system. Parameters are considered below the resolution limit of the system when the relevant marginal posterior is approximately equal to the input marginal prior. In these cases, the measurement uncertainty is sufficiently large that the measurement adds virtually no information about the value of the parameter. Additionally, measurements at the 80 W, 1 Torr condition were highly contaminated by noise due to physical vibrations in our laboratory, producing low SNR. The resulting inference exhibits a high degree of uncertainty that, while plausibly consistent with the observed trends, clearly deviates from the expected value. Our Bayesian framework indicated considerable uncertainty about ν in this case, which illustrates the utility of the UQ procedure. We also note that our reported electron collision frequency values are larger than those predicted by first-order consideration of argon electron-neutral momentum transfer cross sections [34], which almost ubiquitously predicts values less than 10 10 Hz across our operating conditions. We have identified three possible explanations for this discrepancy: 1) our refraction artifact correction scheme is insufficient, 2) reflections at the plasma-vacuum interfaces (which are not considered in our model equations) artificially decrease the measured spectrum magnitude, and/or 3) first-order consideration of electron-neutral momentum transfer cross sections is not sufficient to accurately predict the effective electron collision frequency [35].
In an effort to address the first possibility, we performed a measurement with the plasma discharge intentionally skewed with respect to the THz radiation path to increase the impact of refraction on our electron collision frequency measurements. The resulting spread in uncorrected electron collision frequencies increased by more than an order of magnitude compared to the data shown in Fig. 2. Application of the refraction artifact correction via the Bayesian framework produced an electron collision frequency within 10% of that measured with the plasma discharge operating at the same conditions but without intentional skewing. Though this result does not entirely vindicate the refraction artifact correction, it does not indicate any obvious problems with the correction scheme.
As noted in Sect. 3.1, reflections at the plasma-vacuum interfaces are not included in the model equations due to the small mismatch between the real components of the plasma and vacuum refractive indices. Across all plasma conditions and relevant THz frequency components in this work, the real component of the refractive index is very near unity. As an example, the real component of the plasma refractive index is approximately 0.999 for the 120 W, 1 Torr operating condition. Application of the Fresnel equations at the plasma-vacuum interfaces reveal that reflection-induced changes in the transfer function magnitude are well below the detection capability of the THz-TDS system. Reflections not accounted for by the model equations are therefore not sufficient to explain the discrepancy in electron collision frequency values [8].
Finally, we must consider the possibility that first-order consideration of tabulated argon electron-neutral momentum transfer cross sections is insufficient to accurately predict the effective electron collision frequency. As discussed by Lafleur, et al. [35], simulations with high-fidelity particle-in-cell (PIC) discharge models sometimes find effective electron collision frequencies significantly larger than those predicted by first-order kinetic modeling. Such model results are supported by experimental evidence; unexpectedly large effective electron collision frequencies have been observed in ICP [36], helicon [37], Hall thruster [38,39], planar magnetron [40], and hollow cathode [41] discharges. The discrepancy between our data and the expected values, therefore, is not without precedent and may simply be due to the inadequacies of a simplified model.
Bearing in mind that the electron collision frequency data is unexpected, we still acknowledge the trends shown in Fig. 3. For a given pressure condition, the electron density increases linearly with input discharge power. Moreover, holding power constant, the electron density increases with pressure. The collision frequency data exhibit broadly-similar trends with power and pressure, but these are less pronounced and, given the reported CIs, less conclusive. Collectively, these results may indicate that the marginal energy provided to electrons by additional input power raises the number of ionizing collisions and slightly increases the overall electron collision frequency. Increasing the pressure provides more neutral species, thereby increasing the number of ionizing collisions and, apparently, the total electron collision frequency. Figure 4 shows measured THz pulses alongside the joint posterior of n e and ν, marginalized from p(x|b), for the 120 W, 1 Torr measurement condition. Both the joint and marginal PDFs are Gaussian, and the bivariate PDF does not exhibit any correlation between estimates of the electron density and collision frequency. Based on the Gaussian nature of p(n e , ν|b), Bayesian THz-TDS is a good candidate for the Laplace approximation to the posterior, in which one simply solves for the QoI and nuisance parameters that maximize the posterior, called maximum a posteriori (MAP) estimation, and then fits a Gaussian posterior using the local curvature at the MAP estimate [42]. The resulting Gaussian density can be used to accurately estimate uncertainties about the QoI, which considerably streamlines Bayesian inference and reduces the computational effort required to conduct Bayesian UQ. Though the density and collision frequency are uncorrelated in the posterior PDF, other parameters do exhibit correlations. Figure 5 shows the bivariate marginal posterior PDFs of the density and plasma length (p(n e , L|b)) and collision frequency and plasma length (p(ν, L|b)). A strong, nearly linear correlation exists between the assumed plasma length and resulting electron density. As mentioned in Sect. 4.4, this observation does not indicate a physical correlation between n e and L. Rather, the correlations in Fig. 5 illustrate how our knowledge of these parameters is mediated (and sometimes conflated) by the measurement equations. Information about the plasma length affects our estimate of the number density, and uncertainty about the length has a considerable direct influence on the uncertainty about n e . This is expected, as simplification of the governing equations to measure the electron density from the phase alone produces a linear relationship between the phase and the product of the electron density and plasma length. By comparison, the electron collision frequency is essentially uncorrelated with the plasma length, which is not immediately obvious from the measurement equations. This characteristic was previously noted by Jamison et al. [8], but Bayesian inference provides a simple framework to assess such variable relationships.

Conclusions and outlook
This paper demonstrates the use of a novel Bayesian framework for analyzing THz-TDS plasma diagnostic data. The framework was found to enable computation of both the electron density and collision frequency, as well as their associated uncertainties. The framework also successfully identified measurement conditions for which plasma parameter values were below the resolution floor of the THz-TDS system. Future improvements to the framework could include the use of a sophisticated prior to improve the accuracy of measurements for a specific plasma discharge, higher-fidelity error modeling in the likelihood function, and the extension to complex measurement scenarios not adequately described by the model equations given here.