\subsection{Related Work}
Researchers have been discussing the potential of designing EV charging scheduling algorithms to alleviate grid congestion~\cite{xiong2018vehicle}, to better integrate increasing renewables~\cite{gao2014integrated}, and to accommodate economical objectives as well as operating constraints~\cite{lee2018large}. Real-time algorithm considering arrival and departure statistics are discussed in \cite{tucker2022real}.
Bi-directional charging strategies and vehicle-to-grid (V2G) technologies have been also discussed in \cite{xu2020greenhouse}.

There are some existing work on approaching load shifting problem via solving optimization problem explicitly \cite{chen2021smoothed,lee2019acn, quiros2015statistical} or via time-of-use pricing. Even though electricity prices follow some characteristic of the the load curve (lower cost during the noon and higher cost during the evening), cost optimization using TOU pricing is an indirect approach for optimizing carbon emissions. Based on a recent EV industry report~\cite{sense_report}, there is unutilized potential of carbon reduction via enabling flexible charging windows

Carbon-aware EV charging scheduling is still novel in the way that it combines grid carbon intensity data, EV statistics and optimization framework. \cite{will2022consumer} conducted a survey on EV owners' preferences regarding carbon-neutral charging services. In \cite{ensslen2017empirical}, empirical studies validate that different geographical regions and charging schedules can lead to different carbon emissions.
\cite{colmenar2019electric} proposed a heuristic method to match EV charging with renewable energy resources on a macro scale without considering individual characteristics of charging sessions and charging station's physical constraint. In \cite{chen2019electric}, a demand response scheme with simplified charging models is proposed for the overall US market.

Our technique bears the closest resemblance to the industrial advances mentioned in \cite{parker2020,Sense2021}, and is partly inspired by the carbon-aware cloud computing considering flexible data center loads~\cite{radovanovic2021carbon}. The key distinction lies in how we formulate and solve the optimization problem by taking both carbon intensity forecasts and engineering constraints into account. We also develop carbon intensity forecasting method to accommodate real-time scheduling.

\section{Introduction}
\input{introduction}

\section{Carbon-Aware Charging Framework}
\input{framework}

\section{Method}
\label{sec:model}
\input{model}

\section{Real-Time Scheduling}
\input{online}

\section{Simulation Results}
\label{sec:simulation}
\input{simulation}

\section{Conclusion and Future Work}
In this paper, we explore how EV charging scheduling will positively affect the carbon emissions of EV charging. More interestingly, in this paper, we use the grid-level average carbon intensity as an indicator for optimizing charging facilities' emission profiles. In the future work, we will investigate the choices of carbon and economical objectives, and explore the planning and operation mechanisms considering multiple EV charging stations. We will also explore how pricing mechanisms can affect the session-specific carbon emissions.

\bibliographystyle{IEEEtran}

\subsection{Carbon Intensity}
We denote the time-varying carbon emission (carbon intensity) as $C(t)$ with a unit of kgCO2/kWh or mTCO2(Metric ton CO2)/MWh. Carbon intensity data is a time-series data that reflects the grid-level carbon emission per energy unit. Carbon intensity depends on the power grid's operational plan and carbon intensity of each energy resource.
In this paper, we specifically examined the carbon intensity for the California grid with data provided by CAISO (California Independent System Operator)\footnote{\url{http://www.caiso.com/todaysoutlook/pages/emissions.html}}. To calculate the carbon intensity, we consulted the carbon emission of the grid resources as well as the supply trend of the grid. Figure 2 describes the feature of the two data sets. In Fig \ref{fig:carbon_caiso} (left) there are two main contributors to the state's grid emission: imports and natural gas.
\begin{equation}
    \label{equ:carbon intensity}
    C(t) = \sum_{i=1}^{N} C_i(t) \cdot \frac{P_i(t)}{P_{total}(t)}, \;  \forall i \in \mbox{grid mix};
\end{equation}
where $C_i(t)$ is a time-series data of the carbon intensity of each energy source, $P_i(t)$ is the power supplied by the energy source and $P_{total}(t)$ is the total power generated on the specific day at time $t$.
From the formula and the data, we can expect the carbon intensity to be high during the night and the early morning when solar energy is less available. On the other hand, the carbon intensity of the grid tends to be low during the day as solar energy and other renewables become a major contributor to the state's power supply.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.95\columnwidth]{figures/caiso_carbon_data.png}
    \caption{(Left) The carbon intensity of each resources. We exclude biogas, biomass, geothermal, and coal since it has relatively low carbon contribution in California. (Right) Daily average power supply per resource. We exclude coal and hydro since it is less significant. }
    \label{fig:carbon_caiso}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.75\columnwidth]{figures/seasonal_daily_carbon_intensity_average.png}
    \caption{Average carbon intensity during different seasons.}
\end{figure}
Carbon intensity is not only grid-mix dependent and time sensitive, but is also seasonally sensitive. Fig 3 shows an example of carbon intensity calculated using the above formula throughout the year. We notice that carbon intensity is lower especially during the spring. On top of that, we also notice a shift in time for the carbon intensity curve in autumn as well as winter. We can attribute this to the daylight saving hours policy in California.

\subsection{Charging Station Operation Model}
As mentioned in Section II, our charging framework takes into account of the user constraints as well as the infrastructure constraints. These constraints are discussed below.
We use the discrete time model, where the time step is indexed in $\mathcal{T} := \{1,2,3,\ldots,T\}$.  $\mathcal{V}$ is the set of all EVs. For each vehicle $i \in \mathcal{V}$, we denote the charging demand information is $(x_{i,arrival}, x_{i,depart}, t_{i,arrival}, t_{i,depart})$, where $x_{i,arrival}$ is the state-of-charge SoC of EV $i$ upon its arrival of the charging station, $x_{i,depart}$ is the desired state-of-charge SOC of EV $i$ upon its departure. and $t_{i,arrival}, t_{i,depart}$ is the time index when the vehicle $i$ arrive in / depart from the charging station, $t_{i,arrival}, t_{i,depart} \in \mathcal{T}$.The above four factors are all externally determined by the user. Arrival time $t_{i,arrival}$ and initial energy level $x_{i,arrival}$ are determined when user plugged in the EV. Departure time $t_{i,depart}$ and charge $x_{i,depart}$ are input by the user when they initiate the charging session.

The model aims at deciding the charging power $u_i(t)$ for vehicle $i$ at time $t$, and the state-of-charge (SoC) $x_i(t)$ for EV $i$ at time $t$ given the specific objective.
We are interested in solving the following optimization problem:
\begin{subequations}
    \label{equ:carbon_aware_EV}
     \begin{align}
    \min_{\mathbf{u}} \quad & \sum_{t \in \mathcal{T}} \sum_{i\in \mathcal{V}}  C(t)u_i(t) \Delta T + \lambda \sum_{i\in \mathcal{V}}|x_i(T)-x_{i,depart}| \label{equ:obj}\\
    s.t. \quad &  x_i(0)=x_{i,arrival}, \label{eq:initial_soc}\\
    &  x_{i}(t)=x_i(t-1)+ \frac{\delta u_i(t)}{E_i}, \; t\geq 1; \label{subequ:step}\\
        & 0\leq x_i(t) \leq \bar{x}_{i}, \; i\in \mathcal{V}; \label{subequ:soc}\\
    & 0\leq u_i(t) \leq \bar{u}_i, \; i\in \mathcal{V}; \label{subequ:action}\\
    & \sum_{i \in \mathcal{V}}  u_i(t) \leq \bar{P}, \quad t\in \mathcal{T};\label{subequ:total_u}\\
    & u_i(t)=0, \quad t \notin [t_{i,arrival},t_{i, depart}), \; i \in \mathcal{V}. \label{subequ:time}
     \end{align}
\end{subequations}
The above optimal charging scheduling problem minimizes the total carbon emissions while satisfying charging sessions' time and energy constraints, defined by the objective function~\eqref{equ:obj}. We firstly consider the case when $C(t)$ and EV charging demand information for next $T$ steps are available and communicated to the charging station operator. We denote $\Delta T$ to be the discretization interval, and in our case, $\Delta T$ is 5 minutes aligned with the carbon intensity measurement interval published by CAISO. $\lambda$ is a balanced factor to trade off between carbon reduction and energy delivery. The operational constraints considered in this framework are described as follows:
\newline \textbf{EV charging constraints}: Constraint \eqref{eq:initial_soc} defines the SoC level of EV upon arrival; and  \eqref{subequ:step} depicts the dynamics between charging power and battery SoC, where $E_i$ is the battery energy capacity (kWh). $\delta$ is an EV-specific charging parameter, which is known before the operation and is computed based on the charging efficiency and discretization, e.g., if charging efficiency is $90\%$, and discretization interval is 15 minutes, $\delta = 0.9 \times 15/60 = 0.225$. \eqref{subequ:soc} and \eqref{subequ:action} constrain the admissible charging power for single EV, where the battery maximum SoC is $\bar{x}_i$ and maximum charging power is $\bar{u}_i$.
\newline \textbf{Charging station constraints}: \eqref{subequ:total_u} limits the total available charging power $\bar{P}$ at the charging station. We note that in this paper, we assume there is a known $\bar{P}$ decided by the charging station's transformer.
\newline \textbf{Temporal logic of charging sessions}: \eqref{subequ:time} describes the temporal logic that the charging session has to satisfy, as no charging power shall be delivered whenever EV $i$ has not arrived ($t < t_{i,arrival}$) or has departed $(t \geq t_{i,depart}).$
\newline\textbf{Other Factors to generalize}: The optimization problem that we introduced in (2a) is a general model that can be further generalized. For example, battery health is another important factor to be considered in EV charging. A more complex version of this optimization problem could potentially include penalty terms in the objective to discourage sudden change in charging power, which prevents battery health from deteriorating.

\subsection{Carbon Intensity Day-Ahead Forecasting}
We use the 2021 CAISO data (shown in Fig~\ref{fig:carbon_caiso}) and CAISO day-ahead load forecast data $Load(t)$ to fit and evaluate the carbon predictor model. The grid's day-ahead carbon intensity $\hat{C}(t)$ is closely related to the future electricity demand, calendar features, and past carbon intensity patterns. We randomly choose $80\%$ data to fit the linear regression model and rest data for testing the model performance. The resulting forecasting model is as follows:
\begin{equation}
\begin{split}
\label{equ:carbon_predict}
    \hat{C}(t)&=\beta_0+\beta_1\times Minute + \beta_2 \times Hour + \beta_3 \times Day \\&+ \beta_4 \times Month + \beta_5 \times Load(t)+ \beta_6 \times MA(C(t))_{24}\\&+\beta_7 \times MA(C(t))_{12}+\beta_8 \times MA(C(t))_{1};
    \end{split}
\end{equation}
where $MA(\cdot)_{a}$ denotes the moving average of past $a$ hours. Fig~\ref{fig:carbon_forecast} illustrates the ground truth carbon intensity and forecasted carbon intensity using the fitted model \eqref{equ:carbon_predict}. We note that other statistical and machine learning models can also be used for carbon intensity prediction. We use the linear regression model since it is simple and provides accurate prediction result. For the CAISO case study, the overall testing mean squared absolute error is $0.0085$.
We also demonstrate such forecasts can be readily integrated into the proposed charging scheduler without affecting the algorithm performance. Power grid CO2 emission intensity forecast is also investigated in \cite{leerbeck2020short} by considering additional features and an autoregressive model.
\begin{figure}[t]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/Carbon_Forecasts.png}
    \caption{The ground truth (blue) and predicted carbon intensity (red) for two-day sample in January and July respectively.}
    \label{fig:carbon_forecast}
\end{figure}
We calculate and forecast the average carbon intensity rather than marginal carbon intensity, as in this study, we focus on the operating strategies for single EV charging station. Such single-station operation scheme shall pose minimal effects to transmission grid carbon profiles.  We refer to \cite{carbon2022online} for a more complete review and open-source data on electricity grid carbon intensities.

\subsection{Online Algorithm}
With the predicted carbon intensity $\hat{C}(t)$, we are ready to implement proposed carbon-aware charging scheduler in an online fashion. Denote the optimization problem \eqref{equ:carbon_aware_EV} as $\textbf{SCH}$ with input information coming from the charging infrastructure and EVs' peak power constraints $\bar{P}$ and $\bar{u}$, predicted carbon intensity, the active charging vehicles $\mathcal{V}_k$ along with their states $x_{\mathcal{V}_k}$ where $k \in \mathcal{K}, \mathcal{K} := \{1,2,3,\ldots,K\}$, the algorithm is summarized in Algorithm \ref{alg:algorithm}.
\revise{}

\begin{algorithm}
\caption{Online Carbon Aware EV Charging}
\label{alg:algorithm}
\begin{algorithmic}
\Require Charging infrastructure and each EV's peak charging power limit $\bar{P}$ and $\bar{u}$, operation time length $K$.
\Ensure Charging power $(u_i(1), u_i(2), \ldots, u_i(K), i \in \mathcal{V})$.
\For{$k = 1, \ldots, K$}
\State  $\mathcal{C}_k: = \{\hat{C}(k), \ldots, \hat{C}(k+T-1)\}$
\State  $\mathcal{V}_k :=  \{i \in \mathcal{V} | x_{i,arrival} \leq k \textbf{ AND }  x_{i,depart} > k\}$
\State  $x_{\mathcal{V}_k} :=  \{x_i(k) | i \in \mathcal{V}_k \}$
\State \textbf{compute}
\State \quad
 $(u^*_i(1), \ldots, u^*_i(T), i \in \mathcal{V}_k) = \textbf{SCH}(\mathcal{C}_k, \mathcal{V}_k, x_{\mathcal{V}_k})$
\State \textbf{update} $u_i(k) = u^*_i(1)\,, \forall i \in \mathcal{V}_k$ and $u_i(k) = 0, i \not\in \mathcal{V}_k$.
\State \textbf{update} $x_i(k+1) = x_i(k) + \frac{\delta u_i(k)}{E_i}, \forall i \in \mathcal{V}$.
\EndFor
\end{algorithmic}
\end{algorithm}

At each time step $k = 1, 2, \ldots, K$, we take the forecasted carbon intensity for the next $T$ steps, $\mathcal{C}_k: = \{\hat{C}(k), \ldots, \hat{C}(k+T-1)\}$, the current available EVs $\mathcal{V}_k$ and their SoC levels $x_{\mathcal{V}_k}$, and solves the optimization~\eqref{equ:carbon_aware_EV}. With the optimized charging schedule $(u^*_i(1), \ldots, u^*_i(T), i \in V_k)$, we assign  $u^*_i(1)$ to the charging power of EV $i$ at time step $k$ to $u_i(k) = u^*_i(1)$ for the active vehicles and update state of charge accordingly for every vehicle $x_i(k+1) = x_i(k) + \frac{\delta u_i(k)}{E_i}, \forall i \in \mathcal{V}$.

\subsubsection{Seasonal Shift}
We first present the seasonal carbon emission shift using our method and the baseline EDF method. We select days from March 1 to May 31 as spring, June 1 to August 31 as summer, September 1 to November 30 as autumn and December 1 to February 28 as winter. Fig~\ref{fig:season} shows the comparison of carbon emission and energy delivery (\%) among four seasons with the balanced factor $\lambda = 0.35$.
From the figure, we can see that carbon intensity is lower at noon than other time. Our method takes advantage of the low carbon intensity pattern and shifts the flexible charging demand to the time when carbon intensity is low, such that the carbon emission is reduced efficiently.

\subsubsection{The impact of the balanced factor}\label{sec:impact}
We use the first 20 days in the dataset to show how balanced factor $\lambda$ in ~\eqref{equ:obj} influences the carbon emission and energy delivery quality. Fig~\ref{fig:diff_factor} shows the change of total carbon emission (red curve) and energy loss (blue curve) using different balanced factors. We can see that with the balanced factor decreasing, the carbon emission will decrease while the energy loss will increase. We also compare the performance with the baseline EDF method. Using a balanced factor $\lambda \geq 0.4$, we can achieve the {same} energy delivery while reducing 5.6\% CO2 for 20 days.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/different_factor_3.pdf}
    \caption{The tradeoff between the total carbon emission and percentages of energy loss among different balanced factor using our method and baseline (EDF). We use the first 20 days to simulate the result.}
    \label{fig:diff_factor}
