preparation.tex

% FIXME: Change variable names according to Noise Sources in Bulk CMOS
\chapter{Preparation}%
\label{sec:preparation}
\begin{chapquote}{Lewis Carroll, \textit{Alice in Wonderland}}
``Begin at the beginning,'' the King said gravely, ``and go on till you come to the end: then stop.''
\end{chapquote}
%\section{Grounding and Shielding}
%Add parts from "references\\Grounding and Shielding.pdf"

\section{Laser System}%
\label{sec:prep_laser_system}
%\subsection{Requirements Laser System}
The ARTEMIS experiment is currently in the process of commissioning. Due to the relative abundance of Argon, the ability to create highly charged ion species using an electron beam along with a scientifically relevant transition at $\lambda = \qty{441.25575(17)}{\nm}$ \cite{ar13+_wavelength} with a lifetime of \qty{9.573(6)}{\ms} \cite{ar13+_lifetime} makes \ce{Ar^{13+}} an ideal candidate for this purpose. Figure \ref{fig:bohr_ar13} shows the simplified electronic configuration of the boron-like \ce{Ar^{13+}} investigated.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.28\linewidth}
        \centering
        \import{figures/}{fig_bohr_Ar13.tex}
        \caption{Shell model of boron-like \ce{Ar^{13+}}.}
        \label{fig:bohr_ar13}
    \end{subfigure}
    \begin{subfigure}{0.6\linewidth}
        \centering
        \scalebox{0.75}{%
            \import{figures/}{fig_level_transitions.tex}
        }% scalebox
        \caption{Level diagram of the \ce{Ar^{13+}} transitions accessible at \qty{441}{\nm}. Frequencies with numbers denote optical transitions, letters denote microwave transitions.}
        \label{fig:level_transitions}
    \end{subfigure}
    \caption{Electronic configuration and optically accessible transitions of \ce{Ar^{13+}}.}
    \label{fig:ar13}
\end{figure}

The optical transitions of highly charged \ce{Ar^{13+}} around \qty{441}{\nm} shown in figure \ref{fig:level_transitions} can be used for laser spectroscopy. The transition of interest for commissioning is the transition from the ground state $[(1s)^2 (2s)^2 2p]^2 \, P_{1/2}$ to the first exited state $[(1s)^2 (2s)^2 2p]^2 \, P_{3/2}$ \cite{optical_transitions_ar13}. The Zeeman splitting introduced by the \qty{7}{\tesla} magnet of the ARTEMIS Penning trap results in a frequency splitting of the $^2P_{1/2}$ ground state by $\nu_d = \qty{65}{\GHz}$ and the exited $^2 P_{3/2}$ state of $\nu_a = \qty{130}{\GHz}$. The natural linewidth of these transitions is $\Gamma \approx 2 \pi \times \qty{16.63(1)}{\Hz}$ which is fairly small, but there is substantial Doppler broadening of
\begin{equation}
    \Delta \nu (\lambda = \qty{441}{\nm}, T=\qty{4}{\K}, m=\qty{39.948}{u}) = \frac{2}{\lambda}\sqrt{2 \ln 2 \frac{k_B T}{m}} \approx 2 \pi \times \qty{150}{\MHz} \, , \label{eqn:doppler_broadening}
\end{equation}
seen in the trap which is kept at a temperature of \qty{4}{\K}.

The simplified laser setup shown in figure \ref{fig:laser_system_old} as a schematic was characterised by \citeauthor{thesis_alex} \cite{thesis_alex} and its transfer accuracy of the the wavelength was calculated to be \qty{1.2}{\MHz}. Considering the additional long-term drift of the system led to an upper limit of the absolute wavelength uncertainty of \qty{2.2}{\MHz}. This can be considered more than adequate given the Doppler broadening of \qty{150}{\MHz} calculated in equation \ref{eqn:doppler_broadening}. While being sufficiently accurate, the system has the significant drawback of being fairly complicated to manage if not maintained on a daily basis. Its performance and uncertainty relies on the exact knowledge of the tellurium spectrum surrounding both wavelengths of the lasers.
\begin{figure}[ht]
    \centering
    %\scalebox{1}{%
        \import{figures/}{laser_system_old.tex}
    %} % scalebox
    \caption{Simplified setup of the laser system in use at ARTEMIS prior to this work. Blue lines are laser beams, black lines are electronic signals delivering feedback.}
    \label{fig:laser_system_old}
\end{figure}

To operate the system, the master laser must first be locked to a tellurium reference transition which has to be searched manually using the tellurium spectrum charted by \cite{thesis_alex}. While not complicated, it requires a trained user nonetheless. The transfer cavity is then locked to the master laser, a straight forward process. The spectroscopy laser now has to be locked to the correct fringe of the transfer cavity. This is done by observing the tellurium background again to adjust the diode laser mode to a frequency close to the desired cavity fringe. This is by no means an automated process and requires a competent operator. Although this locking procedure requires training it does work reliably, but there are more issues that are not readily apparent.

One potential problem lies with the master laser. The whole calibration is geared towards a reference transition at \qty{452.756}{\nm} catalogued by \citeauthor{te_reference_lines} \cite{te_reference_lines}. Replacing the master laser in case of a failure is challenging. Blue laser diodes are less flexible when tuning in comparison to their NIR cousins as discussed in \cite{thesis_baus,thesis_alex}. The diode used in the current laser was handpicked for this wavelength. This limits the availability of replacement parts and increases their replacement value. Using another diode and wavelength along with a different tellurium transition would require creating a new tellurium map which is a laborious process.

Another challenge is created by the locking scheme which becomes more complicated when introducing a third laser to complete the closed transition for the laser-microwave double-resonance spectroscopy. The setup in its present condition was already prepared for the second spectroscopy laser which must also be locked to the transfer cavity. Currently both lasers are sent into the transfer cavity with perpendicular polarisations to separate the beams using a polarizing beam splitter. Since the reference laser is running at \qty{453}{\nm} while the spectroscopy laser is using \qty{441}{\nm}, it is also possible to separate those two via a dichroic mirror. While possible, this scheme requires the close overlap of three beams along with their reflection as required for saturation. Such a setup strongly couples the three beam paths and further complicates future adjustments.

The final issue regarding the transfer cavity concerns the required high voltage for the piezoelectric actuator to adjust its length. This piezo requires up to \qty{1.25}{\kV} to reach the necessary translation required for scanning the tellurium spectrum. Not only does this pose a risk to an untrained operator, but it was also discovered during the commissioning of the current system that high frequency noise from the voltage supply is radiated and can enter the experiment. These issues were minimized by testing several power supplies, choosing the lowest noise device and keeping the supply well separated from the rest of the experiment increasing space required for the whole system.

During the commissioning and testing of the laser system it was also found that the laser drivers had issues with blue laser diodes. Those drivers used a modified in-house current source normally used for NIR diodes based on the design of \citeauthor{libbrecht_hall} \cite{libbrecht_hall}. With these drivers an increasing instability was observed when adjusting the current up to the operating point. The origin of this problem lies in the larger operating voltage of the blue laser diode of up to \qty{7}{\V} compared to the more moderate \qty{2.5}{\V} of NIR diodes. The details are discussed in section \ref{sec:compliance_voltage}. Other commercial drivers tested either had the same problem or were far noisier and therefore harder to work with, due to the modification made by the manufacturers to increase the compliance voltage.

To conclude, the commissioning of the laser system has brought up two issues with the current system that impact its availability and performance. The transfer cavity can be considered the Achilles' heel of the system and replacing it with a more flexible alternative like a high performance wavemeter can greatly improve the usability and flexibility of the whole laser system as it breaks up the tight dependency on the tellurium reference laser. Additionally this opens up alternative wavelengths for the spectroscopy of other highly charges ion species. The second issue was found with the laser driver and the apparent lack of commercial solutions sparked the development of a novel laser driver for the next generation of laser diodes to gain direct access to more wavelengths. This is discussed in the next section.

\clearpage
\section{Laser Current Driver}%
\label{sec:laser_current_driver}
% Include Emission wavelength dependence of characteristic temperature of InGaN laser diodes
% Check Diode Laser Characteristics
% I-lamda in Wavelength Dependence of InGaN Laser Diode Characteristics
% also Determination of piezoelectric fields in strained GaInN quantum wells using the quantum-confined Stark effect
Laser diodes are current driven devices, because
\begin{equation*}
    P_{out} \propto I
\end{equation*}
and the diode current $I$ approximately follows the Shockley equation \cite{shockley_diode}
\begin{equation}
    I = I_0 \left( e^{\frac{qV_d}{k_B T}} - 1\right) \, . \label{eqn:shockley}
\end{equation}

where $k_B$ is Boltzmann constant, $T$ the temperature, $q$ the electron charge and $V_d$ the diode voltage. The exponential dependence of the current on the supply voltage calls for a current source to drive a laser diode safely without risking thermal damage due to excessive injection currents.

The primary function of a laser driver is therefore to provide a stable, but user adjustable, current. This current can typically be modulated at frequencies up to several \unit{\MHz} to shape the frequency and amplitude of the laser output light. Additional features, like current and voltage limits, aid in protecting the expensive laser diodes and it is not uncommon to have additional safeguards inside the laser head that are under control of the laser driver like a shorting relay to ensure the laser diode is shorted when the driver is disconnected or disabled.

The focus of this work lies on two types of laser diodes, indium gallium nitride (InGaN) and aluminium gallium arsenide (AlGaAs), but is not limited to those two types. The former material is, for example, used for blue laser diodes at around \qty{450}{\nm}, discussed in the previous section, and laser diodes up to green wavelengths, the latter is used for near-infrared (NIR) laser diodes such as \qty{780}{\nm} laser diodes. Both wavelengths are used for experiments in the Atoms-Photons-Quanta group (in future referred to as \textit{the group}). The former type is used in the ARTEMIS experiment for the spectroscopy of highly charged ions, the latter is extensively used to manipulate and control the rubidium atoms used by the quantum computing experiments. This section deals with the design challenges of such a device used for high precision laser spectroscopy. First the design requirements are established and, from those conditions, technical specifications are developed.

The design requirements are split into three parts which need to be discussed: The environment includes effects like temperature, humidity and time. That section mostly focuses on the ambient temperature though because its effects are the most pronounced. The current source electrical requirements, like drift, noise, output impedance, and modulation bandwidth are discussed as these have a profound impact on the intended application in experiments. Finally, the user interface including the external communication interfaces are defined.

\subsection{Design Goals: Ambient Environment}%
\label{sec:design_goal_environment}
The lasers and its accompanying driver is to be mostly used in a clean laboratory environment. In this particular use case the air is filtered using H14 HEPA filters, but less rigorously controlled environments must be considered as well, because not all fields of application are in optical labs. A mostly dust free industrial environment is considered acceptable as well. Typical lab temperatures are in the range of \qtyrange[range-phrase=\textup{~to~}]{20}{30}{\celsius}. This temperature range was also encountered in the labs discussed in this work before improvements were implemented as part of this work. The upper end of the range must be considered when operating the devices inside a rack where the temperatures are even higher and the device should therefore be tested for its upper limit. A temperature of \qty{35}{\celsius} is a typical value measured inside the racks used in the lab. Humidity is only controlled with dehumidifiers, limiting only the upper bound, resulting in a range of \qtyrange[range-phrase=\textup{~to~}]{15}{60}{\percent rH}.

Figure \ref{fig:lab_temperature_start_of_project} shows a typical one day span of the lab temperature as it was found at the start of this project, plotted using Matplotlib \cite{matplotlib}.
\begin{figure}[ht]
    \centering
    \input{images/temperature_011_2016.pgf}% data/plot_generic.py
    \caption{Temperature in lab 011 of the APQ group on 2016-11-26. Recorded by the LabKraken monitor. See section \ref{sec:res_labkraken} for details.}
    \label{fig:lab_temperature_start_of_project}
\end{figure}

As it can be seen there are strong oscillations of the temperature around the setpoint of \qty{21}{\celsius} as a result of the on–off air conditioning temperature controller. The commercial controller initially installed was using an IMI Heimeier \device{EMO T} thermoelectric actuator \cite{datasheet_heimeier_emo_t}, which is a two-step valve. Although this solution was later replaced by a custom design described in section \ref{sec:lab_temp_control}, these type of controllers are found in many other labs and temperature swings of \qty{2}{\kelvin} must therefore be expected.

These environmental parameters can now be used to estimate the design requirements for the laser driver. In comparison to the other laser system used in this group, the \qty{450}{\nm} system \cite{thesis_baus} required for the spectroscopy of highly charged ions \cite{thesis_alex} at GSI is the more demanding system. This system was found to be more susceptible to changes of the drive current since the wavelength selective filter element was far broader in comparison to a \qty{780}{\nm} system \cite{two_filter_paper}. This laser is stable over regions of tens of \unit{\uA} and requires a maximum drive current of \qty{145}{\mA} \cite{datasheet_osram_pl450b}.

From these considerations, the requirements for the driver can be inferred. It should be able to supply at least \qty{150}{\mA} and stay well within \qty{10}{\uA} over the whole environmental range. Given a worst-case scenario a tolerance of $3\sigma$ (\qty{99.7}{\percent}) must be met \cite{worst_case_design}.

The environmental parameters that mostly affect current sources are temperature and humidity. Air pressure is typically a matter of concern for high voltage systems \cite{IPC-2221B} and secondary in consideration for this design as it is a low voltage system (\qty{<= 48}{\V}). Air pressure effects are also the most expensive to test for, as a pressure chamber is required. Humidity affects electronics both directly though corrosion and also indirectly because the epoxy resin used in the FR-4 PCBs and component moulding is hygroscopic and the absorbed humidity leads to swelling and mechanical stress. This effect is very slow at ambient temperature and can easily take days to show \cite{epoxy_humidity}. This parameter is therefore handled via the long-term stability and not specified separately.

Given environmental conditions, the relative coefficients can be calculated. This estimation assumes a minimum setpoint resolution of 2 steps within the mode-hop-free region of the laser and calculates the \qty{99.7}{\percent} confidence interval. The steps are given in table \ref{tab:dgdrive_tempco}:
\begin{table}[hb]
    \centering
    \begin{tabular}{llr}
        \toprule
        Property& Value& Result \\
        \midrule
        Stable range & \qty{10}{\uA}& \qty{10}{\uA}\\
        2 steps of resolution  & $\div 2$& \qty{5}{\uA} \\
        $1 \sigma$  & $\div 2$& \qty{2.5}{\uA} \\
        Maximum output& \qty{150}{\mA}& \qty{17}{\uA \per \A}\\
        Temperature range& \qty{5}{\K}& \qty{3}{\uA \per \A \per \K}\\
        Worst case ($3 \sigma$)& $\div 3$& \qty{1}{\uA \per \A \per \K}\\
        \bottomrule
    \end{tabular}
    \caption{Estimated requirement for the temperature coefficient of the laser driver.}
    \label{tab:dgdrive_tempco}
\end{table}

While the requirements look moderate at first sight, tuning a quick estimation shown in table \ref{tab:dgdrive_tempco} leads to a temperature coefficient of \qty[per-mode = symbol]{1}{\uA \per \A \per \K} or even tighter when using a higher output driver -- a rather formidable specification for a current source.

Regarding the long-term stability, a \qty{30}{\day} number can be estimated. One may be inclined to call for a drift which is smaller than the stable range, but this would be short-sighted, as there are other factors to consider. The laser including the external resonator has its own figure of merit regarding the spectral drift rate. \citeauthor{ecdl_stability} \cite{ecdl_stability} reported a drift of \qty{2.9}{\MHz \per \hour}, which was attributed either to the external resonator itself, the piezo or the collimation lens. It is most likely that this drift was caused by mechanical changes of the external resonator as it defines the output mode of the laser. The mechanical drift limits the required stability of the current source considerably, as a typical frequency change of the internal resonator with the current of \qty[per-mode=symbol]{3}{\MHz \per \micro \A} \cite{diodelaser_modulation} can be assumed. The (linear) ageing drift of the external resonator over \qty{30}{\day} is equivalent to a \qty{720}{\uA} drift over the same period. For the electronics, the drift is assumed to follow an Arrhenius-like equation resulting from stress, caused during manufacturing. This may eventually change to a slow linear drift after several months of relaxation. The coefficient can either be a positive or negative and leads to
\begin{table}[hb]
    \centering
    \begin{tabular}{llr}
        \toprule
        Property& Value& Result \\
        \midrule
        Ageing drift limit & \qty{720}{\uA}& \qty{720}{\uA}\\
        $1 \sigma$  & $\div 2$& \qty{360}{\uA} \\
        Maximum output& \qty{500}{\mA}& \qty{720}{\uA \per \A}\\
        Worst case ($3 \sigma$)& $\div 3$& \qty{240}{\uA \per \A}\\
        \bottomrule
    \end{tabular}
    \caption{Estimated requirement for the long-term stability over \qty{30}{\day} of the laser driver.}
    \label{tab:dgdrive_stability}
\end{table}

Based on these numbers, it is straightforward to see that the long-term stability of a laser driver is less important than the short-term temperature coefficient since the limiting factor is the mechanical construction of the laser. This necessitates an atomic reference for long-term stability and to compensate for acoustic resonances of the external resonator. Regarding the choice of suitable devices, the tight specification of the temperature coefficient most likely leads to a choice of components that will pass these long-term criteria as well, alleviating a bit the burden of proof as long-term drift specifications are hard to come by since they need a lot of time to validate and cannot be extrapolated from high temperature burn-in tests \cite{voltage_reference_drift}.

%While \qty{800}{\uA \per \A} over a \qty{30}{\day} period may seem large at first, it is actually very hard to accurately produce a current. To put this number in perspective, a commercial high-end current source like the Keithley \device{2600B} is specified for the \qty{100}{\mA} range at about \qty{170}{\uA \per \A} for a \qty{30}{\day} period when calculated from the 1-year specification \cite{datasheet_keithley2600}, again assuming an Arrhenius-like equation as the basis.

All of this leads to the following design specifications regarding the stability of the current driver:
\begin{center}
    \begin{specifications}[label={lst:dgDrive_specs_environment}]{Current source, environmental}
    \begin{itemize}
        \item Temperature range \qtyrange[text-series-to-math, reset-text-series = false, reset-math-version = false, range-phrase=\textup{~to~}]{20}{35}{\celsius}
        \item \textbf{Temperature coefficient \qty[text-series-to-math, reset-text-series = false, reset-math-version = false]{<= 1}{\uA \per \A \per \K}}
        \item Humidity (non-condensing) \qty{<= 75}{\percent rH}
        \item Humidity coefficient not specified, but included in the long-term drift
        \item Maximum altitude not specified
        \item Long-term drift over \qty{30}{\day} \qty{<= 240}{\uA \per \A}
    \end{itemize}
    \end{specifications}
\end{center}

%A basic laser current driver design that has some of the  can be found in the work of \citeauthor{libbrecht_hall} \cite{libbrecht_hall}. While this design contains all the basic features, like a current source, a modulation input and a voltage limit, there are several shortcomings that have emerged over the years with new generations of laser didoes. The laser driver used by legacy applications in this group is based on the aforementioned paper and has been successfully employed in several projects over the years, but several limitations have come up in recent years. In order to derive the design requirements of a new generation of laser drivers the important design elements need to be identified first. The essential design elements are the bulk current source, the modulation current source, the reference element and output programming. The next 4 sections will deal with each element and outline the design goals that were identified while employing the legacy generation of diode drivers in several experiments.

\subsection{Design Goals: Current Source}
% https://www.laserdiodecontrol.com/laser-diode-parameter-overview
% Diode Lasers and Photonic Integrated Circuits (Characteristic temperature)
% Near Threshold Operation of Semiconductor Lasers and Resonant-Type Laser Amplifiers
The change in output current caused by load impedance should be an order of magnitude less than the drift specification to ensure a negligible effect compared to the drift over time. The load resistance presented by the laser diodes most commonly used in our experiments ranges from \qty{50}{\ohm} \cite{datasheet_osram_pl450b} to \qty{30}{\ohm} \cite{datasheet_adl_785} and \qtyrange{10}{15}{\ohm} for \qty{780}{\nm} laser diode \cite{datasheet_sharp_780nm,datasheet_thorlabs_780nm}. The output impedance requirement can therefore be estimated as
\begin{align}
    \frac{R_{load}}{R_{out}} &= \frac{I_{set}}{I_{out}} - 1 \leq \qty[per-mode = symbol]{6.7}{\uA \per \A} \nonumber\\
    R_{out} &\geq \frac{\qty{50}{\ohm}}{\qty[per-mode = symbol]{6.7}{\uA \per \A}} = \qty{7.5}{\mega \ohm}
\end{align}

An output impedance of more than \qty{7.5}{\mega \ohm} for slowly changing loads is a tough requirement, depending on the type of current source, which requires carefully selected components. A high output impedance is, for example, of importance to suppress radiated noise coming from external sources. Especially low frequency components from the mains supply can magnetically couple into the cables, because they are long enough. This noise can be substantial and a high output impedance at low frequencies is therefore important. Other applications will be discussed throughout this work. While a subpar output impedance is more of a limiting factor, the compliance voltage discussed next is a key requirement.

The compliance voltage is the maximum voltage the current source can apply to the load and is another non-ideal component of a real current source. The required voltage strongly depends on the type of laser diode used. The near-infrared laser diodes discussed above have an operating voltage of \qtyrange{1.5}{3}{\V}, while the Osram \device{PL 450B} blue laser diode is specified for \qtyrange{5.5}{7}{\V}. The \qty{7}{\V} required by the Osram laser diode is fairly high for a Fabry–Pérot laser diode and has proven difficult in the past \cite{thesis_baus} as most laser current drivers available are designed for the much lower forward voltage of the near infrared laser diodes. Even higher voltages of around \qtyrange{12}{15}{\V} are required for quantum cascade lasers, but these are currently neither used nor is their use planned in any experiment in the group.

The maximum output current of the laser driver currently required for laser diodes used in the group is \qty{250}{\mA} for the Thorlabs \device{L785H1} \cite{datasheet_thorlabs_780nm}. Therefore a maximum output current of \qty{300}{\mA} is considered sufficient.

The current noise of the laser driver can be estimated from the laser linewidth sought after as the laser frequency is sensitive to the injection current. At low frequencies, about \qty[per-mode=symbol]{-3}{\MHz \per \micro \A} can be attributed to the thermal expansion of the internal resonator of the diode due to resistive heating \cite{diodelaser_modulation}. Above \qty{1}{\MHz} this effect starts declining and exposes the change of the refractive index due to the presence of charge carriers. This high frequency effect is an order of magnitude weaker. Since the frequency sensitivity to current variations of the laser diode drops with higher frequencies, the most important range is from DC to \qty{100}{\kHz}.

To estimate the linewidth requirement, it is important to look at the experimental setup. While the spectroscopy of \ce{Ar^13+} at \qty{4}{\K} is limited to around \qty{150}{\MHz}  as shown on page \pageref{eqn:doppler_broadening}, the quantum computing experiments in this group have more stringent needs. It was shown in \cite{ecdl_stability, ecdl_silicone_housing,ecdl_linewidth_scholten} that with reasonable expense a passive linewidth of less than \qty{100}{\kHz} can be achieved. Using the relationship of the frequency sensitivity to a current modulation of laser diodes, \qty{100}{\kHz} translates to a current noise of \qty{30}{\nA_{rms}} from \qty{1}{\Hz} to \qty{100}{\kHz}. The lower \qty{1}{\Hz} limit is chosen fairly arbitrary, but the presence of $\frac 1 f$-noise inhibits a definition down to DC. There should be negligible amounts noise below \qty{1}{\Hz} compared to the upper \qty{100}{\kHz} though.

The final aspect of the current source that needs to be specified, is the bandwidth of the current steering input. The bandwidth in these terms defines a reasonably flat (\qty{\leq 3}{\dB}) response. As it was discussed above, beyond a frequency of \qty{1}{\MHz}, the frequency sensitivity of the laser diode to current modulation drops by an order of magnitude, altering the transfer function and introducing new challenges for control loops. Therefore a minimum bandwidth of \qty{1}{\MHz} is considered sufficient.

Above \qty{1}{\MHz} it is recommended to either use more dedicated solutions like the direct modulation at the laser head presented in \cite{current_mod_paper} or switch to acousto-optic modulators (AOMs) or electro-optic modulators (EOMs).

This leads to the following requirements regarding the current source of the laser driver:
\begin{center}
    \begin{specifications}[label={lst:dgDrive_specs_electrical}]{Current source, electrical}
    \begin{itemize}
        \item Maximum output current \qty{300}{\mA}, optionally \qty{500}{\mA}
        \item \textbf{Compliance voltage \qty[text-series-to-math, reset-text-series = false, reset-math-version = false]{\geq 8}{\V}}
        \item Output impedance \qty{\geq 7.5}{\mega\ohm} at low frequencies (close to DC)
        \item \textbf{Current noise \qty[text-series-to-math, reset-text-series = false, reset-math-version = false]{\leq 30}{\nA_{rms}} from DC to \qty[text-series-to-math, reset-text-series = false, reset-math-version = false]{100}{\kHz}}
        \item \qty{3}{\dB}-bandwidth of the modulation source \qty{\geq 1}{\MHz}
    \end{itemize}
    \end{specifications}
\end{center}

\subsection{Design Goals: User Interface and Form Factor}
The user interface must allow repeatability and reproducibility of the outputs. The reason is that the laser system is intended to be portable to be moved from the university where it is performance tested to the GSI facility. Within the labs, systems are usually moved from test stands to the actual experiment as well. Requiring as little setup efforts as possible is a big advantage.

The interface must both be accessible both locally and remotely to allow simple adjustment of the parameters while on the bench and also from within the experimental control software. The local controls must be directly accessible to humans without tools to give a better user experience.

The remote user interface is strictly required because the Penning trap and the laser system are spatially separated with the laser system being located in a special laser lab for environmental as well as safety reasons. This separation is about \qty{30}{\meter}. Ideally this remote interface is computer controlled to give full access to all features of the laser system. USB or Ethernet is preferred as this does not require extra hardware in the lab.

Regarding the application programming interface (API), support for both Python and optionally LabVIEW is favoured, as most of the group has switched from LabVIEW to labscript suite \cite{labscript_2013} on Python to run the experiments.

The form factor should allow integration into standard 19-inch racks to allow simple transportation from the experiment location at GSI to the university for testing and calibration.

\begin{center}
    \begin{specifications}[label={lst:dgDrive_specs_api}]{Current source, user interface}
    \begin{itemize}
        \item Local control via the front end without tools
        \item Remote access via a digital interface
        \item Software API supporting \textbf{Python} and optionally LabVIEW
    \end{itemize}
    \end{specifications}
\end{center}

\clearpage
\section{Laser Temperature Controller}%
\label{sec:laser_temperatrure_controller}
The external cavity diode laser (ECDL) design employed at GSI and in this group, based on \cite{ecdl_paris}, consists of two parts: The laser diode, mounted in an aluminium frame containing a collimator, which is mounted in an external resonator also made of aluminium. The aluminum used for the external resonator is AlZn4.5Mg1, also called alloy 7200 \cite{datasheet_laser_alu}. It has a moderate thermal coefficient of expansion of \qty{23.1}{\micro\meter \per \m \per \K}, which is one order of magnitude larger than that of Invar, but is significantly easier to machine.

In order to derive the required stability criteria the laser diode and the external resonator must both be considered. The influence of external parameters on the laser wavelength were discussed in the work of \citeauthor{thesis_tilman} \cite{thesis_tilman}. The temperature sensitivity of a typical near-infrared laser diode at \qty{780}{\nm} along with the external resonator used in this group were calculated to be
\begin{align*}
    K_{T,diode} &\approx \qty{-3}{\GHz \per \K}\\
    K_{T,resonator} &\approx \qty{-9}{\GHz \per \K}\,.
\end{align*}

From these number it is clear that the resonator marks the lower bound. Going to a blue \qty{441}{\nm} laser this criterion is even more critical, because $K_{T,resonator}$ is proportional to the laser frequency and the frequency almost doubles, this leads to a sensitivity of the resonator on the order of
\begin{equation}
    K_{T,resonator} \approx \qty{-16}{\GHz \per \K}\,.
\end{equation}

This implies that in order to match the stability of the laser current driver, the temperate stability should be far better than \qty{1}{\milli \K}. A temperature stability of better than \qty{100}{\micro \K} has been demonstrated before \cite{tempcontroller_10uK,tempcontroller_10uK_jw,tempcontroller_15uK,tempcontroller_30uK,tempcontroller_40uK,tempcontroller_50uK,tempcontroller_65uK}, but all of these solution have in common that they use either multiple layers of shielding and control or elaborate baths into which the subject is submerged. The controller itself is then typically placed inside the controlled environment to shield it from external effects. This type of setup is not feasible in this situation as it would require a considerable redesign of the laser resonator. The laser resonators in use in this group \cite{thesis_tilman} have been set up over the course of several years and there are dozens of them in use. This existing design must therefore be taken into consideration as well. The resonator in its current state does does not have an airtight seal. The sensitivity of the laser frequency to the barometric pressure can be estimated using the formula developed by \citeauthor{ciddor} \cite{ciddor} to be
\begin{equation}
    K_{baro} =  \qty{-75}{\MHz \per \hecto \Pa}\,.
\end{equation}

