# Symplectic eigenvalues of positive-semidefinite matrices and the trace
minimization theorem

Nguyen Thanh Son Institute of Mathematics $\&$ Centre for Advanced Analytics
and Predictive Sciences (CAAPS), University of Augsburg, Universitätsstraße
12a, 86159 Augsburg, Germany and Thai Nguyen University of Sciences, 24118
Thai Nguyen, Vietnam (ntson@tnus.edu.vn). Tatjana Stykel Institute of
Mathematics $\&$ Centre for Advanced Analytics and Predictive Sciences
(CAAPS), University of Augsburg, Universitätsstraße 12a, 86159 Augsburg,
Germany (tatjana.stykel@math.uni-augsburg.de).

###### Abstract

Symplectic eigenvalues are conventionally defined for symmetric positive-
definite matrices via Williamson’s diagonal form. Many properties of standard
eigenvalues, including the trace minimization theorem, are extended to the
case of symplectic eigenvalues. In this note, we will generalize Williamson’s
diagonal form for symmetric positive-definite matrices to the case of
symmetric positive-semidefinite matrices, which allows us to define symplectic
eigenvalues, and prove the trace minimization theorem in the new setting.

###### keywords:

Symplectic eigenvalues, positive-semidefinite matrices, trace minimization.

15A15, 15A18, 70G45

##  1 Introduction

Let us first recall that a matrix $S\in{\mathbb{R}^{2n\times 2n}}$ is referred
to as a _symplectic matrix_ if it satisfies the relation

$S^{T}J_{2n}S=J_{2n}\quad\mbox{with}\quad J_{2n}=\begin{bmatrix}0&I_{n}\\\
-I_{n}&0\end{bmatrix},$

where $I_{n}$ denotes the $n\times n$ identity matrix. The set of such
matrices forms a multiplicative group [9] and is denoted by
${\mathrm{Sp}(2n)}$. We also employ the set of rectangular symplectic matrices
defined as

${\mathrm{Sp}(2k,2n)}=\\{S\in\mathbb{R}^{2n\times
2k}\,:\,S^{T}J_{2n}S=J_{2k}\\}$

for some $k$ with $1\leq k\leq n$. This set was shown to be an embedded
Riemannian manifold of $\mathbb{R}^{2n\times 2k}$ [8, Proposition 3.1].

Derived from Williamson’s work [20], for a $2n\times 2n$ symmetric positive-
definite (spd) real matrix $A$, there exists a matrix $S\in{\mathrm{Sp}(2n)}$
such that

(1) $S^{T}AS=\begin{bmatrix}D&0\\\ 0&D\end{bmatrix},$

where $D=\mathrm{diag}(d_{1},\ldots,d_{n})$ is a diagonal matrix with positive
diagonal entries. The equality (1) looks very like the eigenvalue
decomposition of spd matrices except for the fact that the transformation
matrix $S$ is symplectic instead of orthogonal. The right-hand side of (1) is
termed as _Williamson’s diagonal form_ of $A$. The positive numbers
$d_{1},\ldots,d_{n}$ are referred to as _symplectic eigenvalues_. They form
the _symplectic spectrum_ of $A$. The symplectic spectrum of an spd matrix is
unique while the so-called _diagonalizing matrix_ $S$ in (1) is not unique.
The set of diagonalizing matrices for a given spd matrix was characterized via
a symplectic version of Autonne’s uniqueness theorem in [19, Theorem 3.5].

Symplectic eigenvalues find applications in quantum mechanics, optics, and
stability analysis of structured mechanical systems [11, 2, 15]. Symplectic
eigenvalues and eigenvectors, defined below, can be numerically computed using
a symplectic Lanczos method via the connection with so-called positive-
definite Hamiltonian matrices [1] or by solving a trace minimization problem
using a Riemannian optimization method [8, 19].

Many properties of symplectic eigenvalues for spd matrices have been
investigated. Besides trace minimization mentioned above, these are Cauchy’s
interlacing theorem [4], variational principles and Weyl inequality [6],
Schur-Horn theorem [5], Szegő-type theorem [3], and Lidskii’s (or the
majorization) theorem [13], just to name a few.

It is noteworthy that in the original work [20], $A$ is not necessarily spd
but solely symmetric. It would be a reasonable question to ask if it is still
possible to define symplectic eigenvalues and investigate their properties for
a more general case in which $A$ is symmetric positive-semidefinite (spsd). It
turns out that not every spsd matrix is appropriate for this purpose and even
in the case that symplectic eigenvalues are definable, some properties that
are already known for symplectic eigenvalues of spd matrices are nontrivially
extendable to the new setting. In this note, we will set the condition under
which an spsd matrix enjoys Williamson’s diagonal form which enables the
definition of symplectic eigenvalues. Moreover, we also prove the trace
minimization theorem for symplectic eigenvalues of such matrices because it
can not only be exploited to compute symplectic eigenpairs [19] but it also
provides a direct way to derive the majorization inequality similar to [11,
Theorem 1].

For ease of presentation, we now collect some necessary materials following
[7, 13, 19]. The Kronecker delta function will be denoted and understood as

$\delta_{ij}=\begin{cases}1,\ \mbox{ if }i=j,\\\ 0,\ \mbox{ if
}i\not=j.\end{cases}$

The superscript $\ \cdot^{T}$ applied to a matrix means its transpose. The set
of all $2n\times 2n$ spsd real matrices will be denoted by
$\mathcal{SPSD}(2n)$. Given a square matrix $M$, $\mathrm{tr}(M)$ denotes the
trace of $M$; the _kernel_ or the _null space_ of $M$ is the set $\ker
M=\\{x:Mx=0\\}$; a  _square root_ of $M$ is a square matrix $K$ of the same
size such that $K^{2}=M$. Furthermore,
$\mathrm{span}\\{v_{1},\ldots,v_{\ell}\\}$ stands for the subspace spanned by
the vectors $v_{1},\ldots,v_{\ell}$; applied to a matrix, it means the
subspace spanned by the columns of this matrix. Given
$A_{j}~{}\in~{}\mathbb{R}^{n_{j}\times n_{j}},j=1,\ldots,\ell,$
$\mathrm{diag}(A_{1},\ldots,A_{\ell})$ denotes the (block) diagonal matrix
with $A_{j}$ on the diagonal. A subset $U\subset\mathbb{R}^{2n}$ is said to be
_symplectically orthogonal to_ a subset $V\subset\mathbb{R}^{2n}$ if
$u^{T}J_{2n}v=0$ for all $u\in U,v\in V$. Two pairs of vectors $(u_{1},v_{1})$
and $(u_{2},v_{2})$ in $\mathbb{R}^{2n}\times\mathbb{R}^{2n}$ are said to be
_symplectically normalized_ 111It is slightly different from the definition in
[19]. if

$u_{i}^{T}J_{2n}v_{j}=\delta_{ij}\quad\mbox{and}\quad
u_{i}^{T}J_{2n}u_{j}=v_{i}^{T}J_{2n}v_{j}=0\quad\mbox{for}\quad i,j=1,2.$

The symplecticity of the matrix $[x_{1},\ \ldots,\ x_{2k}]$ is equivalent to
the fact that each two pairs $(x_{i},x_{k+i})$ and $(x_{j},x_{k+j})$ are
symplectically normalized for $i,j=1,\ldots,k$ and $i~{}\not=~{}j.$ It was
shown in [7, Theorem 1.15] that a symplectically normalized set of vectors are
linear independent. Therefore, when $k=n$, such vectors form a basis of
$\mathbb{R}^{2n}$. In that case, it is called a _symplectic basis_. A subspace
$\mathbb{U}\subset\mathbb{R}^{2n}$ is called a _symplectic subspace_ if for
any $x\in\mathbb{U}\backslash\\{0\\}$, there exists $y\in\mathbb{U}$ such that
$x^{T}J_{2n}y=1$. It turns out that we can always construct a symplectic basis
for a symplectic subspace [7, Theorem 1.15]. Conversely, a given subspace is
symplectic if it owns a symplectic basis as shown in Lemma 2.1 below. Two
nonzero vectors written as a matrix $[u,v]\in\mathbb{R}^{2n\times 2}$ form a
_symplectic eigenvector pair_ corresponding to the _symplectic eigenvalue_ $d$
of $A\in\mathcal{SPSD}(2n)$ if it holds

(2) $Au=dJ_{2n}v,\quad Av=-dJ_{2n}u.$

The set $([u,v],d)$ is then generally referred to as a _symplectic eigenpair_.
The existence of Williamson’s diagonal form for any given $2n\times 2n$ spd
matrix tells us that there exists a symplectic basis for $\mathbb{R}^{2n}$
consisting of its symplectic eigenvectors.

##  2 Williamson’s diagonal form of symmetric positive-semidefinite matrices

In this section, we will extend the definition of symplectic eigenvalues to
the case of spsd matrices by explicitly constructing Williamson’s diagonal
forms. For further discussion, we will need the following technical property.

