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

Quantum key distribution (QKD) relying on single photons is regarded as one of the most mature quantum technologies \cite{Pirandola2020,Ma2020}. However, the impossibility of amplifying single photons sets restrictions to the transmission distance. Entanglement-based QKD schemes are able to overcome these range-limitations when embedded in quantum networks\cite{Gisin2007,Kimble2008}, while also exhibiting a lower vulnerability to eavesdropping attacks\cite{Ekert1991,Bennett1992,Acin2007,Gisin2002,Pirandola2020}. For both fiber-based \cite{Wengerowsky2020} and satellite-based \cite{Yin2020} quantum cryptography, the most prominent sources of entangled photon pairs to date are based on the spontaneous parametric down conversion (SPDC) process. These sources are commercially available and can be operated in a large temperature range\cite{Pan2021SPDC}. As a drawback, SPDC sources exhibit approximately Poissonian photon pair emission statistics \cite{Schneeloch_2019}, which severely limits their brightness when a high degree of entanglement - and thus a low qubit error rate (QBER) - is demanded. The biexciton-exciton (XX-X) spontaneous decay cascade in epitaxially grown semiconductor quantum dots (QDs) has been demonstrated to be a viable alternative to SPDC sources due to the sub-Poissonian entangled photon pair emission statistics\cite{Benson2000}. In particular, GaAs QDs obtained by the Al droplet etching technique\cite{Gurioli2019} are capable of emitting polarization-entangled photon pairs with a fidelity to the $\ket{\phi^+}$ Bell state beyond 0.98 \cite{Huber2018,Schimpf2021}, owing to an intrinsically low exciton fine structure splitting (FSS)\cite{Huo2013} and lifetime and a near-zero multi-photon emission probability even at maximum brightness\cite{Schweickert2018}. This allowed the demonstration of entanglement-based QKD with a QBER as low as \SI{1.9}{\percent} \cite{BassoBasset2021,Schimpf2021}.

Independent of the used QD materials, the best performance in terms of entanglement fidelity and biexciton state-preparation efficiency have been obtained by operating the QD sources at temperatures below \SI{10}{K} and by using the resonant two-photon-excitation (TPE) condition \cite{Hafenbrak2007,Muller2014}. One of the drawbacks of the very low working temperatures is that they are difficult to achieve in satellites, where strong payload restrictions have to be met. Starting from about \SI{30}{K}, cryo-coolers for lower temperatures become exceedingly bulkier than higher temperature models\cite{CryoEng}. Regarding TPE, and in general resonant excitation schemes, their main limitation is represented by random charge capture in the QD, normally resulting in a significant drop of the time fraction $\beta$ in which the QD is optically active (also known as blinking) \cite{Jahn2015,Jons2017}. For the generation of single photons via strictly-resonant excitation, blinking has been successfully suppressed by embedding QDs into charge-tunable devices \cite{Warburton2000}, which allow the charge state of a QD to be deterministically controlled \cite{Zhai2020,Somaschi2016}.
The questions whether blinking can be suppressed under TPE is not trivial, since - different from single photon resonant excitation - TPE requires excitation powers about 50 times larger, which can produce free carriers in the continuum via two-photon absorption in the barrier material.
Although two-photon absorption by a QD in a diode structure has been observed via photo-current measurements \cite{Stufler2006}, the usefulness of charge-tunable devices in the context of entangled-photon-pair generation with the TPE method has not been tested.

Here we investigate the optical properties of GaAs QDs embedded in a p-i-n tunnel diode at a temperature of at least \SI{20}{K}, demonstrate blinking-free emission of entangled photon pairs under TPE, and use these photons to successfully implement the BBM92 QKD protocol\cite{Bennett1992} under these conditions. The key generation happens between two buildings, connected by a \SI{350}{m} long single-mode fiber inside the campus of the Johannes Kepler University.

Varying the voltage applied to the diode within a certain window allows for a fine-tuning of the emission wavelength in a range of about \SI{0.2}{nm}, while keeping the blinking-free emission intact. These structures can therefore be vital for applications in quantum networks\cite{Gisin2007}, where a precise matching of the wavelengths of multiple emitters becomes important\cite{Giesz2015,Reindl2017} and the transmission rate scales with $\beta^2$\cite{Jons2017}.

\section{GaAs Quantum Dots in a Tunnel Diode Structure at a Temperature of \SI{20}{K}}

\begin{figure*}
\begin{center}
\begin{tabular}{c}
\includegraphics[]{pics/F1_Sample.pdf}
\end{tabular}
\end{center}
\caption
{ \label{fig:PL}
Photoluminescence properties of GaAs quantum dots (QDs) in a p-i-n diode structure at a temperature of \SI{20}{K}, excited by resonant two-photon excitation (TPE). (a) p-i-n diode structure with a tunnel barrier between the n-doped and the intrinsic region. The inset shows the principle of TPE, with $E_\text{P}$ the laser energy, $E_\text{B}$ the biexciton (XX) binding energy and $S$ the exciton (X) fine structure splitting. (b) Emission spectra at TPE conditions when sweeping the diode voltage $V_\text{P}$ in forward bias. The inset shows the diode current $I$ over $V_\text{P}$. The white dashed line indicates $V_{\text{P,0}}=\SI{0.3}{V}$, at which the diode is operated during the QKD experiment. (c) Emission at $V_\text{P}=V_{\text{P,0}}$. (d) Second order correlation function $g^{(2)}$ of the X signal with a time-bin of \SI{1}{\micro s} at $V_\text{P}=V_{\text{P,0}}$. The  $g^{(2)}$ is shown for the QD in the diode structure (red), indicating an on-time fraction of $\beta=1.00(2)$ and a QD without diode (black, dashed) with a typical value of $\beta\approx 0.3$. (e) Wavelength shift and $\beta$ for different deviations $\delta V_\text{P}=V_{\text{P}}-V_{\text{P,0}}$. The blue dashed lines indicates a value of $\beta$ corresponding to no blinking.}
\end{figure*}

