# On the solid angle of convex sets

J. Bruna, J. Cufí, E. Gallego and A. Reventós  Departament de Matemàtiques
Universitat Autònoma de Barcelona

08193 Bellaterra, Barcelona
Catalonia<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Here we analyze three dimensional analogues of the classical Crofton’s formula
for planar compact convex sets. In this formula a fundamental role is played
by the visual angle of the convex set from an exterior point. A generalization
of the visual angle to convex sets in euclidian space is the visual solid
angle. This solid angle, being an spherically convex set in the unit sphere,
has length, area and other geometric quantities to be considered. The main
goal of this note is to express invariant quantities of the original convex
set depending on volume, surface area and mean curvature integral by means of
integrals of functions related to the solid angle.

###### Key words and phrases:

Invariant measures, convex set, dihedral angle, solid angle, constant width.

###### 1991 Mathematics Subject Classification:

Primary 52A15, Secondary 53C65.

The authors were partially supported by grants 2021SGR01015 (Generalitat de
Catalunya) and PID2021-125625NB-I00. ( Ministerio de Ciencia e Innovación).

##  1\. Introduction and statement of results

The purpose of this note is to give three dimensional analogues of the
classical Crofton’s formula for a planar compact convex set $K$,

(1.1) $\int_{P\notin K}2(\omega-\sin\omega)\,dP=L^{2}-2\pi F.$

Here $L$ is the length of the boundary of $K$, $F$ its area and
$\omega=\omega(P)$ is the _visual angle_ of $K$ as seen from $P$.

A generalization of the visual angle in the plane to compact convex sets in
the euclidean space $\mathbb{E}^{3}$ is the visual solid angle. The solid
angle $\Omega(P)$ of a compact convex set $K$ from a point $P\notin K$ is the
set of unit directions $u$ such that the ray $P+tu,t\geq 0$ meets $K$. Instead
of an arc in the unit cercle $S^{1}$ and its length $\omega$ we have a
(spherically convex) set in the unit sphere $S^{2}$ with a richer geometry,
having area, length, and other geometric quantities that might be considered.

In our analysis we will be lead to two set functions
$\alpha(\Omega),\beta(\Omega)$ defined for spherically convex sets
$\Omega\subset S^{2}$, both replacing $\omega-\sin\omega$ in different senses.
To introduce them we first recall classical notions about spherical geometry,
the main reference being [5]. The radii perpendicular to the support planes of
the cone spanned by $\Omega$ form another cone whose intersection with $S^{2}$
is the so-called _dual curve_ of $\partial\Omega$. The region in $S^{2}$
bounded by this curve and its symmetrical with respect the origin, which we
denote by $\tilde{\Omega}$, consists in the unit directions $v$ such that the
plane $v^{\bot}$ meets $\Omega$. We may identify it with the set of great
circles in $S^{2}$ meeting $\Omega$. Its complement $\tilde{\Omega}^{c}$ in
$S^{2}$ consists of two symmetrical components. Since the scalar product
$\langle v,u\rangle$ does not vanish for $v\in\tilde{\Omega}^{c},u\in\Omega$,
this quantity has constant sign in each component. The one with positive sign
is called the _dual solid cone_ , denoted $\Omega^{*}$.

The set functions $\alpha,\beta$ are respectively defined as

$\alpha(\Omega)=\frac{1}{2}\int\limits_{u\in\Omega,\,v\in\tilde{\Omega}}|\langle
u,v\rangle|\,du\,dv,$

(1.2)
$\beta(\Omega)=\frac{1}{8}\int\limits_{v_{i}\in\tilde{\Omega}}|\det(v_{1},v_{2},v_{3})|\,dv_{1}\,dv_{2}\,dv_{3},$

where $du$ denotes the Lebesgue measure in $S^{2}$.

The factors in front of the integrals are explained by the fact that the set
$\tilde{\Omega}$ doubly parametrizes planes through the origin meeting
$\Omega$, equivalently great circles meeting $\Omega$, so the $v$-integrals
are in fact over this set of planes, see below.

It is immediate to check that both set functions are invariant by rigid
motions $T$ of the sphere, that is
$\alpha(T(\Omega))=\alpha(\Omega),\beta(T(\Omega))=\beta(\Omega)$. We point
out that a straightforward computation shows that the formal analogues in the
plane

$\alpha(I)=\frac{1}{2}\int\limits_{u\in I,\,v\in\tilde{I}}|\langle
u,v\rangle|\,du\,dv,\quad\beta(I)=\frac{1}{4}\int\limits_{v_{i}\in\tilde{I}}|\det(v_{1},v_{2})|\,dv_{1}\,dv_{2},$

where now $I$ is an interval, both equal $\omega-\sin\omega$ up to a constant,
$\omega$ being the length of $I$.

The Crofton type formulas we obtain are then:

######  Theorem 1.1.

For a compact convex set $K\subset\mathbb{E}^{3}$ with mean curvature integral
$M$, volume $V$ and surface area $F$ one has

(1.3) $\frac{1}{2}\pi MF-2\pi^{2}V=\int_{P\notin K}\alpha(\Omega(P))\,dP.$

######  Theorem 1.2.

For a compact convex set $K$ in $\mathbb{E}^{3}$ with mean curvature integral
$M$ and volume $V$ one has

(1.4) $M^{3}-\pi^{4}V=\int_{P\notin K}\beta(\Omega(P))\,dP.$

######  Theorem 1.3.

With the same notations,

(1.5) $\int L(K\cap E)^{2}\,dE=\int_{P\notin K}|\Omega(P)|^{2}\,dP+4\pi^{2}V,$

where $dE$ is the invariant measure for affine planes $E$ in space, $L$
denotes length and $|\Omega|$ surface measure.111Formula (1.5) was previously
considered in private conversations with E. Teufel in 2003.

As a consequence of Theorem 1.1 and Minkowski’s inequality $12\pi V\leq MF$
(see [2]) one has

$\int_{P\notin K}\alpha(\Omega(P))\,dP\geq 4\pi^{2}V,$

with equality only when $K$ is a ball.

Theorem 1.3 has as a consequence:

######  Theorem 1.4.

For all convex sets

(1.6) $\int_{P\notin K}|\Omega(P)|^{2}\,dP\geq 4\pi^{2}V,$

and equality holds if and only if $K$ is a ball.

Part of our analysis will consist in understanding the set functions
$\alpha,\beta$ in terms of metric properties of $\Omega$. For the set function
$\alpha$ a satisfactory description is easily obtained, namely

(1.7) $\alpha(\Omega)=\pi|\Omega|-\langle c(\Omega),c(\Omega^{*})\rangle,$

where $c(\Omega)=\int_{\Omega}u\,du$ and similarly $c(\Omega^{*})$ are
(unweighted) centroids. Again, the analogue of this expression for an arc $I$
in $S^{1}$ ($I^{*}$ being in this case the concentric arc with length
$\pi-\omega$) equals $\omega-\sin\omega$ up to constants. It is possible to
express everything just in terms of $\Omega$. The explicit computation of
$\alpha(\Omega)$ is possible just for spherical caps and other simple cases.

For the set function $\beta$ the description is not so neat. In fact we show
that the difference $\beta(\Omega)-\frac{\pi^{2}}{2}\alpha(\Omega)$ can be
expressed in terms of the dihedral visual angles ${\mathcal{D}}(\Omega,u)$ of
$\Omega$ from points $u\in S^{2}$ not in $\Omega$, the analogue in spherical
geometry of the visual angle. Note that ${\mathcal{D}}(\Omega,u)$ is the angle
between two planes through the origin and $u$ tangent to $\Omega$. From this
it can be seen that a linear combination of equations (1.3) and (1.4) is the
classical Crofton-Herglotz formula, so Theorem 1.1 and Theorem 1.2 can be seen
as equivalent statements (see section 2.3).

Regarding theorem 1.3, it would be interesting to understand the left-hand
side of (1.5) in terms of the geometry of $K$.

##  2\. Proofs of theorems

###  2.1.  Proof of theorems 1.1, 1.2, 1.3 and 1.4

All of them are obtained by mimicking the integral geometry proof of the plane
Crofton formula (1.1). We denote by $E,G$ affine planes and lines in space,
respectively, and by $dE,dG$ the corresponding canonical invariant measures as
used for instance in [5].

From the classical Crofton’s formulas for intrinsic volumes, mean curvature
integral $M$, area of the boundary $F$ and volume $V$ (see [5]) we have, for a
given compact convex set $K\subset\mathbb{E}^{3}$,

(2.1) $\int_{G\cap K\neq\emptyset}dG=\frac{\pi}{2}F,\quad\int L(K\cap
G)dG=2\pi V,\\\ $

and

(2.2) $\int_{E\cap K\neq\emptyset}dE=M,\int L(K\cap
E)dE=\frac{\pi^{2}}{2}F,\int A(K\cap E)dE=2\pi V.$

We shall consider pairs and triples of linear varieties provided with the
product measure. It will be useful to express these measures using the
parametrizations in next lemma, whose proof can be deduced from section 12 of
[5]:

######  Lemma 2.1.

  1. The following holds:

  2. a)

(2.3) $dG\ dE=|\langle u,v\rangle|dG_{P}\ dE_{P}\ dP$

where $P=G\cap E$, $v$ is a unit normal direction of $E$, $u$ is a unit
direction of $G$, $dP$ is the Lebesgue measure on $\mathbb{E}^{3}$, $dG_{P}$
is the measure of lines in $\mathbb{E}^{3}$ through $P$ and $dE_{P}$ the
measure of planes in $\mathbb{E}^{3}$ through $P$.

  3. b)

(2.4) $dE^{1}\ dE^{2}\ dE^{3}=|\det(v_{1},v_{2},v_{3})|dE^{1}_{P}\ dE^{2}_{P}\
dE^{3}_{P}\ dP.$

where $v_{i}$ are the unit normal directions to $E_{i}$ and $P=E_{1}\cap
E_{2}\cap E_{3}.$

  4. c)

For pairs of lines $G^{1},G^{2}$ intersecting at a point $P\in\mathbb{E}^{3}$
we have

(2.5) $dG_{P}^{1}\,dG_{P}^{2}\,dP=dG^{1}_{E}\,dG^{2}_{E}\,dE$

where $dG_{E}$ is the measure of lines within $E$.

In all statements, we can identify $dG_{P},dE_{P}$ with
$\frac{1}{2}du,\frac{1}{2}dv$ respectively.

We can now proceed to prove the announced theorems.

Proof of Theorem 1.1. We consider pairs $(E,G)$ of planes and lines both
meeting $K$. From (2.1), (2.2) and (2.3),

$\frac{\pi}{2}MF=\int\limits_{G\cap K\neq\emptyset,E\cap
K\neq\emptyset}dG\,dE=\int_{\mathbb{E}^{3}}\int\limits_{G_{P}\cap
K\neq\emptyset,E_{P}\cap K\neq\emptyset}|\langle u,v\rangle|dG_{P}\ dE_{P}\
dP.$

If $P\in K$, there is no restriction on $G_{P},E_{P}$. Since they are doubly
parametrized bu $u,v$ respectively,

$\int|\langle u,v\rangle|dG_{P}\ dE_{P}=\frac{1}{4}\int_{u,v\in S^{2}}|\langle
u,v\rangle|du\,dv.$

Obviously the $v$ integral does not depend on $u$, so the above equals

$\pi\int|\langle u,v\rangle|\,dv.$

Choosing $u=(0,0,1)$ and computing in spherical coordinates one obtains the
value $2\pi^{2}$.

If $P\notin K$ then $G_{P}$ is parametrized by $u\in\Omega(P)$ and $E_{P}$ is
doubly parametrized by $v\in\tilde{\Omega}(P)$ whence

$\int\limits_{G_{P}\cap K\neq\emptyset,E_{P}\cap K\neq\emptyset}|\langle
u,v\rangle|dG_{P}\
dE_{P}=\frac{1}{2}\int_{u\in\Omega(P),v\in\tilde{\Omega}(P)}|\langle
u,v\rangle|\,du\,dv=\alpha(\Omega(P)),$

thus proving (1.3). ∎

Proof of Theorem 1.2. Here we use triples of planes meeting $K$ and proceed
analogously to the above proof. Using (2.2) and (2.4) it follows

$M^{3}=\int_{\mathbb{E}^{3}}\int\limits_{E_{P}^{i}\cap
K\neq\emptyset}|\det(v_{1},v_{2},v_{3})|dE^{1}_{P}\ dE^{2}_{P}\ dE^{3}_{P}\
dP.$

Again, if $P$ is within $K$, there is no restriction on the $E_{P}^{i}$ and

