# Completely hyperexpansive operators with finite rank defect operator and de
Branges-Rovnyak spaces

Saee A. Joshi  Ramnarain Ruia Autonomous College
Mumbai- 400019, India<EMAIL_ADDRESS>and  Vinayak. M. Sholapurkar
Bhaskaracharya Pratishthana
Pune- 411004, India<EMAIL_ADDRESS>
###### Abstract.

The process of identifying a Dirichlet-type space $D(\mu)$ for a positive,
Borel measure $\mu$, supported on the unit circle $\mathbb{T},$ with a de
Branges-Rovnyak space was initiated by Sarason [sar]. A characterization of
the symbol for a de Branges-Rovnyak spaces for which the shift operator is a
$2$-isometry, was provided in [kz]. In this paper, capitalizing on the
Aleman’s model [ale] for the cyclic, analytic, completely hyperexpansive
operators, we provide a characterization of cyclic, analytic, completely
hyperexpansive operator with finite rank defect operator in terms of the
symbol for a de Branges-Rovnyak space.

###### Key words and phrases:

de Branges-Rovnyak space, Dirichlet-type space, completely hyperexpansive
operators

###### 2020 Mathematics Subject Classification:

Primary 47B20, 47B32; Secondary 47A67, 46E22

##  1\. Introduction and Preliminaries

It is well known that the class of $2$-isometric operators is a subclass of
the class of completely hyperexpansive operators. A theme of the present paper
hinges around extending the results in [kz] for $2$-isometric operators to
completely hyperexpansive operators, especially in the case when the rank of
the defect operators is finite.

In this paper, we describe a model for a cyclic, analytic, completely
hyperexpansive operator $T$, if the defect operator $\Delta_{T}=T^{*}T-I$ is
of finite rank. A special case of such a model for $2$-isometric operators,
has been obtained in [sc]. This particular case has a bearing on the classical
model for a cyclic, analytic, $2$-isometric operator described by Richter (see
[rich]). Here, we capitalize on the classical model for a cyclic, analytic,
completely hyperexpansive operator given by Aleman (see [ale]) and impose an
additional condition that the defect operator is of finite rank. In both the
classical representations, the corresponding operator is shown to be unitarily
equivalent to the multiplication operator $M_{z}$ on a Dirichlet space
$D(\mu)$ where the measure is supported on the closed unit disk. In
particular, if the operator is $2$-isometric, then the support of $\mu$ lies
inside the unit circle. In the circumstances where the defect operator is of
finite rank, the de Branges-Rovnyak spaces make their appearance. Indeed, in a
series of papers [sc, cgr, kz, sar], the authors have established interesting
connections between Dirichlet spaces and de Branges-Rovnyak spaces. The most
pertinent connections in the context of present paper are [kz, Theorem 1 and
Corollary 3]. We record these results for a ready reference.

######  Theorem 1.1.

Suppose that $b$ is non-extreme in the unit ball of $\mathbb{H}^{\infty}$.
Then the restriction of the shift operator $S|_{\mathcal{H}(b)}$ is concave if
and only if

$b(z)=\frac{c+\gamma z}{1-\beta z},z\in\mathbb{D}\ \mbox{where}\
c,\gamma,\beta\in\mathbb{C}\ \mbox{with}\ |\beta|<1.$

Furthermore, $S|_{\mathcal{H}(b)}$ is $2$-isometry if and only if
$b(z)=\frac{c+\gamma z}{1-\beta z},z\in\mathbb{D}$ where
$c,\gamma,\beta\in\mathbb{C}$ are such that
$1+|\beta|^{2}-|c|^{2}-|\gamma|^{2}=|\beta+\overline{c}\gamma|$.

######  Corollary 1.2 (Chevrot, Guillot, Ransford [cgr]).

Let $\mu$ be a finite positive measure on $\mathbb{T}$ and let $b$ be in the
unit ball of $\mathbb{H}^{\infty}$. Then $D(\mu)=\mathcal{H}(b)$, with
equality of norms, if and only if

$\mu=c\delta_{\lambda}\ \mbox{and}\ b(z)=\frac{\gamma z}{1-\beta z},$

where $\beta\neq 0,\ |\beta|+|\gamma|=1,\ c=\frac{|\gamma|^{2}}{|\beta|}$ and
$\lambda=\frac{\overline{\beta}}{|\beta|}$.

In this paper, we obtain an analogue of Theorem 1.1 for completely
hyperexpansive operators. Further, we also obtain a version of Corollary 1.2
for the Dirichlet space $D(\mu)$, where the measure $\mu$ is supported on the
closed unit disk. In the process, we also prove a generalization of [sc,
Theorem 6.1], extending the result to a completely hyperexpansive operator.
The main result proved in this paper is stated below:

######  Theorem 1.3.

Let $\mu$ be a finite, positive, Borel measure supported on
$\overline{\mathbb{D}}$ and let $b\in\mathbb{H}^{\infty}$ be nonextreme with
$\|b\|_{\infty}\leq 1$. Then $D(\mu)=\mathcal{H}(b)$ with equality of norms,
if and only if $\mu=|\alpha|^{2}\delta_{\lambda},\ \mbox{for some
}\alpha\in\mathbb{C},\ \lambda\in\overline{\mathbb{D}}$ and
$b(z)=\frac{A_{\lambda}z}{1-B_{\lambda}z}$ where $A_{\lambda}$ and
$B_{\lambda}$ are constants depending on $\alpha$ and $\lambda$ such that

$|A_{\lambda}|^{2}=\begin{cases}\frac{|\alpha|^{2}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})-\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right),&\mbox{
if}\ \lambda\neq 0\\\ \displaystyle\frac{1}{1+|\alpha|^{2}},&\mbox{ if}\
\lambda=0\end{cases}$

and

$B_{\lambda}=\frac{|A_{\lambda}|^{2}\overline{\lambda}}{|\alpha|^{2}}.$

We obtain Corollary 1.2 as a special case of Theorem 1.3.
In the remainder of this section, we set up the notation as well as the
context necessary in the sequel. In section 2, we prove that the defect
operator of a multiplication operator $M_{z}$ on a Dirichlet space $D(\mu)$,
where $\mu$ is supported on the closed unit disk, is of finite rank if and
only if $\mu$ is finitely supported on the closed unit disk. In section 3, we
prove a generalization of Theorem 1.1 and finally, capitalizing on the results
in sections 2 and 3, we shall prove the main theorem in section 4.
For a complex Hilbert space $\cal H$, let $\cal B(\cal H)$ denote the $C^{*}$
algebra of bounded linear operators on $\cal H.$ An operator $T\in\cal B(\cal
H)$ is called cyclic if there exists a vector $f\in\cal H$ such that
$\vee\\{T^{n}f:n\geq 0\\}=\cal H.$ We say that $T$ is analytic if $\cap_{n\geq
0}T^{n}\cal H=\\{0\\}.$ An operator $T\in\cal B(\cal H)$ is a $2$-isometry if

$I-2T^{*}T+T^{*^{2}}T^{2}=0.$

An operator $T\in\cal B(\cal H)$ is a $k$-hyperexpansive if

$\sum_{i=0}^{n}(-1)^{i}{n\choose i}T^{*^{i}}T^{i}\leq 0\ ,\ \forall\ 1\leq
n\leq k.$

An operator $T\in\cal B(\cal H)$ is completely hyperexpansive if

$\sum_{i=0}^{n}(-1)^{i}{n\choose i}T^{*^{i}}T^{i}\leq 0\ ,\ \forall n\geq 1.$

For an excellent account of the basic properties of $2$\- isometries, the
reader is referred to [rich, rich1, rich2, rich3]. The theory of completely
hyperexpansive operators has been developed in [ale, aa1, vms]. In what
follows, we need to recall the definitions of two special types of Hilbert
spaces of holomorphic function viz. the Dirichlet space and the de Branges-
Rovnyak space. We shall denote the open unit disk in the complex plane by
$\mathbb{D}$ and the unit circle by $\mathbb{T}.$ Let $\mathbb{H}^{2}$ denote
the classical Hardy space of analytic functions on $\mathbb{D}$ having square
integrable Taylor coefficients and $\mathbb{H}^{\infty}$ denote the space of
all bounded analytic functions on $\mathbb{D}$.

We shall be dealing with the Dirichlet type space $D(\mu)$ where the measure
$\mu$ is supported on the closed unit disk $\overline{\mathbb{D}}.$ Aleman
points out the delicacy in defining the local Dirichlet integrals for the
points $\zeta\in\mathbb{T}$ (see [ale, Proposition 1.7, Chapter IV]). Indeed,
if $\zeta\in\mathbb{D},$ then the local Dirichlet integral of a function
$f\in\mathcal{H}^{2}$ is defined as

$D_{\zeta}(f)=\left\|\frac{f-f(\zeta)}{z-\zeta}\right\|_{\mathbb{H}^{2}}^{2}=\frac{1}{2\pi}\int_{0}^{2\pi}\left|\frac{f(e^{it})-f(\zeta)}{e^{it}-\zeta}\right|^{2}dt.$

If $\zeta\in\mathbb{T},$ then the local Dirichlet integral is defined as

$D_{\zeta}(f)=\left\|\frac{f-f(\zeta)}{z-\zeta}\right\|_{\mathbb{H}^{2}}^{2}$

if the nontangential limit $f(\zeta)$ of $f$ at $\zeta$ exists and
$D_{\zeta}(f)=\infty$ otherwise.

As pointed out in [ale], the necessity of defining the local Dirichlet
integrals for the boundary points, and thereby extending the support of
measure to the closed unit disk is forced by the operator theoretic
considerations. Indeed, with this type of Dirichlet space, one can incorporate
the case of multiplication operator $M_{z}$ being a $2$-isometry (see [ale,
Theorem 1.10, Chapter IV]). Following [ale, Definition 1.8], we recall the
definition of Dirichlet space under consideration for a ready reference.
Let $\mu$ be a finite, non-negative measure supported on
$\overline{\mathbb{D}}.$ The Dirichlet space $D(\mu)$, consists of those
functions $f\in\mathbb{H}^{2}$ with the property that $\int
D_{\zeta}(f)d\mu(\zeta)<\infty.$ The Dirichlet space $D(\mu)$ is a reproducing
kernel Hilbert space with

$\|f\|^{2}=\|f\|_{\mathbb{H}^{2}}^{2}+\int D_{\zeta}(f)\mathrm{d}\mu(\zeta).$

Further, as $\mu$ is finite, the multiplication operator $M_{z}$ is bounded.
In fact, it turns out that $M_{z}$ on $D(\mu)$ is completely hyperexpansive
and if $T\in\cal B(\cal H)$ is a completely hyperexpansive operator, then
there exists a unique nonnegative measure $\mu$ supported on the closed unit
disk $\overline{\mathbb{D}}$ such that $\cal H=D(\mu)$ with the equality of
norms and $T$ is unitarily equivalent to the operator $M_{z}$ on $D(\mu)$ (see
[ale, Theorem 2.5, Chapter IV]). Further, $M_{z}$ is a $2$-isometry if and
only if $\mu$ is supported on $\mathbb{T}$ (see [ale, Theorem 1.10(ii),
Chapter IV ]).

