text:
\section{Introduction}
The optical properties of noble metal nanoparticles are dominated by their ability to support localized surface plasmon (SP) excitations, which can be described as the collective oscillation of the free-electron gas confined to the metal surface. Coupling between light and SPs gives rise to numerous interesting phenomena, such as squeezing light beyond the diffraction limit\cite{Gramotnev:2010} and large enhancements of the local electric field.\cite{Kelly:2003} These SP-induced effects lead to several technological applications, including e.g. improvement of absorption in solar cells,\cite{Atwater:2010} on-chip routing of electromagnetic energy,\cite{Zia:2006} and sensing of biological molecules.\cite{Anker:2008} The ability to probe the plasmonic response of individual nanoparticles is of crucial importance for designing novel structures that can harvest the full potential of plasmonics. One promise of plasmonics is that functional structures can be made much smaller than the wavelength of light, but subwavelength spatial information cannot be obtained with the usual diffraction-limited optical experiments. In contrast to light-based measurement techniques, electron energy-loss spectroscopy (EELS)\cite{Abajo:2010} performed in a transmission electron microscope (TEM) offers the ability to map the plasmonic resonances of sub-wavelength metal nanostructures\cite{Nelayah:2007,Bosman:2007} and in ultra-confined geometries\cite{Raza:2014a} due to the \AA{}ngstr\"{o}m spatial resolution of the TEM and tightly-confined electromagnetic field generated by the swiftly-moving electrons.\cite{Abajo:2010} In conjunction with an electron monochromator, an energy resolution down to a few hundred meV (and decreasing\cite{Krivanek:2014}) can be routinely achieved, allowing for high spatial- and spectral-resolution studies of individual plasmonic structures.

The localized SP resonance energies of metal nanoparticles can be selectively controlled, such as by changing the size, shape, material or the environment of the nanoparticle.\cite{Kelly:2003} However, even for a fixed system, a metal nanoparticle can according to classical electrodynamics show several localized SP resonances corresponding to excitations of different multipolar order. In particular, a spherical metal nanoparticle of radius $R$ with permittivity $\varepsilon(\omega) = \varepsilon_\text{core}(\omega) - \omega_\text{p}^2/\omega^2$ fully embedded in a dielectric medium $\varepsilon_\textsc{b}$ has the non-retarded resonance energies $\omega_l$ governed by\cite{Sernelius:2001}
\begin{equation}
	\omega_l = \frac{\omega_\text{p}}{\sqrt{\varepsilon_\text{core}(\omega_l) + \frac{l+1}{l}\varepsilon_\textsc{b}}}, \label{eq:res_freq_nr}
\end{equation}
where $\omega_\text{p}$ is the plasma frequency of the bulk metal, $\varepsilon_\text{core}(\omega)$ is the frequency-dependent response due to the bound charges of the metal, and, importantly, $l$ denotes the angular momentum of the SP mode, where $l=1$ is the dipole mode, $l=2$ is the quadrupole mode, and so on. Far-field measurement techniques based on plane-wave scattering mainly probe the dipole mode (see Supplementary Note 1), while near-field techniques such as EELS provide access to modes of larger angular momentum as well.\cite{Abajo:2010,Christensen:2014} This property has been advantageously used in many EELS experiments to map and study the multipolar plasmonic response of different metal nanostructures.\cite{Nelayah:2007,Bosman:2007,NGom:2008,Nelayah:2009,Schaffer:2009,Chu:2009,Koh:2011,Nicoletti:2011,Alber:2011,Guiton:2011,Rossouw:2011,Bigelow:2012,Nicoletti:2013,Rossouw:2013,Wiener:2013,Zhou:2014,Tan:2014,Martin:2014} However, study of the higher-order (HO) modes (i.e., $l>1$) in small noble-metal nanoparticles ($R<25$~nm) is hindered by interband transitions and the screening of bound charges [accounted for through $\varepsilon_\textrm{core}(\omega)$ in Eq.~\eqref{eq:res_freq_nr}], which strongly dampens the modes. Additionally, the energy spacing between subsequent $l$ modes is decreased compared to simple metals where $\varepsilon_\textrm{core}\approx 1$, making it difficult to spectrally resolve the modes. These issues can be overcome by encapsulating the nanoparticles in an insulating medium, which shifts the localized SP resonances to lower energies where interband transitions and losses are less pronounced.\cite{Kreibig:1987} EELS measurements on ensembles of encapsulated potassium\cite{vomFelde:1988} and silver nanoparticles\cite{Kreibig:1970} have shown the presence of HO modes, yet unambiguous measurements on the single-particle scale have to our knowledge not been achieved.

On the single-particle scale, several issues of ensemble measurements such as an inhomogeneous size distribution and the coupling between neighbouring particles are automatically prevented. In addition, deterministic studies of single particles constitute an ideal scenario for direct comparison with theoretical EELS simulations, allowing for thorough examination of electrodynamic theory on the few-nanometre scale, where the validity of the constitutive relations in classical theory is debated.\cite{Duan:2012,Scholl:2012} In previous studies, EELS measurements on small metal nanoparticles\cite{Kreibig:1970,vomFelde:1988,Ouyang:1992,Scholl:2012,Raza:2013} were compared with simulations based on optical scattering, thereby neglecting the importance of the position of the electron beam (i.e., impact parameter), and therefore the excitation of HO modes. Theoretical EELS studies on nanoparticles have shown that positioning the electron beam close to the nanoparticle surface more strongly excites SP modes of all angular momenta including HO modes,\cite{Abajo:1999,Abajo:2010,Christensen:2014} providing theoretical evidence that HO modes should be taken into account (through EELS simulations) in the analysis of EELS data on individual small nanoparticles.

In this paper, we report EELS measurements on multi-twinned silver nanoparticles of icosahedral shape encapsulated in silicon nitride (see Fig.~\ref{fig:fig1} for illustration). We present experimental evidence of the excitation of both the dipole ($l=1$) and HO modes ($l>1$) in single silver nanoparticles in the radius range $1-20$~nm. Remarkably, we find that the HO modes can be identified down to a particle radius of only $4$~nm. To our knowledge, this is the first unambiguous experimental observation of HO modes in single silver nanoparticles of this ultrasmall size. We also provide experimental evidence of the significant dependence on the position of the electron beam when acquiring EELS spectra of small silver nanoparticles. Now, in previous studies it was reported that the relative weight of bulk and SP resonances can be controlled by the impact parameter.\cite{Scholl:2012} The experimental novelty here is that the impact parameter moreover also strongly influences the relative weight of the HO resonances. Even more surprisingly, for impact parameters close to the particle surface the spectral weight of the HO modes even exceeds that of the dipole mode. These results demonstrate that HO modes have unanticipated implications for the interpretation of EELS spectra of small nanoparticles.

Comparison with retarded EELS simulations based on classical electrodynamics, where we also take into account the energy resolution of the experimental setup, shows excellent agreement with measurements of nanoparticles down to $4$~nm in radius. Only for particle radii below $4$~nm do we observe a discrepancy between measurements and simulations, providing experimental evidence for a lower size limit for the application of classical electrodynamic theory. At this size scale the EELS signal from the HO modes disappears and we find only a single resonance, interpreted as the dipole mode, which strongly increases in energy with decreasing particle radius. We measure a blueshift of approximately 0.9~eV (from $2.8$~eV to $3.7$~eV) when the particle radius decreases from 4 to 1~nm, which cannot be explained by our simulations within classical electrodynamics. We explain the absence of HO resonances in the EELS signal as a consequence of size-dependent damping in small metal particles,\cite{Kreibig:1969} an effect recently linked to nonlocal response in metals by accounting for diffusion in the free-electron gas.\cite{Mortensen:2014}
\begin{figure}[!t]
	\centering
	\includegraphics[scale=1]{fig1.eps}
	\caption{\textbf{Embedded silver nanoparticles.} Conceptual sketch of a silver (Ag) nanoparticle with radius $R$, which is embedded in a conformal silicon nitride (SiN$_\text{x}$) coating and excited by an electron beam. The impact parameter $b$ denotes the distance between the electron beam and the center of the nanoparticle. The thickness of the silicon nitride layer $t$ is determined from the EELS data and is approximately 35~nm, with about 20~nm below and 15~nm above the particles.}
	\label{fig:fig1}
\end{figure}

\section{Results}
The TEM samples containing encapsulated silver nanoparticles are prepared by first depositing a thin layer of silicon nitride (approximately $15$~nm) on a commercially available silicon nitride TEM membrane ($5$~nm thickness). Subsequently, silver nanoparticles are deposited through sputtering (see Methods), ensuring a large size-distribution, which, crucially, allows all measurements to be performed under identical conditions on the same sample. Lastly, a final thin layer of silicon nitride (again, approximately 15~nm in thickness) is deposited to completely encapsulate the silver particles. The sample is maintained in vacuum during the entire fabrication process. By inspection in scanning TEM (STEM), we find that the silver nanoparticles are multi-twinned icosahedral in shape. The thickness of the silicon nitride layers is a trade-off between signal-to-noise ratio in the EELS measurements and the increase in permittivity of the background environment of the nanoparticles. Here, we find that a thickness of 15~nm provides a good compromise between encapsulation and membrane-induced loss. EELS measurements on a sample with thicker silicon nitride layers (approximately $30$~nm each, i.e., a total thickness of $t=60$~nm) showed a significant noise level in the EELS signal due to increased energy losses of the electrons in the thicker silicon nitride coating (data not included).

The morphology of the encapsulating silicon nitride layer was investigated by high-angular annular dark field (HAADF) STEM imaging on a TEM lamella prepared from a sample with thicker silicon nitride coating using focused ion beam milling. Complementing the cross-sectional view, we have also used plane-view imaging in HAADF STEM combined with energy-dispersive x-ray spectroscopy (EDS) and EELS to determine the morphology and composition of the sample presented in this work. Both of the different analyses reveal that the silicon nitride coating is conformal, i.e., the coating conforms to the topology of the nanoparticles as illustrated schematically in Fig.~\ref{fig:fig1}. The conformal coating ensures that particles of all sizes are being completely encapsulated (see Supplementary Note 3 for more details).

By encapsulating the silver nanoparticles in silicon nitride, several experimental key advantages are obtained as compared to previous STEM EELS studies of silver nanoparticles on substrates.\cite{Ouyang:1992,Scholl:2012,Raza:2013,Scholl:2013,Kadkhodazadeh:2013} Firstly, the encapsulated system ensures that the nanoparticles will not move along the surface (or even detach from the surface) due to the Coulomb forces from the nearby electron beam (which can occur for small metal particles on substrates\cite{Batson:2008,Batson:2011}). Secondly, the formation of silver sulfide and silver oxide is prevented. Finally, as also mentioned earlier, the large dielectric constant of the encapsulating silicon nitride redshifts the dipole and HO modes into a region of significantly lower loss of silver, allowing the probing of different SP modes almost separately in our EELS setup.