$\int|\det(v_{1},v_{2},v_{3})|dE^{1}_{P}\ dE^{2}_{P}\
dE^{3}_{P}=\frac{1}{8}\int_{v_{i}\in
S^{2}}|\det(v_{1},v_{2},v_{3})|\,dv_{1}\,dv_{2}\,dv_{3}.$

The integral in $v_{1},v_{2}$ is independent of $v_{3}$ so the above integral
equals

$\frac{\pi}{2}\int|\det(v_{1},v_{2},v_{3})|\,dv_{1}\,dv_{2}.$

Choosing $v_{3}=(0,0,1)$ and computing in spherical coordinates we get the
value $\pi^{4}$.

If $P\notin K$ then $E_{P}^{i}$ are doubly parametrized by
$v_{i}\in\tilde{\Omega}(P)$ so that

$\int\limits_{E_{P}^{i}\cap
K\neq\emptyset}|\det(v_{1},v_{2},v_{3})|dE^{1}_{P}\,dE^{2}_{P}\,dE^{3}_{P}=\beta(\Omega(P)),$

and (1.4) is proved. ∎

Proof of Theorem 1.3. We use now pairs of intersecting lines, both meeting
$K$. The measure of this set of lines is

$\int_{P\in\mathbb{E}^{3},G_{P}^{i}\cap
K\neq\emptyset}dG_{P}^{1}\,dG_{P}^{2}\,dP=\int_{\mathbb{E}^{3}}\bigg{(}\int_{G_{P}^{i}\cap
K\neq\emptyset}dG_{P}^{1}\,dG_{P}^{2}\bigg{)}\,dP.$

The contribution of $K$ in the $dP$ integral is $(2\pi)^{2}V$, while that of
$K^{c}$ is $\int_{P\notin K}|\Omega(P)|^{2}\,dP$. On the other hand, by (2.5),
this equals

$\int\bigg{(}\int_{G_{1},G_{2}\subset E,G_{i}\cap
K\neq\emptyset}dG^{1}_{E}\,dG^{2}_{E}\bigg{)}dE.$

Since $\int_{G\subset E,G\cap K\neq\emptyset}\,dG=L(K\cap E)$ (the Cauchy-
Crofton formula in the plane $E$), we are done. ∎

We point out that this argument is equivalent to integration on $E$ of the
planar Crofton formula (1.1) for $E\cap K$ within $E$.

Proof of Theorem 1.4. Using the isoperimetric inequality in the plane $E$ the
left hand side of (1.5) is bigger than

$4\pi\int_{E}A(K\cap E)\,dE,$

which by the last equality in (2.2) equals $8\pi^{2}V$, thus proving (1.6). If
equality holds, then $L(K\cap E)^{2}=4\pi A(K\cap E)$ for all $E$, then all
$K\cap E$ are discs, and this implies easily that $K$ is a ball. ∎

###  2.2. On the set function $\alpha$

To prove the relation (1.7) just notice that

$\int_{v\in\tilde{\Omega}}|\langle u,v\rangle|\,dv=\int_{v\in S^{2}}|\langle
u,v\rangle|\,dv-2\int_{v\in\Omega^{*}}|\langle u,v\rangle|\,dv.$

The first integral does not depend on $u$ and equals $2\pi$, while in the
second one $\langle u,v\rangle$ is positive. Altogether,

$\alpha(\Omega)=\pi\int_{u\in\Omega}du-\int_{u\in\Omega,v\in\Omega^{*}}\langle
u,v\rangle\,du\,dv=\pi|\Omega|-\langle c(\Omega),c(\Omega^{*})\rangle.$

When $\Omega$ is a cap in $S^{2}$ with sperical radius $\omega$ then

$\displaystyle\alpha(\Omega)=2\pi^{2}(1-\cos\omega)-\pi^{2}\cos^{2}\omega\sin^{2}\omega.$

In order to express this function in terms of $\Omega$ we find a relation
between the centroids of $\Omega$ and $\Omega^{*}$ using a parametrization of
the boundary of $\Omega$. We consider the orientation in $\Omega$ given by the
unit outward normal to $S^{2}$ and let $\gamma(t)$, $0\leq t\leq\ell$ be the
arc-length parametrization of $\partial\Omega$ with the induced orientation,
so $T=\gamma^{\prime}$ is the unit tangent.

######  Proposition 2.2.

  1. a)

$c(\Omega)=\frac{1}{2}\int_{0}^{\ell}\gamma(t)\times\gamma^{\prime}(t)\
dt,\quad c(\Omega^{*})=\frac{1}{2}\int_{0}^{\ell}k_{g}(t)\gamma(t)\ dt,$

where $k_{g}(t)$ is the geodesic curvature of $\gamma(t)$.

  2. b)

$c(\Omega)+c(\Omega^{*})=\frac{1}{2}\int_{0}^{\ell}\gamma^{\prime}(t)\times\gamma^{\prime\prime}(t)\
dt=\frac{1}{2}\int_{0}^{\ell}k(t)\vec{B}(t)\ dt.$

where $k(t)$ is the curvature of $\gamma(t)$ and $\vec{B}(t)$ its binormal.

###### Proof.

If $u=(x,y,z)$, the first component of $c(\Omega)$ is $\int_{\Omega}xdu$, the
flow through $\Omega$ of the vector field $X=(1,0,0)$. Since $X=\nabla\times
Y,Y=\frac{1}{2}(0,-z,y)$ this component equals

$\frac{1}{2}\int_{\partial\Omega}\langle T,Y(\gamma(t))\rangle\,dt.$

Now $\langle T,Y\rangle$ equals the first component of
$\gamma\times\gamma^{\prime}(t)$. Similarly the other components, so the first
formula in a) is proved. Next, notice that
$\gamma^{*}(t)=\gamma(t)\times\gamma^{\prime}(t)$ parametrizes the dual curve,
the boundary of $\Omega^{*}$, whence one has as well

$c(\Omega^{*})=\frac{1}{2}\int_{0}^{\ell}\gamma^{*}(t)\times(\gamma^{*})^{\prime}(t)\
dt.$

Now,

$\gamma^{*}\times{\gamma^{*}}^{\prime}=(\gamma\times\gamma^{\prime})\times(\gamma\times\gamma^{\prime})^{\prime}=(\gamma\times\gamma^{\prime})\times(\gamma\times\gamma^{\prime\prime})=\det(\gamma,\gamma^{\prime},\gamma^{\prime\prime})\gamma=k_{g}\gamma$

and a) is proved.

In order to prove b) we simplify the notation writing
$\sigma=\gamma\times\gamma^{\prime}+\gamma^{*}\times{\gamma^{*}}^{\prime}.$
Denote by $\vec{T},\vec{N},\vec{B}$ the Frenet frame of $\gamma$. From the
definitions it is easy to see that $\langle\sigma,\vec{T}\rangle=0$. Also

$\langle\sigma,\vec{N}\rangle=\frac{1}{k}\langle\gamma\times\gamma^{\prime},\gamma^{\prime\prime}\rangle+\frac{1}{k}\langle\det(\gamma,\gamma^{\prime},\gamma^{\prime\prime})\gamma,\gamma^{\prime\prime}\rangle=\frac{1}{k}\langle\gamma\times\gamma^{\prime},\gamma^{\prime\prime}\rangle-\frac{1}{k}\det(\gamma,\gamma^{\prime},\gamma^{\prime\prime})=0$

because $\langle\gamma,\gamma^{\prime}\rangle=0$ and so
$\langle\gamma,\gamma^{\prime\prime}\rangle=-1.$ Now we compute
$\langle\sigma,\vec{B}\rangle,$

$\displaystyle\langle\sigma,\vec{B}\rangle$ $\displaystyle=$
$\displaystyle\langle\sigma,\vec{T}\times\vec{N}\rangle=\frac{1}{k}\langle\sigma,\gamma^{\prime}\times\gamma^{\prime\prime}\rangle=\frac{1}{k}\langle\gamma\times\gamma^{\prime}+\gamma^{*}\times{\gamma^{*}}^{\prime},\gamma^{\prime}\times\gamma^{\prime\prime}\rangle=$
$\displaystyle=$
$\displaystyle\frac{1}{k}\langle\gamma\times\gamma^{\prime},\gamma^{\prime}\times\gamma^{\prime\prime}\rangle+\frac{1}{k}\langle\gamma^{*}\times{\gamma^{*}}^{\prime},\gamma^{\prime}\times\gamma^{\prime\prime}\rangle=$
$\displaystyle=$
$\displaystyle-\frac{1}{k}\langle\gamma,\gamma^{\prime\prime}\rangle\langle\gamma^{\prime},\gamma^{\prime}\rangle+\frac{1}{k}\langle
k_{g}\gamma,\gamma^{\prime}\times\gamma^{\prime\prime}\rangle=\frac{1}{k}(1+k_{g}^{2}).$

The curve $\gamma$ being on the unit sphere we have that $k^{2}=1+k_{g}^{2}$;
therefore $\langle\sigma,\vec{B}\rangle=k$ and we conclude that

$\gamma\times\gamma^{\prime}+\gamma^{*}\times{\gamma^{*}}^{\prime}=k\vec{B}.$

Since $\gamma^{\prime}\times\gamma^{\prime\prime}=\vec{T}\times
k\vec{N}=k\vec{B}$ the proposition is proved.

∎

###  2.3. On the set function $\beta$

Here we wish to obtain an alternative expression for (1.2) which will lead us
to the Crofton-Herglotz formula.

Let $u=v_{2}\times v_{3}/|v_{2}\times v_{3}|$. To specify a basis for
$u^{\bot}$ we write $u$ in spherical coordinates,
$u=(\sin\varphi\cos\theta,\sin\varphi\sin\theta,\cos\varphi)$, and define

$e_{1}=\frac{\partial}{\partial\varphi}=(\cos\varphi\cos\theta,\cos\varphi\sin\theta,-\sin\varphi),\qquad
e_{2}=\frac{1}{\sin\varphi}\frac{\partial}{\partial\theta}=(-\sin\theta,\cos\theta,0),$

so that $\\{e_{1},e_{2},v\\}$ is a positive orthonormal basis. We write
$v_{2},v_{3}$ in this basis

$v_{2}=\cos\theta_{2}\cdot e_{1}+\sin\theta_{2}\cdot e_{2},\quad
v_{3}=\cos\theta_{3}\cdot e_{1}+\sin\theta_{3}\cdot e_{2}.$

Then $v_{2},v_{3}$ are parametrized by $u,\theta_{2},\theta_{3}$. Keeping in
mind that the integral is in fact over the set of triplets of great circles
meeting $\Omega$, so $\pm v_{2},\pm v_{3}$ count as one, we let
$0\leq\theta_{2},\theta_{3}\leq\pi$ and require that the angles
$\theta_{2},\theta_{3}$ be within the dihedral angle ${\mathcal{D}}(\Omega,u)$
determined by $\Omega$ and $u$. Then $u\in S^{2},\theta_{2},\theta_{3}$ is a
double parametrization of the set of pairs of great circles meeting $\Omega$.
In this parametrization, it is immediate to check that (cf. [4, (34.1)])

$dv_{2}dv_{3}=|\sin(\theta_{3}-\theta_{2})|d\theta_{2}d\theta_{3}dv,$

while

$|\det(v_{1},v_{2},v_{3})|=|v_{2}\times v_{3}|\cdot|\langle
v_{1},u\rangle|=|\sin(\theta_{3}-\theta_{2})|\cdot|\langle v_{1},u\rangle|.$

Thus

$\beta{\Omega)=\frac{1}{4}\int_{v_{1}\in\tilde{\Omega},u\in
S^{2},\theta_{i}\in{\mathcal{D}}(\Omega,u)}}\sin^{2}(\theta_{3}-\theta_{2})|\langle
v_{1},u\rangle|dv_{1}\,d\theta_{2}\,d\theta_{3}\,du$

Now the integral with respect to $\theta_{2},\theta_{3}$ is easily computed.
Denoting as well by ${\mathcal{D}}(\Omega,u)$ the measure of the dihedral
angle we get

$\beta(\Omega)=\frac{1}{8}\int_{v_{1}\in\tilde{\Omega},\,u\in
S^{2}}({\mathcal{D}}^{2}(\Omega,u)-\sin^{2}{\mathcal{D}}(\Omega,u))|\langle
v_{1},u\rangle|dv_{1}\,du.$

For $\pm u\in\Omega$ one has ${\mathcal{D}}(\Omega,u)=\pi$ whence the
contribution of this part equals $\pi^{2}\alpha(\Omega)/2$. Thus

