\section{Introduction} \label{intro}

RoboPol  \citep{Ramaprakash2019} is a four-channel optical polarimeter installed on the 1.3m telescope at the Skinakas Observatory in Crete, Greece that is primarily used for polarimetry of point sources in the R band. Its successor the Wide Area Linear Optical Polarimeter (WALOP) is under development at the Inter-University Center for Astronomy and Astrophysics (IUCAA) in Pune, India. Images contain artefacts resulting from dust patterns, cosmic ray hits, satellite trails and pixel bleeding contaminating information from celestial objects. With the increasing number of images taken every night from such instruments, it is necessary to automate the analysis of data. However, with humans taken out of the loop it is possible for artefacts to get misidentified as a source and be used in the analysis. This would lead to erroneous results so the detection of such artefacts is imperative.

Early work on detection of artefacts in astronomical images dates to the early 2000s when \cite{Storkey2004} used computer vision techniques such as the Hough Transform to detect linear artefacts like satellite trails, scratches, and diffraction spikes near bright stars. These methods were concerned with detection of linear features and highlighted some of the difficulties of using the Hough Transform  when dealing with light-density variations.

Later on the focus shifted to object identification followed by extracting features for the objects and then classification using these features to separate out artefacts with methods like decision trees and random forests \citep{donalek2008new}. Recent years have seen the compilation of terrestrial datasets like ImageNet consisting of a million labelled images with a thousand categories  such as human faces, digits, vehicles, flowers, animals etc. \cite{deng2009imagenet} followed by the development of deep learning libraries and models using these datasets e.g. the VGG16 architecture \cite{simonyan2014very}.

With deep learning it is possible to skip the sometimes subjective step of feature extraction and go straight to classification after obtaining a labeled dataset  \citep[see e.g][]{Cabrera2017,duev2019RB,duev2019deepstreaks}. This is at the cost of explainability, but with proper validation and test datasets, the results are still reliable. Additional ways to improve the robustness, and faster convergence using techniques like Mask R-CNN and linear scaling combined with normalization are discussed in some recent papers such as \cite{he2017mask,gonzalez2018supervised,burke2019deblending}.

Our task here is to classify objects in RoboPol images into stars and artefacts. RoboPol images contain reflections of bright stars due to the interface between the 2 Wollaston Prisms used in the instrument and it is this dominant class of artefacts that we target here. The interface between the Wollaston Prisms is shown in the diagram of the optical instrument design described in \cite{Ramaprakash2019}, and an example of the reflection artefact in Figure \ref{fig:mask}. The green box in the upper left quadrant shows two horizontally extended artefacts separated vertically. A few stars in the vicinity also got included in the box.

In this paper, we propose to solve the problem of artefact detection for RoboPol images using an appropriately designed Convolutional Neural Network (CNN).
In Section \ref{approach}, we introduce our approach to the problem of detecting artefacts in RoboPol images. We detail the implementation of our method including pre-processing steps, CNN architecture, and visualization of the output, and
in Section \ref{discussion1} we discuss our findings and future possibilities.

\section{Approach}\label{approach}

 \begin{figure}
    \includegraphics[scale=0.33]{artifacts_borders_copy.jpg}
    \caption{The region within the orange lines describes the restricted area from where we do not extract stars or artefacts. Similarly the regions outside the white lines are restricted areas. Note that the lines are exaggerated only for the purpose of representation.}
    \label{fig:mask}
\end{figure}

 The RoboPol database consists of tens of thousands of images taken between 2013 and 2019. We first generate a data-set containing stars and artefacts and then develop a CNN to perform the classification.
The following is an outline of our method:
\newline 1. Create training data for artefacts and stars from RoboPol images through manual labeling. This includes data for validation and testing. The manual labeling was done by visually inspecting about 100 images and recording the pixel coordinates of the artefacts.
\newline 2. Develop a CNN Architecture tuned through hyperparameter variation.
\newline 3. Train the model using training data obtained in step 1.
\newline 4. Validate the model using validation and testing data.
\newline 5. Implement the model to find artefacts in an arbitrary RoboPol image.

\subsection{Training Data for Artefacts and Stars}

Reflection artefacts and stars in RoboPol images have x and y extents from several pixels to a few tens of pixels. We chose a size of 64x64 pixels for our cutouts, with the artefacts and stars centered. Each star appears at four locations due to splitting of light from a single source within the instrument, with the locations lying at the vertices of a diamond. A detailed design implementation is available in \cite{Ramaprakash2019}.