\end{figure}

\subsubsection{Comparison with baselines}
Here we use the entire dataset in 2021. We compare our method with three baseline methods: ES, EDF and TOU.
We test our method with a fix balanced factor $\lambda = 0.3$ (denoted as \texttt{Offline(0.3)}) and $\lambda = 0.4$ (denoted as \texttt{Offline(0.4)}), respectively. We also test the carbon-aware charging algorithm with a flexible factor scheme (denoted as \texttt{Offline(+)}), where we choose $\lambda = 0.3$ for days with vehicle size $n \leq 20$, $\lambda = 0.25$ for days with $n \in [21,30]$ and $\lambda = 0.35$ for other cases.
Secion~\ref{sec:impact} shows that larger balanced factor is more likely to deliver the required energy while leading to less carbon reduction. Since group one (less than 20 charging sessions per day) has relatively low carbon emission, and group three (more than 30 sessions charging sessions per day) has less flexibility in shifting charging demand, we choose larger balanced factor for group one and three to prioritize the energy delivery.
The results are shown in Figure~\ref{fig:compare_adapted}. With $\lambda = 0.4$, the offline optimization can reduce $3.97\%$ carbon emission while slightly improves the energy delivery quality by $0.03\%$ compared to EDF and ES, and it can reduce $1.50\%$ carbon emission and maintain the same energy delivery quality compared to TOU. This is because EDF and ES will schedule charging instantly once EVs arrive, while TOU and our approach adjust the charging power by solving the current optimization instance. TOU rate and carbon intensity is higher after 4pm. To minimize the energy cost, ToU-aware policy is likely to encourage chargers to finish charging before 4pm, and in this way it can also partially help reduce the carbon emission.
With $\lambda = 0.3$, the carbon aware optimization method \texttt{Offline(0.3)} can save $24.7\%$ and $22.71\%$ carbon emission compared to EDF/ES and TOU respectively,
while also lead to $18.20\%$ in reduction in energy delivery quality.
With flexible balanced factors, \texttt{Offline(+)} saves $3.17\%$ carbon and improves $5.55\%$ in delivery quality compared to \texttt{Offline(0.3)}. An interesting future direction is to develop a more principled method of adaptively choosing the balanced factor according to daily charging patterns, that can both maintain the high energy delivery quality and effectively reduce carbon emission.

\begin{table}[htbp]
    \centering
    \begin{tabular}{c c}
    \hline
      \textbf{Strategy}   &  \textbf{Carbon emission (kg)} \\
      \hline
      Equal Sharing   &  2.913 \\
      Earliest Deadline First & 2.913 \\
      ToU price aware &  2.840 \\
      Offline (0.4) & 2.798 \\
      Offline (+)   & 2.127 \\
      Offline (0.3) & 2.195 \\
      \hline
    \end{tabular}
    \caption{Average carbon emission for each charging session using different approaches.}
    \label{tab:compare}
\end{table}
Table~\ref{tab:compare} shows the average carbon emission for each charging session. Compared to EDF, offline with $\lambda=0.4$ reduces 0.115 kg CO2, which is equivalent to travelling 1.15 miles\footnote{\href{https://www.tesla.com/ns_videos/2021-tesla-impact-report.pdf}{Tesla’s 2021 Impact Report}}; offline with adapted balanced factors saves 0.786 kg CO2, corresponding to the travelling distance of 7.86 miles.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/comparea.pdf}
    \caption{The comparison of average energy deliver quality per session (EDQ-Session) and carbon emissions: offline (0.3), offline (+) (smarter balanced factor strategy), offline (0.4), ToU aware (TOU), earliest deadline first (EDF), equal sharing (ES). }
    \label{fig:compare_adapted}
\end{figure}

\subsection{Online EV Charging Scheduling Results}
We next demonstrate the performance of the online scheduling algorithm. We use the carbon intensity prediction and online EV information to generate the charging schedule, that is, at each time step $t$, the model only knows the information of arrived vehicles, and use predicted carbon intensity data for the optimization.

Fig~\ref{fig:online} shows that our online algorithm does not hurt the performance much compared to the offline optima. With $\lambda = 0.3$, the online algorithm delivers $0.21\%$ less energy and generates $0.08\%$ less carbon emission compared to the offline algorithm.
With $\lambda = 0.4$, the online algorithm delivers $0.002\%$ less energy and generates $0.15\%$ more carbon emission compared to the offline algorithm.
Compared to EDF and ES, the real-time scheduling saves $3.81\%$ of carbon emission and can keep delivering the same energy; compared to TOU, it saves $1.50\%$ carbon and maintain the same energy delivery quality.
With balanced factor $\lambda = 0.3$, the online algorithm can further reduce $27.01\%$ of carbon emission at the cost of delivering $14.21\%$ less energy than EDF and ES; and reduces $25.13\%$ carbon at the cost of delivering $12.61\%$ less energy compared to TOU. All results indicate the online algorithm is highly efficient compared to the offline solutions.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/online.pdf}
    \caption{The comparison of energy delivery quality for station (EDQ-Station) and carbon emission: offline (carbon aware offline), online (real time carbon aware), ToU aware (TOU), EDF, and ES. Proposed online algorithm achieves similar performances compared to offline benchmarks.}
    \label{fig:online}
\end{figure}


\section{ Introduction }
\end{center}}
We denote by $\tn$   the set of trigonometric polynomials of degree at  most $N$
\[
s(x)=\sum_{k=-N}^N c_ke^{ikx}
\]
with a positive integer $N$ and $c_k\in\Co$, $k=-N,\dots, N$. Suppose that   $x_1<x_2<\dots<x_r$, $r>1$,   is any sequence in $\T=(-\pi, \pi]$ such that
 \begin{equation}
\min\Bigl\{x_{j+1}-x_j, j=1,\dots, r-1; \ 2\pi-(x_r-x_1)\Bigr\}=\delta>0.
\end{equation}
The  usual large sieve inequality states that
  \begin{equation}
\sum_{j=1}^r|s(x_j)|^2\le \Bigl(\frac{N}{2\pi}+\frac1{\delta}\Bigr)\int_{\T}|s(x)|^2\,dx.
 \end{equation}
See  Selberg \cite[p. 221]{8}, but note the different notation.   Gallagher \cite{4} has given a simple derivation of the large sieve inequalities. It turned out that the method of \cite{4}  can be applied to  $L^p(\T)$-norm.   More precisely, in \cite[p. 96]{3} was proved that  if $s\in\tn$ and $2\le p<\infty$, then
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \Delta_pp^{1/2} \max\biggl(N, \, \frac{4\pi}{\delta}\biggr)\Bigl(\sum_{k=-N}^N|c_k|^{q}\Bigr)^{p/q},
\end{equation}
where $\Delta_p$ is an absolute constant and $1/p+1/q=1$. Next, in \cite[p. 533]{6} the inequality (3) was extended for all $0<p<\infty$ as follows. Let $\Psi$ be  a convex, nonnegative, and nondecreasing function  in $[0,\infty)$. Then for any $s\in T_N$,
  \begin{equation}
\sum_{j=1}^r\Psi(|s(x_j)|^p)\le \Bigl(\frac{N}{\pi}+\frac1{\delta}\Bigr)\int_{\T}\Psi\Bigl(|s(x)|^p(p+1)e/2\Bigr)\,dx,
 \end{equation}
whenever $-\pi<x_1< x_2<\dots <x_r\le \pi$ and $\delta$ is defined by (1). If $\Psi(t)=t$, $0\le t<\infty$, then (4) implies
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \Bigl(\frac{N}{\pi}+\frac1{\delta}\Bigr)\frac{(p+1)e}{2}\|s\|^p_{L^p(\T)}.
 \end{equation}
 In \cite[p. 164]{5}  the estimate (4) was extended to the case of generalized trigonometric polynomials. In particular, for usual trigonometric polynomials , the inequality (5) was improved as follows:
\begin{equation}
 \sum_{j=1}^r|s(x_j)|^p\le \Bigl(\frac{N+1}{2\pi}+\frac1{\delta}\Bigr)\frac{(p+1)e}{2}\|s\|^p_{L^p(\T)}.
 \end{equation}

Note that  inequalities (2)-(6) are also called   forward Marcinkiewicz-Zygmund inequalities (see, e.g., \cite{2}).

In this paper,   we will develop  an approach to inequalities of the type (3) and (5)-(6)  that not use the usual $L^2$-technique. Our approach is based on the  spectral theory of convolution operators on  $\tn$.

 The main result is given in the following theorem. Note that in the sequel, $[x]$ denotes the integer part of a positive number $x$. Also we use the notation $\Gamma(\cdot)$ for the standard gamma function.
\begin{theorem}.
Let $\{x_j\}_{j=1}^r$ be a sequence in $\T$  that satisfies (1). If  $s\in\tn$ and $1\le p<\infty$, then
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \frac{pN\sigma(\delta; N)}{2\sqrt{\pi}}\cdot \frac{ \Gamma\Bigl(p/2\Bigr)}{\Gamma\Bigl(p/2+1/2\Bigr)}\|s\|^p_{L^p(\T)}
\end{equation}
with
\begin{gather}
\sigma(\delta; N)=\begin{array}\{{rl}.
\frac{\pi}{N\delta} ,\ \ \  \ if  &  \frac{\pi}{N\delta}\in \Z, \\
1+\biggl[\frac{\pi}{N\delta}\biggr], & \mbox{otherwise.}
\end{array}
\end{gather}\
 \end{theorem}
If $p$ is a positive integer, then the quantity
\[
 \frac{ \Gamma\Bigl(p/2\Bigr)}{\Gamma\Bigl(p/2+1/2\Bigr)}
 \]
 can be calculated directly by using the relations between $ \Gamma\Bigl(p/2\Bigr)$ and $\Gamma\Bigl(p/2+1/2\Bigr)$.
 \begin{corollary}.
 Let $ l$ be a positive integer. Then under the conditions of Theorem 1 it follows that:

 (i) \ if $p=2l$, then
 \begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \frac{pN\sigma(\delta; N)}{2\sqrt{\pi}}\cdot \frac{ \Gamma\Bigl(p/2\Bigr)}{\Gamma\Bigl(p/2+1/2\Bigr)}\|s\|^p_{L^p(\T)}=\frac{pN\sigma(\delta; N)\cdot 2^{l-1}(l-1)!}{\pi\Bigl(1\cdot 3\cdot 5\cdot \cdots \cdot (2l-1)\Bigr)}\|s\|^p_{L^p(\T)};
\end{equation}
 (ii) \ \ if $p=2l+1$, then
 \begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \frac{pN\sigma(\delta; N)}{2\sqrt{\pi}}\cdot \frac{ \Gamma\Bigl(p/2\Bigr)}{\Gamma\Bigl(p/2+1/2\Bigr)}\|s\|^p_{L^p(\T)}=\frac{pN\sigma(\delta; N)\cdot \Bigl(1\cdot 3\cdot 5\cdot \cdots \cdot (2l-1)\Bigr)}{2^{l+1}\cdot l!}\|s\|^p_{L^p(\T)}.
\end{equation}
 \end{corollary}
\begin{corollary}.
Assume that $s\in\tn$,  $\|s\|_{L^p(\T)}=1$, $p\ge 1$,  and  $\{x_j\}_{j=1}^r$  satisfies (1).

(i) \  If $\frac{\pi}{N\delta}\in \Z$, then
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p<\frac{ p+1}{\delta}.
\end{equation}
(ii) If $\frac{\pi}{N\delta}\not\in \Z$, then
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p< (p+1)\biggl(\frac{N}{\pi}+\frac1{\delta}\biggr).
\end{equation}
\end{corollary}

\section{  Proofs}
Let $M(\T)$ be the Banach algebra of finite complex-valued regular Borel measures on $\T$. The norm in $M(\T)$ is given by the total variation $\|\mu\|$ of $\mu\in M(\T)$. Therefore, the usual Banach space $L^1(\T)$ can be identified with the closed ideal in $M(\T)$ of all measures which are absolutely continuous with respect to the Lebesgue measure $dt$ on $\T$.

Given $\mu\in M(\T)$ and $f\in L^1(\T)$, we define the Fourier transform of $\mu$ and $f$  by
\[
\widehat{\mu}(x)=\int_{\T}e^{-ixt}\,d\mu(t) \quad{\text{ and}}\quad \widehat{f}(x)=\int_{\T}e^{-ixt}f(t)\,dt,
\]
respectively.
For each $u\in L^r(\T)$, $1\le r\le\infty$, on $T_N$ is well defined  the convolution operator
\begin{gather}
A_u(s)(x)=s\ast u(x)=\int_{\T}s(x-t) u(t)\,dt= \int_{\T}\Bigl(\sum_{k=-N}^N c_ke^{ik(x-t)}\Bigr)u(t)\,dt\nonumber\\
=\sum_{k=-N}^N c_k\widehat{u}(k)e^{ikx},
\end{gather}
for each $x\in \T$, where $s\in\tn$. Note that in this definition and also below we assume  that $s$ is a periodic function on the real line with the period equal to  $2\pi$. Below,  the notation $T^p_N$, $1\le p\le\infty$, means that $T_N$ is equipped with the usual $L^p(\T)$ norm.

Suppose that $S_1$  and $S_2$ are two  measurable subsets of $\T$, $\mu_1$ and $\mu_2$ are two non-negative finite measures on $\T$, and $ F:\T^2\to\R$ is a measurable function.  For $1\le p<\infty$, Minkowski's integral inequality \cite[p. 37]{9} states that
\begin{gather}
\biggl[\int_{S_2}\biggl|\int_{S_1}F(x,y)\,d\mu_1(x)\biggr|^p\,d\mu_2(y)\biggr]^{1/p}\nonumber\\
\le\int_{S_1}\biggl(\int_{S_2}|F(x,y)|^p\,d\mu_2(y)\biggr)^{1/p}\,d\mu_1(x).
\end{gather}

\begin{proposition}.
Let $u\in L^q(\T)$, $1\le q\le \infty$, $\|u\|_{L^q(\T)}\neq 0$. Assume that $u$ is continuous,  non-negative and even on $\T$. If
\begin{equation}
{\text{supp}} \ u\subset\Bigl[-\frac{\pi}{2N}\,, \ \frac{\pi}{2N}\Bigr],
\end{equation}
then there exists a trigonometric sum
\begin{equation}
p_u(x)=\sum_{m=-N+1}^{N}\tau_m e^{-i\pi m x/N}
\end{equation}
such that
\begin{equation}
(-1)^m\tau_m> 0,
\end{equation}
$m=-N+1,\dots ,N$, and
\begin{equation}
p_u(n)=\frac1{\widehat{u}(n)}
\end{equation}
for all $n=-N,\dots,N$.
\end{proposition}
PROOF.\quad We start by examining in more details the Fourier transform of $u$. Under the assumptions on $u$, we see that
\begin{equation}
\widehat{u}(x)=\int_{\T}u(t)e^{-ixt}\,dt=2\int_{0}^{\pi/2N}u(t)\cos xt\,dt>0
\end{equation}
for all $x\in [-N, N]$. From this it follows that
\begin{equation}
(\widehat{u})'(x)=-2\int_{0}^{\pi/2N}t u(t)\sin xt\,dt<0
\end{equation}
and
\begin{equation}
(\widehat{u})''(x)=-2\int_{0}^{\pi/2N}t^2u(t)\cos xt\,dt<0
\end{equation}
for all $x\in [0, N]$.

Let $v=1/\widehat{u}$. By (19),  the function $v$  is  well defined and positive on $[-N, N]$. Moreover, we conclude from (19)-(21) that $v$ is an even function and
\begin{equation}
v'(x)>0 \quad {\text{and}}\quad v''(x)>0
\end{equation}
for all $x\in [0, N]$. Therefore,   $v$ is increasing and   convex on $[0,N]$, in particular $v$ is of bounded variation  on $[0, N]$. In particular, this means that $v$  is a function of bounded variation on $\T$. The following is well known: If $f$ is an  $2\pi$-periodic
continuously differentiable even function on $\T$ such that $f$  is of bounded variation on $\T$, then the Fourier series of $f$  converges absolutely (see, eg. \cite[p. 241]{10}).  Thus,
\begin{equation}
v(x)=\sum_{k\in\Z}a_ke^{ikx\pi/N}=\sum_{k\in\Z}a_k\cos\Bigl(\frac{\pi}{N}kx\Bigr)
\end{equation}
with
\begin{equation}
\sum_{k\in\Z}|a_k|<\infty,
\end{equation}
where
\begin{equation}
a_k=\int_{-N}^{N}v(t) e^{-ikt \pi/N}\,dt=2\int_{0}^{N}v(t)\cos\Bigl(\frac{\pi}{N}kt\Bigr)\,dt.
\end{equation}
We claim that
\begin{equation}
(-1)^ka_k>0
\end{equation}
for all $k\in\Z$. Combining (19) with (25), we see that  $a_0>0$. Let $k\ge 1$. Then using
integration by parts, we conclude from (15) that
\begin{equation}
a_k=-\frac{2N}{\pi k}\int_{0}^{N}v'(t)\sin\Bigl(\frac{\pi}{N}kt\Bigr)\,dt=-\frac{2N}{\pi k}\ \sum_{j=0}^{k-1}I_j,
\end{equation}
where
\begin{equation}
I_j=\int_{E_j}v'(t)\sin\Bigl(\frac{\pi}{N}kt\Bigr)\,dt
\end{equation}
and $E_j=[Nj/k, N(j+1)/k]$. Note that the length of $E_j$  is   exactly  half length of period for $\sin\Bigl(\pi kx/N\Bigr)$. Combining this with (22), we see that
\begin{equation}
(-1)^j I_j>0 \quad {\text{and}}\quad |I_j|<|I_{j+1}|
\end{equation}
for all $j=0,\dots,k-1$.  From this,  it is easily seen  that
\[
(-1)^{k-1}\sum_{j=0}^{k-1} I_j>0.
\]
In light of (27) this  proving the claim (26).

For $m=-N+1,\dots, N$, let $\tau_m$ be defined by
\begin{equation}
\tau_m=\sum_{j\in\Z}a_{m+2jN}.
\end{equation}
From (24), we see that the series in (30) converges absolutely. Next, from (26) it follows that,   for each $m$, all  terms of the sequence $\{a_{m+2jN}\}_{j\in\Z}$ have the same sign. In particular, (26) shows that $(-1)^ma_{m+2jN}>0$. Therefore, the trigonometric sum  (16) is well defined and satisfies (17).

Finally, for any $n\in\{-N,\dots, N\}$, combining (23) with (30), we get
\begin{gather}
\frac1{\widehat{u}(n)}=v(n)=\sum_{k\in\Z}a_ke^{ik\pi n/N}=\sum_{m=-N+1}^N\biggl(\sum_{j\in\Z}a_{m+2jN}e^{i(m+2jN)\pi n/N}\Biggr)
 \nonumber\\
=\sum_{m=-N+1}^N\biggl(\sum_{j\in\Z}a_{m+2jN}e^{im\pi n/N}\biggr)=\sum_{m=-N+1}^N\biggl(\sum_{j\in\Z}a_{m+2jN}\biggr)e^{im\pi n/N}\biggr)\nonumber\\=\sum_{m=-N+1}^N\tau_me^{imn\pi/N}=p_u(-n).\nonumber
\end{gather}
As $u$ and $\hat{u}$ are even functions we have $p_u(-n)=p(n)$ for all $n$.
Proposition 4 is proved.

