# Comment on ”Scalings of Rayleigh-Taylor Instability at Large Viscosity
Contrasts
in Porous Media”

E.B. Soboleva  E-mail address<EMAIL_ADDRESS>Ishlinsky Institute for
Problems in Mechanics RAS,
Prospekt Vernadskogo 101-1, Moscow, 119526, Russia

††preprint: APS/123-QED

In a recent Letter, Sabet et al. [1] presented new numerical results on a
development of Rayleigh-Taylor (RT) instability in porous media. Two miscible
fluids are assumed to be of different viscosity. The main finding is that a
symmetry of gravity-driven convection is broken at large viscosity ratios
$M=\mu_{2}^{*}/\mu_{1}^{*}>20$ and convective fingers extend preferably in the
fluid of lower viscosity. If the upper fluid is much more viscous, as
considered in the Letter, upward fingers extend at a much smaller rate than
downward fingers. The authors try to explain the effect revealed by analyzing
the vorticity field. The purpose of this Comment is to show that (1) analyzing
the vorticity field does not allow to explain this effect, (2) it can be
explained by analyzing the velocity in the equation of motion.

(1) Sabet et al. [1] decompose the vorticity field into density and viscosity
counter-parts $\omega^{(\rho)}$ and $\omega^{(\mu)}$, respectively, and
analyze their contributions. The authors treat the vorticity field as the
driving force of convection that is not correct. They wrote, ”The density
component of vorticity $\omega^{(\rho)}$, which is the main driving force
behind the growth of the fingers, is inversely proportional to the local value
of viscosity. In other words, as the viscosity of the mixing layer increases,
the driving force for buoyancy-driven convection is weakened…” Actually, the
fingers grow driven by the gravitational and buoyant forces.

If the vorticity were able to cause the fingering asymmetry (according to the
statement of Authors), it would influence the vertical translation of fluid.
In this case, the vorticity has to change the vertical velocity of finger
tips, which are the leading points of motion, and, particularly, to generate
an additional negative vertical velocity of tips of upward fingers mitigating
their motion. However, as clearly seen in Figure 1(b) in [1], red and blue
colors disappear at tips of both upward and downward fingers indicating that
the value of vorticity is zero here. We see that the vorticity is absent at
tips of both upward and downward fingers and therefore does not affect the
vertical fluid velocity. One could assume that the additional vertical
velocity of tip motion is generated by the vorticity near this tip. However,
we do not find any explanation, why and in what way local rotations of fluid
near the tip are transformed into an additional vertical translation of this
tip.

FIG. 1. Sketch of the problem.

(2) In [1], fluid motions are described by the dimensionless Darcy equation
[see Eq. (2)]. We transform this equation replacing the Darcy velocity
$\mathbf{v}$ with the fluid velocity $\mathbf{V}=\phi^{-1}\mathbf{v}$
[$\mathbf{V}=(U,W)$] and express the z-component of the velocity, $W$. One can
denote $\partial p/\partial z+\rho(c)$ as $T$ and obtain the following:

$W=-\frac{1}{\phi\mu(c)}T,\ T=\frac{\partial p}{\partial z}+\rho(c)\ .$ (1)

The Cartesian coordinate system with the $z$ axis being vertical and oriented
upward and all notations as in [1] are used. On the tip of downward finger,
for example, at point $A$ in Fig. 1, fluid moves down at $W<0$ that leads to
$T>0$ here, according to Eq. (1); the values of $\phi$ and $\mu(c)$ are
positive. On the tip of upward finger at point $B$, $W>0$, and $T<0$.

We find the partial derivative of $W$ with respect to $z$:

$\frac{\partial W}{\partial z}=-\frac{T}{\phi}\frac{\partial}{\partial
z}\left(\frac{1}{\mu(c)}\right)-\frac{1}{\phi\mu(c)}\frac{\partial T}{\partial
z}\ .$ (2)

Further, the case of extremely large $M$ will be considered. In the mixed
zone, the concentration $c$ changes from $0$ to $1$ leading to a very large
increase in the viscosity $\mu$ from $1$ to $M$. The viscosity changes per the
unit of height are especially large on tips of convective fingers, therefore
one can assume that the velocity variations are determined here mostly by the
viscosity inhomogeneities. So, the second term on the right-hand side of Eq.
(2) contributes less than the first term. It means that

$\arrowvert\frac{T}{\phi}\frac{\partial}{\partial
z}\left(\frac{1}{\mu(c)}\right)\arrowvert>\arrowvert\frac{1}{\phi\mu(c)}\frac{\partial
T}{\partial z}\arrowvert\ .$ (3)

Under the last assumption, the sign of the left-hand side of Eq. (2) coincides
with the sign of the first term on the right-hand side, therefore, to analyze
the sign of the derivative $\partial W/\partial z$, the second term can be
neglected. We also transform the derivative of $1/\mu(c)$ and obtain

$\frac{\partial W}{\partial z}\approx
T\frac{1}{\phi\mu^{2}(c)}\frac{\partial\mu(c)}{\partial c}\frac{\partial
c}{\partial z}\ .$ (4)

In Eq. (4), $\partial\mu(c)/\partial c>0$ because $\mu$ is a growing function
of $c$, and $\partial c/\partial z>0$ because $c$ increases with height.
Consequently, the signs of $\partial W/\partial z$ and $T$ coincide. At point
$A$, $T>0$ leading to $\partial W/\partial z>0$. Since $W<0$ here, $W$ becomes
smaller or the absolute value $\mid W\mid$ becomes larger with decreasing $z$.
In other words, if fluid moves down, its absolute velocity $\mid W\mid$
becomes larger. At point $B$, $T<0$ and $\partial W/\partial z<0$, that is $W$
becomes smaller if fluid moves up.

In the case of decreasing $\mu$ with $c$, $\partial\mu(c)/\partial c<0$ and
the signs of $\partial W/\partial z$ and $T$ are opposite. Further reasoning
shows that, on the contrary, upward fingers have to accelerate and downward
fingers have to decelerate. In summary, due to viscosity inhomogeneities,
convective fingers extend preferably in the fluid of lower viscosity. This
conclusion is consistent with one in [1] and confirmed by our own simulation
[2].

The study was supported by the Government program (contract
AAAA-A20-120011690131-7).

## References

  * [1] N. Sabet, H. Hassanzadeh, A. De Wit, and J. Abedi. Scalings of Rayleigh-Taylor Instability at Large Viscosity Contrasts in Porous Media. Phys. Rev. Lett. 126, 094501 (2021). https://doi.org/10.1103/PhysRevLett.126.094501
  * [2] E. B. Soboleva (unpublished).


# A two-vertex theorem for normal tilings

Gábor Domokos, Ákos G. Horváth and Krisztina Regős  Gábor Domokos, MTA-BME
Morphodynamics Research Group and Dept. of Mechanics, Materials and
Structures, Budapest University of Technology and Economics, Műegyetem rakpart
1-3., Budapest, Hungary, 1111<EMAIL_ADDRESS>Ákos G. Horváth, Dept. of
Geometry, Budapest University of Technology and Economics, Egry József utca
1., Budapest, Hungary, 1111<EMAIL_ADDRESS>Krisztina Regős,MTA-BME
Morphodynamics Research Group, Budapest University of Technology and
Economics, Műegyetem rakpart 1-3., Budapest, Hungary, 1111
<EMAIL_ADDRESS>
###### Abstract.

We regard a smooth, $d=2$-dimensional manifold $\mathcal{M}$ and its normal
tiling $M$, the cells of which may have non-smooth or smooth vertices (at the
latter, two edges meet at 180 degrees.) We denote the average number (per
cell) of non-smooth vertices by $\bar{v}^{\star}$ and we prove that if $M$ is
periodic then $\bar{v}^{\star}\geq 2$. We show the same result for the
monohedral case by an entirely different argument. Our theory also makes a
closely related prediction for non-periodic tilings. In 3 dimensions we show a
monohedral construction with $\bar{v}^{\star}=0$.

###### Key words and phrases:

normal tiling, sharp vertex

###### 2010 Mathematics Subject Classification:

52C20

Support of the NKFIH Hungarian Research Fund grant 134199 and of the NKFIH
Fund TKP2021 BME-NVA, carried out at the Budapest University of Technology and
Economics, is kindly acknowledged.

##  1\. Introduction

Both in nature and among man-made structures we often encounter cellular
patterns where a domain is decomposed into mutually disjoint, non-overlapping
cells. Such patterns, also called _tessellations_ or _tilings,_ arise in
geology (as fracture patterns [1, 3, 5, 4, 14]), in engineering and
architecture (as surface patterns or as urban structures [10, 2]), in
geography [6] or in biology (as cell tissues [8, 12]).

In general, we will assume that the boundaries of cells are piecewise
$C^{1}$-smooth, however, in subsection 3.1 we will require piecewise
$C^{2}$-smoothness. We also assume that the tiling is _normal_[7], i.e. that
every cell is topologically equivalent to a disk, the intersection of any two
cells is a connected set or the empty set, all cells are uniformly bounded. In
addition, for inifinite tilings we will assume that they are also _balanced_
[7], i.e. that the average number of vertices and edges exist independently
from the method of the counting, the detailed description of the latter
process is given in Remark 2. Normal tilings have been described in various
settings [9, 11, 13], however, the emphasis was on their combinatorial
properties. Here we aim to describe a curios _metric_ aspect of these
patterns. Cells may have non-smooth or smooth vertices (at the latter, two
edges meet at 180 degrees). We denote the average number (per cell) of the
former by $\bar{v}^{\star}$ and we prove that for periodic patterns
$\bar{v}^{\star}\geq 2$ (see Theorem 1 and Remark 1). In Theorem 1 we also
provide a general formula for the lower bound on $\bar{v}^{\star}$, including
non-periodic patterns on arbitrary, smooth, 2D manifolds. In addition, we
prove the $\bar{v}^{\star}\geq 2$ statement for monohedral tilings (consisting
of congruent cells) based on a different argument. Finally, we offer a
3-dimensional monohedral construction with $v^{\star}=0$.

###  1.1. A simple example and the main result

We start by giving an informal example, illustrating the basic concepts, based
on which we can formulate our main result and we will provide more formal
definitions later. We regard normal tilings of smooth, 2D manifolds. Figure 1
illustrates three such patterns (illustrated in columns (A),(B),(C),
respectively) , all of which are periodic. In the bottom row we indicate some
real-world examples which could be approximated by these patterns. We will
characterize such patterns by the _combinatorial nodal degree $\bar{n}$ and
combinatorial cell degree $\bar{v}$_ which count the average number of cells
meeting at one node and the average number of nodes along the boundary of one
cell. We could also interpret these degrees as counting the number of vertices
meeting at a node and counting the number of vertices along the perimeter of a
cell. In case of periodic patterns, it is easy to compute these averages; for
the particular patterns illustrated in Fig 1 all nodes and cells are identical
and we can see that the combinatorial degrees, illustrated and counted in the
upper row, are identical for all three patterns.

Figure 1. Combinatorial and corner degrees in periodic patterns. (A) Hexagonal
pattern reminiscent of honeycomb with
$(\bar{n},\bar{v})=(\bar{n}^{\star},\bar{v}^{\star})=(3,6)$. (B) Pattern
reminiscent of plant cell tissue or brick wall, with combinatorial degrees
$(\bar{n},\bar{v})=(3,6)$, corner degrees
$(\bar{n}^{\star},\bar{v}^{\star})=(2,4)$. (C) Pattern reminiscent of roof
tiles with combinatorial degrees $(\bar{n},\bar{v})=(3,6)$, corner degrees
$(\bar{n}^{\star},\bar{v}^{\star})=(1,2)$.

