Model and numerical simulation of MKIDs modelling This chapter describes a
study to model two main figures of merit that can describe the performances of
an array. The fabrication yield is defined as the number of resonators,
expressed as a percentage of the total, with individually-distinguishable
resonance frequencies. This represents the amount of pixels that can be used
to reconstruct an image. Two resonators with the same resonance frequency, or
with frequencies that are too close to each other, are inherently
indistinguishable, therefore it is impossible to trace back which MKID
generated the signal. The second parameter of great interest is the energy
resolution of each single resonator. The energy resolution describes the
smallest possible difference in energy of the impinging photons, $\Delta E$,
the detector can identify. A more in-depth description of these figures of
merit is given in the next sections along with the description of a model that
simulates the behaviour of MKIDs in order to predict, optimise and design the
best array possible. Fabrication Yield The fabrication of the detectors for
such instruments is performed in research-grade clean rooms, through processes
which are almost completely human-controlled and not automated. These
processes have plenty of room for optimisation, although the fraction of non-
responding pixels is less than $1\,\%$. Two or more resonators that resonate
at the same frequency while they were not designed for this purpose, are said
to be clashing or colliding. Currently, it is thought that most of the clashes
are due to random variations of the critical temperature of the
superconductor. Although a systematic radial distribution, over a large area ,
of critical temperature is known to arise during TiN${}_{\text}{x}$ sputtering
deposition Vissers_2013, it only plays a negligible role for detectors of chip
size is $\approx$ $4$ cm2 as its effect can be dramatically mitigated with an
accurate placement of the resonators on the test mask that aims to even out
possible drifts in critical temperature and film inhomogeneity . The current
state of the art is that most arrays exhibit a fabrication yield of $75\,\%$
and only for the best ones it goes as high as $80\,\%$ mazin2020. The model
here discussed is based off the approach described by X. Liu et al in
yield_liu. This simple model relies on a few basic assumptions:

  1. 1.

the array of resonators is infinite on both sides of resonator $\\#0$

  2. 2.

the distance, in frequency space, between two adjacent resonators is designed
to be: $f_{n+1}-f_{n}=\Delta=2$ MHz

  3. 3.

the n-th MKID’s resonance frequency is randomly distributed about $\Delta n$
according to a gaussian distribution
$G(f_{n},\sigma)=\frac{1}{\sqrt{2\pi}\sigma}e{{-\frac{1}{2}\left(\frac{(f_{n}-n\Delta)}{\sigma}\right)^{2}}}$

  4. 4.

given a defined window, a resonator that falls within said window is defined
to be clashing.

modelling/model.png

Figure : Representation of the simple model of the MKIDs array. Resonances
$\\#0$ and $\\#\pm 1$ are shown as well as the collision window. The
probability of two resonators colliding is given by the area underneath the
bell curve within the limits of the collision window.

The probability that the $0$th resonator does not clash with the nth
resonator, Pn0, can be evaluated, according to [yield_liu], as one minus the
overlapping area defined by the gaussian distribution and the collision window
as described in Figure 1. Furthermore, due to the translation symmetry of the
model, it is only necessary to evaluate Pn0 for all values of $\pm n$. The
benefit of an infinite array is that allows the following identification:

$P_{n0}=1-\frac{Erf\left(\frac{n\Delta+\chi
w}{\sqrt{2}\sigma}\right)-Erf\left(\frac{n\Delta-\chi
w}{\sqrt{2}\sigma}\right)}{2}$ (1)

Where the integral of a gaussian is evaluated through the error function
Erf(x) defined as in Equation 2

$Erf(x)=\frac{1}{\sqrt{\pi}}\int\limits_{-x}^{+x}e^{-t^{2}}dt$ (2)

Here, $\Delta$ is the nominal spacing between the resonators, and $\chi w$
identifies the collision window through two parameters: $w$ represents the
resonator’s line-width and a $\chi=5$ is assumed under the hypothesis that the
resonators are well distinguishable when separated by $5$ line-widths. This is
to say that the profile of the $0$-th resonator is well contained in a window
five times larges than its typical line-width. Assuming all the resonators to
be independent, the total yield is derived as:

$Yield=\prod\limits_{n=-n_{max}}^{n_{max}}P_{n0}$ (3)

For the sake of a model, it was assumed that the MKID $\\#0$ resonates at the
central frequency of the working bandwidth, $6$ GHz. And only the first $1000$
resonators left and right have been taken into consideration. The
identification 3 is only valid in the case of an infinite array of evenly
spaced resonators. The analysis of this model will be carried out in the next
few pages and the dependencies of the yield from physical parameters such as
the quality factor of the resonators and the critical temperature of the
superconductor will be investigated. It is possible to observe that while
increasing the spacing $\Delta$ between two consecutive elements of the array,
the yield increases accordingly.
A question left unanswered, this far, is how the frequency of an ideal
resonator changes when the desired Tc is not met. A simple approach follows
the discussion described in [Annunziata] and [Tinkham]; the dependence of the
kinetic inductance as a function of the temperature can be evaluated as:

$L_{k}=\frac{l}{w}\frac{R_{s}h}{2\pi^{2}1.76T_{c_{0}}}\tanh\left(\frac{1.76T_{c_{0}}}{K_{B}T}\right)$
(4)

