\section{Criteria}

\section{Introduction}
The adversarial vulnerability of deep neural networks has attracted significant attention in recent years.
With slightly crafted perturbations, the perturbed natural images, or Adversarial Examples~\cite{DBLP:journals/corr/SzegedyZSBEGF13} can mislead state-of-the-art classifiers to output erroneous predictions.
Other than classification, previous works has also proved the existence of adversarial examples in varieties of computer vision tasks: Semantic Segmentation~\cite{DBLP:conf/iccv/XieWZZXY17}, Object Detection~\cite{ijcai2019-134,jia2020fooling}.   Super-Resolution~\cite{Choi_2019_ICCV} and in other fields~\cite{alzantot2018generating,carlini2018audio,huang2017adversarial,dai2018adversarial}, which arises concerns of the public, especially in safety-critical areas.

As extensive methods of generating adversarial examples appear~\cite{DBLP:journals/corr/GoodfellowSS14,carlini2017towards,DBLP:conf/iclr/MadryMSTV18}, the countermeasures are proposed at the same time~\cite{DBLP:journals/corr/abs-1711-00117,DBLP:conf/cvpr/LiaoLDPH018,DBLP:conf/cvpr/LiaoLDPH018,balunovic2020adversarial}, which seems like an arm race.
It is common to see that an defense is defeated by an newer attack within a short time, and vice versa.
However, this arm race does not help us much to understand adversarial examples.
The reasons for adversarial examples' existence remains unclear and it is still intractable to produce adversarially robust deep neural networks.
Currently adversarial training~\cite{DBLP:journals/corr/GoodfellowSS14,DBLP:conf/iclr/MadryMSTV18} is recognized the most effective way to gain adversarial robustness in practice, where the neural networks are forced to play a \textit{min-max} game.
Despite its effectiveness, it has shortages like computationally intensive, time-consuming, and non-provable.

To improve the reliability and trustworthiness of deep neural networks intrinsically, in the last few years, more and more researchers begin to focus on analyzing and explaining deep neural networks in the adversarial setting from different perspectives~\cite{DBLP:journals/corr/GoodfellowSS14,pang2020rethinking,geirhos2018imagenettrained,DBLP:conf/icml/ZhangZ19,NIPS2019_8307,NIPS2019_9394,Wang_2020_CVPR}.
With extensive efforts of the whole community, it is getting clear that why deep neural networks are sensitive to imperceptible perturbations, and some clues are found related to adversarial robustness.

Some recent surveys on adversarial examples have been proposed~\cite{Akhtar2018,wiyatno2019adversarial}, which give comprehensive discussions on existing adversarial attack and defense methods.
To our best knowledge, surveys focusing on understanding adversarial robustness do not exist even though the past few years have witnessed the great efforts made by this community.
We therefore believe that this survey on adversarial robustness can provide up-to-date findings and developments happening in this field.

Our goal in this paper is to give a brief overview of adversarial robustness.
Particularly, we carefully review and analyse adversarial robustness from different perspectives, uniquely highlight the challenges and present future research directions.

\section{Preliminaries}
\subsection{Adversarial attacks}
Adversarial attacks usually refer to finding adversarial examples for well-trained models.
Formally, finding adversarial examples is like solving a optimization problem.
We denote by $f(x, \theta): \mathbb{R}^{m} \rightarrow\{1 \dots k\}$ an image classifier that maps an input image $x$ to a discrete label set $C$ with $k$ classes, in which $\theta$ indicates the parameters of $f$.
Then the adversary tries to find a perturbation $\delta$ as follows:
\begin{equation}\begin{array}{c}
\Delta^*:=\mathop{\arg\min}\limits_{\delta} \{ \delta\in \mathbb{R}^{m}  \mid |\delta|_p \leq \epsilon, f(x+\delta)\neq f(x) \},
\end{array}\end{equation}
where $\epsilon$ is to bound the magnitude of $\delta$, and $p$ can be 0, 2 and $\infty$.
In most cases, $\epsilon$ tends to be small so that the perturbations are imperceptible to human eyes.

\subsection{Adversarial Robustness}
Currently adversarial robustness is defined as the performances of well-trained models facing adversarial examples.
In practice, the widely accepted way to gain adversarial robustness is adversarial training~\cite{DBLP:conf/iclr/MadryMSTV18}.
The idea of adversarial training is simple: it requires targeted models to train on adversarial examples in each training loops.
Formally, the objective function of adversarial training is
\begin{equation}\begin{array}{c}
\mathcal{L}(\theta, x, y):=\mathbb{E}_{(x, y) \sim D}\left[\max \limits_{\delta \in B(x,\varepsilon)} \mathcal{L}(\theta, x+\delta, y)\right]
\label{eqn:adv train}
\end{array}\end{equation}
where $B(x,\varepsilon) := \left\{x+\delta \in \mathbb{R}^{m}  \mid  \ \mid \delta\ \mid _{p} \leq \varepsilon\right\}$, and $\delta$ is obtained via gradient-based methods discussed above.
For most of time, the method used to maximize $\mathcal{L}(\theta, x+\delta, y)$ is Projected Gradient Descent~(PGD) which calculates the gradients iteratively:
\begin{equation}
\delta^{t+1}:=\Pi_{B(x,\varepsilon)}\left(x + \delta^{t}+\alpha \operatorname{sgn}\left(\nabla_{x} \mathcal{L}(\theta, x+\delta^t, y)\right)\right)-x
\end{equation}
where $\Pi_{B(x,\varepsilon)}(\cdot)$ is a projection operator projecting the input into the feasible space $B(x,\varepsilon)$, $t$ is current step, and $\alpha$ is the step size.
Observed from Eq.~(\ref{eqn:adv train}), adversarial training encourages models to predict correctly in an \mbox{$\epsilon$-ball} surrounding data points. Many variants of adversarial training have been developed from this observation~\cite{Zhang2019TheoreticallyPT,jin2020manifold}.

\section{Understanding Adversarial Robustness}

\subsection{Robustness Evaluation}
In this section, we introduce recent works on evaluation of adversarial robustness from different perspectives.

\subsubsection{Benchmarks}
Model accuracy to adversarial examples generated by attacks (FGSM, I-FGSM, FGS, rFGS, DeepFool, CW \cite{7958570}) is the direct indicator of adversarial robustness. \cite{tramer2020adaptive,croce2020reliable} show the insufficiency of such straightforward adaptive approach. \cite{croce2020robustbench} try to establish a systematic standardized benchmark for adversarial robustness by imposing restrictions on the defenses involved.

Additionally, radius of distortion represents the range of adversarial perturbation for generating successful Adversarial Examples. The magnitude of this variant is also meaningful for evaluation of model robustness. With all other variants unchanged, an adversarially-trained model that is robust to AE of larger perturbation radius usually suggests that it has higher adversarial robustness.

\cite{Mu2019MNISTCAR} introduced the MNIST-C dataset by applying visual corruptions on the MNIST testset for benchmarking out-of-distribution robustness. \cite{vaishnavi2019attention} provided two datasets with foreground attention masks that are generated from The German Traffic Sign Recognition (GTSRB) \cite{germanroadsigndata} and MS-COCO to facilitate adversarial training.

\cite{hendrycks2018benchmarking} introduced a series of widely-used visual corruptions and adversarial perturbations to generate a dataset (ImageNet-C)  for corruption robustness and a dataset (ImageNet-P) for adversarial perturbation robustness. \cite{Hendrycks2019NaturalAE} further proposed a dataset (ImageNet-A) which includes 7500 naturally occurring adversarial examples that are encapsulated as a testset for ImageNet-based classifier. After this, in \cite{Hendrycks2020TheMF}, authors introduced a visual-art-based dataset (ImageNet-R) containing renditions of 200 ImageNet classes.

 Other than static images, \cite{Gu2019UsingVT} pioneered in open-sourcing datasets for evaluating robustness to the minute transformations across video frames. Later in \cite{Shankar2019DoIC}, ImageNet-Vid-Robust and YTBB-Robust were proposed to test the robustness on temporal perturbations derived from animated videos.

\subsubsection{Theoretical Bounding}
Pioneering work \cite{Fawzi2015AnalysisOC} has provided fundamental upper bounds and an evaluation metrics for analyzing the robustness of classifiers to adversarial attacks based on the complexity of classification tasks. \cite{DBLP:journals/corr/abs-1808-01688} introduced Accuracy-Robustness Pareto Frontier (ARPF) for deep vision models to evaluate robustness towards various adversarial attacks. \cite{buzhinsky2020metrics} proposed to evaluate the robustness with the help of generative models. Within this work, authors uses latent space performance metrics, which is based on probabilistic reasoning in latent spaces of generative models, to measures the model's resistance to adversarial examples. \cite{NIPS2019_8968} used optimal transport to study the gap between the possible optimal classification accuracy and the currently achieved SOTA ones by adversarially-trained neural networks.

\cite{weng2018evaluating} proposed Cross-Lipschitz Extreme Value for Network Robustness (CLEVER), which is a Extreme Value Theory (EVT) based metric for measuring the robustness of DNN without evaluating adversarial attacking types. However, \cite{Goodfellow2018GradientMC} pointed out that gradient masking
can misguide CLEVER to overestimate the magnitude of perturbation needed to fool the model. This makes CLEVER an inaccurate evaluation of adversarial robustness. Towards improving CLEVER, \cite{DBLP:journals/corr/abs-1810-08640} provided two directions. In the first part, EVT is applied on a new twice-differentiable robustness guarantee to imply a new estimation on robustness score, namely the second-order CLEVER score. Then authors involved Backward Pass Differentiable Approximation (BPDA) to equip CLEVER, showing its broader application on DNN with non-differentiable input transformations.

\subsection{Correlation Analysis}
Adversarial robustness of DNN can be intrinsically related to many properties of the model and datasets. Models with different levels of robustness can also have different bias towards features. In this section, we list several possible directions of research on relations of adversarial robustness and model property, dataset distribution, and training settings.

\subsubsection{Feature Representation}
\cite{NIPS2019_8307} divided feature representations into two categories, namely Robust Features and Non-Robust Features, by a purely human-centric approach. Robust features are features correlated to the object label that can be recognized by human even after generation of AE. Non-robust features are the brittle and incomprehensible features from patterns that can be flipped with adversarial attack. Authors claim that adversarial vulnerability is a direct result of sensitivity to well-generalizing features in the data. By proposing an theoretical framework and experimenting on adversarial attacks, they defend their hypothesis that AEs are a natural consequence of the presence of highly predictive but non-robust features. \cite{vaishnavi2019attention} showed that utilizing foreground attention masks on training sets, which addresses the object's foreground features and eliminates the effect of noise background, can have positive contribution towards increasing PGD-based adversarial robustness.

\cite{pmlr-v97-etmann19a} quantified the relation between model robustness and interpretability of input's saliency map, finding out that the adversarial robustness and alignment of salience map is positively related when testing on non-complex datasets like MNIST. They proved, both theoretically and experimentally, that DNN models that are more robust can produce saliency maps with clearer indication of outstanding features. \cite{ding2018on} discovered the relationship between adversarial robustness and the intrinsic distribution of training data. They hypothesis that alternation on the data distribution that is semantically stable could have significantly different outcome on robustness of adversarially-trained neural networks. Authors also pioneered in raising attentions on developing proper evaluation of adversarial robustness since existing datasets cannot be suffiently reliable after semantically-loseless shift. \cite{pmlr-v97-simon-gabriel19a} found out that, when training loss converges around its minimum, the test accuracy can still be slightly improved but such increment will severely threaten DNN’s robustness because the last increment on accuracy demands a significant increase on loss gradients.

\subsubsection{Model Sparsity}
DNN are computationally expensive and therefore techniques of model pruning have been widely used to maintain the model sparsity. \cite{NIPS2018_7308} conducted research on classification tasks of DNN model to show that, both theoretically and experimentally, the model sparsity and its adversarial robustness has intrinsic relations. The authors find out that the training process of non-linear DNN model can achieve higher adversarial robustness when trainers can increase model sparsity properly through pruning techniques. However, sparse DNN that are over-pruned can be more vulnerable to adversarial attacks than original dense DNN.

\cite{rice2020overfitting} showed that overfitting training in the adversarial context will harm the model robustness. Their empirical results indicate that early stopping will help to reserve the robust performance. Similarly, \cite{Ma_2020} found the the state-of-the-art DNNs designed for large-scale image datasets are often over-parameterized on medical imaging tasks. Their experiments show that the lacking of model sparsity leads to low adversarial robustness of medical-image-based DNNs towards adversarial examples.

