As described in the introduction, we infer the polynomial decay under $H$ from
exponential decay of some induced return map $H_{A}$, where returns to $A$
experience ‘strong’ hyperbolic behaviour. The natural choice for $A$,
following the work of section sec:Hmixing, is the set $\sigma$. We begin by
proving the Bernoulli property for $H_{\sigma}$.

Bernoulli property

The return map $H_{\sigma}$ is Bernoulli with respect to the probability
measure $\mu_{\sigma}=\mu(\sigma)^{-1}\mu$. We will show the conditions
(KS1-3) and (MR); the result then follows from Theorem LABEL:thm:katok-
strelcyn. The return map $H_{\sigma}$ satisfies (KS1-3). Starting with (KS1)
we follow a similar approach to [springham_polynomial_2014], their Lemma 4.1.
We show that there exists $a,C_{1}>0$ s.t. $\forall\,\epsilon>0$,
$\mu_{\sigma}(B_{\varepsilon}(S))\leq C_{1}\varepsilon^{a}$ for
$S=\mathcal{S}\cap\sigma_{1a}$; the argument for the rest of $\mathcal{S}$ is
similar and the result then follows by taking a larger $C_{1}$. Recall the
line segments $\mathcal{L}_{k}$ from ($\dagger$) which for $k\geq 3$ terminate
on the points $\left(0,(k+1)/(4k+2)\right)$ and
$\left(1/(4k-2),(k-1)/(4k-2)\right)$ on the line $\mathcal{L}:y=1/4-x/2$. Let
$P(\varepsilon)$ denote the parallelogram in $\sigma_{1a}$ of width
$2\sqrt{\varepsilon}$, height $\sqrt{\varepsilon}$, with sides aligned with
$x=0$ and $\mathcal{L}$ (see Figure LABEL:fig:sigma1a). For small
$\varepsilon$, $P(\varepsilon)$ contains all line segments $\mathcal{L}_{k}$
where $2\sqrt{\varepsilon}\geq 1/(4k-2)$, i.e. $k\geq k_{0}=\lceil
1/(8\sqrt{\varepsilon})+1/2\rceil$, with

$\mu(B_{\varepsilon}(P(\varepsilon)))=(2\sqrt{\varepsilon}+2\varepsilon)(\sqrt{\varepsilon}+2\varepsilon)=2\varepsilon+6\varepsilon^{3/2}+4\varepsilon^{2}<12\varepsilon.$

The ball $B_{\varepsilon}(P(\varepsilon))$ then covers all of
$B_{\varepsilon}(S)$ except the collection $\mathcal{L}_{k}$, $4\leq k\leq
k_{0}-1$ and the seven line segments $L_{j}$ which terminate on $y=2x$. The
measure of the ball around these latter line segments satisfies

$\mu(B_{\varepsilon}(\cup_{j}L_{j}))\leq
14\,\varepsilon\left(\max_{j}|L_{j}|+2\varepsilon\right)<c_{1}\varepsilon$

for some finite $c_{1}$, so it remains to estimate
$\sum_{k=4}^{k_{0}-1}\mu(B_{\varepsilon}(\mathcal{L}_{k}))$. We can calculate

$|\mathcal{L}_{k}|=\sqrt{\left(\frac{1}{4k-2}\right)^{2}+\left(\frac{k+1}{4k+2}-\frac{k-1}{4k-2}\right)^{2}}=\sqrt{\frac{8k^{2}+4k+1}{4(4k^{2}-1)^{2}}}<\frac{1}{k}$
(1)

so that

$\sum_{k=4}^{k_{0}-1}\mu(B_{\varepsilon}(\mathcal{L}_{k}))<2\varepsilon\sum_{k=4}^{k_{0}-1}\frac{1}{k}+2\varepsilon<4\varepsilon^{2}k_{0}+2\varepsilon\log
k_{0}<c_{2}\varepsilon^{a}$

for some $0<a<1$, $c_{2}>0$ since $k_{0}<\varepsilon^{-1/2}$ and there exists
finite $c$ such that $c\,\varepsilon^{a}>\varepsilon\log\frac{1}{\varepsilon}$
for any $0<a<1$. Since $H_{\sigma}$ is piecewise linear, condition (KS2)
follows trivially and we move onto (KS3). Existence of Lyapunov exponents
almost everywhere follows from Oseledets’ theorem
[oseledets_multiplicative_1968] provided that
$\max\\{\log\|DH_{\sigma}\|,0\\}$ is integrable. This follows from the fact
that if $z\in\sigma$ has return time $R(z;H,\sigma)=k$, then the Jacobian of
$H_{\sigma}$ at $z$ satisfies $\|DH_{\sigma}\|\leq c_{1}k$ for some finite
$c_{1}>0$, and that the measure of the sets
$\\{z\in\sigma\,|\,R(z;H,\sigma)=k\\}$ are of order $k^{-3}$. That these
Lyapunov exponents are non-zero follows from Lemmas LABEL:lemma:itineraries,
LABEL:lemma:cone and an argument similar to that given for $H$ in section
LABEL:sec:OTMhyp. The return map $H_{\sigma}$ satisfies (MR). For a.e
$z\in\sigma$, local manifolds $\gamma_{u}(z)$, $\gamma_{s}(z)$ under
$H_{\sigma}$ align with those of $H$. Note that $H_{\sigma}$ does not
immediately inherit (MR) from $H$ as while successive images of local
manifolds under $H$ contain $h$-segments and $v$-segments, these segments may
not lie in the successive images under $H_{\sigma}$.

[scale=0.8] tomatoRGB255, 99, 71 tealRGB95, 158, 160 (5,0) – (5,10); [very
thick, tomato] (2,2.8) – (3,2.2); [very thick, teal] (10-2,2.8+5) –
(10-3,2.2+5); [fill=gray!40] (0,8) – (10/8,10) – (2.5,10) – (0,5) – (0,8);
[fill=gray!40] (0,9) – (0,10) – (10/16,10) – (0,9); [fill=gray!40] (0,0) –
(5,10) – (5,7) – (10/16,0) – (0,0); [fill=gray!40] (10/8,0) – (2.5,0) – (5,5)
– (5,6) – (10/8,0) ; [fill=gray!40] (5,2) – (50/8,0) – (7.5,0) – (5,5) –
(5,2); [fill=gray!40] (5,1) – (5,0) – (90/16,0) – (5,1); [fill=gray!40] (5,10)
– (10,0) – (10,3) – (90/16,10) – (5,10); [fill=gray!40] (50/8,10) – (7.5,10) –
(10,5) – (10,4) – (50/8,10); [fill=gray!40] (0,150/16) – (0,10) – (1,10) –
(0,150/16); [fill=gray!40] (0,70/8) – (2,10) – (5,10) – (0,7.5);
[fill=gray!40] (0,5) – (10,10) – (10,150/16) – (3,5) – (0,5); [fill=gray!40]
(10,70/8) – (10,7.5) – (5,5) – (4,5) – (10,70/8) ; [fill=gray!40]
(10,150/16-5) – (10,10-5) – (10-1,10-5) – (10,150/16-5); [fill=gray!40]
(10,70/8-5) – (10-2,10-5) – (10-5,10-5) – (10,7.5-5); [fill=gray!40] (10,5-5)
– (0,10-5) – (0,150/16-5) – (10-3,5-5) – (10,5-5); [fill=gray!40] (0,70/8-5) –
(0,7.5-5) – (10-5,5-5) – (10-4,5-5) – (0,70/8-5) ; (0,5) – (10,5); (0,0)
rectangle (10,10);

Figure : The set $\sigma^{\prime}$ superimposed on $\sigma$, their
intersection left white. A $\mathfrak{h}$-segment (red) and a
$\mathfrak{h}^{\prime}$-segment (blue) are plotted in
$\mathcal{R}=\sigma_{3}\cap\sigma_{2}^{\prime}$ and
$\mathcal{R}^{\prime}=\sigma_{2}\cap\sigma_{3}^{\prime}$ respectively.

Let $\mathcal{R}$ denote the quadrilateral $\sigma_{3}\cap\sigma_{2}^{\prime}$
and $\mathcal{R}^{\prime}=\sigma_{2}\cap\sigma_{3}^{\prime}$. Define a
$\mathfrak{h}$-segment as a line segment spanning $\mathcal{R}$ with endpoints
on $\partial\sigma_{3}$. Similarly define a $\mathfrak{h}^{\prime}$-segment as
a line segment spanning $\mathcal{R}^{\prime}$ with endpoints on
$\partial\sigma_{2}$. Examples are plotted in Figure 1. We will show that
there exists $M,N$ such that for all $m\geq M$, $n\geq N$,
$H_{\sigma}^{m}(\gamma_{u}(z))$ intersects
$H_{\sigma}^{-n}(\gamma_{s}(z^{\prime}))$ in either $\mathcal{R}$ or
$\mathcal{R}^{\prime}$. By the remark after Theorem LABEL:thm:mixingOTM we can
find some $n_{2}$ such that $H^{-n_{2}}(\gamma_{s}(z^{\prime}))$ contains a
$h$-segment in $\sigma_{2}^{\prime}$, which in turn contains a
$\mathfrak{h}$-segment in $\mathcal{R}$. As a line segment in $\sigma$, this
$\mathfrak{h}$-segment lies in $H_{\sigma}^{-n_{1}}(\gamma_{s}(z^{\prime}))$
for some $n_{1}\leq n_{2}$. Note that we have a hyperbolic period 2 orbit
$(1/4,1/4)\leftrightarrow(3/4,3/4)$ under $H_{\sigma}$, alternating between
$\mathcal{R}$ and $\mathcal{R}^{\prime}$. Any $\mathfrak{h}$-segment $\Lambda$
contains a point $\zeta$ on the unstable manifold through $(1/4,1/4)$ so
$H_{\sigma}^{-2}(\Lambda)$ contains a point $\zeta^{\prime}$ on the manifold
closer to $(1/4,1/4)$ and extends beyond the boundaries $\partial\sigma_{3}$
by the expansion of $H_{\sigma}^{-2}$. Hence $H_{\sigma}^{-2}(\Lambda)$
contains a $\mathfrak{h}$-segment and by induction so does
$H_{\sigma}^{-2k}(\Lambda)$ for all $k\geq 1$. The odd iterates
$H^{-2k+1}(\Lambda)$ similarly span $\mathcal{R}^{\prime}$ so that

  1. [label=($\vee$)]

  2. 1.

Given arbitrary $z^{\prime}\in\sigma$, there exists $N$ such that for all
$n\geq N$ the image $H_{\sigma}^{-n}(\gamma_{s}(z^{\prime}))$ contains a
$\mathfrak{h}$-segment or a $\mathfrak{h}^{\prime}$-segment.

Define a $v^{\prime}$-segment as a line segment vertically spanning
$\sigma_{2}\cap S_{4}$. Condition (MR) now follows from establishing

  1. [label=($\wedge$)]

  2. 1.

Given arbitrary $z\in\sigma$, there exists $M$ such that for all $m\geq M$ the
image $H_{\sigma}^{m}(\gamma_{s}(z^{\prime}))$ contains both a $v$-segment
_and_ a $v^{\prime}$-segment.

Recall the quadrilaterals $\varsigma_{j}\subset\sigma_{j}$ of points with
return time 1 (see Figure LABEL:fig:returnTime1). It follows from the
definitions of the $\sigma_{j}$ that $H(\varsigma_{2})\subset\sigma_{3}$ and
$H(\varsigma_{3})\subset\sigma_{2}$. The edges of $\varsigma_{3}$ on the
$A_{1},A_{2}$ boundary map into the lines $x=1/2$, $x=1$, in particular onto
the red dashed lines on the boundary of $S_{2}$ in Figure LABEL:fig:HSigma so
that the image of any line segment in $\varsigma_{3}$ which joins these edges
contains a $v^{\prime}$-segment. A analogous result holds for lines segments
traversing $\varsigma_{2}$ and since this behaviour occurs within the return
set $\sigma$ we have that $v$-segments map into $v^{\prime}$-segments under
$H_{\sigma}$ and vice versa. It follows that $v$(’)-segments map into
$v$(’)-segments under $H_{\sigma}^{2}$. It suffices to break into the odd
iterates to satisfy the ‘and’ condition ($\wedge$). By following the steps
(M1-3) in the proof of Theorem LABEL:thm:mixingOTM we can find $m_{2}$ such
that $H^{m_{2}}(\gamma_{u}(z))$ contains a $v$-segment
$\Gamma\subset\sigma_{3}\cap S_{1}$ with $H^{-1}(\Gamma)\subset A_{3}$,
$H^{-2}(\Gamma)\subset A_{2}$, $H^{-3}(\Gamma)\subset A_{1}$,
$H^{-4}(\Gamma)\subset A_{4}$. In particular $H^{-1}(\Gamma)$ lies in
$H(A_{2}\cap H(A_{1}))$, i.e. outside of $\sigma_{2}$, so that $\Gamma$ lies
in $\sigma_{3}\setminus(H(\sigma_{2})\cap\sigma_{3})$. The set
$H(\sigma_{2})\cap S_{1}$ (shown in blue in Figure 1) is the quadrilateral
with corners $(7/68,0)$, $(3/34,0)$, $(27/68,0)$, $(7/17,0)$, which splits
$\sigma_{3}\cap S_{1}$ into left and right parts. We assume first that
$\Gamma$ lies in the right part, intersecting the line $y=1/2$ at some point
$(x_{1},1/2)$ with $x_{1}\geq 7/17$ and the $A_{1},A_{2}$ boundary $y=1/2-x/2$
at some point $(1-2y_{1},y_{1})$. These intersections define a line segment
$\Gamma_{1}\subset\Gamma$, which lies in $A_{1}$ shown in Figure 1. Applying
$F$ maps $(x_{1},1/2)$ to itself (wrapping horizontally around the torus) and
maps $(1-2y_{1},y_{1})$ to $(0,y_{1})$. Applying $G$ then leaves $(0,y_{1})$
invariant and wraps $(x_{1},1/2)$ vertically around the torus to
$(x_{1},1/2+2x_{1})\text{mod}1\equiv(x_{1},-1/2+2x_{1})$. Since $x_{1}\geq
7/17$ we have that $-1/2+2x_{1}\geq 11/34>10/24\geq 1/2-x_{1}/2$ so that
$(x_{1},-1/2+2x_{1})$ lies above the line $y=1/2-x/2$. We restrict again to
$A_{1}$, giving $\Gamma_{2}\subset H(\Gamma_{1})$ with endpoints on
$(x_{1},-1/2+2x_{1})$ and some point $(1-2y_{2},y_{2})$ on $y=1/2-x/2$. This
line meets $y=-1/2+2x$ at $(2/5,3/10)$ so that $y_{2}\leq 3/10$ (see Figure
1). Now $F(\Gamma_{2})$ joins $(0,y_{2})$ to $(5x_{1}-2,2x_{1}-1/2)$ and so
$H(\Gamma_{2})$ joins $(0,y_{2})$ to $(5x_{1}-2,12x_{1}-9/2)$. The set
$A_{3,2}^{1}$ is bounded by the parallel lines $y=7/16-5x/8$ and $y=3/8-5x/8$.
Since $y_{2}<3/8$ and $12x_{1}-9/2\geq 15/34>109/272\geq 7/16-5(5x_{1}-2)/8$
we have that $H(\Gamma_{2})$ contains a segment $\Gamma_{3}$ which traverses
$A_{3,2}^{1}$. The image $H(\Gamma_{3})$ then traverses $Q_{3}$ so that
$H^{2}(\Gamma_{3})\subset H^{4}(\Gamma)$ contains a $v$-segment in
$\sigma_{3}$. Critically we have that $\Gamma_{2},\Gamma_{3}\subset
H(A_{1})\subset\sigma$ but $\Gamma_{4}\subset H(A_{2}\cap H(A_{1}))$ which is
not in $\sigma$. Hence $H_{\sigma}^{3}(\Gamma)$ contains a $v$-segment and we
can apply the $H_{\sigma}^{2}$ result above to show that
$H_{\sigma}^{k}(\Gamma)$ contains a $v$-segment for all $k\geq 2$. Similar
analysis can be applied to $\Gamma$ in the left portion of
$\sigma_{3}\setminus(H(\sigma_{2})\cap\sigma_{3})$. It follows that
$H_{\sigma}^{k}(\Gamma)$ contains a $v^{\prime}$-segment for all $k\geq 3$,
establishing ($\wedge$) with $M=m_{1}+3$.

