\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.  