$\beta(\Omega)=\frac{\pi^{2}}{2}\alpha(\Omega)+\gamma(\Omega),$

with

$\gamma(\Omega)=\frac{1}{8}\int_{v\in\tilde{\Omega},\,\pm
u\notin\Omega}({\mathcal{D}}^{2}(\Omega,u)-\sin^{2}{\mathcal{D}}(\Omega,u))|\langle
v,u\rangle|dv\,du.$

Thus theorems 1.1 and 1.2 imply

$M^{3}-\frac{1}{4}\pi^{3}MF=\int_{P\notin K}\gamma(\Omega(P))\,dP.$

We now insert the definition of $\gamma(\Omega(P))$ and use (2.3). If $u$ the
unit direction of $G,{\mathcal{D}}(\Omega(P),u)$ is the dihedral angle
${\mathcal{D}}(K,G)$ of $K$ as seen from $G$ so the right-hand side above
equals

$\frac{1}{2}\int_{E\cap K\neq\emptyset,G\cap
K=\emptyset}({\mathcal{D}}^{2}(K,G)-\sin^{2}{\mathcal{D}}(K,G))dE\,dG=$
$=\frac{1}{2}\int_{E\cap K\neq\emptyset}dE\bigg{(}\int_{G\cap
K=\emptyset}({\mathcal{D}}^{2}(K,G)-\sin^{2}{\mathcal{D}}(K,G)\bigg{)}\,dG.$

Using (2.2) we obtain the classical Crofton-Herglotz formula [[1] p.75]

$\int_{G\cap
K=\emptyset}({\mathcal{D}}^{2}(K,G)-\sin^{2}{\mathcal{D}}(K,G))\,dG=2M^{2}-\frac{\pi^{3}F}{2}.$

This shows that in the presence of Crofton-Herglotz formula, theorems 1.1 and
1.2 can be seen as equivalent statements.

##  3\. Some inequalities for convex sets of constant width

In this section we will deal with compact convex sets $K$ of constant width.
For each of these sets we have the relation $R+r=a,$ where $a$ is the width of
$K$, and $r,R$ are respectively the inradius and the circumradius of $K$.
Thus, denoting $c=r/R$, we have

$r=\frac{ac}{1+c},\quad R=\frac{a}{1+c}.$

From Jung’s theorem ([3, 3.4.2]) it follows that $c\geq\sqrt{8/3}-1=0.63...$

######  Theorem 3.1.

Let $K$ be a compact convex set of constant width $a$ and $c=r/R$ the quotient
between the inradius and the cirumradius of $K$. Then

(3.1) $\displaystyle\int L(K\cap E)^{2}dE\leq
8\pi^{3}a^{3}\left(\frac{1}{(1+c)^{2}}-\frac{1}{12}\right)\,,$

which is an equality for spheres.

###### Proof.

First we observe that denoting by $p(u)$, $u\in S^{2}$, the support function
of $K$ and $\eta(u)=p(u)-a/2$ one has

$\eta(u)^{2}=p(u)^{2}+a^{2}/4-ap(u).$

Hence

$\int_{S^{2}}\eta(u)^{2}\,du=\int_{S^{2}}p(u)^{2}\,du+\pi a^{2}-2\pi a^{2}\geq
0,$

and so

(3.2) $\int_{S^{2}}p(u)^{2}\,du\geq\pi a^{2}.$

We have

$\displaystyle\int_{\mathcal{A}_{3,k}}L(K\cap
E)^{2}dE=\frac{1}{2}\int_{S^{2}}\int_{0}^{a}L(K\cap
E)^{2}\,dt\,du\leq\frac{1}{2}\int_{S^{2}}\int_{0}^{a}L(S_{R}\cap
E)^{2}\,dt\,du\leq$
$\displaystyle\frac{1}{2}\int_{S^{2}}\int_{0}^{a}\left(2\pi\sqrt{R^{2}-(p(u)-t)^{2}}\right)^{2}dt\,du=\int_{S^{2}}2\pi^{2}\left(R^{2}a-\frac{a^{3}}{3}+a^{2}p(u)-ap(u)^{2}\right)\,du$

and by (3.2)

$\displaystyle\int L(K\cap E)^{2}dE$ $\displaystyle\leq$ $\displaystyle
8\pi^{3}\left(R^{2}a-\frac{a^{3}}{12}\right)=8\pi^{3}a^{3}\left(\frac{1}{(1+c)^{2}}-\frac{1}{12}\right).$

∎

We note that by Jung’s inequality $c\geq\sqrt{8/3}-1$, the above result
implies

$\displaystyle\int L(K\cap E)^{2}dE\leq\frac{7}{3}\pi^{3}a^{3}.$

######  Proposition 3.2.

Let $K$ be a compact convex set of constant width $a$ with $c=r/R$ the
quotient between the inradius and the cirumradius of $K$. Then

$4\pi^{3}a^{3}\bigg{(}\frac{8}{3}\frac{c^{3}}{(1+c)^{3}}-\frac{1}{6}\bigg{)}\leq\int_{P\notin
K}|\Omega(P)|^{2}dP\leq
4\pi^{3}a^{3}\bigg{(}\frac{11-3c(3c^{2}+c-3)}{6(1+c)^{3}}\bigg{)},$

with equalities for spheres, where the lower bound is non negative for
$c>0.657...$.

###### Proof.

The right hand side inequality comes from (3.1) and (1.5) substituting $V$ by
$V_{r}$, where $V_{r}$ is the volume of the insphere $S_{r}$ of $K$. The left
hand side comes subtracting $4\pi^{2}V$ in the easily checked relations

$8\pi^{2}V_{r}=\int L(S_{r}\cap E)^{2}dE\leq\int L(K\cap E)^{2}dE$

and using the inequality $V\leq V_{a/2}$, where $V_{a/2}$ is the volume of the
sphere of radius $a/2$ (see [3]). ∎

######  Remark 3.3.

We note that in terms of the width only, we have

$\int_{P\notin K}|\Omega(P)|^{2}dP\leq\frac{9}{2}\pi^{3}a^{3}(\sqrt{6}-2).$

######  Remark 3.4.

One can ask if equality

$\int L(K\cap E)^{2}dE=\pi MF-4\pi^{2}V$

that holds for spheres is also true for compact convex sets of constant width.
For this case, with the same kind of arguments used above, we are only able to
prove

$\frac{c^{3}-1}{(1+c)^{3}}\leq\frac{1}{16\pi^{3}a^{3}}\left(\int L(K\cap
E)^{2}dE-(\pi
MF-4\pi^{2}V)\right)\leq\frac{-23c^{3}+3c^{2}+3c+17}{24(1+c)^{3}}$

with equalities for spheres (c=1).

## References

  * [1] W. Blaschke.  Vorlesungen uber Integralgeometrie.  VEB Deutscher Verlag der Wissenschaften, Berlin, 3rd edition, 1955.
  * [2] Minkowski H.  Ueber die Begriffe Länge.  Oberfläche und Volume 9, 1:115–121, 1901.
  * [3] H. Martini, L. Montejano, and D. Oliveros.  Bodies of Constant Width, volume 54.  Springer International Publishing, Cham, 2019.
  * [4] J. Rey Pastor and L. A. Santalo.  Geometría integral.  Espasa Calpe, 1951.
  * [5] L.A. Santaló.  Integral geometry and geometric probability.  Cambridge University Press, Cambridge, second edition, 2004.


# Non-hermitian topology in optomechanics – Supplemental material

Wojciech Brzezicki  Institute of Theoretical Physics, Jagiellonian University,
ulica S. Łojasiewicza 11, PL-30348 Kraków, Poland  International Research
Centre MagTop, Institute of Physics, Polish Academy of Sciences, Aleja
Lotników 32/46, PL-02668 Warsaw, Poland  Timo Hyart  Computational Physics
Laboratory, Physics Unit, Faculty of Engineering and Natural Sciences, Tampere
University, FI-33014 Tampere, Finland  Department of Applied Physics, Aalto
University, 00076 Aalto, Espoo, Finland  Francesco Massel  Department of
Science and Industry Systems, University of South-Eastern Norway, PO Box 235,
Kongsberg, Norway

###### Contents

  1. I Linearized optomechanical Hamiltonian: derivation
  2. II Dynamical system equations
  3. III Stability, covariance matrix, entanglement
  4. IV PBC and OBC spectra
  5. V Non-Hermitian topology
  6. VI Two-site model
  7. VII Topological protection
    1. VII.1 Hopping disorder
    2. VII.2 On-site frequency disorder
    3. VII.3 End-sites mechanical dissipation renormalization
  8. VIII Third quantization approach
    1. VIII.1 Hamiltonian
    2. VIII.2 Liouvillian
    3. VIII.3 Time-evolution
    4. VIII.4 Solution of the Lyapunov equation
    5. VIII.5 Correlation functions
    6. VIII.6 Non-zero temperature
    7. VIII.7 Stationary correlation matrix
    8. VIII.8 Entanglement in the presence of unstable end modes
  9. IX Entanglement measurement

##  I Linearized optomechanical Hamiltonian: derivation

We consider here a 1D lattice of optomechanical systems whose phononic degrees
of freedom are coupled through a hopping Hamiltonian. Namely

$H_{\rm opt}=\sum_{\rm i=0}^{N}H_{\rm i}+H_{J}$ (1)

with

$\displaystyle H_{\rm i}$ $\displaystyle=\omega_{\rm c,i}a_{\rm
i}^{\dagger}a_{\rm i}+\omega_{\rm m,i}b_{\rm i}^{\dagger}b_{\rm i}+g_{\rm
i}a_{\rm i}^{\dagger}a_{\rm i}\left(b_{\rm i}^{\dagger}+b_{\rm i}\right)$ (2)
$\displaystyle H_{J}$ $\displaystyle=-J\sum_{\rm i=0}^{N}\left(b_{\rm
i}^{\dagger}b_{\rm i+1}+\,{\rm h.c.}\right)$

where, for each site $\rm i$ $a_{\rm i}$, $b_{\rm i}$ are the photonic and
mechanical degrees of freedom, with resonant frequencies given by $\omega_{\rm
c,i}$ and $\omega_{\rm m,i}$, respectively, and $J$ the hopping energy. This
setup could be realized, for instance, in an optomechanical crystal. The
experimental realization of optomechanical crystals was first discussed in
Ref. [1], building upon previous work in the context of photonic [2, 3] and
phononic [4, 5] crystals. While the concomitant presence of hopping both for
optical and mechanical degrees of freedom can be included in our analysis,
analogously to what was discussed in Ref. [6], we focus here on the case of
purely mechanical hopping.

In our analysis, we assume that all sites are identical ($\omega_{\rm
c,i}=\omega_{\rm c}$, $\omega_{\rm m,i}=\omega_{\rm m}$, $g_{i}=g$). Along the
lines of [7], where the coupling to the environment of a 1D transmon chain was
modified to obtain a periodic modulation of the dissipation properties, we
consider here a spatial-dependent drive of the optomechanical chain in a BRRB
pattern. Each site of the chain is either driven on the blue (“B”) or the red
(“R”) sideband [8].

The linearization of the equations of motion (EOMs) generated by $H_{\rm
i}+H_{J}$, in the presence of an external drive modeled by $H_{\rm
d}=\alpha_{\rm d}e^{-i\omega_{\rm d}t}$ can be performed with the same
strategy used to linearize the Hamiltonian of a single optomechanical system.
Focusing on site $\rm i$ the EOMs for $a_{\rm i}$ and $b_{\rm i}$ can be
written as

$\displaystyle\begin{split}\dot{a}_{\rm i}=-i\omega_{\rm c}a_{\rm
i}-ig_{0}a_{\rm i}\left(b_{\rm i}^{\dagger}-\frac{\kappa}{2}a_{\rm
i}-\alpha_{\mathrm{d}}e^{-i\omega_{\rm d}t}+b_{\rm
i}\right)+\sqrt{\kappa}\,a_{\rm in,i}\end{split}$ (3a)
$\displaystyle\begin{split}\dot{b}_{\rm i}=-i\omega_{\rm m}b_{\rm
i}-ig_{0}a^{\dagger}_{\rm i}a_{\rm i}+iJ(b_{\rm i+1}+b_{\rm
i-1})-\frac{\gamma}{2}b_{\rm i}+\sqrt{\gamma}\,b_{\rm in,i}\end{split}$ (3b)

Writing the optical and mechanical operators as a combination of a strong
classical tone imposed by the drive $\alpha_{\rm d}$ and of a quantum
fluctuation. Namely,

$a_{\rm i}\to\alpha_{\rm i}+a_{\rm i},\ \ \ \\\ b_{\rm i}\to\beta_{\rm
i}+b_{\rm i}.$

We can solve Eqs. (25,26) order-by-order. The zeroth order is given by

$\displaystyle\begin{split}\dot{\alpha}_{\rm i}=i\Delta\alpha_{\rm
i}-\frac{\kappa}{2}\alpha_{\rm i}-\alpha_{\rm d}-ig_{0}\alpha_{\rm
d}\left(\beta_{\rm i}^{*}+\beta_{\rm i}\right)\end{split}$ (4a)
$\displaystyle\begin{split}\dot{\beta}_{\rm i}=-i\omega_{\rm m}\beta_{\rm
i}-\frac{\kappa}{2}\beta_{\rm i}-ig_{0}\left|\alpha_{\rm
d}\right|^{2}+iJ\left(\beta_{\rm i+1}+\beta_{\rm i-1}\right),\end{split}$ (4b)

where we have moved to a frame rotating at the external drive frequency
$\omega_{\rm d}$ and $\Delta=\Omega_{\rm d}-\Omega_{\rm c}$. Since the drive
imposes a constant displacement to the cavity and the mechanical field (in the
frame rotating at $\omega_{\rm d}$), we have $\dot{\alpha}_{\rm i}=0$ and
$\dot{\beta}_{\rm i}=0$, allowing us to solve Eqs. (4a,4b) as a set of coupled
algebraic equations. We note here that, if $g_{0}$ (and the mechanical
displacement) are sufficiently small, namely $g_{0}\left(\beta_{\rm
i}^{*}+\beta_{\rm i}\right)\ll\Delta$, the equations for the cavity field can
be solved independently from the mechanical degrees of freedom.

The first-order equations can be written as

$\displaystyle\begin{split}\dot{a}_{\rm i}=-\frac{\kappa}{2}a_{\rm
i}-ig_{0}\alpha_{\rm i}e^{-i\bar{\Delta}t}\left(e^{i\omega_{\rm m}t}b_{\rm
i}^{\dagger}+e^{-i\omega_{\rm m}t}b_{\rm i}\right)+\sqrt{\kappa}\,a_{\rm
in,i}\end{split}$ (5a) $\displaystyle\begin{split}\dot{b}_{\rm
i}=-\frac{\gamma}{2}b_{\rm i}-ig_{0}e^{i\omega_{\rm
m}t}\left(\alpha_{i}e^{-i\bar{\Delta}t}a^{\dagger}+\alpha_{\rm
i}^{*}e^{i\bar{\Delta}t}a_{\rm i}\right)+iJ\left(b_{\rm i+1}+b_{\rm
i-1}\right)+\sqrt{\gamma}\,b_{\rm in,i}\end{split}$ (5b)

where we have now moved to a frame rotating at the cavity frequency
$\omega_{\rm c}$ and at the mechanical frequency $\omega_{\rm m}$ for $a_{\rm
i}$ and $b_{\rm i}$, respectively. We now consider two specific values for the
driving frequency $\omega_{\rm d}$. The first one is such as to fulfill the
condition $\bar{\Delta}=-\omega_{\rm m}$ (red-detuned drive). Invoking the
rotating wave approximation, we are allowed to neglect in Eqs. (5a,5b) fast
oscillating terms that oscillate at $\bar{\Delta}-\omega_{\rm m}\simeq
2\omega_{\rm m}$. With this approximation, Eqs. (5a,5b) become

$\displaystyle\begin{split}\dot{a}_{\rm i}=-\frac{\kappa}{2}a_{\rm
i}-iG_{-}b_{\rm i}+\sqrt{\kappa}\,a_{\rm in,i}\end{split}$ (6a)
$\displaystyle\begin{split}\dot{b}_{\rm i}=-\frac{\gamma}{2}b_{\rm
i}-iG_{-}a_{\rm i}+iJ\left(b_{\rm i+1}+b_{\rm
i-1}\right)+\sqrt{\gamma}\,b_{\rm in,i}\end{split}$ (6b)

where and we have used the fact that $\bar{\Delta}+\omega_{\rm m}=0$ and
assumed, without loss of generality that $G_{-}=g_{0}\alpha_{i}$. Conversely,
choosing $\omega_{\rm d}$, in such a way that $\bar{\Delta}=\omega_{\rm m}$
(blue-detuned drive), leads to the following EOMs

$\displaystyle\begin{split}\dot{a}_{\rm i}=-\frac{\kappa}{2}a_{\rm
i}-iG_{+}b^{\dagger}_{\rm i}+\sqrt{\kappa}\,a_{\rm in,i}\end{split}$ (7a)
$\displaystyle\begin{split}\dot{b}_{\rm i}=-\frac{\gamma}{2}b_{\rm
i}-iG_{+}a^{\dagger}_{\rm i}+iJ\left(b_{\rm i+1}+b_{\rm
i-1}\right)+\sqrt{\gamma}\,b_{\rm in,i}.\end{split}$ (7b)

The EOMs given in Eqs. (6a-7b) can be thought of as originating from the
following bilinear Hamlitonians (in the laboratory frame)

$\displaystyle H_{\rm i,R}=\omega_{\rm c}a_{\rm i}^{\dagger}a_{\rm
i}+\omega_{\rm m}b_{\rm i}^{\dagger}b_{\rm i}+G_{\rm-}a_{\rm
i}^{\dagger}b_{\rm i}+J(b^{\dagger}_{i}b_{i+1}+b^{\dagger}_{i}b_{i-1})\,{\rm
h.c.}$ (8) $\displaystyle H_{\rm i,B}=\omega_{\rm c}a_{\rm i}^{\dagger}a_{\rm
i}+\omega_{\rm m}b_{\rm i}^{\dagger}b_{\rm i}+G_{\rm+}a_{\rm
i}^{\dagger}b^{\dagger}_{\rm
i}+J(b^{\dagger}_{i}b_{i+1}+b^{\dagger}_{i}b_{i-1})+\,{\rm h.c.}$

depending on whether site $\rm i$ is driven on the red or the blue sideband,
respectively. These two Hamiltonians are the building blocks for the
Hamiltonian given in Eq.(1) of the main text.

##  II Dynamical system equations

We consider here the dynamics induced by the Hamiltonian (26) in the presence
of a coupling to thermal reservoirs both for the mechanical and the optical
degrees of freedom, introducing dissipation rates $\gamma$ and $\kappa$ and
noise operators $b_{\rm in,i}$ and $a_{\rm in,i}$ for the mechanical and
optical degrees of freedom, respectively. The operators $b_{\rm in,i}$ fulfill
the following relations

$\displaystyle\langle b_{\rm in,i}(t)b^{\dagger}_{\rm
in,i}(t^{\prime})\rangle=\left(n_{\rm m}+1\right)\delta(t-t^{\prime})$ (9a)
$\displaystyle\langle b^{\dagger}_{\rm in,i}(t)b_{\rm
in,i}(t^{\prime})\rangle=n_{\rm m}\delta(t-t^{\prime})$ (9b)

and analogous relations hold for $a_{\rm in,i}$

$\displaystyle\begin{split}\langle a_{\rm in,i}(t)a^{\dagger}_{\rm
in,i}(t^{\prime})\rangle=\left(n_{\rm
c}+1\right)\delta(t-t^{\prime})\end{split}$ (10a)
$\displaystyle\begin{split}\langle a^{\dagger}_{\rm in,i}(t)a_{\rm
in,i}(t^{\prime})\rangle=n_{\rm c}\delta(t-t^{\prime}).\end{split}$ (10b)

The Hamiltonian given Eq. (1) of the main text generates the following
dynamics, described by the quantum Langevin equations (QLEs) [9], for the
optical and mechanical degrees of freedom

$\displaystyle\begin{split}\dot{a}_{\rm i}&=-i\left[a_{\rm
i},\hat{H}\right]-\frac{\kappa}{2}a_{\rm i}+\sqrt{\kappa}a_{\rm
i,in}\end{split}$ (11a) $\displaystyle\begin{split}\dot{b}_{\rm
i}&=-i\left[b_{\rm i},\hat{H}\right]-\frac{\gamma}{2}b_{\rm
i}+\sqrt{\gamma}b_{\rm i,in}.\end{split}$ (11b)

Eqs.( 11a,11b) can be written more compactly in terms of a operators vector
$\mathbf{a}=\left[a_{1},b_{1},\dots,a_{4N},b_{4N};a^{\dagger}_{1},b^{\dagger}_{1},\dots
a^{\dagger}_{4N},b^{\dagger}_{4N}\right]^{T}$ and an input-operators vector
$\mathbf{a}_{\rm in}=\left[a_{\rm in\,1},b_{\rm in\,1},\dots,a_{\rm
in\,4N},b_{\rm in\,4N};a^{\dagger}_{\rm in\,1},b^{\dagger}_{\rm
in\,1}\dots,\dots,a^{\dagger}_{\rm in\,4N},b^{\dagger}_{\rm
in\,4N}\right]^{T}$

$\dot{\mathbf{a}}=-i\mathcal{H}_{\rm
NH}\mathbf{a}+\sqrt{\bar{\eta}}\mathbf{a}_{\rm in}$ (12)

with $\bar{\eta}={\rm
diag}\left(\kappa,\gamma,\dots;\kappa,\gamma,\dots\kappa,\gamma\right)$, and

$\mathcal{H}_{\rm NH}=\bar{\sigma}_{z}\mathcal{H}-\frac{i\bar{\eta}}{2}$ (13)

where $\bar{\sigma}_{z}={\rm diag}\left(1,\dots,1;-1,\dots,-1\right)$. The
solution of Eq. (12) can be written as

$\mathbf{a}=\exp\left[-i\mathcal{H}_{\rm
NH}t\right]\mathbf{a}_{0}+\int_{0}^{t}\exp\left[-i\mathcal{H}_{\rm
NH}(t-\tau)\right]\sqrt{\bar{\eta}}\mathbf{a}_{\rm in}d\tau,$ (14)

which corresponds to Eq. (2) of the main text.

##  III Stability, covariance matrix, entanglement

The first aspect in the study of the (quantum) dynamical system described by
Eq. (12), whose solution is given by Eq. (14) (Eq. (2) of the main text) is
its stability. Eq. (12) can be construed as a quantum stochastic differential
equation [9] (quantum version of a Ornstein-Uhlenbeck process [10]), which can
be written, in general, as

$\dot{\bar{a}}=A\bar{a}+\bar{\xi}_{\rm in}.$ (15)

Its stability is guaranteed by the condition ${\rm Re}\,\lambda_{\rm i}<0$
where $\lambda_{\rm i}$ are the eigenvalues of the matrix associated with the
dynamical system, in our case $A=-i\mathcal{H}_{\rm NH}$. This condition
corresponds to the ${\rm Im}\,E<0$ condition on the eigenenergies of the non-
Hermitian Hamiltonian $\mathcal{H}_{\rm NH}$.

From Eq. (15), it is possible to write the Lyapunov equation for the time
evolution of the covariance matrix $V$ for the operators
$\left[a_{1},b_{1},a_{1}^{\dagger},b_{1}^{\dagger},\dots\right]$, namely

$\frac{dV_{\bar{a},\bar{a}^{\dagger}}}{dt}=AV_{\bar{a},\bar{a}^{\dagger}}+V_{\bar{a},\bar{a}^{\dagger}}A^{\dagger}+D$
(16)

where

$V_{\bar{a},\bar{a}^{\dagger}}=\begin{bmatrix}\langle
a_{1}a_{1}^{\dagger}\rangle&\langle a_{1}b_{1}^{\dagger}\rangle&\langle
a_{1}a_{1}\rangle&\langle a_{1}b_{1}\rangle&\cdots\\\ \langle
b_{1}a_{1}^{\dagger}\rangle&\ddots&&\\\ \vdots&&&\\\ \end{bmatrix}$

and $D={\rm
diag}\left[\kappa\left(n_{a}+1\right),\gamma\left(n_{b}+1\right),\kappa
n_{a},\gamma n_{b},\dots\right]$. In our analysis, we are interested in the
stationary value of the covariance matrix. We therefore seek a numerical
solution to the equation

$AV_{\bar{a},\bar{a}^{\dagger}}+V_{\bar{a},\bar{a}^{\dagger}}A^{\dagger}+D=0.$
(17)

From $V_{\bar{a},\bar{a}^{\dagger}}$ it is possible to obtain the covariance
matrix between the quadratures –either mechanical

$\displaystyle X_{\rm m,i}=1/\sqrt{2}\left(b_{\rm i}^{\dagger}+b_{\rm
i}\right),$ $\displaystyle P_{\rm m,i}=i/\sqrt{2}\left(b_{\rm
i}^{\dagger}-b_{\rm i}\right)$

or optical

$\displaystyle X_{\rm o,i}=1/\sqrt{2}\left(a_{\rm i}^{\dagger}+a_{\rm
i}\right),$ $\displaystyle P_{\rm o,i}=i/\sqrt{2}\left(a_{\rm
i}^{\dagger}-a_{\rm i}\right).$

We evaluate the entanglement properties between the quadratures $\xi,\eta={\rm
m,c}$ on sites $\rm i,j$, considering the two mode covariance matrix

$V_{\rm\xi,i,\eta,j}=\begin{bmatrix}\langle
X_{\rm\xi,i}X_{\rm\xi,i}\rangle&\langle
X_{\rm\xi,i}P_{\rm\xi,i}\rangle&\langle
X_{\rm\xi,i}X_{\rm\eta,j}\rangle&\langle X_{\rm\xi,i}P_{\rm\eta,j}\rangle\\\
\langle P_{\rm\xi,i}X_{\rm\xi,i}\rangle&\langle
P_{\rm\xi,i}P_{\rm\xi,i}\rangle&\langle
P_{\rm\xi,i}X_{\rm\eta,j}\rangle&\langle P_{\rm\xi,i}P_{\rm\eta,j}\rangle\\\
\langle X_{\rm\eta,i}X_{\rm\xi,i}\rangle&\langle
X_{\rm\eta,i}P_{\rm\xi,i}\rangle&\langle
X_{\rm\eta,i}X_{\rm\eta,j}\rangle&\langle X_{\rm\eta,i}P_{\rm\eta,j}\rangle\\\
\langle P_{\rm\eta,i}X_{\rm\xi,i}\rangle&\langle
P_{\rm\eta,i}P_{\rm\xi,i}\rangle&\langle
P_{\rm\eta,i}X_{\rm\eta,j}\rangle&\langle
P_{\rm\eta,i}P_{\rm\eta,j}\rangle\end{bmatrix}.$ (18)

Eq. (18) can be written in terms of three $2\times 2$ matrices $\alpha$,
$\beta$, $\gamma$ as

$V_{\rm i,j}=\left[\begin{array}[]{@{}c|c@{}}\alpha_{\rm i,j}&\gamma_{\rm
i,j}\\\ \hline\cr\gamma_{\rm i,j}^{T}&\beta_{\rm i,j}\end{array}\right]$ (19)

Following [11], we can determine the negativity associated with the
entanglement (or lack thereof).

The positivity of the partial transposed state (PPT criterion) can be recast
as the condition $\tilde{\nu}_{-}\geq 1$ on the smallest symplectic eigenvalue
of the covariance matrix evaluated on the partially transposed state. These
can be expressed in terms of Eq. (19) as

$\tilde{\nu}_{\pm}=\sqrt{\frac{\tilde{\Delta}_{\rm
i,j}\pm\sqrt{\tilde{\Delta}_{\rm i,j}-4{\,\rm Det\,}V_{\rm i,j}}}{2}}$ (20)

where $\tilde{\Delta}_{\rm i,j}={\rm Det\,}\alpha_{\rm i,j}+{\rm
Det\,}\beta_{\rm i,j}-2{\,\rm Det\,}\gamma_{\rm i,j}$.

The (logarithmic) negativity

$E_{\mathcal{N}}={\rm max}\left[0,-\log_{2}\tilde{\nu}_{-}\right]$ (21)

gives a measure of the violation of the condition $\tilde{\nu}_{-}\geq 1$, and
thus quantifies the entanglement between the modes localized on sites $\rm i$
and $\rm j$. Furthermore, the diagonal elements of $V_{\rm\xi,i,\eta,j}$ allow
us to infer the number of mechanical excitations on sites $\rm i$ and $\rm j$.

##  IV PBC and OBC spectra

We depict here the spectra for periodic and open boundary conditions for
different values of $G_{+}$. In Fig. 1, we consider the case where both end
states are stable. Panel 1(a) illustrates the case for which ${\rm Im}E_{\rm
e}$ is within the ${\rm Im}E_{\rm b}$ band ($G_{+}=0.2$), whereas panel 1(b)
refers to the case $G_{+}=0.446$, for which a gap develops between ${\rm
Im}E_{\rm e}\to 0^{-}$ and ${\rm Im}E_{\rm b}<0$ states.

Figure 1: Real and imaginary spectra for periodic (first columm) and open
(second column) boundary conditions (columns three and four represent a zoom-
in of column two). Red dots correspond to ${\rm Im}\,E_{\rm e}<0$, i.e.
(stable) end states, blue dots to ${\rm Im}E_{\rm b}$, i.e. bulk states. (a)
$G_{+}=0.05$, (b) $G_{+}=0.242$. $G_{-}=1,J=0.5$, $\gamma=10^{-4}$ (all
energies expressed in units of $\kappa$). For $G_{+}=0.242$, ${\rm Im}\,E_{\rm
e}\simeq 4\cdot 10^{-4}$.

In Fig. 2, we further increase $G_{+}$. At first ($G_{+}=0.26$, panel 2(a))
end states become unstable, while bulk states are still stable. Further
increasing $G_{+}$ ($G_{+}=0.35$,$G_{+}=1$) leads the instability of the bulk
states as well. The case $G_{+}=0.35$ (panel 2(b)) corresponds to the case for
which both the PBC and the OBC spectra are gapless at ${\rm Im}E=0$, while a
gap develops for $G_{+}=1$ (anel 2(c))

Figure 2: Real and imaginary spectra for periodic (first columm) and open
(second column) boundary conditions (columns three and four represent a zoom-
in of column two). Red dots correspond to ${\rm Im}E_{\rm e}<0$, i.e. stable
end states, black dots correspond to ${\rm Im}E_{\rm e}>0$, i.e. unstable end
states, blue dots to ${\rm Im}E_{\rm b}$, i.e. bulk states. (a) $G_{+}=0.35$,
(b) $G_{+}=1$. $G_{-}=1,J=0.5$, $\gamma=10^{-4}$ (all energies expressed in
units of $\kappa$).

##  V Non-Hermitian topology

In order to establish the topological character of the (complex) spectrum of
the chain, we rely on the analysis performed in Ref. [12], where, as outlined
in the main text, the authors establish an exact mapping between the complex
spectrum at ${\rm Re}\,E=0$ of a non-hermitian Hamiltonian –in our case
$\mathcal{H}_{\rm NH}$– and an effective hermitian Hamiltonian given by

${\mathcal{H}}_{{\rm eff}}(\eta,k)=\eta S-iS{\mathcal{H}}_{{\rm NH}}(k).$

The definition of $S$ originates from the symmetry properties of
$\mathcal{H}_{\rm NH}$: it is possible to show that $\mathcal{H}_{\rm NH}$
obeys a non-hermitian chiral symmetry

$S{\cal H}_{{\rm NH}}S^{-1}=-{\cal H}_{{\rm NH}}^{\dagger}$ (22)

where $S=\Pi\,\Sigma$,

$\Sigma=\begin{pmatrix}\sigma&0&0&0\\\ 0&\sigma&0&0\\\ 0&0&\sigma&0\\\
0&0&0&\sigma\end{pmatrix},\quad\Pi=\begin{pmatrix}\sigma&0&0&0\\\
0&\tau&0&0\\\ 0&0&-\tau&0\\\ 0&0&0&-\sigma\end{pmatrix},$

$\sigma=\sigma_{z}\otimes\mathbb{I}_{2}$,$\tau=\mathbb{I}_{2}\otimes\sigma_{z}$,
and $\sigma_{z}$ is the 2x2 Pauli matrix.

We focus on the eigenstates $\psi$ of $\mathcal{H}_{\rm NH}$ with ${\rm Re}=0$
(for open boundary conditions). The states $\psi$ must therefore satisfy the
following relation

$\mathcal{H}_{\rm NH}\psi=i\eta\,\psi.$ (23)

From the symmetry properties expressed in Eq. (22), it is possible to show
that Eq. (23) can be recast as a zero-energy eigenproblem for the _hermitian_
Hamiltonian $\mathcal{H}_{\rm eff}$. Namely

$\mathcal{H}_{\rm eff}(\psi,k)\psi\doteq S\left(\eta-i\mathcal{H}_{\rm
NH}\right)\psi=0.$ (24)

The potentially nontrivial topological properties of $\mathcal{H}_{\rm
eff}(k,\eta)$ are reflected into the topological protection of the ${\rm
Re}=0$ (end) modes of $\mathcal{H}_{\rm NH}$. To ascertain the topological
protection of the end modes, following [7], we can evaluate the Berry
curvature and the Chern number associated with ${\mathcal{H}}_{{\rm eff}}$ as

$\Omega_{k,\eta}=\sum_{{n\leq 8\atop m>8}}{\rm
Im}\frac{2\left[\psi_{k,\eta}^{n\,*}\,\left(\partial_{k}\mathcal{H}^{{\rm
eff}}\right)\,\psi_{k,\eta}^{m}\right]\left[\psi_{k,\eta}^{m\,*}\,\left(\partial_{\eta}\mathcal{H}^{{\rm
eff}}\right)\,\psi_{k,\eta}^{n}\right]}{\left(E_{k,\eta}^{(n)}-E_{k,\eta}^{(m)}\right)^{2}}$
(25)

and

$C=\frac{1}{2\pi}\int_{-\infty}^{+\infty}d\eta\int_{0}^{2\pi}dk\
\Omega_{k,\eta},$ (26)

where the eigenstates are sorted in ascending order of the eigenenergy. In
order to guarantee the quantization of the Chern number for $\mathcal{H}_{\rm
eff}(k,\eta)$, one has to require its periodicity in $k$ and $\eta$. While the
periodicity of the former is guaranteed by the periodicity of
$\mathcal{H}_{\rm NH}(k)$, $\mathcal{H}_{\rm eff}(k,\eta)$ is not periodic in
$\eta$. This issue can be solved by defining a compactified version of
$\mathcal{H}_{\rm eff}(k,\eta)$, sharing with it the same spectrum, through
the relation

$\mathcal{H}_{\rm eff}^{cmp}(k,\eta)=\mathcal{R}_{\eta}\mathcal{H}_{\rm
eff}(k,\eta)\mathcal{R}_{\eta}^{\dagger}$ (27)

with

$\mathcal{R}_{\eta}=\exp\left[i\frac{\pi}{4}\left(1+\tanh\eta\right)\mathcal{G}\right],\qquad\mathcal{G}=\begin{pmatrix}&&&\mathbb{I}_{4}\\\
&&\mathbb{I}_{4}&\\\ &\mathbb{I}_{4}&&\\\ \mathbb{I}_{4}&&&\\\ \end{pmatrix}$
(28)

.

##  VI Two-site model

The stability and dissipation properties of the end states can be approximated
terms of dynamics of a two-site system, with one site driven on the blue
sideband ($G_{+}$) and the other on the red ($G_{-}$). Intuitively, this setup
corresponds to $G_{-}$-driven site acting as an effective bath for the $G_{+}$
driven one, stabilizing its dynamics. In Fourier space, the dynamics of the
two-site system can be written as

${\bf a}=A\,\xi_{\rm in}$ (29)

where ${\bf a}=\left[a_{1},b^{\dagger}_{1},a_{2},b_{2}\right]$, $\xi_{\rm
in}=\left[\sqrt{\kappa}\,a_{\rm in,1},\sqrt{\gamma}\,b^{\dagger}_{\rm
in,1},\sqrt{\kappa}\,a_{\rm in,2},\sqrt{\gamma}\,b_{\rm in,2}\right]$ and

$A^{-1}=\begin{pmatrix}\chi_{\rm c}^{-1}&iG_{+}&0&0\\\ -iG_{+}&\chi_{\rm
m}^{-1}&0&-iJ\\\ 0&0&\chi_{\rm c}^{-1}&-iG_{-}\\\ 0&-iJ&-iG_{-}&\chi_{\rm
m}^{-1}\\\ \end{pmatrix}$ (30)

with $\chi_{\rm m}=\left(\gamma/2-i\omega\right)^{-1}$ and $\chi_{\rm
c}=\left(\kappa/2-i\omega\right)^{-1}$. With an analysis similar to the one
carried out in Sec. II, it is possible to show that the poles of the matrix
$A$ correspond to the (complex) excitation energies for the two-site system,
and, in particular, their imaginary part characterizes its stability
properties. In Fig. 3 we compare the complex eigenenergies of the two-site
system, with the spectrum of the chain, and we see that the overall behavior
of the end state (in red) is well captured by the two-site description. While
an analytical solution to the two-site problem is possible the expressions for
${\rm Im}E_{\rm 2-site}$ are not particularly informative. We provide here the
expression for ${\rm Im}E_{\rm 2-site}$ corresponding to the largest imaginary
part of the two-site eigenenergies in the limit $G_{+}\to 0$ and
$G_{+}\to\infty$. For $G_{+}\to 0$ we have

${\rm Im}\,E_{\rm
2-site}\simeq-\frac{\Delta^{1/3}+\kappa+2\gamma}{6}+\frac{12(G_{-}^{2}+J^{2})-\left(\kappa-\gamma\right)^{2}}{6\Delta^{1/3}}$
(31)

with

$\displaystyle\Delta$
$\displaystyle=\kappa\left(\kappa^{2}-18G_{-}^{2}+36J^{2}\right)+\sqrt{\kappa^{2}\left(\kappa^{2}-18G_{-}^{2}+36J^{2}\right)^{2}-\left(\kappa^{2}-12G_{-}^{2}-12J^{2}\right)^{3}}$
$\displaystyle-3\left(\kappa^{2}-6G_{-}^{2}+12J^{2}\right)\gamma+O(\gamma^{2})$

which, expanding in powers of $J/G_{-}$, can be written as

${\rm Im}\,E_{\rm
2-site}\simeq-\frac{\gamma}{2}-\left(\kappa-\gamma\right)\frac{J^{2}}{G_{-}^{2}}\left(1-\frac{J^{2}}{G_{-}^{2}}\right)+O\left[\left(\frac{J}{G_{-}}\right)^{5}\right]$
(32)

Conversely, for $G_{+}\to\infty$ we get

$\displaystyle{\rm Im}E_{\rm
2-site}\simeq-\frac{\kappa+\gamma}{2}+\sqrt{\frac{(\kappa-\gamma)^{2}}{4}+4G_{+}^{2}}.$
(33)

Figure 3: Imaginary part of the OBC spectrum ${\rm Im}E$ as a function of
$G_{+}$. Red line: ${\rm Im}\,E_{\rm e}$ ; dashed line: approximate result
given by the solution of the two-site problem. (a)
$G_{+}\in\left[0,0.4\right]$ corresponding to Fig. 3 of the main text. (b)
$G_{+}\in\left[0,2\right]$. $G_{-}$,$J=0.5$, $\gamma<=10^{-4}$ (all energies
expressed in units of $\kappa$).

##  VII Topological protection

###  VII.1 Hopping disorder

We plot here the spectra (for periodic- and open-boundary conditions ) in the
presence of disorder for the hopping parameter $J$ ($\delta J$ unformly
distributed in $\left[-0.1J,0.1J\right]$), for different values of $G_{+}$.
For all values of $G_{+}$ it is possible to see that the ${\rm Re}\,E=0$ end
states survive in the presence of disorder. The degeneracy of ${\rm Im}\,E$,
however, is partially lifted: the imaginary part of the spectrum shows the
appearance of pairs of states instead of the quadruplet associated with the
corresponding spectrum in the absence of disorder.

Figure 4: Real and imaginary spectra for periodic (first columm) and open
(second column) boundary conditions (columns three and four represent a zoom-
in of column two). Red dots correspond to ${\rm Im}E_{\rm e}<0$, i.e. stable
end states, black dots correspond to ${\rm Im}E_{\rm e}>0$, i.e. unstable end
states, blue dots to ${\rm Im}E_{\rm b}$, i.e. bulk states. (a) $G_{+}=0.05$,
(b) $G_{+}=0.242$, (c) $G_{+}=1$. $G_{-}=1,J=0.5$, $\gamma=10^{-4}$, $\delta
J/J=0.1$ (all energies expressed in units of $\kappa$).

###  VII.2 On-site frequency disorder

We take here into account the potential difference of on-site frequencies for
either the mechanics or the cavity. Following the derivation in section I it
is easy to show that a site dependence of $\omega_{\mathrm{c}}$ can be
absorbed into the choice of the driving frequency $\omega_{\mathrm{d}}$.
Conversely, a potential site dependence of the mechanical frequency has to be
explicitly taken into account. In Fig. 5, we have plotted the spectra for
$\omega_{m}$ uniformly distributed (in the rotating frame) in
$\left[-0.01\omega_{\mathrm{m}},0.01\omega_{\mathrm{m}}\right]$. The disorder
amplitude is compatible with the current experimental capabilities discussed
in [13].

Figure 5: Real and imaginary spectra for periodic (first columm) and open
(second column) boundary conditions (columns three and four represent a zoom-
in of column two). Red dots correspond to ${\rm Im}E_{\rm e}<0$, i.e. stable
end states, black dots correspond to ${\rm Im}E_{\rm e}>0$, i.e. unstable end
states, blue dots to ${\rm Im}E_{\rm b}$, i.e. bulk states. (a) $G_{+}=0.05$,
(b) $G_{+}=0.242$, (c) $G_{+}=1$. $G_{-}=1,J=0.5$, $\gamma=10^{-4}$,
$\delta\omega_{\mathrm{m}}/\omega_{\mathrm{m}}=0.01$ (all energies expressed
in units of $\kappa$).

As expected from the breaking of the chiral symmetry, Fig 5 displays a
deviation from the ${\rm Re}\,E=0$ condition for the end states. For the
parameters considered here, however, the real part of the spectrum still
exhibits a significant gap for the ${\rm Re}\,E\neq 0$ states.

###  VII.3 End-sites mechanical dissipation renormalization

We depict here the OBC spectra in the presence of measurement-induced
renormalization of the end-states mechanical dissipation. As expected, in
these figures it is possible to see that even for $\gamma_{\rm
end}=100\gamma$, the spectrum is essentially unchanged, as expected from a
topological protection argument.

Figure 6: Real and imaginary spectra for periodic (first columm) and open
(second column) boundary conditions (columns three and four represent a zoom-
in of column two). Red dots correspond to ${\rm Im}E_{\rm e}<0$, i.e. stable
end states, black dots correspond to ${\rm Im}E_{\rm e}>0$, i.e. unstable end
states, blue dots to ${\rm Im}E_{\rm b}$, i.e. bulk states. (a) $G_{+}=0.05$,
(b) $G_{+}=0.242$, $\gamma=10^{-4}$, $\gamma_{\rm end}=10^{-2}$ (all energies
expressed in units of $\kappa$).

##  VIII Third quantization approach

###  VIII.1 Hamiltonian

The real-space Hamiltonian has a form of:

$\displaystyle{\cal H}$ $\displaystyle=$
$\displaystyle\sum_{i}\left(\vec{d}_{i}^{\,\dagger}\hat{P}\vec{d}_{i}+\vec{d}_{i}\hat{Q}\vec{d}_{i}+\vec{d}_{i}^{\,\dagger}\hat{Q}^{\star}\vec{d}_{i}^{\,\dagger}\right)+\sum_{i}\left(\vec{d}_{i}^{\,\dagger}\hat{R}\vec{d}_{i+1}+\vec{d}_{i+1}^{\,\dagger}\hat{R}^{\dagger}\vec{d}_{i}\right)$
(34)

with $\vec{d}_{i}$ describing bosonic degrees of freedom within unit cell $i$
as a vector:

$\displaystyle\vec{d}_{i}=\begin{pmatrix}a_{1,i}&b_{1,i}&a_{2,i}&b_{2,i}&a_{3,i}&b_{3,i}&a_{4,i}&b_{4,i}\end{pmatrix}$
(35)

with matrices $\hat{P}$, $\hat{Q}$ and $\hat{R}$ describing hopping and
pairing of the bosons inside and between unit cells. We have:

$\hat{P}=\begin{pmatrix}0&0&0&0&0&0&0&0\\\ 0&0&0&-1&0&0&0&0\\\
0&0&0&G_{-}&0&0&0&0\\\ 0&-1&G_{-}&0&0&-1&0&0\\\ 0&0&0&0&0&G_{-}&0&0\\\
0&0&0&-1&G_{-}&0&0&-1\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&0&0&-1&0&0\end{pmatrix}\\!,\,\,\hat{Q}=\begin{pmatrix}0&G_{+}&0&0&0&0&0&0\\\
G_{+}&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&G_{+}\\\
0&0&0&0&0&0&G_{+}&0\end{pmatrix}\\!,\,\,\hat{R}=\begin{pmatrix}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&-1&0&0&0&0&0&0\end{pmatrix}\\!.$ (36)

For a system consisting of $N$ unit cells we can set we define a composite
vector:

$\displaystyle\vec{\underline{d}}=\begin{pmatrix}\vec{d}_{1}&\vec{d}_{2}&\dots&\vec{d}_{N}\end{pmatrix}$
(37)

and the Hamiltonian becomes:

$\displaystyle{\cal H}$ $\displaystyle=$
$\displaystyle\vec{\underline{d}}^{\dagger}\left(\hat{T}_{0}\otimes\hat{P}+\hat{T}_{1}\otimes\hat{R}+\hat{T}_{-1}\otimes\hat{R}^{\dagger}\right)\vec{\underline{d}}+\vec{\underline{d}}\left(\hat{T}_{0}\otimes\hat{Q}\right)\vec{\underline{d}}+\vec{\underline{d}}^{\dagger}\left(\hat{T}_{0}\otimes\hat{Q}^{\star}\right)\vec{\underline{d}}^{\dagger}$
(38) $\displaystyle=$
$\displaystyle\vec{\underline{d}}^{\dagger}\hat{H}\vec{\underline{d}}+\vec{\underline{d}}\hat{K}\vec{\underline{d}}+\vec{\underline{d}}^{\dagger}\hat{K}^{\star}\vec{\underline{d}}^{\dagger}$

with $\hat{T}_{0}$ being $N\times N$ identity matrix and $\hat{T}_{1}$ being a
matrix describing shift by one of all unit cells, i.e.,

$\hat{T}_{1}=\begin{pmatrix}0&1\\\ &0&1\\\ &&\ddots&\ddots\\\ &&&0&1\\\
x&&&&0\end{pmatrix},\quad\hat{T}_{-1}:=\hat{T}_{1}^{T}.$ (39)

Here $x=1$ for periodic and $x=0$ for open boundary conditions.

###  VIII.2 Liouvillian

Following the third quantization procedure by Prosen and Seligman [14, 7] we
find that the Liouvillian can be written in the form of

${\cal
L}=-2\begin{pmatrix}\vec{\underline{d}}_{0}^{\,\prime}&\vec{\underline{d}}_{1}^{\,\prime}\end{pmatrix}\hat{X}^{T}\begin{pmatrix}\vec{\underline{d}}_{0}&\vec{\underline{d}}_{1}\end{pmatrix}+\begin{pmatrix}\vec{\underline{d}}_{0}^{\,\prime}&\vec{\underline{d}}_{1}^{\,\prime}\end{pmatrix}\hat{Y}\begin{pmatrix}\vec{\underline{d}}_{0}^{\,\prime}&\vec{\underline{d}}_{1}^{\,\prime}\end{pmatrix}$
(40)

with operators $d_{i,\mu}$ satisfying almost canonical commutation relations
for bosons:

$[d_{i,\mu},d_{j,\nu}^{\prime}]=\delta_{ij}\delta_{\mu\nu}$ (41)

where $d_{j,\nu}^{\prime}$ is not Hermitian conjugate of $d_{j,\nu}$. In the
Liouvillan the matrices defining Hamiltonian enter as blocks of $\hat{X}$ and
$\hat{Y}$:

$\hat{X}=\frac{1}{2}\begin{pmatrix}i\hat{H}^{\star}+\hat{M}&-2i\hat{K}\\\
2i\hat{K}^{\star}&-i\hat{H}+\hat{M}^{\star}\end{pmatrix},\,\hat{Y}=\begin{pmatrix}-i\hat{K}^{\star}&0\\\
0&i\hat{K}\end{pmatrix},$ (42)

where $\hat{M}$ is matrix describing dissipation. In our case it takes a
diagonal form of:

$\displaystyle\hat{M}=\hat{T}_{0}\otimes\frac{1}{2}{\rm
diag}\begin{pmatrix}\kappa&\gamma&\kappa&\gamma&\kappa&\gamma&\kappa&\gamma\end{pmatrix}.$
(43)

Liouvillan is a quadratic form of boson-like operators and can be brought to a
diagonal form of

${\cal L}=-2\sum_{i=1}^{16N}\beta_{i}\xi_{i}^{\prime}\xi_{i},$ (44)

by a suitable symplectic transformation. One finds that the rapidities
$\beta_{i}$ are the eigenvalues of the matrix $\hat{X}$ so that one can
identify $X$ to be proportional to the non-Hermitian Hamiltonian of the
system. Strictly speaking we can identify:

$\hat{H}_{\rm NH}=-2i\hat{X}.$ (45)

###  VIII.3 Time-evolution

To maximally simplify the notation let’s write the Liouvillan as

${\cal
L}=\vec{\alpha}^{\prime}\hat{A}\vec{\alpha}+\vec{\alpha}^{\prime}\hat{B}\vec{\alpha}^{\prime},$
(46)

with $a_{i}$ and $a_{i}^{\prime}$ being almost canonical bosons and $\hat{A}$,
$\hat{B}$ being $L\times L$ matrices. Now, the third-quantized Lindblad
equation reads:

$\frac{d}{dt}\rho(t)={\cal L}[\rho]$ (47)

thus the time-evolution of the density matrix is given by

$\left|\rho(t)\right\rangle=\exp\left[t{\cal
L}\right]\left|\rho(0)\right\rangle.$ (48)

The solution has a form of

$\left|\rho(t)\right\rangle={\cal T}\exp\left[\int_{0}^{t}dt^{\prime}{\cal
L}(t^{\prime})\right]\rho(0).$ (49)

Taking a coherent state as initial

$\left|\rho(0)\right\rangle=\left|\vec{z}(0)\right\rangle=\exp[\vec{z}(0)\vec{\alpha}^{\,\prime}]\left|0\right\rangle,$
(50)

using Suzuki-Trotter decomposition and resolution of identity in coherent
states we get:

$\left|\rho(t)\right\rangle=\frac{1}{\pi^{2LN_{t}}}\int\left(\prod_{i=1}^{N_{t}}d\Im\vec{z}_{i}d\Re\vec{z}_{i}\right)\exp\left[\sum_{i=1}^{N_{t}}S_{i}\right]\left|\vec{z}_{1}\right\rangle,$
(51)

with

$S_{i}=\vec{z}_{i}^{\,\star}\left[1+\frac{t}{N_{t}}\hat{A}\right]\vec{z}_{i+1}+\vec{z}_{i}^{\,\star}\left[\frac{t}{N_{t}}\hat{B}\right]\vec{z}_{i}^{\,\star}-\vec{z}_{i}^{\,\star}\vec{z}_{i}$
(52)

and $\vec{z}_{N_{t}+1}\equiv\vec{z}(0)$ for number of time-steps
$N_{t}\to\infty$. By calculating multi-dimensional Gaussian integral we can
get:

$\left|\rho(t)\right\rangle\\!=\\!\\!\frac{1}{\pi^{2L}}\\!\int\\!d\Im\vec{z}d\Re\vec{z}\exp\\!\left[-\vec{z}^{\,\star}\vec{z}+\vec{z}^{\\!\star}\exp\left[t\hat{A}\right]\vec{z}(0)+\frac{t}{N_{t}}\vec{z}^{\,\star}\left[\hat{B}+\hat{h}\hat{B}\hat{h}^{T}+\dots+\hat{h}^{N-1}\hat{B}(\hat{h}^{T})^{N-1}\right]\vec{z}^{\,\star}\right]\left|\vec{z}\right\rangle.$
(53)

with

$\hat{h}=1+\frac{t}{N}\hat{A}.$ (54)

The series in the exponent can be summed up in a similar way one is summing up
a finite geometric series. Then we get for $N_{t}\to\infty$:

$\displaystyle\left|\rho(t)\right\rangle=\frac{1}{\pi^{2L}}\int
d\Im\vec{z}d\Re\vec{z}\exp\left[-\vec{z}^{\,\star}\vec{z}+\vec{z}^{\,\star}\exp\left[t\hat{A}\right]\vec{z}(0)+\vec{z}^{\,\star}\hat{C}(t)\vec{z}^{\,\star}\right]\left|\vec{z}\right\rangle$
(55)

with $\hat{C}(t)$ being a matrix that satisfy Lyapunov equation of the form

$\hat{A}\hat{C}(t)+\hat{C}(t)\hat{A}^{T}=\exp\left[t\hat{A}\right]\hat{B}\exp\left[t\hat{A}^{T}\right]-\hat{B}.$
(56)

After performing final Gaussian integration the time-evolved density matrix
becomes:

$\left|\rho(t)\right\rangle=\exp\left[\vec{\alpha}^{\prime}\hat{C}(t)\vec{\alpha}^{\prime}\right]\left|\exp\left[t\hat{A}\right]\vec{z}(0)\right\rangle.$
(57)

###  VIII.4 Solution of the Lyapunov equation

The Lyapunov equation is known to have a solution of the closed form:

$\hat{C}(t)=\int_{0}^{\;t}\exp\left[x\hat{A}\right]\hat{B}\exp\left[x\hat{A}^{T}\right]dx.$
(58)

Assuming that $\hat{A}$ can be diagonalized in the eigenbasis $\hat{\beta}$ we
get:

$\hat{\Lambda}=\hat{\beta}^{-1}\hat{A}\hat{\beta}={\rm
diag}\begin{pmatrix}\lambda_{1}&\lambda_{2}&\dots&\lambda_{2L}\end{pmatrix}=\hat{\beta}^{T}\hat{A}^{T}\hat{\beta}^{-1T}.$
(59)

Thus:

$\hat{\beta}^{-1}\hat{C}(t)\hat{\beta}^{-1T}=\int_{0}^{\;t}\exp\left[x\hat{\Lambda}\right]\hat{\beta}^{-1}\hat{B}\hat{\beta}^{-1T}\exp\left[x\hat{\Lambda}\right]dx$
(60)

or:

$\displaystyle\left(\hat{\beta}^{-1}\hat{C}(t)\hat{\beta}^{-1T}\right)_{ij}=\left(\hat{\beta}^{-1}\hat{B}\hat{\beta}^{-1T}\right)_{ij}$
$\displaystyle\times\int_{0}^{\;t}dx\exp\left[x(\lambda_{i}+\lambda_{j})\right].$
(61)

Then by explicitly calculating the integral we finally get:

$\hat{C}(t)=\hat{\beta}\left[\left(\hat{\beta}^{-1}\hat{B}\hat{\beta}^{-1T}\right)\circ\hat{\Gamma}(t)\right]\hat{\beta}^{T},$
(62)

with

$\Gamma_{ij}(t)=\frac{\exp\left[t\left(\lambda_{i}+\lambda_{j}\right)\right]-1}{\lambda_{i}+\lambda_{j}}$
(63)

and $\circ$ denoting Hadamard product of matrices.

###  VIII.5 Correlation functions

After starting from the Liouvillian of the general third-quantized form of Eq.
(46) and getting the time-evolved density matrix $\left|\rho(t)\right\rangle$
in the form of Eq. (57) we can easily compute single or multi-boson
correlation functions. First, we need to make contact with our model
Liouvillian of Eq. (40). We have:

$\vec{\alpha}=\begin{pmatrix}\vec{\underline{d}}_{0}&\vec{\underline{d}}_{1}\end{pmatrix},\quad\hat{A}=-2\hat{X}^{T},\quad\hat{B}=\hat{Y}.$
(64)

Here the ’flavor’ index $\mu=0,1$ of $\vec{\underline{d}}_{\mu}$ is an effect
of third quantization whereas all other indices of $\vec{\underline{d}}$
describe physical bosonic modes of the system. Therefore it makes sense to
define a subblock structure of matrices $\hat{A}$, $\hat{B}$ and $\hat{C}(t)$,
i.e.:

$\hat{A}=\begin{pmatrix}\hat{A}^{00}&\hat{A}^{01}\\\
\hat{A}^{01}&\hat{A}^{11}\end{pmatrix}\\!,\,\hat{B}=\begin{pmatrix}\hat{B}^{00}&\hat{B}^{01}\\\
\hat{B}^{01}&\hat{B}^{11}\end{pmatrix}\\!,\,\hat{C}=\begin{pmatrix}\hat{C}^{00}&\hat{C}^{01}\\\
\hat{C}^{01}&\hat{C}^{11}\end{pmatrix}\\!.$ (65)

Now we can calculate some correlation function. Start with creator-annihilator
average. Following third quantization we get:

$\displaystyle{\rm Tr}\left(d_{i}^{\dagger}d_{j}\rho(t)\right)=\left\langle
0\right|\left(d^{\prime}_{i,0}+d_{i,1}\right)d_{j,0}\left|\rho(t)\right\rangle=\left\langle
0\right|d_{i,1}d_{j,0}\left|\rho(t)\right\rangle.$ (66)

The last expression can be easily calculate from the integral form of the
density matrix $\left|\rho(t)\right\rangle$ of Eq. (55), namely:

$\displaystyle{\rm
Tr}\left(d_{i}^{\dagger}d_{j}\rho(t)\right)=\frac{1}{\pi^{2L}}\int
d\Im\vec{z}d\Re\vec{z}\,\,\,z_{i,1}z_{j,0}\times\exp\left[-\vec{z}^{\,\star}\vec{z}+\vec{z}^{\,\star}\exp\left[t\hat{A}\right]\vec{z}(0)+\vec{z}^{\,\star}\hat{C}(t)\vec{z}^{\,\star}\right].$
(67)

From the Gaussian integral we get:

$\displaystyle{\rm
Tr}\left(d_{i}^{\dagger}d_{j}\rho(t)\right)=C(t)^{10}_{ij}+C(t)^{01}_{ji}+\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{i,1}\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{j,0}.$
(68)

In analogical way we can get other single-boson correlators, i.e.,

$\displaystyle{\rm
Tr}\left(d_{i}^{\dagger}d_{j}^{\dagger}\rho(t)\right)=C(t)^{11}_{ij}+C(t)^{11}_{ji}+\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{i,1}\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{j,1}$
(69)

and

$\displaystyle{\rm
Tr}\left(d_{i}d_{j}\rho(t)\right)=C(t)^{00}_{ij}+C(t)^{00}_{ji}+\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{i,0}\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{j,0}.$
(70)

If the starting state is not vacuum but a generic coherent state then also the
linear averages of creation/annihilation operators are non-trivial. We get:

$\displaystyle{\rm
Tr}\left(d_{i}^{\dagger}\rho(t)\right)=\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{i,1},$
(71) $\displaystyle{\rm
Tr}\left(d_{i}\rho(t)\right)=\left(\exp\left[t\hat{A}\right]\vec{z}(0)\right)_{i,0}.$
(72)

###  VIII.6 Non-zero temperature

In presence of non-zero temperature the Liouvillan of Eq. (40) gets modified
by new form of the $\hat{X}$ and $\hat{Y}$ matrices, one gets:

$\hat{X}=\frac{1}{2}\begin{pmatrix}i\hat{H}^{\star}+\hat{M}-\hat{N}^{\star}&-2i\hat{K}\\\
2i\hat{K}^{\star}&-i\hat{H}+\hat{M}^{\star}-\hat{N}\end{pmatrix},\,\hat{Y}=\begin{pmatrix}-i\hat{K}^{\star}&\hat{N}\\\
\hat{N}^{\star}&i\hat{K}\end{pmatrix},$ (73)

where the definitions of the $\hat{M}$ and $\hat{N}$ block are given by:

$\displaystyle\hat{N}$ $\displaystyle=$
$\displaystyle\hat{T}_{0}\otimes\frac{1}{2}{\rm diag}\begin{pmatrix}\kappa
n_{a}&\gamma n_{b}&\kappa n_{a}&\gamma n_{b}&\kappa n_{a}&\gamma n_{b}&\kappa
n_{a}&\gamma n_{b}\end{pmatrix},$ (74) $\displaystyle\hat{M}$ $\displaystyle=$
$\displaystyle\hat{T}_{0}\otimes\frac{1}{2}{\rm
diag}\begin{pmatrix}\kappa(n_{a}+1)&\gamma(n_{b}+1)&\kappa(n_{a}+1)&\gamma(n_{b}+1)&\kappa(n_{a}+1)&\gamma(n_{b}+1)&\kappa(n_{a}+1)&\gamma(n_{b}+1)\end{pmatrix}.$
(75)

Here $n_{a}$ and $n_{b}$ are thermal occupations of the $a$ and $b$ bosonic
modes, so for zero temperature $\hat{N}=0$ and the $\hat{X}$ and $\hat{Y}$
matrices take the old form of Eq. (42). With new form of the Liouvillan the
rest of derivation is done in the same manner.

Figure 7:  Stationary-state densities for $\gamma=0.0001$ and $\kappa=1.0$.
Dashed line indicates boundary of end-states hidden in the bulk. Densities
diverge in the white excluded region as there is no stationary state.

###  VIII.7 Stationary correlation matrix

Using equations (68-70) it is simple to relate matrix $\hat{C}$ with the
correlation matrix $\hat{V}$. First one has to notice that from the solution
of the Lyapunov equation it follows that $\hat{C}$ is symmetric. Therefore
taking initial condition as vacuum, i.e., $\vec{z}(0)=0$ we get for stationary
state averages:

$\displaystyle C_{ij}^{00}$ $\displaystyle=$ $\displaystyle\langle
d_{i}d_{j}\rangle/2$ $\displaystyle C_{ij}^{11}$ $\displaystyle=$
$\displaystyle\langle
d_{i}^{\dagger}d_{j}^{\dagger}\rangle/2=C_{ij}^{00\star}$ $\displaystyle
C_{ij}^{10}$ $\displaystyle=$ $\displaystyle\langle
d_{i}^{\dagger}d_{j}\rangle/2$ $\displaystyle C_{ij}^{01}$ $\displaystyle=$
$\displaystyle\langle d_{i}d_{j}^{\dagger}-\delta_{ij}\rangle/2=\left\langle
a_{j}^{\dagger}a_{i}\right\rangle/2=C_{ji}^{10}$ (76)

The correlation matrix can be written blockwise as:

$\hat{V}=\begin{pmatrix}\hat{V}^{00}&\hat{V}^{01}\\\
\hat{V}^{10}&\hat{V}^{11}\end{pmatrix}$ (77)

with blocks given by elements: $V^{00}_{ij}=\langle
d_{i}d_{j}^{\dagger}\rangle$, $V^{01}_{ij}=\langle d_{i}d_{j}\rangle$,
$V^{10}_{ij}=\langle d_{i}^{\dagger}d_{j}^{\dagger}\rangle$ and
$V^{11}_{ij}=\langle d_{i}^{\dagger}d_{j}\rangle$. Thus $\hat{C}$ matrix can
be written as:

$\hat{C}=\frac{1}{2}\begin{pmatrix}\hat{V}^{01}&\hat{V}^{00}-1\\\
\hat{V}^{11}&\hat{V}^{10}\end{pmatrix}=\frac{1}{2}\left[\begin{pmatrix}\hat{V}^{00}&\hat{V}^{01}\\\
\hat{V}^{10}&\hat{V}^{11}\end{pmatrix}-\begin{pmatrix}1&0\\\
0&0\end{pmatrix}\right]\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}.$ (78)