Where $R_{s}$ is the resistivity of the film in its normal state, l and w are
the physical length and width of the film, h is its thickness, $T_{c0}$ is the
critical temperature of the film and $k_{B}$ is Boltzmann’s constant.
Given the change in kinetic inductance, and knowing that the resonance
frequency of an LC circuit $f_{r}=\frac{1}{\sqrt{LC}}$ the fractional
frequency shift in resonance frequency induced by a change in critical
temperature, therefore a local change in kinetic inductance, is easily
obtained provided three assumptions are met :

  1. 1.

The capacitance of the resonator is only due to its geometry

  2. 2.

The inductance is only due to the kinetic inductance. The kinetic inductance
fraction $\alpha=1$

  3. 3.

Any change in sheet resistivity Rs with Tc is negligible.

Under these circumstances, the fractional frequency shift is given by Equation
5

$\frac{\Delta
f}{f_{r}}(Tc)=\frac{\frac{1}{\sqrt{L(T_{c})C}}-\frac{1}{\sqrt{L(T_{c_{0}})C}}}{\frac{1}{\sqrt{L(T_{c_{0}})C}}}$
(5)

Here the previous assumptions take a critical part, as
$L(T_{c})=L_{K}(Tc)+L_{g}$ where the latter term represents the geometric
inductance.The hypothesis that $\alpha\approx 1$ allows to neglect $L_{g}$
with respect to $L_{K}$ and the assumption that $R_{s}$ does not change with
stoichiometry allows to simplify Equation 5 as:

$\frac{\Delta
f}{f_{r}}(Tc)=\frac{\sqrt{\frac{1}{T_{c}}\frac{1}{\tanh{\left(\frac{1.76T_{c_{0}}}{2k_{B}T}\right)}}}^{-1}-\sqrt{\frac{1}{T_{c}{}_{0}}\frac{1}{\tanh{\left(\frac{1.76T_{c_{0}}}{2k_{B}T}\right)}}}^{-1}}{\sqrt{\frac{1}{T_{c}{}_{0}}\frac{1}{\tanh{\left(\frac{1.76T_{c_{0}}}{2k_{B}T}\right)}}}^{-1}}$
(6)

Where $\Delta f=f_{r}(Tc)-f_{r}(T_{c_{0}}),$ T${}_{c}{}_{0}$is the desired
critical temperature, $800$ mK for the sake of this discussion; T is the
equilibrium temperature at which the MKIDs array is operated, $100$ mK. Figure
1 shows the absolute value of Equation 6 as a function of the temperature
difference $T_{c}-T_{c0}$.
The main purpose of this calculation is to try and estimate reasonable
boundaries for the width $\sigma$ of the distribution of the frequencies about
their mean value $f_{n}$. Although, Equation 5 is asymmetrical around Tc0 ,
for small $T_{c}$ variations, the two branches can be considered symmetrical
as a fist order approximation. It is assumed that only $0.1\%$ of the
resonators have a critical temperature that differs by more than $200$ mK,
thus a shift in resonance frequency that is greater than $1.2$ GHz. This
justifies the assumption of an upper limit for $\sigma=0.300$ GHz. Such a
broad distribution would be catastrophic for any real array and is certainly
overestimate of any realistic value of $\sigma$ in any well-behaved array .

[width=]modelling/Fig2.png

[width=]modelling/Figure33.png

Figure : Left: fractional frequency shift as a function of critical
temperature fluctuations. Right: yield as a function of $\sigma$, evaluated
for a collision window of $0.5$ MHz. The inset shows a wider span of values of
$\sigma$, up to $0.6$ GHz.

A first evaluation of the fabrication yield is carried out to understand how
well this simple model describes the collective behaviour of an array of
MKIDs. A particular choice of the free-parameters is made; it is worthwhile
describing the fabrication yield as a function of a varying $\sigma$ in order
to justify how big an effect our previous assumptions have, and instead of
evaluating the line-width of the resonators, the current working assumption is
that the total width of the collision window is $0.5$ MHz.
Figure 1 (Right) shows the yield starting to roll off as soon as the standard
deviation of the distribution of the resonance frequency around its centre
value is comparable to the spacing between two adjacent resonators
$\Delta_{r}$. The yield eventually flattens for values of $\sigma\geq 0.2$ GHz
at about $75\%$. Unfortunately, due to its inherent limitations, it is quite
difficult to compare the results produced by the model with experimental
results. Nonetheless, it is important to also include in the dissertation how
the total quality factor $Q_{tot}$ of the resonators comes into play in this
model. By definition, the higher the quality factor $Q$, the sharper is the
resonator’s dips Q is inversely proportional to the line-width
$w=\frac{f_{0}}{Q}$ where $f_{0}$ is the characteristic resonant frequency.
Assuming a $\Delta_{r}=2$ MHz, and based on formula 3, Figure 1 shows a 2D
plot of the fabrication yield as a function of both $\sigma$ and $Q$. Two
trends are visible in the plot. The yield increases with increasing quality;
this says that the sharper the resonator, hence the smaller their line-width
and the clashing window, the higher the fabrication yield is. Also, as
expected, the smaller the width of the frequency scatter distribution, the
higher the yield per fixed values of $Q$.
Very strongly coupled resonators with low quality factors, such as
$Q_{tot}<15000$, represent the case where the resonator’s rejection window is
intrinsically larger than the spacing $\Delta_{r}$ reducing dramatically the
yield to a mathematical zero based on previous assumptions. For very high
values of $\sigma$, the model is, in fact, describing an almost random
distribution of resonant frequencies. It therefore makes sense that the yield
becomes proportional to the rejection window size, thus to the total quality
factor $Q_{tot}$.