The EELS measurements are performed in a monochromated and aberation-corrected FEI Titan operated in STEM mode at an acceleration voltage of 120~kV, providing a probe size of approximately $0.5$~nm and an energy resolution of $0.15\pm0.05$~eV. In Fig.~\ref{fig:fig2}(a) we show the EELS spectrum at fixed impact for an encapsulated silver nanoparticle of radius $R=20\pm 0.2$~nm, which is the largest particle studied in this work. The post-processing routine of the EELS data along with the determination of the particle radius including its uncertainty are detailed in the Methods section. The impact parameter $b$, which defines the distance from the center of the particle to the position of the electron beam, is $b=19.5$~nm, i.e., only $0.5$~nm from the inner surface of the particle. Figure~\ref{fig:fig2}(a) shows two clear peaks at 2.74~eV and 3.25~eV due to the excitation of the dipole and HO modes of the silver nanoparticle, respectively. The corresponding measured EELS intensity maps are shown in Fig.~\ref{fig:fig2}(b) and (c). By comparison of the two maps, it is evident that the HO modes exhibit a significantly larger degree of surface localization than the dipole mode. The different degrees of spatial localization in the EELS intensity maps is of crucial importance, as it affords a degree of control over the excitation of specific SP modes through the careful variation of the impact parameter, as we discuss in more detail in relation with Fig.~\ref{fig:fig3}.

Accompanying the EELS measurements in Fig.~\ref{fig:fig2}, we also present simulated EELS spectra which have been computed using the boundary-element method (BEM),\cite{Hohenester:2012,Abajo:2002} see the Methods section and Supplementary Note 1 for details on the implementation. The icosahedral shape of the nanoparticles are modelled as spheres, an approximation which is justified by the fact that the optical response of these two shapes is very similar.\cite{Noguez:2007,Yang:2011} We model the precise geometry of silver nanoparticles conformally encapsulated in finite-thickness silicon nitride layers, and use a permittivity of $\varepsilon_{\text{SiN}_x}=3.2$ for the silicon nitride. The Kramers--Kronig procedure\cite{Egerton:2011} did not give a reliable result for the dielectric function of the silicon nitride layer, see Supplementary Note 2 for details. We have therefore chosen a pragmatic approach to determining $\varepsilon_{\text{SiN}_x}$ by achieving the best correspondence between the simulated (classical electrodynamics) and measured dipole resonance energies for all of the particle sizes studied, except for the smallest particles ($R<4$~nm) where a blueshift of the dipole mode is observed. As we will see, this very same value for the permittivity of the silicon nitride ($\varepsilon_{\text{SiN}_x}=3.2$) also captures the HO resonance energies for all of the particle sizes examined. Additionally, the same value for the permittivity describes other EELS measurements made on \emph{gold} nanoparticles encapsulated in silicon nitride (data not included).

The red curve in Fig.~\ref{fig:fig2}(a) shows the result of the simulated EELS spectrum for the same particle radius ($R=20$~nm) and impact parameter ($b=19.5$~nm) as in the measurement. Clear spectrally separated peaks due to the excitation of the dipole, quadrupole, and octopole SP modes along with the bulk plasmon are present. We see that the dipole peak is captured in the EELS measurement, while the individual HO modes are not spectrally separated in the measurement. Instead the HO modes accumulate into a single broad peak due to the energy resolution of the EELS setup. To mimic the experimental energy resolution, we convolute our EELS simulation with a Lorentzian point-spread function (PSF) with a full-width at half maximum (FWHM) of 0.15~eV [grey curve in Fig.~\ref{fig:fig2}(a)]. In the convoluted case, the HO modes merge into a single broader peak with a resonance energy that is only slightly lower than in the EELS measurement. This small discrepancy may be due to the faceted nature of the nanoparticles, which may influence the resonance energy of the HO mode as the HO mode is more localized to the surface (and thereby more sensitive to the particular surface morphology) than the dipole mode. An additional feature of the convolution is the significant decrease in intensity of the bulk plasmon, which we also observe in the measured EELS spectrum of Fig.~\ref{fig:fig2}(a). Finally, Fig.~\ref{fig:fig2}(b) and (c) display theoretical EELS maps of the dipole and quadrupole SP modes, respectively, which confirms the experimentally observed stronger surface localization of the HO modes as compared to the dipole mode. In the following, we will discuss how we utilize this SP mode localization to probe the resonance energies of the dipole and HO modes separately. We start by considering the implications of the impact parameter on the EELS signal from a small silver nanoparticle (i.e., varying the impact parameter $b$ for a fixed particle radius $R$). Hereafter, we consider nanoparticles of different sizes with the electron beam close to the particle surface (i.e., varying the particle radius $R$ for almost fixed distance between impact parameters and particle surface $b-R$).
\begin{figure}[!b]
	\centering
	\includegraphics[scale=1]{fig2.eps}
	\caption{\textbf{EELS intensity maps of higher-order surface plasmons.} (\textbf{a}) Experimental EELS data as a function of energy loss of a SiN$_x$-encapsulated silver nanoparticle with radius $R=20 \pm 0.2$~nm excited by an electron beam positioned in the vicinity of the inner particle surface ($b=19.5 \pm 0.25$~nm). The red filled area shows the corresponding BEM EELS simulation, where the labels refer to dipole (DP), quadrupole (QP), octopole (OP), and bulk plasmon (BP). The gray filled area shows the theoretical curve (red filled area) convoluted with a Lorentzian point-spread function (PSF) with a FWHM of 0.15~eV ($f_\textsc{PSF}$). (\textbf{b-c}) Experimental (above), with zero-loss peak removed, and theoretical (below) EELS intensity maps of the same nanoparticle displaying the dipole and higher-order modes, respectively. The measured and simulated EELS intensity maps have identical color scales, respectively.}
	\label{fig:fig2}
\end{figure}

The relative excitation amplitudes of dipole, HO, and bulk plasmon modes depend strongly on the impact parameter $b$, i.e., the position of the electron beam when acquiring EELS data. This dependence is important to understand, since the resonance energies of the SP modes are spectrally close. The ideal situation is to find impact parameters where we primarily excite the dipole mode with a weak signal from the HO modes and vice versa. While challenging, this is to some extent possible: the dipole mode is primarily excited for the electron beam focused several nanometres outside the nanoparticle ($b>R$), while all SP modes including the HO modes are efficiently excited for the electron beam close to the surface of the particle ($b\approx R$). For electron-beam excitations inside the particle ($b<R$), the bulk plasmon is also excited. In Fig.~\ref{fig:fig3} we show experimental and simulated EELS spectra of an encapsulated silver nanoparticle with a radius of $R=9.2\pm0.2$~nm as a function of the impact parameter $b$. Figures~\ref{fig:fig3}(a) and (b) show a schematic and an experimental STEM image of the nanoparticle, respectively, along with the positions of the electron beam. The corresponding theoretical and experimental EELS spectra are shown in Fig.~\ref{fig:fig3}(c) and (d). The impact parameter is varied systematically through a radial progression from particle center to several nanometers outside its surface. In Fig.~\ref{fig:fig3}(c), the unconvoluted theoretical EELS spectra are shown as line plots, while the PSF-convoluted spectra are shown as filled area plots. Considering first the spectra for the largest impact parameters (red curves), it is clear from Fig.~\ref{fig:fig3}(c) that the simulated EELS signal primarily stems from the dipole mode. In the EELS measurement a single symmetric peak is observed, which we accordingly interpret as the dipole mode. A Gaussian function [black solid line in top panel of Fig.~\ref{fig:fig3}(d)] was fitted to determine its resonance energy. As the electron beam approaches the surface of the nanoparticle, the simulated spectra show that HO modes increasingly contribute to the EELS signal (purple and blue curves). This feature is also observed in our EELS measurements (purple curve), where we observe that the peak is no longer symmetric due to the EELS signal from the HO modes. For $b=7.2$~nm (blue curve) an additional peak (besides the dipole peak) shows up. By fitting the sum of two Gaussian functions [grey dashed lines in third panel from top in Fig.~\ref{fig:fig3}(d)] we determine the resonance energies of both of the peaks as the center of each Gaussian. The resonance energy of the low-energy peak is in good agreement with that measured for the dipole resonance energy for impact parameter outside the particle ($b>R$), and the high-energy peak is then attributed to the accumulated contributions from HO modes. As the electron beam penetrates closer to the center of the particle, the simulated EELS signal from the bulk plasmon increases significantly while the signal from the SP modes decreases. In particular the dipole mode decreases strongly in signal since the penetrating electron beam cannot properly induce the antisymmetric dipolar charge distribution in the particle. Briefly summing up, we have shown that the dipole mode can be probed selectively by positioning the electron beam outside the particle ($b>R$), while the signal from the HO modes can be maximized by considering impact parameters close to the surface of the particle ($b\approx R$).
\begin{figure}[!b]
	\centering
	\includegraphics[scale=1]{fig3.eps}
	\caption{\textbf{EELS versus impact parameter.} (\textbf{a-b}) Theoretical setup and experimental STEM image, depicting also the impact parameters considered in the simulated and measured EELS spectra of (\textbf{c-d}), respectively, acquired from a silver nanoparticle with radius $R=9.2 \pm 0.2$~nm. The theoretical spectra (line plots) have been convoluted with a Lorentzian PSF with a FWHM of 0.15~eV (area plots). The convoluted spectra are normalized to unity area in the considered energy range. For clarity, the unconvoluted and convoluted spectra have the same maximal value. The measured spectra have been analyzed according to the procedure described in the Methods section and afterwards normalized to unity area in the same energy range. The fit of a single or two Gaussian functions to the peaks of interest in the experimental spectra is indicated in solid and dashed gray lines. In the measurements, the STEM probe diameter of $0.5$~nm is used as an estimate for the error in the impact parameter $b$.}
	\label{fig:fig3}
\end{figure}

Next, we study the dipole and HO mode resonance energies of encapsulated silver nanoparticles as a function of particle radius $R$, while keeping the electron beam positioned close to the particle surface. To determine the \emph{dipole} resonance energy, we consider only impact parameters several nanometers outside the nanoparticle. In this case, the EELS spectra show a single symmetric peak [as that shown in the red curve of Fig.~\ref{fig:fig3}(d)], which we fit with a single Gaussian function to determine the resonance energy. However, the resonance energy of the \emph{HO modes} is determined by considering EELS spectra acquired by positioning the electron beam close to the surface of the nanoparticle. In these measurements as shown in Fig.~\ref{fig:fig4}(a-g), we consistently observe two peaks on silver nanoparticles of various radii (for fixed distance between electron beam and particle surface $b-R$), in agreement with what we found in the blue spectrum in Fig.~\ref{fig:fig3}(d) where the electron beam was also close to the surface. We stress that the clear and repeated observations in Fig.~\ref{fig:fig4}(a-g) of two SP peaks rather than one in various nanoparticles, and the low noise level, are aided by the fact that the strongest EELS signal is observed when exciting particles close to their surfaces. While individual silver nanoparticles have been studied before with EELS,\cite{Ouyang:1992,Scholl:2012,Raza:2013} the measurements presented here constitute, to our knowledge, the first unambiguous observation of HO modes on this ultra-small size and their controlled excitation.

For every panel in Fig.~\ref{fig:fig4}(a-g) we determine the resonance energy of the HO modes by fitting the spectrum to the sum of two Gaussians (one for the dipole peak and one for the HO peak) and extracting the peak position of the Gaussian function fitted to the HO peak. The two-Gaussian fits accurately capture the features of the EELS data (see also Supplementary Figure~S1 for a residual plot of Fig.~\ref{fig:fig4}), which provides support that the HO modes peak can to a good approximation be described by a single (symmetric) Gaussian function. The reason for this is two-fold. Firstly, the individual HO modes are only closely spaced in energy (with decreasing energy spacing for increasing angular momentum $l$) as a consequence of the maximal allowed resonance energy of the HO modes in classical electrodynamic theory, given by the $l\rightarrow \infty$ limit of Eq.~\eqref{eq:res_freq_nr}. Secondly, the experimental PSF is a spectrally broad (compared to the individual HO modes) and symmetric function, which merges the individual spectrally-close HO modes into a single almost symmetric peak [see also Fig.~\ref{fig:fig2}(a)]. We can provide an accurate estimate for the energy of the $l\rightarrow \infty$ limit, which provides an upper energy limit for the HO modes, by considering the classical $l$-dependent resonance condition for a spherical particle in a homogeneous background medium, governed by\cite{Christensen:2014}
\begin{equation}
	l\varepsilon(\omega) + (l+1) \varepsilon_\textsc{b} = 0, \label{eq:res_cond}