######  Lemma 2.1.

If a subspace $\mathbb{W}\subset\mathbb{R}^{2n}$ owns a symplectic basis, then
it is a symplectic subspace.

######  Proof 2.2.

Assume that the columns of $[w_{1},\ldots,w_{m},w_{m+1},\ldots,w_{2m}]$ with
$1\leq m\leq n$ form a symplectic basis of $\mathbb{W}$. Then any
$x\in\mathbb{W}\backslash\\{0\\}$ can be represented via this basis as
$x=\sum_{j=1}^{m}(\lambda_{j}w_{j}+\lambda_{m+j}w_{m+j})$ with
$\sum_{j=1}^{2m}\lambda_{j}^{2}\not=0$. Setting

$y=\frac{1}{\sum_{j=1}^{2m}\lambda_{j}^{2}}\sum_{j=1}^{m}(\lambda_{j}w_{m+j}-\lambda_{m+j}w_{j}),$

we can verify that $x^{T}J_{2n}y=1$. This implies that the subspace
$\mathbb{W}$ is symplectic.

For a moment, we consider an spd matrix $A$ which has Williamson’s diagonal
form (1). Multiplying (1) with $S^{-T}=J_{2n}SJ_{2n}^{T}$ from the left and
with $S^{-1}=J_{2n}^{T}S^{T}J_{2n}$ from the right, we derive

$A=J_{2n}S\begin{bmatrix}D&0\\\ 0&D\end{bmatrix}(J_{2n}S)^{T}.$

Then, for $0\leq m<n$, we set

$\tilde{A}=J_{2n}S\begin{bmatrix}\tilde{D}_{m+1:n}&0\\\
0&\tilde{D}_{m+1:n}\end{bmatrix}(J_{2n}S)^{T},$

where
$\tilde{D}_{m+1:n}=\mathrm{diag}(0,\ldots,0,d_{m+1},\ldots,d_{n})\in\mathbb{R}^{n\times
n}$. We reverse the above multiplications to obtain

$\tilde{A}S=J_{2n}S\begin{bmatrix}0&-\tilde{D}_{m+1:n}\\\
\tilde{D}_{m+1:n}&0\end{bmatrix}$

and

(3) $S^{T}\tilde{A}S=\begin{bmatrix}\tilde{D}_{m+1:n}&0\\\
0&\tilde{D}_{m+1:n}\end{bmatrix}.$

Several facts can be drawn from the above construction. First, $\tilde{A}$ is
an spsd matrix of rank $2n-2m$. Second, if we write $S$ by its columns as
$S=[s_{1},\ldots,s_{2n}]$, then the kernel of $\tilde{A}$ is
$\mathrm{span}\\{s_{1},\ldots,s_{m},s_{n+1},\ldots,s_{n+m}\\}$. Moreover, it
is a symplectic subspace of dimension $2m$ due to Lemma 2.1. And third, it is
symplectically orthogonal to the subspace
$\mathrm{span}\\{s_{m+1},\ldots,s_{n},s_{n+m+1},\ldots,s_{2n}\\}$ which is
nothing else but the range of $\tilde{A}$.

Obviously, the decomposition (3) is like (1) and therefore it can be thought
of as Williamson’s diagonal form of the spsd matrix $\tilde{A}$. This gives us
a clear vision of how to construct Williamson’s diagonal form of spsd matrices
and under which condition this is possible. In practice, we will proceed
somewhat in the reverse direction. Indeed, given $A\in\mathcal{SPSD}(2n)$
whose kernel is a symplectic space. As any symplectic subspace must have even
dimension [7, Sect. 1.1], we assume that $\dim(\ker A)=2m$ with $0\leq m<n$.
Let $W=[W_{1}\ W_{2}]$ with $W_{1},\ W_{2}\in\mathbb{R}^{2n\times m}$ be a
symplectic matrix whose columns form a symplectic basis of $\ker A$. Note that
such a basis can be computed using the symplectic Gram-Schmidt process [7,
Theorem. 1.15]. We adapt the procedure proposed in [17] for constructing
Williamson’s diagonal form of spd matrices to the more general case where the
matrix $A$ is spsd as follows.

First, it is worth to recall that any spsd matrix has a unique spsd square
root, denoted by $A^{1/2}$, and that $\ker A^{1/2}=\ker A$ [12, Theorem
7.2.6]. Next, we will show that $\ker(A^{1/2}J_{2n}A^{1/2})=\ker A^{1/2}$. It
is trivial that $\ker(A^{1/2}J_{2n}A^{1/2})\supseteq\ker A^{1/2}$. Now,
consider any $y\in\mathbb{R}^{2n}$ such that $A^{1/2}J_{2n}A^{1/2}y=0$. Assume
that $J_{2n}A^{1/2}y\not=0$. Because $\ker A^{1/2}$ is a symplectic subspace,
there exists $z\in\ker A^{1/2}$ such that $z^{T}J_{2n}J_{2n}A^{1/2}y=1$ which
leads to $y^{T}A^{1/2}z=-1$. This contradicts to the fact that $z\in\ker
A^{1/2}$ and hence $J_{2n}A^{1/2}y=0$. As $J_{2n}$ is nonsingular, we deduce
that $A^{1/2}y=0$. This implies the desired equality. Observe moreover that
the matrix $A^{1/2}J_{2n}A^{1/2}$ is skew-symmetric. Then, this matrix can be
transformed into the real Schur form

(4) $Q^{T}A^{1/2}J_{2n}A^{1/2}Q=\mathrm{diag}\left(0_{2\times
2},\ldots,0_{2\times 2},\begin{bmatrix}0&d_{m+1}\\\
-d_{m+1}&0\end{bmatrix},\ldots,\begin{bmatrix}0&d_{n}\\\
-d_{n}&0\end{bmatrix}\right),$

where $0<d_{m+1}\leq\cdots\leq d_{n}$ and $Q$ is a real $2n\times 2n$
orthogonal matrix, see, e.g., [10, Theorem 7.4.1]. After that, we multiply
both sides of (4) with the permutation matrix
$P~{}=~{}[\mathrm{e}_{1},\mathrm{e}_{3},\ldots,\mathrm{e}_{2n-1},\mathrm{e}_{2},\mathrm{e}_{4},\ldots,\mathrm{e}_{2n}]$
from the right and with $P^{T}$ from the left, where $\mathrm{e}_{j}$ denotes
the $j$-th canonical unit vector in $\mathbb{R}^{2n}$. As a result, we obtain

$P^{T}Q^{T}A^{1/2}J_{2n}A^{1/2}QP=\begin{bmatrix}0&D\\\ -D&0\end{bmatrix},$

where $D=\mathrm{diag}(0,\ldots,0,d_{m+1},\ldots,d_{n})$. Finally, setting

(5)
$\tilde{S}=J_{2n}A^{1/2}QP\left[\begin{array}[]{ll}0_{m\times(n-m)}&0_{m\times(n-m)}\\\
0_{(n-m)\times(n-m)}&-\tilde{D}_{m+1:n}^{-1/2}\\\
0_{m\times(n-m)}&0_{m\times(n-m)}\\\
\tilde{D}_{m+1:n}^{-1/2}&0_{(n-m)\times(n-m)}\end{array}\right]=:[\tilde{S}_{1}\
\tilde{S}_{2}]\in\mathbb{R}^{2n\times(2n-2m)},$

where
$\tilde{D}^{-1/2}_{m+1:n}=\mathrm{diag}(1/\sqrt{d_{m+1}},\ldots,1/\sqrt{d_{n}})$
and $\tilde{S}_{1},\tilde{S}_{2}\in\mathbb{R}^{2n\times(n-m)}$, we define the
matrix $S=[W_{1}\ \tilde{S}_{1}\ W_{2}\ \tilde{S}_{2}]$. The following
proposition shows that $S$ is a symplectic matrix which diagonalizes $A$.

######  Proposition 2.3.

Given a matrix $A\in\mathcal{SPSD}(2n)$ whose kernel is a symplectic subspace
of $\mathbb{R}^{2n}$ of dimension $2m$. Then the matrix $S$ constructed as
above is symplectic and it holds

(6) $S^{T}AS=\begin{bmatrix}D&0\\\ 0&D\end{bmatrix},$

where $D=\mathrm{diag}(d_{1},\ldots,d_{n})$ with
$0=d_{1}=\cdots=d_{m}<d_{m+1}\leq\cdots\leq d_{n}$.

######  Proof 2.4.

It follows from $\ker A^{1/2}=\ker A=\mathrm{span}(W)$ that $A^{1/2}W=0$.
Thus, in view of (5), $W^{T}J_{2n}\tilde{S}=0$. Furthermore, taking into
account the fact that the matrices $W$ and $\tilde{S}$ are both symplectic, we
deduce that $S$ is symplectic.

The diagonal form (6) follows directly from the equalities $AW_{j}=0$,
$W_{j}^{T}A=0$, and
$\tilde{S}^{T}_{i}A\tilde{S}_{j}=\delta_{ij}\mathrm{diag}(d_{m+1},\ldots,d_{n})$
for $i,j=1,2$.