We now turn our attention to the definition of a de Branges Rovnyak space.

Following [vol2], recall that the de Branges-Rovnyak space $\mathcal{H}(b)$
associated with a symbol $b\in\mathbb{H}^{\infty},\|b\|_{\infty}\leq 1$ is the
image of $\mathbb{H}^{2}$ under the operator
$(I-T_{b}T_{\overline{b}})^{\frac{1}{2}}$ endowed with the inner product

$\langle(I-T_{b}T_{\bar{b}})^{\frac{1}{2}}f\ ,\
(I-T_{b}T_{\bar{b}})^{\frac{1}{2}}g\ \rangle_{b}\ =\ \langle
f,g\rangle_{\mathbb{H}^{2}}\ ,\ f,g\in
ker((I-T_{b}T_{\bar{b}})^{\frac{1}{2}})^{\perp}$

where $\langle\cdot\rangle_{\mathbb{H}^{2}}$ is the usual inner product on
$\mathbb{H}^{2}$ and $T_{b}$ denotes the Toeplitz operator with respect to
$b$.

For a detailed survey of de Branges-Rovnyak spaces, the reader is referred to
[ball, vol1, vol2, sar1].

Throughout this paper, we choose a symbol $b$ to be a nonextreme point in the
unit ball of $\mathbb{H}^{\infty}$. The choice of a nonextreme symbol $b$
ensures following important results which are useful in the sequel:

  1. (1)

The polynomials are dense in $\mathcal{H}(b)$ (see [vol2, Section 23.4])

  2. (2)

The multiplication operator $M_{z}$ on $\mathcal{H}^{2}$ leaves
$\mathcal{H}(b)$ invariant (see [kz, Proposition 4]).

  3. (3)

There exists a unique outer function $a$ such that $a(0)>0$ and
$|a(\zeta)|^{2}+|b(\zeta)|^{2}=1,\ \mbox{a.e on }\mathbb{T}$. Such a function
is called the pythagorean mate of $b$ and $(b,a)$ are said to form a pair.
Further,

$f\in\mathcal{H}(b)\iff T_{\overline{b}}f\in T_{\overline{a}}\mathbb{H}^{2}$

and $T_{\overline{a}}\mathbb{H}^{2}$ is dense in $\mathcal{H}(b)$ (see [sar1,
pg 24-25]).

##  2\. Cyclic, Analytic, Completely Hyperexpansive Operators With Finite Rank
Defect Operator

In this section, we provide a characterization for the measure $\mu$, so that
the multiplication operator $M_{z}$ on the Dirichlet space $D(\mu)$, has
defect operator of finite rank. At this stage, we refer the readers to the
analogous result [sc, Theorem 6.1], proved in the context of a $2$-isometry.
It turns out that a matrix decomposition of a completely hyperexpansive
operator, as described in [vms], turns out to be useful in proving the
characterization.

######  Proposition 2.1.

Let $\mu$ be a finite positive Borel measure supported on
$\overline{\mathbb{D}}$ and let $\Delta_{M_{z}}=M_{z}^{*}M_{z}-I$ be the
defect operator of the multiplication operator on the Dirichlet space
$D(\mu)$. Then the following are equivalent:

  1. (i)

$\Delta_{M_{z}}$ is of rank $k$

  2. (ii)

There exist constants $c_{i}>0,1\leq i\leq k$ and
$\zeta_{i}\in\overline{\mathbb{D}},1\leq i\leq k$ such that
$\mu=\sum_{i=0}^{k}c_{i}\delta_{\zeta_{i}}$

###### Proof.

Suppose that $\Delta_{M_{z}}$ is of rank k.
It is known that $M_{z}$ on $D(\mu)$ is cyclic, analytic and completely
hyperexpansive ([ale]). Using the decomposition theorem for completely
hyperexpansive operators ([vms, Proposition 3.1]), we have a decomposition for
$M_{z}$ as follows:
$M_{z}=\begin{bmatrix}V&E\\\ 0&W\\\ \end{bmatrix}$ defined on
$\ker(\Delta_{M_{z}})\oplus\overline{ran(\Delta_{M_{z}})}$
where

  1. (a)

$V$ is an isometry

  2. (b)

$V^{*}E=0$

  3. (c)

If $D=E^{*}E+W^{*}W-I$, then $D\geq 0$ and $\ker(D)=\\{0\\}$

  4. (d)

Letting $M_{W}(T)=W^{*}TW$ we have $(I-M_{W})^{m-1}(D)\geq 0,\\\ 1\leq m\leq
k$

We note that $\Delta_{M_{z}}=\begin{bmatrix}0&0\\\ 0&D\\\ \end{bmatrix}$ where
$D$ is invertible and self- adjoint.
Now, letting $m=2$ in (d) above, we have

$D\geq W^{*}DW$ (2.1)

Further, as $D$ is invertible, the inequality (2.1) yields,