Now we can insert such form of matrix $\hat{C}$ into the stationary-state
Lyapunov equation:

$\hat{A}\hat{C}+\hat{C}\hat{A}^{T}=-\hat{B}$ (79)

we get an equation for the $\hat{V}$ matrix:

$\displaystyle\begin{pmatrix}2\hat{M}&0\\\ 0&2\hat{N}\end{pmatrix}=$
$\displaystyle\left[i\begin{pmatrix}\hat{H}-i(\hat{M}-\hat{N})&2\hat{K}\\\
-2\hat{K}&-\hat{H}-i(\hat{M}-\hat{N})\end{pmatrix}\right]\begin{pmatrix}\hat{V}^{00}&\hat{V}^{01}\\\
\hat{V}^{10}&\hat{V}^{11}\end{pmatrix}$ $\displaystyle+$
$\displaystyle\begin{pmatrix}\hat{V}^{00}&\hat{V}^{01}\\\
\hat{V}^{10}&\hat{V}^{11}\end{pmatrix}\left[i\begin{pmatrix}\hat{H}-i(\hat{M}-\hat{N})&2\hat{K}\\\
-2\hat{K}&-\hat{H}-i(\hat{M}-\hat{N})\end{pmatrix}\right]^{\dagger},$ (80)

or simply:

$(i\hat{H}_{\rm NH})\hat{V}+\hat{V}(i\hat{H}_{\rm
NH})^{\dagger}=\begin{pmatrix}2\hat{M}&0\\\ 0&2\hat{N}\end{pmatrix}.$ (81)

For this derivation we have taken advantage of the present current form of the
Hamiltonian with $\hat{H}=\hat{H}^{T}$, $\hat{K}=\hat{K}^{T}=\hat{K}^{\star}$.
Note that the non-Hermitian Hamiltonian $\hat{H}_{\rm NH}$ does not depend on
temperature.

Figure 8:  Entanglement patterns for $\gamma=0.0001$ and $\kappa=1.0$.
Dark/bright color indicates non-vanishing/vanishing stationary state- or
saturation-negativity between two modes denoted in cyan in bottom right corner
of each panel. Red and green lines indicates boundaries of: end-states hidden
in the bulk and end of the system having stationary state. In the excluded
white region the bulk has no stationary state. In the region between green
line and white area the populations are diverging due to the instability of
the end-states but negativity remains finite.

###  VIII.8 Entanglement in the presence of unstable end modes

Surprisingly, the entanglement can be calculated also when there is no
stationary state because the end modes are unstable. In this case there are
four degenerate end-modes which become unstable so the $\hat{C}(t)$ matrix can
be evaluated from Eq. (62) in the limit of $t\to\infty$ by setting:

$\Gamma_{ij}(t\to\infty)\to\frac{1}{\lambda_{i}+\lambda_{j}}\left[\begin{pmatrix}-1&\cdots&-1\\\
\vdots&\ddots&\vdots\\\
-1&\cdots&-1\end{pmatrix}\oplus\begin{pmatrix}\xi&\xi&\xi&\xi\\\
\xi&\xi&\xi&\xi\\\ \xi&\xi&\xi&\xi\\\
\xi&\xi&\xi&\xi\end{pmatrix}\right]_{ij}$ (82)

with $\xi=\exp\left[2t\lambda_{{\rm end}}\right]$. Clearly $\xi$ diverges when
$t\to\infty$. Using this form of $\hat{\Gamma}$ matrix and keeping $\xi$ as
analytical variable we can express $\hat{C}(\xi\to\infty)$ and any two-point
correlation function as analytical expression in $\xi$. Then we can quantify
the entanglement for $t\to\infty$ by computing $\tilde{\nu}_{-}$ and taking
analytical limit of $\tilde{\nu}_{-}(\xi\to\infty)$. Quite non-trivially this
limit is well-defined and finite despite diverging populations of bosons.

##  IX Entanglement measurement

We sketch here a possible detection scheme for the optomechanical entanglement
within the end sites of the chain. Following a known strategy [15], we suggest
that the inclusion of an extra cavity coupled to the mechanical resonators at
each end of the chain can be employed to monitor the mechanical motion of the
extremal sites of the superlattice. To this end, each detection cavity is
weakly driven on the red sideband, inducing an effective beam-splitter
interaction between mechanical ($b_{\rm 1}$, $b_{\rm 4N}$) and readout ($c$,
$d$) degrees of freedom (see Fig. 9).