[width=]modelling/YIELD_25032020.png

Figure : Fabrication yield as a function of the overall quality factor $Q$ and
the width of the frequency scatter distribution, $\sigma$.

It is crucial to state the flaws and the limits of this model and discuss the
possible improvements of said picture. There are two main issues to discuss:
first of all, a gargantuan value of $\sigma$ would produce a $100\%$ yield as
it would scatter most of the resonators outside the working-bandwidth,$4-8$
GHz. It was observed that the yield starts increasing from $\sigma=1$ GHz
onward. Another limit of the model is that, when evaluating the probability of
a collision, only one of the two resonator’s dips is broadened, while the
other is still regarded as if it was infinitely sharp. This effect is only
relevant for very strongly coupled resonators, $Q_{tot}>\frac{1}{\sigma}$,
therefore it is not significant for the plot in Figure 1. Further improvement,
which is a higher-order approximation to what has been discussed so far
resides in the evaluation of the line-width as a function of the specific
resonator’s resonance frequency. In order to summarise the results described
in the previous pages, Table 1 represents some numerical results for
interesting parameters.

σ (MHz) | Q_tot | Yield (%)
---|---|---
0.8 | 20000 | 79.8
0.6 | 20000 | 89.9
1.8 | 30000 | 73.8
1.4 | 30000 | 85.5
1.0 | 30000 | 95.7
2.4 | 50000 | 75.8
1.8 | 50000 | 87.8
1.4 | 50000 | 95.3
4.2 | 100000 | 74.4
2.8 | 100000 | 84.0
1.8 | 100000 | 94.8
Table : Table representing the yield of a 2000 pixel array as a function of
$\sigma$ and the total quality factor $Q_{tot}$

##  1 Resolving Power

For detectors, such as MKIDs, that can measure the energy of the incident
radiation, energy resolution is a critical figure of merit. The best way to
describe it is, borrowing Dr. W. R. Leo’s words, ” […] the extent to which a
detector can distinguish two close lying energies” [WRLeo]. An ideal detector,
when illuminated by monochromatic light would produce a spectrum that is also
monochromatic, and hence described by a Dirac’s $\delta$ function. In reality,
the width of the spectral lines is finite, broadened by several factors which
lead to what is experimentally observed. For the sake of clarity, the energy
resolution of a detector is intended as $\frac{\Delta E}{E}$, The quantity
that will be evaluated most of the time is its inverse $\frac{E}{\Delta E}$ is
referred to as the Resolving Power. At present, deVisser et al.
[deVissersphononiccrystal] claim that the best resolving power achieved
corresponds to $\frac{E}{\Delta E}\sim 55$ for wavelengths between $402-1545$
nm. The next few page will describe how the following effects come into play
in broadening the spectral peak:

  1. 1.

Statistical processes in the MKID

  2. 2.

Amplifier noise - High Electron Mobility Transistor

  3. 3.

Current density inhomogeneity

  4. 4.

Readout frequency

  5. 5.

Two-level systems

###  Energetic processes in the MKID

As stated in [Guo], the best energy resolution an MKID can achieve as a pair-
breaking detector is given by $R=\frac{\Delta
E}{E}={2.355}\sqrt{\frac{F\Delta_{s}}{\eta h\nu}}$, where $\eta\approx 0.57$
is the pair-breaking efficiency [EnergyDownconversion], $h\nu$ is the energy
of the impinging photon, $\Delta_{s}=1.72k_{B}T_{c}$ is the superconducting
gap energy of the adsorbing material. $F$ is the Fano factor, the effect of
which encompasses all the statistical processes that do not lead to pair-
breaking. A typical value of the Fano factor for superconductors is $F=0.2$.
These parameters yield a maximum energy resolution of $0.011$ for
$\lambda_{0}=650$ nm. Unfortunately, there is little room for improvement. For
other kind of energy resolving detectors, such as gas and silicon based
ionisation detectors, the Fano factor is known to vary, quite significantly,
with the material according to Fano_material,fano_gas1. The Fano factor is an
inherent property of the material and perhaps it is worth investigating
whether different superconductors exhibit different Fano factors. the pair-
breaking efficiency, which compares the average number of quasiparticles
produced per photon with the maximum number of quasiparticles that can be
produced, can be marginally improved by changing the design of the resonator
[pairbreakingefficiency]. Studies carried out in SRON by deVisser et al.
[deVissersphononiccrystal] suggest that hot-phonon loss accounts for some of
the losses. In particular, hot-phonon loss contributes to increasing the Fano
factor by a factor $J$.

$R=\frac{\Delta E}{E}={2.355}\sqrt{\frac{(F+J)\Delta_{s}}{\eta h\nu}}$ (7)

An analytical model and experimental measurements demonstrated that for both
Si and SiN substrates, $J\approx 3$, whereas it is higher for Sapphire
substrates. The same paper [deVissersphononiccrystal] proposed the
implementation of MKIDs where the inductors are fabricated on a membrane so to
reduce phonon losses. This should reduce the energy loss through phonon
mediated processes to the silicon substrate, effectively eliminating phonon
loss and thus maximising the pair breaking efficiency. Unfortunately, a value
of $\eta=1$ cannot be achieved since any photon of energy greater than
$2\Delta$ but smaller than $4\Delta$ will only create $2$ quasi particles and
phonon excitations that can not break up any further amount of Cooper Pairs.
For higher energies, the efficiency depends on phonon trapping factor
[pairbreakingefficiency]where improvement is possible.

