probability going to $1$ as $K/n\rightarrow 0$—by a quantity which is of order
$p(K/n)^{1/3}$. Since $p<K$, then
$|\Gamma_{\theta}-\hat{\Gamma}_{\theta}|\rightarrow_{p}0$ when
$K^{4}/n\rightarrow 0$.

Lemma 4 implies that the oracle procedure for each individual regression
yields p-values which are independent uniform distributions so that
$\Gamma_{\theta}$ is a uniform random variable, and subsequently,
$\hat{\Gamma}_{\theta}$ is asymptotically a uniform random variable. ∎

##  Appendix F Proofs for Section 5.2

We first give an intermediate lemma, Lemma 7, which is a deterministic
statement about the power of the test. The lemma is then used to prove Theorem
2.

######  Lemma 7.

Fix an ordering $\theta$, $v\in V$, and $\mathbf{Y}\in\mathbb{R}^{n\times p}$.
Let $\Sigma={\mathbb{E}}([Z_{i},Y_{v,i}][Z_{i},Y_{v,i}]^{T})$ where
$Z_{i}\in\mathbb{R}^{K}$ represents the $K$ basis functions evaluated on
$Y_{{\rm pr}_{\theta}(v),i}$ and let $\lambda_{\min}$ and $\lambda_{\max}$
denote the smallest/largest eigenvalues of $\Sigma$. Let
$\Sigma_{Z,Z}={\mathbb{E}}(Z_{i}Z_{i}^{T})$ and
$\Sigma_{Z,v}={\mathbb{E}}(Z_{i}Y_{v,i})$ denote the sub-matrices of $\Sigma$.
Furthermore, let $W_{i}=\\{\\{h_{j}(Y_{u,i})\\}_{j\in[J],u\in{\rm
pr}_{\theta}(v)},Z_{i},Y_{v,i}\\}\in\mathbb{R}^{K+|{\rm pr}_{\theta}(v)|J+1}$,
$\delta_{1}=\|\frac{1}{n}\sum_{i}W_{i}W_{i}^{T}-{\mathbb{E}}(W_{i}W_{i}^{T})\|_{2}$,
and $g^{2}_{\max}=\max_{j\in[J],u\in{\rm
pr}_{\theta}(v)}\overline{h_{j}(Y_{u})^{2}}$. Suppose
${\delta_{1}}/{\lambda_{\min}}<1/2$ and Alg. 1 is used to test $H_{0}:{\rm
pa}(v)\subseteq{\rm pr}_{\theta}(v)\subseteq{\rm nd}(v)$. Then the resulting
p-value (when sampling an infinite number of bootstrap samples) is bounded
above by

(F.1) $\footnotesize 3\left(\frac{n}{|{\rm
pr}_{\theta}(v)|}\tau^{\star}_{v,{\rm
pr}_{\theta}(v)}\right)^{-1}\left\\{J\overline{Y_{v}^{2}}\times
g^{2}_{\max}+\frac{nJ}{|{\rm pr}_{\theta}(v)|}\delta_{1}^{2}\left[3+9\
\frac{\lambda^{2}_{\max}}{\lambda^{2}_{\min}}+(\delta_{1}^{2}+\lambda_{\max}^{2})\left(\frac{2(\lambda_{\max}/\lambda_{\min}+1)}{\lambda_{\min}}\right)^{2}\right]\right\\}.$

###### Proof.

Let $\hat{\gamma}=P_{\tilde{\eta}}\left(T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})\geq T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})\mid\mathbf{Y}\right)$ be the p-value resulting
from Alg. 1, where $P_{\tilde{\eta}}$ denotes the probability under the
bootstrap distribution conditional on $\mathbf{Y}$. By Markov’s inequality,
for any $d>0$,

(F.2) $\displaystyle P_{\tilde{\eta}}$
$\displaystyle\left(T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})\geq T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})\mid\mathbf{Y}\right)$ $\displaystyle\leq
P_{\tilde{\eta}}\left(\left|T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})-T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})+d\right|\geq d\mid\mathbf{Y}\right)$
$\displaystyle\leq\frac{{\mathbb{E}}_{\tilde{\eta}}\left[\left(T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})-d-T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})\right)^{2}\mid\mathbf{Y}\right]}{d^{2}}$
$\displaystyle\leq\frac{3(T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})-d)^{2}+3{\mathbb{E}}_{\tilde{\eta}}\left(T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})^{2}\mid\mathbf{Y}\right)}{d^{2}}.$

In the remainder of the proof, we let $d^{2}=\frac{n}{|{\rm
pr}_{\theta}(v)|}\tau^{\star}_{v,{\rm pr}_{\theta}(v)}$ and provide
appropriate bounds on
${\mathbb{E}}_{\tilde{\eta}}\left(T(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})^{2}\mid\mathbf{Y}\right)$ and
$(T(\mathbf{Y}_{v},{\rm pr}_{\theta}(v);\mathbf{Y})-d)^{2}$.

Conditioning on the observed data $\mathbf{Y}$ and computing the expectation
over the bootstrap draws, we have

$\displaystyle{\mathbb{E}}_{\tilde{\eta}}\left(T_{2}(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})^{2}\mid\mathbf{Y}\right)$
$\displaystyle=\frac{1}{|{\rm pr}_{\theta}(v)|}\sum_{u\in{\rm
pr}_{\theta}(v)}{\mathbb{E}}_{\tilde{\eta}}\left(\tau(\mathbf{\tilde{Y}}_{v},u,{\rm
pr}_{\theta}(v);\mathbf{Y})^{2}\right)$ $\displaystyle=\frac{1}{|{\rm
pr}_{\theta}(v)|}\sum_{u\in{\rm pr}_{\theta}(v)}\sum_{j\leq
J}{\mathbb{E}}_{\tilde{\eta}}\left(\frac{1}{n}(h_{j}(Y_{u})^{T}\tilde{\eta})^{2}\mid\mathbf{Y}\right)$
$\displaystyle\leq\frac{1}{|{\rm pr}_{\theta}(v)|}\sum_{u\in{\rm
pr}_{\theta}(v)}J\max_{j}{\mathbb{E}}_{\tilde{\eta}}\left(\left(\frac{1}{\sqrt{n}}h_{j}(Y_{u})^{T}\tilde{\eta}\right)^{2}\mid\mathbf{Y}\right).$

For each term, we have

$\displaystyle{\mathbb{E}}_{\tilde{\eta}}\left(\left|\frac{1}{\sqrt{n}}h_{j}(Y_{u})^{T}\tilde{\eta}\right|^{2}\right)$
$\displaystyle={\mathbb{E}}_{\tilde{\eta}}\left(\left|\frac{1}{\sqrt{n}}h_{j}(Y_{u})^{T}(\tilde{Y}_{v}-\mathbf{Z}_{v}\tilde{b})\right|^{2}\right)$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\frac{{\mathbb{E}}_{\tilde{\eta}}\left(\tilde{\eta}_{i}^{2}\right)}{n}h_{j}(Y_{u})^{T}(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})^{T}h_{j}(Y_{u})$
$\displaystyle\leq\frac{{\mathbb{E}}_{\tilde{\eta}}\left(\tilde{\eta}_{i}^{2}\right)}{n}h_{j}(Y_{u})^{T}(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})^{T}h_{j}(Y_{u})$
$\displaystyle\leq\frac{\frac{1}{n}\sum_{i}\hat{\eta}_{i}^{2}}{n}h_{j}(Y_{u})^{T}(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})(I-\mathbf{Z}_{v}[\mathbf{Z}_{v}^{T}\mathbf{Z}_{v}]^{-1}\mathbf{Z}_{v}^{T})^{T}h_{j}(Y_{u})$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\overline{Y_{v}^{2}}\times\overline{h_{j}(Y_{u})^{2}},$

where (i) holds since each element of $\tilde{\eta}$ is independent when
conditioning on $\mathbf{Y}$ and ${\mathbb{E}}(\tilde{\eta}_{i})=0$, so that
${\mathbb{E}}_{\tilde{\eta}}\left(\tilde{\eta}\tilde{\eta}^{T}\right)=\frac{1}{n}\sum_{i}\hat{\eta}_{i}^{2}I$,
and (ii) holds because $I-\mathbf{Z}[\mathbf{Z}^{T}\mathbf{Z}]^{-1}\mathbf{Z}$
is idempotent and
$\frac{1}{n}\sum_{i}\eta_{i}^{2}\leq\frac{1}{n}\sum_{i}Y_{v,i}^{2}$, since
$\eta$ are the residuals after regressing $Y_{v}$ onto $\mathbf{Z}$. Thus,

(F.3) ${\mathbb{E}}_{\tilde{\eta}}\left(T_{2}(\mathbf{\tilde{Y}}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})^{2}\right)\leq J\bar{Y_{v}^{2}}g^{2}_{\max}.$

For the term $(T(\mathbf{Y}_{v},{\rm pr}_{\theta}(v);\mathbf{Y})-d)^{2}$, we
have