This leads to a frequency drift of several \qty{100}{\MHz} due to a pressure drift of around \qty{\pm 10}{\hecto \Pa} observed in the lab over a typical day. A long-term drift of \qty{55}{\hecto \Pa} over the year 2022 was also recorded by the monitoring software LabKraken. On shorter time scales the air pressure varies on the order of tens of \unit{\Pa}. This must be matched by the temperature controller. It is therefore sufficient to call for a stability of \qty{<1}{\milli \K} when using an unsealed resonator. To guarantee such a stability, the resolution of the driver should be at least \qty{200}{\micro \K}, preferably \qty{100}{\micro \K}.

The type of temperature transducer used in the laser design is a \qty{10}{\kilo \ohm} thermistor, so the design must work with this type of sensor, while the support of other sensors like a PT100 is optional.

Finally, a problem often encountered with analog proportional–integral–derivative (PID) controllers when temperature controlling large, well isolated bodies is the long time scales involved. One example found in the lab are high finesse cavities mounted in vacuum enclosures. These extremely stable cavities are extensively used to reduce the linewidth of lasers to a few \unit{\Hz}. The time scales involved necessitate very long integration times $T_i$ or a rather small gain of the integral term $k_i$ of the PID controller. See section \ref{sec:pid_controller_basics} for details on PID controllers and the terminology. An illustration of the problem encountered with an analog controller can be seen in figure \ref{fig:stability_cavity}. It shows the temperature of a Stable Laser Systems \device{VH 6020} cavity housing used for a high finesse Fabry-Pérot cavity which has a time constant of \qtyrange[range-units = single, range-phrase={~to~}]{4}{7}{\hour} \cite{datasheet_vh6020}.
\begin{figure}[ht]
    \centering
    \input{images/stability_cavity.pgf}% data/plot_generic.py
    \caption{Temperature of a Stable Laser Systems \device{VH 6020} controlled by a Team Wavelength \device{HTC1500} temperature controller.}
    \label{fig:stability_cavity}
\end{figure}

The Team Wavelength \device{HTC1500} used in this example is an analog PID controller configured for the maximum specified integration time of \qty{10}{\s} using a \qty{10}{\uF} capacitor. As can be seen by the oscillatory behaviour, this time constant is far too short. Longer time scales can easily be reached by digital controllers which allow very long integration times limited only by the numerical resolution. A digital system can extend the scope of application from lasers to many other systems like those cavities to improve their stability. In addition, a digital system gives more control over the PID tuning parameters, also ensuring repeatability which greatly simplifies setting up new laser systems because a common set of PID parameters can be used a starting point before tuning the controller. Another benefit is the possibility to implement a modified algorithm and additional filters as detailed in section \ref{sec:pid_controller_basics}. This versatility to quickly adapt the programming again increases the number of applications. Integrating autotuning algorithms to help the user find a usable set of parameters reduces setup times of new systems. For all of those reasons, the new controller should be based on a digital design.

The final aspect to be considered is the output driver of the controller. The laser design used in this group uses two Peltier elements to cool both the resonator and the laser diode independently. The driver must therefore integrate two channels. While the biggest TECs currently in use are Laird \device{CP14,127,06,L1,W4.5} which can draw up to \qty{6}{\A} at \qty{15.4}{\V} \cite{datasheet_tec} their optimal coefficient of power is between \qtyrange[range-units = single]{1}{2}{\A} and \qtyrange[range-units = single]{5}{8}{\V}. Having a driver that can output \qty{4}{\A} at \qty{12}{\V} is considered more than sufficient, even for larger TECs and future projects.

Commercial temperature controllers specifying a stability of better than \qty{1}{\milli \kelvin} are hard to come by, especially with multiple channels. Two units were tested for this laser setup. The Vescent \device{SLICE-QTC} and an ILX Lightwave \device{LDT-5948}. The latter is specified for a stability of \qty{5}{\milli \kelvin}, but their application note claims a better performance \cite{appnote_ilx_millikelvin}.

The requirements for the temperature controller can be summarised as:
\begin{center}
    \begin{specifications}[label={lst:dgTemp_requirements}]{Temperature controller}
    \begin{itemize}
        \item Stability: \qty{<1}{\milli \K}
        \item Resolution: \qty{<200}{\micro \K}, \qty{<100}{\micro \K} preferred
        \item Temperature sensor: \qty{10}{\kilo \ohm} thermistor
        \item Two or more channels
        \item Output power: \qty{4}{\A} at \qty{12}{\V}
        \item Digital interface to work with long time scales and reproducible PID parameters
    \end{itemize}
    \end{specifications}
\end{center}

\clearpage
\section{LabKraken}%
\label{sec:prep_labkraken}
\subsection{Design Goals}
LabKraken is designed to be an asynchronous, resilient data acquisition suite that scales to thousands of sensors and across different networks to serve the need for monitoring and automation required for large scale experiments spanning multiple sites. It is written in Python and supports many sensors and instruments found in a scientific environment. Such sensors include Standard Commands for Programmable Instruments (SCPI) capable devices accessible via Ethernet or GPIB or sensors using a serial protocol. Many other Ethernet capable devices are also supported via a simple driver interface.

\subsection{Software Architecture}
To meet the increasing demand for high quality data, LabKraken needs to scale to thousands of sensors which must to be served concurrently. This problem is commonly referred to as the C10K problem as dubbed by \citeauthor{10kProblem} back in 1999 \cite{10kProblem} and refers to handling \num{10000} concurrent connections via network sockets. While today millions of concurrent connections can be handled by servers, handling \num{10000} can still be challenging, especially if the data sources are heterogeneous as is typical for sensor networks of diverse sensors from different manufacturers.

In order to meet the design goals, an asynchronous architecture was chosen and several different approaches were implemented over time. All in all, four complete rewrites of the software were made to arrive at the architecture introduced here. The reason for the rewrites is mostly historical and can be explained by the development of the Python programming language which was used to write the code. The first version was written using Python 2.6 and exclusively supported sensors made by Tinkerforge. With the release of Python 3.5 which supported a new syntax for asynchronous coroutines, the software was rewritten from scratch to support this new syntax, because it made the code a lot more verbose and easier to follow. When Python 3.7 was released asynchronous generator expressions where mature enough to be used in productions and the program was again rewritten to use the new syntax. In 2021 a new approach was taken and the program was once more rewritten with a functional programming style. Some of those approaches will be discussed in the next sections to highlight limits of the programming style used and the improvements made to overcome them. This development underlines important steps in the progress of asynchronous programming made by the Python language in recent years that can be applied to many other problems, for example process control. Specifically so since Python is a very popular language among scientists and used in many experiments. Each of the following sections discusses the same program, but written in different programming styles to show the differences. Especially the last example presenting a function programming style is interesting for experimental control as it gives a clean representation of the data flow from the producer to the consumer \cite{concurrent_programming}.

The example program that will be discussed does the following job. It opens a network connection to a remote (Tinkerforge) sensor platform, then queries the other side for its sensors. When the sensors are returned it looks for a specific sensor, then starts reading data from that sensor to finally print it. The example itself is designed around the Tinkerforge sensors to present a working example instead of the typically used pseudocode. It does represent a common program flow in a sensor application though and the concept is not limited to the Tinkerforge programming API.

\subsubsection{Threaded Design}
The first version of LabKraken used a threaded design approach, because the original libraries of the Tinkerforge sensors are built around threads. Most threaded programs make extensive use of callbacks. These are functions that are passed from the main thread to the worker, typically on creation, and are then called by the worker to inform the main thread of its activity. The downside is that main thread has no knowledge about the caller, the callback might have even been passed on by the worker to another thread.

\lstinputlisting[firstline=1, firstnumber=1, frame=single, breakindent=.5\textwidth, frame=single, breaklines=true, numbers=left, xleftmargin=2em, numberstyle=\tiny, style=mypython]{source/lab_kraken_threads.py}

The program starts at line \num{24} by making a connection to the host, then the first callback functions are registered with the connection object. These callbacks allow the connection thread to signal the program when the connection has been established (\textit{cb\_connected}) and when new sensors are found (\textit{cb\_enumerate}). The main program is now finished and it waits until terminated by the user. All the work is done inside the thread and the program flow unfortunately looses itself in the callbacks which get called by the connection object and their order can only be guessed from the documentation as the main program has no control over it. As the program continues it first enters the \textit{cb\_connected} callback where it will query the host for its sensors in line \num{9}. The answer will be returned through the \textit{cb\_enumerate} callback. This function filters the sensor id for a known sensor and then attaches another callback for the sensor to return data. It then configures the sensor. This program flow is typical for a callback driven design and the reader may imagine how more complex tasks are implemented. As the program grows, more and more layers of callbacks will be added and in the end, the code will be impossible to read without intimate knowledge. The effort of maintaining the callback driven code resulted in the decision to redesign the program when moving to Python 3.

To untangle this problem, Python 3.7 introduced so-called generators. This is a type of expression that will produce values from an iterator. An iterator is an (infinite) ordered series of values or events which can be processed by requesting the next value until the series is exhausted, if it is finite. The main advantage is that the logic of the program stays within the main part and only the data gathering is done outside of this scope. A generator based program is shown next.

\subsubsection{Generator Design}
In addition to a different coding style the code base is moved from a multithreaded to a asynchronous approach using Python asyncio. The difference is that asyncio uses a single thread as opposed to multithreaded code. Multiple threads can run concurrently on multiple processor cores, so different cores can process data at the same time. Asynchronous programs must pause the execution of code paths because they run within a single thread on a single core. This type of programming works best when the tasks are not computationally intensive but input/output bound by external peripherals like a network. While the processor is waiting for the slower network it can work on other tasks. The advantage is that access to shared resources is greatly simplified as these resources will never be accessed at the same time. The code is shown below.
\lstinputlisting[firstline=4, firstnumber=1, frame=single, breakindent=.5\textwidth, frame=single, breaklines=true, numbers=left, xleftmargin=2em, numberstyle=\tiny, style=mypython]{source/lab_kraken_async.py}

The first impression that can be gather from the new design is that the code has become more concise. To understand it, a few Python language keywords used must be introduced. In order to yield control to the next task, the keyword \textit{await} is used which is put in front of a function call. This will pause the current execution and wait until the function has returned with a result. Another important language feature used is a so-called context. A context is created using the \textit{async with} command and it makes sure that after leaving the context certain commands are executed. This can be used to clean up after the the creation and use of certain objects like the Ethernet connection. The connection context will make sure that the Ethernet connection will be properly terminated no matter whether enclosed content was shut down gracefully or not. The iterator uses the \textit{async for} keyword and works asynchronously as well. It pauses the code until a new event can be produced.

The code starts at line \num{26} and runs the \textit{main()} function. This function first connects to the host by creating a context in line \num{15} where the \textit{connection} variable is usable. Using this connection the sensor platform is queried in line \num{16}. In comparison to the precious example it is now far easier to follow the program flow because the context and generator reveal what is happening next. Unfortunately, reading the sensor still requires passing it to a new task because the generator will keep generating more sensors in the meantime. Although the code is split into multiple tasks, the nesting of callbacks as in the previous example is resolved and the readability of the code has improved tremendously.

The only problem is the error handling, because these worker tasks do not communicate with the original task that created them. This can be solved using the so-called observer pattern where an observer tasks watches the workers and handles such event. Directly implementing this pattern creates a myriad of events and event handler registrations. Missing one such event can break the whole program and leads to bugs that are hard diagnose and fix. To simplify this pattern a stream based approach can be applied. The observable is treated as a stream of events which is being processed by the observer using a chain of operators and actions executed in a certain order. This is much like an assembly line where different tasks are executed as the product passes each station.

\subsubsection{Stream Design}
The previously mentioned observer pattern is often implemented using data streams representing the subjects observed while the consumers are the observers. Using functional programming style these data streams can be written in a very concise form as shown in the following version of the example program.
\lstinputlisting[firstline=5, firstnumber=1, frame=single, breakindent=.5\textwidth, frame=single, breaklines=true, numbers=left, xleftmargin=2em, numberstyle=\tiny, style=mypython]{source/lab_kraken_stream.py}

The program starts in line \num{14} and enters the \textit{main()} function. Here, a context is used again to open the connection to the sensor platform. The sensors are queried next and a stream is created to read the reply, then filter for the specified sensor which is then read and the result is printed.

Using this programming style the intent of the program is revealed immediately, even before starting the stream. The syntax uses was borrowed from the Python library \textit{aiostreams} \cite{aiostreams}, which is similar to ReactiveX \cite{reactivex}, a library developed by Microsoft to operate on data streams. The interesting aspect of this code is the use of the pipe operator which inject the result of one function into the next function as its parameters. This way a chain of function calls is created. The \textit{lambda} keyword denotes a small anonymous function, but regular functions can also be used. In combination with an operator like \textit{filter}, \textit{switchmap}, or \textit{print}, they dictate the program flow, hence the name functional programming. These operators need some introduction though. The \textit{filter} operator is simple to understand as it will only pass on the input when the function, to which the input is passed as well, returns true. The \textit{switchmap} operator is more interesting. It is a combination of a \text{map} operator and a \text{switch} operator. The former applies a function to the input and then passes on the output of the function, in this case \textit{read\_temperature} which creates an iterator. The latter operator will take its most recent input and iterate, producing temperature values. This operator will terminate the iteration when a new input is passed and then iterate the new input. This is handy as it automatically makes sure that there is only one reader per sensor and for example new sensor configurations can be injected into the data stream above the \textit{switchmap} which automatically replace the old sensor reader.

This style of programming was found to be ideal for real-time data processing, as it allows to continuously update configurations or add and remove sensors, or even hosts, without having to worry about what happens along the pipeline.

\subsubsection{Device Identifiers}
Every sensor network needs device identifiers. Preferably those identifiers should be unique. Typically a device has some kind of internal identifier. Here are a few examples of the sensors used in the authors network:
\begin{table}[ht]
    \centering
    \begin{tabularx}{0.95\textwidth}{|l|p{3cm}|X|}
        \hline
        Device Type& Identifiers& Example\\
        \hline
        GPIB (SCPI)& \textit{*IDN?}& \small{Keysight Technologies,34470A,MYXXXXXXXX,A.03.03-02.40-03.03-00.52-02-01} or\newline\small{Agilent Technologies,34410A,MYXXXXXXXX,A.03.03-02.40-03.03-00.52-02-01}\\
        \hline
        Tinkerforge& A base58 encoded integer device id& QE9 (163684)\\
        \hline
        LabNode& UUID & cc2f2159-e2fb-4ed9-8021-7771890b37ad\\
        \hline
    \end{tabularx}
    \caption{Device identifiers used by common devices found in the lab. The serial number of the Keysight \device{34470A} DMM was obscured on purpose.}
    \label{tab:common_device_ids}
\end{table}

As it can be seen in table \ref{tab:common_device_ids}, most of these identifiers do not guarantee to uniquely identify a device within a network. The Tinkerforge id is the weakest, as it is a \qty{32}{\bit} integer (\num{4294967295} options), which might easily collide with another id from a different manufacturer. For better readability the id is typically presented as a base58 encoded string. An encoder/decoder example can be found in the TinkerforgeAsync library \cite{TinkerforgeAsync}.

The id string returned by a SCPI device is slightly more useful, but again does not guarantee uniqueness. As per the SCPI specification it returns a string containing \textit{\$manufacturer,\$name,\$serial,\$revision}. Even when ignoring the software revision part which might change on update, the same device might return a different id depending on its settings. The id string shown in table \ref{tab:common_device_ids} relate to the same device, but the latter uses a compatibility flag in the settings.

The only reasonably unique id is used by the LabNodes. The universal unique identifier (UUID) or globally unique identifier (GUID), as dubbed by Microsoft, can be used for networks with participant numbers going into the millions. There are several versions defined in RFC 4122 \cite{rfc_uuid} and the LabNodes use version 4, which is a random \qty{128}{\bit} identifier with \qty{122}{\bit} of entropy. Of the remaining \qty{6}{\bit}, \qty{4}{\bit} are reserved for the UUID version and \qty{2}{\bit} for the variant. This allows to prove the usefulness as a unique id as below.

Calculating the chance of a collision between two random UUIDs is called the birthday problem \cite{BirthdayProblem} in probability theory. The probability of at least one collision in $n$ devices out of $M = 2^{122}$ possibilities can be calculated as follows:
\begin{align}
    p(n) &= 1 - 1 \cdot \left(1 - \frac{1}{M}\right) \cdot \left(1 - \frac{2}{M}\right) \dots \left(1 - \frac{n-1}{M}\right) \nonumber\\
    &= 1 - \prod_{k=1}^{n-1} \left(1 - \frac{k}{M} \right)
\end{align}
Using the Taylor series $e^x = 1+x \dots$, assuming $n \ll M$ and approximating we can simplify this to:
\begin{align}
    p(n) &\approx 1 - \left(e^\frac{-1}{M} \cdot e^\frac{-2}{M} \dots e^\frac{-(n-1)}{M} \right) \nonumber\\
    &\approx 1 - \left(e^\frac{-n(n-1)/2}{M} \right) \nonumber\\
    &\approx 1 - \left(1 - \frac{n^2}{2 M} \right) = \frac{n^2}{2 M}
\end{align}
For one million devices using random UUIDs, this gives a probability of about \num{2e-25}, which is negligible.

In the LabKraken implementation, all devices, except for the LabNodes which already have a UUID, will be mapped to UUIDs using the underlying configuration database. It is up to the user to ensure the uniqueness of the non-UUID ids reported by the devices to ensure proper mapping. These UUIDs can then be used to address and configure each device on the sensor network.

%\subsubsection{Ethernet Bus and Synchronous Buses}
%There are inherent challenges involved with the Ethernet bus for instrumentation. The Ethernet bus is intrinsically asynchronous and multiple controllers can talk to the device at the same time. Not only that, but different processes within the same controller can talk to the same device. This makes deterministic statements about the device state challenging. A device that is not designed to work asynchronously in the first place may have trouble with multiple requests coming in from different clients. This must be kept in mind when using serial adapters like USB or GPIB to Ethernet.

%While it is impossible to rule out the possibility of multiple controllers on a network, care was taken to synchronize the workers within Kraken.
%\subsection{Databases}
%\subsubsection{Cardinality}
%\begin{itemize}
% \item TimescaleDB vs Influx
% \item Example Sensors vs. Experiment
%\end{itemize}

\clearpage
\section{Short Introduction to Control Theory}
This section will give a very brief introduction to some basic concepts of control theory. Many systems require control over one or more process variables. For example, temperature control of a room or a device, or even creating a programmable current from a voltage is one such problem. All of this requires control over a process and is established through feedback, which allows a controller to be aware of the state of the system.

The focus of this section is narrowed down to the concept of feedback and control with regard to developing and understanding PID controllers for temperature control. Simpler feedback loops like those typically used around op-amps will not be primarily considered in this section and are discussed in the relevant part of the documentation. In the following sections, first general properties of the Laplace transform and useful relationships are introduced, then, a model for the system and its controller will be developed, finally, using the model, tuning of the control parameters using different tuning algorithms will be discussed.

\subsection{Introduction to the Transfer Function and the Laplace Domain}%
\label{sec:transfer_function}
There are two types of control systems: open- and closed-loop systems. A system is called open loop, if the output of a system does not feed back to its input as in figure \ref{fig:open_loop}. On the other hand, if the output influences the input of the system via feedback, it is called a closed-loop system, as shown in figure \ref{fig:closed_loop}. Although feedback can be treated in static systems, it is more useful to treat it in dynamic systems, either in the time-domain or the frequency-domain. To discuss these systems, the terminology used in the following section needs to be defined. $G(s)$ is called the transfer function of the system, while $U(s)$ is the input, $Y(s)$ is the output, $s$ is a complex frequency domain variable, $\beta$ is the feedback parameter, also called feedback fraction, as shown in figure \ref{fig:closed_loop}. In this section, upper case letters are used to denote functions in the Laplace domain, while lower case letters are referring to functions in the time domain. Normally, the transfer function is denoted $H(s)$ but to prevent confusion with the Heaviside function $H(t)$, the letter $G$ is used here for the transfer function. In later chapters the common form $H(s)$ is used.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{open_loop.tex}
        \caption{Open-loop system.}
        \label{fig:open_loop}
    \end{subfigure}
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{closed_loop.tex}
        \caption{Closed-loop system.}
        \label{fig:closed_loop}
    \end{subfigure}
    \caption{Block diagram of a closed- and an open-loop control system.}
    \label{fig:feedback_systems}
\end{figure}

It is convenient to express the transfer function in its Laplace transform for several reasons that will be explained below. The systems to be discussed are physical system and hence are causal. That means the output only depends on past and present inputs, but not future inputs. For this reason, only the one-sided or unilateral Laplace transform needs to be considered. It is defined as:
\begin{equation}
    \mathscr{L}\left( f(t) \right) = F(s) = \int_0^\infty f(t) e^{-st}\,dt.
\end{equation}

with $f: \mathbb{R}^+ \to \mathbb{R}$ that is integrable and grows no faster than $c \cdot e^{s_0t}$ for $s_0, c \in \mathbb{R}$. The latter attribute is important for deriving the rules of differentiation and integration.

To understand the benefits of using the Laplace representation of the transfer function, a few useful properties should be discussed with regard to the PID controller. First of all, the Laplace transform is linear:
\begin{align}
    \mathscr{L}\left(a \cdot f(t) + b \cdot g(t) \right) &= \int_0^\infty (a \cdot f(t) + b \cdot g(t)) e^{-st}\,dt \nonumber\\
    &= a \int_0^\infty f(t) e^{-st}\,dt + b \int_0^\infty g(t) e^{-st}\,dt \nonumber\\
    &= a \mathscr{L}\left(f(t)\right) + b \mathscr{L}\left(g(t)\right)
\end{align}