\subsubsection{Model Bias}
\cite{DBLP:conf/icml/ZhangZ19} attempted to discover the inherent bias of Adversarial-Trained Convolutional Neural Networks (AT-CNN) during the task of image recognition. By experimenting on salience maps and alternations of ImageNet destroying either texture or shape, authors find that AT-CNN, which has higher robustness on adversarial attacks, are more biased towards learning global structures (shapes, edges) when recognizing an object. Oppositely, normally-trained CNN with lower robustness demonstrate different bias such that it is more favorable of recognizing the object's texture information. \cite{geirhos2018imagenettrained} drew a similar conclusion stating that, although CNN are widely considered to be shape-biased when recognizing objects, their experiment results on ImageNet-trained CNN show the phenomenon that CNN models are instead strongly biased towards learning the representation of textures rather than shapes. They also find that CNN trained on ImageNet with texture information removed are of higher robustness on classification accuracy when adding various noises to the input image.

\subsubsection{Local Smoothness}
\cite{Yu2019InterpretingAE} found that adversarially-trained DNN models are shown to have more smooth and stable geometry at the peak of loss surface (the visualization of multi-dimension loss function). Such smoothness relieves the under-fitting phenomena on adversarial examples lying in the neighborhood of inputs. \cite{Fawzi2015AnalysisOC} suggested that the phenomenon of adversarial vulnerability is partially because of the low flexibility of classifiers around the varying decision boundary when the difficulty of the classification task changes. \cite{Yu2019InterpretingAE} appliec a visualization approach named decision surface to investigate adversarial robustness. Their approach confirms the intuition that adversarial examples are naturally-existed data points lying in the close neighborhood of the given inputs.

\subsubsection{Corruption Robustness}
\cite{pmlr-v97-gilmer19a} stated that adversarial robustness and corruption robustness are two manifestations of the same underlying phenomenon. Their approach builds up close connections between two concepts by both theoretical proof and empirical experiments. The result suggests that the training process of adversarially-robust DNN model should take improving performance on realistic image corruptions (e.g. ImageNet-C) as an important factor.

\subsection{Underlying Cost}
Enhancing adversarial robustness is not always favorable when taking more performance indicators into consideration. In this section, we will introduce a few widely-concerned dilemma relevant to DNN robustness.
\subsubsection{Adversarial Robustness against Accuracy to Natural Examples}
Despite the success in improving robustness of neural networks to adversarial attacks, adversarially trained models frequently fails to generalize on natural test samples. In this section, we summarize discussions on trade-off between natural accuracy and adversarial robustness.

Early works \cite{Rozsa_2016} showed that normal models with higher accuracy appears to be more robust on adversarial examples. However, authors of \cite{hendrycks2018benchmarking} discovered that more accurate models are actually less robust since their seemingly higher robustness is because of the increased natural accuracy. \cite{Balaji2019InstanceAA} proposed instance-adaptive adversarial training that manipulates the perturbation radius of every training sample to fit a more generalized model with higher natural accuracy. \cite{buzhinsky2020metrics} showed that there is a correlation between the adversarially trained DNN accuracy on clean images and its adversarial robustness in latent spaces. \cite{DBLP:journals/corr/abs-1808-01688} evaluated the recent state-of-art ImageNet-based DNN models on multiple robustness metrics. Authors figured out the linearly negative correlation between the logarithm of model classification accuracy and model robustness. Theoretical analysis on their correlation was also provided in \cite{Zhang2019TheoreticallyPT,tsipras2018robustness}.

More recently, the long standing trade-off between adversarial robustness and natural accuracy is shown to be not necessary unavoidable \cite{Stutz2019DisentanglingAR,yang2020closer}. \citeauthor{yang2020closer} argued that image datasets including CIFAR-10 and MNIST are shown to be distributionally separated and this separation implies that the existence of a robust and accurate classifier can be obtained by rounding a locally-Lipschitz function. Authors also demonstrated that using dropout can facilitate the generalization of adversarial training. Similarly, \cite{song2018improving} showed that such generalization gap of robust training can be narrowed by domain adaptation.

\subsubsection{Adversarial Robustness against Training Resources}
The SOTA adversarial training approaches boast for their high robustness. However, such training process on huge datasets like ImageNet has been only accessible in limited laboratories around the world due to high computing cost, which is partially caused by the high sample complexity during adversarial training \cite{NIPS2018_7749}.

In \cite{Shafahi2019AdversarialTF} authors proposed FREE, an efficient adversarial training algorithm that enables a much cheaper cost of time on adversarial training compared to traditional approaches. Authors argue that it is more efficient to update both the model parameters and image perturbations simultaneously in one backward pass rather than updating them separately. Further after FREE, in \cite{Wong2020FastIB}, authors proposed FAST that introduces random initialization points to make FGSM-based adversarial training as effective as the projected gradient descent based training. They also argue that most of the popular adversarial training approaches can be accelerated by applying standard techniques of efficient training for normal deep networks. However, \cite{Andriushchenko2020UnderstandingAI} argued that adding such unnecessary randomness in FAST does not prevent catastrophic over-fitting where the model drastically loses robustness on a few epoch of adversarial training. Alternatively, authors proposed GradAlign, a new regularization method that prevents catastrophic overfitting by explicitly maximizing the gradient alignment inside the perturbation set. From the perspective of Pontryagin’s Maximum Principle (PMP), authors of \cite{Zhang2019YouOP} discovered that the adversarial perturbation is mostly coupled with the weights of the first layer of the neural networks. This work introduced YOPO, an Neural Network-oriented adversarial training algorithm that splits up the adversary computation and weight updating, where the adversary computation is only focused on the first layer of the proposed network.

\section{Conclusion}
Training robust models that are secure and reliable has become a necessary step towards trustworthy applications of deep learning. In this short survey, we provide an overview on adversarial robustness and present recent approaches of robustness evaluation. After that, we study recent advances in how to analyze and understand underlying properties of adversarial robustness. Finally, we summarize the existing discussions on potential costs of improving adversarial robustness.
\bibliographystyle{named}


\section{Introduction}

The incorporation of Fe ions into the technologically relevant GaN fabricated by metal-organic vapour phase epitaxy above a concentration of 0.4\% cations leads to the self-assembly of embedded Fe$_y$N nanocrystals (NCs) with either ferromagnetic  or antiferromagnetic properties imposed by the growth conditions and by the actual incorporation of the magnetic ions into the wurtzite lattice\,\cite{Bonanni:2007_PRB, Navarro:2010_PRB}. These Fe$_y$N nanocrystals are expected to provide an efficient, robust and energy-saving platform in particular for magnetic data storage\,\cite{Coey:1999_JMMM} and for the generation of polarized spin currents. Further applications include spin current injection into the semiconductor host crystal through spin pumping, and spin current detection by inverse spin Hall effect\,\cite{Chen:2013_NatCom}. Moreover, these semiconductor/NCs hybrid structures offer the possibility to study frustrated magnetic systems and spin glass states\,\cite{Tupan:2020_materials}.\\

Particularly relevant is the $\gamma$’-Fe$_4$N phase, which has a nearly perfect minority spin polarization and is ferromagnetic up to temperatures as high as 765\,K\,\cite{Shirane:1962_PR, Kokado:2006_PRB}, making it suitable for high-density magnetic storage\,\cite{Tsunoda:2010_APEX}. Investigations on the effect of the presence of Ga in $\gamma$’-Fe$_4$N powders are also available\,\cite{Houben:2009_ChemMat, Burghaus:2010_IC, Burghaus:2011_JSSC}. When Ga is incorporated into the $\gamma$’-Fe$_4$N perovskite Pm$\bar{3}$m face-centered cubic (fcc) structure, the Ga atoms can occupy the 1$a$ (corner) or the 3$c$ (faces) Wyckoff positions. It is reported that, Ga tends to occupy the energetically favorable 1$a$ Wyckoff positions, resulting in a maximum occupancy of 90\% Ga and 10\% Fe\,\cite{Burghaus:2010_IC}. The magnetic response of the compound is largely modified through the incorporation of Ga and it ranges from ferromagnetic for $0\,<\,y\,<\,0.25$ to weakly antiferromagnetic for $y\,>0.25$, accompanied by a reduction of the magnetic moment and of the coercive field\,\cite{Burghaus:2011_JSSC, Houben:2009_ChemMat}.\\

In this frame, the growth of ordered arrays of magnetic $\gamma$’-Ga$_y$Fe$_{4-y}$N NCs embedded in GaN\,\cite{Navarro:2012_APL, Grois:2014_Nanotech, Navarro:2019_crystals}, whose size, shape and density can be tuned on demand through the amount of Ga provided during deposition, becomes particularly appealing.  By controlling the growth conditions, the incorporation of Ga atoms into the 1$a$ Wyckoff positions of the fcc $\gamma$’-Fe$_4$N lattice is expected to be promoted, altering the magnetic characteristics of the NCs and therefore, the overall magnetic properties of the layers. In order to implement these nano-objects in ferromagnetic\,\cite{Ando:2011_NatMat, Chen:2013_NatCom}, as well as in antiferromagnetic spintronic devices\,\cite{Wadley:2016_Science, Gomonay:2014_LTP}, it is imperative to identify and quantify the local arrangement of the Ga atoms in the $\gamma$’-Ga$_y$Fe$_{4-y}$N NCs embedded in the GaN host. However, the surrounding GaN matrix challenges the detection -- \textit{e.g.} by conventional element-selective x-ray absorption spectroscopy -- of Ga atoms incorporated into the NCs.\\

Here, through a combination of anomalous x-ray diffraction (AXD) and diffraction anomalous fine structure (DAFS)\,\cite{Maurizio:2020_ASS, Stragier:1992_PRL, Pickering:1993_JACS}, the local-structure of self-assembled $\gamma$’-Ga$_y$Fe$_{4-y}$N nanocrystals embedded in GaN thin layers is investigated, in order to shed light onto the correlation between fabrication parameters, local crystallographic arrangement and overall magnetic properties of the material system. It is found, that due to the conditions imposed by the growth parameters and by the resulting crystallographic environment, the Ga atoms tend to occupy the 3$c$ sites of the fcc structure, while the magnetic characteristics of the NCs -- magnetization and Curie temperature ($T_\mathrm{C}$) -- display an inverse correlation with the amount of Ga atoms occupying 3$c$ sites in the NCs.
The present results demonstrate the outstanding potential of the DAFS method for unravelling the local structure of similar magnetic phase-separated systems even when the NCs and the surrounding matrix contain the same absorber under investigation, though in a different crystallographic landscape.\\