(F.4) $\displaystyle\left(T(\mathbf{Y}_{v},{\rm
pr}_{\theta}(v);\mathbf{Y})-d\right)^{2}$ $\displaystyle=\frac{1}{|{\rm
pr}_{\theta}(v)|}\left(\|\\{\tau(\mathbf{Y}_{v},u,{\rm
pr}_{\theta}(v);\mathbf{Y}\\}_{u\in{\rm
pr}_{\theta}(v)}\|-\|\sqrt{n}\\{\|\tau^{\star}_{v,u,{\rm
pr}_{\theta}(v)}\|\\}_{u\in{\rm pr}_{\theta}(v)}\|\right)^{2}$
$\displaystyle\leq\frac{1}{|{\rm
pr}_{\theta}(v)|}\left\|\\{\tau(\mathbf{Y}_{v},u,{\rm
pr}_{\theta}(v);\mathbf{Y})-\sqrt{n}\|\tau^{\star}_{v,u,{\rm
pr}_{\theta}(v)}\|\\}_{u\in{\rm pr}_{\theta}(v)}\right\|_{2}^{2}$
$\displaystyle\leq\frac{1}{|{\rm pr}_{\theta}(v)|}\sum_{u\in{\rm
pr}_{\theta}(v)}\left(\tau(\mathbf{Y}_{v},u,{\rm
pr}_{\theta}(v);\mathbf{Y})-\sqrt{n}\|\tau^{\star}_{v,u,{\rm
pr}_{\theta}(v)}\|\right)^{2}$ $\displaystyle=\frac{n}{|{\rm
pr}_{\theta}(v)|}\sum_{u\in{\rm
pr}_{\theta}(v)}\left(\|\\{\frac{1}{n}h_{j}(\mathbf{Y}_{u})^{T}\hat{\eta}\\}_{j\in[J]}\|-\|\\{{\mathbb{E}}\left[h_{j}(Y_{u,i})(Y_{v,i}-Z_{i}^{T}b)\\}_{j\in[J]}\right]\|\right)^{2}$
$\displaystyle\leq\frac{n}{|{\rm pr}_{\theta}(v)|}\sum_{u\in{\rm
pr}_{\theta}(v)}\left\|\left\\{\frac{1}{n}h_{j}(\mathbf{Y}_{u})^{T}\hat{\eta}-{\mathbb{E}}\left[h_{j}(Y_{u,i})(Y_{v,i}-Z_{i}^{T}b)\right]\right\\}_{j\in[J]}\right\|^{2}$
$\displaystyle=\frac{n}{|{\rm
pr}_{\theta}(v)|}\sum_{j\in[J]}\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\hat{\eta}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})(Y_{v,i}-Z_{i}^{T}b)\right]\right\|^{2}$
$\displaystyle\leq\frac{n}{|{\rm pr}_{\theta}(v)|}J\max_{j\in
J}\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\hat{\eta}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})(Y_{v,i}-Z_{i}^{T}b)\right]\right\|^{2}.$

For any $j\in J$, we have

(F.5) $\displaystyle\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\hat{\eta}\vphantom{\frac{1}{n}}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})(Y_{v,i}-Z_{i}^{T}b)\right]\right\|^{2}$
$\displaystyle=\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}(\mathbf{Y}_{v}-\mathbf{Z}\hat{b})-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})(Y_{vi}-Z_{i}^{T}b)\right]\right\|^{2}$ $\displaystyle\leq
3\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Y}_{v}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{vi}\right]\right\|^{2}$
$\displaystyle\quad+3\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}\hat{b}-{\mathbb{E}}\left[h_{j}(Y_{ui})Z_{i}^{T}b)\right]\right\|^{2}$
$\displaystyle\leq 3\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Y}_{v}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{vi}\right]\right\|^{2}$
$\displaystyle\quad+9\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}b-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Z_{i}^{T}b)\right]\right\|^{2}+9\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}(\hat{b}-b)\right\|^{2}$ $\displaystyle\leq
3\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Y}_{v}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{vi}\right]\right\|^{2}$
$\displaystyle\quad+9\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Z_{i}^{T}\right]\right\|_{2}^{2}\|b\|^{2}+9\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}\right\|^{2}\left\|\hat{b}-b\right\|^{2}.$

Let $E=\frac{1}{n}\mathbf{Z}^{T}\mathbf{Z}-\Sigma_{\mathbf{Z},\mathbf{Z}}$ and
$e=\frac{1}{n}\mathbf{Z}^{T}Y_{v}-\Sigma_{Z,v}$. By assumption, $\|e\|$ and
$\|E\|$ are both bounded above by $\delta_{1}$. Following the proof of Wang
and Drton, 2020b (, Lemma 1), we have
$\|\Sigma_{Z,Z}^{-1}E\|_{2}<\|\Sigma_{Z,Z}^{-1}\|_{2}\|E\|_{2}\leq\delta_{1}/\lambda_{\min}<1/2$.
Furthermore, we have

$\displaystyle\|\hat{b}-b\|$
$\displaystyle=|(\Sigma_{Z,Z}+E)^{-1}(\Sigma_{Z,v}+e)-\Sigma_{Z,Z}^{-1}\Sigma_{Z,v}|_{2}$
$\displaystyle\leq|(\Sigma_{Z,Z}+E)^{-1}\Sigma_{Z,v}-\Sigma_{Z,Z}^{-1}\Sigma_{Z,v}|_{2}+|(\Sigma_{Z,Z}+E)^{-1}e|_{2}$
$\displaystyle\leq\frac{\|\Sigma_{Z,Z}^{-1}E\|_{2}}{1-\|\Sigma_{Z,Z}^{-1}E\|_{2}}|b|_{2}+\frac{1/\lambda_{\min}}{{1-\|\Sigma_{Z,Z}^{-1}E\|_{2}}}|e|_{2}$
$\displaystyle\leq
2\|\Sigma_{Z,Z}^{-1}\|_{2}\|E\|_{2}|b|_{2}+\frac{2}{\lambda_{\min}}|e|_{2}$
$\displaystyle\leq
2\delta_{1}\|\Sigma_{Z,Z}^{-1}\|_{2}\|b\|_{2}+\frac{2}{\lambda_{\min}}\delta_{1}$
$\displaystyle\leq\frac{2\delta_{1}}{\lambda_{\min}}(\|b\|_{2}+1).$

Since $\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Y}_{v}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{vi}\right]$ is a sub-matrix of
$\frac{1}{n}\sum_{i}W_{i}W_{i}^{T}-{\mathbb{E}}(W_{i}W_{i}^{T})$, we have
$\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Y}_{v}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{vi}\right]\|_{2}\leq\|\frac{1}{n}\sum_{i}W_{i}W_{i}^{T}-{\mathbb{E}}(W_{i}W_{i}^{T})\|\leq\delta_{1}$.
Similarly $\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Z_{i}^{T}\right]\|_{2}\leq\delta_{1}$. In addition,
$\|b\|_{2}=\|\Sigma_{Z,Z}^{-1}\Sigma_{Z,v}\|_{2}\leq\lambda_{\max}/\lambda_{\min}$.
Putting these bounds together, we have

(F.6) $\displaystyle\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\hat{\eta}\vphantom{\frac{1}{n}}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})(Y_{v,i}-Z_{i}^{T}b)\right]\right\|^{2}$ $\displaystyle\leq
3\delta_{1}^{2}+9\delta_{1}^{2}\|b\|^{2}_{2}+9\left\|\hat{b}-b\right\|_{2}^{2}\left(\left\|\frac{1}{n}h_{j}(\mathbf{Y}_{{\rm
pr}_{\theta}(v)})^{T}\mathbf{Z}-{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Z_{i}^{T}\right]\right\|_{2}^{2}+\left\|{\mathbb{E}}\left[h_{j}(Y_{{\rm
pr}_{\theta}(v),i})Y_{{\rm pr}_{\theta}(v),i}\right]\right\|^{2}\right)$
$\displaystyle\leq
3\delta_{1}^{2}+9\delta_{1}^{2}\frac{\lambda^{2}_{\max}}{\lambda^{2}_{\min}}+9(\delta_{1}^{2}+\lambda_{\max}^{2})\left(\frac{2(\lambda_{\max}/\lambda_{\min}+1)}{\lambda_{\min}}\delta_{1}\right)^{2}.$

The result follows by combining (F.3), (F.4), and
(LABEL:eq:boundOnAltDeviation). ∎

### Theorem 2

Suppose $K=O(p)$ and consider a sequence of models where $p,n$ grow, but
$\lambda_{\min},\lambda_{\max},M_{1}$, and $J$ are fixed. Suppose ${\rm
pa}(v)\subseteq{\rm pr}_{\theta}(v)\subseteq{\rm nd}(v)$ does not hold, then
for any fixed level $\alpha$, the hypothesis $H_{0,\theta,v}$ will be rejected
with probability going to $1$ if

  1. (1)

$W_{i}$ is jointly sub-Gaussian with Orlicz-2 norm $M_{1}$,
$\max_{v,j}{\mathbb{E}}(g_{j}(Y_{v})^{2})\leq M_{1}$,
$\max_{v,j}{\mathbb{E}}(Y_{v}^{2})\leq M_{1}$, and $p/n=o(\tau^{\star}_{v,{\rm
pr}_{\theta}(v)})$ and $p/n\rightarrow 0$; or

  2. (2)