tomatoRGB255, 99, 71 tealRGB95, 158, 160 at (-3.5,0) [scale=1.3] at
(-0.15,350/104) $y_{1}$; (0,0) rectangle (5,5); [fill=white] (0,0) – (2.5,5) –
(5,5) – (2.5,0) – (0,0); (0,5) – (10,0); (0,2.5) – (5,0); [fill=gray!80] (0,0)
– (5,10) – (5,7) – (10/16,0) – (0,0); [fill=gray!80] (10/8,0) – (2.5,0) –
(5,5) – (5,6) – (10/8,0) ; [teal] (70/17,5) – (70/68,0) – (30/34,0) –
(270/68,5) – (70/17,5); [very thick] (170/52,350/104) – (4.25,5); [dashed,
thick] (0,350/104) – (4.25,5); [very thick] (0,350/104) –
(0.8724691416999097,5); [very thick] (2.782261032261033,0) – (4.25,3.5); (0,5)
– (5,2.5) – (5,5) – (0,5); [very thick, tomato] (2.782261032261033,0) –
(4.25,3.5); at (4.15,4.3) $\Gamma_{1}$; at (1.3,3.55) $F(\Gamma_{1})$; at
(4.35,3.2) tomato$\Gamma_{2}$black; (0,0) rectangle (5,5); ; at (4,0)
[scale=1.3] at (-0.15,2.98323) tomato$y_{2}$black; (0,0) rectangle (5,5);
[fill=gray!40] (0,150/16-5) – (7,0) – (6,0) – (0,70/8-5) –(0,150/16-5);
[fill=white] (0,0) – (2.5,5) – (5,5) – (2.5,0) – (0,0); (0,5) – (10,0);
(0,2.5) – (5,0); [fill=gray!80] (0,0) – (5,10) – (5,7) – (10/16,0) – (0,0);
[fill=gray!80] (10/8,0) – (2.5,0) – (5,5) – (5,6) – (10/8,0) ; (0,5) – (5,2.5)
– (5,5) – (0,5); [very thick, tomato] (2.782261032261033,0) – (4.25,3.5);
[dashed, thick, tomato] (0,2.98323) – (1.25,3.5); [very thick, tomato]
(0,2.98323) – (1.25,6); [very thick] (0.25249,3.59197) – (0.45807,4.08851); at
(0.55,2.95) tomato$F(\Gamma_{2})$black; at (0.2,4) $\Gamma_{3}$; at (0.65,3.6)
$A_{3,2}^{1}$; at (4.35,3.2) tomato$\Gamma_{2}$black; (0,0) rectangle (5,5); ;
at (-1.02,3.45) $x_{1}$;

Figure : Left: The upper part $\Gamma_{1}\subset A_{1}$ of a $v$-segment in
$\sigma_{3}\setminus(H(\sigma_{2})\cap\sigma_{3})$ and its images
$F(\Gamma_{1})$ (dashed), $H(\Gamma_{1})\cap S_{1}$. Right: A right part
$\Gamma_{2}\subset A_{1}$ of $H(\Gamma_{1})$ and its images $F(\Gamma_{2})$,
$H(\Gamma_{2})\cap S_{1}$. This image necessarily contains a line segment
$\Gamma_{3}$ traversing $A_{3,2}^{1}$.

##  Invariant cones

We now derive specific unstable and stable cone fields for the return map
$H_{\sigma}$, wide enough to ensure invariance (H1.1) yet fine enough to
produce tight bounds on expansion factors, vital for verifying (H5). Define
the cones $\mathcal{C}_{1},\dots,\mathcal{C}_{4}$ by

  1. [label=($\mathcal{C}_{0}$)]

  2. 1.

$3|v_{1}|\geq|v_{2}|\geq 7|v_{1}|/3$, $v_{1}v_{2}>0$,

  3. 2.

$5|v_{1}|/3\geq|v_{2}|\geq\varphi|v_{1}|$, $v_{1}v_{2}<0$,

  4. 3.

$5|v_{1}|/3\geq|v_{2}|\geq\varphi|v_{1}|$, $v_{1}v_{2}>0$,

  5. 4.

$3|v_{1}|\geq|v_{2}|\geq 7|v_{1}|/3$, $v_{1}v_{2}<0$,

and the following stable cones

  1. [label=($\mathcal{C}_{0}^{s}$)]

  2. 1.

$|v_{1}|\geq|v_{2}|$,

  3. 2.

$9/10\geq v_{2}/v_{1}\geq-8/10$,

  4. 3.

$8/10\geq v_{2}/v_{1}\geq-9/10$,

  5. 4.

$|v_{1}|\geq|v_{2}|$.

[scale=1.5] (0,0) circle (3cm); [gray] (0,0) – (10,210/13) – (10,50/3) –
(0,0); [gray] (0,0) – (-10,-210/13) – (-10,-50/3) – (0,0); [gray] (0,0) –
(-10,210/13) – (-10,50/3) – (0,0); [gray] (0,0) – (10,-210/13) – (10,-50/3) –
(0,0); [gray] (0,0) – (10,70/3) – (10,30) – (0,0); [gray] (0,0) – (10,-70/3) –
(10,-30) – (0,0); [gray] (0,0) – (-10,70/3) – (-10,30) – (0,0); [gray] (0,0) –
(-10,-70/3) – (-10,-30) – (0,0); [gray,opacity=0.5] (0,0) – (-10,8) – (-10,-9)
– (0,0); [gray,opacity=0.5] (0,0) – (10,-8) – (10,9) – (0,0); [red]
(-10,-80/5) – (10,80/5); [red] (10,-80/5) – (-10,80/5); [red] (-10,-20) –
(10,20); [red] (10,-20) – (-10,20); [red] (0,-10) – (0,10); (0,0) circle
(2.5cm); [gray,opacity=0.5] (0,0) – (-10,9) – (-10,-8) – (0,0);
[gray,opacity=0.5] (0,0) – (10,-9) – (10,8) – (0,0); (0,0) circle (2cm);
[gray,opacity=0.5] (0,0) – (-10,10) – (-10,-10) – (0,0); [gray,opacity=0.5]
(0,0) – (10,-10) – (10,10) – (0,0); at (-1.38,1.4+0.1) $\mathcal{C}_{1}^{s}$;
at (-1.76,1.67+0.1) $\mathcal{C}_{3}^{s}$; at (-2.28,1.85+0.1)
$\mathcal{C}_{2}^{s}$; at (-1.5,2.68) $\mathcal{C}_{2}$; at (1.55,2.68)
$\mathcal{C}_{3}$; at (-1.1,2.93) $\mathcal{C}_{4}$; at (1.1,2.92)
$\mathcal{C}_{1}$;

Figure : Unstable and stable cone fields $\mathcal{C}_{j}$,
$\mathcal{C}_{j}^{s}$ over the subsets $\sigma_{j}\subset\sigma$ for the
return map $H_{\sigma}$. Also shown in red are the gradients of the line
segments which make up the boundary $\partial\sigma$ which lie outside of all
cone fields.

In the notation of section LABEL:sec:outline, for general $z\in\sigma$ we take
$C_{z}^{u}=\mathcal{C}_{j}$ and $C_{z}^{s}=\mathcal{C}_{j}^{s}$ for
$z\in\sigma_{j}$. These cone fields are plotted in Figure 1.

$DH_{\sigma}$ | $\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{3}$ | $\sigma_{4}$
---|---|---|---|---
$\sigma_{1}$ | $M_{1}$ | - | $M_{3}M_{2}^{k}$ |  $M_{4}$
$M_{4}M_{2}^{k}$ |  |  |  |
$M_{4}M_{3}^{k}$ |  |  |  |
$\sigma_{2}$ | - | - |  $M_{3}$ |
$M_{3}M_{2}$ |  |  |  |
$M_{3}M_{2}^{2}$ |  $M_{4}$ |  |  |
$M_{4}M_{2}$ |  |  |  |
$M_{4}M_{2}^{2}$ |  |  |  |
$\sigma_{3}$ |  $M_{1}$ |  |  |
$M_{1}M_{3}$ |  |  |  |
$M_{1}M_{3}^{2}$ |  $M_{2}$ |  |  |
$M_{2}M_{3}$ |  |  |  |
$M_{2}M_{3}^{2}$ | - | - |  |
$\sigma_{4}$ |  $M_{1}$ |  |  |
$M_{1}M_{2}^{k}$ |  |  |  |
$M_{1}M_{3}^{k}$ | $M_{2}M_{3}^{k}$ | - | $M_{4}$ |
Table : Possible values of the Jacobian $DH_{\sigma}$ at $z$ if
$z\in\sigma_{i}$ (rows) and $z^{\prime}=H_{\sigma}(z)\in\sigma_{j}$ (columns).
Exponent $k$ takes values in $\mathbb{N}$, dashes are shown if no transition
is possible, e.g $H_{\sigma}(\sigma_{1})\cap\sigma_{2}=\varnothing$.

The above cones satisfy $DH_{\sigma}\,C_{z}^{u}\subset C_{z^{\prime}}^{u}$ and
$DH_{\sigma}\,C_{z}^{s}\supset C_{z^{\prime}}^{s}$ for all $z\in\sigma$ where
$DH_{\sigma}$ exists, $z^{\prime}=H_{\sigma}(z)$. We begin with the unstable
cones. Table 1 shows the possible values of $DH_{\sigma}$ at $z$ if
$z\in\sigma_{i}$ and $z^{\prime}\in\sigma_{j}$. The calculations for
$z^{\prime}\in\sigma_{1},\sigma_{4}$ are similar to those made in the proof of
Lemma LABEL:lemma:mappedOTMalignment, noting that each $\mathcal{C}_{j}$ is
contained within $\mathcal{C}$ and, for example,
$M_{1}M_{j}^{k}\mathcal{C}\subset\mathcal{C}_{1}$ for $j=2,3$, and $k\geq 0$.
For $z^{\prime}\in\sigma_{2}$ we verify that
$M_{2}M_{3}^{k}(-1,3)^{T}=(-1)^{k}(-24k+5,-40k-7)^{T}\in\mathcal{C}_{2}$ and
$M_{2}M_{3}^{k}(-3,7)^{T}=(-1)^{k}(60k+11,-100k-15)^{T}\in\mathcal{C}_{2}$ for
all $k\geq 1$ so that $DH_{\sigma}\,C_{z}^{u}\subset C_{z^{\prime}}^{u}$ for
$z\in\sigma_{4}$. For $z\in\sigma_{3}$ we have
$M_{2}(3,5)^{T}=(13,-21)^{T}\in\mathcal{C}_{2}$ and
$M_{2}(13,21)^{T}=(55,-89)^{T}\in\mathcal{C}_{2}$, ensuring invariance in this
particular case also, despite $M_{2}$ being non-hyperbolic. Entirely symmetric
calculations can be made for $z^{\prime}\in\sigma_{3}$, verifying the result
for all unstable cones. For the stable cones, we remark that taking
$C_{z}^{s}=\mathcal{C}_{1}^{s}$ for all $z\in\sigma$ would satisfy
$DH_{\sigma}\,C_{z}^{s}\supset C_{z^{\prime}}^{s}$ but since
$M_{2}^{-1}(1,-1)^{T}=(-1,1)^{T}$ we would be unable to derive sufficient
uniform bounds on expansion factors (H1.2). The matrix $M_{3}^{-1}$ exhibits a
similar problem so we must slim down the cones $C_{z}^{s}$ when
$DH_{\sigma}^{-1}\in\\{M_{2}^{-1},M_{3}^{-1}\\}$ which, observing Table 1, is
for $z\in\sigma_{2},\sigma_{3}$. To remedy this, for such $z$ we slim down the
cones $C_{z}^{s}$ to $\mathcal{C}_{2}^{s},\mathcal{C}_{3}^{s}$ above. As these
cones lie in the wider invariant cone $|v_{1}|\geq|v_{2}|$, the lemma follows
from checking that $DH_{\sigma}\,C_{z}^{s}\supset C_{z^{\prime}}^{s}$ for
$z\in\sigma_{2},\sigma_{3}$. This can be verified via direct calculations.

##  Structure of the singularity set

Using the notation of (H2) in section LABEL:sec:outline, let
$\mathcal{S}_{0}=\cup_{j}\partial\sigma_{j}$, the union of $\partial\sigma$
and the red dashed lines in Figure LABEL:fig:HSigma. The set
$M=\Omega\setminus\mathcal{S}_{0}$ is clearly dense in $\Omega$ and
$H_{\sigma}$ is a $C^{2}$ diffeomorphism from $M\setminus\mathcal{S}_{1}$ onto
$M\setminus\mathcal{S}_{-1}$, being linear on each component. The set
$\mathcal{S}_{0}\cup\mathcal{S}_{1}$ is the countable union of bounded line
segments with the endpoints of each segment terminating on another segment,
giving (H2.2). The gradients of the segments in $\mathcal{S}_{0}$ take values
in $\\{\pm 8/5,\pm 2,\infty\\}$ which avoid unstable and stable cones
$C_{z}^{u}$, $C_{z}^{s}$ (see Figure 1). The gradients of singularity curves
in $\sigma_{1}$ and $\sigma_{4}$ are bounded between -1 and 1 (approaching
these limits as we approach the accumulation points) so lie in
$\mathcal{C}_{1}^{s},\mathcal{C}_{4}^{s}$. The gradients of singularity curves
in $\sigma_{2}$ and $\sigma_{3}$ are bounded between -11/14 and 11/14 so lie
in $\mathcal{C}_{2}^{s},\mathcal{C}_{3}^{s}$ since $11/14<8/10$. Similar
calculations show that the gradients of segments in $\mathcal{S}_{-1}$ lie in
unstable cones. We conclude this section with showing (H2.4). Condition
(LABEL:eq:DxbCondition) can only fail when $\|DH_{\sigma}\|$ becomes
unbounded, i.e. at points $z$ approaching the accumulation points. We consider
the case with $z\in A_{4,2}^{k}$ near $(0,1/4)$, the other cases are similar.
Recall Figure LABEL:fig:sigma1a and the lines $\mathcal{L}_{k}$ from
($\dagger$). We note that $d(z,\mathcal{S}_{1})$ is bounded above by the
length of the segment joining $z=(x,y)$ to $(x,y_{k}(x))$ on
$\mathcal{L}_{k}$, which in turn is bounded above by the height of the segment
joining $(x,y_{k}(x))$ to $(x,y_{k-1}(x))$ on $\mathcal{L}_{k-1}$. This height
is

$y_{k-1}(x)-y_{k}(x)\leq
y_{k-1}\left(\frac{1}{4k-2}\right)-y_{k}\left(\frac{1}{4k-2}\right)=\frac{1}{2(2k-1)^{2}}\leq
c_{1}/k^{2}$

for some constant $c_{1}>0$. The operator norm of $DH_{\sigma}$ over
$A_{4,2}^{k}$ satisfies $\|DH_{\sigma}\|\leq c_{2}k$ for some $c_{2}>0$ so
that (LABEL:eq:DxbCondition) holds for some $c>0$ whenever we choose $b>1/2$.

##  One-step expansion

We will verify (LABEL:eq:oneStep) for the map $f=H_{\sigma}^{2}$, $q=1/2$. We
begin with a basic statement on expansion over unstable curves. Let $M$ be the
constant Jacobian of $f$ over $W_{i}$, $V_{i}=f(W_{i})$. Then

$\lambda^{-}:=\inf_{v\in\mathcal{C}}\frac{\|Mv\|}{\|v\|}\leq\frac{|V_{i}|}{|W_{i}|}\leq\sup_{v\in\mathcal{C}}\frac{\|Mv\|}{\|v\|}=:\lambda^{+}.$

Given any $\varepsilon>0$, consider a piecewise linear approximation
$\hat{W_{i}}$ to $W_{i}$ such that

$\left|\frac{|V_{i}|}{|W_{i}|}-\frac{|\hat{V_{i}}|}{|\hat{W_{i}}|}\right|<\varepsilon$

where $\hat{V_{i}}=f(\hat{W_{i}})$ gives a piecewise approximation for
$V_{i}$. Each of the piecewise components will be line segments aligned with
vectors in $\mathcal{C}$ so that their expansion factors will be bounded by
$\lambda^{\pm}$, giving the result. We also derive basic inequalities on the
length of a given $W_{i}$. Let $\mathcal{L}_{0},\mathcal{L}_{1}$ be the
singularity curves on which $W_{i}$ terminates, write these intersections as
$(x_{0},y_{0})$ and $(x_{1},y_{1})$. Then