In view of the connection of $d_{j},j=1,\ldots,n,$ determined in (6) with the
standard eigenvalues of Hamiltonian matrix $J_{2n}A$, see [19, Proposition
3.2] and references therein, we can conclude that $d_{j},j=1,\ldots,n,$ are
unique. This fact gives rise to the following definition.

######  Definition 2.5.

For $A\in\mathcal{SPSD}(2n)$ with a symplectic kernel, the right-hand side of
(6) is called _Williamson’s diagonal form_ of $A$. The nonnegative diagonal
elements in this form $0\leq d_{1}\leq\cdots\leq d_{n}$ are called the
_symplectic eigenvalues_ of $A$.

It immediately follows from (6) that the symplecticity of the kernel of $A$ is
also a necessary condition for the existence of Williamson’s diagonal form of
$A$. The spsd matrices that do not satisfy this condition are not difficult to
find. Consider, for instance, $A=\mathrm{diag}(I_{n},0_{n\times n})$. Then
$\ker A=\mathrm{span}\\{\mathrm{e}_{n+1},\ldots,\mathrm{e}_{2n}\\}$. One can
directly verify that $x^{T}J_{2n}y=0$ for all $x,y\in\ker A$, which means that
$\ker A$ is an isotropic subspace [7, Section 1.2]. In fact, all spsd matrices
whose kernel is an isotropic subspace do not have Williamson’s diagonal form.

The observation that owning a symplectic kernel is a necessary and sufficient
condition for the existence of Williamson’s diagonal form of an spsd matrix
was previously made in [13, Remark 2.6]. A justification for this observation
can be drawn from [13, Proposition 5.2] where one of the symplectic
projections is uniquely determined by the columns of $W$. Our proof here is
however more constructive.

For the sake of convenience, in the rest of this paper, we restrict ourselves
to the case of spsd matrices that satisfy the symplectic kernel condition
mentioned in Proposition 2.3. Certainly, $\\{0\\}$ is itself a (trivial)
symplectic subspace. Therefore, the derived results hold for spd matrices too.

##  3 Trace minimization theorem

The trace minimization theorem for the symplectic eigenvalues of spd matrices
was established in [11, 4]. To the best of our knowledge, these are the only
sources that can be found in the literature. This theorem was exploited in
[19] to compute symplectic eigenpairs of spd matrices. A popular way of
justifying such a statement is to use minimax principles; see [18, 16] for
standard eigenvalues of Hermitian matrix pencils and [4] for symplectic
eigenvalues of spd matrices. In [4], the derived minimax principles however
invoked the reciprocals of symplectic eigenvalues. This approach is certainly
inapplicable to the case of spsd matrices. In this section, we will prove the
trace minimization theorem for symplectic eigenvalues in a more general
setting where $A$ can be spsd. Note that restricted to the case of spd
matrices, our proof differs from the ones presented in [11, 4]. For
completeness, we restate the trace minimization theorem in a more general form
in which the symplectic eigenvalues can be zero.

######  Theorem 3.1.

Let $A\in\mathcal{SPSD}(2n)$ have symplectic eigenvalues $0\leq
d_{1}\leq\cdots\leq d_{n}$. Then for any $k$ with $1\leq k\leq n$, it holds
that

$\min_{X\in{\mathrm{Sp}(2k,2n)}}\mathrm{tr}(X^{T}AX)=2\sum_{j=1}^{k}d_{j}.$

For a proof of this result, we need some preparation. We start with
reformulating the definition of symplectic eigenvalues and eigenvectors.
Recalling (2), $[u,v]\in\mathbb{R}^{2n\times 2}$ is a symplectic eigenvector
pair corresponding to the symplectic eigenvalue $d$ of $A$ if and only if one
of the following four equivalent conditions holds:

(7) $\displaystyle\begin{split}&\begin{bmatrix}A&0\\\
0&A\end{bmatrix}\begin{bmatrix}u\\\ v\end{bmatrix}=d\begin{bmatrix}\ \
0&J_{2n}\\\ -J_{2n}&0\end{bmatrix}\begin{bmatrix}u\\\ v\end{bmatrix},\\\
&\begin{bmatrix}A&0\\\ 0&A\end{bmatrix}\begin{bmatrix}v\\\
u\end{bmatrix}=-d\begin{bmatrix}\ \ 0&J_{2n}\\\
-J_{2n}&0\end{bmatrix}\begin{bmatrix}v\\\ u\end{bmatrix},\\\
&\begin{bmatrix}A&0\\\ 0&A\end{bmatrix}\begin{bmatrix}\ \ v\\\
-u\end{bmatrix}=d\begin{bmatrix}\ \ 0&J_{2n}\\\
-J_{2n}&0\end{bmatrix}\begin{bmatrix}\ \ v\\\ -u\end{bmatrix},\\\
&\begin{bmatrix}A&0\\\ 0&A\end{bmatrix}\begin{bmatrix}-u\\\ \ \
v\end{bmatrix}=-d\begin{bmatrix}\ \ 0&J_{2n}\\\
-J_{2n}&0\end{bmatrix}\begin{bmatrix}-u\\\ \ \ v\end{bmatrix}.\end{split}$

This fact suggests constructing the matrix pencil

(8) $\mathcal{A}-\lambda\mathcal{J}:=\begin{bmatrix}A&0\\\
0&A\end{bmatrix}-\lambda\begin{bmatrix}\ \ 0&J_{2n}\\\ -J_{2n}&0\end{bmatrix}$

with the following properties.

######  Lemma 3.2.

Let $A\in\mathcal{SPSD}(2n)$ and $\mathcal{A}-\lambda\mathcal{J}$ be the
matrix pencil as in (8). Then the following statements hold:

  1. 1.

The real number $d$ is a symplectic eigenvalue of $A$ if and only if $d$ and
$-d$ are eigenvalues of multiplicities two of the matrix pencil
$\mathcal{A}-\lambda\mathcal{J}$.

  2. 2.

The matrix $\mathcal{A}$ is spsd and the matrix $\mathcal{J}$ is symmetric and
nonsingular.

  3. 3.

The matrix $\mathcal{J}$ possesses only $1$ and $-1$ as its eigenvalues.
Moreover, both of them are of multiplicity $2n$.

######  Proof 3.3.

The first statement follows immediately from (7) while the second one is
straightforward. To justify the last one, let $X_{1}$ and $X_{2}$ be any two
basis matrices of $\mathbb{R}^{n}$. Then it is direct to show that the
matrices

$\begin{bmatrix}X_{1}&\ \ 0\\\ 0&\ \ X_{2}\\\ 0&-X_{2}\\\ X_{1}&\ \
0\end{bmatrix},\quad\begin{bmatrix}\ \ X_{1}&0\\\ \ \ 0&X_{2}\\\ \ \
0&X_{2}\\\ -X_{1}&0\end{bmatrix}$

consist of eigenvectors of $\mathcal{J}$ associated with eigenvalues $1$ and
$-1$, respectively.

The relations in (7) together with Lemma 3.2 establish a connection between
symplectic eigenpairs of $A$ and eigenpairs of the matrix pencil (8). Hence,
the result on trace minimization of positive-semidefinite matrix pencils
derived in [14] is here helpful. Let us first collect the necessary facts from
there. Given two $\ell\times\ell$ symmetric real matrices $A$ and $B$. The
matrix pencil $A-\lambda B$ is said to be _positive-semidefinite_ if there is
a number $\lambda_{0}\in\mathbb{R}$ such that $A-\lambda_{0}B$ is positive-
semidefinite. Such constant $\lambda_{0}$ is called a _definitizing shift_.
Obviously, the matrix pencil $\mathcal{A}-\lambda\mathcal{J}$ in (8) is
positive-semidefinite with a shift $\lambda_{0}=0$.

Furthermore, it was shown in [14, Proposition 4.1] that the eigenvalues of the
positive-semidefinite matrix pencil $A-\lambda B$ with nonsingluar $B$ are
real and the number of non-positive (resp. non-negative) eigenvalues of
$A-\lambda B$ is the same as the number of negative (resp. positive)
eigenvalues of $B$. In other words, if $q$ and $p$ with $q+p=\ell$ is the
number of negative and positive eigenvalues of $B$, respectively, then there
are $q$ non-positive eigenvalues of $A-\lambda B$, namely,
$\alpha_{q}^{-}\leq\cdots\leq\alpha_{1}^{-}\leq 0$, and $p$ non-negative
eigenvalues of $A-\lambda B$, namely,
$0\leq\alpha_{1}^{+}\leq\cdots\leq\alpha_{p}^{+}$. Thus, they can overall be
arranged as

$\alpha_{q}^{-}\leq\cdots\leq\alpha_{1}^{-}\leq\alpha_{1}^{+}\leq\cdots\leq\alpha_{p}^{+}.$