$\displaystyle(D^{\frac{1}{2}}WD^{-\frac{1}{2}})^{*}(D^{\frac{1}{2}}WD^{-\frac{1}{2}})$
$\displaystyle=$ $\displaystyle D^{-\frac{1}{2}}W^{*}DWD^{-\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle D^{-\frac{1}{2}}DD^{-\frac{1}{2}}\ =\ I$

Therefore , $(D^{\frac{1}{2}}WD^{-\frac{1}{2}})$ is a contraction. Now, let
$B=\\{v_{1},v_{2},\cdots v_{k}\\}$ be an orthonormal basis of
${ran(\Delta_{M_{z}})}$ comprising of eigen vectors of
$D^{\frac{1}{2}}WD^{-\frac{1}{2}}$ corresponding to the eigen values
$\zeta_{i},\ 1\leq i\leq k$, where each $\zeta_{i}\in\overline{\mathbb{D}}$.
Now referring to ([ale, Proposition (1.6)]), we have

$\|zf\|^{2}_{D(\mu)}-\|f\|^{2}_{D(\mu)}=\int_{\overline{\mathbb{D}}}|f|^{2}\
\mathrm{d}\mu\ ,\ \forall f\in D(\mu).$ (2.2)

By polarization identity, one obtains,

$\langle zf,zg\rangle-\langle f,g\rangle=\int_{\overline{\mathbb{D}}}f\bar{g}\
\mathrm{d}\mu\ ,\ \forall f,g\in D(\mu).$ (2.3)

So, for $m,n\geq 0$,

$\displaystyle\langle\Delta_{M_{z}}M_{z}^{n}1,z^{m}\rangle$ $\displaystyle=$
$\displaystyle\langle z^{n+1},z^{m+1}\rangle-\langle z^{n},z^{m}\rangle$
$\displaystyle=$ $\displaystyle\int_{\overline{\mathbb{D}}}z^{n}\bar{z}^{m}\
\mathrm{d}\mu.$

That is

$\langle\Delta_{M_{z}}M_{z}^{n}1,M_{z}^{m}1\rangle=\int_{\overline{\mathbb{D}}}z^{n}\bar{z}^{m}\
\mathrm{d}\mu.$ (2.4)

Now note that,
$\\{D^{-\frac{1}{2}}v_{1},D^{-\frac{1}{2}}v_{2},\cdots,D^{-\frac{1}{2}}v_{k}\\}$
also forms a basis for ${ran(\Delta_{M_{z}})}$.
Let $P$ be the projection onto ${ran(\Delta_{M_{z}})}$. Therefore,
$P1=\sum_{i=1}^{k}\beta_{i}D^{-\frac{1}{2}}v_{i}$ for some
$\beta_{i}\in\mathbb{C},1\leq i\leq k$. Therefore, for $m,n\geq 0$,

$\displaystyle\langle\Delta_{M_{z}}M_{z}^{n}1,z^{m}\rangle$ $\displaystyle=$
$\displaystyle\langle\Delta_{M_{z}}z^{n},M_{z}^{m}1\rangle.$ $\displaystyle=$
$\displaystyle\langle DW^{n}P1,W^{m}P1\rangle$

which implies

$\langle\Delta_{M_{z}}M_{z}^{n}1,M_{z}^{m}1\rangle=\sum_{i=1}^{k}\zeta_{i}^{n}\bar{\zeta_{i}}^{m}|\beta_{i}|^{2}.$
(2.5)

Equations (2.4) and (2.5) imply that

$\int_{\overline{\mathbb{D}}}z^{n}\bar{z}^{m}\
\mathrm{d}\mu=\sum_{i=1}^{k}\zeta_{i}^{n}\bar{\zeta_{i}}^{m}|\beta_{i}|^{2}.$

Since polynomials are dense in $D(\mu)$ ([primer, Corollary 7.3.4]), by Stone-
Weiertrass’s theorem, we have

$\mu=\sum_{i=1}^{k}|\beta_{i}|^{2}\delta_{\zeta_{i}}.$

Conversely, suppose that there exist $c_{i}>0$ and
$\zeta_{i}\in\overline{\mathbb{D}},1\leq i\leq k$ such that
$\mu=\sum_{i=1}^{k}c_{i}\delta_{\zeta_{i}}$.
As $\mu$ is supported at $\zeta_{i}\in\overline{\mathbb{D}},\ 1\leq i\leq k$,
for every $f\in D(\mu)$,

$\|f\|^{2}_{D(\mu)}=\|f\|^{2}_{\mathbb{H}^{2}}+\int_{\overline{\mathbb{D}}}\left\|\
\frac{f-f(\zeta)}{z-\zeta}\right\|_{\mathbb{H}^{2}}^{2}\mathrm{d}\mu(z)=\|f\|^{2}_{\mathbb{H}^{2}}+\sum_{i=1}^{k}c_{i}\left\|\
\frac{f-f(\zeta_{i})}{z-\zeta_{i}}\right\|_{\mathbb{H}^{2}}^{2}<\infty.$

This implies that

$f\in D(\mu)\iff\left\|\
\frac{f-f(\zeta_{i})}{z-\zeta_{i}}\right\|_{\mathbb{H}^{2}}^{2}<\infty,\ 1\leq
i\leq k.$

Since $D(\mu)$ is a reproducing kernel Hilbert space, point evaluations must
be bounded. Thus by the reproducing property, there exist linearly independent
vectors $\kappa_{\zeta_{1}},\kappa_{\zeta_{2}},\cdots\kappa_{\zeta_{k}}$in
$D(\mu)$ such that

$\langle f,\kappa_{\zeta_{i}}\rangle=f(\zeta_{i})\ ,\ \mbox{for every}\ f\in
D(\mu),1\leq i\leq k$

We now show that $\Delta_{M_{z}}$ is of rank $k$.
Let $F=\bigvee\\{\kappa_{\zeta_{i}}|1\leq i\leq k\\}$.
Claim: $\overline{ran(\Delta_{M_{z}})}=F$
Now, for $f\in D(\mu)$ consider,

$\displaystyle\langle\Delta_{M_{z}}f,f\rangle$ $\displaystyle=$
$\displaystyle\|M_{z}f\|^{2}_{D(\mu)}-\|f\|^{2}_{D(\mu)}$ $\displaystyle=$
$\displaystyle\int_{\overline{\mathbb{D}}}|f|^{2}\ \mathrm{d}\mu$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{k}c_{i}|f(\zeta_{i})|^{2}$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{k}c_{i}|\langle
f,\kappa_{\zeta_{i}}\rangle|^{2}.$

If $f\in F^{\perp}$ then $\langle\Delta_{M_{z}}f,f\rangle=0$. Since the
operator $\Delta_{M_{z}}\geq 0$, $\Delta_{M_{z}}f=0$ i.e
$f\in\ker(\Delta_{M_{z}})$. Therefore,
$F^{\perp}\subseteq\ker(\Delta_{M_{z}})$ or
$\overline{ran(\Delta_{M_{z}})}\subseteq F$.
Similarly, if $f\in\ker(\Delta_{M_{z}})$ then $\Delta_{M_{z}}f=0$ implies
$\sum_{i=1}^{k}c_{i}|\langle f,\kappa_{\zeta_{i}}\rangle|^{2}=0$. That is,
$\langle f,\kappa_{\zeta_{i}}\rangle=0$ for every $1\leq i\leq k\implies f\in
F^{\perp}$. In other words, $\ker(\Delta_{M_{z}})\subseteq F^{\perp}$ or
$F\subseteq\overline{ran(\Delta_{M_{z}})}$. Therefore,
$\overline{ran(\Delta_{M_{z}})}=F$ and hence $\Delta_{M_{z}}$ must be of rank
$k$. ∎

##  3\.  Completely Hyperexpansive Shift on a $\mathcal{H}(b)$ space

Let $b$ be a non-extreme point of the unit ball of $\mathbb{H}^{\infty}$ and
let $S$ denote the operator of multiplication by $z$ on $\mathbb{H}^{2}$. The
space $\mathcal{H}(b)$ is invariant under $S$. Let $T$ denote the restriction
of $S$ to ${\mathcal{H}(b)}.$ In other words, $T$ is the multiplication by $z$
on the space $\mathcal{H}(b)$.
We now obtain the necessary and sufficient conditions on $b$, so that $T$ is a
completely hyperexpansive operator.

######  Proposition 3.1.

Let $b$ be a nonextreme function in the unit ball of $H^{\infty}$. Then the
restriction of the shift operator to the de Branges Rovnyak space $H(b)$ is
completely hyperexpansive if and only if $b(z)=\frac{c+\gamma z}{1-\beta z}$
with $c,\gamma,\beta\in\mathbb{C}$ and $|\beta|<1$.

###### Proof.

Suppose that $T$ is a complete hyperexpansion. Then, by definition, $T$ is a
$k$-hyperexpansion for every $k\geq 1$. In particular, it is a
$2$-hyperexpansion. Using the necessary and sufficient conditions for $T$ to
be a $2$-hyperexpansion (Theorem 1.1) , $T$ is $2$-hyperexpansive (or concave)
if and only if $b$ has the form $b(z)=\frac{c+\gamma z}{1-\beta z}$ where
$c,\gamma,\beta\in\mathbb{C}$ and $|\beta|<1$.
Conversely, if

$b(z)=\frac{c+\gamma z}{1-\beta z}$

where $c,\gamma,\beta\in\mathbb{C}$ and $|\beta|<1$, then its pythagorean
mate, $a(z)=\frac{\rho-\sigma z}{1-\beta z}$ where $\rho$ and $\sigma$ follow
the conditions

  1. (i)

$\rho\geq|\sigma|$

  2. (ii)

$\rho^{2}+|\sigma|^{2}=1+|\beta|^{2}-|c|^{2}-|\gamma|^{2}\ \mbox{and }\
\rho^{2}|\sigma|^{2}=|\beta+\overline{c}\gamma|^{2}$

  3. (iii)

$\mbox{arg}(\rho)=\mbox{arg}(\beta+\overline{c}\gamma)$.

We observe that for any $h\in\mathcal{H}(b)$,

$\displaystyle\langle(I-3T^{*}T+3{T^{*}}^{2}T^{2}-{T^{*}}^{3}T^{3})h,h\rangle_{b}$
$\displaystyle=$
$\displaystyle(1+\|b\|_{b}^{2})\left(1-\left|\beta-\frac{a^{\prime}(0)}{a(0)}\right|^{2}\right)$
$\displaystyle\left\\{|\langle
h,{S^{*}b}^{2}\rangle_{b}|^{2}-2\mathrm{Re}\left(\bar{\frac{{a^{\prime}(0)}}{a(0)}}\langle
h,S^{*}b\rangle_{b}\overline{\langle h,{S^{*}}^{2}b\rangle_{b}}\right)\right.$
$\displaystyle\left.+\left[\left|\frac{a^{\prime}(0)}{a(0)}\right|^{2}-1\right]|\langle
h,S^{*}b\rangle_{b}|^{2}\right\\}$ $\displaystyle=$
$\displaystyle\left(1-\left|\beta-\frac{a^{\prime}(0)}{a(0)}\right|^{2}\right)\langle(I-2T^{*}T+{T^{*}}^{2}T^{2})h,h\rangle_{b}$

Inductively, it follows that, for every $n\geq 2$ and every
$h\in\mathcal{H}(b)$,

$\left\langle\left(\sum_{i=0}^{n}(-1)^{i}{n\choose
i}{T^{*}}^{i}T^{i}\right)h,h\right\rangle_{b}=\left(1-\left|\beta-\frac{a^{\prime}(0)}{a(0)}\right|^{2}\right)^{n-2}\langle(I-2T^{*}T+{T^{*}}^{2}T^{2})h,h\rangle_{b}$

Therefore, $T$ is completely hyperexpansive if

$\langle(I-2T^{*}T+{T^{*}}^{2}T^{2})h,h\rangle_{b}\leq 0\ \mbox{or}\
1-\left|\beta-\frac{a^{\prime}(0)}{a(0)}\right|^{2}=0.$

In either case, by appealing to Theorem 1.1, we conclude that $T$ must be
$2$-hyperexpansive and thus it follows that $T$ is completely hyperexpansive.

∎

We conclude this section by proving a result required in the proof of the main
theorem of this paper.

######  Proposition 3.2.

Let $b$ be a non-extreme point of unit ball of $\mathbb{H}^{\infty}$ and $a$
be its pythagorean mate. Then the defect operator $\Delta_{T}=T^{*}T-I$ is of
rank 1 and its range is spanned by $S^{*}b$ and its non-zero eigen value
equals $|a(0)|^{-2}\|S^{*}b\|_{b}^{2}$.

###### Proof.

Since $(b,a)$ forms a pair, we have

$\displaystyle T^{*}T$ $\displaystyle=$ $\displaystyle
I+|a(0)|^{-2}S^{*}b\otimes S^{*}b$
$\displaystyle\therefore\Delta_{T}=T^{*}T-I$ $\displaystyle=$
$\displaystyle|a(0)|^{-2}S^{*}b\otimes S^{*}b$

where $(x\otimes y)$ denotes the operator $(x\otimes y)(z)=\langle
z,y\rangle_{b}x\ ,\ z\in\mathcal{H}(b).$
Hence the matrix of the operator $\Delta_{T}$ relative to the decomposition
$\mathcal{H}(b)=\mathbb{C}S^{*}b\oplus(\mathbb{C}S^{*}b)^{\perp}$ is given by

$\begin{bmatrix}|a(0)|^{-2}\|S^{*}b\|_{b}^{2}&0\\\ 0&0\end{bmatrix}$

which implies the required result. ∎

##  4\. Proof of Theorem 1.3

This section is devoted to the proof of the main result of the paper viz
Theorem 1.3.
Let $b\in\mathbb{H}^{\infty}$ be nonextreme with $\|b\|_{\infty}\leq 1$ and
let $\mu$ be a finite positive Borel measure supported on
$\overline{\mathbb{D}}$. Let $T=M_{z}|_{\mathcal{H}(b)}$ denote the
restriction of multiplication by $z$ operator to $\mathcal{H}(b)$ space.
Proof of Theorem 1.3:

###### Proof.

Suppose that $\mathcal{H}(b)=D(\mu)$ with equality of norms.
As $M_{z}$ on $D(\mu)$ is cyclic, analytic and completely hyperexpansive, $T$
must be cyclic, analytic and completely hyperexpansive. Therefore, Theorem 3.1
implies that $b$ must be of the form

$b(z)=\frac{c+\gamma z}{1-\beta z}$

where $c,\gamma,\beta\in\mathbb{C}$ and $|\beta|<1$.
Similarly, Proposition 3.2 implies that the corresponding defect operator
$\Delta_{T}$ has rank 1. Therefore, the corresponding defect operator
$\Delta_{M_{z}}$ on $D(\mu)$ must also have rank 1. So, using Proposition 2.1
we conclude that

$\mu=|\alpha|^{2}\delta_{\lambda}\ ,\ \mbox{for some}\
\alpha\in\mathbb{C},\lambda\in\overline{\mathbb{D}}.$

Now, for Cauchy kernels $k_{w}(z)=\frac{1}{1-\overline{w}z}$ we have,

$D_{\lambda}(k_{w})=\frac{1}{2\pi}\int_{\mathbb{T}}\left|\frac{k_{w}(z)-k_{w}(\lambda)}{z-\lambda}\right||\mathrm{d}z|=\left\|\frac{\overline{w}}{1-\overline{w}\lambda}k_{w}\right\|_{\mathbb{H}^{2}}^{2}$

Since,

$D_{\mu}(f)=\int_{\overline{\mathbb{D}}}D_{\zeta}(f)\mathrm{d}\mu(\zeta)$

therefore we have,

$D_{\mu}(k_{w})=|\alpha|^{2}\left\|\frac{\overline{w}}{1-\overline{w}\lambda}k_{w}\right\|_{\mathbb{H}^{2}}^{2}=\frac{|\alpha|^{2}|w|^{2}}{|1-\overline{\lambda}w|^{2}(1-|w|^{2})}$

Similarly,

$\|k_{w}\|_{b}^{2}=\|k_{w}\|_{\mathbb{H}^{2}}^{2}+\|k_{w}^{+}\|_{\mathbb{H}^{2}}^{2}=\|k_{w}\|_{\mathbb{H}^{2}}^{2}+\frac{1}{(1-|w|^{2})}\frac{|b(w)|^{2}}{|a(w)|^{2}}$

But as $\mathcal{H}(b)=D(\mu)$ with equality of norms we have,

$\frac{1}{(1-|w|^{2})}\frac{|b(w)|^{2}}{|a(w)|^{2}}=\|k_{w}^{+}\|_{H^{2}}^{2}=D_{\mu}(k_{w})=\frac{|\alpha|^{2}|w|^{2}}{|1-\overline{\lambda}w|^{2}(1-|w|^{2})}$
(4.6)

In particular, letting $w=0$ we obtain $b(0)=0$. That is, $b(z)=\frac{\gamma
z}{1-\beta z}$.

Further , if $\phi=\frac{b}{a}$ , then equation (4.6) gives us

$|\phi(w)|^{2}=\frac{|\alpha|^{2}|w|^{2}}{|1-\overline{\lambda}w|^{2}}$

As the result is true for every $w$ , we obtain

$\phi(z)=\frac{\alpha z}{1-\overline{\lambda}z}.$

Let $p(z)=\alpha z,\ q(z)=1-\overline{\lambda}z$. We observe that $\phi$ is a
Smirnov class function. Then, the pair $(b,a)$ can be obtained by the
formulation,

$|a|^{2}=\frac{|q|^{2}}{|p|^{2}+|q|^{2}}\mbox{ and}\
|b|^{2}=\frac{|p|^{2}}{|p|^{2}+|q|^{2}}\ \mbox{a.e on}\ \mathbb{T}\ \mbox{(see
\cite[cite]{[\@@bibref{}{vol2}{}{}, Lemma 23.27]})}.$

Therefore,

$\displaystyle\frac{|\gamma|^{2}}{|1-\beta z|^{2}}$ $\displaystyle=$
$\displaystyle\frac{|\alpha|^{2}}{|1-\overline{\lambda}z|^{2}+|\alpha|^{2}|z|^{2}}\
,\ \mathrm{a.e\ on}\ \mathbb{T}$ $\displaystyle\therefore\
|\gamma|^{2}\left((1+|\lambda|^{2}+|\alpha|^{2})+2\mathrm{Re}(\overline{\lambda}z)\right)$
$\displaystyle=$
$\displaystyle|\alpha|^{2}\left((1+|\beta|^{2})+2\mathrm{Re}(\beta z)\right),\
\mathrm{a.e\ on}\ \mathbb{T}$

Substituting $z=1,\ -1,\ i,\ -i$ we obtain the following:
If $\lambda=0$ then

$|\gamma|^{2}=\frac{1}{1+|\alpha|^{2}}$

and therefore $\beta=0$
Similarly, if $\lambda\neq 0$, then

$|\gamma|^{2}=\frac{|\alpha|^{2}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})\pm\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right).$