\section{Results and Discussion}
The investigated series of GaN thin layers with embedded $\gamma$'-Ga$_y$Fe$_{4-y}$N NCs are fabricated by metal organic vapour phase epitaxy according to the procedure detailed in Refs.\,\cite{Navarro:2019_crystals, Navarro:2012_APL} by varying the Ga supply over the samples' series in order to adjust the incorporation of Ga atoms into the fcc $\gamma$’-Fe$_4$N structure of the NCs. The Fe-doped GaN thin layers, denoted from here on as Ga$\delta$FeN (Section\,\ref{exp}), are grown on a 1\,$\mu$m thick GaN buffer, deposited onto a $c$-plane sapphire (Al$_2$O$_3$) substrate. The samples are covered with a GaN capping layer with thickness ranging between 40\,nm to 300\,nm to prevent the formation of $\alpha$-Fe NCs at the sample surface\,\cite{Li:2008_JCG, Navarro:2019_PRB}. A schematic representation of the sample structure is presented in Fig.\,\ref{fig:gafen}\,(a). The nanocrystals’ size varies depending on the amount of Ga provided during fabrication and is limited to (17$\pm$8)\,nm along the $c$-growth direction, while it ranges between (50$\pm$20)\,nm and (20$\pm$8)\,nm in the layer plane\,\cite{Navarro:2019_crystals}, as exemplified in the plan-view and cross section high-resolution transmission electron microscopy (HRTEM) images of a nanocrystal shown in Figure\,\ref{fig:gafen}\,(b). The selective area diffraction pattern (SADP) of the NCs in the plan-view image allows identifying the GaN matrix and the fcc $\gamma$'-Ga$_y$Fe$_{4-y}$N phase of the NCs from the diffraction spots. Moreover, in the SADP taken along the [001] zone-axis (ZA), the in-plane epitaxial relation $\langle 110 \rangle_\mathrm{NC}$$\parallel$$\langle 1100\rangle_\mathrm{GaN}$ between the nanocrystal and the GaN matrix is evidenced. The nanocrystals density is also tuneable by the amount of Ga provided during growth and is listed in Table\,\ref{tab:param} as a function of the normalized Ga-flow. This parameter is obtained by multiplying the trimethyl-gallium (TMGa) source flow by the time over which the TMGa source is kept open during the growth of the Ga$\delta$FeN layer and of the subsequent GaN capping film, as the incorporation by diffusion of Ga into the NCs during the fabrication of the capping layer cannot be ruled out.
\begin{figure}
	\centering
	\includegraphics[width=13 cm]{Fig1.pdf}
	\caption{(a) Schematic structure of the phase-separated Ga$\delta$FeN layers. (b) Detail of a nanocrystal and its surroundings with corresponding plan-view TEM, SADP and cross-section TEM (XTEM) image. (c) Radial 2$\theta$-$\omega$ HRXRD scans of the investigated samples. (d) Detail of the scan in (c) showing the shift of the (200) $\gamma$'-Ga$_y$Fe$_{4-y}$N diffraction peak to larger angles, \textit{i.e} smaller lattice parameter along the $c$-growth direction, with decreasing normalized Ga flow.}
	\label{fig:gafen}
\end{figure}
\begin{table}
	\centering
	\caption{Investigated samples: growth parameters, lattice parameter along the $c$-growth direction calculated from HRXRD, NCs density\,\cite{Navarro:2019_crystals}, and fit parameters extracted from DAFS.}
	\begin{tabular}{ccccccc}
		\hline
		Sample & Normalized Ga-flow& Cap thickness & Lattice parameter & NC density & $(A-B)/A$ & $y^{3c}$\\
		Number & (arb. units) & (nm) & $\pm 0.005$\,($\AA$) & 10$^{9}$ (NC/cm$^2$) & & \\
		\hline
		S1 & 104  & 40 & 3.789 & 3.9 &  0.07 & 0.10 \\
		S2 & 126  & 50 & 3.786 & 7.1 & 0.13 & 0.18 \\
		S3 & 218  & 140 & 3.797 & 9.8 & 0.19 & 0.26 \\
		S4 & 312  & 230 & 3.798 & 10.2 & 0.22 & 0.30 \\
		S5 & 350  & 300 & 3.806 & 15.0 & 0.23 & 0.31 \\
		\hline
	\end{tabular}
	\label{tab:param}
\end{table}

\subsection{Structure}
The 2$\theta$-$\omega$ high-resolution x-ray diffraction (HRXRD) scans depicted in Fig.\,\ref{fig:gafen}\,(c) for the Ga$\delta$FeN layers show that besides the diffraction peaks of the crystal planes for GaN (0002) and for the underlying Al$_2$O$_3$ (0006) substrate, only an additional diffraction peak assigned to the (200) diffraction planes of the $\gamma$’-Ga$_y$Fe$_{4-y}$N embedded nanocrystals emerges. The $\langle100\rangle$ directions of the NCs lie parallel to the [0001] $c$-growth direction of the wurtzite matrix. The shift of the diffraction peak of the NCs towards larger diffraction angles with decreasing normalized Ga-flow observed in Fig.\,\ref{fig:gafen}\,(d)\,\cite{Navarro:2012_APL, Navarro:2019_crystals}, points at an increase of the lattice parameter of the NCs along the $c$-growth direction that is tentatively assigned to the inclusion of Ga into the $\gamma$’-Fe$_4$N crystal structure\,\cite{Burghaus:2010_IC}. The calculated lattice parameters from the HRXRD are given in Table\,\ref{tab:param}.

In order to confirm the inclusion of Ga atoms into the NCs, AXD and DAFS are carried out at the SIRIUS beamline at SOLEIL at the Ga $K$-edge in the energy range from 10.2\,keV to 11\,keV on the (200) diffraction plane of the NCs. In particular, the combination of the two techniques allows addressing this issue directly since, after the raw signal has been properly corrected according to the procedure detailed here below, no anomalous variation at the Ga absorption edge is expected for a phase such as the Ga-free ferromagnetic $\gamma$’Fe$_4$N one. The main challenge, which can be overcome by employing DAFS, is the presence of Ga atoms in the surrounding matrix and in the capping layer. In order to subtract the substantial x-ray fluorescence background generated by the host crystal, two regions of interest (ROIs) have been selected during the DAFS measurements shown in the inset to Fig.\,\ref{fig:dafs}\,(a): a ROI around the (002) NC diffraction peak (solid line square), and another one (of the same area and close to it) containing solely the fluorescence of the Ga in the surrounding GaN host crystal (dashed line square). The fluorescence from the GaN matrix is then subtracted from the spectra obtained from the (002) NC diffraction, leaving the contribution of Ga to the NCs (002) diffraction only. The raw Ga $K$-edge DAFS signal of the NCs including the fluorescence background due to the GaN surroundings in sample S5 is exemplarily shown in Fig.\,\ref{fig:dafs}\,(a) as a solid line, while the dotted line spectrum in the same figure represents the fluorescence background only. The differential signal presented in Fig.\,\ref{fig:dafs}\,(a) points at the inclusion of Ga atoms into the NCs.
\begin{figure}
	\centering
	\includegraphics[width=15 cm]{Fig2.pdf}
	\caption{(a) Top: Ga $K$-edge DAFS for sample S5 acquired at the (002) diffraction plane of the NCs including the fluorescence from the matrix (---) and the pure matrix fluorescence signal (- -). Inset: intensity map and ROIs from where the two spectra have been obtained. Bottom: differential spectrum obtained from the two spectra shown in (a). (b) Corrected, normalized DAFS for all samples. The Ga-flow during the fabrication increases from top to bottom. (c) Comparison between DAFS simulation (top) and experimental data of S5 (bottom) with Ga occupying 3$c$ sites. (d) Amount of Ga $y$ in 3$c$ sites in the $\gamma$’-Ga$^{3c}_y$Fe$_{4-y}$N NCs embedded in GaN plotted as a function of the normalized Ga-flow, obtained from the fabrication parameters.}
	\label{fig:dafs}
\end{figure}

A detailed analysis of the DAFS spectra has been performed by applying the following corrections to the raw signal in order to recover a signal proportional to the square of the structure factor: absorption correction from air, by the capping layer and by the Ga$\delta$FeN layer, Lorentz factor, polarization, filters, detector efficiency, and geometrical factor\,\cite{Prioetti:1999_PRB}. The data are normalized to the same pre-edge intensity for all samples. In Fig.\,\ref{fig:dafs}\,(b), the corrected, normalized differential DAFS spectra for all samples with increasing Ga-flow from top to bottom are presented. The spectra are shifted vertically for clarity. The variable $(A-B)/A$, with the diffracted intensities $A$ and $B$ marked in the spectra in Fig.\,\ref{fig:dafs}\,(b), quantifies the absorption edge and, therefore, the contribution of the Ga atoms to the structure factor of the reflection. The parameter $A$ corresponds to the normalized and corrected diffracted intensity after the edge, lying around an energy value of 10.4\,keV for all samples, while $B$ is the diffracted intensity at the edge, \textit{i.e.} at the minimum of the DAFS spectrum, located around 10.39\,keV for all samples. Therefore, the difference $A-B$ gives the size of the anomalous diffraction at the Ga $K$-edge. Since the diffracted intensity may vary from sample to sample, this difference is normalized for comparison. In this way, $(A-B)/A$ is proportional to the number of Ga atoms which participate in producing the observed reflection. For Ga-free nanocrystals, $(A-B)/A$ would be 0. Thus, this parameter quantifies how many of the Wyckoff sites -- either 1$a$ or 3$c$ -- in the Pm$\bar{3}$m fcc structure of $\gamma$'-Fe$_4$N are occupied by Ga atoms.

Based on previous literature reports\,\cite{Houben:2009_ChemMat, Burghaus:2010_IC, Burghaus:2011_JSSC}, the DAFS spectra have been first simulated in the region of the Ga $K$ absorption edge assuming that Ga occupies exclusively the 1$a$ sites in the fcc structure. However, with this approach a maximum occupancy of 100\% of 1$a$ sites by Ga in samples S5 and S4 is obtained, which is in contrast to the magnetic properties of the samples discussed in the following section and to the maximal occupation of 90\% achieved experimentally by Burghaus \textit{et al.} for $\gamma$’-Ga$_y$Fe$_{4-y}$N powders\,\cite{Houben:2009_ChemMat,
	Burghaus:2010_IC, Burghaus:2011_JSSC}.

Therefore, the occupation of 3$c$ sites by Ga atoms has been also considered. Assuming that all Ga atoms occupy the 3$c$ sites of the fcc structure, the simulated spectra shown in Fig.\,\ref{fig:dafs}\,(c) is obtained,  yielding $(A-B)/A=$\,0.72. For comparison, the experimental spectrum of S5, also presented in Fig.\,\ref{fig:dafs}\,(c), has a value of $(A-B)/A=$\,0.23. The spectral features (marked with arrows) are well reproduced by the simulation, supporting the assumption of the incorporation of Ga into 3$c$ sites. The values of $(A-B)/A\leq0.23$ indicate that a significant amount of 3$c$ sites is not
occupied by Ga atoms or that there is a number of Ga-free NCs, \textit{i.e.} having the pure $\gamma$-Fe$_4$N phase. The magnitudes of the $(A-B)/A$ parameter obtained for each sample and summarized in Table\,\ref{tab:param}, allow estimating the amount $y$ of Ga occupying 3$c$ sites in the NCs. This is plotted as a function of the normalized Ga-flow in Fig.\,\ref{fig:dafs}\,(d), showing a linear correlation between the two.
\begin{figure}
	\centering
	\includegraphics[width= 15 cm]{Fig3.pdf}
	\caption{(a) Fourier transform of the DAFS spectra and relative fits. (b) and (c) Pm$\bar{3}$m unit cell of $\gamma$’-Fe$_4$N and $\gamma$'-Ga$^{3c}_{0.3}$Fe$_{3.7}$N (samples S5 and S4), respectively, showing the different occupation sites. (d) Next nearest neighbors ($nnn$) of a Ga atom in a 3$c$ occupation site: four Fe atoms located at the vertices (1$a$) and eight Ga atoms located at the faces (3$c$). The nearest neighbors are the N atoms at the center of the unit cell.(e) Average lattice parameter along the $c$-growth direction of the NCs obtained from the (200) peak position in the 2$\theta$-$\omega$ HRXRD scan as a function of the calculated $y$. The values reported by Burghaus \textit{et al.}\,\cite{Burghaus:2010_IC} for $\gamma$’-Ga$_y$Fe$_{4-y}$N, where Ga occupies exclusively 1$a$ sites are included for comparison.}
	\label{fig:dafs2}
\end{figure}

The occupation of the 3$c$ sites by Ga atoms is further confirmed by the detailed analysis of the Fourier transform (FT) of the corrected DAFS oscillations presented in Fig.\,\ref{fig:dafs2}\,(a). A peak at the apparent distance $R$\,$\approx$\,1.5\,$\AA$ in the first Ga-N correlation shell, corresponding to about $\approx$\,2.0\,$\AA$, is observed for all samples, pointing at Ga in the 3$c$ positions. No N first neighbors at a similar distance are expected if Ga occupies 1$a$ sites. Hence, there is no evidence of a significant amount of Ga atoms occupying 1$a$ positions. The quantitative analysis of the data, which has been performed following the approach described in Refs.\,\cite{Prioetti:1999_PRB, Ciatto:2005_PRB, Ciatto:2018_JAP}, and the resulting structural parameters are collected in Table\,\ref{tab2}. For sample S1, the weakness of the fine structure signal does not allow a reliable quantitative analysis.

The unit cells of $\gamma$'-Fe$_4$N and $\gamma$'-Ga$_{0.3}^{3c}$Fe$_{3.7}$N are presented in Figs.\,\ref{fig:dafs2}\,(b) and (c) and the occupation sites and coordination distances in the Pm$\bar{3}$m crystal structure of $\gamma$'-Ga$_y$Fe$_{4-y}$N compounds are listed in Table\,\ref{tab3}. When a Ga atom occupies the 3$c$ Wyckoff position in the Pm$\bar{3}$m crystal sructure of $\gamma$'-Fe$_4$N, the nearest neighbors will be the N atoms at the 1$b$ position separated by 1.89\,$\AA$ from the Ga reference position. Further, in this case, the 3$c$ Ga atom has 12 next nearest neighbors ($nnn$) at a distance of 2.68\,$\AA$: four Fe atoms occupying 1$a$ positions (at the vertices of the structure), and eight Fe(Ga) atoms in 3$c$ positions, as schematically shown in Fig.\,\ref{fig:dafs2}\,(d) for a $\gamma$'-Ga$_{0.3}^{3c}$Fe$_{3.7}$N sample. Therefore, in the first atomic shell of Ga, an unique Ga-N distance $R_{\rm Ga-N}$ has been introduced as a variable in the fits, together with the relative Debye-Waller factor $\sigma^2_{\rm Ga-N}$. In the second shell, $R_{\rm Ga-Fe}$ is the distance between Ga at the 3$c$ position and the four Fe $nnn$ at the vertices of the cell, while $R_{\rm Ga-Ga}$ is the one to the eight Fe(Ga) $nnn$ on the faces, considering Fe and Ga as scattering centers when constructing the theoretical paths used in the fit.

The results reported in Table\,\ref{tab2} show that, for most samples, the ideal crystallographic structure is distorted and the two groups of $nnn$ are found at positions slightly different from the natural ones. The analysis of the shape of the FT and the value of the Debye-Waller factors hint at a global reduction in the amplitude of the signal with reducing the normalized Ga flow. The lineshapes of the different FT spectra are similar with the only exception of sample S2,  where a redistribution of the distances in the second shell is found. The table includes also the extracted
value for the second shell Debye-Waller factors $\sigma^2_{\rm Ga-Fe}$ and $\sigma^2_{\rm Ga-Ga}$.

\begin{table}
	\centering
	\caption{Structural parameters extracted from the quantitative DAFS analysis. $R$-range of the fits = [0.6 - 4.3] $\rm \AA$.}
	\begin{tabular}{ccccccc}
		\hline
		Sample & $R_{\rm Ga-N}$ (${\rm \AA}$) & $R_{\rm Ga-Fe}$ (${\rm \AA}$) & $R_{\rm Ga-Ga}$ (${\rm \AA}$) & $\sigma^2_{\rm Ga-N}$ ($10^{-3}{\rm \AA}^2$) & $\sigma^2_{\rm Ga-Fe}$ ($10^{-3}{\rm \AA}^2)$ & $\sigma^2_{\rm Ga-Ga}$ ($10^{-3}{\rm \AA}^2$)\\
		\hline
		S2  & 1.983 $\pm$ 0.033 & 2.856 $\pm$ 0.127 & 2.799 $\pm$ 0.083 & 0.0115 $\pm$ 0.0000 & 0.0151 $\pm$ 0.0027 & 0.0135 $\pm$ 0.0024\\
		S3  & 1.956 $\pm$ 0.026 & 2.718 $\pm$ 0.025 & 2.841 $\pm$ 0.015 & 0.0085 $\pm$ 0.0057 & 0.0067 $\pm$ 0.0029 & 0.0060 $\pm$ 0.0026\\
		S4  & 1.954 $\pm$ 0.015 & 2.711 $\pm$ 0.020 & 2.838 $\pm$ 0.013 & 0.0036 $\pm$ 0.0031 & 0.0041 $\pm$ 0.0020 & 0.0036 $\pm$ 0.0018\\
		S5  & 1.961 $\pm$ 0.013 & 2.714 $\pm$ 0.025 & 2.843 $\pm$ 0.013 & 0.0019 $\pm$ 0.0022 & 0.0037 $\pm$ 0.0017 & 0.0033 $\pm$ 0.0015\\
		\hline
	\end{tabular}
	\label{tab2}
\end{table}

\begin{table}
	\centering
	\caption{Literature values of the occupation sites and coordination distances in the Pm$\bar{3}$m fcc crystal structure of $\gamma$'-Ga$_y$Fe$_{4-y}$N compounds.}
	\begin{tabular}{cccccc}
		\hline
		Compound & Wyckoff position & Occupation &\multicolumn{3}{c}{Distance ($\AA$)}\\
		& & & 1$a$ & 1$b$ & 3$c$\\
		\hline
		$\gamma$'-Fe$_4$N\,\cite{Jacobs:1995_JAC} & 1$a$  & Fe & 3.79 & 3.28 & 2.68\\
		& 1$b$  & N &  3.28 & 3.79 & 1.89\\
		& 3$c$ & Fe & 2.68 & 1.89 & 2.68\\
		$\gamma$'-GaFe$_3$N\,\cite{Burghaus:2010_IC} & 1$a$  & Ga(0.9),Fe(0.1) & 3.80 & 3.28 & 2.68\\
		& 1$b$  & N &  3.29 & 3.80 & 1.90\\
		& 3$c$ & Fe & 2.68 & 1.90 & 2.68\\
		\hline
	\end{tabular}
	\label{tab3}
\end{table}
The change in the lattice parameter $a$ along the $c$-growth direction of the NCs, as obtained from the HRXRD spectra as a function of the Ga content $y$ is shown in Fig.\,\ref{fig:dafs2}\,(e), where the NCs in the layers fabricated with low normalized Ga-flow, \textit{i.e.} lowest $y$, display also the smallest lattice parameter. The obtained lattice parameters differ significantly from those obtained by Burghaus \textit{et al.} for $\gamma$'-Ga$_y$Fe$_{4-y}$N powders, where Ga occupies exclusively the 1$a$ sites. It was calculated, that the lattice parameter of $\gamma$'-GaFe$_3$N should be larger by about 2\% -- if Ga occupies the 3$c$ sites -- with respect to the one found if Ga occupies only 1$a$ sites\,\cite{Burghaus:2010_IC}. In the NCs considered here, the lattice parameter for the layers with $y^{3c} > 20$\% is comparable or larger than the one reported by Burghaus \textit{et al.}, as observed in Fig.\,\ref{fig:dafs2}\,(e).

\subsection{Magnetic properties}
The $\gamma$'-GaFe$_3$N ternary system is paramagnetic in the temperature range between 20\,K and 300\,K and weakly antiferromagnetic below 20\,K. Interestingly, a local ferromagnetic ordering of Fe atoms building a 13 atoms cluster when 90\% of the 1$a$ sites are occupied by Ga atoms has been reported for $\gamma$’-GaFe$_3$N powders\,\cite{Burghaus:2011_JSSC}. The magnetic response of phase-separated GaN thin layers containing Fe$_y$N embedded nanocrystals consist of two contributions: (i) a paramagnetic Brillouin-like one dominant below 50\,K, attributed to Fe$^{3+}$ ions with $S = 5/2$ substituting for Ga cations in the matrix, and (ii) a component relevant above 50\,K and associated with the ferromagnetic NCs embedded in the layers\,\cite{Navarro:2010_PRB, Gas:2019_MST, Navarro:2020_materials}. In films where antiferromagnetic $\zeta$-Fe$_2$N NCs are formed, the presence of an additional linear non-saturating component was observed\,\cite{Navarro:2010_PRB}.
\begin{figure}
	\centering
	\includegraphics[width=17 cm]{Fig4.pdf}
	\caption{(a) Normalized $M_\mathrm{NC}$ for all samples as a function of magnetic field $H$. (b) Temperature $T$ dependence of the normalized $M_{NC}$ acquired up to 390\,K (symbols). Lines: Langevin $M(T)$ fits. (c) Normalized $M_{NC}^{\mathrm{sat}}$ and $T_\mathrm{C}$ as a function of the amount of Ga $y$ in 3$c$ sites.}
	\label{fig:mag}
\end{figure}

In the samples considered here, the total magnetic moment of the Ga$\delta$FeN layers consists of the two components mentioned above, which can be separated by employing the procedure detailed in Refs.\,\cite{Navarro:2010_PRB,Navarro:2020_materials}. The amount $x_{\mathrm{Fe3+}}$ of Fe$^{3+}$ ions contributing to the Brillouin-like paramagnetism is estimated by considering that the volume of the film nominally containg Fe ions consists of the Ga$\delta$FeN and of the capping layer, where Fe diffusion is expected to take place. For all samples, $x_{\mathrm{Fe3+}}$ \,=\,(0.1$\pm$0.02)\% is obtained, a value well below the established solubility limit of $x_{\mathrm{Fe}}$\,=\,0.4\% for Fe in GaN at the growth conditions considered here\,\cite{Bonanni:2007_PRB, Navarro:2010_PRB}. Since the integrated intensity from the diffraction peak obtained in DAFS is proportional to the volume fraction of the NCs in the samples, the normalized magnetization $M_\mathrm{NC}$ of the nanocrystals is obtained by dividing the magnetic moment of the NCs at 300\,K, after the substraction of the paramagnetic component, by the diffracted intensity. The normalized $M_\mathrm{NC}$ as a function of the field applied in the in-plane configuration ($H\perp$\,[0001], $c$-growth direction) for all samples is presented in Fig.\,\ref{fig:mag}\,(a). The fast saturating behaviour of $M_\mathrm{NC}$, occuring for all samples for magnetic fields $H\geq10$\,kOe, is characteristic of phase-separated Ga$\delta$FeN and resembles a Langevin function\,\cite{Navarro:2020_materials}. However, the presence of a weak magnetic hysteresis indicates that the majority of the NCs is not in thermal equilibrium.

The temperature dependence of the normalized $M_{\mathrm{NC}}$ is shown in Fig.\,\ref{fig:mag}(b). The lines correspond to Langevin $M(T)$ fits, from where the $T_\mathrm{C}$ of the NCs is extrapolated, yielding values between (450$\pm$10)\,K and (500$\pm$10)\,K for the different samples. As an effect of the inclusion of Ga into the crystal lattice, these $T_\mathrm{c}$ values are lower than the 716\,K\,\cite{Dirba:2015_JMMM} and 767\,K\,\cite{Leinweber:1999_JAlloyComp} reported for $\gamma$'-Fe$_4$N. However, the $T_\mathrm{C}$, as well as the saturation value of $M_\mathrm{NC}^{\mathrm{sat}}$ at 300\,K, decrease with $y$, as shown in Fig.\,\ref{fig:mag}(c).
According to Burghaus \textit{et al.}, the magnetic moment of $\gamma$'-Fe$_4$N diminishes with the inclusion of Ga into 1$a$ sites. In contrast, in the nanocrystals considered here, where Ga occupies the 3$c$ sites in the lattice, the magnetization increases with increasing Ga. It is worth to underline, that the values of $M_\mathrm{NC}^{\mathrm{sat}}$ and $T_\mathrm{C}$ reported here originate from all the NCs in the samples, including those where no Ga is incorporated.

The above results show, that the ferromagnetic character and magnetization values of the $\gamma$'-Ga$^{3c}_{y}$Fe$_{4-y}$N NCs embedded in GaN are stable up to 30\% of Ga in the lattice. This suggests, that the inclusion of Ga in 3$c$ sites does not affect significantly the magnetic properties of the system, in contrast to the case of Ga in 1$a$ sites reported for $\gamma$'-Ga$^{1a}_{y}$Fe$_{4-y}$N powders\,\cite{Burghaus:2010_IC,Burghaus:2011_JSSC}.

\section{Conclusion}
It has been shown that the combination of ADX and DAFS is a powerful tool to address the amount and local structure of Ga in fcc $\gamma$’-Ga$_y$Fe$_{4-y}$N nanocrystals embedded in a wurtzite GaN host. The inclusion of Ga into the embedded NCs has been confirmed for all samples and the amount of Ga in the NCs correlates with the normalized Ga-flow employed during the fabrication of the layers: the lower the Ga-flow provided, the less Ga is included in the NCs. A detailed analysis of the FT of the oscillations in the DAFS spectra shows, that Ga preferentially occupies the 3$c$ Wyckoff sites in the $\gamma$’-Ga$_y$Fe$_{4-y}$N lattice, in contrast to previously reported results for the same material system in the powder phase\,\cite{Houben:2009_ChemMat, Burghaus:2010_IC, Burghaus:2011_JSSC}. This is attributed to the arrangement of the surrounding GaN matrix, where Ga tends to occupy lattice sites with N as nearest neighbor, \textit{i.e.} occupying the 3$c$ sites with a reduced Ga-N distance.

The normalized magnetization of the $\gamma$’-Ga$^{3c}_y$Fe$_{4-y}$N nanocrystals embedded in GaN shows an inverse dependence with the amount $y$ of Ga obtained from the DAFS measurements: the more Ga in the NCs, the higher the magnetization. The same behaviour is observed for the $T_\mathrm{C}$ of the NCs, which lies in the range between 450\,K and 500\,K. These results show, that the $\gamma$’-Ga$^{3c}_y$Fe$_{4-y}$N NCs embedded in GaN are ferromagnetic with $T_\mathrm{C}$ values well above room temperature up to a Ga content of 30\% cations, suggesting that the incorporation of Ga into the 3$c$ sites does not affect significantly the magnetic properties of the system.

The present results shed light onto the long-standing issue of the unintentional incorporation of Ga atoms into $\gamma$’-Fe$_4$N nanocrystals embedded in a GaN matrix. The challenge of distinguishing between the Ga atoms in the host matrix and the ones in the embedded NCs -- a task beyond conventional element-sensitive spectroscopies -- has been here overcome with the combination of ADX and DAFS at the Ga $K$-edge, with an approach that can be extended to a large variety of magnetic phase-separated material systems\,\cite{Dietl:2015_RMP}.

\section{Experimental Section}
\label{exp}
\subsection{Sample fabrication}
The investigated samples have been fabricated by metalorganic vapor phase epitaxy by employing trimethyl-gallium (TMGa), ammonia and cyclopentadienyl-iron (CpFe$_2$) as sources for Ga, N and Fe, respectively. The TMGa source is kept at 0$^\circ$C and the CpFe$_2$ source at 17$^\circ$C. The array of NCs is grown on a 1\,$\mu$m thick GaN buffer layer in a $\delta$-growth mode, consisting of alternated growth periods of GaN:Fe and FeN\,\cite{Navarro:2019_crystals, Navarro:2012_APL}, on a $c$-plane sapphire (Al$_2$O$_3$) substrate. A 30\,nm GaN nucleation layer is deposited to compensate the lattice mismatch between substrate and buffer.

\subsection{Characterizations}
The fabricated samples have been pre-characterized by 2$\theta$-$\omega$ radial scans taken along the [0001] growth direction by high-resolution x-ray diffraction (HRXRD) with a 8\,keV Cu $K$-$\alpha$ radiation source with 8\,keV in a Panalytical Material Research Diffractometer. The measurements are performed employing a hybrid monochromator with a 0.25$^\circ$ slit for the incoming beam and a 5\,mm slit in front of a PixCel detector working with 19 detection channels. From the 2$\theta$-$\omega$ radial scans, the NCs phase and lattice parameter along the growth direction are obtained.

Further structural characterization is carried out \textit{via} transmission electron microscopy (TEM) imaging using a JEOL JEM-2200FS TEM microscope operated at 200\,kV in high-resolution imaging (HRTEM) mode. The TEM specimens are prepared in cross-section and in plan-view by a conventional procedure including mechanical polishing followed by Ar$^+$ milling.

In order to quantify the amount of Ga and to obtain the local structure of the Ga atoms in the NCs, anomalous x-ray diffraction (AXD) and diffraction anomalous fine structure (DAFS) at the Ga $K$-edge have been carried out at the SIRIUS beamline of the SOLEIL Synchrotron\,\cite{Ciatto:2016_TSF} selecting the symmetric (002) diffraction plane, expected to contain both Fe and Ga atoms. The setup exploits a double crystal monochromator equipped with Si(111) crystals to scan the energy of the incident photons. In order to reject high-order harmonics and to focus the beam, Pt-coated mirrors are employed, and a high vacuum diffractometer\,\cite{Ciatto:2019_JSR} serves as end-station. The diffraction signal and fluorescence background are collected by using a high vacuum PILATUS3 100K-M two-dimensional detector provided by Dectris.

\emph{Ab initio} DAFS simulations in the proximity of the Ga $K$ absorption edge have been performed using the FDMNES code\,\cite{Bunau:2009_JPCM} in dipolar approximation within a non-muffin-tin finite difference method (FDM) approach\,\cite{Smolentsev:2006_RPC}; simulations are based on clusters with the atoms occupying 3$c$ and 1$a$ sites. Quantitative analysis and fitting of the Ga $K$-edge DAFS oscillations, after correction of the DAFS amplitudes and phase shifts \cite{Prioetti:1999_PRB}, are carried out by using the FEFF package\,\cite{Newville:2001_JSR}, assuming the occupation of 3$c$ sites by Ga.

The magnetic properties are investigated in a Quantum Design superconducting quantum interference device (SQUID) MPMS-XL magnetometer equipped with a low field option at magnetic fields $H$ up to 70\,kOe in the temperature range between 2\,K and 400\,K following an experimental protocol for small signals elaborated to eliminate artifacts and to overcome limitations associated with integral SQUID magnetometry\,\cite{Gas:2019_MST, Sawicki:2011_SST}.

\section*{Conflicts of interest}
The authors declare no conflicts of interest.

\section*{Acknowledgements}
The authors acknowledge financial support by the Austrian Science Fund FWF project numbers V-478 and P31423 and by the CALIPSO project of the European Union. The experiment at the SIRIUS beamline benefited from the SOLEIL beam time allocation No. 20180261.


\section{Introduction}
Restoring the image quality degraded by  blur and noise is an active research topic with applications in astronomy~\cite{roddier1999}, microscopy~\cite{Booth5788}, medical imaging~\cite{FOSTER}, machine vision~\cite{sonka2014}, \emph{etc}.
The blur can be caused by the optical aberration due to light propagation through inhomogeneous medium, by imperfection of the optics of imaging system, by the motion of the camera and/or the object.
In the simplest cases, the effect of the blur can be represented by convolution of the object with a spatially invariant PSF, and a number of approaches dealing with isoplanatic deconvolution with or without \emph{a priori} knowledge of PSF from one or multiple images is known~\cite{Cannon,AyersDainty,Schulz:93,Yaroslavsky:94,Wilding2017}.
In the anisoplanatic case, \emph{e.g.} in deep microscopy imaging~\cite{NaJi} or images deteriorated by motion blur~\cite{Dai2008,Hirsch2010,Sroubek2016}, the shape of the PSF is not constant over the image, the isoplanatic algorithms fail, and anisoplanatic deconvolution algorithms able to evaluate the PSF at every point of the image are required.

However, it is computationally inefficient to implement a method that calculates a unique correction for every point~\cite{Seitz2009,Hirsch2010,Denis2011}. A group of methods uses a global decomposition of the spatially varying~(SV) PSF in order to work with a limited number of spatially static PSFs. Flicker and Rigaut~\cite{Flicker2005} approximate the SV PSF by decomposition into a number of spatially invariant modes. Miraut and Portilla~\cite{Miraut2012} and later Sroubek \textit{et al.}~\cite{Sroubek2016} decompose the the SV PSF into spatially invariant modes by means of the singular value decomposition. An advantage of decomposition of the SV PSF by the singular value decomposition is that it is globally optimal with respect to the approximation error. On the other hand, the use of global decomposition methods comes with large computational load~\cite{Denis2015}.

Another approach to limit the SV PSF to a number of spatially static PSFs is division of the image into isoplanatic subsections. This method was first proposed by Trussel and Hunt~\cite{TrusselHunt}; they divided the image into disjoint subsections with local spatially static PSFs. This method was later used by Costello and Mikhael~\cite{CostelloMikhael} in combination with modified Wiener filtering in order to correct for anisoplanatic aberrations. The implementation of disjoint isoplanatic subsections is known to lead to the appearance of boundary artifacts~\cite{Denis2015}. Nagy and O'Leary~\cite{Nagy97} implemented the overlap-add method, where overlapping isoplanatic subsections were used for deconvolution. The use of this overlap-add method reduces the effects of boundary artifacts in the restored image. This method was later improved by Denis \textit{et al.}~\cite{Denis2015}.

A popular trend in image reconstruction is the use of neural networks. Shajkofci and Liebling~\cite{Shajkofci2018} used a Convolutional Neural Network that requires an optical parametric model. The optical model is used as prior knowledge of the PSF. Another approach is the use image specific \textit{a priori}. Sureau \textit{et al.}~\cite{sureau2019} used a Deep Neural Network architecture to learn supervised setting parameters adapted for galaxy image processing.

The problem becomes even more challenging when no prior knowledge of the PSF or object is available, as is generally the case with medium induced aberrations in microscopy. Hirsch \textit{et al.}~\cite{Hirsch2010} used a similar overlap-add method to perform blind anisoplanatic deconvolution with multiple images of the same object. In blind deconvolution it is the goal to estimate both the object and the aberration or (local) PSF. Many of the blind anisoplanatic correction methods found in literature~\cite{TrusselHunt,LagendijkBiemond1991,GuoLeeTeo1997} use the assumption that the image can be divided into isoplanatic subsections. This group of methods divides the image(s) into isoplanatic subsections and apply conventional isoplanatic blind deconvolution methods in order to estimate the object locally. The performance of methods in this group depends on the speed at which the SV PSF is changing in the image. The assumption that the PSF is spatially static on a small subsection does not always hold. We propose weighted multi-frame deconvolution, a method that prioritizes subsections with highly isoplanatic aberrations for deconvolution. A small weight is assigned when large variation occurs in the PSF locally and parts with small PSF variation are assigned a large weight. These weights determine the contribution to deconvolution of an image frame.

Using multiple images of the same object with different aberrations comes with another problem. The difficulty is that translations in the PSFs cause morph of the images, \textit{i.e.} parts of the object appear at different distances in the images. For local multi-frame deconvolution it is required that the isoplanatic subsections in the images correspond to the same part of the object. This paper demonstrates an algorithm, Tangential Iterative Projections Algorithm for Anisoplanatic Aberrations with Adaptive finite support and Weighted multi-frame deconvolution (TIP4AW), that is specifically designed to deal with this problem.

TIP4AW processes multiple images of the same object in batch mode. All images are considered at once in an iterative process of four steps. In the first step the object is locally estimated by weighted multi-frame deconvolution and the overlap-add method is use to form the estimated object. The second step is the projection of the object on the feasible set of all possible objects. In the third step all local PSFs are estimated by single-frame deconvolution. Isolation of the corresponding subsection in this step is performed by multiplication with an apodization kernel. The fourth and last step is the projection of the estimated local PSFs on the feasible set of all possible PSFs, using adaptive PSF support. Numerical simulations are conducted to show the robustness to noise of this algorithm. Moreover the performance of the algorithm is compared to a state-of-the art algorithm found in literature.

\section{Local blind multi-frame deconvolution}
As described above, a common approach to blind anisoplanatic deconvolution is the use of isoplanatic subsections. On isoplanatic subsections conventional blind deconvolution is performed and the local results are combined to form a corrected image. This section describes anisoplanatic and isoplanatic image formation and elaborates on deconvolution by linear filtering. Furthermore, it is explained how isoplanatic deconvolution can be used in order to deconvolve the image from anisoplanatic aberrations using overlap-add.

\subsection{Image formation}
In an incoherent imaging system the formation of an image can be seen as the the summation of the contribution of all point sources in the object, $o_{k,l}$, where $k,l$ are the discrete spatial coordinates in the object plain. Discrete coordinates are used, because an image is in general spatially samples by a CCD. Let $\tilde{h}_{m,n,k,l}$ be the normalized contribution, or PSF, of point source $o_{k,l}$, where $m,n$ are the spatial coordinates in the image plain. The formation of an image, $i$, is given by
\begin{equation}
    i_{m,n} = \sum_{k=1}^M\sum_{l=1}^N o_{k,l} \tilde{h}_{m,n,k,l} + e_{n,m}, \qquad m=1,...,M, \quad n=1,...,N ,
    \label{eq:DiscAnisoImageAndnoise}
\end{equation}
where $e_{n,m}$ is measurement noise.

The image model of Eq.~\eqref{eq:DiscAnisoImageAndnoise} describes anisoplanatic image formation where all PSFs can be different. In the case of (almost) spatially invariant PSF, the isoplanatic imaging model is used:
\begin{equation}
    i_{m,n} = \sum_{k=1}^M\sum_{l=1}^N o_{k,l} h_{m-k,n-l} + e_{m,n},
    \label{eq:DiscIsoImageAndnoise}
\end{equation}
Equation~\eqref{eq:DiscIsoImageAndnoise} can be recognized as convolution of the object and the PSF, and is denoted in short notation by $i_{m,n}=o_{m,n}\ast h_{m,n} + e_{m,n}$. For the isoplanatic model the convolution theorem can be used and in the spatial frequency domain the convolution in Eq.~\eqref{eq:DiscIsoImageAndnoise} becomes the point-wise multiplication,
\begin{equation}
    I_{m,n} = O_{m,n}\times H_{m,n} + E_{m,n},
    \label{eq:OTF}
\end{equation}
where the capital letter describes a two dimensional Discrete Fourier Transform. The PSF spectrum, $H_{m,n}$, in Eq.~\eqref{eq:OTF} is referred to as the optical transfer function ~(OTF).

\subsection{Deconvolution}
Assuming that Eqs.~\eqref{eq:DiscIsoImageAndnoise} and \eqref{eq:OTF} hold and that the PSF is known, a na\"{i}ve approach to deconvolution is point-wise division of the image spectrum and OTF, $\hat{O}_{m,n}=I_{m,n}/H_{m,n}$. This approach, however, is known to amplify noise and usually leads to a decrease in the signal-to-noise ratio (SNR)~\cite{DecLab2}. A more optimal method is given by the Wiener filter~\cite{Helstrom67},
\begin{equation}
    \hat{O}_{m,n} = \frac{(H_{m,n})^*I_{m,n}}{|H_{m,n}|^2+\frac{1}{\text{SNR}}},
    \label{eq:Wiener}
\end{equation}
where $H^*$ denotes the complex conjugate of $H$.
Equation~\eqref{eq:Wiener} was later modified by Gonzalez and Woods~\cite{GonWoods} for situations where the SNR is unknown.

\subsection{Blind multi-frame deconvolution}
In some applications the exact shape of the PSF is not known, and both the object and PSF should be estimated. Blind deconvolution can be posed as an ill-conditioned minimization problem~\cite{Chaudhuri}. Ayers and Dainty~\cite{AyersDainty} were one of the first to propose a practical algorithm that can solve this problem by alternating between solutions for the object and the PSF and applying constraints in the Fourier domain and spatial domain. Other blind deconvolution methods make use of the temporal variation of the PSF. More than one image of the same object, but with different aberration are used to reconstruct the object and PSF. Such methods are referred to as blind multi-frame deconvolution (BMFD) or blind multi-channel deconvolution. In BMFD a set of $S$ images $\{i_{m,n|s}\}$, with corresponding unknown PSFs $\{h_{m,n|s}\}$, are deconvolved. All images are obtained from the same object,
\begin{equation}
    i_{m,n|s} = \sum_{k=1}^M\sum_{l=1}^N o_{k,l} h_{m-k,n-l|s} + e_{m,n|s}.
    \label{eq:Multi-FrameFormation}
\end{equation}
Different images from the same object can be obtained if the aberration is time-variant. This is often the case in \textit{in vivo} microscopy~\cite{Wilding2018} and in imaging trough turbulence~\cite{Roggemann95}. Schulz~\cite{Schulz:93} demonstrated in 1993 a practical implementation of an algorithm that uses more than one image. Schulz's algorithm is iterative in nature and uses penalized maximum likelihood estimation and a parameterization  of the PSF in order to jointly estimate the object and PSF. A year later, Yaroslavsky \textit{et al.}~\cite{Yaroslavsky:94} showed that a Wiener-like filter,
\begin{equation}
    \hat{O}_{m,n} =
    \frac{\sum_{s=1}^S (H_{m,n|s})^*I_{m,n|s}}{\sum_{s=1}^S |H_{m,n|s}|^2},
    \label{eq:Multi-Frame-Wiener}
\end{equation}
can be used to estimate the object from more than one image when the PSFs are known. In 2017, Wilding \textit{et al.}~\cite{Wilding2017} demonstrate an Ayers and Dainty-like algorithm that uses the filter in Eq.~\eqref{eq:Multi-Frame-Wiener} in order to estimate the object from multiple images. Equation~\eqref{eq:Multi-Frame-Wiener} fails at locations where common zeros in the PSF spectrum occur. Therefore, it is best to record as much images with different PSFs as possible. Since common zeros in the PSF spectrum are more likely to occur in one dimensional~(1D) signals, this algorithm is not optimal in a 1D analysis.

\subsection{Spatial variation of the PSF}
The assumption that the PSF is spatially static is a necessary condition for Eq.~\eqref{eq:OTF} to hold. For anisoplanatic aberrations it can be assumed that Eq.~\eqref{eq:OTF} holds for small enough regions in the image, isoplanatic subsections. This allows the conventional deconvolution techniques to be applied locally. Local results are combined to form the estimated object.

Trussel and Hunt~\cite{TrusselHunt} combined the local results by adding disjoint isoplanatic sections together. This method results in boundary artefacts near the edges of isoplanatic subsections. Nagy and O'Leary~\cite{Nagy97} used the overlap-add method where the isoplanatic subsections overlap each other. An isoplanatic subsection was defined by windows, $W_{\tilde{p}}[m,n]$, where the index $\tilde{p}$ is the window number and $m,n$ are discrete spatial coordinates of the window. These windows must have the following property:
\begin{equation}
    \sum_{\tilde{p}=1}^{\tilde{p}_t} W_{\tilde{p}}[m,n] = 1,
\end{equation}
where $\tilde{p}_t$ is the total number of windows. In overlap-add, the image is locally (de)convolved with a local PSF by means of applying a window $W_{\tilde{p}}$. The global result is found after summation of local results.

\section{Algorithm description}
TIP4AW uses tangential iterative projections~\cite{Wilding2017} and the overlap-add method~\cite{Nagy97,Denis2015} in conjunction with novelties weighted multi-frame deconvolution, Section \ref{:ObjectFourier}, and adaptive PSF support, Section \ref{:PSFSpatial}, to handle anisoplanatic aberrations in multiple images of the same object. It is an algorithm that is demonstrated to be robust to noise and morph.

The method to correct for anisoplanatic aberrations in blurred images requires multiple images of the same object. These images are divided into overlapping subsections such that overlap-add can be applied after deconvolution to avoid the appearance of boundary artifacts. Let $S$ be the number of available images, and let $i_{p,q|s}[m,n]$ denote subsection $p,q$ of image $i_s$. The discrete spatial coordinates $m,n$ are not denoted anymore for ease of reading. Section~\ref{::Discussion} provides more details on the relation between the number of subsections and the quality of the estimated object. Figure~\ref{fig:ImageSubsections}(a) shows an image that is divided into subsections. The subsection $i_{p,q}$ is shaded in red. Every subsection has four quadrants. Subsections $i_{p,q}$ and $i_{p,q+1}$ have two overlapping quadrants. The local PSF for subsection $i_{p,q}$ in image $s$ is denoted by $h_{p,q|s}$. TIP4AW proceeds by iteratively improving estimates of the object and PSF.

Each iteration consists of four steps. In the first step weighted multi-frame deconvolution is used in order to estimate the object locally. Overlap-add is then used to combine local results and find the estimated object. In the second step, the estimated object is improved by applying spatial domain constraints. During the third step, the local PSF is estimated after applying local apodization on the object. In the forth and final step, PSF constraints in the spatial domain are applied to improve the estimated local PSFs.

\subsection{Object step in Fourier domain}
\label{:ObjectFourier}
The algorithm is initialized by assuming that all local PSFs are delta functions, $h_{p,q|s} = \delta[\floor*{m/2},\floor*{n/2}]$, where $\floor*{\cdot}$ is the floor function. As a result, the first estimated object is a pixel-wise average of all images. During the $k$th iteration, weighted multi-frame deconvolution is used to estimate the local object
\begin{equation}
    \bar{o}_{p,q}^{(k)} = \mathcal{F}^{-1} \left\{ \frac{\sum_{s=1}^S a^{(k)}_{p,q|s} \left(H_{p,q|s}^{(k)}\right)^*I_{s}}{\left[ \sum_{s=1}^S a^{(k)}_{p,q|s} |H_{p,q|s}^{(k)}|^2\right]_{>\epsilon\bar{a}_{p,q}}} \right\},
    \label{eq:WMFD}
\end{equation}
where $\mathcal{F}^{-1}\{\cdot\}$ denotes the inverse Discrete Fourier Transform and multiplication and division are element-wise. The notation $N/[D]_{>\alpha}$ means that division is performed if the absolute value of the denominator is larger than $\alpha$, else the denominator is set to infinity, such that the result becomes zero. In Section \ref{sec:weights} is explained how the weights, $a^{(k)}_{p,q|s}$, are determined. $\bar{a}_{p,q}$ is the average $\frac{1}{S}\sum_{s=1}^Sa^{(k)}_{p,q|s}$. The parameter $\epsilon>0$ is related to the SNR, if the SNR is unknown, a small $\epsilon$ can be chosen. Every local object is multiplied element-wise by a window, $W_{p,q}$. The windows are chosen such that $\sum_{p=1}^P\sum_{q=1}^Q W_{p,q}=1$ for every pixel. Here, the windows are congruent to a bi-linear interpolation mask:
\begin{equation}
   W_{p,q}[m,n] = \max\left\{1-\frac{|c_p-m|}{l_m/2},0 \right\}\cdot\max\left\{1-\frac{|c_q-n|}{l_n/2},0\right\},
   \label{eq:mask}
\end{equation}
where $C[c_p,c_q]$ is the centre of subsection $p,q$ and $l_m$ and $l_n$ are the vertical and horizontal lengths of the subsection respectively, and $\max\{\cdot\}$ is a pixel-wise operation. Figure~\ref{fig:ImageSubsections}(b) shows the bi-linear interpolation mask in subsection $p,q$. Extra care must be taken when dealing with exterior subsections, subsections where $p$ or $q$ is equal to $1,P$ or $Q$. All interior subsections are overlapped by four neighbouring subsections and therefore the sum of bi-linear interpolations masks is one uniformly. This is, however, not the case for exterior subsections. It is possible to use local PSFs on the edge of the images. However, it is more difficult to calculate local PSFs on the edge. Instead of calculating the local PSF on the edge, the closest interior PSF is chosen. This is similar to normalizing the windows in Eq.~\eqref{eq:mask}, $W_{p,q}[m,n] = {W_{p,q}[m,n]}/{\sum_{p=1}^P\sum_{q=1}^QW_{p,q}[m,n]}$.
After multiplication of the local object with the corresponding window, the object subsection is formed, Fig.~\ref{fig:ImageSubsections}(c). Summation of these object subsection creates the estimated object
\begin{equation}
    \bar{o}^{(k)} = \sum_{p=1}^P\sum_{q=1}^Q \bar{o}^{(k)}_{p,q}W_{p,q}.
    \label{eq:objectsynthesis}
\end{equation}

\begin{figure}[!ht]%
\centering
\includegraphics[width=\textwidth]{ImageSectionsCombi.pdf}
    \caption{(a)~The deteriorated image is divided into $P$ by $Q$ overlapping subsections. Subsection $i_{p,q}$ is highlighted in red. Every subsection has four quadrants. (b)~The bi-linear interpolation masks, Eq.~\ref{eq:mask}, are used to create one object from all overlapping subsections. Every local object is multiplied with a corresponding mask or window. The sum of these masks is one uniformly. The values of these masks represent the contribution to a point as in bi-linear interpolation. (c)~The object subsection is the local object multiplied with the interpolation mask. The estimated object is the sum of all object subsections.}%
    \label{fig:ImageSubsections}%
\end{figure}

\subsection{Object step in spatial domain}
\label{:ObjectSpatial}
In the second step, the object is projected on the feasible set of all possible object. This is done by applying the following object constraints. The intensity of the estimated object must be non-negative and the object is normalized such that $\left|\left|o^{(k)}\right|\right|_{1} = 1$.
\begin{equation}
    o^{(k)} = \frac{\max\left\{ \text{real}\left(\bar{o}^{(k)}\right),0 \right\}}{\left|\left|\max\left\{ \text{real}\left(\bar{o}^{(k)}\right) ,0\right\}\right|\right|_{1}}.
\end{equation}

\subsection{PSF step in Fourier domain}
\label{:PSFFourier}
In this step, the local PSFs are estimated by local single-frame deconvolution. The image can be deconvolved with the object locally by isolating the correct part in the object. A Gaussian apodization kernel is used to isolate the corresponding part of the object. Apodization also reduces edge effects, such as ringing and ghosts~\cite{DecLab2}. Apodization centred in subsection $p,q$ with center $C[c_p,c_q]$ is performed by multiplication of the object with kernel
\begin{equation}
    K^{w}_{p,q}[m,n] = e^{-\frac{(m-c_p)^2+(n-c_q)^2}{w^2}},
\end{equation}
where $w$ is the width parameter of the Gaussian apodization kernel. $K^{w}_{p,q}$ is the apodization kernel corresponding to local object, $o_{p,q}$, and local PSFs, $h_{p,q|s}$. Then the estimated local PSF is given by single-frame deconvolution of the image and apodized object.
\begin{equation}
   \bar{h}_{p,q|s}^{(k+1)} = \mathcal{F}^{-1} \left\{  \frac{I_s}{ \left[ O_{K^{w}_{p,q}}^{(k)} \right]_{>\epsilon}} \right\},
   \label{eq:PSFest}
\end{equation}
where $O_{K^{w}_{p,q}}^{(k)} = \mathcal{F}\{o^{(k)} K^{w}_{p,q}\}$.

\subsection{PSF step in the spatial domain}
\label{:PSFSpatial}
The final step is the projection of the estimated local PSFs $\bar{h}_{p,q|s}^{(k+1)}$ on the feasible set of PSFs. All PSFs are constraint to be limited in size. Since morph occurs in the images, the local PSFs are unlikely to be centralized, therefore adaptive PSF support is used. Let $\mathbb{X}_{p,q|s}^{(k)}$ denote the adaptive PSF support for local PSF $h_{p,q|s}^{(k+1)}$. The adaptive PSF support is given by a circular area with radius, $r_{\mathbb{X}}$, and its centre is adaptive to the centre of mass of the local PSF. The set of all possible PSFs, $\mathcal{H} = \left\{h|h_{p,q|s}[m,n]=0 \text{ for } [m,n]\not\in\mathbb{X}_{p,q|s}^{(k)}\right\}$, is limited to the size of the support constraint. Furthermore, the PSF is an intensity distribution that contains only non-negative real values and is normalized:
\begin{equation}
    h_{p,q|s}^{(k+1)} = \frac{\max\left\{ \text{real}\left(\bar{h}_{p,q|s}^{(k+1)}\right),0 \right\}}{\left|\left|\max\left\{ \text{real}\left(\bar{h}_{p,q|s}^{(k+1)}\right) ,0\right\}\right|\right|_{1}}.
\end{equation}
As a result the PSFs can be written as the sum of two components. The tangential component of $h_{p,q|s}^{(k+1)}$ denoted as $\left(h_{p,q|s}^{(k+1)}\right)_{\parallel}$, is contained in $\mathcal{H}$. The normal component $\left(h_{p,q|s}^{(k+1)}\right)_{\perp}$ consists of the pixels of the PSF that are not contained in $\mathcal{H}$.
\begin{equation}
    h_{p,q|s}^{(k+1)} = \left(h_{p,q|s}^{(k+1)}\right)_{\parallel} + \left(h_{p,q|s}^{(k+1)}\right)_{\perp}, \quad \left(h_{p,q|s}^{(k+1)}\right)_{\parallel}\in\mathcal{H},\text{ }\left(h_{p,q|s}^{(k+1)}\right)_{\perp} = 0.
\end{equation}
Finally, the local PSF spectrum is given by the Discrete Fourier Transform of the estimated local PSF
\begin{equation}
    H_{p,q|s}^{(k+1)} = \mathcal{F} \left\{ h_{p,q|s}^{(k+1)} \right\}.
\end{equation}
In TIP4AW, these four steps are repeated consecutively. Algorithm \ref{alg1} shows an outline of TIP4AW.

\begin{algorithm}[tbh]
	\caption{TIP4AW}
	\label{alg1}
	\begin{algorithmic}
		\REQUIRE $\{i_1,...,i_S\}$, $k_{\text{max}}$, $r_{\mathbb{X}}$, $P$, $Q$, $w$, $\Delta w$, $\epsilon$, $p_s$
		\RETURN $\hat{o}$, $\{\hat{h}_{p,q|s}\}$
		\STATE initialize $h_{p,q|s}=\delta[\floor*{m/2},\floor*{n/2}]$
		\FOR{$k = 1,...,k_{\text{max}}$}
		\FOR{$p=1,...,P$\tikzmark{top1}}
		\FOR{$q=1,...,Q$}
		\STATE $\bar{o}_{p,q}^{(k)} = \mathcal{F}^{-1} \left\{ \frac{\sum_{s=1}^S a^{(k)}_{p,q|s} \left(H_{p,q|s}^{(k)}\right)^*I_{s}}{\left[\sum_{s=1}^S a^{(k)}_{p,q|s} |H_{p,q|s}^{(k)}|^2\right]_{>\epsilon\bar{a}_{p,q}}} \right\}$
		\ENDFOR
		\ENDFOR            \tikzmark{bottom1}
		\STATE $\bar{o}^{(k)} = \sum_{p=1}^P\sum_{q=1}^Q \bar{o}^{(k)}_{p,q}W_{p,q}$ \tikzmark{top2}
		\STATE $o^{(k)} = \frac{\max\left\{ \text{real}\left(\bar{o}^{(k)}\right),0 \right\}}{\left|\left|\max\left\{ \text{real}\left(\bar{o}^{(k)}\right) ,0\right\}\right|\right|_{1}}$    \tikzmark{bottom2}
		\STATE $O_{K^{w}_{p,q}}^{(k)} = \mathcal{F}\{o^{(k)} K^{w}_{p,q}\}$
		\STATE $O_{K^{w+\Delta w}_{p,q}}^{(k)} = \mathcal{F}\{o^{(k)} K^{w+\Delta w}_{p,q}\}$
		\FOR{$p=1,...,P$\tikzmark{top3}}
		\FOR{$q=1,...,Q$}
		\FOR{$s=1,...,S$}
		\STATE $\bar{h}_{p,q|s}^{(k+1)} = \mathcal{F}^{-1} \left\{  \frac{I_s}{ \left[ O_{K^{w}_{p,q}}^{(k)} \right]_{>\epsilon}} \right\}$
		\STATE $\tilde{h}_{p,q|s}^{(k+1)} = \mathcal{F}^{-1} \left\{  \frac{I_s}{ \left[ O_{K^{w+\Delta w}_{p,q}}^{(k)} \right]_{>\epsilon}} \right\}$
		\ENDFOR
		\ENDFOR
		\ENDFOR \tikzmark{bottom3}
		\STATE $ h_{p,q|s}^{(k+1)} = \frac{\max\left\{ \text{real}\left(\bar{h}_{p,q|s}^{(k+1)}\right),0 \right\}}{\left|\left|\max\left\{ \text{real}\left(\bar{h}_{p,q|s}^{(k+1)}\right) ,0\right\}\right|\right|_{1}}$ \tikzmark{top4}
		\STATE $h_{p,q|s}^{(k+1)}\in\mathcal{H} = \left\{h|h_{p,q|s}[m,n]=0 \text{ for } [m,n]\not\in\mathbb{X}_{p,q|s}^{(k)}\right\}\qquad$  \tikzmark{right}\tikzmark{bottom4}
		\STATE $ a^{(k)}_{p,q|s} = ||h_{p,q|s}-\tilde{h}_{p,q|s}||_F^{-2p_s}$
		\STATE $H_{p,q|s}^{(k+1)} = \mathcal{F} \left\{ h_{p,q|s}^{(k+1)} \right\}$
		\ENDFOR
		\STATE $\hat{o} = o$ and $\hat{h}_{p,q|s} = h_{p,q|s}$
	\end{algorithmic}
	\AddNote{top1}{bottom1}{right}{Object step in Fourier domain (1)}
	\AddNote{top2}{bottom2}{right}{Object step in spatial domain (2)}
	\AddNote{top3}{bottom3}{right}{PSF step in Fourier domain (3)}
	\AddNote{top4}{bottom4}{right}{PSF step in spatial domain, for $h$ and $\tilde{h}$ (4)}
\end{algorithm}

\subsection{Weight determination}
\label{sec:weights}
They key novelty of the algorithm is weighted multi-frame deconvolution. For weighted multi-frame deconvolution the weights, $a_{p,q|s}$, represent the contribution of a frame, $i_s$, that is used in Eq.~\eqref{eq:WMFD} in order to estimate the object. If all weights are chosen $a_{p,q|s}=1$ and $\epsilon=0$ then Eq.~\eqref{eq:WMFD} is exactly the filter in the Fourier domain given in Eq.~\eqref{eq:Multi-Frame-Wiener}. Deconvolution in the Fourier domain requires the assumption that Eqs.~\eqref{eq:DiscIsoImageAndnoise} and \eqref{eq:OTF} are valid. This assumption only holds locally and it is possible that isoplanatic subsection do not coincide, especially when there is a lot of morph in the images. For optimal performance of the algorithms it is favorable that the image subsections are highly isoplanatic, that is with less variant blur. Weighted multi-frame deconvolution makes that highly isoplanatic image subsections contribute more to the deconvolution than image subsections that contain large anisoplanatic aberrations.

When an image frame is locally highly isoplanatic the corresponding weight will be large and likewise a small weight will be assigned to a frame when it is locally anisoplanatic. The metric that is associated with the scale of the variation of the SV PSF is called \textit{isoplanatism}, or PSF variation. The process of finding optimal values for these weights has not been researched yet.

In this paper a method is proposed to demonstrate the effect weighted multi-frame deconvolution has on the image quality of the estimated object. Isoplanatism is measured by comparing the similarly between the estimated local PSF, $h_{p,q|s}$, and complementary estimated local PSF, $\tilde{h}_{p,q|s}$. The local PSFs, as calculated in Eq.~\eqref{eq:PSFest}, depend on the width of the apodization kernel, $K^w$. For the calculation of the estimated local PSFs an apodization kernel, $K^w$, is used and for the complementary estimated local PSFs a wider apodization kernel, $K^{w+\Delta w}$, is used. Since the PSF is spatially variant, the complementary local PSF will be different than the local PSF. The difference is a measure for the local isoplanatism. The weights are given by the following formula:
\begin{equation}
    a_{p,q|s} = ||h_{p,q|s}-\tilde{h}_{p,q|s}||_F^{-2P_s},
    \label{eq:isoplanatism}
\end{equation}
where $||\cdot||_F$ denotes the Frobenius norm and $p_s$ is the isoplanatism sensitivity. This metric was first reported in the Master of Science thesis of one of the authors~\cite{Ketterij2019}.

\section{Results}
The performance of TIP4AW is compared with the space-variant Online Blind Deconvolution with Efficient Filter Flow (OBD/EFF) algorithm in~\cite{Hirsch2010}. This algorithm was implemented in the best of the author's ability. Thereafter the robustness to noise of both algorithms is compared. A MATLAB implementation of TIP4AW and a set of test images can be found here \cite{TIP4AW}.

\subsection{Algorithm comparison}
Before both algorithms can be compared, it is necessary to define an adequate quality metric. As for simulated data the estimated object can be compared to the real object or ground truth. Due to large morph in the images the authors do not recommended to use an error based metric such as the PSNR or mean squared error.

Instead a metric called Fourier Ring Correlation (FRC) which measures the normalized cross-correlation between the real object and the estimated object on rings, $r_n$, in the Fourier domain was used. FRC is originally a method to determine the resolution of reconstructed images in electron microscopy~\cite{frc} and later used in~\cite{banterle} as a resolution criterion in super-resolution microscopy. In this thesis FRC will be used to measure the similarity between the real object and the estimated object,
\begin{equation}
    FRC[r] = \frac{\sum_{r_n\in r}O[r_n]\cdot \hat{O}^*[r_n] }{\sqrt{\sum_{r_n\in r}\left|O[r_n]\right|^2\cdot\sum_{r_n\in r}\left|\hat{O}[r_n]\right|^2}}.
\end{equation}
An FRC curve close to one corresponds to a close match of the real object and estimated object. Since the algorithms are independent of the imaging system, the resolution of the imaging sensor is unknown. The quality of the estimated object can be expressed by a single value that is the largest ring, $r_{n,\text{max}}$, where the FRC was larger than a reference line, the ${2\sigma}$ curve, used in both~\cite{frc,banterle}.

Two numerical experiments were conducted wherein the object was estimated by TIP4AW and OBD/EFF. In each experiment 30 images of the same object were provided to the algorithms. For the experiments the size of the PSF was estimated at 13 pixels and the images are divided into 49 overlapping subsection. For TIP4AW the apodization kernels, $K^{35}$ and $K^{49}$, noise parameter, $\epsilon = 10^{-4.4}$, and isoplanatism sensitivity, $P_s=1.5$, were used in order to maximize the performance. For OBD/EFF a super resolution factor~(SRF) of 2 was used. OBD/EFF uses one image in every iteration, therefore 30 iterations are chosen for TIP4AW. Figure~\ref{fig:CompareHirsch}(a)-(h) shows the ground truth, a typical image frame, and the estimated object for both OBD/EFF and TIP4AW.

\begin{figure}[tbh]
    \centering
    \includegraphics[width=\textwidth]{CompareHirschCombi.pdf}
    \caption{(a)~The ground truth \textit{Puppet}, (b)~a typical image frame of Puppet, (c)~estimated object (Puppet) by OBD/EFF and (d)~estimated object (Puppet) by TIP4AW. (e)~The ground truth \textit{Cameraman}, (f)~a typical image frame of Cameraman, (g)~estimated object (Cameraman) by OBD/EFF and (h)~estimated object (Cameraman) by TIP4AW. OBD/EFF parameters are: $P=Q=7$, PSF size: 13 pixels, SRF: 2. TIP4AW parameters are: $P=Q=7$, $r_{\mathbb{X}}=6$ (13 pixels), $w=35$, $\Delta w=14$, $\epsilon=10^{-4.4}$, $p_s=1.5$. Quantitative results are shown in Fig.~\ref{fig:FRCCompareHirsch}.}
    \label{fig:CompareHirsch}
\end{figure}

Figure~\ref{fig:FRCCompareHirsch} shows a quantitative comparison by the FRC corresponding to the estimated object for both OBD/EFF and TIP4AW. For reference the FRC of the image frames in Fig.~\ref{fig:CompareHirsch}, and the $2\sigma$ curve are shown.

\begin{figure}[!ht]%
\centering
    \subfloat[]{{\includegraphics[width=\textwidth,valign=c]{FRCCompareHirsch186.pdf} }}%
    \vspace{0.1cm}
    \subfloat[]{{\includegraphics[width=\textwidth,valign=c]{FRCCompareHirschCAM186.pdf} }}%
    \caption{The Fourier Ring Correlation provides a quantitative comparison of the results depicted in Fig.~\ref{fig:CompareHirsch}. The quality of the estimated object is compared to the ground truth and can be expressed as the crossing of the FRC and the $2\sigma$ curve. This crossing is denoted by the largest Fourier ring, $r_{n,\text{max}}$. (a) The quality of the image frame, estimated object with OBD/EFF and TIP4AW in Fig.~\ref{fig:CompareHirsch}(b)-(d) are $r_{n,\text{max}}=4$, $7$ and $18$ respectively.
    (b) The quality of the image frame, estimated object with OBD/EFF and TIP4AW in Fig.~\ref{fig:CompareHirsch}(f)-(h) are $r_{n,\text{max}}=4$, $5$ and $18$ respectively.
    }%
    \label{fig:FRCCompareHirsch}%
\end{figure}

\subsection{Effect of additive noise}
In the experiments with TIP4AW and OBD/EFF in this section, Gaussian white noise was added to the images. The pixels in the images have a range from $0$ to $1$ and the Gaussian noise has varying standard deviation. For every noise level the experiments were repeated 50 times. The quality of the estimated object is measured with the Structural Similarity Index (SSIM). SSIM is a metric inspired be the human visual system~\cite{SSIM}. SSIM is a method designed to measure the visual similarity of two images. Here SSIM is used to measure the influence of additive noise. SSIM is a measurement between $0$ and $1$, where values close to $1$ stand for high similarity.
\begin{equation}
SSIM = \frac{ (2\mu_o\mu_{\hat{o}}+c_1) (2\sigma_{o\hat{o}}+c_2) }
{ (\mu_o^2+\mu_{\hat{o}}^2+c_1) (\sigma_{o}^2+\sigma_{\hat{o}}^2+c_2) },
\end{equation}
where $\mu_o$, $\mu_{\hat{o}}$, $\sigma_o$, $\sigma_{\hat{o}}$ and $\sigma_{o\hat{o}}$ are the local means, standard deviations and cross-covariance for the object and estimated object respectively. $c_1$ and $c_2$ are constants.

Figure~\ref{fig:NoiseBox} shows the image quality of the estimated object for different noise levels. The standard deviations of the noise was varying form $\sigma=10^{-5}$ to $\sigma=10^{-1}$.

\begin{figure}[tbh]
    \centering
    \includegraphics[width=\textwidth]{NoiseBox168.pdf}
    \caption{Image quality of the estimated object with TIP4AW and OBD/EFF. White Gaussian noise with different standard deviations (st. dev.) was added to the set of images. The image quality is measured with SSIM. Qualitative results are shown in Fig.~\ref{fig:Noise}.}
    \label{fig:NoiseBox}
\end{figure}

The estimated objects with TIP4AW for a set of input images with different standard deviation are shown in Fig.~\ref{fig:Noise}. The images from left to right are the result of TIP4AW with input images with addition white Gaussian noise with increasing standard deviation.

\begin{figure}[tbh]
    \centering
    \includegraphics[width=\textwidth]{Noise.pdf}
    \caption{In an experiment TIP4AW is used to estimate the original object from sets of input images with different noise levels. White Gaussian noise with varying standard deviations, $\sigma=10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}$, is added to the images. The estimated object corresponding to the sets of input images with different noise levels are depicted in this figure respectively. The quality in SSIM of the estimated objects can be found in Fig.~\ref{fig:NoiseBox}.}
    \label{fig:Noise}
\end{figure}

\section{Discussion and conclusion}
\label{::Discussion}
The term deconvolution refers the inverse process of convolution. Since convolution for image formulation can only be assumed for isoplanatic aberrations, the term anisoplanatic deconvolution is abused in literature when anisoplanatic aberrations are considered. Yet, this paper refers to deconvolution. In this paper weighted multi-frame deconvolution is used to deconvolve multiple frames in order to estimate the object from blurred images. If local isoplanatism in images with anisoplanatic aberrations is sufficient, the performance of TIP4AW can reach what was possible with a batch of images with only isoplanatic aberrations.

TIP4AW was compared to OBD/EFF~\cite{Hirsch2010}. Two data sets, \textit{Puppet} and \textit{Cameraman}, were used and both algorithms were used to calculate the underlying object. Figure.~\ref{fig:CompareHirsch} shows the qualitative results of the experiments and Fig.~\ref{fig:FRCCompareHirsch} show that the estimated object with TIP4AW is of higher quality than OBD/EFF for these two data sets and in terms of FRC. Especially in the Fourier rings corresponding to low spatial frequencies TIP4AW is superior. It is hypothesised that online algorithms that process one image per iteration are not capable of handling large morph in images. TIP4AW is specifically designed for situations where large morph in the images occurs. For this the estimated object with TIP4AW will show a good performance for low spatial frequencies. For the best performance of TIP4AW it is recommended to obtain as much images of the same object as possible. It is not required to align the images, as the algorithm can correct for translations in the object by itself. Obtaining more images will result in a better image quality because it is more likely for isoplanatic subsections to occur and to measure more frequencies within the bandlimit of the system. As demonstrated in~\cite{Wilding2017} and~\cite{Ketterij2019}, TIP-based algorithm have limited performance when the object is sparse, for example in stellar imaging.

In sparse parts of the object the frequency content is not sufficient for accurate estimation of the local PSF. Consequently ringing occurs at the edge of a restored image. In Fig.~2(h) the effect of ringing is noticeable along and normal to the edges of the image. Since the local PSF depends on the image and the object as in Eq.~\eqref{eq:PSFest}, and due to the anisoplanatic nature of the aberration, it is that the effect of ringing occurs non-uniformly. A possible solution for this problem is the calculation of local PSFs at the edges of the images or the use of edge enhancement software.

Besides ringing, shadowing occurs in the restored images. Comparing Fig.~2(d) with ground truth, Fig.~2(a), it can be seen that a shadow has appeared left of the puppet's eye. Also, a shadow occurs at the left of the cameraman's coat, Fig.~2(h), that was not visible in the ground truth, Fig.~2(e). The shadow effect occurs when the estimated local PSF is larger than the adaptive support area. The size of the adaptive support area is a trade-off between accuracy of the estimated local PSF and the convergence of the algorithm.

TIP is known for its robustness to noise~\cite{Wilding2017} and the same holds for TIP4AW. In Figs.~\ref{fig:NoiseBox} and~\ref{fig:Noise} is demonstrated that both TIP4AW and OBD/EFF are robust to noise. The quality of the estimated object after 50 repetitions of both algorithms for different noise levels was measured by the SSIM. As can be seen in Fig.~\ref{fig:NoiseBox}, for low noise levels ($10^{-5}$) the quality of the estimated object with TIP4AW is higher than OBD/EFF in terms of SSIM. As the noise standard deviation increases the quality of the estimated object decreases. Figure~\ref{fig:NoiseBox} shows that the quality with TIP4AW decreases more with increasing noise than OBD/EFF. Therefore, it can be concluded that OBD/EFF is more robust to noise than TIP4AW. Nevertheless, even in noisy situations the quality of the estimated object with TIP4AW will be higher in terms of FRC.

The performance of TIP4AW does not only depend on the number of images and the level of noise. Also the parameters of the algorithm have effect. In general the algorithm requires some iterations to converge to a solution, typically 20 to 30 iterations is sufficient. Further, the user should have a guess about the size of the PSF in order to set the size of the PSF support. If the user has no information or insight in the size of the PSF the algorithm can be run multiple times for different support constraints. The size of the Gaussian apodization kernel is such that its full width at half maximum is about 2 to 4 times larger than the PSF support. The apodization size should be as small as possible while there is still sufficient information for the local PSF to be estimated, therefore the width of the Gaussian apodization kernel depends on the object and aberration. The parameter $\epsilon$ is related to the SNR and is usually small. Moreover, the quality of the estimated object will increase if more subsections are used. The computational effort of TIP4AW is mainly determined by the number of subsections, $PQ$. Figure~\ref{fig:PerformanceVsCost} shows the relation between the performance of TIP4AW measured by the crossing of the FRC and the $2\sigma$ curve and the computational efficiency, $1/(PQ)$.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{PerformanceVsCost1862.pdf}
    \caption{The performance of TIP4AW given by the crossing of the FRC with the $2\sigma$ curve increases with the the number of subsections, $PQ$. In this example the horizontal and vertical number of subsections is equal, $P=Q$. The computational efficiency is given by $1/(PQ)$ and decreases with an increasing number of subsections. As expected the performance of TIP4AW increases as the number of subsections increases.}
    \label{fig:PerformanceVsCost}
\end{figure}

Furthermore, an influence on the quality of the estimated object is the initial aberration in the recorded images. The anisoplanatic effects will always remain in the corrected image. This also holds for other image processing algorithms, unless a correction is made for every pixel independently. It is the strength of TIP4AW, that when some (nearly) isoplanatic regions occurs, this desired region will have a larger contribution to the final result. Therefore, TIP4AW actively suppresses the anisoplanatic effects.

In this paper was chosen for an approach where the least amount of prior knowledge of image formation was assumed. Note that TIP4AW works exclusively with local PSFs and not with information about the unknown phase aberration itself. Only two \textit{a priori} are assumed for the local PSFs being non-negative and real-valued, and limited in size. Both assumptions are applied in the form of a spatial domain constraint to the PSF. The advantage of this approach is that the shape of the PSF is not limited to a parameterized model. The object has a non-negative and real-valued constraint.

This paper demonstrates the potential of weighted multi-frame deconvolution. This novel approach to deconvolution was performed in conjunction with PSF variation or isoplanatism. The calculation of the weights $a_{p,q|s}$ as conducted in this paper is not optimal. Further research should be done about optimal use of weighted multi-frame deconvolution. A possible use for weighted multi-frame deconvolution is known as lucky imaging~\cite{Fried78}. When not the isoplanatism but the size of the PSF is measured, for example by the second moment of the estimated PSF, it is possible to use TIP4AW as a lucky imaging algorithm. In lucky imaging non degraded regions in the image are used for reconstruction of the object.

Another approach in this research could make explicit use of the phase aberrations. Assuming a single conjugation plane outside the pupil plane leads to a phase retrieval problem. Explicit use of the anisoplanatic phase aberration might lead to more coherence between the local PSFs.

To conclude, TIP4AW was presented as an algorithm that can handle with anisoplanatic aberrations and is robust to noise and morph. This algorithm requires multiple images of the same object and only a guess about the size of the PSF. TIP4AW is a practical algorithm for situations where space-temporal variations of the aberration is present. In general TIP4AW is an algorithm that works in batch mode, where all images are processed at the same time. Most of the computational time is spend during the PSF step in the Fourier Domain that is step (3) in algorithm~\ref{alg1}. Most of this computational time can be reduced by calculating all of the local PSFs in parallel. Besides, TIP4AW can also be implemented as an online algorithm. For online computation of TIP4AW, $S'$ images should be processed at the same time, where $S'<S$. In the first iteration images $\{i_1,...,i_{s'}\}$ are processed, in the second iteration the images $\{i_2,...,i_{s'+1}\}$, and in iteration $k$ the images $\{i_k,...,i_{s'+k-1}\}$ are processed. The use of a moving window of $S'$ images allows online implementation of TIP4AW.

\section*{Funding}
This project has received funding from the ECSEL Joint Undertaking (JU) under grant agreement No. 826589.
The JU receives support from the European Union's Horizon 2020 research and innovation programme and Netherlands, Belgium, Germany, France, Italy, Austria, Hungary, Romania, Sweden and Israel.

\section*{Disclosures}

The authors declare no conflicts of interest.