\section{Introduction}

Humans are lifelong learners: in order to function effectively day-to-day, we acquire and accumulate knowledge throughout our lives. The accumulated knowledge makes us efficient and versatile when we encounter new tasks. In contrast to human learning, modern neural network based learning algorithms usually fail to remember knowledge acquired from previous tasks when adapting to a new task (see Figure~\ref{fig:intro} (b)). In other words, it is difficult to generalize once a model is trained on a task. This is the well-known phenomenon of \textit{catastrophic forgetting}~\cite{mccloskey1989catastrophic}. Recent efforts~\cite{shmelkov2017incremental,li2017learning,chaudhry2018efficient,aljundi2018memory} have demonstrated how discriminative models can continually learn a sequence of tasks. Despite the success of these efforts for discriminative models, lifelong learning for generative models remains a challenging and under-explored area.

Lifelong learning methods for discriminative models cannot be directly applied to generative models due to their intrinsic differences. \textit{First}, it is well known that in classification tasks,  the intermediate convolutional layers in deep neural networks are capable of providing generic features.
These features can easily be reused by other models with varying and different classification goals. For generative models, the possibility of such reuse for new generative tasks has not been previously explored, to the best of our knowledge. \textit{Second}, different from discriminative models, the output space of generative models is usually continuous, making it more challenging for the model to maintain the generation quality along the training of a sequence of tasks. \textit{Third}, there could be conflicts between tasks under the generative setting. For discriminative models, it rarely happens that one image has different labels (appears in different tasks). However, for generative models, it is quite common that we want to translate the same image to different domains for different tasks~\cite{isola2016pix2pix,gatys2016image,zhu2017unpaired,choi2018stargan}.

One of the naive approaches to address catastrophic forgetting for generative models is to train a model for each task separately (see Figure~\ref{fig:intro} (c)). Unfortunately, this approach is not scalable in general: as new tasks are added, the storage requirements grow drastically. More importantly, setting the trained model aside without exploiting the benefit it can potentially provide facing a new task would be an inefficient use of resources.
Recent works~\cite{zhai2019lifelong,wu2018memory} have shown promising results on lifelong learning for generative models, but it is also revealed that image quality degradation and artifacts transfer from old to new tasks are inevitable. Therefore, it is valuable to have a continual learning framework designed that (1) is more parameter efficient, (2) preserves the generation quality of both current and previously learned tasks, and (3) can enable various conditional generation tasks across different domains.

In this paper, we introduce a generic continual learning framework \textit{Piggyback GAN} that can perform various conditional image generation tasks across different domains (see Figure~\ref{fig:intro} (d)). The proposed approach addresses the catastrophic forgetting problem in generative continual learning models. Specifically, \textit{Piggyback GAN} maintains a \textit{filter bank}, the filters in which come from convolution and deconvolution layers of models trained on previous tasks. Facing a new task, Piggyback GAN learns to perform this task by reusing the filters in the filter bank by building a set of \textit{piggyback filters} which are factorized into filters from the filter bank. Piggyback GAN also maintains a small portion of unconstrained filters for each task to ensure high-quality generation. Once the new task is learned, the unconstrained filters are appended to the filter bank to facilitate the training of subsequent tasks. Since the filters for the old task have not been altered, Piggyback GAN is capable of keeping the exact knowledge of previous tasks without forgetting any details.

To summarize, our contributions are as follows. We propose a continual learning framework for conditional image generation that is (1) \textit{efficient}, learning to perform a new task by ``piggybacking" on models trained on previous tasks and reusing the filters from previous models, (2) \textit{quality-preserving}, maintaining the generation quality of current task and ensuring no quality degradation for previous tasks while saving more parameters, and (3) \textit{generic}, enabling various conditional generation tasks across different domains. To the best of our knowledge, we are the first to make these contributions for generative continual learning models. We validate the effectiveness of our approach under two settings: (1) paired image-conditioned generation, and (2) unpaired image-conditioned generation.
Extensive qualitative and quantitative comparisons with state-of-the-art models are carried out
to illustrate the capability and efficiency of our proposed framework to learn new generation tasks without the catastrophic forgetting of previous tasks.

\section{Method}

\begin{figure*}[ht]
\begin{center}
  \includegraphics[width=\textwidth]{figures/model/overview-factorization.pdf}
\end{center}
\caption{\textbf{Overview of Piggyback GAN}. Piggyback GAN maintains a \textit{filter bank} for composing convolutional filters to be used for the next task. Given a new task, most filters in the model are trained as \textit{piggyback filters} that are factorized into filters in the filter bank. The remaining small portion of filters are learned without constraints and appended to the filter bank to benefit the training of subsequent tasks.}
\label{fig:overview}
\end{figure*}

In this paper, we study the problem of lifelong learning for generative models. The overall goal is to learn the model for a given sequence of generation tasks $\{T_i\}_{i=1}^N$ under the assumption that the learning process has access to only one task at a time and training data for each task is accessible only once.

Assume that we have a learnt model $M$ trained to perform the first task $T_1$. To proceed with the new task $T_2$, one naive approach is to adapt the learnt model $M$ to the new task by fine-tuning its parameters. This approach suffers from the problem of catastrophic forgetting of task $T_1$. Another naive approach is to learn a separate standalone model $M_i$ for each new task $T_i$ and store all learnt models $\{M_i\}_{i=1}^N$, which drastically scales up total number of parameters without utilizing previous models.

The key idea of this paper is to build a parameter efficient lifelong generative model by ``piggybacking" on existing models and reusing the filters from models trained on previous tasks while maintaining the generation quality to be similar to a single standalone model for all tasks. We achieve this by maintaining a filter bank, and compose the filters for new tasks mostly by factorizing it into filters in the filter bank, with a small portion of filters learnt without constraints. Once the new task is learnt, the filter bank is expanded by adding those unconstrained filters.

\subsection{Piggyback Filter Learning}
\label{sec:piggyback-filter}
For lifelong learning of generative models, given an initial task $T_1$ and an upcoming new task $T_2$, there could be cases where $T_1$ and $T_2$ share similar inputs but different outputs, e.g. $T_1$ is \textit{Photo $\rightarrow$ Monet Paintings} and $T_2$ is \textit{Photo $\rightarrow$ Ukiyo-e Paintings}. These conflicts in input-output space do not occur when each task is trained separately. However, as already mentioned, learning a separate model for each task and storing all learnt models is memory consuming and does not utilize previous models.

For this scenario, the knowledge in the filters from the model trained for $T_1$ could provide valuable information for the training of $T_2$.
For images from domains that are less visually similar, it is well known that in discriminative tasks, the generaliziblity of convolutional neural networks helps maintain its capability in previously unseen domains~\cite{mallyaetal2018,rebuffietal2018}. Such evidence suggests that despite the difference in fine details, certain general patterns that convey semantics could still be covered in previous filter space.
Inspired by these observations, we propose to utilize previously trained filters to facilitate the training of new tasks. In our experiments, we do find that previous filters make steady contribution to learning new tasks regardless of the domain difference across tasks.

Consider a generator network $G_1$ learnt for initial task $T_1$ consisting of $L$ layers of filters $\{F_1^{\ell}\}_{\ell=1}^L$ where $\ell$ denotes the index of layers. For the ${\ell}^{th}$ layer, let the kernel size be $s_w^{\ell} \times s_h^{\ell}$, number of input channels be $c_{in}^{\ell}$ and number of output channels be $c_{out}^{\ell}$. Then $F^{\ell}_1$
is reshaped from 4D into 2D of size $s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell}$ by $c_{out}^{\ell}$, which we denote by $\mathcal{R}(F^{\ell}_1)$. For the upcoming task $T_2$, we learn the filters $F^{\ell}_2$ through factorization operation over filters of the corresponding layers in task $T_1$, namely

\begin{equation}
\begin{split}
    F^{\ell}_{2} \, = \, \mathcal{R}^{-1}(\mathcal{R}(F^{\ell}_{1}) \otimes W^{\ell}_2),
\end{split}
\end{equation}
where $\otimes$ denotes the standard matrix multiplication operation and $\mathcal{R}^{-1}$ denotes the inverse of reshape operation from 2D back to 4D. The derived filters $F^{\ell}_{2}$ is denoted as \textit{piggyback filters} and the corresponding learnable weights $W^{\ell}_2 \in \mathbb{R}^{c_{out}^{\ell} \times c_{out}^{\ell}}$ is denoted as \textit{piggyback weight matrix}. The resulting $F^{\ell}_{2} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times c_{out}^{\ell}}$ is of the same size as $F^{\ell}_{1}$. A bias term could be added to the factorized filters to adjust the final output.

Therefore, by constructing factorized filters, the number of parameters required to be learnt for this layer is reduced from $s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times c_{out}^{\ell}$ to only $c_{out}^{\ell} \times c_{out}^{\ell}$.

\subsection{Unconstrained Filter Learning}
\label{sec:unconstrained-filter}
Parameterizing a new task completely by making use of filters in the initial tasks, though saving substantial storage, cannot capture the key differences between tasks. The filters derived from previous tasks may only characterize certain attributes and may not generalize well to the new task, resulting in poor quality of generation details for the new model. Learning a small number of unconstrained filters gives the model greater flexibility in adapting to the new task. Moreover, learning these unconstrained filters could learn to generate patterns that do not exist in previous tasks, thereby, increasing the power of the model and helping the training of later tasks.