$W_{i}$ is jointly log-concave with Orlicz-1 norm $M_{1}$ and
$p\log(n)^{3}/n=o(\tau^{\star}_{v,{\rm pr}_{\theta}(v)})$ and
$p\log(n)^{3}/n\rightarrow 0$,

where $W_{i}=\left[Y_{v,i},\;\left(\phi_{k}(Y_{{\rm
pr}_{\theta}(v)}):k\in[K]\right),\;\left(h_{j}(Y_{u,i}):j\in[J],u\in{\rm
pr}_{\theta}(v)\right)\right]\in\mathbb{R}^{K+|{\rm pr}_{\theta}(v)|J+1}$.

###### Proof.

The proof simply requires showing that the numerator in (F.3) can be bounded
above by the appropriate quantity with probability going to $1$.

Sub-Gaussian case: Let $t=\log(p)/\sqrt{n}$. Since $Y_{v,i}^{2}$ and
$g_{j}(Y_{u,i})^{2}$ are sub-exponential with Orlicz-$1$ norm bounded by
$M_{1}^{2}$, we have

(F.7) $\displaystyle P\left(\overline{Y_{v}^{2}}\times
g^{2}_{\max}>(M_{1}+t)^{2}\right)$ $\displaystyle\leq
P\left(\overline{Y_{v}^{2}}\geq(M_{1}+t)\cup g^{2}_{\max}\geq(M_{1}+t)\right)$
$\displaystyle\leq
P\left(\overline{Y_{v}^{2}}\geq(M_{1}+t)\right)+pJ\max_{j,u}P\left(\overline{g^{2}_{j}}\geq(M_{1}+t)\right)$
$\displaystyle\leq
2(1+pJ)\exp\left(-\frac{1}{C}\min(\frac{t^{2}}{M_{1}^{4}(1/n)},\left(\frac{t}{M_{1}^{2}(1/n)})\right)\right)$
$\displaystyle\leq
2\exp\left(\log(1+pJ)-\frac{1}{C}\frac{\log(p)^{2}}{M_{1}^{4}}\right).$

Additionally, by Exercise 4.7.3 of Vershynin, (2018),

$P\left[\delta_{1}>cM_{1}^{2}\lambda_{\max}\left(\sqrt{\frac{K+|{\rm
pr}_{\theta}(v)|J+1+t}{n}}+\frac{K+|{\rm
pr}_{\theta}(v)|J+1+t}{n}\right)\right]\leq 2\exp(-t).$

Letting $t=p$ and assuming $n$ is large enough so that $(K+p(J+1)+1)/n<1$,

(F.8) $\displaystyle P$
$\displaystyle\left[\delta_{1}>cM_{1}^{2}\lambda_{\max}2\sqrt{\frac{(K+p(J+1)+1)}{n}}\right]$
$\displaystyle\leq
P\left[\delta_{1}>cM_{1}^{2}\lambda_{\max}\left(\sqrt{\frac{(K+p(J+1)+1)}{n}}+\frac{(K+p(J+1)+1)}{n}\right)\right]\leq
2\exp(-p).$

Therefore, with probability at least
$1-2\exp(-p)-2\exp\left(\log((pJ+1))-\frac{1}{C}\frac{\log(p)^{2}}{M_{1}^{4}}\right)$,
the numerator in (F.1) is bounded above by

$3(M_{1}+\log(p)/\sqrt{n})^{2}+\frac{3nJ}{|{\rm
pr}_{\theta}(v)|}\left(cM_{1}^{2}\lambda_{\max}2\sqrt{\frac{K+p(J+1)+1}{n}}\right)^{2}\\\
\times\left[3+9\frac{\lambda_{\max}^{2}}{\lambda_{\min}^{2}}+(\left(cM_{1}^{2}\lambda_{\max}\right)^{2}+\lambda_{\max}^{2})\left(\frac{2(\lambda_{\max}/\lambda_{\min}+1)}{\lambda_{\min}}\right)^{2}\right].$

Thus, when $K=O(p)$ and $p/n=o(\tau^{\star}_{v,{\rm pr}_{\theta}(v)})$, the
entire quantity in (F.1) goes to $0$.

Log-concave case: Let $t=\log(p)/\sqrt{n}$. Since $Y_{v,i}^{2}$ and
$g_{j}(Y_{u,i})^{2}$ have Orlicz-$1/2$ norm bounded by $M_{1}^{2}$, we have

(F.9) $\displaystyle P\left(\overline{Y_{v}^{2}}\times
g^{2}_{\max}>(M_{1}+t)^{2}\right)$ $\displaystyle\leq
P\left(\overline{Y_{v}^{2}}\geq(M_{1}+t)\cup g^{2}_{\max}\geq(M_{1}+t)\right)$
$\displaystyle\leq
P\left(\overline{Y_{v}^{2}}\geq(M_{1}+t)\right)+pJ\max_{j,u}P\left(\overline{g^{2}_{j}}\geq(M_{1}+t)\right)$
$\displaystyle\leq
2(1+pJ)\exp\left(-\frac{1}{C}\min(\frac{t^{2}}{M_{1}^{4}(1/n)},\left(\frac{t}{M_{1}^{2}(1/n)})\right)^{1/2}\right)$
$\displaystyle\leq
2\exp\left(\log(1+pJ)-\frac{1}{C}\frac{\log(p)^{2}}{M_{1}^{4}}\right).$

To control $\delta_{1}$, we apply Theorem 4.1 of Adamczak et al., (2010),
which states that when $X_{i}$ are i.i.d. log-concave random vectors in
$\mathbb{R}^{p}$, then for any $\varepsilon\in(0,1)$ and constants $c,C$, when
$n\geq C\varepsilon^{-2}\log(2\varepsilon^{-2})^{2}$, then with probability at
least $1-\exp(-c\sqrt{p})$, we have for the sample covariance $S$,

$\|S-\Sigma_{Z,Z}\|\leq\varepsilon\|\Sigma\|.$

Let $\varepsilon=\sqrt{\frac{(K+p(J+1)+1)\log(n)^{3}}{n}}$. The conditions are
satisfied when $p,n$ are large enough so that $(K+p(J+1)+1)\log(n)^{3}/n<1$
and $1\geq
C\log\left[\frac{2n}{(K+p(J+1)+1)\log(n)^{3}}\right]^{2}/\log(n)^{3}$ for the
universal constant $C$. Thus,