The GaAs QDs in a p-i-n tunnel diode structure, as depicted in Fig. \ref{fig:PL}(a), were grown by Molecular Beam Epitaxy (MBE). The first functional layers are six pairs of Distributed Bragg Reflectors (DBRs), grown by alternating Al\ts{0.95}Ga\ts{0.05}As and Al\ts{0.20}Ga\ts{0.80}As layers with \SI{65.2}{nm} and \SI{56.6}{nm} thickness, respectively. The n-doped region is formed by a \SI{95}{nm} thick Al\ts{0.15}Ga\ts{0.85}As layer with an Si concentration of \SI{E18}{cm^{-3}}. A combination of layers (\SI{15}{nm} Al\ts{0.15}Ga\ts{0.85}As plus \SI{8}{nm} Al\ts{0.33}Ga\ts{0.67}As) acts as tunnel barrier between the n-layer and the QDs, which are obtained via the local Al droplet etching technique\cite{Gurioli2019} and covered by the \SI{268}{nm} thick Al\ts{0.33}Ga\ts{0.67}As intrinsic region. A \SI{65}{nm} layer of Al\ts{0.15}Ga\ts{0.85}As, doped with \SI{5E18}{cm^{-3}} carbon atoms forms the p-region. Layers of \SI{5}{nm} Al\ts{0.15}Ga\ts{0.85}As and \SI{10}{nm} GaAs, each doped with \SI{E19}{cm^{-3}} carbon, cap the structure in order to form a conductive surface and to protect the active region from oxidation. The DC voltage $V_\text{P}$ is applied to the p-contact with respect to the grounded n-contact (see supplementary for details about the electrical contacts). A solid immersion lens on top of the structure enhances the overall extraction efficiency to about \SI{3}{\percent}.

The sample containing the GaAs QDs is fixed with silver-glue on the cold-finger of a He flow cryostat and cooled to a temperature of \SI{20}{K}. The temperature is measured with a calibrated silicon diode placed under the cold-finger, so we estimate that the actual sample temperature could be up to \SI{5}{K} higher. One individual QD is optically excited by a focused  pulsed laser with a repetition rate of \SI{80}{MHz} and the laser energy $E_\text{P}$ tuned to half of the biexciton (XX) energy, as depicted in the inset of Fig. \ref{fig:PL}(a). This resonant two-photon excitation (TPE) scheme is the same as used in previous experiments\cite{Schweickert2018,Huber2018,Schimpf2021}.
Figure \ref{fig:PL}(b) depicts the micro-photoluminescence spectra of a QD when adjusting $E_\text{P}$ at the working point $V_\text{P,0}=\SI{0.3}{V}$, used for the further measurements, and then sweeping $V_\text{P}$ from 0 to \SI{1.5}{V}. The inset shows the corresponding diode current. In the voltage range of about \SI{0.15}{V} to \SI{0.35} the QD is in its charge neutral configuration and only the biexciton (XX) and the exciton (X) transition lines are visible. For higher voltages, a single electron tunnels into the QD so that the negative trion (X\textsuperscript{-}) is addressed via phonon-assisted excitation. Figure \ref{fig:PL}(b) shows the spectrum at $V_\text{P,0}$, with the XX excited at pi-pulse. At this voltage, the auto-correlation function $g^{(2)}$ for a time span of \SI{100}{\micro s} and a time-bin of \SI{1}{\micro s} was recorded, depicted as the red data points in Fig. \ref{fig:PL}(d) (In this case, the time-bin is much larger than the excitation period of \SI{12.5}{ns}, so that the anti-bunching from the single-photon emission becomes invisible). The $g^{(2)}$ is one for all times delays, which indicates the complete absence of blinking, i.e., an on-time fraction of $\beta=1/g^{(2)}(0)=1.00(2)$\cite{Jahn2015}. The black dashed line indicates the $g^{(2)}$ typically measured for previously used GaAs QDs without diode, resulting in $\beta\approx0.3$.

It is interesting to note that the tunnel provides a window for $V_\text{P}$ of about \SI{0.1}{V} in which the QDs can be operated without blinking. This offers a tuning range for the central emission wavelengths of a about \SI{0.2}{nm} as depicted for a different, but representative QD in Fig. \ref{fig:PL}(e). The blue dashed line indicates the case of $\beta=1$, corresponding to no blinking. The same tuning-range was observed also for temperatures well below \SI{20}{K}.

\section{Characterization of the Entangled Photon Pairs}

\begin{figure*}
\begin{center}
\begin{tabular}{c}
\includegraphics[]{pics/F2_Entanglement.pdf}
\end{tabular}
\end{center}
\caption
{ \label{fig:Entanglement}
Emission properties relevant for the polarization entanglement, measured at a temperature of $\SI{20}{K}$. (a) Spectra of the individually filtered XX and X emission lines, combined at a 50:50 fiber beam splitter. (b) Single-photon emission characteristics of the XX and X signals, observed by detecting coincidences in a Hanbury-Brown-Twiss arrangement. (c) Lifetime trace of the XX and X signals. The X signal exhibits a slow secondary decay channel, which is not present at temperatures lower than \SI{10}{K}. (d) Example among the 36 recorded coincidence histograms between the XX and X detection, corresponding to a measurement in the HV basis. The red dashed lines indicate the time-bin of \SI{2}{ns} in which the coincidences are summed up to calculate the peak areas. (e) Unpolarized coincidence measurement between the XX and X photons. The excess coincidences at zero time delay stem from a non-unity photon-pair generation probability. (f) Density matrix of the two-photon polarization entangled state of the XX and X photons, recorded by full state tomography.}
\end{figure*}

The QDs of the here used diode sample exhibit an average FSS of about \SI{6}{\micro eV}. Average values below \SI{4}{\micro eV} are regularly achieved for nano-holes created at a substrate temperature of \SI{600}{\degree C} \cite{Huo2013,Keil2017}, while a slightly higher temperature (about \SI{610}{\degree C}) was used here to possibly improve the crystal quality. Instead of resorting to strain tuning to bring the FSS of one individual QD to zero\cite{Huber2018}, the stochastic distribution of the FSS among all QDs on the sample can be used to find a QD with a suitably low FSS for the given use case. In order to estimate the FSS required for a serviceable QBER in a BBM92 QKD arrangement we first calculate the 2-qubit density matrix in polarization space\cite{Hudson2007}, considering only the FSS and the typically observed X lifetime of about \SI{230}{ps}. From the density matrix $\rho$, we calculate the expected QBER via
\begin{equation}
    q=\frac{1}{2}\sum \limits_{i=1}^{4} \braket{O_i|\rho|O_i},
\label{eq:qber}
\end{equation}
where $O_i \in \{H_\text{A} V_\text{B}, V_\text{A} H_\text{B}, D_\text{A} A_\text{B}, A_\text{A} D_\text{B}\}$ are the 2-qubit measurement bases between Alice (A) and Bob (B), in which a measured coincidence corresponds to a false key bit. After a brief scanning we chose a QD with an FSS of \SI{0.96(9)}{\micro eV}, for which we calculate a minimum QBER of \SI{2.7}{\percent}.