###  High Electron Mobility Transistor noise

High Electron Mobility Transistors are likely one of the main limiting factors
when it comes to Energy Resolution. The HEMT discussed here is the Low Noise
Factory LNF-LNC4$\\_$8C which is a ultra-low noise cryogenic amplifier
operating in the $4$-$8$ GHz bandwidth. Its nominal noise temperature, as
stated by the manufacturer, is $2.1$ K.
The HEMT random noise affects the phase signal. In particular, when a photon
strikes the MKID, the phase signal generated is added to this noise and, of
course, the height of the phase-jump has an uncertainty given by the noise
contribution of the HEMT itself. Referencing Figure 1, the quantity that has
to be evaluated is $\mathcal{R}=\frac{\Delta\alpha}{\alpha}$. Assuming that
the MKID taken into consideration is optimised to respond with an $\alpha=\pi$
phase shift when struck by a photon of wavelength $\lambda_{0}=650$ nm. Figure
1 shows the angles just discussed and is of critical importance for the
understanding of the model.

modelling/model_hemt.png

Figure : Schematic of a resonance loop. The large circle describes the cable
delay, the smaller circle is that characteristic of a resonator. Two angles
are highlighted: $\pi$ and $\Delta\alpha$. The latter describes random phase
fluctuations due to the HEMT noise. For the sake of clarity, the centre of the
IQ plane was shifted away from the centre of the loop.

HEMT noise induces random fluctuations in both I and Q. Therefore, a non
illuminated MKID produces a signal in the I-Q plane that is represented as a
scattered cloud. When a photon is detected, the cloud is added to the photon
induced phase signal generating inaccuracy regarding the photon’s phase signal
$\Delta\alpha$ . To clarify, the quantity to evaluate is $\Delta\alpha$, the
angle defined by the radius of the resonant loop and the amplitude of the
phase noise. In particular, it is not trivial to understand how the noise
temperature of the amplifier is then translated into phase noise.
Jonas Zmuidzinas’ 2012 review [Zmuidzinas] helps; assuming that the amplifier
noise is equally distributed between the $I$ and the $Q$ component the noise
component can be described by the following equation:

$HEMT\\_noise=\sqrt{\frac{k_{B}T_{A}G}{2P_{MKID}}BW}$ (8)

A thorough derivation of the equation is given in the Appendix
LABEL:appendix:HEMTnoise where $k_{B}$ is Boltzmann’s constant, TA is the
noise temperature of the amplifier, $G$ is the gain of the amplifier,
$P_{MKID}$ is the power at the HEMT and $BW$ is the bandwidth given by the
Fast Fourier Transform sampling rate, $1$ MHz. Equation 8 describes both the I
and Q components of the HEMT noise in I and Q units. Now, geometrical
considerations, based on maximum and minimum transmission (S21), lead to a
simple formula for $\Delta\alpha$

$\Delta\alpha=\frac{HEMT\\_noise}{radius}$ (9)

Where $HEMT_{noise}$ is the length of the arc described by the cloud, and
$radius$ is the radius of the resonant loop. Assuming perfect transmission in
the cable delay loop, which is described by a circle of unit radius. Knowing
the distance from the origin of the resonance frequency of a non illuminated
MKID is $min|S_{21}|$, and by definition of $Q_{c}$ [Zmuidzinas],
$Q_{C}=\frac{Q}{1-min|S_{21}|}$ , it is possible to obtain

$\Delta\alpha=6\frac{HEMT\\_noise\,\,Q_{c}}{Q}$ (10)

Where $Q_{c}$ is the coupling quality factor, $Q$ is the total quality factor
of the MKID. The value of $\Delta\alpha$ would only describe the rms of said
quantity. There is an additional factor of $3$ to account for the conversion
from a rms quantity to a peak-to-peak quantity.

[scale=0.25]modelling/HEMTNOISE_Power.png

Figure : HEMT noise, $\frac{\Delta\alpha}{\alpha}$ as a function of the power
at the resonator. Standard operating powers are in the $-100$dB range.

While it may seem that the noise component is very low, it is important to
observe that usual readout power, before entering the cryostat, would be
$\approx-50$ dBm, in addition to that, two $20$ dB attenuators have to be
taken into account as well as the attenuation of the coaxial cables going into
the cryostat. Assuming a total power at the HEMT of about $-90$dBm,
$Q_{i}=100000$ and $Q_{c}=45000$, $\frac{\Delta\alpha}{\alpha}\approx 0.037$.
It corresponds to a noise fluctuation of $\pm 3.3^{\circ}$ in good agreement
with common measurements of $\pm 5^{\circ}$ white noise in baseline.

[scale= 0.3]modelling/HEMTNOISE_mammeta.png

Figure : 2D graph showing the dependency of the HEMT noise in degrees,
assuming that a photon produces a $180^{\circ}$ signal, as a function of the
coupling quality factor $Q_{c}$ and the internal quality factor $Q_{i}$. It is
important to observe that the area shaded in yellow,represents the parameter
space where the HEMT noise is equal o greater than $36^{\circ}$.

