identity on $K$ and $K_{j}=\operatorname{span}(\sum_{i\in I_{j}}Y_{i})$ for
all $j\in\\{1,\ldots k\\}$. Thus $\tilde{D}_{P}f(x)\in H^{\prime}$ for a.e.
$x\in U$. If $H^{\prime}=\\{\operatorname{id}\\}$ we are done. So we assume
that $H^{\prime}\neq\\{\operatorname{id}\\}$ and will show that
$\tilde{D}_{P}f(x)$ acts as a constant permutation on each $H^{\prime}$-orbit
$I_{j}$, for $j\in\\{1,\ldots,k\\}$. Then we can argue as in the case $\dim
K=1$ to deduce constancy of the permutation.

We provide some details for the convenience of the reader. Let $I$ be an orbit
and $\tilde{V}_{2,I}=\oplus_{i\in I}\tilde{V}_{2,i}$. The map
$\tilde{D}_{P}f(x)$ preserves $\tilde{V}_{2,I}$ and acts as the identity on
the subspace $K_{I}=\operatorname{span}\\{\sum_{i\in I}\tilde{Y}_{2,i}\\}$.
Thus the argument in the proof of Proposition 12.7 shows that the identities
in that proposition hold with $\lambda=1$. Recall that
$P_{mi}(x)=\delta_{m\sigma_{x}^{-1}(i)}$. Applying the pull-back theorem to
the forms $f_{P}^{*}\omega_{ij}\wedge d(\varphi i_{X}\beta)$ where
$\omega_{ij}$ and $\beta$ are as in (12.10) and (12.13) with $i,j,m\in I$ and
where $X\in V_{1,l}$ with $l\neq m$ we get as in the proof of Lemma (12.18)
the distributional identity $XP_{mi}=0$ for all $m,i\in I$ and all $X\in
V_{1,l}$ with $l\neq m$. Since $\sigma_{x}(I)=I$ we have $\sum_{m\in
I}P_{mi}=1$ for all $i\in I$. Now for $X\in V_{1,m}$ we get
$XP_{mi}=X(1-\sum_{m^{\prime}\in I\setminus\\{m\\}}P_{m^{\prime}i})=0$. Thus
$P_{mi}$ is constant for $m,i\in I$ and hence $\sigma_{x}|_{I}$ is constant. ∎

##  14\. Product quotients of complex Heisenberg groups and higher real
Heisenberg groups

In this section we prove Theorem 11.1 for product quotients
$G=\tilde{G}/\exp(K)$ where $\tilde{G}$ is a product of copies of a higher
real Heisenberg group or complex Heisenberg group. We retain the notation for
product quotients from Section 11. The approach used in this section is
different, and conceptually simpler, than the ones used in Sections 12 and 13.
The main difference between these two situations is that for the higher real
and the complex Heisenberg group there exist linear independent vectors $X,Y$
in $V_{1,i}$ with $[X,Y]=0$ while this is not the case for the first
Heisenberg group. More systematically consider the map

(14.1) $L_{[\,]}:\Lambda_{2}V_{1}\to V_{2},\qquad L_{[\,]}(X\wedge Y)=[X,Y].$

For the complex Heisenberg groups and higher real Heisenberg groups

(14.2) $\ker_{i}:=\ker L_{[\,]}\cap V_{1,i}$

is non-trivial (as follows from a dimension count) while for the first real
Heisenberg group $L_{[\,]}:\Lambda_{2}V_{1,i}\to V_{2}$ is injective. Any
graded automorphism – in particular the Pansu differential – permutes the
subspaces $\ker_{i}$. Now for complex or higher real Heisenberg product
quotients we can pull back a collection of forms $\alpha\in\Lambda^{2}V_{1,i}$
which restrict to a basis of the dual space $\ker_{i}^{\prime}$ and take
$\beta=i_{X}i_{Z}\omega$ where $X\in\ker_{k}$ and $Z\in\oplus_{k^{\prime}\neq
k}V_{1,k^{\prime}}$. The choice of this collection of forms $\alpha$ may look
somewhat adhoc at this point, but it is in fact related to our general
strategy to look for closed forms modulo exact forms. The point is that the
map $L_{[]}$ is the dual of the exterior differential
$d:\Lambda_{1}V_{2}\to\Lambda^{2}V_{1}\simeq(\Lambda_{2}V_{1})^{\prime}$, see
Remark 14.14 after the end of the proof. With the above choice of $\alpha$ and
$\beta$ a short argument similar to the one used for products shows that the
permutation of the subspaces induced by $D_{P}f$ must be locally constant. The
subspace $K$ plays no role in this argument. In fact the argument can also be
used to obtain an alternative proof for products of the complex or the higher
real Heisenberg groups.

We now turn to the formal proof. Note first that (3.12), (3.13), the relation
$[V_{1,i},V_{1,j}]=0$ for $i\neq j$ and (3.3) imply that

(14.3)
$d(i_{Y}i_{Y^{\prime}}i_{Z}\omega)=-i_{Z}i_{[Y,Y^{\prime}]}\omega\quad\hbox{$\forall\,\,Y,Y^{\prime}\in
V_{1,i},\,\,Z\in V_{1,j}$ with $i\neq j$.}$

######  Lemma 14.4.

Let $k\neq k^{\prime}$, $Y,Y^{\prime}\in V_{1,k}$, $Z\in V_{1,k^{\prime}}$ and
$\gamma\in\oplus_{j=1}^{n}\Lambda^{2}(V_{1,j})$. For $U\subset G$ open let
$\varphi\in C^{1}(U)$. Then the following identities hold in $U$.

(14.5) $\displaystyle d(\varphi\,i_{Y}i_{Y^{\prime}}i_{Z}\omega)=$
$\displaystyle\,(Y\varphi)\,i_{Y^{\prime}}i_{Z}\omega-\,(Y^{\prime}\varphi)i_{Y}i_{Z}\omega+(Z\varphi)\,i_{Y}i_{Y^{\prime}}\omega$
$\displaystyle+$ $\displaystyle\,\varphi\,d(i_{Y}i_{Y^{\prime}}i_{Z}\omega).$
(14.6) $\gamma\wedge
d(\varphi\,i_{Y}i_{Y^{\prime}}i_{Z}\omega)=-\gamma(Y,Y^{\prime})\,(Z\varphi)\,\omega+\varphi\,\gamma\wedge
d(i_{Y}i_{Y^{\prime}}i_{Z}\omega).$

If $X\in\Lambda^{2}V_{1,k}$ then

(14.7) $\gamma\wedge
d(\varphi\,i_{X}i_{Z}\omega)=\gamma(X)\,(Z\varphi)\,\omega+\varphi\,\gamma\wedge
d(i_{X}i_{Z}\omega).$

###### Proof.

We have

$d(\varphi\,i_{Y}i_{Y^{\prime}}i_{Z}\omega)-\varphi\,d(i_{Y}i_{Y^{\prime}}i_{Z}\omega)=d\varphi\wedge
i_{Y}i_{Y^{\prime}}i_{Z}\omega.$

Now the graded Leibniz rule (3.4) for the interior derivative gives

$\displaystyle d\varphi\wedge i_{Y}i_{Y^{\prime}}i_{Z}\omega$ $\displaystyle=$
$\displaystyle(Y\varphi)i_{Y^{\prime}}i_{Z}\omega-i_{Y}(d\varphi\wedge
i_{Y^{\prime}}i_{Z}\omega)$ $\displaystyle=$
$\displaystyle(Y\varphi)i_{Y^{\prime}}i_{Z}\omega-
i_{Y}((Y^{\prime}\varphi)i_{Z}\omega)+i_{Y}i_{Y^{\prime}}(d\varphi\wedge
i_{Z}\omega)$ $\displaystyle=$
$\displaystyle(Y\varphi)i_{Y^{\prime}}i_{Z}\omega-(Y^{\prime}\varphi)i_{Y}i_{Z}\omega+(Z\varphi)i_{Y}i_{Y^{\prime}}\omega.$

In the last step we used (3.5). This proves (14.5).

The assertion (14.6) now follows from (3.5) and the fact that
$\gamma(Y^{\prime},Z)=\gamma(Y,Z)=0$ since $Z$ lies in a different first layer
subspace than $Y$ and $Y^{\prime}$. Note that (3.3) gives the identity
$i_{Y}\circ i_{Y^{\prime}}=-i_{Y\wedge Y^{\prime}}$ and thus
$i_{Y}i_{Y^{\prime}}\gamma=-\gamma(Y,Y^{\prime})$.

Finally for $X=-Y\wedge Y^{\prime}$ (14.7) is the same as (14.6) and for
general $X$ the identity (14.7) follows by linearity. ∎

We now assume that $f:G\supset U\rightarrow G$ is a
$W^{1,q}_{\operatorname{loc}}$-mapping for some $q>\nu$ such that $D_{P}f(x)$
is an isomorphism for a.e. $x\in U$.

Recall that $\ker_{i}$ is defined in (14.2). The main point in the proof of
Theorem 11.1 for product quotients of the higher groups is to show the
following identity

(14.8)
$\displaystyle\int_{U}\alpha((D_{P}f)_{*}X)\,(Z\varphi)=0\qquad\forall\varphi\in
C_{c}^{\infty}(U),$
$\displaystyle\forall\alpha\in\oplus_{j=1}^{n}\Lambda^{2}V_{1,j},\,\,X\in\ker_{k},\,\,Z\in
V_{1,k^{\prime}}\quad\hbox{with $k^{\prime}\neq k$.}$

Then we can deduce easily that the permutation induced by $D_{P}f$ is locally
constant.

To prove (14.8) fix $\alpha$, $X$, $Z$ and $\varphi$ and consider the codegree
$3$ forms

$\beta=i_{X\wedge Z}\omega\quad\hbox{and}\quad\eta=\varphi\beta.$

It follows from (14.3) and linearity that $\beta$ is closed. Moreover
$\operatorname{wt}(\beta)=-\nu+3$. Thus Corollary 4.2 implies that

(14.9) $\int_{U}f_{P}^{*}\alpha\wedge d\eta=0.$

Since $D_{P}f$ permutes the first layer subspaces $V_{1,j}$ it follows that
$f_{P}^{*}\alpha\in\oplus_{j=1}^{n}\Lambda^{2}V_{1,j}$. Thus we get from
(14.7) and (14.3)

$f_{P}^{*}\alpha\wedge d\eta=f_{P}^{*}\alpha(X)\,(Z\varphi)\,\omega.$

This yields (14.8).

To deduce properties of the permutation $\sigma_{x}$ induced by $D_{P}f(x)$
recall that $D_{P}f(x)$ maps $\ker_{k}\subset\Lambda_{2}V_{1,k}$ to
$\ker_{\sigma_{x}(k)}$. Let

$P_{kl}(x)=\delta_{\sigma_{x}(k)l}.$

If we choose basis vectors $X_{k,m}$ of $\ker_{k}$ then there exist measurable
functions $G_{mm^{\prime}}:U\to\mathbb{R}$ such that

$(D_{P}f)_{*}X_{k,m}=\sum_{l=1}^{n}P_{kl}\sum_{m^{\prime}}G_{m^{\prime}m}X_{l,m^{\prime}}$

and the matrix $G(x)$ is invertible for a.e. $x$. Let
$\alpha_{k,m^{\prime}}\in\Lambda^{2}V_{1,k}$ be chosen such that the
restrictions to $\ker_{k}$ yield a basis of the dual space $\ker_{k}^{\prime}$
which is dual to the basis $X_{k,m}$, i.e.
$\alpha_{k,m^{\prime}}(X_{k,m})=\delta_{mm^{\prime}}$ for each $k=1,\ldots,n$.
Applying (14.8) with $X=X_{k,m}$ and $\alpha=\alpha_{l,m^{\prime}}$ we get for
each $k$ and $l$