Specifically, for each layer in the generator, we allocate a small portion of unconstrained filters which are learned freely without the constraint to be constructed from previous filters. With a total of $c^{\ell}_{out}$ filters in $F^{\ell}_2$, consider the proportion of unconstrained filters to be $\lambda \, (0 < \lambda \leq 1)$\footnote{$\lambda c^{\ell}_{out}$ could be rounded to the nearest integer. Or $\lambda$ could be chosen to make $\lambda c^{\ell}_{out}$ an integer.}, then we have $\lambda c^{\ell}_{out}$ unconstrained filters and $(1-\lambda) c^{\ell}_{out}$ piggyback filters, and the size of \textit{piggyback weight matrix} $W^{\ell}_2$ has changed to $c_{out}^{\ell} \times (1-\lambda)c_{out}^{\ell}$. Let the unconstrained filters at the $l^{th}$ layer of the $n^{th}$ task $T_n$ be $F_n^{u,\ell}$, and piggyback filters at the $l^{th}$ layer of the $n^{th}$ task $T_n$ be $F_n^{p,\ell}$.

For task $T_1$, all filters are unconstrained filters, namely
\begin{equation}
    F^{\ell}_{1} \, = \, F^{u,\ell}_{1}.
\end{equation}

For task $T_2$, we re-define the filters in $F^2_l$ to be the concatenations of unconstrained filters $F_2^{u,\ell}$ and piggyback filters $F_2^{p,\ell}$, namely $F_2^{\ell}$ is formulated as

\begin{equation}
\begin{split}
    F^{\ell}_{2} \, & = \, [F_2^{u,\ell}, F_2^{p,\ell}] \\
              & = \, [F_2^{u,\ell}, \mathcal{R}^{-1}(\mathcal{R}(F^{u,\ell}_{1}) \otimes W^{\ell}_2)],
\end{split}
\end{equation}
where $F_2^{u,\ell} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times \lambda c_{out}^{\ell}}$, $F_2^{p,\ell} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times (1-\lambda) c_{out}^{\ell}}$ and the resulting $F^{\ell}_{2} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times c_{out}^{\ell}}$ is of the same size as $F^{\ell}_{1}$.

\subsection{Expanding Filter Bank}
\label{sec:filter-bank}

The introduction of unconstrained filters learnt from task $T_2$ brings extra storage requirement and it would be a waste if we set it aside when learning the filters for the upcoming new tasks. Therefore, instead of setting it aside, we construct the piggyback filters for task $T_3$ by making use of the unconstrained filters from both $T_1$ and  $T_2$.

We refer to the full set of unconstrained filters as the \textit{filter bank} which will be expanded every time new unconstrained filters are learnt for a new task. Expanding the \textit{filter bank} could encode more diverse patterns that were not captured or did not exist in previous tasks. By doing so, the model increases the modality of the learned filters and provides more useful information for the subsequent tasks.

When learning the task $T_3$, the filter bank comprises filters corresponding to the task $T_1$ and $T_2$ and is given by $F^{u,\ell}_{1}$ and $F^{u,\ell}_{2}$, resulting in the size of the filter bank $s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times (1+\lambda) c_{out}^{\ell}$. Similarly, for task $T_3$ consider the proportion of unconstrained filters to be same as $T_2$, i.e. $\lambda$, the size of \textit{piggyback weight matrix} $W^{\ell}_3$ would be $(1+\lambda)c_{out}^{\ell} \times (1-\lambda)c_{out}^{\ell}$. We construct $F^{\ell}_{3}$ as the concatenation of unconstrained filters $F^{u,\ell}_{3}$ and piggyback filters $F^{p,\ell}_{3}$ that are constructed from all filters in the \textit{filter bank}, namely $F_3^{\ell}$ is formulated as

\begin{equation}
\begin{split}
    F^{\ell}_{3} \, & = \, [F^{u,\ell}_{3}, F^{p,\ell}_{3}] \\
               & = \, [F^{u,\ell}_{3}, \mathcal{R}^{-1}(\mathcal{R}([F^{u,\ell}_{1}, F^{u,\ell}_{2}]) \otimes W^{\ell}_3) \,],
\end{split}
\end{equation}
where $F_3^{u,\ell} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times \lambda c_{out}^{\ell}}$,  $F_3^{p,\ell} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times (1-\lambda) c_{out}^{\ell}}$ and the resulting $F^{\ell}_{3} \in \mathbb{R}^{s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times c_{out}^{\ell}}$ is of the same size as $F^{\ell}_{1}$ and $F^{\ell}_{2}$.

After task $T_3$ is learnt, the \textit{filter bank} is expanded as
$[F^{u,\ell}_{1}, F^{u,\ell}_{2}, F^{u,\ell}_{3}]$, whose size is $s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times (1+2\lambda) c_{out}^{\ell}$.

To summarize, when learning task $T_{n}$, $F^{\ell}_{n}$ could be written as
\begin{equation}
\begin{split}
    F^{\ell}_{n} \, & = \, [F^{u,\ell}_{n}, F^{p,\ell}_{n}] \\
               & = \, [F^{u,\ell}_{n}, \mathcal{R}^{-1}(\mathcal{R}([F^{u,\ell}_{1}, F^{u,\ell}_{2},...,F^{u,\ell}_{n-1}]) \otimes W^{\ell}_n) \,].
\end{split}
\end{equation}

After task $T_{n}$ is learnt, the \textit{filter bank} is expanded as $[F^{u,\ell}_{1}, F^{u,\ell}_{2}, ..., F^{u,\ell}_{n}]$, whose size is $s_w^{\ell} \times s_h^{\ell} \times c_{in}^{\ell} \times (1+(n-1)\lambda) c_{out}^{\ell}$, and the size of \textit{piggyback weight matrix} $W^{\ell}_{n+1}$ for task $T_{n+1}$ is $(1+(n-1)\lambda)c_{out}^{\ell} \times (1-\lambda)c_{out}^{\ell}$. It should be noted that the weights of filters in the \textit{filter bank} remain fixed along the whole learning process.

\subsection{Learning Piggyback GAN}
We explore two conditional generation scenarios in this paper:
(1) \textit{paired image-conditioned generation}, in which training data contains $M$ pairs of samples \textit{$\{(a_i, b_i)\}_{i=1}^{M}$}, where \textit{$\{a_i\}_{i=1}^{M}$} represent conditional images, \textit{$\{b_i\}_{i=1}^{M}$} represent target images, and correspondence between $a_i$ and $b_i$ exists, namely for any conditional image $a_i$ the corresponding target image $b_i$ is also provided; (2) \textit{unpaired image-conditioned generation}, in which training data contains images from two domains $A$ and $B$, namely images \textit{$\{a_i\}_{i=1}^{M_a} \in A$} and images \textit{$\{b_i\}_{i=1}^{M_b} \in B$}, and correspondence between $a_i$ and $b_i$ does not exist, namely for any conditional image $a_i$ the corresponding target image $b_i$ is not provided.

For both conditional generation scenarios, given a sequence of tasks and a state-of-the-art GAN model, we construct the convolutional and deconvolutional filters in the generator as described in Sec.~\ref{sec:piggyback-filter}, Sec.~\ref{sec:unconstrained-filter} and Sec.~\ref{sec:filter-bank}. The derived filters is updated by the standard gradient descent algorithm and the overall Piggyback GAN model is trained in the same way as any other existing generative models by adopting the desired learning objective for each task.

\section{Conclusion}
We proposed Piggyback GAN for lifelong learning of generative networks, which can handle various conditional generation tasks across different domains. Compared to the naive approach of training a separate standalone model for each task, our approach is parameter efficient since it learns to perform new tasks by making use of the model trained on previous tasks and combining with a small portion of unconstrained filters. At the same time, our model is able to maintain image generation quality comparable to the single standalone model for each task. Since the filters learned for previous tasks are preserved, our model is capable of preserving the exact generation quality of previous tasks. We validated our approach on various image-conditioned generation tasks across different domains, and the qualitative and quantitative results show our model addresses catastrophic forgetting effectively and efficiently.

\section{Related Work}

Our work intersects two previously disjoint lines of research: lifelong learning of generative models, and parameter-efficient network learning.

\textbf{Lifelong learning of generative models.}
For discriminative tasks, e.g.\ classification, many works have been proposed recently for solving the problem of catastrophic forgetting. Knowledge distillation based approaches \cite{shmelkov2017incremental,li2017learning,castro2018end} work by minimizing the discrepancy between the output of the old and new network. Regularization-based approaches \cite{kirkpatrick2017overcoming,chaudhry2018riemannian,aljundi2018memory,chaudhry2018efficient,lopez2017gradient} regularize the network parameters when learning new tasks. Task-based approaches \cite{Serr2018OvercomingCF,aljundi2017expert,rusuetal2016} adopt task-specific modules to learn each task.

For generative tasks, lifelong learning is an underexplored area and relatively less work studies the problem of catastrophic forgetting. Continual generative modeling was first introduced by Seff et al.~ \cite{seff2017continual}, which incorporated the idea of \textit{Elastic Weight Consolidation (EWC)}~\cite{kirkpatrick2017overcoming} into the lifelong learning for GANs. The idea of \textit{memory replay}, in which images generated from a model trained on previous tasks are combined with the training images for the current task to allow for training of multiple tasks, is well explored by Wu et al.~ \cite{wu2018memory} for label-conditioned image generation. However, this approach is limited to label-conditioned image generation and is not applicable
for image-conditioned image generation since no ground-truth conditional inputs of previous tasks are provided. \textit{EWC}~\cite{kirkpatrick2017overcoming} has been adapted from classification tasks to generative tasks of label-conditioned image generation~\cite{seff2017continual,wu2018memory}, but they present limited capability in both remembering previous tasks and generating realistic images. \textit{Lifelong GAN}~\cite{zhai2019lifelong} is a generic knowledge-distillation based approach for conditioned image generation in lifelong learning. However, the image quality of generated images of previous tasks keeps decreasing when new tasks are learned. Moreover, all approaches mentioned above fail in the scenario when there are conflicts in the input-output space across tasks. To the best of our knowledge, Piggyback GAN is the first method that can preserve the generation quality of previous tasks, and enables various conditional generation tasks across different domains.