\end{equation}
from which Eq.~\eqref{eq:res_freq_nr} is derived. The nonretarded approximation used in the derivation of Eq.~\eqref{eq:res_cond} is valid and accurate for silver nanoparticles with radius below 10~nm. To determine the energy of the $l\rightarrow\infty$ limit of Eq.~\eqref{eq:res_cond}, we require first an effective value for the background permittivity $\varepsilon_\textsc{b}$. We estimate the effective background permittivity by rewriting Eq.~\eqref{eq:res_cond} in the case of the dipole mode ($l=1$) as $\varepsilon_\textsc{b} = -\textrm{Re}[\varepsilon(\omega=\omega_\textsc{dp})] / 2$, where $\hbar\omega_\textsc{dp}$ is the resonance energy of the dipole mode. Here, we use the average value for the dipole resonance energies measured from nanoparticles in the radius range $4~\textrm{nm}<R<10~\textrm{nm}$ (Fig.~\ref{fig:fig5}), where the lower radius limit is needed since we measure an increase in the resonance energy for particles smaller than $R<4$~nm (we discuss this observation in more detail in relation to Fig.~\ref{fig:fig5}). From this procedure we find $\varepsilon_\textsc{b} = 3.3$, which is very close to $\varepsilon_{\text{SiN}_x}=3.2$ used for the silicon nitride layers in our BEM simulations, supporting also the interpretation that the smallest nanoparticles ($R<10$~nm) behave as if in a fully embedded homogeneous background environment. With $\varepsilon_\textsc{b} = 3.3$ as the value for the background permittivity, we find from Eq.~\eqref{eq:res_cond} the energy for the $l\rightarrow\infty$ limit to be $\hbar \omega_{l \rightarrow \infty} = 3.27$~eV, which, in the framework of classical electrodynamic theory, provides an upper limit for the resonance energy of the HO modes. Interestingly, the experimental EELS spectra shown in Fig.~\ref{fig:fig4}(a-g) (and the corresponding resonance energies shown in Fig.~\ref{fig:fig5}) indicate that the HO modes are bounded by this classical upper energy limit (shown as dotted lines in Figs.~\ref{fig:fig4} and \ref{fig:fig5}), since we do not measure any of the HO peaks to have a resonance energy larger than $\hbar \omega_{l \rightarrow \infty}$, although we do find some of the resonance energies to be very close to the limit.

As a quantitative measure of the strength of the HO modes, we additionally determine the relative spectral weight of the HO modes $S_\textsc{ho}$ as the ratio of the areas of the HO Gaussian fit to the dipole Gaussian fit. The values for $S_\textsc{ho}$ are shown in each panel in Fig.~\ref{fig:fig4}(a-g). We see that the relative spectral weight of the HO modes often exceeds 50\% and in one instance even surpasses 100\% [Fig.~\ref{fig:fig4}(b)], indicating that the EELS signal from the HO modes dominates that of the dipole mode. These large values for $S_\textsc{ho}$ suggest that EELS data acquired by positioning the electron beam at the surface of the nanoparticle cannot simply be interpreted solely as the dipole mode,\cite{Ouyang:1992,Scholl:2012,Raza:2013} but require careful consideration of the HO modes.
\begin{figure}[!b]
	\centering
	\includegraphics[scale=1]{fig4.eps}
	\caption{\textbf{EELS versus particle sizes.} (\textbf{a}-\textbf{g}) EELS spectra of silver nanoparticles encapsulated in silicon nitride used for determining the resonance energy of the HO surface plasmon modes. The least-squared fit of the sum of two Gaussian functions is indicated in solid gray lines, with each Gaussian indicated in dashed gray lines. The parameter $S_\textsc{ho}$ denotes the spectral weight of the HO modes, calculated as the ratio of the areas of the HO and dipole Gaussian functions. The dotted lines indicate the classical $l\rightarrow\infty$ energy limit. (\textbf{h}) EELS spectrum of encapsulated small silver nanoparticle (typical for all spectra for $R < 4$~nm) showing only a single blueshifted plasmon resonance, attributed to the dipole mode. The text inset of each plot designates the particle radius $R$ and the impact parameter $b$.}
	\label{fig:fig4}
\end{figure}

The experimentally-observed dipole and HO mode resonance energies of all of the nanoparticles studied in this work are summarized in Fig.~\ref{fig:fig5}. Additionally, the superimposed color plot shows the simulated EELS spectra, where an impact parameter set to 1~nm from the inner particle surface is chosen (i.e., $b-R=-1$~nm) to strongly excite the HO modes (as also done in the EELS measurements). The simulated EELS spectra have been convoluted with a Lorentzian PSF and show a strong low-energy peak due to the dipole mode and an only slightly weaker high-energy peak due to the excitation of HO modes. Furthermore, an EELS signal from the bulk plasmon with an energy of approximately 3.8~eV is also present and increases in strength for decreasing particle size (since the impact parameter gets closer to the center of the particle). In the simulations, the dipole mode redshifts slightly with increasing particle radius due to retardation effects, while the HO modes remain almost constant in energy. For particle radii above 4~nm, we find excellent agreement between the observed and simulated resonance energies for both the dipole and HO modes, thereby remarkably validating classical electromagnetic theory down to the few-nanometer scale.

For particle radii below 4~nm, we experimentally observe two distinct features that are not explained by our classical simulations. Firstly, the EELS signal from the HO mode decreases significantly (and basically disappears), since we only observe one single peak in all of our EELS measurements regardless of the impact parameter. We interpret this peak as the dipole mode, since a strong EELS signal is always present for impact parameters outside the nanoparticles, see Fig.~\ref{fig:fig4}(h) for a typical example of such an EELS spectrum. Figure~\ref{fig:fig5} shows that the EELS simulation based on classical electrodynamics does not reproduce this effect. Interestingly, the disappearance of the HO mode can be explained by taking into account nonlocal response in the silver nanoparticles,\cite{Mortensen:2014} which gives rise to size-dependent damping of SPs,\cite{Kreibig:1969} for one thing. Nonlocal response is a consequence of accounting for the wave-vector dependence of the dielectric function (i.e., spatial dispersion), which, in a hydrodynamic approach,\cite{Raza:2015b} originates from convective and diffusive currents of the free-electron gas in metals.\cite{Mortensen:2014} In Supplementary Note 1 we explain in detail the absence of the HO resonance in the EELS spectra of the smallest particles, based on two assumptions: First, the silver nanoparticles are fully embedded in a homogeneous background environment, a valid assumption for the smallest nanoparticles ($R<10$~nm). The second assumption is that the non-classical behavior of the nanoparticles is well described by the generalized nonlocal optical response (GNOR) model,\cite{Mortensen:2014} which shows that size-dependent damping decreases the EELS signal from the HO modes significantly for particle radii smaller than approximately 4~nm. The size-dependent damping is theoretically expected to increase with multipolar order $l$ and is therefore anticipated to influence the HO modes more than the dipole mode.\cite{Raza:2015b} The EELS simulations based on the GNOR model accurately explains our experimental observation of the disappearance of HO modes for $R<4$~nm, which provides strong support for the interpretation that this effect is due to nonlocal response, i.e., a fundamental physical property of silver nanoparticles, and not, e.g., due to the energy resolution of EELS or other instrument-related effects.
\begin{figure}[!t]
	\centering
	\includegraphics[scale=1]{fig5.eps}
	\caption{\textbf{Dipole and higher-order mode resonance energies.} Measured resonance energies as a function of particle radius for the dipole (green dots) and higher-order (light blue dots) surface plasmon modes in silver nanoparticles encapsulated in silicon nitride. The color plot shows the PSF-convoluted EELS simulations. For each particle size in the simulations, the electron beam is positioned 1~nm from the inner particle surface (i.e., $b-R=-1$~nm). Each simulated EELS spectrum has been normalized to unity area.}
	\label{fig:fig5}
\end{figure}

The second measured feature is a huge, nearly discontinuous, increase in dipole resonance energy of approximately 0.9~eV when the particle radius decreases from 4~nm to 1~nm (see Fig.~\ref{fig:fig5}), which is not reproduced in our classical simulations, see also Supplementary Note 1 for additional simulations with a different impact parameter. We are unable to explain the observed blueshift (shown in Figure 5) as a consequence of the position of the electron beam, since the measured resonance energies are still larger than the HO mode energy but lower than the bulk plasmon energy. Here, we note that the blueshift of the SP in similar-sized silver nanoparticles has been observed before.\cite{Kreibig:1985,Kreibig:1995,Ouyang:1992,Tiggesbaumker:1993,Charle:1998,Scholl:2012,Raza:2013} Previous EELS measurements of silver nanoparticles on thin ($<10$~nm) TEM substrates\cite{Scholl:2012,Raza:2013} have shown a (continuous) blueshift of the SP of approximately 0.5~eV in the same size range, which suggests that the non-classical blueshift depends on the effective background permittivity, with larger permittivity giving rise to larger blueshift. However, these previous EELS measurements were performed by acquiring spectra with impact parameters at the particle surface, suggesting, in the light of the results shown in this work, that HO modes may have come into play as well, thereby smoothening the blueshift of the SP to appear more continuous. Calculations based on the nonlocal GNOR model provide qualitative agreement with the EELS measurements shown in Fig.~\ref{fig:fig5}, however, the experimentally-measured blueshift is larger than predicted by this nonlocal theory (see Supplementary Note 1). The origin of the blueshift has been attributed to several other effects besides nonlocal response,\cite{Ruppin:1973,Raza:2013,Mortensen:2014} including quantum-size effects,\cite{Scholl:2012,Genzel:1975} screening from the \textit{d}-band electrons in silver,\cite{Liebsch:1993} and the combination of the latter two.\cite{Monreal:2013} In this work we have shown that nonlocal response indeed affects the SP modes of silver nanoparticles, but as the particle radii decreases, it becomes inherently more difficult to distinguish the relative importance of each of these effects as they all contribute non-negligibly on this size scale. This suggests that the origin of the observed blueshift likely derives from a concerto of these effects.

\section{Discussion and conclusion}
Using EELS, we have experimentally studied the plasmonic response of silver nanoparticles with radii $1-20$~nm encapsulated in a conformal silicon nitride coating. In our EELS spectra, we find, depending on the impact parameter, up to two SP-related resonance peaks due to the excitation of the dipole and higher-order (HO) surface plasmon modes, providing the first unambiguous observation of HO modes in nanoparticles with radius down to only $4$~nm. We find that impact parameters close to the surface of the nanoparticles exhibit the strongest EELS signal from the HO modes, while impact parameters several nanometers outside of the surface mainly probe the dipole mode. For impact parameters close to the inner surface of  nanoparticles with certain radii, the contribution of HO modes to the EELS signal even exceeds that of the dipole mode. Comparison with EELS simulations based on classical electrodynamic theory provides excellent agreement down to particle radii of only $4$~nm. The continued agreement even down to these size-scales is an impressive testament to the qualities of classical electrodynamic theory.