The XX and the X photons generated by TPE are filtered out individually and coupled into single mode fibers. Figure \ref{fig:Entanglement}(a) shows the spectra of the XX and X emission lines merged at a 50:50 fiber beam splitter, of which one output is sent to the spectrometer (The slightly lower XX signal intensity stems from the higher distance from the objective to the fiber collimator compared to the X signal, resulting in a lower coupling efficiency).

Before performing the QKD experiment, a characterization of the single-photon emission characteristics and the polarization entanglement between the XX and the X photons is performed, as those properties primarily determine the QBER during the key generation process. The most important quantities are summarized in Table \ref{tab:performance} and compared to the values for a different QD (in the same diode structure), measured at $\SI{5}{K}$. Figure \ref{fig:Entanglement}(b) depicts a coincidence histogram of an auto-correlation measurement for both the XX and the X signals. Using a time-bin of $\SI{2}{ns}$, the results are $g^{(2)}_\text{XX}(0)=0.034(4)$ and  $g^{(2)}_\text{X}(0)=0.020(3)$, which are comparable to previously measured values using the same optical arrangement at a temperature of \SI{5}{K} (see Table \ref{tab:performance}), but with a slightly higher value for the XX signal. Compared to data acquired at \SI{5}{K}, one can observe a broadening of the peaks belonging to the X signal. Inspecting the lifetime traces, which are shown for both the XX and the X in Fig. \ref{fig:PL}(c), it becomes evident that this broadening stems from a slow decay channel of the X, which overlays with the mono-exponential decay from the bright X to ground state usually observed at \SI{5}{K} (see supplementary for lifetime- and cross-correlation measurements at \SI{5}{K}). A convoluted fit results in an X lifetime of $T_{1,\text{X}}=\SI{252(9)}{ps}$, with the caveat that this value is slightly overestimated due to the presence of the slow decay channel. The XX lifetime is measured as $T_{1,\text{XX}}=\SI{72(3)}{ps}$, which is significantly lower than the \SI{120}{ps} typically observed at $\SI{5}{K}$. These temperature-dependent phenomena are visible in the lifetime traces are subject of further investigation, as they could shed light on the changing dynamics of quasi particles in GaAs QDs at higher temperatures.

A full state tomography\cite{James2001} is performed to determine the degree of entanglement between the XX and the X photons. Figure \ref{fig:Entanglement}(d) shows an example among the 36 recorded XX/X coincidence histograms, corresponding to a measurement in the HV basis. The red dashed lines indicate the time-bin of \SI{2}{ns} in which the coincidences are summed up and compared with the average side peak area. The coincidences stemming from the slow X decay channel (the "tail" on the right side of the peaks), which are partially excluded by this time-bin, account for about \SI{9}{\percent} of the total coincidences per excitation cycle. Figure \ref{fig:Entanglement}(e) depicts the weighted sum of the histograms corresponding to the measurement bases HH, VV, HV and VH (with H the horizontal and V the vertical polarization), normalized by their respective average side peak areas. From this histogram a photon pair generation probability per excitation pulse of $\epsilon=0.87(2)$ is calculated (see supplementary for details), which is marginally lower than the value of $0.91(2)$ observed at \SI{5}{K}.

From the 36 correlation histograms the 2-qubit density matrix in polarization space is calculated and depicted in Fig. \ref{fig:Entanglement}(f). The maximum likelihood estimator used during this process is adapted for the dynamics of the QD light emission (see supplementary). The derived concurrence is $0.713(8)$ and the maximum fidelity to a Bell state is $0.925(5)$, which show a significant drop compared to the values obtained at \SI{5}{K}. However, the expected QBER of \SI{7.45}{\percent} calculated by Eq. \eqref{eq:qber} is still well below the $\SI{11}{\percent}$ allowed for the BBM92 protocol\cite{Bennett1992,Gisin2002}, and therefore suitable for performing QKD.

\begin{table*}[ht]
\caption{Summarized emitter performance for two representative QDs in a tunnel diode structure excited by TPE, measured at temperatures of \SI{5}{K} and \SI{20}{K}, respectively.}
\label{tab:performance}
\begin{center}
\begin{tabular}{|l||c|c||c|c|}
\hline
\rule[-1ex]{0pt}{3.5ex}  Temperature &  \multicolumn{2}{c||}{\SI{5}{K}} & \multicolumn{2}{c|}{\SI{20}{K}} \\
\hline
\rule[-1ex]{0pt}{3.5ex}  & X & XX & X & XX \\
\hline\hline
\rule[-1ex]{0pt}{3.5ex}  $g^{(2)}(0)$ & 0.017(4) & 0.011(3) & 0.020(3) & 0.034(4) \\
\hline
\rule[-1ex]{0pt}{3.5ex}  Lifetime (ps)  & 238(3) & 116(2) & 252(9) & 72(3)\\
\hline
\rule[-1ex]{0pt}{3.5ex}  Pair generation efficiency & \multicolumn{2}{c||}{0.91(2)} &  \multicolumn{2}{c|}{0.87(2)} \\
\hline
\rule[-1ex]{0pt}{3.5ex}  FSS (\SI{}{\micro eV}) & \multicolumn{2}{c||}{1.13(7)} &  \multicolumn{2}{c|}{0.96(9)}\\
\hline
\rule[-1ex]{0pt}{3.5ex}  Calculated Concurrence\textsuperscript{*} & \multicolumn{2}{c||}{0.905} &  \multicolumn{2}{c|}{0.900}\\
\hline
\rule[-1ex]{0pt}{3.5ex}  Measured Concurrence & \multicolumn{2}{c||}{0.904(3)} &  \multicolumn{2}{c|}{0.713(8)}\\
\hline
\rule[-1ex]{0pt}{3.5ex}  Calculated Fidelity to $\ket{\phi^+}$\textsuperscript{*} & \multicolumn{2}{c||}{0.959}&  \multicolumn{2}{c|}{0.960} \\
\hline
\rule[-1ex]{0pt}{3.5ex}  Measured Fidelity to $\ket{\phi^+}$ & \multicolumn{2}{c||}{0.975(1)}&  \multicolumn{2}{c|}{0.925(2)} \\
\hline
\multicolumn{5}{l}{\footnotesize{\textsuperscript{*}Only considering expectation values of measured $g^{(2)}$, X lifetime and FSS\cite{Hudson2007}}}
\end{tabular}
\end{center}
\end{table*}

\section{Quantum Key Distribution with Entangled Photons}