The following theorem provides an important result on trace minimization for
positive-semidefinite matrix pencils.

######  Theorem 3.4.

[14, Theorem 3.1] Let $A,\,B\in\mathbb{R}^{\ell\times\ell}$ be symmetric and
let $B$ be nonsingular with $p$ positive and $q$ negative eigenvalues. Assume
moreover that the matrix pencil $A-\lambda B$ is positive-semidefinite. Let
$p_{1}$ and $q_{1}$ be two integers such that $0\leq p_{1}\leq p$ and $0\leq
q_{1}\leq q$. Then, the function

(9) $\mathbb{R}^{\ell\times(p_{1}+q_{1})}\ni\
X\mapsto\mathrm{tr}(X^{T}AX)\in\mathbb{R},$

subjected to the condition

(10) $X^{T}BX=\mathrm{diag}(I_{p_{1}},-I_{q_{1}}),$

is bounded from below by

(11) $\sum_{j=1}^{p_{1}}\alpha_{j}^{+}-\sum_{j=1}^{q_{1}}\alpha_{j}^{-}.$

If, additionally, there exists a matrix
$X_{0}\,\in\mathbb{R}^{\ell\times(p_{1}+q_{1})}$ satisfying (10) and
consisting of eigenvectors of the matrix pencil $A-\lambda B$ associated with
the eigenvalues
$\alpha_{1}^{+},\ldots,\alpha_{p_{1}}^{+},\alpha_{1}^{-},\ldots,\alpha_{q_{1}}^{-}$,
then the lower bound (11) becomes the minimal value of the function (9) which
is reached at $X_{0}$.

Now, we are ready to prove the trace minimization Theorem 3.1 for symplectic
eigenvalues of spsd matrices.

Proof of Theorem 3.1 Consider the matrix pencil
$\mathcal{A}-\lambda\mathcal{J}$ as in (8). By Proposition 2.3 and Lemma 3.2,
the eigenvalues of $\mathcal{A}-\lambda\mathcal{J}$, sorted in the
nondecreasing order, are

$-d_{n}=-d_{n}\leq\cdots\leq-d_{1}=-d_{1}\leq d_{1}=d_{1}\leq\cdots\leq
d_{n}=d_{n},$

where $d_{1},\ldots,d_{n}$ are the symplectic eigenvalues of $A$. Thus, an
application of Theorem 3.4 to the matrix pencil
$\mathcal{A}-\lambda\mathcal{J}$ for $p_{1}=q_{1}=2k$ yields

(12) $\mathrm{tr}(\mathcal{X}^{T}\mathcal{A}\mathcal{X})\geq
4\sum_{j=1}^{k}d_{j}$

for all $\mathcal{X}\in\mathbb{R}^{4n\times 4k}$ satisfying the constraint

(13) $\mathcal{X}^{T}\mathcal{J}\mathcal{X}=\mathrm{diag}(I_{2k},-I_{2k}).$

Next, for any $S\in{\mathbb{R}^{2n\times 2k}}$, we construct the matrix

(14) $\mathcal{X}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{ll}S&SJ_{2k}^{T}\\\
SJ_{2k}^{T}&S\end{array}\right]\in\mathbb{R}^{4n\times 4k}.$

This matrix fulfills

$\displaystyle\mathrm{tr}(\mathcal{X}^{T}\mathcal{A}\mathcal{X})$
$\displaystyle=\mathrm{tr}(S^{T}AS+J_{2k}S^{T}ASJ_{2k}^{T})$
$\displaystyle=\mathrm{tr}(S^{T}AS+S^{T}ASJ_{2k}^{T}J_{2k})$ (15)
$\displaystyle=2\,\mathrm{tr}(S^{T}AS).$

Moreover, taking $S\in{\mathrm{Sp}(2k,2n)}$, one can verify by a direct
calculation that the matrix $\mathcal{X}$ in (14) satisfies (13). Then, using
(15) and (12), we obtain that

(18) $\displaystyle\min_{S\in{\mathrm{Sp}(2k,2n)}}\mathrm{tr}(S^{T}AS)$
$\displaystyle=\frac{1}{2}\min_{\footnotesize{\begin{array}[]{ll}&\mathcal{X}\mbox{
as in \eqref{eq:proof3}},\\\
&S\in{\mathrm{Sp}(2k,2n)}\end{array}}}\mathrm{tr}(\mathcal{X}^{T}\mathcal{A}\mathcal{X})$
(19)
$\displaystyle\geq\frac{1}{2}\min_{\tilde{\mathcal{X}}^{T}\mathcal{J}\tilde{\mathcal{X}}=\mathrm{diag}(I_{2k},-I_{2k})}\mathrm{tr}(\tilde{\mathcal{X}}^{T}\mathcal{A}\tilde{\mathcal{X}})$
$\displaystyle\geq 2\sum_{j=1}^{k}d_{j},$

in which (19) is due to the fact that the constraint set in the right-hand
side of (18) is a subset of that in (19).

Finally, we choose $S=[u_{1},\ldots,u_{k},v_{1},\ldots,v_{k}]$ consisting of
(symplectically) normalized symplectic eigenvector pairs $[u_{j},v_{j}]$
associated with the symplectic eigenvalue $d_{j}$ of $A$ for $j=1,\ldots,k$.
In this case, we have $\mathrm{tr}(S^{T}AS)=2\sum_{j=1}^{k}d_{j}$ which
completes the proof. $\square$

######  Remark 3.5.

For the computation of the smallest symplectic eigenvalues of the spsd matrix
$A$, a direct use of the symplectic Lanczos procedure [1] will not work
because it involves the inverse of $A$ which does not exist. On the other
side, as the theoretical results in [19] still hold for spsd matrices, Theorem
3.1 plays a key role in enabling the use of the numerical method proposed
there for computing the required symplectic eigenpairs.

## References

  * [1] P. Amodio, On the computation of few eigenvalues of positive definite Hamiltonian matrices, Future Generation Comput. Syst., 22 (2006), pp. 403–411, https://doi.org/10.1016/j.future.2004.11.027.
  * [2] P. Benner, H. Fassbender, and M. Stoll, Solving large-scale quadratic eigenvalue problems with Hamiltonian eigenstructure using a structure-preserving Krylov subspace method, Electron. Trans. Numer. Anal., 29 (2008), pp. 212–229.
  * [3] R. Bhatia and T. Jain, A Szegő type theorem and distribution of symplectic eigenvalues, J. Spectr. Theory, 11 (2001), pp. 1369–1389, https://doi.org/10.4171/JST/377.
  * [4] R. Bhatia and T. Jain, On the symplectic eigenvalues of positive definite matrices, J. Math. Phys., 56 (2015), p. 112201, https://doi.org/10.1063/1.4935852.
  * [5] R. Bhatia and T. Jain, A Schur-Horn theorem for symplectic eigenvalues, Linear Algebra Appl., 599 (2020), pp. 133–139, https://doi.org/10.1016/j.laa.2020.04.005.
  * [6] R. Bhatia and T. Jain, Variational principles for symplectic eigenvalues, Canad. Math. Bull., 64 (2020), pp. 553–559, https://doi.org/10.4153/S0008439520000648.
  * [7] M. de Gosson, Symplectic Geometry and Quantum Mechanics, Advances in Partial Differential Equations, Birkhäuser, Basel, 2006.
  * [8] B. Gao, N. T. Son, P.-A. Absil, and T. Stykel, Riemannian optimization on the symplectic Stiefel manifold, SIAM J. Optim., 31 (2021), pp. 1546–1575, https://doi.org/10.1137/20M1348522.
  * [9] S. Fiori, Solving minimal-distance problems over the manifold of real-symplectic matrices, SIAM J. Matrix Anal. Appl., 32 (2011), pp. 938–968, https://doi.org/10.1137/100817115.
  * [10] G. Golub and C. V. Loan, Matrix Computations. 4th ed, The Johns Hopkins University Press, Baltimore, London, 2013.
  * [11] T. Hiroshima, Additivity and multiplicativity properties of some Gaussian channels for Gaussian inputs, Phys. Rev. A, 73 (2006), p. 012330, https://doi.org/10.1103/PhysRevA.73.012330.
  * [12] R. A. Horn and C. R. Johnson, Matrix Analysis. 2th ed, Cambridge University Press, New York, NY, 2013.
  * [13] T. Jain and H. K. Mishra, Derivatives of symplectic eigenvalues and a Lidskii type theorem, Canad. J. Math., 74 (2022), pp. 457–485, https://doi.org/10.4153/S0008414X2000084X.
  * [14] J. Kovač-Striko and K. Veselić, Trace minimization and definiteness of symmetric pencils, Linear Algebra Appl., 216 (1995), pp. 139–158, https://doi.org/10.1016/0024-3795(93)00126-K.
  * [15] M. Krbek, T. Tyc, and J. Vlach, Inequalities for quantum marginal problems with continuous variables, J. Math. Phys., 55 (2014), p. 062201, https://doi.org/10.1063/1.4880198.
  * [16] I. Nakić and K. Veselić, Wielandt and Ky-Fan theorem for matrix pairs, Linear Algebra Appl., 369 (2003), pp. 77–93, https://doi.org/10.1016/S0024-3795(02)00733-4.
  * [17] K. Parthasarathy, The symplectry group of Gaussian states in $L^{2}(\mathbb{R}^{n})$. In Shiryaev, A., Varadhan, S., Presman, E. (eds), Prokhorov and Contemporary Probability Theory, Springer Proceedings in Mathematics and Statistics 33, Berlin Heidelberg, 2013, pp. 349–369, https://doi.org/10.1007/978-3-642-33549-5_21.
  * [18] A. H. Sameh, and J. A. Wisniewski, A trace minimization algorithm for generalized eigenvalue problem, SIAM J. Numer. Anal., 19 (1982), pp. 1243–1259, https://doi.org/10.1137/0719089.
  * [19] N. T. Son, P.-A. Absil, B. Gao, and T. Stykel, Computing symplectic eigenpairs of symmetric positive-definite matrices via trace minimization and Riemannian optimization, SIAM J. Matrix Anal. Appl., 42 (2021), pp. 1732–1757, https://doi.org/10.1137/21M1390621.
  * [20] J. Williamson, On the algebraic problem concerning the normal forms of linear dynamical systems, Amer. J. Math., 58 (1936), pp. 141–163, https://doi.org/10.2307/2371062.