As a final remark, figure 1 shows a 2D plot which clarifies the dependence of
$\Delta\alpha/\alpha$ on both the internal quality factor $Q_{i}$ and the
coupling quality factor $Q_{c}$. The plot here shown refers to a power value
$P_{g}=-90$ dBm. The area shaded in yellow depicts a portion of the parameter
space where $\Delta\alpha/\alpha\geq 0.2$, or $\pm 36^{\circ}$.

###  Current density inhomogeneity

Another source of uncertainty in the response of an MKID is inhomogeneity of
the current density across the inductor. If the current density varies
spatially across the inductor, the resulting signal from an absorbed photon
will depend on absorption location as variations in current density lead to
changes in responsivity $d\theta/dN_{qp}$; where $\theta$ is the phase angle,
the measured response of the MKID to a photon strike, and $N_{qp}$ is the
number of quasi-particles produced by the photon. Current density uniformity
strongly depends on the resonator’s geometry. Current crowding at right angles
is a well known effect for both normal metals and superconductors [current-
crowding]. The current density is affected by the capacitance of the
interdigitated capacitor. The charge carriers tend to slow down and create
high electric field as they drift towards the two plates of the capacitor
resulting in a small, but non negligible gradient in current. The current
density is highest at the top of the inductor (far from the capacitor) and
sees a typical $5\%$ reduction in proximity of the capacitor. This is because
the change in density of the super current induced by a photon depends non-
linearly on the intensity of the current. This reflects its effect on both the
shift in resonance frequency and the total quality factor $Q_{tot}$. The
change in $Q_{tot}$ accounts for an internal change in $Q_{c}$ which depends
on the resonant frequency of the resonator, while a more significant change in
$Q_{i}$ can be attributed to losses through normal-state electrons produced by
pair breaking. In order to evaluate the effects of current density
inhomogeneities, the current density data, simulated with SONNET with a grid
size in both $x$ and $y$ direction of $0.125\,\mu$m (Figure 1) is imported to
a Python script I developed.The data is read in by the programme and it is
then normalised to the average value over the whole inductor. This creates a
current density map which is, for the most part, in the range [$0$,$1$] with
values greater than $1$ only in the inner part of right-corners, where there
is a significant current crowding. The working assumption for the next steps
is that each single element of the inductor, as defined by the simulation
through the grid size, contributes to the total kinetic inductance of the
resonator with its elemental contribution and the total kinetic inductance is
given by the the sum of the contributions of all the single elements. The
resonator is designed in such a way that its resonant frequency is $f_{0}=6$
GHz and $Q_{c}=30000$.

[scale=0.2]modelling/Figure_1.png

[scale=0.2]modelling/Figure_6.png

Figure : Left: Current density map of the inductor. Right: Phase shift when
struck by a photon as a function of the ($x$,$y$) position.

The Python programme, was fed the current density data produced by SONNET,
translates the current density 2D-map first into a Kinetic Inductance 2D-map,
by knowing the inductance depends on the current density at the location of
each of the cells. Once the kinetic inductance of the whole circuit is known,
its resonant frequency is immediately given by $\tfrac{1}{\sqrt{L_{K}C}}$.
When a photon strikes on a specific point of the inductor, it depletes the
region around it of charge carrying Cooper Pairs increasing the kinetic
inductance of the resonator. As a consequence of this increase of kinetic
inductance, the aforementioned shift in resonance frequency , and phase,
occurs. Figure 1 (Left) shows the expected resonance frequency shift according
to where the photon strikes on the inductor. Here, for the sake of simplicity
the diffusion of the quasi-particles in the superconductor was not taken into
account.
It is possible to evaluate the phase shift registered at the readout
frequency, $6$ GHz for the sake of this simulation, and it can be iterated for
all the possible grid elements produced by SONNET. Manipulating equation
LABEL:eq:s21 to obtain the phase component, it is possible to produce the 2D
map shown in Figure 1 (Right) where it is shown the phase shift, measured at
$f_{0}=6$ GHz, due to an impinging photon as a function of the position where
the photon hits the inductor. In order to simulate a well optimised resonator
at the impinging wavelength, the depletion considered was that which produces
a $\geq 0.95\,\pi$ phase shift.

[scale=0.6]modelling/profile1.png

Figure : Histogram of Figure 1. Where $\alpha$ is $\approx 0.965\,\pi$
radians. The designed resonant frequency, in dark, is $6$ GHz and the total
quality factor is $Q_{tot}=30000$.

As previously discussed in Section 1, the Energy Resolution $R=\frac{\Delta
E}{E}$, can be described in terms of $\frac{\Delta\alpha}{\alpha}$. Figure 1
(Left) represents the frequency histogram of the phase values in Figure 1
(Right). The full width of the profile represents $\Delta\alpha$, whereas the
peak value is representative for $\alpha$. Under these circumstances and,
based on these assumptions, for $f_{0}=6$ GHz and $Q_{tot}=30000$,
$\frac{\Delta\alpha}{\alpha}=0.004$. As a final remark, it is interesting to
evaluate if and how $\frac{\Delta\alpha}{\alpha}$ depends on the total quality
factor $Q_{tot}$. A plot of such dependence is shown in Figure 1. This model
hints that a small improvement can be obtained by increasing the resonator’s
total quality factor.
In reality, the model does not include the effects of dissipative effects
induced by the excited quasi-particles, which behave like normal-state
electrons. These losses account for a reduction in the internal quality factor
$Q_{i}$ and have not been included in the discussion.

[scale=0.35]modelling/Da_a vs Q.png

Figure : Dependence of $\frac{\Delta\alpha}{\alpha}$ as a function of $Q$