Despite having identical combinatorial degrees, the three patterns look
radically different. One way to capture these differences is to compute
_corner degrees_ : the _nodal corner degree_ $\bar{n}^{\star}$ and the _cell
corner degree_ $\bar{v}^{\star}$ count only vertices where the angle of the
half-tangents of the meeting arcs is not 180 degrees (and so, obviously, we
have $\bar{n}^{\star}\leq\bar{n}$, $\bar{v}^{\star}\leq\bar{v}$). In case of
the hexagonal honeycomb in column (A) of Fig. 1, the combinatorial and corner
degrees coincide since none of the examined angles is 180 degrees. In case of
the brick wall pattern the corner degrees are smaller because we do have
vertices with edge-angles 180 degrees. The roof tile pattern has even smaller
corner degrees.

While it is clear that corner degrees may be smaller than combinatorial
degrees, it is not obvious by how much these quantities may differ. In
particular, it is not clear what is the minimum for the cell corner degree
which counts the average number _actual, non-smooth_ corners of a cell. Our
goal is to prove

######  Theorem 1.

Let $\mathcal{M}$ be a smooth, 2D manifold with Euler characteristic $\chi$
and let $M$ be a normal tiling of $\mathcal{M}$ with $E$ edges and corner
degrees $\bar{n}^{\star},\bar{v}^{\star}$. Then we have

$\bar{v}^{\star}\geq 2-\frac{2\chi\bar{n}}{E(\bar{n}-2)+\chi\bar{n}}.$

As we can see, both in the $E\to\infty$ limit and for $\chi=0$ Theorem 1
yields

(1) $\bar{v}^{\star}\geq 2.$

######  Remark 1.

As periodic patterns can be treated as infinite, equation (1) holds for such
patterns.

###  1.2. Structure of paper

In Section 2 we first define the basic concepts and prove Theorem 1. In
Section 3 we offer an alternative proof for the 2D monohedral case as well as
a construction of a 3D monohedral tiling with no vertices.

##  2\. Normal tiling of 2D manifolds

###  2.1. Combinatorial properties

######  Definition 1.

Let $M$ be a normal [7] tiling of a smooth, $d=2$-dimensional manifold
$\mathcal{M}$ and let $M$ have $V$ vertices (0-dimensional cell or nodes), $E$
edges (1-dimensional cells) and $F$ faces (2-dimensional cells or simply
cells). We call the number of edges meeting at the $i$th node ($i=1,2,\dots
V$) the _combinatorial nodal degree_ and denote it by $n_{i}$. We call the
number of edges on the boundary of the $j$th cell ($j=1,2,\dots F$) the
_combinatorial cell degree_ and denote it by $v_{j}$. For finite $V,E,F$ we
denote the average values of the combinatorial nodal and cell degrees by
$\bar{n},\bar{v}$, respectively. We also allow that $V,E,F$ will be infinite,
however, in that case we also assume that $M$ is _balanced_ [7], guaranteeing
that a suitably chosen limit process on a sequence of finite tilings provides
the corresponding averages $\bar{n},\bar{v}$ for the nodal and cell degrees
(for details see Remark 2). We call the $[\bar{n},\bar{v}]$ plane as the
(combinatorial) _symbolic plane_.

######  Remark 2.

On compact manifolds and for finite tilings, counting averages is trivial. On
infinite domains averages are again trivial if $M$ is periodic. For general
normal tilings we compute the averages by counting the nodal and cell degrees
in a finite ball $B(X,\rho)$ with center $X$ radius $\rho$ to obtain
$\bar{n}(X,\rho),\bar{v}(X,\rho)$ associated with the finite normal tiling
$M(X,\rho)$. Then we let the radius $\rho$ approach infinity. We say that the
averages $\bar{n},\bar{v}$ of the degrees exists if the limits
$\lim_{\rho\to\infty}\bar{n}(X,\rho),\lim_{\rho\to\infty}\bar{v}(X,\rho)$
exist and are independent of $X$. Henceforth we will assume that the manifold
$\mathcal{M}$ is either compact and $M$ is finite, or, $M$ is infinite and the
averages converge in the sense described above and boundary effects decay.

######  Remark 3.

Note that on the Euclidean plane the concepts of a tiling being ”normal” and
”balanced” are well-investigated (see Section 3.3 in [7]). Clearly, the
definition of ”normality” formally can be applied on arbitrary manifolds, but
its precise connection with the property ”balanced” is far from obvious,
except on the Euclidean plane (see Theorem 3.3.2 of [7]. In this paper we rely
on the strong conditions on averaging, given in Remark 2).

######  Remark 4.

From the fact that the tiling $M$ is normal, it follows that each 2-cell is
_contractible_ , and the boundary of a 2-cell is a simple closed curve. We
will later discuss $d$-dimensional manifolds and their tilings where we also
use that the $d$-dimensional cells are contractible and all $k<d$ dimensional
faces are also contractible.

######  Remark 5.

Since all edges are assumed to be $C^{1}$-smooth, they have half-tangents at
both endpoints. We also note that the angle of the two half-tangent s can be
calculated from the globally defined differential structure of the smooth
manifold.

######  Definition 2.

We call

$\bar{h}=\frac{\bar{v}\bar{n}}{\bar{v}+\bar{n}}$

the _harmonic degree_ of $M$.

For finite $M$, the Euler formula can be written as

(2) $F-E+V=\chi.$

According to Assumption (1) in Definition 1, and the fact that every edge has
two vertices and belongs to two cells, the combinatorial nodal and cell
degrees can be expressed as

(3) $\bar{n}=\frac{2E}{V},\quad\bar{v}=\frac{2E}{F}.$

We can rewrite equation (2) as

(4) $\frac{F}{2E}+\frac{V}{2E}=\frac{E+\chi}{2E},$

which, using equation (3) and Definition 2 translates into

(5) $\bar{h}=\frac{2E}{E+\chi}.$

In the infinite case we regard an averaging process analogous to the one given
in Remark 2 where the number $E(\rho)$ of edges (and the harmonic degree
$\bar{h}(\rho)$) are defined as the function of the radius $\rho$ of a finite
ball. As $\rho\to\infty$, we have

(6) $\bar{h}:=\lim_{E(\rho)\to\infty}\bar{h}(\rho)=2,$

regardless of the Euler characteristic of $M=\lim_{\rho\to\infty}M(\rho)$.
Also, if $\chi=0$ then we have $\bar{h}=2$, one example is the torus $T^{2}$.

######  Remark 6.

Observe, that the finiteness of the normal tiling $M$ is equivalent to the
compactness of the carrying manifold $\mathcal{M}$.

Figure 2. Illustration of infinite regular patterns on the symbolic plane
$[\bar{n},\bar{v}]$. All patterns lie on the $\bar{h}=2$ hyperbola.

###  2.2. Geometric properties: sharp corners

######  Definition 3.

Let $M$ be a normal and balanced tiling of a smooth, $d=2$-dimensional
manifold $\mathcal{M}$ and let $M$ have $V$ vertices (nodes), $E$ edges and
$F$ faces (cells) and let $n_{i}$ denote the combinatorial degree of the $i$th
node, $(i=1,2,\dots V)$. We number the edges $e_{i,j}$ ($i=1,2\dots
V,j=1,2,\dots n_{i}$) meeting at the $i$th node clockwise consecutively using
their second subscript, so the _pair_ $e_{i,j},e_{i,j+1}$ is not separated by
any other edge. We also regard $e_{i,n_{i}},e_{i,1}$ as one pair, so we
associate $n_{i}$ pairs with the $i$th node. We call a pair _degenerate_ if
the two edges meet at 180 degrees, otherwise we call the pair _generic_. (The
angle between two the edges of an edge-pair is identical to the alternate
interior angle their respective half-tangents at the common vertex, thus we
can distinguish between the angle zero (when the two half-tangents have
opposite direction on their line) and the angle $\pi$ ( when the two half-
tangents have the same direction).) We associate with the $i$th node the
quantity $r_{i}\in\\{0,1,2\\}$, counting the number of degenerate pairs. We
define the _corner degree_ of the $i$th node as

$n^{\star}_{i}=n_{i}-r_{i}$

.

######  Definition 4.

Let $M$ be a normal and balanced tiling of a smooth, $d=2$-dimensional
manifold $\mathcal{M}$ and let $M$ have $V$ vertices (nodes), $E$ edges and
$F$ faces (cells) and let $v_{i}$ denote the combinatorial degree of the $i$th
cell, $(i=1,2,\dots F)$. The edges $e_{i,j}$ ($i=1,2\dots F,j=1,2,\dots
n_{i}$) surrounding the $i$th cell are numbered clockwise consecutively, so
the _pair_ $e_{i,j},e_{i,j+1}$ overlaps at a vertex. We also regard
$e_{i,v_{i}},e_{i,1}$ as one pair, so we associate $v_{i}$ pairs with the
$i$th cell. We call a pair _degenerate_ if the two edges meet at 180 degrees,
otherwise we call the pair _generic_. We associate with the $i$th cell the
quantity $q_{i}\in\\{0,1,2,\dots v_{i}\\}$, counting the number of degenerate
pairs. We define the _corner degree_ of the $i$th cell as

$v^{\star}_{i}=v_{i}-q_{i}$

.

######  Definition 5.

We call the $[\bar{n}^{\star},\bar{v}^{\star}]$ plane the (geometric)
_symbolic plane._

######  Definition 6.

Let

$\rho(M)=1-\frac{\sum_{i=1}^{V}r_{i}}{2V}$

and we call $M$ $\rho$-regular. In particular, we call $M$ regular, semi-
regular or irregular if $\rho(M)=1,1/2,0$, respectively.

###  2.3. Trivial lower bounds on the cell and nodal degrees for large cell
decompositions

Here we regard the $E\to\infty$ limit of a balanced normal tiling . If $M$ is
regular ($\rho=1$) then we have

(7) $\bar{v}=\bar{v}^{\star},\quad\bar{n}=\bar{n}^{\star}.$

Equation (6) implies

(8) $\bar{n},\bar{v}\geq 2,$

which, via (7) is equivalent to

(9) $\bar{n}^{\star},\bar{v}^{\star}\geq 2.$

If, in addition to being regular, $M$ is also combinatorially equivalent to a
convex mosaic then we have

(10) $\bar{n}^{\star},\bar{v}^{\star}\geq 3.$

All the listed bounds are trivial. Also, because of the inverse functional
relationship (6) between $\bar{n}$ and $\bar{v}$, the _lower_ bound of one of
the combinatorial averages is tied to the _upper_ bound of the other
combinatorial average.

If $M$ is not regular $(\rho<1)$ then there is no general functional
relationship between $\bar{n}$ and $\bar{v}$. Based on Definitions 3 and 4 for
non-regular tilings we have

(11) $\bar{v}>\bar{v}^{\star},\quad\bar{n}>\bar{n}^{\star}$

and we can also see that the corner degrees attain their respective minima if
$M$ is irregular, i.e. we have $\rho=0$. Our next goal is to find the lower
bound for $\bar{v}^{\star}$ for irregular normal tilings.

###  2.4. Proof of Theorem 1

###### Proof.

Let $M$ be a regular normal tiling with $F$ faces (cells), $V$ vertices
(nodes) and $E$ edges and respective cell and nodal combinatorial degrees
given in equation (3). Let us _deregularize_ $M$ by transforming generic edge-
pairs into degenerate ones. We perform this by local, continuous
transformation in the vicinity of the nodes. Assume that by such local
transformations we create $2V(1-\rho)$ degenerate pairs to obtain the
$\rho$-regular tiling $M(\rho)$. While the combinatorial degrees remain
constant, the corner degrees can be computed as a function of $\rho$:

(12)
$\bar{v}^{\star}(\rho)=\frac{2E-2(1-\rho)V}{F}=\bar{v}\left(\frac{\bar{n}-2(1-\rho)}{\bar{n}}\right),$
$\bar{n}^{\star}(\rho)=\frac{2E-2(1-\rho)V}{V}=\bar{n}-2(1-\rho).$

We can immediately see that

(13)
$\frac{\bar{v}}{\bar{n}}=\frac{\bar{v}^{\star}(\rho)}{\bar{n}^{\star}(\rho)},$

so the evolution path of topological deregularization on the
$[\bar{n},\bar{v}]$ symbolic plane is a straight line passing through the
origin. By subsitution $\rho=0$ into (12) we obtain

(14) $\bar{n}^{\star}\geq\frac{2E-2V}{V}=\bar{n}-2,$ (15)
$\bar{v}^{\star}\geq\frac{2E-2V}{F}=\bar{v}\left(\frac{\bar{n}-2}{\bar{n}}\right).$

Equation (5) is equivalent to

(16) $\bar{v}=\frac{2E\bar{n}}{\bar{n}(E+\chi)-2E}.$

Substituting (16) into (15) yields the formula of Theorem 1. ∎

Figure 3. Illustration of infinite patterns on the symbolic plane
$[\bar{n}^{\star},\bar{v}^{\star}]$. Observe that combinatorially equivalent
tilings are connected by rays passing through the origin, illustrating formula
(13).

##  3\. Related issues and discussion

###  3.1. The monohedral case in $d=2$ dimensions

Here we discuss the special case where the tiling of the Euclidean plane
consists of congruent cells $\mathcal{C}$. Our next argument requires that at
every smooth point of the boundary of $\mathcal{C}$ the signed curvature
$\kappa$ exists. To guarantee this property we assume here that the smooth
part of the boundary is finite union of two-times continuously differentiable
arcs and for brevity, we call such tilings $C^{2}$-tilings. Theorem 1 implies

######  Corollary 1.

In case of monohedral, normal $C^{2}$-tilings of the Euclidean plane, for the
corner degree of the cell we have $v^{\star}\geq 2$.

Below we give a proof for Corollary 1 which is independent of the proof of
Theorem 1.

###### Proof.

Let us parametrize the circumference of the cell $\mathcal{C}$ by the
arclength $s$ on the unit interval $I\equiv[0,1]$, so we have $s\in[0,1]$. The
boundary of the cell $\mathcal{C}$ is piecewise $C^{2}$-smooth. We will
consider $I$ as the union of the (open) set of smooth points $I_{S}$ and the
(discrete) set of combinatorial nodes $N_{i}$, $i=1,2,\dots v$:

(17) $I=I_{S}\bigcup\\{N_{i}\\}.$

We denote the angle of the tangent by $\alpha(s)$ and at the nodes we will
denote the finite non-negative angles by $\Delta\alpha_{i},(i=1,2,\dots v)$.
We consider every combinatorial node to be a non-smooth point of the boundary
of $C$. The integral $K$ of the (signed) scalar curvature $\kappa$ can be
written as

(18) $K=\int_{I_{S}}\kappa(s)ds+\sum_{i=1}^{v}\Delta\alpha_{i}=2\pi.$

We observe that the discrete angle increments $\Delta\alpha_{i}$ differ from
zero _only_ at corners. To simplify our computation, we re-assign indices
$i\in[1,2,\dots v^{\star}]$ to these non-zero relative angles and so, instead
of (18) we may write:

(19) $K=\int_{I_{S}}\kappa(s)ds+\sum_{i=1}^{v^{\star}}\Delta\alpha_{i}=2\pi.$

Since the tiling is monohedral, the neighbor cells are congruent to the
investigated cell $C$. Now we observe that for every smooth part $\gamma$ of
the boundary there is at least one orientation preserving congruence
$f_{\gamma}:\mathcal{M}\to\mathcal{M}$ which sends the cell $C_{\gamma}$,
overlapping with $C$ along the arc $\gamma$, to $C$. Let $\gamma^{\prime}$ be
the image of $\gamma$ by $f_{\gamma}$. If $\gamma^{\prime}=\gamma$ then the
value of the integral (19) on $\gamma$ is zero. If this is not the case then
the value of the integral on the union $\gamma\cup\gamma^{\prime}$ is zero.
Since $f_{\gamma}(C)$ is the neighbour of $C$ along the arc $\gamma^{\prime}$,
if there exists another smooth part $\gamma^{\star}$ for which
$f_{\gamma^{\star}}(\gamma^{\star})=\gamma^{\prime}$ then we have
$(f_{\gamma}^{-1}\circ f_{\gamma^{\star}})(C)=C$ and so $f_{\gamma}^{-1}\circ
f_{\gamma^{\star}}$ is an orientation preserving symmetry transformation of
$C$. Thus $\gamma^{\prime\prime}:=f_{\gamma}^{-1}\circ
f_{\gamma^{\star}}(\gamma^{\prime})$ is congruent to $\gamma^{\prime}$. If
$\gamma^{\prime\prime}\neq\gamma^{\prime}$ then on the union of $\gamma$,
$\gamma^{\star}$, $\gamma^{\prime}$ and $\gamma^{\prime\prime}$ the value of
the integral is zero. Let us regard the arc $\gamma^{\prime}$ with endpoints
$P$ and $Q$. If $\gamma^{\prime\prime}=\gamma^{\prime}$ then we have two
cases: either the endpoints $P$, $Q$ of the arc $\gamma^{\prime}$ remain fixed
or they are interchanged. In both cases, the orientation preserving congruence
$f_{\gamma}^{-1}\circ f_{\gamma^{\star}}$ fixes the segment $PQ$ and either
$f_{\gamma}^{-1}\circ f_{\gamma^{\star}}$ is the identity (which is a
contradiction) or it is a reflection at the center of $PQ$. This latter
possibility means that the value of the integral is zero in $\gamma^{\prime}$
and also in $\gamma$. It is easy to see by an inductive argument that on the
whole smooth part the integral is zero and so

(20)
$K=\int_{I_{S}}\kappa(s)ds+\sum_{i=1}^{v^{\star}}\Delta\alpha_{i}={{\sum_{i=1}^{v^{\star}}\Delta\alpha_{i}=}}2\pi.$

As $\Delta\alpha_{i}\leq\pi$, this implies that for the number $v^{\star}$ of
corners in $C$ we have $v^{\star}\geq 2$ and thus it proves the Corollary. ∎

######  Remark 7.

It is easy to see that the proof of Corollary 1 also implies that for the
normal, periodic tilings of the Euclidean plane we have $\bar{v}^{\star}\geq
2.$

###  3.2. Construction of a monohedral tiling with no vertices in 3 dimensions

In 3 dimensions the deregularization algorithm, presented in the proof of
Theorem 1, would work in an analogous manner as in 2 dimensions, since in any
dimension $d$, at any point, the number of hypercells (of dimension $d-1$)
tangent to each other is exactly 2. This algorithm reduces the nodal degree by
2, and thus we expect (15) to remain valid. However, unlike in 2 dimensions,
in the 3D case this equation will not provide the actual lower bound on
$\bar{v}^{\star}$: we have to consider that in 3 dimensions one can use a
different deregularization algorithm as well, where not vertices but edges are
being smoothed. The general formulae for the second algorithm appear to be
less trivial as smoothing of edges is not a local operation. Nevertheless, it
is not difficult to establish a sharp lower bound in 3 dimensions:

Figure 4. Construction of a $v^{\star}=0$ monohedral tiling in $d=3$
dimensions. The rectangular prism has its vertical edges at $(x,y)=\pm 1$ and
the saddle-like surface at the top and the bottom are given as
$z=\sqrt{1-x^{2}}-\sqrt{1-y^{2}}\pm 1$.

######  Lemma 1.

In $d=3$ dimensions $\bar{v}^{\star}\geq 0.$

###### Proof.

We prove Lemma 1 by showing a monohedral tiling in 3D with $v^{\star}=0$,
illustrated in Figure 4. ∎

## References

  * [1] P. M. Adler and J-F. Thovert.  Fractures and Fracture Networks.  Springer, Dordrecht, 1999.

  * [2] W. Chang.  Application of tessellation in architectural geometry design.  E3S Web of Conferences, ICEMEE 2018, 38 (03015), 2018, DOI https://doi.org/10.1051/e3sconf/20183803015.
  * [3] Gábor Domokos, Douglas J. Jerolmack, Ferenc Kun, and János Török.  Plato’s cube and the natural geometry of fragmentation.  Proceedings of the National Academy of Sciences, 117(31):18178–18185, 2020.
  * [4] L. Goehring and S. W. Morris.  Scaling of columnar joints in basalt.  Journal of Geophysical Research: Solid Earth, 113(B10), 2008.
  * [5] Lucas Goehring.  Evolving fracture patterns: columnar joints, mud cracks and polygonal terrain.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(2004):20120353, 2013.
  * [6] C. Gold.  Tessellations in gis: Part ii–making changes.  Geo-spatial Information Science, 19:157 – 167, 2016.
  * [7] B. Grünbaum and G.C. Shepard.  Tilings and Patterns.  Freeman and Co., New York, 1987.
  * [8] E. Heller and E. Fuchs.  Tissue patterning and cellular mechanics.  Journal of Cell Biology, 211:219 – 231, 2015.
  * [9] C.E. Horne.  Geometric Symmetry in Patterns and Tilings.  Woodhead Publishing, 2000.
  * [10] B. Jiang.  A topological pattern of urban street networks: Universality and peculiarity.  Physica A, 384:647 – 655, 2007.
  * [11] D. Kazanci and A. Vince.  A property of normal tilings.  Am. Math. Monthly, 111:813 – 816, 2004.
  * [12] E. Laruelle, N. Spassky, and A. Genovesio.  Unraveling spatial cellular pattern by computational tissue shuffling.  Communication Biology, 3, 2020.
  * [13] D. Schattschneider and M. Senechal.  Tilings. Chapter 3 in: Discrete and Computational Geometry.  CRC Press, 2004.
  * [14] S. Steacy and C. Sammis.  An automaton for fractal patterns of fragmentation.  Nature, 353:250–252, 1991.


# Multi-label classification with Transformers for Action Unit Detection

Gauthier Tallec
ISIR
Edouard Yvinec
ISIR/Datakalab
Arnaud Dapogny
Datakalab
Kevin Bailly
ISIR/Datakalab

###### Abstract

Action Unit (AU) Detection is the branch of affective computing that aims at
recognizing unitary facial muscular movements. It is key to unlock unbiased
computational face representations and has therefore aroused great interest in
the past few years. One of the main obstacles toward building efficient deep
learning based AU detection system is the lack of wide facial image databases
annotated by AU experts. In that extent the ABAW challenge paves the way
toward better AU detection as it involves a 2M frames AU annotated dataset. In
this paper, we present our submission to the ABAW3 challenge. In a nutshell,
we applied a multi-label detection transformer that leverage multi-head
attention to learn which part of the face image is the most relevant to
predict each AU.

ImageImage EmbeddingPositional EncodingMulti-Head AttentionAdd & NormMLP
BlockAdd & Norm$\mathbf{K}$$\mathbf{V}$$\mathbf{Q}$$N_{x}\times$Label IDsLabel
EmbeddingMulti-Head AttentionAdd & NormMulti-Head AttentionAdd & NormMLP
BlockAdd & NormLinear$N_{l}\times$SigmoidOutput
probability$\mathbf{K}$$\mathbf{V}$$\mathbf{Q}$$\mathbf{K}$$\mathbf{V}$$\mathbf{Q}$
Figure 1: Overview of the Multi-Label Transformer

##  1 Introduction