# On the measurement of the dispersion relation by a radar and the implication
on the current retrieval

Susanne Støle-Hentschel 111School of Engineering, Tel Aviv University, Israel
222Universit’e Paris-Saclay, ENS Paris-Saclay, Centre Borelli, Gif-sur-
Yvettes, France Benjamin Keeler Smeltzer 333Department of Energy and Process
Engineering, Norwegian University of Science and Technology, N-7491,
Trondheim, Norway 444SINTEF Ocean, Marinteknisk senter, N-7052 Trondheim,
Norway and Yaron Toledo††footnotemark:

###### Abstract

The work analyses the error in current retrievals from images of marine
radars. The study is based on simulations of waves interacting with a shear
current. The measured dispersion is related to the underlying wavenumber-
dependent effective current. The highest tested radar antenna (H=45 m) with
vertical polarization performed the best. For that case the root mean square
error was at most 0.05 ms-1 above the one for the simulated wave field without
imaging mechanism. The observation time of 20 minutes was compared to shorter
windows. Depending on the needed accuracy, the time may be reduced to five
minutes, associated with an loss of accuracy below 12%. The study shows the
error of the current reconstruction depends on the shape of the profile and
varies considerably from realization to realization.

##  1 Introduction

Currents transport sediments as well as nutrients and pollutants. For studying
their distribution in the ocean, current profile along the vertical axis is
needed. While it can be measured well at some depth, e.g. by ADCPs, it is
difficult to obtain the current profile accurately in the uppermost layer
where the water particles are affected by the currents and waves
simultaneously (Davis et al.,, 1981). High frequency (HF) radars are commonly
used to get current estimates close to the surface. Given that the current has
a boundary layer underneath the surface, multiple measurements at different
wavelengths would be required to define the profile. Although Ivonin et al.,
(2004) have shown that very high frequency (VHF) radars may be used to
retrieve currents at three different depth more points would be necessary for
an accurate description of boundary layer.

Marine X-band radars are capable of observing ocean waves of a range of
wavelengths, typically from 30 meters up to hundred meters and more. Young et
al., (1985) was among the first to describe how these images could be used to
extract information on waves and currents by applying a multidimensional
Fourier analysis. Although Young et al., (1985) acknowledged the work of
Stewart and Joy, (1974) resulting in a wavenumber-dependent distortion of the
dispersion relation, for many years the radar community only considered a bulk
current. In recent years there have been multiple attempts for obtaining more
detailed measurements of the current in the upper ocean layer from X-band
radar images, amongst others, Lund et al., (2018) and Campana et al., (2016).
In their work current profiles are inverted from radar images and compared to
ADCPS measurements. As pointed out by Smeltzer and Ellingsen, (2021), a
validation of the current estimates of radar images by ADCPs is questionable.
The deviation between the ADCP and the radar estimate is in the order of
magnitude of the measured current. A slightly similar shear profile comparison
presented in Lund et al., (2020). It did not only use ADCPs for comparison but
a combination of multiple instruments but the problem persists with a ground
truth current given for a time period of 3.5 hours. In addition to the
different time scales, the reference systems for both instruments differ.
While the Doppler-like shift measured by a radar measures the Lagrangian
current (Pizzo et al.,, 2023), ADCPs measure the Eulerian current.

inverting the current from the measured Doppler shifts is an ill posed problem
(Smeltzer and Ellingsen,, 2021). The suggested solutions depend on tunable
parameters like the smoothness of the final profile Campana et al., (2016). As
will be shown below, the selection valid Doppler estimates of the different
wavenumber components has an impact that can not be neglected.

This work aims to quantify errors in shear current profiles estimates by
marine radars. The error results from two sources: 1. The accuracy by which
the radar measures the dispersion relation. 2. The reconstruction of the
underlying shear current. In this study simulations are used to construct
artificial radar images. The a priori knowledge of the exact shear profiles
provides the errors resulting from both sources separately. A range of
different conditions of the radar settings and sea states are compared.
Finally the implications of the measured error on a current inversion are
illustrated by employing the inversion method presented by Smeltzer et al.,
(2019).

##  2 Simulated radar images of waves on a shear current

###  2.1 Simulated wave fields

Table 1: Description of the sea state Name | Parameter | value
---|---|---
Significant wave height | $H_{s}$ | 2.0 m
Peak wave number | $k_{p}$ | 0.073 rad m-1
Peak wave length | $\lambda_{p}$ | 87 m
Peak period | $T_{p}$ | 7.5 s
Directional spreading555according to Mitsuyasu et al., (1975) | $s_{\max}$ | 10, 30, 70
Spatial grid size | $\Delta x=\Delta y$ | 7.5 m
Temporal grid size | $d_{t}$ | 1 s
Covered time period | $T$ | 20 min
Covered patch size | $X=Y$ | 500 m
Peak enhancement factor | $\gamma$ | 3.3
Shear current | $U(z)$ | $\exp(mz)+0.05$
Exponential factor | $C$ | 0.5, 0.2
Current direction | $\psi$ | $30^{\circ}$
Wave direction | $\theta$ | $90^{\circ}$
Water depth | $h$ | 1000 m

Simulated wave fields are based on directional JONSWAP spectra and a shear
current according to the description provided in Table 1. It is assumed that
non-linear features of waves have no major effect on the radar images and the
wave simulations are purely linear. The directional spreading was implemented
following Mitsuyasu et al., (1975). For most investigations the current
profile $U(z)=\exp(0.5z)+0.05$ was employed. Part of the results are also
shown for a exponential profile with a milder shear ($U(z)=\exp(0.2z)+0.05$).

###  2.2 Incorporation of the shear current

The interaction of wave motion and shear currents can be expressed by the
Rayleigh equation, here using the formulation by Ellingsen and Li, (2017)

$w^{\prime\prime}(z)-k^{2}w(z)=\frac{\mathbf{k}\cdot\mathbf{U}^{\prime\prime}(z)}{\mathbf{k}\cdot\mathbf{U}(z)-kc}w(z)$
(1)

where $w$ is the vertical velocity , z, the vertical axis with 0 in the mean
water level and the negative axis towards the sea bottom. $c$ denotes the wave
celerity and $\mathbf{U}$, the two-dimensional shear current.

The boundary condition at $z=\eta(0)$ is

$\left(\mathbf{k\cdot U}_{0}-kc\right)^{2}w^{\prime}(0)-\left(\mathbf{k\cdot
U}_{0}^{\prime}\left(\mathbf{k\cdot
U}_{0}-kc\right)+\mathrm{g}k^{2}\right)w(0)=0,$ (2)

At the sea bottom the vertical velocity vanishes $w(-h)=0$.

The celerity can be retrieved by the shooting method and the solution is very
sensitive to boundary values at the sea surface. For a wide range of profiles,
it is possible to use an approximate solution following the perturbation
approach by Stewart and Joy, (1974) or higher order versions, e.g. Kirby and
Chen, (1989). To first order, the celerity is that of waves without current.
The second order solution adds the effective current,

$U_{\mathrm{eff}}(k)=2k\int_{-\infty}^{0}U(z)\mathrm{e}^{(2kz)}.$ (3)

that can be applied to adapt the dispersion relation

$\omega=\sqrt{k\mathrm{g}\tanh(kh)}+kU_{\mathrm{eff}}(k)\cos(\phi)$ (4)

where $\omega$ is the angular frequency, $k$ the wavenumber, $\mathrm{g}$ the
gravitational constant, h the water depth and $\phi=\theta-\psi$ the angle
between the current $\psi$ and the wavenumber vector $\theta$ (compare Ivonin
et al., (2004)).