(14.10) $\int_{U}P_{kl}G\,(Z\varphi)\,\omega=0\quad\forall\varphi\in
C_{c}^{\infty}(U),\,\,Z\in V_{1,k^{\prime}}\quad\hbox{with $k^{\prime}\neq
k$}.$

Define a function $h$ on the space of matrices by $h(G)=1$ if $\det G\neq 0$
and $h(G)=0$ else. Then $h$ is a Borel function. Since the matrix $G$ in
(14.10) is invertible a.e. in $U$, we have $h(P_{kl}G)=P_{kl}$ a.e. and it
follows from Lemma 2.11 that

(14.11) $\int_{U}P_{kl}\,(Z\varphi)\,\omega=0\quad\forall\varphi\in
C_{c}^{\infty}(U),\,\,Z\in V_{1,k^{\prime}}\quad\hbox{with $k^{\prime}\neq
k$}.$

In other words,

(14.12) $ZP_{kl}=0\quad\forall l,\quad\forall Z\in
V_{1,k^{\prime}},\,\,k^{\prime}\neq k\,\,$

in the sense of distributions. By exchanging $k$ and $k^{\prime}$ we get

(14.13) $ZP_{k^{\prime}l}=0\quad\forall l,\quad\forall Z\in
V_{1,k},\,\,k^{\prime}\neq k.\,\,$

Since $P_{kl}=1-\sum_{k^{\prime}\neq k}P_{k^{\prime}l}$ we deduce that
$ZP_{kl}=0$ for all $Z\in V_{1}$ and all $k,l$. Thus $P_{kl}$ is locally
constant and hence the permutation of the first layer subspaces induced by
$D_{P}f(x)$ is locally constant. This completes the proof.

######  Remark 14.14.

The choice of the forms $\alpha_{k,l}$ can be motivated by our guiding
principle to look for closed left-invariant forms modulo exact left-invariant
forms (see Lemma 4.5). First note by (3.19) we have
$d\alpha(X,Y)=-\alpha([X,Y])$ for $\alpha\in\Lambda^{1}V_{2}$. If we identify
two-forms on $V_{1}$ with the dual space of two-vectors by setting
$\alpha(X\wedge Y)=\alpha(X,Y)$, then we see that exterior differentiation
$d:\Lambda_{1}V_{2}\to(\Lambda_{2}V_{1})^{\prime}$ is dual map of $-L_{[]}$.
Now consider first the case that $\mathfrak{g}=\oplus_{i}\mathfrak{h}_{i}$ is
a direct sum of $n$ copies of the same (higher) Heisenberg group. A natural
space of closed left-invariant forms to detect the permutation $\sigma_{x}$
induced by the Pansu differential is the space $\oplus_{i}\Lambda^{2}V_{1,i}$.
In this case the relevant space of exact forms is given by
$d\Lambda^{1}V_{2}=\oplus_{i}d\Lambda^{1}V_{2,i}$. By the duality between $d$
and $-L_{[]}$ all elements of $d\Lambda^{1}V_{2}$ vanish on all the space
$\ker_{i}$ and the sets $\ker_{i}$ are characterized by this condition. Thus
to find a basis for the quotient space
$\oplus_{i}\Lambda^{2}V_{2,i}/d\Lambda^{1}V_{2}$ it is natural to look for
elements of $\Lambda^{2}V_{2,i}$ whose restriction to $\ker_{i}$ vanishes.
This leads to the basis $\alpha_{k,l}$ used in the proof. The fact that there
exist non-zero $X\in\ker_{i}$ is crucial for finding closed codegree forms
$\beta=i_{Z}i_{X}\omega$ which interact well with the forms $\alpha_{k,l}$.

If $\mathfrak{g}$ is a product quotient $\oplus_{i}\mathfrak{h}_{i}/K$ then it
is still true that each element of $d\Lambda^{1}V_{2}$ vanishes on each space
$\ker_{i}$. Hence the collection of cosets $\alpha_{k,l}+d\Lambda^{1}V_{2}$ is
still linearly independent in the quotient space and the proof shows that it
is sufficiently rich to detect the permutation $\sigma_{x}$.

## References

  * [Ast98] Kari Astala, _Analytic aspects of quasiconformality_ , Proceedings of the International Congress of Mathematicians, Vol. II (Berlin, 1998), no. Extra Vol. II, 1998, pp. 617–626. MR 1648110
  * [BK81] Peter Buser and Hermann Karcher, _Gromov’s almost flat manifolds_ , Astérisque, vol. 81, Société Mathématique de France, Paris, 1981\. MR 619537
  * [Bon06] Mario Bonk, _Quasiconformal geometry of fractals_ , International Congress of Mathematicians. Vol. II, Eur. Math. Soc., Zürich, 2006, pp. 1349–1373. MR 2275649
  * [BV19] Tristan Buckmaster and Vlad Vicol, _Nonuniqueness of weak solutions to the Navier-Stokes equation_ , Ann. of Math. (2) 189 (2019), no. 1, 101–144. MR 3898708
  * [Cap99] Luca Capogna, _Regularity for quasilinear equations and $1$-quasiconformal maps in Carnot groups_, Math. Ann. 313 (1999), no. 2, 263–295. MR 1679786
  * [CC06] Luca Capogna and Michael Cowling, _Conformality and $Q$-harmonicity in Carnot groups_, Duke Math. J. 135 (2006), no. 3, 455–479. MR 2272973
  * [CO15] Michael G. Cowling and Alessandro Ottazzi, _Global contact and quasiconformal mappings of Carnot groups_ , Conform. Geom. Dyn. 19 (2015), 221–239. MR 3402499
  * [CR03] Michael Cowling and Hans G. Reimann, _Quasiconformal mappings on Carnot groups: three examples_ , Harmonic analysis at Mount Holyoke (South Hadley, MA, 2001), Contemp. Math., vol. 320, Amer. Math. Soc., Providence, RI, 2003, pp. 111–118. MR 1979935
  * [Dai99] N. S. Dairbekov, _The morphism property for mappings with bounded distortion on the Heisenberg group_ , Sibirsk. Mat. Zh. 40 (1999), no. 4, 811–823, ii. MR 1721676
  * [DLS09] Camillo De Lellis and László Székelyhidi, Jr., _The Euler equations as a differential inclusion_ , Ann. of Math. (2) 170 (2009), no. 3, 1417–1436. MR 2600877
  * [DLS17] by same author, _High dimensionality and h-principle in PDE_ , Bull. Amer. Math. Soc. (N.S.) 54 (2017), no. 2, 247–282. MR 3619726
  * [DM99] Bernard Dacorogna and Paolo Marcellini, _Implicit partial differential equations_ , Progress in Nonlinear Differential Equations and their Applications, vol. 37, Birkhäuser Boston, Inc., Boston, MA, 1999. MR 1702252
  * [Fou71] Gilles Fournier, _On a problem of S. Ulam_ , Proc. Amer. Math. Soc. 29 (1971), 622. MR 278262
  * [Fox47] R. H. Fox, _On a problem of S. Ulam concerning Cartesian products_ , Fund. Math. 34 (1947), 278–287. MR 27502
  * [GQS67] Victor W. Guillemin, Daniel Quillen, and Shlomo Sternberg, _The classification of the irreducible complex algebras of infinite type_ , J. Analyse Math. 18 (1967), 107–112. MR 217130
  * [Gro86] Mikhael Gromov, _Partial differential relations_ , Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)], vol. 9, Springer-Verlag, Berlin, 1986. MR 864505
  * [Gro87] M. Gromov, _Hyperbolic groups_ , Essays in group theory, Math. Sci. Res. Inst. Publ., vol. 8, Springer, New York, 1987, pp. 75–263. MR 919829
  * [HaK00] Piotr Hajł asz and Pekka Koskela, _Sobolev met Poincaré_ , Mem. Amer. Math. Soc. 145 (2000), no. 688, x+101. MR 1683160
  * [Hei95] Juha Heinonen, _Calculus on Carnot groups_ , Fall School in Analysis (Jyväskylä, 1994), Report, vol. 68, Univ. Jyväskylä, Jyväskylä, 1995, pp. 1–31. MR 1351042
  * [Hei01] by same author, _Lectures on analysis on metric spaces_ , Universitext, Springer-Verlag, New York, 2001. MR 1800917
  * [Hei07] by same author, _Nonsmooth calculus_ , Bull. Amer. Math. Soc. (N.S.) 44 (2007), no. 2, 163–232. MR 2291675
  * [HK95] Juha Heinonen and Pekka Koskela, _Definitions of quasiconformality_ , Invent. Math. 120 (1995), no. 1, 61–79. MR 1323982
  * [HK98] by same author, _Quasiconformal maps in metric spaces with controlled geometry_ , Acta Math. 181 (1998), no. 1, 1–61. MR 1654771
  * [IM93] Tadeusz Iwaniec and Gaven Martin, _Quasiregular mappings in even dimensions_ , Acta Math. 170 (1993), no. 1, 29–81. MR 1208562
  * [Ise18] Philip Isett, _A proof of Onsager’s conjecture_ , Ann. of Math. (2) 188 (2018), no. 3, 871–963. MR 3866888
  * [J5́7] Bjarni Jónsson, _On direct decompositions of torsion-free abelian groups_ , Math. Scand. 5 (1957), 230–235. MR 108454
  * [KKL98] Michael Kapovich, Bruce Kleiner, and Bernhard Leeb, _Quasi-isometries and the de Rham decomposition_ , Topology 37 (1998), no. 6, 1193–1211. MR 1632904
  * [Kle06] Bruce Kleiner, _The asymptotic geometry of negatively curved spaces: uniformization, geometrization and rigidity_ , International Congress of Mathematicians. Vol. II, Eur. Math. Soc., Zürich, 2006, pp. 743–768. MR 2275621
  * [KMX] B. Kleiner, S. Müller, and X. Xie, In preparation.
  * [KR95] A. Korányi and H. M. Reimann, _Foundations for the theory of quasiconformal mappings on the Heisenberg group_ , Adv. Math. 111 (1995), no. 1, 1–87. MR 1317384
  * [LD17] Enrico Le Donne, _A primer on Carnot groups: homogenous groups, Carnot-Carathéodory spaces, and regularity of their isometries_ , Anal. Geom. Metr. Spaces 5 (2017), no. 1, 116–137. MR 3742567
  * [LDX16] Enrico Le Donne and Xiangdong Xie, _Rigidity of fiber-preserving quasisymmetric maps_ , Rev. Mat. Iberoam. 32 (2016), no. 4, 1407–1422. MR 3593530
  * [Lee09] Jeffrey M. Lee, _Manifolds and differential geometry_ , Graduate Studies in Mathematics, vol. 107, American Mathematical Society, Providence, RI, 2009. MR 2572292
  * [M9̈9] Stefan Müller, _Variational models for microstructure and phase transitions_ , Calculus of variations and geometric evolution problems (Cetraro, 1996), Lecture Notes in Math., vol. 1713, Springer, Berlin, 1999, pp. 85–210. MR 1731640
  * [Mau81] R. Daniel Mauldin (ed.), _The Scottish Book_ , Birkhäuser, Boston, Mass., 1981, Mathematics from the Scottish Café, Including selected papers presented at the Scottish Book Conference held at North Texas State University, Denton, Tex., May 1979. MR 666400
  * [Mic08] Peter W. Michor, _Topics in differential geometry_ , Graduate Studies in Mathematics, vol. 93, American Mathematical Society, Providence, RI, 2008. MR 2428390
  * [MM95] G. A. Margulis and G. D. Mostow, _The differential of a quasi-conformal mapping of a Carnot-Carathéodory space_ , Geom. Funct. Anal. 5 (1995), no. 2, 402–433. MR 1334873
  * [Mur81] François Murat, _Compacité par compensation: condition nécessaire et suffisante de continuité faible sous une hypothèse de rang constant_ , Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 8 (1981), no. 1, 69–102. MR 616901
  * [Mv03] S. Müller and V. Šverák, _Convex integration for Lipschitz mappings and counterexamples to regularity_ , Ann. of Math. (2) 157 (2003), no. 3, 715–742. MR 1983780
  * [Nas54] John Nash, _$C^{1}$ isometric imbeddings_, Ann. of Math. (2) 60 (1954), 383–396. MR 65993
  * [OW11] Alessandro Ottazzi and Ben Warhurst, _Contact and 1-quasiconformal maps on Carnot groups_ , J. Lie Theory 21 (2011), no. 4, 787–811. MR 2917692
  * [Pan83] Pierre Pansu, _Croissance des boules et des géodésiques fermées dans les nilvariétés_ , Ergodic Theory Dynam. Systems 3 (1983), no. 3, 415–445. MR 741395
  * [Pan89] by same author, _Métriques de Carnot-Carathéodory et quasiisométries des espaces symétriques de rang un_ , Ann. of Math. (2) 129 (1989), no. 1, 1–60. MR 979599
  * [Rei01] Hans Martin Reimann, _Rigidity of $H$-type groups_, Math. Z. 237 (2001), no. 4, 697–725. MR 1854087
  * [Res89] Yu. G. Reshetnyak, _Space mappings with bounded distortion_ , Translations of Mathematical Monographs, vol. 73, American Mathematical Society, Providence, RI, 1989, Translated from the Russian by H. H. McFaden. MR 994644
  * [RR00] H. M. Reimann and F. Ricci, _The complexified Heisenberg group_ , Proceedings on Analysis and Geometry (Russian) (Novosibirsk Akademgorodok, 1999), Izdat. Ross. Akad. Nauk Sib. Otd. Inst. Mat., Novosibirsk, 2000, pp. 465–480. MR 1847532
  * [Sch93] Vladimir Scheffer, _An inviscid flow with compact support in space-time_ , J. Geom. Anal. 3 (1993), no. 4, 343–401. MR 1231007
  * [Sha04] Yehuda Shalom, _Harmonic analysis, cohomology, and the large-scale geometry of amenable groups_ , Acta Math. 192 (2004), no. 2, 119–185. MR 2096453
  * [Spe69] D. C. Spencer, _Overdetermined systems of linear partial differential equations_ , Bull. Amer. Math. Soc. 75 (1969), 179–239. MR 242200
  * [SX12] Nageswari Shanmugalingam and Xiangdong Xie, _A rigidity property of some negatively curved solvable Lie groups_ , Comment. Math. Helv. 87 (2012), no. 4, 805–823. MR 2984572
  * [Tar79] L. Tartar, _Compensated compactness and applications to partial differential equations_ , Nonlinear analysis and mechanics: Heriot-Watt Symposium, Vol. IV, Res. Notes in Math., vol. 39, Pitman, Boston, Mass.-London, 1979, pp. 136–212. MR 584398
  * [Vod07] S. K. Vodopyanov, _Foundations of the theory of mappings with bounded distortion on Carnot groups_ , The interaction of analysis and geometry, Contemp. Math., vol. 424, Amer. Math. Soc., Providence, RI, 2007, pp. 303–344. MR 2316342
  * [War07] Ben Warhurst, _Tanaka prolongation of free Lie algebras_ , Geom. Dedicata 130 (2007), 59–69. MR 2365778
  * [Whi40] J. H. C. Whitehead, _On the homotopy type of manifolds_ , Ann. of Math. (2) 41 (1940), 825–832. MR 2546
  * [Xie13a] X. Xie, _Quasiconformal maps on non rigid carnot groups_ , arxiv1308.3031, 2013\.
  * [Xie13b] Xiangdong Xie, _Quasisymmetric homeomorphisms on reducible Carnot groups_ , Pacific J. Math. 265 (2013), no. 1, 113–122. MR 3095115
  * [Xie15] by same author, _Quasi-conformal maps on model filiform groups_ , Michigan Math. J. 64 (2015), no. 1, 169–202. MR 3326585