For each image we generate a catalog of sources (including stars, reflection artefacts, and any other connected brightness peaks) using Sextractor\footnote{https:///sextractor.readthedocs.io}. We have roughly 10 artifacts and about 250 stars per image. The catalog comes with flags indicating various conditions such as saturation, proximity to another source, proximity to edge of image etc.\footnote{https://sextractor.readthedocs.io/en/latest/Flagging.html}.
About 70\% of the visually inspected artefacts had no error. This indicated that relying on just flags is not sufficient to separate artefacts.  To obtain training data of stars we make sure that from every image we extract stars of varying brightness and not just from a narrow brightness range. In each image, we chose this range to be 1 star per magnitude-bin for up to 5 magnitudes in each image. Likewise, our training data would contain an uniform distribution of magnitudes of brightness and ensure that we aren't biasing our neural network by providing training images from a limited magnitude range. We do not use all sources so that the sets of stars and artefacts can stay roughly equal, and hence balanced for the classification process.

Unlike typical astronomy images, RoboPol images contain a mask (Figure \ref{fig:mask}). More details about the mask can be found in \cite{Ramaprakash2019}.
The code repository and documentation are available at  https://github.com/delta-papa/Robopol-artifacts

\subsection{CNN Architecture}

\begin{figure*}
    \centering
    \includegraphics[width=0.99\textwidth]{model.jpg}
    \caption{The CNN model consists of 3 convolutional layers, 3 max-pooling layers, 2 fully connected layers and a dropout layer for regularization. In case of the convolutional and max-pooling layers, the size of the layer and the filter size used are also mentioned. For the dense layers, the total number of nodes are shown. The input is a 64x64 pixel PNG image.}
    \label{fig:architecture}
\end{figure*}

We follow the now standard image classification model developed by the Visual Graphics Group (VGG) at Oxford, UK \citep{simonyan2014very}. Our implementation uses 3 convolution layers, 3 max-pooling layers and 2 fully connected layers (see Figure \ref{fig:architecture}). The hidden layers are activated using a ReLU (Rectified Linear Unit) activation. Finally we use a sigmoid activation at the output layer. The first, second and third convolution layer consist of 32, 64 and 128 filters respectively each with a kernel size of 3x3 and stride length of 1. The max-pooling layers use a kernel size of 2x2 pixels. At the end of the 3rd max-pooling layer we use a dropout layer with a probability of dropping a node as 0.4 for regularization ensuring no single parameter of the neural network has a very high coefficient \citep{srivastava2014dropout}. The total number of trainable parameters in our configuraton are 2,452,993 and we use an Adam optimizer to perform back-propagation \citep{kingma2014adam}. The loss function used is a binary cross-entropy loss.

\subsection{Data Augmentation and Training}

While going through the images we saw that most of the artefacts are due to internal reflections and had a horizontal streak-like shape. For proper training of the neural network we need to generate a large data-set of training images. Therefore, to augment the number of images for training, we rotated the cutouts by 180 degrees so that the horizontal nature of the artefacts is preserved.

\begin{figure*}

    \centering
    \begin{subfigure}{0.45\textwidth}
        \centering
        \includegraphics[width=0.9\textwidth]{trainingstartile.jpg}
        \caption{Cutouts of Stars used as training images.}
    \end{subfigure}\hfill
    \begin{subfigure}{0.45\textwidth}
        \centering
        \includegraphics[width=0.9\textwidth]{trainingartifacttile.jpg}
        \caption{Cutouts of Artefacts used as training images.}
    \end{subfigure}
    \caption{Training images for stars (a) and artefacts (b).}
\end{figure*}

To split the data into training and validation we used a 80-20 ratio and shuffled the data-set randomly while splitting to avoid bias. For training, we had a total of  836 images of stars and 925 images of artefacts. The training data was augmented using the ImageDataGenerator class in the high-level Keras\footnote{www.keras.io} API of Python. We performed horizontal and vertical flipping, width and height shifts and shearing.
The shifts were applied to account for possible inaccuracies in centering of the samples.
The validation and training images were both normalized to [0,1] by dividing by 255, the maximum value of an 8 bit image.

\begin{figure}
    \includegraphics[width=0.45\textwidth]{training_validation_plots_model9.jpg}
    \caption{Training and Validation Accuracy for the final model.}
    \label{fig:accuracy}
\end{figure}

The hardware used for performing training was a 2.3 GHz Intel Core i5 processor and the total training time was 40 minutes.

\subsection{Training Performance}
The total number of images chosen for training the model was 1408 (80\% of the 1761 images), and the remaining 353 images were reserved for validation. Training data was used to update the parameters of the model while validation data was used to only evaluate the model's performance after every update. The batch size used was 4. A total of 100 epochs were used in the training and the steps per epoch was set to 1408/4 i.e. 352. We used a learning rate of 0.001. The training accuracy reached about 95\% while the validation accuracy was close to 96\% at the end of 100 epochs as seen in Figure \ref{fig:accuracy}.

\begin{figure}

    \includegraphics[scale=0.17]{confusion_matrix_3conv.jpg}
    \caption{Confusion Matrix for the Binary Classification performed using our model. The normalized numbers are given in brackets. Out of 353 images, 170 were True Positives, 175 were True Negatives, 5 False Positives and 3 False Negatives. The False Positive Rate is 3\% and False Negative Rate is 2\%.}
    \label{fig:confusionmatrix}
\end{figure}

Besides making small changes to the hyperparameters above, we also implemented a network with 2 and 4 convolution layers to see whether there was any advantage in using shallower (2 layers) or deeper (4 layers) CNNs. The training accuracy and validation accuracy reached about 90\% in the shallower network while it reached about 96\% in the deeper network. Although training and validation accuracy may be good indicators of the proper working of a CNN in a binary classification problem there are other important parameters we need to consider when the costs of misclassification are high. For example, we are interested in knowing the false positive rate (sources wrongly classified as stars), the false negative rate (sources falsely classified as artefacts), as also the precision (fraction of sources correctly classified) and recall (fraction of stars correctly classified).  These numbers are summarized through two metrics viz. F1 score and Matthew's Correlation coefficient, and indicate whether our model is working as expected or not.

We need our system to have a high precision and recall score and the F1 score summarizes the 2 scores by taking their harmonic mean.
The Matthew's Correlation Coefficient (MCC) is akin to a correlation coefficient measure between the predicted labels and true labels. A value of +1 indicates perfect positive correlation between the 2 quantities. Equations
\ref{prec}, \ref{recall}, \ref{f1} and \ref{mcc} give the formulae for Precision, Recall, F1 score and MCC respectively. Note that TP, FP, TN, FN stand for True Positives, False Positives, True Negatives and False Negatives respectively.

\begin{equation}\label{prec}
    Precision = \frac{TP}{TP + FP}
\end{equation}

\begin{equation}\label{recall}
    Recall = \frac{TP}{TP + FN}
\end{equation}

\begin{equation}\label{f1}
    F1\ score = \frac{2 \cdot precision\cdot recall}{precision+ recall}
\end{equation}

\begin{equation}\label{mcc}
    MCC = \frac{TP*TN - FP*FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}
\end{equation}

Figures \ref{fig:confusionmatrix} and \ref{fig:roc} show the confusion matrix and ROC curve respectively. The Confusion Matrix tells us the number of true positives, true negatives, false positives and false negatives. The Receiver Operating Characteristics (ROC) curve is a measure of the trade-off between the true positive rate and false positive rate for different values of the threshold used in the classifier. The threshold is a value between 0 and 1. If the probability of the source being a star is greater than the value of the threshold we classify the source as a star else as an artefact. Our threshold is set to 0.5. Ideally, we want the false positive rate to be 0 and true positive rate to be 1. Ideally we expect our ROC curve to have an area (under the curve) to be equal to 1. The Area Under the Curve (AUC) for our model is 0.996 while that for a random classifier is 0.5 as shown by the red dotted line. A zoom-in of Figure \ref{fig:roc} is shown in Figure \ref{fig:roczoom} which shows the values True Positive Rates and False Positive Rates at different thresholds. It shows that our threshold of 0.5 coincides with 0.75 indicating that our classifier can confidently achieve the same True Positive Rate and False Positive Rate at a higher threshold. Table \ref{table:comparison} compares the performance of our model with a shallower and deeper neural network. Based on this comparison, we chose the model with 3 convolution layers as the model with 4 layers actually shows a reduction in F1 score and MCC with a deeper and hence computationally more expensive network, possibly due to over-fitting. Figures \ref{fig:fn} and \ref{fig:fp} show some of the false negatives and false positives.

\begin{figure}

    \includegraphics[scale=0.2]{roc_conv3.jpg}
    \caption{ROC Curve for our final model.}
    \label{fig:roc}
\end{figure}

\begin{figure}
    \includegraphics[scale=0.2]{roc_conv3_zoom.jpg}
    \caption{Zoom in of ROC (Figure \ref{fig:roc}) to highlight the non-ideal area.}
    \label{fig:roczoom}
\end{figure}

\begin{figure}

    \includegraphics[scale=0.16]{fn_invert.jpg}
    \caption{Some of the False Negative classifications where stars were classified as artefacts.}
    \label{fig:fn}
\end{figure}

\begin{figure}
    \includegraphics[scale=0.16]{fp_invert.jpg}
    \caption{Some of the False Positive Classifications where artefacts were classified as stars. We see that the artefacts here didn't have a prominent streak-like shape. Also, the first artefact has a star near the center and the artefact is slightly away from the center.}
    \label{fig:fp}
\end{figure}

\begin{table}
\centering
\caption{Performance comparison of CNN models with different layers. Here FPR is False Positive Rate, FNR is False Negative Rate, MCC is Matthew's Correlation Coefficient.}
\label{table:comparison}
\begin{tabular}{lccc}
 \hline
Parameter &2 Layers &3 Layers& 4 Layers \\
 \hline
 Training Accuracy & 0.90 & 0.95 & 0.90 \\
 Validation Accuracy & 0.91 & 0.96 & 0.95 \\
 FPR & 0.07 & 0.03 & 0.02 \\
 FNR & 0.13 & 0.02 & 0.05 \\
 Precision & 0.88 & 0.98 & 0.95 \\
 Recall & 0.93 & 0.97 & 0.98 \\
 F1 Score & 0.90 & 0.98 & 0.96 \\
 MCC & 0.80 & 0.95 & 0.93\\
 \hline
\end{tabular}
\end{table}

\subsection{Testing and Implementation}

\begin{figure}
    \includegraphics[scale=0.3]{artifact_inver.jpg}
    \caption{Result of testing the classifier on images of artefacts. All artefacts have been classified correctly with high classification probabilities. The classification probability of a source being a star is shown on the top of each cutout.}
    \label{fig:arti_result}
\end{figure}

\begin{figure}
    \includegraphics[scale=0.30]{star_new.jpg}
    \caption{Result of testing the classifier on images of stars. The classification probability of a source being a star is shown on the top of each cutout.}
    \label{fig:star_result}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{flowchart.jpg}
    \caption{A flowchart showing the decision logic of the pipeline.}
    \label{fig:flowchart}
\end{figure}

 We used our model to test 100 randomly chosen images - distinct from the training set - from the 40,000 images of RoboPol taken during the years 2013 and 2014. Our goal was to classify all the sources in each of the 100 images into stars and artefacts and analyze the results. For each image, we obtained  a list of sources, their positions, instrumental magnitudes and extraction flag error by using Sextractor. A histogram of the predicted probability of the sources being stars is shown in Figure \ref{fig:hist1channel}. Out of ~91,000 sources, we have ~88,000 sources classified as stars (90 to 100\% probability of being a star) and 2500 sources classified as artefacts (0 to 10\% probability of being a star). The inset plot shows that in the remaining prediction probability range there are fewer than 10 objects in each bin of size 10\%. This means that $\sim0.05\%$ of the sources had probabilities in the range between 10 to 90\%.

 Figures \ref{fig:arti_result} and \ref{fig:star_result} show the results of classification on test images. Each image contains a single source with known label. At the top of each image is the probability of the source being a star. Sources in Figure \ref{fig:arti_result} are artefacts while those in Figure \ref{fig:star_result} are stars. In both categories, our classification rates are almost always above 90\%.

The implementation pipeline takes an input image and the corresponding Sextractor file as its arguments and produces a list of locations of the detected artefacts along with their location marked in the original image with an associated probability. A decision logic diagram is shown in Figure \ref{fig:flowchart}.

 \begin{figure}
     \includegraphics[scale=0.18]{histogram_random_test.jpg}
     \caption{Histogram of prediction probabilities for ~91,000 sources.}
     \label{fig:hist1channel}
 \end{figure}

\subsection{Visualization with Saliency Maps}\label{viz}
A saliency map helps us find the locations of the pixels in  input images which need to be changed the least to activate the output filter. This means we find the gradient of the input image with respect to the output score. To visualize a saliency map, positive gradients are chosen that would give us the locations of the pixels activating the output filter. In other words this gives us the location of the object of the relevant class in the input image. A saliency map thus gives us the salient features of the class-specific input image that maximize the class score. A detailed mathematical treatment can be found in \cite{journals/corr/SimonyanVZ13}. Figure \ref{fig:saliency} shows the saliency maps for 9 different input images. Observe that for the image in row 3, column 3, only the artefact is visualized in the saliency and not the star at the top right corner of the image. Sextractor's ellipticity measure alone is not sufficient to separate the artefacts.

\begin{figure}
    \includegraphics[scale=0.22]{saliency_combine.jpg}
    \caption{Saliency Map visualization (left) for 9 different input images (right) of artefacts. The Saliency map shows the positive gradients of the image with respect to the artefact class score and thus the locations of the artefacts in the image. These maps show that indeed the classifier is activated by the artefacts themselves and not their background. Figure plotted using Keras Visualization Toolkit \citep{raghakotkerasvis}.}
    \label{fig:saliency}
\end{figure}

\section{Discussion}\label{discussion1}
We used Convolutional Neural Networks to solve the problem of detecting artefacts in polarimetric images. Although the use of CNNs in astronomical image classification is not new, this is the first time that they have been used for detecting artefacts in polarimetric images. The efficiency of the method shows its suitability for use in upcoming polarimetry surveys such as the Polar Areas Stellar Imaging in Polarization High Accuracy Experiment \citep{pasiphae}, which will use the novel Wide Area Linear Optical Polarimeter (WALOP). Our implementation suggests that this method can be reliably used for detecting other kinds of artefacts as well given enough training data. The RoboPol instrument operates down to 16th magnitude in the R1 band. Figure \ref{fig:magscatter} shows that our Deep Learning model can classify stars down to 15.9 magnitude with a prediction probability better than 0.9. We have also plotted the signal to noise ratio (SNR) of the stars on a separate axis. We see that our model works up to SNRs of 15. Thus, our implementation works well with objects within the magnitude range RoboPol observes.

In this paper, we do not use the spatial correlation for stars appearing as a diamond pattern in RoboPol images. That is because the diamond structure in each image of RoboPol is specific to the RoboPol polarimeter design and wouldn't be present in a single image of future polarimeters such as WALOP.

In the RoboPol data-set, we had majority of artefacts due to scattering of light from off-axis stars at the interface of the Wollaston prisms \citep{Ramaprakash2019}. Our method demonstrates that a binary classifier trained on images of stars and artefacts can successfully differentiate between them.  Our training data does not contain enough examples for artefacts such as satellite trails or bleeding pixels and as a result deep learning them is non-trivial without aggressive data augmentation. There already exist methods to remove such artifacts. Out-of-Distribution detection networks \citep{huang2019out} can also be used to detect such infrequent outliers. The final pipeline can incorporate such methods to deliver artefact-free products.

\begin{figure}
    \includegraphics[scale=0.36]{snr.jpg}
    \caption{Scatter plot of prediction probability and signal to noise ratio (SNR) of visually inspected stars from a single image versus the apparent R1 magnitude.  The 42 stars span a magnitude range from 12 to 16.}
    \label{fig:magscatter}
\end{figure}

\section*{Acknowledgements}
The work has been funded by the National Science Foundation under the NSF grant (161547). AM acknowledges support from the NSF (1640818, AST-1815034) and IUSSTF (JC-001/2017). KT acknowledges support  from  the  European Research Council  under  the  European  Union’s Horizon  2020 research  and  innovation  program,  under  grant  agreement  No771282.

\bibliographystyle{mnras}


\section*{ACKNOWLEDGMENTS}
This research was sponsored by the Army Research Office (ARO) under grant numbers W911NF-17-1-0274 and W911NF-12-1-0607. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the ARO or the US Government. The US Government is authorized to reproduce and distribute reprints for government purposes notwithstanding any copyright notation herein. Development and maintenance of the growth facilities used for fabricating samples is supported by DOE (DE- FG02-03ER46028). B.D. acknowledges funding from NSERC. W.A.C. acknowledges funding from NSERC, CIFAR, and FRQNT. We acknowledge the use of facilities supported by NSF through the University of Wisconsin-Madison MRSEC (DMR-1121288). We acknowledge useful discussions with Matthew G. House, and the members of the Vandersypen group, and technical assistance by R. Schouten and R. Vermeulen.

\section{Bayesian inference \label{sec:bayesianInference}}

In this section, we give an algorithm to efficiently perform optimal Bayesian inference of the logical qubit state using the empirically determined statistics of the repetitive readout.

\subsection{Repetitive readout}

Suppose that the logical qubit is repetitively read out $N$ times, with each repetition having a duration $\delta t_{\textrm{rep}}$. We may consider the state of the logical qubit at the discrete times:
\begin{align}
	t_k = k\,\delta t_{\textrm{rep}}, \;\;\; k=0,1,...,{N},
\end{align}
with the $k^{\textrm{th}}$ repetition taking place in the interval $\left[t_{k-1},t_k\right]$. Here, the coherence of the logical qubit in the computational basis plays no role in the statistics of the measurement. For the present purposes, we may therefore model the time evolution of the logical qubit classically. The classical state of the logical qubit at time $t_k$ is labeled $x_k$. The state evolution of the logical qubit up to $k \leq N$ repetitions follows the stochastic time series:
\begin{align}
	\mathbf{X}_{k} = \left\{ x_0, x_1, \dots, x_k \right\}.
\end{align}
In reality, each individual readout is noisy. Thus, the state $x_k$ at each repetition is not directly recorded. Instead, a noisy observation $\mathcal{O}_k$ is recorded. This gives the observation sequence:
\begin{align}
	\mathbf{O}_k = \left\{\mathcal{O}_0, \mathcal{O}_1, \dots ,\mathcal{O}_{k-1}\right\}.
\end{align}
Note that in general, the observations $\mathcal{O}_k$ need not be scalar.

\subsection{Maximum-likelihood decision}

Our goal is to infer the most likely initial state of the logical qubit from the sequence $\mathbf{O}_k$ of noisy repetitive readout outcomes. This is most easily done by calculating the posterior probability ratio:
\begin{align}
	\frac{P(x_0 = 1|\mathbf{O}_k)}{P(x_0 = 0|\mathbf{O}_k)} = \frac{P(\mathbf{O}_k|x_0 = 1)}{P(\mathbf{O}_k|x_0 = 0)} \frac{P(x_0 = 1)}{P(x_0 = 0)}.
\end{align}
The initial state is most likely $1$ if the ratio is larger than unity, and it is most likely $0$ if the ratio is smaller than unity. In the absence of prior information on the logical qubit state, $P(x_0 = 0) = P(x_0 = 1) = 1/2$, the above {\it maximum a posteriori} decision~\cite{kay1998} reduces to calculating the log-likelihood ratio:
\begin{align}
	\lambda^{\textrm{log}}_k = \ln \frac{P(\mathbf{O}_k|x_0 = 1)}{P(\mathbf{O}_k|x_0 = 0)}. \label{eq:logicalLogLikelihoodRatio}
\end{align}
The initial state is now most likely $1$ if $\lambda^{\textrm{log}}_k > 0$, and it is most likely $0$ if $\lambda^{\textrm{log}}_k < 0$. This results in a maximum-likelihood decision~\cite{kay1998}.

\subsection{Logical readout error rate \label{eq:logicalErrorRate}}

The average logical readout error rate $\epsilon$ is given by
\begin{align}
	\epsilon^{\textrm{log}} = \frac{1}{2}\left(\epsilon^{\textrm{log}}_1 + \epsilon^{\textrm{log}}_0\right).
\end{align}
Here, $\epsilon^{\textrm{log}}_1$ and $\epsilon^{\textrm{log}}_0$ are the error rates conditioned on preparation of the logical qubit in states $1$ and $0$ at time $t_0$, respectively. These are given by:
\begin{align}
\begin{split}
	\epsilon^{\textrm{log}}_1 = P(\lambda^{\textrm{log}}_k < 0| x_0 = 1), \\
	\epsilon^{\textrm{log}}_0 = P(\lambda^{\textrm{log}}_k > 0| x_0 = 0).
\end{split}
\end{align}
An experimental estimate of the logical readout error rate $\epsilon^{\textrm{log}}_1$ ($\epsilon^{\textrm{log}}_0$) is obtained by preparing the logical qubit in state $\ket{1}$ ($\ket{0}$) $10^4$ times, measuring the record $\mathbf{O}_k$ and calculating $\lambda^{\textrm{log}}_k$ for each attempt, and counting the number of times where $\lambda_k^{\textrm{log}} < 0$ ($\lambda_k^{\textrm{log}} > 0$). Finally, we note that the average readout fidelity quoted in the main text is defined as $F^{\textrm{log}} = 1 - \epsilon^{\textrm{log}}$, while the logical visibility is $\mathcal{V}^{\textrm{log}} = 1 - 2\epsilon^{\textrm{log}}$.

\subsection{Calculating the log-likelihood ratio}

To calculate $\lambda^{\textrm{log}}_k$, the statistics of the observations $\mathbf{O}_k$ given the initial state of the logical qubit must be known. More precisely, the probability distributions $P(\mathbf{O}_k|x_0 = 1)$ and $P(\mathbf{O}_k|x_0 = 0)$ appearing in Eq.~\eqref{eq:logicalLogLikelihoodRatio} must be calculated. In the following, we show how to model and calculate these distributions for the repetitive readout using the theory of hidden Markov models~\cite{hann2018,nakajima2019,elder2019,bultink2019,yoneda2019}. The hidden Markov model provides a direct connection between the single-repetition probability distributions, $P(\mathcal{O}_k|x_k)$, and the probability distributions for the full measurement record, $P(\mathbf{O}_k|x_0)$, accounting for the dynamics of the logical qubit.  From $P(\mathbf{O}_k|x_0)$, the log-likelihood ratio can be evaluated directly using Eq.~\eqref{eq:logicalLogLikelihoodRatio}.  As described in Sec.~\ref{sec:repetitiveStatistics}, following Eq.~\eqref{eq:softNoiseVector}, here we have determined $P(\mathcal{O}_k|x_k)$ empirically. However, note that these distributions can also be found from an appropriate dynamical model of the readout. See for instance Refs.~\cite{gambetta2007,gammelmark2014,ng2014,danjou2014,danjou2016} and others, where hidden Markov models are used to determine $P(\mathcal{O}_k|x_k)$ at the level of a single repetition.

\subsubsection{Hidden Markov models}

As discussed in the main text, the logical qubit state evolves during the measurement via spin relaxation on a time scale $T_1$. Such a process is Markovian in the sense that the statistics of the state at time $t_{k+1}$ are fully determined by the state at time $t_{k}$:
\begin{align}
	P(x_{k+1}|\mathbf{X}_k) = P(x_{k+1}|x_k). \label{eq:markovCondition}
\end{align}
Because the ancilla qubit is reinitialized after each repetition, the noisy observations $\mathcal{O}_k$ are independent from each other and depend only on the state of the logical qubit at time $t_k$:
\begin{align}
	P(\mathbf{O}_k|\mathbf{X}_k) =  \prod_{l=0}^{k-1} P(\mathcal{O}_l|x_l). \label{eq:whiteNoise}
\end{align}
In other words, the observation noise is white. Finally, prior knowledge about the system state at each time is specified by the prior probability distribution $P(x_k)$ for the state $x_k$ at each time $t_k$. The set $\left\{P(x_{k+1}|x_k), P(\mathcal{O}_k|x_k), P(x_k)\right\}$ defines a `Hidden Markov Model'.

\subsubsection{Forward filtering}

For the hidden Markov models discussed above, the log-likelihood ratio may be calculated with the help of an iterative forward filtering algorithm for the logical qubit state~\cite{gambetta2007,gammelmark2013,gammelmark2014,danjou2014,danjou2016,hann2018,nakajima2019,elder2019,bultink2019,yoneda2019}. Forward filtering consists in calculating the probability distribution of the logical qubit state $x_k$ at time $t_k$ given all previous observations. We denote this probability distribution as
\begin{align}
	\varrho_k(x_k) = P(x_k|\mathbf{O}_{k}).
\end{align}
Using Bayes' rule, the distribution $\varrho_k(x_k)$ may be rewritten as
\begin{align}
	\varrho_k(x_k) = \frac{P(x_k,\mathbf{O}_{k})}{P(\mathbf{O}_{k})} = \frac{P(x_k,\mathbf{O}_{k})}{\sum_y P(y,\mathbf{O}_{k})} = \frac{\ell_k(x_k)}{\sum_y \ell_k(y)}. \label{eq:likelihoodToState}
\end{align}
Here,
\begin{align}
	\ell_k(x_k) = P(x_k,\mathbf{O}_k)
\end{align}
is the joint probability distribution of the state $x_k$ and of the previous observations. One advantage of calculating $\ell_k(x_k)$ instead of $\varrho_k(x_k)$ is that $\ell_k(x_k)$ obeys a linear recurrence relation while $\varrho_k(x_k)$ obeys a non-linear recurrence relation (see Sec.~\ref{sec:filteringEquations}).

We note that the denominator in Eq.~\eqref{eq:likelihoodToState} is the likelihood function
\begin{align}
	\mathcal{L}_k = P(\mathbf{O}_k) = \sum_{x_k} \ell_k(x_k). \label{eq:likelihoodFunction}
\end{align}
Eq.~\eqref{eq:likelihoodToState} now takes the form
\begin{align}
	\ell_k(x_k) = \varrho_k(x_k) \times \mathcal{L}_k. \label{eq:likelihoodToState2}
\end{align}
It is convenient to introduce the vector notation
\begin{align}
\begin{split}
	\ketR{\varrho_k} = \sum_x \varrho_k(x) \ketR{x}, \;\;\;\ketR{\ell_k} = \sum_x \ell_k(x) \ketR{x},
\end{split}
\end{align}
where $\left\{\ketR{x}\right\}$ is a set of basis vectors representing the classical logical qubit states $x$. Eq.~\eqref{eq:likelihoodToState2} is then rewritten as
\begin{align}
	\ketR{\ell_k} = \mathcal{L}_k  \times \ketR{\varrho_k}.\label{eq:likelihoodToStateVec}
\end{align}
The likelihood function, Eq.~\eqref{eq:likelihoodFunction}, may compactly be written as
\begin{align}
	\mathcal{L}_k = \Tr \ketR{\ell_k}. \label{eq:likelihoodFunctionVec}
\end{align}
Here, the trace $\Tr$ of a vector is the sum of its elements. Choosing the basis $\left\{\ketR{1},\ketR{0}\right\}$, the likelihood function may be calculated for the initial states $x_0 = 1$ and $x_0 = 0$ of the logical qubit by setting $\ketR{\varrho_0} = \ketR{\ell_0} = (1,0)^\textrm{T}$ and $\ketR{\varrho_0} = \ketR{\ell_0} = (0,1)^\textrm{T}$, respectively. This enables the maximum-likelihood decision discussed following Eq.~\eqref{eq:logicalLogLikelihoodRatio}.

\subsubsection{The filtering equations \label{sec:filteringEquations}}

For completeness, we now derive the recurrence relation for forward filtering of hidden Markov models. We note that
\begin{align}
\begin{split}
	\ell_{k+1}(x_{k+1}) = P(x_{k+1},\mathbf{O}_{k+1}) &= \sum_{x_k} P(x_{k+1},\mathbf{O}_{k+1}|x_k)P(x_k) \\
	&= \sum_{x_k} P(\mathbf{O}_{k+1}|x_{k+1},x_k)P(x_{k+1}|x_k)P(x_k).
\end{split}
\end{align}
Next we recall that 1) $\mathbf{O}_{k+1} = \left\{\mathcal{O}_1,\mathcal{O}_2,\dots,\mathcal{O}_k \right\}$ does not depend on $x_{k+1}$ by causality and that 2) the observation noise is white, Eq.~\eqref{eq:whiteNoise}. Thus, we have
\begin{align}
\begin{split}
	\ell_{k+1}(x_{k+1}) &= \sum_{x_k} P(\mathbf{O}_{k+1}|x_k)P(x_{k+1}|x_k)P(x_k) \\
	&= \sum_{x_k} P(\mathcal{O}_k|x_k)P(\mathbf{O}_k|x_k)P(x_{k+1}|x_k)P(x_k) \\
	&= \sum_{x_k} P(\mathcal{O}_k|x_k)P(x_{k+1}|x_k)\ell_k(x_k).
\end{split}
\end{align}
This recurrence relation may be written in vector form
\begin{align}
   \ketR{\ell_{k+1}} = V_k(\mathcal{O}_k) \ketR{\ell_k}.
\end{align}
Here, we have introduced the matrix $V_k(\mathcal{O}_k)$ with elements
\begin{align}
   &\braR{x}V_k(\mathcal{O}_k)\ketR{y} = P(x_{k+1}|y_k) P(\mathcal{O}_k|y_k) = w_k^{x y}\, \mathcal{P}_k^y(\mathcal{O}_k),
\end{align}
and we have defined
\begin{align}
\begin{split}
   &w_k^{x y} = P(x_{k+1}|y_k), \\
   &\mathcal{P}_k^x(\mathcal{O}_k) = P(\mathcal{O}_k|x_k). \label{eq:modelMatrices}
\end{split}
\end{align}
The matrix $w_k^{x y}$ describes the transition probabilities for the evolution of the logical qubit state and the vector $\mathcal{P}_k^x(\mathcal{O}_k)$ describes the observation noise for each logical qubit state. Note that there exist corresponding recurrence relations for the state vector $\ketR{\varrho_k}$ and for the likelihood function $\mathcal{L}_k$. They take the form
\begin{align}
	\ketR{\varrho_{k+1}} = \frac{1}{\mathcal{N}_{k+1}} \ketR{\tilde{\varrho}_{k+1}}, \;\;\; \mathcal{L}_{k+1} =  \mathcal{N}_{k+1} \times \mathcal{L}_k,
\end{align}
where
\begin{align}
	\ketR{\tilde{\varrho}_{k+1}} = V_k(\mathcal{O}_k) \ketR{\varrho_k}, \;\;\; \mathcal{N}_{k+1} = \Tr\ketR{\tilde{\varrho}_{k+1}}.
\end{align}

\subsubsection{Numerical algorithm \label{sec:numericalAlgorithm}}

We now provide an efficient numerical algorithm for forward filtering. The algorithm simultaneously calculates the probability vector $\ketR{\varrho_k}$ and the log-likelihood function $\ln \mathcal{L}_k$ as follows:
\begin{enumerate}
	\item Set $k=0$.
	\item Calculate the matrix $V_k(\mathcal{O}_k)$.
	\item Calculate $\ketR{\tilde{\varrho}_{k+1}} = V_k(\mathcal{O}_k) \ketR{\varrho_k}$.
	\item Calculate the norm $\mathcal{N}_{k+1} = \Tr \ketR{\tilde{\varrho}_{k+1}}$.
	\item Update the probability distribution $\ketR{\varrho_{k+1}} = \ketR{\tilde{\varrho}_{k+1}}/\mathcal{N}_{k+1}$.
	\item Update the log-likelihood ratio with $\ln\mathcal{L}_{k+1} = \ln\mathcal{L}_{k} + \ln \mathcal{N}_{k+1}$.
	\item Increase $k$ by one and start again.
\end{enumerate}
Note that we are only interested in estimating the initial state of the logical qubit. Therefore, the matrix $V_k(\mathcal{O}_k)$ may be normalized at each step by any constant factor independent of the qubit state without affecting the maximum-likelihood estimate. In some cases, this may prevent the values of the log-likelihood function from becoming too large.

\section{Repetitive readout statistics \label{sec:repetitiveStatistics}}

The hidden Markov model relevant for the repetitive readout discussed in the main text is obtained by specifying the transition matrix $w_k$ and the noise vector $\mathcal{P}_k$ appearing in Eq.~\eqref{eq:modelMatrices}.

As discussed in the main text, the qubit undergoes relaxation to its spin ground state on a time scale $T_1$. We obtain the value of $T_1$ by simultaneously fitting the measured single-repetition probabilities $P_1$ and $P_0$ shown in Fig.~2b to an expression of the form
\begin{align}
\begin{split}
P_1(t) &= A \epsilon^{-t/T_1} + B,
\end{split}
\end{align}
for the case when Q1 is prepared in $\ket{1}$ before the first measurement cycle, and of the form
\begin{align}
P_0(t) &= B.
\end{align}
for the case when Q1 is initially prepared in $\ket{0}$. Here, $t$ is the total readout time after the beginning of the repetitive readout and $A$ and $B$ are constants. The fit is shown in Fig.~2b and yields $T_1 = 1.8 \pm 0.6 s$. For a relaxation process, the transition matrix $w_k^{x y}$ in the basis $\left\{\ketR{1},\ketR{0}\right\}$ takes the form:
\begin{align}
	w_k =
	\textrm{exp}
	\left[
	\left(
	\begin{array}{cc}
	-1 & 0\\
	 1 & 0
	\end{array}
	\right)
	\frac{\delta t_{\textrm{rep}}}{T_1}
	\right]. \label{eq:lindbladian}
\end{align}
Here, $\delta t_{\textrm{rep}} = 3.263\,\textrm{ms}$.

The noise vector $\mathcal{P}_k^y(\mathcal{O}_k)$ describing the distribution of the outcome $\mathcal{O}_k$ for each logical qubit state is given by
\begin{align}
	\mathcal{P}_k(\mathcal{O}_k) =
	\left(
	\begin{array}{c}
	P(\mathcal{O}_k|x_k = 1) \\
	P(\mathcal{O}_k|x_k = 0)
	\end{array}
	\right). \label{eq:softNoiseVector}
\end{align}
Here, $P(\mathcal{O}_k|x_k = 1)$ and $P(\mathcal{O}_k|x_k = 0)$ are the probability distributions of the readout outcome $\mathcal{O}_k$ for preparation of the logical qubit in states $1$ and $0$, respectively. For the soft decoding procedure discussed in the main text, the readout outcome $\mathcal{O}_k$ is taken to be the peak signal $I_p$ (see Sec.~\ref{sec:singleRepetition}). The distributions of outcomes conditioned on the logical qubit state are then simply the empirically observed distributions of the peak signal $P(I_p|x_k = 1)$ and $P(I_p|x_k = 0)$ displayed in Figs.~3a-b of the main text. For the hard decoding procedure, the readout outcome $\mathcal{O}_k$ is taken to be the binary value $1$ or $0$ obtained by thresholding the peak signal at each repetition. The distributions of outcomes conditioned on the logical qubit state are then given by the conditional single-repetition readout error rates $\epsilon_1$ and $\epsilon_0$:
\begin{align}
\begin{array}{ll}
	P(1|x_k = 1) = 1 - \epsilon_1,  &P(1|x_k = 0) = \epsilon_0, \\
	P(0|x_k = 1) = \epsilon_1,  &P(0|x_k = 0) = 1-\epsilon_0 . \label{eq:hardNoiseVector}
\end{array}
\end{align}
The procedure used to obtain the distributions $P(I_p|x_k = 1)$ and $P(I_p|x_k = 0)$ and the conditional single-repetition readout error rates $\epsilon_1$ and $\epsilon_0$ is detailed in Sec.~\ref{sec:singleRepetition}.

\section{Single-repetition readout calibration \label{sec:singleRepetition}}

Since all repetitions are identical, we use the first repetition to calibrate the single-repetition readout. At $t_0 = 0$, the logical qubit is prepared $10^4$ times in state $x_0 = 1$ and $10^4$ times in state $x_0 = 0$. For each preparation, a readout trace $I(t)$ is recorded for a total time of $2.01\,\textrm{ms}$ in steps of $\delta t = 16.38\,\mu\textrm{s}$. For all readout times $t_R \leq 2.01\,\textrm{ms}$, we extract the peak signal $I_p = \max_t I(t)$~\cite{elzerman2004} and construct the probability distributions of $I_p$ conditioned on the logical qubit state, $P(I_p|x_0=1)$ and $P(I_p|x_0=0)$. For a given readout time $t_R$, single-repetition readout is performed by calculating the single-repetition log-likelihood ratio
\begin{align}
	\lambda(I_p) = \ln \frac{P(I_p|x_0=1)}{P(I_p|x_0=0)}.
\end{align}
If $\lambda(I_p) > 0$ [$\lambda(I_p) < 0$], we decide that the qubit state is most likely $\ket{1}$ [$\ket{0}$]. Therefore, the single-repetition readout error rates $\epsilon_1$ and $\epsilon_0$ for each state are obtained from marginals of the distributions $P(I_p|x_0 = 1)$ and $P(I_p|x_0 = 0)$ as follows:
\begin{align}
\begin{split}
	\epsilon_1 = P(\lambda^{\textrm{log}}<0|x_0 = 1) = \sum_{\left\{I_p | \lambda(I_p)<0\right\}} P(I_p|x_0 = 1), \\
	\epsilon_0 = P(\lambda^{\textrm{log}}>0|x_0 = 0)= \sum_{\left\{I_p|\lambda(I_p)>0\right\}} P(I_p|x_0 = 0), \label{eq:conditionalSingleShot}
\end{split}
\end{align}
and the average single-repetition readout error rate is given by:
\begin{align}
	\epsilon = \frac{1}{2}\left(\epsilon_1 + \epsilon_0\right)
\end{align}
The value of $\epsilon$ is plotted in Fig.~\ref{fig:figS1} as a function of readout time in the absence and in the presence of artificially added noise. We choose the value of $t_R$ that minimizes $\epsilon$. The histograms $P(I_p|x_0 = 1)$ and $P(I_p|x_0 = 0)$ and the error rates $\epsilon_1$ and $\epsilon_0$ corresponding to that optimal value of $t_R$ are those used throughout the main text. In particular, they are used to perform the Bayesian analysis detailed in Sec.~\ref{sec:repetitiveStatistics}. Without adding artifical noise to the readout traces, we find $\epsilon_1 = 32.9 \pm 0.5\%$ and $\epsilon_0 = 16.2 \pm 0.4\%$ at an optimal readout time $t_R = 623\,\mu\textrm{s}$. After adding noise, we find $\epsilon_1 = 41.1\pm 0.5 \%$ and $\epsilon_0 = 42.3\pm 0.5\%$ at an optimal readout time $t_R = 475\,\mu\textrm{s}$. The error bars on $\epsilon_1$ and $\epsilon_0$ are given by the standard deviation of the corresponding binomial error process with $N = 10^4$ samples, $\delta \epsilon_{1} = \sqrt{\epsilon_{1}(1-\epsilon_{1})/N}$ and $\delta \epsilon_{0} = \sqrt{\epsilon_{0}(1-\epsilon_{0})/N}$.
\begin{figure}
\center{\includegraphics[width=\columnwidth]{figS1}}
\caption{Average single-repetition readout error rate $\epsilon$ as a function of readout time $t_R$ in the absence of added noise (solid blue line) and in the presence of added noise (dashed yellow line). The optimal readout time in the absence (presence) of added noise is indicated by the dot-dashed (dotted) line. In both cases, the shaded areas give the statistical uncertainty (see text). Finally, we note that the small discontinuity in the blue curve at $t_R \approx 1.7\,\textrm{ms}$ is a consequence of the finite histogram bin size and has no particular physical significance.}
\label{fig:figS1}
\end{figure}

Note that we have assumed perfect preparation of the logical qubit state throughout. For this device, preparation errors have been reported~\cite{watson2018} to be $\sim 1\%$, smaller than the single-repetition readout errors reported in Sec.~\ref{sec:singleRepetition}. In Sec.~\ref{sec:theoreticalSimulation}, we argue that the average preparation error is in fact around $2-3\%$ ($4-5\%$ for state $\ket{1}$ and $< 1\%$ for state $\ket{0}$) for this run of the experiment. Therefore, the empirically measured single-repetition readout distributions in Figs.~3a-b are close to what they would be without preparation errors and lead to decoding procedures (hard and soft) that are close to optimal. In Sec.~\ref{sec:theoreticalSimulation}, we give rough estimates for the expected logical readout fidelity that would be measured in the absence of preparation errors.

\section{Theoretical simulation \label{sec:theoreticalSimulation}}

To explore the full potential of the repetitive readout discussed in the main text, we perform a numerical simulation. For each qubit state, we randomly sample $10^4$ outcome sequences according to the hidden Markov model described in Sec.~\ref{sec:repetitiveStatistics}. We use the experimentally extracted value of $T_1$ and the empirically measured peak-signal distributions of Fig.~3a as input parameters. A Monte-Carlo estimate of the average logical readout error rate $\epsilon^{\textrm{log}}$ is obtained by thresholding the simulated outcome sequences as described in Sec.~\ref{sec:singleRepetition} and then applying the hard decoding procedure detailed in Sec.~\ref{sec:bayesianInference}.
\begin{figure}
\center{\includegraphics[width=0.8\columnwidth]{figS2}}
\caption{Simulated logical readout error rate (solid blue line) as a function of the number of repetitions for (a) state $\ket{1}$ and (b) state $\ket{0}$. In both cases, the experimentally measured logical readout error rate is shown for comparison (dashed yellow line). The bands give the statistical error associated with the binomial statistics of the sampled errors as explained in Sec.~\ref{sec:singleRepetition}. Moreover, the dotted black curves give the simulated curves corrected to fit the experimental curves accounting for finite preparation errors using Eq.~\eqref{eq:errorComposition}. The best fit values for the preparation errors are $\eta_1 = 4\%$ for state $\ket{1}$ and $\eta_0 = 0.44\%$ for state $\ket{0}$, resulting in an average preparation error of $\eta = 2.2\%$. }
\label{fig:figS2}
\end{figure}

The conditional logical readout error rates $\epsilon_1^{\textrm{log}}$ and $\epsilon_0^{\textrm{log}}$ are plotted in Fig.~\ref{fig:figS2} as a function of the number of repetitions. The experimental results of Fig.~3c are shown for comparison. It is clear from Fig.~\ref{fig:figS2} that the simulated logical readout error rate underestimates the experimentally measured readout error rate as the number of repetitions increases. We believe that this is due to the logical readout error rate becoming comparable to preparation errors for Q1. This assumption is consistent with the simulation and experiment agreeing at low repetition number, where the logical readout error rate is much larger than the expected preparation error (of the order of $1\%$~\cite{watson2018}). In fact, it is possible to estimate the preparation error $\eta$ by fitting the simulated error rate $\epsilon^{\textrm{sim}}$ to the experimentally measured error rate $\epsilon^{\textrm{log}}$ using the error composition relation
\begin{align}
  \epsilon^{\textrm{log}} = (1-2\eta)\epsilon^{\textrm{sim}} + \eta. \label{eq:errorComposition}
\end{align}
We find that the simulation and experiment agree for a preparation error of $\eta_1 \approx 4-5\%$ for state $\ket{1}$ and of $\eta_0 < 1\%$ for state $\ket{0}$, giving an average preparation error of approximately $\eta \approx 2-3\%$. The simulation results suggest that in the absence of preparation errors, the measured fidelity could reach $\sim 98\%$ after 15 repetitions and saturate at $>99\%$ for more than 30 repetitions.

\section{Engineering Gaussian distributions}

It is known that if $P(I_p|1)$ and $P(I_p|0)$ are unimodal Gaussian distributions, soft decoding can reduce the number of readout cycles by up to a factor of two for arbitrarily large readout SNR. For this advantage to exist, the `bit-flip' conversion errors creating the bimodal features of the measured probability distributions must be small enough that the errors are dominated by the Gaussian noise~\cite{danjou2014-2}. Understanding the origin of the conversion errors is thus of great importance. Here, the `bit-flip' conversion errors may arise from an imperfect CROT gate, imperfect ancilla preparation, or imperfect spin-to-charge conversion in the readout of Q2. The CROT and preparation errors will be generically suppressed as the control and  preparation fidelities are improved. The imperfect spin-to-charge conversion in the readout of Q2 may be the most challenging to overcome. In what follows, we discuss avenues to suppress the spin-to-charge conversion errors.

For initialization of Q2 in state $\ket{1}$, the dot should ideally remain empty at all times. Thus, imperfect spin-to-charge conversion arises either from
\begin{enumerate}
\item[1)] the finite time scales for an electron to leave or re-enter the dot. The conversion errors caused by these charge transitions could be suppressed by, e.g., engineering the density of states~\cite{mottonen2010} of the reservoir so that a spin-up electron tunnels out of the quantum dot very rapidly, and so that a spin-down electron tunnels back in the quantum dot very slowly~\cite{danjou2014}. If that were the case, the dot would remain empty at nearly all times and the bimodality of $P(I_p|1)$ would be suppressed.

\item[2)] spin relaxation of Q2 before the electron is able to leave the dot. Indeed, the spin relaxation time of Q2 is of order $1\,\textrm{ms}$, which is only an order of magnitude longer than the observed time scale $\sim 100\,\mu\textrm{s}$ for an electron to escape to the reservoir. The strategy to suppress spin relaxation of Q2 during readout depends strongly its physical origin. However, there is \emph{a priori} no fundamental reason why the spin relaxation time of Q2 cannot be as large as that of Q1 (although it must be noted that increasing the dot-reservoir coupling as suggested in point 1) may negatively affect the spin relaxation time via cotunneling processes).
\end{enumerate}