Ellingsen and Li, (2017) introduced the dimensionless depth-averaged shear

$\delta(\mathbf{k})=\int_{-h}^{0}\frac{\mathbf{k\cdot
U}^{\prime}(z)\sinh\left(2k(z+h)\right)}{kc_{0}\sinh(2kh)}dz$ (5)

that can be used as measure to judge the applicability of the effective
current according to Stewart and Joy, (1974). In all simulations used
$|\delta(\mathbf{k})|<0.08$. The largest $\delta$ for the presented cases
occurs for the shear current $U(z)=\exp(0.5z)+0.05$. The approximated
$\omega(k)$ based on the effective current of Stewart and Joy, (1974) deviates
at most by 0.0022 $\mathrm{rads}^{-1}$ from the solution of the Rayleigh
equation. The error from the approximation is hence acceptable.

###  2.3 Simulated radar images

The radar image is based on a simplified version of the Radar equation (e.g.
Rees, (2013))

$P_{r}=\frac{P_{t}G_{t}\sigma A_{e}}{(4\pi)^{2}r^{4}}$ (6)

where $P_{r}$ and $P_{t}$ are the received and transmitted power respectively,
$G_{t}$ is the gain, $A_{e}$ the effective area (only depending on the radar),
$\sigma$ the radar cross section and $r$ the range. If an uncalibrated radar
is used and we do not take into consideration the variation in areal coverage
we are left with

$P_{r}\propto\frac{\sigma^{0}}{r^{4}}$ (7)

According to Rees, (2013), Eq. 3.51, the backscattering coefficient
$\sigma^{0}$ for the different polarizations $pp$ is

$\sigma_{pp}^{0}=4k^{4}L^{2}(\Delta
h)^{2}\cos^{4}(\theta_{l})\left|f_{pp}(\theta)\right|^{2}\exp\left(-k^{2}L^{2}\sin^{2}(\theta_{l})\right)$
(8)

where $k$ denotes the wavenumber of the electromagnetic wave while $\Delta h$
and $L$ denote characteristic scales to define the surface roughness. Herein,
it is assumed that these will be constant for a given sea state. Hence the
backscattering coefficient only depends on the local incidence angle
$\theta_{l}$. The latter is defined as the angle between the surface normal
and the beam towards the radar (given below).

$\displaystyle f_{HH}(\theta)$
$\displaystyle=\frac{\cos(\theta_{l})-\sqrt{\varepsilon_{r}-\sin^{2}(\theta_{l})}}{\cos(\theta_{l})+\sqrt{\varepsilon_{r}-\sin^{2}(\theta_{l})}},$
(9) $\displaystyle f_{VV}(\theta)$
$\displaystyle=(\varepsilon_{r}-1)\frac{\sin^{2}(\theta_{l})-\varepsilon_{r}\left(1+\sin^{2}(\theta_{l})\right)}{\left(\varepsilon_{r}\cos(\theta_{l})+\sqrt{\varepsilon_{r}-\sin^{2}(\theta_{l})}\right)^{2}},$
(10)

and $k$ is the wave number of the electromagnetic wave, $\Delta h$ the height
variation of the roughness and $L$ the correlation length of the roughness of
the observed surface.

Depending on the water temperature, $\varepsilon_{r}$ is in the range of
$[70^{\circ};80^{\circ}]$. For high grazing angles,

$\displaystyle\left|f_{HH}\right|\approx 1$ (11)

and hence $\sigma_{pp}^{0}\propto\cos^{4}(\theta_{l})$. For vertical
polarization the local incidence angle changes the backscattering coefficient.
The above equation may be approximated by the following:

$\displaystyle\left|f_{VV}\right|\approx\left|\frac{-1-\sin^{2}(\theta_{l})}{\cos(\theta_{l})}\right|$
(12)

The final simplistic radar image is hence given by

$P_{r}\propto\left\\{\begin{matrix}\frac{\cos^{4}(\theta_{l})}{r^{4}}&\mbox{for~{}HH}\\\
\frac{\cos^{2}(\theta_{l})(-1-\sin^{2}(\theta_{l}))^{2}}{r^{4}}&\mbox{for~{}VV}\end{matrix}\right.$
(13)

The amplitude of the signal is the square root of the received power.

The local incidence angle is calculated from the tilt of the surface:

$\cos\left(\theta_{l}\right)=\frac{r\frac{\partial\eta}{\partial
x}+H-\eta}{\sqrt{1+\left(\frac{\partial\eta}{\partial
x}\right)^{2}+\left(\frac{\partial\eta}{\partial
y}\right)^{2}}\sqrt{r^{2}+\left(H-\eta\right)^{2}}}$ (14)

In addition, geometric shadowing as formulated in Nieto Borge et al., (2004)
is applied.

##  3 Extraction of wavenumber dependent Doppler shift

We assume the linear dispersion relation for wave frequency
$\omega_{\mathrm{DR}}$ as a function of wavenumber $\mathbf{k}$ may be
expressed as:

$\omega_{\mathrm{DR}}(\mathbf{k})=\omega_{0}(k)+\mathbf{k}\cdot\mathbf{U}_{\mathrm{eff}}(k),$
(15)

where the $\omega_{0}(k)$ is the dispersion relation in quiescent waters,
$k=|\mathbf{k}|$, $\mathbf{U}_{\mathrm{eff}}(k)$ is the wavenumber-dependent
Doppler shift velocity due to the background current profile Stewart and Joy,
(1974).

The input data is a three-dimensional array or “cube” of sea surface images as
a function of time, denoted as $I(x,y,t)$, where $x,y$ correspond to the
spatial dimensions of the images, and $t$ is time. The physical extent of the
images (assumed to be square) is length $\Delta x$, and the duration of the
time series $T$. The images may in this context refer to either optical
imagery obtained from an aerial platform for example, or radar backscatter
images. The first step in extracting the Doppler shift velocities involves
computing the frequency-wavenumber spectrum from a fast Fourier transform:

$P(\mathbf{k},\omega)=|\mathrm{FFT}\\{I(x,y,t)\\}|^{2}.$ (16)

The spectrum is peaked along the dispersion relation shell. Due to the
symmetries of the Fourier transform of a signal comprised of real numbers,
there are two wave frequencies $\omega_{\pm}$ associated with each unique
wavenumber: $\omega_{+}(\mathbf{k})=\omega_{\mathrm{DR}}(\mathbf{k})$, and
$\omega_{-}(\mathbf{k})=-\omega_{\mathrm{DR}}(-\mathbf{k})$.

In the following steps, we wish to extract $\mathbf{U}_{\mathrm{eff}}(k)$ from
the measured spectrum $P$ over a range of wavenumbers. To do so, we consider
wavenumber bins with center values $k_{i}$ and half-width $\delta k$. A lower
limit on $k_{i}$ of 6-8 multiples of the wavenumber resolution $2\pi/\Delta x$
is typically chosen. On the high wavenumber end, the signal-to-noise ratio
typically limits the Doppler shift extraction, and some examination of the
data is required in choosing the upper wavenumber limit. For each center
wavenumber value $k_{i}$, we define a masked spectrum $F_{i}$ where energy
outside the wavenumber bin is set to zero, i.e.

$F_{i}(\mathbf{k},\omega)=\begin{cases}\sqrt{P(\mathbf{k},\omega)},&\text{if
}|k-k_{i}|\leq\delta k\\\ 0,&\text{otherwise}\end{cases}.$ (17)

We define a characteristic function as

$G(\mathbf{k},\omega,\mathbf{U}_{\mathrm{eff}})=\mathrm{exp}\left[-2\left(\frac{\omega-\omega_{\mathrm{DR}}(\mathbf{k},\mathbf{U}_{\mathrm{eff}})}{a}\right)^{2}\right]+\mathrm{exp}\left[-2\left(\frac{\omega+\omega_{\mathrm{DR}}(-\mathbf{k},\mathbf{U}_{\mathrm{eff}})}{a}\right)^{2}\right],$
(18)

where $a$ is a parameter controlling the spectral width in frequency of the
dispersion shell, and we have included $\mathbf{U}_{\mathrm{eff}}$ as
independent variable affecting the dispersion relation for
$\omega_{\mathrm{DR}}$. The two terms on the right hand side reflect the
contributions from $\omega_{\pm}$ respectively. Previous implementations in
the literature have used a different form of (18) where $G$ is equal to unity
within a narrow frequency interval about the dispersion relation shell, and
equal to zero elsewhere (Serafino et al.,, 2009; Huang and Gill,, 2012; Huang
et al.,, 2016). The scalar product is defined as:

$\langle
A(\mathbf{k},\omega),B(\mathbf{k},\omega)\rangle\equiv\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}A(\mathbf{k},\omega)B(\mathbf{k},\omega)\mathrm{d}\mathbf{k}\mathrm{d}\omega.$
(19)

In practice, as $A,B$ herein correspond to discrete three-dimensional arrays,
the integrals in (19) are be evaluated as a summation of element-wise
multiplications.