But,

$|\gamma|^{2}=\frac{|\alpha|^{2}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})+\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right)$

implies $|\beta|>1$.
Therefore,

$|\gamma|^{2}=\frac{|\alpha|^{2}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})-\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right).$

Hence

$\beta=\frac{\overline{\lambda}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})-\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right).$

Letting $A_{\lambda}=\gamma$ and $B_{\lambda}=\beta$ we obtain the required
result.
Conversely, suppose that $\mu=|\alpha|^{2}\delta_{\lambda},\ \mbox{for some
}\alpha\in\mathbb{C},\ \lambda\in\overline{\mathbb{D}}$ and
$b(z)=\frac{A_{\lambda}z}{1-B_{\lambda}z}$ where $A_{\lambda}$ and
$B_{\lambda}$ are constants depending on $\alpha$ and $\lambda$ such that

$|A_{\lambda}|^{2}=\begin{cases}\left(\frac{|\alpha|^{2}}{2|\lambda|^{2}}\left((1+|\alpha|^{2}+|\lambda|^{2})-\sqrt{(1+|\alpha|^{2}+|\lambda|^{2})^{2}-4|\lambda|^{2}}\right)\right)&\mbox{,
if}\ \lambda\neq 0\\\ \left(\dfrac{1}{1+|\alpha|^{2}}\right)&\mbox{, if}\
\lambda=0\end{cases}$

and

$B_{\lambda}=\frac{|A_{\lambda}|^{2}\overline{\lambda}}{|\alpha|^{2}}.$

Then by Proposition 2.1, the defect operator of $M_{z}$ on $D(\mu)$,
$\Delta_{M_{z}}$ must be of rank 1.
Similarly, applying Proposition 3.1 and Proposition 3.2, we conclude that the
operator $T$ must be cyclic, analytic, completely hyperexpansive and its
defect operator, $\Delta_{T}$ must be of rank 1.

Recall that, by Aleman’s result [ale, Theorem 2.5], T must be unitarily
equivalent to the multiplication by $z$ operator on $D(\mu^{\prime})$ for some
non-negative, Borel measure $\mu^{\prime}$ supported on
$\overline{\mathbb{D}}$ and $\mathcal{H}(b)=D(\mu^{\prime})$ with equality of
norms. But, previous part of this theorem implies that $\mu$ and
$\mu^{\prime}$ must be equal.∎

The above result provides a formulation for computing the non-extreme symbol
$b$, given a positive, point mass measure $\mu$ on $\overline{\mathbb{D}}$, so
that $D(\mu)=\mathcal{H}(b)$ with equality of norms. We use the following
examples to illustrate the construction.

###### Example 4.1 :

Suppose $\mu=\delta_{0}$ i.e $|\alpha|^{2}=1$ and $\lambda=0$. Then by the
formula prescribed in the above theorem, one obtains
$A_{\lambda}=\frac{1}{\sqrt{2}}$ and $B_{\lambda}=0$.

###### Example 4.2 :

Suppose $\mu=\delta_{\frac{1}{2}}$ i.e $|\alpha|^{2}=1$ and
$\lambda=\frac{1}{2}$. Then,
$A_{\lambda}=\left(\dfrac{{9-\sqrt{65}}}{2}\right)^{\frac{1}{2}}$ and
$B_{\lambda}=\dfrac{(9-\sqrt{65})}{4}$.

We conclude this paper by deducing Corollary 1.2 from Theorem 1.3.
Proof of Corollary 1.2:

###### Proof.

Letting $|\lambda|=1$ in Theorem 1.3, we obtain

$|A_{\lambda}|^{2}=\frac{|\alpha|^{2}}{2}\left(2+|\alpha|^{2}-|\alpha|\sqrt{(4+|\alpha|^{2})}\right),$

Similarly,
$B_{\lambda}=\frac{|A_{\lambda}|^{2}\overline{\lambda}}{|\alpha|^{2}}$ implies

$|B_{\lambda}|=\frac{|A_{\lambda}|^{2}}{|\alpha|^{2}}\neq 0$

Now ,

$\displaystyle|B_{\lambda}|+|A_{\lambda}|=1$ $\displaystyle\iff$
$\displaystyle\frac{|A_{\lambda}|^{2}}{|\alpha|^{2}}+|A_{\lambda}|=1$
$\displaystyle\iff$
$\displaystyle|A_{\lambda}|^{2}=\frac{|\alpha|^{2}}{2}\left(2+|\alpha|^{2}\pm|\alpha|^{2}\sqrt{(4+|\alpha|^{2})}\right)$

Therefore,

$|\alpha|^{2}=\frac{|A_{\lambda}|^{2}}{|B_{\lambda}|},\
\lambda=\frac{\overline{B_{\lambda}}}{|B_{\lambda}|}.$

This completes the proof. ∎

### Acknowledgments:

  1. (1)

The authors would like to thank Prof. Sameer Chavan (IIT, Kanpur) for several
useful discussions throughout the preparation of this paper.

  2. (2)

The present work is carried out at the research centre at the Department of
Mathematics, S. P. College, Pune, India (autonomous).

## References


# Optimal Ensemble Control of Matter-Wave Splitting
in Bose-Einstein Condensates