The ABAW3 Challenge [6] is the third edition [7, 13, 6] of Affect Analysis in
the wild [26]. We participated in the AU detection part of the challenge. In
previous editions, top ranked methods took benefit from the heterogeneous
annotations of ABAW3 [9, 12, 11, 10, 8] and used multi-task learning to
jointly learn the different affective computing tasks corresponding to the
different challenges. This year, the use of such strategy was restricted to
the multi-task challenge competitors. We therefore focused on the AU detection
problem specifities.

Action Units (AUs) are a dictionary of unitary muscular activations that was
designed by psychologists to anatomically describe the mechanics of facial
expressions. Theoretically [4], accurate detection of such activations enables
unbiased computational description of human faces and could therefore help
improve face analysis applications altogether.

From a machine learning standpoint, AU detection is a multi-label detection
problem that comes with three main specifities: First, AUs activations are
very local changes of skin texture. Consequently, several methods tried to
design locally adapted features [29] with a notable focus on facial landmark-
based AU attention [14, 20]. Second, AUs are heavily dependant from one
another (mostly for physiological reasons), hence multiple approaches tried to
model those dependencies using label prior dependency learning methods [24,
22]. Finally, AU are events of short duration, therefore video-based datasets,
which include all publicly available AU detection datasets, are bound to be
heavily imbalanced toward the absence of activation. This problem has mostly
been tackled by weighting the binary cross-entropy with AU frequency based
coefficients [21] and adding a dice score contribution to the loss [20].

In this paper, we design a transformer-based architecture for multi-label
detection that we apply to the action unit detection problem. More precisely,
our method learns local features adapted to each AU using a) self-attention on
encoded patches of human faces to produce the features and b) cross-attention
between those features and learned AU token for selecting which feature are
most informative about each AU. Figure 1 summarizes the architecture of our
model.

##  2 Related Work

###  2.1 Action Unit Detection.

The research work in action unit detection is centered around two main
streams, that corresponds to the two main specificities of the problem, namely
the local aspect of AUs and the strong dependencies that AUs displays between
one another.

Exploiting the local aspect of action units require the extraction of local
features. For that purpose, seminal work [29] proposed convolutional layers
with region-wise shared filters to extract features that are adapted to each
face region. Nevertheless, for such method to be efficient face parts need to
always fall in the same rectangular region. To overcome this limitation,
several works [14, 19, 20] guided AU related feature using facial landmarks,
either by learning from predefined landmark neighborhood crops [14], or by
using landmark based AU attention maps [21].

To take advantage of the strong relationships between AU, several works
focused on explicitly modeling those dependencies. In particular,
backpropagation through a probabilistic graphical model (PGM) was adopted in
[1], and an hybrid message passing strategy was used in [22]. Others leveraged
the local nature of AU to assume that label dependencies imply dependencies
between local face zones. Attempts at capturing such spatial dependencies
include attention map learning [21], LSTM-based spatial pooling [18] and more
recently transformer-like architectures [5].

###  2.2 Transformers for Multi-label Detection.

The arrival of transformers [25] has significantly modified the landscape of
state of the art deep learning architectures, first in natural language
processing [2] and more recently in computer vision [3]. So far, transformer-
based multi-label leveraged cross-attention between label-wise learnable
tokens [15] and encoded input parts to learn which portions of the inputs are
useful for predicting each label. In this work, we propose to apply this very
architecture to AU detection.

##  3 Methodology

###  3.1 Background

Most transformer architectures consists in a combination of Multi-Head
Attention modules (MHA) and Multi-Layer Perceptron (MLP) blocks interspersed
by layer normalization (missing quote) and dropout (missing quote) denoted LN
and $\mathbf{D}$ respectively.

The Multi-Head Attention module consists in an ensemble of $N_{h}$ attention
head. Each of these heads computes attention between $n_{q}$ queries
$\mathbf{q}$ stored in $\mathbf{Q}\in\mathbb{R}^{n_{q}\times d}$ and $n$ keys
$\mathbf{K}\in\mathbb{R}^{n\times d}$ for combining $n$ values
$\mathbf{V}\in\mathbb{R}^{n\times d}$. For that purpose, for each head $h$,
query, key and values are encoded using dense layers
$\mathbf{W}_{\mathbf{Q}}^{(h)}\in\mathbb{R}^{d\times
d_{k}},\mathbf{W}_{\mathbf{K}}^{(h)}\in\mathbb{R}^{d\times d_{k}}$ and
$\mathbf{W}_{\mathbf{V}}^{(h)}\in\mathbb{R}^{d\times d_{v}}$:

$\tilde{\mathbf{Q}}^{(h)}=\mathbf{Q}\mathbf{W}_{\mathbf{Q}}^{(h)},\tilde{\mathbf{K}}^{(h)}=\mathbf{\mathbf{K}}\mathbf{W}_{\mathbf{K}}^{(h)},\tilde{\mathbf{V}}^{(h)}=\mathbf{V}\mathbf{W}_{\mathbf{V}}^{(h)}.$
(1)

Then each head performs encoded values combination based on encoded key/query
comparison:

$\mathbf{H}^{(h)}=\text{softmax}\left(\frac{\tilde{\mathbf{Q}}^{(h)}(\tilde{\mathbf{K}}^{(h)})^{T}}{\sqrt{d_{k}}}\right)\tilde{\mathbf{V}}^{(h)}.$
(2)

The $N_{h}$ resulting vectors are concatenated and projected using dense
layers $\mathbf{W}_{O}\in\mathbb{R}^{N_{h}d_{v}\times d}$:

$\text{MHA}(\mathbf{Q},\mathbf{K},\mathbf{V};\mathbf{M})=\text{Concat}(\mathbf{H}^{(1)},\dots,\mathbf{H}^{(N_{h})})\mathbf{W}_{O}.$
(3)

To simplify the notations, we use $d_{k}=d_{v}=d/N_{h}$ (missing quote) and
further define self and cross attention layers as follows:

$\displaystyle\text{SA}(\mathbf{Q};\mathbf{M})$
$\displaystyle=\text{MHA}(\mathbf{Q},\mathbf{Q},\mathbf{Q})$
$\displaystyle\text{CA}(\mathbf{Q};\mathbf{K};\mathbf{M})$
$\displaystyle=\text{MHA}(\mathbf{Q},\mathbf{K},\mathbf{K})$

The MLP Block is composed of two dense layers with dropout on top of each of
them: $\mathbf{W}_{g}\in\mathbb{R}^{d\times d_{mlp}}$ with gelu activation
(missing quote) and $\mathbf{W}_{l}\in\mathbb{R}^{d_{mlp}\times d}$ with
linear activation. Formally for input queries $\mathbf{Q}$:

$\mathbf{Q}_{mlp}=\mathbf{D}(\text{GELU}(\mathbf{Q}\mathbf{W}_{g})),\text{MLP}(\mathbf{Q})=\mathbf{D}(\mathbf{Q}_{mlp}\mathbf{W}_{l}).$
(4)

###  3.2 Multi-Label Transformer

####  3.2.1 Self Attention on input image

Recent work [3] on transformer architectures unlocked the use of self
attention mechanism on images. However it has been shown that applying such
mechanism on raw images requires large training datasets with a lot of
variability which is only partially the case for ABAW3 database. For that
purpose we use a pretrained encoder $g_{\mathbf{W}_{I}}$ to first encode the
raw images $\mathbf{X}$ in a patch-based fashion:

$\tilde{\mathbf{X}}^{(0)}=g_{\mathbf{W}_{I}}(\tilde{\mathbf{X}})\in\mathbb{R}^{n_{x}\times
d},$ (5)

Second we encode positions by adding a different token to each encoded patch:

$\mathbf{X}^{(0)}=\tilde{\mathbf{X}}^{(0)}+\mathbf{I}_{n_{x}}\mathbf{W}_{n_{x}}$
(6)

where $\mathbf{W}_{n_{x}}\in\mathbb{R}^{n_{x}\times d}$ stores each patch
position embedding. Finally, we apply $N_{x}$ layers of self-attention on the
extracted patches. Layer $l$ consists in first, input self attention:

$\mathbf{X}^{(l)}_{xx}=\text{SA}_{xx}^{(l)}(\mathbf{X}^{(l)}),\tilde{\mathbf{X}}^{(l)}_{xx}=\text{LN}[\mathbf{X}^{(l)}+\mathbf{D}(\mathbf{X}^{(l)}_{xx})],$
(7)

followed by MLP block based encoding:

$\mathbf{X}^{(l)}_{mlp}=\text{MLP}^{(l)}(\tilde{\mathbf{X}}^{(l)}_{xx}),\mathbf{X}^{(l+1)}=\text{LN}[\tilde{\mathbf{X}}^{(l)}_{xx}+\mathbf{D}(\mathbf{X}^{(l)}_{mlp})].$
(8)

####  3.2.2 Cross Attention on label tokens

Previous works [15] on multi-label transformers leverage $N_{l}$ layers of
multi-head attention between learnable label tokens and input parts in order
to learn which portion of the input is the more informative about each label.
Formally, we denote
$\mathbf{T}^{(0)}=\mathbf{I}_{L}\mathbf{W}_{t}\in\mathbb{R}^{L\times d}$ the
$L$ labels ids encoded with dense layer $\mathbf{W}_{t}\in\mathbb{R}^{L\times
d}$ and $\mathbf{X}^{(N_{x})}$ the $n_{x}$ considered parts extracted from
input. Layer $l$ consists in, first, token self attention:

$\mathbf{T}^{(l)}_{tt}=\text{SA}_{tt}^{(l)}(\mathbf{T}^{(l)}),\tilde{\mathbf{T}}^{(l)}_{tt}=\text{LN}[\mathbf{T}^{(l)}+\mathbf{D}(\mathbf{T}^{(l)}_{tt})].$
(9)

Second, token image cross attention:

$\mathbf{T}^{(l)}_{tx}=\text{CA}_{tx}^{(l)}(\tilde{\mathbf{T}}_{tt}^{(l)},\mathbf{X}^{(N_{x})}),\tilde{\mathbf{T}}^{(l)}_{tx}=\text{LN}[\tilde{\mathbf{T}}^{(l)}_{tt}+\mathbf{D}(\mathbf{T}^{(l)}_{tx})],$
(10)

and finally, MLP block based encoding:

$\mathbf{T}^{(l)}_{mlp}=\text{MLP}^{(l)}(\tilde{\mathbf{T}}^{(l)}_{tx}),\mathbf{T}^{(l+1)}=\text{LN}[\tilde{\mathbf{T}}^{(l)}_{tx}+\mathbf{D}(\mathbf{T}^{(l)}_{mlp})].$
(11)

At the end of the $N_{l}$-th layer, the vector of predictions is computed by
projecting the $L$ encoded tokens using a dense layer
$\mathbf{W}_{p}\in\mathbb{R}^{d\times 1}$ with sigmoid activation $\sigma$:

$\mathbf{p}=S(\mathbf{T}^{(N_{l})}\mathbf{W}_{p})\in\mathbb{R}^{T}.$ (12)

Task $t$ distribution is then estimated as follows:

$\displaystyle\log p(y^{(t)}\mid\mathbf{x},\mathbf{W})$
$\displaystyle=-\text{BCE}(y^{(t)},p^{(t)}),$ (13) $\displaystyle=y^{(t)}\log
p^{(t)}+(1-y^{(t)})\log(1-p^{(t)}),$

where BCE stands for binary cross entropy and $\mathbf{W}$ englobes all the
network parameters.

Training is done under the assumption that tasks are independent given the
input. Therefore it consists in the minimization of the following maximum
likelihood based loss:

$\mathcal{L}(\mathbf{W})=\sum_{t=1}^{T}\text{BCE}(y^{(t)},p^{(t)}).$ (14)