Another interesting property is the derivative and integral of a function $f$. The function $f$ must, of course, be differentiable and grow no faster than the exponential function as defined above:
\begin{align}
    \mathscr{L}\left(\frac{df}{dt}\right) &= \int_0^\infty \underbracket{f'(t)}_{v'(t)} \underbracket{\vphantom{f'(t)}e^{-st}}_{u(t)}\,dt \nonumber\\
    &= \left[e^{-st} f(t) \right]_0^\infty - \int_0^\infty (-s)f'(t)\,dt \nonumber\\
    &= -f(0) + s \int_0^\infty f'(t)\,dt \nonumber\\
    &= s F(s) - f(0)
\end{align}

\begin{align}
    \mathscr{L} \left( \int_0^t f(\tau)\,d\tau \right) &= \int_0^\infty \left(\int_0^t f(\tau)\,d\tau e^{-st} \right)\,dt \nonumber\\
    &= \int_0^\infty \underbracket{e^{-st}\vphantom{\int_0^t}}_{v'(t)} \underbracket{\int_0^t f(t)\,d\tau}_{u(t)}\,dt \nonumber\\
    &= \left[\frac{-1}{s} e^{-st} \int_0^t f(t)\,d\tau \right]_0^\infty - \int_0^\infty \frac{-1}{s} e^{-s\tau} f(\tau)\,d\tau \nonumber\\
    &= 0 + \frac{1}{s} \int_0^\infty e^{-s\tau} f(\tau)\,d\tau \nonumber\\
    &= \frac{1}{s} F(s) \label{eqn:lapace_integration}
\end{align}

If the initial state $f(0)$ can be chosen to be $0$, the differentiation becomes a simple multiplication by $s$, while the integration becomes a division by $s$. These three properties greatly simplify the calculations required for studying a proportional–integral–derivative controller in section \ref{sec:pid_controller_basics}.

Finally, the most important aspect is the possibility to give a simple relation between the input $u(t)$ and the output $y(t)$ of a system. This relationship between input and output of a system as shown in figure \ref{fig:open_loop} is given by the convolution, see e.g. \cite{pid_basics}. Assuming the system has an initial state of $0$ for $t<0$, hence $u(t<0) = 0$ and $g(t<0) = 0$, one can calculate:
\begin{equation}
    y(t) = (u \ast g)(t) = \int_0^\infty u(\tau) g(t-\tau)\,d\tau
    \label{eqn:convolution}
\end{equation}

Applying the Laplace transform, greatly simplifies this:
\begin{align}
    Y(s) &= \int_0^\infty e^{-st} y(t)\,dt \nonumber\\
    \overset{\ref{eqn:convolution}}&{=} \int_0^\infty \underbrace{e^{-st}}_{e^{-s(t-\tau)}e^{-s\tau}} \int_0^\infty u(\tau) g(t-\tau)\,d\tau\,dt \nonumber\\
    &= \int_0^\infty \int_0^t e^{-s(t-\tau)} e^{-s\tau} g(t-\tau) u(\tau)\,d\tau\,dt \nonumber\\
    &= \int_0^\infty e^{-s\tau} u(\tau)\,d\tau \int_0^\infty e^{-st} g(t)\,dt \nonumber\\
    &= U(s) \cdot G(s)
\end{align}

This formula is a lot simpler than the convolution of $u(t)$ and $g(t)$, therefore the use of the Laplace transform has become very popular in control theory.

Having derived some of the most useful properties, it is interesting to look at a few functions, which are heavily used in control theory, like a function delayed by the time interval $\theta$. To demonstrate its properties, let $f(t-\theta)$ be
\begin{equation}
    g(t) \coloneqq \begin{cases} f(t-\theta), & t \geq \theta \\ 0, & t < \theta \end{cases} \,. \label{eqn:delayed_f}
\end{equation}

The reason for this definition is that it is mandatory for the system to be causal. This means, it is impossible to get information from the future ($t<\theta$). To satisfy this requirement, any constant other than \num{0} may be chosen, as is done later in section \ref{sec:pid_tuning_rules}, when determining tuning parameters and fitting experimental data to a model. An example of such a time delayed function $g(t)$ is shown in figure \ref{fig:heaviside_delayed}.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.4\linewidth}
        \centering
        \scalebox{0.75}{%
            \import{figures/}{laplace_no_delay.tex}
        } % scalebox
        \caption{Original signal $f(t)$.}
        \label{fig:heaviside}
    \end{subfigure}
    \begin{subfigure}{0.4\linewidth}
        \centering
        \scalebox{0.75}{%
            \import{figures/}{laplace_time_delay.tex}
        } % scalebox
        \caption{Delayed signal $f(t-2)$.}
        \label{fig:heaviside_delayed}
    \end{subfigure}
\end{figure}

The Laplace transform of a delayed signal $g(t)$ can be calculated as follows:
\begin{align}
    \mathscr{L}\left( g(t) \right) &= \int_0^\infty g(t) e^{-st}\,dt \nonumber\\
    \overset{\ref{eqn:delayed_f}}&{=} \int_\theta^\infty f(t-\theta) e^{-st}\,dt \nonumber\\
    \overset{\tau \coloneqq t-\theta}&{=} \int_0^\infty f(\tau) e^{-s(\tau+\theta)}\,d\tau \nonumber\\
    &= e^{-s\theta} \int_0^\infty f(\tau) e^{-s\tau} \,d\tau \nonumber\\
    &= e^{-s\theta} F(s) \label{eqn:laplace_delayed}
\end{align}

To satisfy the causality requirement in the time domain, the Heaviside function $H(t)$ can be used to give a more concise representation of $g(t)$:
\begin{align}
    \mathscr{L}\left( f(t-\theta) H(t-\theta) \right) = e^{-s\theta} F(s) \label{eqn:laplace_causality}
\end{align}

Lastly, the Laplace transform of $e^{at}$ is given, which is commonly used in differential equations:
\begin{align}
    \mathscr{L}\left(e^{at} \right) &= \int_0^\infty e^{(a-s)t}\,dt = \frac{1}{a-s} \left[e^{(a-s)t} \right]_0^\infty = \frac{1}{s-a} \label{eqn:laplace_exponential}
\end{align}

Using these tools, it is possible calculate the transfer function of a closed-loop temperature controller, which will be done in the next section.

\subsection{A Model for Temperature Control}%
\label{sec:temperature_control_model}
\begin{figure}[ht]
    \centering
    %\scalebox{1}{%
        \import{figures/}{first-order_model.tex}
    %} % scalebox
    \caption{Simple temperature model of a generic system.}
    \label{fig:first-order_model_room}
\end{figure}

In order to describe a closed-loop system using a transfer function $G(s)$, one has to first create a model for the process and the controller involved. This section will derive the simple, but very useful first order model with dead-time. This model can be derived from the idea that the system at temperature $T_{system}$ has a thermal capacitance $C_{system}$, an influx of heat $\dot Q_{load}$ from a thermal load and a controller removing heat from the system through a heat exchanger with a resistance of $R_{force}$. Additionally, there is some leakage through the walls of the system to the ambient environment via $R_{leakage}$. This analogy of thermodynamics with electrodynamics allows to create the model shown in figure \ref{fig:first-order_model_room}. Since this model is to be used for a temperature controller, more simplifications can be made and a so-called small-signal model can be developed as opposed to the large-signal model shown above. The small-signal model is an approximation around a working point that is valid for small deviations around it, similar to a Taylor approximation. The small-signal model can be used to calculate the system response to small changes of the controller output in order to estimate the control parameters at that working point.

Using the small-signal approach, the system response can be split into a constant and a dynamic part: the 0\textsuperscript{th} and 1\textsuperscript{st} order of the Taylor approximation. In order to simplify the system shown in figure \ref{fig:first-order_model_room}, the assumption can be made that the system load $\dot Q_{load}$ and the flux through $R_{leakage}$ is \textit{reasonably stable}. \textit{Reasonably stable} means that it can be treated as small deviations and additionally any changes are within the bandwidth of the controller and well suppressed. This allows to treat them as (almost) constant effects. Constant effect can be neglected in the small-signal model because they only manifest as a constant offset applied to the output of the controller and have no dynamics. This leaves only the room with its heat capacity and the heat exchanger as dominant factors in the small-signal model shown in figure \ref{fig:first-order_model}. Here $T_{force}$ and $T_{system}$ were replaced by $T_{in}$ and $T_{out}$, $R_{force}$ and $C_{system}$ were replaced by $R$ and $C$ for better readability.
\begin{figure}[htb]
    \centering
    \scalebox{1}{%
        \import{figures/}{first-order_model_kirchhoff.tex}
    } % scalebox
    \caption{Simplifications of the temperature model of a room lead to this first order model.}
    \label{fig:first-order_model}
\end{figure}

This is the classic $RC$ circuit. To calculate the transfer function, a relationship between $T_{in}$ and $T_{out}$ is required and exploiting the analogy of thermodynamics and electrodynamics again, using Kirchhoff's second law, following the arrow in figure \ref{fig:first-order_model} one finds:
\begin{alignat}{1}
    \sum T_i &= 0 \nonumber\\
    T_{in}(t) - \dot{Q}(t) R - \frac 1 C \int \dot{Q}(t)\,dt &= 0 \label{eqn:first-order_model_kirchhoff}
\end{alignat}

Taking the Laplace transform, applying equation \ref{eqn:lapace_integration}, solving for $ \dot Q(s)$ and using $T_{out} = \frac{1}{sC} \dot Q(s)$ to replace $\dot Q$, equation \ref{eqn:first-order_model_kirchhoff} can be written as:
\begin{align*}
    T_{in}(s) - \dot{Q}(s) R - \frac{1}{sC} \dot{Q}(s) &= 0\\
    \dot{Q}(s) = \frac{T_{in}(s)}{R-\frac{1}{sC}} &= \frac{T_{out}}{\frac{1}{sC}}
\end{align*}

This allows to calculate the transfer function of the process $P$ as
\begin{align}
    P(s) &= \frac{T_{out}}{T_{in}} = \frac{\frac{1}{sC}}{R-\frac{1}{sC}} \nonumber\\
    &= \frac{1}{sRC + 1} \nonumber\\
    &= \frac{1}{1 + s\tau} = \frac{K}{1 + s\tau}\,. \label{eqn:first-order_model}
\end{align}
with the system gain $K$ and the time constant $\tau$. In case of the $RC$ circuit, the gain is $1$, but other systems may have a gain factor of $K \neq 1$. This is generally the case when using any type of sensor that converts the measurand into the input signal. $K$ is therefore included here for the sake of generality.

Equation \ref{eqn:first-order_model} is called the transfer function of a first order model, because its origin is a differential equation of first order. This model describes homogeneous systems like a room very well, as can be seen in section \ref{sec:pid_controller_tuning}, but in order to derive the transfer function including the controller and the sensor some more work is required to derive the sensor transfer function.

Expanding on figure \ref{fig:open_loop} and equation \ref{eqn:convolution} the open-loop transfer function of the process and its sensor becomes:
\begin{equation}
    G(s) = P(s) \cdot S(s)
\end{equation}
and the block diagram changes to
\begin{figure}[htb]
    \centering
    %\scalebox{1}{%
        \import{figures/}{open_loop_full.tex}
    %}% scalebox
    \caption{Open-loop system with sensor.}
\end{figure}

The transfer function of the sensor, given an ideal linear transducer, can be modeled as a delay line with delay $\theta$ and $f(t-\theta) = H(t-\theta)$. A sensor gain of $1$ is assumed here, because any system gain is assumed be included in the parameter $K$ of the process transfer funciton. Using equation \ref{eqn:laplace_delayed}, $S(s)$ can be written as
\begin{equation}
    S(s) = e^{-\theta s} .
\end{equation}

The full system model including the time delay can now be written as:
\begin{equation}
    G(s) = \frac{K}{1 + s\tau} e^{-\theta s} \label{eqn:first-order_plus_dead_time_model}
\end{equation}

This is called a first order plus dead-time model (FOPDT) or first order plus time-delay model (FOPTD). While the Laplace representation is useful to mathematically explore the mode, in order to fit experimental data to this model it is more convenient to transform the transfer function \ref{eqn:first-order_plus_dead_time_model} into the time domain. To have a meaningful result, an input $U(s)$ is required, because $G(s)$ is only a transformation. In principal, any function can serve this purpose, but a step function is typically used, for example by \citeauthor{ziegler_nichols} \cite{ziegler_nichols} and many others \cite{tuning_rules,pessen_integral,simc,simc_paper,pid_controllers_for_time_delay_systems,pi_stabilization_of_fopdt_systems, pid_basics}. The step function is both simple to calculate and to apply to a real system in form of a controller output change. This technique is also explored in more detail in section \ref{sec:pid_controller_tuning}. Using equations \ref{eqn:laplace_delayed} and \ref{eqn:laplace_exponential}, the Heaviside $H(t)$ step function transforms into
\begin{equation}
    \mathscr{L} \left(u(t) \right) = U(s) = \mathscr{L} \left( \Delta u H(t) \right) = \frac{\Delta u}{s}
\end{equation}
with the step size $\Delta u$. The output response $Y(s)$ of the system to the step can then be calculated analytically.
\begin{align}
    Y(s) &= U(s) \cdot G(s)\nonumber\\
    &= \frac{\Delta u}{s} \frac{K}{1 + s\tau} e^{-\theta s} \nonumber\\
    &=  K \Delta u \frac{1}{s (1 + s\tau)} e^{-\theta s} \nonumber\\
    &= K \Delta u \left(\frac{1}{s} - \frac{\tau}{s\tau+1} \right) e^{-\theta s} \nonumber\\
    &= K \Delta u \left(\frac{1}{s} - \frac{1}{s+\frac{1}{\tau}} \right) e^{-\theta s}
\end{align}

To derive $y(t)$, the inverse Laplace transform of $Y(s)$ is required. Unfortunately, this is not as simple as the Laplace transform. Fortunately though the required equations were already derived in equations \ref{eqn:lapace_integration} and \ref{eqn:laplace_exponential}. Making sure causality is guaranteed as shown in equation \ref{eqn:laplace_causality}, the simple first order model can be transformed back into the time domain.
\begin{align}
     y(t) &= \mathscr{L}^{-1} \left(Y(s)\right) \nonumber\\
     &= K \Delta u \mathscr{L}^{-1} \left(\frac{1}{s} e^{-\theta s} \right)  - K \mathscr{L}^{-1} \left( \frac{1}{s+\frac{1}{\tau}} e^{-\theta s} \right) \nonumber\\
    \overset{\ref{eqn:laplace_exponential}}&{=} K \Delta u \cdot 1 \cdot H(t-\theta) - \left(e^{-\frac{t-\theta}{\tau}} \right) H(t-\theta) \nonumber\\
    &= K \Delta u \left(1-e^{-\frac{t-\theta}{\tau}} \right) H(t-\theta) \label{eqn:first-order_plus_dead_time_model_time-domain}
\end{align}

The time domain solution of the FOPDT model can now be used to extract the parameters $\tau$, $\theta$ and $K$ from a real physical system.

The procedure can be summarised from the above as follows. The controller must be set to a constant output and the room must be given time to reach equilibrium. Once the temperature has settled, an output step of $\Delta u$ is applied. The system will respond after a time delay and then follow an exponential function. A simulation of the step response applied to a first order model with time delay is shown in figure \ref{fig:fopdt}. The gain is $K=1$. The solid black line shows the response of the transfer function, including the system and the sensor. The dashed lines show the individual components, the Heaviside function governing the delay and the exponential term of the system. The controller output step $\Delta u = 1$ is applied at $t=0$ and not shown explicitly. From figure \ref{fig:fopdt} it can be clearly seen that the sensor does not register a change until the time delay $\theta$ has passed and the Heaviside function changes from $0$ to $1$. Then the system responds with an exponential decay towards \num{1}.
\begin{figure}[ht]
    \centering
    \input{images/FOPDT_theory.pgf}% data/simulations/sim_laplace_fopdt.py
    \caption{Time domain plot of a first order plus dead time model showing individual components of the model and the composite function $y(t)$. Model parameters used: $K= \Delta u = 1$, $\tau=2$, $\theta=4$.}
    \label{fig:fopdt}
\end{figure}

So far, only open-loop systems were discussed. Using the FOPDT model, the system parameters can now be extracted from an existing system using a fit to the time domain reaction of such a system to a step input. Having extracted the system parameters, the next step is to design a controller around the system and close the loop to realize a controlled system. This is shown in the next section.

\subsection{PID Controller Basics}%
\label{sec:pid_controller_basics}
While there are many different types of controllers, like the bang–bang controller utilized in many temperature controller, which turns on at a certain threshold and turns off at another one, producing the saw-tooth shaped room temperature curve shown in figure \ref{fig:lab_temperature_start_of_project}, a continuous control system is desired to keep fluctuations to a minimum. The most commonly used controller type for non-integrating systems is the proportional–integral–derivative (PID) controller \cite{pid_in_industry}. A non-integrating system is a system without memory whose steady state does not depend on previous inputs. The advantage of applying a PID controller is that the controller does not need any special knowledge about the system model. A universal PID is simple to implement and can be tuned to control a wide range of different systems. While there are many variations of the PID algorithm \cite{pid_controller}, this section only introduces the basic, parallel, PID controller commonly used in digital implementation and deals with some of the shortcomings in practical applications.

In order to extend the FOPDT system derived in the previous section \ref{sec:temperature_control_model}, with the PID controller one has to move to a closed-loop system. Adding to figure \ref{fig:closed_loop} and inserting a new control block into the transfer function yields figure \ref{fig:closed_loop_pid}.
\begin{figure}[ht]
    \centering
    \scalebox{1}{%
        \import{figures/}{closed_loop_pid.tex}
    }% scalebox
    \caption{Closed-loop system with a PID controller.}
    \label{fig:closed_loop_pid}
\end{figure}

The error signal $E(s)$ used by the PID controller is the difference between the setpoint and the control parameter, in this case the room temperature. The transfer function of the PID controller can be split into three parts. A proportional part that is proportional to the error representing the present state, an integral part that is proportional to the accumulated error, representing the past state, and a derivative part that is proportional to the change in the error, extrapolating into the future. Analytically, it can be written as
\begin{align}
    c(t) &= k_p e(t) + k_i \int_0^t e(\tau) \,d\tau + k_d \frac{\mathrm{d}e(t)}{\mathrm{d}t} \label{eqn:pid_controller}\\
    C(s) &= k_p + k_i \frac{1}{s} + k_d s \,. \label{eqn:pid_controller_laplace}
\end{align}

The following discussion will mostly focus on equation \ref{eqn:pid_controller}, because, the time-domain equation is the one that can be implemented in software. As hinted above, there are a few shortcomings with the classic PID equation, when used in a real system which requires dynamic changes of the the setpoint or $k_i$.

The first problem that needs addressing is occurring when changing the PID parameter $k_i$, because equation \ref{eqn:pid_controller} is given for a time-independent $k_i$. Assuming a settled system without external disturbances, the output is fully determined by the integrator value because the error is zero. Now, when $k_i$ is changed, the output immediately changes, due to the change of the integral term. This is unintended. To fix it, the integral term must be changed to
\begin{equation}
    k_i \int_0^t e(\tau) \,d\tau \Rightarrow \int_0^t k_i(\tau) e(\tau) \,d\tau \,.
\end{equation}

This way, when adjusting $k_i$, its new value is applied to future error values only and there is no sudden kick.

The next issue is called \text{derivative kick}. When looking at the derivative part of equation \ref{eqn:pid_controller}, it can be seen that when instantly changing the setpoint, as in a step function, $\frac{\mathrm{d}e(t)}{\mathrm{d}t} \to \infty$. This behaviour is not intended and to fix this, the derivative part can be modified as follows.
\begin{align}
    \frac{\mathrm{d}e(t)}{\mathrm{d}t} &= \frac{\mathrm{d}\left(u(t) - y(t)\right)}{\mathrm{d}t} = \underbrace{\cancel{\frac{\mathrm{d}u(t)}{\mathrm{d}t}}}_{\to \infty} - \frac{\mathrm{d}y(t)}{\mathrm{d}t} \nonumber\\
    &=- \frac{\mathrm{d}y(t)}{\mathrm{d}t}
\end{align}

The new derivative term is equal to the unmodified one, except in the case of setpoint changes. Removing the setpoint from the equation, the controller behaves as intended. This solution is sometimes called \textit{derivative on measurement} as opposed to \textit{derivative on error}.
\begin{figure}[hb]
    \centering
    \input{images/sim_pid_controller_bode.pgf}% data/simulations/sim_pid_controller_bode.py
    \caption{Magnitude plot over frequency of the PID controller transfer function. Both the ideal PID controller and the PID controller with a filtered derivative are shown.}
    \label{fig:sim_pid_controller}
\end{figure}

Another issue is caused by the derivative term with noisy inputs. Assuming there is a very short input spike due to noise, the differential of the derivative term will again be sent to very high values, pushing the output away from the correct value and forcing the controller to slowly rebalance.

To further discuss the problem and its solution it is best to visit the frequency domain and visualize the transfer function of the PID controller as shown in figure \ref{fig:sim_pid_controller}. The ideal PID controller without filtering of the derivative displays a very strong response to low frequency inputs. This is due to the integral action, which removes any (constant) offset. It needs to have infinite gain at DC to push the offset to zero. In reality this is limited by input noise. Then follows a plateau with a magnitude of $k_p$ for the proportional term and finally the differential gain starts growing in magnitude and keeps steadily growing with rising frequency, just as expected.

With some knowledge about the process or the sensor it is possible to define an upper frequency, above which inputs become unrealistic and must therefore be unwanted noise. By filtering the derivative term with a first order filter causes it to roll off and its gain becomes constant as shown in figure \ref{fig:sim_pid_controller}. By adding the filter the PID controller transfer function changes to

\begin{equation}
    C(s) = k_p + k_i \frac{1}{s} + \frac{k_d s}{1 + s \alpha k_d} \,. \label{eqn:pid_controller_filtered}
\end{equation}

Typically $\alpha$ is in the range of \numrange{0.05}{0.2} \citep[p. 129]{pid_controller}.

An alternative is to filter the whole input. Depending on the filter cutoff, there is not much difference to equation \ref{eqn:pid_controller_filtered}, because the filter will not touch the proportional and integral part of the transfer function if both are well within its passband.

From figure \ref{fig:sim_pid_controller} it can also be seen, why in some publications, the gain $k_p$ is applied to all three terms and $k_i$ and $k_d$ are replaced with $T_i$ and $T_d$ to accommodate for that.
\begin{equation}
    C(s) = k_p \left(1 + \frac{1}{T_i s} + \frac{T_d s}{1 + s \alpha T_d} \right) \label{eqn:pid_controller_series}
\end{equation}

Using this form allows to shift the overall gain up and down keeping its shape instead of just the $k_p$ part, thus changing the corner frequencies. The alternative form is only given here for the sake of completeness. The author uses the ideal form shown in equation \ref{eqn:pid_controller_laplace} with the parameters $k_p$, $k_i$, and $k_d$ wherever possible.

This concludes the discussion of the PID controller and the introduction of the basic terms. It now begs the question how the controller interacts with the system and how to derive the optimal PID parameters from a given system or model. Thus, the next section discusses controller tuning rules and their effect on the system performance.

\subsection{PID Tuning Rules}%
\label{sec:pid_tuning_rules}
While many PID tuning rules can be found in the literature, their application depends on the underlying system and the desired system response. This section will discuss several proposed solutions and compare them to the authors use case. The section aims to give a simple method to determine decent PI/PID parameters for the applications found in the lab. Among the methods discussed are the most classic set of tuning rules developed by \citeauthor{ziegler_nichols} \cite{ziegler_nichols}, and an improved version called Skogestad Internal Model Control (SIMC) presented by \citeauthor{simc_paper} \cite{simc_paper} which promises better performance for non-integrating systems. These rules all include simple instructions to extract the necessary parameters using pen and paper. Using a computer and fitting algorithms, the bar for \textit{simple} has been raised considerably, so more complex approaches can be undertaken which extract more parameters from the system. Using these additional parameters, more precise control is promised by \citeauthor{pid_basics} \cite{pid_basics, advanced_pid_control} with a method called AMIGO. Finally, it is possible to shape the control loop to result in a desired transfer function. This technique is mostly used in motor control \cite{pid_controller,advanced_pid_control} and also requires the model parameters.

All of these rules will be compared against a demo model of a room to explain the details. It is the first order model with delay which was derived in equation \ref{eqn:first-order_plus_dead_time_model}. The discussion is limited to the FOPDT model, because the systems treated in this work could be modelled very well using this equation. Higher-order models are discussed in more details for example in \cite{advanced_pid_control,pid_controller,simc_paper}, in case the reader encounters such a system and feels the need to extract the model parameters.
\begin{equation}
    G(s) = \frac{K e^{-\theta s}}{1 + s \tau} \label{eqn:demo_process_model}
\end{equation}

The following parameters were extracted from lab 011 of the APQ group, using the techniques shown in section \ref{sec:temperature_control_model} using equation \ref{eqn:first-order_plus_dead_time_model_time-domain}. The details are discussed in section \ref{sec:pid_controller_tuning}. The system gain $K$ was scaled to the full scale output (\qty{4095}{\bit}) of the controller, hence the somewhat strange unit \unit[per-mode=power]{\K \bit\per\bit}.
\begin{table}[hb]
    \centering
    \begin{tabular}{ccc}
        \toprule
        Gain K& Lag $\tau$& Delay $\theta$ \\
        \midrule
        \qty[per-mode=power]{13.07}{\K \bit\per\bit}& \qty{395}{\s}& \qty{187}{\s}\\
        \bottomrule
    \end{tabular}
    \caption{Example paramters extracted from lab 011 using the techniques shown here and as applied in section \ref{sec:pid_controller_tuning}.}.
    \label{tab:pid_example_model}
\end{table}

Before detailing the tuning parameters, the loop shaping method will be explained first, because it cannot only be used to derive custom rules but was also used to create the SIMC rules proposed by \citeauthor{simc_paper} \cite{simc_paper}. The aim of this method is to derive a controller that shapes the model in such a way that a desired system response to setpoint changes is achieved. A general closed-loop system with a controller $C$ and a system $G$ is shown in figure \ref{fig:closed_loop_controller}. This will be used as a basis to find the required controller for a desired transfer function $\frac{Y(s)}{U(s)}$.
\begin{figure}[ht]
    \centering
    \scalebox{1}{%
        \import{figures/}{closed_loop_controller.tex}
    }% scalebox
    \caption{Closed-loop system $G$ with a controller $C$.}
    \label{fig:closed_loop_controller}
\end{figure}

Starting with the transfer function of the controlled system, made up of the controller and the system, most experimenters would, at least in a feverish dream, prefer a transfer function of the following divine form
\begin{equation*}
    \frac{Y(s)}{U(s)} = 1 \,,
\end{equation*}
but unfortunately life is more profane and there is no controller that will always (and with warp speed) force a system to a certain setpoint. One may therefore settle for the second-best choice, a first order low pass with a slow roll-off and a small delay, which must be added to ensure causality. One therefore arrives at
\begin{equation}
    \frac{Y(s)}{U(s)} = \frac{e^{-\theta s}}{1 + s \tau_c}\,, \label{eqn:desired_transfer_function}
\end{equation}
where $\tau_c$ is the closed-loop time constant and a measure for the aggressiveness of the controller. A small $\tau_c$ results in a more aggressive controller with faster response.

For the system shown figure \ref{fig:closed_loop_controller} the closed-loop transfer function is found to be
\begin{align*}
    \frac{Y(s)}{U(s)} &= \frac{C(s) G(s)}{C(s) G(s) + 1} \\
    \Rightarrow C(s) &= \frac{1}{G(s)} \frac{1}{\frac{Y(s)}{U(s)} -1}
\end{align*}

This loop now needs to be shaped into the desired transfer function given in equation \ref{eqn:desired_transfer_function}, so substituting $\frac{Y(s)}{U(s)}$ yields
\begin{align}
    C(s) &= \frac{1}{G(s)} \frac{e^{-\theta s}}{s \tau_c +1 - \underbrace{e^{-\theta s}}_{\approx 1 - \theta s}}\\
    &\approx \frac{1}{G(s)} \frac{e^{-\theta s}}{s (\tau_c + \theta)} \,.
\end{align}

$e^{-\theta s}$ was approximated using a first order Taylor expansion. The desired controller response now only depends on the system (including the sensor) to be controlled. So, substituting the system equation \ref{eqn:demo_process_model} results in
\begin{align}
    C(s) &= \frac{1}{K} \frac{s \tau + 1}{(\tau_c + \theta) s} \nonumber\\
    &= \underbrace{\frac{1}{K} \frac{\tau}{\tau_c + \theta}}_{k_p} + \underbrace{\frac{1}{K} \frac{1}{\tau_c + \theta}}_{k_i} \frac{1}{s}\,.
\end{align}

This is a PI controller with $k_p = \frac{1}{K} \frac{\tau}{\tau_c + \theta}$ and $k_i = \frac{1}{K} \frac{1}{\tau_c + \theta}$. From these calculations, it can be seen that a first order model can be fully treated using a PI controller. Second order (and higher order) models typically necessitate a PID or more sophisticated controller for optimal control. The problems discussed in this work mainly focus temperature control of (mostly) homogeneous objects, so the focus lies on the PI controller for most of the remaining section but the ideas and simulations can similarly be applied to the PID controller as well. Any caveats to be expected when treating a PID instead of a PI controller will be discussed.

Using the loop shaping technique, it is fairly easy to derive custom rules in case the model parameters can be extracted. As mentioned above, one such loop-shaped tuning rule is the SIMC rule set and the authors of those rules give advice for an ample variety of different models and also investigate the parameter choice regarding stability, load, and setpoint disturbances. Before attempting a custom approach, it is therefore recommended to check \cite{simc_paper} for an appropriate set of rules for more complex models in order to save time and effort.
\begin{table}
    \centering
    \begin{tabular}{lcccc}
        \toprule
        Tuning Rule& $k_p$& $T_i$ & $T_d$ & Source \\
        \midrule
        Z-N PI & $\frac{0.9 \tau}{K \theta}$ & $\frac{\theta}{0.3}$ & -- & \cite{ziegler_nichols}\\
        Z-N PID & $\frac{1.2 \tau}{K \theta}$ & $2 \theta$ & $\frac{\theta}{2}$ & \cite{ziegler_nichols}\\
        SIMC PI & $\frac{\tau}{K (\tau_c + \theta)}$ & $\min\left(\tau, 4 (\tau_c+\theta)\right)$ & -- & \cite{simc_paper}\\
        SIMC PID & $\frac{\tau_1}{K (\tau_c + \theta)}$ & $\min\left(\tau_1, 4 (\tau_c+\theta)\right)$ & $\tau_2$ & \cite{simc_paper}\\
        AMIGO PI & $\frac{0.15}{K} + \left(0.35 - \frac{\tau \theta}{\left(\tau + \theta\right)^2}\right) \frac{\tau}{K \theta}$ & $0.35 \theta + \frac{13 \tau^2 \theta}{\tau^2 + 12 \tau \theta + 7 \theta^2}$ & -- & \citep[p. 228]{advanced_pid_control}\\
        AMIGO PID & $\frac{1}{K} \left(0.2 + 0.45 \frac{\tau}{\theta}\right)$ & $\frac{0.4 \theta + 0.8 \tau}{\theta + 0.1 \tau} \theta$ & $\frac{0.5 \tau \theta}{0.3 \theta + \tau}$ & \citep[p. 233]{advanced_pid_control}\\
        \bottomrule
    \end{tabular}
    \caption{PI/PID parameters for different tuning rules. The PI controllers assume a first order model, the PID rules are required when dealing with a second order model.}
    \label{tab:pid_tuning_parameters}
\end{table}

For reasons of brevity, in table \ref{tab:pid_tuning_parameters}, the PID parameters are given as $k_p$, $T_i$ and $T_d$ as introduced in equation \ref{eqn:pid_controller_series}. $k_i$ and $k_d$ can be calculated from
\begin{align*}
    k_i &= \frac{k_p}{T_i}\\
    k_d &= k_p T_d\,.
\end{align*}

Regarding the SIMC PI/PID algorithm, \citeauthor{simc_paper} \cite{simc_paper} and \citep[ch. 5]{simc} suggests using $\tau_c = \theta$ for “\textit{tightest possible subject to maintaining smooth control}“. Following this recommendation, the minimum can be calculated from the parameters given in table \ref{tab:pid_example_model} on page \pageref{tab:pid_example_model} as $\min\left(\tau, 4 (\tau_c+\theta)\right) = \min\left(\tau, 8 \theta\right) = \tau$.

Using the rules above, the full system can be simulated now. This was done using Python. The simulation source code can be found in \external{data/simulations/sim\_pid\_controller.py} as part of the online supplemental material \cite{supplemental_material}. The simulation can be used to model arbitrary PI(D) controller and arbitrary models can be used as well. It allows to compare different settings before applying them to a real system. It also considerably shortens deployment times because especially for systems with long time scales, it becomes difficult to test several parameter sets on the fly, thus a simulation can reduce deployment time to a few minutes instead of hours.

The simulation emulates the PID controller developed for the lab temperature controller. By default it has a sampling rate of \qty{1}{\Hz}. The simulation  will apply a setpoint change of \qty{+1}{\K} \qty{10}{\s} into the simulation. After the simulation, it will plot the time domain response of the controlled system. The setpoint change in this scenario is very similar to the load disturbances that are expected. Typically a noise source is used here instead, but in contrast to the statistical noise, which could be used to test for disturbance rejection, the situation in labs are different and cannot be modelled with stationary noise. While there is some noise coming from the sensor and the lab, the major disturbances are usually caused by the experimenters instead of the lab itself. These are events like a device being switched on or off for an extended period of time, longer than the controller needs to settle. This is equivalent to a setpoint change in terms of the error term in equation \ref{eqn:pid_controller}, since there is no difference in the error term between a setpoint and a process variable change. Do note, that this is not true for the PID controller, whose derivative term directly works on the measurement (or process variable) as this was explicitly implemented above. For PID controllers, there is a difference between the setpoint change behaviour and system noise rejection. This must be kept in mind and tested accordingly.

Simulating the model above and using the PI parameters derived from table \ref{tab:pid_tuning_parameters}, gives the plot shown in figure \ref{fig:pid_controller_comparison}.
\begin{figure}[ht]
    \centering
    \input{images/sim_pid_controller_comparison.pgf}% data/simulations/sim_pid_controller.py
    \caption{Different PI Controllers tuned with parameter derived using the following methods: Ziegler-Nichols, SIMC and AMIGO. The system model is the FOPTD model for room 011.}
    \label{fig:pid_controller_comparison}
\end{figure}

As it can be seen in figure \ref{fig:pid_controller_comparison}, the Ziegler-Nichols tuning rule produces a very aggressive PI controller that shows quite a bit ringing, which is undesired for this application. The AMIGO rules are rather conservative, but do not produce any overshoot. The SIMC rules have proven the most useful for this application so far. This experience is in line with the results from \citeauthor{thesis_liebmann} \cite{thesis_liebmann}, who tested different PID tuning algorithms for their viability for temperature control in the labs discussed here.

To conclude, several PID tuning rules were presented and using a Python simulation tool it is possible to test a set of PID parameters before implementation. Using an example based on parameters extracted from a real environment, the different tuning rules were applied to a model for a real lab and the SIMC tuning rules were found to give the best results for this application. The reader should now be able to extract the model parameters from physical systems and have the tools to choose an optimal set of tuning parameters for the PID controller. Further reading recommendations are for a broad overview \cite{pid_controller}, and for more details \cite{advanced_pid_control}.

\clearpage
\section{Noise and Allan Deviation}%
\label{sec:allan_deviation}
The Allan variance \cite{adev} $\sigma_A^2(\tau)$ is a two-sample variance and used as a measure of stability. The Allan deviation $\sigma_A(\tau)$ is the square root of the variance. Originally, the Allan variance was used to quantify the performance of oscillators, namely the frequency stability, but it can be used to evaluate any quantity. In order to define the Allan variance, a few terms need to be defined first. A single measurement value of a time series $y(t)$ can be written as
\begin{equation}
    \bar y_k(t) = \frac{1}{\tau} \int_{t_{k}}^{t_{k}+\tau} y(t)\,dt . \label{eqn:allan_variance_measurement}
\end{equation}
This is the $k$-th measurement with a measurement time or integration time $\tau$. The latter term is frequently used for digital multimeters (DMM). $t_k$ is the start of the $k$-th sampling interval including the dead time $\theta$
\begin{equation}
    t_{k+1} = t_k + T
\end{equation}
with
\begin{equation}
    T \coloneqq \tau + \theta .
\end{equation}

\begin{figure}[hb]
    \centering
    \scalebox{1}{%
        \import{figures/}{allan_variance_definitions.tex}
    }% scalebox
    \caption{Measurement interval according to equation \ref{eqn:allan_variance_measurement}. The shaded region is the signal acquisition period.}
    \label{fig:allan_variance_definitions}
\end{figure}

Using this, the deviation over $N$ samples is defined as \cite{adev,psd_to_adev}
\begin{equation}
    \sigma_y^2(N,T,\tau) = \left\langle \frac{1}{N-1} \left(\sum _{k=0}^{N-1}\bar y_k^2(t)-\frac{1}{N}\left(\sum _{k=0}^{N-1} \bar y_k(t)\right)^2\right)\right\rangle
\end{equation}
The $\langle \; \rangle$ denotes the (infinite time) average over all measurands $y_k$ or, simply put, the expected value.

The Allan variance is a special case of this definition with zero dead-time ($\theta=0$) and only 2 samples:
\begin{align}
    \sigma_A^2(\tau) &= \sigma_A^2(N=2,T=\tau,\tau) \label{eqn:allan_coefficients}\\
    &= \left\langle \frac{\left(\bar y_{k+1} - \bar y_k \right)^2}{2} \right\rangle
\end{align}
It can be shown \cite{psd_to_adev} that \ref{eqn:adev_estimator} is indeed more useful than $\sigma_A^2(N\to\infty,T=\tau,\tau)$, because $\sigma_A^2(N=2,T=\tau,\tau)$ converges for processes that do not have a convergent $\sigma_A^2(N\to\infty,T=\tau,\tau)$.

In practice, no experiment can take an infinite number of samples, so typically the Allan variance must be estimated using a number of samples $m$:
\begin{equation}
    \sigma_A^2(\tau) \approx \frac1 m \sum_{k=1}^m \frac{\left(\bar y_{k+1} - \bar y_{k} \right)^2}{2} \label{eqn:adev_estimator}
\end{equation}
This estimation can lead to artifacts in the results as discussed later. In order to derive the Allan variance from a set of data points, the different values of $\tau$ are usually obtained by averaging over a number of samples as there is no dead time (by definition of the Allan variance).

Additionally, the Allan variance is mathematically related to the two-sided power spectral density $S_y(f)$ \cite{psd_to_adev}:
\begin{equation}
    \sigma_A^2(\tau) = 2 \int_0^\infty S_y(f) \frac{\sin^4\left( \pi f \tau \right)}{(\pi f \tau)^2}\,df \label{eqn:psd_to_adev}
\end{equation}

and therefore all processes that can be observed in the power spectral density plot can also be seen in the Allan deviation. The inverse transform however, is not always possible as shown by \citeauthor{inverse_adev} \cite{inverse_adev}.

Distinguishing different noise processes using the Allan deviation will be elaborated in the next section.

%TODO: Add Shot noise
\subsection{Identifying Noise in Allan Deviation Plots}
It was already mentioned by \citeauthor{adev} in \cite{adev} that types of noise, whose spectral density follows a power law
\begin{equation}
    S(f) = h_{\alpha} \cdot f^\alpha \label{eqn:power_law}
\end{equation}
can be easily identified in the Allan deviation plot. The constant $h_\alpha$ is called the power (intensity) coefficient. The most common types of noise encountered in experimental data and their representations can be found in table \ref{tab:adev_alpha}, which serves as a summary of this section. Since those types of noise are present in any measurement or electronic device, it warrants a further discussion to understand their root causes and ideas to minimize them. While not a type of noise, linear drift can also be easily identified in the Allan deviation plot. It is therefore included in table \ref{tab:adev_alpha} as well.
\begin{table}[ht]
    \centering
    \begin{tabular}{lcc}
        \toprule
        Amplitude noise type& Power-law coefficient $\alpha$& Allan variance $\sigma_A^2$\\
        \midrule
            White noise & $0$& $\frac 1 2 h_0 \tau^{-1}$ \cite{adev_noise_types}\\
            Flicker noise& $-1$& $2 \ln 2 \, h_{-1} \tau^0$ \cite{adev_noise_types}\\
            Random walk noise& $-2$& $\frac 3 2 \pi^2 h_{-2} \tau^{1}$ \cite{adev_noise_types}\\
            Burst noise& $0 \textrm{ and } -\!2$& $y_{rms}^2\frac{\bar \tau^2}{\tau^2} \left(4 e^{-\frac{\tau}{\bar \tau}} - e^{-\frac{2 \tau}{\bar \tau}} + 2 \frac{\tau}{\bar \tau} - 3 \right)$\\
            Drift & --& $\frac 1 2 D^2 \tau^2$ \cite{adev_drift}\\
        \bottomrule
    \end{tabular}
    \caption{Power law representations of different noise types using the Allan variance.}
    \label{tab:adev_alpha}
\end{table}

In order to arrive at a good understanding of the features seen in an Allan deviation plot, this section will provide the reader with examples of each type of noise and the corresponding time domain, power spectral density and Allan deviation plot. Since a complete overview is not available in current literature, all required mathematical descriptions and simulation tools will be discussed here. The simulations were done using Python and the source code is linked to in the discussions. The files are found in the online supplemental material found at \cite{supplemental_material}. Using these scripts, all the graphs shown can be recreated and explored further.

\subsubsection{White Noise}%
\label{sec:white_noise}
White noise is probably the most common type of noise found in measurement data. Johnson noise found in resistors, caused by the random fluctuation of the charge carriers, is one example of mostly white noise up to a bandwidth of \qty{100}{\MHz}, from where on quantum corrections are required \cite{nist_johnson_noise}. Amplifiers also tend to have a white noise spectrum at higher frequencies.

For the latter reason, white noise typically makes up for a considerable amount of noise in measurements, unless one works at very low frequencies. White noise is a series of uncorrelated random events and therefore characterised by a uniform power spectral density, which means there is the same power in a given bandwidth at all frequencies up to infinity. White noise therefore has infinite power (variance). In reality a measurement is always limited in bandwidth and hence the above property of a constant power spectral density only holds within that bandwidth. Those bandlimited samples of white noise thus have a finite variance.
Since white noise is so common, a few of its properties should be mentioned. One such property is that the variance $\sigma_{x+y}^2$ of two uncorrelated variables $x$ and $y$ adds as:
\begin{equation}
    \sigma_{x+y}^2  = \sigma_x^2 + \sigma_y^2 + \underbrace{2\,\mathrm{Cov}(x,y)}_{\text{uncorrelated}\, =\, 0}\ = \sigma_x^2 + \sigma_y^2 \label{eqn:adding_white_noise}
\end{equation}

This results in simple addition rules for variances from different sources, but it must be stressed here that this property is only valid for uncorrelated sources like white noise, although it is usually incorrectly applied to all measurements which unfortunately obscures rather than clarifies the uncertainties involved.

In order to demonstrate the effect of white noise in Allan deviation plots, it was simulated using the \textit{AllanTools} library written by \citeauthor{allantools} \cite{allantools}. The noise generator is based on the work of \citeauthor{noise_generation} \cite{noise_generation}. The full Python program code is published online \cite{supplemental_material} and found in \external{data/simulations/sim\_allan\_variance.py}. To allow better comparison, all noise densities are normalised to give an Allan deviation of $\sigma_A(\tau_0)=1$, with $\tau_0$ being the smallest time interval.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/white_noise_amplitude.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Time domain}
        \label{fig:white_noise_time}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/white_noise_psd.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Power spectral density}
        \label{fig:white_noise_psd}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/white_noise_adev.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Allan deviation}
        \label{fig:white_noise_adev}
    \end{subfigure}
    \caption{Different representations of white noise.}
    \label{fig:white_noise_simulated}