\textbf{Parameter-efficient network learning.}
In recent years, there has been increasing interest in better aligning the resource consumption of neural networks with the computational constraints of the platforms on which they will be deployed \cite{chenetal2018,wuetal2019,yangetal2018}. One important avenue of active research is parameter efficiency. Trained neural networks are typically highly over-parameterized \cite{deniletal2013}. There are many effective strategies for improving the parameter efficiency of trained networks: for example, pruning algorithms learn to remove redundant or less important parameters \cite{heetal2019,yuetal2018}, and quantization algorithms learn to represent parameters using fewer bits while maintaining task performance \cite{jacobetal2018,jungetal2019,zhangetal2018}. Parameter efficiency can also be encouraged during training via architecture design \cite{maetal2018,sandleretal2018}, sparsity inducing priors \cite{tartaglioneetal2018}, resource re-allocation \cite{qiaoetal2019}, or knowledge distillation from a teacher network \cite{tungmori2019,zagoruykokomodakis2017}.
Our work is orthogonal to these methods, which may be applied on top of Piggyback GAN.

The focus of this work is parameter-efficient lifelong learning: we aim to leverage weight reuse and adaptation to improve parameter efficiency when extending a trained network to new tasks, while avoiding the catastrophic forgetting of previous tasks. Several approaches have been explored for parameter-efficient lifelong learning in a discriminative setting. For example, progressive networks \cite{rusuetal2016} train a new network for each task and transfer knowledge from previous tasks by learning lateral connections, a linear combination of layer activations from previous tasks is computed and added to the inputs to the corresponding layer of the new task. While our approach computes a linear combination of filters to construct new filters and uses them the same as normal filters. Rebuffi et al.~ \cite{rebuffietal2018} aims to train a universal vector that are shared among all domains providing generic filters for all domains, and task-specific adapter modules are added to adjust the trained network to new tasks. While we are working on solving a sequence of tasks, the filters trained for the initial task are hardly generic enough for later tasks. Inspired by binary neural networks, Mallya et al.~ \cite{mallyaetal2018} introduced the piggyback masking method (which inspired the naming of Piggyback GAN). This method learns a binary mask over the network parameters for each new task. Tasks share the same backbone parameters but differ in the parameters that are enabled. Multi-task attention networks
\cite{liuetal2019} learn soft attention modules for each new task. Task-specific features are obtained by scaling features with an input-dependent soft attention mask.
To the best of our knowledge, Piggyback GAN is the first method for parameter-efficient lifelong learning in a generative setting.

\section{Experiments}
We evaluate Piggyback GAN under two settings: (1) paired image-conditioned generation, and (2) unpaired image-conditioned generation. We first conduct an ablation study on the piggyback filters and unconstrained filters. We also demonstrate the generalization ability of our model by having same $T_2$
following different $T_1$. Finally, we compare our model with state-of-the-art approach~\cite{mallyaetal2018} proposed for discriminative models (i.e., classification models), which also shares the idea of ``piggybacking" on previously trained models by reusing the filters. Different from our approach, \cite{mallyaetal2018} learns a binary mask applying on the filters of a base network for each new task.

\noindent \textbf{Training Details.} All the sequential generation models are trained on images of size $256 \times 256$. We use the Tensorflow~\cite{abadi2016tensorflow} framework with Adam optimizer~\cite{kingma2014adam}. We set the parameters $\lambda=\frac{1}{4}$ for all experiments. For paired image-conditioned generation, we use UNet architecture~\cite{isola2016pix2pix} and the last layer is set to be task-specific (a task-specific layer contains only unconstrained filters) for our approach and all baselines. For unpaired image-conditioned generation, we adopt the architecture~\cite{johnson2016perceptual,zhu2017unpaired} which have shown impressive results for neural style transfer, and last two layers are set to be task-specific for our approach and all baselines. For both conditional generative scenarios, bias terms are used to adjust the output after factorization.

\noindent \textbf{Baseline Models.} We compare Piggyback GAN to the following baseline models: (a) \textit{Full}: The model is trained on single task, which could be treated as the ``upper bound" for all approaches. (b) \textit{$\frac{1}{4}$Full} and \textit{$\frac{1}{2}$Full}: The model is trained on single task, and $\frac{1}{4}$ or $\frac{1}{2}$ number of filters used in baseline \textit{Full} is used. (c) \textit{Pure Factorization (PF)}: Model is trained with piggyback filters that are purely constructed from previously trained filters. (d) \textit{Sequential Fine-tuning (SFT)}: The model is fine-tuned in a sequential manner, with parameters initialized from the model trained/fine-tuned on the previous task. (e)\textit{~\cite{mallyaetal2018}}: a state-of-the-art approach for discriminative model (i.e. classification model) which reuses the filters by learning and applying a binary mask on the filters of a base network for each new task.

\noindent \textbf{Quantitative Metrics.} In this work, we use two metrics \textit{Acc} and \textit{Fréchet Inception Distance (FID)~\cite{heusel2017gans}} to validate the quality of the generated data. \textit{Acc} is the accuracy of the classifier network trained on real images and evaluated on generated images (higher Acc indicates better generation results). \textit{FID} is an extensively used metric to compare the statistics of generated samples to samples from a real dataset. We use this metric to quantitatively evaluate the quality of generated images (lower FID indicates higher generation quality).

\subsection{Paired Image-conditioned Generation}
We first demonstrate the effectiveness of Piggyback GAN on 4 tasks of paired image-conditioned generation, which are all image-to-image translations, on challenging domains and datasets with large variation across different modalities~\cite{cordts2016cityscapes,isola2016pix2pix,tylevcek2013spatial,yu2014fine}.
The first task is \textit{semantic labels} $\rightarrow$ \textit{street photos}, the second task is  \textit{maps} $\rightarrow$ \textit{aerial photos}, the third task is \textit{segmentations} $\rightarrow$ \textit{facades}, and the fourth task is \textit{edges} $\rightarrow$ \textit{handbag photos}.

\noindent \textbf{Ablation study on choice of $\lambda$.} First we conduct an ablation study on the choice of different values of $\lambda$. For each upcoming task, we explore a set of 5 values of $\lambda$: $0$, $\frac{1}{8}$, $\frac{1}{4}$, $\frac{1}{2}$ and $1$. $\lambda=0$ corresponds to the special case where all filters are piggyback filters, and $\lambda=1$ corresponds to the special case where all filters are unconstrained filters. Figure~\ref{fig:lambda} illustrates the FID for different $\lambda$. For example, $\lambda=\frac{1}{4}$ balances well the performance and number of parameters. While parameter size scales almost linearly with $\lambda$, the improvement in performance (lower FID is better) slows down gradually.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.65\textwidth]{figures/experiments/lambda_2.png}
    \caption{FID for different $\lambda$. Ablation study of choice of different values of $\lambda$. For all upcoming new tasks, the performance improves as $\lambda$ increases. However, while parameter size scales almost linearly with $\lambda$, the improvement of performance slows down gradually.}
    \label{fig:lambda}
\end{figure}

\noindent \textbf{Ablation study on model components.} We also conduct an ablation study on the piggyback filters and unconstrained filters by comparing Piggyback GAN with the baseline models \textit{Full}, \textit{$\frac{1}{4}$Full} and \textit{Pure Factorization (PF)} on Task 2 (\textit{maps} $\rightarrow$ \textit{aerial photos}) given the same model trained on Task 1 (\textit{semantic labels} $\rightarrow$ \textit{street photos}).

The quantitative evaluations of all approaches are summarized in Table~\ref{table:ablation-paired}. The baseline \textit{$\frac{1}{4}$Full} produces the worst classification accuracy. Since the details like building blocks are hardly visible, the classifier sometimes mistakes category \textit{facades} for category \textit{maps}.  The baseline \textit{PF} generates images that contain more details as compared to the images generated using \textit{$\frac{1}{4}$Full}. This suggests that \textit{piggyback filters} can provide valuable information of the patterns in the training data. High \textit{FID} scores for baselines \textit{$\frac{1}{4}$Full} and \textit{PF} indicate that generation qualities for both approaches are poor, and it is observed that the generated images are blurry, resulting in lots of missing details, and also contain lots of artifacts for both approaches. Our Piggyback GAN is parameter efficient and produces images having similar quality as the \textit{Full} model.

\begin{table}[h!]
\centering
\begin{tabular}{p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{2.5cm}<{\centering}}
\toprule
 & Full & $\frac{1}{2}$Full  & $\frac{1}{4}$Full & PF & Piggyback GAN \\
\cmidrule[0.04em](lr{.2em}){1-6}
Acc & 97.90 & 88.10 & 84.87 & 99.82 & 97.99 \\
FID & 156.23 & 189.08 & 285.71 & 303.90 & 171.04 \\
\bottomrule
\end{tabular}
\caption{Model components. Ablation study of model components on paired image-conditioned generation tasks. Different models are trained and evaluated on task 2, FID and classification score on task 2 is reported.}
\label{table:ablation-paired}
\end{table}

\noindent \textbf{Generalization ability.} To demonstrate the generalization ability of our model, we learn the same $T_2$: \textit{maps} $\rightarrow$ \textit{aerial photos} with three different initial tasks $T_1$: \textit{semantic labels} $\rightarrow$ \textit{street photos}, \textit{segmentations} $\rightarrow$ \textit{facades}, and \textit{edges} $\rightarrow$ \textit{handbag photos} to evaluate whether Piggyback GAN could generalize well. The quantitative evaluation of all approaches are summarized in Table~\ref{table:generality}. The experimental results indicate that Piggyback GAN performs stable. The consistent Acc and FID scores indicate that the generated images on task $T_2$ from all three initial tasks $T_1$ have similar high image quality.