###  Sampling frequency

The last effect that is taken into consideration is the interplay between the
sampling rate of the ADCs, that also affects the amplification noise, see
Section 1 in the readout electronics, and the quasi-particle lifetime. The
cardinal theorem of interpolation, as formulated by Nyquist [Nyquist] and
Shannon [Shannon], describes a classical problem in sampling. When the ADC
digitises the analogue signal at its input, it creates a form of uncertainty
in the maximum of the phase signal. According to the aforementioned theorem,
during sampling, the frequencies higher than half of the sampling frequency,
cannot be properly sampled and reconstructed, therefore both the rising and
the trailing edge of the phase signal are subject to uncertainties. It is
impossible to identify with absolute certainty where the peak lays and what
its real value is. The ADC used on ROACH boards is commonly operated with a
$1$ MHz sampling rate, therefore, it is possible to constrain the location of
the maximum within a $1\,\mu$s time frame. In particular, $1\,\mu$s is enough
time to take account of the ring-down time of the resonator. Given this
uncertainty, and knowing the quasi-particle recombination time, it is easy to
evaluate how a $1\,\mu$s interval is reflected into an error on the phase
value. A schematic picture is shown in Figure 1. The recombination time for
this model is $100\,\mu$s; which is a typical value for TiNx resonators.

modelling/sampling.png

Figure : The uncertainty on phase angle depends on the sampling rate and the
quasi-particle recombination time.

$\frac{\Delta\alpha}{\alpha}=1-e^{-\frac{\Delta t}{\tau_{qp}}}$ (11)

Where $\Delta t$ is the time interval between two points being sampled and
$\tau_{qp}$ is the quasi-particle recombination time. For a sampling rate of
$1$ MHz and a recombination time of $100\,\mu$s,

$\frac{\Delta\alpha}{\alpha}=0.00995$ (12)

Obviously, a faster sampling rate coupled to a faster FFT would help, as one
order of magnitude alone would result in an improvement in
$\frac{\Delta\alpha}{\alpha}$ from $\sim 0.01$ to $\sim 0.001$. The other
option, would require the use of a superconducting layer with an intrinsic
recombination time that is longer than $100\,\mu$s. Although this is possible
in theory, depending on the application, it may not be possible for practical
reasons as a longer recombination time would imply a longer dead time of the
detector.

###  Two-Level System Noise

MKIDs exhibit an ”excess” phase noise that is not fully described by the
components previously discussed which is commonly referred as two-level
systems (TLS) noise. The internal quality factor $Q_{i}$ of a resonator is
often limited not by the losses internal to the superconductor, rather by
dissipation through TLS. TLS noise is primarily due to the coupling of the
electric dipole moment of two-level systems with the electric field of the
resonator. In particular, since the interdigitated capacitors tend, by their
own nature, to sustain high electric fields the capacitors are considered the
TLS-noise producing element of an MKID. Even though a semi-empirical model for
TLS noise has been developed, I will not delve into it. Instead it is
interesting to evaluate the contribution to the phase noise solely due to TLS.
When probing MKIDs arrays with tones at -$90$ dBm, a typical phase noise of
$\pm 5^{\circ}$ has been observed. According to the calculations carried out
in the previous section, one might assume the HEMT phase noise of $\pm
3.3^{\circ}$; with the remainder, $\pm 1.7^{\circ}$ can be addressed as TLS
noise. TLS noise is well known to saturate as the power of the probe tones
increases resulting in an increase in the signal to noise ratio. The
description of this behaviour is beyond the scope of this work and will not be
discussed. Instead, Table 1, shows the value of $\frac{\Delta\alpha}{\alpha}$
for different values of the noise component due to TLS.

TLS noise | $\pm 0.5^{\circ}$ | $\pm 0.7^{\circ}$ | $\pm 2.5^{\circ}$ | $\pm 10^{\circ}$
---|---|---|---|---
$\frac{\Delta\alpha}{\alpha}$ | 0.0056 | 0.0189 | 0.0278 | 0.111
Table : Contributions to $\frac{\Delta\alpha}{\alpha}$ for different values of
TLS phase noise, in degrees.

###  Energy resolution

The connection between the energy resolution $R=\frac{\Delta E}{E}$ and the
parameter here discussed $\frac{\Delta\alpha}{\alpha}$ still has to be
established. First of all, a cumulative value for
$\frac{\Delta\alpha}{\alpha}$ can be obtained on the hypothesis that the four
contributions are independent and their quadratures can be added up to an
overall distribution width.

$\frac{\Delta\alpha}{\alpha}=\sqrt{\left(\frac{\Delta\alpha}{\alpha}\right)_{HEMT}^{2}+\left(\frac{\Delta\alpha}{\alpha}\right)_{Currents}^{2}+\left(\frac{\Delta\alpha}{\alpha}\right)_{Readout}^{2}+\left(\frac{\Delta\alpha}{\alpha}\right)_{TLS}^{2}}=0.0429$
(13)

Where the HEMT noise accounts for the most part, and the contribution due to
the current inhomogeneity is about an order of magnitude smaller than that of
the HEMT. The three quantities are represented in Table 1.

$\left(\frac{\Delta\alpha}{\alpha}\right)_{HEMT}$ | $\left(\frac{\Delta\alpha}{\alpha}\right)_{Currents}$ |  | $\left(\frac{\Delta\alpha}{\alpha}\right)_{TLS}$
---|---|---|---
0.037 | 0.004 | 0.0095 | 0.0189