Lastly, the ABAW3 challenge is based on videos, i.e. on continuous frame
sequences. As AUs are segment-level events as opposed to frame-level events
(i.e. activation of a specific AU usually lasts for at least a couple of
frames): to impose such structure in the transformer prediction, we simply
smoothed these predictions using a moving mean of window size $w_{s}=20$
frames.

##  4 Experiments

Submission ID | moving mean | frequency weights | Dice | DISFA/BP4D pretraining | mixup | Valid F1 | Test Mean F1
---|---|---|---|---|---|---|---
1 | ✗ | ✓ | ✓ | ✗ | ✓ | 52.7 | -
2 | ✓ | ✓ | ✓ | ✗ | ✓ | 53.8 | -
3 | ✓ | ✗ | ✗ | ✓ | ✗ | 50.3 | -
4 | ✓ | ✗ | ✓ | ✓ | ✗ | 51.0 | -
5 | ✓ | ✗ | ✓ | ✓ | ✓ | 51.7 | -
Table 1: Summary of submissions to the ABAW3 Challenge

###  4.1 Network Architecture

For all our experiments we use an (Imagenet-pretrained) InceptionV3 backbone
from which the pooling layer is removed and replaced by a 1D convolution with
$128$ output channels, on top of which we add a single self-attention layer
with $N_{h}=8$. As far as cross attention layers are concerned we use
$N_{l}=2$ as well as a dropout rate of $0.1$.

###  4.2 Training Strategy

Actions units detection is a heavily imbalanced problem, therefore we use a
number of solutions to mitigate this problem while learning: we weight each
example (term of the sum in Equation 14) with frequency based coefficients.
Furthermore, we add a Dice score contribution [20] to the final loss in order
to adapt to the F1Score based evaluation.

For loss optimization we use AdamW optimizer [16]. For the Inception part of
the network we use an exponentially decaying learning rate with initial value
$5e-4$ and decay rate $0.99$. For the transformer part of our network we use
the learning schedule in [25]. We scale it with respect to the number of
patches $n_{x}=64$ in the input image self attention part, and with respect to
the number of labels $L$, in the label tokens cross-attention part.

###  4.3 Pretraining and Data Augmentation

We investigate pretraining the network on already available data, e.g. BP4D
[28] and DISFA [17] databases. However, these datasets do not have exatcly the
same AU annotation as in ABAW3 dataset: Hence, we handle this by pretraining
on the reunion of the 3 different AU sets and masking the bce loss for the AU
that are not annotated for each example of each dataset.

Furthermore, we apply data-augmentation by combining geometrical (random
rotation, horizontal flips, random zoom) and color-based (channel drop, random
brightness) augmentations. Furthermore we make use of label smoothing [23] and
image mixup [27].

Table 1 reports the results of our different submissions to the challenge.

##  5 Conclusion

Conclude with results on the challenge.

## References

  * [1] Ciprian Corneanu, Meysam Madadi, and Sergio Escalera.  Deep structure inference network for facial action unit recognition.  In Proceedings of the European Conference on Computer Vision (ECCV), September 2018.
  * [2] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [3] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929, 2020.
  * [4] Rosenberg Ekman.  What the face reveals: Basic and applied studies of spontaneous expression using the Facial Action Coding System (FACS).  Oxford University Press, USA, 1997.
  * [5] Geethu Miriam Jacob and Bjorn Stenger.  Facial action unit detection with transformers.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7680–7689, 2021.
  * [6] Dimitrios Kollias.  Abaw: Valence-arousal estimation, expression recognition, action unit detection & multi-task learning challenges.  arXiv preprint arXiv:2202.10659, 2022.
  * [7] D Kollias, A Schulc, E Hajiyev, and S Zafeiriou.  Analysing affective behavior in the first abaw 2020 competition.  In 2020 15th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2020)(FG), pages 794–800.
  * [8] Dimitrios Kollias, Viktoriia Sharmanska, and Stefanos Zafeiriou.  Face behavior a la carte: Expressions, affect and action units in a single network.  arXiv preprint arXiv:1910.11111, 2019.
  * [9] Dimitrios Kollias, Viktoriia Sharmanska, and Stefanos Zafeiriou.  Distribution matching for heterogeneous multi-task learning: a large-scale face study.  arXiv preprint arXiv:2105.03790, 2021.
  * [10] Dimitrios Kollias, Panagiotis Tzirakis, Mihalis A Nicolaou, Athanasios Papaioannou, Guoying Zhao, Björn Schuller, Irene Kotsia, and Stefanos Zafeiriou.  Deep affect prediction in-the-wild: Aff-wild database and challenge, deep architectures, and beyond.  International Journal of Computer Vision, pages 1–23, 2019.
  * [11] Dimitrios Kollias and Stefanos Zafeiriou.  Expression, affect, action unit recognition: Aff-wild2, multi-task learning and arcface.  arXiv preprint arXiv:1910.04855, 2019.
  * [12] Dimitrios Kollias and Stefanos Zafeiriou.  Affect analysis in-the-wild: Valence-arousal, expressions, action units and a unified framework.  arXiv preprint arXiv:2103.15792, 2021.
  * [13] Dimitrios Kollias and Stefanos Zafeiriou.  Analysing affective behavior in the second abaw2 competition.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 3652–3660, 2021.
  * [14] Wei Li, Farnaz Abtahi, Zhigang Zhu, and Lijun Yin.  Eac-net: A region-based deep enhancing and cropping approach for facial action unit detection.  In 2017 12th IEEE International Conference on Automatic Face & Gesture Recognition (FG 2017), pages 103–110. IEEE, 2017.
  * [15] Shilong Liu, Lei Zhang, Xiao Yang, Hang Su, and Jun Zhu.  Query2label: A simple transformer way to multi-label classification.  arXiv preprint arXiv:2107.10834, 2021.
  * [16] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  In International Conference on Learning Representations, 2019.
  * [17] S Mohammad Mavadati, Mohammad H Mahoor, Kevin Bartlett, Philip Trinh, and Jeffrey F Cohn.  Disfa: A spontaneous facial action intensity database.  IEEE Transactions on Affective Computing, 4(2):151–160, 2013.
  * [18] Xuesong Niu, Hu Han, Songfan Yang, Yan Huang, and Shiguang Shan.  Local relationship learning with person-specific shape regularization for facial action unit detection.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 11917–11926, 2019.
  * [19] Zhiwen Shao, Zhilei Liu, Jianfei Cai, and Lizhuang Ma.  Deep adaptive attention for joint facial action unit detection and face alignment.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 705–720, 2018.
  * [20] Zhiwen Shao, Zhilei Liu, Jianfei Cai, and Lizhuang Ma.  Jaa-net: Joint facial action unit detection and face alignment via adaptive attention.  International Journal of Computer Vision, 129(2):321–340, 2021\.
  * [21] Zhiwen Shao, Zhilei Liu, Jianfei Cai, Yunsheng Wu, and Lizhuang Ma.  Facial action unit detection using attention and relation learning.  IEEE Transactions on Affective Computing, 2019.
  * [22] Tengfei Song, Zijun Cui, Wenming Zheng, and Qiang Ji.  Hybrid message passing with performance-driven structures for facial action unit detection.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6267–6276, 2021.
  * [23] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2818–2826, 2016.
  * [24] Gauthier Tallec, Arnaud Dapogny, and Kevin Bailly.  Multi-order networks for action unit detection.  arXiv preprint arXiv:2202.00446, 2022.
  * [25] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [26] Stefanos Zafeiriou, Dimitrios Kollias, Mihalis A Nicolaou, Athanasios Papaioannou, Guoying Zhao, and Irene Kotsia.  Aff-wild: Valence and arousal ‘in-the-wild’challenge.  In Computer Vision and Pattern Recognition Workshops (CVPRW), 2017 IEEE Conference on, pages 1980–1987. IEEE, 2017.
  * [27] Hongyi Zhang, Moustapha Cisse, Yann N. Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  In International Conference on Learning Representations, 2018.
  * [28] Xing Zhang, Lijun Yin, Jeffrey F Cohn, Shaun Canavan, Michael Reale, Andy Horowitz, and Peng Liu.  A high-resolution spontaneous 3d dynamic facial expression database.  In 2013 10th IEEE international conference and workshops on automatic face and gesture recognition (FG), pages 1–6. IEEE, 2013.
  * [29] Kaili Zhao, Wen-Sheng Chu, and Honggang Zhang.  Deep region and multi-label learning for facial action unit detection.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3391–3399, 2016.


# Generalized thermo vacuum state derived by the partial trace method
††thanks: Worked supported by the National Natural Science Foundation of China
under Grant 10775097.

Li-yun Hu1,2 and Hong-yi Fan1 Corresponding author. E-mail<EMAIL_ADDRESS>1Department of Physics, Shanghai Jiao Tong University, Shanghai 200030, China
2College of Physics & Communication Electronics, Jiangxi Normal University,
Nanchang 330022, China

###### Abstract

By virtue of the technique of integration within an ordered product (IWOP) of
operators we present a new approach for deriving generalized thermo vacuum
state which is simpler in form that the result by using the Umezawa-Takahashi
approach, in this way the thermo field dynamics can be developed. Applications
of the new state are discussed.

Keywords: the partial trace method, generalized thermal vacuum state, thermo
field dynamics, the IWOP technique

PACS numbers: 03.65.-w, 42.50.-p

##  I Introduction

In nature every system is immersed in an environment, the problem about the
system interacting with the environment is a hot topic in quantum information
and quantum optics. To describe system-environment evolution Takahashi-Umezawa
introduced thermo field dynamics (TFD) 1 ; 2 , with which one may convert the
calculations of ensemble averages at finite temperature

$\left\langle
A\right\rangle=\mathtt{tr}\left(A\rho\right)/Z\left(\beta\right),\text{
}\rho=e^{-\beta H},$ (1)

to equivalent expectation values with a pure state
$\left|0(\beta)\right\rangle$, i.e.,

$\left\langle A\right\rangle=\left\langle
0(\beta)\right|A\left|0(\beta)\right\rangle,$ (2)

where $\beta=1/kT$, $k$ is the Boltzmann constant, and
$Z\left(\beta\right)=\mathtt{tr}\rho=\mathtt{tr}e^{-\beta H}$ is the partition
function; $H$ is the system’s Hamiltonian. Then how to find the explicit
expression of $\left|0(\beta)\right\rangle?$ If one expands
$\left|0(\beta)\right\rangle$ in terms of the energy eigenvector set of $H$,
$\left|0(\beta)\right\rangle=\sum_{n}\left|n\right\rangle f_{n}(\beta),$ and
then substituting it into Eq.(2), which results in
$f_{n}^{\ast}(\beta)f_{m}(\beta)=Z^{-1}\left(\beta\right)e^{-\beta
E_{n}}\delta_{n,m}$ (after comparing with Eq.(1)). By introducing a fictitious
mode, $\left\langle\tilde{n}\right.\left|\tilde{m}\right\rangle=\delta_{n,m},$
then Takahashi-Umezawa obtained

$\left|0(\beta)\right\rangle=Z^{-1/2}\left(\beta\right)\sum_{n}e^{-\beta
E_{n}/2}\left|n,\tilde{n}\right\rangle.$ (3)

Thus the worthwhile convenience in Eq.(2) is at the expense of introducing a
fictitious field (or called a tilde-conjugate field, denoted as operator
$\tilde{a}^{\dagger}$) in the extended Hilbert space, i.e., the original
optical field state $\left|n\right\rangle$ in the Hilbert space $\mathcal{H}$
is accompanied by a tilde state $\left|\tilde{n}\right\rangle$ in
$\mathcal{\tilde{H}}$. A similar rule holds for operators: every Bose
annihilation operator $a$ acting on $\mathcal{H}$ has an image $\tilde{a}$
acting on $\mathcal{\tilde{H}}$,
$\left[\tilde{a},\tilde{a}^{\dagger}\right]=1$. These operators in
$\mathcal{H}$ are commutative with those in $\mathcal{\tilde{H}}$.