For initialization of Q2 in state $\ket{0}$, the dot should ideally remain occupied at all times. The probability of an electron leaving the dot should in principle be exponentially suppressed in the ratio of the thermal energy to the Zeeman splitting. However, the observed probability of transitioning to the reservoir is too large to be explained by such thermal suppression. Instead, a measurement of the transition probability as a function of the plunger gate voltage of Q2 suggests that the electron escapes the $\ket{0}$ state via an excited quantum dot state. The value of the plunger voltage at which the transitions are suppressed is consistent with an excited state with energy tens of $\mu\textrm{eVs}$ above the ground state. We conjecture that this state is an excited valley state with the same spin (spin-down) as the ground state. Therefore, a likely explanation for the dominant transition mechanism is excitation to a higher valley state (via, e.g., absorption of energy from the biased charge sensor) followed by a transition from the excited valley state to the reservoir. The bimodality of $P(I_p|0)$ could therefore be suppressed by engineering valley splitting much larger than the charge sensor bias.

The above discussion highlights the importance of understanding the underlying physics for optimization of qubit readout.


\section{Introduction} \label{sec:intro}

The Transiting Exoplanet Survey Satellite (TESS) searches for extrasolar
planets by observing a sector of the sky for 27.4 days and using high-quality
photometric measurements of the stars in the field to look for transits
\citep{RickerEtal:TESS}.  A TESS sector covers 84\degr{} of ecliptic
\replaced{longitude}{latitude} (\replaced{$|l|>6$\degr{}}{$|\beta|>6$\degr{}}
) with four cameras, each with a 24\degr{} field of view (21-arcsec pixels).
Images are obtained every 2~minutes and then coadded into 30-minute exposures
that are saved as Full-Frame Images (FFI).  During its 27.4-day window, each
sector may produce as many as \replaced{1348}{1315} observations.  TESS will observe 26 of
these sectors (13 South and 13 North of the ecliptic plane) covering 90\% of
the sky during its 2-year primary mission.