# Chaos and Entanglement in Non-Markovian Optomechanical Systems

Pengju Chen  Nan Yang  Austen Couvertier  Quanzhen Ding  Rupak Chatterjee
Ting Yu<EMAIL_ADDRESS>Center for Quantum Science and Engineering
Department of Physics, Stevens Institute of Technology
Hoboken, New Jersey 07030, USA.

###### Abstract

We study the chaotic motion of an optomechanical system coupled to a non-
Markovian environment. We show that the environmental memory time can
significantly affect chaos in an enhancing way. In addition to classical
chaotic motion, the quantum entanglement in the presence of chaos is
investigated. It is found that both the environmental memory and chaos can
lift up bipartite entanglement in a non-linear optomechanical system. These
observations may help expand our understanding of the transition from
classical to quantum dynamics.

††preprint: APS/123-QED

##  I Introduction

Optomechanics has provided a powerful tool to study both classical and quantum
dynamics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. By using non-
Markovian style “baths” interacting with such systems of interest, one is able
to fine tune a memory time parameter found in these non-Markovian reservoirs
and observe their effects on the system in question. In particular, we focus
on the effects of this memory time parameter on the signature of classical
chaos and quantum entanglement within our optomechanical system. While many
have studied Markovian open systems using the widely used ansatz of the Born-
Markov approximation, which assumes a memory-less environment and a weak
coupling between system and environment [16, 17, 18, 19, 20], the theory on
non-Markovian quantum open systems is less widely used but more realistic for
many true experimental systems. With the development of non-Markovian quantum-
state diffusion (NMQSD) [21, 22, 23, 24, 25, 26, 27, 28, 29], we are able to
study the non-Markovian effects of open optomechanical systems where the bath
of the system has a ‘memory’ such that it may feed information back to the
system in a non-trivial way.

Classical chaos has exhibit fascinating features in connections with various
topics in math and physics, such as fractal dimension [30], dynamical systems
[31], and of course, quantum systems [32], while showing numerous applications
[33, 34]. Many attempts have been made to study chaotic dynamics of quantum
systems including semi-classical orbit methods [35, 36, 37, 38] and
statistical and spectral chaos using random matrix methods [39, 40]. The route
to chaos generated by a optomechanical system within a Markov environment has
been studied in [6]. This naturally leads to the question of what would happen
to such a system embedded in non-Markovian environments. This is the first
part of our investigation below.

As entangled states may be generated between these coupled optomechanical
quantum systems, it is expected that the qualitative nature of the dynamics
can affect entanglement [41]. This will be the second part of our study
investigating how signatures of classical chaos can affect quantum
entanglement. Different from a linearized system, our system has a non-linear
Hamiltonian resulting in non-Gaussian states. Recent methods [42, 43] are
employed to determine the negativity of our entangled state and a time
averaged entanglement measure is used to provide a comparison with the
prototypical quantitative measure of classical chaos, the maximal Lyapunov
exponent ($\bm{LE}$).

Additionally, it has been shown that non-Markovian environment plays a pivotal
role in entanglement dynamics [44, 45, 46], which has shown usefulness in
preverving and enhancing entanglement [47, 48] in optomechanial systems. It
would be interesting to see whether it’s still the case with the presence of
classical chaos.

The paper is organized as follows. Section II will give a detailed description
of our optomechanical system coupled to a non-Markovian environment. Section
III and IV will provide results of our numerical simulations, with section III
focusing on the non-Markovian reservoir’s influence on classical chaos.
Section IV explores the relation between entanglement and chaos within the
optomechanical non-Markovian open system. Finally, a conclusion of our results
is given.

Figure 1: Schematic diagram of the optomechanical system, where the light
field is considered to be in non-Markovian environment. The mechanical bath is
set to be Markov to account for mechanical loss.

##  II Optomechanical model coupled to Non-Markovian environment

In this section, we introduce the optomechanical system and the non-Markovian
environment that will be used to investigate chaos and entanglement. The
optomechanical system is where a mechanical mode couples with an optical mode,
as exhibited in Figure 1.

Assuming $\hbar=1$, the system is described by the Hamiltonian [1] given by

$H_{s}=\big{[}-\Delta+g_{0}(b+b^{\dagger})\big{]}a^{\dagger}a+\Omega
b^{\dagger}b+\alpha_{\textit{{L}}}(a^{\dagger}+a),$ (1)

where $a$ and $b$ are the annihilation operators of the optical mode and
mechanical mode which satisfy the commutation relations $[a,a^{\dagger}]=1$
and $[b,b^{\dagger}]=1$. $\Delta$ is the detuning parameter and $\Omega$ is
the cantilever frequency. The detuning is defined as
$\Delta\equiv\omega_{L}-\omega_{\mathrm{cav}}$, where $\omega_{L}$ is the
pumping laser frequency and $\omega_{\mathrm{cav}}$ is the resonance frequency
of the cavity mode $a$. $\alpha_{L}$ is the pumping laser amplitude and
$g_{0}$ is the coupling constant between the two modes.

Since the optomechanical system is modelled as an open system, the cavity has
radiative loss and cantilever has mechanical damping with respect to their
local environments. The two baths are Bosonic, with the optical bath is set to
be non-Markovian, which is the main focus of our study and can be adjusted
experimentally [49]. While the mechanical bath is considered Markov to account
for mechanical loss. The two baths are separated.

Using the first order approximation of the NMQSD equation [22], the master
equation takes the form (details can be seen in Appendix A)

$\dot{\rho}=-i[H_{s},\rho]+\Gamma
D\big{[}b,\rho\big{]}+\bigg{\\{}f_{0}(t)[a,\rho
a^{\dagger}]+if_{1}(t)[a^{\dagger},[H_{s},a]\rho]+f_{2}(t)[a^{\dagger},[a^{\dagger},a]a\rho]+H.C.\bigg{\\}},$
(2)

where $H.C.$ stands for Hermitian conjugate. The Lindblad term $\Gamma
D\big{[}b,\rho\big{]}=\Gamma\big{\\{}b\rho
b^{\dagger}-\frac{1}{2}(b^{\dagger}b\rho+\rho b^{\dagger}b)\big{\\}}$
represents the Markov mechanical bath. $f_{0}$, $f_{1}$ and $f_{2}$ are time-
dependent coefficients which are given by

$\begin{split}&f_{0}(t)=\int_{0}^{t}\alpha(t,s)\;ds,\\\
&f_{1}(t)=\int_{0}^{t}\alpha(t,s)(t-s)\;ds,\\\
&f_{2}(t)=\int_{0}^{t}\int_{0}^{s}\alpha(t,s)\alpha(s,u)(t-s)\;du\;ds.\end{split}$
(3)

$\alpha(t,s)$ is the Ornstein-Uhlenbeck (O-U) correlation function

$\alpha(t,s)=\frac{\kappa\gamma}{2}e^{-\gamma|t-s|},$ (4)

where $1/\gamma$ is the memory time, and $\kappa/\Omega=1$ is the optical
damping rate.

We express the system parameters in units of $\Omega$ and introduce the
dimensionless time parameter $\tau=\Omega t$. Additionally, we introduce two
dimensionless variables [50, 51, 6]

