\section{Introduction}
A strong classical electromagnetic background can lead to vacuum instability and produce particle-antiparticle pairs. This process of particle creation from the quantum vacuum was first studied in 1951 by Schwinger under a constant electric field, and this phenomenon is known as the Schwinger effect \cite{PhysRev.82.664}.This particle creation paradigm has crucial importance for nonequilibrium processes in heavy-ion collisions\cite{24,25,26} as well as astrophysical phenomena\cite{27}and the search for nonlinear and nonperturbative effects in ultraintense laser systems\cite{12,13,14}.
\newline
Particle production is the process of evolving a quantum system from an initial equilibrium configuration to a new final equilibrium configuration via an intermediate  non-equilibrium evolution caused by a strong field background.
Such strong field QED pair production process has been investigated using scattering calculations \cite{Harvey:2009ry}, exact solutions \cite{Dunne:1998ni}, semi-classical techniques \cite{Kim:2000un,DiPiazza:2004lsj}, Monte Carlo simulations \cite{Gies:2005bz} and quantum kinetic equations \cite{Alkofer:2001ik,Blaschke:2005hs}
Quantitative description of particle production at all times in time-dependent electromagnetic field is not possible due to  the absence of unique separation into positive and negative energy states at intermediate times and these positive and negative states   are well-defined  only at asymptotically early and late times where  the field vanishes. A common approach is to define particle numbers in terms of an adiabatic basis using Bogoliubov transformation \cite{2,3,4,5}.
In the adiabatic basis, we examine the problem of pair production in a time-varying spatially uniform Sauter electric field  which has been studied by various authors \cite{6,44}. who calculated  the  number of particles created  at  the asymptotic time but the problem of particle  production at the finite time is not studied for Sauter-Pulse electric field.
We looking for the evolution of the quantum system at some initial time $t_0$ in the vacuum state but now what will be the properties of the quantum system at finite time $t$ ?  What happens to the system properties at that non-asymptotic time?
By finding the exact analytic solution of the mode function, which allows us to study the behaviour of the quantum vacuum during pair production at any instant of time for time-dependent Sauter pulsed electric fields in spinor quantum electrodynamics (QED).
We use  natural units and set $ \hslash = c = m = 1 $ and express all variables in terms of the electron mass unit.

\section{Theory}
Consider the creation of the electron-positron pair from vacuum  by a linearly polarized  time-dependent spatially uniform  electric field along the $z$-axis which is characterized by the four-vector potential $A_\mu(t) = (0, 0, 0, A(t))$ with $E(t) = - \dot{A} (t).$
In our calculations,  we choose the Sauter field in the direction of $z-$axis, which is given by
\begin{equation}
    E(t) = - \dot{A} (t)= E_0 sech^2(t/\tau)
\end{equation}
and corresponding vector-potential
\begin{equation}
    A(t) = -E_0  \tau tanh(t/\tau).
\end{equation}
where $E_0$ amplitude of the electric field and $\tau$ pulse duration of the electric field.
The one-particle momentum distribution function $f(\textbf{p},t)$ is an important quantity in the description of the particle production process in the time-dependent electric field.
As it was demonstrated in \cite{41,42}, f(\textbf{p},t) can be written in terms of the mode solution as
\begin{equation}
 f(\textbf{p},t) = \frac{{ (\omega + \Pi) |\dot\Psi^{(\pm)} \pm i \omega \Psi^{(\pm)}|^2}}{2 \omega}
 \label{f}
\end{equation}
where the mode equation \cite{40,43}
\begin{equation}\label{dir_eq1}
\ddot\Psi^{(\pm)}+[\omega^2(t)\pm ie\dot A(t)]\Psi^{(\pm)}=0
\end{equation}
with
\begin{equation}\label{omega}
\omega(t)^2=\Pi^2+\epsilon_\perp ^2,\quad \Pi=p_z+eA(t),\quad \epsilon_\perp^2=\textbf{p}_\perp^2+m^2.
\end{equation}
For Sauter-Field, we are able to exactly solve the mode equation \ref{dir_eq1}, which becomes the hyper-geometric differential equation and its solution is given in terms of hypergeometric functions. Using this solution, we calculate the one-particle distribution function in terms of hypergeometric functions. The time evolution of the particle distribution function $f(\textbf{p},t)$ in momentum space is studied for $E_0 = 0.2  $ and $\tau = 10.$
\begin{figure}[t]
\begin{center}
{
\includegraphics[width = 2.5in]{t_neg7_p3.eps}
\includegraphics[width = 2.5in]{t_0_p3.eps}
\includegraphics[width = 2.5in]{t_11_p3.eps}
\includegraphics[width = 2.5in]{t_24.5_p3.eps}
\includegraphics[width = 2.5in]{t_30_p3.eps}
\includegraphics[width = 2.5in]{t_36_p3.eps}
\includegraphics[width = 2.5in]{t_45_p3.eps}
\includegraphics[width = 2.5in]{t_100_p3.eps}
}
\caption{Time evolution of quasi-particle distribution function in   the longitudinal momentum space $f(p_\parallel, t)$  at different times. The transverse momentum is considered to zero, and all the units are taken in the electron mass units.The field parameters are  $E_0=0.2 E_c$ and $ \tau =10.$}
   	\label{fig:1}
\end{center}
\end{figure}
\section{Result}\label{Result}
In this section, we discuss about both the longitudinal and transverse  momentum spectrum of the created particle using the solution of the mode equation in \ref{f} for $E_0 = 0.2$ and $ \tau = 10 $.

\subsection{Longitudinal momentum spectrum}
In Figure \ref{fig:1}, we show how the particle-creation proceeds. The produced electron is accelerated to the negative z-direction by the electric field. At early times (QEPP region \cite{46}) when electric field$E(t)$ is increasing, momentum distribution shows the smooth Gaussian-like structure and after reaching  a maximum value at time $t=0$ electric field  starts decaying,  and we see the peak of the momentum distribution function  is shifted as shown in figure \ref{fig:1}(a),(b),(c).
As time progresses, the distribution function tends to be peaked around the zero value of the longitudinal momentum as seen in figure \ref{fig:1}(d). In this process, it shows a complex behavior of splitting and manifests oscillation arising at $t = 3 \tau $ (near REPP region ) where the electric field nearly vanished see figure \ref{fig:1}(e),(f). These complex behavior disappear at the asymptotic time ( $t= 10 \tau$), and  we see a smooth Gaussian-like structure \ref{fig:1}(h).
\begin{figure}[t]
\begin{center}
{
\includegraphics[width = 2.4755in]{t_p1_0.eps}
\includegraphics[width = 2.55in]{t_p1_7.eps}
\includegraphics[width = 2.55in]{t_p1_14.eps}
\includegraphics[width = 2.55in]{t_p1_22.eps}
\includegraphics[width = 2.55in]{t_p1_28.eps}
\includegraphics[width = 2.55in]{t_p1_35.eps}
\includegraphics[width = 2.55in]{t_p1_50.eps}
\includegraphics[width = 2.55in]{t_p1_100.eps}
}
\caption{Time evolution of quasi-particle distribution function in   the transverse  momentum space $f(p_\perp, t)$  at different times. The longitudinal momentum is considered to be zero, and all the units are taken in electron mass units.The field parameters are  $E_0=0.2$ and $ \tau =10.$}
   	\label{FG:2}
