A new method to measure the unloaded quality factor of superconducting cavities

We describe an method that measures the unloaded quality factor $Q_0$, the external quality factor $Q_E$, and the cavity detuning $\Delta\omega$ with a recursive least-squares algorithm. It combines a large number of consecutive measurements to successively improve an estimate of fit parameters that asymptotically converges to the ``real'' values. Exploiting the large amount of data acquired by a digital low-level RF system permits us to reach this asymptotic regime in a moderate time frame of seconds to minutes. Simulations show that the method works both for critically coupled and over-coupled cavities. A new calibration method adresses very tight tolerances of the method on system parameters.


Introduction
The unloaded quality factor of superconducting cavities, commonly referred to as Q 0 , parameterizes the losses due to the energy dissipated in the walls of the cavity.These dynamic losses constitute a major heat load for the cryogenic system, especially for accelerators operating in continuous-wave mode [1,2,3], because it limits their energy efficiency.Measuring Q 0 is therefore of paramount importance for the successful operation of superconducting accelerating structures and the accelerators that use them.
Measuring Q 0 in critically coupled cavities is usually done in a vertical cryostat and is based on comparing the incident power to the power escaping through the power coupler, the measuring antenna, and via losses in the cavity walls [4].A careful analysis of the associated uncertainties is documented in [5].Those were found to be in the several-percent range.Additional systematic errors, especially due to re-reflections from the circulator are discussed in [6,7].Refined methods, exploiting both amplitude and phase information, are discussed in [8,9].
Measuring Q 0 in over-coupled cavities equipped with high-power couplers, usually found in cryomodules, is more difficult.In particular, because the power escaping through the power couplers is orders of magnitude larger than the power dissipated in the cavity walls or escaping through the measuring antenna.In the past this made it practically impossible to measure Q 0 by analyzing radiofrequency (RF) signals, where [10] is an exception.Therefore, calorimetric measurements are normally used.Either by measuring the increase of the helium pressure [11], if the enclosing helium vessel is closed, or by measuring the flow of evaporated helium [12], if helium cooling is operated in steady state.The situations is aggravated, if multiple cavities are jointly connected to a single helium vessel, which only allows to measure their combined power dissipation.The cited accuracies for the determination of Q 0 are in the 10 % to 15 % level, reached after averaging the helium flow for times measured in hours.
In this report, we extend the analysis from [13] and describe a new method that measures the unloaded quality factor Q 0 , the external quality factor Q E , and the cavity detuning ∆ω for critically coupled and for over-coupled cavities.It uses a recursive least-squares algorithm and combines a large number of successive measurements to improve an estimate of the fit parameters.Affected by process noise only, this type of measurement can be shown to asymptotically converge to the "real" values [14], whereas measurement noise can slightly bias the asymptotic value.We will discuss this further below.Nevertheless, exploiting the large amount of data acquired by a digital low-level RF system permits us to reach the asymptotic regime in a moderate time frame of seconds or minutes.
This report is organized as follows.We first discuss the model and the system identification process in Section 2, before simulating a generic system in Section 3, where we find that, especially for over-coupled systems, the sensitivity to systematic errors is extreme.We address this problem in Section 4 with a new calibration method, specific for this type of measurement, before concluding.