Figure 9: Cartoon depiction of the left end in the presence of an extra
detection cavity weakly driven on the red sideband (intracavity mode $c$,
input/output modes $c_{\rm in}$, $c_{\rm out}$).

In the following we focus, without loss of generality on the first site of the
chain. The detection cavity described by bosonic operators $c$ and
$c^{\dagger}$, is weakly driven on the red sideband, giving rise (in the
appropriate frame) to the following EOM for $c$ in Fourier space

$-i\omega c=\frac{\kappa}{2}c-iG_{\rm R}b_{\rm 1}+\sqrt{\kappa}\,c$ (83)

where, $G_{\rm R}$ and $\kappa_{\rm R}$ are the coupling induced by the
detection drive and the detection cavity linewidth, respectively. Hereafter,
we will assume $\kappa_{\rm R}=\kappa$. Taking into account the usual input-
output relations

$\displaystyle c_{\rm out}$ $\displaystyle=\sqrt{\kappa}\,c-c_{\rm in}$ (84a)
$\displaystyle a_{\rm out}$ $\displaystyle=\sqrt{\kappa}\,a-a_{\rm in}$ (84b)

[16], the output field for the readout cavity can be written as

$\displaystyle c_{\rm out}$ $\displaystyle=\left(\chi_{\rm
c}\kappa-1\right)c_{\rm in}-i\chi_{\rm c}G_{\rm R}\sqrt{\kappa}\,b_{\rm 1}.$
(85)