\end{center}
\end{figure}
\subsection{Transverse momentum spectrum }
Time evolution of the transverse momentum spectrum of created quasi-particle pairs are depicted in the figure \ref{FG:2}.The shape of the transverse
momentum distribution is changed when the electric field is switched on. As we see in the figure \ref{FG:2}(c),(d) at the transient stage the peak at $p_\perp = 0$ splits into weakly
pronounced peaks at $p_\perp \approx  \pm 0.67m$ and $p_\perp \approx  \pm 0.84m.$
At the asymptotic time (REPP stage \cite{46}), that splitting disappears and distribution function $f(p_\perp)$ shows the smooth structure with a peak at zero transverse momentum as shown in figure \ref{FG:2}(h).
\section{Conclusion}
We have analyzed the electron-positron pair creation from  the vacuum by a time-like Sauter pulse electric field. We have seen that the particle Longitudinal momentum distributions exhibit the oscillating structure at  a finite time where the electric field is nearly zero and this oscillating structure can be understood in the Dynamical Tunneling picture\cite{45}.
Dynamical tunneling can be understood as the process that develops in time inter-band dynamics in momentum representation for particles. There are different possibilities of particle tunneling through different channels in momentum space( or momentum states).
As time progresses, different possibilities can occur: (a) particle can tunnel at time $t$ with  $p_\parallel'-$value of momentum.
(b)At early stage,  particle with a lower momentum value ($p_\parallel < p_\parallel'$) can tunnel and then be accelerated to momentum $p_\parallel'$,
(c)Particle with a higher momentum value ($p_\parallel > p_\parallel'$) can tunnel at  time, $t_1$, and is then decelerated to get that momentum $p_\parallel'$ value.
Finally, when  the individual probability amplitudes of this process are added together, they result in quantum interference at time  $t$.
At the asymptotic time, those processes do not share the same phase information and become random due to averaging over particle paths that do not show quantum interference effects, which leads to a smooth distribution function.
However, multiphoton absorption is not seen for the value of the Keldysh parameter $\gamma < 1,$ but it is believed to give a smaller contribution in comparison to what we get from tunneling. Since, in our calculation, it is not much less than one ($\gamma = 0.5$). As a result, another process of particle creation with momentum $p_\parallel$, namely the multiphoton process, cannot be ruled out. These are possibilities, and extensive quantitative  analyses of   these processes  are in progress\cite{DEEPAK:2023}.
The transverse momentum distribution shows only the splitting of smooth structure and the absence of quantum interference, which is an obvious interference effect that occurs only in the direction of the electric field.
\newline
A more detailed investigation of the  momentum distribution of created particles at the finite time is currently underway \cite{DEEPAK:2023}.

\section{Acknowledgments}
DEEPAK gratefully acknowledge the financial support from Homi Bhabha National Institute (HBNI) for carrying out this research work. Deepak also  thanks the organizers of the QED Laser Plasmas International Workshop 2022, Dresden for the possibility to present our work.


\section{Introduction}
Surface water maps provide important knowledge to help understand the impacts and relationships between the environment and human activities. Knowing where water is and how it changes and relates to the surrounding environment enables successful responses that are closely related to water resources \citep{poff2016sustainable,viviroli2011climate,florke2018water}, biodiversity \citep{cawse2021nasa}, and human well-being \citep{sanders2022large}.

Optical remote sensing technology has become the most widely adopted and powerful tool for surface water mapping \citep{huang2018detecting}. Standard practice for extracting water using optical imaging is to exploit the spectral characteristics of water that the reflectance of water is significantly lower in the infrared channel than other channels compared to other objects \citep{xu2006modification, ji2009analysis}. Large bodies of water, such as oceans, rivers, and lakes, have been the focus of water mappings since early Earth observation satellite missions \citep{pekel2016high,khandelwal2017approach}, and reflectance-based water mapping has shown successful performance on medium- and low-resolution ($>$100$m$) satellite imagery \citep{feyisa2014automated,yamazaki2015development}.

With the development of high-resolution imaging technology, mapping of very small water bodies ($<$1$ha$ or 10,000$m^2$) also become possible. Particularly, as small water bodies, such as small lakes, ponds, low-level streams, and ditches, are closely related to people's daily lives and are highly affected by anthropogenic footprints, the demand for a very high-resolution map of water bodies has been rapidly increasing these days \citep{kelly2022opportunities,xu2020towards,hoekstra2018urban,biggs2017importance,riley2018small}. However, in high-resolution images, conventional reflectance-based water mapping does not show satisfactory performance \citep{ogilvie2018surface} and often underestimates the extent of small water bodies \citep{pickens2020mapping,mao2022high}. This is mainly because the reflectance of small water bodies is highly varying in high-resolution images while high-resolution images require multiple image acquisitions to cover large-area which will again significantly increase the spectral variability.

To address this problem, efforts have been made to develop more sophisticated methods. Integrating multi-temporal \citep{pickens2020mapping} or different types of sensors (e.g., radar) \citep{ahmad2019fusion} or developing data-driven methods \citep{ko2015classification,isikdogan2017surface, chen2020novel, wang2020urban} are examples.
However, fundamentally, algorithms based on optical images cannot be free from spectral characteristics of water whose values vary depending on image acquisition conditions e.g., atmospheric conditions and sensor specifications) \citep{martins2017assessment} or surface water conditions (e.g., turbidity and surface roughness) \citep{odermatt2012review,kim2016evaluation}.

In addition, current optical image-based water mapping algorithms are limited to providing only 2D information. Meanwhile, even simply for identifying the direction of water flow, 3D topographic information is needed. Obviously, 3D water information can provide better information for understanding the dynamics of water, and it enables in-depth analysis of how water relates to its surrounding environments and improves the management of water resources and water-related disasters \citep{musa2015review,jarihani2015satellite,arrighi2019effects}.
Thus, the authors note that airborne laser scanning will eventually be required for surface water maps to be truly useful as the LiDAR data is needed to obtain high-resolution 3D topographic information. Nevertheless, as airborne laser scanning has not been effectively used for surface water mapping, separate resources and efforts have been required for completing a 3D hydrography.

This paper presents a scalable water mapping method using topographic airborne LiDAR data. The scalable method refers to the method that yields consistently satisfactory results across wide and varied landscapes. Our method operates in a fully automated manner, only requires topographic airborne LiDAR data, and extracts accurate and high-resolution water bodies including the elevation of water bodies. Extensive experiments with large ($\approx$ 2,500$km^{2}$) and diverse landscapes (urban, coastal, and mountainous areas) confirmed that our algorithm can generate more accurate results than NDWI-based methods even without any additional parameter tunings. Moreover, the proposed method enables high-resolution mapping of a full 3D topographic map (water and terrain) with only airborne LiDAR data.

The remainder of this paper is organized as follows. Section 2 illustrates related works. Section 3 describes the proposed surface water mapping methods. Section 4 illustrates datasets and experimental design. Section 5 provides experimental results. Section 6 concludes the paper.

\section{Related works}
\subsection{Optical imagery-based water mapping}

Optical imagery has been dominantly used for surface water mapping. With multiple satellite missions at different spatial resolutions ranging from 1 km to 1-m resolution, water bodies of the Earth have been mapped for more than several decades \citep{pekel2016high,khandelwal2017approach}. In particular, with recent advancements in satellite, airborne, and UAV imaging technology, a large-area mapping of small water bodies has attracted more attention \citep{becker2019unmanned, straffelini2021mapping,qayyum2020glacial, sogno2022remote}.

Although the higher spatial resolution imagery enables observing small water bodies, large-area mapping of small water bodies is more complicated than that with the lower resolution imageries for large water bodies. This is because, first, small water bodies have more wide variety of spectral reflectance \citep{ogilvie2018surface}. Second, shadow, occlusion, and relief displacement of urban structures make water mapping more challenging \citep{yang2018urban}
. Third, large-area mapping with high-resolution imagery necessitates multiple different image acquisitions, while different acquisitions may require different threshold values for water index-based mapping \citep{feyisa2014automated}.

To improve the accuracy, some recent studies have developed more sophisticated decision rules.
\citep{yang2018urban} alleviated errors from built-up areas by utilizing multiple water and shadow indices.
\citep{dong2022mapping} integrated the spatial information into the classification rule by introducing the roughness of the water index.
In addition, there have been studies that used data-driven approaches such as support vector machine \citep{sun2015soft}, random forest \citep{ko2015classification}, and deep learning methods \citep{chen2020novel, isikdogan2017surface, wang2020urban}.

More often, water has been mapped with other multiple classes including ice and cloud \citep{wu2018scribble} or more comprehensive land covers \citep{song2019patch, robinson2019large,zhang2022urbanwatch}. Most methods obtained satisfactory results under the condition that sufficient high-quality training samples are available.
However, acquiring high-quality training samples for large areas is a challenging task, particularly in high-resolution image, for the same reason that surface water mapping using spectral characteristics is difficult \cite{pickens2020mapping, liu2021production,murray2022high}.
Moreover, the decision rule of data-driven approaches is generally not clear and their results are biased towards the training data distribution \citep{moreno2012unifying}, which makes quality control of their output difficult.

\subsection{Airborne LiDAR based water mapping}

Relatively very few studies have used topography airborne laser scanning for surface water mapping. Topographic airborne laser scanning refers to the system operated with the near infrared laser. We exclude bathymetry LiDAR in our discussion. The main idea of surface water mapping with LiDAR is to take advantage of the intensity values of airborne LiDAR and the property that the point density of LiDAR point clouds within water areas is much lower than in non-water areas as water hardly reflects the laser point.

\citep{brzank2008aspects} performed a supervised fuzzy classification using using height, intensity, and point density as input. \citep{hofle2009water} modeled both intensity and point density based on the sensing geometry and performed a seeded region growing segmentation and object-based classification for surface water mapping. \citep{smeeckaert2013large} used a support vector machine for water extraction based on historical coastline data. This study refined the coastline data and conducted a supervised classification using 3D point-based features. \citep{crasto2015lidar} performed a decision tree-based data-driven approach using point density, elevation, and intensity of LiDAR points as inputs. More recently, \citep{malinowski2016local} input radiometric and geometric laser point features of full-waveform LiDAR into supervised classifiers to conduct flood mapping. \citep{shaker2019automatic} input multispectral LiDAR-based features into the log-likelihood classification model and performed supervised classification for water mapping. Also, \citep{yan2019scan} performed water body mapping by developing a method that collects training samples based on the ratio of intensity and elevation for each scan line and performs a classification using the maximum likelihood model.

To sum up, state-of-the-arts LiDAR-based water mapping methods exploit the properties of water which are characterized by the low intensity and low point density over water bodies. Also, many studies tried to calibrate the intensity values. However, intensity values are easily confounded by many other factors such as sensor orientation, scan angles, and turbidity and roughness of water surfaces, which require even more concrete calibration than that for the optical images. In the end, LiDAR-based water mapping methods ends up relying on data-driven approaches which require human interventions and high-quality training data; while the performance of data-driven approaches highly depends on the quality of training data and is susceptible to data distribution shifts, just like the case of optical image-based methods.

\subsection{Limitations and opportunities}
A major barrier to developing scalable water body mapping in both optical and LiDAR-based methods has been an inconsistent spectral characteristic of water in remote sensing data (reflectance in optical imagery; intensity in airborne LiDAR).
This problem becomes more challenging in high-resolution image, it will involve many more small water bodies whose spectral properties are more diverse and more influenced by the surrounding environment.
In addition, large-area mapping requires multiple image acquisitions because individual swath width is limited, especially in high spatial resolution imaging; while multiple image acquisitions increase the likelihood of involving different atmospheric conditions or may involve different types of sensors, which will eventually make scalable water mapping much more difficult.

In the meantime, topography is defined by terrain and water, and they interact and determine each other's shape. For the complete topography mapping, both terrain and water maps are needed. In addition, as discussed, we need a three-dimensional terrain to make surface water map truly impactful. This is why airborne laser scanning is necessary for the high-resolution 3D hydrography mapping. National-scale airborne LiDAR data has been widely collected nowadays. It gives us a great opportunity for large-area mapping at high-resolution. However, since the performance of water mapping using LiDAR is not as effective as using optical imagery, airborne LiDAR data has been used only for terrain mapping and has not been effectively used for water mapping.

In this paper, we developed a scalable water mapping workflow that only requires topographic airborne LiDAR data. The workflow uses only the 3D coordinates of the earth and leverages the geometric and signal properties of water from topographic airborne LiDAR data. As the intensity data is not used, any sensor-specific calibration is not needed, and uncertainties from spectral variability can be avoided. Our workflow extracts water segments based on the LiDAR's point density and extends the extent of each segment based on the robust assumption that the local surfacd water is flat. With the proposed method, it is now possible to obtain full 3D topographic map in a fully automated manner.

\section{Proposed water mapping method}
\subsection{Overall strategy}
\label{sec:Finely rasterized DSM generation}
A standard LiDAR sensor for topographic scanning emits near-infrared laser pulses whose most of the energy is absorbed by the water or specularly reflected away from the receiver's field of view. Therefore, water bodies can be characterized by a large number of laser point dropouts in topographic airborne LiDAR data. In other words, local point density of LiDAR point clouds can be a variable that distinguishes water from non-water. However, the tricky part is that even within the same water body, local point densities can vary greatly, and there are non-water regions where the local point density is lower than that of a water body due to occlusion. The main reason for the point density variability within water body is varying incidence angle. Theoretically, the number of recorded laser points should gradually decrease as incidence angle increase. However, not only the surface conditions on the ground and water, but also the sensor orientation due to the physical movement of the airborne platform greatly influences the local point density over the surface. To overcome the limitation of solely relying on point density, previous studies have utilized intensity information. However, intensity values are also greatly affected by sensing conditions, which could make scalable water body mapping more difficult.

The proposed method simply uses only 3D coordinate information of point clouds. It does not used the intensity information.
Instead, we add an assumption that local surface water is flat. This assumption is robust as the way of water is always downwards and locally connected.
Figure 1 shows a pictorial example of typical airborne laser scanning.
\begin{figure}[b]
	\centering

	\includegraphics[width=0.45\textwidth]{figs/water_profile.jpg}
	\caption{Profile view of an airborne laser scanning}
	\label{fig:fig1}
\end{figure}
The red dot represents points that are reflected and observed by the receiver. Only red dots' 3D coordinates are the input of the proposed method. As shown, low point density can either indicates a water body or an occluded area by buildings. Therefore, the proposed method first considers the area with low point density as the initial water segment but excludes the false water segment caused by occlusion of high-rise buildings. On the other hand, some regions of the water body can have as a high point density as the ground due to the low incidence angle, which makes the initial water segment not comprehensive. Thus, the proposed method extends the initial water segment using the water elevation-based region merging. Specifically, the proposed method extends the initial water segment by merging other water segments whose elevations are the same as the initial water segment and that are connected to the initial segment at the same time. To sum up, the proposed method consists of the following two steps (1) initial water segment extraction, and (2) water elevation-based region merging (WERM).

\subsection{Initial water segments extraction}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/water_workflow_15.jpg}
	\caption{A graphical illustration of the proposed method}
	\label{fig:fig2}