The strong presence of HO modes in small silver nanoparticles has profound influence on the interpretation of previous EELS measurements on nanoparticles on substrates.\cite{Ouyang:1992,Scholl:2012,Raza:2013} In particular, we have provided thorough experimental evidence that the EELS signal depends strongly on the impact parameter. Positioning the electron beam as close as possible to the nanoparticle surface will excite the HO modes as well as the dipole mode, so that the single or second observed surface-plasmon peaks in the EELS spectra actually are to be understood as the combined effect of many multipole surface plasmon resonances. To ensure that the EELS signal stems only from the dipole mode, we propose to position the electron beam several nanometers from the particle surface. Hereby, signal from HO modes and the bulk plasmon may be minimized.

For particle radii below 4~nm, we observe the disappearance of the HO modes and an increase in dipole resonance energy with decreasing particle size. The former is a consequence of size- and $l$-dependent surface plasmon damping, which are both effects due to nonlocal response and thereby a fundamental physical property of silver nanoparticles. Regardless of the impact parameter, we find only a single resonance for radii below 4~nm, which increases in energy with decreasing particle size. The lack of dependence on impact parameter suggests that these very small particles no longer support separate bulk and surface modes, but only a single plasmonic mode. We interpret the nature of this mode as a surface dipole mode, since the mode can be excited with impact parameters outside the particle, a feature not possible for bulk plasmons. The dipolar nature of the plasmon mode is further substantiated by earlier optical far-field measurements on similar-sized nanoparticles,\cite{Kreibig:1995,Tiggesbaumker:1993,Charle:1998} which only probe the dipole mode.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Expeditionary science is the act of collecting \emph{in situ} samples and observations of the natural world for the purposes of scientific discovery and model development. Transient, dynamic phenomena --- hydrothermal plumes, algal blooms, lava flows --- are of interest in many disciplines of observational science. It is often impossible for a single, static sensor to capture a comprehensive picture of these mesoscale spatiotemporal phenomena due to their large spatial scales and dynamic temporal nature. Large networks of static sensor nodes could address this sensing challenge, but the transience and sparsity of the target phenomena necessitates logistically impractical sensor density. Mobile robots are uniquely well-positioned to tackle this problem via extended deployments carrying heterogeneous sensor payloads that actively seek and perform targeted surveys of ephemeral environments. However, robots deployed for expeditionary science today currently execute open-loop, preset trajectories, like ``lawnmowers'' (back-and-forth grid-based pattern) hand-designed by human scientists. In dynamic environments, this open-loop execution often results in sparse measurements of the target phenomena or misses a short-lived target entirely. Given the cost of expeditionary field operations and the value of the data collected, it is critical to improve the efficiency and efficacy of robots as scientific tools.

Building robotic platforms for expeditionary science requires an autonomy stack that can integrate observations from heterogeneous sensors into a model of a spatiotemporal system and use this model to plan informative trajectories that target a specific scientific objective. This kind of autonomy poses many challenges for integrating probabilistic modeling and decision-making. Previous work in informative path planning (IPP) \cite{Hitz2017}, adaptive sampling/experimental design \cite{krause2008near}, and decision-making under uncertainty \cite{sunberg2018online} has tackled aspects of the expeditionary science problem, especially in static environments using data-driven models and information-based rewards. However, existing methodologies do not trivially extend to spatiotemporal environments, leaving key challenges such as model and dynamics learning and decision-making in large-scale environments unaddressed. In this abstract, we discuss these and other key challenges, and present some preliminary work conducted for hydrothermal plume discovery and mapping in the deep sea.

\section{Expeditionary Science as a Decision-Making Problem}
In an expeditionary science mission, a robot is tasked with collecting scientifically useful measurements of an unknown, partially-observable spatiotemporal environment under time, energy, and dynamical constraints. We formulate this sequential decision-making problem as a partially observable Markov decision-process (POMDP). Let $\Pi(\cdot)$ denote the space of probability distributions over the argument. A finite horizon POMDP can be represented as tuple: $(\mathcal{S}, \mathcal{A}, T, R, \mathcal{Z}, O, b_0, H, \gamma)$, where $\mathcal{S}$ are the states, $\mathcal{A}$ are the actions, and $\mathcal{Z}$ are the observations. At planning iteration $t$, the agent selects an action $a \in \mathcal{A}$ and the transition function $T: \mathcal{S} \times \mathcal{A} \to \Pi(\mathcal{S})$ defines the probability of transitioning between states in the world, given the current state $s$ and control action $a$. After the state transition, the agent receives an observation according to the observation function $O: \mathcal{S} \times \mathcal{A} \to \Pi(\mathcal{Z})$, which defines the probability of receiving an observation, given the current state $s$ and previous control action $a$. The reward function $R: \mathcal{S} \times \mathcal{A} \to \mathbb{R}$ serves as a specification of the task. A POMDP is initialized with belief $b_0$ and plans over horizon $H$ with discount factor $\gamma$.

Due to the stochastic, partially observable nature of current and future states, the realized reward in a POMDP is a random variable. Optimal planning is defined as finding a policy $\{\pi_t^*: \Pi(\mathcal{S}) \to \mathcal{A}\}_{t=0}^{H-1}$ that maximizes expected reward: $\mathbb{E} \Big[ \sum_{t=0}^{H-1} \gamma^t R\big(S_t, \pi_t(b_t)\big) \mid b_0 \Big]$, where $b_t$ is the updated belief at time $t$, conditioned on the history of actions and observations. The recursively defined horizon-$h$ optimal value function $V^*_h$ quantifies, for any belief $b$, the expected cumulative reward of following an optimal policy over the remaining planning iterations: $V_0^{*}(b) = \truemax_{a \in \mathcal{A}} \mathbb{E}_{s \sim b}[R(s, a)]$ and
\begin{align}
     V_h^{*}(b) &=  \max_{a \in \mathcal{A}} \mathbb{E}_{s \sim b}[R(s, a)] + \gamma \int_{z \in \mathcal{Z}} P(z \mid b, a) V_{h-1}^*(b^{a,z}) \text{d}z \hspace{0.6cm} h \in [1, H-1],  \label{eq:value}
\end{align}
where $b^{a,z}$ is the updated belief after taking control action $a$ and receiving observation $z$, computed via Bayes rule using the transition $T$ and observation $O$ functions. The optimal policy at horizon $h$ is to act greedily according to a one-step look ahead of the horizon-$h$ value function. However, Eq.~\ref{eq:value} is intractable for large or continuous state, action, or observation spaces and thus the optimal policy must be approximated. Well-designed algorithmic and heuristic choices are necessary to develop efficient, practical, and robust planning algorithms for expeditionary sciences.

\section{Advances to Expeditionary Science}
\label{sec:fieldops}
As a case study we highlight a scientific mission to the Gulf of California in November 2021, in which a research vessel deployed an autonomous underwater vehicle (AUV) for deep sea discovery and mapping of hydrothermal vent plumes (see Fig.~\ref{fig:diagram}). Characterizing hydrothermalism has implications for global nutrient and energy budgets, but nearly two-thirds of all hypothesized hydrothermal vents are yet undiscovered \cite{beaulieu2015undiscovered}. AUVs with point sensors and the ability to execute trajectories under operational constraints can ``sniff'' for vents in the water column and use observations to estimate seafloor vent characteristics, but tidal advection, diffusion, and turbulent mixing all contribute to complicated spatiotemporal observations. Additionally, accurate vent inference requires the disambiguation of the spatially-small, chemically-rich buoyant stem from the spatially-vast neutrally-buoyant layer that compose plumes \cite{morton1956turbulent}. Here, we present some of our work tackling this expeditionary robotics problem, and highlight key challenges we found in this context.

\begin{figure}[h!]
    \centering
    \includegraphics[width=\columnwidth]{figures/plume_model_concept_2.png}
    \caption{In November 2021, AUV Sentry was deployed in the Gulf of California to detect and map hydrothermal plume structures in the deep ocean. Plumes are composed of a buoyant stem and neutrally-buoyant layer. These plumes are dynamic, subject to tidal advection, diffusion, and turbulent mixing. We infer the spatiotemporal structure of plumes from sparse partial-observations and plan informative trajectories that respect operational constraints of AUV Sentry. Field data are shown on the right from two onboard science point-sensors. The dORP/dt sensor shows areas of reactive water, likely indicating buoyant stem detections. The OBS sensor shows areas of high turbidity, which can persist throughout a plume. Combining these sensor signals, we can infer the plume structure and dynamics.}
    \label{fig:diagram}
\end{figure}

\subsection{Spatiotemporal IPP for the Maximum Seek-and-Sample Problem}
We formulated the problem of finding the buoyant stem within the multi-modal non-buoyant plume layer using a robot as a POMDP and developed \textsc{Plumes }: \textbf{P}lume \textbf{L}ocalization under \textbf{U}ncertainty using \textbf{M}aximum-valu\textbf{E} information and \textbf{S}earch \cite{flaspohler2019information}, an autonomy stack for target-seeking in multi-modal environments. \textsc{Plumes } utilized a Gaussian Process (GP) belief to model a single time-varying chemical signal, a progressively-widened Monte Carlo tree search, and a novel task-driven information-theoretic reward function. In spatiotemporal systems, \textsc{Plumes } was able to find and track moving buoyant-stem targets and repeatedly explore an environment as uncertainty grew under stochastic plume dynamics \cite{preston2019adaptive}. \textsc{Plumes } addressed challenges in continuous search for planning and balancing exploration-exploitation behavior under sparse reward signals, and demonstrated how the sum of algorithmic components can lead to a potent system. However, \textsc{Plumes } also revealed the challenges of expeditionary science, particularly the difficulty of simultaneously training and utilizing a spatiotemporal belief model, and developing computationally tractable real-time nonmyopic planners.

\subsection{Abstraction via Task-driven Macro-action Discovery}
One key barrier to tractable, real-time nonmyopic planning in \textsc{Plumes } was the need to search over long planning horizons, chaining together low-level action primitives to explore the environment. The development of abstract, high-level action representations would enable planners to search over shorter planning horizons, saving compute effort during real-time planning. Instead of using human-designed abstractions of the robot's continuous action space, our recent work learns a set of high-level actions that are useful for a given scientific task directly from a low-level specification of the POMDP model \cite{flaspohler2020belief}. The algorithm takes into account the robot's current uncertainty about the state of the world, the stochasticity of the underlying environmental dynamics, and the current task to develop a set of high-level actions that can be used online while providing formal performance guarantees about the resulting ``abstract'' policies.

\subsection{Iterative Mission Structures and Physically-Informed Uncertainty Models}
\textsc{Plumes } used a general, data-driven model of the environment and assumed that the robot's action space consisted of flexible, low-level motion primitives. To improve the operational effectiveness of \textsc{Plumes } at sea, we proposed to: 1) make use of data-efficient probabilistic models that explicitly use physical equations that govern plume expression in the water column, and 2) develop planners that directly encode the physical and operational constraints of the AUV. We developed and tested this new algorithmic expeditionary system during field trials in November 2021. The system made use of a belief model called \textsc{Phumes }: \textbf{PH}ysically-informed \textbf{U}ncertainty \textbf{M}odel of \textbf{E}nvironment \textbf{S}patiotemporality, which sampled uncertain parameters of a simplified model of buoyancy-driven hydrothermalism (e.g., vent temperature, tidal crossflow) to generate multiple ``envelopes'' of possible plume-derived waters. These were aggregated into a probabilistic forecast over the spatiotemporal volume to be sampled by the AUV. A trajectory optimizer with knowledge of vehicle and ship operational constraints then scaled lawnmower primitives to maximize the number of ``in plume'' measurements that were to be collected at the continuum between buoyant stem and neutrally-buoyant layer. \textsc{Phumes } addressed open challenges in robotic decision-making related to injection of prior knowledge or scientific inductive bias for data efficient training in complex domains and scalable planning with expensive belief models.