$\sigma=g_{0}/\kappa,\ P=\frac{8\alpha_{L}^{2}g_{0}^{2}}{\Omega^{4}}.$ (5)

The pumping parameter $P$ gives the strength of the laser input of the cavity.
The quantum-classical scaling parameter $\sigma$ is set to be
$\sigma=x_{\mathrm{zpt}}/x_{\mathrm{res}}=g_{0}/\kappa=0.1$, where
$x_{\mathrm{zpt}}=\sqrt{\hbar/(2m\Omega)}$ are the zero-point fluctuations of
the cantilever (with mass $m$). Furthermore, $x_{\mathrm{res}}$ is the
resonance width of the cavity, and $g_{0}$ is the optomechanical single-photon
coupling strength [51, 6]. Since $x_{\mathrm{res}}$ is a classical quantity,
$\sigma$ shows how close the quantum dynamics of the optomechanical system is
to the classical limit, where $\sigma\to 0$ provides the classical limit.

We use the re-scaled creation and annihilation operators $\langle
a_{1}\rangle=[\Omega/(2\alpha_{L})]\langle a\rangle$, $\langle
b_{1}\rangle=(g_{0}/\Omega)\langle b\rangle$. Applying the trace expectation
$\langle\dot{A}\rangle=tr(A\dot{\rho})$, the complete set of dynamical
equations can be seen in Appendix B Eq. (27), here we give the first two
equations

$\begin{split}\frac{d}{d\tau}\langle
a_{1}\rangle=&-i(1+f_{1})\bigg{\\{}\langle a_{1}b_{1}\rangle+\langle
a_{1}b_{1}^{\dagger}\rangle-\frac{\Delta}{\Omega}\langle
a_{1}\rangle+\frac{1}{2}\bigg{\\}}-\frac{f_{0}^{\ast}+f_{2}}{\Omega}\langle
a_{1}\rangle,\\\ \frac{d}{d\tau}\langle
b_{1}\rangle=&-i\bigg{\\{}\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle+\langle
b_{1}\rangle\bigg{\\}}-\frac{\Gamma}{2\Omega}\langle b_{1}\rangle.\end{split}$
(6)

Eq. (27) forms a set of closed dynamical equations which will be used to
investigate both chaos and entanglement of our system. Semi-classical (SC)
approximation is made since the coupling $g_{0}$ is weak. Since the
Hamiltonian is non-linear, an exact set of equations will have too many terms.
Therefore, we use a SC approximation such that terms higher than second order
are split.

(a) $\gamma=10$

(b) $\gamma=2$

(c) $\gamma=0.8$

Figure 2: Phase diagrams with different correlation frequencies $\gamma$ at
$\Delta=-0.70$, $P=1.25$. The coordinates $(Z_{1},Z_{2},Z_{3},Z_{4})$ are the
real and imaginary parts of $\langle a_{1}\rangle$ and $\langle b_{1}\rangle$.
The 4th coordinate $Z_{4}$ is represented by scaled colours. As memory time
increases ($\gamma$ decreases), the dynamics goes from regular to chaotic.
Fig. 2(c) exhibits the emerging of chaos as $\gamma$ is lowered to 0.8.

(a) $\gamma=10$

(b) $\gamma=2$

Figure 3: Bifurcation diagrams and corresponding Lyapunov exponents
($\bm{LE}$), with the laser pumping parameter $P=1.37$ and detuning varying
from -1.2 to -0.4. In Fig. 3(a), the chaotic region is bounded within a small
segment $\Delta\in[-1.03,-0.92]$. In Fig. 3(b), the chaotic region expanding
to $\Delta\in[-1.13,-0.99]$ and a new chaotic region $\Delta\in[-0.83,-0.52]$
emerges, with some inter-adjacent regular regions inside.

(a) $\gamma=10$

(b) $\gamma=1$

Figure 4: Pictures of chaotic regions of the optomechanical systems with
different memory times plotted in the $P$-$\Delta$ plane. The colour scale
shows the value of maximal Lyapunov exponent ($\bm{LE}$) of every data point.
Comparing Fig. 4(a) and 4(b), the chaotic area expands as the memory time gets
longer ($\gamma$ decreases).

##  III Chaos and Non-Markovian Environments

The mechanism of how a non-Markovian environment affects a chaotic system is
the focus of this section. In general, chaotic systems tend to be very
sensitive to parameter changes. The simulation methods used below allow one to
manually adjust such a parameter, the non-Markovian memory time, which can
have a delayed effect on the dissipative process of our open system. Such an
effect may significantly influence chaos generation and therefore, we examine
various memory times and their effects of chaos generation.

Our simulation is mainly based on observing the optomechanical system while
changing the memory time of the optical bath represented by the parameter
$\gamma$ (the inverse of memory time). Furthermore, we vary the system
parameters $\Delta$ and $P$ to get a comprehensive picture of the chaos
distribution of our system. The initial states are set to be the vacuum
states. We use the maximal Lyapunov exponent ($\bm{LE}$) as the indicator of
chaos, which is calculated using the Wolf’s method of phase reconstruction
[52]. If the $\bm{LE}$ is positive, the system is considered to be chaotic.

By taking the real and imaginary parts of the time series of $\langle
a_{1}\rangle$ and $\langle b_{1}\rangle$ from Eq. (27), we can generate the
four-dimensional phase diagrams based on different memory times (the inverse
of $\gamma$). Figure (2) below illustrates the system dynamics under different
memory times for a fixed point ($P=1.25$ and $\Delta=-0.70$). The increase of
$\gamma$ results in the system changing from regular to chaotic, which serves
as a key example of memory induced chaos.

Next, we plot the bifurcation diagrams which mark the stable points of the
cantilever oscillation. We fix $P=1.37$. At $\gamma=10$, the period-doubling
bifurcation takes place, as is exhibited in Fig. (3(a)), where the chaotic
region is bounded within a small segment $\Delta\in[-1.03,-0.92]$. For the
$\gamma=2$ case, Fig. (3(b)) shows the chaotic region expanding to
$\Delta\in[-1.13,-0.99]$ while a new chaotic region at
$\Delta\in[-0.83,-0.52]$ emerges, with some inter-adjacent regular regions
inside. The comparison between Fig. (3(a)) and (3(b)) clearly indicates that
the longer memory time can, not only expand chaotic regions, but also induce
chaos from the previously non-chaotic areas. Note that $\gamma=2$ is not far
away from the Markov limit, but it still causes significant change to chaotic
dynamics, further indicating that the chaos is very sensitive to parameter
changes, including the memory time.

To summarize the appearances of chaos, the $\bm{LE}$ of every data point is
plotted to form global chaos landscapes, as is shown in Fig. (4(a)) for
$\gamma=10$ and (4(b)) for $\gamma=1$. The parameter ranges are set to be
$P\in[0.8,1.6]$ and $\Delta\in[-1.4,-0.4]$. As the pumping increases, we see
more chaotic motion [6, 51, 50]. In both cases, the chaotic regions have some
fine structures of inter-adjacent regular regions. For our model, we have
shown that the non-Markovian condition expands the chaotic regions while
lowering the pumping bar for chaos generation. For the the case of $\gamma=10$
(Fig. 4(a)), the chaos emerging level is $P=1.37$, while for $\gamma=1$ (Fig.
4(b)), the level is lowered to $P=1.08$.

Comparing Fig. (4(a)) and (4(b)), we have seen that the increase of memory
time can expand chaotic areas in the parameter plane and decrease the pumping
threshold for chaos generation. While the amplitude damping is the major
dissipating channel, the environmental memory can slow down the dissipation of
energy. The dissipative dynamics will undergo temporal revivals due to the
memory effect. This results in the emerging of chaos that requires less
pumping energy, thus lowering the pumping threshold while expanding chaotic
regions across the map. A similar phenomena has been discussed in [48] where
the memory enhanced entanglement generation was also due to the back-flow in
the dissipation channel in the form of information. Since energy carries
information, our model can be understood in a similar fashion.

(a) $\gamma=10$

(b) $\gamma=1$

Figure 5: The entanglement strength $\bm{En}$ on the $P$-$\Delta$ plane. The
colour scale measures the value of $\bm{En}$. The entanglement shows similar
fine structures as the chaos.

(a) $P=1.3$

(b) $P=1.4$

Figure 6: Bifurcation diagram, maximal Lyapunov exponent ($\bm{LE}$),
comparing with average negativity of the partial transposition ($\bm{En}$),
when $\gamma=10$.  Figure 7: Phase diagram and the corresponding long term
dynamical evolution of the NPT $\mathcal{N}$, with $P=1.4$. The first roll
shows $\Delta=-1.1$ (non-chaotic) and the second roll $\Delta=-1.0$ (chaotic).
The dynamical evolution of $\mathcal{N}$ appears solid simply because it is
oscillating very quickly.

(a) $P=1.3$

(b) $P=1.4$

Figure 8: Bifurcation diagram, maximal Lyapunov exponent ($\bm{LE}$),
comparing with average negativity of the partial transposition ($\bm{En}$),
when $\gamma=1$.

##  IV Entanglement and Chaos

We now turn to the quantum aspects of our optomechanical system to investigate
entanglement under the presence of chaos and memory time. The optomechanical
system we are studying is a weakly coupled bipartite system. There have been
several studies focusing on chaos and its effects on entanglement [53, 54, 55,
41, 56, 57]. One previous study looks at the effects of environment memory on
entanglement [48].

Since the non-linear Hamiltonian of our system (1) causes our physical states
to be non-Gaussian, the conventional ways of calculating inseparability [58,
59, 60] would cause imaginary terms that cannot be ignored. Therefore, a non-
conventional method is needed to calculate the entanglement of the two modes
of our system. Here, we will use a hierarchy of necessary and sufficient
conditions for the negativity of the partial transposition (NPT) in terms of
observable moments [42]. Though this leads only to sufficient conditions for
entanglement, it can be applied to a variety of quantum states (including non-
Gaussian states).

###  IV.1 Entanglement Measurement

A brief introduction of Shchukin and Vogel’s (SV) criterion is provided here
[42]. It is known that a bipartite quantum state is entangled if NPT is
achieved [43, 61]. We start with a separable bipartite state $\hat{\rho}$. The
condition for its separability is that its partial transposed density operator
$\hat{\rho}^{\text{PT}}$ can be written in the following form (assuming that
we are transposing the second Hilbert space state, which won’t affect our
result)

$\hat{\rho}^{\text{PT}}=\sum^{\infty}_{n=0}p_{n}\hat{\rho}^{(n)}_{1}\otimes\hat{\rho}^{(n)\text{T}}_{2}.$
(7)

Based on this property of separable states, we can employ NPT as a sufficient
condition for inseparability, which is referred to as the Peres-Horodecki
condition [61, 62, 63]. The separability is determined through a matrix of
moments of the partial transposition, which is given by

$M_{pqrs,nmkl}=\langle a^{\dagger q}a^{p}a^{\dagger n}a^{m}b^{\dagger
s}b^{r}b^{\dagger k}b^{l}\rangle^{\text{PT}}.$ (8)

Shchukin and Vogel [42] introduced a way of ordering the moments
$M_{pqrs,nmkl}$ as follows

$\begin{split}&1,\ \langle a\rangle,\ \langle a^{\dagger}\rangle,\ \langle
b\rangle,\ \langle b^{\dagger}\rangle,\ \langle a^{2}\rangle,\ \langle
a^{\dagger}a\rangle,\ \langle a^{\dagger 2}\rangle,\ \langle ab\rangle,\\\
&\langle a^{\dagger}b\rangle,\ \langle b^{2}\rangle,\ \langle
ab^{\dagger}\rangle,\ \langle a^{\dagger}b^{\dagger}\rangle,\ \langle
b^{\dagger}b\rangle,\ \langle b^{\dagger 2}\rangle,\ ...\end{split}$ (9)

