\section{Introduction}
Several advances in scattering techniques inspired the interest in neutron science, which provides different applications in areas such as chemistry, physics, biology, medicine and engineering research. These advances are strongly dependent on detector developments, which has also grown in the last years, along with the capability to obtain bright monochromatic neutron sources, whether in large or small facilities \cite{PIETROPAOLO20201, HALLWILTON2014, Ledoux2021, DIGIOVINE2021}.

The $^{3}$He  proportional counter is one of the most widely used type of neutron detectors, still widely used. Beyond its historical tradition, it presents several advantages, such as the large neutron capture cross-section ($\sigma =$ \SI{5325}{\barn} for \SI{25}{\milli\electronvolt} \cite{tendl}), low sensitivity to gamma-rays, and the possibility of working on a wide neutron energy range in the low energy region \cite{PIETROPAOLO20201}. However, due to a shortage of this gas since the early 2000s \cite{He3shortage}, the scientific community keeps searching for alternative solutions with $\ce{^{3}_{}He}$-free detectors \cite{Sacchetti2015}, employing \mbox{$^{113}$Cd ,} $^{157}$Gd , $^{10}$B\mbox{}, and $^{6}$Li  \cite{MIYAKE2011, PIETROPAOLO2013, Schooneveld_2016, COMPAAN2020, ALBANI2020} as neutron converters. Nevertheless, substituting $^{3}$He  has the price of reducing the neutron detection efficiency or increasing the detection complexity.

One versatile alternative for neutron detection consists of using one of the referred neutron converter isotopes with the Gas Electron Multiplier (GEM) \cite{Sauli1997}, a microstructure widely used to detect charged particles that can cover large sensitive areas,  present fair energy and position resolutions as well as robustness. These aspects consolidated these microstructures, which are already in use in several applications, such as high energy physics \cite{ALICE, COMPASS, ABBANEO2017298}, muon tomography \cite{Muon_tomography}, X-rays fluorescence imaging \cite{geovoXRF}, and, more recently, neutron detection \cite{KLEIN2011, KOHLI2016, ZHOU2020, ZHOU2021}.

This work presents a double-GEM detector which uses $^{10}$B$_4$C  as a neutron converter, deposited on its aluminum cathode. The ionizing products generated by the $^{10}$B(n,$\alpha$)$^{7}$Li  capture reaction, shown in Eq.~\ref{eq:decays}, are a critical factor in detector's design. Since \mbox{$^{10}$B } has a relatively large neutron capture cross section of \SI{3870}{b} for \SI{25}{\milli\electronvolt} \cite{tendl} and about \SI{19}{\percent} natural abundance \cite{elements-abundance}, one of its compounds, $^{10}$B$_4$C, is particularly interesting for our purposes, given that it also presents good wear resistance besides chemical and thermal stability \cite{ESSref3}.
\begin{equation}
    \text{n} + \ce{^{10}_{}B} \longrightarrow
    \begin{cases}
    \ce{^{4}_{}He} + \ce{^{7}_{}Li}: \quad \text{ 2.792 MeV (ground state, about 6\%)}\\
    \ce{^{4}_{}He} + \ce{^{7}_{}Li^{*}}: \quad \text{2.310 MeV (excited state, about 94\%).}
\end{cases}
\label{eq:decays}
\end{equation}
This prototype can be easily assembled and offers the possibility to cover large areas, an essential requirement for some applications. The detector was tested in the IEA-R1 nuclear research reactor \cite{PERROTTA1998}, at the Nuclear and Energy Research Institute (IPEN) in S\~ao Paulo, Brazil, using a monochromatic thermal neutron beam. Its versatility resulted from quickly obtaining neutron position sensitivity using a few electronics channels and presenting gain stability.

\section{Experimental setup}
\label{sec:experimental}
Our double-GEM detector is composed of a stack with a \SI{0.5}{\milli\meter} thick aluminum cathode coated with amorphous enriched boron carbide ($^{10}$B$_4$C ), followed by two GEM foils, provided by CERN.

The deposition of the $^{10}$B$_4$C   film on the Al cathode was done by direct-current magnetron sputtering in ESS Linköping Detector Coatings Workshop, located at Linköping, Sweden. The facility is equipped with an industrial deposition unit (CC800/9, CemeCon AG) and has been dedicated to the development and production of B$_4$C coatings as the solid converting layers in neutron detectors. Details about the deposition machine and characterization of the B$_4$C thin films for neutron detectors can be found in \cite{ESSref1,ESSref2,ESSref3}.

In this work, two $^{10}$B$_4$C  targets ($^{10}$B\mbox{ }enrichment $>$ 95 wt.\%) were installed on sputtering cathodes in the deposition machine. The Al plate was attached to another Al backing plate to achieve $^{10}$B$_4$C  coating on only one surface, where the whole assembly was mounted to an electrically grounded sample table (see Fig. \ref{fig:ESS}). The system was first pumped down to a primary vacuum pressure of $\sim$\SI{1.5e-4}{Pa} ($\sim$\SI{1.1e-6}{Torr}) while the Al plate was heated up to $\sim$\SI{280}{\degreeCelsius} for \SI{3}{\hour} mainly to remove the water residuals on the surfaces. After an additional \SI{2}{\hour} cooling, the plate was then treated with radio frequency plasma etching in $\sim$\SI{0.35}{Pa} of pure Ar to improve the coating’s adhesion by removing native oxides on the substrate surface and increasing surface roughness. Thereafter, the system was filled up with $\sim$\SI{0.8}{Pa} ($\sim$\SI{6}{mTorr}) of Ar before the beginning of the $^{10}$B$_4$C  deposition. During the deposition, the sample table was continuously rotated to obtain a better thickness uniformity. A polished Si reference was mounted at an equivalent position of the Al plate in the same deposition run for thickness measurement of the coating and the thickness of $^{10}$B$_4$C  coating was measured to be $\approx$\SI{2.2}{\micro\metre} by a profilometer (DektakXT, Bruker).

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{Figures/Figure_ESS.pdf}
    \caption{The assembly of the Al cathode plate (highlighted with the dotted box) on the Al backing plate for deposition.}
    \label{fig:ESS}
\end{figure}

In addition to the deposition of the converter layer, the geometrical disposition of the elements within the detector is a key factor regarding the signal we want to measure. The drift, transfer, and induction regions of our prototype are represented in Fig. \ref{fig:scheme}. The drift region has \SI{2}{\milli\meter} thickness and is biased with \SI{100}{V}, while the transfer and induction regions have both \SI{1}{\milli\meter} thickness and are biased with \SI{300}{V} and \SI{400}{V}, respectively. The detector works under an open flow of an Ar/CO$_2$ (90/10) gas mixture, with a flow of \SI{6}{L/h} through the $\sim$\SI{0.85}{\liter} volume aluminum chamber. The readout plane, produced at CERN, consists of \num{256}+\num{256} strips in X and Y directions, covering a \SI{100}{\centi\meter\squared} sensitive area.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.7\textwidth]{Figures/ndetector_electronics6.pdf}
    \caption{Detector scheme, showing the drift, transfer and induction regions. Not to scale. Figure adpated from \cite{LucasMaster}.}
    \label{fig:scheme}
\end{figure}

The thickness of the drift region determines the sensitivity of the detector to gamma rays background. In this context, a common requirement for thermal neutron detectors is to be insensitive to the gamma flux created by radiation that escapes the reactor core. To evaluate the reduction of such background, two thicknesses were tested for the drift region: \SI{20}{mm} and \SI{2}{mm}. As the reduction of the gamma radiation background was considerable for the \SI{2}{mm}, we operated the detector using this drift region configuration. This is possible because the reaction products' high ionization power generates a large amount of primary electrons. That way, even a small segment of their tracks is enough to provide a large signal.

The signal multiplied by the GEM stack is collected by \num{256}+\num{256} strips, which are interconnected by a resistive chain, as shown in Fig. \ref{fig:rchain_scheme}, composed of \SI{60.00(6)}{\ohm{}} SMD resistors. There are five signals collected to reconstruct a neutron image: two signals for the X direction, two for the Y direction and the charge signal, collected from the bottom electrode of the last GEM (see Fig. \ref{fig:scheme}). The charge signal is used as the trigger for the data acquisition, as shown in Fig. \ref{fig:eletronic_scheme}, since it is position-independent, coming directly from the multiplication. Considering the expected extraction efficiency $\epsilon_{\text{ext}}$ for the last GEM \cite{geovo}, the charge signal collected at its bottom is given by $1-\epsilon_{\text{ext}}$, which is $\sim$ \SI{64}{\percent} of the multiplied charge for our configuration.

\begin{figure}[H]
    \centering
    \includegraphics[height=0.5\textwidth]{Figures/resistivechain_GEovane.pdf}
    \caption{PCB resistive chain to fit the detector signal output. Each board has 128 SMD resistors and a PANASONIC\textregistered 130-pin connector. Since the readout has 256 strips for each direction, two plates are used for each direction. Figure from \cite{geovo}.}
    \label{fig:rchain_scheme}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{Figures/electronics_model15.pdf}
    \caption{Signal acquisition setup scheme. The trigger circuit is pictured in blue. Figure from \cite{LucasMaster}.}
    \label{fig:eletronic_scheme}
\end{figure}

In order to reconstruct the event's position, we use Eq.~\ref{eq:charge_dist} \cite{geovoXRF}, where $Q$ is the charge information collected from the bottom of the last GEM, and $L_x$ and $L_y$ are the lengths of the sensitive area for each direction. As indicated in Fig. \ref{fig:eletronic_scheme}, $X_1$ and $X_2$ are the signals collected at each end of the resistive chain in the X direction, while $Y_1$ and $Y_2$ are the equivalent signals in the Y direction.
\begin{equation}
\label{eq:charge_dist}
    X_{*} = L_x\frac{X_2 - X_1}{Q} \qquad Y_{*} = L_y\frac{Y_2 - Y_1}{Q}