\section{Core Challenges and Opportunities in Expeditionary Robotics}
\label{sec:challenges}
As identified in Sec.~\ref{sec:fieldops}, there are a number of open challenges in reinforcement learning and decision-making in the context of expeditionary robotics. Here we present a brief survey of these open challenges and opportunities informed by years of collaborating with scientists to develop expedition frameworks in diverse application areas.

\subsection{Belief Representations}

\textbf{Epistemic and aleatoric uncertainty:}
Reducing epistemic uncertainty of a spatiotemporal environment requires access to a model of the underlying dynamical system, or a data-driven technique that can uncover it. Extracting physically-meaningful quantities from observational data is typically performed post-expedition using computationally expensive numerical models ``tuned'' by observations. While this lends itself well to Bayesian inference formulations, it is intractable for practical decision-making. Data-driven techniques for model discovery \cite{raissi2019physics} may be arguably more tractable, but generally suffer small-data challenges. Developing models that overcome the challenges of efficiently characterizing spatiotemporal dynamics from streaming, sparse observations would generally improve expeditionary robotics. Additionally, there is a unique opportunity to enable computation of proxies for aleatoric uncertainty, which are well-described in spatiotemporal environments with measures of chaotic motion (e.g., Lyapunov exponents) inferred from data \cite{blanchard2019analytical}. The implication that aleatoric uncertainty can be estimated has yet to be utilized to, e.g., assess the attainable resolution of a model or set planning horizons.

\textbf{Scientific knowledge as inductive bias:}
The kernel of a GP, the loss function in a neural network, or the activation functions between layers in a deep network can all be viewed as forms of inductive bias in a learning problem. For data-driven discovery of spatiotemporal dynamics, improving sample efficiency by leveraging opportunities to inject scientific knowledge to alleviate the learning burden is an open problem. While canonical numerical models of spatiotemporal phenomena are too computationally expensive to directly incorporate into e.g., GP kernels, the physical principles that underlie these models can be more easily summarized. ``Physically-informed'' data-driven probabilistic representations have been demonstrated outside of expeditionary robotics \cite{raissi2019physics} and some work within IPP \cite{salam2019adaptive} shows rich opportunities for analyzing and extending these methods for larger environments and longer planning horizons.

\textbf{Low-dimensional state embeddings:}
Expressing a spatiotemporal environment completely would require an exceedingly large, high-dimensional representation. Model order reduction (MOR) techniques reduce the dimensionality of spatiotemporal systems to a set of weights and vectors that sufficiently describe patterns in the dynamics. Uncovering low-dimensional state embeddings from partially-observed expedition data is a general challenge; uncovering a \emph{useful} embedding for a specific decision-making problem is additionally challenging. Access to such an embedding would reduce the computational burden of representing belief in large environments for planning.

\subsection{Decision-Making}

\textbf{Rollout-based planning with expensive belief models:}
State-of-the-art planners for POMDP problems often make use of rollout-based planning in tree search frameworks; continuous search variables are handled using strategies such as progressive widening or scenario sampling \cite{sunberg2018online}. However, these planners require extensive online simulations for each rollout performed. Forward-simulating the dynamics and observational models for complex, spatiotemporal phenomena can be computationally intensive, which often limits the feasible look-ahead horizon in real-time operations on computationally-limited robotic platforms.  Planners that selectively or adaptively perform expensive rollouts, automatically adjust the planning horizon based on the dynamics of the environmental system, or make use of continuous, offline planners would enable improved decision-making for expeditionary science.

\textbf{Information rewards and task-driven exploration:} Due to partial observability and stochastic dynamics in spatiotemporal contexts, a decision-maker must operate with significant and often growing state uncertainty. However, not all state uncertainty impacts task performance and uniform information gathering strategies can be inefficient. Understanding the value of information for accomplishing a task is a known challenge for planning under uncertainty and this is particularly true for expeditionary robotics. Recent works that develop heuristic information rewards \cite{flaspohler2019information} or task-driven value of information metrics \cite{flaspohler2020belief} begin to build the tools necessary for expeditionary robotic planning.

\textbf{Robust planning under model mismatch and uncertainty:}
Scientific models, whether data-driven or based on physical principles, are always imperfect representations of a robot's environment. Model mismatch or uncertainty in key model parameters leads to discrepancies between the environmental predictions that a robot uses during planning and its real-time observations. Planning robot trajectories that entirely miss a phenomenon due to overconfidence in an incorrect model is detrimental to scientific objectives. Planners must develop policies or trajectories that are robust to model mismatch and uncertainty, or are guaranteed to perform as well as a simple, naive data collection strategy.

\textbf{Interpretable and operational decision-making:}
Decision-making algorithms must interface with and are constrained by a variety of stakeholders, including scientists, robot operators, and engineers. For example, when deploying an AUV from an oceanographic research vessel, the decision-making algorithm must account for ship scheduling, timing delays, weather, and multi-vehicle operations. This requires developing flexible planners that can understand and account for these complex constraints. Additionally, stakeholders are often concerned with robot safety and data quality. Producing plans that are interpretable for scientific and operational stakeholders is key for building trust and confidence in scientific autonomy.

\section{Final Thoughts}
Expeditionary science and robotics motivates a set of interesting reinforcement learning and decision-making problems that are not well-addressed by current state-of-the-art methods. There is increasing need to enable better \emph{in situ} methodologies for monitoring and characterizing large spatiotemporal environments as accelerated climate change transforms delicate carbon budgets, ecosystem networks, and weather patterns. Robotic technologies can play a transformative role in understanding these evolving trends and assessing policy interventions. Novel developments towards the challenges listed here may have constructive and complementary effects in other application spaces, as many domains, such as material discovery and robotic manipulation, have overlapping challenges and opportunities. Our contributions, including \textsc{Plumes } \cite{flaspohler2019information,preston2019adaptive}, macro-action discovery \cite{flaspohler2020belief}, the \textsc{Phumes } model and trajectory optimizer for operational missions, and ongoing work in physically-informed deep kernel learning, represent algorithmic and systems contributions towards more effective autonomous systems in expeditionary science and begin to address some of the key challenges that we propose.

\bibliographystyle{IEEEtran}


text:
\section*{Acknowledgments} located above the reference section.
\thanks{The first and last authors are supported by [NIHR01 EB02875 and NSF MDS-2010778]. The second author is supported by KAUST.}

\thanks{$^*$Corresponding author: Moo K. Chung}

\begin{document}
\maketitle

\centerline{\scshape
Soumya Das$^{{\href{mailto:sundew122436@gmail.com}{\textrm{\Letter}}}*1}$
and Hernando Ombao$^{{\href{mailto:hernando.ombao@kaust.edu.sa}{\textrm{\Letter}}}2}$
and Moo K. Chung$^{{\href{mailto:mkchung@wisc.edu}{\textrm{\Letter}}}1}$
}

\medskip

{\footnotesize
 \centerline{$^1$University of Wisconsin-Madison, USA}
}

\medskip

{\footnotesize

 \centerline{$^2$King Abdullah University of Science and Technology, Saudi Arabia}
}

\bigskip

\begin{abstract}
Developing reliable methods to discriminate different transient brain states that change over time is a key neuroscientific challenge in brain imaging studies. Topological data analysis (TDA), a novel framework based on algebraic topology, can handle such a challenge. However, existing TDA has been somewhat limited to capturing the static summary of dynamically changing brain networks. We propose a novel dynamic-TDA framework that builds persistent homology over a time series of brain networks. We construct a Wasserstein distance based inference procedure to discriminate between time series of networks. The method is applied to the resting-state functional magnetic resonance images of human brain. We demonstrate that the proposed dynamic-TDA approach can distinctly discriminate topological patterns between male and female brain networks.
\end{abstract}

\section{Introduction}

The human brain and its dynamical changes is an archetypal example of complex systems. Within this system, brain functions depend on constant interplays between local information processing and efficient global integration of information \cite{sporns.2005, bullmore.2009, siemon.2014}. Analyzing such interactions in brain pathology is crucial to reveal the underpinnings of neurodegenerative diseases \cite{uhlhaas.2006, uhlhaas.2010, hu.2015}. In brain network analysis, the neural interactions are encoded as a graph consisting of nodes and edges. Often, the whole brain is parcellated into several disjoint regions, which serve network nodes \cite{tzourio.2002, desikan.2006, hagmann.2007, fornito.2016, arslan.2018} whereas the correlations between different parcellations serve as edge weights. In most graph theory based analysis, features such as node degrees and clustering coefficients are acquired from adjacency matrices after thresholding edge weights \cite{sporns.2003, vanwijk.2010}. However, this leads the final statistical results and interpretations that rely on the choice of the threshold \cite{lee.2012.TMI, chung.2013.MICCAI}. Due to this limitation, there arises a need to develop a multiscale network model that can produce dependable outcomes regardless of the choice of threshold. Recently popular Topological data analysis (TDA) can fill this gap providing a topologically consistent solution across varying thresholds \cite{edelsbrunner.2000, zomorodian.2005, singh.2008, ghrist.2008, carlsson.2008, wasserman.2018}. Within TDA, persistent homology based approaches have become increasingly popular as a tool for analyzing different brain imaging data because they can capture the persistences \cite{oballe2021bayesian, song.2020.ISBI} of different topological features that are robust under different scales \cite{bubenik.2007, ghrist.2008, lee.2012.TMI, chung.2019.NN}. The persistences are usually summarized and expressed using barcodes \cite{song.2020.arXiv}. Although such approaches have been applied to increasingly diverse biomedical problems, they are mostly limited to investigating the static summary of dynamically changing data such as functional magnetic resonance images (fMRI) and electroencephalography \cite{chung.2019.ISBI, bassett.2017, wang.2018.annals}. However, a few recent studies applied TDA to capture dynamic patterns of data including the applications to financial data \cite{gidea.2018} and gene expression data \cite{perea.2015, song.2020.ISBI}.

Motivated by these studies, we develop a novel dynamic-TDA based approach that can statistically discriminate two groups of brain networks. This dynamic-TDA framework builds persistent homology over a multivariate time series \cite{song.2020.ISBI} and constructs a test statistic using the Wasserstein distance between barcodes. The proposed Dynamic-TDA is applied on the resting-state fMRI (rs-fMRI) human and uncovers interesting topological differences in the temporal dynamics between males and females. Unlike many existing approaches limited to  static summaries of dynamically changing images \cite{chung.2019.ISBI, bassett.2017, wang.2018.annals}, the proposed approach helps us investigate the dynamic pattern of the images \cite{gidea.2018, perea.2015, song.2020.ISBI}. Motivated by time-delay embedding (TDE), we propose new dynamic embedding that can extract temporal profile of topological changes in multivariate time series  \cite{von.2010.tde}.

\section{Methods}
The human brain networks can be represented as graphs, where a node represents a parcellated brain region and an edge represents the interaction between two such brain regions \cite{chung.2019.NN}. Then we build a filtration on brain graphs.