The moments $M_{pqrs,nmkl}$ in Eq. (8) can be expressed in terms of the
moments of the original state

$\langle a^{\dagger q}a^{p}a^{\dagger n}a^{m}b^{\dagger s}b^{r}b^{\dagger
k}b^{l}\rangle^{\text{PT}}=\langle a^{\dagger q}a^{p}a^{\dagger
n}a^{m}b^{\dagger l}b^{k}b^{\dagger r}b^{s}\rangle.$ (10)

For our quantum state, it is convenient to use a higher-order test involving a
small number of minors, which was introduced in [43]. The matrix determinant
is given as follows (using the re-scaled $a_{1}$ and $b_{1}$)

$D_{HO}=\begin{vmatrix}1&\langle a_{1}b_{1}^{\dagger}\rangle\\\ \langle
a_{1}^{\dagger}b_{1}\rangle&\langle
a_{1}^{\dagger}a_{1}b_{1}^{\dagger}b_{1}\rangle\end{vmatrix}.$ (11)

If there exist a negative determinant

$D_{HO}<0,$ (12)

then the NPT has been demonstrated, which provides a sufficient condition for
entanglement. Conveniently, the set of mean values are already provided by the
semi-classical equations of motion (27), and the SC approximation is used such
that $\langle a_{1}^{\dagger}a_{1}b_{1}^{\dagger}b_{1}\rangle=\langle
a_{1}^{\dagger}a_{1}\rangle\langle b_{1}^{\dagger}b_{1}\rangle$.

We define $\mathcal{N}$ to denote the negativity of the partial transposition
(NPT)

$\mathcal{N}=\mathrm{max}[0,-D_{HO}].$ (13)

A positive value for $\mathcal{N}$ indicates that the states are entangled. We
then define $\bm{En}$ as the long-time average of $\mathcal{N}(t)$ [56]. It
provides an overview of the entanglement intensity as is given by

$\bm{En}=\frac{1}{T}\int_{0}^{T}\mathcal{N}(t)dt,$ (14)

where $T\to\infty$. Since $\bm{En}$ is only determined by parameters $P$,
$\Delta$ and $\gamma$, we can calculate the averge NPT ($\bm{En}$) of every
data point analogous to the $\bm{LE}$ for chaos.

###  IV.2 Simulation Results

The following simulations show the quantum entanglement under the influence of
chaotic motion and non-Markovian environment. We consider two memory
parameters, namely the close to Markov case and the non-Markovian case.

First, we focus on the close to Markov case, with $\gamma=10$. Overall, the
graph (Fig. 5(a)) of average NPT ($\bm{En}$) on the parameter plane reflects
almost the same fine structure as the chaos graph (Fig. 4(a)). It is worth
noting that the very small $\bm{LE}$ spikes in Fig. (4(a)) coincide with sharp
$\bm{En}$ spikes. To provide a closer look, we show the bifurcation,
$\bm{LE}$, and $\bm{En}$ diagrams, fixing $\gamma=10$. With $P=1.3$, Fig.
(6(a)) shows two very small $\bm{LE}$ spikes coincide with two huge $\bm{En}$
spikes. Fig. (6(b)) demonstrates that the chaotic area
$\Delta\in[-1.066,-0.94]$ corresponding to an area of entanglement increase.
As shown in Fig. (7), the chaos can enhance the bipartite entanglement for the
non-linear optomechanical system, our results are in tune with the
observations for a quantum kicked top system [56]. In summary, the bipartite
entanglement generation can be increased by classical chaotic dynamics.
However, it should be noted that the $\bm{En}$ level (as shown in Fig. 4(a))
can be lower than some non-chaotic areas on the parameter plane meaning that
the enhancement effect of entanglement by chaos is localized (on the parameter
plane).

Furthermore, in the non-Markovian $\gamma=1$ case (Fig. 5(b)), we can still
see some fine structures related to the chaos graph (Fig. 4(b)), where the
$\bm{En}$ level around chaotic areas is larger than the $\gamma=10$ case (Fig.
4(a)). In Fig. (8(b)), we can observe the slight enhancement of entanglement
by chaos. Yet, that effect dies down at $\Delta>-0.78$ while in Fig. (8(a)),
it terminates at $\Delta>-0.74$. Overall, the $\bm{En}$ level is higher than
the $\gamma=10$ case and the distinction of $\bm{En}$ between chaotic and non-
chaotic areas is much less, which is similar to the results demonstrated in a
linear optomechanical system [48], whose explanation is that the memory effect
creates information back flow such that the dissipative dynamics will
experience temporal revivals. The dissipation and the back flow from the
environment may reach a new balance point so that the steady state
entanglement ($t\to\infty$) may be dependent on its memory. The result of this
mechanics is the enhancement of entanglement by the environmental memory time,
in our non-linear optomechanical system. Though this effect is still
localized, as one cannot tell whether all areas on Fig. (5(b)) have higher
$\bm{En}$ than Fig. (5(a)). In summary, both the environmental memory and the
chaotic dynamics can significantly affect entanglement generation.

##  V Conclusion

We have analyzed classical chaos in a non-Markovian optomechanical system and
measured the corresponding quantum entanglement. We have shown that
environmental memory can increase chaos generation in the optomechanical
system. For the non-linear optomechanical system, it is shown that the memory
effect causes the energy back flow resulting in the chaos region expanding
while lowering the pumping bar for chaos generation. As for the entanglement,
we observed two effects. The first one is that bipartite entanglement is
sensitive to chaos and can be enhanced by it. The second effect is that the
environmental memory time can increase bipartite entanglement generation. In
conclusion, in the non-linear optomechanical system, non-Markovian memory time
can be an enhancing factor for classical chaos generation, while chaos and
environmental memory both can increase entanglement generation.

## Acknowledgement

This work is partially supported by the ART020-Quantum Technologies Project.
We thank Dr. Mengdi Sun, Kenneth Mui, and Yifan Shi for helpful discussion on
optomechanics.

##  Appendix A The non-Markovian quantum state diffusion (NMQSD) equation

The following section provides a detailed derivation of the non-Markovian
master equation.

Both environments are set to be at zero temperature, and separate from each
other. Since the Markov environment is well-understood and its formation can
be added to the master equation easily, we can focus on the optical bath and
use the non-Markovian quantum state diffusion (QSD) equation to derive the
master equation.

The non-Markovian bath for the optical mode can be described by a set of
harmonic oscillators [48]

$H_{B}=\sum_{j}\omega_{j}c^{\dagger}_{j}c_{j},$ (15)

where $c_{j}$ and $c_{j}^{\dagger}$ are annihilation and creation operators
satisfying the commutation relation $[c_{i},c_{j}^{\dagger}]=\delta_{ij}$. In
general, the interaction Hamiltonian between the system and its Bosonic bath
is described by

$H_{I}=\sum_{j}g_{j}(Lc^{\dagger}_{j}+L^{\dagger}c_{j}),$ (16)

where $L=a$ is the Lindblad operator representing the optical damping and
$g_{j}$ are the system-bath coupling intensities. The interaction Hamiltonian
is in the form of a rotating wave approximation which assumes the condition of
weak coupling strength ($g_{j}\ll\Omega$).

Assuming the optomechanical system and the environment are not correlated
initially, the state of the optomechanical system connecting to a non-
Markovian bath can be represented by the non-Markovian quantum state diffusion
(NMQSD) equation [23, 21] governing the time evolution of the a stochastic
pure state $|\psi_{t}(z^{\ast})\rangle$, which is given by

$\partial_{t}|\psi_{t}(z^{\ast})\rangle=\bigg{[}-iH_{s}+az^{\ast}_{t}-a^{\dagger}\overline{O}(t,z^{\ast})\bigg{]}|\psi_{t}(z^{\ast})\rangle,$
(17)

where $O(t,s,z^{\ast})\psi_{t}\equiv\frac{\delta\psi_{t}}{\delta z_{s}}$ and
$\overline{O}(t,z^{\ast})\equiv\int^{t}_{0}ds\alpha(t,s)O(t,s,z^{\ast})$ with
the initial condition $O(t,s=t,z^{\ast})=a$.
$\alpha(t,s)=\frac{\kappa\gamma}{2}e^{-\gamma|t-s|}$ is the Ornstein-Uhlenbeck
(O-U) correlation function.
$z_{t}^{\ast}=-i\sum_{j}g_{j}z^{\ast}_{j}e^{i\omega_{j}t}$ is a colored
complex Gaussian process satisfying

$M[z_{t}^{\ast}z_{s}]=\alpha(t,s).\ M[z_{t}z_{s}]=0,$ (18)

where $M[\cdot]\equiv\int\frac{dz^{2}}{\pi}e^{-|z|^{2}}[\cdot]$ denotes the
ensemble average over the noise $z_{t}$.

To solve (17), one needs to find the operator $O(t,s,z^{\ast})$. Under the
condition that the memory time is not too long, we consider the expansion of
$O(t,s,z^{\ast})$ in powers of $(t-s)$ [22]

$O(t,s,z^{\ast})=O(s,s,z^{\ast})+\frac{\partial O(t,s,z^{\ast})}{\partial
t}\bigg{|}_{t=s}(t-s)+\cdots,$ (19)

which is a systematic expansion of the non-Markovian QSD. The zeroth-order
term corresponds to the standard Markov QSD when $1/\gamma\to 0$. We choose
the first-order approximation of the operator $O(t,s,z^{\ast})$ since the
first-order term is the most important correction to Markovian dynamics given
that the memory time is assumed to be short.

At the time point $t=s$, the expression of the operator $O(t,s,z^{\ast})$ is
given by [22]

$\displaystyle O(s,s,z^{\ast})$ $\displaystyle=a,$ (20)
$\displaystyle\frac{\partial O(t,s,z^{\ast})}{\partial t}\bigg{|}_{t=s}$
$\displaystyle=-i[H_{s},a]-\int^{s}_{0}\alpha(s,u)du[a^{\dagger},a]a,$ (21)

where $H_{s}$ is the system Hamiltonian and $a$ is the Lindblad operator. At
this time point, the operator $\overline{O}$ becomes

$\overline{O}=f_{0}(t)a-f_{1}(t)i[H_{s},a]-f_{2}(t)[a^{\dagger},a]a,$ (22)

where

$\begin{split}&f_{0}(t)=\int_{0}^{t}\alpha(t,s)\;ds,\\\
&f_{1}(t)=\int_{0}^{t}\alpha(t,s)(t-s)\;ds,\\\
&f_{2}(t)=\int_{0}^{t}\int_{0}^{s}\alpha(t,s)\alpha(s,u)(t-s)\;du\;ds.\end{split}$
(23)

One can turn (17) into a master equation by taking the ensemble mean over the
noise $z_{t}$ and introducing the reduced density matrix
$\rho_{t}=M[|\psi_{t}(z^{\ast})\rangle\langle\psi_{t}(z^{\ast})|]$. When the
environment is not far away from Markov, the dependence of the operator
$\overline{O}(t,z^{\ast})$ on the noise $z_{t}$ is negligible. Under this
approximation, the master equation [22] takes the form

$\frac{d}{dt}\rho_{t}=-i[H_{s},\rho_{t}]+[a,\rho_{t}\overline{O}(t)^{\dagger}]-[a^{\dagger},\overline{O}(t)\rho_{t}].$
(24)

For further details about the derivation of the master equation, see [24, 28].

Since the mechanical bath is separate and Markov, it can be added to the
master equation by simply including the appropriate Lindblad term which is
given by