Andre Luiz P. de Lima1, Andrew K. Harter2, Michael J. Martin2, and Anatoly
Zlotnik3 This project was supported by the LDRD program and the Center for
Nonlinear Studies at Los Alamos National Laboratory. Research conducted at Los
Alamos National Laboratory is done under the auspices of the National Nuclear
Security Administration of the U.S. Department of Energy under Contract No.
89233218CNA000001.1Andre Luiz P. de Lima is a Ph.D. student in the Electrical
& Systems Engineering Department at Washington University in St. Louis, St.
Louis, MO, USA<EMAIL_ADDRESS>Andrew K. Harter and Michael J. Martin are
with the Quantum Technologies Team at Los Alamos National Laboratory, Los
Alamos, NM, USA {harter<EMAIL_ADDRESS>Anatoly Zlotnik is in the Applied
Mathematics & Plasma Physics Group, Los Alamos National Laboratory, Los
Alamos, NM, USA<EMAIL_ADDRESS>
###### Abstract

We present a framework for designing optimal optical pulses for the matter-
wave splitting of a Bose-Einstein Condensate (BEC) under the influence of
experimental inhomogeneities, so that the sample is transferred from an
initial rest position into a singular higher diffraction order. To represent
the evolution of the population of atoms, the Schrödinger’s equation is
reinterpreted as a parameterized ensemble of dynamical units that are
disparately impacted by the beam light-shift potential in a continuous manner.
The derived infinite-dimensional coupled Raman-Nath equations are truncated to
a finite system of diffraction levels, and we suppose that the parameter that
defines the inhomogeneity in the control applied to the ensemble system is
restricted to a compact interval. We first design baseline square pulse
sequences for the excitation of BEC beam-splitter states following a previous
study, subject to dynamic constraints for either a nominal system assuming no
inhomogeneity or for several samples of the uncertain parameter. We then
approximate the continuum state-space of the ensemble of dynamics using a
spectral approach based on Legendre moments, which is truncated at a finite
order. Control functions that steer the BEC system from an equivalent rest
position to a desired final excitation are designed using a constrained
optimal control approach developed for handling nonlinear dynamics. This
representation results in a minimal dimension of the computational problem and
is shown to be highly robust to inhomogeneity in comparison to the baseline
approach. Our method accomplishes the BEC-splitting state transfer for each
subsystem in the ensemble, and is promising for precise excitation in
experimental settings where robustness to environmental and intrinsic noise is
paramount.

##  I INTRODUCTION

Numerous challenges in quantum science and technology involve design,
observation, or control of bilinear dynamical systems [1]. A compelling
application in metrology involves cold atom interferometry, in which optical
standing waves are used to split and recombine matter waves [2]. In this
setting, a time-varying laser pulse is applied to steer a quantum system
between states of interest. This experimental approach has demonstrated
significant sensitivity to accelerations and rotations, and has the potential
to provide real improvement in fundamental physics experiments [3]. However,
experimental uncertainties are inherent in these systems, and this reduces the
precision in the splitting of matter waves, which can lead to inconsistent
diffraction fidelity [4].

Robustness to experimental inhomogeneities, such as laser intensity, has been
investigated in atom interferometry applications. Previous approaches have
attempted to minimize or nullify the influence of these dynamical
disturbances, for example adiabatic transfer techniques, which do not depend
on exact laser intensity [5]. Other methods seek to optimally design light
pulses to improve tolerance to systematic inhomogeneities by using specific
control profile shapes, such as composite pulses [6], or large-scale
optimization, such as GRadient Ascent Pulse Engineering (GRAPE) [7]. A
practical means to represent inhomogeneities is through the light amplitude
parameterized description of the system time-evolution dynamics. This
mathematical description enables the state evolution of the sample to be
represented by a collection of ordinary differential equations, each of which
describes the possible evolution taken by an individual atom in the experiment
subject to the inexact effect of the applied laser pulse.

Ensemble control theory was developed to address the need to understand the
evolution of systems comprised by a very large (potentially infinite) number
of similar dynamical units, such as molecular, atomic and quantum systems,
when manipulated by a single universally applied control field [8]. This area
of mathematical systems theory focuses on defining controllability conditions
[9], measurement processing [10], and control design approaches [11] for
parameterized dynamical systems. Ensemble control theory has been used to
address challenging control tasks for large-scale dynamical problems,
including brain dynamics [12] and Magnetic Resonance Imaging (MRI) [13].
Furthermore, ensemble control was shown to be effective in pulse design for
Nuclear Magnetic Resonance (NMR) accounting for the influence of systematic
disturbances [14, 13]. The mathematical setting used for control synthesis for
bilinear ensemble systems [13] can be applied to develop laser pulses for the
atom interferometry procedure, and indeed the simplification of the Raman-Nath
equations truncated to a two-level system results in the Bloch equations that
appear in NMR [15, 16].

Ensemble systems are challenging because of infinite-dimensionality that
results from a parameter that varies continuously on a compact interval.
Whereas the state of each individual ensemble element is finite-dimensional,
the state of an ensemble system lies on a Hilbert space. In this context,
methods designed for unified representations have been studied with the goal
of reducing computational complexity. Moment representations have been shown
to be effective in significantly reducing the size of ensemble system dynamics
[17]. Moment dynamics are essentially a spectral approximation of the ensemble
in Hilbert space, with the spectral accuracy property that ensures exponential
convergence to the continuum state with increasing approximation order. This
representation enables design of control sequences for ensemble systems in a
space of greatly reduced dimension without significant loss of precision [18].

In this paper, we revisit the problem of matter-wave splitting involving a
Bose-Einstein condensate (BEC) in a standing light wave potential [15]. We
consider optical standing-wave beam splitters for which the effect of an
optical pulse is subject to experimental inhomogeneities, which we represent
by an uncertain factor applied to the light wave amplitude envelope. This
results in an ensemble system, whose evolution we describe using the Legendre
moment dynamics. In this reduced space, we design pulse sequences that
transfer the state of the ensemble to a desired set-point, and evaluate the
improved fidelity of the BEC momentum transfer within the quantum state space.

This project is organized as follows. In Section II, we derive the
mathematical representation of the BEC-splitting system from the Schrödinger’s
equation, and characterize the source of inhomogeneity. In Section III, we
then reproduce a state-of-the-art method of using square pulse sequences for
optical splitting of matter-waves without compensating for inherent
inhomogeneity [2]. In Section IV, we describe the Legendre moment
representation of the BEC splitting dynamics, our proposed method for light
pulse design within the moment space, and the resulting optimal control
formulation. We present the computed optimal controls and simulations of the
resulting momentum transfers in Section V for several diffraction orders,
which show a significant improvement in fidelity with respect to the benchmark
method. Finally, we conclude in Section VI.

##  II ENSEMBLE SYSTEM FOR BEC SPLITTING

We consider a dilute BEC composed of an initially stationary population of
atoms, which is exposed to a standing-wave light field that excites the
cluster of particles with similar, yet distinct intensities. The governing
dynamics for this system are described by a one-dimensional Schrödinger’s
equation for a single atom. Inhomogeneity in the light beam is represented by
a multiplicative parameter $\varepsilon$ that scales the light shift potential
amplitude $\Omega(t)$, which leads to a wave function

$\\!\\!i\dot{\psi}(x,t,\varepsilon)\\!=\\!\left(-\frac{\hbar}{2m}\frac{d^{2}}{dx^{2}}\\!+\\!\varepsilon\Omega(t)cos\left(2k_{0}x\right)\right)\psi(x,t,\varepsilon).$
(1)

The time-dependent parameter $\Omega(t)$ serves as the control input to the
system available to the experimenter, and $k_{0}$ is the wave vector of the
light field. We suppose that $\varepsilon\in K$ is the physical representation
of the unequal influence of the control pulse on individual atoms in the BEC
sample, where $K\equiv[1-\delta,1+\delta]$ is a compact interval. Each atomic
particle can be viewed as corresponding to a unique $\varepsilon$, which
permits us to use this parameter as an identifier to index the individual
systems in the sample.

Equation (1) can be further developed into a sequence of coupled Raman-Nath
equations by describing the wave function within the Bloch basis. The
resulting differential equations dictate the time evolution of the beam-
splitter state populations (represented by ground state $C_{0}$, and symmetric
and anti-symmetric superpositions of momentum states $C_{2n}^{+}$ and
$C_{2n}^{-}$, for $n=0,1,2,\ldots$) and the manner in which these interact.
This infinite-dimensional system of ordinary differential equations (ODEs) is
[15]

$i\dot{C_{2n}}(k,t)=\frac{\hbar}{2m}(2nk_{0}+k)^{2}C_{2n}(k,t)\\\
+\frac{\varepsilon\Omega(t)}{2}[C_{2n-2}(k,t)+C_{2n+2}(k,t)].$ (2)

We make several assumptions in order to further simplify these infinite-
dimensional Raman-Nath equations:

  * •

The dependence on parameter $k$ is dropped by assuming the system to be
narrowly distributed around the value of $k=0$.

  * •

The coupling of the states $C_{0}$ and $C_{2n}^{+}$ with their anti-symmetric
counterparts $C_{2n}^{-}$ is considered negligible by assuming that $k\ll
k_{0}$.

  * •

A value $N^{+}$ is defined such that all levels above $C_{2N^{+}}^{+}$ are
permanently unpopulated. This is achievable by assuming that all atoms are at
rest at the beginning and that $\Omega(t)/2\ll(2N^{+})^{2}\hbar
k_{0}^{2}/2m=(2N^{+})^{2}\omega_{r}$, where $\omega_{r}$ is the photon recoil
frequency.

This mathematical setting has been developed in a variety of previous studies
on matter-wave splitting techniques [19, 15, 2]. The resulting state-space
system is described by

$\frac{d}{dt}\begin{bmatrix}C_{0}(t,\varepsilon)\\\
C_{2}^{+}(t,\varepsilon)\\\ \vdots\\\
C_{2N}^{+}(t,\varepsilon)\end{bmatrix}=-i\omega_{r}A(\varepsilon,\Omega(t))\begin{bmatrix}C_{0}(t,\varepsilon)\\\
C_{2}^{+}(t,\varepsilon)\\\ \vdots\\\ C_{2N}^{+}(t,\varepsilon)\end{bmatrix},$
(3)

where $A(t,\varepsilon)$ is a real $(N^{+}+1\times N^{+}+1)$ symmetric matrix
defined by

$A(\varepsilon,\Omega(t))=\begin{bmatrix}0&0&0&\ldots&0\\\ 0&4&0&\ldots&0\\\
0&0&16&\ldots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\ldots&(2N^{+})^{2}\end{bmatrix}\\\
+\frac{\varepsilon\Omega(t)}{2\omega_{r}}\begin{bmatrix}0&\sqrt{2}&0&0&\ldots&0\\\
\sqrt{2}&0&1&0&\ldots&0\\\ 0&1&0&1&\ldots&0\\\ 0&0&1&0&\ddots&\vdots\\\
\vdots&\vdots&\vdots&\ddots&\ddots&1\\\ 0&0&0&\ldots&1&0\end{bmatrix}.$ (4)

For the system in Equations (3)-(4), we define the state vector
$C^{+}=[C_{0},C_{2}^{+},...,C_{2N^{+}}^{+}]\in\mathbb{C}^{N^{+}+1}$. A goal of
this study is to enable constrained optimization approaches for pulse design,
for which the complex dynamics in Equations (3)-(4) must be converted to a
real-valued system. We thus define
$C^{+}_{\mathbb{R}}\in\mathbb{R}^{2(N^{+}+1)}$ such that
$C^{+}_{\mathbb{R}}=[C_{0,Re},C_{0,Im},C^{+}_{2,Re},C^{+}_{2,Im},...,C^{+}_{2N^{+},Re},C^{+}_{2N^{+},Im}]$,
in terms of which the time evolution dynamics are

$\frac{d}{dt}C^{+}_{\mathbb{R}}=\omega_{r}A(\varepsilon,\Omega(t))\otimes\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}C^{+}_{\mathbb{R}},$ (5)

where $\otimes$ is the Kronecker product.

In the above representation, the momentum dynamics for each atom are truncated
to a finite-dimensional quantum system of $N^{+}+1$ levels. Nonetheless, due
to the presence of parameter $\varepsilon$, there is a collection of
infinitely many systems that must be steered simultaneously by a common
control input. Next, we examine methods to design control pulses for momentum
transfer of the system (5) that are invariant to a continuum of values of the
uncertain parameter $\varepsilon$.

##  III SPLITTING MATTER USING SQUARE PULSES

A previously proposed approach for splitting matter waves uses a sequence of
two square pulses for the envelope of an optical input that is applied to a
dilute BEC in a standing light wave potential. The approach has been shown to
have a high degree of precision in simulation, and can be applied in practice
[15, 20]. The drawback of this method is its use of a single nominal system to
define the dynamic response of each atom in the BEC to the optical pulse, so
there is no compensation for inherent inhomogeneities. We describe the square-
pulse design procedure below, and extend the original methodology with pulses
of equal amplitude to the design of sequences of pulses with different
amplitudes, which has been shown to enable higher fidelity [2].

###  III-A Square optical pulse design for a single system

The square pulse optimal design problem involves the optimization of five
parameters, which are illustrated in Figure 1. The decision variables are the
time durations $\tau_{1}$ and $\tau_{3}$ of the pulses, the time interval
$\tau_{2}$ between the pulses, and the amplitudes $\Omega_{1}$ and
$\Omega_{2}$ of the first and second square pulses, respectively. These
parameters determine the final state of the system after the pulse is
employed, which can be described by continuously applying the variation of
parameters formula to Equation (3). The state at the terminal time
$\tau_{ps}=\tau_{1}+\tau_{2}+\tau_{3}$ of the pulse sequence is

$\begin{gathered}C^{+}(\tau_{ps},\varepsilon)=e^{-i\omega_{r}\tau_{3}A_{3}}e^{-i\omega_{r}\tau_{2}A_{2}}e^{-i\omega_{r}\tau_{1}A_{1}}C^{+}(0,\varepsilon)\\\
\text{s.t. }A_{1}=A(\varepsilon,\Omega_{1}),\,A_{2}=A(\varepsilon,0),\\\
A_{3}=A(\varepsilon,\Omega_{2}),\,\tau_{ps}=\tau_{1}+\tau_{2}+\tau_{3}.\end{gathered}$
(6)