Model and system identification
We describe the superconducting cavity as resonator with resistor R, capacitor C, and inductor L coupled in parallel that is driven by a current I.After averaging over the RF oscillations, the envelope of the voltage across the components is described by V = V r + iV i , where V r is the real part of the voltage (in-phase, I) and V i is the imaginary part (quadrature phase, Q).This dynamical system can be written as [15] where ∆ω is the cavity detuning, ω is the resonance frequency of the cavity, and ω 12 = ω/2Q L is the cavity bandwidth with the loaded quality factor 1 Here Q E is the external quality factor, Q 0 the unloaded quality factor, and Q t ≫ Q 0 is the quality factor of the antenna that measures the field in the cavity.I = I r + iI i is the current that drives the cavity.It is measured by directional couplers that determine the forward current I + = I + r + iI + i and the reflected current , which leads us to [13] dVr dt dV i dt Converting the system to discrete time with time step ∆t results in in the system identification process.We assume it is uncorrelated, has magnitude σ m , and is characterized by Similar to [13] we isolate the dependence on the fit parameters by first writing Equation 3 as with and That ω 0 depends on both Q 0 and Q t indicates that we cannot disentangle their respective contributions using this method.Since normally Q t ≫ Q 0 this should not be a big problem.Otherwise, β t = Q 0 /Q t must be determined by conventional methods [4,5] of the critically coupled cavity.Continuing with the discussion of the system identification process, we write the right-hand side of Equation 5as This allows us to write Equation 5 in the compact form with which has the same structure already encountered in [13].Following in the same vein we write successive versions of Equation 8 in a vectorized form and solve it with the Moore-Penrose pseudoinverse [16] to arrive at where we introduce the vector ⃗ q T of the estimate of fit parameters after T iterations.
The matrix U T rapidly grows to unmanageable proportions such that we resort to a recursive algorithm that solves it one step at a time.We therefore introduce P −1 T = U ⊤ T U T and its initial value P 0 = p 0 1 and express P T +1 through P T in the following way And here the previously mentioned bias originates.Since the matrices G t depend linearly on the voltages ⃗ 4therefore also appears quadratically and on average causes a bias proportional to σ 2 m , which turned out to be insignificant in all simulations reported below.
It turns out that retrieving P T +1 from inverting P −1 T + G ⊤ T +1 G T +1 can be efficiently done with the Woodbury matrix identity [17] With the substitutions A −1 = P T , V = G ⊤ T +1 , and W ⊤ = G T +1 , we obtain We now continue to find ⃗ q T +1 by writing Equation 11for T + 1 Equation 14 and 15 constitute the algorithm to continuously update estimates for the three components of ⃗ q, the bandwidth q(1) = ω E ∆t, the detuning q(2) = ∆ω∆t, and q(3) = ω 0 ∆t as new voltage and current measurements become available.
Like in [13] we introduce a "forgetting factor" α = 1 − 1/N f , where N f is the time horizon over which old measurements are discounted.This is of limited practical value here, because it limits the achievable accuracy of the method.It might, however, become useful to investigate the stability of the measurements and whether some unaccounted external factors cause the fit parameters to vary over time.To avoid the limitation of the accuracy in this report, we choose N f = 10 6 in all simulations, a value much larger than the number of iterations used.With α included, Equation 12 reads P −1 and Equation 14 becomes whereas Equation 15turns into In these equations the calculations are moderately time-consuming, because the involved quantities are matrices.In particular, P T is a 3 × 3 matrix, G T +1 is a 2 × 3 matrix, and the inversion involves the 2 × 2 matrix α1 + G T +1 P T G ⊤ T +1 .