$\sqrt{(x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2}}\leq|W_{i}|\leq|y_{1}-y_{0}|\sqrt{1+\frac{1}{g^{2}}}$

where $g=\inf|v_{2}/v_{1}|$ over $(v_{1},v_{2})^{T}\in\mathcal{C}$. Noting
that the lower bound is trivial, we focus on the upper bound. Since $g>0$ for
all unstable cones $\mathcal{C}$, the projection of $W_{i}$ to the $y$-axis is
injective. Without loss of generality suppose $y_{1}>y_{0}$, then we can
parameterise $W_{i}$ as a curve $(x(y),y)$ for $y_{0}\leq y\leq y_{1}$. Now

$\split|W_{i}|&=\int_{y_{0}}^{y_{1}}\sqrt{\left(\frac{\mathrm{d}x}{\mathrm{d}y}\right)^{2}+\left(\frac{\mathrm{d}y}{\mathrm{d}y}\right)^{2}}\d{y}\\\
\leq(y_{1}-y_{0})\sup_{y_{0}\leq y\leq
y_{1}}\sqrt{1+\left(\frac{\mathrm{d}x}{\mathrm{d}y}\right)^{2}}\\\
\leq(y_{1}-y_{0})\sqrt{1+\frac{1}{g^{2}}}$

as tangent vectors $(x^{\prime}(y),1)^{T}$ to $W_{i}$ lie in $\mathcal{C}$.
Let $P_{1}=\\{(0,1/4),(1/2,1/4),(1/2,3/4),(1,3/4)\\}$ denote the accumulation
points similar to that of $\sigma_{1a}$,
$P_{2}=\\{(1/4,1/2),(1/4,1),(3/4,0),(3/4,1/2)\\}$ the accumulation points
similar to that of $\sigma_{1b}$. Let $\varepsilon$ be small. Given a set $P$,
let $B_{\varepsilon}(P)$ denote the union of the balls
$B_{\varepsilon}(p)\cap\sigma$, centred at $p\in P$ of radius $\varepsilon$.
The following describes the images of balls about $P_{1}\cup P_{2}$ under
$H_{\sigma}$. Given small $\varepsilon>0$, there exists some
$\varepsilon^{\prime}>0$ such that $H_{\sigma}(B_{\varepsilon}(P_{1}\cup
P_{2}))$ covers $B_{\varepsilon^{\prime}}(P_{1}\cup P_{2})$. We describe the
covering of $B_{\varepsilon^{\prime}}((1/2,3/4))$, analysis for the other
points in $P_{1}\cup P_{2}$ is analogous. For any $\varepsilon>0$,
$B_{\varepsilon}(P_{1}\cup P_{2})$ contains the sets $A_{4,3}^{k}$ for all
$k\geq k_{0}$ where $k_{0}\in\mathbb{N}$ depends on $\varepsilon$. Each
$A_{4,3}^{k}$ consists of two quadrilaterals, one in the ball around $(1/4,1)$
and the other in the ball around $(1/4,1/2)$. Figure 1 shows this latter
quadrilateral, with corners on the points

$r_{1}=\left(\frac{k+1}{4k+2},\frac{k+1}{2k+1}\right),\quad
r_{2}=\left(\frac{k-1}{4k-2},\frac{1}{2}\right),\quad
r_{3}=\left(\frac{k}{4k+2},\frac{1}{2}\right),\quad
r_{4}=\left(\frac{k+2}{4k+6},\frac{k+2}{2k+3}\right).$

Since $DH_{\sigma}$ is constant on $A_{4,3}^{k}$, given by the integer valued
matrix $M_{4}M_{3}^{k}=(-1)^{k}\big{(}\smallmatrix 2k+1&-2k-2\\\
-6k-26k+5\big{)}$, its image $H_{\sigma}\left(A_{4,3}^{k}\right)$ is given by
the quadrilateral with corners given by
$M_{4}M_{3}^{k}\,r_{j}^{T}\text{mod}1$. For odd $k$ we can calculate these
corners as

$r_{1}^{\prime}(k)=\left(\frac{1}{2}+\frac{1}{4k+2},\frac{3}{4}-\frac{5}{8k+4}\right),\quad
r_{2}^{\prime}(k)=\left(\frac{1}{2}+\frac{1}{4k-2},\frac{3}{4}-\frac{5}{8k-4}\right),$
$r_{3}^{\prime}(k)=\left(\frac{1}{2},\frac{3}{4}+\frac{1}{8k+4}\right),\quad
r_{4}^{\prime}(k)=\left(\frac{1}{2},\frac{3}{4}+\frac{1}{8k+12}\right),$

shown in Figure 1. For even $k$ the corners of $A_{4,3}^{k}$ in the ball
around $(1/4,1)$ map into the $r_{j}^{\prime}$. Writing this quadrilateral as
$Q(k)$, since $r_{2}^{\prime}(k+1)$ = $r_{1}^{\prime}(k)$ and
$r_{3}^{\prime}(k+1)=r_{4}^{\prime}(k)$ we have that $\cup_{k\geq k_{0}}Q(k)$
is the polygon with corners $r_{2}^{\prime}(k_{0})$, $r_{3}^{\prime}(k_{0})$,
and
$\lim_{k\to\infty}r_{1}^{\prime}(k)=\lim_{k\to\infty}r_{4}^{\prime}(k)=(1/2,3/4)$.
Noting that $r_{3}^{\prime}(k_{0})>3/4$ and $r_{2}^{\prime}(k_{0})$ lies on
the line $y-\frac{3}{4}=-\frac{5}{2}(x-\frac{1}{2})$, there exists
$\varepsilon^{\prime}$ such that $\cup_{k\geq k_{0}}Q(k)$ covers all points
$(x,y)\in B_{\varepsilon^{\prime}}((1/2,3/4))$ with
$y\geq\frac{3}{4}-\frac{5}{2}(x-\frac{1}{2})$. The image
$H_{\sigma}(B_{\varepsilon}((1,3/4))\cap A_{4})$ fills the remaining portion
of $B_{\varepsilon^{\prime}}((1/2,3/4))$, since $H_{\sigma}=H$ on $A_{4}$,
$H(1,3/4)=(1/2,3/4)$, $DH\,(-2,1)^{T}=(0,-1)^{T}$, and
$DH\,(0,-1)^{T}=(2,-5)^{T}$. Condition (LABEL:eq:oneStep) holds for
$H_{\sigma}^{2}$ when there exists $\varepsilon>0$ such that $W\cap
B_{\varepsilon}(P_{1}\cup P_{2})=\varnothing$. We claim that an unstable curve
$W$ of vanishing length, bounded away from the accumulation points, is split
into at most 9 components $W_{i}$ by the singularity set for $H_{\sigma}^{2}$.
The upper bound follows from analysis of the original singularity set for
$H_{\sigma}$. Let $P_{F}$ denote the set of fixed points under $H$,
$P_{F}=\\{(0,1/2),(1/2,0),(1/2,1/2),(1,1)\\}$. Observing Figure
LABEL:fig:sigmaSingSet, if $W\cap B_{\varepsilon}(P_{F})\neq\varnothing$ then
$W$ is split by $\mathcal{S}$ into at most 5 components $W_{j}$, and if $W\cap
B_{\varepsilon}(P_{F})=\varnothing$ then the upper bound is 3. We consider
these cases separately. Take, for example, $W\cap
B_{\varepsilon}((0,1/2))\neq\varnothing$. Observing Figure LABEL:fig:sigma1b,
four of the components $W_{j}$ map into $A_{4}^{\prime}$ under $H_{\sigma}$,
and their images lie in some sector $B_{\varepsilon^{\prime}}((1,1/2))\cap
A_{4}^{\prime}$. We can take $\varepsilon$ small enough that this sector lies
entirely in $A_{4}$, so that no further splitting occurs during the next
iterate of $H_{\sigma}$. The other component $W\cap A_{1}$ maps into some
sector $B_{\varepsilon^{\prime}}((0,1/2))\cap A_{1}^{\prime}$ and is split
into at most 5 components, giving at most $N=9$ components in total. The other
cases $W\cap B_{\varepsilon}(p)\neq\varnothing$, $p\in P_{F}$, are analogous.
Now suppose $W\cap B_{\varepsilon}(P_{F})=\varnothing$. $\mathcal{S}$ splits
$W$ into at most 3 components $W_{j}$ and, by Lemma  and the above, each
$H_{\sigma}(W_{j})$ is bounded away from the accumulation points $P_{1}\cup
P_{2}$ and the fixed points $P_{F}$. Hence each $H_{\sigma}(W_{j})$ is split
into at most 3 components during the next iterate of $H_{\sigma}$, again
giving at most $N=9$ components in total. The weakest expansion of
$DH_{\sigma}^{2}$ over cones $\mathcal{C}_{j}$ on $\sigma_{j}$ using the
euclidean norm is that of $M_{1}M_{4}=\big{(}\smallmatrix-3&8\\\ -821\big{)}$
on $\sigma_{1}$ (or equivalently $M_{4}M_{1}$ on $\sigma_{4}$), and is given
by

$c=\frac{\|M_{1}M_{4}(3,7)^{T}\|}{\|(3,7)^{T}\|}\sqrt{\frac{(-9+56)^{2}+(-24+147)^{2}}{3^{2}+7^{2}}}=\sqrt{\frac{8669}{29}}\approx
17.29$

so that, by Lemma 1, $|V_{i}|\geq c\,|W_{i}|$ for each component $W_{i}$. Now
for $q=1/2$ we have

$\split\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}&=\sum_{i}\sqrt{\frac{|W|}{|V_{i}|}}\frac{|W_{i}|}{|W|}\\\
=\sum_{i}\sqrt{\frac{|W_{i}|}{|V_{i}|}}\sqrt{\frac{|W_{i}|}{|W|}}\\\
\leq\frac{1}{\sqrt{c}}\sum_{i=1}^{N}\sqrt{\frac{|W_{i}|}{|W|}}.$

Letting $x_{i}=|W_{i}|/|W|$ and taking vectors
$u=\left(\sqrt{x_{1}},\dots,\sqrt{x_{N}}\right)$, $v=(1,\dots,1)$ we have that
$\sum_{i=1}^{N}x_{i}=1$ and so
$\left(\sum_{i=1}^{N}\sqrt{x_{i}}\right)^{2}=\left(u\cdot
v\right)^{2}\leq(u\cdot u)\,(v\cdot v)=\left(\sum_{i=1}^{N}x_{i}\right)N=N$ by
the Cauchy-Schwarz inequality. Hence

$\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}\leq\frac{\sqrt{N}}{\sqrt{c}}<1$

since $N\leq 9<c$.

[][] [scale=0.85*7] (1.95,4.5) rectangle (3.2,6); in 1,…,60 [very thin] ( 5+
(2*+2)*(5-7.5)/(2*+1) , 5 )–( 5+ (2*+2)*(7.5-7.5)/(2*+1) , 7.5 ); [black]
(2.47,5) – (2.53,5.06) – (2.53,5) – (2.47,5) ; (0,5) – (5,5);
[fill=white](0,0) – (5,10) – (5,0) – (0,0); [blue, thick] (1713/700,4.9)
–(2.49,5)– (2.8,5.9); [red, thick] (6639/2515,5.4348) – (739/275,613/110);
[red,fill=red,opacity=0.2] (2,5) – (30/14,5) – (50/18,50/9) – (40/14,40/7) –
(2,5); at (2.13, 5.05) $A_{4,3}^{k}$; at (6639/2515-0.05,5.4348+0.05) $W_{k}$;
at (2.4,4.95) $W_{\star}$; at (2.75,5.9) blue$W$black; [scale=0.5] at (2,5)
$\bullet$; at (2,4.95) $r_{2}$; [scale=0.5] at (30/14,5) $\bullet$; at
(30/14,4.95) $r_{3}$; [scale=0.5] at (50/18,50/9) $\bullet$; at
(50/18+0.05,50/9) $r_{4}$; [scale=0.5] at (40/14,40/7) $\bullet$; at
(40/14+0.05,40/7) $r_{1}$; at (30/14,4.7) $A_{1}\cap\sigma$; [][] at (0,0)
[scale=0.85*12] (4.96,6.6) rectangle (5.35,7.7); in 1,…,80 ( 2.5 +
(2*+1)*(5-5)/(2*) , 5 )–( 2.5 + (2*+1)*(10-5)/(2*) , 10 ); [fill=white] (5,10)
– (5,7.5) – (10,10) – (5,10); at (5.20,7.36) $A_{1,3}^{l}$; (5,7.5) – (10,10);
[fill=white] (5,10) rectangle (4.8,5); [red,fill=red,opacity=0.3]
(5+10/(4*8+2),7.5 - 50/(8*8+4) ) – (5+10/(4*8-2),7.5 - 50/(8*8-4) ) – (5, 7.5
+ 10/(8*8+4) ) – (5, 7.5 + 10/(8*8+12) ) – (5+10/(4*8+2),7.5 - 50/(8*8+4) );
[scale=0.5] at (5+10/(4*8+2),7.5 - 50/(8*8+4) ) $\bullet$; at (5+10/(4*8+2) -
0.02,7.5 - 50/(8*8+4) ) $r_{1}^{\prime}$; [scale=0.5] at (5+10/(4*8-2),7.5 -
50/(8*8-4) ) $\bullet$; at (5+10/(4*8-2)-0.02,7.5 - 50/(8*8-4) )
$r_{2}^{\prime}$; [scale=0.5] at (5, 7.5 + 10/(8*8+4) ) $\bullet$; at (5+0.02,
7.5 + 10/(8*8+4) ) $r_{3}^{\prime}$; [scale=0.5] at (5, 7.5 + 10/(8*8+12) )
$\bullet$; at (5-0.02, 7.5 + 10/(8*8+12) ) $r_{4}^{\prime}$; [fill=black]
(5,7.48) – (5,7.5) – (5+0.03,7.5+0.015) – (5,7.48); ; [dashed] (-0.4,1.4)
circle (0.6cm); [dashed] (-0.4+0.6*cos(315),1.4+0.6*sin(315)) –
(2.6+2*cos(160),0.1+2*sin(160)); [dashed, fill=white] (2.6,0.1) circle (2cm);
[scale=1.5] at (2.6+1.05,-1.2) red$U_{k}$black; (2.6,0.1) circle (2cm);
[thick] (2.6+2*cos(170),0.1+2*sin(170)) – (2.6+2*cos(80),0.1+2*sin(80));
[thick] (2.6+2*cos(250),0.1+2*sin(250)) – (2.6+2*cos(-5),0.1+2*sin(-5));
[red,fill=red,opacity=0.2] (1,2) – (2.5,2.6) – (4.5,-2.4) – (3,-3) – (1,2);
[red,very thick] (2.6+2*cos(110),0.1+2*sin(110)) –
(2.6+2*cos(295),0.1+2*sin(295)); [blue,very thick]
(2.6+1.48*cos(109),0.1+1.48*sin(109)) –
(2.6+1.24*cos(296.5),0.1+1.24*sin(296.5)); [scale=1.5] at (2.6+0.6,0.1+0.2)
blue$U_{k,l}$black;

Figure : Part (a) shows an unstable curve $W$ passing near to the accumulation
point $(1/4,1/2)$, split into $W_{\star}$ below $y=1/2$ and the collection
$W_{k}\subset A_{4,3}^{k}$. Part (b) shows the image
$U_{k}=H_{\sigma}(W_{k})\subset H_{\sigma}\left(A_{4,3}^{k}\right)$, which for
odd $k$ lies near the accumulation point $(1/2,3/4)$ and contains subcurves
$U_{k,l}\subset A_{1,3}^{l}$.