In the baseline scenario, we suppose that $\varepsilon=1$. The control design
goal is to transfer the momentum state of the sample from an initial rest
state $C^{+}(0)=[1,0,0,\ldots,0]^{T}$ to a final desired state $C_{f}^{+}$.
This desired state is typically defined with respect to the $L^{2}$-norm (in
quantum momentum space), in which an energy level $2n$ is reached at the end
of the sequence of pulses, so that $\lVert C_{f,2n}^{+}\rVert_{2}=\lVert
C_{f}^{+}\rVert_{2}=1$. For this control task, the optimization problem can be
formulated as

$\displaystyle\min_{\Omega_{1},\Omega_{2},\tau_{1},\tau_{2},\tau_{3}}$
$\displaystyle\lVert|C_{f}^{+}|-|C^{+}(\tau_{ps},1)|\rVert_{2}\text{ as in
Eq.}\eqref{eq:time_evolution}$ (7) $\displaystyle s.t.$
$\displaystyle\textrm{Dynamics in Equations
\eqref{eq:BEC_ensemble_A}-\eqref{eq:BEC_ensemble_Real}}.$

The procedure described above follows the original methodology for the design
of square pulse sequences for matter wave splitting. Because a nominal
parameter $\varepsilon$ is used, the approach does not account for non-uniform
influence of the pulse on the dynamical behavior of individual atoms in the
BEC, and we expect that control performance is not robust to values of
$\varepsilon$ away from unity. We extend the formulation below by adding
dynamic constraints for samples of $\varepsilon\subset K$ as an attempt to
account for the dynamics of an ensemble of atoms that are affected
inhomogeneously by the applied pulse sequence.

Figure 1: Illustration of a square pulse sequence control $\Omega(t)$ and
parameters.

###  III-B Square optical pulse design for ensemble system

We extend the optimization problem (7) defined in the previous section to
account for multiple values of $\varepsilon$. Because it is intractable to
define dynamic constraints for each system in an infinite-dimensional
ensemble, we define an objective function as a summation of minimum $L_{2}$
norm objectives for $m$ values of $\varepsilon$ sampled directly from the
parameter domain $K$. The optimization is conducted subject to dynamic
constraints for each of the same $m$ values of $\varepsilon$.

We define the finite collection of parameters
$\varepsilon_{1:m}=\\{\varepsilon_{1},\varepsilon_{2},\ldots,\varepsilon_{m}\\}\subset
K$, such that $\varepsilon_{1}=\mathop{\arg\min}_{K}\varepsilon$,
$\varepsilon_{m}=\mathop{\arg\max}_{K}\varepsilon$, with intermediate samples
uniformly spaced on $[\varepsilon_{1},\varepsilon_{m}]$. The simultaneous
optimization for all systems indexed by $i\in\\{1,2,\ldots,m\\}$ is formulated
as

$\displaystyle\min_{\Omega_{1},\Omega_{2},\tau_{1},\tau_{2},\tau_{3}}$
$\displaystyle\sum_{i=1}^{m}\lVert|C_{f}^{+}|-|C^{+}(\tau_{ps},\varepsilon_{i})|\rVert_{2}$
(8) $\displaystyle s.t.$ $\displaystyle\textrm{Dynamics in Equations
\eqref{eq:BEC_ensemble_A}-\eqref{eq:BEC_ensemble_Real}
$\varepsilon=\varepsilon_{i}$},$ $\displaystyle\qquad\forall i=1,\ldots,m.$

The control pulse sequences that are obtained by solving the optimization
problems (7) and (8) will be used as benchmarks in Section V to demonstrate
the performance improvement gained by way of the proposed moment dynamics
method.

##  IV CONSTRAINED ENSEMBLE PROBLEM USING MOMENT DYNAMICS

We present a novel approach to improve light pulse design methods for matter-
wave splitting by compensating for the inherent inhomogeneities that are
present in physical systems. In this section, we derive a polynomial moment
based representation of the ensemble dynamics of the matter wave splitting
process, and formulate an ensemble control problem for optimal state transfer
of a BEC sample in momentum space.

The polynomial moment ensemble approach has been applied recently to design
open-loop controls to transfer ensemble systems between initial and target
states in Hilbert space, demonstrating low terminal state error [17, 21, 9].
The method enables a control problem for an infinite-dimensional continuum to
be reduced to a finite-dimensional control problem by transforming the state
space using a polynomial approximation that can then be truncated at a
sufficiently high order. The transformation facilitates the design of a single
control input that steers the moments to a desired representation in the
spectral moment space, which is equivalent to steering each individual element
of the ensemble to its target in the state space [14].

The moment-based approach for ensemble control systems can be defined using
any basis from a diverse collection of orthogonal polynomials. For this
project, we have chosen to use Legendre polynomials as the basis of the moment
system, because this particular basis has desirable properties for
approximation of continuous functions on compact intervals [18]. This
representation has also been used previously for controlling the Bloch
equations that represent a two-level quantum system in the context of NMR
[14].

###  IV-A Moment Systems using Legendre Polynomials

The concept of moment space representation begins with a function defined
within a separable Hilbert space $\mathcal{H}$. If a basis
{$|\psi_{k}\rangle:k\in\mathbb{N}^{d}$} can be defined (for simplicity, we
will use $d=1$), then the moment of order $k$ for a function $x(t)$ can be
defined by

$m_{k}(t)=\langle x(t)|\psi_{k}\rangle.$ (9)

For an ensemble of systems indexed by a parameter varying in the compact
interval $K=[-1,1]$, it is possible to obtain the equivalent Legendre moments
by using the respective orthogonal polynomials $P_{k}(\varepsilon)$, obtained
using the recursive relation given, after normalization, by

$\varepsilon P_{k}=c_{k-1}P_{k-1}+c_{k}P_{k+1},$ (10)

where $P_{0}(\varepsilon)=1/\sqrt{2}$,
$P_{1}(\varepsilon)=\sqrt{3/2}\varepsilon$, $\ldots$, are the normalized
Legendre polynomials and $c_{k}=(k+1)/\sqrt{(2k+1)(2k+3)}$. In this setting,
the Legendre moment of order $k$ for the vector of quantum states
$C^{+}(t,\varepsilon)$ is defined by

$m_{k}(t)=\langle
C^{+}(t),P_{k}\rangle=\int_{-1}^{1}C^{+}(t,\varepsilon)P_{k}(\varepsilon)d\varepsilon.$
(11)

Taking the time derivative of Equation (11) yields the equivalent moment
dynamics of the ensemble system. The use of the Legendre polynomial basis to
define moment spaces was shown to possess a variety of advantages for
representing ensemble systems. For instance, the Legendre polynomials form an
orthogonal basis on the Hilbert space $\mathcal{H}$, which induces an isometry
between the moment space and the quantum ensemble system. The preservation of
metrics through this transformation has advantages for defining tractable
optimization formulations to represent optimal control problems for ensemble
systems. The recursive relation in Equation (10) also enables the description
of moment dynamics defined by bounded operators. Finally, the control profile
is preserved with transformation between the moment space and state space, so
that practical control design is possible using the moment dynamics method.

###  IV-B BEC Splitting in the Moment Space

We apply the transformation in Equation (11) to the time-evolution model for
the quantum system described by Equations (4)-(5) to obtain the dynamics of
the corresponding Legendre moments. The quantum ensemble system is defined for
the ensemble parameter $\varepsilon\in K=[1-\delta,1+\delta]$, such that
$\delta\in(0,1)$, which is a different compact interval than the domain of the
Legendre polynomial basis. Therefore, we define a parameter $\varepsilon^{*}$
such that $\varepsilon=1+\delta\varepsilon^{*}$, meaning that
$\varepsilon^{*}\in[-1,1]$, and the moments are calculated using the Hilbert
space related to parameter $\varepsilon^{*}$. In this setting, the moment
dynamics for $m(t)=[m_{0}(t),m_{1}(t),\ldots,m_{N-1}(t)]$ are described by

$\displaystyle\dot{m(t)}=\omega_{r}\left[I(N)\otimes
A(1,\Omega(t))\otimes\begin{bmatrix}0&-1\\\ 1&0\end{bmatrix}\right.$ (12)
$\displaystyle\quad\quad\left.+\mathcal{C}(N)\otimes\left(A(\delta,\Omega(t))-A(\delta,0)\right)\otimes\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}\right]\\!\cdot\\!m(t),$

where $I(N)$ is an identity matrix of dimension $N\times N$ and
$\mathcal{C}(N)$ is defined by

$\mathcal{C}(N)=\begin{bmatrix}0&c_{0}&&&\\\ c_{0}&0&c_{1}&&\\\
&c_{1}&\ddots&\ddots&&\\\ &&\ddots&0&c_{N-2}\\\ &&&c_{N-2}&0\end{bmatrix}.$
(13)

With the Legendre moment dynamics as defined in Equation (12), we can state an
optimization problem that approximates the optimal control problem for the
ensemble in state-space. The orthogonality property of the Legendre polynomial
basis defines an isometry between the two studied spaces, which enables a
metric to be used to specify an objective for the optimization problem.
Moreover, it is also possible to include inequality constraints in the
optimization problem, which we will apply in the control design. We are
interested specifically in enforcing limits on the absolute amplitude and the
rate of change of the control function.

Specifying the final state in moment space requires an additional
reformulation. In the state-space, the optimization objective is defined in
terms of the absolute value of the final state, as in Equation (7). To express
an objective function that depends on the state of an ensemble system, we
formulate the objective in moment space. To that end, we minimize the
magnitude of the undesired energy levels by aiming to nullify their states,
i.e., $\lVert
C_{f}^{+}\circ(\mathbf{1}_{2(N^{+}+1)}-e_{n+1}\otimes\mathbf{1}_{2})\rVert_{2}=0$.
The resulting optimization problem in moment space is

$\displaystyle\min_{\Omega(t)}$ $\displaystyle\lVert
m(T)\circ[\mathbf{1}_{N}\otimes(\mathbf{1}_{2(N^{+}+1)}-e_{n+1}\otimes\mathbf{1}_{2})]\rVert_{2}$
(14) $\displaystyle s.t.$ $\displaystyle\textrm{Dynamics in Equations
\eqref{eq:BEC_moment}-\eqref{eq:BEC_moment_C}},$
$\displaystyle\Omega_{min}\leq\Omega(t)\leq\Omega_{max},$
$\displaystyle\Delta\Omega_{min}\leq\dot{\Omega}(t)\leq\Delta\Omega_{max},$

where $T$ is a predefined time horizon, similar to the parameter $\tau_{ps}$
used in problems (7) and (8).

##  V CONTROL SYNTHESIS RESULTS