\begin{table}[h!]
\centering
\begin{tabular}{p{1.2cm}<{\centering} p{2.5cm}<{\centering} p{2.2cm}<{\centering} p{1.5cm}<{\centering}}
\toprule
$T_1$ & \makecell{\textit{semantic labels} \\ $\downarrow$ \\ \textit{street photos}} & \makecell{\textit{segmentations} \\ $\downarrow$ \\ \textit{facades}} & \makecell{\textit{edges} \\ $\downarrow$ \\ \textit{handbags}} \\
\cmidrule[0.04em](lr{.2em}){1-4}
Acc & 97.99 & 97.08 & 98.09  \\
FID & 171.04 & 174.67 & 166.20  \\
\bottomrule
\end{tabular}
\caption{Generalizability. Given different models trained for different initial tasks $T_1$, same new task $T_2$ is learnt. The consistent Acc and FID scores indicate that our model generalize well and perform consistently well given different set of filter banks achieved from different tasks.}
\label{table:generality}
\end{table}

\noindent \textbf{Comparison with SOTA method and baselines.} We compare Piggyback GAN with the state-of-the-art approach~\cite{mallyaetal2018}, \textit{Sequential Fine-tuning (SFT)} and baseline \textit{Full}. Baseline \textit{Full} can be considered as the ``upper bound" approach, which provides the best performance that Piggyback GAN and \cite{mallyaetal2018} could achieve. For both Piggyback GAN and sequential fine-tuning, the model of \textit{Task2} is initialized from the same model trained on \textit{Task1}. The same model also serves as the backbone network for approach~\cite{mallyaetal2018}.

Generated images of each task for all approaches are shown in Figure~\ref{fig:paired} and the quantitative evaluations of all approaches are summarized in Table~\ref{table:comparisons-paired}. It is clear that the sequentially fine-tuned model completely forgets the previous task and can only generate incoherent edges2handbags-like (\textit{edges} $\rightarrow$ \textit{handbag photos})-like patterns. The classification accuracy suggests that state-of-the-art approach~\cite{mallyaetal2018} is able to adapt to the new task. However, it fails to capture the fine details of the new task, thus resulting in large artifacts in the generated images, especially for task 2 (\textit{maps} $\rightarrow$ \textit{aerial photos}) and task 4 (\textit{edges} $\rightarrow$ \textit{handbag photos}). In contrast, Piggyback GAN learns the current generative task while remembering the previous task and can better preserve the quality of generated images while shrinking the number of parameters by reusing the filters from previous tasks.

\begin{table}[h!]
\centering
\begin{tabular}{p{1.5cm}<{\centering} p{2.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering}}
\toprule
 & Piggyback GAN & \cite{mallyaetal2018} & SFT & Full \\
\cmidrule[0.04em](lr{.2em}){1-5}
Acc & 89.80 & 88.09 & 24.92 & 91.10\\
FID & 137.87 & 178.16 & 259.76 & 130.71 \\
\bottomrule
\end{tabular}
\caption{Quantitative evaluation among different approaches for continual learning of paired image-conditioned generation tasks. Different models are trained and evaluated on tasks 1-4 based on the same model trained on task 1. The average FID score over 4 tasks is reported.}
\label{table:comparisons-paired}
\end{table}

\begin{figure*}[h!]
\begin{center}
  \includegraphics[width=\textwidth]{figures/experiments/comp-paired.png}
