\begin{array}{c} J_{X,\lambda}(s_1)\\J_{X,\lambda}(s_2)\\ \vdots\\ J_{X,\lambda}(s_k) \end{array}\right) 
\stackrel{d}{\rightarrow} \mathcal N_C\Big( \left( \begin{array}{c} 0\\0\\ \vdots \\ 0 \end{array}\right), \underbrace{\left( \begin{array}{cccc} f_{X,\lambda}(s_1,s_1) &    f_{X,\lambda}(s_1,s_2) &  \ldots &  f_{X,\lambda}(s_1,s_k) \\
 f_{X,\lambda}(s_2,s_1) &    f_{X,\lambda}(s_2,s_2) & \ldots &  f_{X,\lambda}(s_2,s_k) \\ 
  \vdots & \vdots & \ldots & \vdots \\ f_{X,\lambda}(s_k,s_1) &    f_{X,\lambda}(s_k,s_2)& \ldots &  f_{X,\lambda}(s_k,s_k)\end{array} \right)}_{\displaystyle =\Sigma_\lambda} \Big),$$ where $\mathcal N_C$ denotes a circularly-symmetric complex Gaussian distribution with mean zero and complex-valued covariance matrix $ \Sigma_\lambda$. Furthermore, for two different frequencies $0 < \lambda_j \neq \lambda_k <\pi$, the corresponding vectors of finite Fourier transforms $J_{X,\lambda_j}$ and $J_{X,\lambda_k} $ are asymptotically independent; see, e.g., Theorem 5 in Cerovecki and Hörmann (2017). These properties of the finite Fourier transform as well as the fact that $ \widehat{p}_{X,\lambda} (\sigma,\tau) = J_{X,\lambda}(\sigma)\overline{J}_{X,\lambda}(\tau)$, for $\sigma, \tau \in [0,1]$, is the periodogram kernel, motivate the following bootstrap procedure to approximate the distribution of the test statistic $\mathcal U_T$ defined in (\[eq:test-stat-main\]).

1.   For $ \lambda_t=2\pi t/n$, $t=1,2, \ldots, N$, estimate the pooled spectral density operator $ \widehat{\mathcal F}_{\lambda_t}$ by $$\label{eq.sp-pooled}
      \widehat{\mathcal F}_{\lambda_t} = \frac{1}{2}\widehat{\mathcal F}_{X,\lambda_t} +  \frac{1}{2}\widehat{\mathcal F}_{Y,\lambda_t}$$ and denote by $ \widehat{f}_{\lambda_t}(\sigma,\tau)$, for $ \sigma, \tau \in \{s_1, s_2, \ldots, s_k\}$, the corresponding estimated pooled spectral density kernel.

2.   Generate two independent vectors 