\end{equation}

As mentioned before, the detector was tested at the IEA-R1 nuclear reactor using the same neutron beam of the AURORA diffractometer \cite{PARENTE2010678} (see Fig. \ref{fig:beam_detector}). This monochromatic neutron beam is selected out of a collimated radial neutron beam from the core of the reactor, which operates at \SI{4.5}{\mega\watt}. A multi-blade focusing silicon monochromator is used to select these neutrons with a fixed wavelength of \SI{1.399}{\angstrom} (\SI{41.8}{\milli \electronvolt}) providing a \SI{6.22(19)e4}{n\, \centi\metre^{-2}\second^{-1}} flux ($\phi$), measured by the Nuclear Metrology Laboratory at IPEN with the $^{197}\mathrm{Au}(n,\gamma)^{198}\mathrm{Au}$ reaction, using a \SI{12}{\milli\metre} diameter thin gold foil, centered at the beam outlet (letter E of Fig. \ref{fig:beam_detector}).

The detector (letter A of Fig. \ref{fig:beam_detector}) and the preamplifiers' box (letter C of Fig. \ref{fig:beam_detector}) were mounted over an aluminum plate. This setup was fixed to a X-Y-Z precision table, that allowed the fine tuning of the detector's position with respect to the beam.

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{Figures/BEAM_compressed.pdf}
    \caption{Position of the detector (A) relatively to the beam. The following elements are also indicated: ground shielding of the y-direction resistive chains (B), preamplifiers box (C), high voltage inputs (D), and the beam outlet (E).}
    \label{fig:beam_detector}
\end{figure}

\section{Results and discussion}
In this section, we present the main experimental results concerning the characterization of the detector: the spatial resolution obtained by two different methods, the detector efficiency, and the gain stability. We also discuss how these results relate to the previous expectations obtained with simulations.

\subsection{Position calibration and resolution}
To perform the spatial calibration of the detector, we used handmade cadmium masks which has a high neutron absorption cross section. Holes of different diameters were drilled in three masks. The positions and areas of the holes were measured using an image obtained with a table scanner. One of the produced masks in shown in Fig. \ref{fig:calibration}, with its \num{15} holes, and its neutron image, obtained with the mask centered with respect to the detector's sensitive area.

\begin{figure}[H]
\begin{subfigure}[T]{.5\textwidth}
  \centering
  \includegraphics[width=0.8\linewidth]{Figures/Mascara_furos_03new_compressed.pdf}
  \label{fig:calib1}