Simulations
In the simulations we assume a generic superconducting cavity with a resonance frequency ω/2π = 10 9 Hz and Q ′ 0 = 10 9 .As in [13] we scale the voltages and currents to their maximum values and denote the corresponding real and imaginary values by v r , v i , i r , and i i respectively.For the process noise we assume σ p = 10 −4 × V max and for the measurement noise σ m = 10 −3 × V max .Once the cavity reaches a steady state, the first and last column of G t in Equation 9 become linearly dependent.We therefore have to pulse the generator that excites the cavity.In all simulations we turn it on and off every 1000 iterations.We refer to the MATLAB [18] code on github [19] for the details of the implementation.We first consider a close-to-critically-coupled system and assume Q E = 8×10 8 , which results in a bandwidth of f 12 = ω 12 /2π = 1.125 Hz.For the detuning we assume ∆ω = ω 12 /2.Since both ω E and ω 0 are very small we chose a moderately slow sampling rate of 1 kSample/s or ∆t = 10 −3 s.If the LLRF system operates at a higher rate, low-pass filtering and subsequent undersampling the output is beneficial because it reduces the measurement noise and reduces the bias mentioned before.
The top-left part in Figure 1 shows the normalized voltages and currents for 10000 iterations, which corresponds to 10 s in real time.The pulsing cur- rents (lower panel) and the responding voltages (upper panel) are clearly visible.The top-right plot shows f E = ω E /2π, ∆f = ∆ω/2π, and f Q 0 = ω 0 /2π as the simulation progresses.We observe that these fit parameters approach constant values after a few thousand iterations.The lower-left part shows Q E = ω/ω E and Q ′ 0 = ω/ω 0 .Both quantities approach the correct values on the same time scale as the frequencies.The part on the lower right shows the relative uncertainties σ(Q)/Q (the "relative error bars") derived from the diagonal elements of the empirical covariance matrix P T .For example the error bar of Q ′ 0 is approximately given by σ(Q ′ 0 ) = P T (3, 3)σ m .We find that these statistical uncertainties quickly approach the percent level.
We evaluate the sensitivity to systematic uncertainties by either changing the resistance R or the calibration scale factor to determine ⃗ V ′ T by 5 % when producing the data, thus causing a systematic mismatch between the assumed and the "real" model.We find that these changes also affect the fitted values of Q E and Q ′ 0 by about 5 %.We conclude that the method is reasonably robust in this regime of operation.
Next we consider an over-coupled cavity with Q E = 10 6 , which results in a bandwidth of f 12 = ω 12 /2π = 500 Hz.Again we use ∆ω = ω 12 /2 in the simulations.Since Q E is much lower, the system reacts much more speedily and we choose a smaller value of ∆t = 10 −4 s.The simulation extends over 6 × 10 4 iterations which corresponds to six seconds real time.
The four parts in Figure 2 follow the same order as in Figure 1 where the top left shows the pulsing voltages and currents and the top right shows the fit parameters with ω 0 /2π ≈ 1 Hz and ω E /2π ≈ 1000 Hz due to the large difference between Q E and Q ′ 0 .These frequencies lead to the fit values for Q ′ 0 and Q E shown on the bottom left.Within about 2 × 10 4 iterations they approach the values used to generate the data.The corresponding relative uncertainties, shown on the bottom right, are 2 % and less than 10 −3 , respectively.We conclude that the algorithm actually recovers the quality factors with good resolution, despite their three-order-of-magnitude difference in magnitude.
The sensitivity to calibration uncertainties of the resistance R and voltage measurement, on the other hand, is extreme.Scale errors, off by 5 × 10 −4 , either of the shunt impedance or the voltage measurements scale almost double the fitted value of Q ′ 0 .Repeating the sensitivity analysis for a cavity with Q E = 10 7 we find that ten times larger scale errors of 5 × 10 −3 double the fitted value of Q ′ 0 , whereas in a cavity with Q E = 10 8 scale errors of 5 × 10 −2 double Q ′ 0 .This extreme sensitivity to systematic errors warrants a discussion of methods to calibrate the hardware.