$\Gamma D\big{[}b,\rho\big{]}=\Gamma\big{\\{}b\rho
b^{\dagger}-\frac{1}{2}(b^{\dagger}b\rho+\rho b^{\dagger}b)\big{\\}},$ (25)

where $\Gamma/\Omega=10^{-3}$ is the mechanical damping.

Using the first order approximation of $\overline{O}(t)$ (22), the master
equation takes the final form

$\dot{\rho}=-i[H_{s},\rho]+\Gamma
D\big{[}b,\rho\big{]}+\bigg{\\{}f_{0}(t)[a,\rho
a^{\dagger}]+if_{1}(t)[a^{\dagger},[H_{s},a]\rho]+f_{2}(t)[a^{\dagger},[a^{\dagger},a]a\rho]+H.C.\bigg{\\}},$
(26)

where $H.C.$ stands for Hermitian conjugate.

##  Appendix B The complete set of dynamical equations

$\begin{split}\frac{d}{d\tau}\langle
a_{1}\rangle=&-i(1+f_{1})\bigg{\\{}\langle a_{1}b_{1}\rangle+\langle
a_{1}b_{1}^{\dagger}\rangle-\frac{\Delta}{\Omega}\langle
a_{1}\rangle+\frac{1}{2}\bigg{\\}}-\frac{f_{0}^{\ast}+f_{2}}{\Omega}\langle
a_{1}\rangle,\\\ \frac{d}{d\tau}\langle
b_{1}\rangle=&-i\bigg{\\{}\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle+\langle
b_{1}\rangle\bigg{\\}}-\frac{\Gamma}{2\Omega}\langle b_{1}\rangle,\\\
\frac{d}{d\tau}\langle a_{1}^{\dagger}a_{1}\rangle=&-\frac{i}{2}(\langle
a_{1}^{\dagger}\rangle-\langle
a_{1}\rangle)-\frac{1}{\Omega}(f_{0}+f_{0}^{\ast}+f_{2}+f_{2}^{\ast})\langle
a_{1}^{\dagger}a_{1}\rangle\\\
&+i(f_{1}^{\ast}-f_{1})(-\frac{\Delta}{\Omega}+\langle
b_{1}^{\dagger}\rangle+\langle b_{1}\rangle)\langle
a_{1}^{\dagger}a_{1}\rangle+\frac{i}{2}(f_{1}^{\ast}\langle a_{1}\rangle-
f_{1}\langle a_{1}^{\dagger}\rangle),\\\ \frac{d}{d\tau}\langle
b_{1}^{\dagger}b_{1}\rangle=&-i\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle(\langle b^{\dagger}_{1}\rangle-\langle
b_{1}\rangle)-\frac{\Gamma}{\Omega}\langle b_{1}^{\dagger}b_{1}\rangle,\\\
\frac{d}{d\tau}\langle
a_{1}^{2}\rangle=&-2i(1+f_{1})\bigg{\\{}(-\frac{\Delta}{\Omega}+\langle
b_{1}^{\dagger}\rangle+\langle b_{1}\rangle)\langle
a_{1}^{2}\rangle+\frac{\langle
a_{1}\rangle}{2}\bigg{\\}}-\frac{2}{\Omega}(f_{0}^{\ast}+f_{2})\langle
a_{1}^{2}\rangle,\\\ \frac{d}{d\tau}\langle
b_{1}^{2}\rangle=&-2i\bigg{\\{}\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle\langle b_{1}\rangle+\langle
b_{1}^{2}\rangle\bigg{\\}}-\frac{\Gamma}{\Omega}\langle b_{1}^{2}\rangle,\\\
\frac{d}{d\tau}\langle a_{1}b_{1}\rangle=&-i(1+f_{1})\bigg{\\{}(\langle
b_{1}^{2}\rangle+\langle
b_{1}^{\dagger}b_{1}\rangle)a_{1}-\frac{\Delta}{\Omega}\langle
a_{1}b_{1}\rangle+\frac{1}{2}\langle b_{1}\rangle\bigg{\\}}\\\
&-i\bigg{\\{}\langle a_{1}b_{1}\rangle+\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle\langle
a_{1}\rangle+(\frac{g_{0}}{\Omega})^{2}\langle
a_{1}\rangle\bigg{\\}}-\frac{1}{\Omega}(f_{0}^{\ast}+f_{2})\langle
a_{1}b_{1}\rangle,\\\ \frac{d}{d\tau}\langle
a_{1}b_{1}^{\dagger}\rangle=&-i(1+f_{1})\bigg{\\{}(\langle b_{1}^{\dagger
2}\rangle+\langle
b_{1}^{\dagger}b_{1}\rangle)a_{1}+(\frac{g_{0}}{\Omega})^{2}\langle
a_{1}\rangle-\frac{\Delta}{\Omega}\langle
a_{1}b_{1}^{\dagger}\rangle+\frac{1}{2}\langle
b_{1}^{\dagger}\rangle\bigg{\\}}\\\ &+i\bigg{\\{}\langle
a_{1}b_{1}^{\dagger}\rangle+\frac{P}{2}\langle
a_{1}^{\dagger}a_{1}\rangle\langle
a_{1}\rangle+(\frac{g_{0}}{\Omega})^{2}\langle
a_{1}\rangle\bigg{\\}}-\frac{1}{\Omega}(f_{0}^{\ast}+f_{2})\langle
a_{1}b_{1}^{\dagger}\rangle,\end{split}$ (27)

where we have used the relations (in terms of $a$ and $b$)

$\begin{split}&\langle a^{\dagger}\rangle=\langle a\rangle^{\ast},\ \langle
b^{\dagger}\rangle=\langle b\rangle^{\ast},\ \langle
aa^{\dagger}\rangle=\langle a^{\dagger}a\rangle+1,\\\ &\langle
bb^{\dagger}\rangle=\langle b^{\dagger}b\rangle+1,\ \langle a^{\dagger
2}\rangle=\langle a^{2}\rangle^{\ast},\ \langle b^{\dagger 2}\rangle=\langle
b^{2}\rangle^{\ast},\\\ &\langle a^{\dagger}b^{\dagger}\rangle=\langle
ab\rangle^{\ast},\ \langle a^{\dagger}b\rangle=\langle
ab^{\dagger}\rangle^{\ast}.\end{split}$ (28)

## References

  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * Jiang  _et al._ [2017] X. F. Jiang, L. B. Shao, S. X. Zhang, X. Yi, J. Wiersig, L. Wang, Q. H. Gong, M. Loncar, L. Yang, and Y. F. Xiao, Science 358, 344 (2017).
  * Sciamanna [2016] M. Sciamanna, Nat. Photon. 10, 366 (2016).
  * Monifi  _et al._ [2016] F. Monifi, J. Zhang, S. K. Özdemir, B. Peng, Y. X. Liu, F. Bo, F. Nori, and L. Yang, Nat. Photon. 10, 399 (2016).
  * Navarro-Urrios  _et al._ [2017] D. Navarro-Urrios, N. E. Capuj, M. F. Colombano, P. D. Garcia, M. Sledzinska, F. Alzina, A. Griol, A. Martinez, and C. M. Sotomayor-Torres, Nat. Commun. 8, 14965 (2017).
  * Bakemeier  _et al._ [2015] L. Bakemeier, A. Alvermann, and H. Fehske, Phys. Rev. Lett. 114, 013601 (2015).
  * Carmon  _et al._ [2007] T. Carmon, M. C. Cross, and K. J. Vahala, Phys. Rev. Lett. 98, 167203 (2007).
  * Larson and Horsdal [2011] J. Larson and M. Horsdal, Phys. Rev. A 84, 021804(R) (2011).
  * Lee  _et al._ [2009] S. B. Lee, J. Yang, S. Moon, S. Y. Lee, J. B. Shim, S. W. Kim, J. H. Lee, and K. An, Phys. Rev. Lett. 103, 134101 (2009).
  * Lü  _et al._ [2015] X. Y. Lü, H. Jing, J. Y. Ma, and Y. Wu, Phys. Rev. Lett. 114, 253601 (2015).
  * Piazza and Ritsch [2015] F. Piazza and H. Ritsch, Phys. Rev. Lett. 115, 163601 (2015).
  * Sun and Sukhorukov [2014] Y. Sun and A. A. Sukhorukov, Opt. Lett. 39, 3543 (2014).
  * Walter and Marquardt [2016] S. Walter and F. Marquardt, New J. Phys. 18, 113029 (2016).
  * Wang  _et al._ [2014] G. L. Wang, L. Huang, Y. C. Lai, and C. Grebogi, Phys. Rev. Lett. 112, 110406 (2014).
  * Yang  _et al._ [2015] N. Yang, J. Zhang, H. Wang, Y. X. Liu, R. B. Wu, L. Q. Liu, C. W. Li, and F. Nori, Phys. Rev. A 92, 033812 (2015).
  * Dalibard  _et al._ [1992] J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * Gisin and Percival [1992] N. Gisin and I. C. Percival, J. Phys. A: Math. Gen. 25, 5677 (1992).
  * Carmichael [1993] H. Carmichael,  _An Open Systems Approach to Quantum Optics Springer-Verlag_ (Springer Berlin, Heidelberg, Berlin, 1993).
  * Wiseman and Milburn [1993] H. M. Wiseman and G. J. Milburn, Phys. Rev. A 47, 1652 (1993).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Strunz  _et al._ [1999] W. T. Strunz, L. Diósi, and N. Gisin, Phys. Rev. Lett. 82, 1801 (1999).
  * Yu  _et al._ [1999] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A 60, 91 (1999).
  * Diósi  _et al._ [1998] L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A 58, 1699 (1998).
  * Strunz and Yu [2004] W. T. Strunz and T. Yu, Phys. Rev. A 69, 052115 (2004).
  * Yu [2004] T. Yu, Phys. Rev. A 69, 062107 (2004).
  * Jing and Yu [2010] J. Jing and T. Yu, Phys. Rev. Lett. 105, 240403 (2010).
  * Yang  _et al._ [2012] H. Yang, H. Miao, and Y. Chen, Phys. Rev. A 85, 040101 (2012).
  * Chen  _et al._ [2014] Y. Chen, J. Q. You, and T. Yu, Phys. Rev. A 90, 052104 (2014).
  * Xu  _et al._ [2014] J. Xu, X. Zhao, J. Jing, L. A. Wu, and T. Yu, J. Phys. A: Math. Theor. 47, 435301 (2014).
  * Moon [1992] F. C. Moon,  _Chaotic and Fractal Dynamics: An Introduction for Applied Scientists and Engineers_ (Wiley, Hoboken, 1992).
  * Strogatz [2001] S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering_ (Westview Press, Boulder, 2001).
  * Poulin [2002] D. Poulin, Physics Department and IQC, University of Waterloo (2002).
  * Field  _et al._ [1995] S. Field, N. Venturi, and F. Nori, Phys. Rev. Lett. 74, 74 (1995).
  * Field  _et al._ [1997] S. B. Field, M. Klaus, M. G. Moore, and F. Nori, Nature 388, 252 (1997).
  * Nakamura [1993] K. Nakamura,  _Quantum Chaos: A New Paradigm of Nonlinear Dynamics_ (Cambridge University Press, Cambridge, 1993).
  * Heller [1984] E. J. Heller, Phys. Rev. Lett. 53, 1515 (1984).
  * Heller [2018] E. J. Heller,  _The semiclassical way to dynamics and spectroscopy_ (Princeton University Press, 2018).
  * Gutzwiller [1990] M. C. Gutzwiller,  _Chaos in classical and quantum mechanics_ (Springer Verlag, New York, 1990).
  * Ullmo [2008] D. Ullmo, Rep. Prog. Phys. 71, 026001 (2008).
  * Wright and Weaver [2010] M. Wright and R. Weaver,  _New directions in linear acoustics and vibration: quantum chaos, random matrix theory and complexity_ (Cambridge University Press, 2010).
  * Fujisaki  _et al._ [2003] H. Fujisaki, T. Miyadera, and A. Tanaka, Phys. Rev. E 67, 066201 (2003).
  * Shchukin and Vogel [2005] E. Shchukin and W. Vogel, Phys. Rev. Lett. 95, 230502 (2005).
  * Gomes  _et al._ [2009] R. M. Gomes, A. Salles, F. Toscano, P. H. Souto Ribeiro, and S. P. Walborn, Proc. Natl. Acad. Sci. 106, 21517 (2009).
  * Chou  _et al._ [2008] C. H. Chou, T. Yu, and B. L. Hu, Phys. Rev. E 77, 011112 (2008).
  * Paz and Roncaglia [2008] J. P. Paz and A. J. Roncaglia, Phys. Rev. Lett. 100, 220401 (2008).
  * Wang and Clerk [2013] Y. D. Wang and A. A. Clerk, Phys. Rev. Lett. 110, 253601 (2013).
  * Cheng  _et al._ [2016] J. Cheng, W. Z. Zhang, L. Zhou, and W. Zhang, Sci. Rep. 6, 1 (2016).
  * Mu  _et al._ [2016] Q. Mu, X. Zhao, and T. Yu, Phys. Rev. A 94, 012334 (2016).
  * Liu  _et al._ [2011] B. H. Liu, L. Li, Y. F. Huang, C. F. Li, G. C. Guo, E. M. Laine, H. P. Breuer, and J. Piilo, Nat. Phys. 7, 931 (2011).
  * Marquardt  _et al._ [2006] F. Marquardt, J. G. E. Harris, and S. M. Girvin, Phys. Rev. Lett. 96, 103901 (2006).
  * Ludwig  _et al._ [2008] M. Ludwig, B. Kubala, and F. Marquardt, New J. Phys. 10, 095013 (2008).
  * Wolf  _et al._ [1985] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, Physica D 16, 285 (1985).
  * Furuya  _et al._ [1998] K. Furuya, M. C. Nemes, and G. Q. Pellegrino, Phys. Rev. Lett. 80, 5524 (1998).
  * Lakshminarayan [2001] A. Lakshminarayan, Phys. Rev. E 64, 036207 (2001).
  * Bandyopadhyay and Lakshminarayan [2002] J. N. Bandyopadhyay and A. Lakshminarayan, Phys. Rev. Lett. 89, 060402 (2002).
  * Wang  _et al._ [2004] X. Wang, S. Ghose, B. C. Sanders, and B. Hu, Phys. Rev. E 70, 016217 (2004).
  * Ghose and Sanders [2004] S. Ghose and B. C. Sanders, Phys. Rev. A 70, 062315 (2004).
  * Simon [2000] R. Simon, Phys. Rev. Lett. 84, 2726 (2000).
  * Duan  _et al._ [2000] L. M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 84, 2722 (2000).
  * Adesso  _et al._ [2004] G. Adesso, A. Serafini, and F. Illuminati, Phys. Rev. A 70, 022318 (2004).
  * Peres [1996] A. Peres, Phys. Rev. Lett. 77, 1413 (1996).
  * Horodecki  _et al._ [1996] R. Horodecki, M. Horodecki, and P. Horodecki, Phys. Rev. A 222, 21 (1996).
  * Horodecki  _et al._ [2001] M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Rev. A 283, 1 (2001).