\begin{figure}
\begin{center}
\begin{tabular}{c}
\includegraphics[]{pics/F3_KeyGen.pdf}
\end{tabular}
\end{center}
\caption
{ \label{fig:KeyGen}
Key generation in the BBM92 protocol over a time span of about 8 hours. (a) QKD arangement. Alice and the photon source are situated on an optical table, Bob is placed in a movable box on a table in another building and connected with the source via a \SI{350}{m} long single mode fiber. (b) Qubit error rate (QBER) during the key generation with an average of \SI{8.42}{\percent}. (c) Raw (black) and secure (red) key rate during the key generation with average values of \SI{54.8}{bits/s} and \SI{14.2}{bits/s}, respectively. (d) Encryption of a bitmap with the dimensions of $144\times 150$ pixels and a color-depth of 24 bit, resulting in a total size of \SI{64.8}{kbyte}. The encryption with Alice's key yields a scrambled message ready to be sent over a public channel. (e) Decryption at Bob's site when using an uncorrected key (left) and a corrected key (right).}
\end{figure}

The XX and X photons in their individual single mode fibers are distributed to Alice and Bob, as sketched in Fig. \ref{fig:KeyGen}(a). The infrastructure is identical to the one used in a previous QKD experiment\cite{Schimpf2021}. After an initial synchronization and polarization correction routine, the key generation is performed overnight for a total duration of about 8 hours. The observed QBER, depicted in Fig. \ref{fig:KeyGen}(b), has an average value of \SI{8.42}{\percent}.

The QBER shows a minimum of \SI{6.15}{\percent}, which is even lower than the \SI{7.45}{\percent} estimated from $\rho$, shown in Fig. \ref{fig:Entanglement}(f). This discrepancy probably arises from the way the time-synchronization between Alice and Bob is maintained: In order to dynamically adjust the relative time delay between the arrival times of the XX and X photon at Alice and Bob, respectively, a peak in the continuously measured cross-correlation function is tracked. All photons around the peak-maximum in a time window of \SI{2}{ns} are used for key generation. We found that already a shift of the center of this time-bin by \SI{0.5}{ns} (about the time-resolution of our single photon detectors) can change the QBER by around \SI{2}{\percent}, because a higher contribution of photons stemming from the slow X decay channel, shown in Fig. \ref{fig:Entanglement}(c), decreases the entanglement (see supplementary for details). The time-window during the key generation was therefore probably shifted with respect to the one used for calculating $\rho$. These findings are not only interesting for a better understanding of the underlying mechanisms degrading the entanglement at higher temperatures, but also indicate that with careful choice of the time-window (within the limits of the detector resolution) a compromise between QBER and efficiency during the key generation process can be made.

We attribute the changing QBER over the course of the measurement to the following factors: The FSS changed from the initial \SI{0.96(9)}{\micro eV} to \SI{1.35(11)}{\micro eV} over the course of the QKD measurements, which can add about \SI{2}{\percent} to the QBER due to an accelerated phase rotation process leading to a decreased fidelity to the ideal Bell state on average\cite{Hudson2007}. A changing FSS was never noticed before in these kinds of samples. A possible origin could be a changing strain state in the silver glue used to stick the sample on the chip carrier. The second contribution to a rising QBER could be a varying ambient temperature, that affects the polarization rotation exerted by the fibers, which then leads to a larger deviation from the $\ket{\phi^+}$ Bell state, for which our QKD setup is designed. After the initial polarization correction no active correction was performed during the 8 hours of key generation.

Figure \ref{fig:KeyGen}(c) presents the raw (black) and secure (red) key generation rate. The secure key rate depends on the raw key rate and the QBER during the respective time window, as the QBER largely determines the efficiency of the key reconciliation process as described in Ref.\cite{Schimpf2021}. A change in raw key rate over time occurs due to a slight drift of the cryostat relative to the objective, leading to a decreasing excitation- and collection-efficiency. We counter this effect by an automatized movement of the X/Y position of the cryostat via linear stages to optimize the average detector count rates, should they fall below a certain threshold. The average raw key rate over the full time span is found to be \SI{54.8}{bits/s}, with a resulting average secure key rate of \SI{14.2}{bits/s}.

The generated (raw) key is used to encrypt a bitmap with a size of \SI{64.8}{kbyte}, see Fig. \ref{fig:KeyGen}(d), using a one-time-pad procedure. The decryption is depicted in Fig. \ref{fig:KeyGen}(e), for the cases when the raw key pair was used (left) and the corrected one (right).

\section{Discussion and Conclusion}

In this work, we have demonstrated a blinking-free source of polarization-entangled photon pairs, based on a GaAs quantum dot operated at a temperature of at least \SI{20}{K}. The intrinsically low fine structure splitting, owing to the local Al droplet etching technique\cite{Gurioli2019}, together with the employed p-i-n diode allows us to generate an uninterrupted stream of photon pairs with a fidelity to the $\ket{\phi^+}$ Bell state of \SI{0.925(2)} when using the pulsed two-photon-excitation scheme\cite{Hafenbrak2007,Muller2014}. The device also allows the fine-tuning of the emission wavelength within a range of \SI{0.2}{nm} while keeping the blinking-free operation intact, which is favourable for interconnecting multiple sources to quantum networks\cite{Reindl2017,Jons2017,Kimble2008,Gisin2007}.

The source was used to demonstrate quantum key distribution (QKD) via the BBM92 protocol\cite{Bennett1992} between two parties in two different buildings of the Johannes Kepler University, connected via a \SI{350}{m} long underground single mode fiber. The average qubit error rate (QBER) was \SI{8.42}{\percent}, which is well below the threshold of \SI{11}{\percent} required for generating secure keys when using the BBM92 protocol\cite{Bennett1992,Gisin2002}. From an initial key rate of \SI{54.8}{bits/s} error-free and privacy-amplified\cite{Schimpf2021} keys could be distilled with a rate of \SI{14.2}{bits/s}.

Comparing the decay dynamics of the biexciton and exciton states at temperatures of \SI{20}{K} and \SI{5}{K} allows us to identify a secondary slow decay channel of the exciton as the major entanglement degrading mechanism. While the physical origin of this observation and the elaboration of possible solutions will require further investigations, we find that the QBER during QKD can be optimized by a mild time filtering (see supplementary). This work makes us  optimistic that by combining electrical control with advanced photonic processing \cite{Liu2019,Wang2019} and strain-tuning platforms \cite{Trotta2015,Huber2018} will lead to nearly ideal sources of entangled photon pairs that can be operated in demanding environments.