For a harmonic oscillator the Hamiltonian is $\hbar\omega a^{\dagger}a,$
$\left|n\right\rangle=a^{\dagger n}/\sqrt{n!}\left|0\right\rangle,$Takahashi-
Umezawa obtained the explicit expression of $\left|0(\beta)\right\rangle$ in
doubled Fock space,

$\left|0(\beta)\right\rangle=\text{sech}\theta\exp\left[a^{\dagger}\tilde{a}^{\dagger}\tanh\theta\right]\left|0\tilde{0}\right\rangle=S\left(\theta\right)\left|0\tilde{0}\right\rangle,$
(4)

which is named thermo vacuum state, and $S\left(\theta\right)$ thermo
operator,

$S\left(\theta\right)\equiv\exp\left[\theta\left(a^{\dagger}\tilde{a}^{\dagger}-a\tilde{a}\right)\right],$
(5)

which is similar in form to the a two-mode squeezing operator except for the
tilde mode. $\theta$ is a parameter related to the temperature by
$\tanh\theta=\exp\left(-\frac{\hbar\omega}{2kT}\right).$

An interesting question thus challenges us: For the Hamiltonian being
$H=\omega a^{\dagger}a+\kappa^{\ast}a^{\dagger 2}+\kappa a^{2},$ then what is
the corresponding thermo vacuum state? One may wonder if this question is
worth of paying attention since this $H$ can be diagonalized by the Bogoliubov
transformation as a new harmonic oscillator, correspondingly, the thermo
vacuum state for $H$ can be obtained by acting the same transformation on
$\left|0(\beta)\right\rangle$ in (4) (see Eq. (A9) in the Appendix). To make
this issue worthwhile, we emphasize that we shall adopt a completely new
approach to construct thermo vacuum state and our result is simpler in form
than that in Eq. (A10). Our work is arranged as follows. In Sec. 2 by re-
analyzing Eqs. (1)-(2) we shall introduce a new method (the partial trace
method) to find the explicit expression of $\left|0(\beta)\right\rangle$ in
(4). Then using this method, we obtain the expression of
$\left|0(\beta)\right\rangle$ in Eq. (4) in Sec. 3. For the degenerate
parametric amplifier, we derive a generalized thermal vacuum state
$\left|\phi(\beta)\right\rangle$ in Sec. 4. Section 5 is devoting to
presenting some applications of $\left|\phi(\beta)\right\rangle$.

##  II The partial trace method

Following the spirit of TFD, for a density operator $\rho=e^{-\beta
H}/Z\left(\beta\right)$ with Hamiltonian $H$, we can suppose that the ensemble
averages of a system operator $A$ may be calculated as

$A=\mathtt{tr}\left(\rho
A\right)=\left\langle\psi(\beta)\right|A\left|\psi(\beta)\right\rangle,$ (6)

where $\left|\psi(\beta)\right\rangle$ corresponds to the pure state in the
extended Hilbert space.

Let $\mathtt{Tr}$ denote the trace operation over both the system freedom
(expressed by $\mathtt{tr}$) and the environment freedom by
$\widetilde{\mathtt{tr}}$, i.e.,
$\mathtt{Tr}=\mathtt{tr}\widetilde{\mathtt{tr}}$, then we have

$\displaystyle\left\langle\psi(\beta)\right|A\left|\psi(\beta)\right\rangle$
$\displaystyle=$
$\displaystyle\mathtt{Tr}\left[A\left|\psi(\beta)\right\rangle\left\langle\psi(\beta)\right|\right]$
(7) $\displaystyle=$
$\displaystyle\mathtt{tr}\left[A\widetilde{\mathtt{tr}}\left|\psi(\beta)\right\rangle\left\langle\psi(\beta)\right|\right].$

Note that

$\widetilde{\mathtt{tr}}\left|\psi(\beta)\right\rangle\left\langle\psi(\beta)\right|\neq\left\langle\psi(\beta)\right|\left.\psi(\beta)\right\rangle,$
(8)

since $\left|\psi(\beta)\right\rangle$ involves both real mode and fictitious
mode. Comparing Eq.(7) with Eq.(1) we see

$\widetilde{\mathtt{tr}}\left|\psi(\beta)\right\rangle\left\langle\psi(\beta)\right|=e^{-\beta
H}/Z\left(\beta\right).$ (9)

Eq.(9) indicates that, for a given Hamiltonian $H$, if we can find a density
operator of pure state $\left|\psi(\beta)\right\rangle$ in doubled Hilbert
space, whose partial trace over the tilde freedom may lead to density operator
$e^{-\beta H}/Z\left(\beta\right)$ of the system, then the average value of
operator $A$ can be calculated as an equivalent expectation value with a pure
state $\left|\psi(\beta)\right\rangle,$ i.e., $\left\langle
A\right\rangle=\mathtt{tr}\left(Ae^{-\beta
H}/Z\left(\beta\right)\right)=\left\langle\psi(\beta)\right|A\left|\psi(\beta)\right\rangle$.

In particular, when $H=\hbar\omega a^{\dagger}a,$ a free Bose system, Eq. (9)
becomes

$\widetilde{\mathtt{tr}}\left|0(\beta)\right\rangle\left\langle
0(\beta)\right|=\left(1-e^{-\beta\hbar\omega}\right)e^{-\beta\hbar\omega
a^{\dagger}a}\equiv\rho_{c},$ (10)

$\rho_{c}$ is the density operator of chaotic field. This equation enlightens
us to have a new approach for deriving $\left|0(\beta)\right\rangle\colon$
$\left|0(\beta)\right\rangle\left\langle 0(\beta)\right|$ in doubled Hilbert
space should be such constructed that its partial trace over the tilde freedom
may lead to density operator $\rho_{c}$ of the system. In the following we
shall employ the technique of integration within an ordered product (IWOP) of
operators 3 ; 4 ; 5  to realize this goal.

##  III Derivation of $\left|0(\beta)\right\rangle$ in Eq.(4) via the new
approach

Using the normally ordered expansion formula 6

$e^{-\beta\hbar\omega
a^{\dagger}a}=\colon\exp\left\\{\left(e^{-\beta\hbar\omega}-1\right)a^{\dagger}a\right\\}\colon,$
(11)

(where the symbol $\colon\colon$ denotes the normal ordering form of
operator), and the IWOP technique we have

$\displaystyle\colon\exp\\{\left(e^{-\beta\hbar\omega}-1\right)a^{\dagger}a\\}\colon$
(12) $\displaystyle=$ $\displaystyle\int\frac{d^{2}z}{\pi}\colon
e^{-\left|z\right|^{2}+z^{\ast}a^{{\dagger}}e^{-\beta\hbar\omega/2}+zae^{-\beta\hbar\omega/2}-a^{{\dagger}}a}\colon.$

Remembering the ordering form of vacuum projector operator
$\left|0\right\rangle\left\langle 0\right|=\colon e^{-a^{{\dagger}}a}\colon$,
we can rewrite Eq.(12) as

$\displaystyle\colon\exp\\{\left(e^{-\beta\hbar\omega}-1\right)a^{\dagger}a\\}\colon$
(13) $\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{z^{\ast}a^{{\dagger}}e^{-\beta\hbar\omega/2}}\left|0\right\rangle\left\langle
0\right|e^{zae^{-\beta\hbar\omega/2}}\left\langle\tilde{z}\right.\left|\tilde{0}\right\rangle\left\langle\tilde{0}\right.\left|\tilde{z}\right\rangle,$

where $\left|\tilde{z}\right\rangle$ is the coherent state 7 ; 8  in fictitous
mode

$\text{ \ \
}\left|\tilde{z}\right\rangle=\exp\left(z\tilde{a}^{\dagger}-z^{\ast}\tilde{a}\right)\left|\tilde{0}\right\rangle,\tilde{a}\left|\tilde{z}\right\rangle=z\left|\tilde{z}\right\rangle,\left\langle\tilde{0}\right.\left|\tilde{z}\right\rangle=e^{-\left|z\right|^{2}/2\text{\
}}.$ (14)

Further, multipling the factor $\left(1-e^{-\beta\hbar\omega}\right)$ to both
sides of Eq.(13) and using the completeness of coherent state
$\int\frac{d^{2}z}{\pi}\left|\tilde{z}\right\rangle\left\langle\tilde{z}\right|=1$,
we have

$\displaystyle\left(1-e^{-\beta\hbar\omega}\right)\times\text{ Eq}.(\ref{13})$
(15) $\displaystyle=$
$\displaystyle\left(1-e^{-\beta\hbar\omega}\right)\int\frac{d^{2}z}{\pi}\text{
}\left\langle\tilde{z}\right|e^{z^{\ast}a^{{\dagger}}e^{-\beta\hbar\omega/2}}\left|0\tilde{0}\right\rangle\left\langle
0\tilde{0}\right|e^{zae^{-\beta\hbar\omega/2}}\left|\tilde{z}\right\rangle$
$\displaystyle=$
$\displaystyle\left(1-e^{-\beta\hbar\omega}\right)\int\frac{d^{2}z}{\pi}\text{
}\left\langle\tilde{z}\right|e^{a^{{\dagger}}\tilde{a}^{{\dagger}}e^{-\beta\hbar\omega/2}}\left|0\tilde{0}\right\rangle\left\langle
0\tilde{0}\right|e^{a\tilde{a}e^{-\beta\hbar\omega/2}}\left|\tilde{z}\right\rangle$
$\displaystyle=$
$\displaystyle\widetilde{\mathtt{tr}}\left[0(\beta)\left\langle
0(\beta)\right|\right],$

where

$\left|0(\beta)\right\rangle=\sqrt{1-e^{-\beta\hbar\omega}}\exp\left[a^{\dagger}\tilde{a}^{\dagger}e^{-\beta\hbar\omega/2}\right]\left|0\tilde{0}\right\rangle,$
(16)

which is the same as Eq.(4). Thus, according to Eq.(10), from the chaotic
field operator we have derived the thermo vacuum state, this is a new
approach, which has been overlooked in the literature before.

##  IV Generalized thermo vacuum state
$\left|\phi\left(\beta\right)\right\rangle$

Now, we consider a degenerate parametric amplifier whose Hamiltonian is

$H=\omega a^{\dagger}a+\kappa^{\ast}a^{\dagger 2}+\kappa a^{2},$ (17)

whose normalized density operator $\rho$ is

$\rho\left(\mathtt{tr}e^{-\beta H}\right)=e^{-\beta H}=e^{-\beta\left(\omega
a^{\dagger}a+\kappa^{\ast}a^{\dagger 2}+\kappa a^{2}\right)}.$ (18)

Recalling that
$\frac{1}{2}\left(a^{\dagger}a+\frac{1}{2}\right),\frac{1}{2}a^{\dagger 2}$
and $\frac{1}{2}a^{2}$ obey the SU(1,1) Lie algebra, thus we can derive a
generalized identity of operator 9 ; 10  as follows:

$\displaystyle\exp\left[fa^{\dagger}a+ga^{\dagger 2}+ka^{2}\right]$ (19)
$\displaystyle=$ $\displaystyle e^{-f/2}e^{\frac{ga^{\dagger
2}}{\mathcal{D}\coth\mathcal{D}-f}}e^{\left(a^{\dagger}a+\frac{1}{2}\right)\ln\frac{\mathcal{D}\text{sech}\mathcal{D}}{\mathcal{D}-f\tanh\mathcal{D}}}e^{\frac{ka^{2}}{\mathcal{D}\coth\mathcal{D}-f}},$