\end{figure}

Figure \ref{fig:white_noise_simulated} shows a sample of white noise in its three different forms. Figure \ref{fig:white_noise_time} is the time series representation from which the power spectral density was calculated and is shown in figure \ref{fig:white_noise_psd}. The dashed line shows the expectation value of the power spectral density and the Allan deviation.

From this simulation, several features can be observed. First of all, the power spectral density is flat and constant with $h_0 = 2$, which is in accordance with table \ref{tab:adev_alpha} and the normalisation mentioned earlier. Figure \ref{fig:white_noise_adev} shows the typical $\tau^{-\frac 1 2}$ dependence of white noise in the Allan deviation plot. This immediately explains, why filtering white noise scales with $\frac{1}{\sqrt{n}}$ with $n$ being the number of samples averaged.

\subsubsection{Burst Noise}%
\label{sec:theory_burst_noise}
Burst noise, popcorn noise, or sometimes referred to as random telegraph signal is a random bi-stable change in a signal and is caused by generation-recombination processes. This happens, for example, in semiconductors if there is a site that can trap an electron for a prolonged period of time and then randomly release it. Impurities causing lattice defects are discussed in this context \cite{kay2012operational,burst_noise_psd,popcorn_noise_orgin,technote_ti_popcorn_noise}. Such lattice defects can also be introduced by ion implantation during doping. Fortunately, this type of noise has become less prevalent in modern manufacturing processes, because the quality of the semiconductors has improved. But if a trap site is located very close to an important structure, for example a high precision Zener diode, its effect might be so strong that it can be clearly seen.

The discussion is split into two parts. First the power spectral density is calculated and then the Allan variance is calculated using that result.

The spectral density of burst noise caused by a single trap site was derived in \cite{burst_noise_wiener_khinchin} by \citeauthor{burst_noise_wiener_khinchin}. \citeauthor{burst_noise_wiener_khinchin} used the autocorrelation function of the burst noise signal and applied the Wiener-Khinchin (Wiener-Хи́нчин) theorem, which connects the autocorrelation function with the power spectral density. A more detailed derivation can be found in \cite{fundamentals_of_noise_processes}, in this paper the preconditions like stationarity of the process, are also discussed. The burst noise signal consists of two energy levels, called $0$ and $1$, split by $\Delta y$. Multiple burst noise signals can be superimposed in a real device. This would then result in multiple levels, but they can be treated separately. The measurement interval over an even number of transitions, so that one ends in the same state as the measurement has started, is the time $T$. The mean lifetime of the levels is called $\bar \tau_0$ and $\bar \tau_1$:
\begin{equation}
    \bar \tau_{0} \approx \frac 1 N \sum_{i}^N \tau_{0,i} \qquad \bar \tau_{1} \approx \frac 1 N \sum_{i}^N \tau_{1,i}
\end{equation}

Figure \ref{fig:burst_noise} shows a burst noise signal along with the definitions above.
\begin{figure}[hb]
    \centering
    %\scalebox{1}{%
        \import{figures/}{burst_noise.tex}
    %} % scalebox
    \caption{A random burst noise signal.}
    \label{fig:burst_noise}
\end{figure}

Using these definitions, one can then derive \cite{burst_noise_wiener_khinchin}:
\begin{align}
    R_{xx}(T) &= \Delta y^2 \cdot \frac{\bar \tau_1 \bar \tau_0 e^{-\left(\frac{1}{\bar \tau_1}+\frac{1}{\bar \tau_0}\right)T}}{\left(\bar \tau_1 + \bar \tau_0\right)^2} \quad \text{and} \label{eqn:burst_noise_correlation}\\
    S(\omega) &= 4 R_{xx}(0) \frac{\frac{1}{\bar \tau_1} + \frac{1}{\bar \tau_0}}{\left(\frac{1}{\bar \tau_1} + \frac{1}{\bar \tau_0}\right)^2 + \omega^2} \qquad \omega > 0 . \label{eqn:burst_noise_psd}
\end{align}
Note, that the power spectral density is the one-sided version, hence an additional factor of $2$ is included. The constant term was omitted here and can usually be neglected, because it is not relevant for calculating the power spectral density as it only contributes a single peak at $\omega=0$. Using the following definitions of the average time constant and the duty cycle
\begin{align}
    \frac{1}{\bar \tau} &= \frac{1}{\bar \tau_1} + \frac{1}{\bar \tau_0} \quad \mathrm{and} \label{eqn:definition_bar_tau}\\
    D_i &= \frac{\bar \tau_i}{\bar \tau_1 + \bar \tau_0} \quad i \in \{0 ; 1\}
\end{align}

equations \ref{eqn:burst_noise_correlation} and \ref{eqn:burst_noise_psd} can be rewritten to give a more intuitive form.
\begin{align}
    R_{xx}(T) &= \Delta y^2 D_1 D_0 \, e^{-\left(\frac{1}{\bar \tau_1}+\frac{1}{\bar \tau_0}\right)T}\\
    S(\omega) &= 4 R_{xx}(0) \frac{\bar \tau}{1 + \omega^2 \bar \tau^2} \label{eqn:burst_noise_lorentzian}
\end{align}

The special case $\bar \tau_0 = \bar \tau_1$ with $D_i=\frac 1 2$ is the previously mentioned case of random telegraph noise.

$R_{xx}(0)$ can be identified as the mean squared value of $y$:
\begin{equation}
    y_{rms} = \sqrt{R_{xx}(0)} \,.
\end{equation}

Equation \ref{eqn:burst_noise_lorentzian} is a Lorentzian function and from this, it can be easily seen that a single trap site has a power spectral density that is proportional to $\frac{1}{f^2}$ at high frequencies and is flat at low frequencies.

With the spectral density in hand, it is now possible to calculate the Allan variance as it was done by \citeauthor{allen_dev_flicker} in \cite{allen_dev_flicker} for the classic example of random telegraph noise where $\bar \tau_1 = \bar \tau_0$. Do note that table I given by \citeauthor{allen_dev_flicker} shows the total number of events instead of the instantaneous number of events typically given. Hence their notation must be multiplied by $\frac{1}{\tau^2}$ (or $\frac{1}{T^2}$ in their notation). For the generic case with $\bar \tau_1$, $\bar \tau_0$ and the definition of $\bar \tau$ given in equation \ref{eqn:definition_bar_tau} one finds for the Allan variance of burst noise:
\begin{equation}
    \sigma^2_A(\tau) = R_{xx}(0) \frac{\bar \tau^2}{\tau^2} \left(4 e^{-\frac{\tau}{\bar \tau}} - e^{-\frac{2 \tau}{\bar \tau}} + 2 \frac{\tau}{\bar \tau} - 3 \right) \label{eqn:burst_noise_avar}
\end{equation}

Having arrived at equations \ref{eqn:burst_noise_lorentzian} and \ref{eqn:burst_noise_avar} of the power spectral density and Allan variance, it it now possible to model it. For this purpose, parts of the Python library \textit{qtt} \cite{qtt} was used. This algorithm written by \citeauthor{qtt} implements continuous-time Markov chains to simulate the burst noise signal. The result can be see in figure \ref{fig:burst_noise_simulated}. For these simulations one of the time constants, namely the lifetime of the lower state $\bar \tau_0$ was held constant, while the lifetime of the upper state was varied to show the effect of different $\bar \tau$. By looking at the time domain in figure \ref{fig:burst_noise_time} it can be seen that the maximum average number of state changes can be observed, when $\bar \tau_1 = \bar \tau_0$. If $\bar \tau_1 > \bar \tau_0$ the system will favour the upper, while if $\bar \tau_1 < \bar \tau_0$ it will favour the lower state instead. This explains why the noise is strongest for random telegraph noise when $\bar \tau_1 = \bar \tau_0$, which can also be seen in the power spectral density plot in figure \ref{fig:burst_noise_psd}. Looking at the Allan deviation in figure \ref{fig:burst_noise_adev} confirms this, but also shows another interesting implication as it shows an obvious maximum. If the application allows a choice over the sampling interval $\tau$, the effect of the burst noise can be mitigated by staying well clear of the maximum.

The small deviation from the analytical solution in figure \ref{fig:burst_noise_adev} suggesting an upwards trend at large $\tau$ is a typical so-called end-of-data error. As it was discussed above, the Allan deviation can only be estimated given a limited number of samples using equation \ref{eqn:adev_estimator} and going to longer $\tau$ means there are fewer samples to average over.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.8\linewidth}
        \centering
        \scalebox{1}{%
            \input{images/burst_noise_amplitude.pgf}% data/simulations/sim_burst_noise.py
        } % scalebox
        \caption{Time domain}
        \label{fig:burst_noise_time}
    \end{subfigure}
    \begin{subfigure}{0.8\linewidth}
        \centering
        \scalebox{1}{%
            \input{images/burst_noise_psd.pgf}% data/simulations/sim_burst_noise.py
        } % scalebox
        \caption{Power spectral density}
        \label{fig:burst_noise_psd}
    \end{subfigure}
    \begin{subfigure}{0.8\linewidth}
        \centering
        \scalebox{1}{%
            \input{images/burst_noise_adev.pgf}% data/simulations/sim_burst_noise.py
        } % scalebox
        \caption{Allan deviation}
        \label{fig:burst_noise_adev}
    \end{subfigure}
    \caption{Different representations of burst noise for different $\bar \tau_1$ and fixed $\bar \tau_0 = \qty{1}{\s}$.}
    \label{fig:burst_noise_simulated}
\end{figure}

The burst noise equations can be used to gain further insight into other types of noise. The first one is Shot noise, which is commonly found in photodetectors and lasers. Here, electrons or photons are created at discrete intervals resulting in an instantaneous signal. This means that the lifetime of the upper level is very short in comparison to the lower level ($\tau_1 \ll \tau_0$) equation \ref{eqn:burst_noise_psd} becomes:
\begin{align}
    S_{Shot}(\omega) = S_{\tau_1 \ll \tau_0}(\omega) &= 4 \Delta y^2 \frac{\tau_1}{\tau_0} \frac{\frac{1}{\bar \tau_1}}{\left(\frac{1}{\bar \tau_1}\right)^2 + \omega^2}\nonumber\\
    &= 4 \Delta y^2 \frac{1}{\tau_0} \frac{1}{\frac{1}{\tau_1^2}+\omega^2}\\
    \overset{\omega \ll 1/\tau_0}&{\approx} 4 \Delta y^2 \frac{\tau_1^2}{\tau_0} = \text{const.}
\end{align}

Typically, a very large number of such events happen. When not counting single events, but rather a stream, the relation $\omega \ll 1/\tau_0$ is valid and hence the result is a white spectrum as $S_{Shot}(\omega)$ is constant with respect to $\omega$ --- just as observed in photodetectors and lasers.

The other interesting occurrence is a case where many trap sites with different time constants are contributing to the noise. This can change the shape of the spectrum from $f^{-2}$ to $f^{-1}$ and is discussed in the next section.

\clearpage
\subsubsection{Flicker Noise}%
\label{sec:flicker_noise}
Flicker noise is also called $\frac 1 f$-noise and it can be observed in many naturally occurring phenomena. Its origin is not clear, even though there have been many explanations. An overview can be found in \cite{flicker_noise_overview, flicker_noise_overview2, origins_1_f_noise}. This section concentrates on flicker noise in electronic devices. In thick-film resistors, for example, it was shown to extend over at least 6 decades without any visible flattening \cite{1_f_noise_thick_film}. In transistors, flicker noise is caused by the existence of generation-recombination noise or burst noise discussed in the previous section \cite{origins_1_f_noise}. If there are many uncorrelated trap sites which contribute to the total noise, the envelope of the noise spectral density changes from $\frac{1}{f^2}$ to $\frac{1}{f^1}$ as shown in figure \ref{fig:flicker_noise_evelope}
\begin{figure}[hb]
    \centering
    \input{images/flicker_noise_envelope.pgf}% data/simulations/sim_flicker_noise_envelope.py
    \caption{Multiple overlapping Lorentzian noise sources forming a $\frac 1 f$-like shape.}
    \label{fig:flicker_noise_evelope}
\end{figure}

Given that no trap site can store an electron indefinitely, the number of trap sites $N$ with a certain time constant $\frac 1 2 \bar \tau = \bar \tau_0 = \bar \tau_1$ must decline when going to longer time scales. Assuming $N$ is inversely proportional to the time constant $\bar \tau$
\begin{equation}
    N(\tau) \propto \frac{1}{\bar \tau}\,, \label{eqn:flicker_noise_weight_function}
\end{equation}
which can be motivated if the trapping process is thermally activated \cite{1_f_noise_motivation} and using equation \ref{eqn:burst_noise_lorentzian} from the previous section, multiplying the weight function \ref{eqn:flicker_noise_weight_function} and integrating over all possible storage times gives:
\begin{align}
    S(\omega) &= \lim_{t \to \infty} \int_0^t N(\bar \tau) \, 4 R_{xx}(0) \frac{\bar \tau}{1 + \omega^2 \bar \tau^2} \, d\bar\tau \nonumber\\
    \overset{\bar \tau_0 = \bar \tau_1}&{=} 4 R_{xx}(0)\, C_N \lim_{t \to \infty} \int_0^t \frac{1}{1 + \omega^2 \bar\tau^2} \, d\bar\tau \nonumber\\
    &= \frac{4 R_{xx}(0)\, C_N}{\omega} \lim_{t \to \infty}  \arctan{\bar\tau \omega} \Big|_{\bar\tau=0}^t \nonumber\\
    &= \frac{4 R_{xx}(0)\, C_N}{\omega} \cdot \frac{\pi}{2} \nonumber\\
    &= \frac{2 \pi R_{xx}(0)\, C_N}{\omega}\\
    S(f) &= h_{-1} f^{-1}
\end{align}

$C_N$ is the proportionality constant of \ref{eqn:flicker_noise_weight_function} and $h_{-1}$ is the power coefficient introduced in \ref{eqn:power_law}. This shows that for a large number of distributed trap sites, a noise spectrum of $f^{-1}$ is found.

Using equation \ref{eqn:psd_to_adev}, the Allan variance can be calculated from the power spectral density:
\begin{align}
    \sigma_A^2(\tau) &= 2 h_{-1} \int_0^\infty \frac{1}{f} \frac{\sin^4\left( \pi f \tau \right)}{(\pi f \tau)^2}\,df \nonumber\\
    &=2 \ln 2 \, h_{-1}
\end{align}

Again, using the \textit{AllanTools} library \cite{allantools}, flicker noise was simulated to give an impression of its properties.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/flicker_noise_amplitude.pgf}% data/simulations/sim_allan_variance.py
            %FIXME: Correct height
        } % scalebox
        \caption{Time domain}
        \label{fig:flicker_noise_time}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/flicker_noise_psd.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Power spectral density}
        \label{fig:flicker_noise_psd}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/flicker_noise_adev.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Allan deviation}
        \label{fig:flicker_noise_adev}
    \end{subfigure}
    \caption{Different representations of flicker noise.}
    \label{fig:flicker_noise_simulated}
\end{figure}

While it is not immediately evident from the power spectral density, the Allan deviation plot explains very well, why additional filtering does not affect flicker noise. No matter how long the integration time, the variance will still be the same.

The small wiggles at longer $\tau$ are typical end-of-data errors caused by spectral leakage, because there are insufficient samples to average over \cite{adev_long_tau}. As it was discussed above, the Allan deviation can only be estimated using equation \ref{eqn:adev_estimator} given a limited number of samples. Therefore, at $\frac{\tau}{2}$ there are only $2$ samples left, so there is no averaging possible to improve the estimate of the Allan deviation, which causes the oscillations at low frequencies or large $\tau$.

As a last remark, a commonly used definition in combination with flicker noise is the corner frequency $f_c$. The corner frequency appears in situations where there is both flicker and white noise present. It is the crossover point in frequency, where the flicker noise is equal compared to the white noise.
\begin{equation}
    f_c = \frac{h_{-1}}{h_0} \label{eqn:corner_frequency}
\end{equation}
It can be graphically extracted from the power spectral density plot by drawing a line trough the flicker noise and the white noise and finding the intersection. This can be seen in figure \ref{fig:adev_example_psd} on page \pageref{fig:adev_example_psd}. The corner frequency can be found where the horizontal dashed blue and green line meet.

\subsubsection{Random Walk}%
\label{sec:random_walk}
Random walk noise can be attributed to environmental factors such as temperature \cite{random_walk_fm} and diffusion processes, the latter contributing to the ageing effect seen in semiconductors.
It is a process, where in each time step the change is randomly determined to be either a positve or negative step with equal probability and a fixed step size. Its mean is
\begin{equation}
    \langle y_n \rangle = \langle e_1 + e_2 + \dots e_n \rangle = \underbrace{\langle e_1 \rangle}_{=\,0} + \langle e_2 \rangle + \dots + \langle e_n \rangle = 0 \, ,
\end{equation}
but its variance
\begin{equation}
    \sigma_y^2 = \langle y_n^2 \rangle - \underbrace{\langle y_n \rangle}_{=\,0} = \sigma_{e_1}^2 + \sigma_{e_2}^2 + \dots \sigma_{e_n}^2 = n \sigma_e^2
\end{equation}
goes with $n$ (or $t$). It therefore is not a stationary process as can also be seen in figure \ref{fig:random_walk_adev}.

The power spectral density can be calculated \cite{psd_to_adev,noise_generation} to be
\begin{equation}
    S(f) = h_{-2} \frac{1}{f^2}
\end{equation}
and the Allan deviation can again be calculated from the spectral density
\begin{align}
    \sigma_A^2(\tau) &= 2 h_{-2} \int_0^\infty \frac{1}{f^2} \frac{\sin^4\left( \pi f \tau \right)}{(\pi f \tau)^2}\,df \nonumber\\
    &=\frac{2}{3} \pi^2 h_{-2}\, \tau
\end{align}

The \textit{AllanTools} library \cite{allantools} can then be used to simulate the random walk.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/random_walk_noise_amplitude.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Time domain}
        \label{fig:random_walk_time}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/random_walk_noise_psd.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Power spectral density}
        \label{fig:random_walk_psd}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/random_walk_noise_adev.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Allan deviation}
        \label{fig:random_walk_adev}
    \end{subfigure}
    \caption{Different representations of random walk noise.}
    \label{fig:random_walk_noise_simulated}
\end{figure}


\clearpage
\subsubsection{Drift}
Finally, the last feature of the Allan deviation plot that needs to be discussed is drift. Drift happens at very long time scales and describes a linear dependence of the measurand on time. This is also part of the ageing effect. \citeauthor{adev_drift} discussed the effect of drift \cite{adev_drift} on the Allan variance and found the following relationship:
\begin{align}
    \sigma_A^2(\tau) = \frac{D^2}{2} \tau^2
\end{align}
with slope of the drift $D$.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/drift_amplitude.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Time domain}
        \label{fig:drift_time}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/drift_adev.pgf}% data/simulations/sim_allan_variance.py
        } % scalebox
        \caption{Allan deviation}
        \label{fig:drift_adev}
    \end{subfigure}
    \caption{Different representations of linear drift.}
    \label{fig:drift_noise_simulated}
\end{figure}

\subsubsection{Dead Time}%
\label{sec:dead_time}
The coefficients given in the previous examples were derived using the assumption that all samples in a measurement are continuous with a dead time $\theta = 0$. Unfortunately, measurements sometimes have a dead time that is non negligible. This problem was extensively discussed by \citeauthor{psd_to_adev} \cite{psd_to_adev}. \citeauthor{adev_frequency_counter} even developed special models to account for the algorithms of modern frequency counters \cite{adev_frequency_counter}. While some frequency counters support gapless measurements, the situation is entirely different for digitizers and digital multimeters. Several settings commonly used affect the dead time, which can be considerable. It is therefore important to discuss typical measurement settings for voltmeters to estimate the errors that arise from those settings. The focus of this discussion lies on the dead time introduced by digital multimeters, but the application is not limited to this field.

The most commonly used settings that affect the dead time of a voltmeter are autozeroing and line synchronization. Autozeroing is done by adding additional measurements to the normal input integration cycle. To correct for the zero offset drift a zero measurement is added where the ADC is switched to the low terminal. Additionally, some devices add a reading of the reference voltage to correct for gain errors. The implementation details and type of measurements are manufacturer dependent.

The other setting, which can be enabled in voltmeters, is the line synchronization to increase the noise rejection of the instrument. This setting synchronizes the start of a measurement to the zero crossing of the power line. Depending on the instrument, this might cause a delay of one power line cycle (PLC) after each measurement if the instrument is not capable of processing the previous measurement while at the same time recording another one.

A simple measurement with dead time is shown in figure \ref{fig:allan_variance_definitions} on page \pageref{fig:allan_variance_definitions}. That model assumes that the dead time is constant and is always added after the actual integration time $\tau$. This is rarely true for real measurement data as many devices and even ADCs use internal averaging and autozeroing to produce a measurement. The actual dead time is therefore spread over the whole measurement and not limited to the end of the measurement. An example is the Keysight \device{3458A} DMM that automatically switches to averaging when selecting integration times greater than \qty{10}{\plc}. The reason is simple as for longer integration times, more and more flicker noise starts contributing to the measurement. The measurement is therefore split into single measurements of \qty{10}{\plc} and, using autozeroing, the flicker noise is suppressed. This is discussed in more detail as an example in section \ref{sec:autozero}. The mathematical problem of a distributed dead time was already noted by \citeauthor{adev_noise_types} \cite{adev_noise_types} and it is distinctively different from the calculations made by \citeauthor{psd_to_adev} \cite{psd_to_adev} for a single dead time at the end of the measurement. The exact mathematical treatment is complex and is beyond the scope of this work, especially considering that autozeroing does a lot more than just adding dead time at the end of the measurement. Fortunately, using a few assumptions, the problem can be greatly simplified.

An interesting observation can be made for white noise. Since it is uncorrelated, it makes no difference whether it is sampled in full, or only partially, and therefore the Allan deviation for a white noise process with or without dead time is the same:
\begin{equation}
    \sigma^2(N,T, \tau) = \sigma^2(N=2,T=\tau, \tau) = \sigma_A^2(\tau) \frac 1 2 h_0 \tau^{-1}