\end{figure*}

In the proposed method, initial water segments are extracted based on the assumption that the point density over water bodies is much lower than in non-water areas as water bodies hardly reflect laser points. Please refer to “1. Water segments extraction” of Figure 2.

First, the proposed method rasterized the raw airborne LiDAR data into a high-resolution 2D gridded space, digital surface model (DSM). If the point clouds are rasterized into high-resolution 2D gridded space, there will be pixels where any laser point is not registered, and areas with more non-registered pixels are more likely to be water bodies. Our method calculates the average point density ($P$) of a DSM by dividing the number of non-empty cells by the number of total cells. Then, the number of non-empty cells in a certain size of a sliding window ($N pixels$) will have a binomial distribution $B(N, P)$. Specifically, $B(N, P/2)$ was used to compensate for the imbalance of point density due to scanning overlap and to avoid overly detecting water. Based on the binomial distribution, a lower confidence bound was used for the decision boundary of the water classification. Simply put, the center of sliding window will be classified as a water pixel if the sliding window contains relatively fewer numbers of laser points than other windows.

As a default, we rasterized the point clouds into a 0.5-meter resolution grid, and we used a window of 9 x 9 and a confidence level value of 2 (the critical z-score) as a default. In addition, as shown in Figure 1, as lower local point density can be due to an occlusion from adjacent high-rise buildings, the proposed method adopts a 3D building mapping algorithm \citep{song2022towards} as an optional operation to exclude false positives by creating a building buffer. This operation is only for urban areas with high-rise buildings, as non-urban areas rarely have occlusion large enough to be extracted by water.

\subsection{Water Elevation-based Region Merging (WERM)}

As the initial water segment is not comprehensive, the proposed method extends the initial water segment based on the assumption that the elevation of the connected water segment is the same. Please refer to “2. Water Elevation-based Region Merging (WERM)” in Figure 2.

First, for each initial water segment, regions of the same elevation are extracted by slicing the DSM based on the elevation. Then, if any extracted region is connected to the initial water segment, the region is merged into the initial water body. In this way, the algorithm WERM expands initial water segments to neighboring regions whose elevations are very close to each initial water segment. Finally, all initial water segments and the merged water segment will be the final surface water map.

Specifically, in WERM, we picked the $10th$ percentile of elevations of each initial water segment as the elevation of the water segment.
As a default, WERM finds regions of the same elevation by slicing the DSM within $\pm 0.1m$ elevation range considering the vertical precision of airborne LiDAR data. Of all initial segments of water, WERM only works for bodies of water larger than 500$m^{2}$ to exclude insignificant water bodies such as small pools and puddles. Lastly, WERM runs twice as the extended water segment after the first run can give a more reliable elevation value of each segment. The default numerical values ($10th$ percentile, $\pm 0.1m$, and 500$m^{2}$) used in WERM are empirically determined and can be modified. A more detailed discussion in parameter values can be found in Section 6.

\section{Datasets and experimental design}
Extensive and various landscapes were used for the evaluation as the development of a scalable water mapping algorithm is our goal. Five different large-scale datasets were employed for the evaluation (Figure 3). The five datasets include three metropolitan areas, one mountainous area, and one coastal area.  Seas, lakes, rivers, streams, ponds, wetlands, ditches, and snows were included in the dataset. Approximately 2,500$km^{2}$ area, 10 billion pixels evaluated at 0.5-meter resolution.

\begin{figure*}
	\centering

	\includegraphics[width=0.8\textwidth]{figs/water_fivedatasets_new.jpg}
	\caption{Five large-area datasets over the US: Denver, Wind River, Dallas, Orlando, and Hollywood}
	\label{fig:fig3}
\end{figure*}

For the quantitative evaluation, we employed the water body layer of the USGS’s High Resolution National Hydrography Dataset Plus (NHDPlus HR) \cite{moore2016road, moore2019user}. NHDPlus HR is an integrated set of geospatial data layers, including the best available National Hydrography Dataset (NHD) and the National Watershed Boundary Dataset (WBD). Also, we compared the proposed method to NDWI-based thresholding methods.

Specifically, we used two different NDWI-based methods that have different levels of privileges. For the first method, we found a single optimal threshold value that makes the method obtain the highest possible overall accuracy compared to the ground truth (NHDPlus HR) for each experimental area. For the second method, we divided each dataset into small-sized tiles and found the optimal threshold value for every individual small-sized tile. Thus, the first method (“NDWI-G-Opt”) creates a surface water map using the global optimal threshold value; while the second method (“NDWI-L-Opt”) creates a surface water map that is stitched with small-sized tiles that are locally optimized to obtain the highest performance for each tile. Although the two NDVI-based methods are not feasible in real practice as the groundtruth must be unknown, we used them to evaluate the proposed method’s performance with more competitive criteria. On the other hand, for the proposed method, we used the same parameter values for the entire dataset. No site-specific optimization was carried out. In this way, we investigate the scalability and the performance of the proposed method in comparison with highly competitive NDVI-based methods to which we have given special benefits.

\begin{table*}[width=2\linewidth
\caption{Summary of the experimental datasets}\label{tbl1}
\begin{tabular*}{\tblwidth}{LLLLLL}
\toprule
Name& Denver& Wind River& Dallas & Orlando & Hollywood\\
\midrule
Location& Denver, CO & Wind River, WY & Dallas, TX & Orlando, FL & Hollywood, FL\\
\midrule
Dimension & 25-km by 25-km & 13-km by 13-km & 35-km by 35-km & 20-km by 20-km &  25-km by 2.5-km\\
\midrule
Geography & Metropolitan area & Mountainous area & Metropolitan area & Metropolitan area & Coastal area\\
\hline
& Leica TerrianMapper & Leica TerrianMapper & Leica ALS 80 & Leica ALS 80 & Riegl VG-1560i\\
& 2.4 points$/m^2$ & 2.0 points$/m^2$ & 3.0 points$/m^2$ & 9.8 points$/m^2$ & 8.2 points$/m^2$\\
Scanning details & AGL: 2.7-3.1km & AGL: 3.0-3.2km & AGL: 1.8km & AGL: 1.4km & AGL: 1.3km\\
& scan angle: 40$^{\circ}$ & scan angle: 40$^{\circ}$ & scan angle: 35$^{\circ}$ & scan angle: 40$^{\circ}$ & scan angle: 60$^{\circ}$ \\
& overlap: 20\% & overlap: 25\% & overlap: 30\% & overlap: 60\% & overlap: 30\%\\
\hline
Terrain elevation & 25-km by 25-km & 13-km by 13-km & 35-km by 35-km & 20-km by 20-km &  25-km by 2.5-km\\
\midrule
LiDAR acquisition
&
2020.05.-2020.06.
&
2019.08.-2019.09.
&
2019.03.-2019.07.
&
2018.12.-2019.03.
&
2018.06.
\\
\midrule
NAIP acquisition
&
2019.08.-2019.09.
&
2019.08
&
2020.10.-2020.11.
&
2019.11.
&
2019.11.
\\

\bottomrule
\end{tabular*}
\end{table*}

Table 1 summarizes the five datasets. Not only geographic diversity but also diversity of laser scanning conditions are taken into account for the assessment. All airborne LiDAR data were from the USGS's 3D Elevation Program (3DEP) (https://www.usgs.gov/3d-elevation-program)\cite{stoker2022accuracy}. For the NDWI-based water mapping, near infrared and green bands of the National Agriculture Imagery Program’s optical imagery (NAIP) was used.

The experimental dataset is very large. Providing averaged quantitative results for all datasets might be not sufficient to investigate the performance in-depth. Also, NHD plus datasets include snows, wetlands, and some mangroves as water class while the proposed method and NDWI-based method are intended to extract only surface water bodies.

To address these issues, tile-based evaluation was introduced in addition to the conventional quantitative evaluation to conduct more detailed evaluations in various ways. The tile-based evaluation method proceeds as follows. First, the entire experimental area was divided into small-sized tiles. Then, quantitative evaluation was conducted after excluding tiles that have significant amounts of wetlands and snow cover. In addition, we performed a separate evaluation for only selected challenging tiles based on the magnitude of disagreements among surface water maps. To be specific, we calculated a tile-based difference in intersection over unions (IoU) of the proposed method and the NDWI-based method. Then, we selected challenging areas based on the tile-based IoU difference. Lastly, we selected three sets of tiles for each dataset and investigated the difference both quantitatively and qualitatively. In this way, we tried to evaluate the performance of the proposed algorithm effectively and objectively in extensive datasets.

\section{Experimental results}
\label{sec:Experimental results}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/denver.jpg}
	\caption{Results of the Denver dataset}
	\label{fig:fig4}
\end{figure*}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/windriver.jpg}
	\caption{Results of the Wind River dataset}
	\label{fig:fig5}
\end{figure*}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/dallas.jpg}
	\caption{Results of the Dallas dataset}
	\label{fig:fig6}
\end{figure*}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/orlando.jpg}
	\caption{Results of the Orlando dataset}
	\label{fig:fig7}