\subsection{Graph filtration}
Let $G = (V, \bm{w})$ be a weighted graph, where $V$ is a set of nodes and $\bm{w}$ is a vector of edge weights. A binary graph $G_\epsilon = (V, \bm{w}_\epsilon)$ of $G$ is a graph with binary edge weights $\bm{w}_\epsilon = (w_{\epsilon,i})$:
\begin{align*}
w_{\epsilon, i} = \begin{cases}
1, \text{~~~if~} w_{i}>\epsilon,\\
0,  \text{~~~otherwise.}
\end{cases}
\end{align*}
A graph filtration of $G$ is defined as a collection of nested binary networks:
\begin{align*}
G_{\epsilon_0} \supset G_{\epsilon_1} \supset \cdots \supset G_{\epsilon_k},
\end{align*} where $\epsilon_0 <\epsilon_1<\cdots<\epsilon_k$ are filtration values \cite{2021.vijay,chung.2019.NN, song.2020.arXiv}. Usually edge weights are chosen for filtration values. The graph filtration of a graph with $4$ nodes is illustrated in Figure~\ref{fig:bd_decomposition}.

Unlike the Vietoris-Rips filtration \cite{carlsson2022applying, zhangpersistent}, the filtration does not produce more than 1D homology.   In a binary network $G_\epsilon$, which is a simplicial complex consisting of only nodes and edges, $0$-dimensional ($0$D) holes are \textit{connected components} and $1$-dimensional ($1$D) holes are \textit{cycles}. The number of connected components and the number of independent cycles in a graph are referred to as $0$th Betti number $(\beta_0)$ and $1$st Betti number $(\beta_1)$ respectively \cite{carlsson.2009.bulletin}.

\begin{figure}[!h]
	\centering
	{\includegraphics[width=13cm]{graph_filannoted}}
	\caption{Illustration of graph filtration and birth-death decomposition of a graph with $4$ nodes and $4$ edges. The birth-death decomposition partitions the edge set into disjoint birth and death sets. The numbers are edge weights.}
\label{fig:bd_decomposition}
\end{figure}

In a graph filtration, the sequential removal of edges leads to either birth of a connected component or death of a cycle \cite{chung.2019.NN}. Once a connected component is born, it never dies. Therefore its death values are assumed to be at $\infty$ and ignored. On the other hand, all the cycles are considered to be born at $-\infty$ and their birth values  ignored. Thus we can partition edge weights into the birth set  $B(G)$ and death set $D(G)$. The birth set  $B(G)$ consists of the sequence of increasing birth values
\begin{align*}
B(G): \epsilon_{b_1}<\cdots<\epsilon_{b_{m_0}},
\end{align*} where $m_0 = |V|-1$ \cite{song.2020.arXiv}. $B(G)$ characterizes the $0$D homology of a graph filtration and is  referred as $0$D barcode. The death set  $D(G)$  consists of he sequence of increasing death values
\begin{align*}
D(G): \epsilon_{d_1}<\cdots<\epsilon_{d_{m_1}},
\end{align*}  where $m_1 =  {(|V|-1)\left(|V|-2\right)}/{2}$ \cite{song.2020.arXiv}. $G(G)$ characterizes the $1$D homology of a graph filtration and is  referred as $1$D barcode. Every edge of a graph must be in either $0$D barcode or $1$D barcode and thus have the \textit{birth-death decomposition} of a graph. We illustrate birth-death decomposition of a graph in Figure~\ref{fig:bd_decomposition}. If $0$D and $1$D barcodes are visualized  using the persistent diagrams, they form 1D scatter points vertically (1D barcode) or horizontally (0D barcode). This also diverges from the Vietoris-Rips filtration, which produces 2D scatter points as the 1D persistent digram \cite{ berry.2020,bubenik.2015, chazal.2014,edelsbrunner.2008}.

\subsection{Wasserstein distance}
Barcodes characterize the topology of a network \cite{ghrist.2008}. Therefore, the topological similarity or dissimilarity between two networks can be quantified using the \textit{Wasserstein distance} between barcodes  or persistent diagrams \cite{liang.2018, mi.2020,1974.vallender, song.2021.MICCAI}.

Let $G_1 = \left(V^{\Phi}, \bm{w}^{\Phi}\right)$ and $G_2 = \left(V^{\Psi}, \bm{w}^{\Psi}\right)$ be two given networks with $|V^{\Phi}| = |V^{\Psi}|=p$ nodes. Let $P^\Phi$ and $P^\Psi$ be the corresponding persistent diagrams. The $r$-Wasserstein distance between $P^\Phi$ and $P^\Psi$ is then given by
$$\mathfrak{L}_r(P^\Phi, P^\Psi) = \inf_{\tau: P^\Phi \to P^\Psi} \Big( \sum_{x \in P^\Phi} \| x - \tau(x) \|^r \Big)^{1/r}$$
over every possible permutation or bijection $\tau$ between $P^\Phi$ and $P^\Psi$ \cite{berwald2018computing, panaretos.2019}. When $r \to \infty$, the $\infty$-Wasserstein distance is given by $$\mathfrak{L}_\infty(P^\Phi, P^\Psi) = \inf_{\tau: P^\Phi \to P^\Psi} \max_{x \in P^\Phi} \| x - \tau(x) \|.$$ The $\infty$-Wasserstein distance is also known as \textit{Bottleneck distance} \cite{cohen.2007.stability}. The standard methods for computing the Wasserstein distance is to use the Kuhn-Munkres approach or the Hungarian algorithm that takes $\mathcal{O}(p^3)$ run time \cite{edmonds.1972,munkres1957algorithms, songdechakraiwut2020topological, su2015optimal, su2015shape}. However, for graph filtrations which produces 1D scatter points as barcodes or persistent diagrams, the Wasserstein distances can be computed in  $\mathcal{O}( p \log p)$ run time by simply matching in sorted order in  $\mathcal{O}(p \log p)$ run time.

\begin{theorem}\label{thm1} For 1D scatter points as persistent diagrams $P^\Phi$ and $P^\Psi$, we have
$$\mathfrak{L}_r(P^\Phi, P^\Psi) = \inf_{\tau: P^\Phi \to P^\Psi} \Big( \sum_{x \in P^\Phi} |x - \tau(x) |^r \Big)^{1/r} = \Big( \sum_{x \in P^\Phi} | x - \tau^*(x) |^r \Big)^{1/r},$$ where $\tau^*$ maps the $i$-th smallest value in $P^\Phi$ to the $i$-th smallest value in $P^\Psi$ for all $i$. For $r\to \infty$,
$$\mathfrak{L}_\infty(P^\Phi, P^\Psi) = \inf_{\tau: P^\Phi \to P^\Psi} \max_{x \in P^\Phi} | x - \tau(x)| = \max_{x \in P^\Phi} |x - \tau^*(x)|.$$
\end{theorem}

The first statement in Theorem \ref{thm1} is due to \textit{majorization} \cite{marshall.1979.inequalities}. Vector $\bm{a}$ majorizes another vector $\bm{b}$ if the components of $\bm{a}$ are ``more spread out" than are the components of $\bm{b}$. This is exactly our case if we let $$\bm{a} = \{x-\tau(x): {x} \in P^\Phi\} \text{~~and~~} \bm{b} = \{x-\tau^*(x): {x} \in P^\Phi\}.$$
In such cases, we have
\begin{align*}
\sum_{{x} \in P^\Phi} f(x-\tau(x)) \geq \sum_{{x} \in P^\Phi} f(x-\tau^*(x))
\end{align*} for any convex function $f$.  For $f(x)=|x|^r$, we have Theorem \ref{thm1}. The second statement then follows
\begin{align*}
\mathfrak{L}_\infty(P^\Phi, P^\Psi) =\lim_{r\to\infty} \mathfrak{L}_r(P^\Phi, P^\Psi) =  \lim_{r\to\infty}   \Big( \sum_{x \in P^\Phi} | x - \tau^*(x) |^r \Big)^{1/r} =\max_{x \in P^\Phi} |x - \tau^*(x)|.
\end{align*}

The Wasserstein distance satisfies the following inequality
\begin{align*}
\mathfrak{L}_p(P^\Phi, P^\Psi) \leq \mathfrak{L}_r(P^\Phi, P^\Psi) \leq \mathfrak{L}_1(P^\Phi, P^\Psi)
\end{align*}
for $r<p<\infty$. Since the $1$-Wasserstein distance is the upper bound, it is expected not perform well against other  $r$-Wasserstein distances especially in the presence of outliers. We illustrate this in a simulation study  (Figure \ref{fig:out}). We simulated two edge weight vectors of length $100$ from $\mathcal{N}(0, 0.01^2)$. We introduce outliers into one of the edge weight vectors by replacing entries by values simulated from $\mathcal{N}(1.5, 0.01^2)$. Then computed $1$-,  $2$-, and  $\infty$-Wasserstein distances between them and plotted the resulting distances in Figure~\ref{fig:out}. We observe that the $1$-Wasserstein distance is severely inflated when more outliers are introduced. On the other hand, $2$- and $\infty$-Wasserstein distances are significantly less impacted by the presence of outlying data.

\begin{figure}[!h]
	\centering
	{\includegraphics[width=13cm]{outliers.pdf}}
	\caption{Plot of 1-,  2-, and $\infty$-Wasserstein distances for varying number of outliers in the data.}
\label{fig:out}
\end{figure}

We use Theorem~\ref{thm1} to find the exact form of the {\em Wasserstein distance} between PDs of two networks as follows.
\begin{theorem}\label{theorem:thm2}
 The $r$-Wasserstein distance between the $0$D barcodes for graph filtration is given by
$$\mathfrak{L}^{b}_{r}(P^\Phi, P^\Psi) = \left( \sum_{k=1}^{m_0} \left|\epsilon_{b_k}^{\Phi} - \epsilon_{b_k}^{\Psi}\right|^r\right)^{1/r},$$
where $\epsilon_{b_k}^{\Phi}$ and $\epsilon_{b_k}^{\Psi}$ are the $k$-th smallest birth values of the connected components. The $r$-Wasserstein distance between the $1$D barcodes for graph filtration is given by
$$\mathfrak{L}^{d}_{r}(P^\Phi, P^\Psi)=  \left( \sum_{k=1}^{m_1} \left|\epsilon_{d_k}^{\Phi} - \epsilon_{d_k}^{\Psi}\right|^r\right)^{1/r},$$
where $\epsilon_{d_k}^{\Phi}$ and $\epsilon_{d_k}^{\Psi}$ are the $k$-th smallest death values of the cycles. For $\infty$-Wasserstein distance, the distqances are given by
\begin{align*}
\mathfrak{L}^{b}_{\infty}(P^\Phi, P^\Psi) =\max_{1\leq k \leq m_0}\left|\epsilon_{b_k}^{\Phi} - \epsilon_{b_k}^{\Psi}\right|\text{~~~and~~~}\mathfrak{L}^{d}_{\infty}(P^\Phi, P^\Psi) =\max_{1\leq k \leq m_1}\left|\epsilon_{d_k}^{\Phi} - \epsilon_{d_k}^{\Psi}\right|.
\end{align*}
\end{theorem}

\subsection{Dynamic embedding}