\end{equation}

Consequently, if the dead time is added at a frequency high enough, so that the input amplifier output is dominated by white noise, the dead time will have no influence on the Allan variance.

Finally, \citeauthor{psd_to_adev} \cite{psd_to_adev} notes that for measurement durations or averaging times $T \gg T_0$, the Allan variance with respect to $T$ shows an asymptotic behaviour of $\sigma_A^2(T) \to \sigma_A^2(\tau)$.

\subsection{Example}%
\label{sec:noise_example}
Using the results from the previous sections, it is possible to simulate a typical measurement sample containing white noise, flicker noise and random walk behaviour. The simulation was written in Python using the \textit{AllanTools} library \cite{allantools} to generate the time domain data, which was then converted to a power spectrum using the algorithm developed by \citeauthor{welch} \cite{welch} to estimate the power spectral density. The Allan deviation was calculated using the \textit{AllanTools} library. The full Python source is available at \cite{supplemental_material} and found in \external{data/simulations/sim\_allan\_variance\_example.py}. The time domain data shown here were downsampled from $2^{25}$ data points to \num{2000} points for faster plotting using the Largest-Triangle-Three-Buckets (LTTB) algorithm created by \citeauthor{lttb} \cite{lttb} and also available as a Python library. The downsampling algorithm chosen is optimal for this application because it aims to visually keep the result the same by favouring parts of the data where there is more dynamics. The only difference noticeable to the author is that the edges of the white noise plot are a slightly rougher. The full data set can be obtained using the source code given above if one desires. The power spectrum and the Allan deviation were always calculated from the full dataset. The data of the power spectrum were additionally binned to be evenly spaced on a logarithmic scale. This considerably reduced the high frequency noise and made the plot easier to read while not negatively impacting its shape.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/example_white_noise_amplitude.pgf}% data/simulations/sim_allan_variance_example.py
        } % scalebox
        \caption{White noise}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/example_flicker_noise_amplitude.pgf}% data/simulations/sim_allan_variance_example.py
        } % scalebox
        \caption{Flicker noise}
    \end{subfigure}
    \begin{subfigure}{0.32\linewidth}
        \centering
        \scalebox{0.75}{%
            \input{images/example_random_walk_noise_amplitude.pgf}% data/simulations/sim_allan_variance_example.py
        } % scalebox
        \caption{Random walk}
    \end{subfigure}
    \caption{Three separate noise components that were summed together to simulate a typical noise source.}
    \label{fig:adev_example_noise_types}
\end{figure}

The three time series shown in figure \ref{fig:adev_example_noise_types} were sequentially generated using a fixed seed for the random number generator to ensure repeatability as long as the order of creation is kept the same. For generating the noise, the algorithm presented by \citeauthor{noise_generation} \cite{noise_generation}, implemented in the \textit{AllanTools} Python library was used. The noise strength parameters were deliberately chosen in such a way that both the white noise and the random walk part have more noise power than the flicker noise. This allows to distinguish them in the following plots at both extremes of the frequency scale and time scale. Finally, the three types of noise data were summed together to give the combined signal, which is shown in figure \ref{fig:adev_example_time}, again downsampled using LTTB. The summed series clearly shows the white noise content and it is possible to deduce either flicker or random walk noise, but it is highly obscured due to the amount of white noise. Using only the time domain plot makes it very hard to clearly distinguish the type of noise present, let alone estimate the individual noise power of the three components. Therefore, a different analysis tool is called for.
\begin{figure}[ht]
    \centering
    \input{images/example_combined_noise_amplitude.pgf}% data/simulations/sim_allan_variance_example.py
    \caption{A simulated time series containing white noise, flicker noise and random walk behaviour.}
    \label{fig:adev_example_time}
\end{figure}

A common approach to identify noise sources is the power spectrum. It is easily accessible, even in real-time using spectrum analysers and, utilizing the computational power of modern computers, large time-domain data sets can be converted using the fast Fourier transform, making this the method of choice in the lab. The power spectral density  of figure \ref{fig:adev_example_time} is shown in figure \ref{fig:adev_example_psd}. It allows to clearly separate the white noise part from the other $f^{\alpha}$ components. The dashed lines representing the individual components were plotted using the $h_\alpha$ values calculated from the input parameters of the simulation. The noise spectral density $h_0$ of the white noise signal can be easily extracted even by hand without resorting to a fit. This yields $h_{0} = \qty{2e-3}{\per \Hz} $. $h_{-1}$ and $h_{-2}$ can be extracted as well using a fit to
\begin{equation}
    S(f) = \sum_{\alpha = -2}^0 h_\alpha f^\alpha \, .
\end{equation}
The noise corner frequency $f_c$ can either be calculated from $h_0$ and $h_{-1}$ using equation \ref{eqn:corner_frequency} or determined graphically by constructing a tangent with a slope of $-1$ to the spectral density. From the intersection of the blue $h_0$ line and the green $h_{-1}$ line the corner frequency is found to be $f_c \approx \qty{1.8}{\kHz}$.
\begin{figure}[hb]
    \centering
    \input{images/example_combined_noise_psd.pgf}% data/simulations/sim_allan_variance_example.py
    \caption{A simulated power spectrum containing white noise, flicker noise and random walk behaviour.}
    \label{fig:adev_example_psd}
\end{figure}

To get an even better representation of the individual noise contributions, the Allan variance or Allan deviation can be used. The Allan deviation plot shown in figure \ref{fig:adev_example_adev} gives very clean results and all noise components can be clearly identified. The individual components were plotted using dashed lines as well.
\begin{figure}[ht]
    \centering
    \input{images/example_combined_noise_adev.pgf}% data/simulations/sim_allan_variance_example.py
    \caption{A simulated Allan deviation containing white noise, flicker noise and random walk behaviour.}
    \label{fig:adev_example_adev}
\end{figure}

The Allan variance was calculated using the overlapping Allan variance algorithm \cite{oadev_definition} and only Allan deviation values for frequency values of $(1, 2, 4)$ per decade were plotted. The overlapping Allan variance gives a better confidence at longer intervals or lower frequencies, allowing to identify very low frequency noise like the random walk shown here. Reference \cite{oadev_definition} also gives a very good comparison of other algorithms to identify even more noise types in data sets like phase noise. Plotting only three values per decade improves the clarity of the plot, because at longer $\tau$s, even though the overlapping Allan variance is used, some oscillations inevitably show up. Using fewer values of $\tau$ causes less distractions in this case.
From the figure \ref{fig:adev_example_adev}, the Allan deviation of the flicker noise can be estimated from the flat minimum to be around $2.3$ or $\sqrt{5}$. Using table \ref{tab:adev_alpha} the Allan variance can be converted to
\begin{equation*}
    h_{-1} = \frac{5}{2 \ln 2} \approx 3.6
\end{equation*}

Using the previously found $h_0$, this corner frequency is calculated using equation \ref{eqn:corner_frequency} to be:
\begin{equation*}
    f_c = \frac{5}{\qty{2e-3}{\per \Hz} \cdot 2 \ln 2} \approx \qty{1.8}{\kHz}
\end{equation*}

This is obviously the same result as the one from the geometric approach above.

This concludes the examples section for different noise types. The reader should now be able to identify different types of noise in measurement data and have learnt to appreciate the value that the Allen variance brings to the table. An example was presented that applied all techniques shown in this section to extract information about the noise sources in a dataset. Additionally, Python source code is provided to further explore the topic.

\clearpage
\section{Autozeroing}%
\label{sec:autozero}
Autozeroing (AZ), sometimes called zero-drift or dynamic offset compensation, is such an important concept that it must be discussed in its own right. The need for autozeroing comes from the typical behaviour of amplifiers. Every amplifier has some offset, be it small or large, and especially at high gains, this offset becomes a problem for high precision measurements. To make matters worse, this offset is not stable over time and drifts with both time and temperature. It can therefore not be calibrated out once, it must be permanently adjusted during operation, depending on environmental conditions. This procedure is called autozeroing.

There are many different ways to implement autozeroing and regarding operational amplifiers a good overview can be found in \cite{horowitz1989}. As an example, the autozero cycle for the Keithley \device{Model 2002} and the Keysight \device{3458A} Multimeter is shown in figure \ref{fig:dmm_autozero_comparison}. Keithley uses a more complex and slower algorithm, while HP implemented a simpler but faster algorithm. The most simple (digital) approach is to regularly switch the input from the signal to zero, take a reading, then subtract this reading from all subsequent readings until a new zero reading is taken. An alternative approach adds another measurement of the reference voltage to apply a gain correction as well. This is done by the Keithley \device{Model 2002} and works very well to suppress gain drift in the input amplifier due to temperature changes but increases the time between samples by another \qty{50}{\percent}. The Keysight/HP \device{3458A} in the other hand calculates those gain corrections only during the manual auto-calibration (ACAL) routine to maintain a higher throughput.
\begin{figure}[hb]
    \centering
    %\resizebox {0.8\textwidth} {!} {
        \import{figures/}{dmm_autozero.tex}
    %} % resizebox
    \caption{Auto-zero phases of the Keysight \device{3458A} and Keithley \device{Model 2002}.}
    \label{fig:dmm_autozero_comparison}
\end{figure}

\subsection{Offset Nulling}%
\label{sec:autozero_offset_nulling}
Offset nulling is the most basic approach to autozeroing. It aims to remove the offset drift of an amplifier. Especially at high gains, the offset, which is multiplied by the gain, can be substantial. In order to explain how offset nulling works and how it shapes the spectrum, it is best to discuss it based on an example. While this technique can also be found in many integrated circuits, it is more noticeable in DMMs, because it is a switchable option. Therefore, the example data set simulated is based on the parameters of the aforementioned Keysight \device{3458A} multimeter. The corner frequency and the white noise floor is modeled after the \qty{10}{\V} range of the \device{3458A} \cite{3458A_noise_floor, sampling_with_3458A} with the values given below. Do note that both references \cite{3458A_noise_floor, sampling_with_3458A} contain a typographical error. The corner frequency of the noise floor is erroneously given as \qty{0.5}{\Hz}, but should be \qty{1.5}{\Hz}. This can be seen in figure 2.35 in \citep[p. 116]{sampling_with_3458A}, where the noise spectral density is plotted and it was also confirmed with the author \cite{lapuh_email_corner_frequency}. The data used in this section is generated using the Python \textit{AllanTools} library \cite{allantools} and the simulation source code can be found in \external{data/simulations/sim\_autozero.py} as part of the online supplemental material \cite{supplemental_material}.
\begin{figure}[ht]
    \centering
    \scalebox{1}{%
        \import{figures/}{offset_nulling_definitions.tex}
    } % scalebox
    \caption{Integration sequences of the offset nulling algorithm. Solid lines denote sampled data. Red is the input signal, green is the zero reading and blue is the dead time required for switching inputs.}
    \label{fig:dmm_autozer_offset_nulling}
\end{figure}

For this simulation, a noise-free and arbitrarily chosen \qty{10}{\V} input is assumed to be sampled by the device at a sampling rate of \qty{10}{\plc} at \qty{50}{\Hz}, the same rate discussed previously on page \pageref{sec:dead_time}. As it will be shown, the actual mean value of the input signal has no bearing on the outcome of the calculation when considering offset nulling, but its value must be considered for other types of autozeroing as discussed in section \ref{sec:autozero_gain} and is included here only for the sake of completeness.

Figure \ref{fig:dmm_autozer_offset_nulling} shows the individual sequences of the offset nulling algorithm. First, the source is sampled for $\tau_s = \qty{10}{\plc}$, then the input is switched to the LO terminal. While this operation is very fast and takes less than \qty{1}{\ms} \cite{article_3458A_input_impedance}, if the instrument is synchronized to the line frequency the zero measurement will nonetheless be delayed until the next zero crossing, hence the dead-time $\theta = \qty{1}{\plc}$. Finally, the zero reference is measured for another $\tau_r = \qty{10}{\plc}$ and then the instrument switches back to the HI terminal.

The data is simulated in the following way: First, two sets of noise data are generated, a white noise spectrum with a noise spectral density of \qty[power-half-as-sqrt, per-mode=symbol]{165}{\nV \Hz\tothe{-0.5}} and a flicker noise spectrum with an intensity scaled to result in a final spectrum with a corner frequency of \qty{1.5}{\Hz}. The required flicker noise intensity is calculated using equation \ref{eqn:corner_frequency}. To get a good low frequency estimate, $2^{20} \approx 10^{6}$ values were generated. Finally, the two noise data sets are summed with the noise-free input source to give the final result. Other effects, such as power-line hum are neglected in this simple simulation because it would needlessly overcomplicate the example and limit the educational value. The same goes for higher order random-walk $f^{-2}$ noise components, which can be introduced by temperature fluctuations and other environmental effects and would be present in a real measurement.
\begin{figure}[ht]
    \centering
    \input{images/autozero_raw_amplitude.pgf}% data/simulations/sim_autozero.py
    \caption{Time series data with white noise and flicker noise.}
    \label{fig:autozero_raw_time}
\end{figure}

The time domain plot of the simulation is shown in figure \ref{fig:autozero_raw_time}. The white noise component is clearly visible, while the $f^{-1}$ flicker noise can be recognized, but its strength can hardly be estimated. It was already shown in section \ref{sec:noise_example} that different types of noise have different frequency components and those can be distinguished in the frequency domain, which leads to the next approach.
\begin{figure}[hb]
    \centering
    \input{images/autozero_raw_psd.pgf}% data/simulations/sim_autozero.py
    \caption{Simulated power spectrum of a Keysight \device{3458A} containing white noise and flicker noise. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_raw_psd}
\end{figure}

The noise power spectral density shown in figure \ref{fig:autozero_raw_psd} is calculated from the time series given in figure \ref{fig:autozero_raw_time} and confirms the flicker and white noise content. The theoretical white noise floor is shown as a horizontal dashed blue line and the flicker noise as a dashed green line. The \qty{1.5}{\Hz} corner frequency, which is defined as the intersection between the $f^{-1}$ noise and the white noise floor can be easily identified using those lines. It is evident that the \qty{5}{\Hz} sampling frequency with a \qty{2.5}{\Hz} bandwidth does not allow the spectral density to fully settle to the white noise floor.

From the power spectral density is can be seen that higher frequencies have a significantly lower noise spectral density than low frequencies. It is therefore most beneficial to do measurements at higher frequencies. To discuss the optimal measurement interval, the Allan deviation is an excellent tool.
\begin{figure}[ht]
    \centering
    \input{images/autozero_raw_adev.pgf}% data/simulations/sim_autozero.py
    \caption{Simulated Allan deviation of the input amplifier of a Keysight \device{3458A} containing white noise and flicker noise. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_raw_adev}
\end{figure}

The Allan deviation it plotted in figure \ref{fig:autozero_raw_adev} and shows two distinct regions. Short $\tau$ display an asymptotic behaviour towards white noise with a $\tau^{−0.5}$ dependence and at longer $\tau$ the constant flicker noise region can be identified. At very long $\tau$ typical end-of-data oscillations can be seen, which are the result of the limited confidence of the Allan deviation estimator as previously discussed and can therefore be safely ignored. The Allan deviation clearly demonstrates the performance of the device at longer integration times and it is obvious that beyond an integration time of about \qty{1}{\second} or \qty{50}{\plc} no additional information can be extracted from the measurement and the variance is constant. This leads to the need for autozeroing to remove the flicker noise. It can be shown \cite{autozero_with_dead_time} that subtracting a reference measurement from the actual measurement data removes all correlated effects. Since flicker noise is autocorrelated, it can be removed by subtracting a zero measurement.

To demonstrate autozeroing, two cases will be discussed. Going back to figure \ref{fig:dmm_autozer_offset_nulling} it can be seen that between switching inputs, a dead time $\theta$ is added. For a first discussion, this dead time is neglected and then the effect of adding a dead time is discussed in a next step.

Using figure \ref{fig:autozero_raw_adev} it was shown that integrating over flicker noise, does not reduce the variance. In order to have as little flicker noise content in the final measurement values as possible, it is clear that the autozeroing should be done as fast as feasible to keep the flicker noise content out. This allows to calculate the expected variance of the autozeroed measurement. The noise of the input measurement $x$ and the reference measurement $y$ are the same, because in this model the only noise source comes from the input amplifier, as the input signal is assumed to be noise-free. The zero level is, by definition, noise-free. As discussed above, the autozero interval is chosen, in such a way that its variance is dominated by white noise. The variance $\sigma^2$ of the combined measurement of $x-y$ can then be calculated using equation \ref{eqn:adding_white_noise}:
\begin{equation}
    \sigma_{x-y}^2 = \sigma_x^2 + \sigma_y^2 \label{eqn:autozeroing}
\end{equation}

By subtracting the zero reading the amplifier noise is effectively added twice to the final result, once for the input measurement and once for the zero measurement. Additional noise from the input signal noise would simply be added to this as it is uncorrelated as well.

Do note, that the number of samples is now half the number before applying autozeroing. This leads to an interesting effect. Taking for example a data set containing only white noise with a variance $\sigma^2$ and removing half the samples obviously does not change the variance as white noise is not correlated, but subtracting the samples is effectively decimating the data set and since the sampling rate is halved, the Nyquist band is halved as well. Unfortunately the input noise bandwidth stays the same. The second Nyquist band is then folded back into the first, thus doubling the noise power density.

To conclude, it is expected that the variance doubles and the power spectral density quadruples. These considerations can be compared to the simulated data. In order to apply the autozeroing algorithm to the simulated data set, the constant noise-free part, namely the \qty{10}{\V} signal, of the noisy input signal was nulled for every odd value and the residual noise was subtracted from the signal value. The result in the time domain is shown in figure \ref{fig:autozero_time}.
\begin{figure}[hb]
    \centering
    \input{images/autozero_amplitude.pgf}% data/simulations/sim_autozero.py
    \caption{Simulated measurement with autozeroing applied.}
    \label{fig:autozero_time}
\end{figure}

When comparing figure \ref{fig:autozero_time} to figure \ref{fig:autozero_raw_time} on page \pageref{fig:autozero_raw_time} it is immediately evident that the $f^{-1}$ flicker noise component is no longer present. The difference in white noise strength is difficult to compare and it must be turned to the power spectral density again. When calculating the spectral density it is important to remember that the sampling rate is now halved because the odd samples were subtracted. The result is shown in figure \ref{fig:autozero_psd} along with dashed lines showing the noise content prior to applying the autozero algorithm as was done in figure \ref{fig:autozero_raw_psd} on page \pageref{fig:autozero_raw_psd}.

The power spectral density in figure \ref{fig:autozero_psd} confirms an increase in the white noise power as discussed above and using the graph it can be worked out that the white noise power $\sqrt{h_{-1}}$ has increased from \qty[power-half-as-sqrt, per-mode=symbol]{165}{\nV \Hz\tothe{-0.5}} to \qty[power-half-as-sqrt, per-mode=symbol]{489}{\nV \Hz\tothe{-0.5}}, an increase by a factor of $\sqrt{8.8}$, which is more than estimated from equation \ref{eqn:autozeroing}. Including the factor of $2$ introduced by the decimation, the increase of $\sqrt{h_{-1}}$ was gauged to be by a factor of $\sqrt{4}$ . The cause of the additional noise was already mentioned above. There is still some substantial $f^{-1}$ noise present at the autozero frequency of \qty{5}{\Hz}. This type of noise is not uncorrelated and therefore the covariance is not zero, hence equation \ref{eqn:adding_white_noise} does not strictly hold and additional correlated noise is leaking into the result. This hypothesis can be confirmed by increasing the sampling frequency by an order magnitude. Doing this, the white noise floor of the autozero measurement now only increases by a factor of $\sqrt{4.5}$, which is close to the expected factor of $\sqrt{4}$. This means that the autozeroing frequency should be chosen to be at least a decade above the noise corner frequency to be most effective.
\begin{figure}[ht]
    \centering
    \input{images/autozero_psd.pgf}% data/simulations/sim_autozero.py
    \caption{Simulated power spectrum of a Keysight \device{3458A} with autozeroing applied. The dashed lines denote the noise present prior to applying the autozero algorithm. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_psd}
\end{figure}

Nonetheless, down to very low frequencies, $f^{-1}$ noise is effectively suppressed and the spectral density is almost perfectly flat. For reference, the dashed lines show the noise content that was present in the dataset prior to autozeroing, which is less white noise, but far more flicker noise.
\begin{figure}[hb]
    \centering
    \input{images/autozero_adev.pgf}% data/simulations/sim_autozero.py
    \caption{Simulated Allan deviation of a Keysight \device{3458A} with autozeroing applied. The dashed lines denote the deviation prior to applying the autozero algorithm. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_adev}
\end{figure}

The Allan deviation plot in figure \ref{fig:autozero_adev} also confirms that white noise is the only component and shows a $\tau^{-\frac 1 2}$ dependence for the full range of integration times.

From this plot it can be seen that for measurement times longer than about \qty{2}{\s} or \qty{100}{\plc}, autozeroing has a clear benefit over a measurement without autozeroing. It must be noted though that, judging from this simulation, the device would reach a noise floor of \qty[per-mode = symbol]{0.01}{\uV \per \V} only at integration times of slightly more than \qty{10}{\s}, while the datasheet claims \qty{2}{\s}. Do note, that this simulation is for the \qty{10}{\V} range of the DMM and therefore \qty[per-mode = symbol]{0.01}{\uV \per \V} is \qty{0.1}{\uV_{rms}}. It is therefore likely that the noise parameters of a real device are better than the numbers used in the simulation. Additionally, the datasheet likely refers to an instrument that is synced to a \qty{60}{\Hz} power line frequency which shifts the sampling frequency up by \qty{20}{\percent} and, as discussed, reduces the noise floor because more noise content is white noise at the autozero interval. In this simulation the \qty{0.01}{\uV \per \V} (\qty{0.1}{\uV_{rms}}) noise level would be reached at exactly \qty{10}{\s} when using a line frequency of \qty{60}{\Hz}. For the purpose of demonstrating the autozeroing algorithms these subtleties are irrelevant.

For the comparison of different ADC integration intervals before applying autozeroing figure \ref{fig:autozero_nplcs_adev} can be consulted. Using the Allan deviation makes it very simple to compare noise figures for identical measurement times $\tau$, yet different integration times, before autozeroing is applied. The simulation source code can be found in \external{data/simulations/sim\_optimal\_autozero.py} as part of the online supplemental material \cite{supplemental_material}.
\begin{figure}[ht]
    \centering
    \input{images/autozero_nplcs_adev.pgf}% data/simulations/sim_optimal_autozero_v2.py
    \caption{Allan deviation for different ADC integration intervals before applying the AZ algorithm. Dead time $\theta = \qty{0}{\s}$. The dashed line denotes the Allan variance without autozeroing. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_nplcs_adev}
\end{figure}

It can be seen that with an increasing integration time before applying the AZ algorithm more uncertainty is accumulated due to the $f^{-1}$ content which cannot be filtered. As a result, after removing the $f^{-1}$ content using autozeroing more time is required for filtering until the same Allan deviation can be reached. From these simulations it can be concluded that if there is only a negligible dead time $\theta$ involved when switching the inputs, it is advantageous to switch early, while white noise is still dominating the noise content.% The two lines for \qty{1}{\PLC} and \qty{2}{\PLC} are clearly the ones with the lowest uncertainty

Finally, the case of a non-negligible dead time shall be treated. When the dead time has to be considered, it is clear that the autozero frequency cannot be arbitrarily increased, because an increasing proportion of sampling time is lost to the dead time. This effective loss in sampling time then increases the noise spectral density due to aliasing as discussed above. To show this effect, the simulation above is modified to include a dead time of \qty{1}{\plc} as detailed in figure \ref{fig:dmm_autozer_offset_nulling}. The dead time is added once after each measurement because the input is switched after each measurement. There are also alternative switching patterns like the one proposed by \citeauthor{autozero_with_dead_time} \cite{autozero_with_dead_time} splitting the measurement interval in two and instead of measuring HI-LO-HI-LO-HI-LO, to measure HI-LO-LO-HI-HI-LO. This scheme has both advantages and disadvantages, because $f^{-1}$ flicker noise is correlated and its autocorrelation function decays with $\text{const.} - \ln(\tau)$ \cite{flicker_noise_autocorrelation,flicker_noise_autocorrelation2}. Therefore constantly changing the order of subtracted samples is not as efficient in removing the noise as the normal autozero procedure because neighbouring samples are highly correlated. Only when the dead time is large in comparison to the measurement time, this method yields an advantage. Some measurements also allow for another scheme. If the measurement is differential, the HI and LO input can be inverted without incurring the noise penalty of equation \ref{eqn:autozeroing} because both measurements taken contain the desired data. This puts the autozeroing closer to a synchronous detection scheme, but this is outside the scope of this discussion. For the sake of simplicity, only the case of a HI-LO-HI-LO measurement mentioned first is treated here. To compare the zero dead time case with the non-negligible dead time case, the Allan deviation for different integration times is again evaluated in the same way as it was in figure \ref{fig:autozero_nplcs_adev}. The results are shown in figure \ref{fig:autozero_deadtime_nplcs_adev}.
\begin{figure}[ht]
    \centering
    \input{images/autozero_deadtime_nplcs_adev.pgf}% data/simulations/sim_optimal_autozero_v2.py
    \caption{Allan deviation for different ADC integration intervals before applying the AZ algorithm. Dead time $\theta = \qty{1}{\plc}$. The dashed line denotes the Allan variance without AZ. The line frequency is \qty{50}{\Hz}.}
    \label{fig:autozero_deadtime_nplcs_adev}
\end{figure}

Figure \ref{fig:autozero_deadtime_nplcs_adev} demonstrates that the effectiveness of the AZ scheme no longer keeps increasing with an ever rising switching frequency. Instead, there is an optimal autozero interval. Above this optimal frequency, the portion of time spent with dead time is getting too large and too little information is collected. For the parameters chosen for this simulation ($f_c = \qty{1.5}{\Hz}$ and \qty[power-half-as-sqrt, per-mode=symbol]{165}{\nV \Hz\tothe{-0.5}}), \qty{5}{\plc} at \qty{50}{\Hz} is the optimal interval. If the corner frequency is shifted to a lower frequency, the optimum shifts more towards \qty{10}{\plc}. The same goes for a higher line frequency of \qty{60}{\Hz}. This explains, why HP chose \qty{10}{\plc} as the maximum integration time. For integration times higher than that, software averaging is used delivering the performance shown in figure \ref{fig:autozero_deadtime_nplcs_adev} along the \qty{10}{\plc} line.

It should be stressed here that the dead time is not the only factor to consider when choosing the autozero interval. For example, in case of an amplifier, switching the input also adds an error current due to the charge injection of the switching transistors. This may negatively impact the measurement of a high impedance source. These additional drawbacks are implementation specific and must already be considered during the design phase.

\subsection{Gain Correction}%
\label{sec:autozero_gain}
The effect of the gain correction, where the input value $x$ is scaled by a scaling factor $y$ to adjust the gain error, can be calculated, assuming white noise, as follows:
\begin{align}
    \sigma_{x \cdot y}^2 &= \langle x^2 y^2 \rangle - \langle x y \rangle^2 \nonumber\\
    &= \langle x^2 \rangle \langle y^2 \rangle + \underbrace{2\,\mathrm{Cov}\left(x^2,y^2\right)}_{\text{uncorrelated} \, = \, 0} - \left( \langle x \rangle \langle y \rangle + \underbrace{2\,\mathrm{Cov}\left(x,y\right)}_{=\, 0} \right)^2 \nonumber\\
    &= \left(\sigma_x^2 + \langle x \rangle^2\right) \cdot \left(\sigma_y^2 + \langle y \rangle^2\right) - \langle x \rangle \langle y \rangle \nonumber\\
    &= \sigma_x^2 \sigma_y^2 + \sigma_x^2 \langle y \rangle^2 + \sigma_y^2 \langle x \rangle^2 \label{eqn:variance_multiplied}
\end{align}

With respect to the gain correction, equation \ref{eqn:variance_multiplied} can be further simplified. The scaling factor is derived from the reference voltage $V_{ref}$ and normalised using $\frac{V_{ref, measured}}{V_{ref}}$. The expected value, therefore is $\langle y \rangle \approx 1$, as the ADC full scale gain should not drift much. Furthermore, $\sigma_y^2$ is scaled by the constant $1/V_{ref}$ and $\sigma_x^2 \sigma_y^2 \ll \sigma_x^2$. The latter should be true for any measurement of significance.
\begin{equation}
    \sigma_{x \cdot y}^2 \approx \sigma_x^2 + \sigma_y^2 \langle x \rangle^2
\end{equation}

The gain correction noise therefore behaves similar to the offset correction case, except that it scales with the input voltage $x$ and has no effect with a shorted input, while fully introducing its additional noise when a full scale input is applied.


% check \cite{psd_to_adev} Appendix II for details on dead time
% Compare PSD in Generation-Recombination Noise, Allan Variance, and Low-Frequency Gain Instabilities in Microwave Amplifiers to our controller. The hump look similar. Due to popcorn noise