The normalized scalar product between the spectrum and characteristic function
$G$ is then:

$V(\mathbf{U}_{\mathrm{eff}})=\frac{\langle
F_{i}(\mathbf{k},\omega),G(\mathbf{k},\omega,\mathbf{U}_{\mathrm{eff}})\rangle}{\sqrt{\langle
F_{i},F_{i}\rangle\langle G,G\rangle}}.$ (20)

Maximizing $V$ with respect to the Doppler shift velocity
$\mathbf{U}_{\mathrm{eff}}$ implies the greatest overlap between the
dispersion relation shell defined by the characteristic function and that of
the true spectrum.

The maximization is performed in two steps. The first involves evaluating $V$
over a wide search grid range covering the expected possible values of
$\mathbf{U}_{\mathrm{eff}}$. The value of the current with maximum $V$ over
this first search range is $\mathbf{U}_{\mathrm{eff}}^{\mathrm{initial}}$. The
second step further maximizes $V$ using
$\mathbf{U}_{\mathrm{eff}}^{\mathrm{initial}}$ as an initial guess to a
function maximizing routine, which finds a more precise value
$\mathbf{U}_{\mathrm{eff}}^{\mathrm{final}}$. In a MATLAB implementation, the
solver fminsearch is used to minimize $-V(\mathbf{U}_{\mathrm{eff}})$.

##  4 Results

###  4.1 Current estimates from simulated radar images and simulated waves

Figure 1: Root mean square error and bias of the effective current (left) and
direction (right) based on 50 Simulations. Simulated wave field and 6
different radar settings for the current $U(z)=\exp(0.5z)+0.05$. Figure 2:
Root mean square error and bias of the effective current and direction based
on 50 Simulations. Simulated wave field and 6 different radar settings for the
current$U(z)=\exp(0.2z)+0.05$.

The effective current is best estimated for the simulated wave field (compare
figures 1 and 2). The simulated radar images deviate with different magnitude
depending on the elevation of the radar antenna and the polarization. A height
of 45 m and horizontal polarization give the best results with a maximum RMS
deviation of about 0.05 m/s. With these radar conditions the estimate of the
current direction outperforms even the pure simulated waves up to a wavenumber
of 0.27 rad/m. The bias has a different sign for simulated waves and simulated
radar images. The effect of the bias for the current inversion will be
discussed below. It is worth noticing that the bias for the simulated waves
and the simulated radar images have opposite signs. The bias has an order of
magnitude similar to the RMS error. Hence, the imaging mechanism
systematically shifts the dispersion relation depending on the elevation of
the antenna height and the polarization. The results for the simulations for a
current with less shear in the upper layer ($U(z)=\exp(0.2z)+0.05$) were
compared.

####  4.1.1 Wave directional bandwidth and wave-current angle

Figure 3: The influence of the directional spreading for
$\mathrm{s}_{\max}=10,30,70$ (following Mitsuyasu et al., (1975) where large
spreading is associated with a low value). Statistics based on 10 simulated
wave fields for each directional spreading.

In an ideal setting for remote sensing of currents, the wave field would
consist of waves traveling in all directions. Realistic spectra are finite in
directional bandwidth however, and in this section we analyze the effect of
the directional bandwidth on the extracted Doppler shifts. The Doppler shifts
of the values $s_{\mathrm{max}}=10,30,70$ were compared for 10 cases of
simulated wave fields each (compare Fig. 3). As the directional bandwidth
increases, the Doppler shifts is improved especially in the lower wavenumbers.
It should be noted that the bias changes sign for the estimate of the
effective current for different cases.

Figure 4: Root mean square error and bias for different angles between
simulated waves and current (magnitude:left and angle: right).

The angle between the wave system and the current has a stronger effect on the
results than the configurations of the radar system in the low wavenumber
regime (compare Figure 4). For higher wavenumbers, the difference is still
noticable. The more current and waves are parallel, the lower the error of the
magnitude of the current. For the angle, it is opposite but the differences
are marginal.

####  4.1.2 Temporal and spatial window size

The finite dimensions of the sea surface image cube result in corresponding
spectral resolutions $\Delta\omega$ and $\Delta k$. The spectral resolution in
turn limits the accuracy with which the dispersion relation and Doppler shift
velocities can be determined. In this section we explore the effect of the
temporal and spatial window sizes on the accuracy of the Doppler shifts. In
figure 5 we show the Doppler shift velocities as a function of wavenumber for
different temporal durations $T$ and window lengths $L$ as indicated in the
respective legends. The results are not particularly sensitive to $T$ for
values greater than $300s$ (five minutes). The governing spectral resolution
is imposed by the spatial window which should not be reduced below five peak
wave lengths.

Figure 5: Dispersion retrieval for different temporal and spatial windows.

The sensitivity to the window size $L$ is strong, especially for larger
wavenumbers. The window size limits the minimum wavenumbers where Doppler
shift velocities may be extracted, thus having consequences for the range of
depths for which the current profile can be reconstructed.

To interpret the results of Fig. 5, we estimate the uncertainty scalings of
the dispersion relation due to the spectral resolution (Smeltzer et al.,,
2019):

$\Delta c_{\Delta k}=\frac{\mathrm{d}\omega_{0}}{\mathrm{d}k}\Delta k/k$ (21)

$\Delta c_{\Delta\omega}=\Delta\omega/k.$ (22)

Figure 6: Uncertainty scaling of the Doppler shift velocities based on
spectral resolution $\Delta k$ and $\Delta\omega$.

The parameters $\Delta c_{\Delta k}$ and $\Delta c_{\Delta\omega}$ may be
interpreted as the spectral resolution in units of velocity. We show in figure
6 $\Delta c_{\Delta k}$ and $\Delta c_{\Delta\omega}$ for $T=1200$ s and
$\Delta x=500$ m. Over all relevant wavenumbers, $\Delta c_{\Delta k}$ is
greater than $\Delta c_{\Delta\omega}$, meaning that uncertainty in the
dispersion relation is dominated by the wavenumber resolution as opposed to
frequency resolution. We now can interpret the results shown in Fig. 5: the
errors in the Doppler shift velocities are more sensitive to decreasing
$\Delta x$ due to the increase in $\Delta c_{\Delta k}$, whereas decreasing
$T$ up to point has little effect as $\Delta c_{\Delta k}$ still is dominant.
It is thus only when $T$ is decreased such that $\Delta
c_{\Delta\omega}\sim\Delta c_{\Delta k}$ that the uncertainty of the
dispersion relation becomes affected.

##  5 Profile reconstruction

This section evaluates how the error from the measurement of the dispersion
relation may affect the reconstruction of the current profile. Multiple
reconstruction methods have been suggested. A summary may be found in Smeltzer
and Ellingsen, (2021). Details on the inversion method employed herein can be
found in Smeltzer et al., (2019). It is not the purpose of this work to
provide the best current reconstructions but to illustrate how errors in the
measured dispersion relation change the result of the inversion.

In the course of this work it has been noticed that the shear profile with the
stronger shear studied above is not well suited for the reconstruction of the
current. Even when employing input of the correct dispersion relation. The
analysis of the effective current is therefore restricted to the case with
lower shear, $U(z)=\exp(0.2z)+0.05$.

The result for low wavenumbers is strongly affected by spectral leakage.

Figure 7: Current inversions for simulated waves ($\eta$) and simulated radar
images (polarization and antenna height as provided in the title) for 50
examples each. The true shear current is drawn as black dashed line. It is
$U(z)=\exp(0.5z)+0.05$ in the upper row and $U(z)=\exp(0.2z)+0.05$ in the
lower row.

Figure 8: Overall performance of current inversions for $U(z)=\exp(0.5z)+0.05$
(left) and $U(z)=\exp(0.2z)+0.05$ (right) as shown in Fig. 7.

The shear currents are retrieved by the PEDM algorithm (Smeltzer et al.,,
2019). Each realization is represented by single line in Fig. 7. Superimposed
lines result in a a darker colour in the plot. The true profiles are shown as
dashed black lines. Due to the strong bias that some of the dispersion
measurements have shown, a pre-processing step was needed. The effective
current is a smooth quantity and sudden jumps in the data can easily be
detected and eliminated. Further, the sampling variability leads to a non-
smooth curve. By fitting the remaining points to a polynomial of degree five
reduces the unwanted noisiness in the measured dispersion relation and
improves the shear current reconstruction.