Condition (LABEL:eq:oneStep) holds for $H_{\sigma}^{2}$ when $W\cap
B_{\varepsilon}(P_{2})\neq\varnothing$ for all $\varepsilon>0$. We begin with
the case $W\cap B_{\varepsilon}((1/4,1/2))\neq\varnothing$ and let
$\varepsilon\to 0$. We may choose $\delta$ sufficiently small so that $W$
intersects $A_{1}\cap\sigma$ and some collection of sets $A_{4,3}^{k}$,
$k_{0}\leq k\leq k_{1}$, where $k_{1}\to\infty$ as $\varepsilon\to 0$,
$k_{0}\to\infty$ as $\delta\to 0$. Therefore, $\mathcal{S}$ splits $W$ into a
lower component $W_{\star}\subset A_{1}\cap\sigma$ and upper components
$W_{k}\subset A_{4,3}^{k}$, illustrated in Figure 1. We study how the images
of these components under $H_{\sigma}$ are split up by $\mathcal{S}$. Recall
the corners $r_{j}(k)$ which define $A_{4,3}^{k}$ near $(1/4,1/2)$. The curve
$W_{k}$ has endpoints on $r_{1}r_{2}$ and $r_{3}r_{4}$, and all tangent
vectors to $W$ lie in $\mathcal{C}_{1}$. For odd $k$ the image
$U_{k}=H_{\sigma}(W_{k})$ is a curve joining $r_{1}^{\prime}r_{2}^{\prime}$ to
$r_{3}^{\prime}r_{4}^{\prime}$, with tangent vectors aligned in
$M_{4}M_{3}^{k}\,\mathcal{C}_{1}$. This curve is split by $\mathcal{S}$ into
an upper portion $U_{k,\star}\subset A_{4}\cap\sigma$, and a collection
$U_{k,l}\subset A_{1,3}^{l}$ for some consecutive range $l_{0}\leq l\leq
l_{1}$ which depends on $k$. Each $A_{1,3}^{l}$ is bounded by the lines

$\mathcal{L}_{l}:y-\frac{1}{2}=\frac{2l}{2l+1}(x-1/4)$ (2)

and $\mathcal{L}_{l-1}$, hence a lower bound on $l_{0}(k)$ is given by the
largest $l$ such that $r_{2}^{\prime}(k)$ lies on or above
$\mathcal{L}_{l-1}$. One can verify that $r_{2}^{\prime}(k)$ lies on
$\mathcal{L}_{l-1}$ when $k=7l-4$ and approaches $(1/2,3/4)$ monotonically in
$x$ and $y$ so that $l_{0}(k)\geq\lfloor\frac{k+4}{7}\rfloor$. To determine an
upper bound on $l_{1}$, note that $r_{4}^{\prime}r_{1}^{\prime}$ lies on the
line

$y-\frac{3}{4}-\frac{1}{8k+12}=-\frac{6k+8}{2k+3}\left(x-\frac{1}{2}\right),$
(3)

meeting the $A_{4}$ boundary $\mathcal{L}:y=1/2+x/2$ at the point

$(x_{k},y_{k})=\left(\frac{7k+10}{14k+19},\frac{21k+29}{28k+38}\right).$ (4)

We similarly calculate that the line $\mathcal{L}_{l}$ meets $y=1/2+x/2$ at
the point

$(X_{l},Y_{l})=\left(\frac{l}{2l-1},\frac{1-3l}{2-4l}\right).$

The intersection of $U_{k}$ with $y=1/2+x/2$ must be some point $(x,1/2+x/2)$
with $x\geq x_{k}$ so that an upper bound on $l_{1}(k)$ is the smallest $l$
such that $x_{k}\geq X_{l}$, which reduces to $l\geq 7k+10$, hence
$l_{1}(k)\leq\lceil 7k+10\rceil=7k+10$. For even $k$ the splitting behaviour
is entirely analogous, with $H_{\sigma}(W_{k})$ intersecting $\mathcal{S}$ in
the neighbourhood of $(1,1/4)$. For the lower component $W_{\star}$, the image
$U_{\star}$ = $H_{\sigma}(W_{\star})=H(W_{\star})$ lies in a neighbourhood of
$H(1/4,1/2)=(1/4,1)$ and is split by $\mathcal{S}$ into a collection
$U_{\star,j}\subset A_{4,3}^{j}$, $j_{0}\leq j\leq j_{1}$, where
$j_{1}\to\infty$ as $\varepsilon\to 0$, $j_{0}\to\infty$ as $\delta\to 0$.
Write $W_{\star,j}=H_{\sigma}^{-1}(U_{\star,j})$,
$W_{k,\star}=H_{\sigma}^{-1}(U_{k,\star})$,
$W_{k,l}=H_{\sigma}^{-1}(U_{k,l})$, then $W$ splits into components

$W=\left(\bigcup_{j\geq j_{0}}W_{\star,j}\right)\cup\left(\bigcup_{k\geq
k_{0}}W_{k,\star}\right)\cup\left(\bigcup_{k\geq
k_{0}}\bigcup_{l=l_{0}}^{l_{1}}W_{k,l}\right)$ (5)

on which $DH_{\sigma}^{2}$ is constant. Let
$V_{i}=H_{\sigma}(U_{i})=H_{\sigma}^{2}(W_{i})$, then for $q=1/2$:

$\split\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}&=\sum_{i}\sqrt{\frac{|W_{i}|}{|V_{i}|}}\sqrt{\frac{|W_{i}|}{|W|}}\\\
=\sum_{j\geq
j_{0}}\sqrt{\frac{|W_{\star,j}|}{|V_{\star,j}|}}\sqrt{\frac{|W_{\star,j}|}{|W|}}+\sum_{k\geq
k_{0}}\sqrt{\frac{|W_{k,\star}|}{|V_{k,\star}|}}\sqrt{\frac{|W_{k,\star}|}{|W|}}+\sum_{k\geq
k_{0}}\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{|W_{k,l}|}{|V_{k,l}|}}\sqrt{\frac{|W_{k,l}|}{|W|}}\\\
\leq\sum_{j\geq
j_{0}}\sqrt{\frac{1}{\Lambda_{\star,j}}}\sqrt{\frac{|W_{\star,j}|}{|W|}}+\sum_{k\geq
k_{0}}\sqrt{\frac{1}{\Lambda_{k,\star}}}\sqrt{\frac{|W_{k,\star}|}{|W|}}+\sum_{k\geq
k_{0}}\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{1}{\Lambda_{k,l}}}\sqrt{\frac{|W_{k,l}|}{|W|}}$

by Lemma 1, where $\Lambda_{i}$ is the minimum expansion factor of
$DH_{\sigma}^{2}$ on $W_{i}$ over the cone $\mathcal{C}_{1}$. Define
$W_{\diamond}=W\setminus W_{\star}$ and let $0\leq p\leq 1$ denote the
proportion $|W_{\diamond}|=p\,|W|$, then

$\split\liminf_{\delta\to
0}\sup_{W:|W|<\delta}\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}\leq\sup_{0\leq
p\leq 1}\Bigg{(}&\lim_{j_{0}\to\infty}\sum_{j\geq
j_{0}}\sqrt{\frac{1}{\Lambda_{\star,j}}}\sqrt{\frac{(1-p)|W_{\star,j}|}{|W_{\star}|}}\\\
+\lim_{k_{0}\to\infty}\sum_{k\geq
k_{0}}\sqrt{\frac{1}{\Lambda_{k,\star}}}\sqrt{\frac{p\,|W_{k,\star}|}{|W_{\diamond}|}}\\\
+\lim_{k_{0}\to\infty}\sum_{k\geq
k_{0}}\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{1}{\Lambda_{k,l}}}\sqrt{\frac{p\,|W_{k,l}|}{|W_{\diamond}|}}\,\Bigg{)}.$

We put upper bounds on each of these sums using lower bounds on the expansion
factors $\Lambda_{i}$ and geometric bounds on the curves $U_{i}$ terminating
on $\mathcal{S}$. We use asymptotic notation $f\sim g$ for functions $f,g$ if
$f/g\to 1$, and write $f\lesssim g$ if there is some function $h$ such that
$f\leq h\sim g$. Starting with the first sum, $DH_{\sigma}^{2}$ is given by
$M_{4}M_{3}^{j}M_{1}=(-1)^{j}\left(\smallmatrix-2j-3&-6j-8\\\
6j+818j+21\right)$ on each component $W_{\star,j}$ with minimum expansion
factors given by

$\split\Lambda_{\star,j}&=\inf_{7/3\leq m\leq
3}\sqrt{\frac{\left(-2j-3-6jm-8m\right)^{2}+(6j+8+18jm+21)^{2}}{1+m^{2}}}\\\
\sim\inf_{7/3\leq m\leq
3}\sqrt{\frac{\left(2+6m\right)^{2}+(6+18m)^{2}}{1+m^{2}}}\,j=\frac{48\sqrt{145}}{29}\,j:=c_{\star}\,j.$
(6)

Each curve $U_{\star,j}$ has tangent vectors in $M_{1}\mathcal{C}_{1}$
satisfying $41/17\leq|v_{2}|/|v_{1}|\leq 17/7$, $v_{1}v_{2}\geq 0$. For each
$j>j_{0}$, $U_{\star,j}$ traverses $A_{4,3}^{j}$ so that (making a similar
calculation to equation LABEL:eq:sigma1bhk) Lemma  gives

$\frac{a_{\star}}{j^{2}}\lesssim|U_{\star,j}|\lesssim\frac{b_{\star}}{j^{2}}$
(7)

for $a_{\star}=\frac{13}{80}\sqrt{2}$,
$b_{\star}=\frac{41}{192}\sqrt{1+17^{2}/41^{2}}$ (calculated in the appendix,
section LABEL:sec:upperCalculationsAppendix). The upper bound also trivially
holds for $j=j_{0}$. Let $\Lambda_{1}^{+}$, $\Lambda_{1}^{-}$ denote the
maximum and minimum expansion factors of $M_{1}$ over $\mathcal{C}_{1}$, then

$|W_{\star}|=\sum_{j\geq j_{0}}|W_{\star,j}|\geq\sum_{j\geq
j_{0}}\frac{|U_{\star,j}|}{\Lambda_{1}^{+}}\gtrsim\frac{a_{\star}}{\Lambda_{1}^{+}}\sum_{j\geq
j_{0}+1}\frac{1}{j^{2}}\geq\frac{a_{\star}}{\Lambda_{1}^{+}(j_{0}+1)}\sim\frac{a_{\star}}{\Lambda_{1}^{+}j_{0}}$

where we have used the fact that

$\frac{1}{j^{2}}\geq\frac{1}{j(j+1)}=\frac{1}{j}-\frac{1}{j+1}$

and considered the telescoping sum. Similarly $|W_{\star,j}|\lesssim
b_{\star}/\left(\Lambda_{1}^{-}j^{2}\right)$ so that

$\frac{|W_{\star,j}|}{|W_{\star}|}\lesssim\frac{b_{\star}\Lambda_{1}^{+}j_{0}}{a_{\star}\Lambda_{1}^{-}}\,\frac{1}{j^{2}}.$

Hence

$\split\sum_{j\geq
j_{0}}\sqrt{\frac{1}{\Lambda_{\star,j}}}\sqrt{\frac{(1-p)|W_{\star,j}|}{|W_{\star}|}}&\lesssim\sum_{j\geq
j_{0}}\sqrt{\frac{1}{c_{\star}j}}\sqrt{\frac{(1-p)b_{\star}\Lambda_{1}^{+}j_{0}}{a_{\star}\Lambda_{1}^{-}}\,\frac{1}{j^{2}}}\\\
=\sqrt{\frac{(1-p)b_{\star}\Lambda_{1}^{+}j_{0}}{c_{\star}a_{\star}\Lambda_{1}^{-}}}\sum_{j\geq
j_{0}}j^{-3/2}\\\
\leq\sqrt{\frac{(1-p)b_{\star}\Lambda_{1}^{+}}{c_{\star}a_{\star}\Lambda_{1}^{-}}}\,2\sqrt{\frac{j_{0}}{j_{0}-1}}\\\
\to
2\sqrt{\frac{(1-p)b_{\star}\Lambda_{1}^{+}}{c_{\star}a_{\star}\Lambda_{1}^{-}}}$
(8)

as $j_{0}\to\infty$, where we have used $\sum_{j\geq
j_{0}}j^{-3/2}\leq\int_{j_{0}-1}^{\infty}x^{-3/2}\d{x}$. Moving onto the next
summation, $\Lambda_{k,\star}$ is determined by
$M_{4}^{2}M_{3}^{k}=(-1)^{k}\left(\smallmatrix 14k+5&-14k-12\\\
-34k-1234k+29\right)$ and satisfies

$\Lambda_{k,\star}\sim\inf_{7/3\leq m\leq
3}\sqrt{\frac{(14-14m)^{2}+(34-34m)^{2}}{1+m^{2}}}\,k=\frac{104}{\sqrt{29}}\,k=:c_{\diamond}\,k.$

For odd $k$ the curve $U_{k,\star}$ has endpoints $(1/2,y_{0})$ on
$r_{3}^{\prime}r_{4}^{\prime}$ and $(2y_{1}-1,y_{1})$ on $\mathcal{L}$, where
$y_{1}$ is bounded by the intersections of $r_{4}^{\prime}r_{1}^{\prime}$ and
$r_{3}^{\prime}r_{2}^{\prime}$ with $\mathcal{L}$ (see Figure 1). An upper
bound on $|y_{1}-y_{0}|$ is given by taking $(1/2,y_{0})=r_{3}^{\prime}$ and
$(2y_{1}-1,y_{1})=r_{4}^{\prime}r_{1}^{\prime}\cap\mathcal{L}$. Noting (4),
this gives

$|y_{1}-y_{0}|\leq\frac{3}{4}+\frac{1}{8k+4}-\frac{21k+29}{28k+38}=\frac{6k+9}{56k^{2}+104k+38}\sim\frac{6}{56}k^{-1}.$

$U_{k,\star}$ has tangent vectors $(v_{1},v_{2})^{T}$ in
$M_{4,3}^{k}\,\mathcal{C}_{1}\subset\mathcal{C}_{4}$, so that
$|v_{2}/v_{1}|\geq 7/3$. By Lemma  we then have $|U_{k,\star}|\lesssim
b_{\diamond}/k$, where $b_{\diamond}=\frac{6}{56}\sqrt{1+\frac{9}{49}}\approx
0.117$. The minimum expansion factor of
$M_{4}M_{3}^{k}=(-1)^{k}\left(\smallmatrix 2k+1&-2k-2\\\ -6k-26k+5\right)$
over $\mathcal{C}_{1}$ is given by

$\inf_{7/3\leq m\leq
3}\sqrt{\frac{(2-2m)^{2}+(6-6m)^{2}}{1+m^{2}}}\,k=\frac{8\sqrt{145}}{29}\,k=:\gamma\,k$

which gives $|W_{k,\star}|\lesssim b_{\diamond}/\left(\gamma k^{2}\right)$.
The analysis for even $k$ is analogous and gives the same upper bound. We next
require a lower bound on $|W_{\diamond}|$. For $k>k_{0}$, $W_{k}$ is a curve
with tangent vectors in $\mathcal{C}_{1}$ which traverses $A_{4,3}^{k}$.
Making the same calculation as (LABEL:eq:sigma1bhk), $|W_{k}|$ is bounded
below by the shortest path across $A_{4,3}^{k}$, the line segment passing
through $r_{4}$ with gradient $3$. That’s

$|W_{k}|\geq\sqrt{\left(\frac{1}{16k^{2}+28k+6}\right)^{2}+\left(\frac{3}{16k^{2}+28k+6}\right)^{2}}\sim\frac{\sqrt{10}}{16k^{2}}$
(9)

so that

$|W_{\diamond}|\geq\sum_{k\geq k_{0}+1}|W_{k}|\gtrsim\frac{a}{k_{0}}$ (10)

with $a:=\sqrt{10}/16$. Hence

$\sum_{k\geq
k_{0}}\sqrt{\frac{1}{\Lambda_{k,\star}}}\sqrt{\frac{p\,|W_{k,\star}|}{|W_{\diamond}|}}\lesssim\sum_{k\geq
k_{0}}\sqrt{\frac{1}{c_{\diamond}k}}\sqrt{\frac{pb_{\diamond}k_{0}}{a\gamma}\,\frac{1}{k^{2}}}\to
2\sqrt{\frac{pb_{\diamond}}{c_{\diamond}a\gamma}}$

as $k_{0}\to\infty$, by a similar argument to (8). For the third summation,
$\Lambda_{k,l}$ is determined by the matrix

$M_{1}M_{3}^{l}M_{4}M_{3}^{k}=(-1)^{k+l}\pmatrix{-}48kl-10k-18l-3&48kl+10k+42l+8\\\
-112kl-26k-42l-8112kl+26k+98l+21$

and satisfies (for large $k,l$)