Table : Typical contributions to $\frac{\Delta\alpha}{\alpha}$ from the main
sources of noise: HEMT, Current inhomogeneity, Readout frequency and TLS.

Even though it is intuitive that $\frac{\Delta\alpha}{\alpha}$ and
$R=\frac{\Delta E}{E}$ are related, it is not straightforward to compare
$R_{Fano}=0.01116$ and $\frac{\Delta\alpha}{\alpha}=0.0385$. While deriving an
analytic function that describes the angular displacement in the IQ plane as a
function of the photon energy is possible, and an attempt is described in
Appendix LABEL:appendix:eq-DeltaE, it is perhaps not the optimal choice for
the sake of this work. A more practical approach when working with detectors
is to use a number of calibration wavelengths. By illuminating the detector
with monochromatic light of selected wavelengths ($\lambda_{i}...\lambda_{n})$
and based on the response of the detector at each wavelength, it is possible
to trace a calibration curve. Often times, this curve is simply given by a
linear interpolation between every two adjacent calibration points. In each
segment of the calibration curve the identification described by Equation 14
is correct:

$\frac{\Delta E}{E}=\frac{\Delta\alpha}{\alpha}$ (14)

The only assumption that is necessary to explicit is that at every wavelength
the width of the angular distribution $\Delta\alpha$ caused by a a photon of
wavelength $\lambda$ is either confined within the individual segment defined
by $\lambda_{j}<lambda\lambda_{j+1}$ or on its extrapolation. Based on these
assumptions, the energy resolution corresponding to a value of
$\frac{\Delta\alpha}{\alpha}=0.0429$ is $R_{\alpha}=0.0752$. The two
contributions to the energy resolution can now be added in quadrature:

$R=\sqrt{R_{Fano}^{2}+R_{\alpha}^{2}}=0.0748$ (15)

This means that when illuminated by a red photon, $\lambda=650$ nm, the
wavelength used as an example for the entire length of this chapter, the MKID
taken into consideration will detect its energy to be $1.907\pm 0.142$ eV
corresponding to a resolving power, $P=\frac{E}{\Delta E}=11.52$. It is worth
taking a moment to break down the different contributions that limit the
detector’s performance. In order to do so, it is possible to individually
neglect the elements in Equations 13 and 15 and then list them from the most
to the least detrimental: The noise component that limits the performance of
the detector the most is that induced by the HEMT. If one could, ideally,
reduce it to zero, the resolving power would increase to $P=17.35$. The second
element that contributes is given by hot-phonon losses. If one could reduce
phonon losses to a negligible amount, the resolving power would increase to
$25.75$. Then comes two level systems (TLS) noise, neglecting which, would
bring the resolving power up to $P=47.24$. The readout electronics reduces the
resolving power from $P=76.65$ and finally the contributions given by the Fano
noise and current inhomogeneities, which, when the contribution due to all the
other noise sources have been neglected, alone limit the resolving power to
$P=142.50$. One final comment is on the wavelength dependencies of some the
main contributors to the energy resolution limits. In particular, it is
immediate that the readout sampling frequency has no dependence on the
wavelength at which the resonator is being illuminated. It is intuitive that
the current inhomogeneity is present in the superconducting inductor
independently on the impinging photon’s wavelength. Obviously, a resonator
optimised to produce a $180^{\circ}$ phase shift for $650$ nm will produce a
smaller signal when detecting less energetic photons, and the ratio
$\frac{\Delta\alpha}{\alpha}$ due to it is dependent on the wavelength. Things
get more complicated when considering the effects of the HEMT, and its noise
temperature, on the signal. In this case, $\Delta\alpha$ can be assumed to
remain constant independently on the wavelength on the impinging photon,
whereas the absolute value of $\alpha$ changes with wavelength, effectively
decreasing the signal to noise ratio, and thus reducing the energy resolution
of the detector. Furthermore, this case is not easy to discuss as the
identification in Equation LABEL:eq:daa_1 and Equation 14 are not valid
anymore. Under this circumstance it is not correct to claim that $E_{max}$
describes the value of $E$ if the resonator produces a phase response that is
significantly different from $\pi$ radians. Finally, the Fano contribution
deserves some attention as it may have gone unstated that the value of
$R_{Fano}$ is wavelength dependent, and it is ranges from $\approx 0.075$ for
highly energetic blue photons ($300$ nm) to $\approx 0.017$ for infrared
photons ($1500$ nm).

###  Energy dependence