Because of the large field-of-view, comets and asteroids serendipitously
appear in the TESS data.  Approximately 50~comets will be bright enough to be
detected in the two years of TESS observations, and we are using the
long-duration, high-cadence observations to monitor these comets for
rotational variability, outburst activity and coma studies.  The data are
ideal for exploring lightcurve behavior, because the 30-minute cadence
permits studies of short-term variability, while the 27-day sequences allow
measurements of periodicities $>$24~hrs (which are difficult to obtain from
ground-based observations) and provide the opportunity to search for changes
in the rotation rate \citep[e.g.,][]{BodewitsEtal:tgk}.  The TESS
observations are also well-suited for outburst studies, as the continuous
coverage and high-precision photometric capabilities allow searches for even
small events.  For any outbursts that occur during a 27-day sector window,
the 30-minute cadence will accurately document the start-time, brightening
profile, and peak brightness---characteristics that are rarely captured for
typical outburst detections---allowing for detailed comparisons that could
provide clues to the mechanisms at work in \replaced{these}{this} phenomenon.  Finally, the
wealth of observations allow for studies of the spatial structures in
cometary comae, including the ability to coadd large amounts of data (phased
to rotation periods if they are known) to investigate faint features that
might exist.

In its approach to perihelion, comet 46P/Wirtanen appeared within the TESS
field of view.  Wirtanen is a particularly interesting comet for
demonstrating TESS's potential contributions to cometary science. It is
a small hyperactive comet that made an historically close approach to the
Earth (0.077 AU) in December 2018, and is a potential target for future
spacecraft missions. We report here on first results from the observations of
a comet with TESS.