If we assume that the detection tone $G_{\rm R}$ is much weaker than the
driving tone $G_{+}$, the dynamics of the optomechanical system is essentially
unchanged by the detection process, simply leading to the renormalization of
$\gamma$ and $n_{\rm eff}$. As far as the superlattice is concerned, we have
demonstrated that our protocol is robust with respect to changes in the
thermal population and, in Section VII.3 of the SI, we have numerically
demonstrated the fact that, owing to the topological features of our system,
enlargement is robust also for large variations of the mechanical linewidth on
the end sites.

Assuming that the cavities are coupled to a zero-temperature bath, we can
relate the (normal-ordered) power spectrum of the correlations between the
output fields

$\displaystyle S^{\rm out}_{\rm
1,2}\left(\omega\right)=\int\frac{d\omega^{\prime}}{2\pi}\braket{O^{\dagger}_{\rm
1,out}(\omega)O_{\rm 2,out}(\omega^{\prime})}$ (86)

(where $O_{\rm 1,out},O_{\rm 2,out}=a_{\rm out},c_{\rm out}$) to the power
spectrum of the correlations between $a$ and $c$ (and, consequently, $b$) [9,
16], given that

$\displaystyle S^{\rm out}_{\rm 1,2}\left(\omega\right)=\kappa S_{\rm
1,2}\left(\omega\right).$ (87)

## References

  * Safavi-Naeini and Painter [2010] A. H. Safavi-Naeini and O. Painter, Design of optomechanical cavities and waveguides on a simultaneous bandgap phononic-photonic crystal slab, Optics Express 18, 14926 (2010).
  * Painter  _et al._ [1999] O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, Two-Dimensional Photonic Band-Gap Defect Mode Laser, Science 284, 1819 (1999).
  * Vučković  _et al._ [2002] J. Vučković, M. Lončar, H. Mabuchi, and A. Scherer, Design of photonic crystal microcavities for cavity QED, Physical Review E 65, 016608 (2002).
  * Vasseur  _et al._ [2007] J. O. Vasseur, A.-C. Hladky-Hennion, B. Djafari-Rouhani, F. Duval, B. Dubus, Y. Pennec, and P. A. Deymier, Waveguiding in two-dimensional piezoelectric phononic crystal plates, Journal of Applied Physics 101, 114904 (2007).
  * III and El-Kady [2009] R. H. O. III and I. El-Kady, Microfabricated phononic crystal devices and applications, Measurement Science and Technology 20, 012002 (2009).
  * Ludwig and Marquardt [2013] M. Ludwig and F. Marquardt, Quantum Many-Body Dynamics in Optomechanical Arrays, Physical Review Letters 111, 073603 (2013).
  * Brzezicki  _et al._ [2023] W. Brzezicki, M. Silveri, M. Płodzień, F. Massel, and T. Hyart, Non-Hermitian topological quantum states in a reservoir-engineered transmon chain, Physical Review B 107, 115146 (2023).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Gardiner and Zoller [2004] C. W. Gardiner and P. Zoller,  _Quantum noise_ , Springer (Springer, 2004).
  * Gardiner [2004] C. W. Gardiner,  _Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences_ , Springer Series in Synergetics (Springer series in synergetics, 2004).
  * Adesso and Illuminati [2007] G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, Journal of Physics A: Mathematical and Theoretical 40, 7821 (2007).
  * Brzezicki and Hyart [2019] W. Brzezicki and T. Hyart, Hidden Chern number in one-dimensional non-Hermitian chiral-symmetric systems, Physical Review B 100, 161105 (2019).
  * Youssefi  _et al._ [2022] A. Youssefi, S. Kono, A. Bancora, M. Chegnizadeh, J. Pan, T. Vovk, and T. J. Kippenberg, Topological lattices realized in superconducting circuit optomechanics, Nature 612, 666 (2022).
  * Prosen and Seligman [2010] T. Prosen and T. H. Seligman, Quantization over boson operator spaces, Journal of Physics A: Mathematical and Theoretical 43, 392004 (2010).
  * Vitali  _et al._ [2007] D. Vitali, S. Gigan, A. Ferreira, H. R. Böhm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Optomechanical Entanglement between a Movable Mirror and a Cavity Field, Physical Review Letters 98, 030405 (2007).
  * Walls and Milburn [2008] D. F. Walls and G. J. Milburn,  _Quantum Optics_, Quantum Optics (Springer Berlin Heidelberg, 2008).