For reference, I am reading Lemma 2.2 here.
Suppose that $D\subseteq \mathbb{R}^n$ is non-empty bounded domain with piece-wise smooth boundary. Let $(u_j)$ be a sequence of Dirichlet eigenfunctions on $D$ with corresponding eigenvalues $$ 0 < \lambda_1^2\leq \lambda_2^2 \leq \dots. $$ For each $j\in\mathbb{N}$, define $$ h_j = \frac{1}{\lambda_j}. $$ Fix a function $a \in S(T^\ast D)$ and let $A = Op_{h_j}(a)$ be the quantization of $a$. We denote by $U(t;h_j)$ the Schrödinger propagator. That is, $$ U(t;h_j) = \exp\left({-\frac{itP_j}{h_j}}\right) $$ where $P_j = \frac{1}{2}\left(h_j^2\Delta -1\right)$.
As in the document linked above, we are able to make the following observation $$ \langle Au_j, u_j \rangle_{L^2} =\langle U(-t;h_j) A U(t;h_j) u_j, u_j \rangle_{L^2} =\langle Op_{h_j}(a\circ \varphi_t)u_j, u_j \rangle_{L^2} + \color{red}{O_t(h_j)} $$ where the last step follows from Egorov's theorem and $\varphi_t$ denotes the geodesic billiard flow.
The big-O notation means that our hidden term is bounded above by $C_t h_j$ for all $h_j$ sufficiently small. Can the constants $C_t$ be bounded uniformly in $t$? Also, since the bound only holds for $h_j$ small, does the 'how small' depend on $t$?
We then "average both sides over $t\in [0,T]$". I'm assuming this means integrating over $t\in [0,T]$ and dividing by $T$. This should yield $$ \langle Au_j, u_j \rangle_{L^2} = \color{red}{\langle Op_{h_j}\left(a_T\right) u_j, u_j \rangle_{L^2}} + {O_T(h_j)} $$ where $$ a_T(x, \omega) = \frac{1}{T}\int_0^T a\circ\varphi_t(x, \omega) \,\mathrm{d}{t}. $$
Is the term $\langle (a_T u_j, u_j \rangle_{L^2}$ obtained after changing the order of integration? If so, what integral representation are we using and how does one justify this application of the Fubini-Tonelli theorem?
Am I possibly misunderstanding anything? Any help is appreciated!