\section{TESS Observations} \label{sec:observations}

During the 2018 TESS Sector~3 observations, Wirtanen appeared in the field of
view of Camera~2, CCD~\#3. From September~20 to October~17, 1288 FFI images
of the comet were acquired.  During this time, the comet's heliocentric
distance dropped from 1.51 to 1.30~AU and the geocentric distance from 0.58
to 0.36~AU, while the solar phase angle changed from 23.2\degr{} to
29.8\degr{}.  During the first and last three days of the
sector observations, the spacecraft was testing new pointing software,
increasing the jitter \citep{FausnaughEtal:TESS_sec3} and causing $\sim$340
observations to be considered ``non-science-quality''.  Other suboptimal
images include a sequence in which the scattered light from Earth was not
well removed and occasional frames that are smeared due to TESS momentum
dumps that occur every $\sim$2.5 days.  Although some of these data may be
recoverable, we have removed them for the initial analyses reported here.  In
total, we used 935 images spanning 20.4~days.

\section{Data Reduction and Analyses} \label{sec:reduction}

The FFI data are crowded with stars that must be removed to minimize their
interference with measurements of the moving comet and we adopted the TESS
user-provided {\em Difference Image Analysis} software (DIA;
\citealt{OelkersStassun:DIA}) to perform
this task. DIA is designed to highlight objects whose brightnesses vary with
time, and as an intermediate step in the reduction procedures the scattered
light (from the Earth, Moon, etc.) and background stars are removed from each
image.  We use these intermediate, cleaned images, an example of which is
shown in Figure~\ref{fig:cleaned_image}, for our comet analyses.  The
procedure does a good job removing the fainter stars and the wings of
brighter stars, but due to the undersampled PSF
\citep{VanderEtal:TESS_instrum}, uncorrelated residuals remain at the centers
of bright stars.  These residuals can still interfere with the comet
photometry, but at a significantly lower level than the stars in the original
images.