We will denote by $\delta_a$ the usual Dirac measure supported on $a\in\T$.

\begin{proposition}.
Let  positive numbers $p$ and $q$ satisfy $1/p+1/q=1$. Under the conditions of Proposition 4 on  $u\in L^q(\T)$ it follows that the operator (13) possesses on $\tnp$ a bounded inverse of the type
\begin{equation}
A^{-1}_u(s)(x)=\int_{\T}s(x-t)\,d\mu(t),
\end{equation}
where
\begin{equation}
\mu=\sum_{m=-N+1}^{N}\tau_m \delta_{\pi m/N}
\end{equation}
and  $\{\tau_m\}_{-N+1}^{N}$ are defined by (30), using (23) for $v=1/\hat{u}$.  Furthermore,
\begin{equation}
\|A^{-1}_u\|_{\tnp}=\|\mu\|_{M(\T)}=\sum_{m=-N+1}^N|\tau_m|=\frac1{\widehat{u}(N)}.
\end{equation}
\end{proposition}
PROOF.\quad Let $s(x)=\sum_{k=-N}^Nc_ke^{ikx}\in \tn$. From (16) it follows that
\begin{equation}
s\ast\mu(x)=\sum _{k=-N}^Nc_k\biggl(\int_{\T}e^{ik(x-t)}\,d\mu(t)\biggr)=\sum _{k=-N}^Nc_k\widehat{\mu}(k)e^{ikx},
\end{equation}
where
\[
\widehat{\mu}(k)=\int_{\T}e^{-ikt}\,d\mu(t)=\sum_{m=-N+1}^N\tau_me^{-i\pi km/N}=p_u(k),
\]
for all $k=-N,\dots, N$. Now, taking into account (18), we conclude from from (13) that (31) defines  the inverse of the operator $A_u$.

Now,  according to  (34), we see that  the set  $\{\widehat{\mu}(k)=p_u(k): \ k=-N,\dots,N\}$ coincides with the spectrum of $A^{-1}_u$. Therefore, if $| A^{-1}_u|_{\tnp}$ denotes the spectral radius of $A^{-1}_u$, then
\[
| A^{-1}_u|_{\tnp}= \max\{|\widehat{\mu}(k)|=|p_u(k)|: \ k=-N,\dots,N\}.
\]
Let us recall that  $|H|_X\le \|H\|_X$, i.e. the spectral radius is not greater than the operator norm,  for any bounded linear operator $H$ on a normed space  $X$.
Combining this with (18), we get
\begin{gather}
\|A^{-1}_u\|_{\tnp}\ge \max_{-N\le k\le N}|\widehat{\mu}(k)|= \max_{-N\le k\le N}|p_u(k)|=\max_{-N\le k\le N}\biggl\{\frac1{|\widehat{u}(k)|}\biggr\}\nonumber\\
=\frac1{\min_{-N\le k\le N}|\widehat{u}(k)|}
\ge\frac1{\widehat{u}(N)}=p_u(N)=\sum_{M=-N+1}^N \tau_me^{i\pi mN/N}=\nonumber\\
\sum_{M=-N+1}^N (-1)^m\tau_m=\sum_{M=-N+1}^N |\tau_m|=\|\mu\|.
\end{gather}
Other hands, from Minkowski's inequality (14) it is easily to see that
\begin{gather}
\|A^{-1}_u(s)\|_{\tnp}=\biggl(\int_{\T}\biggl|\int_{\T}s(x-t)\,d\mu(t)\biggr|^p\,dx\biggr)^{1/p}\nonumber\\
\le \biggl(\int_{\T}\biggl|\int_{\T}|s(x-t)|\,d|\mu|(t)\biggr|^p\,dx\biggr)^{1/p}
\nonumber\\
\le \int_{\T}\biggl(\int_{\T}|s(x-t)|^p\,dx\biggr)^{1/p}\,d|\mu|(t)= \|s\|_{\tnp}\|\mu\|,
\end{gather}
where $|\mu|$ denote the variation of $\mu$. Thus, (35) with (36) show   (33), and the Proposition 5  is proved.

PROOF OF THEOREM 1.\quad
Assume that $u\in L^q(\T)$, $\|u\|_{L^q(\T)}=1$ and $u$ satisfies the conditions of Proposition 4.  If $s\in\tn$, then by H\"{o}lder's inequality, we get
\begin{gather}
|s(x)|^p= \biggl|\int_{\T}A^{-1}_u(s)(x-t)u(t)\,dt\biggr|^p= \biggl|\int_{-\pi/2N}^{\pi/2N}A_u^{-1}(s)(x-t)u(t)\,dt\biggr|^p\nonumber\\
\le \biggl(\int_{-\pi/2N}^{\pi/2N}\Bigr|A_u^{-1}(s)(x-t)\biggl|^p\,dt\biggr) \biggl(\int_{-\pi/2N}^{\pi/2N}|u(t)|^q\,dt\biggr)^{p/q}\nonumber\\
= \int_{-\pi/2N}^{\pi/2N}\Bigr|A^{-1}_u(s)(x-t)\biggl|^p\,dt\nonumber
\end{gather}
for each $x\in\T$. Therefore,
\begin{gather}
\sum_{j=1}^r|s(x_j)|^p\le \sum_{j=1}^r\int_{-\pi/2N}^{\pi/2N}\biggl|A_u^{-1}(s)(x_j-t)\biggr|^p\,dt\nonumber\\ =\sum_{j=1}^r\int_{x_j-\pi/2N}^{x_j+\pi/2N}\biggl|A_u^{-1}(s)(y)\biggr|^p\,dy.
\end{gather}
We will denote by $E_j$ the set
\[
E_j=\Bigl(x_j-\frac{\pi}{2N} , x_j+\frac{\pi}{2N}\Bigr],
\]
$ j=1,\dots, r$. Using the fact that $A_u^{-1}(s)$ is a trigonometric polynomial, i.e., a continuous and periodical function on $\R$, we conclude that
\begin{equation}
\int_{x_j-\pi/2N}^{x_j+\pi/2N}\Bigl|A_u^{-1}(s)(y)\Bigr|^p\,dy  = \int_{E_j}\Bigl|A_u^{-1}(s)(y)\Bigr|^p\,dy.
\end{equation}
 Let $x\in\R$ and assume that $x\in E_j$, for $j=i, i+1,\dots, i+k$ with some  non-negative integer $k$.  Then we claim that:

 (i) \ if $\pi/N\delta\in\Z$, then
 \begin{equation}
k+1\le \frac{\pi}{N\delta};
\end{equation}
 (ii) \ if $\pi/N\delta\not\in\Z$, then
\begin{equation}
k\le  \Bigl[\frac{\pi}{N\delta}\Bigr],
\end{equation}
where $[\cdot]$ is the usual integer part of a  real number.

Indeed, since $x\in E_i\cap E_{i+k}$, it follows that $x_i+\pi/2N> x_{i+k}- \pi/2N$. Then
\begin{equation}
x_{i+k}-x_i< \frac{\pi}N.
\end{equation}
Other hands, by (1)  we get
\begin{equation}
x_{i+k}-x_i\ge k\delta.
\end{equation}
Combining (42) with (41), we see that
\begin{equation}
k< \frac{\pi}{N\delta}.
\end{equation}
Now note  that $k$ is an integer. Hence if $\pi/N\delta\in\Z$, then (43) implies (39). For $\pi/N\delta\not\in\Z$, we get (40), which yields our  claim.

Thus, each $x\in\T$ can  belong to at most  $\sigma_p$ intervals $E_j$, $j=1,\dots, p$,  where $\sigma_p$ was defined by (8).

Now combining (35), (36) with (37), (38), (31) and (8), we get
 \begin{gather}
\sum_{j=1}^r|s(x_j)|^p\le\sum_{j=1}^r\int_{x_j-\pi/2N}^{x_j+\pi/2N}\biggl|A_u^{-1}(s)(y)\biggr|^p\,dy\le \sigma(\delta; N)\int_{\T}\Bigl|A^{-1}_u(s)(y)\Bigr|^p\,dy=\nonumber\\
= \sigma(\delta; N) \|A^{-1}_u(s)\|^p_{L^p(\T)}\le \sigma(\delta; N)\|A_u^{-1}\|^p_{T^p_N }\|s\|^p_{L^p(\T)}=
 \frac{\sigma_p}{\bigl(\hat{u}(N)\bigr)^p}\|s\|^p_{L^p(\T)}.
\end{gather}
Next, the estimate (44) can be improved to
\begin{equation}
\sum_{j=1}^r|s(x_j)|^p\le \frac{\sigma(\delta; N)}{\sup_{u}\Bigl(\hat{u}(N)\Bigr)^p}\ \|s\|^p_{L^p(\T)},
\end{equation}
where the supremum extends over all admissible $u$ as described in the statement of Proposition 4. We claim that for such an $u$ we have
\begin{equation}
\sup_{u}(\widehat{u}(N))^p=\|\cos Nx\|^p_{L^p[-\pi/2N, \pi/2N]}=\frac{\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}}{N}.
\end{equation}
Indeed, H\"{o}lder's inequality implies that
\begin{gather}
|\widehat{u}(N)|^p= \Bigl|\int_{\T}u(x)\cos Nx\  dx\Bigr|^p= \Bigl|\int_{-\pi/2N}^{\pi/2N}u(x)\cos Nx\  dx\Bigr|^p\nonumber\\
\le \int_{-\pi/2N}^{\pi/2N}|\cos Nx|^p\  dx \cdot \|u\|_{L^q[-\pi/2N, \pi/2N]}^p=\int_{-\pi/2N}^{\pi/2N}|\cos Nx|^p\,dx\nonumber\\
=\frac1{N}\int_{-\pi/2}^{\pi/2}|\cos t|^p\, dt=\frac{\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}}{N}.
\end{gather}
 Moreover, since we used H\"{o}lder's inequality, it follows that the estimate (47) is exact and the equality is attained if
\[
u(x)= \theta\cos ^{p-1} Nx \cdot\chi_{[-\pi/2N, \pi/2N]}(x),
\]
 where $\chi_{[-\pi/2N, \pi/2N]}$ is the indicator function of the interval $[-\pi/2N, \pi/2N]$  and $\theta \in\Co$ is such that   $\|u\|_{L^q[-\pi/2N, \pi/2N]}=1$. Therefore, our claim (46) is proved.

 Combining (45) with (46), we get
 \begin{equation}
 \sum_{j=1}^r|s(x_j)|^p\le\frac{\sigma(\delta; N)\cdot N}{\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}}\|s\|^p_{L^p(\T)}.
 \end{equation}
 Next,
  \begin{equation}
\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}=\int_{-\pi/2}^{\pi/2} \cos^p t\, dt=B\Bigl(\frac12; \frac{p+1}{2}\Bigr),
 \end{equation}
 (see, e.g., \cite[p. 142]{7}), where $B$ is Euler's beta function defined by
 \[
 B(a;b)=2\int_0^{\pi/2}\sin^{2a-1}\theta\cos^{2b-1}\theta\, d\theta
 \]
 for $\Re a$, $\Re b>0$.  Applying the following connection between the beta and the usual gamma function  (\cite[p. 142]{7}
 \[
 B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)},
 \]
 we conclude from (49) that
  \begin{equation}
\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}=\frac{\Gamma\Bigl(\frac12\Bigr)\Gamma\Bigl(\frac{p}2+\frac12\Bigr)}{\Gamma\Bigl(\frac{p}2 +1\Bigr)}
 \end{equation}
Since $\Gamma\Bigl(p/2+1\Bigr)=\Gamma\Bigl(p/2\Bigr)\cdot p/2$ and $\Gamma\Bigl(1/2\Bigr)=\sqrt{\pi}$ (see, e.g., \cite[p.p.  137-138]{7}), it follows from (49) and (50) that
\[
\|\cos t\|^p_{L^p[-\pi/2, \pi/2]}=\frac{2\sqrt{\pi}\cdot \Gamma\Bigl(\frac{p}2+\frac12\Bigr)}{\Gamma\Bigl(\frac{p}2\Bigr)}.
 \]