# An Embarrassingly Simple Consistency Regularization Method for Semi-
Supervised Medical Image Segmentation

###### Abstract

The scarcity of pixel-level annotation is a prevalent problem in medical image
segmentation tasks. In this paper, we introduce a novel regularization
strategy involving interpolation-based mixing for semi-supervised medical
image segmentation. The proposed method is a new consistency regularization
strategy that encourages segmentation of interpolation of two unlabelled data
to be consistent with the interpolation of segmentation maps of those data.
This method represents a specific type of data-adaptive regularization
paradigm which aids to minimize the overfitting of labelled data under high
confidence values. The proposed method is advantageous over adversarial and
generative models as it requires no additional computation. Upon evaluation on
two publicly available MRI datasets: ACDC and MMWHS, experimental results
demonstrate the superiority of the proposed method in comparison to existing
semi-supervised models. Code is available at:
https://github.com/hritam-98/ICT-MedSeg

Index Terms—  Semi-supervised Learning, Medical Image Segmentation,
Interpolation, Consistency Regularization

##  1 Introduction

Supervised medical image segmentation is a widely explored problem in computer
vision, achieving exponential growth recently. These methods mostly rely upon
deep learning, requiring large-scale pixel-wise annotation data. However,
acquiring huge amounts of labelled medical data is tedious and expensive, thus
methods alleviating this requirement are highly expedient. Semi-Supervised
Learning (SSL) is a promising direction in this regard, requiring only a few
labelled data and compensating for the large portion of unlabelled data by
generating pseudo labels. Recently, SSL-based methods have been widely
recognized for their superior performance in medical image segmentation. They
not only eliminate the necessity of large-scale annotations but also produce
accurate segmentation results that are very close to those obtained from
supervised models.

Recently, Bai _et al._ [1] proposed an iterative SSL framework where pseudo
labels are generated after every iteration and refined using conditional
random field (CRF). Adversarial learning is another popular method for
utilizing unlabelled data. Zhang _et al._ [2] proposes a deep adversarial
network (DAN) that encourages the segmentation of unlabelled data to be
similar to those of labelled data. The state-of-the-art SSL segmentation model
mean-teacher network is recently extended by Yu _et al._ [3] for uncertainty-
aware left atrium segmentation.

Consistency regularization is another plausible solution in this direction,
that encourages realistic perturbations of an unlabelled image to produce
consistent segmentation maps. For example, Bortsova _et al._ [4] proposes a
transformation-consistent segmentation network, capable of exploring the
equivariance of elastic perturbations for precise lung X-ray segmentation and
Li _et al._ [5] proposes a semi-supervised segmentation framework using
transformation consistency, encouraging consistent predictions of the network-
in-training for different perturbation of the same input. Unlike these
methods, which rely upon the low-density region assumption, our proposed
method chooses perturbation directed towards another unlabelled sample,
thereby reducing the necessity of expensive gradient calculation.

On the other hand, interpolation-based regularizers have been achieving state-
of-the-art performance in various tasks and across multiple architectures.
Recently, this idea has been extended to an unsupervised setting by Berthelot
_et al._ [6] where the authors propose that utilizing realism of the latent
space interpolation from autoencoder can improve model learning. Driven by
this speculation and the aforementioned success of consistency regularization
in SSL methods, in this paper we investigate an interpolation-based mixing
technique in a semi-supervised setting and its utility in medical image
segmentation.

##  2 Methodology

Fig. 1: Overall framework of our proposed architecture.

###  2.1 Consistency Regularization with Interpolation

Consistency Regularization has been used in existing literature as a means to
enhance the robustness of a network pipeline by enforcing the network against
various perturbations on the unlabelled data, to increase the generalizability
of these new data points. Among them, the most effective perturbation would be
in the adversarial direction - in the direction almost perpendicular to the
decision boundary between the positive and negative examples, i.e., the
direction in which the network is most liable to misclassify the pixels.
However, most existing literature incorporates perturbations that may or may
not be in the adversarial direction, and thus, results in loss of
generalizability. Some other methods do perturb the input in the adversarial
direction but require a large amount of unlabelled data and thus, might not be
feasible in the biomedical domain.