The Legendre moments method for matter-wave splitting ensemble control design
is evaluated by comparing its performance with that of the square pulse
sequence method. We first compare the Legendre moments approach with the
controls obtained by solving the problem in Equation (7). A similar comparison
is then done with respect to controls obtained by solving the problem in
Equation (8), which seeks to account for variation in $\varepsilon$ by
sampling values directly in the ensemble space $K$. To compare the performance
of the various control designs in achieving the goal of transferring the
ensemble to the target in Hilbert space on $K$, we define an index $I_{e}$
based on the objective function for the pulse design in the quantum state
space. This performance index, which we seek to minimize, is defined as

$I_{e}=\int_{1-\delta}^{1+\delta}\lVert|C_{f}^{+}|-|C^{+}(T,\varepsilon)|\rVert_{2}d\varepsilon.$
(15)

###  V-A Computational Approach

We solve the optimal control problem (14) using a general iterative scheme for
computing optimal control inputs for nonlinear systems [22]. For a general
continuous-time nonlinear control system of the form
$\dot{x}(t)=f(x(t),u(t))$, where $x(t)$ and $u(t)$ are state and control
vectors, respectively, a zero-order hold assumption is used to define the
evolution of the system using a piece-wise constant control input. Beyond
using a traditional first-order discretization and linearization, we use an
approach that takes advantage of a higher-order Taylor series expansion to
accurately represent linearization of nonlinear dynamics and their Jacobian at
each iteration [23]. Beginning with an initial guess of the control function,
the algorithm solves an approximation of the optimization problem in (14) to
yield the optimal vector $\Delta U$ of time-discretized variations in the
control input that improves the objective value. This results in a sequence of
quadratic programs, which, though each can be efficiently solved, require
highly complex algebraic formulation of the dynamics in Equations (12)-(13)
when truncated at a high momentum level $N^{+}$ and sufficiently high
polynomial moment order to represent the quantum and ensemble dynamics in
their respective Hilbert spaces. In contrast to previous applications of this
approach to the control of quantum ensemble systems [18], the size and
complexity of the dynamics examined here requires the iterative use of
symbolic algebra to compute the Taylor expansion of the nonlinear dynamics,
which becomes extremely complex.

Once the moment dynamics and the respective symbolic expressions are obtained,
the pulse is designed by the iterative application of quadratic programming to
solve the optimization problem in Equation (14). To apply the methodology
defined in [23], the time domain $[0,T]$ is finely discretized by using a
small time step. We use a nominal value of $\Delta t=0.001$ over a total time
horizon of $T=3$, and this results in a total of 3000 variables and 12000
constraints for optimizing the variation $\Delta U\in\mathbb{R}^{3000}$ that
is added to the control function at each iteration. The same momentum level
truncation at $N^{+}=9$ and moment order $N=20$ are used in all of the
computations.

Algorithm 1 Iterative Moment Pulse Design Algorithm

1:$n$, Initial pulse guess $U$, Moment dynamics $M(t,m(t),u(t))$, number of
time steps $N_{steps}$, $\lambda$,
$f(m(t))=m(t)\circ[\mathbf{1}_{N}\otimes(\mathbf{1}_{2(N^{+}+1)}-e_{n+1}\otimes\mathbf{1}_{2})]$

2:$\text{Cost Function in Equation
\ref{eq:moment_optimization}}\leq\textbf{Tolerance}$

3:while $\lVert f(m(T))\rVert\geq\textbf{Tolerance}$ do

4: $m(0),m(\Delta t),\ldots,m(T)\leftarrow F(t,m(t),U)$

5: $i\leftarrow 1$

6: while $i\leq N_{steps}$ do

7: $A_{i}\leftarrow\frac{\partial}{\partial m}F((i)\Delta t,m((i)\Delta
t),U_{i+1})$

8: $B_{i}\leftarrow\frac{\partial}{\partial U}F((i)\Delta t,m((i)\Delta
t),U_{i+1})$

9: $i\leftarrow i+1$

10: end while

11: $H\leftarrow[A_{N_{steps}}\ldots A_{2}B_{1}|A_{N_{steps-1}}\ldots
A_{3}B_{2}|\ldots$ $|A_{N_{steps}}B_{N_{steps}-1}|B_{N_{steps}}]$

12: $\Delta U\leftarrow\min_{\Delta U}\Delta U^{T}(H^{T}H+\lambda I)\Delta U$
$+f(m(T))^{T}H\Delta U$

13: $U\leftarrow U+\Delta U$

14:end while

The computational procedure used to synthesize ensemble controls is described
in Algorithm 1, in which a cost function is iteratively reduced until a
desired error tolerance is reached. The methodology can be described as
consisting of three stages: (1) propagation of moment states using the current
control profile; (2) linearization of time-localized dynamics; (3) and an
update of the control profile. Both the propagation of moment states and
respective linearizations are performed using a Taylor series expansion of
higher order over the symbolic expression that represents the dynamics in
Equation (12). This approach is adopted for improved precision in the
linearized approximation defined by the set of matrices $A_{k}$ and $B_{k}$,
which are coupled in a matrix $H$, used to estimate the variation of the
moment state at the final time step. The optimization problem can be
represented using a quadratic program with a quadratic objective function and
linear constraints, which we solve using a quadratic programming solver. This
function has a penalty term scaled using a parameter $\lambda$, which is used
in order to avoid high variations in the defined control variation $\Delta U$,
which would invalidate previously obtained linearized transformations. The
related code is developed and executed on MATLAB™, version R2023a.

In the following computational studies, the constraint bound values defined in
Equation (14) are set to $\Omega_{max}=100$, $\Delta\Omega_{min}=-500s^{-1}$,
and $\Delta\Omega_{min}=500s^{-1}$. We examine control designs for various
target states in momentum space with real controls where $\Omega_{min}=-100$,
as well as with strictly non-negative real-valued control functions where
$\Omega_{min}=0$. Our computational results are described in the following
section. Table I contains a list of control synthesis problems specified by
different values of parameters $n$ and $\Delta\Omega_{min}$, and the resulting
time required to compute the ensemble controls. We observe that the
computation time is affected by the increase in complexity required to achieve
a higher momentum level $n$. All simulations were performed on a 2023 MacBook
Pro with Apple(R) M2 Pro Processor at 3.4 GHz and 16GB of memory.

TABLE I: Computation times for moment pulse design. $n$ in $C_{f}^{+}$ | $\Delta\Omega_{min}$ | Computation time (s)
---|---|---
1 | 0 | 866.4169
2 | 0 | 827.9643
3 | 0 | 1079.8969
4 | 0 | 1039.5165
1 | -100 | 700.5696
2 | -100 | 845.3960
3 | -100 | 1001.5283
4 | -100 | 1088.8991
Figure 2: Comparison of the performance of optimal square pulse sequences with
moment ensemble controls. For all target states, we design for $\delta=0.1$,
and use a time horizon of $T=3$. The desired momentum levels for these
computations are (from left to right) $C_{2}^{+}$, $C_{4}^{+}$, $C_{6}^{+}$,
and $C_{8}^{+}$. The plots show (from top to bottom): the final error achieved
for the square pulse and the moment ensemble method as a function of the
ensemble parameter $\varepsilon$; the optimal control function for strictly
non-negative values of $\Omega(t)$; and the optimal control when allowing
real-valued inputs. Note that the skew-rate limit is binding in the control
solutions computed for $n=4$ in the right most column, and this is what limits
the terminal state fidelity.

###  V-B Moment Ensemble and Square Pulse Design

To demonstrate the performance improvement of the Legendre moment dynamics
method with respect to the square pulse sequence obtained by solving problem
(7), we compare controls computed for an ensemble defined by $\delta=0.1$. The
square pulses are designed using a maximum momentum level $N^{+}=24$, and the
parameters that define the optimal pulses for target momentum states
$n=1,2,3,4$ are given in Table II. For control synthesis using the proposed
moment ensemble method, we use a moment order $N=20$ and truncate the momentum
level at $N^{+}=9$. The resulting controls are validated by applying them to
simulations of the entire ensemble with $\varepsilon\in[1-\delta,1+\delta]$
with the momentum levels truncated above $N^{+}=24$.

TABLE II: Optimal square pulse design parameters. $n$ | $\Omega_{1}/\omega_{r}$ | $\Omega_{2}/\omega_{r}$ | $\omega_{r}\tau_{1}$ | $\omega_{r}\tau_{2}$ | $\omega_{r}\tau_{3}$
---|---|---|---|---|---
1 | 3.9865 | 2.2849 | 0.4744 | 0.9427 | 0.4181
2 | 13.0036 | 9.5440 | 1.000 | 0.7190 | 1.000
3 | 32.4012 | 34.7591 | 0.1905 | 0.5523 | 0.1913
4 | 41.4215 | 41.4263 | 1.2653 | 0.8002 | 1.9952

The real-valued and positive ensemble controls and terminal state error
functions for the ensemble controls and square pulse sequences are compared in
Figure 2 for target momentum levels $n=1$, $2$, $3$, and $4$. The square pulse
sequence performs as expected in each case, with error that is negligible at
the nominal value of $\varepsilon=1$, but which quickly increases as
$\varepsilon$ diverges from unity. In contrast, we see that the terminal state
error for the simulation in which the controls obtained by solving problem
(14) are applied to the ensemble remains quite low for all values of the
ensemble parameter $\varepsilon\in[0.9,1.1]$ for target momentum levels $n=1$,
$2$, $3$, and still shows a significant performance advantage in the case of a
target state $C_{f}^{+}=C_{f,2n}^{+}$ with $n=4$.

The proposed method shows clear improvement with respect to all tested square
pulse sequences in overall performance for momentum transfer of the ensemble
of systems for positively-constrained controls as well as real-valued ones.
Observe that the amplitude and variation of the controls obtained for $n=4$
are greater than for the lower valued momentum state targets, which can be
explained by the need to manipulate higher frequency dynamics. Indeed, the
constraints on amplitude and variation for control function are binding at
many times during the optimization horizon, which necessarily limits the
degree to which the algorithm can meet the optimization objective. There are
inherent trade-offs between the constraint bound values, the optimization
horizon, and the complexity of the control task, which we have illustrated
here.

Figure 3: Comparison of control performance, including the expanded square
pulse design method for a collection of systems. The plots represent the
design parameters (from left to right): $\delta=0.1$ and $n=1$; $\delta=0.1$
and $n=4$; and$\delta=0.4$ and $n=1$.

The performance achieved by the proposed method is quantified in Table III, in
which the index $I_{e}$ defined in Equation (15) is given for the square pulse
controls and the positive and real valued ensemble controls as shown in Figure
2, for which the index values are denoted by $I_{e,sp}$, $I_{e,+}$ and
$I_{e,\mathbb{R}}$, respectively. As seen in the top row of Figure 2, the
results in Table III show that the positive valued ensemble controls result in
a significant decrease in terminal error with respect to the square pulse
sequence, and the real-valued ensemble controls are even better. We note that
in practice, the control represents the amplitude envelope of an optical
pulse, so $\Omega(t)$ must remain positive in the present setting.
Reformulation of the Raman-Nath equations and adjustment of the experimental
setting may enable real-valued control inputs, however.