Dynamically changing functional brain imaging data can be viewed as a $d$-variate time series, where $d$ is the number of brain regions we are  investigating. The traditional approach for transforming a univariate time series data into graph data is {\em time delay embedding} (TDE), which is  motivated by Taken's theorem \cite{seversky.2016}. Initially, TDE was used to recover the underlying dynamics of time series \cite{takens.1981}. Recently, it has been used in analyze human speech \cite{brown.2009} and classification of time series data \cite{perea.2015, emrani.2014}. TDE transforms each of the $d$ univariate time series {\em separately} over sliding window of size $M$ to construct $d$ different point clouds in $\mathbb{R}^{M}$ (Figure \ref{fig:tde_schematic}-left). However, it is not immediately clear  how to extend TDE to functional brain imaging data.
TDE also requires an additional step of estimating the dimension $M$ somehow from data \cite{song.2020.ISBI}.
Motivated by TDE, we propose a new {\em dynamic embedding} method, where we converts a $d$-variate time series $\bm{x}_t$ into a time-varying sequence of a point cloud in $\mathbb{R}^{d}$.

\begin{figure}[!t]
	\centering
	{\includegraphics[width=1\linewidth]{dynamic-embdding}}
	\caption{Left: time delay embedding of univariate time series. Sliding window of size 3 is used to embed the time series a static cloud point in $\mathbb{R}^3$.  The proposed dynamic embedding a trivariate time series into a time-varying graph. We used the sliding window of size 5 that produces graphs with $p=5$ nodes in $\mathbb{R}^3$.}
\label{fig:tde_schematic}
\end{figure}

Let $\bm{x}_t=(x_{1t}, \ldots, x_{dt})^\top$ be a $d$-variate time series. For each time point $t$, $\bm{x}_t$ will be projected as a scatter point in $\mathbb{R}^d$ by taking $x_{it}$ as the $i$-th coordinate of the scatter point. Time series $\bm{x}_t, \ldots, \bm{x}_{t+p-1}$ in sliding window $[t, t+p-1]$ will be embedded into $p$ scatter points in $\mathbb{R}^d$ \cite{song.2020.ISBI}. The Euclidean distance between scatter points will be used as edge weights in building time varying graph $G_t = (V_t, \bm{w}_t)$, where vertex set $V_t = \{ \bm{x}_t, \ldots, \bm{x}_{t+p-1} \}$
and the edge weights $ \bm{w}_t$ are Euclidean distances between the corresponding scatter points. Figure \ref{fig:tde_schematic}-right illustrates an example with a trivariate time series that is turned into a time-varying graph in $\mathbb{R}^3$. Subsequently, we build graph filtrations to obtain the dynamically changing $0$D and $1$D barcodes (Figure \ref{fig:simu_beta_bd}):
$$B(G_t) = \{{\epsilon}_{b_{1t}}, \ldots, {\epsilon}_{b_{m_0t}}\} \text{~~and~~} D(G_t) = \{{\epsilon}_{d_{1t}}, \ldots, {\epsilon}_{d_{m_1t}}\}.$$
While TDE implicitly fuses the dynamic aspect of the time series into a cloud point, one distinct advantage of the proposed Dynamic-TDA method is that the dynamic topological aspect of the time series can be easily extracted.

Consider  two $d$-variate time series
$\bm{x}_t=(x_{1t}, \ldots, x_{dt})^\top$ and $\bm{y}_t=(y_{1t}, \ldots, y_{dt})^\top$ \cite{kakizawa1998discrimination, shumway2003time, fryzlewicz2009consistent, bohm2010classification}. We are interested in determining the topological difference between the time series. Following the embedding process, we build time varying graphs
$X_t$ and $Y_t$ with $p$ nodes corresponding to time series $\bm{x}_t$ and $\bm{y}_t$. The total number of finite birth values of connected components and the total number of death values of cycles are \cite{songdechakraiwut2020topological}
$$m_0 = p-1 \text{~~and~~} m_1 = 1+ \frac{p\left(p-3\right)}{2}.$$

We can compute the  $0$D and $1$D barcodes at time point $t$ corresponding to time series $\bm{x}_t$ and $\bm{y}_t$.
The $r$-Wasserstein distance between the $0$D barcodes at time point $t$ is the given by
\begin{align*}
\mathfrak{L}_{r,t}^b(\bm{x}, \bm{y}) = \left(\sum \limits_{k=1}^{m_0} \left|{\epsilon}_{b_{kt}}^{\bm{x}} - {\epsilon}_{b_{kt}}^{\bm{y}}\right|^r\right)^{1/r},
\end{align*}
where ${\epsilon}_{b_{kt}}^{\bm{x}}$ and ${\epsilon}_{b_{kt}}^{\bm{y}}$ are the $k$-th smallest birth values associated with graphs $X_t$ and $Y_t$ respectively. Similarly, the $r$-Wasserstein distance between the $1$D barcodes at time point $t$ is
\begin{align*}
\mathfrak{L}_{r,t}^d(\bm{x}, \bm{y}) = \left(\sum \limits_{k=1}^{m_1} \left|{\epsilon}_{d_{kt}}^{\bm{x}} - {\epsilon}_{d_{kt}}^{\bm{y}}\right|^r\right)^{1/r},
\end{align*}
where ${\epsilon}_{d_{kt}}^{\bm{x}}$ and ${\epsilon}_{d_{kt}}^{\bm{y}}$ are the $k$-th smallest death values associated with graphs $X_t$ and $Y_t$ respectively. Similarly, $\infty$-Wasserstein distance is given as
$$\mathfrak{L}_{\infty,t}^b(\bm{x}, \bm{y})  =  \max_{1\leq k \leq m_0} \left|{\epsilon}_{b_{kt}}^{\bm{x}} - {\epsilon}_{b_{kt}}^{\bm{y}}\right|,$$
$$\mathfrak{L}_{\infty,t}^d(\bm{x}, \bm{y})  =  \max_{1\leq k \leq m_0} \left|{\epsilon}_{d_{kt}}^{\bm{x}} - {\epsilon}_{d_{kt}}^{\bm{y}}\right|.$$

\begin{figure}[t]
	\centering
	\includegraphics[width=0.6\linewidth]{simu_beta}
	\caption{Plot of time series of $0$D and $1$D barcodes obtained using the dynamic embedding method on a simulated 4-variate time series with Beta$(5,2)$ noises in the simulation study 2. This produces time varying 4-nodes graphs. We observe periodic behaviors in both the $0$D and $1$D barcodes, which are consequences of the sinusoidal variation in the simulation model.}
	\label{fig:simu_beta_bd}
\end{figure}

Subsequently, we propose a statistic combining the $r$-Wasserstein distances over all time points as
\begin{align*}
\mathcal{L}_r(\bm{x}, \bm{y})
= \sum_t \mathfrak{L}_{r,t}^b(\bm{x}, \bm{y})  + \sum_t \mathfrak{L}_{r,t}^d(\bm{x}, \bm{y}).
\end{align*}
For $r\to \infty$, we have
\begin{align*}
\mathcal{L}_\infty(\bm{x}, \bm{y})
=\sum \limits_{t} \mathfrak{L}_{\infty,t}^b(\bm{x}, \bm{y})    + \sum \limits_{t} \mathfrak{L}_{\infty,t}^d(\bm{x}, \bm{y})  .
\end{align*}

A large value of $\mathcal{L}_r(\bm{x}, \bm{y})$ or  $\mathcal{L}_\infty(\bm{x}, \bm{y})$ implies that the dynamic barcodes produced from  the two time series $\{\bm{x}_t\}$ and $\{\bm{y}_t\}$ are  {\em topologically distant}  and there exists a large difference in the $0$D and $1$D topological features between the two time series $\bm{x}_t$ and $\bm{y}_t$. Since the distribution of the proposed test statistic $\mathcal{L}_r(\bm{x}, \bm{y})$ or $\mathcal{L}_\infty(\bm{x}, \bm{y})$ is not known, we apply the permutation test, which do not assume any probability distribution for inference \cite{thompson.2001, zalesky.2010, nichols.2002, winkler.2016}.

\section{Simulation study}
Since there is no ground truth in real data, we performed the simulation studies where the ground truths are specified.
This controlled experiments allows us to determine the performance of the proposed method.

\subsection{Study 1} Consider a trivariate time series from a vector autoregressive model of order one \cite{zivot2006vector, lutkepohl2013vector, haslbeck2021tutorial}:
\begin{align*}
\bm{x}_t = \bm{a} + \bm{b}t + \bm{\Theta}\bm{x}_{t-1}  + \bm{\xi}_t,
\end{align*}
where noise $\bm{\xi}_t \thicksim \mathcal{N}_3 (\bm{0}, \bm{\Sigma})$ at each time $t$.
Also consider the other trivariate time series that is simply Gaussian white noise at each time point:
\begin{align*}
\bm{y}_t  \thicksim \mathcal{N}_3 (\bm{0}, \bm{\Sigma}).
\end{align*}
We chose the trend parameters $\bm{a} = (1, 1, 0)^\top$ and $\bm{b} = (1.5, 2, 0)^\top$ for $1 \leq t \leq 30$.
The coefficient and covariance matrices were chosen as
\begin{align*}
\bm{\Theta} =  \begin{bmatrix}
0.2 &-0.1 &0.5\\
-0.4 &0.2 &0 \\
-0.1 &0.2 &0.3
\end{bmatrix}, \quad  \bm{\Sigma} = \begin{bmatrix}
0.1 & 0.01 & 0.3\\
0.01 & 0.5 & 0 \\
0.3 & 0 & 1
\end{bmatrix}.
\end{align*}

We used $\mathcal{L}_2({\bf x}, {\bf y})$ as the test statistic and applied the permutation test with 100,000 permutations. The simulation is independently performed 10 times and we averaged the results.
The obtained $p$-values are $10^{-6} \pm 3.16 \times 10^{-6}$ indicating strong topological difference between two trivariate time series.

We also checked if the method is not incorrectly producing signals when there is expected to be no topological differences. We generated another time series
\begin{align*}
\bm{z}_t = \bm{a} + \bm{b}t + \bm{\Theta}\bm{x}_{t-1}  + \bm{\zeta}_t,
\end{align*}
where noise $\bm{\zeta}_t \thicksim \mathcal{N}_3 (\bm{0}, \bm{\Sigma})$ at each time $t$. Since $\bm{x}_t$ and $\bm{z}_t$ are identical time series model with independently generated noises, we should not expect the method to detect any topological differences. The simulation is independently performed 10 times and we averaged the results.
The obtained $p$-values are $0.42 \pm 0.30$ indicating the method is not able to detect any topological differences.

\begin{figure}[!h]
	\centering
	{\includegraphics[width=8cm]{simulation1.pdf}}
	\caption{The plots of three independently simulated trivariate time series study 1. Each trivariate time series is represented as a random walk in $\mathbb{R}^3$. They are generated from the same vector autoregressive model and should not show topological differences. The proposed Wasserstein distance  is not detecting any difference.}
\label{fig:simulation1}
\end{figure}

\subsection{Study 2} We conducted another simulation study to determine the performance of the proposed method in
 discriminating two $d$-variate time series with Beta noises. We considered a time series model with a trend, periodicity, and perturbation