Thus, we have proposed a pixel-wise data perturbation strategy as consistency
regularization in our work, which operates as described. The network is
trained in such a way so as to ensure stable and accurate segmentation of
image points interpolated from existing points. Considering two unlabelled
image data-points ${u_{1}}$ and ${u_{2}}$, we interpolate another unlabelled
image point $\mathcal{M}_{\alpha}(u_{1},u_{2})$, where
${\mathcal{M}_{\alpha}(u_{1},u_{2})=\alpha u_{1}+(1-\alpha)u_{2}}$, for some
hyperparameter $\alpha$. Now, consistency regularization is applied between
the output of the interpolated image data-point
$f(\mathcal{M}_{\alpha}(u_{1},u_{2})$ and the interpolation of the outputs of
the original unlabelled points $\mathcal{M}_{\alpha}(f(u_{1}),f(u_{2}))=\alpha
f(u_{1})+(1-\alpha)f(u_{2})$. This exploits the fact that the network learns
to predict a pixel-level segmentation mask of the input images, and further,
consistency is maintained between the outputs of the interpolated inputs and
the interpolated outputs of the original inputs.

Thus, the unlabelled samples in the datasets are used to generate the new
interpolated images and the corresponding pseudo-labels. It takes two
unlabelled images as input and returns the interpolated image and the
corresponding pseudo-label, which is used by the network pipeline. Therefore,
the Consistency Regularization technique can be summarized as:

$\mathcal{M}_{\alpha}(f_{\theta^{\prime}}(u_{1}),f_{\theta^{\prime}}(u_{2}))\simeq
f_{\theta}(\mathcal{M}_{\alpha}(u_{1},u_{2}))$ (1)

This data mixing technique would help the model learn more robust features
improving the semi-supervised learning on subsequent (target) tasks since
random perturbations do not guarantee adversarial perturbation.

###  2.2 Objective Function

Consider labelled samples $(x_{i},y_{i})\sim L_{l}$ from joint distribution
$P(X,Y)$ and unlabelled samples $(u_{i},u_{j})\sim L_{ul}$ from borderline
distribution $P(X)=\frac{P(X,Y)}{P(X|Y)}$. Using SGD for every iteration $t$,
the encoder-decoder parameter $\theta$ is updated minimising the objective
function:

$\mathcal{L}=\mathcal{L_{CE}}+r(t).\mathcal{L_{U}}$ (2)

where $\mathcal{L_{CE}}$ is the cross entropy loss applied over the labelled
data $L_{l}$ and $\mathcal{L_{U}}$ is the interpolation consistency
regularization loss applied over the unlabelled data $L_{ul}$, $r(t)$ is the
ramp function adjusting the weight of $\mathcal{L_{U}}$ after every iteration.
$\mathcal{L_{U}}$ is calculated over $(u_{i},u_{j})$ of sampled minibatches
and the pseudo labels $\tilde{y}_{i}=\mathcal{F}_{\theta^{\prime}}(u_{i})$ and
$\tilde{y}_{j}=\mathcal{F}_{\theta^{\prime}}(u_{j})$ ($\theta^{\prime}$ is the
exponential moving average of $\theta$). Next, interpolation
$u_{m}={\mathcal{M}}_{\alpha}(u_{i},u_{j})$ and model prediction
$\tilde{y}_{m}=\mathcal{F}_{\theta}(u_{m})$ are computed updating $\theta$ to
bring $\tilde{y}_{m}$ closer to the interpolation of the pseudo labels,
$\mathcal{M}_{\alpha}(\tilde{y}_{i},\tilde{y}_{j})$. The deviation in
$\tilde{y}_{m}$ and ${\mathcal{M}}_{\alpha}(\tilde{y}_{i},\tilde{y}_{j})$ is
penalised using the mean squared loss. Therefore, $\mathcal{L_{U}}$ can be
expressed as:

$\displaystyle\mathcal{L_{U}}=E_{u_{i},u_{j}\sim
L_{ul}}l(\mathcal{F}_{\theta}({\mathcal{M}}_{\alpha}(u_{i},u_{j})),$ (3)
$\displaystyle{\mathcal{M}}_{\alpha}(\mathcal{F}_{\theta^{\prime}}(u_{i}),\mathcal{F}_{\theta^{\prime}}(u_{j})))$

Our overall approach is depicted in Figure 1 and algorithm 1.

Input:

$L_{l}$: Distribution of labelled samples; $L_{ul}$: Distribution of
unlabelled samples

Define:

$f_{\theta}(\cdot)$: Segmentation network with trainable parameter $\theta$

$f_{\theta^{\prime}(\cdot)}$: Segmentation network with parameter
$\theta^{\prime}$\- exponential moving average of $\theta$

$T$: Total number of iterations; $r(t)$: ramp function

$\lambda$: exponential moving average change rate

$\mathcal{M}_{\alpha}(u_{1},u_{2})=\alpha u_{1}+(1-\alpha)u_{2}$

while (t $\leq$ T) do

Sample labelled mini-batch $\rhd\>\\{(x_{p},y_{p})\\}_{p=1}^{P}\sim L_{l}$

Supervised CE loss
$\rhd\>\mathcal{L}_{\mathcal{CE}}(\\{(f_{\theta}(x_{p}),y_{p})\\}_{p=1}^{P})$

Sample two unlabelled batches $\rhd\>\\{(u_{i},u_{j})\\}_{k=1}^{U}\sim L_{ul}$

Generate pseudo-labels
$\rhd\>\\{\tilde{y_{i}},\tilde{y_{j}}\\}_{k=1}^{U}=\\{f_{\theta^{\prime}}(u_{i},u_{j})\\}_{k=1}^{U}$

Interpolation
$\rhd\>u_{m}=\mathcal{M}_{\alpha}(u_{i},u_{j}),y_{m}=\mathcal{M}_{\alpha}(\tilde{y_{i}},\tilde{y_{j}})$

Interpolated pseudo-label $\rhd\>\tilde{y}_{m}=f_{\theta}(u_{m})$

Unsupervised loss
$\rhd\>\mathcal{L}_{\mathcal{U}}=MSE(\\{y_{m},\tilde{y}_{m}\\}_{m=1}^{U})$

Overall model loss
$\rhd\>\mathcal{L}=\mathcal{L}_{\mathcal{CE}}+r(t)\mathcal{L}_{\mathcal{U}}$

Gradient computation
$\rhd\>\mathcal{G}_{\theta}\longleftarrow\mathcal{L}\cdot\nabla_{\theta}$

Update parameter
$\rhd\>\theta^{\prime}\longleftarrow(1-\lambda)\theta+\lambda\theta^{\prime}$

$\theta\longleftarrow step(\mathcal{G}_{\theta},\theta)$

end while

return $\theta$

Algorithm 1 Pseudo-code of our proposed method.

##  3 Experiments and Results

Table 1: Comparison of the proposed method with State-of-the-art methods on
the ACDC 2017 Dataset and MMWHS dataset.

Methods | ACDC 2017 | MMWHS
---|---|---
| labelled data
---
(as % of the training set)
DSC(%) | HD (mm) |  | labelled data
---
(as % of the training set)
DSC(%) | HD(mm)
Peng et al.[7] | 5% | 85.76 | - | 13.30% | 55.75 | -
Mix-Up[8] | 10% | 86.3 | - | 50% | 79.6 | -
Chaitanya et al[9] | 10% | 88.6 | - | 50% | 79.4 | -
Adversarial Training[2] | 20% | 79.1 | 5.16 | 50% | 77.9 | 3.20
Ours | 10% | 89.8 | 4.47 | 40% | 79.83 | 3.05

###  3.1 Dataset and Experimental Setup

The experiments for evaluating the model performance was carried out on two
public datasets: the ACDC 2017 [10] dataset consisting of 100 short-axis
cardiac MRI volumes with expert annotations for three structures: left and
right ventricle and myocardium and the MMWHS [11] dataset, hosted in MICCAI
2017, consisting of 20 cardiac MRI samples with expert annotations for seven
structures: left and right ventricle, left and right atrium, pulmonary artery,
myocardium, and ascending aorta. The datasets were split into train,
validation, and test sets, where 2 volumes were considered for validation
purposes in both the datasets. The test set consists of 20 and 5 volumes for
ACDC 2017 and MMWHS, respectively. The rest of the data was used for training.
To study the segmentation performance with respect to the labelled data, we
have experimented with 1.25%, 2.5%, and 10% labelled volumes for ACDC 2017 and
10%, 20%, and 40% labelled volumes for the MMWHS dataset, following the
existing works in literature [9]. We have used ResNet-50 as the encoder
backbone of our architecture, using an ADAM optimizer with an initial learning
rate of 1e-5. For evaluation purposes, we have used three widely used metrics:
Dice Similarity Score (DSC), Average Symmetric Distance (ASD), and Hausdorff
Distance (HD). Average of all the metric scores over all the classes and
reported in this paper. Mixing parameter $\alpha$ was set experimentally.

###  3.2 Results on the ACDC 2017 Dataset

We achieve DSC of 73.56%, 79.05%, and 89.80% on the ACDC 2017 dataset for
1.25%, 2.5%, and 10% labelled volumes respectively. For the same set of
experiments, ASD of 1.981, 1.546, and 1.229 and HD of 6.719, 5.820, and 4.473
were obtained respectively. The merits of the model can be better appreciated
by referring to the comparison with other state-of-the-art methods, shown in
Table 1, where we have compared our obtained results using 10% labelled
volumes only. A higher Dice score at a lower percentage of labelled data
indicates that the model is capable of generalizing the data efficiently by
making use of a lesser number of labelled samples. As seen in the table, our
approach uses a low number of labelled samples and produces the best results
in terms of DSC and HD, whereas methods like [2] use double the amount of
labelled data, yet achieve inferior performance as compared to our results,
proving the efficiency and robustness of our model. Only Peng _et al._ [7]
uses a lower number of labelled data, but produces poorer results as well. Our
method produces a reasonably better result by making use of only 10% of the
labelled data available for training.

Fig. 2: Comparison of results obtained by using different percentages of
labelled data in training set on ACDC 2017 (left) and MMWHS (right) dataset.

###  3.3 Results on the MMWHS Dataset

Upon evaluation on the MMWHS dataset, DSC of 66.86%, 71.24%, and 79.83% have
been achieved using 10%, 20%, and 40% labelled volumes, respectively. ASD
scores of 2.671, 2.207, and 1.819 and HD values of 5.073, 4.288, and 3.047
were achieved for these experiments, respectively. The efficacy of the
proposed technique can be verified by comparison with existing state-of-the-
art techniques as shown in Table 1, where we have compared our results
obtained using 40% annotated data only, yet achieving superior performance as
compared to methods with $50\%$ labelled data. The results were reported from
the papers directly. All the other parameters were the same as those methods,
to maintain fair comparison.

Fig. 3: Comparative image of the segmentation masks of our model, along with
the other semi-supervised models

Peng _et al._ [7] use the least number of labelled samples in their method,
but they have produced vastly substandard results. Whereas, our method obtains
better results as compared to most state-of-the-art methods, by using 40%
labelled data in the entire training set. The results achieved by Zhang _et
al._ [8, 9] are very close to ours, however, they have used a larger
percentage of labelled data for training.

Figure 2 represents the results obtained using different percentages of
labelled data from the training set. It is obvious that employing more
labelled data while training will produce superior performance. However, we
can observe from the figure that our model can produce results very close to
the supervised method (that uses 100% labelled train data) by employing only a
small percentage (10% for ACDC 2017 and 40% for MMWHS) of annotations. Figure
3 offers a visual comparison of the proposed method with a few of the state-
of-the-art methods existing in contemporary literature. A careful inspection
of the image shows that our method generates labels that are the closest to
the ground truth (GT) when compared to the other methods. ASD scores were not
reported in most of the methods in the literature for both the datasets and
hence not considered for comparison.

##  4 Conclusion and Future Work

Recently, machine learning has had a transformative impact on various
applications, but its effectiveness is impacted due to the lack of labelled
data, especially in medical imaging, where availing of those annotations is
extremely difficult. Here we have tried to address this shortcoming by
developing a semi-supervised learning strategy that encourages consistency
regularization by interpolation. It is advantageous over the previous SSL
models in multiple aspects: unlike adversarial perturbations or generative
models, it requires almost no additional computations. Moreover, our proposed
model outperforms several state-of-the-art SSL methods by utilizing fewer
annotations than these methods. In future, we plan to extend the work by
interpolating at the hidden representations to boost the performance and
evaluate the performance in general computer vision problems to assess its
robustness. We further plan to extend the experimentation by adding more
variability between the two interpolating images.

##  5 Compliance with Ethical Standards

This research study was conducted retrospectively using human subject data
made available in open access. Ethical approval was not required as confirmed
by the license attached with the open-access data.

## References

  * [1] W. Bai, O. Oktay, M. Sinclair, H. Suzuki, M. Rajchl, G. Tarroni, B. Glocker, A. King, P. Matthews, and D. Rueckert,  “Semi-supervised learning for network-based cardiac mr image segmentation,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2017, pp. 253–260.
  * [2] Y. Zhang, L. Yang, J. Chen, M. Fredericksen, D. Hughes, and D. Chen,  “Deep adversarial networks for biomedical image segmentation utilizing unannotated images,”  in International conference on medical image computing and computer-assisted intervention. Springer, 2017, pp. 408–416.
  * [3] L. Yu, S. Wang, X. Li, C. Fu, and P. Heng,  “Uncertainty-aware self-ensembling model for semi-supervised 3d left atrium segmentation,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2019, pp. 605–613.
  * [4] G. Bortsova, F. Dubost, L. Hogeweg, I. Katramados, and M. de Bruijne,  “Semi-supervised medical image segmentation via learning consistency under transformations,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2019, pp. 810–818.
  * [5] X. Li, L. Yu, H. Chen, C. Fu, L. Xing, and P. Heng,  “Transformation-consistent self-ensembling model for semisupervised medical image segmentation,”  IEEE Transactions on Neural Networks and Learning Systems, vol. 32, no. 2, pp. 523–534, 2020.
  * [6] D. Berthelot, C. Raffel, A. Roy, and I. Goodfellow,  “Understanding and improving interpolation in autoencoders via an adversarial regularizer,”  arXiv preprint arXiv:1807.07543, 2018.
  * [7] Jizong Peng, Marco Pedersoli, and Christian Desrosiers,  “Boosting semi-supervised image segmentation with global and local mutual information regularization,”  arXiv preprint arXiv:2103.04813, 2021.
  * [8] H. Zhang, M. Cisse, Y. D., and D. Lopez-Paz,  “mixup: Beyond empirical risk minimization,”  arXiv preprint arXiv:1710.09412, 2017.
  * [9] K. Chaitanya, E. Erdil, N. Karani, and E. Konukoglu,  “Contrastive learning of global and local features for medical image segmentation with limited annotations,”  Advances in Neural Information Processing Systems, 2020.
  * [10] O. Bernard, A. Lalande, C. Zotti, F. Cervenansky, X. Yang, et al.,  “Deep learning techniques for automatic mri cardiac multi-structures segmentation and diagnosis: Is the problem solved?,”  IEEE transactions on medical imaging, vol. 37, no. 11, pp. 2514–2525, 2018.
  * [11] X. Zhuang and J. Shen,  “Multi-scale patch and multi-modality atlases for whole heart segmentation of mri,”  Medical image analysis, vol. 31, pp. 77–87, 2016.