$\|S-\Sigma_{Z,Z}\|\leq\sqrt{\frac{((K+p(J+1)+1)\log(n)^{3}}{n}}\lambda_{\max}$

and $\delta_{1}^{2}\leq\frac{(K+pJ+1)\log(n)^{3}}{n}\lambda_{\max}^{2}$ with
probability at least $1-\exp(-c\sqrt{(K+|{\rm pr}_{\theta}(v)|J+1})$. Thus,
with probability going to $1$, we have that the numerator in (F.1) is bounded
above by

$\frac{nJ}{|{\rm
pr}_{\theta}(v)|}\frac{(K+pJ+1)\log(n)^{3}\lambda_{\max}^{2}}{n}\left[3+9\frac{\lambda_{\max}^{2}}{\lambda_{\min}^{2}}+(2\lambda_{\max}^{2})\left(\frac{2(\lambda_{\max}/\lambda_{\min}+1)}{\lambda_{\min}}\right)^{2}\right]\\\
+(M_{1}+\log(p)/\sqrt{n})^{2}.$

Since the denominator is $\frac{n}{|{\rm pr}_{\theta}(v)|}\tau^{\star}_{v,{\rm
pr}_{\theta}(v)}$, when $K=O(p)$ and $p\log(n)^{3}/n=o(\tau^{\star}_{v,{\rm
pr}_{\theta}(v)})$, the entire quantity in (F.1) goes to $0$. ∎

##  Appendix G Full Simulation Details

We provide the results for the comparison of goodness-of-fit procedures for
$p=10,15,20,30,45$.

Table 5. Comparison of testing procedures when $n\approx p^{5/4}$. Each proportion in the table has been multiplied by 1000. In Cols 3-8 Here, the value in the table is bold if it exceeds the nominal size of $.1$ by 2 standard deviations. In Cols 9-14, the value in the table is bold if the procedure has the largest power (or is within 2 standard deviations of the procedure with the largest power) for that particular setting. For each setting, if the procedure’s empirical size is significantly above the nominal level, then we do not display it’s empirical power. In the columns, ‘W’ is our proposed procedure, ‘S’ indicates Sen and Sen, (2014), ‘RO’ and ‘RL’ indicate the OLS and lasso variants of Shah and Bühlmann, (2018), ‘M’ indicates Berrett and Samworth, (2019), and ‘H’ indicates Schultheiss et al., (2023). | Size | Power
---|---|---
dist | p | W | S | RO | RL | M | H | W | S | RO | RL | M | H
gamma | 10 | 42 | 50 | 92 | 84 | 114 | 0 | 144 | 68 | 182 | 248 | 144 | 0
15 | 72 | 88 | 76 | 100 | 116 | 0 | 268 | 100 | 262 | 340 | 174 | 0
20 | 88 | 90 | 88 | 78 | 114 | 0 | 264 | 156 | 242 | 346 | 200 | 0
30 | 98 | 110 | 40 | 86 | 122 | 0 | 322 | 180 | 276 | 490 | 196 | 0
45 | 104 | 192 | 36 | 204 | 186 | 10 | 332 |  | 330 |  |  | 14
laplace | 10 | 52 | 38 | 88 | 82 | 76 | 0 | 52 | 60 | 126 | 90 | 122 | 0
15 | 74 | 60 | 100 | 86 | 80 | 0 | 122 | 98 | 106 | 78 | 130 | 0
20 | 80 | 94 | 78 | 78 | 108 | 0 | 112 | 96 | 118 | 98 | 132 | 0
30 | 90 | 110 | 76 | 84 | 112 | 0 | 162 | 86 | 108 | 104 | 178 | 0
45 | 76 | 142 | 38 | 164 | 106 | 8 | 136 |  | 100 |  | 158 | 10
lognorm | 10 | 60 | 70 | 88 | 80 | 88 | 0 | 234 | 88 | 248 | 236 | 202 | 0
15 | 64 | 100 | 70 | 80 | 144 | 0 | 418 | 172 | 368 | 390 |  | 0
20 | 116 | 130 | 42 | 76 | 220 | 0 | 478 |  | 358 | 410 |  | 0
30 | 76 | 198 | 32 | 106 | 312 | 0 | 558 |  | 412 | 562 |  | 0
45 | 106 | 198 | 16 | 236 | 506 | 36 | 578 |  | 486 |  |  | 62
mixed | 10 | 58 | 62 | 108 | 66 | 94 | 0 | 112 | 60 | 146 | 188 | 126 | 0
15 | 72 | 98 | 98 | 80 | 108 | 0 | 248 | 104 | 244 | 296 | 208 | 0
20 | 74 | 102 | 60 | 66 | 112 | 0 | 324 | 140 | 238 | 372 | 188 | 0
30 | 68 | 122 | 56 | 84 | 146 | 0 | 344 | 184 | 274 | 484 |  | 0
45 | 98 | 156 | 58 | 180 | 186 | 12 | 438 |  | 358 |  |  | 42
unif | 10 | 28 | 34 | 98 | 90 | 80 | 0 | 38 | 28 | 90 | 140 | 56 | 0
15 | 72 | 38 | 122 | 90 | 62 | 0 | 44 | 70 | 72 | 134 | 66 | 0
20 | 88 | 74 | 98 | 72 | 78 | 0 | 48 | 58 | 82 | 118 | 48 | 0
30 | 78 | 66 | 108 | 100 | 70 | 0 | 78 | 102 | 106 | 160 | 48 | 0
45 | 112 | 112 | 104 | 200 | 48 | 12 | 48 | 110 | 96 |  | 62 | 8
weibull | 10 | 54 | 54 | 96 | 108 | 104 | 0 | 244 | 76 | 240 | 298 | 216 | 0
15 | 82 | 102 | 72 | 98 | 124 | 0 | 346 | 154 | 330 | 384 | 272 | 0
20 | 76 | 132 | 50 | 70 | 146 | 0 | 442 |  | 386 | 480 |  | 0
30 | 86 | 188 | 44 | 78 | 222 | 0 | 490 |  | 420 | 612 |  | 0
45 | 98 | 174 | 18 | 168 | 258 | 14 | 516 |  | 472 |  |  | 20
Table 6. Comparison of testing procedures when $n\approx p^{2}$. Each proportion in the table has been multiplied by 1000. In Cols 3-8 Here, the value in the table is bold if it exceeds the nominal size of $.1$ by 2 standard deviations. In Cols 9-14, the value in the table is bold if the procedure has the largest power (or is within 2 standard deviations of the procedure with the largest power) for that particular setting. For each setting, if the procedure’s empirical size is significantly above the nominal level, then we do not display it’s empirical power. In the columns, ‘W’ is our proposed procedure, ‘S’ indicates Sen and Sen, (2014), ‘RO’ and ‘RL’ indicate the OLS and lasso variants of Shah and Bühlmann, (2018), ‘M’ indicates Berrett and Samworth, (2019), and ‘H’ indicates Schultheiss et al., (2023). | Size | Power
---|---|---
dist | p | W | S | RO | RL | M | H | W | S | RO | RL | M | H
gamma | 10 | 66 | 160 | 14 | 76 | 106 | 94 | 860 |  | 806 | 792 | 922 | 458
15 | 98 | 174 | 2 | 104 | 106 | 116 | 934 |  | 950 | 946 | 994 | 730
20 | 102 | 172 | 0 | 62 | 60 | 118 | 980 |  | 994 | 996 | 998 | 896
30 | 90 | 146 | 0 | 102 | 66 | 120 | 988 |  | 1000 | 1000 | 1000 | 962
45 | 116 | 144 | 0 | 324 | 74 | 136 | 992 |  | 998 |  | 1000 |
laplace | 10 | 90 | 126 | 38 | 106 | 120 | 88 | 342 | 304 | 182 | 180 | 496 | 318
15 | 92 | 114 | 10 | 82 | 238 | 130 | 384 | 376 | 236 | 266 |  |
20 | 90 | 104 | 4 | 92 | 340 | 90 | 412 | 438 | 312 | 350 |  | 652
30 | 98 | 116 | 0 | 124 | 436 | 118 | 466 | 492 | 458 | 620 |  | 846
45 | 124 | 102 | 2 | 280 | 522 | 98 | 550 | 554 | 660 |  |  | 958
lognorm | 10 | 100 | 258 | 12 | 72 | 440 | 140 | 948 |  | 880 | 798 |  |
15 | 98 | 212 | 0 | 74 | 602 | 148 | 982 |  | 970 | 944 |  |
20 | 78 | 184 | 0 | 72 | 742 | 168 | 994 |  | 998 | 994 |  |
30 | 94 | 176 | 0 | 110 | 892 | 162 | 996 |  | 1000 | 1000 |  |
45 | 102 | 180 | 0 | 322 | 984 | 192 | 1000 |  | 1000 |  |  |
mixed | 10 | 92 | 204 | 44 | 102 | 144 | 108 | 834 |  | 774 | 752 |  | 606
15 | 102 | 172 | 40 | 72 | 136 | 116 | 932 |  | 944 | 950 |  | 884
20 | 120 | 186 | 34 | 76 | 166 | 102 | 958 |  | 970 | 978 |  | 966
30 | 82 | 162 | 40 | 116 | 244 | 102 | 980 |  | 994 | 986 |  | 998
45 | 92 | 112 | 70 | 338 | 302 | 128 | 986 | 978 | 994 |  |  |
unif | 10 | 96 | 100 | 122 | 100 | 12 | 62 | 38 | 154 | 200 | 278 | 56 | 212
15 | 114 | 102 | 170 | 106 | 0 | 86 | 82 | 194 |  | 426 | 6 | 438
20 | 114 | 112 | 176 | 104 | 0 | 84 | 136 | 252 |  | 640 | 2 | 616
30 | 114 | 102 | 248 | 148 | 0 | 100 | 172 | 232 |  |  | 2 | 800
45 | 108 | 118 | 334 | 276 | 0 | 84 | 260 | 280 |  |  | 0 | 926
weibull | 10 | 82 | 212 | 22 | 96 | 306 | 114 | 918 |  | 906 | 830 |  | 754
15 | 114 | 224 | 0 | 92 | 310 | 158 | 974 |  | 974 | 972 |  |
20 | 108 | 186 | 0 | 58 | 272 | 122 | 992 |  | 1000 | 998 |  | 990
30 | 90 | 214 | 0 | 124 | 286 | 110 | 992 |  | 1000 | 1000 |  | 1000
45 | 94 | 176 | 0 | 348 | 380 | 148 | 998 |  | 1000 |  |  |


# Shaping High-Order Diffraction-Free Beams Through Continuous Superposition
of Bessel Beams

Michel Zamboni-Rached1, Jéssyca Nobre-Pereira1 and João Quaglio1 E-mail:
<EMAIL_ADDRESS>
(1Universidade Estadual de Campinas, Faculdade de Engenharia Elétrica e
Computação, Campinas, SP, Brasil.)

###### Abstract

Recognized for their non-diffracting properties, Bessel beams can be
conveniently combined to generate the so-called Frozen Waves, which are
monochromatic beams endowed with topological charge and whose longitudinal
intensity pattern can be shaped according to a previously chosen function.
Continuous superposition of Bessel beams is specially suitable for micrometer-
scale domains, being highly relevant for applications in optical tweezers,
particle trapping and atom guidance. Previous studies have successfully
constructed micrometer Frozen Wave solutions with null topological charge;
nevertheless, some challenges persist in obtaining exact solutions when
dealing with higher topological charge values (i.e., higher-order micrometer
Frozen Waves). Typically, a topological charge raising operator is used to
elevate the order; however, such solutions face significant issues when the
orders become excessively high. In this paper, based in continuous
superposition of higher order Bessel beams, we develop a novel analytic and
exact solution for higher-order micrometer Frozen Waves, so obtaining a method
for modeling the intensity of beams with arbitrary topological charge within
micrometer spatial domains. We also investigate the behavior of the electric
field’s longitudinal component in such highly non-paraxial regimes through a
vectorial treatment.

##  1 Introduction

It is well known that solutions of the cylindrical Helmholtz equation given by
$\psi(\rho,\phi,z)=J_{\nu}(\eta\rho)\exp(i\beta z)\exp(i\nu\phi)$, and known
as Bessel beams [1], are resistant to diffraction in non-guided media. Here,
$\eta=\sqrt{k^{2}-\beta^{2}}$ and $\beta$ are the transverse and longitudinal
wavenumbers respectively, $k=\omega/c$ is the wavenumber ($\omega$ is the
angular frequency and $c$ the light speed), and $\nu$ is the beam’s
topological charge, also referred as the beam’s order. Alongside plane waves,
which are also non-diffracting, ideal Bessel beams are not physically
realizable since their transverse fields are not square integrable, i.e., they
require an infinite energy flux through any transverse plane to maintain their
self-healing indefinitely along the propagation. In practice, we can relay on
truncated versions [2] or gaussian-apodized versions, as the Bessel-Gauss
beams [3, 4]. Such type of beams is also referred to as Localized Waves or
non-diffracting waves[5, 6], and significant progress has been made in this
research field, whether in theory, experiments, or applications[7, 8, 9, 10].

Over recent years, superposition of Bessel beams has been used for modelling
the longitudinal intensity pattern (LIP) of diffraction-resistant beams. Such
method is called Frozen Wave (FW), as an allusion to the static envelope [11,
12, 13, 14] of the resulting beams. The Frozen Wave (FW) method has proven
versatile, enabling the creation of structured light beams endowed with
intriguing properties [15, 16], even in challenging scenarios like absorbing
media [17, 18] and stratified dielectric structures [19, 20]. For instance,
numerous studies have explored the theoretical and experimental control of
polarization [21, 22] and orbital angular momentum density [23, 24] along the
propagation of optical beams.

The FW beams beams may have applications across various fields, as in optical
tweezers [14, 25, 26], atom or particle guidance [27], acoustics [24],
holography [28] and are also promising in high-resolution microscopy [29],
microfabrication and micromanipulation [30, 31].

In this paper, we initially introduce the well-established mathematical method
[32, 33] for deriving a solution for micrometer zero-order FWs through the
continuous superposition of zero-order Bessel beams, resulting in a
superposition of Mackinnon-type beams [34]. However, generating a higher-order
FW, i.e., a FW of order $\nu\neq 0$, represents a significant challenge due to
the lack of analytical solutions for continuous superposition of higher-order
Bessel beams. Given this difficulty, the commonly used approach involves
applying an operator, here referred to as the Topological Charge Raising
Operator (TCRO), $\nu$ times to the zero-order FW solution to achieve the
desired $\nu$-order. Essentially, this method constitutes a recursive
differentiation process, and in [25], we were able to systematize it by
presenting the resulting field as a closed mathematical solution in terms of
spherical Bessel functions, so enabling the generation of $\nu$ order FWs.
However, the TCRO approach proved to work for a limited range of orders, as it
modifies the original spectrum (i.e., the spectrum of the original zero-order
FW) and thus can affect the resulting Longitudinal Intensity Pattern (LIP).

To address this issue, we introduce a novel method, that we call the Expansion
Method, for the generation of FWs of arbitrary order. It employs calculations
based on Fourier series expansions, resulting in promising advancements for
generating higher-order FWs beams in micrometer spatial regions. Due to the
highly non-paraxial behavior of these beams, we also conducted vectorial
treatment based on Maxwell’s equations to account for the longitudinal
component of the electric field.

The paper is structured as follows: Section 2.1 details the mathematical
construction of micrometer zero-order Frozen Waves (FWs). Subsequently,
Section 2.2 elucidates the TCRO method, followed by the expansion method and
its adaptation for a simple vectorial case in Section 2.3. Both methods are
illustrated with examples. Section 3 provides our general remarks and
conclusion.

##  2 Mathematical Methodology

###  2.1 The FW method based on continuous superposition of zero-order Bessel
beams

An efficient way to obtain optical beams with null topological charge and
structured in micrometer regions is through the FW method based on a
superposition of forward-propagating and equal frequency Bessel beams with
distinct longitudinal (and consequently transverse) wavenumbers $\beta$, that
is111For simplicity, in this paper we will omit the time factor $\exp(-i\omega
t)$. [32],

$\Psi_{0}(\rho,z)=\int_{-\omega/c}^{\omega/c}S(\beta)J_{0}(\eta(\beta)\rho)e^{i\beta
z}d\beta,$ (1)

where

$\eta(\beta)=\sqrt{k^{2}-\beta^{2}}$ (2)

is the transverse wave number.

The aim of the method is to obtain, through the superposition (1), a beam
whose LIP along the $z$ axis is given by the modulus square of a function
$F(z)$ of the form

$F(z)=f(z)\mathrm{e}^{iQz}\,\,,$ (3)

where $f(z)$, referred to as the morphological function, is a function of our
choice, and the parameter $Q$ is a positive constant, specifically a fraction
of $k=\omega/c$, and its role is to ensure that $S(\beta)$ is negligible for
$\beta<0$, thereby ensuring that the superposition is primarily composed of
forward-traveling waves. Mathematically, we demand that
$|\Psi_{0}(\rho=0,z,t)|^{2}\approx|F(z)|^{2}$, and the FW method demonstrates
that this can be achieved by considering the spectral function $S(\beta)$
represented by the following Fourier series.

$S(\beta)=\sum_{n=-\infty}^{\infty}A_{n}\exp\left(in\frac{2\pi}{K}\beta\right)\,\,,$
(4)

with

$A_{n}=\frac{1}{K}F({-2\pi n}/{K})\,\,,$ (5)

where $K=2\omega/c$. The solution of 1 with the spectrum 4 is given the
following superposition of Mackinnon-type beams [34]

$\Psi_{0}\left(\rho,z,t\right)=\mathrm{e}^{-i\omega
t}\sum_{n=-\infty}^{\infty}F\left(-\frac{2n\pi}{K}\right)\mathrm{sinc}\left[\xi(\rho,z)\right]\text{,}$
(6)

with:

$\xi(\rho,z)=\sqrt{\frac{\omega^{2}}{c^{2}}\rho^{2}+\left(n\pi+\frac{\omega}{c}z\right)^{2}}\text{.}$
(7)

Solution 6 represents a zero-order Frozen Wave (i.e., $\nu=0$, and thus null
topological charge), whose Longitudinal Intensity Pattern (LIP) can be easily
modeled in micrometer-scale regions using a function $F(z)$ of our choice. The
spot radius of the resulting beam is approximately given by $\rho_{0}\approx
2.4/\left(k^{2}-Q^{2}\right)^{1/2}$.

###  2.2 The Topological Charge Raising Operator Method

In some contexts, such as atom guidance and occasionally in optical tweezers,
obtaining solutions for higher-order Frozen Waves ($\nu=1,2,3,...$) becomes
essential. Due to the lack of a known method for solving the integral in
Eq.(1) with (4) and the integrand given by a Bessel beam of order $\nu\neq 0$,
the procedure used to obtain such higher-order beams involves the utilization
of the following operator:

$U=e^{i\phi}\Big{[}\frac{\partial}{\partial\rho}+\frac{i}{\rho}\frac{\partial}{\partial\phi}\Big{]}\text{,}$
(8)

which is referred to here as the Topological Charge Raising Operator (TCRO)
due to its ability to increase the order of a beam by one unit if it possesses
a well defined topological charge. This method relies on applying the TCRO
$\nu$ times to a zero-order Frozen Wave (FW) to obtain a $\nu$-order solution
(symbolically, we write $U\Psi_{\nu}=\Psi_{\nu+1}$). Unlike the case of zero
order, the $\nu$-order FW has its intensity pattern distributed over a
cylindrical surface with a radius given by

$\rho_{\nu}\approx\frac{\zeta_{\nu}}{\left(k^{2}-Q^{2}\right)^{1/2}}\text{,}$
(9)

where $\zeta_{\nu}$ represents the maximum value of $J_{\nu}(\zeta)$ when
$\zeta=\zeta_{\nu}$. This is illustrated in Fig. (1b), which shows an example
of a second-order FW.

It has been demonstrated in Ref. [25] that applying the Topological Charge
Raising Operator (TCRO) $\nu$ times to the zero-order Frozen Wave, as given by
Eq.6, yields a FW beam of order $\nu$, whose solution is described by the
following exact and analytical expression:

$\begin{split}\Psi_{\nu}=(-1)^{\nu}N_{\nu}k^{\nu}e^{i\nu\phi}\sum_{n=-\infty}^{\infty}F\left(-\frac{2n\pi}{K}\right)\Bigg{(}\frac{\rho}{\sqrt{\rho^{2}+(z+\frac{n\pi}{k})^{2}}}\Bigg{)}^{\nu}j_{\nu}(\xi)\,\,,\end{split}$
(10)

where $N_{\nu}=-(k^{2}-Q^{2})^{-\nu/2}/(J_{\nu})_{max}$ is a normalization
constant222The normalization constant $N_{\nu}$ does not emerge naturally in
the TCRO application process. We can understand the necessity of this constant
as a multiplicative factor in the final solution when we realize that by
applying the TCRO $\nu$ times to the integral solution (1) of the zero-order
beam, two important things occur. The first is the emergence of the factor
$(-1)^{\nu}(k^{2}-\beta^{2})^{\nu/2}$ multiplying the original spectrum
$S(\beta)$, and the second is the substitution of the zero-order Bessel
function, $J_{0}(\sqrt{k^{2}-\beta^{2}}\rho)$, with the $\nu$-order Bessel
function, $J_{\nu}(\sqrt{k^{2}-\beta^{2}}\rho)$ (along with the appearance of
the term $\exp(i\nu\phi)$). This causes the resulting beam, although generally
having a longitudinal pattern with a shape similar to that of the zero-order
beam, to have amplitude values differing by a constant. It can be shown that
this constant is approximately given by $N_{\nu}$, with $(J_{\nu})_{max}$
being the maximum value of $J_{\nu}(\cdot)$, and $j_{\nu}(\cdot)$ is the
spherical Bessel function of the first kind.

This technique, of obtaining a higher-order Frozen Wave from a zero-order one
through successive applications of the TCRO, has enabled the creation of a
broader range of structured micrometer beams endowed with orbital angular
momentum, which can be utilized in a variety of applications.

As previously mentioned, in general, the longitudinal intensity pattern of the
$\nu$-order Frozen Wave, given by Eq. (10), will exhibit the same Longitudinal
Intensity Pattern (LIP) as the original zero-order Frozen Wave, with the
intensity shifted from $\rho=0$ to a cylindrical surface of radius
$\rho_{\nu}$. However, it is important to highlight that the spectrum of the
resulting $\nu$-order beam will not be identical to the spectrum $S(\beta)$ of
the original zero-order beam. Instead, it will be given by
$(-1)^{\nu}\left(k^{2}-\beta^{2}\right)^{\nu/2}S(\beta)$ (for more details,
refer to Ref. [25]) and, consequently, can be considerably different from
$S(\beta)$, potentially resulting in undesirable distortions in the
Longitudinal Intensity Pattern (LIP) of the higher-order beam. When the
original Frozen Wave (zero-order Frozen Wave) has a wide spectrum, a
characteristic feature for beams structured in micrometer-scale regions, this
distortion effect in the higher-order Frozen Wave given by Eq.10 can occur
even for small values of $\nu$.

To address the aforementioned limitations associated with the TCRO method, we
will introduce in Section 2.3 an exact formulation, referred to here as the
expansion method, designed to provide a higher-order Frozen Wave (FW) with the
same spectrum as the original zero-order one, making it immune to distortions
that may occur when successively applying the TCRO to the initial zero-order
beam.

####  2.2.1 Example

In all examples we consider the vacuum and used the following parameters:
angular frequency $\omega=3.54\times 10^{-15}$ rad/s, implying
$\lambda_{0}=0.532\,\mu\mathrm{m}$ for the wavelength, and
$Q=0.85\frac{\omega}{c}$. Considering a morphological function $f(z)$ given by
a supergaussian of order 14 centered at $z_{1}=0\,\mu\mathrm{m}$ and width
$2*Z_{1}=20\,\mathrm{\mu m}$, the chosen function $F(z)$ is:

$F(z)=\mathrm{exp}\left(-\left(\frac{z-z_{1}}{Z_{1}}\right)^{14}\right)\mathrm{exp}(iQz).$
(11)

This example illustrates the application of the TCRO to obtain a second-order
micrometer Frozen Wave (FW). Using Eq. (6), it is possible to construct a
zero-order FW with the desired Longitudinal Intensity Pattern (LIP) given by
the modulus square of 11 and a spot radius of approximately $\rho_{0}\approx
0.18,\mu\mathrm{m}$, as shown in Fig. (1a). Applying the TCRO twice to this
zero-order FW yields the solution 10 with $\nu=2$, and its intensity is
presented in Fig. (1b). It is evident that the LIP is preserved and
concentrated over a cylindrical surface with a radius of $\rho_{2}\approx
0.49\,\mu\mathrm{m}$.

Figure 1: Longitudinal intensity patterns of (a) a zero order FW and (b)
second order FW obtained from the TCRO method.

###  2.3 The Expansion Method

####  2.3.1 Scalar case

Our main goal here is to obtain the analytical and exact solution for the
following continuous superposition of Bessel beams of order $\nu\neq 0$:

$\Psi(\rho,\phi,z)=\emph{N}\,e^{i\nu\phi}\int_{-\omega/c}^{+\omega/c}S(\beta)J_{\nu}(\eta(\beta)\rho)e^{i\beta
z}d\beta,$ (12)

where $\emph{N}=1/(J_{\nu})_{max}$ and $S(\beta)$ is given by (4,5). With
this, we will have a beam of topological charge $\nu$, endowed with the
property $|\Psi(\rho=\rho_{\nu},\phi,z)|\approx|F(z)|^{2}$, where $F(z)$ is
given by Eq.(3), and $\rho_{\nu}$ by Eq.(9). In other words, we will obtain a
structured hollow beam with a transverse radius of $\rho_{\nu}$ and a LIP
chosen at will, without any risk of obtaining a beam with distortions in the
chosen spatial pattern and alterations in the spectrum, as observed in the
TCRO method.

Unfortunately, a solution for the integral in (12) with $S(\beta)$ given by
(4) is not known. Now, using an ingenious mathematical trick, we will obtain
the desired solution. To do so, let’s write the integral solution (12) as

$\Psi(\rho,\phi,z)=\emph{N}\,e^{i\nu\phi}\int_{-\omega/c}^{+\omega/c}\frac{\Lambda(\beta)}{\eta(\beta)}J_{\nu}(\eta(\beta)\rho)e^{i\beta
z}d\beta,$ (13)

where

$\Lambda(\beta)=\eta(\beta)S(\beta)\,\,,$ (14)

with $\eta(\beta)$ given by Eq.(2).

Expanding $\Lambda(\beta)$ into a Fourier series, we have

$\Lambda(\beta)=\sum_{n=-\infty}^{\infty}B_{n}\exp\left(i\frac{2\pi
n}{K}\beta\right),$ (15)

with the coefficients $B_{n}$ given as

$B_{n}=\frac{1}{K}\int_{-\omega/c}^{+\omega/c}\eta(\beta)S(\beta)e^{-i\frac{2\pi
n}{K}\beta}d\beta.$ (16)

Now, by using Eq.(15) into Eq.(13) and using the result [35]

$\begin{split}&\int_{0}^{a}\frac{\cos(cx)}{\sqrt{a^{2}-x^{2}}}J_{\nu}(b\sqrt{a^{2}-x^{2}})dx\\\
&=\frac{\pi}{2}J_{\frac{\nu}{2}}\left[\frac{a}{2}(\sqrt{b^{2}+c^{2}}-c)\right]J_{\frac{\nu}{2}}\left[\frac{a}{2}(\sqrt{b^{2}+c^{2}}+c)\right],\end{split}$
(17)

we find that the exact analytical solution for Eq.(12) is given by

$\begin{split}&\Psi(\rho,\phi,z)=\emph{N}\,\pi e^{i\nu\phi}\\\
&\times\sum_{n}B_{n}J_{\frac{\nu}{2}}\left[\frac{\omega}{2c}\left(\sqrt{\rho^{2}+\left(z+\frac{2\pi
n}{K}\right)^{2}}-\left(z+\frac{2\pi n}{K}\right)\right)\right]\\\ &\times
J_{\frac{\nu}{2}}\left[\frac{\omega}{2c}\left(\sqrt{\rho^{2}+\left(z+\frac{2\pi
n}{K}\right)^{2}}+\left(z+\frac{2\pi
n}{K}\right)\right)\right]\,\,,\end{split}$ (18)

In order for the solution (18) to be fully defined, it is necessary to know
the values of the coefficients $B_{n}$ given Eq.(16), whose integral could
naturally be solved numerically. However, here we will obtain the analytical
solution for it.

We proceed by expanding $S(\beta)$ and $\eta(\beta)$ into Fourier series as
follows:

$S(\beta)=\sum_{m=-\infty}^{\infty}A_{m}\exp\left(i\frac{2\pi
m}{K}\beta\right),$ (19)

where, as before,

$A_{m}=\frac{1}{K}F\left(\frac{-2\pi m}{K}\right),$ (20)

and

$\eta(\beta)=\sum_{u=-\infty}^{\infty}C_{u}\exp\left(i\frac{2\pi
u}{K}\beta\right),$ (21)

with the coefficients $C_{u}$ determined by

$C_{u}=\frac{1}{K}\int_{-\omega/c}^{+\omega/c}\eta(\beta)e^{-i\frac{2\pi
u}{K}\beta}d\beta.$ (22)

By writing $\eta(\beta)$ as $\frac{\omega}{c}\sqrt{1-x^{2}}$ where
$x^{2}=\frac{\beta^{2}}{(\omega^{2}/c^{2})}$, we have that

$\begin{split}C_{u}=\frac{k}{2}\int_{-1}^{1}\sqrt{1-x^{2}}e^{-iu\pi x}dx\
&=-k\sqrt{\pi}\left(\frac{1}{u\pi}\right)\Gamma\left(\frac{3}{2}\right)J_{1}(-u\pi)\,\,,\end{split}$
(23)

where $\Gamma(\cdot)$ is the gamma function.

Finally, by using Eq. (19) into Eq. (16), we have

$B_{n}=\sum_{m=-\infty}^{\infty}A_{m}\frac{1}{K}\int_{-\omega/c}^{+\omega/c}\eta(\beta)e^{-i\frac{2\pi(n-m)}{K}\beta}d\beta,$
(24)

and so

$B_{n}=\sum_{m}A_{m}C_{n-m}.$ (25)

With this, using Eqs.(20,23), the coefficients $B_{n}$ are given by

$B_{n}=\sum_{m=-\infty}^{\infty}\,-\frac{1}{2\sqrt{\pi}(n-m)}\Gamma\left(\frac{3}{2}\right)F\left(\frac{-2\pi
m}{K}\right)J_{1}(-(n-m)\pi)\,.$ (26)

To summarize, a Frozen Wave of order $\nu\neq 0$, obtained by the continuous
superposition of Bessel beams through Eq.(12) with $S(\beta)$ given by Eqs.
(4,5), is described by the exact analytical solution provided in Eq. (18),
where the coefficients $B_{n}$ are determined by Eq.(26).

####  2.3.2 Vectorial case

Micrometer optical Frozen Waves are shaped in such small dimensions
(nonparaxial regime) that a scalar analysis of the phenomenon is not
appropriate. Taking this into account, we will adopt a vectorial approach
here, considering an optical beam whose electric field is given by

$\textbf{E}=E_{0}\Psi(\rho,\phi,z)\mathbf{\hat{x}}+E_{z}\mathbf{\hat{z}},$
(27)

Here, the spatial modeling is performed on the transverse component of the
field, $E_{x}=E_{0}\Psi(\rho,\phi,z)$, with $\Psi(\rho,\phi,z)$ given by the
solution (18) from Section 2.3.1, and $E_{0}$ is a constant with a value of
one and units of electric field. The axial component $E_{z}$ is then obtained
from Gauss’s law:

$E_{z}=-E_{0}\frac{\partial}{\partial x}\int\Psi(\rho,\phi,z)dz,$ (28)

where $\frac{\partial}{\partial x}$ will be written later in cylindrical
coordinates. After inserting Eq. (18) into Eq. (28), it is evident that the
integral with respect to $z$ is challenging to perform analytically. To
address this issue, we use a strategy involving the spectral representation of
$\Psi(\rho,\phi,z)$ provided by Eq. (13). This approach simplifies the
integration, as the only term depending on $z$ in this expression is
$\mathrm{e}^{i\beta z}$. Consequently, after using (13) into (28), we get

$E_{z}=-\emph{N}\,E_{0}\frac{\partial}{\partial
x}\left[e^{i\nu\phi}\int_{-\omega/c}^{+\omega/c}\frac{\Lambda(\beta)}{i\beta\eta(\beta)}J_{\nu}(\eta(\beta)\rho)e^{i\beta
z}d\beta\right]$ (29)

with $\Lambda(\beta)$ given by Eq. (14) and $\eta(\beta)$ by Eq.(2).

To start with the evaluation of the integral appearing in Eq.(29), let us
express $\eta(\beta)S(\beta)/i\beta$ in Fourier series

$\eta(\beta)\frac{S(\beta)}{i\beta}=\sum_{n=-\infty}^{\infty}B^{\prime}_{n}\exp\left(i\frac{2\pi
n}{K}\beta\right)\,\,,$ (30)

with the coefficients $B^{\prime}_{n}$ given by

$B_{n}^{\prime}=\frac{1}{K}\int_{-\omega/c}^{+\omega/c}\eta(\beta)\frac{S(\beta)}{i\beta}e^{-i\frac{2\pi
n}{K}\beta}d\beta.$ (31)

By using Eq.(30) into Eq.(29), and using the result given by Eq.(17), we
obtain

$\begin{split}&E_{z}(\rho,\phi,z)=\left(\cos\phi\frac{\partial}{\partial\rho}-\frac{\sin\phi}{\rho}\frac{\partial}{\partial\phi}\right)\left\\{N\,\pi\mathrm{e}^{i\nu\phi}\times\right.\\\
&\left.\sum_{n}B_{n}^{\prime}J_{\frac{\nu}{2}}\left[\frac{\omega}{2c}\left(\sqrt{\rho^{2}+\left(z+\frac{2\pi
n}{K}\right)^{2}}-\left(z+\frac{2\pi n}{K}\right)\right)\right]\right.\\\
&\left.\times
J_{\frac{\nu}{2}}\left[\frac{\omega}{2c}\left(\sqrt{\rho^{2}+\left(z+\frac{2\pi
n}{K}\right)^{2}}+\left(z+\frac{2\pi
n}{K}\right)\right)\right]\right\\}\,\,,\end{split}$ (32)

where we have written $\frac{\partial}{\partial x}$ in cylindrical
coordinates.

Equation (32) is the exact analytical solution for the axial component of the
electric field of the FW beam, but for it to be complete, it is necessary to
know the values of the coefficients $B^{\prime}_{n}$ given by Eq.(31), whose
integral has to be solved. To do this, we will use the that

$\frac{S(\beta)}{i\beta}=\sum_{m=-\infty}^{\infty}A_{m}^{\prime}\exp\left(i\frac{2\pi
m}{K}\beta\right),$ (33)

with

$A_{m}^{\prime}=\frac{1}{K}g\left(\frac{-2\pi m}{K}\right),$ (34)

and

$g\left(\frac{-2\pi m}{K}\right)=\left[\int
F(z)dz\right]_{z=-\frac{2n\pi}{K}}.$ (35)

Now, by using Eq. (33) into Eq. (31) we have

$B_{n}^{\prime}=\sum_{m=-\infty}^{\infty}A_{m}^{\prime}\frac{1}{K}\int_{-\omega/c}^{+\omega/c}\eta(\beta)e^{-i\frac{2\pi(n-m)}{K}\beta}d\beta,$
(36)

then

$B_{n}^{\prime}=\sum_{m}A_{m}^{\prime}C_{n-m}.$ (37)

where $C_{n-m}$ can be obtained from Eq. (23). Therefore, the final expression
of $B_{n}^{\prime}$ is given by

$B_{n}^{\prime}=\sum_{m}\,-\frac{1}{2\sqrt{\pi}(n-m)}\Gamma\left(\frac{3}{2}\right)g\left(\frac{-2\pi
m}{K}\right)J_{1}(-(n-m)\pi)\,\,,$ (38)

with $g\left(\frac{-2\pi m}{K}\right)$ given by Eq.(35).

####  2.3.3 Example

Example 1. Here, we construct two scalar sixth-order FWs with the LIP dictated
by the function given by Eq.(11), where the first one, Fig. (2a), is obtained
using the TCRO method, and the second one, Fig.(2b), using the expansion
method. In both cases, we use the same wavelength and parameter $Q$ values as
in the example from Section 2.3.1.

It is evident that the FW obtained with the TCRO method exhibits distortions
at the edges. This occurs because, as previously mentioned, the TCRO method
can reproduce the intended LIP within certain limitations. In this case, the
six successive applications (since $\nu=6$) of the TCRO on the zero-order FW
are sufficient to considerably alter the spectrum of the resulting beam,
leading to the distortion in the desired LIP.

On the other hand, we observe that the FW obtained by the expansion method is
immune to such distortions, providing the desired LIP. This is due to the fact
that the spectrum of the resulting sixth-order beam is, in fact, the spectrum
$S(\beta)$.

Figure 2: Intensities of two 6th-order scalar FWs designed to have the same
LIP, with: (a) FW obtined from the TCRO method; (b) FW obtained from the
expansion method.

Example 2. In this example, we apply the expansion method considering the
vectorial nature of the beam, more specifically considering a beam with an
electric field given by Eq. (27), with $\Psi$ and $E_{z}$ given by Eqs. (18)
and (32), respectively, where the transverse component $E_{x}=\Psi$ is the
same sixth-order FW from the previous example obtained by the expansion
method. Figures (3a) and (3b) show, respectively, the 3D intensities of the
transverse and axial components of the electric field, with the inset figures
displaying the intensity fields on the $z=0$ plane. As mentioned earlier, for
highly non-paraxial beams, a longitudinal component emerges, as observed in
Figure (3b) and, in some cases, it can present a non-negligible part in the
composition of the resulting beam.

Figure 3: a) 3D intensity of the transverse electric field of the FW beam
obtained from the expansion method; b) 3D intensity of the axial component
$E_{z}$ of the electric field of the FW beam obtained from the expansion
method. Insets: intensity fields on the $z=0$ plane.