\clearpage
\section{Current Sources}
% TODO: The FET Constant-Current Source/Limiter
Throughout this work the concept of current sources is widely used, for example section \ref{sec:laser_current_driver} discusses a current source to drive laser diodes and the temperature controller discussed in section \ref{sec:temperature_controller} uses a current source to measure the resistance of a temperature sensitive resistor. While there are many more use cases, this section will limit the discussion to a few examples used by the devices presented in this work. Namely, a unidirectional transconductance amplifier with an operational-amplifier in conjunction with a field-effect transistor and a bidirectional Howland current pump invented by Bradford Howland in 1962 and first published in 1964 by \citeauthor{howland_current_source} \cite{howland_current_source}. The discussion will start with the properties of the ideal current source and, based on this, develop a more accurate model. The models developed typically represent the static, time-independent case unless explicitly stated. First, the unidirectional current source is treated, then the bidirectional Howland current pump is discussed.

\subsection{Current Sink and Current Source}%
\label{sec:current_sink_current_source}
The question whether to use a current source or a current sink is elemental for the design of a laser driver. Figure \ref{fig:current_sink_source} shows different configurations of current sinks and sources with respect to the laser diode.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.225\linewidth}
        \centering
        \import{figures/}{current_source_high.tex}
        \caption{Source with\protect\\grounded LD.}
        \label{fig:current_source_high}
    \end{subfigure}
    \begin{subfigure}{0.225\linewidth}
        \centering
        \import{figures/}{current_source_low.tex}
        \caption{Source with\protect\\floating LD.}
        \label{fig:current_source_low}
    \end{subfigure}
    \begin{subfigure}{0.225\linewidth}
        \centering
        \import{figures/}{current_sink_high.tex}
        \caption{Sink with\protect\\grounded LD.}
        \label{fig:current_sink_high}
    \end{subfigure}
    \begin{subfigure}{0.225\linewidth}
        \centering
        \import{figures/}{current_sink_low.tex}
        \caption{Sink with\protect\\floating LD.}
        \label{fig:current_sink_low}
    \end{subfigure}
    \caption{Different configurations of current sinks and sources with respect to the laser diode. A green check mark denotes a fail-safe configuration when accidentally shorting one or more pins of the diode to the laser chassis, illustrated by a dashed connection.}
    \label{fig:current_sink_source}
\end{figure}

The optimal configuration depends on the laser diode and safety aspects in terms of protecting the laser diode. The protection of the laser diode is discussed first. The laser resonator is assumed grounded in the setup. This is not the design case, but incorrect assembly can facilitate this condition. While not intended, there are numerous ways to also accidentally short-the diode to ground and since there are no immediate consequences arising from it, when the controller is disconnected, it might easily be overlooked. This blunder should not bear the risk of destroying an expensive laser diode. To ensure this, a configuration where the laser diode is shorted out, instead of the current source or sink, must be chosen. That way, the laser diode is automatically removed from the circuit in case of an error condition.
Choosing between a current sink and a current source is more subtle. If the other shell of the laser diode is connected to the anode, a current sink can be considered to keep the diode can at ground potential. This is not an issue with the laser design in this group though, because the laser diode mount is floating. Another aspect is the electronics. A current source is typically implemented using p-channel field-effect transistors, while current sinks are using n-channel transistors and additionally the input of a current source is referenced to the positive supply, while the sink is referenced to the negative supply. Using the negative supply as a reference for control signals brings more challenges than vice versa, because typically integrated components like digital-to-analog converters prefer working with positive voltages and would need additional support to be floated to a negative reference. This makes a current source simpler to implement in this scenario and this work focuses on the current source. In principle all methods that will be discussed can be applied to a current sink as well.

\subsection{Ideal Current Source}
\label{sec:ideal_current_source}
The ideal current source as shown in figure \ref{fig:ideal_current_source} has two major properties besides the output current $I_{out}$, the output impedance $R_{out}$ and the compliance voltage, which are best understood when looking at the two equivalent representations of a current source separately. On the left in figure \ref{fig:ideal_current_source_norton}, the Norton representation can be seen. Norton's theorem reduces any linear circuit to a current source, shown in green, with a parallel resistance $R_{out}$, usually called output resistance or impedance. On the right, the Thévenin representation can be seen, which simplifies a circuit as a voltage source, also shown in green, with a series resistance.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{current_source_norton.tex}
        \caption{Norton representation.}
        \label{fig:ideal_current_source_norton}
    \end{subfigure}
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{current_source_thevenin.tex}
        \caption{Thévenin representation.}
        \label{fig:ideal_current_source_thevenin}
    \end{subfigure}
    \caption{An ideal current source with output impedance $R_{out}$ and noise $e_n$.}
    \label{fig:ideal_current_source}
\end{figure}

First, the output impedance is discussed. Ideally, $R_{out}$ is infinite and all current is forced to flow through the load. Given a finite output impedance leads to a decreased accuracy of $I_{out}$, because it is influenced by the load impedance as
\begin{equation}
    I_{out} = I_{set} \cdot \frac{R_{out}}{R_{load} + R_{out}} \, .
\end{equation}

In addition to a decreased accuracy, inserting a noise voltage source between the current source and the load as shown in figure \ref{fig:ideal_current_source} in orange, has the same effect as a changing load resistance and due to the finite output impedance $R_{out}$, any voltage noise $e_n$ translates to current noise $i_n$ through the load as
\begin{equation}
    i_n = \frac{e_n}{R_{load} + R_{out}} \approx \frac{e_n}{R_{out}} \, ,
\end{equation}

again making a high output impedance desirable to suppress noise sources between the current source and the load.

Going to figure \ref{fig:ideal_current_source_thevenin} of a current source in the Thévenin representation allows discussing the compliance voltage property. As it was said above, the output impedance of an ideal current source is infinite and so is the maximum output voltage of said current source. A finite output impedance immediately implies a finite supply voltage to keep the current to a finite limit, which dictates a maximum output voltage. This is called the compliance voltage.

\subsection{The Field-Effect Transistor Current Source}%
\label{sec:mosfet_current_source}
% Good slides can be found here: https://www.ittc.ku.edu/~jstiles/312/handouts/
Given the limited supply voltage of a real current source drives the need for a resistive element that has a finite resistance and infinite, or very high, frequency dependent dynamic impedance to react to load changes. One such pass element, having these properties, is a field-effect-transistor (FET). A junction-gate field-effect transistor (JFET) or metal–oxide–semiconductor field-effect transistor (MOSFET) can be used either as a current source or sink, depending on its doping. A p-channel FET, which uses a positive doping of the channel, is a current source, while an n-channel FET works as a current sink. This discussion is focussing on the p-channel FET with MOSFETs at its centre, because it covers the bulk of the laser current driver design in section \ref{sec:laser_current_driver}.
\begin{figure}[hb]
    \centering
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{p-channel_jfet.tex}
        \caption{P-Channel JFET.}
        \label{fig:pjfet}
    \end{subfigure}
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{p-channel_mosfet.tex}
        \caption{P-Channel MOSFET.}
        \label{fig:pmos}
    \end{subfigure}
    \caption{The simplified semiconductor structure of a JFET and a MOSFET.}
    \label{fig:FETs}
\end{figure}

The difference between a JFET and a MOSFET is the gate structure as illustrated in figure \ref{fig:FETs}. While a MOSFET has an insulated gate, the JFET does not. This reduces the gate leakage current, typically by about three orders of magnitude and allows to forward bias the device since there is no diode, resulting in larger current handling capacity. So for low currents up to a few \unit{\mA} or low noise applications, JFETs are preferred, while MOSFETs can handle several hundred ampere. The same mathematical approach can be applied to both types of FETs though. The other difference between a JFET and a MOSFET is the fact that JFETs are only available as depletion-mode (normally-on) devices, while MOSFETs are available as both depletion and enhancement (normally-off) devices. The reason is the gate structure as mentioned above. An enhancement-mode device does not conduct when the gate-to-source voltage $V_{GS} = \qty{0}{\V}$, so $V_{GS}$ must be decreased or the junction enhanced for the device to allow conduction. This is not possible with an uninsulated gate like a simple n-p junction of a JFET, which would then start conducting. A p-channel depletion-mode device on the other hand conducts at $V_{GS} = \qty{0}{\V}$ and $V_{GS}$ must be increased and the junction depleted to reduce the current, which is possible with the uninsulated gate, because the n-p junction is reverse biased. The annotated circuit symbol and the quantities used to discuss the device properties are shown in figure \ref{fig:fet_symbols}.
\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{jfet_pins.tex}
        \caption{P-channel JFET.}
        \label{fig:fet_symbols_jfet}
    \end{subfigure}
    \begin{subfigure}{0.4\linewidth}
        \centering
        \import{figures/}{pmos_pins.tex}
        \caption{P-channel MOSFET.}
        \label{fig:fet_symbols_mosfet}
    \end{subfigure}
    \caption{Basic p-channel FET circuit.}
    \label{fig:fet_symbols}
\end{figure}

A p-channel FET has its source (S) connected to the positive supply and the drain (D) is connected to a more negative voltage, typically the load. For the MOSFET the gate (G) is biased below the source to allow conduction. The source is usually connected to the substrate for solitary devices as shown in figure \ref{fig:pmos}. This will be assumed in all further discussions and the consequences of a substrate that is biased differently are omitted here. The interested reader may look up these details in \cite{mosfet_details}.

As it was hinted above, if appropriately biased, a FET can be considered a voltage controlled current source. This property can be seen in figure \ref{fig:fet_curret_gate_bias}.
\begin{figure}[hb]
    \centering
    \input{images/mosfet_current_gate_bias.pgf}% data/plot_generic.py
    \caption{Simulated drain current for different gate bias voltages of an \device{IRF9610} p-channel MOSFET.}
    \label{fig:fet_curret_gate_bias}
\end{figure}

Figure \ref{fig:fet_curret_gate_bias} shows the current $I_D$ flowing out of the drain of a p-channel MOSFET over the drain-to-source voltage $V_{DS}$ which is applied across the FET. For illustrative purposes an example p-channel MOSFET was chosen and its \textit{Simulation Program with Integrated Circuit Emphasis} (SPICE) model \cite{irf9610_spice,irf9610_spice_better} was used to generate the data, yet the overall shape is the same for all FETs. For more information on modelling MOSFETs in SPICE, \citep[p. 442]{spice_mosfets} can be consulted. There are two regions, the first region, where $V_{DS} > V_{GS} - V_{th}$, demonstrates an almost linear correlation of the channel current and the voltage across the device. This is called the ohmic region, where the MOSFET behaves much like a (gate-) voltage controlled resistor and can be described \cite{shockley_fet_equations} as
\begin{equation}
    I_{D,ohmic} = \underbrace{\kappa (V_{GS} - V_{th}) V_{DS}}_{\text{ohmic}} - \underbrace{\frac 1 2 \kappa V_{DS}^2}_{\text{pinch off}} \, .
\end{equation}
For small voltages $V_{DS}$ the output current is proportional to the applied voltage $V_{DS}$ across the channel just like a normal resistor, giving rise to its name ohmic region. As the voltage increases further $I_D$ starts leveling off because $V_{DS}$ starts affecting the channel conductivity. The channel is slowly getting pinched off at one end and becomes tapered. The reason is that the voltage $V_{DS}$ is dropped across the length of the channel. This voltage drop is linear with $V_{DS}$, resulting in a $-V_{DS}^2$ dependency of the current, reducing the conductivity of the channel. $V_{th}$ is called the threshold voltage of a MOSFET or pinch-off voltage $V_p$ in case of a JFET and is the voltage at which a current starts flowing.

The parameter $\kappa$ is a device specific parameter and depends on process parameters and the geometry of the device.
\begin{equation}
    \kappa = \kappa' \frac W L = \mu C_{ox} \frac W L
\end{equation}
$\mu$ is the electron mobility, which is about \qty{1350}{\square \cm \per \V} for n-channel MOSFETs and about \qty{540}{\square \cm \per \V} for p-channel MOSFETs \cite{fet_equations}. $C_{ox}$ is the gate-oxide capacitance per unit area and determined by the thickness $t_{ox}$ of the silicon dioxide layer of the gate
\begin{equation}
    C_{ox} = \frac{\epsilon_{ox}}{t_{ox}} \approx \frac{3.9 \cdot \epsilon_{0}}{t_{ox}} \approx \frac{\qty{3.45e-11}{\F \per \m}}{t_{ox}}\,,
\end{equation}
$W$ is the width of the channel, and $L$ is the length of the channel.

The letter $\kappa$ is used here instead of the usual $k$ as it is used by \citeauthor{fet_equations} \cite{fet_equations} to avoid confusion with the Boltzmann constant $k_B$. Unfortunately, $\kappa$ is not well controlled \cite{horowitz1989}, because it is not just determined by the size, but also the doping of the material. While the size of the structure can be well controlled to within a few \unit{\nm} using lithography masks, the doping is a matter of temperature and time in a diffusion furnace. The ohmic mode of operation is, for example, used in switches or linear voltage regulators to control the output voltage of the regulator, forming a low impedance voltage source and not the desired current source. This brings up the next region to discuss.

Once the voltage $V_{DS}$ has reached $V_{GS} - V_{th}$, the channel is fully pinched off, any further increase in $V_{DS}$ will not lead to an increase in $I_D$, in other words the output resistance becomes infinite. The MOSFET is said to be pinched-off or in saturation. In practice there still is a small influence of $V_{DS}$ on the channel. While the depth can no longer decrease as its length is \num{0} at one end already, the channel will retract a small amount in length with increasing $V_{DS}$. This is taken into account by the factor $\lambda$, called channel-length modulation. The drain current in saturation can now be described \cite{shockley_fet_equations} as
\begin{equation}
    I_{D,sat} = \underbrace{\frac 1 2 \kappa \left(V_{GS} - V_{th} \right)^2}_{\text{ideal FET}} (1 + \lambda V_{DS}) \, . \label{eqn:mosfet_saturation}
\end{equation}

The parameter $\lambda$ is the first order Taylor expansion of the length dependence of $\kappa$ and typically is small and on the order of \qtyrange[range-units = single, per-mode=power]{0.01}{0.05}{\per \volt} for p-channel MOSFETs \citep[p. 23]{mosfet_flicker_noise}. It mainly depends on the length of the channel to which it is inversely proportional, since the channel length defines the slope of the tapered channel. Sometimes the value $\frac{1}{\lambda}$ is also referred to as the Early voltage $V_A$. It is noteworthy that more modern processes choose a smaller channel length to reduce the on-state resistance of the MOSFET because the main application of a MOSFET nowadays is as a switch. The reduced channel length makes the MOSFET more susceptible to the channel length modulation effect. This will be discussed in more detail in section \ref{sec:component_selection}, when choosing a suitable MOSFET.

Going back to figure \ref{fig:fet_curret_gate_bias} the effect of the channel-length modulation can be seen as a small slope of $I_D$ in the saturation region.

Combining the previous equations, the FET drain current behaviour can be summed up as
\begin{equation}
    I_D = \begin{cases}
        0 & \text{if } V_{GS} - V_{th} < 0\\
        \kappa (V_{GS} - V_{th}) V_{DS} - \frac 1 2 \kappa V_{DS}^2 & \text{if } V_{GS} - V_{th} >= 0 \text{ and } V_{DS} < V_{GS} - V_{th}\\
        \frac 1 2 \kappa \left(V_{GS} - V_{th} \right)^2 (1 + \lambda V_{DS}) & \text{if } V_{GS} - V_{th} >= 0 \text{ and } V_{DS} \geq V_{GS} - V_{th}
    \end{cases}
    \label{eqn:mosfet_id_large_signal}
\end{equation}

The saturation region is the region of interest for building a high output impedance current source, because for a wide range of $V_{DS}$ the current remains almost constant and can be adjusted using the gate voltage $V_{GS}$. As a reminder, for the p-channel MOSFET, all voltages are reversed. $V_{GS}$, $V_{th}$, $V_{DS}$, $\kappa$ and $I_D$ are negative. Some datasheets therefore only give the magnitude of those quantities. The important aspect to remember is that for the p-channel enhancement-mode MOSFET the gate must be biased negative with respect to the source pin by a least the threshold voltage ($V_{GS} < V_{th}$ or $|V_{GS}| > |V_{th}|$) to turn the transistor on and allow current to flow.

Before proceeding to the precision current source in section \ref{sec:precision_current_source}, the concept of conductance and transconductance must be explored. The transconductance describes the relationship of the input voltage with the output current. The conductance is a measure for how well current flows from input to output. The transconductance $g_m$ and the channel conductance $g_{DS}$ are defined as
\begin{align}
    g_{m, sat} &\coloneqq \left. \frac{\partial I_{D,sat}}{\partial V_{GS}} \right|_{V_{DS} = const} = \kappa \left(V_{GS} - V_{th} \right) (1 + \lambda V_{DS}) \, , \label{eqn:mosfet_gm}\\
    &= \sqrt{2 \kappa I_D \left(1+ \lambda V_{DS}\right)} \approx \sqrt{2 \kappa I_D} \label{eqn:mosfet_gm_approximation} \\
    g_{DS, sat} &\coloneqq \left. \frac{\partial I_{D,sat}}{\partial V_{DS}} \right|_{V_{GS} = const} = \frac{1}{2} \kappa \left(V_{GS} - V_{th} \right)^2 \lambda\\
    &= \frac{I_D}{\frac{1}{\lambda} + V_{DS}} = \frac{1}{R_o} \approx I_D \lambda \label{eqn:mosfet_gds}\,.
\end{align}
The transconductance $g_m$, as a measure of the current gain with respect to the gate-source voltage of the MOSFET, is proportional to the square root of the drain current $I_D$. The inverse of the channel conductance $g_{DS}$ is called output resistance $R_o$ and discussed below. Typically the $V_{DS}$ term in the denominator of the output resistance in equation \ref{eqn:mosfet_gds} can be neglected.

The meaning of $g_{m}$ and $g_{GS}$ can be best understood when looking at a mathematical model of the MOSFET. These models come in varying complexity and either as a large-signal or small-signal model. Only the latter is used here. The small-signal model, is a first order Taylor approximation around the working point, for a constant gate-source voltage $V_{GS}$ and constant drain-source $V_{DS}$, hence both $g_{m}$ and $g_{GS}$ are constants.
\begin{align}
    I_D &\approx \frac{\partial I_D}{\partial V_{GS}} \Delta V_{GS} + \frac{\partial I_D}{\partial V_{DS}} \Delta V_{DS}\\
    &= g_{m} \Delta V_{GS} + g_{DS} \Delta V_{DS}\\
    &= g_{m} v_{GS} + \frac{1}{R_o} v_{DS} = i_D \label{eqn:mosfet_id_small_signal}
\end{align}
The lower case letters denote the variables of the small-signal model as they only change very little compared to the working point parameters.
From \ref{eqn:mosfet_id_small_signal} it can be seen that the $g_{DS}$ term adds to the output current and is proportional to $v_{DS}$. Comparing this with figure \ref{fig:ideal_current_source_norton}, the proportionality constant can be identified as $\frac{1}{R_o}$ like proposed above. Just like the ideal current source in figure \ref{fig:ideal_current_source}, the model can be given in the Norton or Thévenin representation, both shown in figure \ref{fig:mostfet_small_signa_model}.
\begin{figure}[hb]
    \centering
    \begin{subfigure}{0.43\linewidth}
        \centering
        \import{figures/}{mosfet_small_signal.tex}
        \caption{Small-signal model of a saturated MOSFET including the output resistance. The output resistance models the channel-length modulation as given by equation \ref{eqn:mosfet_id_small_signal}.}
        \label{fig:mostfet_small_signa_model_model_norton}
    \end{subfigure}
    \begin{subfigure}{0.43\linewidth}
        \centering
        \import{figures/}{mosfet_small_signal_t-model.tex}
        \caption{MOSFET model in Thévenin representation.}
        \label{fig:mostfet_small_signa_model_thevenin}
    \end{subfigure}
    \caption{Equivalent MOSFET models in Norton and Thévenin representations.}
    \label{fig:mostfet_small_signa_model}
\end{figure}

A detailed graphic derivation of the Thévenin representation can be found in \cite{fet_equations}. The Thévenin representation will prove especially valuable when treating circuits with a resistance in the source leg.
The small-signal model now shows that the output impedance is dependent on the channel-length modulation $\lambda$ and $v_{DS}$. Typically, $\frac{1}{\lambda} \gg v_{DS}$, so $\lambda$ is the most important factor governing the output impedance of a MOSFET.

To give an example of the output impedance of a MOSFET, parameters were taken from the aforementioned SPICE model of the \device{IRF9610}. Do note that these parameters of the model are tuned to match certain operating conditions by their creators and only present an estimation of the real MOSFET. Using the example parameters from table \ref{tab:current_source_parameters}, $I_D=\qty{250}{\mA}$, $\lambda = \qty[per-mode=power]{4}{\per \milli \volt}$, $V_{DS}=\qty{3.5}{\V}$ equation \ref{eqn:mosfet_gds} yields
\begin{equation}
    R_{out} = R_{o}\left(I_D=\qty{250}{\mA}, \lambda = \qty[per-mode=power]{4}{\per \milli \volt}\right) = \qty{1014}{\ohm} \overset{V_{DS} = 0}{\approx} \qty{1}{\kilo \ohm} \, , \label{eqn:mosfet_rout_irf9610}
\end{equation}
which is not very convincing as a current source. The insignificant impact of $V_{DS}$ on the output impedance can be seen when dropping the $V_{DS}$ term, which leads to an output impedance of \qty{1}{\kilo \ohm}. In textbooks this dependence is therefore usually neglected. To improve $R_{out}$, the focus thus lies on the $\lambda$ dependence. The model derived from equation \ref{eqn:mosfet_id_small_signal} can be used to do so, leading to the precision current source presented next.

% This can be demonstrated building a simple current source and then cascoding it. A very simple current source can be built using a JFET. As mentioned above, a JFET is a depletion-mode device and is already turned on at $V_{GS} = \qty{0}{\V}$. To turn it off the gate voltage must be increased above the source leg. For an illustration, refer to figure \ref{fig:jfet_curret_gate_bias}, which is very similar to the MOSFET behaviour.
%
% \begin{figure}[ht]
%     \centering
%     \input{images/jfet_current_gate_bias.pgf}
%     \caption{Simulated drain current for different gate bias voltages of a \device{2N5460} p-channel JSFET.}
%     \label{fig:jfet_curret_gate_bias}
% \end{figure}
%
% The topmost curve of figure \ref{fig:jfet_curret_gate_bias} is the case with a direct connection of the gate to the source. Above about $V_{DS}=\qty{4}{\V}$, the JFET works as a current source, although the effect of the Early voltage given in equation \ref{eqn:mosfet_saturation} can be clearly seen with a slight dependence of the output current on $V_{DS}$. With increasing $V_{GS}$, it can be observed that the slope of $I_D$ flattens.
%
%
% An example for a cascode is shown in figure \ref{fig:current_source_jfet_cascode}.
%
% \begin{figure}[ht]
%     \centering
%     \begin{subfigure}[t]{0.3\linewidth}
%         \centering
%         \import{figures/}{current_source_fet_no_bias.tex}
%         \caption{JFET current source.}
%         \label{fig:current_source_jfet_no_bias}
%     \end{subfigure}%
%     %\hfill%
%     \begin{subfigure}[t]{0.3\linewidth}
%         \centering
%         \import{figures/}{current_source_fet_bias.tex}
%         \caption{Self-biased JFET current source.}
%         \label{fig:current_source_jfet_bias}
%     \end{subfigure}%
%     %\hfill%
%     \begin{subfigure}[t]{0.3\linewidth}
%         \centering
%         \import{figures/}{current_source_fet_cascode.tex}
%         \caption{Cascoded JFET current source.}
%         \label{fig:current_source_jfet_cascode}
%     \end{subfigure}
%     \caption{Different types of JFET current sources with increasing output impedance.}
%     \label{fig:current_source_jfet}
% \end{figure}

\subsection{Precision Current Source}%
\label{sec:precision_current_source}
In the previous section \ref{sec:mosfet_current_source} it was shown in equation \ref{eqn:mosfet_id_small_signal} that the output impedance of a MOSFET depends on the channel-length modulation $\lambda$ and is too low for practical purposes. On the quest to improve the output impedance of the MOSFET circuit in figure \ref{fig:mostfet_small_signa_model_model_norton}, the most obvious solution would be to simply add a source resistor $R_s$ into the circuit as shown in in figure \ref{fig:pmos_current_source_resistor}. At first glance this may seem to only add a series resistance to $R_o$, but the attempt is more intriguing and will lead to an even better solution.
\begin{figure}[ht]
    \centering
    \begin{subfigure}[t]{0.45\linewidth}
        \centering
        \import{figures/}{current_source_resistor.tex}
        \caption{MOSFET with source resistor $R_s$ to improve the output impedance $R_{out}$.}
        \label{fig:pmos_current_source_resistor}
    \end{subfigure}%
    \begin{subfigure}[t]{0.45\linewidth}
         \centering
         \import{figures/}{pmos_small_signal_resistor.tex}
         \caption{Small-signal Thévenin model.}
         \label{fig:pmos_current_source_resistor_small_signal}
     \end{subfigure}%
     \caption{Circuit of a MOSFET with source degeneration resistor and equivalent Thévenin model.}
\end{figure}

Before calculating the output impedance, we shall have a look at $v_{GS}$ and the input signal $v_i$ derived from it. With the introduction of the source resistor $R_s$, $v_i$ no longer equals $v_{GS}$, because $\frac{1}{g_m}$ now forms a voltage divider with $R_s$ and it follows
\begin{equation}
    v_{GS} = v_i \frac{\frac{1}{g_m}}{R_s + \frac{1}{g_m}} = v_i \frac{1}{1 + g_m R_s} \,.
\end{equation}
This implies a reduction in gain by the factor $\frac{1}{1 + R_s g_m}$ compared to the previously discussed approach. The cause of this reduction is negative feedback. To understand this, imagine that with a constant $v_i$ and hence a constant current $I_D$ flowing, a changing load resistance is trying to modulate $I_D$. Any increase in $I_D$ will cause the voltage across $R_s$ to rise, reducing $v_{GS}$, because $v_i$ is still constant. The decreasing $v_{GS}$ will then reduce $I_D$, thus introducing negative feedback. Having realized there is negative feedback present, it can be postulated that the reduction in input sensitivity, or effective transconductance, will be passed on to the output impedance. This very interesting relationship will now be derived.

To calculate the output impedance, figure \ref{fig:pmos_current_source_resistor_small_signal} can be simplified by grounding $v_i$, because there is no AC component as there is no current flowing through the insulated MOSFET gate and is not modulated. The load $R_{load}$ resistance must is replaced by an AC test voltage $v_{load}$ to modulate $I_D$. These changes result in the small-signal model shown in figure \ref{fig:pmos_common_gate_amplifier}. This configuration is also called a common-gate amplifier.
\begin{figure}[ht]
    \centering
    \import{figures/}{mosfet_small_signal_cg.tex}
    \caption{Small-signal model of the common-gate amplifier with source resistance $R_s$.}
    \label{fig:pmos_common_gate_amplifier}
\end{figure}

The (dynamic) output impedance is given by
\begin{equation}
    R_{out,cg} = \frac{v_{load}}{i_D}\,, \label{eqn:mosfet_rout}
\end{equation}
with $i_D = i_S$, since there is no gate current. $v_{load}$ can easily be calculated by looking at figure \ref{fig:pmos_common_gate_amplifier} and equals the total voltage across $R_o$ and $R_s$. $v_{GS}$ can also be found, because the gate is grounded. With the resistance $\frac{1}{g_m}$ at one end, the voltage at the source pin must be $-v_{GS}$.
\begin{align}
    v_{load} &= \left(i_D - i\right) R_o + i_S R_s \nonumber\\
    &= \left(i_D - g_m v_{gs}\right) R_o + i_D R_s \nonumber\\
    &= \left(i_D + g_m i_D R_s\right) R_o + i_D R_s \label{eqn:mosfet_cg_vout}
\end{align}
Using equations \ref{eqn:mosfet_rout} and \ref{eqn:mosfet_cg_vout} gives
\begin{equation}
    R_{out,cg} = \left(1 + g_m R_s\right) R_o + R_s \label{eqn:mosfet_cg_rout}
\end{equation}
for the output impedance.

This result is interesting, as it can be be immediately seen that the output impedance scales very quickly with the transconductance $g_m$ and $R_s$. As it was already speculated above, the reduction in the transconductance $\frac{1}{1 + g_m R_s}$ of the MOSFET is transferred to the output impedance, which is increasing by the inverse of the loss in transconductance.

Going back to the quest for an increased output impedance, it is apparent that increasing $R_s$ quickly raises the output impedance, as it scales with $gm_m R_o$, but it would come at the cost of a significantly reduced compliance voltage. Therefore, other means need to be explored. As we have seen, the scale factor $gm_m R_o$ is explained by feedback and this leads to another solution. The amount of feedback can be increased further using an operational amplifier (op-amp) as shown in figure \ref{fig:precision_current_source}.
\begin{figure}[ht]
    \centering
    \import{figures/}{precision_current_source.tex}
    \caption{Transconductance amplifier with a p-channel MOSFET.}
    \label{fig:precision_current_source}
\end{figure}

The output impedance of this transconductance amplifier is amplified by the open-loop gain of the op-amp as shown in appendix \ref{sec:transfer_function_transconductance}, while the transfer function greatly simplifies to
\begin{align}
    R_{out} &\approx A_{ol} \left(g_m R_o R_s + R_o + R_s \right) \nonumber\\
    I_{out} &\approx \frac{V_{ref}}{R_s} \label{eqn:current_source_transfer_function}
\end{align}