Substituting this into (48), we obtain (7). Theorem 1 is proved.

POOF OF COROLLARY 2. \quad It is known that for a nonnegative integer $n$,
\[
\Gamma(n+1)=n! \qquad {\text{and}} \qquad \Gamma\Bigl(n+\frac12\Bigr)=\sqrt{\pi}\cdot\frac{1\cdot 3\cdot 5\cdot \cdots \cdot(2-1)}{2^n}
\]
 (see \cite[p.  139]{7}). Using this, the representations (9) and (10) can be verified by straightforward calculation.


\section{1. Experimental Methodology} \label{sec:experimental_methodology}

\subsection{Preparation of lipid-coated particles}

Fluorescently labeled, lipid-coated particles were created by coating silica micro-beads with a supported lipid bilayer (SLB) containing a minority fraction of fluorescently tagged lipid.
1,2-dioleoyl-sn-glycero-3-phos-phocholine (DOPC) and 1,2-dioleoyl-sn-glycero-3-phospho-L-serine (DOPS) were purchased from Avanti Polar Lipids.
Atto 647-1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE-Atto 647) was purchased from ATTO-TEC GmbH.
Silica microspheres (diameter 2.5 $\mu$m; catalog code: SS05000) were purchased from Bangs Laboratories.
Small unilamellar vesicles (SUVs) were formed using an established sonication method \cite{bakalar2018size}.
In brief, a lipid film containing DOPC, 5\% DOPS, and 0.5\% DOPE-Atto 647 was dried under nitrogen and then under vacuum for 30 minutes.
The film was rehydrated in Milli-Q (MQ) water to 0.2 mg/mL lipids, sonicated at low power using a tip sonicator (Branson SFX250 Sonifier) at 20\% of maximum, 1 s/2 s on/off, for three minutes.
MOPS buffer was added at a final concentration of 50 mM MOPS, pH 7.4, 100 mM NaCl to the resulting SUV mixture.

Silica microspheres were cleaned using a 3:2 mixture of sulfuric acid:hydrogen peroxide (Piranha) for 30 minutes in a bath sonicator, spun at 1000 g, and washed 3 times before being resuspended in MQ water.
To form SLBs on the beads, 50 $\mu$L of SUV solution was mixed with 10 $\mu$L of the cleaned bead suspension.
The bead/SUV mixture was incubated for 15 minutes at room temperature while allowing the beads to sediment to the bottom of the centrifuge tube.
Beads were washed 5 times with MQ water by gently adding/removing the liquid without resuspending the beads into solution.
The fluidity of the SLB was verified by imaging beads on a glass coverslip at high laser intensity, where the diffusion of labeled lipids was visible after photo-bleaching a small region.
Lipid-coated beads were deposited into a chamber containing MQ water and sealed off to eliminate drift.
The beads settled down to the bottom of the chamber and all experiments were conducted in 2D.

\subsection{Optical tweezer setup and calibration}

An array of moving harmonic traps was generated using optical tweezers (Tweez 305, Aresis Ltd; Ljubljana, Slovenia), using an IR laser (1064 nm) with a maximum power of 5 W continuous wave (CW).
We selected a trap-to-trap switching rate of 100 kHz to ensure that the particles will effectively feel a continuous harmonic potential.
We used a 16 $\times$ 16 array of traps, which results in $\approx 2.5$ ms time delay to illuminate all trap positions.
This time delay is significantly smaller than the Brownian and oscillatory convection timescales in our system, ensuring that the particles experience a continuous harmonic potential.
A custom MATLAB script was written to construct a time trajectory of oscillatory trap positions for each cell lattice position and incorporated into the tweezer software.
The trap focus was adjusted to the mid-plane of the colloids sitting above the substrate.
Laser powers were adjusted from 0.05-0.5 W to vary the trap stiffness from $\kappa = 0.5$-6 $kT/\mu\mathrm{m}^2$.

The trap stiffness $\kappa$ was calibrated by measuring the equilibrium probability distribution of the particles in a stationary array of traps.
For each laser power, $\kappa$ was obtained by binning particles by their radial position $r$ from the center of the trap and fitting the binned data to a Boltzmann distribution, $P(r) = (\kappa/2\pi)\mathrm{e}^{-\kappa r^{2}/(2kT)}$.
An example of a distribution and fit is shown in Fig.~\ref{Fig:SI1}.
We verified that there are no variations in trip stiffness between different lattice positions in the array.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.7\linewidth]{Figs/SI/SIFig_Trap_Stiffness.jpg}
		\caption{%
			Measurement of trap stiffness $\kappa$ from the equilibrium probability distribution of particles diffusing in a harmonic well generated by optical tweezers.
			Data are fit to a Boltzmann distribution to obtain $\kappa$ ($\kappa = 4~kT /\mu\mathrm{m}^2$ in the case shown).
			This measurement was averaged over all 16 $\times$ 16 trap positions in the lattice array and repeated for every laser power used in this study.
		}
		\label{Fig:SI1}
	\end{center}
	\vspace{-18pt}
\end{figure}

The trap width $W_{\text{trap}}$ was determined from a separate set of experiments.
Two traps were placed side-by-side with center-to-center separation distance $W$.
The first trap, containing a trapped particle, was held fixed while the position of the second trap was varied; the average position $\langle x_{i}(t) \rangle$ of the particle was measured as a function of the separation distance $W$ (Fig.~\ref{Fig:SI2}).
When the second trap is placed far away, no interference is observed on the average position of the particle.
However, as the second trap is moved closer, $W < 3$ $\mu$m for a particle of radius $a = 1.25$ $\mu$m, the average position drifts towards the second trap.
We found that the average particle position remains approximately constant within the range of separation distances of $W =$ 3-3.5 $\mu$m, giving an approximate trap width $W_{\text{trap}} \approx 3.2$ $\mu$m.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.7\linewidth]{Figs/SI/SIFig_Trap_Width.jpg}
		\caption{%
			Measurement of trap width $W_{\text{trap}}$.
			A second trap was placed at varying separation distances from the first trap containing a trapped bead.
			We measured the time-averaged position of the trapped bead, $\langle x_i(t) \rangle$, for varying separation distances at fixed trap stiffness.
			We found that the average position is pulled towards the second trap at distances $W < 3$ $\mu$m and is approximately constant in the range $W = 3$-3.5 $\mu$m. This gives an average trap width $W_{\text{trap}} \approx 3.2$ $\mu$m.
		}
		\label{Fig:SI2}
	\end{center}
	\vspace{-18pt}
\end{figure}

\subsection{Measurement of diffusivity}

The long-time self diffusivity was determined by particle tracking. All imaging was carried out on an inverted Nikon Ti2-Eclipse microscope (Nikon Instruments) using a water-immersion objective (Plan Apochromat VC 60x, numerical aperture 1.2, water).
Lumencor SpectraX Multi-Line LED Light Source was used for excitation (Lumencor, Inc).
Fluorescent light was spectrally filtered with an emission filter (680/42; Semrock, IDEX Health and Science) and imaged on a Photometrics Prime 95 CMOS Camera (Teledyne Photometrics).
In order to achieve satisfactory long-time statistics, particle trajectories were measured for times much larger than all other timescales in the system (including the diffusive timescale $\gamma L^{2}/kT$, oscillation period $2\uppi/\omega$, and trapping timescale $\gamma/\kappa$). A modified MATLAB script, based on the IDL code by Crocker and Grier \cite{crocker1996methods,crockerweeksIDL,blairdufresneMATLAB}, was used to track the individual particles by identifying each particle center and tracking its trajectory over time using an image stack with one frame taken every 1-2 s.
Particles that were immobile (due to defects) were filtered out so as not to be considered during image post-processing.

The average diffusivity tensor is classically defined in terms of the long-time derivative of the mean squared displacements (MSD) of the particles:
\begin{equation}
    \overline{\tens{D}} = \lim_{t \rightarrow \infty} \frac{1}{2}
    \frac{\mathrm{d}}{\mathrm{d}t}
    \langle \Delta \bm{R}(t) \Delta \bm{R}(t) \rangle
    ,
    \label{eq:SI_diffusivity_tensor}
\end{equation}
where $\bm{R}$ denotes the {\it global} position vector [related to the {\it local} position vector $\bm{r}$ by Eq. \eqref{eq:SI_coordinate_conversion}, below] and the angle brackets $\langle\,\cdot\,\rangle$ denote an {\it ensemble} average (not to be confused with the {\it cell} average defined in the main text). The MSD tensor over a time interval $t$ is computed from the formula,
\begin{equation}
	\braket{
	\Delta \bm{R} (t)
	\Delta \bm{R} (t)
	}
	=
	\frac{1}{N_{\text{p}}}
	\sum_{i=1}^{N_{\text{p}}}
	\lim_{\tau\rightarrow\infty}
	\frac{1}{\tau-t}
	\int_{0}^{\tau-t}
	\left[ \bm{R}_{i}(s + t) - \bm{R}_{i}(s)  \right]
	\left[ \bm{R}_{i}(s + t) - \bm{R}_{i}(s)  \right]
	\, \mathrm{d} s
	,
	\label{eq:SI_msd_experiments}
\end{equation}
where $\bm{R}_{i}(t)$ denotes the global position of the $i$th particle at time $t$.
In Eq.~\eqref{eq:SI_msd_experiments}, the squared displacement of a particle with index $i$ is first averaged over all time windows of duration $t$ within the interval $\tau$ of the particle's trajectory. This ``time average'' for each $i$th particle, evaluated in the limit as $\tau \rightarrow \infty$, is subsequently averaged over all particles $i = 1, 2, \dots, N_{\text{p}}$ to approximate the ensemble average of all squared displacements with satisfactory statistics.
At long times, the MSD tensor $\langle \Delta \bm{R}(t) \Delta \bm{R}(t) \rangle$ oscillates with fixed amplitude about a steady, linear growth. Thus, the long-time derivative of the MSD can be measured by simply dividing by time, leading to the relation,
\begin{equation}
    \overline{\tens{D}} = \lim_{t \rightarrow \infty} \frac{1}{2 t}
     \langle \Delta \bm{R}(t) \Delta \bm{R}(t) \rangle
     .
     \label{eq:SI_diffusivity_experiments}
\end{equation}
Equation \eqref{eq:SI_diffusivity_experiments} was used to measure the diffusivity from the measured particle trajectories (see Fig.~\ref{Fig:SI3}). Trajectories were averaged over a sufficiently long time interval $\tau$ to ensure linear growth, and the time integral in Eq.~\eqref{eq:SI_msd_experiments} was discretized using the left Riemann sum.
Statiscal errors in the MSD were calculated using a bootstrap algorithm \cite{ross2020introduction}.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.5\linewidth]{Figs/SI/SIFig_MSD_Example.jpg}
		\caption{%
			Representative mean squared displacements $\langle \Delta x(t) \Delta x(t) \rangle / (2t)$ (black symbols) and $\langle \Delta y(t) \Delta y(t) \rangle / (2t)$ (blue symbols) measured using Eq.~\eqref{eq:SI_msd_experiments} for Brownian particles diffusing through an oscillating array of harmonic traps.
			Diffusivities reported in the main text were computed from the long-time plateaus of these curves, using Eq.~\eqref{eq:SI_diffusivity_experiments}. Statistical errors were calculated using the bootstrap algorithm \cite{ross2020introduction} over the entire observation time window.
		}
		\label{Fig:SI3}
	\end{center}
	\vspace{-18pt}
\end{figure}

The particle resistivity $\gamma$ used in all theoretical calculations was calibrated by measuring the Stokes-Einstein-Sutherland diffusivity $D_{0} = kT/\gamma \approx 0.105$ $\mu$m$^2$/s of particles diffusing in the absence of a harmonic potential. For a spherical particle of radius $a$ in a fluid of viscosity $\eta$, the particle resistivity is given by $\gamma = 6 \uppi \eta a K_{D}$, where $K_{D}$ is a drag-correction factor to account for the hydrodynamic interaction with a nearby wall (in our case, the substrate floor). For our system with $a = 1.25$ $\mu$m and $\eta = 1$ cP, we estimate the drag-correction factor to be $K_{D} = kT/(6\uppi \eta a D_{0}) \approx 1.63$, corresponding to a particle-to-wall spacing of about 0.5 $\mu$m according to Fax\'en's formula \cite{happel2012low}. This gives a particle resistivity of $\gamma \approx 9.49$ $kT\cdot\text{s}/\mu\text{m}^2$.

\section{2. Taylor-Dispersion Theory} \label{sec:dispersion_theory}

\subsection{Derivation of Eqs. (3)-(5): governing equations for the probability density and displacement}

The starting point for deriving the basic equations in the main text is the single-particle Smoluchowski equation,
\begin{equation}
	\frac{\partial P(\bm{R}, t)}{\partial t}
	=
	-
	\bm{\nabla}_{\bm{R}} \cdot \bm{J} (\bm{R}, t)
	,
	\label{eq:SI_smoluchowski_eqn}
\end{equation}
where $P(\bm{R}, t)$ is the probability density of finding a Brownian particle at a {\it global} position $\bm{R}$ and time $t$ and
\begin{equation}
	\bm{J} (\bm{R}, t)
	=
	\bm{u} (t) P
	-
	\frac{1}{\gamma}
	[
	kT \bm{\nabla}_{\bm{R}} P
	+
	P \bm{\nabla}_{\bm{R}} V (\bm{R})
	]
	\label{eq:SI_probability_flux}
\end{equation}
is the probability flux. The spatial periodicity of the potential-energy field allows us to convert the ``global'' position $\bm{R}$ to the ``local'' position $\bm{r}$ via the transformation,
\begin{equation}
	\bm{R} = \bm{n} L + \bm{r},
	\label{eq:SI_coordinate_conversion}
\end{equation}
where $\bm{n}$ contains the lattice indices of a given periodic cell. In terms of lattice and local coordinates, $V(\bm{R}) \equiv V(\bm{r})$, $P(\bm{R},t) \equiv P_{\bm{n}} (\bm{r}, t)$, and $\bm{J} (\bm{R}, t) \equiv \bm{J}_{\bm{n}} (\bm{r},t)$.

In the following, we employ the ``flux-averaging'' approach of Brady and coworkers \cite{morris1996self,Zia2010,Takatori2014,burkholder2017tracer,burkholder2019fluctuation,peng2020upstream}.
First, we define the continuous wavevector $\bm{k}$ and apply the discrete Fourier transform $\hat{(\,\cdot\,)} \equiv \sum_{\bm{n}} (\, \cdot\,) \mathrm{e}^{\mathrm{i} \bm{k} \cdot \bm{n} L}$ to Eqs.~\eqref{eq:SI_smoluchowski_eqn}-\eqref{eq:SI_probability_flux}, obtaining
\begin{equation}
	\frac{\partial \hat{P}(\bm{k}, \bm{r}, t)}{\partial t}
	=
	-
	(\mathrm{i} \bm{k} + \bm{\nabla}_{\bm{r}})  \cdot \hat{\bm{J}} (\bm{k}, \bm{r}, t)
	,
	\label{eq:SI_smoluchowski_eqn_ft}
\end{equation}
\begin{equation}
	\hat{\bm{J}} (\bm{k}, \bm{r}, t)
	=
	\bm{u} (t) \hat{P}
	-
	\frac{1}{\gamma}
	[
	kT (\mathrm{i} \bm{k} + \bm{\nabla}_{\bm{r}}) \hat{P}
	+
	\hat{P} \bm{\nabla}_{\bm{r}} V (\bm{r})
	]
	.
	\label{eq:SI_probability_flux_ft}