\section{Introduction}
Colloidal dispersions exist in many shapes and forms, and understanding their behavior is relevant to biological,\cite{Don2003} ecological,\cite{Jackson2009} and industrial\cite{Alexopoulos2007} processes. The structures observed in such dispersions include for example amorphous colloidal glasses\cite{Weeks2012} and arrested colloidal gel networks.\cite{Pusey1995,Sandkuhler2004} The macroscopic behavior of colloidal dispersions is in large part determined by this spatial arrangement, and material properties such as the thermal conductivity and rheology vary strongly with the micro- and mesoscopic organization of the particles.\cite{Aksay1990,Lin2003,Egelhaaf2009,Burgi2010,Schurtenberger2013a} Thus, determining this spatial organization and understanding its link to the macroscopic behavior is a major point of interest within soft matter science.\cite{Dickinson2013,Wagner2017,Furst2019} For example, knowledge about the local structure of colloidal particles in gel networks\cite{Tanaka2008,Sprakel2017,Hermes2019} has led to the development of models that couple the network structure to macroscopic parameters such as the storage modulus.\cite{Morbidelli2001} Such models often require structural data that describe the gel network over a broad range of length scales, and microscopy and numerical simulations both offer powerful tools for probing the formation kinetics and structural properties of particle gels on this broad range of length scales.\cite{Solomon2003,Bartlett2005,Helgeson2015,Vermant2020,Walstra1995,Cates2002,Sandy2017,Zia2018} These methods are often combined with computational analyses that allow for the extraction of structural parameters,\cite{Grier1996,Weitz2001b} such as pair correlation functions, structure factors, and coordination numbers, that are particularly useful for describing the local network structure.\cite{Vermant2011,Solomon2012,Cohen2018,Simonsen2019} \\
Applying these analyses to experimental and computational datasets,\cite{Sciortino2008,Schmiedeberg2016,Furst2019} valuable insight regarding the local structure of such particle networks can be gained. However, such datasets also contain structural information on longer length scales, although a limited number of tools exist that can analyze the datasets reliably and efficiently.\cite{Dolors2007,Eiser2016} Existing analysis tools primarily focus on quantifying the material thickness and porosity of the network.
These analyses however require well-resolved individual network segments, which is not always available when analyzing microscopy data due to limitations in image resolution, especially with small particles.\cite{Egelhaaf2008,Xiaowei2009,Stenhammar2020}\\
In this Article, we will present a versatile computational method called ArGSLab, Arrested and Gelated Structures Laboratory, developed as a MATLAB package providing a set of tools to quantitatively analyze network structures in particle aggregates such as colloidal gels. ArGSLab performs a mesoscopic structural analysis that uses experimental or computational data to reconstruct three-dimensional gel networks\cite{Tanaka2019} as graphs of nodes and links.\cite{Stenhammar2020} ArGSLab furthermore enables quantification of the network structure using mesoscopic structural parameters such as node densities, branching densities and branch length statistics. As input, ArGSLab handles both image stacks from microscopy and sets of particle coordinates from either experiment or simulation. Crucially, the software does not require microscopy data containing well-resolved particle positions, enabling the analysis of structures formed by particles down to the same size range as the point-spread function. The software then converts the input data into single voxel-thin backbone structures, using routines based on established algorithms.\cite{Groen1980,Chu1994,Fratzl2017} The generated network structures are then transformed into a set of nodes and links for further statistical, graph theory-based analysis and visualization. ArGSLab can be accessed either via the provided graphical user interface (GUI) or directly through the MATLAB script. \\
In the following, we will first summarize the main features of ArGSLab and briefly discuss the overall workflow. Then, we will discuss the algorithms in some detail, and thereafter discuss the different output options. Finally, we will illustrate the capabilities of the software using three different data sets: (i) a set of microscopy images of colloidal gels with well-resolved particles, (ii) similar data but with poorly resolved individual particles, and (iii) coordinate sets from Brownian dynamics simulations.

\section{Software description}

\begin{figure}[ht]
	\centering
	\includegraphics[width=\linewidth]{Pipeline.png}
	\caption{Overview of the ArGSLab workflow. As input, the software uses either image stacks or particle coordinate sets, which are preprocessed into binarized images and then transformed into ``skeletons''. These go through a cleaning step, and are then transformed into a set of nodes and links for quantitative analysis. ArGSLab provides a histogram of the link lengths, the density of nodes and links in the system, and the tortuosity of the network.}
	\label{Pipeline}
\end{figure}

An overview of the general ArGSLab workflow is shown in Fig.~\ref{Pipeline}. The four major stages of the analysis are: (i) preprocessing, (ii) skeletonization, (iii) cleaning and (iv) quantitative analysis. In the preprocessing stage, coordinate sets are first transformed into image stacks, which are then successively Gaussian blurred, binarized and morphologically closed. Small unconnected structures are then removed, and small holes are filled. The skeletonization stage consists of the iterative erosion of voxels that belong to the material network, until a single-voxel thick backbone skeleton remains. This is then transformed into a graph of nodes and links. The cleaning stage removes artefacts from the skeletonization process, and the fourth and final stage contains a thorough analysis of the graph, producing statistics such as link length distribution, tortuosity and node density. In the following, we will describe these features.\\