Calibration
We emphasize that all hardware must be meticulously calibrated following the discussion in [4,5] as well as calibrating the propagation delays by de-embedding [20] the ancillary hardware.
Once this is accomplished we turn our attention to the new algorithm and note that it depends entirely on the definition of the matrix G t from Equation 9. We observe that any scale factors of currents or voltages can be absorbed in the resistance R, which serves as a factor to make currents and voltages commensurate.The systematic errors discussed at the end of the previous section were caused by such an imbalance between the currents and voltages that we can compensate by adjusting R. For example, simultaneously increasing all voltages ten-fold and, at the same time, increasing R ten-fold, leaves the fit-parameters unaffected.This is easy to see, because both G t and ⃗ y t+1 from Equation 9 are multiplied by ten, but this common factor cancels, once ⃗ q T is calculated in Equation 11.Being the only free parameter in the algorithm, we therefore have to determine R to the 10 −4 level.
We determine R from Equation 1 in steady state (dV /dt = 0) and onresonance (∆ω = 0).With input currents at a constant value, the steady state is typically reached after a few seconds.If the data acquisition system runs with a high sampling rate, the algorithm from [13] can validate that ∆ω is indeed zero.Under these conditions Equation 1 simplifies to Note that here we use the total current ⃗ I = ⃗ I + + ⃗ I − instead of the forward current ⃗ I + .We thus need to add both channels from the directional coupler normally used to measure the forward current only.By inspecting Equation 18 we find that the bandwidth ω 12 cancels, leaving us with In the latter equation we added the prime to the voltages to indicate that their measurement is affected by measurement noise of magnitude σ m .Moreover, we explicitly added the subscript t to denote that even in steady state the noise causes the values to vary with time.Thus we can also repeat measurements to improve the accuracy of determining R. Equation 19 has indeed the same structure as Equation 8with ⃗ V ′ t taking the role of ⃗ y t+1 , the resistance R taking the role of ⃗ q T , and ⃗ I t taking the role of G t .We can therefore use the same recursive least-squares algorithm already described in Section 2, especially Equations 16 and 17, to determine R. The corresponding equations now read and where q T is the steadily improving estimate of the resistance R. Again, we refer to [19] for the details of the implementation in MATLAB.
The left-hand plot in Figure 3 shows the difference between the estimated value of R and the true value, often referred to as estimation error, as the calibration progresses for 10 6 iterations after steady state was reached.Here we do not use any forgetting factor and set α = 1.We observe that the estimation error steadily decreases below the 10 −4 level that is required to determine the unloaded quality factor in over-coupled cavities.The right-hand plot shows the empirical covariance matrix P T as the calibration progresses.It steadily decreases as well.From P T we can obtain the approximate error bars of R by evaluating √ P T σ m .We thus find that the described calibration procedure indeed allows us to determine R with high precision.All systematic scale factors are then incorporated in R as long as the hardware used in the calibration is the same one used to measure Q ′ 0 .

Conclusions
We described a new method to determine the unloaded quality factor of superconducting cavities.The method works both for critically coupled and over-coupled cavities, but is very sensitive on systematic calibration errors.This sensitivity can be absorbed into the resistance R that can be calibrated with high precision.
Both the identification of the quality factors and the calibration of R employ recursive least-squares algorithms, which are known to converge towards the "real" values.
We have, however, to keep in mind that the method relies on the fit parameters being constant while the measurement is ongoing.This can be verified, provided the data acquisition system and the analog-to-digital converters sample at a high rate.Low-pass filtered and undersampled signals are passed to the algorithm described in this paper, while the high-speed signals are passed to the algorithm described in [13].The latter can then be operated in parallel to determine the detuning ∆ω and the bandwidth ω 12 at high speed to ensure that the parameters are indeed constant.A rapid increase in ω 12 would indicate a quench, while microphonics or Lorentz force detuning show up as variations of ∆ω.
In particular at very high power levels, the Lorentz force detuning could become a problem, especially due to pulsing the input current.This might limit the new method to lower power levels.Optionally, other, more gentle excitation patterns, for example by amplitude-modulating the input power with a sinusoidal variation, can be explored.This would still ensure the algorithm to work, though possibly at the expense of slower convergence.Optimizing these methods is, however, outside the scope of this report.
The need to pulse or modulate the input power interferes with continuouswave operation of an accelerator, but only very briefly, because the time needed to achieve reasonable accuracies is very short, only a few seconds or minutes are needed.This could either be scheduled during dedicated machine development shifts or during brief interruptions of regular operations.

Figure 2 :
Figure 2: Simulation of an over-coupled cavity: voltages and currents (top left); fit parameters (top right); Q E and Q ′ 0 (bottom left); relative uncertainties (bottom right).

Figure 3 :
Figure 3: The estimation error (left) during the calibration of R and the empirical covariance matrix P T (right) as a function of the number of iterations.