where we have set $\mathcal{D}^{2}=f^{2}-4kg.$ Thus Comparing Eq.(18) with
Eq.(19) we can recast Eq.(18) into the following form

$\left(\mathtt{tr}e^{-\beta H}\right)\rho=\sqrt{\lambda
e^{\beta\omega}}\exp\left[E^{\ast}a^{\dagger
2}\right]\exp\left[a^{\dagger}a\ln\lambda\right]\exp\left[Ea^{2}\right],$ (20)

where we have set

$\displaystyle D^{2}$ $\displaystyle=$
$\displaystyle\omega^{2}-4\left|\kappa\right|^{2},$ $\displaystyle\lambda$
$\displaystyle=$ $\displaystyle\frac{D}{\omega\sinh\beta D+D\cosh\beta D},$
(21) $\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{-\lambda}{D}\kappa\sinh\beta D.$

Further, using the formula in Eqs.(11) and (13), we have

$\displaystyle\left(\mathtt{tr}e^{-\beta H}\right)\rho$ (22) $\displaystyle=$
$\displaystyle\sqrt{\lambda e^{\beta\omega}}\exp\left[E^{\ast}a^{\dagger
2}\right]\colon\exp\left\\{\left(\lambda-1\right)a^{\dagger}a\right\\}\colon\exp\left[Ea^{2}\right]$
$\displaystyle=$ $\displaystyle\sqrt{\lambda
e^{\beta\omega}}\int\frac{d^{2}z}{\pi}e^{E^{\ast}a^{\dagger
2}+\sqrt{\lambda}z^{\ast}a^{{\dagger}}}\left|0\right\rangle\left\langle
0\right|e^{Ea^{2}+\sqrt{\lambda}za}\left\langle\tilde{z}\right.\left|\tilde{0}\right\rangle\left\langle\tilde{0}\right.\left|\tilde{z}\right\rangle$
$\displaystyle=$ $\displaystyle\sqrt{\lambda
e^{\beta\omega}}\int\frac{d^{2}z}{\pi}\left\langle\tilde{z}\right|e^{E^{\ast}a^{\dagger
2}+\sqrt{\lambda}z^{\ast}a^{{\dagger}}}\left|0\tilde{0}\right\rangle\left\langle
0\tilde{0}\right.e^{Ea^{2}+\sqrt{\lambda}za}\left|\tilde{z}\right\rangle$
$\displaystyle=$ $\displaystyle\sqrt{\lambda
e^{\beta\omega}}\widetilde{\mathtt{tr}}\left[e^{E^{\ast}a^{\dagger
2}+\sqrt{\lambda}a^{{\dagger}}\tilde{a}^{{\dagger}}}\left|0\tilde{0}\right\rangle\left\langle
0\tilde{0}\right.e^{Ea^{2}+\sqrt{\lambda}a\tilde{a}}\right]$
$\displaystyle\equiv$ $\displaystyle\left(\mathtt{tr}e^{-\beta
H}\right)\widetilde{\mathtt{tr}}\left[\left|\phi\left(\beta\right)\right\rangle\left\langle\phi\left(\beta\right)\right|\right],$

which indicates that the pure state in doubled Fock space for Hamiltonian in
(17) can be considered as

$\left|\phi\left(\beta\right)\right\rangle=\sqrt{\frac{\lambda^{1/2}e^{\beta\omega/2}}{Z\left(\beta\right)}}e^{E^{\ast}a^{\dagger
2}+\sqrt{\lambda}a^{{\dagger}}\tilde{a}^{{\dagger}}}\left|0\tilde{0}\right\rangle,$
(23)

where the partition function $Z\left(\beta\right)$ is determined by

$Z\left(\beta\right)=\mathtt{tr}e^{-\beta H}=\mathtt{tr}\left\\{\sqrt{\lambda
e^{\beta\omega}}e^{E^{\ast}a^{\dagger 2}}\colon
e^{\left(\lambda-1\right)a^{\dagger}a}\colon e^{Ea^{2}}\right\\}.$ (24)

Using $\int\frac{d^{2}z}{\pi}\left|z\right\rangle\left\langle z\right|=1$ and
the integral formula 11

$\displaystyle\int\frac{d^{2}z}{\pi}\exp\left(\zeta\left|z\right|^{2}+\xi
z+\eta z^{\ast}+fz^{2}+gz^{\ast 2}\right)$ (25) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\zeta^{2}-4fg}}\exp\left[\frac{-\zeta\xi\eta+\xi^{2}g+\eta^{2}f}{\zeta^{2}-4fg}\right],$

whose convergent condition is Re$\left(\zeta\pm f\pm g\right)<0,\
$Re$\left(\frac{\zeta^{2}-4fg}{\zeta\pm f\pm g}\right)<0,$ we can get

$Z\left(\beta\right)=\sqrt{\frac{\lambda
e^{\beta\omega}}{\left(1-\lambda\right)^{2}-4\left|E\right|^{2}}}=\frac{e^{\beta\omega/2}}{2\sinh\left(\beta
D/2\right)}.$ (26)

Thus the normalized state for Eq.(18) in doubled Fock space is given by

$\left|\phi\left(\beta\right)\right\rangle=\sqrt{2\lambda^{1/2}\sinh\left(\beta
D/2\right)}e^{E^{\ast}a^{\dagger
2}+\sqrt{\lambda}a^{{\dagger}}\tilde{a}^{{\dagger}}}\left|0\tilde{0}\right\rangle,$
(27)

and the internal energy of system is

$\left\langle H\right\rangle_{e}=-\frac{\partial}{\partial\beta}\ln
Z\left(\beta\right)=\frac{D\coth\left(\beta D/2\right)-\omega}{2},$ (28)

which leads to the distribution of entropy

$\displaystyle S$ $\displaystyle=$
$\displaystyle-k\mathtt{tr}\left(\rho\ln\rho\right)=\frac{1}{T}\left\langle
H\right\rangle_{e}+k\ln Z\left(\beta\right)$ (29) $\displaystyle=$
$\displaystyle\frac{D}{2T}\coth\left(\beta
D/2\right)-k\ln\left[2\sinh\left(\beta D/2\right)\right].$

In particular, when $\kappa=0,$ leading to $D=\omega,$ so Eq.(27) just reduces
to $\left|0(\beta)\right\rangle$ with $\omega\rightarrow\hbar\omega$, and
Eqs.(28) and (29) become
$\frac{\omega}{2}\left(\coth\left(\beta\omega/2\right)-1\right)$ and
$\frac{\omega}{2T}\coth\left(\beta\omega/2\right)-k\ln\left[2\sinh\left(\beta\omega/2\right)\right],$
respectively, as expected 12 . Thus by virtue of the technique of IWOP we can
display the partial trace method to deduce the pure state representation for
some new density operators of light field at finite temperature.

##  V Applications of generalized thermo vacuum state

###  V.1 Internal energy distribution of the system

As an application of Eq.(27), we can evaluate the each term’s contribution to
energy in Hamiltonian. Based on the idea from Eq.(1) to (2), the system
operator $A$ can be calculated as $\left\langle
A\right\rangle_{e}=\left\langle\phi\left(\beta\right)\right|A\left|\phi\left(\beta\right)\right\rangle.$
Thus uisng the completeness of coherent state and the integral formula Eq.(25)
as well as noticing
$\left\langle\phi\left(\beta\right)\right.\left|\phi\left(\beta\right)\right\rangle=1$,
$\left(1-\lambda\right)^{2}-4\left|E\right|^{2}=4\lambda\sinh^{2}\left(\beta
D/2\right),$ then we have

$\displaystyle\left\langle\omega a^{\dagger}a\right\rangle_{e}$
$\displaystyle=$
$\displaystyle\omega\left\langle\phi\left(\beta\right)\right|\left(aa^{\dagger}-1\right)\left|\phi\left(\beta\right)\right\rangle$
(30) $\displaystyle=$ $\displaystyle 2\omega\lambda^{1/2}\sinh\left(\beta
D/2\right)\frac{\partial}{\partial\lambda}$
$\displaystyle\times\int\frac{d^{2}z}{\pi}e^{-\left(1-\lambda\right)\left|z\right|^{2}+Ez^{2}+E^{\ast}z^{\ast
2}}-\omega$ $\displaystyle=$
$\displaystyle\lambda^{1/2}\frac{\partial}{\partial\lambda}\frac{2\omega\sinh\left(\beta
D/2\right)}{\sqrt{\left(1-\lambda\right)^{2}-4\left|E\right|^{2}}}-\omega$
$\displaystyle=$
$\displaystyle\frac{\omega}{2}\left(\frac{\omega}{D}\coth\beta D/2-1\right),$

and

$\displaystyle\left\langle\kappa^{\ast}a^{\dagger 2}\right\rangle_{e}$
$\displaystyle=$ $\displaystyle 2\kappa^{\ast}\lambda^{1/2}\sinh\left(\beta
D/2\right)\frac{\partial}{\partial
E^{\ast}}\frac{1}{\sqrt{\left(1-\lambda\right)^{2}-4EE^{\ast}}}$ (31)
$\displaystyle=$
$\displaystyle-\frac{\left|\kappa\right|^{2}}{D}\allowbreak\coth\beta D/2,$

as well as

$\left\langle\kappa
a^{2}\right\rangle_{e}=-\frac{\left|\kappa\right|^{2}}{D}\allowbreak\coth\beta
D/2.$ (32)

From Eqs.(31) and (32) we see that the two items ($\kappa^{\ast}a^{\dagger 2}$
and $\kappa a^{2})$ have the same energy contributions to the system, as
expected. Combing Eqs.(30)-(32) we can also check Eq.(28).

###  V.2 Wigner function and quantum tomogram

The Wigner function plays an important role in studying quantum optics and
quantum statistics 13 ; 14 . It gives the most analogous description of
quantum mechanics in the phase space to classical statistical mechanics of
Hamilton systems and is also a useful measure for studying the nonclassical
features of quantum states. In addition, the Wigner function can be
reconstructed by measuring several quadratures
$P\left(\hat{x}_{\theta}=\hat{x}\cos\theta+\hat{p}\sin\theta\right)$ with a
homodyne detection and then applying an inverse Radon transform—quantum
homodyne tomography 15 . Using Eq.(26) one can calculate conveniently the
Wigner function and quantum tomogram. Recalling that the single-mode Wigner
operator $\Delta\left(z\right)$ in coherent state representation is given by
16 ; 17

$\Delta\left(\alpha\right)=e^{2\left|\alpha\right|^{2}}\int\frac{d^{2}z}{\pi^{2}}\left|z\right\rangle\left\langle-z\right|e^{-2\left(z\alpha^{\ast}-z^{\ast}\alpha\right)},$
(33)

thus the Wigner function is

$\displaystyle W\left(\alpha\right)$ $\displaystyle=$
$\displaystyle\left\langle\phi\left(\beta\right)\right|\Delta\left(\alpha\right)\left|\phi\left(\beta\right)\right\rangle$
(34) $\displaystyle=$ $\displaystyle
e^{2\left|\alpha\right|^{2}}\int\frac{d^{2}z}{\pi^{2}}\left\langle\phi\left(\beta\right)\right.\left|z\right\rangle\left\langle-z\right.\left|\phi\left(\beta\right)\right\rangle
e^{-2\left(z\alpha^{\ast}-z^{\ast}\alpha\right)}$ $\displaystyle=$
$\displaystyle\frac{\tanh\left(\beta
D/2\right)}{\pi}e^{-\frac{2}{D}\left[\omega\left|\alpha\right|^{2}+\left(\kappa\alpha^{2}+\kappa^{\ast}\alpha^{\ast
2}\right)\right]\tanh\left(\beta D/2\right)},$