$\Lambda_{k,l}\sim\inf_{7/3\leq m\leq
3}\sqrt{\frac{(48-48m)^{2}+(112-112m)^{2}}{1+m^{2}}}\,kl=64\,kl=:c\,kl.$

We can show an upper bound $|U_{k,l}|\lesssim b/l^{2}$ where
$b=\frac{3}{32}\sqrt{1+\frac{9}{49}}\approx 0.102$ (see section
LABEL:sec:upperCalculationsAppendix) so that $|W_{k,l}|\lesssim b/(\gamma
kl^{2})$. Now by (10),

$\split\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{1}{\Lambda_{k,l}}}\sqrt{\frac{p|W_{k,l}|}{|W_{\diamond}|}}&\lesssim\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{1}{c\,kl}}\sqrt{\frac{pbk_{0}}{a\gamma
kl^{2}}}\\\ \leq\sqrt{\frac{1}{c\,k}}\sqrt{\frac{pbk_{0}}{a\gamma
k}}\sum_{l=\lfloor\frac{k+4}{7}\rfloor}^{7k+10}l^{-3/2}\\\ \leq
2\sqrt{\frac{1}{c\,k}}\sqrt{\frac{pbk_{0}}{a\gamma
k}}\left(\frac{1}{\sqrt{\lfloor\frac{k+4}{7}\rfloor-1}}-\frac{1}{\sqrt{7k+10}}\right)\\\
\sim 2\sqrt{\frac{1}{c\,k}}\sqrt{\frac{pbk_{0}}{a\gamma
k}}\left(\sqrt{7}-\frac{1}{\sqrt{7}}\right)\frac{1}{\sqrt{k}}.$

Letting $h=(\sqrt{7}-1/\sqrt{7})^{2}=36/7$, we have that

$\split\sum_{k\geq
k_{0}}\sum_{l=l_{0}}^{l_{1}}\sqrt{\frac{1}{\Lambda_{k,l}}}\sqrt{\frac{p|W_{k,l}|}{|W_{\diamond}|}}&\lesssim
2\sqrt{\frac{pbhk_{0}}{ca\gamma}}\sum_{k\geq k_{0}}k^{-3/2}\\\ \to
4\sqrt{\frac{pbh}{ca\gamma}}$

as $k_{0}\to\infty$. Hence for $q=1/2$

$\liminf_{\delta\to
0}\sup_{W:|W|<\delta}\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}\leq\sup_{0\leq
p\leq
1}\left(2\sqrt{\frac{(1-p)b_{\star}\Lambda_{1}^{+}}{c_{\star}a_{\star}\Lambda_{1}^{-}}}+2\sqrt{\frac{pb_{\diamond}}{c_{\diamond}a\gamma}}+4\sqrt{\frac{pbh}{ca\gamma}}\right).$
(11)

It is simple to show that for $s,t>0$ the function
$f(p)=s\sqrt{1-p}+t\sqrt{p}$ always attains its maximum value at
$p=t^{2}/(s^{2}+t^{2})$. Hence letting

$s=2\sqrt{\frac{b_{\star}\Lambda_{1}^{+}}{c_{\star}a_{\star}\Lambda_{1}^{-}}}\approx
0.450,\quad
t=2\sqrt{\frac{b_{\diamond}}{c_{\diamond}a\gamma}}+4\sqrt{\frac{bh}{ca\gamma}}\approx
0.639$

gives

$\liminf_{\delta\to
0}\sup_{W:|W|<\delta}\sum_{i}\left(\frac{|W|}{|V_{i}|}\right)^{q}\frac{|W_{i}|}{|W|}\leq
s\sqrt{\frac{s^{2}}{s^{2}+t^{2}}}+t\sqrt{\frac{t^{2}}{s^{2}+t^{2}}}\approx
0.781<1$

as required. The analysis is analogous for $W$ near $(1/4,1)$ and extends to
$W$ near $(3/4,0)$ and $(3/4,1/2)$ using the symmetry $T(x,y)=(1-x,y+1/2)$
which commutes with $H_{\sigma}$ (as seen in the proof of Lemma
LABEL:lemma:unstableGrowth). Equivalent analysis verifies the two step
expansion for curves near the other accumulation points $p\in P_{1}$. We
provide the relevant calculations to the appendix, Proposition
LABEL:prop:lowerTwoStep. We are now ready to apply Theorem
LABEL:thm:chernovZhang.

##  Decay of correlations

The return map $H_{\sigma}:\sigma\to\sigma$ enjoys exponential decay of
correlations. In particular it admits a Young tower with base $\Delta_{0}$
satisfying the exponential tail bound

$\mu\left(\\{z\in\sigma\,|\,R(z,H_{\sigma},\Delta_{0})>n\\}\right)\leq\mathrm{const}\,\theta^{n}$
(12)

for all $n\geq 1$ where $\theta<1$ is some constant. We run through the
conditions for applying Theorem LABEL:thm:chernovZhang. Invariance of the
unstable and stable cone fields $C_{z}^{u}$, $C_{z}^{s}$ was the subject of
section 1, satisfying (H1.1). Condition (H1.2) follows by taking
$\Lambda=\sqrt{85/41}$, with this lower bound attained by considering the
expansion of $M_{2}^{-1}$ over the cone boundary of $\mathcal{C}_{2}^{s}$ with
gradient $-8/10$. Noting Remark LABEL:remark:firstRemark, we next show
(H1.3’). The cone fields are continuous over the components $\sigma_{j}$ of
$\Omega\setminus\mathcal{S}_{0}$, indeed they are constant. Noting that all of
the stable cones lie within $\mathcal{C}_{1}^{s}$ and all of the unstable
cones lie in $\mathcal{C}$, a positive angle between stable and unstable cone
fields follows from $\mathcal{C}_{1}^{s}\cap\mathcal{C}=\varnothing$. (H2) was
the subject of section 1. Unstable manifolds provide a class of $H_{\sigma}$
invariant unstable curves which satisfy the regularity conditions listed in
[chernov_statistical_2009]. Piecewise linearity of the map trivially implies
their bounded curvature and bounds on distortion; absolute continuity follows
from Lemma 1. (H4) follows from Proposition 1, where we showed that
$H_{\sigma}$ is Bernoulli with respect to the normalised Lebesgue measure on
$\sigma$. Noting Remark LABEL:remark:firstRemark, (H5) follows for the map
$H_{\sigma}^{2}$ by Propositions , 1, LABEL:prop:lowerTwoStep.


# Continuous-variable assisted thermal quantum simulation

Dan-Bo Zhang  Guangdong Provincial Key Laboratory of Quantum Engineering and
Quantum Materials, GPETR Center for Quantum Precision Measurement and SPTE,
South China Normal University, Guangzhou 510006, China  Frontier Research
Institute for Physics, South China Normal University, Guangzhou 510006, China
Guo-Qing Zhang  Guangdong Provincial Key Laboratory of Quantum Engineering and
Quantum Materials, GPETR Center for Quantum Precision Measurement and SPTE,
South China Normal University, Guangzhou 510006, China  Frontier Research
Institute for Physics, South China Normal University, Guangzhou 510006, China
Zheng-Yuan Xue  Guangdong Provincial Key Laboratory of Quantum Engineering and
Quantum Materials, GPETR Center for Quantum Precision Measurement and SPTE,
South China Normal University, Guangzhou 510006, China  Frontier Research
Institute for Physics, South China Normal University, Guangzhou 510006, China
Shi-Liang Zhu<EMAIL_ADDRESS>Guangdong Provincial Key Laboratory of Quantum
Engineering and Quantum Materials, GPETR Center for Quantum Precision
Measurement and SPTE, South China Normal University, Guangzhou 510006, China
Frontier Research Institute for Physics, South China Normal University,
Guangzhou 510006, China  Z. D. Wang<EMAIL_ADDRESS>Department of Physics and
HKU-UCAS Joint Institute for Theoretical and Computational Physics at Hong
Kong, The University of Hong Kong, Pokfulam Road, Hong Kong, China  Frontier
Research Institute for Physics, South China Normal University, Guangzhou
510006, China

###### Abstract

Simulation of a quantum many-body system at finite temperatures is crucially
important but quite challenging. Here we present an experimentally feasible
quantum algorithm assisted with continuous-variable for simulating quantum
systems at finite temperatures. Our algorithm has a time complexity scaling
polynomially with the inverse temperature and the desired accuracy. We
demonstrate the quantum algorithm by simulating finite temperature phase
diagram of the Kitaev model. It is found that the important crossover phase
diagram of the Kitaev ring can be accurately simulated by a quantum computer
with only a few qubits and thus the algorithm may be readily implemented on
current quantum processors. We further propose a protocol implementable with
superconducting or trapped ion quantum computers.

_Introduction.–_ Simulation of quantum many-body systems has been an incentive
for building quantum computers Feynman (1982). Recent advances of scaling-up
quantum processors enable us to simulate steady states and dynamics of a
quantum system at zero temprature to a larger size Barends  _et al._ (2015);
Bernien  _et al._ (2017); Kandala  _et al._ (2017); Zhang  _et al._ (2017),
and remarkably, to make a simulation of quantum phase transition to be a
reality Islam  _et al._ (2011); Zhang  _et al._ (2017). However, simulation of
a quantum many-body system at finite temperatures is even more challenging as
the state in a quantum computer is usually a pure quantum state, while
simulating a quantum system at finite temperatures requires to simulate a kind
of mixed states, namely quantum thermal (Gibbs) states. The relevance of
thermal states is ubiquitous, and its simulaiton is not only important for
physics itself, such as understanding high-Tc superconductivity Lee  _et al._
(2006) , but also can provide quantum speedup for optimization Van Apeldoorn
_et al._ (2017).

Thermal quantum simulation (TQS) requires a well control of both quantum
coherence and temperature, which challenges the current quantum platforms. A
quantum algorithm based on quantum phase estimation may involve a large number
of auxiliary qubits and complicated quantum circuits, which is not suitable
for near-term quantum simulators Terhal and DiVincenzo (2000); Poulin and
Wocjan (2009); Temme  _et al._ (2011); Riera  _et al._ (2012). Recent hybrid
quantum-classical variational algorithms require less quantum resource and are
feasible in implementation Wu and Hsieh (2019); Verdon  _et al._ (2019); Liu
_et al._ (2019); Chowdhury  _et al._ (2020); Wang  _et al._ (2020), but it
should be trained for each Hamiltonian at every temperature and thus is not a
general solution. Moreover, it waits for a guarantee of quantum advantages in
time complexity.

An alternative way is to use continuous variables (CV and also called qumode)
for encoding and processing high-density information by exploiting its
infinite dimensional Hilbert space Weedbrook  _et al._ (2012); Lau  _et al._
(2017). Notably, a hybrid approach of incorporating both qubits and qumodes
has been shown to have a potential advantage to make the best-of-both-worlds
Furusawa and Van Loock (2011); Andersen  _et al._ (2015); Liu  _et al._
(2016); Gan  _et al._ . Moreover, the mainstream platforms of quantum
computers often have naturally existing continuous variables, such as motional
modes of trapped ions Leibfried  _et al._ (2003); Monroe and Kim (2013) and
cavity modes of superconducting circuits Wallraff  _et al._ (2004); Paik  _et
al._ (2011); Devoret and Schoelkopf (2013), making a hybrid variable approach
physical realizable. This may enable us to design hybrid-variable quantum
algorithms for TQS with both quantum advantage and feasible physical
implementation.

In this paper, we present a quantum algorithm for thermal quantum simulation
assisted with an auxiliary qumode, which has a time complexity scaling
polynomially with the inverse temperature $\beta$ and the desired accuracy
$\epsilon$. The algorithm converts thermal information, encoded in the CV
resource state parameterized with $\beta$, into temperature of the quantum
system. Moreover, by revealing an equivalence relation of the quantum
algorithm, we further propose adaptive TQS, allowing TQS at varied $\beta$
with a proper-chosen resource state, which enables flexibility of algorithm
design in practice. To show the power of adaptive TQS, we consider thermal
states of two textbook models, the quantum Ising model Sachdev  and Kitaev
ring Kitaev (2001), in the quantum critical regime, and demonstrate our
algorithm can accurately determine the crossover temperature. This indicates
an interplay between quantum and thermal fluctuations, which underlines the
quantum criticality at finite temperatures, can be faithfully captured. An
intriguing result here is that the important crossover phase diagram of the
Kitaev model in periodic condition (Kitaev ring) can be accurately simulated
by a quantum computer with only a few qubits and thus the algorithm may be
readily implemented on current quantum processors. We also propose an
experimental protocol implementable with a superconducting or trapped ion
quantum computer. Our work opens a new avenue for simulating finite
temperature quantum systems by exploiting the power of continuous variables.

_Thermal quantum simulation._ Consider a quantum system described by a
Hamiltonian $H$, which can be mapped to a quantum computer with $N$ qubits.
The energies and eigenstates of the Hamiltonian governed by the Schrodinger
equation $H|u_{n}\rangle=E_{n}|u_{n}\rangle$ ($n=0,1,...,D-1$,where
$D=2^{N}$). At a finite temperature $T=1/\beta$, the system in equilibrium is
in a quantum thermal state $\rho(\beta)=e^{-\beta H}/Z(\beta)$, where
$Z(\beta)=\text{Tr}e^{-\beta H}$ is the partition function. Our goal is to
prepare a pure quantum state $|\psi(\beta)\rangle$ in an enlarger Hilbert
space which has the property
$\rho(\beta)=\text{Tr}_{A}|\psi(\beta)\rangle\langle\psi(\beta)|$, where
$\text{Tr}_{A}$ denotes the partial trace of some ancillary degrees of freedom
addressed later. One can verify that, partial trace of the second partite of
the thermofield double (TFD) state defined as

$|\psi(\beta)\rangle=\sum_{n}\frac{e^{-\beta
E_{n}/2}}{\sqrt{\mathcal{Z}(\beta)}}|u_{n}\rangle\otimes|u^{*}_{n}\rangle$ (1)

is just $\rho(\beta)$. We propose a quantum process as follows,

$|\psi(\beta)\rangle=\sqrt{\mathcal{C}}e^{-\beta H/2}\otimes
I|\psi(0)\rangle,$ (2)

where $\mathcal{C}$ is a normalization factor. Here
$|\psi(0)\rangle=\frac{1}{\sqrt{D}}\sum_{n}|n\rangle\otimes|n\rangle$ is an
infinite temperature TFD which can be written as a product of N copies of Bell
state (see Supplementary Material (SM) for a derivation), and thus is easy to
prepare. The central task is then to construct $e^{-\beta H/2}$.

We introduce an auxiliary qumode to represent a nonunitary $e^{-\beta H/2}$ as
an integral of unitaries Lau  _et al._ (2017); Arrazola  _et al._ (2019);
Zhang  _et al._ (2019, 2020), which extends the linear-combination-of-
unitaries to the case of CV Gui-Lu (2006); Childs  _et al._ (2017); Van
Apeldoorn  _et al._ (2017); Gilyen  _et al._ (2018); Arrazola  _et al._
(2019). Note that

$e^{-\beta
h/2}=\int_{-\infty}^{\infty}dpR(\beta,p)e^{-ihp},~{}~{}R(\beta,p)=\frac{2}{\pi}\frac{\beta}{\beta^{2}+4p^{2}}$
(3)

holds for $h\geq 0$. It implies that we shall add a constant to $H$ to
guarantee a nonnegative ground state energy in the algorithm. Then in the
basis $\\{|u_{n}\rangle\\}$ that $H$ is diagonal, one can verify that

$\displaystyle e^{-\beta
H/2}=\int_{-\infty}^{\infty}dpR(\beta,p)e^{-iHp}\propto\langle
0_{q}|e^{-iH\hat{p}}|R(\beta)\rangle,$ (4)

where
$|R(\beta)\rangle=\sqrt{\beta\pi}\int_{-\infty}^{\infty}dpR(\beta,p)|p\rangle_{p}$.
Here we denote $|p\rangle_{p}$ ($|q\rangle_{q}$) as basis of continuous
variable quadrature $\hat{p}$ ($\hat{q}$). Eq. (4) shows a scheme that the
nonunitary operator $e^{-\beta H/2}$ can be implemented on a quantum computer,
assisted with an ancillary qumode. The qumode is prepared at
$|R(\beta)\rangle$, evolves jointly with the system by $e^{-iH\hat{p}}$, and
is finally projected onto $|0\rangle_{q}$. As the information related to
temperature is encoded in $|R(\beta)\rangle$, we may call it as a resouce
state.