\end{figure*}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/hollywood.jpg}
	\caption{Results of the Hollywood dataset}
	\label{fig:fig8}
\end{figure*}

Figure 4 shows the results for the Denver dataset. Figure 4-A shows the CIR image and extracted surface water map from the proposed method for the entire Denver dataset. OUR represents surface water maps generated from the proposed method. Figure 4-B highlighted three subset regions that showed significantly different results. CIR, NDWI, and four different surface water maps were visualized. Overall accuracy (OA) and IoU compared to USGS's surface water map were also noted. Figure 4-C shows the statistics of the three methods. C1 boxplot displays the overall accuracies of the three methods. C2, C3, and C4 show the comparison results for three different sampled sets. C2 shows the results from the top 10 challenging tiles that showed the largest IoU difference among different methods. C3 shows the results after excluding tiles whose groundtruths contain significant amounts of snow or wetland cover. Lastly, C4 shows the results from the entire tiles of the Denver dataset.

As shown in Figure 4-A, it is confirmed that the visually identifiable large water bodies are generally well-extracted. The entire 25-km by 25-km test area was divided into 100 square tiles, of which the tile showing the largest difference were enlarged and displayed in Figure 4-B. As shown in the first row of Figure 4-B, the proposed method hardly generated false positives in the urban area compared to the other two NDWI-based methods. NDWI-L-opt used the locally optimized threshold to obtain the highest possible OA in the given tile. Even though it is tuned using the USGS's water map, there are a large amount of spurious false positives in the urban area. NDWI-G-opt was able to detect some water bodies that were not extracted in NDWI-L-opt, but more false positives were detected in the urban area, resulting in lower OA and IoU.

Similarly, in the second row of Figure 4-B, the proposed method accurately extracted water bodies. However, NDWI-L-opt failed to detect some water bodies; while NDWI-G-opt generated too many false positives. Lastly, in the third row, the proposed method showed very similar results to the USGS water body. A closer look at the stream in the upper left shows that the proposed method can depict the detailed branches of the stream in more detail than the body of water in the USGS. On the other hand, the two NDWI-based methods were not able to accurately detect streams. If the NDWI threshold had been set to a lower value, the NDWI-based method would have been able to detect the stream, but false positives would have increased significantly.

As shown in boxplot C1 in Figure 4, overall accuracies are not significantly different although the proposed method has the least deviation overall. C2 shows the quantitative result by selecting only the top 10 tiles with the largest difference in IoU value among all 100 tiles. The proposed method recorded the highest accuracy in all metrics. The recall was particularly higher than other methods, indicating that the proposed method can achieve a high detection rate in small water bodies such as streams, as confirmed in the third row of Figure 4-B. The following C3 is the evaluation result using the filtered groundtruth, excluding 20 tiles that have a fairly large error in the water map of the USGS. Lastly, C4 is the evaluation result using the entire 100 tiles without filtering any tiles. In both C3 and C4, the proposed method recorded higher accuracies than NDWI-based methods overall. NDWI-L-opt obtained the highest precision. This is partly because NDWI-L-opt is tuned to get the highest OA for each tile, which makes NDWI-L-opt tend to map water bodies passively and avoid false positives.

Figure 5 shows the experimental results of the Wind River dataset. Similar to the layout in Figure 4, Figure 5-A shows the CIR image and the surface water map generated by the proposed method for the entire area, Figure 5-B shows three subset regions, and Figure 5-C shows the statistics of the experiment results.

As can be seen in Figure 5-A, the Wind River dataset is a completely mountainous terrain with no man-made structures, and steep and rugged terrains cast many shadows. Seeing the surface water map generated by the proposed method, it can be confirmed that the proposed method distinguishes between water areas and shadow areas that are difficult to distinguish visually.

Figure 5-B selects and shows three areas where there is a large difference in accuracy between water extraction methods. As shown in the first row, the proposed method clearly distinguishes the water body and the shadow caused by the terrain and records very high accuracy, but the method using the NDWI fails to distinguish between water and shadow despite applying an optimal threshold. In the NDWI-based methods, there were noises in the areas shaded by rugged terrain and trees as shown in the second and third rows. On the other hand, the proposed method generally showed very similar results to USGS and recorded very high accuracy. One error in the upper left corner of the third row was due to the snow cover.

As shown in boxplot C1 in Figure 5, the proposed method recorded higher OAs in general, and its variation was also very low. On the other hand, the NDWI-based method recorded significantly lower OA. Moreover, for all quantitative metrics in C2–C4, our proposed algorithm recorded significantly higher accuracy than NDWI-based methods.

The Dallas dataset is a 35-km by 35-km metropolitan area covering various types of water bodies. Including several large lakes, the Trinity River and streams branching off from it are widely distributed. There are also many small wetlands, streams, and ponds throughout.

Figure 6-B shows three sampled areas that show significant differences among surface water maps. The first row shows the reservoirs, streams, and lakes of the urban area. The proposed method recorded the highest values in both OA and IoU, and when visually compared with the USGS water body, the water body layer generated by the proposed method showed little error. Rather, small creeks that were not extracted in USGS were well detected by the proposed method. Meanwhile, NDWI-L-opt and NDWI-G-opt were not able to extract water bodies as well as the proposed method. In particular, NDWI-G-opt had lots of omission errors. Considering that NDWI-G-opt used a single threshold for the entire area to yield the highest OA, the result indicates that a single fixed threshold for large areas cannot guarantee satisfactory performance. In the second row of Figure 6-B, NDWI-L-opt recorded the highest overall accuracy. However, it was found that the USGS water layer has some errors probably due to the temporal discrepancy. When compared to the CIR image, the proposed method extracted small creeks well; while the two NDWI-based methods failed to detect creeks successfully. This sample shows an example that differences in quantitative metrics do not accurately represent actual performance. A similar result was also confirmed in the third row of Figure 6-B.

As shown in Figure 6-C, the overall performance of the proposed method was quite similar to that of NDWI-L-opt. In general, the proposed method tends to detect more water, considering the high recall and low precision. This might be attributed to the fact that the proposed method detects small water bodies more comprehensively as shown in Figure 6-B. Considering that NDWI-based methods find the optimal threshold using groundtruth, the result confirms the proposed method could be a better alternative way for surface water mapping.

The Orlando dataset covers a 20-km by 20-km metropolitan area. The area is dotted with hundreds of lakes and ponds over very flat terrain. As shown in Figure 7-A, lakes and ponds show various optical characteristics, which makes it challenging for image-based methods to extract accurate water bodies. However, since the proposed method exploits laser property and operates with the robust assumption that the connected water has the same elevation, it successfully extracted lakes and ponds.

Figure 7-B highlighted the three sampled regions. The first row contains several small lakes and ponds in the urban area. The proposed method showed satisfactory results when compared with CIR images and USGS's water layer, and recorded high accuracy in OA and IoU. Meanwhile, two NDWI-based methods failed to successfully extract some lakes and ponds.
This is because some non-water areas have higher NDWI values than water bodies. Lowering the threshold would have allowed all the water bodies to be extracted, but that would have generated too many false positives. NDWI-L-opt detected a larger amount of water as it used a lower threshold than NDWI-G-opt, but it costs more noise in non-water areas. The second row shows a unique characteristic of the proposed method in an algae-covered lake.
Unlike NDWI-based methods, the proposed method extracted the complete water body even where algae were covered. This is because algae-covered lakes still have the same level as other algae-free portions.

The third row shows the results of tiles where high-rise buildings and lakes coexist. As shown in the NDWI image, the shadows caused by high-rise buildings sometimes have higher NDWI values than the water body, resulting in many false positives in NDWI-based methods. On the other hand, since the proposed method uses a building buffer, most of the false positives could be prevented.

The proposed method also outperformed NDWI-based methods in quantitative metrics as shown in Figure 7-C. In challenging tiles, the performance gap was distinctive, especially in IoU, recall, and F1. This is because NDWI-based methods kept their thresholds low to retain overall accuracy high. A similar pattern is shown in C3 and C4 as well.

The Hollywood dataset covers a coastal urban area facing the ocean and is composed of various shapes of creeks, bays, and artificial waterways. As shown in Figure 8-A, the proposed method extracted water bodies well overall.

Figure 8-B shows a total of 3 tiles out of 10 tiles in the Hollywood dataset. Overall, both the proposed method and the two NDWI-based methods showed satisfactory performance. In all cases, OA and IoU were higher than 0.9 and 0.8, respectively. Compared to the USGS water body layer, the main errors were associated with the shadow that created by tall buildings along the coast. Otherwise, errors were either small-sized, insignificant inland water bodies or false negatives due to water transport vehicles parked along the waterway. Among different methods, the proposed method obtained the highest accuracy in both OA and IoU as it hardly generates errors from shadows or occlusions. One unique error observed in this dataset was a false negative due to the rapid current from the boat in the middle of a large water body. This is because the spectral characteristics of the water changed temporarily due to the rapids. Although such a small size error can be easily removed by post-processing, this case confirms the advantage of the proposed method that it is hardly affected by water surface conditions.

Figure 8-C provides quantitative results. Three excluded tiles of C3 are the bottom 3 tiles as the USGS's water map classified mangroves as water bodies. Except for those three tiles, all methods showed satisfactory results, and the proposed method obtained the highest accuracy in all metrics in all cases.

\section{Discussion}
\subsection{The impact of parameters}

\begin{figure*}
	\centering

	\includegraphics[width=1\textwidth]{figs/water_parameter.jpg}
	\caption{The impact of parameters: z-score (Z), elevation range (ER), and minimum size (MS)}
	\label{fig:fig9}
\end{figure*}

As confirmed by experimental results with extensive datasets with different laser scanning conditions, the proposed method can generate highly accurate surface water maps without additional parameter tuning. Additional experiments found that tuning the parameters within reasonable ranges had little effect on the overall statistics. Nevertheless, we have found that additional parameter tuning sometimes can produce more accurate results.

In this section, we provided experimental results after changing three key parameters: z-score (Z), elevational range (ER), and minimum size (MS). Z determines the confidence interval for finding the initial water segment. For example, a Z of 1.96 represents a 95\% confidence level from a standard normal distribution. The default value for Z is 2 in our method. We conducted additional experiments by changing Z to either 1.5 or 2.5. Thus, if Z is 1.5, it is expected that the method will detect water more sensitively. When Z is 2.5, the method will detect water more conservatively. ER is a small elevation interval for finding areas with the same elevation in the algorithm WERM. The default of ER is 0.1 (a range of $\pm$ 0.1m) considering the typical vertical precision of airborne laser scanning systems. We altered the default value to either 0.05 or 0.15 so that each case extends its initial water segment more conservatively and more aggressively. Finally, the parameter MS determines the minimum size of the initial water segment to be extended to the adjacent area. To prevent insignificant bodies of water, such as puddles, from extending to the surrounding terrain, the default is set to 500 square meters. In this parameter analysis, we experimented with 100$m^{2}$ and 1000$m^{2}$, making the algorithm expand the initial water segment more conservatively and aggressively, respectively.