\end{center}
\caption{Comparison among different approaches for continual learning of paired image-conditioned generation tasks. Piggyback GAN can learn the current task without forgetting the previous ones, and preserve the generation quality while shrinking the number of parameters by ``piggybacking" on an existing model.}
\label{fig:paired}
\end{figure*}

For paired image generation, we also compare our approach with \textit{Lifelong GAN}~\cite{zhai2019lifelong} and \textit{Progressive Network (PN)}~\cite{rusuetal2016}. The FID of \textit{PN} over all tasks is $141.29$. \textit{Lifelong GAN} suffers from quality degradation: the FID on cityscapes increases from 122.52 to 160.28 then to 221.07 as tasks are added, while Piggyback GAN ensures no quality degradation since the filters for previous tasks are not altered. These results show the effectiveness and the flexibility of combining the two types of filters of Piggyback GAN.

\subsection{Unpaired Image-conditioned Generation}

We also apply Piggyback GAN to another challenging scenario: unpaired image-conditioned generation $\mathrm{domain} \, A \rightarrow \mathrm{domain} \, B$. The model is trained on unpaired data, where correspondence between domain $A$ and domain $B$ does not exist, namely for each image in domain $A$ there is no corresponding ground-truth image in domain $B$. We apply our model on 2 sequences of tasks: tasks of image-to-image translation and tasks of style transfer.

\noindent \textbf{Tasks of image-to-image translation.} We convert the same 4 tasks from the paired scenario to the unpaired scenario, and compared our approach with \textit{Full}, the state-of-the-art approach~\cite{mallyaetal2018} and \textit{Sequential Fine-tuning (SFT)}. The results are shown in Table~\ref{table:comparisons-unpaired}. We observed that the Sequential Fine-tuning (SFT) cannot remember previous tasks and suffers catastrophic forgetting. Our approach produces images with high quality on par with Full model, while [24] is capable of learning each new task, however the generation quality is poor.

\noindent \textbf{Tasks of style transfer.} Two tasks in a given sequence may share the same input domain but have different output domains, e.g. $T_1$ is \textit{Photo $\rightarrow$ Monet Paintings} and $T_2$ is \textit{Photo $\rightarrow$ Ukiyo-e Paintings}. While this is not a problem when training each task separately, it does cause problems for lifelong learning. We explore a sequence of tasks of unpaired image-conditioned generation as described above. The first task is unpaired image-to-image translation of \textit{photos} $\rightarrow$ \textit{Monet paintings}, the second task is unpaired image-to-image translation of \textit{photos} $\rightarrow$ \textit{Ukiyo-e paintings}.

We compare Piggyback GAN against the state-of-the-art approach~\cite{mallyaetal2018}, \textit{Sequential Fine-tuning (SFT)} and baseline \textit{Full} which serves as the ``upper bound" approach. For both Piggyback GAN and sequential fine-tuning, the model of Task 2 is initialized from the same model trained on Task 1, which also serves as the backbone network for approach~\cite{mallyaetal2018} to allow for fair comparison.

The quantitative and qualitative evaluations for comparison among different approaches are shown in Table~\ref{table:comparisons-unpaired} and Figure~\ref{fig:unpaired}, respectively. The baseline \textit{SFT} completely forgets the previous \textit{Monet} style and can only produce images in the \textit{Ukiyo-e} style. Both the state-of-the-art approach~\cite{mallyaetal2018} and Piggyback GAN are able to adapt to the new style.  The classification accuracy and FID scores indicate that Piggyback GAN can better preserve the quality of generated images, and produce images with styles most closely resembling \textit{Ukiyo-e} paintings.

\begin{table}[h!]
\centering
\begin{tabular}{p{1.5cm}<{\centering}| p{2.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering} p{1.5cm}<{\centering}}
\toprule
 Tasks & \multicolumn{4}{c}{4 Image-to-image Translation Tasks}\\
\cmidrule[0.04em](lr{.2em}){1-5}
 Methods & Piggyback GAN & \cite{mallyaetal2018} & SFT & Full \\
\cmidrule[0.04em](lr{.2em}){1-5}
Acc	& 90.30	& 88.29	& 24.85	& 91.04\\
FID	& 113.89 & 149.51& 262.37 & 109.46 \\
\bottomrule
\multicolumn{5}{c}{}\\
\toprule
 Tasks & \multicolumn{4}{c}{2 Style Transfer Tasks}\\
\cmidrule[0.04em](lr{.2em}){1-5}
 Methods & Piggyback GAN & \cite{mallyaetal2018} & SFT & Full \\
\cmidrule[0.04em](lr{.2em}){1-5}
Acc & 77.97 & 70.58 & 50.00 & 78.97 \\
FID & 106.95 & 118.26 & 135.05 & 101.09 \\
\bottomrule
\end{tabular}
\caption{Unpaired image conditioned generation tasks. Different models are trained on all tasks based on the same model trained on task 1. The average score over all tasks is reported.}
\label{table:comparisons-unpaired}
\end{table}

\begin{figure*}[h!]
\begin{center}
  \includegraphics[width=\textwidth]{figures/experiments/comp-unpaired.png}
\end{center}
\caption{Comparison among different approaches for continual learning of unpaired style transfer tasks. Piggyback GAN can preserve the generation quality the current task most without forgetting the previous ones.}
\label{fig:unpaired}
\end{figure*}


\section{Introduction}

Due to the fact that it is unpredictable, violent, and often hidden in clouds, the common phenomenon of lightning still harbors many secrets such as the conditions for initiation and the processes that allow the formation of conducting ionized channels called leaders, just to name a few. During any part of its development a lightning flash emits copious amounts of electromagnetic radiation over a very broad frequency spectrum, ranging from static electric fields to gamma-ray emission. For this work we use the emission in the Very-High Frequency (VHF) band to image the three dimensional (3D) lightning development in time. The major advantage of VHF-imaging over, for example, optical imaging, is that VHF is not obscured by the presence of clouds allowing for large distance observations.

The imaging of lightning in three dimensions using VHF-interferometric-based methods has a history dating back to the early early work of \citep{Mardiana:2002} where 3D is reached by combining angular measurements from two different stations. This work was followed up by others, such as \cite{Akita:2014,Stock:2014,LiuHY:2018} and culminated recently in the work of \cite{Jensen:2021} reaching a resolution of 200~m using 6 antennas.

We have developed similar techniques for imaging lightning with LOFAR \cite{Haarlem:2013}, which has resulted in much improved resolution. In this previous work we had used the signals of several hundred antennas for 3D imaging \cite{Hare:2019}. In this impulsive-imaging approach the time of the peak in the cross correlation between a pulse on the reference station and pulses in other stations is used in a chi-square fitting procedure to find the most likely location for the source. This impulsive imager was optimized a step further in \cite{Scholten:2021-init} and can find over 200 source locations per millisecond where meter-scale accuracy is reached.

Despite the history in the lightning field of naming these previous techniques as "interferometry", extremely few of them \cite{Stock:2014-PhD, Tilles:2019} use actually true interferometry in the sense of coherently summing the electric fields measured by many antennas.
In this work we show the results of using a true interferometric beamforming technique for imaging lightning in 3D using hundreds of LOFAR antennas with maximal baselines of 100~km.
In beamforming the signals from hundreds antennas are coherently added and the imaged volume is rastered into voxels, volumetric pixels.
Where most interferometric approaches integrate over considerable time periods, we divide the beamformed time trace for each voxel in narrow slices. A slice may be as small as 100~ns, the impulse-response time of the LOFAR antennas. For each time slice and each voxel the coherent intensity is calculated, thus forming a 4D (three space and one time dimension) intensity image. For each time slice the position of maximal coherent intensity is taken as the position of the source. In spite of the high time resolution we still resolve sources that emit with intensity well below the noise limit, where the noise level of a single antenna is dominated by the galactic background \cite{Mulrey:2019}. Since we can reach a time resolution of the order of the impulse-response time of our system we name this procedure Time Resolved Interferometric 3-Dimensional (TRID) imaging. In \secref{Method} the TRID imaging technique is discussed in detail.

Lightning, being an electric discharge phenomenon, has positively and negatively charged propagating structures, called leaders. While the positive leaders propagate in a rather continuous way, the negative leaders, in contrast, propagate discontinuously through a stepping process. In this stepping process copious amounts of VHF are emitted making them ideal test cases for our imaging method. As a first application we thus have used TRID imaging to investigate negative stepped leaders in \secref{NegLead}. The main draw-back of interferometric 3D imaging is that it is compute-intensive, thus without large dedicated compute resources only a relatively small region of space can be imaged. In \secref{LargeScale} we explain our two-step approach where the entire lightning flash is imaged with the impulsive imager \cite{Scholten:2021-init} and from this image the interesting parts are selected that are subsequently imaged using the TRID imager.
We show in \secref{Dynamics} that the TRID imager allows to resolve the dynamics of a corona flash for some negative leaders. As a second application we show in \secref{PowerLaw} that the intensity spectrum of VHF-pulses emitted from negative leaders follows an almost perfect inverse-square power-law over four orders of magnitude (with a scaling exponent of $-1.8$).
We end with a short discussion in \secref{summ}.

\section{Methods}\seclab{Method}

The LOFAR infrastructure consists of several thousand antennas distributed in stations over much of Europe which a dense core in the Netherlands. For the present analysis we use the antennas installed in the Netherlands. The LOFAR antennas are arranged in stations each containing 96 Low Band Antennas (LBA, inverted v-shaped dipoles operating in the 30 -- 80~MHz band) and a similar number of High Band Antennas (100-200~MHz, not used here). For practical reasons based in the observation mode the data of 12 dipoles each from all 37 Dutch stations is used where the geographic spread is shown in \figref{LOFAR-NL} covering a logarithmic baseline distribution up to 100~km.  At the core of LOFAR, called the Superterp, the antenna stations are more closely packed. Upon a trigger the RAM buffers containing raw data from every antenna are frozen and read out over glass fiber to a central storage. This allows for off-line processing of the time traces (sampled at 200~MHz) of all antennas. The atomic clocks in the stations allow for a timing stability of a fraction of a nanosecond over the recording time span of maximally 5 seconds.

\begin{figure}[th]
\centering{\includegraphics[bb=3.5cm 3.5cm 25.5cm 25.cm,clip, width=.47\textwidth]{LOFARmapFigure} }

	\caption{Layout of the Dutch LOFAR stations. The core of LOFAR is indicated by the red $\bigoplus$\ sign while the yellow stars show the general location of flashes A and B that are  discussed in this work. The black frame indicates the general area ($140\times 140$~km$^2$) where flashes can be mapped accurately and is centered at the LOFAR core.  Stations marked in red are not included in TRID imaging (see text).}	 \figlab{LOFAR-NL}
\end{figure}

The basic implementation of beamforming imaging is relatively straightforward. The part of the atmosphere where the lightning is to be imaged is rastered into voxels, see \secref{Vox}. For each voxel the time traces of all antennas are summed while accounting for the difference in travel time to each antenna. This yields for each voxel a time trace which is the coherent sum of all antennas, as discussed in \secref{Vox}. Adding the traces for the different voxels is always performed for a fixed time span in one particular `reference' antenna, usually an antenna from the core of LOFAR. Thus we can find the source location of a particular structure seen in the trace of the reference antenna. The choice of reference antenna is not important as long as it is in the dense core.

The time trace of each voxel is cut in slices of a fixed duration. The slicing is done such that they are synchronous when correcting for the signal travel time from the reference antenna to the voxel. For all slices the coherent intensity is determined resulting in a voxelated intensity profile for each time slice, similar to what is shown in \figref{B-NLa-166}. For each time slice the maximum of this voxelated intensity profile is determined and used as the source location, as discussed in more detail in \secref{Max}.
Plotting the position of all sources will result in the beam-formed (or interferometric) image of a segment of the flash in space and time where the length of the time slices determines our final time resolution since any detail within a slice is summed over. In the following the different steps involved in this procedure are discussed in more detail.

\subsection{Space \& time grids and summing time traces}\seclab{Vox}

The part of the atmosphere that is of interest is rastered into voxels. We have found that this is most conveniently done in spherical coordinates with the center at the reference antenna in the dense LOFAR-core. In this way we can more easily account for the fact that our resolution in the two transverse directions is much better than in the radial direction. The optimal grid-spacing is dependent on the areal spread of the antenna locations. For the example discussed in this work we include antennas within 50~km distance from the core where the antennas are distributed on an irregular logarithmic grid with a dense core that is designed to be optimal for astrophysical applications, see \cite{Haarlem:2013} and \figref{LOFAR-NL}. The voxel grid is chosen such that the maximal time shift for any antenna for two bordering voxels is about 1~ns, our time calibration accuracy. This implies for the examples discussed in  \secref{NegLead} a grid spacing of 0.003$^\circ$ in the azimuth angle ($\hat\phi$), 0.01$^\circ$ in elevation angle ($\hat\theta_e$, upwards from the horizon), and 10~m in the radial direction ($\hat{R}$). For a source at 50~km distance this implies a grid spacing of about 1~m transverse to the line of sight (from the reference antenna) and 10~m along the line of sight. A typical grid used in  \secref{NegLead} spans over $(60 \times 30 \times 20)$ voxels in $(\hat\phi,\,\hat\theta_e,\,\hat{R})$.

A section of the time trace in the reference antenna is selected with typically a length of 0.3~ms. For each voxel the
time traces of all antennas are added after time shifting them to account for the travel time differences of a signal from the voxel. This yields the coherent (or beam formed) time trace for that voxel. All coherent voxel time traces are evaluated at identical times for the reference antenna.

The actual time shift of the time traces is done in Fourier space by applying a frequency-dependent phase shift to the trace for each antenna. After summation, the traces are transformed back to the time domain. The actual antenna time traces are taken longer than the required 0.3~ms to be able to account for the maximal relative time shift over the full voxelated volume.
We thus obtain for each voxel a coherent time trace, with all antennas contributing, of 0.3~ms length with additional trailing ends where there will be partial coherence. These trailing end are excluded from the analysis.

\begin{figure}[th]
\centering{\includegraphics[ width=0.49\textwidth]{InterfContourFig1} }
	\caption{A typical interferometer image for Flash B. The axes show the grid in voxel number where $(0,0,0)$ is the center voxel in the image cube.  Both images are made for a thin (one voxel thick) section through the center voxel. To show the general features of the intensity plot, the hyper cube is chosen such that the maximum intensity is at the center voxel, which generally is not the case. The typical length scale for $d\phi$ and $d\theta_e$ is about 1~m. The intensity unit [gb] is explained in \secref{Cal}. }	 \figlab{B-NLa-166}
\end{figure}

The intensity may be integrated over the full time trace of 0.3~ms which washes out all time dependence over the integration range.
However, a much better time resolution is generally required in order to image the dynamics of the flash.
To obtain this, we cut the time trace from each voxel into narrow slices and we integrate the power over each slice for every voxel to obtain the voxelated intensity profile for each slice, giving an image like \figref{B-NLa-166}. The optimal slice length depends on the research question of interest, taking it short resolves  finer time structures but makes the imaging more noisy, taking it long  improves the signal-to-noise ratio, and thus the sensitivity, but increases the chance that multiple sources appear simultaneously and are confused.
The smallest time step where subsequent time-frames are considered to be independent is the width of the impulse response function.
The impulse-response time is determined by the frequency dependence of the LOFAR antennas, amplifiers, and filters. An easy estimate of this is obtained from the full width at half maximum of the narrow peaks in the spectrum, about 50~ns.
As shown in \secref{NegLead} very stable results are obtained even for a time slice as short as 100~ns, which is a little longer than the impulse response time.
Since we have shown that we can reach down to close to the impulse-response time of the system we name this Time-Resolved Interferometric 3D (TRID) imaging.

\subsection{Antenna weighting}\seclab{AntWeight}

Since the array of antennas has an extent that is comparable to the distance to the source region, the received signal strength varies greatly over the array. If this is not taken into account in the analysis one would effectively use the antennas that are near the source region only and thus potentially use only a fraction of the imaging accuracy of the array as a whole. To improve this aspect of the performance we include weighting factors for the antennas to compensate for the imbalance in signal strength. In interferometry applications for astronomy (2D, not 3D imaging) there is a considerable amount of literature regarding the advantages and disadvantages of different weighting schemes, see \cite{Briggs:1999, Yatawatta:2014} for comprehensive reviews.

For this work we have chosen for a weighting scheme that compensates in lowest order the signal strength variations in the antennas due to their distance from the source, where the weight is capped for antennas far from the source for which the signal to noise ratio becomes unfavorable. In addition it takes into account the fact that the largest density of antennas is at the core.

The amplitude of an emitted signal drops inversely proportional to distance $R_{as}=\sqrt{D_{as}^2+h_s^2}$, where $D_{as}$ is the distance from a certain antenna to the source in the horizontal plane and $h_s$ denotes the altitude of the source. For this weighting scheme  it is assumed that the antennas are all in a horizontal plane since that part of the Netherlands is flat to a good approximation, in all other interferometry calculations the proper 3-dimensional locations of the antennas are used. Another important factor to account for is the antenna gain that depends on the azimuth $\phi$ and elevation $\theta_e$ angles of the source with respect to the antenna. The antenna gain vanishes for sources at the horizon and thus (since it is an analytic function of the angles) depends on $\theta_e$ as $\sin{(\theta_e)}$. For almost all cases of interest the sources are at small elevation angles for which we can assume $\sin{(\theta_e)}\approx h_s/R_{as}$ where the antenna position enters in $R_{as}$ only.

The antenna gain depends also on $\phi$ as well as on the polarization of the signal. Since an analysis of the polarization of the radiation falls outside the scope of this work as it would make the analysis considerably more complicated, we have ignored this dependence. Even though the antennas are not in a small azimuthal angle range from the source (full opening angle about 60$^\circ$, see \figref{LOFAR-NL}) we still achieve high-quality results.

Combining the $1/R_{as}$ drop of signal strength with the proportionality of the antenna gain to $\sin{(\theta_e)}$, we obtain an antenna weighting factor
\begin{equation}
W_a=\left\{
  \begin{array}{lr}
    R_{as}^2/R_{rs}^2 & \rm{if} \quad R_{as}^2/R_{rs}^2 < W_{\rm max}\\
     W_{\rm max} & \rm{if} \quad  R_{as}^2/R_{rs}^2 \ge  W_{\rm max}
  \end{array} \;. \eqlab{W}
\right.
\end{equation}
where the weight is normalized to unity for the reference antenna at a distance of $R_{rs}$. In addition the weight is capped to a maximum of  $W_{\rm max}=1.2$ for distant stations where the signal to noise ratio is getting worse. It is conceivable that a different choice for $ W_{\rm max}$, for example by making it dependent of the signal-to-noise ratio, will improve the image accuracy.

\subsection{Antenna calibration and Intensity}\seclab{Cal}

The antenna timings are calibrated for each flash following the procedure outlined in \cite{Hare:2019,Scholten:2021-init}. Per flash 20  -- 30 bright stand-alone pulses are selected from the whole flash where care is taken that their source locations roughly cover the extent of the flash. For all sources the location as well as the antenna timings are searched in a simultaneous fit using the search algorithm discussed in \cite{Scholten:2021-init}.

Since for interferometry the stability of the phase of the signal over all antennas is important, we perform an additional check of this phase for a single distinct pulse in the spectrum selected for TRID imaging. If this phase is off by more than 90$^\circ$ the antenna will not be used. Usually this eliminates less than 1\% of all antennas. This may also result in eliminating a whole station from the analysis for a particular flash. The latter could be due to the fact that we have ignored the azimuthal angle dependence of the antenna function.

The gain of each antenna is calibrated by normalizing the noise level to unity. The noise level is determined from the parts of the recorded trace for which there is no lightning activity detected within 0.3~ms.  Since this noise level is largely due to non-terrestrial sources, i.e.\ radiation produced by the Milky Way, we name this the Galactic Background [gb] level.

The coherent or interferometric intensity is calculated as the intensity of the signal from the source as received by the reference antenna and is expressed in units of [gb]). It is tacitly assumed that the azimuth-angle dependence of the antenna gains is limited and is effectively averaged-out. The elevation-angle dependence is to a large extent accounted for by the weighting factors (see \secref{AntWeight}) when calculating the coherent intensity.

All LOFAR antennas used in this work are inverted v-shape dipoles where X- (Y-) dipoles, corresponding to odd (even) numbered antennas, are oriented in the NE-SW (NW-SE) plane. Our analysis is performed separately for X- and Y-dipoles since they differ significantly in their sensitivity to polarized radiation. The antenna function (the Jones matrix, specifying for each dipole and polarization the gain depending on angle and frequency) has to be used \cite{Haarlem:2013} to convert the measured intensity to the absolute intensity of the source, which we have not done in this work. The absolute calibration of the LOFAR antennas is performed in \cite{Mulrey:2019}.
Extracting intensities (and even polarization) of a source is considerably simplified if it is assumed that the source is point-like emitting a delta-pulse (flat frequency spectrum over our frequency range) and linearly polarized, however such an analysis falls outside the scope of the present work.

\subsection{Determining the location of maximal intensity}\seclab{Max}

\figref{B-NLa-166} shows a typical intensity profile for a time slice of 10~$\mu$s. For this particular one the coordinates of the hypercube have been chosen such that the peak intensity is at the center. Generally this is of course not the case. The procedure used to determine the position of the maximum is discussed in this section. Besides the main maximum there are several secondary maxima, due to so-called side-beams, a well known feature in interferometric imaging. The positions and intensities of these side beams is determined by the antenna configuration. The LOFAR irregular lattice of antennas suppresses the intensities of these side beams to a large extent. Secondary maxima could also be due to background noise or due to multiple sources in the same time slice.
It is straightforward to determine the voxel with maximal intensity. To reach an inter-voxel accuracy we have implemented two interpolation procedures, quadratic and barycentric interpolation.
\begin{description}
\item[Quadratic] The intensity of the voxels around the voxel with the maximum intensity are fit with a paraboloid,
    \begin{equation}
    I_p(\vec{x})=I_0 + \sum_i\left[\frac{1}{2} A_i x_i^2 + B_i x_i\right] + \sum_{i<j} R_{i,j} x_i x_j \;,
    \end{equation}
    where $x_i$ is the $i^{th}$ grid coordinate, $i=1,2,3$. The coefficients $I_0$, $A_i$, and $R_{i,j}$ are fitted to the grid points bordering the maximum. The inter-voxel maximum $\vec{x_p}$ is taken at the point where the paraboloid reaches its maximum.
\item[Barycentric] The barycentric maximum is calculated as
    \begin{equation}
    \vec{x}_b= \sum_{\vec{x}} \left[(I(\vec{x})-I_{\rm Th})\, \vec{x} \right] /  \sum_{\vec{x}} \left[(I(\vec{x})-I_{\rm Th}) \right] \;,
    \end{equation}
    where the sum runs over all voxels with an intensity exceeding the threshold value $I_{\rm Th}$. This threshold is taken as $I_0/1.2$ (where $I_0$ is the maximum voxel intensity) or the largest value of a voxel on the outer surface of the grid, whichever is larger.
\end{description}
When there are many active sources in small area, the barycentric interpolation will yield some weighted average position, while the quadratic interpolation yields the position of the strongest source.
It is observed that the details of intensity surface are complicated with small ripples (order 10\% in intensity over distances of 20~m depending on the source location) likely due to side beams. These ripples make the quadratic interpolation unstable for a coarse grid. Using barycentric interpolation these ripples are efficiently averaged yielding a properly interpolated maximum. The examples shown in \secref{NegLead} use a fine grid and for this reason the quadratic interpolation is used.

\section{Negative Leaders}\seclab{NegLead}

As an example of the possibilities offered by the new TRID imager we apply it to some normal negative leaders. These are selected as their structure is assumed to be relatively simple and one thus obtains a good insight in the intrinsic accuracy of the TRID imager. We use two negative leaders from each of two flashes, Flash A occurring on April 24, 2019 at 21:30:56.221 UTC and Flash B on August 14, 2020 at 14:14:58.669 UTC. These times are taken as $t=0$ for each flash. We show a zoomed-in view on the fine structure of one negative leader of Flash B and statistics for all four negative leaders.

The location of the two flashes with respect to the LOFAR core is shown in \figref{LOFAR-NL}. Both flashes have been imaged using the impulsive imager, where Flash A was reported on in \cite{Scholten:2021-init,Scholten:2021-RNL} and Flash B in \cite{Scholten:2021-HANL}. For the impulsive imager all stations are used, including those marked in red in \figref{LOFAR-NL}, while for TRID imaging we have excluded these stations because they are at a distance larger than 50~km from the LOFAR core. This limit on distance is taken since for further distances the signal becomes weak, or the spread in angles to the source too large. For Flash B also station RS208, RS508, and RS307 were excluded since for these the phase-stability for selected pulses could not be shown, see \secref{Cal}.

The sections of the negative leaders which we will use for the discussion in this work are given in \tabref{NegLead}.
\begin{table}[ht]
\centering
\begin{tabular}{||c|c||r@{ -- }l|r@{ -- }l|r@{ -- }l|r@{ -- }l||}
\hline
Flash & Leader & \multicolumn{2}{c|}{time [ms]} & \multicolumn{2}{c|}{North [km]} & \multicolumn{2}{c|}{East [km]} & \multicolumn{2}{c||}{height [km]} \\
\hline \hline
A & NL-A1 & 170& 178 &    21 & 23  &  18&20 &   2.5&3.5 \\
A & NL-A2 & 177&182  &    21&23    &  21&22 &   3. & 3.5  \\
B & NL-B1 & 238& 245 & $-25$&$-25$ &  37&38 &  5.3&5.8 \\
B & NL-B2 & 540& 552 & $-34$&$-34$ &  36&37 &  4.9& 5.5 \\
\hline
\end{tabular}
\caption{The negative leader sections that are used for further analysis. \tablab{NegLead}}
\end{table}
In \secref{LargeScale} we show that while the image quality on a large scale may be somewhat better for the TRID imager as compared to the impulsive imager, this does not outweigh the additional effort in CPU power. As is shown in \secref{Dynamics} the true strength of the TRID imager lies in resolving fainter sources keeping good time resolution. Another strong point of the TRID imager is that it is complete in the sense that it will locate all sources with a strength exceeding background in the image-cube with only minor caveats. The impulsive imager is in this respect much less complete since the located sources  are not restricted by the small hypercube (which can also be an advantage) and furthermore need to be well separated from earlier and later sources. The completeness of the TRID-imaged sources allows to analyze the VHF-intensity spectrum, as discussed in \secref{PowerLaw}.

\subsection{Larger scale structures}\seclab{LargeScale}

\begin{figure*}[th]
\centering{\includegraphics[ bb=1.0cm 2.4cm 24.5cm 25.7cm,clip, width=0.49\textwidth]{Fig3A} }
\centering{\includegraphics[ bb=1.0cm 2.4cm 24.5cm 25.7cm,clip, width=0.49\textwidth]{Fig3B}}
	\caption{A comparison of the images for section `NL-B1' from a negative leader of Flash B (see \tabref{NegLead}) as obtained by the impulsive imager (left panels) and the TRID imager (right panels) using the X-dipoles from LOFAR. The boxes on the right indicate the zoom-in volume for \figref{Intf-NLaZM}, left. }	 \figlab{Intf-NLaH}
\end{figure*}

Interferometric imaging is rather expensive. For the present implementation it typically takes half a day on a single CPU to produce an TRID image of a volume of about 1~km$^3$ for a trace of 0.3~ms while it takes about 1 hour to image a section of 100~ms of a flash with the impulsive imager, however, since the interferometric imaging code can easily be paralleled, this problem may not be as constraining at it appears. Both estimates are for a typical LOFAR configuration with a few hundred antennas.   In \figref{Intf-NLaH} the results of the two imaging procedures are compared for a section of about 1~km of a negative leader from Flash B. This represents only a small section of the full flash that covers an area of the order $25 \times 10$~km$^2$ and is active for 2~s. For the TRID image (\figref{Intf-NLaH}, right) only strong sources ($I>30$~gb) are plotted resulting in 399 sources while the quality cuts on the impulsive image (left) leaves 159 sources. The two imaging methods are subject to different systematic and statistical uncertainties and have a different resolution. One should therefore focus on general structures and not on individual points and as is seen from \figref{Intf-NLaH}, there is no significant difference between the two. The sources in the TRID image show a bit less scatter than those from the impulsive imager but the differences are marginal.

We observe that the received coherent power integrated over 0.3~ms for a section of a  negative leader of Flash A, located at the NE of the array, is about twice as large for the Y-dipoles as it is for the X-Dipoles. For Flash B (at the SE of the core) the intensity ratio is reversed. This is due to the azimuth-angle dependence of the antenna gain. For sources in the NE horizontally polarized radiation is recorded almost exclusively in the Y-dipoles with a gain (in power) that is about a factor two larger than that for vertical polarization recorded exclusively by the X-dipoles in this configuration while for sources in the SE the ratio is reversed. From the observation of the intensity ratio for X- and Y-dipoles one thus can conclude that on average there is no preferred polarization direction for the VHF emission from a negative leader.

\subsection{Detailed view of a negative leader}\seclab{Dynamics}

\begin{figure*}[th]
\centering{\includegraphics[ bb=1.0cm 2.4cm 24.5cm 30.4cm,clip, width=0.49\textwidth]{SourcesNLaZ} }
\centering{\includegraphics[ bb=1.0cm 2.4cm 24.5cm 30.4cm,clip, width=0.49\textwidth]{SourcesCorFlash} }
	\caption{Left panels show a zoom in on a section of the negative leader shown in \tabref{NegLead}. The sources found by the impulsive imager (44 sources) and in color those from TRID with a $I>30$~gb (99 sources) are indicated by black open circles. The arrow indicates the radial direction and points to the core of LOFAR. The right shows an even further zoom-in on one particular burst as imaged using TRID with an intensity threshold of $I>0.3$~gb showing 153 sources. In light (darker) grey the interferometric more intense sources are shown that fall after (before) the time frame of the plot. The area of the dots is proportional to the intensity for the colored strong sources. \figlab{Intf-NLaZM}}
\end{figure*}

The strength of TRID-imaging is shown in \figref{Intf-NLaZM} where zoom-ins are made at two different levels on the leader track shown in \figref{Intf-NLaH}. The left side of \figref{Intf-NLaZM} shows the stronger sources ($I>30$~gb) only, while for the zoom-in on the right all sources with an intensity $I>0.3$~gb (i.e.\ with a strength below the noise level of a single antenna) are shown.

One aspect evident from \figref{Intf-NLaZM} is that our imaging accuracy is very directional dependent, as was remarked already in \secref{Vox}. Sources occurring in a narrow time window, i.e. have the same color, show a much larger spread in the ($\hat{R}$) direction (radially outward from the core at (0,0), indicated by the black arrow on the left plan view of \figref{Intf-NLaZM}) than in the transverse direction, about a factor 10 difference.

\figref{Intf-NLaZM} shows clearly that the real strength of TRID imaging lies in its ability to image the finer details of a flash.
The time-colored sources on the left of \figref{Intf-NLaZM} are obtained through TRID imaging while the sources marked by black circles have been found using the impulsive imager, using the same quality and interferometric-intensity cuts as in \figref{Intf-NLaH} for both. The top left panel clearly shows that also this leader has a distinct burst structure that is typical for negative leader propagation and was analyzed extensively in \cite{Hare:2020}.
The right panels of \figref{Intf-NLaZM} show a zoom in on the  burst indicated by the black boxes on the left side. The size of the colored circles reflect source intensity. The limit on the intensity has been decreased to 0.3~gb, well below the noise level in the reference antenna, which is the reason that many more sources show (a total of 254 for the zoomed-in section). As is to be expected, the positioning accuracy of the weaker sources is less and they thus show more scatter, in particular in the radial direction. The grey sources in the right panels mark the positions of intense TRID-imaged sources that fall outside the imaged time span of 241.65 -- 241.75~ms used for the colored sources. The earlier sources are darker grey (at time $t < 241.31$~ms) and the lighter ones occur later than the colored sources (at $t>241.89$~ms). For all grey sources the dot-size is constant.

In the plan view (\figref{Intf-NLaZM}, bottom right) the leader propagates from the lower left corner to the middle right, i.e.\ from the dark grey dots to the light grey ones. The height v.s.\ time panel, the top of the right side of \figref{Intf-NLaZM}, shows that the burst itself has a very distinct structure, starting with some intense sources, the corona flash itself, followed by weaker sources with a decreasing density. From the plan view, the bottom panel, it is seen that they follow the propagation direction of the leader, as expected for streamers ensuing from a corona flash.
Unfortunately the image quality is worse in the radial direction (almost perpendicular to the propagation direction of the negative leader) and it thus cannot be seen if the streamer activity occurs in a narrow forward cone or in a broad fanned-out structure. The propagating ionization front has a width (in the propagation direction) of less than 5~m which is about the locating accuracy we reach for these small intensity sources. The streamers from the flash took 100~$\mu$s to propagate over a distance of about 20~m. The average velocity is about $2\times 10^5$~m/s, but it can be seen that the velocity of the streamers has dropped considerably towards the end of their reach. At this distance the following corona flash starts, as indicated by the light grey dots, some 160~$\mu$s after this flash ceases to show VHF-emission.

The structure shown at the right of  \figref{Intf-NLaZM} is general for all bursts we see on the two negative leader sections we have investigated for Flash B.
All burst start with some strong sources and many weaker sources where the number decreases in time. After 100~$\mu$s this has decreased to about the background level and the next burst occurs not more than 0.5~ms later at the furthest distance where streamer activity was observed from the previous flash. For the time period between the end of one burst and the start of the next the area is VHF-quiet and we detect only background fluctuations.

The observed structure is very reminiscent of what is observes for High Altitude Negative Leaders (HANL) in \cite{Scholten:2021-HANL} using the impulsive imager. The main difference is that a HANL step length covers a few 100~m while it is only a few 10~m in the present case. In the case of the HANL a filament structure could be distinguished. At the end of the reach of this filaments new coronal bursts were seen to occur, a considerable time after the filament activity had died away.

A possible mechanism is  that in the corona flash some poorly conductive channels (streamers) are formed that stop propagating when a maximal range is reached. Charge, however, continues to move along these streamers to accumulate at the end point. After a certain time, probably depending on conductivity and steamer length, the accumulated charge is sufficient for a break down creating the following corona flash. This would involve a process that may be similar to the stepping seen in laboratory when voltage changes linearly \cite{Kochkin:2014}.
In the current standard picture for the stepping process a space stem grows towards the head of the negative leader where upon contact charge flows towards the outer end of the space stem, heating the space stem to become a leader and initiating a following corona flash. In such a picture one would expect some VHF emission from the space stem when the sudden increase in (negative) charge causes a sudden drop in its potential. In this initial work we see absolutely no evidence for this process, which needs to be investigated further in future work.

The spatial development of corona flashes could not be distinguished as well on the negative leaders for Flash A, although it is seen for some.  This might be related to the fact that the investigated negative leaders for Flash A are at a lower altitude than for Flash B. Combined with the observations in \cite{Scholten:2021-HANL} this would support a strong altitude dependence of the stepping process, as was already suggested in \cite{WuT:2015, Lyu:2016}. We will investigate the occurrence of corona flashes in a future, more extensive study of negative leaders using  TRID imaging.

\subsection{Intensity spectrum}\seclab{PowerLaw}

\begin{figure}[th]
\centering{\includegraphics[ bb=0.2cm 2.0cm 22cm 15cm,clip, width=0.49\textwidth]{IntensityDistr} }
	\caption{Distribution of source intensities for the four negative leaders sections given in \tabref{NegLead}. The analysis shown is using the Y-dipoles only. A power law dependence is shown to guide the eye.}	 \figlab{IntfIntS}
\end{figure}

The impulsive imager can reconstruct the source location of a pulse only when it is clearly distinct in the time trace. In practice this means that only for the strongest pulses the source locations can be reconstructed and the reconstruction efficiency drops rapidly with pulse density. This imaging efficiency is difficult to quantify. Since the TRID imager does not rely on pulse finding, the spectrum slicing is independent of the pulse structure, the efficiency for reconstructing the source locations of the received intensity is almost 100\% for sources in the image cube, with the provision  that at most a single source position can be reconstructed within the time frame of a slice.
In practice this means that (almost) all sources are reconstructed when the density of sources in time is less than the inverse of the slicing time. Another provision is that there should  be no strong sources in the vicinity of the image cube since otherwise the images are blurred by the side beams from these strong sources.
This completeness of the spectrum allows for an analysis of the intensity spectrum.

The distribution of the interferometric source intensities (see \secref{Cal}) is given in \figref{IntfIntS} using a slicing time of 100~ns. The distributions are calculated for the four negative leader sections listed in \tabref{NegLead}. The interferometry analysis is performed separately for the X- and Y-dipoles that have different orientations and thus different polarization sensitivities. The extracted distributions are very similar for the two dipoles and only the results for the Y-dipole analysis is given in \figref{IntfIntS}.
The plotted distributions have been normalized (by dividing by a factor $F=N_t \, I_{\rm Th}$ where $N_t$ is the total number of sources in the distribution and $I_{\rm Th}$ is the threshold for including sources) to express the similarity for the different leaders more clearly. Additionally the plotted distributions have been scaled by the square of the intensity in order to emphasize smaller details.
The factor $F$ is chosen such that a perfect inverse square power law dependence will show as a horizontal line at unity.

The intensity distributions for the four negative leader sections are all very similar as is shown in \figref{IntfIntS}. The distributions follow a power law, and to guide the eye an $I^{-1.8}$ dependence is shown by the black line.
At smaller intensities the distributions saturate and even drop off to zero because at most a single source, the strongest, is located per time slice of 0.1~$\mu$s. The weaker sources are therefore `masked' by the stronger ones and the power spectrum is no longer complete. The main differences between the distributions show in the saturation at small intensities and is related to the measured power. The measured power depends on the altitude of the source as well as the distance, see \secref{AntWeight} for a discussion, and differs for X- and Y-dipoles.
At the largest intensities there are signs of a drop-off in the power law, however this may also be due to limited statistics.

\begin{table}[ht]
\centering
\begin{tabular}{||c||c|c|c||}
\hline
Leader & ${\cal N}$ & $\alpha$ & $\gamma$ \\
\hline
NL-A1-Y&   0.88 &  1.76 &    1.49 \\%th=2
NL-A1-X &  0.89 &  1.79 &    1.04 \\
\hline
NL-A2-Y &  0.85 &  1.75 &    1.24 \\%th=2.
NL-A2-X &  0.87 &  1.78 &    1.09 \\
\hline
NL-B1-Y &  1.3 &  1.85 &     0.110 \\
NL-B1-X &  1.0 &  1.81 &     0.135 \\
\hline
NL-B2-Y &  1.3 &  1.88 &     0.102 \\
NL-B2-X &  1.5 &  1.90 &    0.220 \\
\hline
\end{tabular}
\caption{The values for the normalization ${\cal N}$, the power $\alpha$, and the small-intensity suppression coefficient $\gamma$ in \eqref{PowLaw} as extracted for fits to the intensity distributions for the for negative leader sections listed in \tabref{NegLead}. Separate analyses are made for the X- and Y-dipoles. \tablab{IntAmp}}
\end{table}

To account for the suppression at small intensities we fit the normalized distributions $N(I)$ with a modified power law,
\eqref{PowLaw}.
\begin{equation}
N(I)= {\cal N} \,I^{-\alpha}  \,e^{-\gamma/I}\;, \eqlab{PowLaw}
\end{equation}
where $I$ is expressed in units of [gb]. The last factor, dependent on $\gamma$, suppresses the distribution at small amplitudes in good agreement with the data. The values for the fitted normalization ${\cal N}$, the power $\alpha$, and the small-intensity suppression factor $\gamma$ are given in \tabref{IntAmp} for four negative leader sections given in \tabref{NegLead} and the two different dipoles.

The fitted powers $\alpha$ are almost the same for all eight cases. The main differences are seen in the values for $\gamma$ indicating that for Flash A the Y-dipoles have a larger gain than the X-dipoles while for Flash B the situation is reversed, see the discussion in \secref{LargeScale}. Since Flash A occurs closer to the core, the $\gamma$ coefficients for Flash A are all larger than those for B. In a future work we will investigate the extent of the power law dependence in more detail, in particular if this applies to all parts of the flash and the cut-off in the power law which one would expect at the highest intensities since the total radiated energy should not diverge.

It is interesting to note that power-law dependencies are usually related to disruptive processes, \cite{MengFZ:2019,Stumpf:2012,Virkar:2014}, and in this respect it should have been no surprise to find this in association with a corona flash. On the other hand it is surprising to observe a power-law, which is scale invariant, as one is dealing with a process where many scales enter such as the length scale that is given by the stepping length, the time scale set by the stepping time, or the scale set by the total charge in the leader tip.

A power-law dependence was observed already in \cite{Machado:2021} using the data from the impulsive imager where an unbiased intensity spectrum could be obtained for the 10\% most intense sources.

\section{Summary and outlook}\seclab{summ}

Time Resolved Interferometric 3D-imaging gives an unprecedented insight in the finer dynamics of lightning propagation. Using hundreds of antennas with baselines of up to 100~km allows to image the more weakly emitting structures with (a better than) 10~m scale precision. The TRID imaging procedure should be considered complementary to the impulsive imager presented in \cite{Scholten:2021-init} that builds an image by locating the sources of single, reasonably well isolated, pulses in the spectrum. The impulsive imager can efficiently image the complete flash while the TRID-imager allows to image much fainter aspects of the flash as well as a somewhat improved locating accuracy, however  at the expense of being much more CPU-intensive.

As a first, exploratory, application of TRID imaging we have applied it to observe the finer details of negative leader bursts and imaged the temporal and spatial development of a corona flash. We show that it expands up to the  distance where the flash of the following step starts. In addition we show that the intensities of VHF pulses emitted by negative leaders follow an almost perfect inverse-square power law extending to the largest intensities. Each of these subjects requires a much more in-depth analysis which we plan for forthcoming publications.

In spite of the unprecedented possibilities of the present method to study the fine dynamics of lightning, there are still several aspects where the method should be improved. A minor improvement would be to optimize the antenna-weighting scheme discussed in \secref{AntWeight} or implement a more sophisticated construction of the position of the maximum (or maxima), see \secref{Max}. An important improvement is foreseen through the implementation of polarization, angle, and frequency dependence of the antenna gain when performing the coherent sum over all antennas. Apart from lifting the requirement that the flash is seen at small elevation angles and in a not too large azimuth angle range, see \secref{Cal}, this has the potential of improving the imaging quality, giving access to the direction of the electrical current moment in the VHF-source, and allowing the implementation of a true de-convolution of the beam structure giving the possibility to locate multiple sources per window.

\section{acknowledgements}
The project has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme [grant number  640130];
We furthermore acknowledge financial support from FOM [FOM-project 12PR304];
BMH is supported by the NWO [grant number VI.VENI.192.071];
AN is supported by the DFG [grant number NE 2031/2-1];
\\
This paper is based on data obtained with the International LOFAR Telescope (ILT). LOFAR~\cite{Haarlem:2013} is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Universit\'{e} d'Orl\'{e}ans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK.

\section{Data statement}\seclab{data}

The data are available from the LOFAR Long Term Archive (for access see \citep{LOFAR-LTA}).
To download this data, please create an account and follow the instructions for ``Staging Transient Buffer Board data'' at \citep{LOFAR-LTA}. In particular, the utility ``wget'' should be used as follows:
\\ {\footnotesize \verb!wget https://lofar-download.grid.surfsara.nl/lofigrid!}
\\\hspace*{7em} {\footnotesize \verb!/SRMFifoGet.py?surl="location"!}
\\where ``location'' should be specified as:
\\ {\small \verb!srm://srm.grid.sara.nl/pnfs/grid.sara.nl/data/lofar!}
\\\hspace*{7em} {\small \verb!/ops/TBB/lightning/!} \\ followed by
\\ {\small \verb!L668464_D20180921T194259.023Z_"stat"_R000_tbb.h5!} (for Flash A)
\\ {\small \verb!L792864_D20200814T143241.768Z_"stat"_R000_tbb.h5!} (for Flash B)
\\ {\small \verb!L792864_D20200814T143241.768Z_"stat"_R000_tbb.h5!} (for Flash C)
and where
``stat'' should be replaced by the name of the station, CS001, CS002, CS003, CS004, CS005, CS006, CS007, CS011, CS013, CS017, CS021, CS024, CS026, CS028, CS030, CS031, CS032, CS101, CS103, RS106, CS201, RS205, RS208, RS210, CS301, CS302, RS305, RS306, RS307, RS310, CS401, RS406, RS407, RS409, CS501, RS503, RS508, or RS509.
The source code used for TRID imaging can be found at \citep{Scholten_LofarImaging:2020, Scholten-v19:2021}.
All figures in this work have been made using the Graphics Layout Engine (GLE) \cite{GLE} plotting package.