TESS pipeline products are calibrated to $\mathrm{e^\textrm{-}~s^{-1}}$,
which can be converted to TESS magnitudes (15,000
$\mathrm{e^\textrm{-}~s^{-1}}$ corresponds to $T=10.0$;
\citealt{VanderEtal:TESS_instrum}).  Using the Web TESS Viewing
Tool\protect{\footnote{\url{https://heasarc.gsfc.nasa.gov/cgi-bin/tess/webtess/wtv.py}}},
which requires an object's brightness in multiple bandpasses (we used
typical comet colors $V-J=1.47\pm0.17$ and $J-H=0.52\pm0.12$,
\citealt{HartmannEtal:cometcolor}; and $V-R=0.50\pm0.03$,
\citealt{Jewitt:colors}) we produced a rough conversion to $R$ magnitudes,.
The systematic uncertainty in $R$ is dominated by the uncertainties in the
comet colors, with a conservative estimate of $\pm0.3$~mag.  Relative
photometry from image to image, on the other hand, is more tightly
constrained, with errors $<0.017$~mag.  Flux units are computed using the
conversion coefficient for a $R=0$ star of
$2.18\times10^8~\mathrm{W~m^{-2}~\mu m^{-1}}$ \citep{BesselEtal:zeropoints}.

\begin{figure}[hbt]
\plotone{Fig_image_clean.pdf}
\caption{Example of the star and background removal, showing the region of a
  sample frame (tess2018284152940-s0003-2-3-0123-s\_ffic\_sa) surrounding
  comet Wirtanen before (left) and after (right) the DIA cleaning
  process. The images \added{are 1.22$\times10^6$~km across and} are
    displayed with the same logarithmic scale.\label{fig:cleaned_image}}
\end{figure}

\subsection{Lightcurve and Outburst} \label{sec:outburst}

Centering on the optocenter of the coma, we measured the brightness of the
comet in multiple apertures with fixed radii ranging from 15,000 to 40,000~km
at the distance of the comet.  To produce our lightcurve, we ultimately
selected the 25,000~km aperture, which is large enough to contain at least
3~PSFs at all comet distances, but minimizes the number of residual stars
crossing the aperture.  The resulting lightcurve, shown in
Figure~\ref{fig:lightcurve}, shows the comet continuously brightening over
the observation window, but more importantly, it reveals that TESS captured a
moderate-sized (0.5~mag) outburst---the only one to be reported during the
2018 apparition.  This event is likely to be the best \added{temporally}
characterized natural
outburst ever recorded, with the nominal pre-outburst behavior, the onset and
rise to the peak, and ultimately the falloff over the course of several
weeks, all being observed at a 30-minute cadence.  The observations show a
two-phase brightening, with a rapid $\sim$1~hr jump commencing on
September~26.12\added{$\pm$0.01}, followed by a more gradual increase that continued for
another $\sim$8~hr.  After peaking, the coma began a roughly exponential fading
that lasted 15-20~days.  Because the comet's pre-outburst brightening rate is
not well constrained, it is not clear whether or not the comet had returned
to its baseline level by the end of the observation window.

\begin{figure}[hbt]
\plotone{Fig_outburst_photom.pdf}
\caption{Comet Wirtanen's lightcurve as measured in a 25,000~km
  \added{radius} aperture from the TESS images. The top panel shows the
  comet's secular brightening (black points) and reveals an outburst that
  begins at September~26.12. Removing an assumed linear baseline (dotted
  line) highlights the outburst behavior (red points).  The lower panel
  expands the section around the onset of the outburst, showing the rapid,
  hour-long jump (0.21~$\mathrm{mag~hr^{-1}}$), followed by a more gradual
  brightening (0.034~$\mathrm{mag~hr^{-1}}$) that peaks $\sim$8~hrs later.
  See Section~\ref{sec:reduction} for a discussion of the
  uncertainties. \label{fig:lightcurve}}
\end{figure}

Following up on this discovery, we inspected the images for additional
information about the outburst.  We registered and coadded the data in 3~hour
blocks to improve the S/N, and then enhanced the post-outburst images by
dividing out a pre-outburst frame.  The resulting sequence (e.g., Figure
\ref{fig:outburst_image}) shows the rapid brightening of the central coma,
followed by an outflow of material that is roughly symmetric around the
nucleus. The leading edge of this outflow had a projected velocity
$\sim$800~$\mathrm{m~s^{-1}}$, suggesting that it was gas related to, and
possibly driving, the outburst.  Given the detector bandpass (600--1000~nm)
the gas was most likely dominated by CN, with emission bands at 914.1 and
787.3~nm \citep{schleicher:thesis}.  The gas contributes only a few percent
of the total brightness \replaced{increase, and shows no signature in
the lightcurve as it exits the photometric aperture}{within the aperture}.

\begin{figure}[hbt]
\plotone{Fig_outburst_image.pdf}
\caption{Image sequence showing the outburst's effect on Wirtanen's coma.
  The first panel shows the pre-outburst morphology, enhanced by dividing out
  a 1/$\rho$ profile.  To show changes in the coma with time, panels 2--6 are
  enhanced by \replaced{dividing}{subtracting} out the (unenhanced) image
  from panel~1, \added{scaled by the linear trend shown in
    Figure~\ref{fig:lightcurve}}.  Panels~2 and~3 bracket the onset of the
  outburst (September~26.12) and~\replaced{3--5}{4--6} show the bright
  central condensation and the rapidly expanding gas cloud.  \deleted{The gap
    on the West side of the gas cloud is an artifact resulting from dividing
    out the bright dust tail.} Each panel is 400,000~km across, with North
  up and East to the left. \added{The light blue circle denotes a
    25,000~km radius aperture.}  \label{fig:outburst_image}}
\end{figure}

In images over the few days following the outburst, the rapidly moving gas
diffuses away, while, in contrast, the bright central peak (presumably dust)
remained highly concentrated, showing little radial expansion over the
following 15-20 days. This material left our photometric aperture primarily
when it dispersed down tail under the influence of radiation pressure
\added{(see Figure~\ref{fig:dusttail_image})}.  Thus, the optically dominant
grains ejected in this event must have been moving slowly, with expansion
velocities of a few tens of~m/s, (\added{consistent with those seen in
  outbursts in comet 49P/Arend-Rigaux and 67P/Churyumov-Gerasimenko}
\citep{EisnerEtal:arend,LinEtal:67Poutburst}) or else they would have left
the photometric aperture more rapidly than was observed.

\begin{figure}[hbt]
\plotone{Fig_dusttail.pdf}
\caption{\added{Image sequence showing the effects of Wirtanen's outburst on
    its dust tail.  The first panel highlights the outburst, followed by
    images that show a narrow tail forming in the anti-solar direction over
    subsequent days. Each panel consists of 8~hrs worth of images, centered
    on the stated date and rotated to align the sunward direction with that
    of the Sep~27 panel before coadding.  Images were enhanced by combining
    all images between September~28 and October~9 (again rotated to align the
    sunward direction) into a single temporally averaged template, which was
    then divided out of each of the eight coadded panels.  Images are 100
    pixels across, which ranges from 780,000~km on September~27 to 640,000~km
    on October~8.}  \label{fig:dusttail_image}}
\end{figure}

\subsection{Dust Trail} \label{sec:trail}

Examination of the TESS images also revealed that Wirtanen has a previously
unreported dust trail.  Due to the presence of the comet itself, the DIA
software oversubtracts the background in the vicinity of the coma and trail.
To better investigate the trail, we grouped images into 5-day bins,
re-projected the raw frames into the comet's rest frame with the velocity
vector aligned, masked the bright stars, and then combined the images via the
median value at each pixel.  The resultant images were additionally filtered
to remove the CCD banding artifacts.  A sample mosaic is presented in
Figure~\ref{fig:trail}.

The dust trail extends in both directions beyond the~24\degr{} TESS
field-of-view, which covers a mean anomaly range ($\Delta M$)
of~$-1.6$\degr{} to~0.7\degr{}.  A Gaussian function fit to the width of the
trail at $\Delta M=-0.16$\degr{} (-1.7\degr{} on the sky) has a peak of
29.3~$\mathrm{mag~arcsec^{-2}}$ and a 1$\sigma$ width of
$(5.5\pm1)\times10^4$~km.  The width and optical depth
($\tau=7.0\times10^{-12}$) are consistent with the
faint end of the range of trails observed at infrared wavelengths
\citep{IshiguroEtal:DustTrails}, assuming a 4\% geometric albedo.

\begin{figure}[hbt]
\plotone{46p-tess-trail.pdf}
\caption{Sample mosaic showing Wirtanen's dust trail, produced by combining
  240~images from the 5-day period from September~28 to October~3. The top
  panel shows the full mosaic, spanning CCDs~3 and~4 in Camera~2; sections at
  the nucleus and outer ends are expanded below.  In the lower frame of each
  image pair, a cross is located at the nucleus' position and the comet's
  orbit is overlain, showing that the trail is centered on the comet's path.
  The direction arrows in the inset are shown for the nucleus' position at
  the midtime of the five-day window.  \label{fig:trail}}
\end{figure}

\section{Discussion} \label{sec:discussion}

Outbursts have been observed in comets for over a century
\citep[e.g.,][]{Hughes:outbursts,FiloChury:outburst}.  A number of proposals
have been put forth as trigger mechanisms
\citep{Hughes:outburstmech,GronkWeso:outburstrev}, though it is generally
believed that a single explanation is insufficient to account for all cases.
Similarly, little is known about the physical processes that govern the
behavior of the ejecta during the initial stages of the event.  The main
problem is that the majority of outbursts are detected only after they have
peaked in intensity and the quasi-exponential fading stage provides few
observational constraints on the conditions present during the explosive
phase.  The excellent coverage of the Wirtanen event allows us to use its
characteristics as a baseline model for the behavior of a moderate sized
outburst that can be compared and contrasted with other outbursts.  With this
information, we can begin exploring the physical processes at work.

\subsection{Characterization of Wirtanen's Outburst} \label{sec:character}

To begin this investigation, we can evaluate some potential trigger
mechanisms by looking at a comet's history.  Wirtanen is not known for
exhibiting outburst activity, and only three previous events were reported
since its 1948 discovery: 1991~October~7 (+16~days from perihelion, $<$1~mag,
\citealt{Kronk:cometog6}); 2002~September~25 (+29~days, $\sim$2~mag,
\citealt{Yoshida:website2002}, \citealt{Kidger:CGdust2}); and 2008~May~16
(+86~days, 2~mag, \citealt{Kidger:astroncirc2446}).  Including the 2018
outburst on September~26 (--81~days), there is no correlation in the orbital
position for these outbursts.  Furthermore, given the similar viewing
geometry on alternate apparitions, there are often observations from $\sim$11
years before or after the observed outbursts that rule out any events
occurring at those points on other orbits.  Thus, the lack of correlation
between outbursts and the comet's true anomaly or heliocentric distance,
combined with the general rarity of these events, suggests that a
persistently volatile region, as described by \cite{Miles:sw1_cryovol},
is not the cause of the Wirtanen events.

Next, we can characterize the scale of the event to provide context for
comparisons to other events.  The outburst produced a 0.5~mag brightness
increase.  For a more quantitative measurement, we convert this brightening
to an increase in reflective surface area, and then to a mass of the ejected
material.  We use the photometric analysis outlined by \cite{Jewitt:photom},
and assume that the brightness increase is dominated by sunlight reflected
off dust grains ejected into the coma. For an initial brightness $R=13.1$ and
a total brightening of 0.5~mag, we find an increase in the optical cross
section of ${\sim}2{\times}10^8$~m$^2$ (using the geometric parameters at the
start of the outburst, the Halley-Marcus composite dust phase function
compiled by
D.~Schleicher\protect\footnote{\url{https://asteroid.lowell.edu/comet/dustphase.html}},
and an assumed 4\% albedo).

Adopting a grain population with density 2000~$\mathrm{kg~m^{-3}}$, radii
ranging from 0.1~$\mu$m to 1~mm and a \added{differential size distribution
with a} power-law index of $-3.5$, we find that the above cross section
converts to a total dust mass of ${\sim}3{\times}10^6$~kg.  Unfortunately,
the mass tends to be dominated by the largest particles and the dust
properties in the outburst are not well constrained, so changing the assumed
upper limit of the grain size and/or the power law index can alter the
computed mass by as much as two orders of magnitude for plausible
populations.  Thus, our best estimate is that the outburst likely produced
$10^5-10^8$~kg of material.  This is equivalent to a crater a few 10s of
meters in radius (for a bulk nucleus density of 500~$\mathrm{kg~m^{-3}}$),
which is consistent with the outburst-associated features seen on comet
9P/Tempel~1 \citep{BeltonEtal:cryovol}.  Pits known to be associated with
outbursts are also seen on comet 67P/Churyumov-Gerasimenko
\citep{VincentEtal:cometactivity}.  Although these features are smaller
($<$10~m) than that computed for Wirtanen, the related C-G outbursts were
also smaller, and were not detected from Earth.

An important factor in the Wirtanen outburst is that it is one of the first
outbursts whose rise in intensity is well documented, showing a two-phase
increase in brightening.  Given that the gas was measured to expand outward
at $\sim$800~$\mathrm{m~s^{-1}}$, while the dust exhibited significantly lower speeds, we
know that during the first eight hours following the initiation of the
outburst, even the leading edge of the gas cloud would not have escaped our
25,000~km aperture.  This means that the two-phase brightness profile must
reflect temporal characteristics of the outburst (e.g., the changes in slope
are not due to material leaving the aperture).  We suggest that the initial,
hour-long period of brightening represents the energetic phase in which gas
and dust were rapidly ejected into the coma.  The second, more gradual phase
of brightening is probably due to the continued expansion of the dust from an
initially dense state to an optically thin regime.  However, it is also
possible that this phase arises from other causes: a continued excavation, at
a much slower rate as the outburst subsides; temporary enhanced emission of
gas and dust produced by sublimation of newly exposed ices; or the increase
in reflective surface area as a small number of grains in the coma gradually
fragment into smaller particles.

\subsection{Comparison to Other Outbursts} \label{sec:compare}

Our literature search revealed only three other occasions in which
high-cadence coverage was obtained of the brightening phase of an outburst: a
2017 event in 29P/Schwassmann-Wachmann~1 (SW1), the 2007 outburst of comet
17P/Holmes, and the DI experiment in comet Tempel~1 in 2005.  As discussed
below, there are significant differences in the basic nature of each of these
events, but we can compare and contrast the observed behavior in the early
stages of each case, to explore whether similar mechanisms may be at work.

SW1 is known for frequent outbursts.  About half seem to be periodic
\citep{MilesEtal:sw1_anatomy}, suggesting that the nucleus has volatile-rich
hot spots that are triggered by diurnal heating cycles.  On 2017~July~2, SW1
experienced a $\sim$2~mag outburst and \cite{MilesEtal:sw1outburst} reported
on high-cadence photometry of this event.  Few results from these
measurements have been reported, but \replaced{for comparison, we used the
  same computation described above to estimate an}{if we adopt the same
  assumptions as used in the Wirtanen estimate above, we compute an} increase
in the SW1 optical cross section of ${\sim}10^{11}$~m$^2$.  \added{Because
  conditions are dramatically different in the two comets, this is likely a
  poor comparison, but it suggests that the SW1 outburst was} significantly
larger than \added{that} seen in Wirtanen.  Even so, the behavior \added{of the
  lightcurve} during the initial stages appears to be nearly identical in
character to the Wirtanen event, with a rapid rise for $\sim$0.5~hr, followed
by a more gradual increase that continued for at least another hour, when the
observations ended.  The fact that both events exhibit similar behavior, even
in their timescales, suggests that the \added{physcial} mechanisms
\replaced{at work}{governing the ejecta} are the same, even though the
magnitude of the events were dramatically different.

Comet Holmes experienced a 15~mag brightening---the largest outburst ever
recorded---around 2007~October~23.7 \citep{MontaltoEtal:holmes,
  SteveJewitt:holmes}.  The onset of the outburst was not captured, but the
42-hour duration of the rise allowed observers to capture the later stages of
the brightness increase \citep{TrigoEtal:holmes, HsiehEtal:holmes,
  LiEtal:holmes}.  This lightcurve shows that there was a change in the rate
of brightening around October~24.1, around half a day into the outburst,
where the rate became {\em steeper}, before inflecting and then peaking a day
later.  This steepening in slope was shown to be the result of a rapid
cascade of large dust grains fragmenting into successively smaller particles,
which dramatically increased the reflective surface area in the coma
\citep{HsiehEtal:holmes, SteveJewitt:holmes}.  In the days following the
outburst, the leading edge of the ejected dust cloud was measured to have a
velocity $\sim$550~$\mathrm{m~s^{-1}}$ \citep[e.g.,][]{TrigoEtal:holmes}, an
order of magnitude higher than the dust velocities in Wirtanen, which likely
reflects the enormous amount of energy released in the event.

Our final comparative example is the Deep Impact experiment at comet Tempel~1
(T1).  Although this was a man-made outburst, it was known to be the result
of an impact and should exhibit the same phenomena as a naturally-triggered
impact event.  Furthermore, because it was planned, it was intensely observed
from both DI and from the ground, allowing remote observations to be
connected to specific events seen by the spacecraft.  High cadence photometry
of the event shows a three-phase brightening \citep{MeechEtal:DI_campaign,
  FernandezEtal:Tempel1, MitchellEtal:tempel1, KuppersEtal:Rosetta_DI}.  The
T1 lightcurve starts with a very sharp increase for the first $\sim$1~min,
followed by 6~min of gradual brightening, and then another 10-15 minutes with
a somewhat steeper slope before it flattened out at $\sim$2~mag above its
starting brightness.  The first two phases mimic the behavior of Wirtanen and
SW1, while the later steepening is comparable to that seen in Holmes.  The
outburst started to fade $\sim$45 mins after impact.  The peak velocity of
the ejecta was ~$\sim$200~$\mathrm{m~s^{-1}}$.  It is interesting to note that details in the
observed phenomena, including changes in the brightness slope, peak
intensity, and timing of the observed features, vary somewhat depending on
the aperture size, wavelength, etc., and exploring these differences should
provide additional information for exploring the behavior of the ejecta.

As was suggested in Section~\ref{sec:character} the two-phase brightness
increase is related to the initial ejection of highly dense material that
then becomes optically thin as it expands outward.  The T1 data provide a
test of this conjecture.  \cite{KoloEtal:tempel1} used DI observations of the
ejecta in front of the comet's limb to measure the optical thickness of the
dust. They found that the cloud started out optically thick, but tended to
thin within a few seconds in most areas.  However, a few bands of material
remained optically thick for as long as a minute.  This time
\replaced{scale}{frame} corresponds to that of the steep segment of the
photometric profile, suggesting that the change in slope may indeed be due to
the expansion of the ejecta into an optically thin regime.  \added{The
difference in timescales between Wirtanen and T1 are likely to be due to
differences in velocity as well as the properties of the ejecta.}

The comparisons between the initial stages of these four outbursts reveal
both similarities and differences.  The similarities suggest that common
processes govern the physics, regardless of the cause or size of the event,
while the differences, such as the shorter timescale seen in T1 or the long
duration of the Holmes brightness increase, could represent characteristic
signatures that provide clues to the mechanisms involved.  Results for comet
Wirtanen also show the promise that TESS observations hold for cometary
science.

\acknowledgments
This paper includes data collected by the TESS mission, which are publicly
available from the Mikulski Archive for Space Telescopes (MAST). Funding for
the TESS mission is provided by NASA's Science Mission directorate.

\vspace{5mm}
\facility{TESS}

\software{DIA \citep{OelkersEtal:DIA_CSTAR,OelkersStassun:DIA}, Astropy
  \citep{AstropyEtal:astropy, AstropyEtal:opensci}, SEP \citep{Barbary:sep}
}