\end{equation}
Next, we spatially average Eqs.~\eqref{eq:SI_smoluchowski_eqn_ft}-\eqref{eq:SI_probability_flux_ft} over one periodic cell according to $\braket{\,\cdot\,} \equiv L^{-2}\int_{L^{2}} (\,\cdot\,) \, \mathrm{d} \bm{r}$, apply the divergence theorem, and invoke periodic boundary conditions to obtain the continuity equation,
\begin{equation}
	\frac{\partial \hat{\rho}(\bm{k}, t)}{\partial t}
	=
	-
	\mathrm{i} \bm{k} \cdot \braket{\hat{\bm{J}}} (\bm{k}, t)
	,
	\label{eq:SI_smoluchowski_eqn_ft_avg}
\end{equation}
\begin{equation}
	\braket{\hat{\bm{J}}} (\bm{k}, t)
	=
	\bm{u} (t) \hat{\rho}
	-
	\frac{1}{\gamma}
	[
	kT \mathrm{i} \bm{k} \hat{\rho}
	+
	\braket{\hat{P} \bm{\nabla}_{\bm{r}}V}
	]
	,
	\label{eq:SI_probability_flux_ft_avg}
\end{equation}
where $\hat{\rho} (\bm{k}, t) \equiv \braket{\hat{P}} (\bm{k}, t)$ is the Fourier-transformed number density. Eqs.~\eqref{eq:SI_smoluchowski_eqn_ft_avg}-\eqref{eq:SI_probability_flux_ft_avg} represent the macroscopic transport equations for the periodic lattice.

Next, we define the structure function $\hat{G}(\bm{k}, \bm{r}, t)$ as
\begin{equation}
	\hat{P}(\bm{k}, \bm{r}, t)
	=
	\hat{\rho} (\bm{k}, t)
	\hat{G}(\bm{k}, \bm{r}, t)
	.
	\label{eq:SI_conditional_probability}
\end{equation}
Multiplying Eq.~\eqref{eq:SI_smoluchowski_eqn_ft_avg} by $\hat{G}$, subtracting from Eq.~\eqref{eq:SI_smoluchowski_eqn_ft}, and dividing through by $\hat{\rho}$ then gives
\begin{flalign}
	\frac{\partial \hat{G}(\bm{k}, \bm{r}, t)}{\partial t}
	&=
	-
	\hat{\rho}^{-1} [
	\mathrm{i} \bm{k} \cdot ( \hat{\bm{J}} - \braket{\hat{\bm{J}}} \hat{G} )
	+
	\bm{\nabla}_{\bm{r}} \cdot \hat{\bm{J}}
	]
	\nonumber\\
	&=
	- \bm{u} (t) \cdot \bm{\nabla}_{\bm{r}} \hat{G}
	+
	\frac{kT}{\gamma}
	\nabla_{\bm{r}}^{2} \hat{G}
	+
	\frac{1}{\gamma}
	\bm{\nabla}_{\bm{r}} \cdot
	[
	\hat{G} \bm{\nabla}_{\bm{r}} V (\bm{r})
	]
	+
	\mathrm{i} \bm{k} \cdot
	\left(
	\frac{2kT}{\gamma} \bm{\nabla}_{\bm{r}} \hat{G}
	+
	\frac{1}{\gamma}
	[
	\hat{G} \bm{\nabla}_{\bm{r}} V (\bm{r})
	-
	\braket{\hat{G} \bm{\nabla}_{\bm{r}}V} \hat{G}
	]
	\right)
	,
	\label{eq:SI_structure_field_eqn}
\end{flalign}
where in the last line we have substituted Eqs.~\eqref{eq:SI_probability_flux_ft}, \eqref{eq:SI_probability_flux_ft_avg}, and \eqref{eq:SI_conditional_probability}. Taylor-expanding $\hat{G}$ about $\bm{k} = \bm{0}$,
\begin{equation}
	\hat{G}(\bm{k}, \bm{r}, t)
	=
	g(\bm{r}, t)
	+
	\mathrm{i} \bm{k} \cdot \bm{d} (\bm{r}, t)
	+
	\cdots
	,
	\label{eq:SI_taylor_series}
\end{equation}
substituting the expansion into Eq.~\eqref{eq:SI_structure_field_eqn}, and collecting terms of like order in $\mathrm{i} \bm{k}$ yields the ordered set of equations,
\begin{flalign}
	\frac{\partial g(\bm{r}, t)}{\partial t}
	+
	\bm{u} (t) \cdot \bm{\nabla}_{\bm{r}} g
	-
	\frac{kT}{\gamma}
	\nabla_{\bm{r}}^{2} g
	-
	\frac{1}{\gamma}
	\bm{\nabla}_{\bm{r}} \cdot
	[
	g \bm{\nabla}_{\bm{r}} V (\bm{r})
	]
	=
	0
	,
\end{flalign}
\begin{flalign}
	\frac{\partial \bm{d}(\bm{r}, t)}{\partial t}
	+
	\bm{u} (t) \cdot \bm{\nabla}_{\bm{r}} \bm{d}
	-
	\frac{kT}{\gamma}
	\nabla_{\bm{r}}^{2} \bm{d}
	-
	\frac{1}{\gamma}
	\bm{\nabla}_{\bm{r}} \cdot
	[
	\bm{d}\bm{\nabla}_{\bm{r}} V (\bm{r})
	]^{\dag}
	=
	\frac{2kT}{\gamma} \bm{\nabla}_{\bm{r}} g
	+
	\frac{1}{\gamma}
	[
	g\bm{\nabla}_{\bm{r}} V (\bm{r})
	-
	\braket{g \bm{\nabla}_{\bm{r}}V} g
	]
	.
\end{flalign}
The last two equations are exactly Eqs.~\eqref{eq:g_eqn} and \eqref{eq:d_eqn} from the main text. Conservation of probability requires the $g$- and $\bm{d}$-fields to satisfy periodic boundary conditions as well as the normalization conditions $\braket{g} = 1$ and $\braket{\bm{d}} = \bm{0}$.

\subsection{Derivation of Eqs. (6)-(7): effective drift velocity and diffusivity}

The effective drift velocity $\bm{U}(t)$ and diffusivity $\tens{D} (t)$ of the Brownian particle are related to the Fourier-transformed, average flux $\braket{\hat{\bm{J}}}$ via the large-wavelength expansion,
\begin{equation}
	\braket{\hat{\bm{J}}} (\bm{k}, t)
	=
	\hat{\rho}
	\left[
		\bm{U} (t)
		-
		\mathrm{i} \bm{k} \cdot \tens{D} (t)
		+
		\cdots
	\right]
	.
	\label{eq:SI_average_flux_ft_eqn1}
\end{equation}
In order to derive expressions for $\bm{U}$ and $\tens{D}$, we insert Eqs.~\eqref{eq:SI_conditional_probability} and \eqref{eq:SI_taylor_series} into \eqref{eq:SI_probability_flux_ft_avg}, obtaining
\begin{flalign}
	\braket{\hat{\bm{J}}} (\bm{k}, t)
	&=
	\hat{\rho}
	\left(
	\bm{u} (t)
	-
	\frac{1}{\gamma}
	[
	kT \mathrm{i} \bm{k}
	+
	\braket{\hat{G} \bm{\nabla}_{\bm{r}}V}
	]
	\right)
	\nonumber\\
	&=
	\hat{\rho}
	\left[
	\bm{u} (t) - \frac{1}{\gamma} \braket{g \bm{\nabla}_{\bm{r}}V}
	-
	\mathrm{i} \bm{k}
	\left(
	\frac{kT}{\gamma} \tens{I}
	+
	\frac{1}{\gamma}\braket{\bm{d} \bm{\nabla}_{\bm{r}}V}
	\right)
	+
	\cdots
	\right]
	.
	\label{eq:SI_average_flux_ft_eqn2}
\end{flalign}
Equating terms of like order in $\mathrm{i} \bm{k}$ in Eqs.~\eqref{eq:SI_average_flux_ft_eqn1} and \eqref{eq:SI_average_flux_ft_eqn2} furnishes the expressions,
\begin{equation}
	\bm{U}(t)
	=
	\bm{u} (t)
	-
	\frac{1}{\gamma} \braket{g \bm{\nabla}_{\bm{r}} V} (t)
	,
\end{equation}
\begin{equation}
	\tens{D}(t)
	=
	\frac{kT}{\gamma} \tens{I}
	+
	\frac{1}{\gamma}\braket{\bm{d} \bm{\nabla}_{\bm{r}}V} (t)
	,
\end{equation}
which are exactly Eqs.~\eqref{eq:effective_drift}-\eqref{eq:effective_diffusivity} in the main text.

\section{3. Numerical Method} \label{sec:numerical_method}

Eqs.~\eqref{eq:g_eqn} and \eqref{eq:d_eqn} were solved using the finite-element method in COMSOL Multiphysics$^\text{\textregistered}$ (Version 5.5) with the ``Coefficient Form PDE'' physics interface. An $L \times L$ square cell was set up and discretized into triangular elements (Fig.~\ref{Fig:SIMeshes}). Periodic boundary conditions were applied to the $g$- and $\bm{d}$-fields at the edges of the cell. Studies were run using both time-dependent ($\bm{u}\not=\bm{0}$) and stationary ($\bm{u}=\bm{0}$) solvers. For the time-dependent studies, the $g$- and $\bm{d}$-fields were initialized to uniform values $1$ and $\bm{0}$, respectively, and time-advanced using the backward differentiation formula with a timestep $\Delta t = 0.001(2\uppi/\omega)$ until a periodic steady state was achieved. The number of periods needed to reach steady state generally increased with the oscillation frequency. For the stationary studies, the equations were solved iteratively using Newton's method and the normalization conditions $\braket{g} = 1$ and $\braket{\bm{d}} = \bm{0}$ were implemented as weak-form constraints. Upon solving for the $g$- and $\bm{d}$-fields, Eqs.~\eqref{eq:effective_drift} and \eqref{eq:effective_diffusivity} were evaluated using a fourth-order domain integration method and (in the time-dependent studies) subsequently time-averaged over the final oscillation period.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.7\linewidth]{Figs/SI/SIFig_Meshes.png}
		\caption{%
			Triangular meshes used for the finite-element calculations. Meshes containing (a) 1132 elements (for the time-dependent studies) and (b) 29,018 elements (for the stationary studies) were used. Coarser meshes were used in the time-dependent calculations to save computational time.
		}
		\label{Fig:SIMeshes}
	\end{center}
	\vspace{-18pt}
\end{figure}

\section{4. Asymptotic Limits} \label{sec:asymptotic_limits}

\subsection{Derivation of Eq. (8): stationary traps with shallow potential wells}

If the harmonic traps held in a fixed configuration, $\bm{u} = \bm{0}$ and the $g$- and $\bm{d}$-fields achieve a steady state. Equations \eqref{eq:g_eqn} and \eqref{eq:d_eqn} then simplify to
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} g (\bm{r})
	+
	\bm{\nabla}_{\bm{r}} \cdot
	[
	g(\bm{r}) \bm{\nabla}_{\bm{r}} V (\bm{r})
	]
	=
	0
	,
	\label{eq:SI_g_eqn_steady}
\end{flalign}
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} \bm{d}(\bm{r})
	+
	\bm{\nabla}_{\bm{r}} \cdot
	[
	\bm{d}(\bm{r})\bm{\nabla}_{\bm{r}} V (\bm{r})
	]^{\dag}
	=
	-
	2kT \bm{\nabla}_{\bm{r}} g
	-
	g\bm{\nabla}_{\bm{r}} V
	+
	\braket{g \bm{\nabla}_{\bm{r}}V} g
	.
	\label{eq:SI_d_eqn_steady}
\end{flalign}
Eq.~\eqref{eq:SI_g_eqn_steady} may be solved subject to the constraint $\braket{g} = 1$ to get the Boltzmann distribution,
\begin{equation}
	g (\bm{r})
	=
	\frac{\mathrm{e}^{-V(\bm{r})/kT}}{\braket{\mathrm{e}^{-V/kT}}}
	.
	\label{eq:SI_g_field_steady}
\end{equation}
The governing equation for the $\bm{d}$-field, Eq.~\eqref{eq:SI_d_eqn_steady}, then simplifies to
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} \bm{d}
	+
	\bm{\nabla}_{\bm{r}} \cdot
	(
	\bm{d}\bm{\nabla}_{\bm{r}} V
	)^{\dag}
	=
	-
	kT \bm{\nabla}_{\bm{r}} g
	.
	\label{eq:SI_d_eqn_steady_2}
\end{flalign}

Eq.~\eqref{eq:SI_d_eqn_steady_2} cannot be solved analytically in general. However, for ``shallow'' potential wells, $\Delta V \ll kT$, we may Taylor-expand Eq.~\eqref{eq:SI_g_field_steady} as
\begin{equation}
	g
	=
	1
	-
	\frac{V - \braket{V}}{kT}
	+
	\frac{V^{2} - \braket{V^{2}} - 2 \braket{V} (V - \braket{V})}{2(kT)^{2}}
	+
	\cdots
	,
\end{equation}
so that Eq.~\eqref{eq:SI_d_eqn_steady_2} becomes
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} \bm{d}
	+
	\bm{\nabla}_{\bm{r}} \cdot
	(
	\bm{d}\bm{\nabla}_{\bm{r}} V
	)^{\dag}
	=
	\left( 1 - \frac{V - \braket{V}}{kT} + \cdots \right)\bm{\nabla}_{\bm{r}} V
	.
	\label{eq:SI_d_eqn_steady_3}
\end{flalign}
To solve Eq.~\eqref{eq:SI_d_eqn_steady_3}, we expand the $\bm{d}$-field in a perturbation series,
\begin{equation}
	\bm{d} (\bm{r})
	=
	\bm{d}^{(0)} (\bm{r})
	+
	\bm{d}^{(1)} (\bm{r})
	+
	\cdots
	,
	\label{eq:SI_d_field_perturbation_series}
\end{equation}
where $\bm{d}^{(0)} = O(\Delta V/kT)$, $\bm{d}^{(1)} = O[\Delta V^{2}/(kT)^{2}]$, and so on. Inserting Eq.~\eqref{eq:SI_d_field_perturbation_series} into \eqref{eq:SI_d_eqn_steady_3} and collecting terms of like order in $\Delta V / kT$ yields the ordered set of equations,
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} \bm{d}^{(0)}
	=
	\bm{\nabla}_{\bm{r}} V
	,
	\label{eq:SI_d_eqn_steady_4a}
\end{flalign}
\begin{flalign}
	kT
	\nabla_{\bm{r}}^{2} \bm{d}^{(1)}
	=
	-
	\bm{\nabla}_{\bm{r}} \cdot
	(
	\bm{d}^{(0)}\bm{\nabla}_{\bm{r}} V
	)^{\dag}
	-
	\frac{V - \braket{V}}{kT} \bm{\nabla}_{\bm{r}} V
	,
	\label{eq:SI_d_eqn_steady_4b}
\end{flalign}
subject to the constraints $\braket{\bm{d}^{(0)}} = \bm{0}$, $\braket{\bm{d}^{(1)}} = \bm{0}$, etc. Since $V$ and $\bm{d}$ are spatially periodic, Eqs.~\eqref{eq:SI_d_eqn_steady_4a}-\eqref{eq:SI_d_eqn_steady_4b} may be sequentially solved by means of Fourier series:
\begin{equation}
	\bm{d}^{(0)} (\bm{r})
	=
	-
	\frac{1}{kT}
	\sum_{\bm{q} \not= \bm{0}}
	\frac{\mathrm{i} \bm{q}}{q^{2}}
	V_{\bm{q}} \mathrm{e}^{\mathrm{i} \bm{q} \cdot \bm{r}}
	,
	\label{eq:SI_d0_soln}