\end{subfigure}%
\begin{subfigure}[T]{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{Figures/MASKmatrix2.pdf}
  \label{fig:calib2}
\end{subfigure}
\caption{Left: digital image of the mask obtained with a table scanner. Right: calibrated neutron image.}
\label{fig:calibration}
\end{figure}

The calibration is given by a linear dependency between the coordinates of certain holes in arbitrary units and their coordinates in known distances. The chosen holes are colinear in the X and Y directions, normally the ones in the central column for Y direction and the central row for X direction. We measured the detector spatial resolution by two methods: edge and point spread functions. The first one consists of using edges to fit step functions, as shown in Fig. \ref{fig:edges}.

The data was fitted with the complementary error function
\begin{equation}
    \erfc(p) = 1 - \frac{2}{\sqrt{\pi}} \int_{0}^{p} e^{-t^{2}}dt,
    \label{eq:erf_function}
\end{equation}
where
\begin{equation}
    p(x) = \frac{x-\mu}{\sqrt{2}\sigma}.
    \label{eq:PvarErf_function}
\end{equation}

This function provides the same parameters of interest ($\sigma$ and $\mu$) as the edge spread function (ESF). In this work, we define the position resolution as the FWHM of the Gaussian (or the 10-to-90\,\% distance of the ESF), i. e., 2.35$\sigma$.

\begin{figure}[H]
    \begin{subfigure}{.5\textwidth}
      \centering
      \includegraphics[width=1\linewidth]{Figures/mask_edges.pdf}
    \end{subfigure}
    \begin{subfigure}{.5\textwidth}
      \centering
      \includegraphics[width=1\linewidth]{Figures/edge_fit.pdf}
    \end{subfigure}
\caption{Left: neutron image of one of the masks, with a selected region of an edge marked as the red rectangle. Right: Complementary error function ($\erfc$) adjusted in the y-projection of the chosen area, the spatial resolution obtained is \SI{3.0(1)}{\milli\metre}. Figure adapted from \cite{LucasMaster}.}
\label{fig:edges}
\end{figure}

We measured the position resolution with the edge spread function in different regions of the detector's sensitive area and along X and Y coordinates. The results are shown in Fig. \ref{fig:edgesmap}.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.55\textwidth]{Figures/edges6.pdf}
    \caption{The points where the spatial resolution was measured using the edge spread function. The set of spatial resolution values resulted in an average value of \SI{2.8(3)}{\milli\metre} for the X direction and \SI{3.0(1)}{\milli\metre} for the Y direction, compatible with each other, where the uncertainties were given by the standard deviation of each set. Consequently, it is reasonable to consider the detector's resolution as the average of the resolutions for both directions, which results \SI{2.9(3)}{\milli\metre}.}
    \label{fig:edgesmap}
\end{figure}

The second method relies on the fact that the FWHM of an ideally punctual hole converges to the spatial resolution of the detector, which is the width of the point spread function (PSF). The projection $P_{meas}$ obtained for a hole with a given diameter is the convolution of the PSF, for which it is safe to assume a 2D Gaussian function, with the projection of a circular area. This projection over the X direction can be written as
\begin{equation}
    P_{meas}(x) = (P_x \circledast g)(x):=\int_{-\infty}^{\infty} f(\tau) g(x-\tau,\sigma)d\tau,
    \label{eq:convolution}
\end{equation}
where $g$ is the Gaussian function and $P_x$ the projection of a circular area over the X direction. Fitting this convolution for holes of several diameters allows us to estimate the spatial resolution. The convolution shown in Eq.~\ref{eq:convolution} is computed numerically and the fit\footnote{the fitted parameter is the resolution itself, since it is linearly dependent on $\sigma$.} is done using the Least Square Method. The result is shown in Fig. \ref{fig:convFits}.

\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{Figures/Ajuste_resolucaoX_noChi2.pdf}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{Figures/Ajuste_resolucaoY_noChi2.pdf}
\end{subfigure}
\caption{Left: Convolution fit for holes with different diameters, resulting in the spatial resolution in the X (left) and Y (right) directions. Figures adapted from \cite{LucasMaster}.}
\label{fig:convFits}
\end{figure}

Although compatible with the \SI{3}{\milli\metre} obtained with the ESF, the values obtained with this method are slightly better. This can be explained by the fact that the FWHM assigned to each diameter is the average of the FWHM for several runs using holes of that diameter.

The values obtained with the two different techniques result in a spatial resolution compatible with at least \SI{3}{\milli\metre}.

\subsection{Detector efficiency and gain stability}
\subsubsection{Analytic Calculation}
It is possible to obtain the expected detection efficiency $\varepsilon_d$ from analytical calculations. The first step is to compute the capture efficiency $\varepsilon_c$, which depends on the macroscopic cross-section $\Sigma$ for the converter material and its thickness $d =$ \SI{2.2}{\micro\metre}. Taking the converter $(^{10}\mathrm{B}_4\mathrm{C})$ density $\rho=$  \SI{2.3}{\gram\per\centi\meter\cubed} \cite{Hoglund_depo}, its neutron capture cross section for \SI{41.8}{\milli\electronvolt} neutrons \mbox{$\sigma =$ \SI{2994(9)}{b}} \cite{tendl}, its molar mass $\mathrm{M}_{\mathrm{B}_4\mathrm{C}} = $ \SI{52}{\gram\per\mol}, and
the stoichiometry $a$ of $^{10}$B in a $^{10}$B$_4$C  converter layer, which is \num{4}, one can obtain the number of target nuclei of the converter:
\begin{equation}
    N = a\frac{\rho}{\text{M}}N_A \approx 1.065 \times 10^{23} \, \text{cm}^{-3}.
    \label{eq:N_target_nuclei}
\end{equation}

The result of Eq.~\ref{eq:N_target_nuclei} is used to compute the   macroscopic cross-section \cite{HUSSEIN20071} for the converter layer:
\begin{equation}
    \Sigma = N\sigma = 3.19(1)\times10^{-2}\, \mu\text{m}^{-1}.
\end{equation}

The neutron beam intensity after crossing our converter layer is
\begin{equation}
    I_{\text{beam}} = e^{-\Sigma d}\approx 93.23(12)\%.
    \label{eq:beamintensity}
\end{equation}
Therefore, the capture efficiency is given by the percentage of absorbed neutrons:
\begin{equation}
    \varepsilon_c = 1 - e^{-\Sigma d}\approx 6.77(58)\%.
    \label{eq:conversionefficiency}
\end{equation}

When the neutron is captured too far from the lower side of the converter layer (at the boundary with the gaseous region), it might not be detected because the products of the capture reaction do not have enough energy to exit the layer, as pictured in Fig. \ref{fig:layer_geometricC} (case A). Therefore, such products have to fulfill a set of geometrical requisites that depend on the capture's position, which are shown in Fig. \ref{fig:layer_geometricC} (case B and C). Finally, the products that enter the drift region with too low energy are also not detected since there is an energy threshold to handle electrical noise and other acquisition system specifications.

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{Figures/layer_angle3-jinstNEW.pdf}
    \caption{The spheres of radius $\tau_{\alpha}$ and $\tau_{Li}$ define the accessible region for the $\alpha$ and $^{6}$Li  nuclei, respectively, inside the converter layer (A). Given certain capture point, only the products within the solid angle $\Omega$ reaches the gas (B). The solid angle $\Omega$ is a function of the angle $\text{cos}^{-1}(x/\tau)$, where $x$ is the smallest distance between the capture point and the drift region (C). Figure adapted from \cite{LucasMaster}.}
    \label{fig:layer_geometricC}
\end{figure}

In order to account for the constrains discussed above, the values for $\tau_{\alpha}$ and $\tau_{\text{Li}}$, for both possible decay channels of Eq.~\ref{eq:decays}, were evaluated from direct integration of the average dE/dx plot generated through SRIM \cite{ZIEGLER20101818}.  It is also necessary to consider the energy threshold, which was done by finding $\tau$ for which the average remaining energy for the given particle is equal to the threshold used in the experiment, of \SI{100(10)}{\kilo\electronvolt}. This process is shown in Fig. \ref{fig:Eremainingplot}.

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{Figures/plotalldepintR3new_v5.pdf}
    \caption{Remaining energy for each $^{10}$B(n,$\alpha$)$^{6}$Li  reaction product as a function of depth for a theoretically infinite $^{10}$B$_4$C  converter layer. The $\tau$ values used for the calculations are the depth at which the particle's remaining energy equals the threshold. This condition gives: $\tau_{\alpha}$(\SI{1775}{\kilo\electronvolt}) $=$ \SI{3.97}{\micro\metre}, $\tau_{\alpha}$(\SI{1470}{\kilo\electronvolt}) $=$ \SI{3.13}{\micro\metre}, $\tau_{\mathrm{Li}}$(\SI{1015}{\kilo\electronvolt}) $=$ \SI{1.54}{\micro\metre} and $\tau_{\mathrm{Li}}$(\SI{840}{\kilo\electronvolt}) $=$ \SI{1.31}{\micro\metre}.}
    \label{fig:Eremainingplot}
\end{figure}

The considerations above are translated to the calculation for the detection efficiency $\varepsilon_d$ as follows:

\begin{equation}
    \varepsilon_d = F\varepsilon_c = F(1 - e^{-\Sigma d}),
    \label{eq:deteff1}
\end{equation}
where $F$ is the share of capture products emitted within the solid angle $\Omega$ of Fig. \ref{fig:layer_geometricC}B, and is computed as

\begin{equation}
    F = \sum_{\text{all products}}F(d,\tau_i) = \sum_{\text{all products}}\int_{0}^{d}\int_{0}^{\text{cos}^{-1}(x/\tau_i)}d\Omega dx \approx 0.477,
    \label{eq:detectablefraction}
\end{equation}
for a boron carbide layer of thickness \SI{2.2}{\micro\metre}. As result, we obtain the neutron detection efficiency
\begin{equation}
    \varepsilon_d = 3.22(18)\,\%.
    \label{eq:deteff_final}
\end{equation}

\subsubsection{Simulations}
Monte Carlo simulations play an important role in the development and optimization of this kind of detector. For this purpose, GEANT4 is a powerful tool to describe the detector and for the simulation of the interaction of particles with matter \cite{Agostinelli2003}, allowing the simulation of the nuclear reaction, like Eq.~\ref{eq:decays}, as well as the transport of the reaction products inside the detector.

The detector was simulated with GEANT4 using the physics list \textit{QGSP\_BERT\_HP} (which has high precision models for low energy neutrons \cite{geant4_physicsReference}), according to the device described in section \ref{sec:experimental}. The simulation consisted in a parallel uniform beam with \num{e6} neutrons of \SI{1.399}{\angstrom} ($\approx41.8\, \mathrm{meV}$), hitting the detector in a direction perpendicular to the converter layer.

Taking 14 different thicknesses ranging from $\SI{0.25}{\micro\meter}$ to $\SI{5}{\micro\meter}$ of the boron carbide converter layer, we evaluate the thermal neutrons detection efficiency by counting the charged particles that leave the boron converter and reach the gaseous volume. For a nominal thickness of $\SI{2.2}{\micro\meter}$, an efficiency of $\SI{3.14(10)}{\percent}$ was obtained. While the maximum efficiency was achieved with a thickness of $\SI{2.5}{\micro\meter}$ resulting in $\SI{3.17(10)} {\percent}$, as shown in Fig. \ref{fig:sim_eff}.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.575\textwidth]{Figures/efficiency_mona3.pdf}
    \caption{Detection efficiency for one layer of boron carbide, simulated with GEANT4. The fit for maximum efficiency is an arbitrary function adjusted just to guide the eye within the interval that comprises the data around the peak (continuous line).}
    \label{fig:sim_eff}
\end{figure}

Using the low-energy electromagnetic model Livermore \cite{geant4_physicsReference}, we simulated the deposited energy by the reaction products inside the $\SI{2}{\milli\meter}$ of the drift region presented in Fig. \ref{fig:sim_spectra} together with the experimental data. This simulation does not include the energy threshold of our system nor the charge multiplication process, i.e., the acquired energy is equal to the energy deposited by the reaction products in the absorption region of the detector.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.7\textwidth]{Figures/h3new4v3.pdf}
    \caption{Experimental energy spectrum compared with the GEANT4 simulation. The simulation does not include the effect of the detector finite energy resolution and was normalized to the maximum value of the experimental energy distribution. Figure from \cite{LucasMaster}.}
    \label{fig:sim_spectra}
\end{figure}

\subsubsection{Experimental measurements}
The detection efficiency is evaluated by dividing the number of detected neutrons $N_{det}$ by the total number of neutrons $N_{tot}$ hitting the detector. We calculated the total number of neutrons from the neutron flux $\phi$, measured as described in section \ref{sec:experimental}, the area $A$ selected from the beam and the time length $\Delta t$ of the measurement. The neutron flux was corrected using the counts of a fission chamber installed at the neutron beam outlet using as reference the number of counts registered by this same chamber during the neutron flux calibration process. With these values, the experimental neutron detection efficiency is given by

\begin{equation}
\label{eq:deteff}
    \varepsilon = \frac{N_{det}}{N_{tot}} = \frac{N_{det}}{{\phi A \Delta t}}.
\end{equation}

Using cadmium masks with squared holes, we selected different areas from the beam: \SI{0.36}{\centi\metre\squared} and \SI{1.44}{\centi\metre\squared}, at the same position where the flux measurement was taken. The total number of detected neutrons $N_{det}$ was counted with negligible uncertainties at these rates. We obtained \SI{2.75(26)}{\percent} detection efficiency for the \SI{0.36}{\centi\metre\squared} mask and \SI{2.57(15)}{\percent} for the \SI{1.44}{\centi\metre\squared}, which are compatible within $1\,\sigma$ with each other.

The gain stability was monitored in real-time, in order to provide necessary corrections to the image reconstruction in offline data analysis. A gain shift would imply changes in efficiency since the energy spectrum would be displaced from its original position with respect to the energy threshold. Besides that, changes in the gain over time may bring the detector to an unstable operation in terms of discharge probability.

The data shown in Fig. \ref{fig:TimeEvo} present a method to analyze the gain stability: the energy channel of the run was divided into batches of \SI{4}{\minute} each (separated by the dashed vertical lines) and plotted over time. Then, we locate the spectrum peak using a peak identification algorithm \cite{MORHAC2000108} and fit a Gaussian to it, registering its center as the cross marker in Fig. \ref{fig:TimeEvo}. Slow gain changes due to environmental changes are easily seen as a tendency in the sequence of the points, whose energy spectrum can be drifted by a reasonable amount as it happens in the case reported by \cite{geovoXRF}. The tendency analysis can be carried out by fitting the behaviour of the difference between each peak position and the average position of all peaks, which we labelled as $\Delta_{\text{peak}}$ in the Fig. \ref{fig:TimeEvo}.

\begin{figure}[H]
  \centering
  \includegraphics[width=1\linewidth]{Figures/plotFig15_final3.pdf}
  \caption{Temporal stability analysis with the tendency analysis. This run has 4793145 events and was divided into batches of \SI{4}{\minute} (\num{123403} entries each), generating 39 peaks. The average of the peak position was \num{897(2)} ADC channels. The fit for $\Delta_{\text{peak}}$ resulted an angular coefficient $-0.028(9)$ ADC channel/min and a linear coefficient equals to  $2.2(8)$ ADC channels.}
  \label{fig:TimeEvo}
\end{figure

There is a slight drift of the spectrum's peak within about \SI{2.5}{\hour}, as the fit result of Fig. \ref{fig:TimeEvo} suggests. However, this drift is smaller than \num{10} ADC channels, corresponding to about \SI{0.2}{\percent\per\hour} change, which is not statistically sufficient to be interpreted as a reasonable change. Also, at the beginning of the operation, some slow gain stabilization of these detectors are expected, as reported in other works \cite{CHATTERJEE2021165749, Chatterjee_2020, HAUER2020164205, ALFONSI20126, CHERNYSHOVA2020111755}. The monitoring remains essential since it can be used to correct the detector's efficiency if significant changes happen, using the detector's current rate and the neutron monitor fission chamber.

\section{Final Considerations}
A neutron detector was built with a single boron carbide layer and a resistive charge division readout. This detector shows the feasibility of a position-sensitive detector that can be easily mounted in a neutron beam line based on simple standard electronics. The detector also presented negligible sensitivity to gamma rays (with no events when the neutron beam is off), an essential feature in higher radiation background environments like nuclear reactor research facilities. If needed, the gamma sensitivity can be controlled by changing the drift region, and a new energy pre-calibration can be quickly obtained with new simulation results.

The average of the obtained thermal neutron detection efficiency is \SI{2.66(30)}{\percent}, which agrees with the prediction of \SI{3.22(18)}{\percent} and with the preliminary simulation of \SI{3.14(10)}{\percent}, and can be reduced by altering the layer thickness allowing different applications, such as beam monitors.
At the current configuration, one can use the detector to measure the beam profile given that the position resolution obtained was at least \SI{3}{mm}. There is also the possibility to raise the efficiency by adding new converter layers over the GEMs, as seen in other works \cite{MBGEM,KLEIN2011, KOHLI2016}. For future works and following that same idea, it is possible to use boron coated Thick-GEMs, operating with unitary effective gain, raising the thermal neutron detection efficiency.

\acknowledgments
This work was supported by CNPq grant 156767/2019-8.\\
H. Natal da Luz acknowledges project GAČR GA21-21801S (Czech Science Foundation).


\section*{CVPRW review}

	%

	%

	\title{An Architecture for the detection of GAN-generated Flood Images with Localization Capabilities  \vspace{-0.7cm}}
	%

	\author{\IEEEauthorblockN{Jun Wang, Omran Alamayreh, Benedetta Tondi and Mauro Barni}
		\IEEEauthorblockA{Department of Information Engineering and Mathematics, University of Siena, Via Roma, 53100 Siena, Italy\\

		}

	}

	%

	\maketitle

	\begin{abstract}

		In this paper, we address a new image forensics task, namely the detection of fake flood images generated by ClimateGAN architecture. We do so by proposing a hybrid deep learning architecture including both a detection and a localization branch, the latter being devoted to the identification of the image regions manipulated by ClimateGAN. Even if our goal is the detection of fake flood images, in fact, we found that adding a localization branch helps the network to focus on the most relevant image regions with significant improvements in terms of generalization capabilities and robustness against image processing operations. The good performance of the proposed architecture is validated on two datasets of pristine flood images downloaded from the internet and three datasets of fake flood images generated by ClimateGAN starting from a large set of diverse street images.
	\end{abstract}

	%

	\IEEEpeerreviewmaketitle

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

	\enlargethispage{\baselineskip}

	Nowadays, Generative Adversarial Networks (GANs)  can synthesise artificial images with a quality that can easily fool humans. As a consequence,  Multimedia Forensics (MF) researchers have increased their efforts to develop tools capable of detecting GAN-generated contents.
	%
	Climate change is almost invariably considered as one of the most serious threats to humankind and one of the hardest challenges our society will have to face with in the coming years.
	The visualization of the effects of climate change can play a major role to develop new weather analysis tools and to raise public awareness. For this reason, several works have been proposed to generate weather-sensitive images using GAN networks \cite{kawakami2021semi,karacan2016learning,rothmeier2021let,lin2019weagan,schmidt2021climategan}.

	Though the purpose of these works is to promote research and raise the awareness on climate change, the risk exists that climate-sensitive synthetic images are used maliciously or exploited within organised disinformation campaigns. It is, then, of paramount importance to develop tools capable of distinguishing natural and synthetic weather images.

	In this paper, we present the first method for the detection of GAN-generated weather images. Specifically, we focus on the detection of flood images generated by the \textit{ClimateGAN} network, which can be used to transform general street images into images showing extreme flood conditions.

	%
	The proposed network has a hybrid architecture based on ResNet50 \cite{he2016deep}. The basic ResNet50 architecture is used to extract the features allowing the distinction between fake and natural images. Even if our goal is to distinguish natural and fake flood images, however, the features extracted by the convolutional layers of ResNet50, are also used to drive a fully convolutional architecture in charge of localizing the manipulated image region. One of the main findings of our work, in fact, is that adding a localization task on top of the feature extraction layers, has a positive effect on detection accuracy, even if localization is not very accurate. The reason for such a counterintuitive behaviour is that by asking the network to localize the manipulation, we somehow force the feature extraction layers to focus on the most significant parts of the image.
	%

	%
	The effectiveness of the proposed method was validated on the collected flood images showing a very good detection accuracy. The performance of the localisation module
	is also generally good, even if we remark that the final goal of our architecture is the detection of GAN generated images and the addition of a localisation module is mainly an expedient to force the network to look at the most significant image regions.
	Among the strengths of the proposed solution we stress that: i) the detector can be trained with a reduced number of images (400 natural and 400 fake images in our experiments); ii) the detector exhibits very good generalisation capabilities, as demonstrated by the experiments we run on natural images gathered from various sources across the Internet; iii) the proposed method shows good robustness against common image processing methods such as image resizing, addition of Gaussian noise, JPEG compression, and filtering.

	The rest of the paper is organized as follows. In Section \ref{sec:related}, we briefly

	overview methods for the detection of synthetic images. In Section \ref{sec:data}, we describe the datasets we have used in our experiments. The proposed forgery detector is detailed in Section \ref{sec:methodology}. The results of the experiments we carried out are reported and discussed in Section \ref{sec:experiments}. In Section \ref{sec:conclusion}, we conclude the paper with some final remarks and clues for future research.

	\section{Related work}
	\label{sec:related}

	Traditionally, the pipeline for the design of an effective image forgery detector consists in
	extracting a set of hand-crafted features \cite{li2020identification,mccloskey2019detecting,hu2021exposing}, and using them to detect the presence of forgery artifacts.
	For the detection of GAN-generated images and forgeries, some works have proposed to analyze the frequency artefacts introduced by GANs due to the upsampling operation applied by such networks \cite{yu2019attributing,liu2021spatial,luo2021generalizing}.
	Recently, data-driven detectors based on various CNN architectures, e.g.,\cite{rossler2019faceforensics++, chai2020makes,zhao2021multi},
	are gaining more and more popularity, due to their superior performance, at least when the test conditions match those used during training.
	%

	%
	Forgery localization has received relatively less attention, yet the possibility to localize where an image has been tampered with would be particularly useful to explain the detection results and improve the generalization capability of the detector.
	%