where we have noticed
$\left(1+\lambda\right)^{2}-4\left|E\right|^{2}=\frac{4D\cosh^{2}\beta
D/2}{\omega\sinh\beta D+D\cosh\beta D}$ and used the integral formula (25). In
particular, when $\kappa=0,$ Eq.(34) reduces to

$W\left(\alpha\right)=\frac{\tanh\left(\beta
D/2\right)}{\pi}\exp\left\\{-2\left|\alpha\right|^{2}\tanh\left(\beta
D/2\right)\right\\},$ (35)

which is just the Wigner function of thermo vacuum state
$\left|0\left(\beta\right)\right\rangle$.

On the other hand, we can derive the tomography (Radon transform of Wigner
function) of the system by using Eq.(26). Recalling that, for single-mode
case, the Radon transform of the Wigner operator is just a pure state density
operator 18 ,

$\int\delta\left(q-fq^{\prime}-gp^{\prime}\right)\Delta\left(\alpha\right)\mathtt{d}q^{\prime}\mathtt{d}p^{\prime}=\left|q\right\rangle_{f,g\text{
}f,g}\left\langle q\right|,$ (36)

where $\alpha=(q+\mathtt{i}p)/\sqrt{2},$ and $\left(f,g\right)$ are real,

$\left|q\right\rangle_{f,g}=C\exp\left[\frac{\sqrt{2}}{A}qa^{{\dagger}}-\frac{e^{\mathtt{i}2\varphi}}{2}a^{{\dagger}2}\right]\left|0\right\rangle,$
(37)

and
$C=\left[\pi\left(f^{2}+g^{2}\right)\right]^{-1/4}\exp\\{-q^{2}/[2(f^{2}+g^{2})]\\},A=f-\mathtt{i}g=\sqrt{f^{2}+g^{2}}e^{-\mathtt{i}\varphi}.$
Eq.(37) is named as the intermediate coordinate-momentum representation 18 .
From Eq.(36) and Eq.(26) it then follows that the tomogram can be calculated
as

$\mathcal{R}\left(q\right)_{f,g}\equiv\left\langle\phi\left(\beta\right)\right.\left|q\right\rangle_{f,g\text{
}f,g}\left\langle
q\right.\left|\phi\left(\beta\right)\right\rangle=\int\frac{d^{2}z}{\pi}\left|{}_{f,g}\left\langle
q\right|\left\langle\tilde{z}\right.\left|\phi\left(\beta\right)\right\rangle\right|^{2}.$
(38)

Then submitting Eqs.(37) and (26) into Eq.(38), we obtain

$\displaystyle\mathcal{R}\left(q\right)_{f,g}$ $\displaystyle=$
$\displaystyle\frac{2\sinh\left(\beta
D/2\right)}{C^{-2}\sqrt{\lambda+\left|G\right|^{2}/\lambda}}$ (39)
$\displaystyle\times\exp\left\\{2q^{2}\left[\frac{1-\lambda\text{Re}G^{-1}}{\left|A\right|^{2}\left(\lambda+\left|G\right|^{2}/\lambda\right)}+\text{Re}\frac{2E}{A^{2}G}\right]\right\\},$

where we have used Eq.(25) and set $G=1+2e^{\mathtt{i}2\varphi}E$. Eq.(39) is
the positive-definite tomogram, as expected. As far as we are concerned, this
result has not been reported in the literature before.

In sum, by virtue of the technique of integration within an ordered product
(IWOP) of operators we have presented a new approach for deriving generalized
thermo vacuum state which is simpler in form that the result by using the
Umezawa-Takahashi approach, in this way the thermo field dynamics can be
developed.

Appendix:

As a comparison of our new approach with the usual way of deriving thermo
vacuum state in TFD theory, in this appendix, we shall derive the explicit
expression of $\left|\phi\left(\beta\right)\right\rangle$ by diagonalizing
Hamiltonian (17). For this purpose, we introduce two unitary operators: one is
a single mode squeezing operator,

$S=\exp\left(\frac{\nu}{\mu}\frac{a^{{\dagger}2}}{2}\right)\exp\left[\left(a^{{\dagger}}a+\frac{1}{2}\right)\ln\frac{1}{\mu}\right]\exp\left(-\frac{\nu}{\mu}\frac{a^{2}}{2}\right),$
(A1)

where $\mu$ and $\nu$ are squeezing parameters satisfying the unitary-modulate
condition $\mu^{2}-\nu^{2}=1$; and the other is a rotational operator,
$R=\exp\left(\frac{i\phi}{2}a^{\dagger}a\right),$ which lead to the following
transformations,

$\displaystyle SaS^{{\dagger}}$ $\displaystyle=\mu a-\nu a^{{\dagger}},\text{
}Sa^{{\dagger}}S^{{\dagger}}=\mu a^{{\dagger}}-\nu a,$ (A2) $\displaystyle
S^{{\dagger}}aS$ $\displaystyle=\mu a+\nu a^{{\dagger}},\text{
}S^{{\dagger}}a^{{\dagger}}S=\mu a^{{\dagger}}+\nu a,$ (A3)

and

$RaR^{{\dagger}}=ae^{-\frac{i\phi}{2}},\text{
}Ra^{{\dagger}}R^{{\dagger}}=a^{{\dagger}}e^{\frac{i\phi}{2}}.$ (A4)

Thus, under the unitary transform $SR$, we have (setting
$\kappa=\left|\kappa\right|e^{i\phi}$)

$\displaystyle H^{\prime}$ $\displaystyle=SRHR^{{\dagger}}S^{{\dagger}}=\omega
Sa^{\dagger}aS^{{\dagger}}+\left|\kappa\right|Sa^{\dagger
2}S^{{\dagger}}+\left|\kappa\right|Sa^{2}S^{{\dagger}}$
$\displaystyle=\left(\omega\mu^{2}+\omega\nu^{2}-4\left|\kappa\right|\mu\nu\right)a^{{\dagger}}a+\left(\omega\nu-2\left|\kappa\right|\mu\right)\nu$
$\displaystyle+\left(\left|\kappa\right|\mu^{2}+\left|\kappa\right|\nu^{2}-\omega\mu\nu\right)\left(a^{{\dagger}2}+a^{2}\right).$
(A5)

In order to diagonalize Eq.(A5), noticing $\mu^{2}-\nu^{2}=1$ and making
$\left|\kappa\right|\left(\mu^{2}+\nu^{2}\right)-\omega\mu\nu=0,$ whose
solution is given by

$\mu^{2}=\frac{\omega}{2\omega^{\prime}}+\frac{1}{2},\nu^{2}=\allowbreak\frac{\omega}{2\omega^{\prime}}-\frac{1}{2},\omega^{\prime}=\sqrt{\omega^{2}-4\left|\kappa\right|^{2}}.$
(A6)

then Eq.(A5) becomes

$H^{\prime}=\omega^{\prime}\left(a^{\dagger}a+\frac{1}{2}\right)-\frac{1}{2}\omega.$
(A7)

i.e., the diagonalization of Hamiltonian is completed.

According to Eq.(16), the thermal vacuum state corresponding to density
operator $\rho^{\prime}=e^{-\beta H^{\prime}}/\mathtt{tr}\left(e^{-\beta
H^{\prime}}\right)=e^{-\beta\omega^{\prime}a^{\dagger}a}/\mathtt{tr}\left(e^{-\beta\omega^{\prime}a^{\dagger}a}\right)$
is given by

$\left|0(\beta)\right\rangle=\sqrt{1-e^{-\beta\omega^{\prime}}}\exp\left[a^{\dagger}\tilde{a}^{\dagger}e^{-\beta\omega^{\prime}/2}\right]\left|0\tilde{0}\right\rangle.$
(A8)

Thus the generalized thermal vacuum state is

$\displaystyle\left|\phi^{\prime}(\beta)\right\rangle$
$\displaystyle=R^{{\dagger}}S^{{\dagger}}\left|0(\beta)\right\rangle$
$\displaystyle=\sqrt{1-e^{-\beta\omega^{\prime}}}R^{{\dagger}}S^{{\dagger}}\exp\left[a^{\dagger}\tilde{a}^{\dagger}e^{-\beta\omega^{\prime}/2}\right]\left|0\tilde{0}\right\rangle.$
(A9)

Using the transformation in (A3), (A4) and noticing Eq.(A1) as well as
$\frac{\nu}{\mu}=\sqrt{\frac{\omega-\omega^{\prime}}{\omega+\omega^{\prime}}}$,
we can finally put Eq.(A9) into the following form

$\displaystyle\left|\phi^{\prime}(\beta)\right\rangle$
$\displaystyle=\sqrt{\left(1-e^{-\beta\omega^{\prime}}\right)/\mu}\exp\left[\frac{1}{\mu}e^{-\left(\beta\omega^{\prime}+i\phi\right)/2}a^{{\dagger}}\tilde{a}^{\dagger}\right.$
$\displaystyle\left.-\frac{\nu e^{-i\phi}}{2\mu}a^{{\dagger}2}+\frac{\nu
e^{-\beta\omega^{\prime}}}{2\mu}\tilde{a}^{\dagger
2}\right]\left|0\tilde{0}\right\rangle.$ (A10)

Comparing Eq.(A10) with Eq.(27), we see that Eq.(27) is simpler in form than
that in Eq.(A10).

ACKNOWLEDGEMENT: Worked supported by the National Natural Science Foundation
of China under Grant 10775097 and 10874174. E-mail<EMAIL_ADDRESS>
## References

  * (1) Y. Takahashi, H. Umezawa, Collective Phenomena 2 (1975) 55\.
  * (2) Memorial Issue for H. Umezawa, Int. J. Mod. Phys. B 10 (1996) 1563.
  * (3) Hong-yi Fan, Hai-liang Lu and Yue Fan, Ann. Phys. 321 (2006) 480.
  * (4) A. Wünsche, J. Opt. B: Quant. Semiclass. Opt. 1 (1999) R11.
  * (5) Hong-yi Fan, J. Opt. B: Quant. Semiclass. Opt. 5 (2003) R147.
  * (6) W. H. Louisell, Quantum Statistical Properties of Radiation (New York: John Wiley 1973).
  * (7) R. J. Glauber, Phys. Rev. 130 (1963) 2529; R. J. Glauber, Phys. Rev. 131 (1963) 2766.
  * (8) J. R. Klauder and B. S. Skargerstam, Coherent States (World Scientific, Singapore 1985).
  * (9) S. M. Barnett and P. M. Radmore, Methods in Theoretical Quantum Ooptics (Clarendon press, Oxford 1997).
  * (10) Li-yun Hu and Hong-yi Fan, Commun. Theor. Phys. (Beijing, China) 51 (2009) 321.
  * (11) R. R. Puri, Mathematical Methods of Quantum Optics (Springer-Verlag, Berlin, 2001), Appendix A.
  * (12) R. K. Pathria, Statistical Mechannics, 2nd ed. Elsevier (World Scientific, Singapore 2001) Pte Ltd.
  * (13) E. Wigner, Phys. Rev. 40 (1932) 749.
  * (14) Wolfgang P. Schleich, Quantum Optics in Phase Space, (Wiley-VCH, Birlin 2001).
  * (15) A. I. Lvovsky and M. G. Raymer, Rev. Mod. Phys. 81 (2009) 299 and references therein.
  * (16) M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge: Cambidge University Press, 1997).
  * (17) Li-yun Hu and Hong-yi Fan, J. Opt. Soc. Am. B 25 (2008) 1955.
  * (18) Hong-yi Fan and Hai-ling Chen, Chin. Phys. Lett. 18 (2001) 850.