\begin{align*}
\bm{x}_t = (a+bt)\bm{1} + c\sin \left(\frac{2\pi t}{P}\right)\bm{1} + \bm{\xi}_t, \end{align*}
where noise follows the Beta distribution with shape parameter parameters $\alpha$ and $\beta$, i.e.,
$\bm{\xi}_t \thicksim \text{Beta}(\alpha, \beta)$. The shape parameters $\alpha$ and $\beta$ allow to take a wide variety of distributional patterns. We simulated two $d$-variate time series using the same model with different values of $\alpha$ and $\beta$.  We used trend parameters $a = c =0.1$ and $b=0.001$, period of the sinusoidal component $P=8$, length of the time series $T=50$, dimension $d=10$. Three different distributions Beta$(1, 1)$, Beta$(5, 2)$ and Beta$(1, 5)$ are used. The above model is used to simulate time series at  $p=5$ nodes. The proposed method is used to compute the Wasserstein distances between them.  The permutation test with 100,000 permutations was used. The simulation is indecently performed 10 times and average $p$-values are tabulated in Table~\ref{Table:simu-dtda}.
\begin{table}
\centering
		\caption{The average $p$-values obtained using the Dynamic-TDA method with 1-, 2-, and $\infty$-Wasserstein distances for various pairs of noise distributions. The $p$-values smaller than 0.05 indicate that our method can identify topological differences in time series at a $95\%$ confidence level.\vspace{0.2cm}}
		\label{Table:simu-dtda}
		\renewcommand{\arraystretch}{1.2}
		\setlength{\tabcolsep}{18pt}
		\resizebox{\columnwidth}{!}{%
		\begin{tabular}{c|ccc}
			\centering
			Distribution & $\mathcal{L}_1(\bm{x}, \bm{y})$  & $\mathcal{L}_2(\bm{x}, \bm{y})$ & $\mathcal{L}_\infty(\bm{x}, \bm{y})$\\
			\hline
			Beta$(1,1)$ vs. Beta$(5,2)$ & $0.0001$ & $0.0000$ & $0.0000$  \\
			Beta$(1,1)$ vs. Beta$(1,5)$ & $0.0000$  & $0.0000$ & $0.0000$ \\
			Beta$(5,2)$ vs. Beta$(1,5)$ &  $0.0230$  & $0.0182$ & $0.0260$ \\
			\hline
			Beta$(1,1)$ vs. Beta$(1,1)$ & $0.6795$ & $0.9563$ & $0.9114$ \\
			Beta$(5,2)$ vs. Beta$(5,2)$ & $0.1726$ & $0.5581$ & $0.9410$ \\
			Beta$(1,5)$ vs. Beta$(1,5)$ & $0.3285$ & $0.5186$ & $0.5105$ \\
			\hline
		\end{tabular}}
\end{table}

\begin{figure}[!h]
	\centering
	{\includegraphics[width=6.2cm]{dtda1}}
	{\includegraphics[width=6.2cm]{dtda2}}
	\caption{Histogram of $2$-Wasserstein distance and the corresponding observed test statistics (in dotted red lines) for the scenarios: Beta$(1,1)$ vs. Beta$(5,2)$ (left) and Beta$(1,1)$ vs. Beta$(1,1)$ (right).}
\label{fig:hist}
\end{figure}

Time series drawn from the same Beta distribution produced large $p$-values and those drawn from different Beta distributions had $p$-values smaller than $0.05$. We conclude that the Dynamic-TDA method can discriminate two multivariate time series whenever there is a difference at a $95\%$ confidence level. Figure~\ref{fig:hist} plots the histograms of $\mathcal{L}_2(\bm{x}, \bm{y})$ for two specific examples Beta$(1,1)$ vs. Beta$(5,2)$ (top left) and Beta$(1,1)$ vs. Beta$(1,1)$ (top right).

\section{Application to functional brain images} \label{real}
\subsection{Data description} We used  a resting-state fMRI dataset collected as part of the Human Connectome Project (HCP) \cite{vanessen.2012, glasser2013minimal}. The dataset comprises fMRI scans of $412$ subjects ($172$ males and $240$ females) measured over $1200$ time points using a gradient-echoplanar imaging sequence \cite{song.2020.arXiv, 2021.vijay}. Since there are more than $300000$ voxels in fMRI, we parcellated the brain into 116 non-overlapping anatomical regions using the Automated Anatomical Labeling (AAL) template and averaged rsfMRI over each parcellation  (Figure \ref{fig:AAL}) \cite{tzourio.2002,desikan.2006, hagmann.2007, arslan.2018}.  The head movements cause serious spatial artifacts in functional connectivity \cite{power.2012, van.2012, satterthwaite.2012, caballero.2017}. We calculated the framewise displacements from the three translational displacements and three rotational displacements at each time point to measure the head movement from one volume to the next. The volumes with framewise displacement larger than $0.5$ mm and their neighbors (one back and two forward time points) were scrubbed \cite{van.2012, power.2012, huang.2020.NM}. More than one third of $1200$ volumes were scrubbed in $12$ subjects. We  removed the 12 subjects ($4$ males and $8$ females).

\begin{figure}[t]
	\centering
	{\includegraphics[width=1\linewidth]{rsfMRI4regions}}
	\caption{Time series of functional brain connectivity between left and right hemispheres in 4 regions. There are total 6670 time series of brain connectivity between 116 region we analyzed in the study.}
\label{fig:AAL}
\end{figure}

\begin{figure}[!h]
	\centering
	{\includegraphics[width=13cm]{female_male}}
	\caption{Static female (left) and male (right) brain networks  obtained by computing the Pearson correlation of rsfMRI in 116 brain regions (black nodes) over whole time point. Correlation networks are too dense to visualize so we only displayed connections above 0.5.}
\label{fig:fm}
\end{figure}

\begin{figure}[!h]
	\centering
	{\includegraphics[width=13.5cm]{avg}}
	{\includegraphics[width=13.5cm]{avg_polish}}
	\caption{Top: The average time series of functional brain connectivity for female (solid black) and male (dotted black) averaged over edges. Bottom: Corresponding detrended female (solid black) and male (dotted black) time series. We removed the first 60 time points in the analysis.}
\label{fig:fm_time}
\end{figure}

Subsequently we computed the Pearson correlation between brain regions over sliding window of size 30 in measuring the strength of time varying brain connectivity.   There are total $(116 \times 115)/2 = 6670$ time series of brain connectivity between 116 regions (Figure \ref{fig:AAL}).  Often we represent the {\em static} summary of brain networks by computing Pearson correlation over the whole time points (Figure \ref{fig:fm}). These time series are feed into the proposed pipeline. The time series are averaged over all connections to determine if there is an overall trend (Figure \ref{fig:fm_time}-top).
We observe a sharp decline in the initial measurements of fMRI, which is usually considered as artifacts included by the large variance at the beginning of each MRI scan \cite{diedrichsen.2005}. Thus, we removed the first $60$ time points and only used the remaining $1140$ time points in the subsequent analysis. There exists a slow increasing trend. We fit a smoothing public spline  $f(t)$ to to time varying functional connectivity $x_t$ at each edge and determined the trend  by minimizing \cite{green1993nonparametric, rj1990generalized}:
$$\sum_{t=1}^T \big| x_t -  f(t) \big|^2 + \lambda \int  f''(t)^2 \, \text{d}t.$$
We then removed the fitted trend and obtained the detrended time series $x_t- f(t)$ (Figure~\ref{fig:fm_time}-bottom).

\subsection{Global topological inference} We determined if there is overall topological difference between resting-state functional brain connectivity between males and females. We constructed time-varying graphs over with $p=20$ nodes, computed dynamic barcodes. We performed the permutation test with half million random permutations using the Wasserstein distance  $\mathcal{L}_2(\bm{x}, \bm{y})$. We obtain the observed test statistic as $3257.2$ and the corresponding p-value $< 0.001$. Thus, we conclude that the proposed method distinctly discriminates the functional connectivity between females and males. Figure~\ref{fig:real_hist} plots the histogram of the test statistic and the corresponding observed test statistic (in dashed red).

\begin{figure}[!t]
	\centering
	{\includegraphics[width=12cm]{real_dtda_1}}
	\caption{The histogram of the $2$-Wasserstein distance obtained using the resting-state fMRI of human. The dashed red line represents the observed value of the distance.}
\label{fig:real_hist}
\end{figure}

\subsection{Local topological inference} It is of interest to quantify the differences between the females and males in a human brain network at the connection level. We used TDE, which can perform localized topological inference at the edge level. We limited the analysis to  connections that give
Pearson correlation above 0.7 over the whole time point using the univariate time series of connectivity. We applied TDE and compute the corresponding $2$-Wasserstein based statistic between females and males at each edge separately. A higher value of the statistic means the female and male are different while a lower value implies no significant difference at the corresponding edge. Any edge that shows statistic value above 4.5783 is significant at a $95\%$ confidence level. We plotted the $2$-Wasserstein based TDE statistic  map in Figure~\ref{fig:whole}.

\begin{figure}[t]
	\centering
	{\includegraphics[width=1\linewidth]{central}}
	\caption{The map of the $2$-Wasserstein based TDE statistic of the brain. A higher value of distance indicates larger topological distance of underlying time series of connectivity at the edge level. We observe the most connected regions are within the occipital lobe, specifically the calcarine fissure and surrounding cortex, cuneus and lingual gyrus. Also, the precentral gyrus shows the difference. Any connection above disatnce 4.5783 is statistically difference at a $95\%$ confidence level.}
\label{fig:whole}
\end{figure}

\subsection{Modeling dynamic barcodes} The proposed dynamic embedding produces time series of $0$D and $1$D barcodes, which can be further modeled to interpolate topology at another unknown time point. Below we illustrate modeling such dynamic aspects of the topological features of the brain data. As an  illustration, we considered $3$ edges in the frontal brain lobes of the both female and male brain networks. We applied our method  to the trivariate time series to construct time-varying graphs with $p=40$ nodes and computed the barcodes. We plotted the time series of one of the $0$D barcodes (in solid red) of the females in Figure~\ref{fig:spline}. The high frequency of the points in the plot readily suggests applying a high degree polynomial interpolation to the data. However, such interpolation often causes Runge's phenomenon \cite{fornberg2007runge, epperson1987runge, boyd1992defeating}. To avoid this, we used the cubic spline to fit the series. We choose the smoothing parameter $\lambda=0.001$. The parameter is chosen empirically so that it captures the overall patterns of the barcodes. The fitted spline (in dashed blue) is superimposed on the original barcode. To determine the goodness of fit of the model, we used the mean squared error (MSE):
\begin{align*}
\text{MSE} = \sum (\text{observed} - \text{fitted})^2.
\end{align*} The MSE of the fitted spline to the considered time series of barcode is $0.0069$ suggesting a reasonable fit to the original data.

\begin{figure}[t]
	\centering
	{\includegraphics[width=1\linewidth]{represent}}
	\caption{Plot of a time series of a representative $0$D barcode (solid red) and the corresponding fitted spline (dashed blue).}
\label{fig:spline}
\end{figure}

\section{Conclusion} TDA based approaches have been previously applied to static brain imaging studies \cite{chung.2019.ISBI, bassett.2017, wang.2018.annals}, where topological features that characterize the brain signals are identified. However, they are limited to investigating static summaries of brain networks. Studies limited to static summaries are bound to give incomplete results. In this paper, we propose new dynamic-TDA in statistically discriminating between two groups of multivariate time series. The proposed framework builds persistent homology over a multivariate time series and constructs a statistic based on $r$-Wasserstein distance between persistences. We applied this method to resting state fMRI data of female and male human subjects and discover a significant topological difference in terms of combined $0$D and $1$D topological distance between them. Considering the success of the proposed approach in discriminating multivariate time series, we believe that the proposed method can also be applied to other type of time series data.

\section*{Acknowledgments}
We would like to thank D Vijay Anand and Sixtus Darkurah for discussion on TDA. This study is funded by NIHR01 EB02875 and NSF MDS-2010778.

\bibliographystyle{aims}