Following the arguments discussed in the previous section, it is established
that the Energy Resolution is wavelength dependent, therefore the empirical
evaluation of said dependency could be of great interest. Investigating in
detail the effect requires great effort, moreover, interlinking the several
physical processes that contribute to the description of the effect is not
easy. An easier approach is to imagine a resonator of which the response is
above the noise level over the whole range of wavelengths of interest: from
$650$ to $1650$ nm. Ideally, the resonator would be optimised to produce a
$180^{\circ}$ phase signal when struck by a $650$ nm photon and would produce
a signal which is still discernible from the noise at $1650$ nm, let us assume
one order of magnitude above the noise threshold, e.g. a $50^{\circ}$ phase
shift. Ideally, one could imagine the phase shift produced by the resonator to
be directly dependent on the wavelength of the impinging photon. Under these
assumptions, the height of the phase pulse would be described by Equation 16
$\alpha=180-0.13(\lambda-650)$ (16) The different components that make up
$\frac{\Delta\alpha}{\alpha}$ might or might not dependent on the value of
$\alpha$: the contributions due to current crowding and inhomogeneity, TLS
noise, and the HEMT noise depend on the pulse height. Given the same amount of
noise produced by the resonator or the HEMT will have a different impact on
the detector’s performance depending on the height of the phase pulse. On the
contrary, given the self-similarity of the exponential decay, the contribution
due to the readout electronics and its sampling frequency does not depend on
the height of the signal. Under these assumptions, Equation 13 becomes:
$\frac{\Delta\alpha}{\alpha}\left(\lambda\right)=\sqrt{\left(\frac{\Delta\alpha}{\alpha}\left(\lambda\right)\right)_{HEMT}^{2}+\left(\frac{\Delta\alpha}{\alpha}\left(\lambda\right)\right)_{Currents}^{2}+\left(\frac{\Delta\alpha}{\alpha}\left(\lambda\right)\right)_{TLS}^{2}+\left(\frac{\Delta\alpha}{\alpha}\right)_{Readout}^{2}}$
(17) where the $\lambda$ dependency of the parameter $\alpha$ is described by
Equation 16. The curve described by Equation 17 is shown in Figure 1. The
equation that links $\frac{\Delta\alpha}{\alpha}$ to the energy resolution
$\frac{\Delta E}{E}$ is only accurate when the resonators are well optimised
for the specific wavelength, which means, when the phase pulse it produces is
$\approx 180^{\circ}$. Keeping in mind that this identification is only exact
for wavelengths close to $650$ nm, it is possible to extend it over the whole
wavelength range of interest in order to gather qualitative information about
the trending of the curve, which is shown in Figure 1. The scatter plot shows
resolving power values measured by different groups for MKIDs made of
different materials: Guo et al., [Guo], (Green) for their TiN/Ti/TiN MKIDs
measured a value of $P=3.7$ at $1550$ nm, in agreement with what predicted by
the model $P=3.77$, Szypryt et al., [Szypryt_2017], (Red) obtained a resolving
power in the wavelength range $808-1310$ nm that ranges between $8.1-5.8$; and
Mazin, [mazin2020], (Cyan) reports a resolving power of $10$ at $808$ nm for
hafnium (Hf) MKIDs. The experiments seem to be in good agreement with the
trending of the modelled curve.

[scale=0.23]modelling/SNR-wavelength.png

[scale=0.23]modelling/P-energy.png

Figure : Left: Signal to noise ratio $\frac{\Delta\alpha}{\alpha}$ as a
function of the wavelength of the impinging photon. Right: Resolving power
$\frac{E}{\Delta E}$ as a function of the energy of the photons that
illuminate the resonator. The scatter plot represents experimentally measured
points as collected from: Green - [Guo], Cyan - [mazin2020], Red -
[Szypryt_2017].

##  2 Possible improvements

A final section is necessary to point out the results described in the
previous pages and what technological improvements can be imagined might be
relevant in future works. When it comes to the yield of an array, there is
very little improvement that can be foreseen for superconductors such as
$TiN_{x}$ for which it is known that $T_{c}$ homogeneity is very difficult.
Other superconductors can exhibit a better uniformity, $TiN-Ti-TiN$
multilayers have not been tested thoroughly yet. The deposition of uniform and
high quality thin films is well optimised for both $Ti$ and $TiN$. On a more
technological stand point, an enhancement in the maximum frequency at which a
low noise-temperature HEMT can operate, might allow the choice of sacrificing
the close packing of an extra amount of resonators to improve the pixel yield.
A further approach to maximise the yield would be the retuning of the
colliding resonators through post-production processes such as trimming and
DC-bias retuning [DeLucia-multiplexing] (See chapter LABEL:retuning). As
previously discussed, the importance of the energy resolution of the detector
is central and it is crucial for most application that it is as optimised as
possible. The last few pages describe a model created $ad\,hoc$ to investigate
the effects of several parameters on the energy resolution of an MKID: the
noise temperature of the HEMT, which appeared to be the most significant
contribution, the inhomogeneity in current density in the inductor of the
resonator, the sampling frequency of the readout electronics and the Fano
noise that dominates the pair breaking process. The dominant contribution is
provided by the HEMT, and very little can be done to minimise it. Of course a
better HEMT has been in the community’s wish list for several years now, and
hopefully the parametric amplifiers currently under development [Zobrist_2019]
might be crucial for the further development of this technology. That said,
Figure 1 shows that a higher internal quality factor $Q_{i}$ can help improve
the signal to noise ratio. Further improvement, on non dominant factors, is
possible. New geometries can be investigated in order to maximise the current
homogeneity, among these, a design with a large inductor could be interesting
in order to minimise the effects of the current crowding at the right corners,
whereas a tapered inductor could reduce the current density gradient. An in
depth study of the physical properties of new materials, such as the
multilayers and the Granular Aluminium, might produce interesting results in
terms of quasi-particle life time, which would impact the contribution to the
energy resolution given by the readout electronics. The Fano component issue
is currently being addressed by the SRON group in the Netherlands, by
maximising the pair breaking efficiency $\eta$. In the case of semiconductors
[Fano_material] and gases [fano_gas1], the Fano factor has been observed to be
material dependent. It might be worth investigating if there is a significant
change in Fano noise across a broad variety of superconductors.