\subsection{Algorithms}
\label{secalgo}
As input, the package accepts either .tiff image stacks from confocal microscopy or 3D particle coordinate sets from, \textit{e.g.}, computer simulations. Coordinate diagrams are first projected onto a .tiff image stack using user-defined voxel dimensions, with an optional convolution using a default or user-defined point spread function. The user-provided or generated image stacks are then treated with a Gaussian blur and binarized (\emph{i.e.}, transformed into images containing only black or white pixels). The purpose of treating images with a Gaussian blur is to reduce noise in the final skeleton, predominantly in the form of unrealistically small side chains, a known artefact of this type of skeletonization.\cite{Chu1994,Tougne2015} The binarization step is performed with a threshold individually calculated for each image slice. Determination of the threshold can be done either manually or using automated routines, where the automatic determination uses Otsu's method,\cite{Otsu1979} whereas the manual method determines, for each slice, the intensities $I_{10}$ and $I_{90}$, which are the pixel intensity values below which 10\% and 90\% of individual pixel intensities fall, respectively. A user-defined value $V_u$ between 0 and 1 then defines the threshold intensity $I_T$, with
\begin{equation}
I_T = I_{10} + V_u(I_{90}-I_{10}).
\end{equation}
The threshold determination is performed on each slice individually in order to compensate for differences in signal intensities between image slices: this is helpful when investigating microscopy data where a scattering length density mismatch leads to the average fluorescence intensity being dependent on the focal depth in the sample, or where fluorescence bleaching effects play a role.\cite{Egelhaaf2008} Images are then treated with a morphological closing step, consisting of subsequent dilation and erosion steps~\cite{Dougherty1992} performed with a user-defined ellipsoidal structural element. This procedure removes noise and effectively fills small gaps between particles, ensuring contiguity where necessary.\cite{Devaux2011} Structures are connected should their binarized shapes be closer than twice the structural element radius. Subsequently, unconnected clusters smaller than a threshold fraction are removed, with larger unconnected clusters being analyzed as separate network structures. Finally, holes in the structure that are enclosed on all sides in 3D and are smaller than a user-defined threshold value are detected and filled in, reducing unphysical holes in the skeleton. The sizes of the Gaussian blur kernel, the morphological closing element and the minimal cluster size take either default or user-defined values.\\
Subsequently, the package performs the skeletonization, using algorithms adapted from \citet{Fratzl2017} This process consecutively removes material voxels bordering non-material voxels, with the condition that they are not end points, non-Euler invariant points, and are Euler simple points. Removal of a voxel considers the 3$\times$3 voxel cube $D_3(i)$ surrounding voxel $i$, \textit{i.e.} a cube of 27 voxels with voxel $i$ at the center. An end point is defined as the case where $D_3(i)$ contains in total two material voxels, of which one is $i$. Euler invariant points are defined as voxels that, upon removal, do not lead to a change of the Euler characteristic of the total structure.\cite{Pudney1998,Hansen2010,Baja2016} We define an Euler simple point as a voxel $i$ whose removal neither changes the connectivity nor the Euler characteristic of $D_3(i)$;\cite{Pudney1998,Omar2013} a more detailed discussion of these definitions can be found in \citet{Chu1994} The transformation of a skeleton into a graph of nodes and links is then performed by assuming that (i) each voxel $i$ whose surrounding $D_3(i)$ contains exactly 2 material voxels is an end node, (ii) each voxel $i$ whose $D_3(i)$ contains exactly 3 material voxels belongs to a link and (iii) each voxel $i$ whose $D_3(i)$ contains more than 3 material voxels is a branching node.\\
While the skeletonization process robustly extracts the relevant backbone structure, several known unrealistic features can arise, often due to noisy or irregular data.\cite{Chu1994} The subsequent cleaning steps remove such unwanted features, and are performed iteratively using the following procedures:
\begin{enumerate}
\item Unphysically small side chains, arising from irregular structure surfaces,\cite{Chu1994,Chen2018} are excluded by removing all side chains smaller than a certain threshold value, by default set to 1.5$\sigma$ with $\sigma$ the particle diameter.
\item The removal of \textit{canal nodes}, \textit{i.e.} nodes that are connected to exactly two links. Such nodes and their connecting links are effectively one single link and therefore merged.
\item Nodes that are in close proximity are collected into a single node, as dense sections of gels with irregular edges often lead to complex and unphysical skeleton sections with many interconnected nodes and links. A set of nodes to be collected contains all nodes that are closer than a threshold value to another node in the same set. The threshold value is by default set to 0.9$\sigma$.
\end{enumerate}
After these procedures, a new skeleton is formed from the cleaned graph.\\

\begin{figure}[ht]
	\centering
	\includegraphics[width=.5\linewidth]{vis_results.png}
	\caption{Visualized output. A) slices of the original image .tiff stack overlaid with the skeleton (green) and branching nodes (blue). Since these are 2D slices of a 3D image, interrupted skeleton strands continue in adjacent $z$-slices. B) 3D rendering of the same skeleton.}
	\label{vis_res}
\end{figure}