Figure 1: An illustration of preparing quantum thermal state $\rho(\beta)$
from N copies of Bell states, through coupling the system with a qumode by
unitary evolution $e^{-iH\hat{p}/2}$. The qumode is initialed at the resource
state $|R(\beta)\rangle$ and is finally projected onto $|0\rangle_{q}$.
Thermal state $\rho(\beta)$ is obtained by discarding (tracing out) the
additional $N$ ancillary qubits.

_Quantum algorithm.–_ We now present the procedure of quantum algorithm for
preparing thermal quantum simulation. The quantum algorithm is probabilistic
since it postselects the quadrature $\hat{q}$ to zero. In practice, the
measurement should have a finite precision, which reduces the accuracy while
raises the success rate. We model this effect by projecting to a squeezing
state $|0,s\rangle=s^{-\frac{1}{2}}\pi^{-\frac{1}{4}}\int
dp|p\rangle_{p}e^{-p^{2}/2s^{2}}$, which squeezes the quadrature $\hat{q}$ by
a factor $s$. The quantum algorithm, as illustrated in Fig. 1, has three
steps:

  1. 1.

_State preparation._ Prepare N copies of Bell states,
$|\psi(0)\rangle=\prod_{i=1}^{N}B_{i,i+N}$, where
$B_{i,i+N}=\frac{1}{\sqrt{2}}(|0_{i}0_{i+N}\rangle+|1_{i}1_{i+N}\rangle)$ is a
Bell state. A qumode is initialized in a resource state $|R(\beta)\rangle$.
The total system is in a state
$|\Psi_{0}\rangle=|\psi(0)\rangle\otimes|R(\beta)\rangle$. Therefore, to
simulate a Hamiltonian $H$ which is encoded in $N$ qubits, our scheme needs
$2N$ qubits and one qumode.

  2. 2.

_Unitary evolution._ A constant is added to $H$ to make positivity of the
spectrum. then the unitary evolution $e^{-iH\hat{p}}\otimes I$ is implemented
and it couples the system with the qumode. The unitary can be decomposed with
Trotter-Suzuki formula, which will be discussed later. The state turns to be

$|\Psi_{1}\rangle=\sum_{n}\int_{-\infty}^{\infty}dp\phi_{n}(\beta,p)|u_{n}\rangle\otimes|u^{*}_{n}\rangle\otimes|p\rangle_{p},$
(5)

where $\phi_{n}(\beta,p)=\sqrt{\frac{\beta\pi}{D}}R(\beta,p)e^{-E_{n}p}$.

  3. 3.

_Projection._ Project the qumode onto a squeezing state $|0,s\rangle$, the
unnormalized final state (discarding the qumode) is

$|\tilde{\psi}(\beta)\rangle=\sum_{n}a(E_{i},\beta,s)|u_{n}\rangle\otimes|u^{*}_{n}\rangle,$
(6)

where
$a(E,\beta,s)=s^{-\frac{1}{2}}\pi^{-\frac{1}{4}}\int_{-\infty}^{\infty}dp\phi_{n}(\beta,p)e^{-\frac{p^{2}}{2s^{2}}}$.
The success rate is $O(\frac{\mathcal{Z}(\beta)}{sD})$.

Note at $s\rightarrow\infty$ limit, $a(E,\beta,s)\propto e^{-\beta E/2}$,
which exactly equals to the quantum thermal state at an inverse temperature
$\beta$.

One issue of the quantum algorithm is that the resource state can not be
produced for free, e.g., $R(\beta\rightarrow 0)$ and
$R(\beta\rightarrow\infty)$ can not be efficiently prepared since they are
infinitely squeezed. To solve this problem, we reveal an equivalent relation
of the quantum algorithm, which says that $\rho(\beta)$ can be simulated with
different pairs of resource states and unitary evolutions, namely

$\left(R(\beta_{0}),e^{-iH\hat{p}\frac{\beta}{\beta_{0}}}\right)\mapsto\rho(\beta),$
(7)

where $\beta_{0}$ is adjustable. The equivalent relation is a direct
consequence of an invariance under $\beta\rightarrow a\beta$ and $H\rightarrow
H/a$ in Eq. (4).

The equivalent relation allows us to use a fixed resource state to simulate
the thermal state at varied $\beta$, which we call as adaptive TQS. From the
perspective of quantum resource theory Chitambar and Gour (2019), the
equivalent relation may reveal a conversion between static resource (preparing
resource state) and dynamical resource (constructing unitary $e^{-iH\hat{p}}$)
and their trade-off for thermal quantum simulation. It enriches the
flexibility and feasibility of algorithmic implementation , as difficulties of
resource state preparation and Hamiltonian evolution vary on different quantum
platforms.

_Time complexity.–_ We give a runtime analysis for a relative error $\epsilon$
in the partition function. As it is not practical to ignore the cost of
preparing resource state (e.g., at $\beta\rightarrow 0,\infty$ limits), we
consider adaptive TQS with a resource state $|R(\beta_{0})\rangle$ that is
easy to prepare. The runtime then relies on the circuit complexity of
constructing the unitary operator $e^{-iH\hat{p}\beta/\beta_{0}}$ and the
success rate.

For a Hamiltonian $H=\sum_{i=1}^{M}c_{i}H_{i}$ with local terms, one can refer
to Trotter-Suzuki formula Lloyd (1996) to decompose
$e^{-iH\hat{p}\beta/\beta_{0}}$. Typically, it includes terms like evolutions
of
$\sigma^{\alpha}_{i}\hat{p},\sigma^{\alpha}_{i}\sigma^{\gamma}_{i+1}\hat{p}$
($\alpha,\gamma=x,y,z$), which can be viewed as parity-dependent displacement
operator and can be decomposed as basic gates of qubits and a hybrid gate
$e^{i\theta\sigma^{x}\hat{p}}$ (see SM). The circuit complexity is
$O(M^{3}\epsilon^{-1}\beta^{2})$, and may be improved with more advanced
techniques Low and Chuang (2017); Campbell (2019); Childs  _et al._ (2019).
For a precision $\epsilon$, the squeezing factor should be
$s=O(\epsilon^{-\frac{1}{2}})$. Using amplitude amplification G. Brassard and
Tapp (2002), the success rate becomes
$O(\sqrt{\frac{\mathcal{Z}(\beta)}{sD}})$, and the algorithm should be run
repeatedly with
$O(\sqrt{\frac{D}{\mathcal{Z}(\beta)}}\epsilon^{-\frac{1}{4}})$ times. In
total, the time complexity is
$O(\sqrt{\frac{D}{\mathcal{Z}(\beta)}}M^{3}\beta^{2}\epsilon^{-\frac{5}{4}})$,
which improves polynomially from methods using quantum phase estimation Poulin
and Wocjan (2009). In addition, we found that the quantum algorithm gives a
bound from below for the partition function (see SM).

Figure 2: Performance of preparing thermal states of single qubit with two
different approaches: $\beta$ ($\tilde{\beta}$) stands for results from TQS
(adaptive TQS).

_Demonstration with single qubit.–_ To warm up, we simulate thermal state of a
single qubit to highlight some features of the quantum algorithm. For the
numeral simulation, we develop a classical simulator for hybrid-variable
quantum computing based on the open-source Qutip Johansson  _et al._ (2013).
The Hamiltonian is $H=g(\sigma_{x}+c)$ (here $g>0$ and $c=1$ enforces all
eigenvalues are nonnegative). After initialing the cavity mode into a resource
state and preparing a Bell state, an unitary $e^{-ig(\sigma^{x}+1)\hat{p}t}$
performs on the cavity mode and the system qubit, and finally the cavity mode
is squeezed and projected onto zero photon. We use two approaches for
preparing thermal states at $\beta=1,2,3$. The first fixes $t=1$ and uses
resource states $|R(\beta)\rangle$ with $\beta=1,2,3$, respectively. The
second is adaptive TQS, fixing the resource state $|R(\beta=1)\rangle$ and
using $t=1,2,3$, respectively. Due to finite squeezing at the projection,
those simulated thermal states can only approximate the exact thermal states.
We use trace distance to measure the precision,
$d_{r}(\beta)=\frac{1}{2}\text{Tr}|\tilde{\rho}(\beta)-\rho(\beta)|$. In Fig.
2, we can see that trace distances decrease rapidly with an increase of
squeezing factor, at the cost of decreasing success rate, which are expected.
Using adaptive TQS, the increasing of precision can be faster, and the trace
distance can be small even at small $s$. Moreover, a moderate truncation of
phonon number (e.g., up to $7$ even-number Fock states) can reach a good
precision (see SM). The above results indicate that we can use an adaptive
approach with approximated resource state for thermal quantum simulation.

_Simulation of the finite-temperature phase diagram.–_ We now consider our TQS
with two textbook models, the Kitaev ring (a spinless p-wave superconductor)
Kitaev (2001) and quantum Ising model Sachdev . The Hamiltonian of the Kitaev
ring reads

$H_{K}=-J\sum_{i=1}^{L}(c_{i}^{\dagger}c_{i+1}+c_{i}^{\dagger}c_{i+1}^{\dagger}+h.c.)-\mu\sum_{i=1}^{L}c_{i}^{\dagger}c_{i},$
(8)

where fermionic operators $c_{L+1}=c_{1}$ as we consider the periodic
condition. The model has a topological phase transition at
$\lambda_{K}\equiv\mu/2J=1$. Using Jordan-Wigner transformation,
$c_{i}=\prod_{j=1}^{i-1}\sigma_{j}^{z}\sigma_{i}^{-},~{}c_{i}^{\dagger}=\sigma_{i}^{+}\prod_{j=1}^{i-1}\sigma_{j}^{z},~{}c_{i}^{\dagger}c_{i}=\frac{1}{2}(\sigma_{i}^{z}-1)$,
the Hamiltonian of the Kitaev ring can be mapped into a spin model

$H_{S}=-h\sum_{i}^{L}\sigma_{i}^{z}-J\sum_{i}^{L-1}\sigma_{i}^{x}\sigma_{i+1}^{x}-J\sigma_{1}^{y}P_{L}\sigma_{L}^{y}+E_{0},$
(9)

where $h=-\frac{\mu}{2}$, $P_{L}=\prod_{i=2}^{L-1}\sigma_{i}^{z}$ is a string
operator, $E_{0}$ is added for assuring nonnegative spectrum (which is
demanded for the quantum algorithm). The spin model $H_{S}$ become the quantum
Ising model if $\sigma_{1}^{y}P_{L}\sigma_{L}^{y}$ is replaced by
$\sigma_{L}^{x}\sigma_{1}^{x}$, and it has a phase transition at
$\lambda_{S}\equiv h/J=1$. This difference of the Ising and Kitaev models has
a big effect on the quantum criticality at small sizes, as addressed below.
The finite-temperature phase diagrams of these two models in the infinite
lattice size are the same, which have an important V-shape crossover structure
Sachdev  as shown in Fig. 3.

From numeral calculation (see SM for details), it is shown that almost
$L\approx 80$ qubits is required for the quantum Ising chain to show well
shaped crossover temperatures, as seen in Fig 3, and $L$ should be larger
closer to the QPT point. For the Kitaev ring, in contrast, the temperature
crossovers are very close in shape for a quite large range even for small $L$.
The subtle lies in that Kitaev ring always has an energy gap
$\Delta=J|1-\lambda_{K}|$, while the quantum Ising model will have low-lying
in-gap states when $\lambda_{S}<1$, leading to small crossover temperatures
for $\lambda_{S}<1$, which is obvious for $L=10,20$. Thus, we refer to a small
size Kitaev ring for simulating finite-temperature phase diagram with quantum
computers.

Figure 3: Finite temperature phase diagram of the one-dimensional Kitaev ring
and quantum Ising model.

To address the feasibility of the physical implementation, we adopt adaptive
TQS for simulating thermal states of the Kitaev ring, using a resource state
$|R(\beta=4)\rangle$ and a squeezing factor $s=10$. The crossover temperature
for each $\lambda_{K}$ is determined by the temperature that maximizes the
magnetic susceptibility $\chi_{m}$ (corresponding to fluctuation of fermion
number in the Kitaev ring), which can be obtained by measuring the
magnetization
$\mathcal{M}=\frac{1}{L}\sum_{i}\text{Tr}\rho(\beta)\sigma_{i}^{z}$ and then
calculating the magnetic susceptibility
$\chi_{m}=\frac{\partial\mathcal{M}}{\partial h}$ using a finite difference.
The result is shown in Fig. 4. For $L=2$, we compare results of different
truncation for the resource state ($N_{c}$ is the number of Fock basis of
even-number photons), and it can be seen that the $V$-shape crossover
temperature approaches the exact one when increasing $N_{c}$. We also simulate
$L=2,3,4,5$, which demonstrates very close crossover temperatures. Therefore,
a remarkable advantage of our quantum algorithm is that observation of the
important crossover temperature phase diagram of the Kitaev ring needs only a
few qubits.

Figure 4: Simulation of crossover temperature for the Kitaev ring in the
quantum critical regime with adaptive TQS. (a). Crossover temperatures that
use different truncating of photon number $N_{c}$ for the resource state
$|R(\beta=4)\rangle$.(b).Crossover temperature for different lattice size,
$L=2,3,4,5$.

_Experimental realization.–_ We now discuss physical implementation of the
quantum algorithm, which relies on hybrid variable quantum computing.
Promising candidates include trapped ions Leibfried  _et al._ (2003); Zhu  _et
al._ (2006); Häffner  _et al._ (2008); Monroe and Kim (2013); Zhang  _et al._
(2018) and superconducting circuits Krastanov  _et al._ (2015); Li  _et al._
(2017); Hofheinz  _et al._ (2008); Heeres  _et al._ (2015); Li  _et al._
(2017); Hu  _et al._ (2019); Song  _et al._ (2019), etc. We take
superconducting circuit system as an example, for their well controllable CV
cavity mode and its coupling to the qubits. The scheme can be
straightforwardly used in trapped ion quantum computer.

Firstly, the resource state can be expanded in a Fock space,
$|R(\beta)\rangle=\sum_{n=0}^{\infty}r_{n}|2n\rangle$, where
$r_{n}=\int_{-\infty}^{\infty}dq\sqrt{\frac{\beta}{2}}e^{-\beta|q|/2}H_{2n}(q)e^{-\frac{q^{2}}{2}}$
and $H_{n}(q)$ is the $n$-th order Hermite function. With a truncation of
photon number, the state can be prepared in a cavity by superposing even-
number photons. This can be achieved with a sequence of qubit rotation and
Jaynes-cumming type qubit-cavity coupling Law and Eberly (1996); Hofheinz  _et
al._ (2009), or using the number-dependent arbitrary phase gate and
displacement operator in the dispersive regime Heeres  _et al._ (2015); Wang
_et al._ (2017). Note that resource states around $\beta=2$ can be
approximated with a few components of small number photons, and are feasible
in both approaches. Thus, we can chose one such resource state and use it for
adaptive TQS. Secondly, construction of $e^{-iH\hat{p}t}$ can be compiled into
basic qubit gates and only one hybrid variable gate
$e^{i\theta\sigma^{x}\hat{p}}$, as discussed before. All are standard quantum
operations in superconducting circuit system, and remarkably, the hybrid gate
$e^{i\theta\sigma^{x}\hat{p}}$ can be readily realized in the strong coupling
limit Niemczyk  _et al._ (2010); Forn-Díaz  _et al._ (2017, 2019). Thirdly,
projection onto a squeezing state $|0,s\rangle$ can be implemented by first
performing a squeezing on the CV mode, and then post-selecting the CV mode to
the vacuum state (zero photon state). Further, we can measure the system to
access the thermal state by quantum state tomography, or studying quantum
statistical mechanism by measuring physical quantities, such as heat capacity,
magnetic susceptibility, etc. As for the Kitaev ring, we note that evolution
of nonlocal term $\sigma_{1}^{y}P_{L}\sigma_{L}^{y}\hat{p}$ can be constructed
efficiently (see SM). Thus, all above well-developed quantum operations can
render a feasible implementing protocol for thermal quantum simulation of the
Kitaev ring, in order to illuminate the novel quantum critical regime on small
quantum processors.