For each case 50 realizations were analyzed. The shear is underestimated in
the uppermost layer. In deeper layers the estimated shear current has an
offset where the direction of the offset depends on the profile. For very
strong shear the method tends to overshoot in lower layers. Overall, the
results for the simulated waves and the different radar configurations lead to
similar results. The differences become apparent by evaluating the root mean
square (RMS) error as shown in Fig. 8. The first bar shows the error when
using the correct effective current as input to the method. It shows that the
method is not able to obtain the correct profile even with the correct input,
the error also depends strongly on the shape of the current profile. The
method is better suited for cases of milder shear as the RMS error for exact
input can be reduced by one half with respect to a stronger current. For the
Doppler shifts extracted from simulated waves the error is in the range of the
exact input or above. Adding the imaging mechanism leads to a wider spreading
of solutions for the different realizations. For the two example profiles,
this effect is stronger in the profile with the lesser shear. The increase in
spread is more pronounced for horizontal polarization. The 95% confidence
interval ranges up to 3 cm/s up to a water depth of 10 m. The size of the
confidence intervals varies up to 1 cm/s between the cases. Overall, the
antenna height of $H=45$ m is favourable and vertical polarization should be
used.

##  6 Discussion and conclusion

The presented work confirms that X-band radars are capable of measuring a part
of the dispersion relation corresponding to wavenumbers in the range from 0.1
to 0.35 rad/m. By the offset from the theoretical dispersion relation the
underlying effective current can be estimated. Depending on the current
profile and the wavenumber the RMS lies below 0.1 m/s and below 0.04 m/s for
$k>0.15$rad/m. The reconstruction of the shear current is less accurate.
Although the average performance over 50 realizations shows reasonable
results, there is a considerable spread of single realizations. It has been
shown by simulations that the imaging mechanisms amplifies both the deviations
from the true profiles, and the spread for different realizations. When
comparing the performance of the inversion method by inserting the exact
effective current, the error depends strongly on the input profile. For the
studied examples the error was doubled.

It should be noted that the study implements the radar measurement as a
Eulerian measurement rather than a Lagrangian one (compare (Pizzo et al.,,
2023)). Therefore the Stokes drift that is normally incorporated in current
retrievals from radar images is not contained in the simulations. To our
understanding it does not affect the findings herein as the Stokes drift is
neither present in the wave model nor in the radar image.

Further improvements of the complete scheme may be expected by resolving the
dispersion by instruments of higher resolution. As the shear has the strongest
effect in the uppermost ocean layer more information in the short wave segment
would improve the results. The observed decrease of the error for increasing
wavenumbers results from a more evenly distribution of energy in the spectral
domain for higher frequencies. For one, the number of waves in relation to the
domain is larger and shorter waves have a higher directional spreading.

As outlined above, the spatial resolution limits the range of suitable
wavenumbers for the current estimates at the lower end. The influence of the
window size may proof challenging in reality. To comply with limitations on
the homogeneity of current fields, the patch size was limited to around 500 m
in each direction. The wave length was deliberately chosen to not fit into the
window. The algorithm of the current retrieval is subject to spectral leakage
that causes larger errors. It may help to use wider patches for the lower
wavenumbers, assuming that the current in deeper layers is less inhomogeneous.
This is however only valid in deep water.

Figure 9: Analysis of influence of random noise on the retrieval of the
current. Left: RMS error based on 1000 realizations. Right: Error distribution
at three meters below the mean water level.

Further improvement of the inversion algorithm with a stronger focus on
possible noise and biases in the input would be advisable. In the presented
study a crucial step was to sort out Doppler estimates in the low wavenumber
regime to obtain good shear current reconstructions. Hence, the method is very
sensitive to errors. In its current form the simulation of radar images is
free of noise and errors are probably underestimated. An additional analysis
of the influence of random noise in the current inversions was conducted. 1000
realizations of effective currents ($U(z)=\exp(0.5z)+0.05,\psi=30^{\circ}$
with random noise were analyzed. It confirms that the errors increase with
depth but there is no bias introduced (compare Fig. 9).

It would be desirable to extend the presented work to intermediate and shallow
waters. However, simulating realistic bathymetries with waves and currents is
not straightforward. The inhomogenity of condtions cause a smearing of the
dispersion cone.

For future campaigns on retrieval of the current based on radar images, a high
antenna (e.g. 45 m) and vertical polarization are preferable (for a distance
of 200 meters between the radar and the observed waves). If the evolution of
the current is of interest 5-10 minute time windows are suggested as trade-off
between temporal resolution and low error. The patch size of the observation
should include more than five peak wavelengths if homogeneity of the current
can be assumed.

## Acknowledgements

This work was funded by Minerva Fund and the European Research Council through
the HIGHWAVE project (grant no. 833125).

## References

  * Campana et al., (2016) Campana, J., Terrill, E. J., and De Paolo, T. (2016).  The development of an inversion technique to extract vertical current profiles from x-band radar observations.  Journal of Atmospheric and Oceanic Technology, 33(9):2015–2028.
  * Davis et al., (1981) Davis, R., deSzoeke, R., and Niiler, P. (1981).  Variability in the upper ocean during mile. part ii: Modeling the mixed layer response.  Deep Sea Research Part A. Oceanographic Research Papers, 28(12):1453–1475.
  * Ellingsen and Li, (2017) Ellingsen, S. Å. and Li, Y. (2017).  Approximate dispersion relations for waves on arbitrary shear flows.  Journal of Geophysical Research: Oceans, 122(12):9889–9905.
  * Huang et al., (2016) Huang, W., Carrasco, R., Shen, C., Gill, E. W., and Horstmann, J. (2016).  Surface current measurements using x-band marine radar with vertical polarization.  IEEE Transactions on Geoscience and Remote Sensing, 54(5):2988–2997.
  * Huang and Gill, (2012) Huang, W. and Gill, E. (2012).  Surface current measurement under low sea state using dual polarized x-band nautical radar.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(6):1868–1873.
  * Ivonin et al., (2004) Ivonin, D. V., Broche, P., Devenon, J.-L., and Shrira, V. I. (2004).  Validation of hf radar probing of the vertical shear of surface currents by acoustic doppler current profiler measurements.  Journal of Geophysical Research: Oceans, 109(C4).
  * Kirby and Chen, (1989) Kirby, J. T. and Chen, T.-M. (1989).  Surface waves on vertically sheared flows: approximate dispersion relations.  Journal of Geophysical Research: Oceans, 94(C1):1013–1027.
  * Lund et al., (2020) Lund, B., Haus, B., Graber, H., Horstmann, J., Carrasco, R., Novelli, G., Guigand, C., Mehta, S., Laxague, N., and Özgökmen, T. (2020).  Marine x-band radar currents and bathymetry: An argument for a wave number-dependent retrieval method.  Journal of Geophysical Research: Oceans, 125(2):e2019JC015618.
  * Lund et al., (2018) Lund, B., Haus, B. K., Horstmann, J., Graber, H. C., Carrasco, R., Laxague, N. J., Novelli, G., Guigand, C. M., and Özgökmen, T. M. (2018).  Near-surface current mapping by shipboard marine x-band radar: A validation.  Journal of Atmospheric and Oceanic Technology, 35(5):1077–1090.
  * Mitsuyasu et al., (1975) Mitsuyasu, H., Tasai, F., Suhara, T., Mizuno, S., Ohkusu, M., Honda, T., and Rikiishi, K. (1975).  Observations of the directional spectrum of ocean wavesusing a cloverleaf buoy.  Journal of Physical Oceanography, 5(4):750 – 760.
  * Nieto Borge et al., (2004) Nieto Borge, J. C., Rodriguez, G. R., Hessner, K., and Gonzalez, P. I. (2004).  Inversion of marine radar images for surface wave analysis.  Journal of Atmospheric and Oceanic Technology, 21.
  * Pizzo et al., (2023) Pizzo, N., Lenain, L., Rømcke, O., Ellingsen, S. Å., and Smeltzer, B. K. (2023).  The role of lagrangian drift in the geometry, kinematics and dynamics of surface waves.  Journal of Fluid Mechanics, 954:R4.
  * Rees, (2013) Rees, W. G. (2013).  Physical Principles of Remote Sensing:.  Cambridge University Press, Cambridge, 2 edition.
  * Serafino et al., (2009) Serafino, F., Lugni, C., and Soldovieri, F. (2009).  A novel strategy for the surface current determination from marine x-band radar data.  IEEE Geoscience and Remote Sensing Letters, 7(2):231–235.
  * Smeltzer et al., (2019) Smeltzer, B. K., Æsøy, E., Ådnøy, A., and Ellingsen, S. Å. (2019).  An improved method for determining near-surface currents from wave dispersion measurements.  Journal of Geophysical Research: Oceans, 124(12):8832–8851.
  * Smeltzer and Ellingsen, (2021) Smeltzer, B. K. and Ellingsen, S. Å. (2021).  Current mapping from the wave spectrum.  In Ocean Remote Sensing Technologies: High frequency, marine and GNSS-based radar, Radar, Sonar and Navigation, pages 357–378. Institution of Engineering and Technology.
  * Stewart and Joy, (1974) Stewart, R. H. and Joy, J. W. (1974).  Hf radio measurements of surface currents.  In Deep sea research and oceanographic abstracts, volume 21, pages 1039–1049. Elsevier.
  * Young et al., (1985) Young, I. R., Rosenthal, W., and Ziemer, F. (1985).  A three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents.  JGR, 90:1049–1059.