\subsection{Visualization}
\label{secvis}
ArGSLab provides output that allows the user to visualize and quantify the analyzed network structures. The software also provides output that helps optimizing the parameters of the skeletonization process. The standard visualizations are: (i) the original image stack overlaid with the calculated skeleton and branching nodes, as shown in Fig.~\ref{vis_res}A (data from \citet{Schmiedeberg2016}), and (ii) a fully rotatable 3D rendering of the skeleton backbone, as shown in Fig.~\ref{vis_res}B.\cite{Law2005,Yeh2020}\\
\begin{figure}[t]
	\centering
	\includegraphics[width=\linewidth]{Optimization.png}
	\caption{Slices from 3D input and a selection of ArGSLab's output image stacks. A) Raw input image. B) Raw input image overlaid with a green and red color filter, produced in the binarization step. Green indicates voxels determined to contain material (material voxels), and red to contain no material (non-material voxels). Black voxels had little to no intensity in the original image. C) Preprocessing effects after binarization: white sections are material voxels, unaltered by the preprocessing, green sections are material voxels added by preprocessing, and red sections are previous material voxels that have been removed. D) Overlay of the raw data with the final skeleton and nodes.}
	\label{Optimization}
\end{figure}As described in Section~\ref{secalgo}, several parameters and thresholds can be changed from their default values to optimize the output. The package provides several types of visualization to help optimize those parameters, a selection of which are shown in Fig.~\ref{Optimization}. Starting from a raw image stack (Fig.~\ref{Optimization}A), a new image stack is output after the binarization step (Fig.~\ref{Optimization}B). This is intended to help setting the binarization threshold, should the automatic method be insufficient. Fig.~\ref{Optimization}C shows the effects of further preprocessing after binarization, \textit{i.e.} the morphological closing, hole filling, and the removal of smaller unconnected structures. The final results are shown in Fig.~\ref{Optimization}D, showing an overlay of the raw image with the final skeleton.\\

\subsection{Statistical analysis}
From the graph of branching nodes and connecting links, a number of statistical quantities can be extracted, such as the number and density of nodes and links. The total number of nodes and links $N_N$ and $N_L$ are calculated together with the average number of links per node $\frac{N_L}{N_N}$. The density of nodes and links $\rho_N$ and $\rho_L$ respectively are given in units of ${\mu}\text{m}^{-3}$ and $\sigma^{-3}$ with $\sigma$ the particle diameter. Furthermore, ArGSLab provides a histogram showing the distribution of link lengths in the network, presenting the normalized link count $\tilde{N}(\Lambda)$ as a function of link length $\Lambda$, which we define as the length of a link between two nodes as measured along the link path. $\tilde{N}(\Lambda)$ is related to the fractal dimension $d_f$, a parameter often used for characterizing aggregated systems:\cite{Aksay1990,Morbidelli2003,Morbidelli2016} a propensity for long links and low $d_f$ are both characteristic of open aggregates. Another output measure is the tortuosity $\xi$, which is defined as the average trajectory length following a backbone strand $\lambda$, divided by the calculated Euclidean distance $\lambda_{Euc}$ between their start and end points A and B:
\begin{equation}
\xi = \left\langle \frac{\lambda(\text{A,B})}{\lambda_{Euc}(\text{A,B})} \right\rangle,
\end{equation}
where A and B are the start and end points of paths along the backbone intersecting respectively the upper and the lower box face in either Cartesian direction, and the angular brackets denote averaging over all such paths for a given configuration. $\xi$ provides a measure of the strand erraticity,\cite{Gladden1992,Liu2006,Harvey2018,Stenhammar2020} and is often used in catalysis,\cite{Hu2017} blood vessel analysis,\cite{Tomaszewski2018} or satellite geography.\cite{Harvey2018} All quantitative analyses are corrected for box edge effects. Furthermore, additional analyses beyond the ones described here are facilitated by ArGSLab as the skeleton backbone and the node and link lists can be exported as .mat or .txt files, after which analyses tailored to the specific system can be applied.\\

\begin{figure}
	\centering
	\includegraphics[width=.6\linewidth]{GUI.png}
	\caption{Screenshot of the GUI, with fine-tuning tools enabled.}
	\label{GUI}
\end{figure}

\subsection{Graphical User Interface}
The package can be accessed either using the provided MATLAB wrapper code, or using the graphical user interface (GUI), shown in Fig.~\ref{GUI}. The GUI gives intuitive access to the full capabilities of the code and requires only a few inputs: the folder path containing the input data, the number of $z$-slices to be analyzed in each file (in case not all slices should be analyzed), the voxel dimensions, and the average particle diameter. Ticking the ``fine-tuning options'' box reveals several further options that allow for optimization of output quality and can typically be left at their default values. These are described in-depth in the GUI tooltip menus and the package manual. The GUI is designed to be an intuitive tool for using ArGSLab, allowing users to run the full analysis without loss of functionality as compared to the wrapper script.\\

\section{Examples}
In the following we will illustrate the use of ArGSLab by showing the processing of three example datasets. Dataset 1 consists of microscopy images of a set of colloidal gels made from PMMA spheres with an added depletant that causes attractive forces leading to gel formation. The gels are imaged using confocal laser scanning microscopy (CLSM) and the large particles allow for individual particles to be well-resolved.\cite{Schmiedeberg2016} The results from analyzing this dataset are compared to the corresponding results from coordinate sets of the particle centers, extracted using routines following Crocker and Grier.\cite{Grier1996} Dataset 2 contains CLSM images of gels of smaller, pNIPAm microgel particles (as described in Section~\ref{secmatmet}) that cannot be individually resolved. Dataset 3 is a set of particle coordinates obtained from Brownian dynamics simulations of colloidal gels formed from colloidal particles with attractive patches icosahedrally arranged on the particle surface.\cite{Stenhammar2020} Further details about the preparation can be found in Section~\ref{secmatmet}, and the results are summarized in Fig.~\ref{HSExample} and Table~\ref{HSExampleTab}.\\
Figs.~\ref{HSExample}A-D show the results for Dataset 1. The quantitative analyses in Fig.~\ref{HSExample}D and Table~\ref{HSExampleTab} reflect the similarity between experimental skeletons and the ones obtained from the corresponding coordinate data. The small discrepancies between the results can be explained as follows: the coordinate extraction routines identify the vast majority of particles, but atypical particles, due to synthetic or imaging anomalies, are sometimes not picked up. While this only leads to negligible errors in tracking experiments, ArGSLab's structural analysis identifies this as broken links. Such effects can also arise from the determination of a particle center slightly off from its real particle center due to experimental noise. ArGSLab can partially accommodate such anomalies by choosing larger structural element sizes. \\
The results from Dataset 2 (Figs.~\ref{HSExample}E-H and Table~\ref{HSExampleTab}) clearly demonstrate the software's ability of extracting quantitative data from structures constructed from non-resolvable particles. Nonetheless, the lower data quality often leads to the need for a more careful selection of optimization parameters. Skeletonization of data with low contrast, low signal-to-noise ratio, or poor $z$-resolution is known to lead to skeletons with small, unphysical side chains, or with a multitude of nodes in dense cluster-like sections.\cite{Chu1994} Tuning the parameters of the cleaning routines can alleviate these types of noise: larger Gaussian blurring kernels or morphological closing structural elements often cause dense clusters to be correctly interpreted as closed entities. Furthermore, the node collection threshold in the second cleaning step can be increased, collecting these virtual nodes into a single node. Furthermore, the terminal branch length threshold (third cleaning step in Section~\ref{secalgo}) can be raised to remove unphysical side chains. Changing these parameters can affect the number of nodes and links found, and should be kept constant between datasets for full comparability. In all our examples we have used the default parameters, which provides optimal comparability between different experiments.\\
One important consideration when setting the analysis parameters is that determining whether two particles are to be considered in contact or not is a non-trivial question\cite{Makse2007,Egelhaaf2011a,Hsiao2020} which will affect the final skeleton. Whether ArGSLab considers two segments in contact or not is implicitly governed by the binarization and the morphological closing steps. In brief, two gel segments will be considered in contact if their material voxels after binarization are closer together than the diameter of the morphological closing element. This choice is rationalized by considering that colloidal gels require a significant attractive potential, and intersegment gaps are often metastable, while noisy images often lead to apparent gaps. Optimal skeleton output therefore requires a careful tuning of the morphological closing element size and the binarization threshold, with the preselected standard parameters as good starting points.\\
The results from Dataset 3 (Figs.~\ref{HSExample}I-L and Table~\ref{HSExampleTab}) further demonstrate that ArGSLab offers a direct comparison between experimental systems and model systems studied by simulations. The large difference in the node and link densities $\rho_N$ and $\rho_L$ between Dataset 1 and Datasets 2 and 3 results from the much higher volume fraction in the former system ($\phi = 0.2$ compared to $\phi = 0.05$) that leads to much higher node and link densities in Dataset 1. The difference in $\phi$ also results in a significantly lower $\tilde{N}(\Lambda)$ at high $\Lambda$ for Dataset 1, and the tortuosity $\xi$ is also lowered by the higher $\rho_N$. It is also apparent that the gels in Datasets 2 and 3 have a stronger propensity for having end nodes rather than branching nodes from the lower $N_L/N_N$ ratios in Table~\ref{HSExampleTab}.\\

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{combined_results.png}
	\caption{Analysed data from the three example datasets: A-D) Dataset 1, corresponding to dataset C5 in Ref.~\citenum{Schmiedeberg2016}. E-H) Dataset 2, containing particle gels formed by attractive microgel particles, synthesized as in Ref.~\citenum{Schurtenberger2020} I-L) Brownian dynamics simulation data of particles with patchy attractions, from Ref.~\citenum{Stenhammar2020}. Panels (A,E,I) show image sections from the raw data, panels (B,F,J) overlays of the same image and the final skeleton, panels (C,G,K) 3D renderings of the full skeleton (green) and  branching nodes (blue), and panels (D,H,L) the normalized link length distribution $\tilde{N}(\Lambda)$. The two datasets in panel (D) corresponds to $\tilde{N}(\Lambda)$ obtained using either the raw image data, or from a particle coordinate set extracted from the same data, demonstrating the equivalence of the two approaches. All 3D skeleton renderings correspond to a volume of approximately 1.7 $\times 10^4 \sigma^3$, and the scale bars represent 10$\sigma$.}
	\label{HSExample}
