# Attainable set for rank $3$ step $2$ free Carnot group
with positive controls 111The work is supported by the Russian Science
Foundation under grant 22-21-00877 (https://rscf.ru/en/project/22-21-00877/)
and performed in Ailamazyan Program Systems Institute of Russian Academy of
Sciences.

A. V. Podobryaev
A. K. Ailamazyan Program Systems Institute of RAS
<EMAIL_ADDRESS>
###### Abstract

We find the attainable set for a control system on the free Carnot group of
rank $3$ and step $2$ with positive controls. This kind of control systems is
connected with the theory of free Lie semigroups; with some estimates for
probabilities of inequalities for independent random variables; with the
nilpotent approximation of robotic control systems and with contour recovering
without cusps in image processing. We investigate the boundary of the
attainable set with the help of the Pontryagin maximum principle for the time-
optimal control problem. We study extremal trajectories that correspond to
bang-bang, singular and mixed controls. We obtain upper bounds for the number
of switchings for optimal controls. This implies a parametrization of the
boundary faces of the attainable set.

Keywords: attainable set, Carnot group, time-optimal control problem, bang-
bang control, nilpotent Lie semigroup

AMS subject classification: 93B03, 49K15, 35R03, 22E25, 20M05.

##  1 Introduction

Consider the vector space $\mathbb{R}^{r}\times\mathbb{R}^{\frac{r(r-1)}{2}}$
of pairs of a vector $x$ and a skew-symmetric matrix $y$. This vector space
equipped with the multiplication rule

$(x,y)\cdot(x^{\prime},y^{\prime})=(x+x^{\prime},y+y^{\prime}+(xx^{\prime
T}-x^{\prime}x^{T}))$ (1)

is called _a free Carnot group $G$ of rank $r$ and step $2$_. The following
control system

$\begin{array}[]{ll}\dot{x}_{i}=u_{i},&i=1,\dots,r,\\\
\dot{y}_{ij}=x_{i}u_{j}-x_{j}u_{i},&i<j\\\ \end{array}$ (2)

is invariant under the left action of the group $G$. Control systems on Carnot
groups are some kind of cornerstones in geometric control theory [1] due to
existence of a nilpotent approximation for general control systems [2]. In
particular, such control systems appear in several robotic systems [3] and in
some models for contour reconstruction without cusps in image processing [4].
We consider system (2) with controls $u_{1},\dots,u_{r}\geqslant 0$. Our goal
is a description of the attainable set $\mathcal{A}$ from the point $(0,0)$
for this control system.

For any set of positive numbers $(c_{1},\dots,c_{r})$ consider _a dilation_
that is a map $D_{c_{1},\dots,c_{r}}:G\rightarrow G$ such that

$D_{c_{1},\dots,c_{r}}(x_{i},y_{ij})=(c_{i}x_{i},c_{i}c_{j}y_{ij}),\qquad
i=1,\dots,r,\qquad i<j.$ (3)

Note that dilations preserve system (2).

Denote by $\mathcal{A}_{1}=\\{(x,y)\in\mathcal{A}\,|\,x=(1,\dots,1)\\}$ the
section of the attainable set. Since $\mathcal{A}$ coincides with the closure
of the orbit of $\mathcal{A}_{1}$ with respect to dilations (3), then it is
sufficient to describe the section $\mathcal{A}_{1}$.

H. Abels and È. B. Vinberg [5] suggested a probability interpretation of the
section $\mathcal{A}_{1}$. It turns out that the coordinates of points of this
set are $y_{ij}=2p_{ij}-1$ for $i<j$, where $p_{ij}=P(\xi_{i}<\xi_{j})$ for
independent random variables $\xi_{1},\dots,\xi_{r}$ such that
$P(\xi_{1}=\dots=\xi_{r})=0$. So, the set
$\mathcal{B}={{1}\over{2}}((1,\dots,1)^{T}+\mathcal{A}_{1})$ has a probability
interpretation.

For $e^{t_{1}X_{i_{1}}}e^{t_{2}X_{i_{2}}}\dots
e^{t_{n}X_{i_{n}}}\in\mathcal{A}_{1}$ from (1) we obtain

$\sum\limits_{l\,|\,i_{l}=i}{t_{i}}=1,\qquad
p_{ij}=\sum\limits_{l<m\,|\,i_{l}=i,i_{m}=j}{t_{l}t_{m}}.$ (4)

In paper [5] the attainable set $\mathcal{A}$ was regarded as _a two-step
nilpotent Lie semigroup of rank $r$_ and the set $\mathcal{B}$ (that is
equivalent to the section $\mathcal{A}_{1}$) was described for $r=3$ with the
help of an algebraic method. The result implies some not obvious bounds for
probabilities $P(\xi_{i}<\xi_{j})$. But it seems that the method is
unapplicable in more complex cases $r>3$.

In this paper we propose a constructive and algorithmic approach for $r=3$.
The conjecture is that this method will be useful also for $r>3$.

The main idea is to consider the time-optimal problem on the group $G$ and
then to describe extremal trajectories that come to the boundary of the
attainable set. They are optimal trajectories for a reduced problem. The
Pontryagin maximum principle implies necessary conditions of optimality. This
allows us to find the extremal trajectories. Then second order optimality
conditions distinguish optimal trajectories. More precisely, we get the number
of control switchings for different kinds of optimal trajectories (bang-bang,
singular and mixed). This implies a parametrization of the boundary faces of
the section $\mathcal{A}_{1}$ of the attainable set. Note that endpoints of
bang-bang trajectories with periodical control sweep non-trivial faces of the
boundary of the attainable set. We give the answer in terms of the set
$\mathcal{B}$ and the coordinates $p=p_{12}$, $q=p_{23}$, $r=p_{31}$ on this
set, see (4).

The paper has the following structure. In Section 2 we introduce a time-
optimal problem and necessary conditions for extremal trajectories coming to
the boundary of the attainable set. Then we describe types of extremal
trajectories in Section 3. We find an upper bound for the number of control
switchings for bang-bang trajectories in Section 4. We study singular and
mixed trajectories in Sections 5 and 6, respectively. Finally, we summarize
obtained results in the main theorem in Section 7.

The author would like to thank prof. A. M. Tsirlin for useful discussions and
E. A. Podobryaeva for the figures.

##  2 Time-optimal problem

In this section we formulate a time-optimal control problem, apply the
Pontryagin maximum principle and derive the adjoint subsystem of ODEs that
describes extremal trajectories. Here we consider the general case of
arbitrary $r\geqslant 2$.

Consider a time-optimal problem

$\begin{array}[]{ll}(x,y)\in\operatorname{Lip}{([0,T],G)},&\\\
x(0)=0,&x(T)=x^{1},\\\ y(0)=0,&y(T)=y^{1},\\\
\dot{x}_{i}=u_{i},&i=1,\dots,r,\\\ \dot{y}_{ij}=x_{i}u_{j}-x_{j}u_{i},&i<j,\\\
(u_{1},\dots,u_{r})\in
L^{\infty}([0,T],U),&T\rightarrow\mathrm{min},\end{array}$ (5)

where the set of controls is a simplex

$U=\\{(u_{1},\dots,u_{r})\,|\,u_{1}+\dots+u_{r}=1,\,u_{1},\dots,u_{r}\geqslant
0\\}.$

######  Proposition 1.

Every admissible trajectory of control system $\eqref{eq-time-optimal}$ is
optimal.

###### Proof.

Using an argument from [6], we obtain that the distance from a point of an
admissible trajectory to the plane $x_{1}+\dots+x_{r}=1$ is proportional to
the time. ∎

######  Proposition 2.

An optimal trajectory of problem $\eqref{eq-time-optimal}$ comes to the
boundary of the attainable set iff it is time-optimal for the problem
$\eqref{eq-time-optimal}$ factorized by the subspace
$\mathbb{R}(1,\dots,1)^{T}$ (so called reduced problem).

###### Proof.

The original problem is an extended problem of the reduced one, see for
example [1]. So, optimal trajectories of the reduced problem lift to the
trajectories coming to the boundary of the attainable set of the original one.
∎

For $i=1,\dots,r$ let $X_{i}$ be a basis vector field corresponding to
coordinate $x_{i}$. Denote by $h_{i}=\langle X_{i},\,\cdot\,\rangle$ linear on
the fibers of the cotangent bundle $T^{*}G$ functions. Introduce a family of
functions on $T^{*}G$ depending on a parameter $u=(u_{1},\dots,u_{r})^{T}$:

$H_{u}(\lambda)=u_{1}h_{1}(\lambda)+\dots+u_{r}h_{r}(\lambda),\qquad\lambda\in
T^{*}G.$ (6)

Necessary conditions of optimality follow from the Pontryagin maximum
principle [7, 1]. These conditions are the same for extremal trajectories of
problem (5) coming to the boundary of the set $\mathcal{A}$ and for optimal
trajectories of the reduced control system.

######  Theorem 1.

If $(\tilde{x},\tilde{y})$ and $\tilde{u}$ are an optimal process for problem
$\eqref{eq-time-optimal}$, then there exist a curve
$\lambda\in\operatorname{Lip}{([0,T],T^{*}G)}$,
$\pi(\lambda(t))=(\tilde{x}(t),\tilde{y}(t))$ and a number $\nu\in\\{0,1\\}$
such that for a.e. $t\in[0,T]$ we have

$\begin{array}[]{ll}\dot{\lambda}(t)=\vec{H}_{\tilde{u}(t)}(\lambda(t)),&\lambda(t)\neq
0,\\\ H_{\tilde{u}(t)}(\lambda(t))=\max\limits_{u\in
U}{H_{u}(\lambda(t))},&H_{\tilde{u}(t)}(\lambda(t))\equiv\nu,\\\ \end{array}$

where $\pi:T^{*}G\rightarrow G$ is the natural projection and $\vec{H}$ is the
Hamiltonian vector field corresponding to a Hamiltonian $H$.

Recall that the curve $\lambda$ is called _an extremal_ and its projection
$\pi(\lambda)$ is called _an extremal trajectory_.

Let $Y_{ij}=[X_{i},X_{j}]$ and $h_{ij}=\langle Y_{ij},\,\cdot\,\rangle$. Then
the adjoint subsystem of the Hamiltonian system of the Pontryagin maximum
principle reads as

$\dot{h}(t)=R\tilde{u}(t),\qquad\dot{R}=0,$ (7)

where $h=(h_{1},\dots,h_{r})^{T}$ and $R=(h_{ij})$ is a skew-symmetric matrix.

It is easy to see that if $\nu=0$, then $u=0$ and the corresponding extremal
trajectory is constant $(0,0)$. So, we assume that $\nu=1$. It follows from
the maximum condition that the adjoint trajectories $h(\,\cdot\,)$ for
extremals lie on the boundary of the quadrant

$Q=\\{(h_{1},\dots,h_{r})^{T}\in(\mathbb{R}^{r})^{*}\,|\,h_{1},\dots,h_{r}\leqslant
1\\}$

that is a dual set to the set of controls $U$.

##  3 Types of extremals

Here and below we fix $r=3$. Introduce some notation. For $i=1,2,3$ let
$F_{i}$ be the face of the quadrant $Q$ defined by the condition $h_{i}=1$.
Next, for $i,j=1,2,3$ denote by $E_{ij}=F_{i}\cap F_{j}$ the edges of $Q$.
Finally, let $E=E_{12}\cup E_{23}\cup E_{31}$ be the union of the edges and
$V=(1,1,1)^{T}$ be the vertex.

Let $h:[0,T]\rightarrow(\mathbb{R}^{3})^{*}$ be a trajectory of the adjoint
subsystem. Note that if $h(t)$ is inside some face of the quadrant $Q$, then
the corresponding control is at the vertex of the triangle $U$. If $h(t)$ is
on some edge of $Q$, then the corresponding control lies on the edge of the
triangle $U$. If $h(t)=V$, then the control may be anywhere at the triangle
$U$.

######  Proposition 3.

For problem $\eqref{eq-time-optimal}$ there are three types of extremals
$\lambda:[0,T]\rightarrow T^{*}G$:
$(1)$ bang-bang, i.e., such that

$\operatorname{card}{\\{t\in[0,T]\,|\,h(t)\in E\\}}<\infty;$

$(2)$ singular, i.e., $h(t)\in E_{i}$ for some $i=1,2,3$ or $h(t)\in V$;
$(3)$ mixed, i.e., a concatenation of a finite number of singular and bang-
bang arcs, where $h:[0,T]\rightarrow(\mathbb{R}^{3})^{*}$ is the corresponding
trajectory of the adjoint subsystem.

###### Proof.

The function

$C=h_{1}h_{23}+h_{2}h_{31}+h_{3}h_{12}$ (8)

is a Casimir function. In particular, it is an integral of system (7). Hence,
the trajectories $h(t)$ lie in the intersections of the boundary of the
quadrant $Q$ with two-dimensional affine subspaces $C=\mathrm{const}$.
Possible situations are presented in Fig. 1 for bang-bang trajectories and the
singular trajectory at the vertex $V$ and Fig. 2 where non-trivial singular
and mixed arcs appear. ∎

(a) $1$ or $2$ switchings

(b) periodical control

Figure 1: Adjoint subsystems for bang-bang trajectories and one singular
trajectory.

(a) concatenation of singular and bang arcs

(b) concatenation of two singular arcs

Figure 2: Adjoint subsystems with singular trajectories.

##  4 Bang-bang trajectories

In this section we get an upper bound for the number of switchings on optimal
bang-bang curves using the following theorem by A. A. Agrachev and R. V.
Gamkrelidze [8]. Earlier this method was successfully applied for analysis of
bang-bang extremals on the Cartan group [9].

######  Theorem 2.

Assume that $q(\,\cdot\,)$ is an extremal curve, $u(\,\cdot\,)$ is an extremal
control of problem $\eqref{eq-time-optimal}$ and $\lambda(\,\cdot\,)$ is the
corresponding extremal that is unique upto a multiplication by a positive
scalar. Let $u(t)=u^{i}$ for $t\in(t_{i},t_{i+1})$, where
$0=t_{0}<t_{1}<\dots<t_{k+1}=T$, and let $V_{i}$ be the velocity corresponding
to the control $u^{i}$. Define recursively the following operators:

$P_{0}=P_{1}=\mathrm{id},\qquad P_{i}=P_{i-1}\circ
e^{(t_{i}-t_{i-1})\operatorname{ad}{V_{i-1}}},\qquad i=2,\dots,k.$

Define the vector fields $Z_{i}=P_{i}V_{i}$. If the quadratic form

$G(\alpha)=\sum\limits_{0\leqslant i<j\leqslant
k}{\alpha_{i}\alpha_{j}\langle\lambda(t_{1}),[Z_{i},Z_{j}](q(t_{1}))\rangle}$

is not negative semidefinite on the space

$W=\\{(\alpha_{0},\dots,\alpha_{k})\in\mathbb{R}^{k+1}\bigm{|}\sum\limits_{i=0}^{k}{\alpha_{i}}=0,\,\sum\limits_{i=0}^{k}{\alpha_{i}Z_{i}(q(t_{1}))}=0\\},$

then the trajectory $q(\,\cdot\,)$ is not optimal.

######  Proposition 4.

If an extremal curve has more than 4 switchings, then it is not optimal.

###### Proof.

Let us apply Theorem 2 to an extremal curve with $k=5$ switchings. For
definiteness let

$\begin{array}[]{lll}V_{0}=X_{1},&V_{1}=X_{2},&V_{2}=X_{3},\\\
V_{3}=X_{1},&V_{4}=X_{2},&V_{5}=X_{3}\\\ \end{array}$

and $\tau_{i}=t_{i}-t_{i-1}$ for $i=1,\dots,5$. We have

$\begin{array}[]{c}Z_{0}=X_{1},\qquad Z_{1}=X_{2},\qquad
Z_{2}=X_{3}+\tau_{2}Y_{23},\\\
Z_{3}=X_{1}-\tau_{2}Y_{12}-\tau_{3}Y_{13},\qquad
Z_{4}=X_{2}-\tau_{3}Y_{23}+\tau_{4}Y_{12},\\\
Z_{5}=X_{3}+(\tau_{2}+\tau_{5})Y_{23}+\tau_{4}Y_{13}.\\\ \end{array}$

Note that $\tau_{5}=\tau_{2}$, since this is the time of passing the same edge
of a triangle by the extremal trajectory, see Fig. 1b. Hence, from the system
of equations for the subspace $W$ we get

$\begin{array}[]{lll}\alpha_{0}=-\frac{\tau_{4}}{\tau_{3}}\alpha_{5},&\alpha_{1}=-\frac{\tau_{2}}{\tau_{3}}\alpha_{5},&\alpha_{2}=-\alpha_{5},\\\
\alpha_{3}=\frac{\tau_{4}}{\tau_{3}}\alpha_{5},&\alpha_{4}=\frac{\tau_{2}}{\tau_{3}}\alpha_{5},&\alpha_{5}\in\mathbb{R}.\\\
\end{array}$

Next, for quadratic form $G$ we obtain

$G(\alpha)=\frac{2\tau_{2}\tau_{4}}{\tau_{3}}\left(\frac{h_{12}}{\tau_{3}}+\frac{h_{23}}{\tau_{4}}-\frac{h_{31}}{\tau_{2}}\right)\alpha_{5}^{2}.$

Note that since we have a triangle as a section of $Q$ and then the
coefficients of $C$ are positive $h_{12},h_{23},h_{31}>0$. Using (8) it is not
difficult to see that

$\tau_{2}=\frac{K}{h_{12}h_{23}},\qquad\tau_{3}=\frac{K}{h_{23}h_{31}},\qquad\tau_{4}=\frac{K}{h_{31}h_{12}},$

where $K=h_{12}+h_{23}+h_{31}-C>0$. Therefore, the sign of $G(\alpha)$
coincide with the sign of
$h_{12}h_{23}h_{31}+h_{23}h_{31}h_{12}-h_{31}h_{12}h_{23}=2h_{12}h_{23}h_{31}>0$.
Thus, the quadratic form $G$ is positive definite on the subspace $W$. So, the
extremal trajectory under consideration can not be optimal. ∎

Now we describe a part of the boundary $\partial\mathcal{A}_{1}$ (or
equivalently $\partial\mathcal{B}$) corresponding to the endpoints of bang-
bang extremal curves with the number of switchings less than $4$.

######  Proposition 5.

Bang-bang trajectories with not more than $3$ switchings correspond to
$\partial\mathcal{B}$.

###### Proof.

Obviously $\mathcal{B}\subset[0,1]^{3}$. Using (4) we conclude, that
trajectories without switchings or with one switching correspond to the
vertices $A_{1}=(1,0,0)$, $B_{2}=(0,1,0)$, $C_{1}=(0,0,1)$ of the cube
$[0,1]^{3}$, see Fig. 3. Next, trajectories with two switchings match to the
vertices $A_{2}=(1,0,1)$, $C_{2}=(0,1,1)$, $D_{1}=(1,1,0)$. Finally,
trajectories with three switchings respond to the edges of the triangles
$A_{1}B_{2}C_{1}$ and $A_{2}C_{2}D_{1}$. These edges lie on the boundary of
the cube. This implies that these edges are in $\partial\mathcal{B}$. ∎

Figure 3: The cube in the $(p,q,r)$-space containing the set $\mathcal{B}$.

######  Proposition 6.

Extremal trajectories with $4$ switchings correspond to quadratic surfaces in
$\mathcal{B}$ of the type

$p+qr=1\qquad\text{or}\qquad(1-p)+(1-q)(1-r)=1$

for any cyclic permutations of the variables

$p=P(\xi_{1}<\xi_{2}),\qquad q=P(\xi_{2}<\xi_{3}),\qquad
r=P(\xi_{3}<\xi_{1}).$

###### Proof.

This follows directly from (4). The sequence of controls defines the type of
the equation. Even permutations of controls correspond to the first type of
equations and odd permutations correspond to the second one. ∎

##  5 Singular trajectories

The next proposition allows us to consider singular trajectories with a
control on the boundary $\partial U$.

######  Proposition 7.

If a control for an extremal trajectory lies in the interior of the control
set $U$, then this trajectory lies in the interior of the attainable set
$\mathcal{A}$.

###### Proof.

Assume by contradiction that this trajectory comes to the boundary
$\partial\mathcal{A}$. It follows that this trajectory is optimal for the
reduced problem. Project the corresponding control from the center of the set
$U$ to its boundary and use this projection as a new control. Thus, this new
control allows to go faster along the same extremal trajectory of the reduced
problem. ∎

Therefore, the study of singular trajectories corresponding to the vertex $V$
is reduced to the study of bang-bang trajectories, singular trajectories
corresponding to an edge of the quadrant $Q$ and mixed trajectories.

Let a singular trajectory corresponds to the edge $E_{ij}$. Consider a free
two-step Carnot subgroup $H_{ij}$ such that its Lie algebra is generated by
the elements $X_{i}$ and $X_{j}$. (The subgroup $H_{ij}$ is isomorphic to _the
Heisenberg group_.) Denote the attainable set for the subgroup $H_{ij}$ by
$\mathcal{A}(H_{ij})$.

######  Proposition 8.

It is sufficient to consider singular curves corresponding to the edge
$E_{ij}$ with two switchings to attain the set $\partial\mathcal{A}(H_{ij})$.
Three switchings are sufficient to attain any point of the set
$\mathcal{A}(H_{ij})$.

###### Proof.

Since the factor of the set $\mathcal{A}(H_{ij})$ by dilations is just an
interval $\mathcal{A}_{1}(H_{ij})=[-1,1]$, then this statement could be
checked directly. ∎

##  6 Mixed trajectories

######  Proposition 9.

The endpoints of mixed trajectories that are concatenation of a singular arc
corresponding to an edge of the quadrant $Q$ and a bang arc _(_ see Fig.
_2a)_ match to the edges of the set $\mathcal{B}$.

###### Proof.

It can be checked directly using (4) that the ends of these trajectories
correspond to the edges $A_{1}A_{2}$, $A_{2}B_{2}$, $B_{2}C_{2}$,
$C_{2}C_{1}$, $C_{1}D_{1}$, $D_{1}A_{1}$ of the cube. Thus, these endpoints
correspond to $\partial\mathcal{B}$. ∎

######  Proposition 10.

The endpoints of mixed trajectories that are concatenation of two singular
arcs corresponding to two edges of the quadrant $Q$ _(_ see Fig.  _2b)_ match
to the triangular faces $A_{1}A_{2}B_{2}$, $B_{2}C_{2}C_{1}$ and
$C_{1}D_{1}A_{1}$ of the set $\mathcal{B}$.

###### Proof.

Consider two edges of the quadrant $Q$, say $E_{12}$ and $E_{23}$. The
attainable set for such mixed trajectories is
$\mathcal{A}(H_{12})\cdot\mathcal{A}(H_{23})=\partial\mathcal{A}(H_{12})\cdot\mathcal{A}(H_{23})$.
The last equality follows from Proposition 8. So, it is sufficient to make one
switching on the first singular arc and two switchings on the second one. The
ends of these trajectories correspond to the triangles $A_{1}A_{2}B_{2}$,
$B_{2}C_{2}C_{1}$ and $C_{1}D_{1}A_{1}$. Since these triangles lie on the
boundary of the cube, then they are parts of the boundary
$\partial\mathcal{B}$. ∎

##  7 Conclusion

We considered all types of extremal trajectories for the time-optimal problem
(5) and now we can summarize the result.

######  Theorem 3.

$(1)$ The attainable set $\mathcal{A}$ for a free two-step Carnot group of
rank $3$ can be obtained applying all dilations $\eqref{eq-dilation}$ to the
section $\mathcal{A}_{1}=2\mathcal{B}-(1,\dots,1)^{T}$ and taking the closure.
$(2)$ The set $\mathcal{B}$ is a curved polyhedron carved from the cube
$[0,1]^{3}$ with $(p,q,r)$-coordinates by the surfaces

$p+qr=1,\qquad(1-p)+(1-q)(1-r)=1$

and surfaces obtained from these surfaces by any permutations of the variables
$p,q,r$.
$(3)$ The vertices, edges and faces of $\mathcal{B}$ correspond to the
endpoints of the extremal trajectories of problem $\eqref{eq-time-optimal}$
with not more than $2$, $3$ and $4$ switchings, respectively.

## References

  * [1] A. A. Agrachev and Yu. L. Sachkov, Control Theory from the Geometric Viewpoint. Encyclopaedia of Mathematical Sciences, vol. 87. Springer-Verlag, 2004.
  * [2] A. A. Agrachev and A. V. Sarychev, “Filtrations of a Lie algebra of vector fields and the nilpotent approximation of controllable systems,” Dokl. Akad. Nauk SSSR, vol. 295, no. 4, pp. 777–781, 1987.
  * [3] A. A. Ardentov and A. P. Mashtakov, “Control of a mobile robot with a trailer based on nilpotent approximation,” Automation and Remote Control, vol. 82, pp. 73–92, 2021.
  * [4] A. P. Mashtakov, “Sub-Riemannian geometry in image processing and modeling of the human visual system,” Rus. J. Nonlin. Dyn., vol. 15, no. 4, pp. 561–568, 2019.
  * [5] H. Abels and È. B. Vinberg, “On free two-step nilpotent Lie semigroups and inequalities between random variables,” J. Lie Theory, vol. 29, no. 1, pp. 79–87, 2019.
  * [6] A. A. Ardentov, E. Le Donne and Yu. L. Sachkov, “A sub-Finsler problem on the Cartan group,” Proceedings of the Steklov Institute of Mathematics, vol. 304, pp. 42–59, 2019.
  * [7] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchenko, The Mathematical Theory of Optimal Processes. Pergamon Press, Oxford, 1964.
  * [8] A. A. Agrachev and R. V. Gamkrelidze, “Symplectic geometry for optimal control,” Nonlinear controllability and optimal control, Monogr. Textbooks Pure Appl. Math., vol. 133, pp. 263–277, Dekker, New York, 1990.
  * [9] A. A. Ardentov, E. Le Donne and Yu. L. Sachkov, “Sub-Finsler geodesics on the Cartan group,” Regular and Chaotic Dynamics, vol. 24, no. 1, pp. 36–60, 2019.


# Casimir force between two substrates within three different layers using the
scattering approach

A. Moradian  Department of Science, Campus of Bijar, University of Kurdistan,
Bijar, Iran  M. R. Setare<EMAIL_ADDRESS>Department of Science, University
of Kurdistan, Sanandaj, Iran  S. A. Seyedzahedi  Department of Science,
University of Kurdistan, Sanandaj, Iran

###### Abstract

We compute the Casimir force for a system composed of two layers as substrates
within three different homogenous layers. We use the scattering approach along
with the Matsubara formalism in order to calculate the Casimir force at finite
temperature. We impose the appropriate boundary condition on the tangential
components of the electric and magnetic fields to construct the reflection
matrices. Our calculation is simple such that one can extend it to systems
containing inhomogeneous layers as good candidates for designing nanomachines.
We find that the Casimir force is proportional to $L^{-5}$ where $L$ is the
thickness of the inner layer.

###### pacs:

03.65.Nk, 42.50.Lc

##  I Introduction

The Casimir effect l1  resulting from modifying of the vacuum fluctuations due
to the insertion of the boundaries has entered a new era of novel accurate
measurements. According to the rapid progress of the nanotechnology and the
introduction of the Casimir force offering new possibility for designing
nanomechanical systems l5 ; l5-2 ; l5-3 ; l5-4 , it will be useful to utilize
multilayers in the investigating systems to evaluate Casimir energy. Having a
repeating arrangement of thin layers of two different materials, this class of
materials plays an important role in different branches of science and
technology. Studying the properties of multilayers continues to grow
exponentially due to their application at the biological interfaces for
controlling drug molecules’ release and permission as an instance l5-1 . The
non-pairwise addition effect in the van der Waals-Lifshitz interaction has
been investigated in a multilayer system in l6 ; l7 . They have explored the
effect of the presence of other layers in the multilamellar geometries
applying an algebra of $2\times 2$ matrices resulting in the interesting
consequence of the nonadditivity of the Casimir effect. Considering two
periodic dielectric gratings, Lambrecht and Marachevsky have presented an
exact calculation to obtain the Casimir energy l9 . It is worth mentioning
that comparing their theoretical calculation and the measurement performed by
Chan et al. l10  for grating-sphere geometry manifests a meaningful agreement.
Evaluating a modal approach in l11  -based on the scattering approach l12 ;
l12-1 ; l12-2 ; l12-3 ; l12-4 \- the researchers have investigated the finite
temperature Casimir interaction force between two periodic nanostructures and
they have also illustrated the flexibility of this formalism. The Casimir
energy between a plate and a nanostructured surface at arbitrary temperature
has been calculated in the framework of the scattering theory in l13  and as a
significant consequence of this investigation it is illustrated that for
grating geometries the contribution of the thermal part of the Casimir energy
is intensified at small separation distances. Modeling grating as a dielectric
function depending on the space and frequency as well as using variable phase
method, Graham in l14  has presented an appropriate approach to investigate
the Casimir effect for gratings with deep corrugations.

In this work considering two half-space mediums as substrates and a symmetric
array composed of three homogenous layers of alternating materials in between
(see Fig.1), we present reflected and transmitted amplitudes with the
efficient approach of l4-1 . With respect to the properties of a dielectric,
we regard the permittivity of the materials to depend on the frequency in
order to assume a real-world system. There are different approaches to obtain
the Casimir energy, but not all efficient in calculating that for multilayers,
gratings and other similar cases. The divergency occurring in calculating the
Casimir energy dose not appear in the scattering formalism and besides with
the advent of this approach, efficient evaluation of the Casimir energy for a
wide variety of geometries, materials, and external conditions has been made
possible l15 . According to these features, we initiate from a scattering
approach with functional determinant representation at finite temperature to
determine the Casimir energy for the configuration of Fig.1. In a qualitative
sense, we obtain the value of the free energy that agrees with the amount
evaluated by Podgornik and his colleagues in l7  using a different approach.
The ability of the scattering approach which allows one to consider non-
trivial geometries at finite temperature together with a realistic description
of the material properties, permits us to calculate the Casimir energy for a
configuration with a periodic structure in its layers (i.e. inhomogeneous
configuration). This approach makes it possible to investigate such a
configuration in future with the purpose of utilizing that in designing
nanomachines.

Figure 1: A system composed of two half spaces of a material characterized
with the permittivity $\varepsilon_{R}(\omega)$ as substrates within three
alternating homogenous layers described by permittivities
$\varepsilon_{B}(\omega)$ and $\varepsilon_{P}(\omega)$.

##  II Constructing reflection matrices

Let us consider two half-space mediums with a symmetric array composed of
three layers of alternating materials in between, as depicted in Fig.1. $B$,
i.e. the inner homogeneous region, is $0<z<L$ with the permittivity
$\varepsilon_{B}(\omega)$. Over this, there exists another homogeneous region,
$L<z<L+a$, labeled $P$ with the permittivity $\varepsilon_{P}(\omega)$ and
above that we have an $R$-labeled homogeneous region characterized by the
permittivity $\varepsilon_{R}(\omega)$. Below region $B$, symmetric to the
upper half of this system, in the region $-a<z<0$ there is another medium $P$
described with the permittivity $\varepsilon_{P}(\omega)$ and under that we
have another $R$-region characterized by $\varepsilon_{R}(\omega)$.

We consider a plane wave propagating along the $z$-axis reflected and
transmitted through the nearby layers. Considering that this configuration is
$t$, $x$ and $y$ invariant, we can extract a factor $e^{i(k_{x}x+k_{y}y-\omega
t)}$ from all of the components of the electromagnetic fields. In the absence
of the upper half of Fig.1, we introduce the $y$-components of the electric
and magnetic fields in the medium $B$ as:

$\left\\{\begin{matrix}E_{y}^{B}=I_{B}^{e}e^{-i\gamma_{B}z}+R_{B}^{e}e^{i\gamma_{B}z}\\\
H_{y}^{B}=I_{B}^{h}e^{-i\gamma_{B}z}+R_{B}^{h}e^{i\gamma_{B}z}\end{matrix}\right.$
(1)

where $I^{e/h}_{B}$ and $R^{e/h}_{B}$ are the incident and reflection
coefficients of the electric (magnetic) waves. $\gamma_{B}$ indicates the $z$
component of the wave vector in the medium $B$, which can be introduced for
the medium $j=B$, $P$ or $R$ as

$\gamma_{j}=i\bigg{(}-\varepsilon_{j}(\omega)\frac{\omega^{2}}{c^{2}}+k_{x}^{2}+k_{y}^{2}\bigg{)}^{1/2}$
(2)

In this relation $\omega$ is the frequency, $c$ denotes the speed of light and
$k_{x}$ and $k_{y}$ are the longitudinal components of the wave vectors.
Regarding the Drude model, we assume the relative permittivity of the region
$\varepsilon_{j}(\omega)$ as function of frequency including dissipation.

Now we determine the reflection coefficients of Eq. (1). In this regard we
introduce the $y$-component of the fields in the $P$ layer as

$\left\\{\begin{matrix}E_{y}^{P}=T_{P}^{e}e^{-i\gamma_{P}z}+R_{P}^{e}e^{i\gamma_{P}z}\\\
H_{y}^{P}=T_{P}^{h}e^{-i\gamma_{P}z}+R_{P}^{h}e^{i\gamma_{P}z}\end{matrix}\right.$
(3)

including a reflection and also a transmission contribution for both of the
fields. Considering that region $R$ is extremely large in comparison with
other regions of the system, we can write the $y$-components of the electric
and magnetic fields in this region as the following simple form

$\left\\{\begin{matrix}E_{y}^{R}=T_{R}^{e}e^{-i\gamma_{R}(z+a)}\\\
H_{y}^{R}=T_{R}^{h}e^{-i\gamma_{R}(z+a)}\end{matrix}\right.$ (4)

where $T^{e}_{R}$ and $T^{h}_{R}$ are the transmitted electric and magnetic
coefficients through the medium $R$, respectively. Due to the thickness of the
$P$ layer and for the sake of convenience of imposing the boundary condition,
a translation of magnitude $a$ becomes obvious in this relation. We treat the
$y$-components of the electromagnetic fields of different mediums in detail.
Now with the aid of the Maxwell’s curl equations, the $x$-components of the
fields in the medium $j$ are given by

$\left\\{\begin{matrix}E_{x}^{j}=\frac{ik_{y}}{\varepsilon_{j}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}\partial_{x}E_{y}^{j}-\frac{i\omega\mu_{0}}{\varepsilon_{j}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}\partial_{z}H_{y}^{j}\\\
H_{x}^{j}=\frac{ik_{y}}{\varepsilon_{j}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}\partial_{x}H_{y}^{j}+\frac{i\omega\varepsilon_{j}(w)\varepsilon_{0}}{\varepsilon_{j}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}\partial_{z}E_{y}^{j}\end{matrix}\right.$
(5)

where $\varepsilon_{0}$ and $\mu_{0}$ are the physical constants corresponding
to the vacuum permeability and permittivity and $\vec{H}$ denotes the magnetic
field strength.

Now we want to obtain the reflection coefficients of medium $B$, while we have
to consider the upper half of Fig.1 as well as the lower half. In this regard,
we made the following simple change in the variables of Eqs. (1), (3) and (4)
analogous to an $L$ translation along the positive direction of the $z$-axis
and an incident onto the upper surface

$\begin{matrix}z\rightarrow L-z\end{matrix}$ (6)

and therefore, at the interface of medium $B$ and its upper $P$ medium, the
$y$-components of the electric and magnetic fields can be written as

$\left\\{\begin{matrix}E_{y}^{B}=I_{B}^{e}e^{-i\gamma_{B}(L-z)}+R_{B}^{e}e^{i\gamma_{B}(L-z)}\\\
H_{y}^{B}=I_{B}^{h}e^{-i\gamma_{B}(L-z)}+R_{B}^{h}e^{i\gamma_{B}(L-z)}\end{matrix}\right.$
(7)

Here $I_{B}^{e}$ and $I_{B}^{h}$ are the incident coefficients of the electric
and magnetic fields in the medium $B$, respectively. $R_{B}^{e}$ and
$R_{B}^{h}$ that are the reflection coefficients of the electric and magnetic
fields play an important role in constructing the reflection matrix of the
upper half of Fig.1. For medium $P$, one can write the $y$-components of the
electric and magnetic fields analogous to Eq. (3) as the following

$\left\\{\begin{matrix}E_{y}^{P}=T_{P}^{e}e^{-i\gamma_{P}(L-z)}+R_{P}^{e}e^{i\gamma_{P}(L-z)}\\\
H_{y}^{P}=T_{P}^{h}e^{-i\gamma_{P}(L-z)}+R_{P}^{h}e^{i\gamma_{P}(L-z)}\end{matrix}\right.$
(8)

where $T_{P}^{e/h}$ denotes transition from medium $P$ to medium $R$ in the
upper half of Fig.1 and $R_{P}^{e/h}$ describes the reflection from the
interface of regions $P$ and $R$. Finally for the $y$-components of the
electric and magnetic fields of medium $R$ we have

$\left\\{\begin{matrix}E_{y}^{R}=T_{R}^{e}e^{-i\gamma_{R}(L+a-z)}\\\
H_{y}^{R}=T_{R}^{h}e^{-i\gamma_{R}(L+a-z)}\end{matrix}\right.$ (9)

where $T_{R}^{e}$ and $T_{R}^{h}$ are the transition coefficients of medium
$R$. The $x$-components of these fields are given by Eq. (5).

Imposing the continuity boundary condition on the introduced longitudinal
components of both electric and magnetic fields on the interface of the layers
at $z=0$, we have

$\left\\{\begin{array}[]{ll}\frac{-k_{x}k_{y}}{\varepsilon_{B}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(I_{B}^{e}+R_{B}^{e})-\frac{\omega\mu_{0}\gamma_{B}}{\varepsilon_{B}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(I_{B}^{h}-R_{B}^{h})\par=\frac{-k_{x}k_{y}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{e}+R_{P}^{e})-\frac{\omega\mu_{0}\gamma_{P}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{h}-R_{P}^{h})\\\
\\\ I_{B}^{e}+R_{B}^{e}=T_{P}^{e}+R_{P}^{e}\\\ \\\
\frac{-k_{x}k_{y}}{\varepsilon_{B}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(I_{B}^{h}+R_{B}^{h})+\frac{\omega\varepsilon_{B}\varepsilon_{0}\gamma_{B}}{\varepsilon_{B}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(I_{B}^{e}-R_{B}^{e})\par=\frac{-k_{x}k_{y}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{h}+R_{P}^{h})+\frac{\omega\varepsilon_{P}\varepsilon_{0}\gamma_{P}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{e}-R_{P}^{e})\\\
\\\ I_{B}^{h}+R_{B}^{h}=T_{P}^{h}+R_{P}^{h}\end{array}\right.$ (10)

We exert the mentioned boundary condition on the surface $z=-a$ to obtain the
following set of equations

$\left\\{\begin{array}[]{ll}\frac{-k_{x}k_{y}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{e}\,e^{i\gamma_{P}a}+R_{P}^{e}\,e^{-i\gamma_{P}a})-\frac{\omega\mu_{0}\gamma_{P}}{\varepsilon_{p}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{h}\,e^{i\gamma_{P}a}-R_{P}^{h}\,e^{-i\gamma_{P}a})\par=\frac{-k_{x}k_{y}}{\varepsilon_{R}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}T_{R}^{e}-\frac{\omega\mu_{0}\gamma_{R}}{\varepsilon_{R}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}T_{R}^{h}\\\
\\\ T_{P}^{e}\,e^{i\gamma_{P}a}+R_{P}^{e}\,e^{-i\gamma_{P}a}=T_{R}^{e}\\\ \\\
\frac{-k_{x}k_{y}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{h}\,e^{i\gamma_{P}a}+R_{P}^{h}\,e^{-i\gamma_{P}a})+\frac{\omega\varepsilon_{P}\varepsilon_{0}\gamma_{P}}{\varepsilon_{P}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}(T_{P}^{e}\,e^{i\gamma_{P}a}-R_{P}^{e}\,e^{i\gamma_{P}a})\par=\frac{-k_{x}k_{y}}{\varepsilon_{R}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}T_{R}^{h}+\frac{\omega\varepsilon_{R}\varepsilon_{0}\gamma_{R}}{\varepsilon_{R}(\omega)\omega^{2}/c^{2}-k_{y}^{2}}T_{R}^{e}\\\
\\\
T_{P}^{h}\,e^{i\gamma_{P}a}+R_{P}^{h}e^{-i\gamma_{P}a}=T_{R}^{h}\end{array}\right.$
(11)

Considering that the incident coefficients of the electric and magnetic fields
are indeed inputs of our problem, we first assume $I_{B}^{e}=1$ and
$I_{B}^{h}=0$ to solve this system of equations with the purpose of obtaining
desired variables $R_{B}^{e}$ and $R_{B}^{h}$ as the elements of the first
column of the reflection matrix of the inner medium analogous to the lower
half of Fig.1, i.e. $R_{l}(\omega)$. Then assuming $I_{B}^{e}=0$ and
$I_{B}^{h}=1$, one is able to obtain the other column of this reflection
matrix. We construct the reflection matrix of region $B$ for the
electromagnetic wave reflected from the lower nearby layers with the following
recipe

$R_{l}(\omega)=\left(\begin{array}[]{cc}R_{B}^{e}(I_{B}^{e}=1)&R_{B}^{e}(I_{B}^{h}=1)\\\
R_{B}^{h}(I_{B}^{e}=1)&R_{B}^{h}(I_{B}^{h}=1)\\\ \end{array}\right)$ (12)

where, as an instance, $R_{B}^{e}(I_{B}^{e}=1)$ refers to assuming
$I_{B}^{e}=1$ along with $I_{B}^{h}=0$ in Eqs. (10) and (11) to obtain
$R_{B}^{e}$ in these system of equations. To construct $R_{u}(\omega)$ that is
the reflection matrix due to the upper half of our system, we have performed
the same procedure for the tangential components of the electromagnetic fields
corresponding to the upper half of the system introduced in Eqs. (7), (8) and
(9).

##  III Calculating Casimir force in region $B$

In this section we will use the scattering approach to obtain the Casimir
force for the system composed of two substrates within three homogeneous
alternating materials (see Fig.1). For this purpose we construct a matrix M
depicting a full-round trip through medium $B$ as

$\textit{M}(\omega)=R_{u}(\omega)e^{i\gamma_{B}L}R_{l}(\omega)e^{i\gamma_{B}L}$
(13)

where $L$ denotes the thickness of medium $B$, $e^{i\gamma_{B}L}$ refers to
the free photon propagation in this medium and $R_{u}(\omega)$ and
$R_{l}(\omega)$ are the upper and lower reflection matrices that we have
presented in the previous section. Assuming $\omega=i\zeta_{n}$ with
$\zeta_{n}=\frac{2\pi nk_{B}T}{\hbar}$, we utilize the Matsubara frequency to
obtain the finite temperature Casimir force per unit area from the following
equation

$\mathcal{F}(L)=2\pi
k_{B}T{\sum_{n=0}^{\infty}}{{}^{\prime}}\int_{-\infty}^{+\infty}\int_{-\infty}^{\infty}tr[(1-\textit{M}_{n})^{-1}\partial_{L}\textit{M}_{n}]dk_{x}dk_{y}$
(14)

where $k_{B}$ is the Boltzmann constant, $T$ denotes the temperature and $n$
is an integer (see l9 ; l13  and references there in). To proceed, we assume
materials of regions $B$, $P$ and $R$ to be water, lipid and water in order to
follow the configuration presented in l7  closely. Performing the calculations
for this special case as it is illustrated in Fig.2, we obtain the Casimir
force in the $SI$ system of units as

Figure 2: Plot of the Casimir interaction force as a function of the
dimensionless parameter $l=L/a$ in units of $F_{0}=10^{-11}\>N$ for the
configuration of Fig.1 at room temperature. Here we have considered water and
lipid for the alternating layers with $a=4\>nm$ as the thickness of medium
$B$. As it is expected for large values of the dimensionless parameter $l$,
this Casimir force converges rapidly due to the feature of the Casimir effect.

$F(l)=-6.176\times 10^{-9}/\,l^{5}$ (15)

in which $l$ is a dimensionless parameter introduced as $l=L/a$ and $a$ is the
thickness of medium $P$. Consequently the Casimir energy of this configuration
is proportional to $L^{-4}$ which is the same as the energy dependence
obtained in l7  for the investigated configuration. We should emphasize that
in the asymptotic limit $L\gg a$, assuming $a=4\>nm$ and $L=100\>nm$ a free
energy equal to $1.5\times 10^{-23}\>j$ has been obtained in l7  which is in a
good agreement with the value of the Casimir energy that we have computed for
the parameter $l=L/a=100/4$.

##  IV Conclusion

In this paper applying the scattering approach along with the Matsubara
formalism, we have presented a calculation of the finite temperature Casimir
force for a system composed of two half space materials with a three-layer
array of alternating materials between them. We have considered Drude model
for dielectric functions in order to take the materials’e dissipation and
frequency dependence into account. Introducing longitudinal components of the
electromagnetic fields in different regions of Fig.1, we have imposed the
continuity boundary condition on these fields to construct the reflection
matrices. We have considered mediums of water and lipid to follow the
configuration presented in l7 , closely. Assuming $l=L/a=100/4$ with $L$ and
$a$ as the thickness of mediums $B$ and $P$, we have computed a Casimir energy
equal to $1.5\times 10^{-23}\>j$ at room temperature. This result fits with
the value of the energy corresponding to $a=4\>nm$ and $L=100\>nm$ in the
$L\gg a$ limit of l7 . This agreement illustrates the validity of our simple
and transparent procedure which is fortunately able to evaluate the Casimir
energy for structures with inhomogeneous layers which are recognized as good
candidates for engineering Casimir effect with the purpose of designing
miniaturized systems.

## References

  * (1) H. B. G. Casimir, Proc. K. Ned. Akad. Wet. 51, 793 (1948).
  * (2) M.F. Miri and R. Golestanian, Appl. Phys. Lett. 92, 113103(2008).
  * (3) M.F. Miri, V. Nekouie and R. Golestanian, Phys. Rev. E 81, 016104 (2010).
  * (4) M. Nasiri, A. Moradian and M.F. Miri, Phys. Rev. E. 82, 037101 (2010).
  * (5) A. Moradian and M.F. Miri, Eur. Phys. J. E. 34, 12 (2011).
  * (6) J. A. Jaber and J. B. Schlenoff, Current Opinion in Colloid and Interface Science 11, 324 (2006).
  * (7) R. Podgornik,P. L. Hansen, and V. A. Parsegian,J. Chem. Phys. 119, 1070 (2003).
  * (8) R. Podgornik,R. H. French, and V. A. Parsegian,J. Chem. Phys. 124, 044709 (2006).
  * (9) A. Lambrecht and V. N. Marachevsky, Phys. Rev. Lett. 101, 160403 (2008).
  * (10) H. B. Chan, Y. Bao, J. Zou, R. A. Cirelli, F. Klemens, W. M. Mansfield, and C. S. Pai, Phys. Rev. Lett. 101, 030401 (2008).
  * (11) P. S. Davids, F. Intravaia, F. S. Rosa, and D. A. Dalvit, Phys. Rev. A 82, 062111 (2010).
  * (12) M.T. Jaekel and S. Reynaud, J. Phys. I France 1, 1395 (1991).
  * (13) E. I. Kats, Sov. Phys. JETP 46, 109 (1977).
  * (14) A. Lambrecht, P. A. M. Neto, and S. Reynaud, New J. Phys. 8, 243 (2006).
  * (15) O. Kenneth and I. Klich, Phys. Rev. Lett. 97, 160401 (2006).
  * (16) T. Emig, N. Graham, R. L. Jaekel, and M. Kardar, Phys. Rev. Lett. 99, 170403 (2007).
  * (17) R. Guérout, J. Lussange, H. B. Chan, A. Lambrecht, and S. Reynaud, Phys. Rev. A 87, 052514 (2013).
  * (18) N. Graham, Phys. Rev. A 90, 032507 (2014).
  * (19) M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, J. Opt. Soc. Am. A 12, 1077 (1995).
  * (20) S. J. Rahi, T. Emig, and R. L. Jaffe, Lecture Notes in Physics, Preprint arXiv:1007.4355 (2010).


# Multiple states of two-dimensional turbulence above topography

Jiyang He1,2  Yan Wang1,2<EMAIL_ADDRESS>1Department of Ocean Science, The
Hong Kong University of Science and Technology, Hong Kong, China2Center for
Ocean Research in Hong Kong and Macau, The Hong Kong University of Science and
Technology, Hong Kong, China.

###### Abstract

The recent work of Siegelman & Young (PNAS, vol. 120(44), 2023, pp.
e2308018120) revealed two extreme states reached by the evolution of unforced
and weakly-damped two-dimensional turbulence above random rough topography,
separated by a critical kinetic energy $E_{\\#}$. The low- and high-energy
solutions correspond to topographically-locked and roaming vortices,
surrounded by non-uniform and homogeneous background potential vorticity (PV),
respectively. However, we found that these phenomena are restricted to the
particular intermediate length scale where the energy was initially injected
into the system. Through simulations initialised by injecting the energy at
larger and smaller length scales, we found that the long-term state of
topographic turbulence is also dependent on the initial length scale. If the
initial length scale is comparable with the domain size, the flow field
evolves to the minimum-enstrophy state proposed by Bretherton and Haidvogel
(J. Fluid Mech., vol.78(1), 1976, pp. 129-154), without the presence of
vortices, especially when the energy is no larger than $E_{\\#}$. As the
initial length scale becomes smaller, more vortices nucleate, and become more
mobile. Simultaneously, the background PV tends to homogenization, even if the
energy is below $E_{\\#}$. These results complement the phenomenology of
topographic turbulence documented by Siegelman & Young, by showing that the
minimum-enstrophy and background PV homogenization states can be reached by
large- and small-scale initial fields, respectively, whatever the energy is.

###### keywords:

##  1 Introduction

Two-dimensional, rapidly-rotating turbulence above topography (topographic
turbulence hereafter) serves as an effective reduced model for studying large-
scale oceanic motions over rough and sloping seafloor. This model has provided
insight into the emergence of prevailing undercurrents over continental
shelves and slopes (Wang & Stewart, 2018) and persistent along-bathymetry
flows above topographic depressions and seamounts (Solodoch et al., 2021).
Topographic turbulence is also relevant to long-lived, anticyclonic vortices
locked to large-scale topographic bowls in the ocean, such as the Lofoten
Basin Eddy and the Mann Eddy in the North Altantic (see Köhl, 2007; LaCasce et
al., 2024). Since current climate models remain to fall short of resolving
turbulent flows in the ocean, employing the topographic turbulence model to
investigate phenomena arising from it will facilitate the parameterization of
ocean turbulence over bumpy seabeds in these models (Holloway, 1992; Radko,
2023; Eaves et al., 2024).

Theories of topographic turbulence aimed to derive the steady states via
variational principles. For an energy-conserving system, the potential
enstrophy is minimized (Bretherton & Haidvogel, 1976), leading to a linear
relation between the potential vorticity (PV hereafter) and geostrophic
streamfunction. The other approach is to maximize the entropy and is thus
related to the equilibrium statistical mechanics (Salmon et al., 1976). The
ensembles of PV and streamfunction also satisfy a linear relation. According
to Carnevale & Frederiksen (1987), the two solutions are actually consistent
with each other, and nonlinearly stable.

Recent numerical experiments show that topographic turbulence shares some
common features with flat-bottom two-dimensional turbulence (2DT). For
instance, condensated states commonly found in 2DT (Boffetta & Ecke, 2012)
have been identified in topographic turbulence over small-scale topography
(Zhang & Xie, 2023; Gallet, 2024). Siegelman & Young (2023) (SY23 hereafter)
showed the coexistence of vortices and a background flow in topographic
turbulence, the former of which are prevalent in 2DT (Mcwilliams, 1984, 1990).
They found that the mobility of vortices and the uniformity of background
flows depend on the energy of the system: on the low-energy branch, the
vortices are locked to topographic features, and the background PV is
spatially non-uniform; on the high-energy branch, the vortices roam throughout
the domain, and the background PV is homogenized (Rhines & Young, 1982). These
two branches are separated by a critical energy.

SY23 examined the minimum-enstrophy hypothesis of Bretherton & Haidvogel
(1976) against their numerical simulations, and identified significant
discrepancies: the high-energy solutions displayed homogeneous background PV
that is not predicted by the minimum-enstrophy hypothesis at the same high
energy; the long-term total enstrophy far exceeded that of the theoretical
prediction for almost all energies. SY23 attributed these discrepancies to the
presence of long-lived vortices: the roaming of vortices mixes the background
PV toward homogeneity, and vortex shielding inhibits enstrophy from cascading
to smaller scales.

In this work, we will show that the emergence of vortices might relate to the
specific intermediate length scale of the initial fields prescribed by SY23,
at which the initial total enstrophy is much higher than the theoretical
minimum value and the large margin may serve as the seeds of vortices. As
opposed to fixing an initially energized wavenumber, we conduct numerical
experiments similar to those of SY23 but initialised by monoscale fields
sweeping a broad range of wavenumbers, from scales comparable to the domain
size to those much smaller than that chosen by SY23. The purpose here is to
investigate how topographic turbulence evolves in response to different scales
of initial conditions: can the emergence of vortices be inhibited and the
minimum-enstrophy state be reached if the initial scale is adequately large;
what would occur if the initial scale is smaller than that of SY23?

##  2 Framework

###  2.1 Theory

As in SY23, we consider an unforced, single-layer quasi-geostrophic (QG) flow
on a $f$-plane in a doubly-periodic domain ($L\times L$, with $L$ denoting the
domain size). The governing equation of the QG PV $q(x,y,t)$ reads

$\frac{\partial q}{\partial t}+\frac{\partial\psi}{\partial x}\frac{\partial
q}{\partial y}-\frac{\partial\psi}{\partial y}\frac{\partial q}{\partial
x}=D\zeta,\quad q=\zeta+\eta=\nabla^{2}\psi+\eta,$ $None$

where $\psi$ and $\zeta$ are the geostrophic streamfunction and relative
vorticity, respectively. The operator $D$ stands for the dissipation of PV
that will be specified later. Effects from bottom topography are encoded in
the topographic PV, $\eta(x,~{}y)=-f_{0}h_{1}/h_{0}$, for a total depth of
$h_{0}+h_{1}(x,y)$ with small fluctuations $h_{1}(x,y)$ around a constant
average $h_{0}$, scaled by the local Coriolis frequency $f_{0}$ (see SY23).

In the absence of dissipation, the system (2.1) conserves the kinetic energy
and total enstrophy, defined respectively as

$E=\frac{1}{2}\left<|\nabla\psi|^{2}\right>=\frac{1}{2L^{2}}\iint|\nabla\psi|^{2}dxdy\quad\mbox{and}\quad
Q=\frac{1}{2}\left<q^{2}\right>=\frac{1}{2L^{2}}\iint q^{2}dxdy.$ $None$

Here $\left<\cdot\right>$ represents an domain-average operation. When the
added dissipation in the governing equation (2.1) operates at high
wavenumbers, the energy is approximately conserved, while the enstrophy drops
by a large amount as a result of forward cascade. This motivated Bretherton &
Haidvogel (1976) to predict that an energy-conserving system would evolve into
a minimum-enstrophy state. The state can be obtained by performing variational
calculation of the functional $\mathcal{L}=Q+\mu E$, with $\mu$ denoting the
Lagrange multiplier. The resulting Euler-Lagrange equation,

$q_{\ast}=\nabla^{2}\psi_{\ast}+\eta=\mu\psi_{\ast},$ (3)

dictates a linear relation between $q_{\ast}$ and $\psi_{\ast}$. Solutions to
(3) are actually inviscid steady solutions to the QG PV equation (2.1).

To solve the equation (3), the topographic PV can be Fourier-transformed into

$\eta(x,~{}y)=\sum_{\boldsymbol{k}}\eta_{\boldsymbol{k}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{x}},$
(4)

where $\boldsymbol{k}$ represents the wavenumber vector. As in SY23, the
topographic PV is chosen to satisfy a power spectrum of $k^{-2}$
($k=|\boldsymbol{k}|$), by which the Fourier amplitudes in the expression (4)
shall take the form

$\eta_{\boldsymbol{k}}=\alpha\mathrm{e}^{\mathrm{i}\phi_{\boldsymbol{k}}}k^{-3/2}.$
(5)

The normalization parameter $\alpha$ is adjusted to produce a specific value
of the root mean square of the topographic PV, defined as
$\eta_{rms}=\sqrt{\sum_{\boldsymbol{k}}|\eta_{\boldsymbol{k}}|^{2}}$.
Following SY23, we prescribe $\eta_{rms}=10^{-6}~{}\mathrm{s}^{-1}$.
Meanwhile, the random phase $\phi_{\boldsymbol{k}}$ in (5) serves to produce a
realization of $\eta$. Throughout this work, the topography remains fixed, and
is shown in figure 1.

Figure 1: Topographic potential vorticity
$\eta(x,y)(\mathrm{unit:}~{}10^{-6}~{}\mathrm{s}^{-1})$ used throughout this
work.

With the prescribed topography, the solution to the Euler-Lagrange equation
(3) can be obtained in spectral space as

$\psi_{\ast}(x,y,\mu)=\sum_{\boldsymbol{k}}\frac{\eta_{\boldsymbol{k}}\mathrm{e}^{\mathrm{i}\boldsymbol{k}\cdot\boldsymbol{x}}}{\mu+k^{2}}.$
(6)

The corresponding energy and minimum enstrophy are then

$E=\frac{1}{2}\sum_{\boldsymbol{k}}\frac{k^{2}|\eta_{\boldsymbol{k}}|^{2}}{\left(\mu+k^{2}\right)^{2}}\quad\mbox{and}\quad
Q_{min}=\frac{1}{2}\sum_{\boldsymbol{k}}\frac{\mu^{2}|\eta_{\boldsymbol{k}}|^{2}}{\left(\mu+k^{2}\right)^{2}},$
$None$

respectively. With a given energy level $E$ for an energy-conserving system,
the Lagrange multiplier $\mu$ and minimum enstrophy $Q_{min}$ can be sought
through the above relations.

Figure 2$(a,b)$ show the resulting $\mu$ and $Q_{min}$ as functions of $E$,
respectively, scaled by $k_{1}^{2}$ ($k_{1}=2\pi/L$ represents the fundamental
wavenumber) and $Q_{\eta}=\eta_{rms}^{2}/2$. In addition, the energy has been
scaled by the critical energy level $E_{\\#}$, corresponding to the particular
case of $\mu=0$ and thus of homogenized PV ($q_{\\#}=0$). At this critical
energy, the minimum enstrophy vanishes ($Q_{min}=0$; see figure 2$(b)$). The
critical energy $E_{\\#}$ separates the low- and high-energy branches, which
correspond to positive and negative values of $\mu$, respectively. These
theoretical solutions are identical to those shown in SY23 (see their figure
2$(a,b)$), who further found that $E_{\\#}$ is a good separator of high- and
low-energy branches of their numerical experiments.

Figure 2: Minimum-enstrophy solutions as in SY23 (see their Figures 2(a)–(b)).

###  2.2 Numerical setups

We solve the QG PV equation (2.1) using the open-source pseudo-spectral
package GeophysicalFlows.jl running on GPUs (Constantinou et al., 2021). Time
is stepped forward by a fourth-order Rounge-Kutta scheme. The dissipation term
($D\zeta$) in (2.1) is implemented via a spectral filter of relative
vorticity, which is applied to high wavenumbers at the end of each time step.
We choose exactly the same parameter setups as in SY23: the domain size is
$10^{6}$ meters; the resolution is $1024\times 1024$; the integration time is
$1.5\times 10^{9}$ seconds (around $47.53$ years); time step is $1500$
seconds, which is reduced further to $500$ seconds for computationally
demanding runs.

Figure 3: Variations of initial enstrophy $Q(0)$ with the initial wavenumber
$k_{ini}$ for different energies, compared with the minimum enstrophy
$Q_{min}$.

We initialise the computations by a series of random monoscale relative
vorticity fields with a wide range of enstrophy, which is achieved by
different combinations of the wavenumber and energy level. For all cases, the
Fourier components of the initial relative vorticity $\zeta(x,y,0)$ are chosen
randomly within the bandwidth $[k_{ini}-0.5k_{1},k_{ini}+0.5k_{1}]$, where
$k_{ini}$ is the wavenumber of the initial monoscale field. In practice, we
choose a series of

$k_{ini}\in\\{1,2,4,8,12,16,24,32,48,64\\}k_{1},$ (8)

corresponding to large-to-small initial monoscales, and

$E\in\\{0.05,0.25,0.50,1.00,2.00\\}E_{\\#},$ (9)

corresponding to low-to-high energy levels. Note that SY23 prescribed the
initial monoscale at approximately $16k_{1}$, which is intermediate in our
selected range of wavenumbers (8). The initial enstrophy for different
combinations of $k_{ini}$ and $E$ is shown in figure 3, estimated as
$Q(0)\approx k_{ini}^{2}E+Q_{\eta}$. That is, with a prescribed energy level,
$Q(0)$ increases at a rate of $k_{ini}^{2}$, which means that the larger
amount of initial enstrophy is contained at the smaller scales. As shown in
figure 3, the initial enstrophy exceeds the theoretical minimum enstrophy for
all cases. However, the former is adequately close to the latter in the
small-$k_{ini}$ regime, as opposed to the case of intermediate wavenumber
(highlighted by the vertical dotted-dashed line in figure 3) considered by
SY23. Owing to the initial/minimum enstrophy proximity, small-$k_{ini}$ runs
may readily converge to the minimum-enstrophy states.

##  3 Results

In this section, we analyse the solutions of our numerical experiments.
Particular attention is paid to the dependence of the long-term state on the
initial wavenumber $k_{ini}$. Results of $k_{ini}=16k_{1}$, consistent with
those in SY23, are highlighted in red color upon comparison with other
wavenumber cases.

###  3.1 Time evolution of energy and enstrophy

The time evolution of energy and enstrophy is shown in figure 4 for
simulations of the critical energy $E/E_{\\#}=1$ (results of other energies
exhibit similar behaviours and are not shown). As shown in figure 4($a$), the
energy loss during run time is larger for larger $k_{ini}$, because the
corresponding initial spectral position of energy injection is closer to the
filtering region at high wavenumbers. However, the energy loss at the final
time (47.53 years) amounts to a negligible fraction of the initial energy,
within $0.3\%$ for $k_{ini}\leq 16k_{1}$ and around $3\%$ for
$k_{ini}=64k_{1}$ (see first column of numbers in table 1). Therefore, energy
is nearly conserved for all runs, which is the necessary condition for
examining the minimum-enstrophy hypothesis.

Figure 4: Time evolution of energy ($a$) and enstrophy ($b$) for the $E_{\\#}$ runs. $k_{ini}$ $E$ | $0.05E_{\\#}$ | $0.25E_{\\#}$ | $0.50E_{\\#}$ | $E_{\\#}$ | $2E_{\\#}$
---|---|---|---|---|---
theory | 1.000; 0.457 | 1.000; 0.174 | 1.000; 0.057 | 1.000; 0.000 | 1.000; 0.108
$1{k}_{1}$ | 0.997; 0.465 | 0.999; 0.177 | 0.999; 0.061 | 0.999; 0.017 | 0.999; 0.193
$2{k}_{1}$ | 0.997; 0.472 | 0.999; 0.184 | 0.999; 0.064 | 0.999; 0.006 | 0.999; 0.299
$4{k}_{1}$ | 0.997; 0.475 | 0.998; 0.181 | 0.999; 0.092 | 0.999; 0.030 | 0.999; 2.105
$8{k}_{1}$ | 0.997; 0.468 | 0.998; 0.224 | 0.998; 0.323 | 0.999; 0.220 | 0.999; 2.420
$12{k}_{1}$ | 0.996; 0.480 | 0.998; 0.605 | 0.998; 1.173 | 0.998; 1.321 | 0.998; 5.238
$16{k}_{1}$ | 0.996; 0.582 | 0.997; 1.305 | 0.997; 2.081 | 0.997; 1.816 | 0.998; 8.734
$32{k}_{1}$ | 0.990; 1.170 | 0.992; 3.835 | 0.992; 5.766 | 0.992; 11.76 | 0.992; 20.98
$48{k}_{1}$ | 0.981; 2.022 | 0.983; 6.945 | 0.983; 9.439 | 0.984; 25.73 | 0.984; 36.76
$64{k}_{1}$ | 0.969; 3.115 | 0.971; 11.12 | 0.971; 16.50 | 0.972; 32.41 | 0.972; 63.38
Table 1: Energy $E(47.53~{}\mbox{years})/E(0)$ (first column of numbers) and
enstrophy $Q(47.53~{}\mbox{years})/Q_{\eta}$ (second column of numbers) for
runs with different $k_{ini}$ and $E$, compared against the theoretical
predictions (first row).

By contrast, the enstrophy decreases drastically during run time, as shown in
figure 4($b$). Consistent with the findings of SY23, there is an initial fast
drop followed by a slow decrease. Generally, for a fixed energy level,
larger-$k_{ini}$ runs, or equivalently, runs with higher initial enstrophy
(see figure 3), maintain higher enstrophy during run time. Runs with small
$k_{ini}$ are intriguing: the enstrophy approaches the minimum value $Q_{min}$
during run time. At the final time (see second column of numbers in table 1),
the enstrophy $Q$ of small-$k_{ini}$ ($1k_{1}$, $2k_{1}$, $4k_{1}$) runs
prescribed with low energy ($E/E_{\\#}\leq 1$) is extremely close to the
corresponding minimum value $Q_{min}$. However, for the high energy case of
$2E_{\\#}$, the final enstrophy of the $1k_{1}$ run is almost twice the
minimum value. For the intermediate case of $k_{ini}=16k_{1}$ considered by
SY23 and cases with larger $k_{ini}$, the final enstrophy is much larger than
the minimum value, except for the lowest-energy runs ($0.05E_{\\#}$). Thus,
the minimum-enstrophy states can actually be reached by runs with large-scale
initial conditions, especially by those with low energy ($E/E_{\\#}\leq 1$).
These findings contrast markedly with the conclusions of SY23.

###  3.2 Long-term states

Next we inspect the long-term snapshots in runs with typical energy levels at
the final time (47.53 years), and compare those with the minimum-enstrophy
states $\zeta_{\ast}$ and $q_{\ast}$ obtained theoretically. The corresponding
animations, showing the time evolution to the long-term states, can be found
in the supplementary materials.

Figure 5: Long-term snapshots of $0.25E_{\\#}$ runs for
$k_{ini}=[1,16,64]k_{1}$ (second-fourth columns) compared with the minimum
enstrophy state at the same energy level (first column): first row, $\zeta$;
second row, $q$; third row, $q_{res}$.

Figure 5 compares the snapshots of relative vorticity (first row) and PV
(second row) from three runs with different $k_{ini}$ against the minimum-
enstrophy state at the low energy of $E=0.25E_{\\#}$. The state of the
smallest $k_{ini}=1k_{1}$ run (figure 5($b,f$)) resembles the minimum-
enstrophy state (figure 5($a,e$)), displaying a non-uniform PV field shaped by
the low-pass-filtered topography throughout the domain. This resemblance is
consistent with the proximity of the corresponding enstrophy (see table 1),
and further demonstrates that the minimum-enstrophy state can actually be
reached with very large-scale initial conditions. However, although no strong
vortices emerge, there exist some weak extrema of PV (see figure 5$(f)$) at
the topographic depressions and elevations (see figure 1 for topography). As
$k_{ini}$ increases to $16k_{1}$, strong vortices emerge at these locations,
with cyclones and anticyclones locked to topographic elevations and
depressions, respectively (figure 5($g$)). The background relative vorticity
becomes stronger (compare figure 5($c$) with ($b$)), leading to the weaker
background PV (compare figure 5($g$) with ($f$)). As $k_{ini}$ increases
further to $64k_{1}$, more vortices show up. Notably, the background relative
vorticity strengthens further and the background PV is nearly homogenized
(figure 5($d,h$)). A larger initial wavenumber provides more seeds of
vortices; more vortices nucleate afterwards and mix the background PV more
efficiently to homogeneity. Note that the vortices in the $64k_{1}$ run remain
locked to fixed locations but exhibit higher mobility compared to those in the
$16k_{1}$ run (see supplementary movies of the corresponding runs). In
summary, at a low energy level, the resemblance with the minimum-enstrophy
state and the near homogenization of background PV are observed in simulations
with two limiting initial length scales, with an emergent transitional state
as the initial length scale varies in between. These phenomena cannot be
observed in low-energy simulations with a single intermediate wavenumber of
$16k_{1}$ considered by SY23.

Now we turn to the runs with the critical energy level of $E_{\\#}$, as shown
in figure 6. Once again, the solutions of the small-$k_{ini}$ runs are similar
to the minimum-enstrophy state, with nearly homogenized PV (see figure
6$(h,i)$). This demonstrates that the minimum-enstrophy state can be reached,
together with the analogous enstrophy produced numerically and theoretically
(table 1). As $k_{ini}$ increases, more vortices nucleate and roam throughout
the domain. The background PV is then completely homogenized.

Figure 6: Long-term snapshots of $E_{\\#}$ runs for
$k_{ini}=[1,2,8,16,64]k_{1}$ (second-sixth columns) compared with the minimum
enstrophy state at the same energy level (first column): first row, $\zeta$;
second row, $q$. Figure 7: Same as figure 6 but for $2E_{\\#}$ runs.

A parallel comparison is made further among simulations prescribed with the
high energy level of $2E_{\\#}$, shown in figure 7. As discussed earlier, the
final-state enstrophy does not approach the minimum value, even for the case
of the smallest initial wavenumber $1k_{1}$ (table 1). However, the solution
of the $1k_{1}$ run is still similar to the minimum enstrophy state to some
extent (figure 7$(g,h)$), both displaying non-uniform PV fields imprinted by
the low-pass-filtered topography. However, as shown in figure 7$(h,i)$ for
runs of $1k_{1}$ and $2k_{1}$, a cyclone and a anticyclone are locked
respectively to the topographic depression and elevation, exactly opposing the
relation between vortex polarity and depth in the low-energy runs. These
gigantic vortices inhibit enstrophy from further decaying toward the
theoretical minimum value. Toward the case of $8k_{1}$, the vortices become
smaller in size and more localised; they orbit around the large-scale features
of topography (see the corresponding supplementary movie). The background PV
is essentially homogenized. As $k_{ini}$ increases further, more vortices
emerge and tend to roam throughout the domain. Note that, although the
background PV of the $64k_{1}$ run is homogenized, it is not identically zero.

###  3.3 Empirical background flow

Figure 8: Variations of the empirical slope $\mu_{emp}$ between $q$ and $\psi$
with the initial wavenumber $k_{ini}/k_{1}$ at four energy levels: $(a)$
$E=0.25E_{\\#}$; $(b)$ $E=0.50E_{\\#}$; $(c)$ $E=E_{\\#}$; $(d)$ $E=2E_{\\#}$.
Horizontal dashed lines represent the theoretical Lagrange multiplier based on
the minimum-enstrophy principle at the same energy levels. Vertical dotted-
dashed lines highlight the wavenumber of $k_{init}/k_{1}=16$ considered in
SY23.

As in SY23, we extract the background flow from the data via linear fitting
between $q$ and $\psi$. The vortices in the flow field are outliers and
removed by a threshold. The points with $|q|/\eta_{rms}$ below the threshold
are retained to extract the background state. Then applying least-squares
fitting to $q$ and $\psi$ of the remaining points yields an approximate linear
relation,

$q\approx\mu_{emp}\psi.$ (10)

When extracting the empirical slope $\mu_{emp}$ for each time instant, we
found that its value does not vary significantly in time after an initial
adjustment phase. The empirical slope $\mu_{emp}$ should depend on the energy
$E$, as illustrated by SY23, and the initial wavenumber $k_{ini}$ based on our
experiments. The results of $\mu_{emp}$ are shown in figure 8. We choose three
thresholds for extracting $\mu_{emp}$, and find only quantitative differences.
The runs with the critical energy $E_{\\#}$ show that $\mu_{emp}$ is always
close to zero, indicating near PV homogenization. For small-$k_{ini}$ runs,
the empirical slope $\mu_{emp}$ is close to the theoretical $\mu(E)$
(horizontal dashed line) shown in figure 2$(a)$. Even at the high energy level
of $2E_{\\#}$, a negative $\mu_{emp}$ is consistent with the theoretical
$\mu(2E_{\\#})$. This is another justification of the existence of the minimum
enstrophy state achievable by numerical simulations but not obtained by SY23.
As the initial wavenumber increases, the magnitude of $\mu_{emp}$ tends to
zero, implying a transition from a non-uniform background PV field imprinted
by the low-pass-filtered topography to PV homogenization, consistent with the
preceding observations from the long-term snapshots. The numerical solutions
of SY23 (see vertical dotted-dashed lines) might capture a confined range of
the transition regime, particularly in the low-energy cases.

As in SY23, the empirical slope $\mu_{emp}$ is used to define a “residual” PV
field, $q_{res}=q-\mu_{emp}\psi$, which is shown in figure 5$(i,j,k)$.
Evidently, the background flow is sufficiently removed, thus confirming the
$q$-$\psi$ linear relation (10).

##  4 Conclusion

In this work, we study the states of unforced, weakly-decaying two-dimensional
turbulence above random rough topography, initialised by a series of monoscale
fields covering a broad range of length scales and energy levels. Our
observations complement those of Siegelman & Young (2023), who took a single,
intermediate monoscale for initialisation into account. Some of the
phenomenological descriptions of topographic turbulnece from Siegelman & Young
(2023) shall be updated. The minimum-enstrophy solutions of Bretherton &
Haidvogel (1976) can be reached by initial fields whose length scales are
comparable with the domain size, especially when the energy is no larger than
$E_{\\#}$. As the initial length scale decreases, the higher initial enstrophy
provides more seeds for vortex nucleation. More emerging vortices shield
enstrophy from cascading toward high wavenumbers, and efficiently mix the
background PV. A non-uniform background PV field imprinted by large-scale
topography, at low or high energy, is thus weakened and tends to
homogenization. Our results highlight the sensitivity of topographic
turbulence to the length scales of initial conditions. Given that the depth-
invariant mode of geostrophic turbulence responsible for driving mass and heat
transports was constantly found to be energized across a range of wavenumbers
(Chen, 2023), this work calls for ongoing efforts for refining
parameterisations of ocean turbulence over bumpy seafloor by taking the
energization length scales into account.

[Acknowledgements]This work is supported by the Research Grants Council (RGC)
of Hong Kong under awards Early Career Scheme 26307720 and General Research
Fund 16305321, and by the Center for Ocean Research (CORE), a joint research
center between Laoshan Laboratory and HKUST.

[Declaration of interests]The authors report no conflict of interest.

## References

  * Boffetta & Ecke (2012) Boffetta, G. & Ecke, R. E. 2012 Two-Dimensional Turbulence. Annu. Rev. Fluid Mech. 44 (1), 427–451.
  * Bretherton & Haidvogel (1976) Bretherton, F. P. & Haidvogel, D. B. 1976 Two-dimensional turbulence above topography. J. Fluid Mech. 78 (1), 129–154.
  * Carnevale & Frederiksen (1987) Carnevale, G. F. & Frederiksen, J. S. 1987 Nonlinear stability and statistical mechanics of flow over topography. J. Fluid Mech. 175 (-1), 157.
  * Chen (2023) Chen, S.-N. 2023 Revisiting the baroclinic eddy scalings in two-layer, quasigeostrophic turbulence: Effects of partial barotropization. J. Phys. Oceanogr. 53 (3), 891 – 913.
  * Constantinou et al. (2021) Constantinou, Navid C., Wagner, Gregory LeClaire, Siegelman, Lia, Pearson, Brodie C. & Palóczy, André 2021 GeophysicalFlows.jl: Solvers for geophysical fluid dynamics problems in periodic domains on CPUs & GPUs. J. Open Source Softw. 6 (60), 3053.
  * Eaves et al. (2024) Eaves, R. E., Maddison, J. R., Marshall, D. P. & Waterman, S. 2024 An energy and enstrophy constrained parameterization of barotropic eddy potential vorticity fluxes. Preprint. Preprints.
  * Gallet (2024) Gallet, B. 2024 Two-dimensional turbulence above topography: Condensation transition and selection of minimum enstrophy solutions, arXiv: 2404.06475.
  * Holloway (1992) Holloway, G. 1992 Representing Topographic Stress for Large-Scale Ocean Models. J. Phys. Oceanogr. 22 (9), 1033–1046.
  * Köhl (2007) Köhl, A. 2007 Generation and Stability of a Quasi-Permanent Vortex in the Lofoten Basin. J. Phys. Oceanogr. 37 (11), 2637–2651.
  * LaCasce et al. (2024) LaCasce, J. H., Palóczy, A. & Trodahl, M. 2024 Vortices over bathymetry. J. Fluid Mech. 979, A32.
  * Mcwilliams (1984) Mcwilliams, J. C. 1984 The emergence of isolated coherent vortices in turbulent flow. J. Fluid Mech. 146, 21–43.
  * Mcwilliams (1990) Mcwilliams, J. C. 1990 The vortices of two-dimensional turbulence. J. Fluid Mech. 219 (-1), 361.
  * Radko (2023) Radko, T. 2023 A generalized theory of flow forcing by rough topography. J. Fluid Mech. 961, A24.
  * Rhines & Young (1982) Rhines, P. B. & Young, W. R. 1982 Homogenization of potential vorticity in planetary gyres. J. Fluid Mech. 122 (-1), 347.
  * Salmon et al. (1976) Salmon, R., Holloway, G. & Hendershott, M. C. 1976 The equilibrium statistical mechanics of simple quasi-geostrophic models. J. Fluid Mech. 75 (4), 691–703.
  * Siegelman & Young (2023) Siegelman, L. & Young, W. R. 2023 Two-dimensional turbulence above topography: Vortices and potential vorticity homogenization. PNAS 120 (44), e2308018120.
  * Solodoch et al. (2021) Solodoch, A., Stewart, A. L. & McWilliams, J. C. 2021 Formation of Anticyclones above Topographic Depressions. J. Phys. Oceanogr. 51 (1), 207–228.
  * Wang & Stewart (2018) Wang, Y. & Stewart, A. L. 2018 Eddy dynamics over continental slopes under retrograde winds: Insights from a model inter-comparison. Ocean Model. 121, 1–18.
  * Zhang & Xie (2023) Zhang, L.-F. & Xie, J.-H. 2023 Spectral condensation in quasi-geostrophic turbulence above small-scale topography, arXiv: 2311.16612.