Among the entire experimental areas, we excerpted one area where the difference according to the parameters was the most prominent.
Figure 9-A shows the CIR image and the surface water map generated from the proposed method with the default parameter. This area is a part of the Dallas dataset. It has the Trinity River meandering on the left and the Dallas Southside Wastewater Treatment Plant on the right. It is a very flat terrain, with a variety of water bodies, including rivers, wetlands, lakes, reservoirs, water tanks, ponds, and puddles.

Figure 9-B shows the resulting surface water map after varying three parameters' values (i.e., Z, ER, and MS). On the left side of the center, parameters are adjusted so that the water body can be extracted more conservatively, and on the right side, parameters are adjusted so that the water body can be extracted more aggressively. First, when Z was changed from the default value to 2.5 and 1.5, respectively, no significant changes occurred overall. Considering that Z works as a reference point for distinguishing water and non-water pixels using local point density, results suggest that the classification of water pixels based on the point density is quite robust. On the other hand, when the ER was changed from the default value to 0.05 and 0.15, respectively, significant changes were observed in the resulting surface water maps. First, when ER was changed to 0.05, the error as if the small pond was flooded to the surrounding terrain on the right side of Figure 9-A was improved. On the other hand, when the ER is changed to 0.15, this "flooded" error becomes wider. A similar pattern was observed when the parameter MS was changed. When MS was set to 1000, most small ponds were not extended as they were smaller than 1,000$m^{2}$. However, when is set to 100, some smaller ponds were extended to adjacent terrains.

Figure 9 shows an example of how parameters can affect the resulting surface water map. The results showed that ER and MS parameters can significantly affect the resulting surface water map. However, since the two parameters are physically meaningful parameters, as they describe water level and water area, respectively, the parameter tuning and its result are quite intuitive and easy to explain. In other words, the error and uncertainty of the proposed method can be managed relatively well compared to that of the optical image-based method.

\subsection{3D surface water map}

The proposed method can generate a 3D surface water map that includes elevation information of each water body as a by-product. 2D surface water map is very useful geospatial data to support numerous studies and decision makings. However, a 2D surface water map alone cannot say where water is flowing from and where it goes. In addition, 3D information of water is helpful to monitor the water cycle or to monitor and prevent water-related pollution and disasters. Until now, most 3D hydrology has been produced by placing a 2D water layer created from optical images on a DEM. For example, the USGS commissions and distributes the production of hydro-flattened DEMs, which requires post-processing of DEM using a separate GIS layer of water. This procedure necessarily entails additional errors due to temporal difference and registration errors.
With the proposed method, a full 3D topography (both terrain and water) can be produced with the same original data \citep{song2022new}.
We expect the proposed method can shorten complicated procedures for hydro-flattened DEM and 3D water mapping while improving their accuracies

\section{Conclusion}

This paper presents a surface water mapping method that provides reliable and accurate results at high resolution across a wide variety of landscapes. The proposed method uses only airborne LiDAR data based on the natural law that surface water is flat as gravity always pulls liquid molecules down. Avoiding the use of optical properties, our method enables fully automated and scalable surface water mapping. Extensive experiments with large ($\approx$ 2,500$km^{2}$) and diverse landscapes (urban, coastal, and mountainous areas) confirmed that our method can generate more accurate results than NDVI-based methods even without any parameter tunings. We confirmed its overwhelming performance of water extraction in varied conditions placed in diverse environments. We confirmed its advantage, especially in steep and rugged mountains with shadows and perennial snows and complex urban environments with very small water bodies. We expect the proposed method will expedite an accurate large-scale full 3D topography mapping.


\section{Introduction}
In the Aharonov-Bohm double-slit experimental setup, a long but thin solenoid of electric current is introduced between the two slits (Fig.~1). At a distance $R$ from the central axis of a long solenoid, comprising circular loops of current,  the magnetic field (in cgs units) is \cite{PU85}
 \begin{eqnarray}
\label{eq:84b.4b1a1}
\mathbf {B} &=&\frac{4\pi I M}{c} \hat z\,,\,\,\,\,\,\,\,\,\,\,\,\, R<r\\
\label{eq:84b.4b1a2}
\mathbf {B} &=&0 \,.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, R>r
\end{eqnarray}
where $I$ is the constant electric current flowing in the solenoid, with $M$ current loops per unit length along the $z$-axis, and $r$ is the radius of each circular loop.

Thus there is no magnetic field outside the solenoid.

Vector potential, defined at a field point by $\nabla \times \mathbf {A}=\mathbf {B}$, is easily determined for such a long solenoid of current, using Stokes' theorem, to get \cite{25}
\begin{eqnarray}
\label{eq:84a.4b2a}
\oint \mathbf {A} \cdot {\rm d}\mathbf {x}&=&\int( \nabla \times \mathbf {A}) \cdot {\rm d}{\bf a}=\int \mathbf {B} \cdot  {\rm d}{\bf a}\,,
\end{eqnarray}
which, from the cylindrical symmetry, possesses a finite component only along the azimuthal direction, with  ${\bf A}=A_\phi \hat\phi$, having a value at $R$ as
\begin{eqnarray}
\label{eq:84b.4b1a}
A_\phi &=&\frac {2\pi R I M}{c} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ R<r\\
\label{eq:84b.4b1b}
A_\phi &=&\frac {2\pi r^2 I M}{cR} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\
 R>r
\end{eqnarray}

An electric charge $Q$, of mass $m$ and moving with a velocity ${\bf v}$, and accordingly having a mechanical momentum ${\bf p}_{\rm m}=m{\bf v}$ (we consider here only non-relativistic cases), when passing through a region with a finite vector potential $\bf A$, gets associated with it, in addition, an electromagnetic (EM) momentum
\begin{eqnarray}
\label{eq:84a.3}
{{\bf p}_{e}}&=&\frac{Q{\bf A}}{c}\,.
\end{eqnarray}
This EM momentum is independent of either mass $m$ or velocity ${\bf v}$ of the charge $Q$, which may thus even be stationary. This form of momentum in classical  electromagnetism,  suggested by Maxwell  \cite{Ma65}, has been discussed at some length in the literature  \cite{Ko78,Se98,Gr12,Mc19a}. From Eq.~(\ref{eq:84a.3}) the momentum would be finite for a non-zero $\bf A$ at the location of $Q$ outside the solenoid of current, even if the magnetic field ${\displaystyle \mathbf {B} =\mathbf {0}}$ there.
Two equal charges,  $Q_1$  and $Q_2$, placed symmetrically on two opposite sides of the solenoid, and thus having equal and opposite vector potentials, ${\bf A}_1$, ${\bf A}_2$ with ${\bf A}_1=-{\bf A}_2$ at their locations, will accordingly, possess equal and opposite EM momentum vectors, $Q_1{\bf A}_1/c=-Q_2{\bf A}_2/c$.

In quantum mechanics, the wave function associated with the electric charge $Q$, because of the additional EM momentum ${{\bf p}_{e}}$, develops an extra phase shift
\begin{eqnarray}
\label{eq:84a.4}
{\displaystyle \varphi}=\frac{1}{\hbar}\int {{\bf p}_{\rm e}} \cdot {\rm d}\mathbf {x}=\frac{Q}{c\hbar}\int \mathbf {A} \cdot {\rm d}\mathbf {x}\,.
\end{eqnarray}
Then the two charge beams moving along two different paths, but each having the same start and end points as the other, will acquire accordingly, a phase difference
\begin{eqnarray}
\label{eq:84a.4a}
\Delta \varphi&=&\frac{Q}{c\hbar}\oint \mathbf {A} \cdot {\rm d}\mathbf {x}
=\frac{Q}{c\hbar}\int \mathbf {B} \cdot  {\rm d}{\bf a}
=\frac{Q\:\Phi}{c\hbar}\,,
\end{eqnarray}
which is thus determined by the magnetic flux $\Phi$ through the area enclosed between the paths, even though at the location of either beam the magnetic field is nil. The choice of gauge for ${\bf A}$ does not affect $\Delta \varphi$ \cite{29}. Once Eq.~(\ref{eq:84a.3}), from classical  electromagnetism, is accepted to be providing a measure of EM momentum in the system, everything else seems to follow from quantum mechanics.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig1.jpg}
\caption{A schematic of the Aharonov-Bohm experimental setup. Coherent beams of electrons, passing  through two slits, $S1$ and $S2$, separated along the y-axis, form an interference pattern on the screen. When a long but thin solenoid of electric current is introduced behind the two slits, midway between the two beams, an extra phase between the two electron beams appears, as shown by a shift in the interference pattern. In the region of electron beams, the magnetic field is nil ($\mathbf {B}=0$) though the vector potential is finite ($\mathbf {A}\ne 0$).}
\end{center}
\end{figure}

Such a phase difference has actually been inferred in the Aharonov-Bohm setup, from an observed shift in the interference pattern \cite{Ch60,To86,Os86}.
From a physical perspective, however, it is not clear from where does such a mysterious EM momentum appear, whose presence has been verified experimentally, implying that there is some interaction between the solenoid and the charge beams even in the absence of any electromagnetic field at the locations of the beams.
This interaction, ostensibly through the vector potential at the location of either beam, is responsible for EM momenta that get reflected in the observed quantum mechanical phase shift between two beams.
But to date it still remains a mystery how the solenoid influences the beams in the absence of  EM fields at their locations, or how do the beams interact with the solenoid.