TABLE III: Performance results for moment ensemble controls for $\delta=0.1$. $n$ | $I_{e,sp}$ | $I_{e,+}$ | $I_{e,\mathbb{R}}$
---|---|---|---
1 | 0.0141 | 0.0021 | 0.0004
2 | 0.0579 | 0.0012 | 0.0004
3 | 0.0515 | 0.0030 | 0.0006
4 | 0.1768 | 0.0608 | 0.0260

###  V-C Moment Ensemble and Ensemble Square Pulse Design

We also compare the ensemble control method to the square pulse controls
obtained by solving the optimization problem in Equation (8). The comparison
is made for three control design cases. First, the target state is $n=1$ with
$\delta=0.1$; the second case has target state $n=4$ with $\delta=0.1$; and
the third aims for $n=1$ with $\delta=0.4$. The parameter values obtained by
solving the square pulse design problem (8) using multiple $\varepsilon_{i}$
samples for $i=1,\ldots,m$ are given in Table IV. The results for these
scenarios are shown in Figure 3, and the values of the performance index
$I_{e}$ for the compared controls are shown in Table V. The table includes
values $I_{e,sp}$, $I_{e,+}$, $I_{e,\mathbb{R}}$, $I_{e,3}$, and $I_{e,10}$
for the nominal parameter square pulse, the positive ensemble control, the
real-valued ensemble control, the square pulse with $m=3$ samples, and the
square pulse with $m=10$ samples, respectively, for all three control design
cases described above.

TABLE IV: Ensemble Square pulse design parameters. $n$ | $\delta$ | $m$ | $\Omega_{1}/\omega_{r}$ | $\Omega_{2}/\omega_{r}$ | $\omega_{r}\tau_{1}$ | $\omega_{r}\tau_{2}$ | $\omega_{r}\tau_{3}$
---|---|---|---|---|---|---|---
1 | 0.1 | 3 | 4.6404 | 4.2278 | 0.1943 | 0.9618 | 0.5286
1 | 0.4 | 3 | 10.6354 | 6.6525 | 1.2038 | 3.5000 | 0.5420
4 | 0.1 | 3 | 64.3286 | 27.5580 | 0.0924 | 3.4724 | 0.6647
1 | 0.1 | 10 | 4.0446 | 4.2140 | 0.3769 | 0.6386 | 0.6767
1 | 0.4 | 10 | 60.0000 | 4.6034 | 0.0159 | 2.6005 | 0.5322
4 | 0.1 | 10 | 60.8008 | 30.2755 | 0.0967 | 3.4694 | 0.6798

The results show the limitations of the square pulse design, even when
multiple samples are used over the ensemble parameter domain. For the first
case shown at left in Figure 3, the ensemble square pulse approach (Problem
(8)) provides little improvement in relation to the nominal square pulse
(Problem (7)). The few degrees of freedom of the method and the small domain
$K$ limits the improvement gained by sampling in the ensemble space. In the
second and third cases, however, the improvement over using only one nominal
value of $\varepsilon$ is noticeable. Nonetheless, the ensemble moment method
exhibits significant improvements over the square pulse designs based on
sampling.

TABLE V: Performance results for the computations including the ensemble square pulse approach. $n$ | $\delta$ | $I_{e,sp}$ | $I_{e,+}$ | $I_{e,\mathbb{R}}$ | $I_{e,3}$ | $I_{e,10}$
---|---|---|---|---|---|---
1 | 0.1 | 0.0141 | 0.0021 | 0.0004 | 0.0127 | 0.0125
4 | 0.1 | 0.1768 | 0.0608 | 0.0260 | 0.0984 | 0.0965
1 | 0.4 | 0.2278 | 0.0030 | 0.0006 | 0.5884 | 0.1818

We conclude that directly sampling the ensemble parameter space $K$ for the
optimal control problem is not an effective means to compensate for
inhomogeneities in the experimental setting. In contrast, the proposed
ensemble control method promises to achieve improved and homogeneous
performance for an entire BEC ensemble in a matter wave interferometry
experiment.

##  VI CONCLUSION

We present a method for designing constrained optical pulses for optimal
matter-wave splitting of a Bose-Einstein condensate in the presence of
experimental inhomogeneities, which induce a parameter uncertainty within the
momentum evolution of the collection of atoms. An ensemble of systems is
parameterized by the factor that modulates the optical pulse, and Legendre
moments are used to represent the optimal control problem for an infinite
dimensional system as an optimization problem in moment space. The results
demonstrate that the fidelity achieved by the proposed method has advantages
over recently developed state-of-the-art optical pulse sequences. In
particular, the method achieves precise momentum transfer of the BEC to high
diffraction orders with inherent robustness to inhomogeneity in the effect of
the optical pulse on atoms in the sample, and yields continuously-varying
pulses that can be tuned to maximize the effectiveness of equipment in
experimental settings.

We expect that future research will extend our results to more detailed models
of BEC-splitting, as well as control protocols for additional experimental
steps for state preparation and measurement. The moment method can be extended
to incorporate more parameters to improve robustness to additional sources of
uncertainty and inhomogeneity. Our computational approach could be improved
and further analyzed to characterize the trade-offs in this setting between
truncations in the moment and quantum Hilbert spaces, constraint bound values,
pulse duration, and control performance. Based on the successful validation of
our approach in simulation, we expect that optimal pulses can be tested in
experimental settings to confirm anticipated improvements in performance of
matter-wave diffraction metrology.

## References

  * [1] Christiane P. Koch, Ugo Boscain, Tommaso Calarco, Gunther Dirr, Stefan Filipp, Steffen J. Glaser, Ronnie Kosloff, Simone Montangero, Thomas Schulte-Herbrüggen, Dominique Sugny, et al.  Quantum optimal control in quantum technologies. strategic report on current status, visions and goals for research in europe.  EPJ Quantum Technology, 9(1):19, 2022.
  * [2] Mary Clare Cassidy, Malcolm G Boshier, and Lee E Harrell.  Improved optical standing-wave beam splitters for dilute Bose–Einstein condensates.  Journal of Applied Physics, 130(19), 2021.
  * [3] B. Canuel, F. Leduc, D. Holleville, Alexandre Gauguet, J. Fils, A. Virdis, A. Clairon, N. Dimarcq, Ch J. Bordé, A. Landragin, et al.  Six-axis inertial sensor using cold-atom interferometry.  Physical Review Letters, 97(1):010402, 2006.
  * [4] P. A. Altin, M. T. Johnsson, Vladimir Negnevitsky, G. R. Dennis, Russell Paul Anderson, J. E. Debs, S. S. Szigeti, K. S. Hardman, Shayne Bennetts, G. D. McDonald, et al.  Precision atomic gravimeter based on bragg diffraction.  New Journal of Physics, 15(2):023009, 2013.
  * [5] Matt Jaffe, Victoria Xu, Philipp Haslinger, Holger Müller, and Paul Hamilton.  Efficient adiabatic spin-dependent kicks in an atom interferometer.  Physical Review Letters, 121(4):040402, 2018.
  * [6] Alexander Dunning, Rachel Gregory, James Bateman, Nathan Cooper, Matthew Himsworth, Jonathan A Jones, and Tim Freegarde.  Composite pulses for interferometry in a thermal cold atom cloud.  Physical Review A, 90(3):033608, 2014.
  * [7] Jack Saywell, Max Carey, Mohammad Belal, Ilya Kuprov, and Tim Freegarde.  Optimal control of Raman pulse sequences for atom interferometry.  Journal of Physics B: Atomic, Molecular and Optical Physics, 53(8):085006, 2020.
  * [8] Jr-Shin Li and Navin Khaneja.  Control of inhomogeneous quantum ensembles.  Physical review A, 73(3):030302, 2006.
  * [9] Shen Zeng and Frank Allgoewer.  A moment-based approach to ensemble controllability of linear systems.  Systems & Control Letters, 98:49–56, 2016.
  * [10] Yao-Chi Yu, Vignesh Narayanan, ShiNung Ching, and Jr-Shin Li.  Learning to control neurons using aggregated measurements.  In 2020 American Control Conference (ACC), pages 4028–4033. IEEE, 2020.
  * [11] Zaki Leghtas, Alain Sarlette, and Pierre Rouchon.  Adiabatic passage and ensemble control of quantum systems.  Journal of Physics B: Atomic, Molecular and Optical Physics, 44(15):154017, 2011.
  * [12] Anatoly Zlotnik, Raphael Nagao, István Z Kiss, and Jr-Shin Li.  Phase-selective entrainment of nonlinear oscillator ensembles.  Nature communications, 7(1):10788, 2016.
  * [13] Jr-Shin Li and Navin Khaneja.  Ensemble control of bloch equations.  IEEE Transactions on Automatic Control, 54(3):528–536, 2009.
  * [14] Xin Ning, André Luiz P. de Lima, and Jr-Shin Li.  NMR pulse design using moment dynamical systems.  In IEEE 61st Conference on Decision and Control (CDC), pages 5167–5172, 2022.
  * [15] Saijun Wu, Ying-Ju Wang, Quentin Diot, and Mara Prentiss.  Splitting matter waves using an optimized standing-wave light-pulse sequence.  Physical Review A, 71(4):043602, 2005.
  * [16] Anatoly Zlotnik and Shin Li.  Iterative ensemble control synthesis for bilinear systems.  In 51st IEEE Conference on Decision and Control (CDC), pages 3484–3489. IEEE, 2012.
  * [17] Vignesh Narayanan, Wei Zhang, and Jr-Shin Li.  Moment-based ensemble control.  arXiv preprint arXiv:2009.02646, 2020.
  * [18] Jr-Shin Li, Wei Zhang, and Yuan-Hung Kuan.  Moment quantization of inhomogeneous spin ensembles.  Annual Reviews in Control, 54:305–313, 2022.
  * [19] Alexander D. Cronin, Jörg Schmiedmayer, and David E. Pritchard.  Optics and interferometry with atoms and molecules.  Reviews of Modern Physics, 81(3):1051, 2009.
  * [20] Wei Xiong, Xuguang Yue, Zhongkai Wang, Xiaoji Zhou, and Xuzong Chen.  Manipulating the momentum state of a condensate by sequences of standing-wave pulses.  Physical Review A, 84(4):043616, 2011.
  * [21] Yao-Chi Yu, Vignesh Narayanan, and Jr-Shin Li.  Moment-based reinforcement learning for ensemble control.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–12, 2023.
  * [22] Shen Zeng.  Iterative optimal control syntheses illustrated on the Brockett integrator.  IFAC-PapersOnLine, 52(16):138–143, 2019.
  * [23] Minh Vu and Shen Zeng.  Iterative optimal control syntheses for nonlinear systems in constrained environments.  In 2020 American Control Conference (ACC), pages 1731–1736. IEEE, 2020.