_Summary.–_ We have proposed a quantum algorithm for thermal quantum
simulation assisted with an auxiliary CV resource state. We have confirmed its
power by simulating finite-temperature phase diagram of the Kitaev ring, and
found that the important crossover phase diagram of the model can be
accurately determined by a quantum computer with only a few qubits. Thus, our
work may pave the way for studying finite temperature quantum systems in
experiments.

###### Acknowledgements.

This work was supported by the Key-Area Research and Development Program of
GuangDong Province (Grant No. 2019B030330001), the National Key Research and
Development Program of China (Grant No. 2016YFA0301800), the National Natural
Science Foundation of China (Grants No. 91636218 and No. U1801661), the Key
Project of Science and Technology of Guangzhou (Grant No. 201804020055), and
the CRF of Hong Kong (No. C6005-17G).

## References

  * Feynman (1982) Richard P. Feynman, “Simulating physics with computers,” Int. J. Theor. Phys. 21, 467–488 (1982).
  * Barends  _et al._ (2015) R. Barends, L. Lamata, J. Kelly, L. García-Álvarez, A. G. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Yu Chen, Z. Chen, B. Chiaro, A. Dunsworth, I. C. Hoi, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, E. Solano, and John M. Martinis, “Digital quantum simulation of fermionic models with a superconducting circuit,” Nature Communications 6, 7654 (2015).
  * Bernien  _et al._ (2017) Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S. Zibrov, Manuel Endres, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin, “Probing many-body dynamics on a 51-atom quantum simulator,” Nature 551, 579–584 (2017).
  * Kandala  _et al._ (2017) Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242–246 (2017).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z. X. Gong, and C. Monroe, “Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator,” Nature 551, 601–604 (2017).
  * Islam  _et al._ (2011) R. Islam, E. E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G. D. Lin, L. M. Duan, C. C. Joseph Wang, J. K. Freericks, and C. Monroe, “Onset of a quantum phase transition with a trapped ion quantum simulator,” Nature Communications 2, 377 (2011).
  * Lee  _et al._ (2006) Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen, “Doping a mott insulator: Physics of high-temperature superconductivity,” Rev. Mod. Phys. 78, 17–85 (2006).
  * Van Apeldoorn  _et al._ (2017) J. Van Apeldoorn, A. Gilyen, S. Gribling, and R. de Wolf, “Quantum sdp-solvers: better upper and lower bounds,” 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS). Proceedings , 403–414 (2017).
  * Terhal and DiVincenzo (2000) Barbara M. Terhal and David P. DiVincenzo, “Problem of equilibration and the computation of correlation functions on a quantum computer,” Phys. Rev. A 61, 022301 (2000).
  * Poulin and Wocjan (2009) David Poulin and Pawel Wocjan, “Sampling from the thermal quantum gibbs state and evaluating partition functions with a quantum computer,” Phys. Rev. Lett. 103 (2009), 10.1103/PhysRevLett.103.220502.
  * Temme  _et al._ (2011) K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, “Quantum metropolis sampling,” Nature 471, 87–90 (2011).
  * Riera  _et al._ (2012) Arnau Riera, Christian Gogolin, and Jens Eisert, “Thermalization in nature and on a quantum computer,” Phys. Rev. Lett. 108, 080402 (2012).
  * Wu and Hsieh (2019) Jingxiang Wu and Timothy H. Hsieh, “Variational thermal quantum simulation via thermofield double states,” Phys. Rev. Lett. 123, 220502 (2019).
  * Verdon  _et al._ (2019) Guillaume Verdon, Jacob Marks, Sasha Nanda, Stefan Leichenauer, and Jack Hidary, “Quantum hamiltonian-based models and the variational quantum thermalizer algorithm,” (2019), arXiv:1910.02071 [quant-ph] .
  * Liu  _et al._ (2019) Jin-Guo Liu, Liang Mao, Pan Zhang, and Lei Wang, “Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits,” (2019), arXiv:1912.11381 [quant-ph] .
  * Chowdhury  _et al._ (2020) Anirban N. Chowdhury, Guang Hao Low, and Nathan Wiebe, “A variational quantum algorithm for preparing quantum gibbs states,” (2020), arXiv:2002.00055 [quant-ph] .
  * Wang  _et al._ (2020) Youle Wang, Guangxi Li, and Xin Wang, “Variational quantum gibbs state preparation with a truncated taylor series,” (2020), arXiv:2005.08797 [quant-ph] .
  * Weedbrook  _et al._ (2012) Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84, 621–669 (2012).
  * Lau  _et al._ (2017) H. K. Lau, R. Pooser, G. Siopsis, and C. Weedbrook, “Quantum machine learning over infinite dimensions,” Phys. Rev. Lett. 118, 080501 (2017).
  * Furusawa and Van Loock (2011) Akira Furusawa and Peter Van Loock,  _Quantum teleportation and entanglement: a hybrid approach to optical quantum information processing_ (John Wiley & Sons, 2011).
  * Andersen  _et al._ (2015) U. L. Andersen, J. S. Neergaard-Nielsen, P. van Loock, and A. Furusawa, “Hybrid discrete- and continuous-variable quantum information,” Nat. Phys. 11, 713–719 (2015).
  * Liu  _et al._ (2016) N. N. Liu, J. Thompson, C. Weedbrook, S. Lloyd, V. Vedral, M. L. Gu, and K. Modi, “Power of one qumode for quantum computation,” Phys. Rev. A 93, 052304 (2016).
  * (23) H. C. J. Gan, Gleb Maslennikov, Ko-Wei Tseng, Chihuan Nguyen, and Dzmitry Matsukevich, “Hybrid quantum computation gate with trapped ion system,”  1908.10117 .
  * Leibfried  _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, “Quantum dynamics of single trapped ions,” Rev. Mod. Phys. 75, 281–324 (2003).
  * Monroe and Kim (2013) C. Monroe and J. Kim, “Scaling the ion trap quantum processor,” Science 339, 1164–1169 (2013).
  * Wallraff  _et al._ (2004) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R. S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, “Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics,” Nature 431, 162–167 (2004).
  * Paik  _et al._ (2011) Hanhee Paik, D. I. Schuster, Lev S. Bishop, G. Kirchmair, G. Catelani, A. P. Sears, B. R. Johnson, M. J. Reagor, L. Frunzio, L. I. Glazman, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, “Observation of high coherence in josephson junction qubits measured in a three-dimensional circuit qed architecture,” Phys. Rev. Lett. 107, 240501 (2011).
  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, “Superconducting circuits for quantum information: An outlook,” Science 339, 1169–1174 (2013).
  * (29) S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, Cambridge, England, 1999).
  * Kitaev (2001) A Yu Kitaev, “Unpaired majorana fermions in quantum wires,” Physics-Uspekhi 44, 131–136 (2001).
  * Arrazola  _et al._ (2019) Juan Miguel Arrazola, Timjan Kalajdzievski, Christian Weedbrook, and Seth Lloyd, “Quantum algorithm for nonhomogeneous linear partial differential equations,” Phys. Rev. A 100, 032306 (2019).
  * Zhang  _et al._ (2019) Dan-Bo Zhang, Zheng-Yuan Xue, Shi-Liang Zhu, and Z. D. Wang, “Realizing quantum linear regression with auxiliary qumodes,” Phys. Rev. A 99, 012331 (2019).
  * Zhang  _et al._ (2020) Dan-Bo Zhang, Shi-Liang Zhu, and Z. D Wang, “Protocol for implementing quantum nonparametric learning with trapped ions,” Phys. Rev. Lett. 124, 010506 (2020).
  * Gui-Lu (2006) Long Gui-Lu, “General quantum interference principle and duality computer,” Commun. Theor. Phys. 45, 825–844 (2006).
  * Childs  _et al._ (2017) Andrew M. Childs, Robin Kothari, and Rolando D. Somma, “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,” Siam. J. Comput. 46, 1920–1950 (2017).
  * Gilyen  _et al._ (2018) A. Gilyen, Su Yuan, Low Guang Hao, and N. Wiebe, “Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics arxiv,” arXiv , 67 pp.–67 pp. (2018).
  * Chitambar and Gour (2019) Eric Chitambar and Gilad Gour, “Quantum resource theories,” Rev. Mod. Phys. 91, 025001 (2019).
  * Lloyd (1996) Seth Lloyd, “Universal quantum simulators,” Science 273, 1073–1078 (1996).
  * Low and Chuang (2017) Guang Hao Low and Isaac L. Chuang, “Optimal hamiltonian simulation by quantum signal processing,” Phys. Rev. Lett. 118, 010501 (2017).
  * Campbell (2019) Earl Campbell, “Random compiler for fast hamiltonian simulation,” Phys. Rev. Lett. 123 (2019), 10.1103/PhysRevLett.123.070503.
  * Childs  _et al._ (2019) Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu, “A theory of trotter error,” (2019), arXiv:1912.08854 [quant-ph] .
  * G. Brassard and Tapp (2002) M. Mosca G. Brassard, P. Høyer and A. Tapp,  _Quantum Amplitude Amplification and Estimation_ (Contemporary Mathematics Series Millenium Volume 305, 2002).
  * Johansson  _et al._ (2013) J. R. Johansson, P. D. Nation, and Franco Nori, “Qutip 2: A python framework for the dynamics of open quantum systems,” Comput. Phys. Commun. 184, 1234–1240 (2013).
  * Zhu  _et al._ (2006) Shi-Liang Zhu, C. Monroe, and L. M. Duan, “Trapped ion quantum computation with transverse phonon modes,” Phys. Rev. Lett. 97, 050505 (2006).
  * Häffner  _et al._ (2008) H. Häffner, C. F. Roos, and R. Blatt, “Quantum computing with trapped ions,” Physics Reports 469, 155–203 (2008).
  * Zhang  _et al._ (2018) Junhua Zhang, Mark Um, Dingshun Lv, Jing-Ning Zhang, Lu-Ming Duan, and Kihwan Kim, “Noon states of nine quantized vibrations in two radial modes of a trapped ion,” Phys. Rev. Lett. 121, 160502 (2018).
  * Krastanov  _et al._ (2015) Stefan Krastanov, Victor V. Albert, Chao Shen, Chang-Ling Zou, Reinier W. Heeres, Brian Vlastakis, Robert J. Schoelkopf, and Liang Jiang, “Universal control of an oscillator with dispersive coupling to a qubit,” Phys. Rev. A 92, 040303 (2015).
  * Li  _et al._ (2017) Linshu Li, Chang-Ling Zou, Victor V. Albert, Sreraman Muralidharan, S. M Girvin, and Liang Jiang, “Cat codes with optimal decoherence suppression for a lossy bosonic channel,” Phys. Rev. Lett. 119, 030502 (2017).
  * Hofheinz  _et al._ (2008) Max Hofheinz, E. M. Weig, M. Ansmann, Radoslaw C. Bialczak, Erik Lucero, M. Neeley, A. D. O’Connell, H. Wang, John M. Martinis, and A. N. Cleland, “Generation of fock states in a superconducting quantum circuit,” Nature 454, 310–314 (2008).
  * Heeres  _et al._ (2015) Reinier W. Heeres, Brian Vlastakis, Eric Holland, Stefan Krastanov, Victor V. Albert, Luigi Frunzio, Liang Jiang, and Robert J. Schoelkopf, “Cavity state manipulation using photon-number selective phase gates,” Phys. Rev. Lett. 115, 137002 (2015).
  * Hu  _et al._ (2019) L. Hu, Y. Ma, W. Cai, X. Mu, Y. Xu, W. Wang, Y. Wu, H. Wang, Y. P. Song, C. L. Zou, S. M. Girvin, L. M. Duan, and L. Sun, “Quantum error correction and universal gate set operation on a binomial bosonic logical qubit,” Nat. Phys. 15, 503–508 (2019).
  * Song  _et al._ (2019) Chao Song, Kai Xu, Hekang Li, Yu-Ran Zhang, Xu Zhang, Wuxin Liu, Qiujiang Guo, Zhen Wang, Wenhui Ren, Jie Hao, Hui Feng, Heng Fan, Dongning Zheng, Da-Wei Wang, H. Wang, and Shi-Yao Zhu, “Generation of multicomponent atomic schrödinger cat states of up to 20 qubits,” Science 365, 574–577 (2019).
  * Law and Eberly (1996) C. K. Law and J. H. Eberly, “Arbitrary control of a quantum electromagnetic field,” Phys. Rev. Lett. 76, 1055–1058 (1996).
  * Hofheinz  _et al._ (2009) Max Hofheinz, H. Wang, M. Ansmann, Radoslaw C. Bialczak, Erik Lucero, M. Neeley, A. D. O’Connell, D. Sank, J. Wenner, John M. Martinis, and A. N. Cleland, “Synthesizing arbitrary quantum states in a superconducting resonator,” Nature 459, 546–549 (2009).
  * Wang  _et al._ (2017) W. Wang, L. Hu, Y. Xu, K. Liu, Y. Ma, Shi-Biao Zheng, R. Vijay, Y. P Song, L. M. Duan, and L. Sun, “Converting quasiclassical states into arbitrary fock state superpositions in a superconducting circuit,” Phys. Rev. Lett. 118, 223604 (2017).
  * Niemczyk  _et al._ (2010) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross, “Circuit quantum electrodynamics in the ultrastrong-coupling regime,” Nat. Phys. 6, 772–776 (2010).
  * Forn-Díaz  _et al._ (2017) P. Forn-Díaz, J. J García-Ripoll, B. Peropadre, J. L. Orgiazzi, M. A Yurtalan, R. Belyansky, C. M Wilson, and A. Lupascu, “Ultrastrong coupling of a single artificial atom to an electromagnetic continuum in the nonperturbative regime,” Nat. Phys. 13, 39–43 (2017).
  * Forn-Díaz  _et al._ (2019) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, “Ultrastrong coupling regimes of light-matter interaction,” Rev. Mod. Phys. 91, 025005 (2019).
  * Nielsen and Chuang (2010) Michael A. Nielsen and Isaac L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).

##  Appendix A Supplementary Material

##  Appendix B Thermofield double state.

We prove that the purified thermal state at finite temperatures can be written
as $N$ pairs of Bell states. Denote $\\{E_{n}\\}$ and $\\{|u_{n}\rangle\\}$
($n=0,1,...,2^{N}-1$) as eigenvalues and eigenstates for a Hamiltonian $H$
which can be encoded with $N$ qubits. The system of $H$ at infinite
temperatures is described by a complete mixed state. To purify this complete
mixed state, we use a thermofield double (TFD) state at infinite temperatures,

$|\psi(\beta=0)\rangle\equiv\frac{1}{\sqrt{D}}\sum_{n}|u_{n}\rangle\otimes|u^{*}_{n}\rangle=\frac{1}{\sqrt{D}}\sum_{n}|n\rangle\otimes|n\rangle,$
(10)

where $D=2^{N}$ is the dimension of Hilbert space. We now prove the equation.
Note that both $\\{|u_{n}\rangle\\}$ and $\\{|n\rangle\\}$ form a complete set
of basis for the Hilbert space of $N$ qubits, we can find an unitary
transformation that $|u_{n}\rangle=\sum_{n}U_{nm}|m\rangle$. Then,

$\displaystyle\sum_{n}|u_{n}\rangle\otimes|u^{*}_{n}\rangle$ $\displaystyle=$
$\displaystyle\sum_{nmm^{\prime}}U_{nm}U^{*}_{nm^{\prime}}|m\rangle\otimes|m^{\prime}\rangle$
(11) $\displaystyle=$
$\displaystyle\sum_{mm^{\prime}}\delta_{mm^{\prime}}|m\rangle\otimes|m^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\sum_{n}|n\rangle\otimes|n\rangle$

We prove that $|\psi(0)\rangle=\prod_{i=1}^{N}B_{i,i+N}$, where
$B_{i,i+N}=\frac{1}{\sqrt{2}}(|0_{i}0_{i+N}\rangle+|1_{i}1_{i+N}\rangle)$.
With a binary expression $n=n_{1}n_{2}..,n_{N}$ ($n_{i}=0,1$), we can write,

$\displaystyle|\psi(\beta=0)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{D}}\sum_{n}|n\rangle\otimes|n\rangle$ (12)
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{D}}\sum_{n_{1}n_{2}...n_{N}}|n_{1}n_{2}...n_{N}\rangle\otimes|n_{1}n_{2}...n_{N}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{D}}\prod_{i=1}^{N}\sum_{n_{i}=0,1}|n_{i}\rangle_{i}\otimes|n_{i}\rangle_{i+N}$
$\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\frac{1}{\sqrt{2}}(|0_{i}0_{i+N}\rangle+|1_{i}1_{i+N}\rangle)$