The presence of an EM momentum to have a classical origin in the system has been attempted for a {\em moving} charge  \cite{Bo08}. However, according to Eq.~({\ref{eq:84a.3}) one should be able to account for the EM momentum in the system {\em even for a stationary charge} in a classical explanation, which somehow has not been successful \cite{Mc19b}.
The literature on both experimental and theoretical fronts is so vast that we make reference to a recent review article \cite{Ba09}.
The failure of an explanation within  classical physics  has led to quantum mechanical topological explanation where the Aharonov-Bohm phase of the wave function of a charged particle depends on the
topology of the space it moves in, assuming that the presence of a solenoid of current makes the configuration space non-simply connected  \cite{Ah83,Ah16,Ah17}. Such non-locality features of quantum mechanics may have deep philosophical implications  \cite{16,He97}; classical explanations of EM momentum are not in vogue in the contemporary literature.

\section{EM momentum of a charge in a vector potential from classical perspective}
Since Eq.~({\ref{eq:84a.3}), relating the momentum to the vector potential, has its genesis in classical electromagnetism, one should be able to resolve the issue of EM momentum, without invoking quantum effects, within the classical physics itself.
A convincing argument that $Q{\bf A}/c$ in some respects does represent momentum within the classical electromagnetism, comes from the following.

Suppose the electric current $I$ in the solenoid is slowly decreased at a constant rate, Then from Eq.~({\ref{eq:84b.4b1b}), the vector potential $\mathbf {A}$ would decrease too,
giving rise, in the absence of a scalar potential, to an electric field $c \mathbf {E}=-\partial\mathbf {A}/\partial t$.
This electric field would exert a force on the stationary $Q$, changing its mechanical momentum
\begin{eqnarray}
\label{eq:84a.5a}
\frac{{\rm d}{\bf p}_{\rm m}}{{\rm d}t}&=&m\frac{{\rm d}{\bf v}}{{\rm d}t}=Q\mathbf {E}\nonumber\\
&=&-\frac{Q}{c}\frac{{\rm d}{\bf A}}{{\rm d}t}
=-\frac {2\pi Q r^2 M}{c^2R}\frac{{\rm d}{I}}{{\rm d}t} \hat\phi
\,,
\end{eqnarray}
implying thereby that the mechanical momentum $m{\bf v}$ of $Q$ increases at the cost of $Q{\bf A}/{c}$, with the latter  decreasing for ${\rm d}{I}/{\rm d}t <0$.
This suggests
$Q{\bf A}/c$ to be a form of momentum, called EM momentum ${\bf p}_{\rm e}$ of the charge $Q$, along with conservation of ${\bf p}_{\rm m}+{\bf p}_{\rm e}=m{\bf v}+{Q{\bf A}}/{c}$, called generalized momentum \cite{Go50}, as
\begin{eqnarray}
\label{eq:84a.6}
\frac{\rm d}{{\rm d}t}\left[{\bf p}_{\rm m}+{\bf p}_{\rm e}\right]&=& 0\,,
\end{eqnarray}
all within classical electromagnetism.

However, to date, it still remains an unresolved enigma where after all could such a `mysteriously hidden' momentum be residing since there is no obvious net linear motion in the system, especially when the charge is considered to be stationary.
As the term ``momentum'' conjures up a vision of some kind of linear motion, a question arises where, after all such motion, if any, is lying in the system? The only non-random motion ostensibly present in the system is in the drift velocities of the current carrying charges in the steady current loop. However, a linear momentum cannot be solely due to the drift velocities, as any such momentum vector integrated over a closed circuit would be zero. Moreover, from Eq.~(\ref{eq:84a.3}), the momentum in question involves, not just the electric current that gives rise to ${\bf A}$, but also the specification of charge $Q$ and {\em its location} (where ${\bf A}$ is to be evaluated), although any movement of $Q$ does not enter into picture.

We explore here accordingly, from a classical physics perspective, the mysterious momentum in the Aharonov-Bohm setup, endeavouring to possibly unravel wherein the momentum lies in the system. We shall demonstrate here explicitly how an electric charge stationary at a location, where there may be a finite vector potential ${\bf A}$, though no magnetic field, does give rise to EM momentum, mysteriously latent in the system and which is reflected in the Aharonov-Bohm experiments.
Moreover, as will be seen, the momentum in question is not confined to and localized at some specific location, like that of charge $Q$; the non-local characteristic of momentum is evident  in such a case even within the classical physics picture itself.

In the case of $N$ discrete charges $q_{\rm j}$, with velocity vectors ${\bf v}_{\rm j}$, the vector potential ${\bf A}$ at the  location of the charge $Q$ is computed from the summation
\begin{eqnarray}
\label{eq:84a.4e1}
{\bf A}({\bf x}_{0})&=&\frac {1}{c}\sum_{j=1}^{N} \frac {q_{\rm j}\,{\bf v}_{\rm j}}{|{\bf x}_{\rm j}-{\bf x}_{0}|}\,,
\end{eqnarray}
where $|{\bf x}_{\rm j}-{\bf x}_{0}|$ is the distance of charge $q_{\rm j}$, moving with velocity ${\bf v}_{\rm j}$, from the location ${\bf x}_{0}$ of the charge $Q$.

Then we have
\begin{eqnarray}
\label{eq:84a.4e2}
\frac {Q {\bf A}}{c}&=&\frac {1}{c^2}\sum_{j=1}^{N} \frac {Q\, q_{\rm j}}{|{\bf x}_{\rm j}-{\bf x}_{0}|}\,{\bf v}_{\rm j}\,.
\end{eqnarray}

In order to comprehend the electromagnetic momentum in the system from a physical perspective, instead of the usual way of looking at it in terms of the vector potential $A$ at the location ${\bf x}_0$ of the charge $Q$, we can interpret Eq.~(\ref{eq:84a.4e2}) in terms of the scalar potential $\phi_{\rm j}=Q/|{\bf x}_{\rm j}-{\bf x}_{0}|$, due to $Q$ at the location ${\bf x}_{\rm j}$ of charge $q_{\rm j}$.
For that we rewrite Eq.(\ref{eq:84a.4e2}) as
\begin{eqnarray}
\label{eq:84a.4e3}
\frac {Q {\bf A}}{c}&=&\sum_{j=1}^{N} \frac {\phi_{\rm j}\:q_{\rm j}}{c^2}\,{\bf v}_{\rm j
\,.
\end{eqnarray}

We can use the energy-mass relation, to express the potential energy $q_{\rm j}\,\phi_{\rm j}$ of a charge $q_{\rm j}$, owing to the presence of charge $Q$ at ${\bf x}_0$, in terms of its mass equivalent
\begin{eqnarray}
\label{eq:84a.4d4}
\Delta m_{\rm j}&=&\frac {q_{\rm j}\,\phi_{\rm j}}{c^2} =\frac {1}{c^2} \frac {Q\, q_{\rm j}}{|{\bf x}_{\rm j}-{\bf x}_0|}\,.
\end{eqnarray}
Then we can write
\begin{eqnarray}
\label{eq:84a.4e4}
\frac {Q {\bf A}}{c}&=&\sum_{j=1}^{N} {\Delta m_{\rm j}\,{\bf v}_{\rm j}}
\,,
\end{eqnarray}
which can now be readily recognized as an electromagnetic momentum, ${\bf p}_{e}$, in the system.
It should be noted that $\Delta m_{\rm j}$ here has nothing to do with mass $m_{\rm j}$ of the $j$th charged particle and that the electromagnetic momentum ${\bf p}_{e}$ in Eq.~(\ref{eq:84a.4e4}) is not sum of the kinetic momentum, $\Sigma m_{\rm j}\,{\bf v}_{\rm j}$, of moving charged particles. Also it is not possible to localize the potential energy $q_{\rm j}\,\phi_{\rm j}$ or the equivalent mass $\Delta m_{\rm j}$ at either of the charge locations, ${\bf x}_{0}$ or ${\bf x}_{\rm j}$. Nor could one pinpoint the electromagnetic momentum ${\bf p}_{e}$ at the location ${\bf x}_{0}$ of charge $Q$, one has to instead take a holistic view that the system comprises an electromagnetic momentum ${\bf p}_{e}$, without localizing it, even from a classical physics perspective.

For a continuous distribution of moving charges or a current density ${\bf j}({\bf x})=\rho({\bf x})\,{\bf v}({\bf x})$, the vector potential $\bf A$ at a field point ${\bf x}_0$ is determined from the volume integral  \cite{1,2,25,PU85}
\begin{eqnarray}
\label{eq:84a.4b}
{\bf A}({\bf x}_0)&=&\int\frac {{\bf j}({\bf x})}{c|{\bf x}-{\bf x}_0|}\,{\rm d}\tau=\frac {1}{c}\int\frac {\rho({\bf x})\,{\bf v}({\bf x})}{|{\bf x}-{\bf x}_0|}\,{\rm d}\tau\,,
\end{eqnarray}
where $|{\bf x}-{\bf x}_0|$ is the distance from the point ${\bf x}_0$ of the charge element $\rho({\bf x})\,{\rm d}\tau$, moving with a velocity $\bf v({\bf x})$, here ${\rm d}\tau$ denotes an element of volume.

Then we can write
\begin{eqnarray}
\label{eq:84a.4d1}
\frac{Q{\bf A}}{c}=\frac {1}{c^2}\int\frac {Q\,\rho({\bf x})\,{\bf v}({\bf x})}{|{\bf x}-{\bf x}_0|}\,{\rm d}\tau\,,
\end{eqnarray}
Because of the scalar potential $\phi({\bf x})=Q/|{\bf x}-{\bf x}_{0}|$ at ${\bf x}$ due to $Q$, the system comprising a charge density $\rho({\bf x})$, possesses $\rho({\bf x})\,\phi({\bf x})$ as potential energy per unit volume. Then in the expression
\begin{eqnarray}
\label{eq:84a.4c1}
\frac{Q{\bf A}}{c}&=&\frac {1}{c^2}\int {\phi({\bf x})\,\rho({\bf x})\,{\bf v}({\bf x})}\,{\rm d}\tau\,,
\end{eqnarray}
we could use the energy-mass relation, to express the potential energy density in terms of its equivalent mass density, $\mu({\bf x})={\rho({\bf x})\phi({\bf x})}/{c^2}$, to write
\begin{eqnarray}
\label{eq:84a.4d2}
\frac{Q{\bf A}}{c}&=&\int{\mu({\bf x})\,{\bf v}({\bf x})}\,{\rm d}\tau={{\bf p}_{e}}\,.
\end{eqnarray}
Here $\mu({\bf x}){\bf v}({\bf x})$ is the momentum density in the system,
whose volume integral yields the EM momentum ${{\bf p}_{\rm e}}$ in the system, owing to the presence of charge $Q$ at ${\bf x}_0$.
\section{Electromagnetic momentum of a charge outside a solenoid of current}
We want now to examine the EM momentum of a charge outside a solenoid and we shall show here that Eqs.~(\ref{eq:84a.4d1}), (\ref{eq:84a.4c1}), or equivalently (\ref{eq:84a.4d2}), lead to an electromagnetic field momentum, latent in the system. We should clarify that this EM momentum is different from what sometimes is referred to as `hidden momentum'  \cite{Sh67} and which actually is a mechanical momentum present in the system.

A long, thin solenoid, carrying a steady electric current $I$, can be considered as a superposition of large number of small planar current loops, each carrying current $I$, with, say, $M$ loops stacked per unit length of the solenoid, plus  a long straight wire carrying current $I$ along the axis of the solenoid. Inside the solenoid the magnetic field would still be $B_{\rm z}=4 \pi I M/c$ (Eq.~({\ref{eq:84b.4b1b})). On the outside, however, there will be an azimuthal field $B_{\phi}= {2I}/{cR}$ due to current $I$ along the axis of the solenoid \cite{PU85}, which nonetheless, would not affect $\Delta \varphi$ (Eq.~({\ref{eq:84a.4a})), since flux $\Phi$ enclosed between the two beams will not change. Therefore we shall henceforth ignore the axial current in the solenoid

A small planar loop carrying a current $I$ around area $\bf a$ of the loop constitutes, irrespective of its shape, a magnetic dipole ${\bf m}=I{\bf a}/c$, giving rise to a vector potential \cite{PU85} ${\bf A}={{\bf m}\times\hat{\bf R}}/{R^2}$ at ${\bf R}$ from the loop. This implies for the charge $Q$ at ${\bf R}$, from Eq.~({\ref{eq:84a.3}), an EM momentum
\begin{eqnarray}
\label{eq:100.4}
\!\!\!\!\!\!\!\!\!{\bf p}_{\rm e} &=&\frac{Q\,{\bf A}}{c}=\frac{Q\,{\bf m}\times\hat{\bf R}}{cR^2}=\frac {({\bf E} \times {\bf m})}{c}=\frac {({\bf E} \times {\bf a})\,I}{c^2},
\end{eqnarray}
where ${\bf E}=-Q\hat{\bf R}/{R^2}$ is the electric field due to the charge $Q$  at the location of the small current loop. Thus Eq.~({\ref{eq:84a.3}), expressing EM momentum due to the vector potential of the current loop at the location of $Q$, represents implicitly a mutual interaction of the charge $Q$ and the current loop since from Eq.~({\ref{eq:100.4}), the EM momentum vector could as well be considered due to the cross product of electric field ${\bf E}$ of $Q$ and magnetic moment ${\bf m}$ of the current loop. Of course this does not in any way  resolve the issue of momentum in this apparently static system.

Since the current loop may consist of a conducting wire, the electric field of the charge $Q$ will not extend inside the loop wire, as the induced surface charge density there would tend to cancel any external static electric field inside the wire, leaving only the perpendicular components at the surface, to make the loop equipotential. In that case the system, in the absence of mutual interaction of charge $Q$ and the current loop, may not possess an EM momentum \cite{Jo94}, contrary to what would have been otherwise expected from Eq.~({\ref{eq:84a.3}). However, the EM momentum inferred for the system from the experimentally  observed shift in the fringe patterns \cite{Ch60,To86,Os86} indicates that there may be  something amiss in the above arguments.

Actually, there is a rather subtle issue involved here as in these experiments there are two coherent charged beams, emerging simultaneously from slits $S1$ and $S2$, assumed to be symmetrically placed on either side of the thin solenoid or its small current loops. Thus one has to consider the EM momentum simultaneously for a pair of equal charges, placed symmetrically, on two opposite sides of the current loop. To be specific, we designate the charges as  $Q_1$ from slit $S1$ and $Q_2$ from $S2$, lying respectively at distances ${\bf R}_1$ and ${\bf R}_2$, with ${\bf R}_2=-{\bf R}_1$, measured from the loop position.
Thus, at least to a first order, the electric fields, say ${\bf E}_2$ and  ${\bf E}_1$ at the loop and the corresponding scalar potentials across the loop will be equal and opposite for the two charges, making the loop effectively equipotential, without the aid of induced surface charges that would otherwise get formed there.

Since the induced surface charge on the conducting coil due to either charge may be nil, the mutual interaction between the current loop and each charge would still be present, with the resulting momentum associated with either charge being equal and opposite to that associated with the other charge,
${{\bf p}_{e}}_2=-{{\bf p}_{e}}_1$, which show up as a shift in the interference pattern.
This can be seen from the conservation of generalized momentum. Suppose the electric current $I$ in the solenoid is slowly reduced at some constant rate. Then the induced electric field vector (Eq.~({\ref{eq:84a.5a}))
at ${\bf R}_1$ would exert force on $Q_1$, changing its mechanical momentum. Now  any such change in the mechanical momentum of a charge is possible, from the conservation of generalized momentum (Eq.~({\ref{eq:84a.6})), only at the cost of its EM momentum, $\Delta {\bf p}_{\rm m1}=-\Delta {\bf p}_{\rm e1}$, implying the presence of ${{\bf p}_{e}}_1$ in the system. Similar is the argument for ${{\bf p}_{e}}_2$ in the case of $Q_2$.

Thus one could proceed with the investigation of EM momentum associated with either charge, without considering any induced surface charges in the current loop that otherwise might have formed to cancel the electric fields, thus leaving the mutual interaction between corresponding charge and the current loop unperturbed.
\subsection{Electromagnetic momentum of a charge outside a small current loop}
We apply our above results to show, from a physical perspective, that in the presence of an external charge, a current loop does give rise to an EM momentum.
For this, we consider a small rectangular loop, carrying a steady current, $I=j \sigma=n e  v_{\rm d} \sigma $, where $n$ is the number density of charges in the circuit,
$e$ is the electric charge of conducting charged particles, $v_{\rm d}$ is their drift velocity and $\sigma$ is the cross-section of the current-carrying wire. Then from Eq.~(\ref{eq:84a.4e1}) or (\ref{eq:84a.4b}), the vector potential $\bf A$ at the location ${\bf x}_0$ of $Q$ is determined from the integral over the circuit length
\begin{eqnarray}
\label{eq:84a.4b2b}
{\bf A}({\bf x})&=&\frac {1}{c}\oint\frac {n e \sigma}{|{\bf x}-{\bf x}_0|}{\bf v}_{\rm d}\,{\rm d}{l}\,,
\end{eqnarray}
where ${\rm d}{l}$ is an element of  the circuit, with ${\bf v}_{\rm d}$ as the drift velocity along  the direction of the current at the location ${\bf x}$ of the circuit.
Then the EM momentum of the system, from Eq.~(\ref{eq:84a.4d1}), can be written as
\begin{eqnarray}
\label{eq:84a.4c}
{{\bf p}_{e}}&=&\oint\frac{{Q}\,{n e \sigma}{\bf v}_{\rm d}}
{ {c^2}|{\bf x}-{\bf x}_0|}{\rm d}{l}=\oint\frac{{\bf v}_{\rm d}}{c^2}{\rm d}\,{\cal E}=\oint{\bf v}_{\rm d}\,{\rm d}{m}.
\end{eqnarray}
Here ${\rm d}{\cal E}={Q}\,{n e \sigma}{\rm d}{l}/|{\bf x}-{\bf x}_0|$ is the electric potential energy of current carrying charges in the volume element $\sigma{\rm d}{l}$ in the presence of charge $Q$ at ${\bf x}_0$, while  ${\rm d} m={\rm d}{\cal E}/(c^2)$ is the equivalent electric mass element. The momentum in Eq.~(\ref{eq:84a.4c}) is not the same as the kinetic momentum, $\oint m n \sigma {\bf v}_{\rm d}\,{\rm d}{l}$, of the electric current carriers, each supposedly of individual mass $m$. Such a kinetic momentum of the electric current carriers for a steady current, in any case, yields a nil value when summed over the whole current loop ($ m n \sigma \oint{\bf v}_{\rm d}\,{\rm d}{l}=0$), however, $\oint {\bf v}_{\rm d}\,{\rm d}{m}$ over the loop does result in a finite value, as will be shown below.

\subsection{Electromagnetic momentum of charges placed on symmetrically opposite sides a small current loop}
We have still to show that the pair of equal charges, $Q_1$ and $Q_2$, placed symmetrically, on two opposite sides of the current carrying solenoid will give rise to equal and opposite EM momentum in the system, something not so readily apparent from Eqs.~(\ref{eq:84a.4e4}), (\ref{eq:84a.4d2}) or (\ref{eq:84a.4c}).

Let us consider charge $Q_1$ at a distance $R_1$ along the $-y$ direction, from a small rectangular current loop (Fig.~2). We can calculate the EM momentum of the system by considering pairs of current elements, each of infinitesimal length ${\rm d}{l}$, placed symmetrically on two opposite sides of the loop.
The drift velocity in the arm $DA$ of the current loop is along $\hat{\bf x}$ direction and the current element is nearer to the charge $Q_1$, thus making a higher contribution to the EM momentum which is along $\hat{\bf x}$ direction, while the drift velocity in the arm $BC$ is along $-\hat{\bf x}$ direction and the current element being farther from $Q_1$, makes a contribution which is relatively lower and is along $-\hat{\bf x}$ direction. The EM momentum of the whole  system, for a small current loop, with $r\ll R_1$ (Fig.~2), is
\begin{eqnarray}
\label{eq:100.1}
{\bf p}_{\rm e1} &=& \frac{{Q_1}\,{n e \sigma}\,s\,{v}_{\rm d}}{{c^2}(R_1-r)} \hat{\bf x}-\frac{{Q_1}\,{n e \sigma}\,s\,{v}_{\rm d}}{{c^2}(R_1+r)} \hat{\bf x}
\,.
\end{eqnarray}
The first term on the right hand side is the contribution of the current loop from the arm DA  while the second term is from the arm BC, each arm of length $s$.
The contributions from the arms $AB$ and $CD$ to the EM momentum, being equal and opposite, cancel.

To a first order (for $r/ R_1\ll 1$), we get
\begin{eqnarray}
\label{eq:100.2}
{\bf p}_{\rm e1} &=& \frac{{2Q_1}\,{n e \sigma}\,s\,{v}_{\rm d}\,r}{{c^2}R_1^2} \hat{\bf x
=\frac {({\bf E}_1 \times {\bf a})\,I}{c^2}\,,
\end{eqnarray}
where ${\bf E}_1=Q_1{\hat{\bf y}}/{R_1^2}$, is the electric field of the charge $Q_1$ at the location of the small current loop and ${\bf a}=2\:{\bf s}\times {\bf r}=2s\,r\, \hat{\bf z}$ is the area vector of the current loop.
The EM momentum in the system, being directly proportional to the drift velocity (${p}_{\rm e}\propto {v}_{\rm d}$) of current carriers, is thus zero to start with for ${v}_{\rm d}=0$ (when $I=0$) and increases as ${v}_{\rm d}$ increases with $I$.
It is, however, interesting that a finite linear momentum exists in the system owing  to the presence of charge $Q_1$ outside the current loop, even when $Q_1$ as well as the current loop are both stationary, and momentum vector, if any, due to current carrying charges adds to zero over the closed loop ($ m n \sigma \oint{\bf v}_{\rm d}{\rm d}{l}=0$).

As for the charge $Q_2$, lying on opposite side of the current loop at distance $R_2$ along the $y$ direction, the contribution of EM momentum from the arm BC of the current loop will be higher than that of the arm DA, as a result the net EM momentum will be along $-x$ direction. It will be so even though the potential energy and the mass equivalent for  $Q_2$ is similar as for  $Q_1$, the opposite directions of drift velocities in arms  DA and BC will make ${\bf p}_{\rm e2}=-{\bf p}_{\rm e1}$.

The electric current in the loop could be due to electrons instead of positive charges, essentially making no difference to any of our arguments. Also, at the locations of charge $Q_1$ or $Q_2$, the electric field as well as the scalar potential of the lattice of positive ions, is equal and opposite to that of negative current carrying electrons for an over all charge neutral current loop, however, that does not cancel the EM momenta ${{\bf p}_{e1}}$ and ${{\bf p}_{e2}}$, arising from the drift velocities of current carrier electrons as the positive ions, fixed in the lattice, do not move with any drift velocity. There is though no net electric potential energy in the system due to $Q_1$ or $Q_2$, nonetheless it gives rise to a net EM momentum.
\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig2.jpg}
\caption{A charge $Q_1$ lies a distance $R_1$ from a small current loop $ABCD$, carrying a uniform current $I$. Two equal current elements, ${\rm d}{l}$,
separated by distance $2r$, on two opposite sides of the current loop, are shown.}
\end{center}
\end{figure}

We could replace the rectangular shape of the loop with a polygon comprising a larger number of sides or in the limit even with a circular loop without changing the final result. In fact, by a  side by side superposition of a sufficiently large number of small rectangular loops, any such shape of the current loop could be realized. Further, the electric field $\bf E$ now need not necessarily be only in the plane of the current loop. Equation~(\ref{eq:100.2}), therefore, is a fairly general result as long as the current loop is small enough for $\bf E$ to be considered uniform over its extent.
\subsection{An alternate computation of electromagnetic momentum for current loop in electric field of charge}
The presence of EM momentum (Eq.~(\ref{eq:100.2}) in the system, can be worked out in an alternate way, without starting from Eq.~(\ref{eq:84a.3}). In fact, along with the computation of the EM momentum, one could even derive Eq.~(\ref{eq:84a.3}) as well, that way.
\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig3.jpg}
\caption{A current loop $ABCD$ is carrying current $I$. There is an electric field ${\bf E}$, uniform over the current loop. The electric field ${\bf E}$ does positive  work, $I {\bf E} \cdot {\rm d}{\bf l}$ on a current element ${\rm d}{\bf l}$ in the arc $ABC$ of the current loop and an equal and opposite work on a similar current element ${\rm d}{\bf l}$ in the arc $CDA$, the two current elements separated by a distance $S$, on two opposite sides of the current loop. This implies transfer of an infinitesimal electromagnetic energy, between two current elements across the shaded area ${\rm d}{\bf a}$ of the circuit, which represents an element of electromagnetic momentum in the system.}
\end{center}
\end{figure}

We consider a small circular loop, $ABCD$, carrying a constant current $I=j \sigma$, with $\bf {j}$ as the current density and $\sigma$ the cross-section of the current-carrying wire. There is an electric field $\bf E$, uniform over the dimensions of the current loop (Fig.~3).
The electric field does work on the current density ${\bf j}$ at a rate ${\bf j}\cdot {\bf E}$ per unit volume  \cite{1,2,25}.
Accordingly, on a current element d$\bf l$ of the loop, the electric field does work
at a rate ${\bf j}\cdot {\bf E}\,\sigma\, {\rm d}{l}= I\, {\rm d}{\bf l \cdot E}$.
In section $ABC$ of the loop,  the electric field
$\bf E$ is doing positive work while in section $CDA$, the work done is negative.

Thus, power is fed through the electric field into the system through the side $ABC$ while a similar amount of power is being drained off from the side $CDA$.  Effectively, a continuous transfer of electric energy per unit time is taking place from arm $CDA$ to $ABC$ across the current loop due to the presence of the electrical field (Fig.~3). Even Though in this process there is no net
change in the total energy content of the system, there is nevertheless a
linear momentum associated with this energy flux across the loop from side $CDA$ to $ABC$.

The electromagnetic momentum can be calculated easily if we consider a pair of current
elements placed symmetrically on two opposite sides of the loop (Fig.~3). The power being fed into the circuit element ${\rm d}{\bf l}$ in the side $ABC$
through the electric field is $I {\bf E} \cdot {\rm d}{\bf l}$, while a similar amount of power is being drained off from a similar circuit element${\rm d}{\bf l}$ in the side $CDA$. Effectively, an energy
$I {\bf E} \cdot {\rm d}{\bf l}$ is being transported per unit time along $\textbf{s}$, the vector joining the  current element in side $CDA$ to that in $ABC$, implying a momentum
\begin{equation}
\label{eq:p31.1a}
{\rm d}{\bf  P}=I({\bf E} \cdot {\rm d}{\bf l}) \frac {\bf s}{c^2}= {\bf E} \times {\rm d}{\bf a}\frac {I}{c^2},
\end{equation}
where d$\bf a={\bf s}\times {\rm d}{\bf l}$ is the element of area vector contained between these two current elements across the loop.
An integration over the entire loop gives the total momentum
\begin{equation}
\label{eq:p31.1b}
{\bf  P}_{\rm e}=\frac {({\bf E} \times {\bf a})\,I}{c^2}= \frac{\bf E \times \bf m}{c}=\frac{Q{\bf m}\times\hat{\bf R}}{R^2c}=\frac{Q {\bf A}}{c}\,,
\end{equation}
Thus in this way we get not only the EM momentum, which is the same as in Eq.~(\ref{eq:100.2}), we also derive Eq.~(\ref{eq:84a.3}).

There are other, similar examples in physics where momentum is present in the system due to equal and opposite work being done on spatially separated parts of the system, implying a flux of energy between these spatially separated parts, implying momentum in the system.
For instance, a similar continuous transport of electromagnetic energy per unit time across an electric circuit between its opposite arms, but with no change in the net energy of the system, has been shown elsewhere \cite{31} to explain the presence of linear electromagnetic momentum in a stationary system comprising a pair of crossed electric and magnetic dipoles, where nothing obviously is moving or no temporal changes are occurring in the system.
A perfect fluid under pressure, having a bulk motion even with non-relativistic velocities has finite momentum proportional to pressure that does work on two opposite ends of a fluid element giving rise to momentum in the system \cite{84}. There is also an opposite example where a charged parallel plate capacitor, moving parallel to the plate separation, has finite electromagnetic energy, but in spite of its motion, has zero electromagnetic momentum in the system \cite{15,84a}.
A moving spherical charge distribution, representing a classical electron model, where an equal and opposite work done by the opposite forces of the leading and trailing hemispheres, gives rise to an energy flux and thereby an electromagnetic momentum  in the system that explains the more than a century-old famous factor of $4/3$ in the electromagnetic momentum of such a system \cite{15,84a}.

One would normally expect the power difference between arms $ABC$ and $CDA$ to be compensated by the agency tending to maintain a uniform and steady electric current in the loop, with an equal amount of mechanical energy transfer rate from arm $ABC$ to $CDA$ in the circuit that itself might  entail a mechanical momentum
(sometimes called hidden momentum \cite{Sh67}).
However, in the present case, with equal charges, $Q_1$ and $Q_2$, on opposite sides of the current loop, a constant current in both arms will be maintained because of their equal and opposite electric fields, ${\bf E}_1$ and ${\bf E}_2$. Moreover, the energy flux from arm $CDA$ to $ABC$ because of $Q_1$ will be compensated by an equal energy flux from arm $ABC$ to $CDA$ due to $Q_2$. However, as was discussed earlier, associated with individual charges there would still be present EM momenta, ${{\bf p}_{e}}_1$, ${{\bf p}_{e}}_2$, source of the phase difference, $\Delta \varphi$, experimentally observed between the two charge beams.

Now, we can compute from Eq.~(\ref{eq:100.2}) the total electromagnetic linear momentum associated with the charge $Q_1$ and the solenoid of current by summing over $M l$ current loops of the solenoid, and using ${B}_{\rm z}= {4\pi I M}/{c}$ inside the solenoid (Eq.~(\ref{eq:84b.4b1a1})), to get
\begin{eqnarray}
\label{eq:84b.6}
{\bf p}_{\rm e1}=\sum_{1}^{Ml} \frac {({\bf E}_1 \times {\bf a})\,I}{c^2}=\int\frac {({\bf E}_1 \times {\bf B})}{4\pi c} {\rm d}\tau\,,
\end{eqnarray}
which is the volume integral of the EM momentum density, $({\bf E} \times {\bf B})/4\pi c$ over the solenoid.
For the charge $Q_2$, with ${\bf E}_2=-{\bf E}_1$, Eq.~(\ref{eq:84b.6}) implies ${\bf p}_{\rm e2}=-{\bf p}_{\rm e1}$.

The seat of the field momentum (Eq.~(\ref{eq:84b.6})) might appear to be within the solenoid, however, one has to take the holistic view that the EM momentum actually lies in the {\em composite system} of the charge plus solenoid, as it has been emphasized  \cite{Va12} that the composite system is represented by one state. Accordingly the system acts as a whole, giving rise to the momentum, that gets reflected in the Aharonov-Bohm quantum interference experiment. This non-localized interaction seems to be the explanation of this intriguing phenomenon from a classical physics perspective. The vector potential ${\bf A}$ may be still considered in this case as a convenient, intermediary mathematical step with $Q{\bf A}/c$ presenting a fa\c{c}ade of the mutual electric interaction of the conducting current carriers in the solenoid and an external charge $Q$, which gives rise to an  EM momentum in the system.

\section*{Declarations}
The author has no conflicts of interest/competing interests to declare that are relevant to the content of this article. No funds, grants, or other support of any kind was received from anywhere for this research.