##  3 Conclusions

The continuous superposition of Bessel beams is a powerful tool for spatial
modeling of the optical field in micrometer regions. This methodology, known
as the Frozen Wave method, is well-developed for the zero-order case, i.e.,
the superposition of zero-order Bessel beams, providing analytical and exact
solutions for structured light beams with zero topological charge. However,
applying the method for higher orders, aiming to obtain analytical and exact
solutions that provide structured beams with orbital angular momentum, is
mathematically more complex. The only known method to mitigate such
difficulties was the use of a Topological Charge Raising Operator (TCRO). By
applying this operator $\nu$ times to a structured zero-order beam, it can
produce a structured beam of order (topological charge) $\nu$. However, the
successive application of the TCRO affects the spectrum of the resulting beam,
leading to distortion in the desired spatial pattern.

Addressing this issue, we have developed what we refer to as the ”expansion
method,” which definitively resolves the challenges of the Frozen Wave method
for higher orders. It provides an analytical and exact approach to obtaining
structured beams of arbitrary order (topological charge) in micrometer
regions. Furthermore, considering that structured beams in micrometer regions
are highly non-paraxial, we extended the expansion method to the vectorial
case, taking into account the axial component of the electric field. This
method is particularly relevant for research in atom guidance, optical
trapping, micro-manipulation of particles, microlithography, imaging, and
medical applications.