\end{figure}

\begin{table*}
\centering
\caption{Values from statistical analysis of the example datasets in Fig.~\ref{HSExample}.}
\label{HSExampleTab}
\begin{tabular}{ l | l | l | l | l }
Type & $\rho_N$ [$\sigma^{-3}$] & $\rho_L$ [$\sigma^{-3}$] & $N_L/N_N$ & $\xi$\\
\hline
\hline
Dataset 1 (images) & $7.86 \times 10^{-2}$ & $1.34 \times 10^{-1}$ & $1.70$ & $2.03 \pm 0.14$ \\
Dataset 1 (coordinates) & $8.36 \times 10^{-2}$ & $1.44 \times 10^{-1}$ & $1.72$ & $2.05 \pm 0.09$ \\
\hline
Dataset 2 & $7.18 \times 10^{-3}$ & $8.67 \times 10^{-3}$ & $1.22$ & $2.23 \pm 0.12$ \\
\hline
Dataset 3 & $9.10 \times 10^{-3}$ & $1.02 \times 10^{-2}$ & $1.13$ & $2.25 \pm 0.13$ \\
\end{tabular}
\end{table*}

\section{Conclusions and Perspectives}
In this Paper, we have demonstrated the capabilities of ArGSLab, a tool for analyzing networks of aggregated particulate matter, such as colloidal gels. We have shown that ArGSLab is able to analyze both microscopy and simulation data, and allows for direct comparison between different datasets. The software provides a quantitative analysis of the mesoscopic network structure, which is more difficult to analyze than the local particle arrangements. Importantly, ArGSLab also provides a method to quantify systems where individual particles are not resolvable and the local structure thus cannot easily be probed. The graph theory approach that our package applies enables comprehensive analyses of aggregated structures, and it yields several parameters allowing for a unique insight into their structure. Moreover, ArGSLab can be used to transform experimental or simulation data into a set of nodes and links for further analysis, and thus provides powerful information for the development of new analyses tailored to user-specific needs. Since the method is developed to analyze any network that is composed of smaller segments, ArGSLab can be applied to networks beyond colloidal gels, such as protein networks,\cite{Becker2011} bacterial aggregates in biofilms\cite{Thomas2018} and biological cell networks.\cite{Hilfiker2019} By providing for an easy means to extract global network parameters, ArGSLab provides a powerful method of extracting mesoscopic structural parameters that can be related to microscopic and macroscopic properties, for instance obtained from light scattering or rheology.

\section{Materials and Methods}
\label{secmatmet}
Dataset 1 is the same as dataset C5 in Ref.~\citenum{Schmiedeberg2016}. Particles are PMMA spheres with a diameter $\sigma = 1.72 \mu$m, aggregated due to depletion forces. We analyze 30 stacks with total imaged volume $8.85\times10^{4} \mu\text{m}^{3}$ or $1.72\times10^4 \sigma^{3}$ per stack, voxel dimensions of 0.061$\times$0.061$\times$0.12 $\sigma^3$, and a volume fraction $\phi$ of approximately $0.2$. Further details can be found in Ref.~\citenum{Schmiedeberg2016}.\\
Dataset 2 consists of polystyrene/poly-N-Isopropylacrylamide core-shell microgels, synthesized using a method described in Ref.~\citenum{Schurtenberger2020}. In collapsed form, $\sigma = 398$nm and particles aggregate irreversibly due to a combination of van der Waals interactions and other short-range attractions. We analyzed 40 image stacks with a total imaged volume per stack of $1.08\times10^{4} \mu\text{m}^{3}$ or $1.71\times10^5 \sigma^{3}$, voxel dimensions of 0.054$\times$0.054$\times$0.12 $\sigma^3$, with $\phi = 0.05$. The samples are prepared and imaged using identical methods as described in Ref.~\citenum{Stenhammar2020}.\\
Dataset 3 was obtained from Brownian dynamics using particles with attractive patches icosahedrally placed on the surface as described in Ref.~\citenum{Stenhammar2020}. We analyze the final structures from four different Brownian dynamics simulations. To compensate for local fluctuations, the results from 10 temporally close coordinate diagrams were averaged for each simulation. The total system volume was $1.0\times10^5 \sigma^{3}$ with $\phi = 0.05$. Further details can be found in Ref.~\citenum{Stenhammar2020}.\\

\section*{Computational Efficiency}
On an Intel\textsuperscript{\textregistered} octacore i7-7700 CPU with 3.60 GHz running Windows 10 and MATLAB 2021a, extracting a skeleton and subsequent analysis from one image from dataset 1 (512x512x151 pixels, $\phi = 0.2$) takes approximately 5 minutes. Starting from the coordinate diagram only adds a few seconds to this calculation time. Additional visualization and optimization visualization adds 6 minutes. The code has been tested on iOS, Windows and Linux platforms, with comparable computation times. Decreasing image size and $\phi$ reduces computational time.

\section*{Resources}
A project page for the package exists at GitHub.com (github.com/jimmink/argslab), where the standalone MATLAB package and GUI for Windows, iOS and Linux are freely available, together with a user manual and example datasets. An idealized, computer-generated image stack is also available, intended to help users familiarize themselves with ArGSLab and its parameters. A walkthrough of the analysis of this image is available as Supplementary Information.

\section*{Acknowledgements}
We gratefully acknowledge financial support from the Alexander von Humboldt Foundation (JNI), the European Research Council (ERC-339678-COMPASS (PS)) and the Swedish Research Council (Grant numbers 2018-04627 (PS) and 2019-03718 (JS)).

\section*{Supplementary Information}
A walkthrough of the ArGSLab analysis of the freely available benchmarking image.

\footnotesize{