\end{equation}
\begin{equation}
	\bm{d}^{(1)} (\bm{r})
	=
	\frac{1}{2(kT)^{2}}
	\sum_{\bm{q} \not= \bm{0}}
	\sum_{\bm{q}' \not= \bm{0}}
	\cdot
	\left( \frac{\mathrm{i} \bm{q}}{q^{2}} + \frac{2 \mathrm{i} \bm{q} \cdot (\bm{q} - \bm{q}')\bm{q}'}{q^{2} q'^{2}} \right)
	V_{\bm{q} - \bm{q}'} V_{\bm{q}'}
	\mathrm{e}^{\mathrm{i} \bm{q} \cdot \bm{r}}
	,
	\label{eq:SI_d1_soln}
\end{equation}
where $\bm{q}$ is the discrete wavevector and $V_{\bm{q}} \equiv L^{-2} \int_{L^{2}} [V(\bm{r}) - \braket{V}] \mathrm{e}^{- \mathrm{i} \bm{q} \cdot \bm{r}} \, \mathrm{d} \bm{r}$ denotes the Fourier integral of $V$.

By use of Eqs.~\eqref{eq:effective_diffusivity} and \eqref{eq:SI_d_field_perturbation_series}, the effective diffusivity of the Brownian particle is given by
\begin{flalign}
	\tens{D}
	&=
	\frac{kT}{\gamma} \tens{I}
	+
	\frac{1}{\gamma} \braket{\bm{d} \bm{\nabla}_{\bm{r}} V}
	\nonumber\\
	&=
	\frac{kT}{\gamma} \tens{I}
	+
	\frac{1}{\gamma} \braket{\bm{d}^{(0)} \bm{\nabla}_{\bm{r}} V}
	+
	\frac{1}{\gamma} \braket{\bm{d}^{(1)} \bm{\nabla}_{\bm{r}} V}
	+
	\cdots
	.
	\label{eq:SI_diffusivity_tensor_smallV_expansion}
\end{flalign}
Multiplying Eqs.~\eqref{eq:SI_d0_soln} by $\bm{\nabla}_{\bm{r}} V = \sum_{\bm{q}\not=\bm{0}} \mathrm{i} \bm{q} V_{\bm{q}} \mathrm{e}^{\mathrm{i} \bm{q} \cdot \bm{r}}$ and averaging over an $L \times L$ cell yields the force-displacement dyads,
\begin{equation}
	\braket{\bm{d}^{(0)} \bm{\nabla}_{\bm{r}} V}
	=
	-
	\frac{1}{kT}
	\sum_{\bm{q} \not= \bm{0}}
	\frac{\bm{q}\bm{q}}{q^{2}}
	|V_{\bm{q}}|^{2}
	,
	\label{eq:SI_d0_gradV_avg}
\end{equation}
\begin{equation}
	\braket{\bm{d}^{(1)} \bm{\nabla}_{\bm{r}} V}
	=
	\frac{1}{2(kT)^{2}}
	\sum_{\bm{q} \not= \bm{0}}
	\sum_{\bm{q}' \not= \bm{0}}
	\left( \frac{\bm{q}\bm{q}}{q^{2}} + \frac{2 \bm{q} \cdot (\bm{q} - \bm{q}')\bm{q}'\bm{q}}{q^{2} q'^{2}} \right)
	V_{\bm{q} - \bm{q}'} V_{\bm{q}'} V_{-\bm{q}}
	,
	\label{eq:SI_d1_gradV_avg}
\end{equation}
where $|V_{\bm{q}}|^{2} \equiv V_{\bm{q}} V_{-\bm{q}}$. Thus, the diffusivity tensor $\tens{D}$ admits the Fourier-series representation,
\begin{flalign}
	\tens{D}
	&=
	\frac{kT}{\gamma} \,
	\left[
	\tens{I}
	-
	\frac{1}{(kT)^{2}}\sum_{\bm{q} \not= \bm{0}}
	\frac{\bm{q}\bm{q}}{q^{2}}
	|V_{\bm{q}}|^{2}
	+
	\frac{1}{2(kT)^{3}}
	\sum_{\bm{q} \not= \bm{0}}
	\sum_{\bm{q}' \not= \bm{0}}
	\left( \frac{\bm{q}\bm{q}}{q^{2}} + \frac{2 \bm{q} \cdot (\bm{q} - \bm{q}')\bm{q}'\bm{q}}{q^{2} q'^{2}} \right)
	V_{\bm{q} - \bm{q}'} V_{\bm{q}'} V_{-\bm{q}}
	+
	\cdots
	\right]
	.
	\label{eq:SI_diffusivity_smallV_1}
\end{flalign}

An alternative expression for $\tens{D}$ can be obtained by writing leading-order displacement field as the negative gradient of a potential,
\begin{equation}
	\bm{d}^{(0)} (\bm{r})
	=
	- \frac{1}{kT} \bm{\nabla}_{\bm{r}} \varPhi(\bm{r})
	,
	\label{eq:SI_diffusivity_as_gradient_of_potential}
\end{equation}
where $\varPhi(\bm{r})$ satisfies the Poisson equation,
\begin{equation}
	\nabla^{2}_{\bm{r}} \varPhi(\bm{r}) = - [V(\bm{r}) - \braket{V}],
	\label{eq:SI_poisson_eqn}
\end{equation}
subject to the closure $\braket{\varPhi} = 0$. The Fourier-series solution of Eq.~\eqref{eq:SI_poisson_eqn} is
\begin{equation}
	\varPhi(\bm{r})
	=
	\sum_{\bm{q} \not= \bm{0}}
	q^{-2}
	V_{\bm{q}} \mathrm{e}^{\mathrm{i} \bm{q} \cdot \bm{r}}
	.
	\label{eq:SI_poisson_field_soln}
\end{equation}
By use of Eqs.~\eqref{eq:SI_d0_soln}, \eqref{eq:SI_d1_gradV_avg}, and the convolution theorem, it can be shown that
\begin{equation}
	\braket{\bm{d}^{(1)} \bm{\nabla}_{\bm{r}} V}
	=
	\frac{1}{2kT}
	\braket{(V - \braket{V})^{2} \bm{\nabla}_{\bm{r}} \bm{d}^{(0)}}
	+
	\braket{(\bm{\nabla}_{\bm{r}} \bm{d}^{(0)} \cdot \bm{\nabla}_{\bm{r}} V)\bm{d}^{(0)}}
	.
	\label{eq:SI_d1_gradV_avg_identity}
\end{equation}
Then, by Eqs.~\eqref{eq:SI_d0_gradV_avg}, \eqref{eq:SI_diffusivity_as_gradient_of_potential}, \eqref{eq:SI_poisson_field_soln}, and \eqref{eq:SI_d1_gradV_avg_identity}, it follows that
\begin{equation}
	\braket{\bm{d}^{(0)} \bm{\nabla}_{\bm{r}} V}
	=
	- \frac{1}{kT} \braket{\bm{\nabla}_{\bm{r}} \varPhi \bm{\nabla}_{\bm{r}} V}
	,
	\label{eq:SI_d0_gradV_avg_2}
\end{equation}
\begin{equation}
	\braket{\bm{d}^{(1)} \bm{\nabla}_{\bm{r}} V}
	=
	\frac{1}{2(kT)^{2}}
	\left[
	-
	\braket{(V - \braket{V})^{2} \bm{\nabla}_{\bm{r}} \bm{\nabla}_{\bm{r}} \varPhi}
	+
	2 \braket{(\bm{\nabla}_{\bm{r}} \bm{\nabla}_{\bm{r}} \varPhi \cdot \bm{\nabla}_{\bm{r}} V)\bm{\nabla}_{\bm{r}} \varPhi}
	\right]
	.
	\label{eq:SI_d1_gradV_avg_2}
\end{equation}
Substituting Eqs.~\eqref{eq:SI_d0_gradV_avg_2}-\eqref{eq:SI_d1_gradV_avg_2} into \eqref{eq:SI_diffusivity_tensor_smallV_expansion} then gives the alternative representation,
\begin{flalign}
	\tens{D}
	&=
	\frac{kT}{\gamma}
	\left(
	\tens{I}
	-
	\frac{1}{(kT)^{2}}
	\braket{\bm{\nabla}_{\bm{r}} \varPhi \bm{\nabla}_{\bm{r}} V}
	+
	\frac{1}{2(kT)^{3}}
	\left[
	-
	\braket{(V - \braket{V})^{2} \bm{\nabla}_{\bm{r}} \bm{\nabla}_{\bm{r}} \varPhi}
	+
	2 \braket{(\bm{\nabla}_{\bm{r}} \bm{\nabla}_{\bm{r}} \varPhi \cdot \bm{\nabla}_{\bm{r}} V)\bm{\nabla}_{\bm{r}} \varPhi}
	\right]
	+
	\cdots
	\right)
	.
	\label{eq:SI_diffusivity_smallV_2}
\end{flalign}

Since $V(\bm{r})$ is isotropic, only the trace of the steady diffusivity tensor need be computed: $D \equiv \tfrac{1}{2} \tens{D} : \tens{I}$. Using Eq.~\eqref{eq:SI_diffusivity_tensor_smallV_expansion}, the scalar diffusivity $D$ is given by
\begin{equation}
	D
	=
	\frac{kT}{\gamma}
	+
	\frac{1}{2\gamma} \braket{\bm{d}^{(0)} \cdot \bm{\nabla}_{\bm{r}} V}
	+
	\frac{1}{2\gamma} \braket{\bm{d}^{(1)} \cdot \bm{\nabla}_{\bm{r}} V}
	+
	\cdots
	.
	\label{eq:SI_diffusivity_scalar_smallV_expansion}
\end{equation}
Taking the trace of Eqs.~\eqref{eq:SI_d0_gradV_avg_2}-\eqref{eq:SI_d1_gradV_avg_2}, integrating by parts, and applying Eq.~\eqref{eq:SI_poisson_eqn} then gives
\begin{flalign}
	\braket{\bm{d}^{(0)} \cdot \bm{\nabla}_{\bm{r}} V}
	&=
	- \frac{1}{kT} \braket{\bm{\nabla}_{\bm{r}} \varPhi \cdot \bm{\nabla}_{\bm{r}} V}
	\nonumber\\
	&=
	-\frac{1}{kT} \braket{(V - \braket{V})^{2}}
	,
	\label{eq:SI_d0_gradV_avg_3}
\end{flalign}
\begin{flalign}
	\braket{\bm{d}^{(1)} \cdot \bm{\nabla}_{\bm{r}} V}
	&=
	\frac{1}{2(kT)^{2}}
	\left[
	-
	\braket{(V - \braket{V})^{2} \nabla^{2}_{\bm{r}} \varPhi}
	+
	2 \braket{(\bm{\nabla}_{\bm{r}} \bm{\nabla}_{\bm{r}} \varPhi \cdot \bm{\nabla}_{\bm{r}} V) \cdot\bm{\nabla}_{\bm{r}} \varPhi}
	\right]
	\nonumber\\
	&=
	\frac{1}{2(kT)^{2}}
	\left[
	\braket{(V - \braket{V})^{3}}
	+
	\braket{\bm{\nabla}_{\bm{r}} (|\bm{\nabla}_{\bm{r}} \varPhi|^{2}) \cdot \bm{\nabla}_{\bm{r}} V}
	\right]
	.
	\label{eq:SI_d1_gradV_avg_3}
\end{flalign}
Inserting Eqs.~\eqref{eq:SI_d0_gradV_avg_3}-\eqref{eq:SI_d1_gradV_avg_3} into \eqref{eq:SI_diffusivity_scalar_smallV_expansion} then gives
\begin{equation}
	D
	=
	\frac{kT}{\gamma}
	\left(
	1
	-
	\frac{1}{2(kT)^{2}}
	\braket{(V - \braket{V})^{2}}
	+
	\frac{1}{4(kT)^{3}}
	\left[
	\braket{(V - \braket{V})^{3}}
	+
	\braket{\bm{\nabla}_{\bm{r}} (|\bm{\nabla}_{\bm{r}} \varPhi|^{2}) \cdot \bm{\nabla}_{\bm{r}} V}
	\right]
	+
	\cdots
	\right)
	.
	\label{eq:SI_diffusivity_scalar_smallV_expansion}
\end{equation}
The last expression is exactly Eq.~\eqref{eq:diffusivity_soft_traps} from the main text.

\subsection{Derivation of Eq. (9): stationary traps with deep potential wells}

For stationary, ``deep'' potential wells, $\Delta V \gg kT$, the small-potential perturbation series \eqref{eq:SI_d_field_perturbation_series} fails to converge. Unfortunately, no exact analytical solution of Eq.~\eqref{eq:SI_d_eqn_steady_2} is readily available. However, one can take advantage of the fact that, for deep potential wells, the probability density is strongly localized near the origin $\bm{r} = \bm{0}$ of the lattice cell where the potential-energy field $V(\bm{r})$ is minimized. Then, a useful {\it approximation} of the $\bm{d}$-field is
\begin{flalign}
	\bm{d} (\bm{r})
	&\approx
	- \bm{r} g (\bm{r})
	\nonumber\\
	&=
	- \frac{\bm{r} \mathrm{e}^{-V(\bm{r})/kT}}{\braket{\mathrm{e}^{-V/kT}}}
	.
	\label{eq:SI_d_field_largeV_approx}
\end{flalign}
Eq.~\eqref{eq:SI_d_field_largeV_approx} is the particular solution of Eq.~\eqref{eq:SI_d_eqn_steady_2} and conserves probability, $\braket{\bm{d}} = \bm{0}$. However, this particular solution clearly violates the periodic boundary conditions at the edges of the lattice cell $x = \pm L/2$, $y=\pm L/2$, incurring an error of $O(L\mathrm{e}^{-\Delta V/kT}/\braket{\mathrm{e}^{-V/kT}})$ that decreases in magnitude with increasing trap stiffness. Fig.~\ref{Fig:SI_StrongTrap_dx_vs_x} compares the approximation, Eq.~\eqref{eq:SI_d_field_largeV_approx}, against the ``exact'' numerical solution for the displacement field, showing very good agreement. The slight error in the approximation is due to the neglect of the homogeneous solution of Eq.~\eqref{eq:SI_d_eqn_steady_2}, which is complicated by the 2D potential-energy field given by Eq.~\eqref{eq:harmonic_potential}. It will be shown that the error in this approximation for the $\bm{d}$-field quantitatively (though not qualitatively) impacts the prediction for the effective diffusivity.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=1\linewidth]{Figs/SI/SIFig_StrongPotentialLimit_Dx_vs_X.png}
		\caption{
		Comparison of numerical solution for the steady displacement field density $d_{x}(x,y)$ against the particular solution [see Eq.~\eqref{eq:SI_d_field_largeV_approx}] for a stiff trap, $\kappa=5$ $kT/\mu$m$^2$. (a) 2D contour plot of $d_{x}$ with line traces at four distinct values of $y$. (b) Plot of $d_{x}$ against $x$ for each line trace shows favorable agreement to Eq.~\eqref{eq:SI_d_field_largeV_approx}.
		}
		\label{Fig:SI_StrongTrap_dx_vs_x}
	\end{center}
	\vspace{-18pt}
\end{figure}

Using Eq.~\eqref{eq:harmonic_potential} for $V(\bm{r})$ and Eq.~\eqref{eq:SI_d_field_largeV_approx} for $\bm{d}(\bm{r})$, the force-displacement dyad that appears in Eq.~\eqref{eq:diffusivity_soft_traps} can now be approximated as
\begin{flalign}
	\braket{\bm{d} \bm{\nabla}_{\bm{r}} V}
	&\approx
	-
	\frac{1}{\kappa}
	\frac{\braket{\mathrm{e}^{-V/kT}\bm{\nabla}_{\bm{r}} V \bm{\nabla}_{\bm{r}} V}}{\braket{\mathrm{e}^{-V/kT}}}
	,
	\label{eq:SI_d_gradV_avg_largeV}
\end{flalign}
where we've used the fact that $\bm{\nabla}_{\bm{r}} V = \kappa \bm{r}$ for $r \le \frac{1}{2} W_{\text{trap}}$ and $=\bm{0}$ otherwise. Defining the well depth as $\Delta V =\tfrac{1}{8} \kappa W_{\text{trap}}^{2}$, the cell averages in Eq.~\eqref{eq:SI_d_gradV_avg_largeV} become
\begin{flalign}
	\braket{\mathrm{e}^{-V/kT}}
	&=
	\frac{2\uppi kT}{\kappa L^{2}}\left( 1 - \mathrm{e}^{-\Delta V/kT} \right)
	+
	\bigg( 1- \frac{2\uppi \Delta V}{\kappa L^{2}} \bigg) \mathrm{e}^{-\Delta V/kT}
	,
\end{flalign}
\begin{flalign}
	\braket{\mathrm{e}^{-V/kT}\bm{\nabla}_{\bm{r}} V \bm{\nabla}_{\bm{r}} V}
	&=
	\frac{2\uppi (kT)^{2}}{L^{2}}
	\left[
		1
		-
		\left( 1 + \frac{\Delta V}{kT} \right)\mathrm{e}^{-\Delta V/kT}
	\right]
	\tens{I}
	.
\end{flalign}
Substitution into Eq.~\eqref{eq:SI_d_gradV_avg_largeV} then gives, upon simplification,
\begin{flalign}
	\braket{\bm{d} \bm{\nabla}_{\bm{r}} V}
	&\approx
	kT
	\left\{
		-1
		+
		\left[
		1
		+
		\frac{2\uppi kT}{\kappa L^{2}}
		\left(
			\mathrm{e}^{\Delta V/kT}
			-
			\frac{\Delta V}{kT}
			-
			1
		\right)
		\right]^{-1}
	\right\}
	\tens{I}
	\nonumber\\
	&\approx
	\left(
	- k T
	+
	\frac{\kappa L^{2}}{2\uppi}
	\mathrm{e}^{-\Delta V / kT}
	\right)
	\tens{I}
	\quad
	\text{for}
	\quad
	\Delta V \gg kT.
\end{flalign}
Substitution into Eq.~\eqref{eq:effective_diffusivity} and replacing $\kappa \tens{I}$ by $(\bm{\nabla}_{\bm{r}}\bm{\nabla}_{\bm{r}}V)|_{\bm{r}=\bm{0}}$ then gives the following approximation for the diffusivity tensor:
\begin{equation}
	\tens{D}
	\approx
	\frac{L^{2}}{2\uppi \gamma}
	\mathrm{e}^{-\Delta V / kT}
	(\bm{\nabla}_{\bm{r}}\bm{\nabla}_{\bm{r}}V)|_{\bm{r}=\bm{0}}
	,
\end{equation}
or, upon taking one-half the trace,
\begin{equation}
	D
	\approx
	\frac{L^{2}}{4\uppi \gamma}
	\mathrm{e}^{-\Delta V / kT}
	(\nabla_{\bm{r}}^{2}V)|_{\bm{r}=\bm{0}}
	.
	\label{eq:SI_diffusivity_largeV_approx}
\end{equation}
This is exactly the form that would be predicted by Kramers' theory for the escape of a Brownian particle from a deep potential well \cite{kramers1940brownian,brinkman1956brownian,brinkman1956brownian2}.
Comparison of Eq.~\eqref{eq:SI_diffusivity_largeV_approx} to numerical calculations of $D$ indicates the qualitatively correct dependence on the trapping strength, but quantitative discrepancies due to errors in the approximation \eqref{eq:SI_d_field_largeV_approx} for the $\bm{d}$-field (see Fig.~\ref{Fig:SI_StrongTrap_D_vs_K_LogLinearPlot}). Quantitative agreement can be obtained by renormalizing the above result by a factor that depends upon the ratio $W_{\text{trap}} / L$. Therefore, we write
\begin{equation}
	D
	\propto
	\frac{L^{2}}{4\uppi \gamma}
	\mathrm{e}^{-\Delta V / kT}
	(\nabla_{\bm{r}}^{2}V)|_{\bm{r}=\bm{0}}
	\label{eq:SI_diffusivity_largeV_approx_proportionality}
\end{equation}
up to a proportionality constant. Eq.~\eqref{eq:SI_diffusivity_largeV_approx_proportionality} is identical to Eq.~\eqref{eq:diffusivity_stiff_traps} from the main text.
For traps of diameter $W_{\text{trap}} = 3.2$ $\mu$m spaced a distance $L=6$ $\mu$m apart, a proportionality constant of 1.5 gives quantitative agreement with the exact dispersion theory (see Fig.~\ref{Fig:SI_StrongTrap_D_vs_K_LogLinearPlot}).

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.5\linewidth]{Figs/SI/SIFig_StrongPotentialLimit_ExponentialTrends.png}
		\caption{%
		Log-linear plot of diffusivity $D$ against trap stiffness $\kappa$. The full numerical solution (solid curve) is compared against Eq.~\eqref{eq:diffusivity_stiff_traps} (dashed curves) using two different constants of proportionality. Irrespective of the numerical prefactor, Eq.~\eqref{eq:diffusivity_stiff_traps} demonstrates the appropriate scaling with the trapping strength and is consistent with Kramers' theory of activated escape. A proportionality constant of 1.5 gives quantitative agreement with the exact solution for the specific geometry considered in this study.
		}
		\label{Fig:SI_StrongTrap_D_vs_K_LogLinearPlot}
	\end{center}
	\vspace{-18pt}
\end{figure}

\subsection{Derivation of Eq. (11): oscillating traps in the high-frequency limit}

In the high-frequency limit, the potential-energy field is cycled in the $x$-direction at a rate much faster than the response time of the Brownian particle. A reasonable model for this system is a quasi-steady, uniform convection in the $x$-direction, for which we make the ansatz $g = g(y)$ and $d_{y} = d_{y} (y)$ (for the time being, we will ignore the $d_{x}$-field). Eqs.~\eqref{eq:g_eqn} and \eqref{eq:d_eqn} then simplify to
\begin{equation}
	kT
	\frac{\mathrm{d}^{2}g}{\mathrm{d} y^{2}}
	+
	\frac{\partial V}{\partial y} \frac{\mathrm{d} g}{\mathrm{d} y}
	+
	\left( \frac{\partial^{2} V}{\partial x^{2}} + \frac{\partial^{2} V}{\partial y^{2}} \right) g
	=
	0
	,
	\label{eq:SI_g_eqn_1d}
\end{equation}
\begin{equation}
	kT
	\frac{\mathrm{d}^{2} d_{y}}{\mathrm{d} y^{2}}
	+
	\frac{\partial V}{\partial y} \frac{\mathrm{d} d_{y}}{\mathrm{d} y}
	+
	\left( \frac{\partial^{2} V}{\partial x^{2}} + \frac{\partial^{2} V}{\partial y^{2}} \right) d_{y}
	=
	-
	2 kT \frac{\mathrm{d} g}{\mathrm{d} y}
	-
	g \frac{\partial V}{\partial y}
	+
	\bigg\langle g \frac{\partial V}{\partial y} \bigg\rangle g
	.
	\label{eq:SI_dy_eqn_1d}
\end{equation}
Averaging Eqs.~\eqref{eq:SI_g_eqn_1d}-\eqref{eq:SI_dy_eqn_1d} over the $x$-direction only and defining the modified potential,
\begin{equation}
	v(y)
	=
	\frac{1}{L} \int_{-L/2}^{L/2} V(x,y) \, \mathrm{d} x
	,
\end{equation}
then gives
\begin{equation}
	kT
	\frac{\mathrm{d}^{2}g}{\mathrm{d} y^{2}}
	+
	\frac{\mathrm{d}}{\mathrm{d} y} \left( g\frac{\partial v}{\partial y} \right)
	=
	0
	,
	\label{eq:SI_g_eqn_1d_2}
\end{equation}
\begin{equation}
	kT
	\frac{\mathrm{d}^{2} d_{y}}{\mathrm{d} y^{2}}
	+
	\frac{\mathrm{d}}{\mathrm{d} y} \left( d_{y} \frac{\partial v}{\partial y} \right)
	=
	-
	2 kT \frac{\mathrm{d} g}{\mathrm{d} y}
	-
	g \frac{\mathrm{d} v}{\mathrm{d} y}
	+
	\bigg\langle g \frac{\mathrm{d} v}{\mathrm{d} y} \bigg\rangle g
	,
	\label{eq:SI_dy_eqn_1d_2}
\end{equation}
where we have applied the conditions $V(L/2,y)=V(-L/2,y)$ and $(\partial V / \partial x)|_{x = \pm L/2} = 0$. Here, it is understood that the cell average of a one-dimensional (1D) function $f(y)$ simplifies to a 1D average in the $y$-direction, $\braket{f} = L^{-1} \int_{-L/2}^{L/2} f(y)\, \mathrm{d}y$.

Eqs.~\eqref{eq:SI_g_eqn_1d_2}-\eqref{eq:SI_dy_eqn_1d_2} are the 1D versions of Eqs.~\eqref{eq:SI_g_eqn_steady}-\eqref{eq:SI_d_eqn_steady}. The solution of Eq.~\eqref{eq:SI_g_eqn_1d_2} for the $g$-field, subject to the constraint $\braket{g} = 1$, is the 1D analog of Eq.~\eqref{eq:SI_g_field_steady}:
\begin{equation}
	g(y)
	=
	\frac{\mathrm{e}^{-v(y)/kT}}{\braket{\mathrm{e}^{-v/kT}}}
	.
\end{equation}
Eq.~\eqref{eq:SI_dy_eqn_1d_2} then simplifies to
\begin{flalign}
	kT
	\frac{\mathrm{d}^{2} d_{y}}{\mathrm{d} y^{2}}
	+
	\frac{\mathrm{d}}{\mathrm{d} y} \left( d_{y} \frac{\partial v}{\partial y} \right)
	&=
	-
	kT \frac{\mathrm{d} g}{\mathrm{d} y}
	\nonumber \\
	&=
	\frac{\mathrm{e}^{-v/kT}}{\braket{\mathrm{e}^{-v/kT}}}
	\frac{\mathrm{d} v}{\mathrm{d} y}
	,
	\label{eq:SI_dy_eqn_1d_3}
\end{flalign}
which is the 1D analog of Eq.~\eqref{eq:SI_d_eqn_steady_2}. Unlike the 2D problem, the 1D problem admits an exact analytical solution:
\begin{flalign}
	d_{y} (y)
	&=
	- y g (y)
	+
	c_{1} \mathrm{e}^{- v(y)/kT}
	\int_{0}^{y}
	\mathrm{e}^{v(\eta)/kT} \, \mathrm{d} \eta
	+
	c_{2} L \mathrm{e}^{-v(y) / kT}
	\nonumber\\
	&=
	- \frac{y \mathrm{e}^{-v(y)/kT}}{\braket{\mathrm{e}^{-v/kT}}}
	+
	c_{1} \mathrm{e}^{- v(y)/kT}
	\int_{0}^{y}
	\mathrm{e}^{v(\eta)/kT} \, \mathrm{d} \eta
	+
	c_{2} L \mathrm{e}^{-v(y) / kT}
	.
	\label{eq:SI_dy_field_1d_general_soln}
\end{flalign}
The first term on the right-hand side of Eq.~\eqref{eq:SI_dy_field_1d_general_soln} is simply the particular solution of Eq.~\eqref{eq:SI_dy_eqn_1d_3}; it is the 1D analog of Eq.~\eqref{eq:SI_d_field_largeV_approx}, which was used to approximate the full solution in the strong-potential limit. The remaining terms in Eq.~\eqref{eq:SI_dy_field_1d_general_soln} are the homogeneous solutions, with constants $c_{1}$, $c_{2}$ that must be determined from the periodicity and normalization conditions,
\begin{subequations}
\begin{gather}
	d_{y} (L/2) - d_{y} (-L/2) = 0,
	\label{eq:SI_dy_1d_bc1}
	\\
	\braket{d_{y}}
	=
	\frac{1}{L} \int_{-L/2}^{L/2} d_{y}(y) \, \mathrm{d}y
	=
	0.
	\label{eq:SI_dy_1d_bc2}
\end{gather}
	\label{eq:SI_dy_1d_bcs}%
\end{subequations}
Inserting Eq.~\eqref{eq:SI_dy_field_1d_general_soln} into \eqref{eq:SI_dy_1d_bcs}, setting $v(L/2) = v(-L/2)$, and solving for the two unknowns $c_{1}$ and $c_{2}$ gives
\begin{subequations}
\begin{flalign}
	c_{1} &= \braket{\mathrm{e}^{-v/kT}}^{-1} \braket{\mathrm{e}^{v/kT}}^{-1}
	,
	\\
	c_{2} &=
	\frac{1}{L} \braket{\mathrm{e}^{-v/kT}}^{-2}
	\left(
		\braket{y\mathrm{e}^{-v/kT}}
		-
		\braket{\mathrm{e}^{v/kT}}^{-1}
		\bigg\langle \mathrm{e}^{-v/kT}\int_{0}^{y} \mathrm{e}^{v(\eta)/kT} \mathrm{d} \eta \bigg\rangle
	\right)
	.
\end{flalign}
	\label{eq:SI_dy_1d_constants}%
\end{subequations}

With the solution for $d_{y}(y)$ fully specified, it remains to compute the effective diffusivity along the $y$-axis. Multiplying Eq.~\eqref{eq:SI_dy_field_1d_general_soln} by $\mathrm{d} v / \mathrm{d}y$, applying the inverse chain rule, and averaging over the $y$-direction gives
\begin{flalign}
	\bigg\langle d_{y} \frac{\mathrm{d} v}{\mathrm{d}y}\bigg\rangle
	&=
	kT
	\left(
	\braket{\mathrm{e}^{-v/kT}}^{-1}
	\bigg\langle y \frac{\mathrm{d}\mathrm{e}^{-v/kT}}{\mathrm{d}y} \bigg\rangle
	-
	c_{1}
	\bigg\langle \frac{\mathrm{d}\mathrm{e}^{-v/kT}}{\mathrm{d}y}
	\int_{0}^{y}
	\mathrm{e}^{v(\eta)/kT} \, \mathrm{d} \eta
	\bigg\rangle
	-
	c_{2} L
	\bigg\langle \frac{\mathrm{d}\mathrm{e}^{-v/kT}}{\mathrm{d}y} \bigg\rangle
	\right).
	\label{eq:SI_dy_gradv_1d_avg}
\end{flalign}
Inserting Eqs.~\eqref{eq:SI_dy_1d_constants} into \eqref{eq:SI_dy_gradv_1d_avg} and integrating by parts then gives, after some simplification,
\begin{flalign}
	\bigg\langle d_{y} \frac{\mathrm{d} v}{\mathrm{d}y}\bigg\rangle
	&=
	k T
	\left(
		-1
		+
		\braket{\mathrm{e}^{-v/kT}}^{-1} \braket{\mathrm{e}^{v/kT}}^{-1}
	\right)
	.
\end{flalign}
Since $d_{y}$ is independent of $x$, $\braket{d_{y} (\partial V / \partial y)} = \braket{d_{y} (\mathrm{d} v / \mathrm{d} y)}$. Thus, the $yy$ component of Eq.~\eqref{eq:effective_diffusivity} simplifies to
\begin{flalign}
	\overline{D}_{yy}
	&=
	\frac{kT}{\gamma}
	+
	\frac{1}{\gamma}
	\bigg\langle d_{y} \frac{\mathrm{d} v}{\mathrm{d}y}\bigg\rangle
	\nonumber\\
	&=
	\frac{kT}{\gamma}
	\braket{\mathrm{e}^{-v/kT}}^{-1} \braket{\mathrm{e}^{v/kT}}^{-1}
	,
	\label{eq:SI_Dyy_1d}
\end{flalign}
where an overbar is used to denote the long-time average over one periodic cycle.
This is the classical result for diffusion of a Brownian particle in a 1D periodic potential \cite{lifson1962self,Festa1978}.

Up until now, we have neglected the $d_{x}$-field, which appears in the $xx$-component of Eq.~\eqref{eq:effective_diffusivity} and, therefore, influences the effective diffusivity along the $x$-axis. To a first approximation, we assume that the gradients in the $x$-direction have been ``smeared out'' so that dispersion in that direction is negligible: $\braket{d_{x} (\partial V / \partial x)} \approx 0$. This approximation is consistent with a model of dispersion in an effectively 1D potential. Therefore, the $xx$-component of Eq.~\eqref{eq:effective_diffusivity} (time-averaged) is simply the Stokes-Einstein-Sutherland diffusivity:
\begin{equation}
	\overline{D}_{xx} = \frac{kT}{\gamma}
	.
	\label{eq:SI_Dxx_1d}
\end{equation}
Eqs.~\eqref{eq:SI_Dyy_1d} and \eqref{eq:SI_Dxx_1d} are exactly the same as Eq.~\eqref{eq:diffusivity_high_frequency} from the main text.

\section{5. Brownian Dynamics Simulations} \label{sec:brownian_dynamics_simulations}

The Langevin equation of motion corresponding to Eqs.~\eqref{eq:SI_smoluchowski_eqn}-\eqref{eq:SI_coordinate_conversion} is given by
\begin{equation}
	\frac{\mathrm{d} \bm{r}_{i} (t)}{\mathrm{d} t}
	=
	-\bm{u}(t)
	-
	\frac{1}{\gamma} \bm{\nabla}_{\bm{r}} V[\bm{r}_{i}(t)]
	+
	\sqrt{\frac{2kT}{\gamma}}\bm{B}_{i}(t)
	,
	\qquad
	i = 1, 2, \dots, N_{\text{p}}
	,
	\label{eq:SI_langevin_eqn}
\end{equation}
where $i$ is the particle index, $N_{\text{p}}$ is the total number of particles in the system, and $\bm{B}_{i} (t)$ is a white-noise source with statistics,
\begin{equation}
	\braket{\bm{B}_{i}(t)}
	=
	\bm{0},
	\qquad
	\braket{\bm{B}_{i}(t)\bm{B}_{i}(t')}
	=
	\delta(t - t')
	\tens{I}
	.
	\label{eq:SI_fluctuation_dissipation_theorem}
\end{equation}
[Note that the angle brackets $\langle\,\cdot\,\rangle$ appearing in Eq.~\eqref{eq:SI_fluctuation_dissipation_theorem} denote {\it ensemble} averages and are not to be confused with the {\it cell} average defined in the main text.]
The potential-energy field $V(\bm{r})$ and convective velocity $\bm{u}(t)$ appearing in Eq.~\eqref{eq:SI_langevin_eqn} are given by Eqs.~\eqref{eq:harmonic_potential} and \eqref{eq:trap_velocity}, respectively. Interactions between particles have been neglected, so the $N_{\text{p}}$ equations of motion are uncoupled.
For the purpose of numerically time-advancing Eq.~\eqref{eq:SI_langevin_eqn}, it is convenient to shift to the laboratory frame in which the position of each particle is measured as $\bar{\bm{r}}_{i}(t) = \bm{r}_{0}(t) + \bm{r}_{i}(t)$, where $\bm{r}_{0}(t) = \int_{0}^{t} \bm{u} (\tau) \, \mathrm{d} \tau = \hat{\bm{e}}_{x} A \sin{(\omega t)}$ denotes the time-dependent position of the moving traps. In this frame, Eq.~\eqref{eq:SI_langevin_eqn} becomes
\begin{equation}
	\frac{\mathrm{d} \bar{\bm{r}}_{i}(t)}{\mathrm{d} t}
	=
	-
	\frac{1}{\gamma} \bm{\nabla}_{\bar{\bm{r}}} V[\bar{\bm{r}}_{i}(t)-\bm{r}_{0}(t)]
	+
	\sqrt{\frac{2kT}{\gamma}}\bm{B}_{i}(t)
	,
	\qquad
	i = 1, 2, \dots, N_{\text{p}}
	.
	\label{eq:SI_langevin_eqn_relative}
\end{equation}
Here, the convective term has been eliminated and the potential-energy field oscillates in time.

In our Brownian dynamics simulations, we numerically advanced Eq.~\eqref{eq:SI_langevin_eqn_relative} using the GPU-enabled HOOMD-blue software package \cite{anderson2020hoomd}. A system of $N_{\text{p}} = 10,000$ particles was initialized at random positions within a periodically replicated $L\times L$ cell and advanced for $\tau = 10,000$ s (2.78 h) using a time step $\Delta t = 1$ ms. Fig.~\ref{Fig:SI_Gfield_Theory_vs_Simulation} shows that the simulated probability density shows excellent agreement with the deterministic solution of the corresponding Smoluchowski equation [Eq.~\eqref{eq:g_eqn}]. The MSD and effective diffusivity of the particles were then computed exactly as in the experiments using Eqs.~\eqref{eq:SI_msd_experiments}-\eqref{eq:SI_diffusivity_experiments}, wherein the time integral was discretized using the left Riemann sum.

\begin{figure}[!h]
	\begin{center}
\includegraphics[width=1\linewidth]{Figs/SI/SIFig_Gfield_TheoryVsSimulation.png}
		\caption{%
		Comparison of the convected probability density $g(x,y,t)$ for a stiff trap near the critical frequency ($\kappa=5$ $kT/\mu$m$^2$, $\omega/2\uppi=18.33$ mHz) from
		(a) deterministic solution of the Smoluchowski equation [Eq.~\eqref{eq:g_eqn}] and (b) stochastic simulation of the Langevin equation [Eq.~\eqref{eq:SI_langevin_eqn}].
		}
		\label{Fig:SI_Gfield_Theory_vs_Simulation}
	\end{center}
	\vspace{-18pt}
\end{figure}

\subsection{Derivation of Eq. (12): convective escape of a Brownian particle from a harmonic well}

We wish to estimate the critical oscillation frequency $\omega_{\text{max}}$ at which a Brownian particle rattling around the bottom of a potential-energy well is convected near the edge of the well with ample probability for escape. To make such an estimate, we start with the Langevin equation, Eq.~\eqref{eq:SI_langevin_eqn}, simplified for a single particle in a harmonic well $V(\bm{r}) = \tfrac{1}{2} \kappa r^{2}$:
\begin{equation}
	\frac{\mathrm{d} \bm{r} (t)}{\mathrm{d} t}
	=
	-
	\frac{\kappa}{\gamma} \bm{r}(t)
	-
	\bm{u}(t)
	+
	\sqrt{\frac{2kT}{\gamma}}\bm{B}(t)
	.
	\label{eq:SI_langevin_eqn_harmonic_well}
\end{equation}
Eq.~\eqref{eq:SI_langevin_eqn_harmonic_well} may be straightforwardly integrated with the initial condition $\bm{r}(0) = \bm{0}$ to give the fluctuating particle position,
\begin{flalign}
	\bm{r} (t)
	&=
	\mathrm{e}^{-\kappa t / \gamma}
	\int_{0}^{t}
	\mathrm{e}^{\kappa s / \gamma}
	\left(
		-\bm{u}(s)
		+
		\sqrt{\frac{2kT}{\gamma}} \bm{B}(s)
	\right)
	\mathrm{d} s
	.
	\label{eq:SI_langevin_eqn_harmonic_well_soln}
\end{flalign}
Substituting Eq.~\eqref{eq:trap_velocity} into \eqref{eq:SI_langevin_eqn_harmonic_well_soln} for the convective velocity then gives, upon integration,
\begin{flalign}
	\bm{r} (t)
	&=
	-\hat{\bm{e}}_{x} A
	\left( \frac{\gamma \omega /\kappa}{1 + (\gamma \omega / \kappa)^{2}} \right)
	\left(
		\cos{(\omega t)}
		+
		\frac{\gamma \omega}{\kappa}
		\sin{(\omega t)}
		-
		\mathrm{e}^{-\kappa t / \gamma}
	\right)
	+
	\sqrt{\frac{2kT}{\gamma}}
	\mathrm{e}^{-\kappa t / \gamma}
	\int_{0}^{t}
	\mathrm{e}^{\kappa s / \gamma}
	\bm{B}(s)
	\, \mathrm{d} s
	.
	\label{eq:SI_langevin_eqn_harmonic_well_soln_2}
\end{flalign}
The first term on the right-hand side of Eq.~\eqref{eq:SI_langevin_eqn_harmonic_well_soln_2} is the deterministic part of the fluctuating particle particle position, which is driven by oscillatory convection and attenuated by the trapping force. The second term is the stochastic part due to Brownian motion. The mean displacement and mean squared displacement of the particle respectively capture strength of these deterministic and stochastic elements:
\begin{flalign}
	\braket{\bm{r}(t)}
	=
	-\hat{\bm{e}}_{x} A
	\left( \frac{\gamma \omega /\kappa}{1 + (\gamma \omega / \kappa)^{2}} \right)
	\left(
		\cos{(\omega t)}
		+
		\frac{\gamma \omega}{\kappa}
		\sin{(\omega t)}
		-
		\mathrm{e}^{-\kappa t / \gamma}
	\right)
	,
	\label{eq:SI_drift_in_harmonic_well}
\end{flalign}
\begin{flalign}
	\braket{(\bm{r}(t) - \braket{\bm{r}(t)})(\bm{r}(t) - \braket{\bm{r}(t)})}
	=
	\frac{kT}{\kappa} \left( 1 - \mathrm{e}^{-2\kappa t / \gamma} \right) \tens{I}
	,
	\label{eq:SI_variance_in_harmonic_well}
\end{flalign}
where we have applied the white-noise statistics, Eq.~\eqref{eq:SI_fluctuation_dissipation_theorem}, of the fluctuating $\bm{B}$-field.

After waiting a long enough time $t \gg \gamma /\kappa$, the exponential terms in Eqs.~\eqref{eq:SI_drift_in_harmonic_well}-\eqref{eq:SI_variance_in_harmonic_well} die off and we are left with an oscillating particle probability with variance $k T / \kappa$ given by Eq.~\eqref{eq:SI_variance_in_harmonic_well}. The amplitude of these oscillations are found from the extrema of the particle drift, Eq.~\eqref{eq:SI_drift_in_harmonic_well}:
\begin{equation}
	\sup_{t\ge 0}
	| \langle \bm{r}(t)\rangle |
	=
	\frac{\gamma \omega A/\kappa}{\sqrt{1+ (\gamma \omega / \kappa)^{2}}}
	\approx
	\frac{\gamma \omega A}{\kappa}
	\quad
	\text{for}
	\quad
	\frac{\gamma \omega}{\kappa} \ll 1
	.
\end{equation}
Thus, the basin of probability of size $\sim \sqrt{kT/\kappa}$ oscillates with amplitude $\sim \gamma \omega A / \kappa$ about the center of the potential-energy well. As the frequency $\omega$ is increased, the oscillations become more pronounced. The particle is expected to escape a well of finite width $W_{\text{trap}}$ when the spatial extent of the particle probability density crosses the edge of the well, at a critical frequency $\omega_{\text{max}}$:
\begin{equation}
	\tfrac{1}{2} W_{\text{trap}}
	\approx
	\frac{\gamma \omega_{\text{max}} A}{\kappa}
	+
	\sqrt{\frac{kT}{\kappa}}
	,
\end{equation}
or, solving for $\omega_{\text{max}}$,
\begin{equation}
	\omega_{\text{max}}
	\approx
	\frac{\kappa}{\gamma A}
	\left(
		\tfrac{1}{2} W_{\text{trap}}
		-
		\sqrt{\frac{kT}{\kappa}}
	\right)
	.
\end{equation}
The last expression is exactly Eq.~\eqref{eq:critical_frequency} from the main text.

\section{6. Additional Data}

In addition to measuring the effective diffusivity $\overline{\tens{\bm{D}}}$ as a function of the oscillation frequency $\omega$, we also varied the amplitude $A$ while holding the frequency fixed.
The strength of the convective velocity $\bm{u}(t) = \hat{\bm{e}}_{x} \omega A \cos{(\omega t)}$ may be modified by varying either the amplitude $A$ or the frequency $\omega$.
Fig.~\ref{Fig:SI4} plots $\overline{D}_{xx}$ and $\overline{D}_{yy}$ against $A$ for a fixed trap stiffness $\kappa = 5~kT / \mu\mathrm{m}^2$ and frequency $\omega/2\uppi = 18.3$ mHz.
This frequency corresponds to the critical frequency $\omega_{\text{max}}$ (for which $\overline{D}_{xx}$ is maximized) for $\kappa = 5~kT / \mu\mathrm{m}^2$ and $A = 5$ $\mu$m, as shown in the main text (see Fig.~\ref{fig:Fig3}).
We find that the $\overline{D}_{xx}$ is non-monotonic and achieves a maximum at $A = 5$ $\mu$m.
For amplitudes $A > 5$ $\mu$m, the convective motion is fast compared to the particle response time. Consequently, the particles sample regions outside of the harmonic well and their average diffusivity along the convection axis is reduced.

\begin{figure}[!h]
	\begin{center}
		\includegraphics[width=0.5\linewidth]{Figs/SI/SiFig_AmplitudeSweep.png}
		\caption{%
			Effective diffusivity as a function of oscillation amplitude at a fixed trap stiffness $\kappa = 5~kT / \mu\mathrm{m}^2$ and frequency $\omega/2\uppi = 18.3~\mathrm{mHz}$.
			Like Fig.~\ref{fig:Fig3} in the main text, $\overline{D}_{xx}$ is non-monotonic and reaches a maximum when the convection strength balances the harmonic trap strength.
			At very large amplitudes, the particle cannot quickly respond to the rapidly oscillating trap and explores regions outside of the harmonic well.
		}
		\label{Fig:SI4}
	\end{center}
	\vspace{-18pt}
\end{figure}

\pagebreak
\section{7. Supplemental Movies}
Below, we describe the Supplemental Movies associated with this manuscript. All time stamps corresponds to hours:minutes:seconds.

\begin{enumerate}
	\item[] \textbf{S1.} Experimental micrographs of silica particles with radius $a = 1.25~\mu$m diffusing through a stationary array of harmonic traps (6$\times$6 grid shown) with varying trap stiffness.

	\item[] \textbf{S2.} Microscopic Brownian dynamics simulations of a small sample of particles diffusing through a stationary array of harmonic traps (6$\times$6 grid shown) with varying trap stiffness (same parameters as in S1).

	\item[] \textbf{S3.} Experimental micrographs of silica particles with radius $a = 1.25~\mu$m diffusing through an oscillating array of stiff traps (6$\times$6 grid shown) with varying oscillation frequency and fixed trap stiffness $\kappa = 5~kT/\mu\mathrm{m}^2$. The second part of the movie shows the trajectories of several tagged particles.

	\item[] \textbf{S4.} Microscopic Brownian dynamics simulations of a small sample of particles diffusing through an oscillating array of stiff traps (6$\times$6 grid shown) with varying oscillation frequency and fixed trap stiffness $\kappa = 5~kT/\mu\mathrm{m}^2$ (same parameters as in S3).

	\item[] \textbf{S5.} Macroscopic Brownian dynamics simulations of 10,000 particles diffusing through a stationary array of harmonic traps (60$\times$60 grid shown) over long length and time scales, varying the trap stiffness.

	\item[] \textbf{S6.} Macroscopic Brownian dynamics simulations of 10,000 particles diffusing through an oscillating array of stiff traps (60$\times$60 grid shown) over long length and time scales, varying the oscillation frequency at a fixed trap stiffness $\kappa = 5~kT/\mu\mathrm{m}^2$.

	\item[] \textbf{S7.}  2D contour plots of the displacement field density $d_{x}(x,y,t)$ in an $L\times L$ periodic cell containing an oscillating harmonic trap, varying the oscillation frequency at a fixed trap stiffness $\kappa = 5~kT/\mu\mathrm{m}^2$ (same parameters as in S6). Bottom row plots the long-time average over one periodic cycle.

\end{enumerate}