## References

  * [1] J. Durnin, “Exact solutions for nondiffracting beams. i. the scalar theory,” J. Opt. Soc. Am. A 4, 651–654 (1987).
  * [2] M. Zamboni-Rached, E. Recami, and M. Balma, “Simple and effective method for the analytic description of important optical beams when truncated by finite apertures,” Appl. Opt. 51, 3370–3379 (2012).
  * [3] F. Gori, G. Guattari, and C. Padovani, “Bessel-gauss beams,” Optics Communications 64, 491–495 (1987).
  * [4] K. Huang, P. Shi, G. W. Cao, K. Li, X. B. Zhang, and Y. P. Li, “Vector-vortex bessel–gauss beams and their tightly focusing properties,” Opt. Lett. 36, 888–890 (2011).
  * [5] H. Hernandez-Figueroa, M. Zamboni-Rached, and E. Recami, _Localized Waves (book of 387 pages)_ (2008).
  * [6] H. Hernandez-Figueroa, E. Recami, and M. Zamboni-Rached, _Non-Diffracting Waves(book of 500 pages)_ (J. Wiley, 2014).
  * [7] P. Saari and K. Reivelt, “Evidence of $\mathit{X}$-shaped propagation-invariant localized light waves,” Phys. Rev. Lett. 79, 4135–4138 (1997).
  * [8] E. Recami, “On localized “x-shaped” superluminal solutions to maxwell equations11work partially supported by infn, murst and cnr (italy), and by capes, cnpq (brazil). this paper first appeared in preliminary form as report infn/fm–96/01 (i.n.f.n.; frascati, 1996), its number in the electronic lanl archives being # physics/9610012; oct. 15, 1996; cf. ref. [24].” Physica A: Statistical Mechanics and its Applications 252, 586–610 (1998).
  * [9] K. Huang, P. Shi, G. W. Cao, K. Li, X. B. Zhang, and Y. P. Li, “Vector-vortex bessel–gauss beams and their tightly focusing properties,” Opt. Lett. 36, 888–890 (2011).
  * [10] G. A. Siviloglou and D. N. Christodoulides, “Accelerating finite energy airy beams,” Opt. Lett. 32, 979–981 (2007).
  * [11] M. Zamboni-Rached, “Stationary optical wave fields with arbitrary longitudinal shape by superposing equal frequency bessel beams: Frozen waves,” Opt. Express 12, 4001–4006 (2004).
  * [12] M. Zamboni-Rached, E. Recami, and H. E. Hernández-Figueroa, “Theory of “frozen waves”: modeling the shape of stationary wave fields,” J. Opt. Soc. Am. A 22, 2465–2475 (2005).
  * [13] T. A. Vieira, M. R. R. Gesualdi, and M. Zamboni-Rached, “Frozen waves: experimental generation,” Opt. Lett. 37, 2034–2036 (2012).
  * [14] R. A. B. Suarez, L. A. Ambrosio, A. A. R. Neves, M. Zamboni-Rached, and M. R. R. Gesualdi, “Experimental optical trapping with frozen waves,” Opt. Lett. 45, 2514–2517 (2020).
  * [15] M. Zamboni-Rached, E. Recami, T. A. Vieira, M. R. Gesualdi, and J. Nobre-Pereira, “Structured light by linking diffraction-resistant spatially shaped beams,” Phys. Rev. Appl. 10, 034023 (2018).
  * [16] J. A. V. Mendonça, M. Zamboni-Rached, and E. Recami, “Arrays of frozen waves: Some theory and experiments,” Optics Communications 482, 126576 (2021).
  * [17] M. Zamboni-Rached, “Diffraction-attenuation resistant beams in absorbing media,” Opt. Express 14, 1804–1809 (2006).
  * [18] M. Zamboni-Rached and M. Mojahedi, “Shaping finite-energy diffraction- and attenuation-resistant beams through bessel-gauss–beam superposition,” Phys. Rev. A 92, 043839 (2015).
  * [19] G. de A. Lourenço-Vittorino and M. Zamboni-Rached, “Modeling the longitudinal intensity pattern of diffraction resistant beams in stratified media,” Appl. Opt. 57, 5643–5651 (2018).
  * [20] J. ao Quaglio and M. Zamboni-Rached, “Transmission of spatial-shaped diffraction-resistant beams through stratified dielectric media: finite energy formulation,” J. Opt. Soc. Am. A 35, 2079–2086 (2018).
  * [21] M. Corato-Zanarella and M. Zamboni-Rached, “Electromagnetic frozen waves with radial, azimuthal, linear, circular, and elliptical polarizations,” Phys. Rev. A 94, 053802 (2016).
  * [22] M. Corato-Zanarella, A. H. Dorrah, M. Zamboni-Rached, and M. Mojahedi, “Arbitrary control of polarization and intensity profiles of diffraction-attenuation-resistant beams along the propagation direction,” Phys. Rev. Appl. 9, 024013 (2018).
  * [23] A. H. Dorrah, M. Zamboni-Rached, and M. Mojahedi, “Controlling the topological charge of twisted light beams with propagation,” Phys. Rev. A 93, 063864 (2016).
  * [24] Q. Liang, D. Yang, Y. Zhang, Y. Zheng, and L. Hu, “Probability of orbital angular momentum modes carried by a finite energy frozen wave in turbulent seawater,” OSA Continuum 3, 2429–2440 (2020).
  * [25] J. Nobre-Pereira, V. de Angelis, L. Ambrosio, and M. Zamboni-Rached, “Exact analytical solutions for micrometer structured vortex beams with applications to optical tweezers,” Optics & Laser Technology 171, 110374 (2024).
  * [26] L. A. Ambrosio and M. Zamboni-Rached, “Analytical approach of ordinary frozen waves for optical trapping and micromanipulation,” Applied Optics 54, 2584 – 2593 (2015). Cited by: 39.
  * [27] E. G. P. Pachon, M. Zamboni-Rached, A. H. Dorrah, M. Mojahedi, M. R. R. Gesualdi, and G. G. Cabrera, “Architecting new diffraction-resistant light structures and their possible applications in atom guidance,” Opt. Express 24, 25403–25408 (2016).
  * [28] A. H. Dorrah, P. Bordoloi, V. S. de Angelis, J. O. de Sarro, M. Zamboni-Rached, L. A. Ambrosio, and F. Capasso, “Longitudinally structured light sheets for 3d holography,” in _2022 Conference on Lasers and Electro-Optics (CLEO),_ (2022), pp. 1–2.
  * [29] R. M. Herman and T. A. Wiggins, “Production and uses of diffractionless beams,” J. Opt. Soc. Am. A 8, 932–942 (1991).
  * [30] A. Marcinkevi\vcius, S. Juodkazis, S. Matsuo, V. Mizeikis, and H. Misawa, “Application of bessel beams for microfabrication of dielectrics by femtosecond laser,” Japanese Journal of Applied Physics 40, L1197–L1199 (2001).
  * [31] J. Arlt, V. Garces-Chavez, W. Sibbett, and K. Dholakia, “Optical micromanipulation using a bessel light beam,” Optics Communications 197, 239–245 (2001).
  * [32] M. Zamboni-Rached, L. A. Ambrosio, A. H. Dorrah, and M. Mojahedi, “Structuring light under different polarization states within micrometer domains: exact analysis from the maxwell equations,” Optics Express 25, 10051–10056 (2017).
  * [33] M. Zamboni-Rached and E. Recami, “Subluminal wave bullets: Exact localized subluminal solutions to the wave equations,” Phys. Rev. A 77, 033824 (2008).
  * [34] L. Mackinnon, “A nondispersive de broglie wave packet,” Foundations of Physics 8, 157–176 (1978).
  * [35] “Table of integrals, series, and products,” in _Table of Integrals, Series, and Products (Seventh Edition),_ D. Zwillinger, V. Moll, I. Gradshteyn, and I. Ryzhik, eds. (Academic Press, Boston, 2007), p. iv, seventh edition ed.