Most of the appraoches dealing with the localization of GAN genereated contents focus on faces, namely, deepfake.
	A first technique focusing on the localization of manipulated image areas was proposed in \cite{songsri2019complement}, by relying on paired facial landmarks information and a face image as input. Huang et al. \cite{huang2020fakelocator} proposes a \textit{FakeLocator} model to localize the texture artefacts introduced by the up-sampling steps usually applied by GAN architectures. By focusing on a specific task, i.e., face swapping, \cite{mazaheri2022detection} and \cite{nguyen2019multi} obtain tampering localization by segmenting the input images into manipulated regions, and then share the features for the detection task.
	Similarly, \textit{Face X-ray} \cite{li2020face} analyzes the artefacts introduced when two images are blended and use this information to guide the training of the detector.
Finally, Dang et al. \cite{dang2020detection} and Zhao et al. \cite{zhao2021multi} use an attention mechanism to improve the detection accuracy and to identify the image region that mostly contributes to the detection.

	In our work, we propose a hybrid framework, wherein localization is used to aid the detection of GAN manipulated flood images

	To the best of our knowledge, this is the first work focusing on the detection of GAN manipulated flood images.

	\enlargethispage{\baselineskip}

	\section{Dataset construction}
	\label{sec:data}

	\subsection{Natural flood images}

	With regard to natural flood images, we relied on two public datasets released recently.

	%
	The \textbf{Roadway Flooding Image (RWFI)} dataset \cite{sazara2019detecting} consists of 441 images with different scenes from urban, suburban and natural settings.

	The images have all the same size equal to 385$\times$512. Mask images of flooded areas are also available. We used this dataset for training.
	%
	The \textbf{WSOC Flood Image} dataset \cite{zaffaroni2020water} contains 490 real flood images gathered from Twitter and online news using "floods" as the keyword for the search query. 439 of such images are also accompanied by a binary mask image identifying the flooded areas.
	We used these 439 flood images to test the generalization capability of the proposed detector.

	\subsection{ClimateGAN Images}
	\label{sec:climate}

	%

	We focused on the detection of synthetic flood images generated by  the \textit{ClimateGAN} network described in \cite{schmidt2021climategan}.

	\textit{ClimateGAN} is an architecture for image-to-image translation whose goal is to produce images with extreme flooding conditions from street images.

	%

	%
	By using the pre-trained \textit{ClimateGAN} model released by the authors of \cite{schmidt2021climategan}\footnote{\url{https://github.com/cc-ai/climategan}}, we built three synthetic flood image datasets

	by using diverse street images. Specifically, we first created 3,900 GAN flood images (hereafter referred to as \textbf{StreetG} dataset), starting from the Cityscapes \cite{cordts2016cityscapes} and Kitti \cite{geiger2013vision} video datasets described in \cite{Alamayreh2021gan}.

	For each video, we took a frame every 60 frames and used the central cropped area of size 512 $\times$ 512 as input of the \textit{ClimateGAN}. We also built two additional small \textit{ClimateGAN} datasets, namely the \textbf{WebG132} and \textbf{WebG504} datasets, to validate the generalization capability of the proposed detector. The street view images used to build these datasets were collected from the internet by using the Mapillary tool (\url{https://www.mapillary.com}). More specifically, WebG132 contains 132 \textit{ClimateGAN} flood images built as  \textbf{StreetG}. WebG504 is an extended version of WebG132 built by resizing the input images rather than cropping. WebG504 consists of 504 GAN flood images.
	%

	Some examples of the fake images contained in the datasets are given in Fig. \ref{fig:gan}. The ClimateGAN image datasets split in training and testing are available at the link \url{https://drive.google.com/drive/folders/12Qd1-T8J0jCwYYqiOoARsahyElB--Jai?usp=sharing}.

	\begin{figure}
		\centering
		\includegraphics[width=0.95\columnwidth]{figures/GANflood.pdf}
		\vspace{-0.3cm}
		\caption{GAN generated flood images. Top: StreetG; Middle: WebG132; Bottom: WebG504.}
		\label{fig:gan}
	\end{figure}

	\begin{figure}[t]
		\centering

		\includegraphics[width=1\linewidth]{figures/ClimateGAN_Figure4.pdf}

		\vspace{-0.3cm}
		\caption{Overview of the proposed detection/localization architecture.}
		\label{fig:framework}
		\vspace{-0.3cm}
	\end{figure}

	\enlargethispage{\baselineskip}

	\section{The Proposed Detector}
	\label{sec:methodology}

	In this section, we present our framework for the detection of GAN-manipulated flood images and the localization of the manipulated flooded areas.

	%
	Given a flood image $x$, with $x \in \mathbb{R}^{H \times W \times 3}$, the goal of our system is to associate to $x$ a label $y$ and a mask $M$, in such a way that $y =1$, if $x$ is a manipulated image and $0$ otherwise, and the mask $M$ indicates the pixels where the image has been manipulated (pixels for which $M = 1$ indicate the manipulated areas).
	The main intuition behind our work is that asking the network to localize the manipulated areas also helps the detection process, since, during training, we can instruct the network to focus on the image areas that have been most likely been manipulated, namely the water regions.

	The overall architecture of the proposed method is shown in Fig. \ref{fig:framework}.

	To start with, given an image $x$, we pass it through a backbone ResNet50 network to extract semantic features $SF$.

	%
	The features extracted in this way are exploited for the detection and localization tasks. Specifically, for the localization task, we consider a fully convolutional network (FCN) \cite{long2015fully},

	where the localization map $M_{loc}$ is obtained by using two convolutional layers and a sigmoid function layer. At the same time, the semantic features $SF$ obtained from the last convolutional layer of ResNet50 are sent to an AvgPool layer and a Fully Connected (FC) layer to distinguish between GAN and real images.
	%
	%
	In order to train the proposed architecture, we rely on two loss functions,

	one associated with the detection task and one with the localization task. The detection loss is defined as:
	%
	\begin{equation}

		L_{det}=\frac{1}{N}\sum_{i}Y_{i} \log(y_i) + (1-Y_{i})\log(1-y_i)
	\end{equation}
	%
	where the sum is extended to all the images in the training set, $y_i$ is the output of the detector in correspondence to the input image $x_i$, and $Y_i$ is the corresponding ground truth,  which determines whether or not an image is manipulated.
	In a similar way, the localization loss is defined as
	%

	\begin{equation}

		L_{loc}=\frac{1}{S}\sum_{i}^{N}\sum_{j}^{H\times W}G_{ij} \log(M_{ij}) + (1-G_{ij})\log(1-M_{ij})
	\end{equation}

	%
	with $S = N\times H\times W$, and where $M_i \in \mathbb{R}^{H \times W}$ is the estimated manipulation map for the test image $i$, $M_{ij}$ the value assumed by the map in correspondence of the $j$-th pixel, and $G_{ij}$ the corresponding ground-truth.

	The total loss used to train the model, then,  is the weighted sum of the two losses above, namely:
	%

	$L_{tot} = \lambda _{det} \cdot L_{det} + \lambda_{loc} \cdot L_{loc}$,

	%
	where $\lambda_{det}$ and $\lambda_{loc}$ must be set in such a way to balance the importance of the two loss terms.

	\section{Experiments}
	\label{sec:experiments}
	\subsection{Setting}
	\textbf{Baseline methods.}

	To prove the effectiveness of the proposed method, we selected two popular deep neural networks, namely Xception and ResNet50, which have been proven to be effective for several forgery detection applications (see, for instance,  \cite{Alamayreh2021gan} and \cite{gragnaniello2021gan}). We considered various training strategies. We trained

	an Xception and ResNet50 model as standard binary detectors, to decide if the input image has been generated by \textit{ClimateGAN} or not. In the ResNet50 case, the architecture coincides with the classification branch of the proposed hybrid architecture. Besides, we also trained four models by using an additional mask image as an input, to guide the attention of the network toward water regions. In detail,  two models (referred to with the label MUL) take as input the image $x$ multiplied by a binary mask indicating the water regions. The other two models (referred to with the label CAT), take as input a 4-channel image, with the fourth band containing a binary mask highlighting the water areas.
	%
	The performance obtained by the above networks are compared to the detection performance of our hybrid architecture, hereafter referred to as ResNet50\_hyb.
	%
	It is worth stressing that, in contrast to the networks taking as an additional input the mask with the water regions, our method does not need such information. In fact, we force the network to look at the water regions indirectly, by asking the network to localise the manipulated areas.

	\textbf{Implementation details}.

	All the models were trained with the same configuration from the scratch for 30 epochs.

	For the optimization, we used the Adam optimizer with a learning rate of 0.0001, a mini-batch size of 16. The input images were resized to 224$\times$224$\times$3 and normalized with mean [0.485, 0.456, 0.406] and variance [0.229, 0.224, 0.225]. The training images were augmented by [saturation, brightness, contrast] adjustment with factors equal to [0.05, 0.05, 0.05].
	%
	We trained the networks by using the RWFI dataset for real images and an equal number of GAN-generated flood images from the StreetG dataset. The training datasets were split into training and validation subsets with percentages 80\% and 20\% respectively. All the other datasets were used for testing.

	For the proposed hybrid network, we experimentally set $\lambda_{det} = 0.4$ (and $\lambda_{loc} = 0.6$), that is the setup that gave the best results.
	%
	To evaluate the detection effectiveness, we consider the true positive (TPR) and true negative rates (TNR) of the decision.

	We also calculate the area under the ROC curve (AUC) performance by pairing each fake image dataset with the real image dataset. To evaluate the localization performance, we use the balanced Pixel Accuracy (bPA) and Intersection over Union (IoU), which for a binary classifier are defined as:
	\begin{equation}

		\text{bPA}=\frac{1}{2N_t}\sum_{t=1}^{N_t}\sum_{i=0}^{1} \frac{p_{ii}}{\sum_{j=0}^{1}p_{ij}},

	\end{equation}
	and
	\begin{equation}

		IoU = \frac{1}{N_t}\sum_{t=1}^{N_t}\frac{\left|\left\{ M_t \equiv 1 \right\} \cap \left\{ G_t \equiv 1 \right\} \right|}{\left|\left\{ M_t \equiv 1 \right\} \cup \left\{ G_t \equiv 1 \right\} \right|},
	\end{equation}
	where $p_{ij}$ indicates the number of pixels in class $i$ classified as class $j$, $\{M_t \equiv 1\}$ is the set of pixels of value 1 in the $M_t$ mask images,

	and $N_t$ indicates the number of test images.

	\subsection{Results}

	\begin{table}
		\centering
		\footnotesize
		\renewcommand\arraystretch{0.8}
		\setlength{\tabcolsep}{0.55mm}

		\caption{Results (\%)  for different datasets.}
		\vspace{-0.2cm}

		\begin{tabular}{@{}lccccccc@{}}
			\toprule
			& WSOC & \multicolumn{2}{c}{StreetG} & \multicolumn{2}{c}{WebG132} & \multicolumn{2}{c}{WebG504} \\
			\midrule
			Methods & TNR\% & TPR\% & AUC\% & TPR\% & AUC\% & TPR\% & AUC\% \\
			\midrule
			Xception & 96.8 & 98.0 & 99.5 & 37.9 & 84.3 & 57.7 & 91.8 \\
			ResNet50 & 95.0 & 98.6 & 99.4& 47.0 & 84.0 & 64.9 & 91.5 \\
			Xception+M (CAT) & 94.5 & 86.4 & 96.6 & 41.7 & 75.5 & 65.9 & 89.2\\
			Xception+M (MUL) & 96.8 & 97.5 & 99.5 & 65.1& 95.9 & 73.4 & 96.9\\
			ResNet50+M (CAT) & 92.0 & 99.6 & 98.7 & 47.0 & 83.0 & 60.5 & 88.6\\
			ResNet50+M (MUL) & 95.2 & 98.3 & 99.3 & 70.5 & 96.1 & 79.2 & 96.7\\

			ResNet50\_hyb & \textbf{98.0} & \textbf{100} & \textbf{100} & \textbf{93.4} & \textbf{99.0} & \textbf{95.4} & \textbf{98.9}\\
			\midrule
			Methods & bPA & bPA & IoU & bPA & IoU & bPA & IoU \\
			\midrule

			ResNet50\_hyb & 98.6 & 96.1 & 92.0 & 84.4 & 62.5 & 87.9& 71.2\\
			\bottomrule
		\end{tabular}

		\label{tab:1}
		\vspace{-0.1cm}
	\end{table}

	\begin{figure}
		\centering
		\includegraphics[width=0.9\columnwidth]{figures/Vis.pdf}
		\caption{Visualization of output masks and CAM maps for ResNet50\_hyb.}
		\label{fig:mask}
		\vspace{-0.2cm}
	\end{figure}

	\enlargethispage{\baselineskip}

	\begin{figure*}[h!]
		\centering

		\includegraphics[width=0.495\linewidth]{figures/JPEGandResize.pdf}

		\includegraphics[width=0.495\linewidth]{figures/BlurAndFilter.pdf}
		\caption{Robustness results (\%) in terms of TNR, TPR, and AUC in the presence of JPEG compression with a quality factor of 50 (top, 1-4), resizing with factor 0.5 (bottom, 1-4), Median filtering, window size = 3$\times$3 (top, 5-8) and Gaussian blur, window right size = 3$\times$3(bottom, 5-8).}
		\label{fig:robust2}
		\vspace{-0.2cm}
	\end{figure*}
	%

	Table \ref{tab:1} shows the results on the four test datasets. The baseline methods achieve good classification performance on the StreetG dataset, with a TPR around 98\% and an AUC equal to 99\%. However, the performance dropped on the datasets that were not used for training. The performance of ResNet50 and Xception are the worst, likely due to the small number of images available for training.

	The way the water mask images are used during training also has a significant impact on the performance. By looking at the results in rows 3 to 6, multiplying the mask and the input image results in better performance on WebG132 and WebG504 datasets, than simply concatenating the mask to the image bands.

	%
	The hybrid detection/localization method proposed in this paper (ResNet50\_hyb) achieves the best performance on all the datasets, with accuracies always well above 90\% and often close to 100\%. In particular, the proposed method shows a good generalization capability, achieving very good performance on WebG132 and WebG504 datasets, with 93.4\% and 95.4\% TPR values, and 99.0\% and 98.9\% AUCs. This is a remarkable result, given that such good performance is obtained without the additional information provided by the water mask.
	%
	The superior performance provided by the hybrid model is even more interesting if we consider the localization accuracy. By looking at the bPA and IoU values reported in

	the last row of Table \ref{tab:1}, we can see that the proposed network provides good results on WSOC and StreetG, even if we can notice a performance drop on WebG132 and WebG504 datasets. The remarkable conclusion we can draw is that asking the network to localize the tampered regions helps to distinguish fake images from real ones, even when localization does not work very well,
	since it forces the detection branch to focus on the manipulated area.
	Such intuition is confirmed by the CAM maps \cite{selvaraju2017grad}, computed with regard to the detection task, shown in the last row of Fig. \ref{fig:mask}.

	%

	Table \ref{tab:3} and Figure \ref{fig:robust2} show the robustness of the methods against various image processing operators, including JPEG compression with quality factor 50, image downsampling with resizing 0.5, Median filtering (3$\times$3), Gaussian blur (3$\times$3), and Gaussian noise addition with zero mean and variance equal to 0.003\footnote{In all our experiments images take values in the [0,1] range.}.
	%

	Overall, all the models are robust to the image processing methods used in our tests, with the exception of Gaussian noise addition (see Table \ref{tab:3}). In the Gaussian noise case, in fact, the performance of the Xception+mask (MUL) model dropped from 97.5\%, 65.1\% and 73.4\% to 1\%, 0\%, and 0\%.
	ResNet50+mask (MUL) also suffers the addition of Gaussian noise. The bad behaviour of the models adopting a multiplicative approach to include the information provided by the water mask indicates that, in the absence of information about the non-tampered regions, the impact of noise is amplified, since the texture information only is not enough to detect the presence of tampering. As further evidence of this intuition, we observe that Xception and ResNet50 show stronger robustness. The positive effect of considering the entire image rather than the water areas only is evident also with regard to the proposed hybrid method, which, as it can be seen from Table \ref{tab:3}, retains good performance also in the presence of noise addition.
	%

	With regard to the other attacks (see Figure

	\ref{fig:robust2})
	the proposed method always achieves the best results with a minor performance drop on the WSOC real dataset only.

	\begin{table}
		\centering
		\footnotesize
		\renewcommand\arraystretch{0.8}
		\setlength{\tabcolsep}{0.6mm}

		\caption{Results (\%) in the presence of Gaussian noise.
		\vspace{-0.2cm}
		\begin{tabular}{@{}lccccccc@{}}
			\toprule
			& WSOC & \multicolumn{2}{c}{StreetG} & \multicolumn{2}{c}{WebG132} & \multicolumn{2}{c}{WebG504} \\
			\midrule
			Methods & TNR\% & TPR\% & AUC\% &TPR\% & AUC\% & TPR\% & AUC\% \\
			\midrule
			Xception & 98.0 & 77.7 & 98.0 & 18.9 & 78.5 & 34.7 & 88.2 \\
			ResNet50 & 98.4 & 94.4 & 99.3 & 35.6 & 85.1 & 47.2 & 92.0 \\
			Xception + M (CAT) & 90.4 & 60.0 & 83.1 & 25.6 & 56.3 & 46.0 & 76.9\\
			Xception + M (MUL) & \textbf{100} & 1.0 & 93.9 & 0 & 91.9 & 0 & 92.5\\
			ResNet50 + M (CAT) & 94.5 & 89.5 & 97.6 & 26.5 & 78.4 & 41.1 & 86.0\\
			ResNet50 + M (MUL) & 97.5 & 34.1 & 89.5 & 15.1 & 90.5 & 15.1 & 90.8\\

			ResNet50\_hyb & 97.0 & \textbf{100} & \textbf{100} & \textbf{84.1} & \textbf{97.6} & \textbf{86.5} & \textbf{97.8}\\
			\bottomrule
		\end{tabular}

		\label{tab:3}
		\vspace{-0.5cm}
	\end{table}

	\section{Conclusion}
	\label{sec:conclusion}

	In this paper, we have presented a ResNet50-based hybrid network for the detection of GAN-generated flood images. In particular, we exploited localization during training to help the network to focus on the most significant parts of the analyzed image, with a positive impact on the detection performance. To train and test the proposed network, we have gathered more than 4,000 fake flood images generated by \textit{ClimateGAN} and 880 genuine flood images. The experiments we carried out show the excellent performance of our network, its good generalization capabilities, and its robustness against various image processing operations.
	As a future work, we are planning to consider different generative models. Being our approach general, we also plan to extend it to the detection of different GAN-generated climate disasters and forgeries.

	Meanwhile, it would be interesting to combine this work with geolocalization, that is, given a flood image, first localize the geographical location the image is referring to, and then verify the plausibility that a flood actually occurred.

	%

	%

	%

	%

	%

	\section*{Acknowledgment}

	This work has been partially supported by the China Scholarship Council(CSC) , file No. 202008370186, by the
	PREMIER project under contract PRIN 2017 2017Z595XS-
	001, funded by the Italian Ministry of University and Research,
	and by the Defense Advanced Research Projects Agency (DARPA) and the
	Air Force Research Laboratory (AFRL) under agreement number FA8750-20-2-1004. The U.S. Government is authorized to
	reproduce and distribute reprints for Governmental purposes
	notwithstanding any copyright notation thereon.

	\bibliographystyle{IEEEtran}


\subsection{Models}

\par{To generate the periodic model, the hybrid functional PBE0 was used, with the "light" basis sets of FHI-aims. The light basis sets were sufficient in this case because they give a sufficiently accurate representation of the geometric and electronic structure and because truly accurate values for energies were unnecessary. A $5\times 5\times 5$ $k$-grid was needed for meV convergence. }

\par{The QM/MM model was built from the relaxed structure of the periodic model following the procedure outlined in previous work \cite{RN8650} and was found to be reliable for studies of other metal oxides. The chosen QM-region was a $3\times 3\times 3$ grid of W\textsuperscript{6+} ions, and all the O\textsuperscript{2$-$} ions between them, totaling 81 atoms. }

\par{The potential energy functions used in the MM region were the Islam 1998 potentials \cite{RN8651}. These potentials involve Buckingham potentials between ions and a harmonic potential between anion cores and their electron shell. These simple potentials give the correct formation energies and dielectric constants for WO\textsubscript{3} but cannot replicate the complex structure of room-temperature monoclinic I ($\gamma$-WO$_3$) form. The partially covalent nature of WO\textsubscript{3} and the different lengths for its six W--O bonds due to pseudo Jahn--Teller distortions cannot be replicated with such simple potentials. As a consequence, all the MM nuclei needed to be fixed during energy minimization to prevent instability; if allowed to relax, the active MM region would shift to a cubic structure, creating spurious interactions at each of its boundaries. In addition, any ions in the QM region bound to multiple ions in the MM region were also fixed. }

\par{Such a strict fixing of ions is necessary for the model to work but, importantly, limits what it is capable of modeling. This model can only determine structural changes within a single WO\textsubscript{6} octahedron. If the correct structural change should be a particularly large-scale one, with multiple WO\textsubscript{6} octahedra distorting simultaneously, a substantially larger---and potentially computationally infeasible---calculation would be required. However, any understanding gained from knowing the electronic structure at the valence and conduction band boundaries, as well as the behavior of a fully localized system, can still be insightful.}

\par{In the excited state characterized by the XAFS studies, an electron is excited from the valence band into the conduction band. This excitation should result in the formation of an electron hole and a charge carrier. To investigate how this excitation process affects the structure and bond length, the neutral structure was energy-minimized in a periodic QM/MM model and then the structure with a hole was energy-minimized using QM/MM, then again for a carrier. This approach should reveal which (if either) state would influence bond lengths.}

\par{The HOMO is nonbonding with respect to any W--O bond; however, the LUMO is bonding with respect to the W-O5 bond, and antibonding with respect to the W-O2 bond.}

\par{Bond lengths for W-O1 through W-O6 were calculated for the relaxed neutral QM/MM models. The electron--hole and charge-carrier models varied from the neutral models by containing one fewer or one more electron in the QM region, respectively. The results are shown in Table \ref{fig:bond_lengths}, and are all expressed in angstroms. Figure \ref{fig:WO3MSPPT} shows the structures of the GS and MS WO$_3$. }

\par{Reassuringly, the results for the neutral QM/MM models well match the experimental values. The maximum variance between models for any one bond length is approximately 0.04 \AA.}

\par{Between the neutral and hole QM/MM models, the bond lengths do not substantially change. This result was expected because removing an electron from an orbital that is neither bonding nor antibonding with respect to any W--O bond should not strongly affect any W--O bond length. The carrier QM/MM models, by contrast, showed a lengthening of the W-O2 bond and a shortening of the W-O5 bond compared with the neutral model. This result is expected because the LUMO orbital is bonding along the W-O2 bond and antibonding along the W-O5 bond. A spin-density analysis of the carrier model showed that the additional electron is mostly located on the central W (spin +0.52), with the rest located on the nearby W (spin +0.1 to +0.2).}

\begin{center}
\begin{table}
\caption{Values for W--O bond length for a range of models.}
\begin{tabular}{c c c c}
 \hline
 \ \ \ & \ \  QM/MM \ \ & \ \  QM/MM\ \  & \ \ QM/MM\ \  \\

 & Neutral & Hole & \underline{\textbf{Carrier}} \\
 \hline
 W-O1 (c)  & 1.78 & 1.77 & 1.78 \\
 W-O2 (b)  & \textbf{1.78} & 1.79 & \textbf{1.85} \\
 W-O3 (a) & 1.87 & 1.87 & 1.88 \\
 W-O4 (a) & 1.88 & 1.88 & 1.89 \\
 W-O5 (b) & \textbf{2.09} & 2.07 & \textbf{1.99} \\

 W-O6 (c) & 2.16 & 2.17 & 2.18 \\
 \hline
\end{tabular}
\label{fig:bond_lengths}
\end{table}
\end{center}

\begin{figure}[hbt!]
\includegraphics[width=0.7
\textwidth]{WO3MSPPT}
\caption{QM/MM results for structures of GS and MS WO$_3$}
\label{fig:WO3MSPPT}
\end{figure}

\section{\label{sec:level1}Introduction}
Semiconductor oxides are used as photocatalysts to harvest sunlight to split water into H$_2$ and O$_2$.\cite{Strunk2021}
Electrons in the valence band, which is mainly composed of O 2$p$ orbitals, are excited by photoirradiation to the conduction band, followed by stabilization of the photoelectrons to form the metastable state (MS), which is often assumed to be associated with surface or bulk defect sites.\cite{Strunk2021}
The local structure and electronic state of the MS need to be elucidate to understand the photocatalysis of water splitting and to improve the performance of water-splitting photocatalysts.

A pump--probe (PP) X-ray absorption fine structure (XAFS) investigation of the photocatalysts using pulsed X-rays emitted from an X-ray free electron laser (XFEL) and synchrotron radiation (SR) has revealed element-specific local electronic and geometrical structures, thereby providing  new scientific knowledge about photocatalysis that cannot be obtained by other techniques.\cite{RN7738}
WO$_3$ is used as a Z-scheme water-splitting photocatalyst that is sensitive to visible light, which is the main component of sunlight.\cite{abe2010}
We have acquired PP X-ray absorption near edge structure (PP-XANES) spectra of WO$_3$ after pulsed laser irradiation and have found that the MS was created within 150 ps, followed by relaxation to the ground state (GS) with a time constant of 2 ns.\cite{RN257,RN25,RN2588,RN7257}
We proposed that the photoelectrons were trapped at distorted stoichiometric W sites in the crystal lattice, similar to the trapping of polarons, instead of at surface or bulk defects.\cite{RN2588}
The formation time for the MS, 150 ps, was longer than that for the other oxides such as TiO$_2$ and Fe$_2$O$_3$, consistent with previously-reported optical measurement results\cite{RN2269}.
It is desirable to acquire the PP extended X-ray absorption fine structure (PP-EXAFS) spectrum to confirm the structure change because EXAFS is sensitive to local structural features such as bond lengths.\cite{RN1856}
However, PP-EXAFS spectral analysis is much more difficult than PP-XANES spectral analysis \cite{RN2588} because of the lower PP-EXAFS signal intensity and because the multishell fitting required for compounds with a complex structure, such as WO$_3$, is inherently challenging.
We addressed these two difficulties by taking the difference spectra between before and after photoabsorption to obtain the EXAFS spectrum of MS WO$_3$.
The difference spectra were  accumulated and then analyzed using the constrained thorough search (CTS) method to elucidate the bond lengths in MS WO$_3$.\cite{RN7005,RN8640}

In the present paper, we report the PP-EXAFS results for the structure of the MS WO$_3$ state to clarify the local structure change and to attempt to explain the long formation time for MS WO$_3$.

The WO$_3$ was purchased from Wako Chemicals, and 50--200 nm WO$_3$ particles were dispersed in ultrapure water.
The concentration of WO$_3$ was 2.5 mmol L$^{-1}$.
The WO$_3$ jet was supplied to the cross point of the laser used for photoexcitation (Ti-sapphire with a 1 ps pulse width, 945 Hz repetition rate, and 400 nm wavelength) and the X-ray beam emitted from the Photon Factory Advanced Ring (PF-AR, 6.5 GeV, 60 mA, single bunch operation with pulse width of 100 ps, and pulse interval of 1.26 $\mu$s).
The fluorescence X-rays were monitored using a scintillation counter with a Cu filter to attenuate the elastic X-rays.
Details of the experimental setup are available elsewhere.\cite{RN257}
The difference spectrum before and after the photoirradiation, $\Delta \mu(=\mu_{API}-\mu_{BPI}) $, was accumulated, thereby reducing the possible fluctuations.
The $\mu_{API}$ and $\mu_{BPI}$ are the spectra recorded after and before the photoirradiation, respectively.
The spectrum of the MS WO$_3$, $\mu_{MS}$, was obtained by assuming that $\mu_{API}$ was a linear combination of the spectrum of the ground state (GS), $\mu_{GS}$, and $\mu_{MS}$.
${\mu}_{API}={{\alpha\mu}_{MS}+{\left(1-\alpha\right)\mu}_{GS}}$, where $\alpha$, the fraction of MS WO$_3$, was obtained from the XANES spectrum, as shown in Fig. \ref{fig:XANES} in the supporting information (\ref{sec:XANES}).
${\mu}_{MS}$ was processed and analyzed using the EXAFS analysis programs REX\cite{RN3616} and LARCH \cite{RN6645}, as described in \ref{sec:Analysis}.
The Fourier transformation ranges of $\Delta k$ and $\Delta r$ were 5.5 ($=8.5-3$) \AA$^{-1}$  and 1 ($=2-1$) \AA \ , respectively, so that the amount of information, $M$, was approximately 5 ($M=2 \Delta k \cdot \Delta r)/\pi + 2)$\cite{PhysRevB.48.9825}.
The inversely Fourier transformed data were first analyzed by a curve-fitting method using the LARCH package.\cite{RN6645}
Further  details of the structure, including bond lengths, were derived from the EXAFS spectra using the CTS method, as explained in \ref{sec:CTS}.\cite{RN7005,RN8640,RN7246}

Fig. \ref{chi_1} shows the EXAFS oscillation of the MS WO$_3$, $\chi(k)_{MS}$, and its Fourier transform.
The main peak at 1.4 \AA \  was inversely Fourier transformed to the $k$-space, and the data were analyzed by one-shell curve fitting.
The one-shell curve fitting results are shown in Fig. \ref{fig:CF} and Table \ref{tab:CF1} in \ref{sec:RCF}.
We conducted the analysis using two methods.
In the first method, the Debye--Waller factor was adjusted freely; in the second method, the Debye--Waller factor was constrained as a function of the bond length given by Eq.\ref{eq:DWANDR} in \ref{sec:CTS}.
Both results indicated that the bond length contracted in the MS WO$_3$; however, details of the structure were unclear.

\begin{figure}
 \centering
  \includegraphics[width=\linewidth,keepaspectratio,clip]{Figure1a.eps}
  \vspace{-2.5 cm}
  \caption{(a)$\Delta\chi\left(k\right)$ for W L$_3$-edge EXAFS spectra between before and 150 ps after the irradiation. (b) $\chi_{MS}\left(k\right)$  for the MS WO$_3$.  (c) The Fourier transform of $\chi_{MS}\left(k\right)$. The broken line shows the range of Fourier filtering. (d) Inversely Fourier transformed data(black line) for the MS WO$_3$ with the results of CTS(red line).  The enlarged one is available in Fig. \ref{fig:TSFITTING}}
  \label{chi_1}
\end{figure}

We carried out the CTS analysis based on the uniform prior probability principle.\cite{RN8640}
We assumed the following:
(1) MS WO$_3$ has uniformly distorted octahedra, similar to GS WO$_3$.
(2) Debye--Waller (DW) factors are functions of the bond length in Eq.\ref{eq:DWANDR} in \ref{sec:RCF}.\cite{RN1924} A greater length corresponds to a larger DW factor.
(3) EXAFS spectroscopy is sensitive to short bonds because of the $\exp (-2r/\lambda -2 \sigma^2 k^2 ) / r^2 $ term.

The structural parameters were limited to the value of $M$ ($\approx$5).
The GS WO$_3$ has a monoclinic structure in which two distorted octahedral sites are present, as shown in Table \ref{tab:CFR}.\cite{RN8155}
The bond lengths in GS WO$_3$ have 12 different values.
Two sites have a similar local structure; the corresponding W--O distances are therefore denoted as $r'''_1$, $r'''_2$, $r'''_3$, $r'''_4$, $r'''_5$, and $r'''_6$ in Table \ref{tab:CFR} in \ref{sec:CTS}.
We intended to demonstrate that $r'''_1 \ne r'''_2$  in the MS WO$_3$, as predicted from the W L$_1$-edge XANES results.\cite{RN2588}
According to assumption (3) (i.e., that shorter bonds ($r'''_1$, $r'''_2$) have greater contributions to the EXAFS spectra), we assumed that $r'''_1$ and $r'''_2$ were independently optimized and we constrained $r'''_3 = r'''_4 = r_3$ and $r'''_5 = r'''_6=r_4$, as shown in Table \ref{tab:CFR}. The CTS was carried out for the four parameters $r_1(=r'''_1)$, $r_2(=r'''_2)$, $r_3$, and $r_4$ with the corresponding coordination numbers, $N_1=N_2=1$ and $N_3=N_4=2$, as shown in Table \ref{tab:tableS1} and Table \ref{tab:CFR} in \ref{sec:CTS}.
Table \ref{tab:cts} shows the CTS results that reproduced the crystal data given in Table \ref{tab:CFR} in \ref{sec:CTS} for the GS WO$_3$.
Note that $r_1 $ and $r_2$ are approximately the same as shown in Table \ref{tab:cts} even if the two parameters are independently changed.
Fig. \ref{fig:histoGS} shows the occurrence histograms.
We found that the histograms for $r_1$ and $r_2$ have similar distributions and the same peak positions.
Consequently, for the GS WO$_3$, we conclude that $r_1 $ and $r_2 $ agreed with each other within the error.

Table \ref{tab:cts} shows the CTS results for the MS WO$_3$.
Fig. \ref{fig:TSFITTING} shows a comparison between the Fourier filtered $\chi_{MS}(k) $ and the calculated one.
The observed data were well reproduced by the parameters determined by the CTS method.
The R-factor was 0.02.
In the case of the MS WO$_3$, we found that $r_1$ and $r_2$ were not equal and that $r_1$ was 0.11 \AA  \ shorter than that for the GS WO$_3$.
The histograms of $r_1$ and $r_2$ show different peak positions (Fig. \ref{fig:HistoMS}).
Bonds other than $r_1$ in the MS WO$_3$ were  longer than the corresponding bonds in the GS WO$_3$, in good agreement with previously reported WO$_3$ L$_1$ XANES results.\cite{RN2588}
The photoexcited MS WO$_3$ exhibited a further distorted structure.

\begin{figure}

 \centering
    \includegraphics[width=80 mm]{HISTOGRAMMS2a.eps}
    \caption{(a)--(d) Distribution of bond lengths that satisfy R-factor $<$ 0.10 for the structure parameters for the MS WO$_03$.  Dashed lines show the fitting results obtained using a Gaussian function. }
    \label{fig:HistoMS}
\end{figure}
\begin{table}[htbp]
    \caption{Results of constrained thorough search analysis for the W L$_3$-edge EXAFS spectra of GS and MS WO$_3$.    "Crystal" indicates crystallographic data.\cite{RN8155}  }
        \label{tab:cts}
    \centering
    \begin{tabular}{c c c c c  c}
    \hline
        sample &$r_1$ & $r_2$ & $r_3$ &$r_4$ & \begin{tabular}{c}
R-\\factor \end{tabular}\\
         \hline
GS &1.77$\pm 0.08 $&1.77$\pm 0.08$ &1.91$\pm 0.04 $ &2.17$\pm 0.09$ &0.049\\
         Crystal & (1.75 & 1.77 & 1.90 & 2.13 )&\\
         \hline
         MS
        & 1.66$\pm 0.06$ &1.81$\pm 0.04$ &1.93$\pm 0.06$ &2.19$\pm 0.07$ & 0.012\\
         \hline
    \end{tabular}
\end{table}

The CTS analysis showed that the MS WO$_3$ exhibited a further distorted structure, where the shortest W--O bond was further shortened from 1.77  \AA \ to 1.66  \AA. \cite{RN2588}
This structure change was interpreted as follows in a previous paper.\cite{RN2588}
The photoelectron excited by the pulsed laser would occupy the d$_{xy}$ orbital of W, which is on a plane perpendicular to the shortest W=O bond.  The orbital is located at the bottom of the conduction band in energy scale, forming a $\pi^*$ antibonding orbital with 2$p_x$ and 2$p_y$ for oxygen atoms in the $xy$ plane.
Thus, the W--O bonds in the $xy$ plane were elongated, pushing up the W atom in the direction of the shortest W--O bond.\cite{RN2588}
As described later and in \ref{sec:MO} and in \ref{sec:qmmm}, this previous interpretation is not consistent with the results of
the density functional theory (DFT) and quantum mechanics/molecular mechanics (QM/MM) calculations.

In the literature, the MS WO$_3$ was assumed to be located in bulk or surface defect sites.\cite{RN1799,RN8356, RN2222}
However, we disproved these possibilities because the amount of MS states is much larger than the amount of defects.
Yamagata {\it et al.} recently reported that introducing defects by reducing WO$_3$ decreased the lifetime of the MS WO$_3$.\cite{RN8614}
Surface defects might be less likely because the surrounding environment did not affect the lifetime of the MS WO$_3$, as shown in \ref{sec:react}.

We proposed the MS WO$_3$ should be located at the normal lattice sites in the bulk stoichiometric WO$_3$ and be distorted from the structure of the GS WO$_3$.
Such distortion would stabilize the MS WO$_3$, similar to the stabilization of polarons.
\begin{figure}[htbp]
 \centering
  \includegraphics[width=\linewidth,keepaspectratio,clip]{Fig3a.eps}
  \caption{The local structure change in the GS WO$_3$ (left panel) and the MS WO$_3$ ( right panel).   Large red and small blue circles represent O and W atoms, respectively.  yellow broken circles indicate the position of W in GS. Thick blue arrows show the movement of W atoms in the MS WO$_3$.  W atoms along the $c$-axis move in the direction corresponding to W=O contraction. Orange filled and empty arrows in the orbitals correspond to electrons with spin and holes, respectively.  Grey filled circles indicate holes.   The unpaired electron spins are all parallel in the polaron cluster.}
  \label{fig:polaron}
\end{figure}

However, we still faced two problems.
First, the DFT and QM/MM calculations for the MS WO$_3$ local structure indicated that the less-distorted structure was stable, with the shortest W=O bond being elongated as shown in Table \ref{tab:GAUSSIAN} in \ref{sec:MO}.
Second, it was unclear why the formation of the MS WO$_3$ structure took 150 ps.
The vibration frequencies of W--O are on the order of 100 cm$^{-1}$ or 3 THz or 0.3 ps.
Fe$_2$O$_3$ hematite and CuWO$_4$ took less than 1 ps to produce a polaronic MS.\cite{RN7556,RN8560}
Thus, the simple isolated distorted MS structure should form much faster.
According to the L$_3$ edge result, the W$^{5+}$ was formed and was gradually changed to the MS structure.\cite{RN25}
Such a long time scale might be related to the collective structure change.
If the single distorted structure of the MS WO$_3$ was located in the GS WO$_3$, it should have a greater distortion energy between the surrounding GS WO$_3$.
Thus, the theoretical calculation suggested a less distorted structure.
If the two distorted MS WO$_3$ structures gathered and were aligned along the $c$-axis, along which the shortest W=O bonds were assumed to be headed, the distortion energy should be released.
Consequently, the formation of the well-ordered distorted MS WO$_3$ cluster shown in Fig. \ref{fig:polaron} should require a long time.

The problem with the well-ordered MS WO$_3$ clusters is the charge balance, which should show a large negative charge.
To compensate for the negative charge, holes should combine with photoexcited electrons in the well-ordered MS WO$_3$ cluster.In this case, the recombination of photoexcited electrons and holes might be enhanced.
If the spins of the hole and the photoexcited electron near the hole are  parallel or in a triplet state, the well-ordered distorted MS cluster should be stabilized and have a certain lifetime; however, we lack direct evidence supporting this scenario.
If we could carry out large-scale DFT calculations of the MS excited state with a large cluster size, we could support our hypothesis for the formation of such a well-defined distorted MS WO$_3$ cluster.
Such calculations will be the topic of future work.

Picosecond time-resolved W L$_3$-edge EXAFS spectroscopy was conducted to characterize the structural change in MS WO$_3$ in detail.
The change of the W L$_3$-edge EXAFS spectra was characterized.
We applied the CTS method to analyze the EXAFS data for the MS WO$_3$, where one of the shortest bonds was further shortened by 0.11 \AA \  from its original bond length, consistent with the structure proposed on the basis of W L$_1$-edge XANES spectroscopy.
This structural change in the MS WO$_3$ led the distorted MS structure, which should be stabilized, to form well-ordered distorted clusters in conjunction with the formation of the triplet state.

The authors would like to express their gratitude to Prof. Masaaki Yoshida at Yamaguchi University for his technical advice and to the Japan Society for the Promotion of Science (JSPS) for their support through a Grant-in-Aid for Exploratory Research (No. 26620110), a Grant-in-Aid for Scientific Research (A) (Nos. 15H02173 and 20H00367), a Grant-in-Aid for JSPS Research Fellow (No. 15J07459), and a Grant for collaborative research in the Institute for Catalysis, Hokkaido University (No. 15A1004).
The EXAFS experiments were carried out at beamline station BL-NW14A with the approval of PF-PAC (Prop. Nos. 2013G166, 2015G541, 2015G542).

See Supplemental Material at [URL will be inserted by publisher] for \ref{sec:XANES}Estimation of the ratio of photoexcitation;\ref{sec:Analysis} EXAFS data processing and the curve-fitting analysis;\ref{sec:RCF} Results of curve fitting;\ref{sec:CTS} Constrained thorough search; \ref{sec:GSWCTS} CTS analysis of the EXAFS for the WO$_3$ in the ground state ;\ref{sec:react} Surroundings-controlled W L$_3$-edge XANES measurements; \ref{sec:MO} Details of density functional theory (DFT) structure optimization for a single WO$_6$ complex in monoclinic and orthorhombic environments; \ref{sec:qmmm} Details of the QM/MM structural optimization for neutral, hole, and carrier states.