In addition to the increased output impedance, the current $I_D = I_{out}$ can now steered by adjusting $V_{ref}$ and is, given sufficient loop gain of the op-amp, no longer dependent on the MOSFET but rather only on the sense resistor $R_s$.

This has the added benefit that it is possible to leverage the tight accuracy and precision of a resistor over the poor specifications of a MOSFET. Resistors can be manufactured with tolerances of less than \qty{100}{\micro \ohm \per \ohm}, which is orders of magnitude better than FETs, which can be matched to low \unit{\percent} values with patience.

Using the example parameters from table \ref{tab:current_source_parameters}, the output impedance in saturation can now be calculated again for $I_{out}=\qty{250}{\mA}$ and the ideal \device{IRF9610} model with the addition of an idealized \device{AD797} op-amp using the worst-case specifications.
\begin{equation}
    R_{out} \approx \qty[per-mode=power]{2}{\volt \per \uV} \left(\qty{0.64}{\siemens}\cdot \qty{1014}{\ohm} \cdot \qty{30}{\ohm} + \qty{1014}{\ohm} + \qty{30}{\ohm} \right) \approx \qty{40}{\giga\ohm} \label{eqn:current_source_output_impedance}
\end{equation}

From these consideration, it can be seen that the open-loop gain and the unity-gain bandwidth of the op-amp essentially determine the properties of the current source, given that $R_{id} \gg R_s$ and $R_o \gg R_s$. This will be important for selecting an operational amplifier later.

The next section will focus on the MOSFET and discuss the compliance voltage of the current source, which was only briefly touched during the introduction. It will give rise to criteria for selecting a MOSFET for the precision current source.

\subsection{Compliance Voltage}%
\label{sec:compliance_voltage}
The compliance voltage of a current source is the maximum voltage it can output to maintain the requested output current. For an ideal current source, the compliance voltage is infinite, but it is obviously limited in the physical world.

The precision current source discussed in section \ref{sec:precision_current_source} has several limiting factors of the compliance voltage, which shall be discussed now. The compliance voltage is taxed most at the maximum output current $I_{out,max}$. Thus for the following discussion, the output is always treated as set to maximum.

Looking at figure \ref{fig:precision_current_source} of the precision current source it is immediately evident that the output voltage can be calculated by subtracting the voltage across the source resistor $V_{R_s}$ and the MOSFET $V_{DS}$ from the supply voltage $V_{sup}$
\begin{equation*}
    V_{out} = V_{sup} - V_{R_s} - V_{DS} = V_{sup} - V_{ref} - V_{DS}\,.
\end{equation*}

The voltage $V_{R_s}$ is given by equation \ref{eqn:current_source_transfer_function} and equal to the setpoint voltage and hence given by the system parameters. This leads to the question of the minimum working point voltage $V_{DS}$ at $I_{out,max}$. As a reminder, from equation \ref{eqn:mosfet_id_large_signal} and figure \ref{fig:fet_curret_gate_bias} one can see that the drain current is almost constant over $V_{DS}$ in the saturation region, and in the ohmic region is proportional to $V_{DS}$. The transition point from the ohmic region to the saturation region is at $V_{DS} = V_{GS} - V_{th}$ and putting this into equation \ref{eqn:mosfet_id_large_signal} yields for the drain current
\begin{align}
    I_D &= \frac{1}{2} \kappa V_{DS}^2 \left(1+ \lambda V_{DS}\right) \nonumber\\
    \Rightarrow V_{DS} &\approx \sqrt{\frac{2 I_D}{\kappa}}\\
    &\approx \qty{784}{\mV}
\end{align}

The latter result was calculated using the example parameters from table \ref{tab:current_source_parameters}. At this point it can already be postulated that the MOSFET will severely change in its function as a current source for $V_{DS} < \qty{0.78}{\V}$. To quantify this, one has to look at the output impedance of the transconductance amplifier once again. In the last section, the output impedance was only treated for the saturation region, but this time, $R_{out}$ must be considered over a wide range of $V_{DS}$, thus not only in the saturation region but also in the ohmic region. Instead of using the small-signal model as before, which assumed only small changes of $V_{DS}$, a large-signal model must be applied, which also includes the non-linear nature of the piece-wise defined equation \ref{eqn:mosfet_id_large_signal} of the drain current.

For the sake of simplicity, a SPICE simulation of figure \ref{fig:precision_current_source} was carried out in LTSpice \cite{ltspice}. Solving this analytically bears no educational value over the numerical solution shown below as will be seen. Additionally, the SPICE simulation also offers the opportunity to add additional, parasitic elements to the model to evaluate their effect, for example, the capacitive nature of the MOSFET gate.

The simulation itself is numerically challenging and the typical approaches will lead to the limits of the numerical precision. To make the simulation feasible, the large-signal model is broken down into several small segments. For each of these segments, the small-signal model at its respective working point is evaluated and then the result joined back together to reconstruct the large-signal model sought. How this is done in detail, is shown in appendix \ref{sec:ltspice_current_source} as it is beyond the scope of this section. The final result was calculated for two different frequencies, one frequency was deliberately chosen so low (\qty{1}{\micro\Hz}) that it is well below the dominant pole of the op-amp, meaning that the full open-loop gain applies and the other frequency chosen was \qty{1}{\MHz}, were the gain had dropped to \qty[per-mode=power]{10}{\V \per V}. This is shown in figure \ref{fig:ltspice_output_impedance_simulation}.
\begin{figure}[ht]
    \centering
    \input{images/ltspice_output_impedance_simulation.pgf}% data/plot_generic.py
    \caption{Simulated output impedance for the precision current source from figure \ref{fig:precision_current_source} at DC and \qty{1}{\MHz} over the drain-source voltage.}
    \label{fig:ltspice_output_impedance_simulation}
\end{figure}

Looking at figure \ref{fig:ltspice_output_impedance_simulation} clearly shows the effect of entering the ohmic region of the MOSFET. Over a range of \qty{100}{\mV} starting at the \qty{0.78}{\V} calculated above, the output impedance drops by two orders of magnitude and then keeps dropping at an exponential rate with decreasing $V_{DS}$. The same effect applies to the output impedance at \qty{1}{\MHz}, although the starting value of the output impedance is around \qty{200}{\kilo\ohm} due to the reduced gain from the op-amp at \qty{1}{\MHz}. It can also be seen that $R_{out}$ levels off at \qty{30}{\ohm}, the value of the sense resistor.

This overall effect of leaving the saturation region is so drastic that the compliance voltage must be defined in such a way that the MOSFET remains in saturation and this leads to
\begin{equation}
    V_{comp} = V_{sup} - V_{ref} - \sqrt{\frac{2 I_D}{\kappa}} \,.
\end{equation}

Now turning to the supply voltage, it is limited by the op-amp which must drive the gate of the MOSFET all the way up to the supply to turn off the current source. The reference voltage is, unless one divides it down dictated by the reference chosen. This, unfortunately, leaves only little room for the MOSFET and it must be carefully chosen not limit the compliance voltage too much.

At this point a fallacy the author has observed multiple times must be addressed. In order to address the limited compliance voltage, one may be tempted to use multiple MOSFETs in parallel to divide the current between the MOSFETs and thereby reduce the voltage that needs to be dropped across the FET proportional to $\frac{1}{\sqrt{N}}$, where $N$ is the number of MOSFETs paralleled.

Imagine the following modified circuit of the precision current source shown in figure \ref{fig:precision_current_source_two_mosfets} with two MOSFETs in parallel. For clarity the gate resistors required are not included.
\begin{figure}[ht]
    \centering
    \import{figures/}{precision_current_source_2fets.tex}
    \caption{Transconductance amplifier with two p-channel MOSFETs in parallel.}
    \label{fig:precision_current_source_two_mosfets}
\end{figure}

While at first this seems like a solution to the limited $V_{DS}$, it is not recommended for a number of reasons given here.

The first reason is, MOSFET specifications are very loose, notably the threshold voltage $V_{th}$, the transconductance $g_m$ and the capacitances, but the latter is of little concern here. These tolerances limit the usefulness of paralleling MOSFETs to certain conditions, for example, when using the MOSFETs as a switch ad not as a current source. The difference between its use as a switch and a current source is the thermal load, which is a lot higher when using the MOSFET as current source. In this respect it seems to be a common misunderstanding that MOSFETs are immune to thermal runaway. This is mostly true when using them as a switch, fully turned on and in the ohmic region. In this case, there are two effects occurring, the first is that the (absolute) value of $V_{th}$ decreases with temperature, thus increasing $I_D$ and the second effect is, $R_{DS,on}$ is rising with temperature \cite{mosfet_thermal_runaway}. The latter effect is, depending on the physical design, stronger, but it depends on the specific MOSFET. A detailed analysis of paralleling MOSFETs as switches can be found in \cite{paralleling_mosfets}. In case the MOSFET is operating in pinch-off and not the ohmic region, $R_{DS,on}$ has no influence on the current, therefore, the only effect at work is the decreasing $V_{th}$. Depending on the difference in $V_{th}$ between the paralleled MOSFETs, one MOSFET will take most of the current and power. Adding source resistors can compensate for this by pushing down the source voltage as the current goes up. This will then reduce $V_{GS}$. The size of the resistor depends on the transconductance $g_m$ and the temperature coefficient of $V_{GS}$, which is around \qtyrange[range-units = single]{1.5}{2}{\mV \per \K} \cite{mosfet_vgs_tempco}. Unfortunately, \qty{1}{\ohm} or \qty{2}{\ohm}, will already use up, most of the benefits gained in compliance voltage as will be shown below.

The second reason why paralleling MOSFETs is not desirable can be seen when
remembering equation \ref{eqn:mosfet_id_large_signal}. It is known that the transition from the unwanted ohmic region to the saturation region is
\begin{equation}
    V_{DS} \geq V_{GS} − V_{th}\,.
\end{equation}

Looking at \ref{fig:precision_current_source_two_mosfets}, it can be seen that $V_{GS}$ is set by the op-amp and is the same for both MOSFETs because their gates and sources are connected. However, $V_{th}$ is device specific and according to the datasheet of the example \device{IRF9610} \cite{datasheet_IRF9610} $V_{th}$ values can show a spread of as much as \qtyrange[range-units = single, range-phrase={~to~}]{-2}{-4}{\V}, although \cite{appnote_mosfet_parameter_spread} suggests that MOSFETs from the same reel show a spread of only \qty{\pm 125}{\mV} of $V_{th}$  within the same batch for consecutive devices. The \qty{125}{\mV} was found for the \device{BUK7S1R5-40H} \cite{datasheet_BUK7S1R5}, which was sampled in this report. The number given in the report is for $3\sigma$ and, assuming the datasheet values for the spread are also referring to $3\sigma$, the spread found in the report is about twice as high as the datasheet value of \qtyrange[range-units = single]{2.4}{3.6}{\V}. Assuming similar numbers for \device{IRF9610} MOSFET used in our examples, this leads to \qty{\pm 208}{\mV} for the \device{IRF9610}, again applying $3\sigma$. Using this number, a Monte Carlo simulation (not quite, because the dice were biased to yield a Gaussian distribution instead of equal probabilities) was run using LTSpice, simulating the circuit shown in figure \ref{fig:precision_current_source_two_mosfets} and also the original circuit using only one MOSFET. For this simulation the current source was set to \qty{250}{\mA} as per table \ref{tab:current_source_parameters}. The load voltage was set to
\begin{equation*}
    V_{DS, parallel} = \sqrt{\frac{2 \frac{I_d}{2}}{\kappa}} \approx \qty{555}{\mV}\,,
    V_{DS, single} = \sqrt{\frac{2 I_d}{\kappa}} \approx \qty{784}{\mV} \,.
\end{equation*}

$\frac{I_D}{2}$ was used to calculate $V_{DS, parallel}$ for the parallel configuration to show the effect assuming perfect current sharing between the MOSFETs. Additionally, a configuration with an increased safety margin of $1 \sigma = \qty{70}{\mV}$ added to $V_{DS, parallel}$ was investigated. \num{4000} samples were drawn and the spread of the output impedance was calculated for each circuit. The results are shown as a histogram in figure \ref{fig:ltpsice_mosfet_mc_output_impedance}. The counts give the number of cases for each bin of the output impedance.
\begin{figure}[ht]
    \centering
    \input{images/ltspice_mosfet_mc_output_impedance.pgf}% data/plot_ltspice_monte-carlo.py
    \caption{Results of a Monte Carlo simulation of the output impedance for different configurations of MOSFETs.}
    \label{fig:ltpsice_mosfet_mc_output_impedance}
\end{figure}

Unsurprisingly, there is no variance of the output impedance in the single MOSFET case in accordance to what can be seen in appendix \ref{sec:transfer_function_transconductance}. The op-amp gain simply suppresses all device properties of the MOSFET. The slight variation of $g_m$ for different samples was not simulated, because this variation stems from the variation of $\kappa$ and goes as $\frac{1}{\sqrt{\kappa}}$, so its effect is not as pronounced as the threshold.

In case of two MOSFETs, the output impedance varies over an order of magnitude from about \qtyrange[range-units = single]{1.8}{52}{\giga \ohm}. Even when increasing the drain-source voltage by $1 \sigma = \qty{70}{\mV}$ to \qty{625}{\mV} on average, the spread is still an order of magnitude. Only when increasing $V_{DS}$ to around \qty{700}{\mV}, the situation stabilises, but then the net gain from this measure has shrunk to a meager \qty{84}{\mV}. It can be seen from this simulation that the system-to-system spread becomes very unstable in tough situations. Such instability can also be brought into the system by temperature effects as $V_{th}$ is temperature dependent as discussed above. Additionally, it may suffer from thermal runaway unless each individual MOSFET is laid out to carry the full current.

\subsection{Noise Sources}%
\label{sec:current_source_noise}
The fundamentals of different types of noise were already introduced in section \ref{sec:allan_deviation}. Here, a subset of those noise types is revisited. It is expected that the dominant noise observed in the current source circuit is $\frac{1}{f}$-noise at low frequencies and white noise towards higher frequencies. All noise components will be converted to the so-called input referred notation to make the noise sources comparable. This can be easily understood when looking at two amplifiers with different gain. If both of them add a fixed amount of noise to the output signal, the absolute amount of noise may be the same, but the signal to noise ratio shows a different picture. To compare these amplifiers it is useful to divide the noise by the transfer function (gain) of the amplifier. This is called the input-referred noise, since it treats the noise in relation to the input signal. Additionally, when calculating noise figures, the noise bandwidth is always considered to be \qty{1}{\Hz} unless specified otherwise.
\begin{figure}[ht]
    \centering
    \import{figures/}{precision_current_source.tex}
    \caption{Transconductance amplifier with a p-channel MOSFET. Copied from page \pageref{fig:precision_current_source}.}
    \label{fig:precision_current_source_noise}
\end{figure}

Noise sources are ubiquitous in the circuit shown in figure \ref{fig:precision_current_source_noise}, which is repeated here from figure \ref{fig:precision_current_source} on page \pageref{fig:precision_current_source} for clarity. The resistor $R_s$, the MOSFET, the op-amp, the setpoint voltage $V_{ref}$ and the supply voltage $V_{sup}$ can all contribute noise to the output current. Fortunately, some of those noise contributions are either very small or well suppressed in this design, so each component must be briefly discussed to see if they can indeed be safely neglected.

Starting with the supply voltage $V_{sup}$, it can be seen that any change of this voltage affects the string $R_s$-$Q$-$R_{load}$. From equation \ref{eqn:transconductance_amplifier_transfer_function}, it is known that if the op-amp gain is high (true within the bandwidth of the op-amp) all disturbances of the voltage across $R_s$ will be suppressed and the output current is only defined by the reference input and $R_s$. Looking closer, the supply noise is present at the inverting and non-inverting input of the op-amp with the same magnitude. If there is no current flowing into the op-amp pins, which is true for low frequencies, the noise is affecting both pins equally and it will be suppressed by the common-mode-rejection ratio (CMRR) of the device. Fortunately, this is a strong quality of precision op-amps and values of more than \qty[per-mode=power]{1}{\uV \per \volt} are not uncommon. The op-amp will therefore take care of the supply noise at low frequencies. At high frequencies the parasitic capacitance of the input pins and the reduced gain and CMRR come into play. To take care of this, it is therefore prudent to filter the supply voltage for high frequency noise.

The next noise source is the reference voltage. The reference is directly connected to the input and its noise dictates most of the circuit noise. While the high-frequency noise can again be filtered to some extend, the low frequency noise, which is mostly $\frac{1}{f}$-noise cannot be filtered as was shown in section \ref{sec:flicker_noise}, so it must be kept low from the start and the reference selected for low flicker noise.

The MOSFET as a noise source is considered in appendix \ref{sec:mosfet_noise} and the interested reader may find the derivation of the MOSFET noise component there. The two types of noise that need to be considered are the flicker noise of the MOSFET and its wideband thermal noise as calculated in equation \ref{eqn:current_noise_mosfet} on page \pageref{eqn:current_noise_mosfet}
\begin{equation*}
    i_{n} = \sqrt{\underbrace{4 k_B T \frac{2}{3} g_m}_{\text{thermal}} + \underbrace{\frac{K_f I_D}{C_{ox} L^2} \frac{1}{f}}_{\text{flicker}}} \,.
\end{equation*}

To calculate the input referred noise in order to show that the MOSFET noise will be suppressed by the op-amp, the current noise needs to be divided by the open-loop gain derived as equation \ref{eqn:transconductance_amplifier_open_loop_gain} on page \pageref{eqn:transconductance_amplifier_open_loop_gain}
\begin{equation}
    e_{n,FET} = \frac{i_n}{A_f} = \frac{\sqrt{4 k_B T \frac{2}{3} g_m + \frac{K_f I_D}{C_{ox} L^2} \frac{1}{f}}}{\frac{A_{op}}{R_s} \frac{g_m \left(R_o || R_s || R_{id}\right)}{g_m \left(R_o || R_s || R_{id}\right) + 1}} \,.
\end{equation}

Looking at the parameters from table \ref{tab:current_source_parameters}, it is found that $\left(R_o || R_s || R_{id}\right) \approx R_s$ and $e_n$ can be simplified to
\begin{align*}
    e_{n,FET} &\approx \frac{\sqrt{4 k_B T \frac{2}{3} g_m + \frac{K_f I_D}{C_{ox} L^2} \frac{1}{f}}}{A_1 \frac{1}{R_s + \frac{1}{g_m}}}\\
    &\approx \frac{R_s + \frac{1}{g_m}}{A_1} \sqrt{4 k_B T \frac{2}{3} g_m + \frac{K_f I_D}{C_{ox} L^2} \frac{1}{f}}\\
    \overset{A_1 \to \infty}&{=} 0
\end{align*}

Unless the MOSFET transconductance $g_m$ or the gain of the op-amp $A_1$ become very small, the noise of the MOSFET is very well suppressed. If the wideband thermal noise contribution is small (which it is, see \ref{sec:mosfet_noise}) and the flicker noise corner frequency is within the bandwidth of the op-amp, the noise contribution from the MOSFET can be neglected.

The noise contribution from the sense resistor $R_s$ is the Johnson–Nyquist noise, see section \ref{sec:white_noise}, which when transformed to its Norton representation can be written as a current noise
\begin{equation}
    i_{n,R} = \sqrt{\frac{4 k_B T}{R_s}} \label{eqn:current_noise_resistor}\,.
\end{equation}

Additionally, it was shown that depending on the material of the resistive element, a flicker noise component can also be present. This is especially prevalent in carbon and thick-film resistors \cite{flicker_noise_carbon_film,1_f_noise_thick_film}. While thin-film resistors are less noisy, their performance varies greatly between different models \cite{resistor_current_noise_ligo}, so their make and model must be carefully selected for the application. Metal foil and wirewound resistors were shown to perform best and have almost no flicker noise \cite{resistor_current_noise_ligo,flicker_noise_foil_resistor_beev}. Using a high quality resistor, the flicker noise can be neglected and only the thermal noise must be taken into account.

The sense resistor is part of the feedback network and therefore it contributes fully to the noise of the transconductance amplifier. Input referred, the current noise must be divided by the closed-loop gain $A_f$ given by \ref{eqn:transfer_function_closed_loop} on page \pageref{eqn:transfer_function_closed_loop}.
\begin{equation}
    e_{n,R} = i_{n,R} \cdot \beta \approx i_n \cdot R_s = \sqrt{4 k_B T R_s} \label{eqn:noise_sense_resistor}
\end{equation}

The final component to be discussed is the operational amplifier. Although the op-amp is a rather complex device, its noise can be modeled with sufficient accuracy by a small number of noise sources. The typical noise model of an op-amp is shown in figure \ref{fig:op-amp_noise_model}.
\begin{figure}[ht]
    \centering
    %\scalebox{1}{%
        \import{figures/}{op-amp_noise_model.tex}
    %} % scalebox
    \caption{Noise model of the operational amplifier.}
    \label{fig:op-amp_noise_model}
\end{figure}

In figure \ref{fig:op-amp_noise_model} one can see that there are three noise sources required to treat the op-amp. The input voltage noise source $e_{n}$ and two input current noise sources $i_n$. The current noise noise sources $i_n$ are assumed to be mostly uncorrelated. This assumption leads to an upper bound as can be seen from figure \ref{fig:op-amp_input_stage}, which shows the input differential amplifier, which is the first stage of a typical bipolar op-amp.
\begin{figure}[hb]
    \centering
    %\scalebox{1}{%
        \import{figures/}{op-amp_input_stage.tex}
    %} % scalebox
    \caption{Bipolar op-amp input stage with noise sources.}
    \label{fig:op-amp_input_stage}
\end{figure}

Of the three noise sources $i_{n,p}$, $i_{n,n}$ and $i_{n,EE}$ only $i_{n,p}$ and $i_{n,n}$ are uncorrelated, because it is the input bias current of the individual transistors, and only the effect of $i_{n,EE}$ is correlated, because the current of the emitter bias current source is equally distributed between the two input transistors. Since effects of equal magnitude and sign cancel out due to the differential nature of the input stage, correlated effects are suppressed. An equal magnitude can be assumed, because the gain of the two transistors is well matched, due to their close proximity on the semiconductor die. Therefore assuming all noise is uncorrelated presents an upper bound for the current noise $i_n$. A more detailed analysis can be found in \cite{op-amp_noise_correlation}. Due to the matching of the transistors, the magnitude $i_{n,p}$ and $i_{n,n}$ are also closely matched, hence, in the model used here, they are assumed equal and referred to as $i_n$ from now on. These two current noise sources cannot be combined with the voltage noise source, because they depend on the external impedance connected to the op-amp, so this final step must be done by the circuit designer as shown next.

Both the voltage noise $e_n$ and the current noise $i_n$ can be found in the datasheet of the op-amp. Typically, these values given are input-referred values. For the complete circuit as in figure \ref{fig:precision_current_source_noise} it is possible to calculate the full noise contribution of the op-amp as
\begin{equation}
    e_{n,op} = \sqrt{e_n^2 + e_{n,+}^2 + e_{n,-}^2}\,,
\end{equation}
assuming the noise sources are uncorrelated. The input referred voltage noise $e_{n,-}$ of the inverting input resulting from the current noise can be calculated in a similar fashion as $e_{n_R}$ in equation \ref{eqn:noise_sense_resistor}. It is likewise part of the feedback network and must therefore be divided by the closed-loop gain $A_f$ as before.
\begin{equation}
    e_{n,-} \approx i_n \cdot R_s
\end{equation}

The current noise of the non-inverting input can be translated by considering its input impedance. This is determined by the filter circuit of the reference voltage which is required to remove the high frequency noise as discussed above. Assuming an RC-filter of first order,
the output impedance can be calculated from the transfer function of the low-pass filter, derived in equation \ref{eqn:first-order_model}
\begin{align}
    R_{out,f} &= R_{f} \cdot A = \frac{R_{f}}{1+sR_{f}C}  \label{eqn:output_impedance_rc_filter}\\
    \lim_{s \to 0} R_{out,f} &= R_{f} \nonumber\\
    \lim_{s \to \infty} R_{out,f} &= 0 \,.\nonumber\\
     e_{n,+} &\approx \frac{i_n R_{f}}{1+sR_{f}C}
\end{align}

Looking at the output impedance of the filter, it can be seen that for high frequencies, the output impedance approaches \num{0}, while for low frequencies it is $R_{f}$. This has the effect that if the filter corner frequency $\omega_0 = \frac{1}{RC}$ is close to the flicker noise corner frequency of the reference voltage there is is almost no wideband current noise contribution from $e_{n,+}$ as well. Only the $\frac{1}{f}$ component of the op-amp current noise multiplied with $R_{f}$ is left. Ideally, this is lower than the reference noise to have negligible impact and must be kept in mind when selecting an op-amp.
This leads to the total noise of the op amp
\begin{equation}
    e_{n,op} = \sqrt{e_n^2 + (i_n R_s)^2 + \left|\frac{i_n R_{f}}{1+sRC}\right|^2} \,.
\end{equation}

To conclude, table \ref{tab:current_source_noise_contributers} is given as a reference for the noise contributions in the low-frequency and also the wideband domain. From this table, it can be seen that the only wideband-noise contributors are the reference resistor and the op-amp. The low-frequency contributors are the voltage reference and the op-amp, since they have a strong flicker noise component. A low-noise, precision op-amp typically has far less low frequency noise than a voltage reference and the dominant low frequency contributor remains the voltage reference.
\begin{table}[ht]
    \centering
    \begin{tabular}{lll}
        \toprule
        Noise component& Low frequency& Wideband \\
        \midrule
        $V_{sup}$ & $\approx 0$ & $\approx 0$\\
        MOSFET & $\approx 0$ & $\approx 0$\\
        $V_{ref}$ & $\sqrt{e_{n,ref}^2 + 4 k_B T R_{f}} $ & $\approx 0$\\
        $R_s$ & $\sqrt{4 k_B T R_s}$ & $\sqrt{4 k_B T R_s}$\\
        Op-amp & $\sqrt{e_n^2 + i_n^2 (R_s^2 + R_{f}^2)}$ & $\sqrt{e_n^2 + i_n^2 R_s^2}$\\
        \bottomrule
    \end{tabular}
    \caption{Input referred noise components of the transconductance amplifier. Multiply by $\frac{1}{R_s}$ to get the output referred current noise.}
    \label{tab:current_source_noise_contributers}
\end{table}

From those findings it is clear that the most important choices regarding the noise contributions are a good quality metal-foil or wirewound sense resistor $R_s$, a low noise voltage reference and a low noise op-amp. Regarding the low-noise op-amp it is critical to decide between low voltage noise or low current noise. This choice depends on the value of $R_s$. For typical values of $R_s$ below \qty{1}{\kilo\ohm}, voltage noise is the dominating effect. The reference should be chosen for low flicker noise.

\subsection{Component Selection}%
\label{sec:component_selection}
This section deals with selecting the right components for the precision current source presented in section \ref{sec:precision_current_source}. The focus lies on the requirements defined in section \ref{sec:laser_current_driver}, notably specifications \ref{lst:dgDrive_specs_environment} and \ref{lst:dgDrive_specs_electrical}. Most attention will be on the MOSFET, the operational amplifier, and the voltage reference. For these components examples from literature are given and are compared to the requirements. Discussed first is the voltage reference, because this will define several parameters down the road. Then the op-amp is considered, for which several examples from scientific publications and other alternatives are shown and the best solution is presented. Finally, the selection parameters for the MOSFET will be elaborated. The reader must be warned though that the lineup of p-channel MOSFETs in production is decreasing, with more and more products being discontinued in favour of n-channel MOSFETs. The component recommendations may therefore already be outdated.

Numerous laser driver designs for different applications and laser diodes can be found in literature \cite{libbrecht_hall, laser_driver_mosfet_noise, laser_driver_digital, laser_driver_digital_update, laser_driver_qcl_space, laser_driver_qcl_taubman, laser_driver_qcl_taubman_multiplexer}. Even though \citeauthor{libbrecht_hall} \cite{libbrecht_hall} were not the first to present their circuit, a similar solution can already be found in \cite{laser_driver_old}, their design stands out for its simplicity. The designs mentioned can be divided into two groups. High power drivers for quantum cascade lasers (QCL) typically featuring a compliance voltage $V_c$ of more than \qty{10}{\V} and output currents of up to several ampere based on the work of \citeauthor{laser_driver_qcl_taubman} \cite{laser_driver_qcl_taubman} and medium power devices for laser diodes having a lower compliance voltage of around \qty{2}{\V} and capable of driving a few hundred \unit{\mA} based on the work of \citeauthor{libbrecht_hall} \cite{libbrecht_hall}. The requirements of this work mostly fall into the latter category, except for the compliance voltage, which is targeted for blue laser diodes and $V_c \ge \qty{8}{\V}$. All these designs share one common aspect though, the type of voltage reference. Most laser drivers in literature and also commercial products are designed around low-noise, low-drift buried Zener diode voltage references, namely the Analog Devices (ADI) \device{LM399} \cite{datasheet_LM399} or ADI \device{LTZ1000} \cite{datasheet_LTZ1000}.

The buried types of voltage references are Zener diodes. that are created within the bulk silicon using ion implantation. This reduces noise caused by surface contamination \cite{zener_diode_stability}. These diodes are not true Zener diodes though, but called Zeners nonetheless, and use a mix of Zener and avalanche breakdown to compensate the temperature coefficient.