It should be reminded that other purification protocols are possible, for
instance,
$|\psi(\beta=0)\rangle\equiv\frac{1}{\sqrt{D}}\sum_{n}|u_{n}\rangle\otimes|n\rangle$.

##  Appendix C Exact free energy

For the Kitaev ring and the quantum Ising model under periodic boundary
condition (PBC), we give methods to exactly calculate their free energy at
different parameters and temperatures. Specified attention should be given to
the quantum Ising model, since its energy spectrum should be calculated in the
even-parity and odd-parity subspaces, respectively. From the free energy the
magnetic susceptibility can be obtained.

Kitaev ring. For the Kitaev ring, the Hamiltonian after diagonalization in
k-space with Bogoliubov transformation can be written as,

$H_{K}=\sum_{i=0}^{L-1}\xi(k_{i})(\eta_{k_{i}}^{\dagger}\eta_{k_{i}}-\frac{1}{2}),$
(13)

where the single-particle spectrum is
$\xi(k)=2J\sqrt{1+\lambda_{K}^{2}-2\lambda_{K}\cos{k}}$, and we set
$\lambda_{K}=2\mu/J$. For this non-interacting fermionic model that follows
Fermi-Dirac statistics, the partition function can be easily calculated as,

$Z(\beta)=\prod_{i=0}^{L-1}(e^{-\beta\xi(k_{i})/2}+e^{\beta\xi(k_{i})/2}).$
(14)

Then, the free energy per site is,

$f(\beta)=-\beta^{-1}\log{Z(\beta)}/L=-\beta^{-1}\left[-\frac{1}{L}\sum_{i=0}^{L-1}\log{\cosh{\frac{\beta\xi(k_{i})}{2}}}-\log{2}\right].$
(15)

At the thermodynamical limit, it turns to be,

$f(\beta)=-\beta^{-1}\log{Z(\beta)}/L=-\beta^{-1}\left[\frac{1}{2\pi}\int_{0}^{2\pi}dk\log{\cosh{\frac{\beta\xi(k)}{2}}}-\log{2}\right].$
(16)

Quantum Ising model. The quantum Ising model under PBC is,

$H_{I}=-h\sum_{i}^{L}\sigma_{i}^{z}-J\sum_{i}^{L-1}\sigma_{i}^{x}\sigma_{i+1}^{x},$
(17)

where $\sigma_{L+1}^{\alpha}=\sigma_{1}^{\alpha}$. Using the Jordan-Wigner
transformation, it becomes a fermionic Hamiltonian,

$H_{I}=-J\sum_{i=1}^{L-1}(c_{i}^{\dagger}c_{i+1}+c_{i}^{\dagger}c_{i+1}^{\dagger}+h.c.)-\frac{h}{2}\sum_{i=1}^{L}c_{i}^{\dagger}c_{i}-J(-1)^{\hat{N}}(c_{L}-c_{L}^{\dagger})(c_{1}+c_{1}^{\dagger}),$
(18)

where $\hat{N}=\sum_{i=1}^{L}c_{i}^{\dagger}c_{i}$. Compared with the Kitaev
ring, there is a global phase $(-1)^{\hat{N}}$ that depends on the parity of
fermion number. Thus, we should discuss even and odd parity subspace
separately.

  1. 1.

Even-parity subspace. We should have $c_{L+1}=-c_{1}$ in accordance with the
third term of Eq. (18). This enforces $k_{i}$ to take values $k_{i}=\frac{2\pi
i}{L}+\pi/L$,($i=0,1,...,L-1$). By Bogoliubov transformation, the Hamiltonian
can be diagonalized as,

$H^{e}_{I}=\sum_{i=0}^{L-1}\xi(k_{i})(\eta_{k_{i}}^{\dagger}\eta_{k_{i}}-\frac{1}{2}),$
(19)

where $\xi(k)=2J\sqrt{1+\lambda_{S}^{2}-2\lambda_{S}\cos{k}}$, and we have
introduced $\lambda_{S}=h/J$. It should be emphasized that only even number of
excitations can allow for this Hamiltonian due to the constraint of even-
parity subspace. This makes that a formula as in Eq. (14) can not apply here.

  2. 2.

Odd-parity subspace. We have $c_{L+1}=c_{1}$, and $k_{i}$ to take values
$k_{i}=\frac{2\pi i}{L}$,($i=0,1,...,L-1$). The Hamiltonian in this subspace
writes,

$H^{o}_{I}=\sum_{i=1}^{L-1}\xi(k_{i})(\eta_{k_{i}}^{\dagger}\eta_{k_{i}}-\frac{1}{2})-2J(1-\lambda_{S})(c_{0}^{\dagger}c_{0}-\frac{1}{2}).$
(20)

The last term comes from the fact that a spinless p-wave pairing can not
happen at $k=0$. Similarly, only odd-number excitations are allowed for the
Hamiltonian $H^{o}_{I}$.

With Hamiltonians in both even-parity and odd-parity subspaces, we can get all
eigenstates for the quantum Ising model, and then calculate the partition
function. However, this can be intractable when $L$ is large, e.g., $L=50$,
because there is an exponential growth number of terms to be added.

We use a recursion relation for the partition function which is scalable for
large $L$. Let $J$ and $G$ are contributions to the partition function from
the even and odd parity subspace, separately. Then, $Z=J+G$. We discuss how to
calculate $J$, and $G$ can be calculated similarly. Denote $J^{e}_{n}$
($J^{o}_{n}$) as a partition function that comes from even (odd) number
single-particle excitations for modes in $\eta_{k_{i}}^{\dagger}$ that $0\leq
i\leq n$. Then, we have a recursion relation,

$J^{e}_{n}=J^{e}_{n-1}+e^{-\beta\xi(k_{n})}J^{o}_{n-1},~{}~{}~{}J^{o}_{n}=e^{-\beta\xi(k_{n})}J^{e}_{n-1}+J^{o}_{n-1},~{}~{}~{}J^{e}_{1}=e^{-\beta
E_{0}},~{}~{}~{}J^{o}_{1}=e^{-\beta E_{0}}e^{-\beta\xi(k_{0})},$ (21)

where $E_{0}=-\frac{1}{2}\sum_{i=0}^{L-1}\xi(k_{i})$ is the ground state
energy. By writing the recursion relation in a matrix form, we can calculate
$J=J^{e}_{L}$ very efficiently. Note that for the odd-parity subspace we
should get $G=G^{o}_{L}$, correspondingly.

With obtained partition function $Z=J+G$, the free energy can be got using
$F=-T\log{Z}$. In practice, we can take $e^{-\beta E_{0}}$ outside the
expression $Z$ since it can be a very large number.

##  Appendix D More simulation results.

Effect of truncation of photon number.For practice, we consider a truncation
of the Fock space and the resource state thus superposes a finite set of
photon numbers. We show that a moderate truncation is enough for thermal
quantum simulation with good accuracy, which is shown in Fig. 5, illustrated
by thermal state of a single qubit.

Figure 5: Effects on truncation of resource state in the Fock space for
preparing thermal states of single qubit with two different approaches:
$\beta$ ($\tilde{\beta}$) stands for results from TQS (adaptive TQS). Here,
$N_{c}$ is the number of Fock basis of even-number photons, for instance, for
$N_{c}=4$, $|0\rangle,|2\rangle,|4\rangle,|6\rangle$ are used.

Free energy. Accuracy simulation of free energy is key for thermal numeral
simulation. Here, we present numeral results for free energies for the Kitaev
ring for different $\lambda$ and $T$, as well as how it depends on the
squeezing factor. As seen in Fig. 6, free energy can be simulated very
accurately. Moreover, free energies at small squeezing factors are larger than
exact results, which is consisted to theoretical analysis that says the
quantum algorithm gives a bound from above for the free energy.

Figure 6: Free energies for the Kitaev ring at $L=2$. Dash lines correspond to
exact results and markers are results from adaptive TQS. The left figure is
obtained for a squeezing factor $s=10$, and the right figure is obtained at
temperature $T=0.4$. All are calculated under a truncation of $N_{c}=35$.

##  Appendix E Decomposition of hybrid-variable quantum operators.

To simulate $e^{-iH\hat{p}t}$, we use a trotterization to decompose it into
short time evolutions, evolving operators $e^{-i\delta\sigma^{i}_{z}\hat{p}}$,
$e^{-i\delta\sigma^{i}_{x}\sigma^{i+1}_{x}\hat{p}}$,
$e^{-i\delta\sigma^{1}_{y}P_{L}\sigma^{L}_{y}\hat{p}}$, and so on. We show
that all those operators can be compiled into basic quantum gates of qubits
and only one qubit-qumode hybrid quantum gate $e^{-i\delta\sigma_{x}\hat{p}}$,
as illustrated in Fig. 7, which is similar to the technique in the textbook
Nielsen and Chuang (2010). For convenience, we also denote
$X=\sigma^{x},Y=\sigma^{y},Z=\sigma^{z}$. Here, we can use some basic single-
qubit transformations,

$HZH=X,~{}~{}~{}HXH=Z,~{}~{}~{}HYH=-Y,~{}~{}~{}SXS^{\dagger}=Y.$ (22)

And those evolve CNOT gates,

$CZ_{2}C=Z_{1}Z_{2},~{}~{}~{}CZ_{1}C=Z_{1},~{}~{}~{},CX_{1}C=X_{1}X_{2},~{}~{}~{}CX_{2}C=X_{2},$
(23)

where we denote $C$ as CNOT gate with qubit 1 as control and qubit 2 as
target. We remark $e^{-i\delta\sigma^{i}_{x}\sigma^{i+1}_{x}\hat{p}}$ can be
seen as parity-dependence displacement operator, where the direction of
displacement depends on the parity $\sigma^{i}_{x}\sigma^{i+1}_{x}$.

For the Kitaev ring, specifically evolutions of
$\sigma_{1}^{y}P_{L}\sigma_{L}^{y}\hat{p}$ can be decomposed, following the
method in Fig. (7), where we have shown the case $L=2,3$. For other $L$, one
can use an identity,

$e^{-i\delta\sigma_{1}^{y}\sigma_{2}^{z}...\sigma_{L}^{y}\hat{p}}=S_{1}H_{1}S_{L}H_{L}He^{-i\delta\sigma_{1}^{z}\sigma_{2}^{z}...\sigma_{L}^{z}\hat{p}}H_{1}H_{L}S_{1}^{\dagger}S_{L}^{\dagger},$
(24)

where $e^{-i\delta\sigma_{1}^{z}\sigma_{2}^{z}...\sigma_{L}^{z}\hat{p}}$ can
be constructed with the same fashion as
$e^{-i\delta\sigma_{1}^{z}\sigma_{2}^{z}}$ and
$e^{-i\delta\sigma_{1}^{z}\sigma_{2}^{z}\sigma_{3}^{z}}$.

Figure 7: Decomposition of hybrid quantum gates with a series of CNOT and a
qubit-qumode gate $e^{-i\delta\sigma_{x}\hat{p}}$.

##  Appendix F Requirement of squeezing factor

The task is to calculate the requirement of $s$ with the required accuracy
$\epsilon$ for simulating quantum thermal state with the inverse temperature
$\beta$. For this we use a relative accuracy $\epsilon$ of both the partition
function and the quantum Gibbs distribution. It is expected that the former is
less demanding for squeezing factor under relative accuracy $\epsilon$. In the
following, we find that both may give the same dependence of the squeezing
factor $s$ on $\epsilon$. We first consider TQS for $\beta$ using
$|R(\beta)\rangle$, assuming that an arbitrary resource state can be prepared.
Then, we consider adaptive TQS with an easy-to-prepare resource state.

A relative accuracy $\epsilon$ of the partition function is defined as

$e(\beta,s):=\frac{|\tilde{\mathcal{Z}}(\beta,s)-\mathcal{Z}(\beta)|}{\mathcal{Z}(\beta)}\leq\epsilon.$
(25)

Here
$\tilde{\mathcal{Z}}(\beta,s)\propto\int_{E_{\text{0}}}^{E_{\text{1}}}a^{2}(E,\beta,s)D(E)dE$,
where $D(E)$ is spectral density of the system. Note that
$\tilde{\mathcal{Z}}(\beta,s)-\mathcal{Z}(\beta)$ can be approximated as

$\displaystyle\int_{E_{\text{1}}}^{E_{\text{0}}}D(E)dE\int_{-\infty}^{\infty}dp\frac{2}{\pi}\frac{\beta}{\beta^{2}+4p^{2}}(e^{-\frac{p^{2}}{2s^{2}}}-1)e^{-iEp}$
$\displaystyle\times\int_{-\infty}^{\infty}dp\frac{2}{\pi}\frac{\beta}{\beta^{2}+4p^{2}}(e^{-\frac{p^{2}}{2s^{2}}}+1)e^{-iEp}$
$\displaystyle=-16\int_{E_{\text{0}}}^{E_{\text{1}}}dED(E)(e^{\frac{\beta^{2}}{8s^{2}}}-1)(e^{\frac{\beta^{2}}{8s^{2}}}+1)e^{-\beta
E},$ (26)

where in the third line we have applied Cauchy’s residue theorem in complex
analysis, with an integral of $p$ in the half complex plane enclosing the pole
$p=-i\beta/2$ (the other pole is not used since the integral must depend on
$E$). We note that Eq. (F) is small than zero, which suggests that the quantum
algorithm would give a bound from below for the partition function, and
consequently a bound from above for the free energy.

For small $s/\beta$, using $e^{x}\approx 1+x$, then

$\displaystyle
e(\beta,s)\sim\frac{2\beta^{2}}{s^{2}}\frac{\int_{E_{\text{0}}}^{E_{\text{1}}}dED(E)e^{-\beta
E}}{\int_{E_{\text{0}}}^{E_{\text{1}}}dED(E)e^{-\beta
E}}=\frac{2\beta^{2}}{s^{2}}$ (27)

It is interesting to note that the final expression is independent of the
details of the spectral function. For TQS of $\beta$ using a resource state
$|R(\beta)\rangle$, to achieve $e(\beta,s)\leq\epsilon$ at inverse temperature
$\beta$, the required squeezing factor is $s\sim\beta\epsilon^{-\frac{1}{2}}$.

We now consider adaptive TQS at inverse temperature $\beta$ with a resource
state $|R(\beta_{0})\rangle$. Eq. (F) then is,

$\displaystyle\int_{E_{\text{1}}}^{E_{\text{0}}}D(E)dE\int_{-\infty}^{\infty}dp\frac{2}{\pi}\frac{\beta_{0}}{\beta_{0}^{2}+4p^{2}}(e^{-\frac{p^{2}}{2s^{2}}}-1)e^{-iEp\frac{\beta}{\beta_{0}}}$
(28)
$\displaystyle\times\int_{-\infty}^{\infty}dp\frac{2}{\pi}\frac{\beta_{0}}{\beta_{0}^{2}+4p^{2}}(e^{-\frac{p^{2}}{2s^{2}}}+1)e^{-iEp\frac{\beta}{\beta_{0}}}$
$\displaystyle=-16\int_{E_{\text{0}}}^{E_{\text{1}}}dED(E)(e^{\frac{\beta_{0}^{2}}{8s^{2}}}-1)(e^{\frac{\beta_{0}^{2}}{8s^{2}}}+1)e^{-\beta
E}.$

One can get $s\sim\epsilon^{-\frac{1}{2}}$, which is independent of $\beta$.
This can explain why adaptive TQS is less demanding for squeezing. The cost,
however, is that the complexity to construct the unitary evolution increases
with a factor $\beta^{2}$.

For a relative error $\epsilon$ of the Gibbs distribution, it requires that

$e(\beta,s):=\frac{|a^{2}(E,\beta,s)/\tilde{\mathcal{Z}}(\beta,s)-e^{-\beta
E}/\mathcal{Z}(\beta,s)|}{e^{-\beta E}/\mathcal{Z}(\beta,s)}\leq\epsilon.$

It can be verified that a scaling of $s\sim\beta\epsilon^{-\frac{1}{2}}$
($s\sim\epsilon^{-\frac{1}{2}}$) can satisfy this relative error for TQS
(adaptive TQS). In fact, following the same technique of relative error of
partition function, the same scaling of squeezing factor with $\beta$ and
$\epsilon$ can be obtained.