The Zener effect is the tunneling of electrons through the barrier formed between the valence band and conduction band. It has a negative temperature coefficient, because an increase in temperature reduces the size of the band gap. This effect dominates at low currents.

Avalanche breakdown, on the other, hand describes a mechanism in which free electrons (due to temperature) are accelerated to such energies that they knock out other electrons, causing an avalanche of electrons. This effect has a positive temperature coefficient \cite{tempco_avalanche_breackdown} and occurs at higher currents. While the zero temperature coefficient point is around \qty{5}{\V}, this operating point implies a high susceptibility to changes in the reverse current. Typically, the Zener voltage is shifted slightly upwards to result in a net positive coefficient, which is then compensated by the negative temperature coefficient of a forward biased diode \cite{zener_diode_stability}. This results in the common Zener diode voltage of around $\qty{6.2}{\V} + \qty{0.7}{\V} = \qty{6.9}{\V}$. In comparison to other types of diodes, buried Zeners have the best stability and lowest noise. $V_{ref} \approx \qty{7}{\V}$ is therefore set in stone for the reasons given above. Table \ref{tab:overview_buried_zener_diodes} lists some commercially available buried Zener diodes. All of these diodes are manufactured by ADI as they are the sole manufacturer left on the market to produce this kind of diodes.
\begin{table}[ht]
    \centering
    \begin{tabular}{lllll}
        \toprule
        Component& Voltage& Temperature coefficient & Stability& Package\\
        \midrule
        \device{LT1021} & \qty{7}{\V} & \qtyrange[range-units = single]{2}{5}{\uV \per \V \per \K} & \qty{15}{\uV \per \V \per \kilo\hour\tothe{0.5}} & SO-8\\
        \device{LT1027} & \qty{5}{\V} & \qtyrange[range-units = single]{1}{2}{\uV \per \V \per \K} & not specified & SO-8\\
        \device{LM399} & \qty{7}{\V} & \qtyrange[range-units = single]{0.3}{1}{\uV \per \V \per \K} & \qty{8}{\uV \per \V \per \kilo\hour\tothe{0.5}} & TO-46\\
        \device{ADR1399} & \qty{7}{\V} & \qtyrange[range-units = single]{0.2}{1}{\uV \per \V \per \K} & \qty{7}{\uV \per \V \per \kilo\hour\tothe{0.5}} & TO-46\\
        \device{LTZ1000} & \qty{7.2}{\V} & \qty{0.05}{\uV \per \V \per \K} & \qty{0.3}{\uV \per \V \per \kilo\hour\tothe{0.5}} & TO-99\\
        \device{ADR1000} & \qty{6.6}{\V} & \qty{<0.2}{\uV \per \V \per \K} & \qty{0.2}{\uV \per \V \per \kilo\hour\tothe{0.5}} & TO-99\\
        \bottomrule
    \end{tabular}
    \caption{List of commercially available buried Zener diodes and selected properties.}
    \label{tab:overview_buried_zener_diodes}
\end{table}

Choosing a voltage reference can be done according to specification \ref{lst:dgDrive_specs_environment}. A temperature coefficient of \qty{<= 1}{\uA \per \A \per \K} rules out any non-hermetic unheated voltage reference. Using a hermetic package also improves the stability against humidity as the epoxy used for an SO-8 package is hydrophilic and swells when exposed to water vapour causing pressure on the die, resulting in a change of the output voltage. The hermetic voltage references can be divided into two groups, the \device{LM399} and the newer \device{ADR1399} in one group and the \device{LTZ1000} and its newer counterpart \device{ADR1000} in another. While the \device{LM399} requires very few external components, the external circuit for the \device{LTZ1000} is far more elaborate requiring more parts and space. Additionally, the \device{LTZ1000} is more than four times the price of the \device{LM399} in quantities of \num{10} at the time of writing. Last but not least, the stability and temperature coefficient of the \device{LTZ1000} cannot be matched by the performance of the sense resistor in the precision current source, so the sense resistor gives a lower bound of about \qty{0.5}{\uA \per \A \per \K}. Unless the better low frequency noise performance is required, the \device{LM399} and \device{ADR1399} are the more economical parts. The performance of the latter two references will be discussed in section \ref{sec:zener_diode_selection} and the sense resistor is considered next.

With the maximum reference voltage of \qty{7}{\V} known, a sense resistor between \qty{14}{\ohm} (\qty{500}{\mA}) and \qty{28}{\ohm} (\qty{250}{\mA}) is required. The Vishay \device{VPR221Z} are high power low drift metal foil resistors in a TO-220 package and a solid choice here. The low value resistors, combined with the requirement for a low current noise output of the source, limits the choice of op-amps to bipolar low-noise devices or discrete implementations. Table \ref{tab:overview_bipolar_op-amps} lists some choices compiled from the literature sources, which will now be discussed.
\begin{table}[ht]
    \centering
    \begin{tabular}{llll}
        \toprule
        Component& Wideband-noise& Low frequency noise & Temperature coefficient \\
        \midrule
        \device{LT1028} & \qty[power-half-as-sqrt]{0.85}{\nV \per \Hz\tothe{0.5}} & \qty{35}{\nV_{p-p}} & \qty{0.2}{\uV \per \K}\\
        \device{AD797} & \qty[power-half-as-sqrt]{0.9}{\nV \per \Hz\tothe{0.5}} & \qty{50}{\nV_{p-p}} & \qty{0.2}{\uV \per \K}\\
        \device{ADA4898} & \qty[power-half-as-sqrt]{0.9}{\nV \per \Hz\tothe{0.5}} & not specified & \qty{1}{\uV \per \K}\\
        \device{ADA4004} & \qty[power-half-as-sqrt]{1.8}{\nV \per \Hz\tothe{0.5}} & \qty{150}{\nV_{p-p}} & \qty{0.7}{\uV \per \K}\\
        \device{AD8671} & \qty[power-half-as-sqrt]{2.8}{\nV \per \Hz\tothe{0.5}} & \qty{77}{\nV_{p-p}} & \qty{0.3}{\uV \per \K}\\
        \bottomrule
    \end{tabular}
    \caption{List of low-noise precision bipolar operational amplifiers with typical performance properties.}
    \label{tab:overview_bipolar_op-amps}
\end{table}

The low value of the sense resistor makes a bipolar op-amp the preferred choice, because they have a very low voltage noise and their current noise and input bias current do not interfere with such a low value resistor. While a discrete solution using matched JFETs or bipolar transistors may push the input noise even lower, the temperature stability, circuit complexity and again the size speaks against this option, so the discussion will be limited to integrated solutions only.

To find a reference point for the choice of op-amp, the thermal noise of the sense resistor must be looked at. A \qty{28}{\ohm} sense resistor has a thermal noise of
\begin{equation*}
    e_n\left(R = \qty{28}{\ohm} T = \qty{23}{\celsius}\right) = \qty[power-half-as-sqrt]{0.67}{\nV \per \Hz\tothe{0.5}} \,.
\end{equation*}

This means that even the lowest noise op-amp from table \ref{tab:overview_bipolar_op-amps} dominates the wideband-noise. With this said, the component choices made in literate can be discussed. The \device{AD8671} chosen by \citeauthor{laser_driver_digital} \cite{laser_driver_digital} only makes sense, because they have chosen a very large filter resistor $R_{f}$ of $2 \times \qty{3}{\kilo\ohm}$. The \device{ADA4004} was used by Moglabs in the \device{DLC-102}, again likely due to the high values of $R_{f}$ used. The \device{ADA4898} might seem like a good choice at first sight but the very limited (in terms of precision op-amps) open-loop gain of \qty{0.14}{\V \per \uV} makes this op-amp an inexpensive, but poor choice. The final choice is between the \device{AD797} and the \device{LT1028}, both op-amps have very similar specifications, but there is a peculiarity in the datasheet of the \device{LT1028} \cite{datasheet_LT1028}. The \textit{High Frequency Voltage Noise
vs Frequency} plot shows a noise bump at around \qty{400}{\kHz}. The original application of the \device{LT1x28} is for audio frequencies, which is well below that bump, but in this case poses a problem. The publication by \citeauthor{libbrecht_hall} \cite{libbrecht_hall} also blames the \device{LT1028} for the noise peak around \qty{400}{\kHz}. Moreover, this peak is the reason why \citeauthor{laser_driver_mosfet_noise} \cite{laser_driver_mosfet_noise} found the \device{LT1028} to have a higher integrated noise than the \device{AD797}. Additionally to the superior noise performance, the \device{AD797} (B-grade) has excellent specifications overall. The open-loop gain is between \qtyrange[range-units = single]{2}{20}{\V \per \uV}, the supply rejection is greater than \qty{1}{\uV \per \V}, the bias current is almost constant between \qtyrange[range-units = single]{20}{100}{\celsius} and the unity gain bandwidth is around \qty{10}{\MHz}. Finally it does have a high output drive capability of \qty{50}{\mA}, which allows to drive large MOSFETs. These features make the \device{AD797} the ideal op-amp for low-value sense resistors, although it puts limits on the maximum filter resistor to limit the low frequency current noise contribution.

Finally, the choice of MOSFETs can be discussed. As it was shown in section \ref{sec:mosfet_current_source} in equation \ref{eqn:mosfet_saturation}, the channel length modulation plays an important role in increasing the channel conductance $g_{DS}$ and limiting the output impedance. To reduce the channel length modulation a longer channel is preferred. Manufacturers do not give these numbers, nor information on the manufacturing process. Older technologies like the planar (lateral) FET are better suited for operating in the saturation region than the modern trench (vertical) FET. Trench MOSFETs are geared towards a low on-state resistance $R_{DS,on}$, which is important for switching applications, but their lower resistance comes from a shorter channel. One of the few planar MOSFETs still available on the market is the HEXFET, which was designed for switching applications, but proves useful nonetheless. High voltage MOSFETs also have longer channels than low voltage MOSFETs, so searching for MOSFETs that are rated for \qtyrange[range-units = single]{60}{100}{\V} or more can narrow down the candidates. While the output impedance is a factor worth keeping in mind, the most important aspect is whether the MOSFET can drive the load regarding the compliance voltage. To outline the problem, one can refer again to the example parameters from table \ref{tab:current_source_parameters}.

Assuming a supply voltage of \qty{15}{\V} and the \device{AD797} op-amp, the current source supply voltage $V_{sup}$ is limited to about \qtyrange[range-units = single]{11}{12}{\V}, because the \device{AD797} is no rail-to-rail op-amp and its output only swings within \qty{3}{\V} of the rail (minimum) and the input is limited to within \qty{2}{\V} of the rail (minimum). Considering the maximum $V_{ref}$ at full output of \qty{7}{\V} and a load voltage of \qty{3}{\V} in case of the \device{L785H1} \cite{datasheet_thorlabs_780nm} laser diode used as an example in this section leaves only
\begin{equation}
    V_{DS,min} = V_{sup} - V_{ref} - V_{load} = \qty{12}{\V} - \qty{7}{\V} - \qty{3}{\V} = \qty{2}{\V} \label{eqn:minimum_mosfet_vds}
\end{equation}
for the MOSFET -- a serious challenge.

To find a suitable MOSFET, one has to consult the \textit{Typical Output Characteristics} graph in the datasheet. By using the maximum output current specification it is possible to estimate the minimum drain-source voltage $V_{DS}$ to keep the MOSFET in saturation at the given maximum output current. This again narrows down the list of candidates.

The final aspect is the capacitive nature of the MOSFET gate. This property was briefly brushed in appendix \ref{sec:mosfet_noise} and the parasitic capacitances can be found in figure \ref{fig:mosfet_parasitic_capacitors}. The \device{AD797} can drive fairly large capacitive loads and several hundred \unit{\pF} are possible. It is best to keep the input capacitance $C_{iss}$ below \qty{500}{\pF}. The output impedance of the \device{AD797}, is about \qty{10}{\ohm} at \qty{1}{\MHz} and rising by an order of magnitude at \qty{10}{\MHz}. \qty{500}{\pF} results in an impedance of around \qty{300}{\ohm} dropping by an order of magnitude at \qty{10}{\MHz}, so keeping the capacitance low, allows for a higher bandwidth of the current source.

Using these guidelines, searching a MOSFET across a lot of manufactures can still be tedious, but the distributor Digikey, for example, allows filtering and sorting by voltage and input capacitance. The following MOSFETs in table are given as an example and can be chosen for their respective current ranges.
\begin{table}[ht]
    \centering
    \begin{tabular}{llll}
        \toprule
        MOSFET& Maximum $V_{DS}$& Input capacitance $C_{iss}$ & Current range \\
        \midrule
        \device{IRF9610} & \qty{200}{\V}& \qty{170}{\pF} & \qtyrange[range-units = single]{100}{250}{\mA}\\
        \device{IRF9Z10} & \qty{50}{\V}& \qty{270}{\pF} & \qtyrange[range-units = single]{250}{500}{\mA}\\
        \device{IRF9Z14} & \qty{60}{\V}& \qty{270}{\pF} & \qtyrange[range-units = single]{250}{500}{\mA}\\
        \bottomrule
    \end{tabular}
    \caption{Example MOSFETs for a current source and recommended current ranges.}
    \label{tab:example_mosfet_selection}
\end{table}

The current range of the MOSFETs in table \ref{tab:example_mosfet_selection} is given based on the datasheet, making sure that the MOSFET can be biased into saturation for the estimated minimum $V_{DS}$ according to \ref{eqn:minimum_mosfet_vds}. The \device{IRF9Z10} is a lower voltage version of the \device{IRF9Z14} and the \device{IRF9Z14} should be preferred if available. Those MOSFETs starting with \textit{IRF} are all HEXFETs formerly made by International Rectifier, whose MOSFET business was bought by Vishay in 2007.

To summarize the component selection. The ADI \device{AD797} op-amp is a superior choice among the op-amps compared and the recommended choice for being low noise with enormous gain, high bandwidth and a strong drive current. A MOSFET based on an older process, to ensure stable performance in saturation, like the Vishay \device{IRF9Z14} is a good choice for medium power applications. The reference of choice is a a buried Zener diode, because of their low noise and high stability. The ADI \device{LM399} or \device{ADR1399} is the most economical choice. Regarding the sense resistor, it must be able to dissipate up to $\qty{500}{\mA} \cdot \qty{7}{\V} = \qty{3.5}{\W}$ with minimal drift, making the Vishay \device{VPR221Z} a very good choice.

\subsection{Current Source Example Parameters}%
\label{sec:current_source_summary}
Throughout this section, example calculations are performed to give the reader an idea of real-life parameters that can be applied with the theoretical models. These parameters are summarised in table \ref{tab:current_source_parameters}, including their origin.
\begin{table}[ht]
    \centering
    \begin{tabular}{lll}
        \toprule
        Parameter& Value& Source \\
        \midrule
        MOSFET drain current $I_D$ & \qty{250}{\mA} & \device{L785H1} \cite{datasheet_thorlabs_780nm}\\
        MOSFET $\kappa$ & \qty[per-mode=power]{0.813}{\ampere \per \square\volt} & \device{IRF9610} SPICE model \cite{irf9610_spice}\\
        MOSFET channel length modulation $\lambda$ & \qty[per-mode=power]{4}{\per \milli \volt} & \device{IRF9610} SPICE model \cite{irf9610_spice}\\
        MOSFET source voltage & \qtyrange{3.5}{4}{\V} & section \ref{sec:component_selection}\\
        Source/Sense Resistor $R_s$ & \qty{14}{\ohm} or \qty{28}{\ohm} & section \ref{sec:component_selection}\\
        Op-amp differential input impedance $R_{id}$ & \qty{7.5}{\kilo\ohm} & \device{AD797} \cite{datasheet_AD797}\\
        Op-amp open-loop gain $A_{ol}$ & \qty[per-mode=power]{2}{\volt \per \uV} & \device{AD797} \cite{datasheet_AD797}\\
        Op-amp gain bandwidth product $GBP$ & \qty{10}{\MHz} & \device{AD797} \cite{datasheet_AD797}\\
        \bottomrule
    \end{tabular}
    \caption{Parameters used throughout this section and their sources.}
    \label{tab:current_source_parameters}
\end{table}

\subsection{Howland Current Pump}%
\label{sec:howland_current_source}
This section will discuss two versions of the popular Howland current source (HCS) \cite{howland_current_source}. Both the traditional Howland current pump and the so-called \textit{improved} Howland current pump which is similar, but changes some of its properties for good and bad. Both are bidirectional current sources and can be used for high frequency current modulation or the generation of a precision floating current source. The Howland current source is most useful for small currents of up to a few \unit{mA}, because its output impedance is limited and mostly depends on the resistors surrounding the op-amp as will be shown below. The design discussed here aims at an output current gain of \qty{1}{\mA \per \V}, which has proven useful for diode lasers over time. The full circuit is shown in figure \ref{fig:howland_current_source}. It is shown as the improved Howland current source but can be configured as the classic Howland current source by shorting out $R_{2a}$.
\begin{figure}[ht]
    \centering
    %\scalebox{0.5}{%
        \import{figures/}{howland_current_source.tex}
    %} % scalebox
    \caption{The Howland current source. Using $R_{2a} = \qty{0}{\ohm}$ is the classic version, while $R_{2a} \neq \qty{0}{\ohm}$ is the \textit{improved} version. $C_c$ is a compensation capacitor to enhance stability.}
    \label{fig:howland_current_source}
\end{figure}

As can be seen, these two current sources designs require a set of either \num{4} or \num{5} resistors depending on the desired configuration as the classic Howland current source does not need $R_{2a}$. The output impedance of the Howland current source is derived in appendix \ref{sec:appendix_howland_current_source} and the result is equation \ref{eqn:appendix_howland_output_impedance_resistors_gain} repeated here,
\begin{equation}
    R_{o,m,A} = \left(\frac{R^2 - R_{2a}^2}{R}\right) \frac{\left(A R + R - \epsilon + 1\right)}{A \left(R + \epsilon - 1\right) + 2 R - 2 \epsilon + 2}\,,\label{eqn:howland_output_impedance_resistors_gain}
\end{equation}
where $A$ is the frequency dependent gain of the op-amp as discussed in equation \ref{eqn:op-amp_gain}, $R = R_1 = R_2 = R_3 = R_4$ and $\epsilon$ is the resistor mismatch factor introduced in \ref{eqn:howland_mismatching_factor}. The mismatch factor can be applied for worst-case calculations for a given resistor tolerance $T$ using equation \ref{eqn:howland_current_source_worst_case_mismatch}.

By using the formula above, it is possible to calculate the output impedance for a modulation current source of a laser driver. Since the the output impedance of the \textit{improved} Howland current source is worse than that of the classic Howland current source according to equation \ref{eqn:howland_output_impedance_resistors_gain}, the discussion will begin here. Experience has shown that an input sensitivity of \qty{1}{\mA \per \V} is a reasonable choice, giving enough headroom to steer the laser and apply a modulation as mentioned above. The resistor value can be calculated from this input sensitivity to be \qty{1}{\kilo\ohm}. To reduce manual labor it is recommended to use an array of \num{4} matched resistors like the Vishay \device{MORN} \cite{datasheet_MORN} or the ADI \device{LT5400} \cite{datasheet_LT5400}, both quad resistor arrays offering a ratio matching between \qty{0.5}{\percent} and \qty{0.01}{\percent}. Commonly available is the \device{MORN} array with \qty{0.05}{\percent} tolerance and the \device{LT5400} with \qty{0.01}{\percent} tolerance and even \qty{0.005}{\percent} between neighbouring resistor pairs. For such arrays and an amplifier with $A=\num{e7}$, the worst case output impedance can be calculated as
\begin{align*}
    R_{out}(R=\qty{1}{\kilo\ohm}, T=\qty{0.05}{\percent}) &\approx \qty{499}{\kilo\ohm} & R_{out}(R=\qty{1}{\kilo\ohm}, T=\qty{0.01}{\percent}) &\approx \qty{2.5}{\mega\ohm} \,.
\end{align*}

$R_{out}$ is proportional to $\frac{1}{T}$, since the amplifier gain is very high and the tolerances are small. The upper limit imposed by the op-amp gain can be calculated according to equation \ref{eqn:howland_output_impedance_loop_gain} as
\begin{equation*}
    R_{out, max} (A=\qty{10}{\V \per \uV}) \approx \qty{2.5}{\giga\ohm}\,.
\end{equation*}

Comparing these values with the requirements from specification \ref{lst:dgDrive_specs_electrical}, $R_{out} \qty{\geq 7.5}{\mega \ohm}$, it is obvious that the Howland current source can meet them, but not without extra trimming or a tighter matched resistor array (T \qty{<= 0.01}{\percent}) which is hard to come by and expensive. In order to get a better understanding of the output impedance and to show the effect of different tolerances, a Monte Carlo simulation was conducted. The parameters assumed a Gaussian distribution of the resistor values with $3 \sigma = T$ due to the tolerance. The Monte Carlo simulation was then run \num{e5} times and the output impedance was extracted for each run. Finally, the absolute value (there are positive and negative values, see below) was binned into \num{100} bins of \qty{500}{\kilo\ohm} width and plotted as a histogram.

The LTSpice simulation file can be found at \external{source/spice/howland\_current\_source\_ideal.asc} and the result, limited to bin values $\qty{\leq 50}{\mega\ohm}$ for better readability, is shown in figure \ref{fig:ltpsice_howland_mc_output_impedance}.
\begin{figure}[ht]
    \centering
    \input{images/ltspice_howland_mc_output_impedance.pgf}% data/plot_ltspice_monte-carlo.py
    \caption{Histogram of the Monte Carlo simulations of the output impedance for the classic Howland current source with different resistor tolerances, assuming $A = \infty$, $3 \sigma = T$, $R=\qty{1}{\kilo\ohm}$. The improved Howland current sources uses $R = \qty{10}{\kilo\ohm}$ and $R_{2b} = \qty{1}{\kilo\ohm}$. The number of simulation runs was \num{e5}.}
    \label{fig:ltpsice_howland_mc_output_impedance}
\end{figure}

Figure \ref{fig:ltpsice_howland_mc_output_impedance} gives a more complete picture of the expected output impedance distribution when implementing an (improved) Howland current source without trimming. Do note, that figure \ref{fig:ltpsice_howland_mc_output_impedance} only shows the absolute value of the output impedance. The impedance can be either positive or negative, depending on the resistor ratios. The probability is evenly distributed between negative and positive impedance, producing a left-handed negative copy of the plot in figure \ref{fig:ltpsice_howland_mc_output_impedance}. For the purpose of improving readability of the plot, only the absolute value is plotted. A negative impedance means that with increasing load voltage, the output current increases as well. For the purpose of of driving laser diodes, this distinction is irrelevant as the output impedance mostly determines the noise immunity.

First, the basic Howland current source is discussed. Using \qty{0.05}{\percent} tolerance resistors the lower limit of \qty{500}{\kilo\ohm} can be easily identified as the leftmost bin of the histogram that contains a non-zero number of counts. The maximum probability is reached at around \qty{2}{\mega\ohm}. Integrating over the probability density of the values gives a \qty{31.5}{\percent} probability to end up with an output impedance of at least \qty{7.5}{\mega\ohm}. This is not nearly enough to meet the specificaction. Going to \qty{0.01}{\percent} tolerance resistors, the lower limit of \qty{2.5}{\mega\ohm} can again be identified by the absence of counts in the lower bins. The maximum probability is around \qtyrange{10}{20}{\mega\ohm} and the chance of getting an output impedance of more than the targeted \qty{7.5}{\mega\ohm} is \qty{95.6}{\percent}, a number far closer to the desired $3\sigma$ value (\qty{99.7}{\percent}) and can be considered acceptable given the high cost of acquiring a resistor array with better specifications.

The improved Howland current source was simulated with $R=\qty{10}{\kilo\ohm}$ and $R_{2b} = \qty{1}{\kilo\ohm}$. It is therefore expected that according to equation \ref{eqn:improved_howland_output_impedance_equal_resistors} the minimum output impedance is about a factor of two higher than compared to the basic Howland current source. Looking at figure \ref{fig:ltpsice_howland_mc_output_impedance}, this assumption holds. For a resistor tolerance of $T=\qty{0.05}{\percent}$, the minimum bin is \qty{1.5}{\mega\ohm}, with the maximum probability at \qty{4.5}{\mega\ohm}. The chance of having an output impedance of more than \qty{7.5}{\mega\ohm} is \qty{79.7}{\percent}, which is quite an improvement over the \qty{31.5}{\percent} of the basic current source, but still not usable. Using $T=\qty{0.01}{\percent}$ resistors, the improved Howland current source has a minimum impedance of \qty{8.5}{\mega\ohm}, which is sufficient to meet the requirements and the maximum is between \qty{20}{\giga\ohm} and \qty{30}{\giga\ohm}. The results are summarised again in table \ref{tab:howland_current_source_summary}.
\begin{table}[hb]
    \centering
    \begin{tabular}{lccc}
        \toprule
        Configuration& Tolerance $T$ & $R_{out, min}$ Bin& P($R_{out} \geq \qty{7.5}{\mega\ohm}$) \\
        \midrule
        \multirow{2}{*}{HCS} & \qty{0.05}{\percent}& \qty{500}{\kilo\ohm}& \qty{31.5}{\percent}\\
        & \qty{0.01}{\percent} & \qty{4}{\mega\ohm}& \qty{95.6}{\percent} \\
        \multirow{2}{*}{Improved HCS} & \qty{0.05}{\percent}& \qty{1.5}{\mega\ohm}& \qty{79.7}{\percent} \\
        & \qty{0.01}{\percent}& \qty{8.5}{\mega\ohm}& $\qty{100}{\percent}$\\
        \bottomrule
    \end{tabular}
    \caption{Summary of the Monte Carlo simulations for different current source configurations.}
    \label{tab:howland_current_source_summary}
\end{table}

In comparison to the MOSFET based precision current source, for which an output impedance of \qty{40}{\giga\ohm} was calculated in equation \ref{eqn:current_source_output_impedance}, the Howland current source is the weak link, when both are combined in a laser driver. To improve this situation the resistors can either be trimmed or a JFET or MOSFET cascode can be added to improve the output impedance. Do note a cascode is not bidirectional though. Several trimming options were discussed by \citeauthor{howland_pease} \cite{howland_pease}, but trimming at this level will prove difficult. The desired trim for the classic Howland current source and $R=\qty{1}{\kilo\ohm}$ can be calculated as
\begin{align*}
    R_{out} &\approx \frac{R}{\epsilon} \approx \frac{R}{4T} \geq \qty{7.5}{\mega\ohm}\\
    \Rightarrow T &\approx \qty{33}{\micro\ohm \per\ohm}\,.
\end{align*}

The final aspect that needs to be discussed is the compliance voltage. The compliance voltage of the HCS and improved HCS was calculated in equations \ref{eqn:howland_current_compliance_voltage} and \ref{eqn:improved_howland_current_compliance_voltage} as
\begin{align*}
    V_{c, HCS} &\leq \frac{1}{2} V_{o,max} & V_{c, iHCS} &\leq V_{o,max} - V_{in}
\end{align*}

From those equations it can be seen that there is a significant difference between the the basic Howland current source and the improved HCS. The compliance voltage of the basic Howland current source only depends on the maximum output voltage of the op amp and hence the supply voltage. The improved HCS depends on both the supply voltage and the input voltage. This makes it unsuitable for a laser driver modulation source, because the maximum modulation input typically is roughly the same voltage as the maximum output of a typical op-amp, since the laser driver is modulated or steered by another box using the same power rail. The improved Howland current source can only by employed in situations, where the input is well defined and a lot lower than the power supply rail of the HCS op-amp. For a laser driver modulation current source, the basic HCS is more suitable since it is independent of the input. To achieve a high output impedance, a high quality matched array is necessary, or alternatively, a difference amplifier like the ADI \device{LT1997} \cite{datasheet_LT1997} can be used. Those amplifiers contain both the amplifier and a resistor network that is tightly matched. In case of the \device{LT1997}, the matching is better than \qty{60}{\micro\ohm \per \ohm}. Unfortunately, the choice of resistor values in this case is even more limited than the range of resistor arrays.

The problem of a limited choice of resistor values as arrays or integrated resistors can be addressed by a circuit that combines both the basic HCS and the improved HCS. This circuit is shown in figure \ref{fig:buffered_howland_current_source}.
\begin{figure}[hb]
    \centering
    %\scalebox{0.5}{%
        \import{figures/}{buffered_howland_current_source.tex}
    %} % scalebox
    \caption{A buffered Howland current source combining the improved HCS and the basic HCS.}
    \label{fig:buffered_howland_current_source}
\end{figure}

By adding another op-amp $R_2a$ is now connected to a low impedance node, just like in case of the basic Howland current source. This leads to the same formulae regarding the output impedance and the compliance voltage as the improved HCS, but removes the matching requirement of $R_2a$ and $R_2b$. The full derivation can be found in the Python notebook at \external{data/simulations/howland\_current\_source.ipynb} as part of the online supplemental material \cite{supplemental_material}.

With the buffered HCS it is possible to use either \num{4} matched resistors or a difference amplifier with integrated resistors. The single resistor $R_5$ is used to set the output current and can be chosen independently to result in
\begin{equation}
    I_{out} = \frac{V_{in}}{R_5}\,.
\end{equation}

% TODO: Investigation of Long-Term Drift of NTC Temperature Sensors with less than 1 mK Uncertainty