Uncertainty relation of mixed states by means of Wigner-Yanase-Dyson
information

D. Lia111email<EMAIL_ADDRESS>X. Lib, F. Wangc, H. Huangd,
X. Lie, L. C. Kwekf

a Dept of mathematical sciences, Tsinghua University, Beijing 100084 CHINA

b Department of Mathematics, University of California, Irvine, CA 92697-3875,
USA

c Insurance Department, Central University of Finance and Economics, Beijing
100081, CHINA

d Electrical Engineering and Computer Science Department

University of Michigan, Ann Arbor, MI 48109, USA

e Dept. of Computer Science, Wayne State University, Detroit, MI 48202, USA

f National Institute of Education, Nanyang Technological University, 1 Nanyang
Walk, Singapore 637616

Centre for Quantum Technologies, National University of Singapore, 3 Science
Drive 2, Singapore 117543

Institute of Advanced Studies (IAS), Nanyang Technological University, 60
Nanyang View Singapore 639673

Abstract

The variance of an observable in a quantum state is usually used to describe
Heisenberg uncertainty relation. For mixed states, the variance includes
quantum uncertainty and classical uncertainty. By means of the skew
information and the decomposition of the variance, a stronger uncertainty
relation was presented by Luo in [Phys. Rev. A 72, 042110 (2005)]. In this
paper, by using Wigner-Yanase-Dyson information which is a generalization of
the skew information, we propose a general uncertainty relation of mixed
states.

PACS 03.65.Ta

Keywords: Heisenberg uncertainty relation, the skew information, Dyson
information

##  1 Introduction

In quantum measurement theory, the Heisenberg uncertainty principle provides a
fundamental limit for the measurements of incompatible observables. On the
other hand, as dictated by Cramer-Rao’s lower bound, there is also an ultimate
limit for the resolution of any unbiased parameter (see for instance, [1]),
and this lower bound is given by a quantity called Fisher information. A long
time ago, Wigner demonstrated that it is more difficult to measure observables
that do not commute with some additive conserved quantity. Thus, observables
not commuting with some conserved quantity cannot be measured exactly and only
approximate measurement is possible. This trade-off in measurement forms the
basis of the well-known Wigner-Araki-Yanase theorem. In their study of quantum
measurement theory, Wigner and Yanase introduced a quantity called the skew
information. As shown in [2], the skew information is essentially a form of
Fisher information.

The skew information for a mixed state $\rho$ relative to a self-adjoint
“observable”, $A$, is defined as $I(\rho$, $A)=$
$-\frac{1}{2}\mbox{Tr}\rho^{1/2}$, $A]^{2}$. This definition was subsequently
generalized by Dyson as
$I_{\alpha}(\rho,A)=-\frac{1}{2}\mbox{Tr}([\rho^{\alpha},A][\rho^{1-\alpha},A])$,
where $0<\alpha<1$ [3]. When $\alpha=1/2$, $I_{\alpha}(\rho$, $X)$ is reduced
to the skew information. The convexity of $I_{\alpha}(\rho,A)$ was finally
resolved by Lieb[4, 5].

The von Neumann entropy of $\rho$, defined as $S(\rho)=-tr\rho\ln\rho$, has
been widely used as a measure of the uncertainty of a mixed state. This
quantity, profoundly rooted in quantum statistical mechanics, possesses
several remarkable and satisfactory properties. Like all measures, the von
Neumann entropy, together with its classical analog called the Shannon
entropy, is not always the best measure under certain contexts. In [6, 7, 2,
8], the skew information was proposed as means to unify the study of
Heisenberg uncertainty relation for mixed states.

It is well know in the standard textbooks that the Heisenberg uncertainty
relation for any two self-adjoint operators $X$ and $Y$ is given by

$V(\rho,X)V(\rho,Y)\geq\frac{1}{4}||\mbox{Tr}(\rho[X,Y]||^{2}.$ (1)

Note that $[$,$]$ is commutator, i.e. $[A$, $B]=AB-BA$ and the variance of the
observable $X$ with respect to $\rho$ is

$V(\rho,X)=\mbox{Tr}(\rho X^{2})-(\mbox{Tr}(\rho X))^{2}.$ (2)

A similar definition applies to $V(\rho,Y)$.

When $\rho$ is a mixed state, Luo showed that the variance comprises of two
terms: a quantum uncertainty term and a classical uncertainty term[6, 7]. He
separated the variance into its quantum and classical part by using the skew
information. He interpreted $I(\rho$, $X)$ as the quantum uncertainty of $X$
in $\rho$ by the Bohr complementary principle and $V(\rho,X)-I(\rho$, $X)$ as
the classical uncertainty of the mixed state. He then considered
$U(\rho,X)=\sqrt{V^{2}(\rho,X)-[V(\rho,X)-I(\rho,X)]^{2}}$ as a measure of
quantum uncertainty. Thus, he obtained the following two inequalities for the
uncertainty relation.

$I(\rho,X)J(\rho,Y)\geq\frac{1}{4}||\mbox{Tr}(\rho[X,Y]||^{2}.$ (3)
$U(\rho,X)U(\rho,Y)\geq\frac{1}{4}||\mbox{Tr}(\rho[X,Y]||^{2}.$ (4)

where $J(\rho$, $Y)=\frac{1}{2}\mbox{Tr}\\{\rho^{1/2}$, $Y_{0}\\}^{2}$, and
$Y_{0}=Y-\mbox{Tr}(\rho Y)$. The notation $\\{$ $\\}$ is the anticommutator,
i.e. $\\{A$, $B\\}=AB+BA$.

This article is organized as follows: In section 2, we discuss various
properties of the Wigner-Yanase-Dyson information. We show using a counter
example that it need not satisfy the uncertainty relation obtained from the
skew information. In section 3, we formulate an uncertainty relation for
Wigner-Yanase-Dyson information. Finally, in section 4, we reiterate our main
results. We have also provided two appendices concerning the proof of the new
uncertainty principle and additivity of the Wigner-Yanase-Dyson information.

##  2 Wigner-Yanase-Dyson information violates Heisenberg uncertainty relation

In this paper, we extend the above discussion to Wigner-Yanase-Dyson
information. The skew information proposed by Dyson can also be written as

$\displaystyle I_{\alpha}(\rho,X)$ $\displaystyle=$
$\displaystyle\mbox{Tr}(\rho X^{2})-\mbox{Tr}(\rho^{\alpha}X\rho^{1-\alpha}X)$
(5) $\displaystyle=$ $\displaystyle\mbox{Tr}(\rho
X_{0}^{2})-\mbox{Tr}(\rho^{\alpha}X_{0}\rho^{1-\alpha}X_{0})\text{, }$

where $X_{0}=X-\mbox{Tr}(\rho X)$. $I_{\alpha}(\rho,X)$ is positive from Eq.
(LABEL:q-info-2). Similarly, we define
$J_{\alpha}(\rho,Y)=\frac{1}{2}tr(\\{\rho^{\alpha}$,
$Y_{0}\\}\\{\rho^{1-\alpha}$, $Y_{0}\\})$. When $\alpha=1/2$,
$J_{\alpha}(\rho$, $Y)$ is reduced to $J(\rho$, $Y)$. As well, we can define
$J_{\alpha}(\rho,X)$, $J_{\alpha}(\rho,A)$, and $J_{\alpha}(\rho,B)$. By
calculating,

$\displaystyle J_{\alpha}(\rho,Y)=$ $\displaystyle\mbox{Tr}(\rho
Y_{0}^{2})+\mbox{Tr}(\rho^{\alpha}Y_{0}\rho^{1-\alpha}Y_{0})=$
$\displaystyle\mbox{Tr}(\rho
Y^{2})+\mbox{Tr}(\rho^{\alpha}Y\rho^{1-\alpha}Y)-2(\mbox{Tr}\rho Y)^{2}.$ (6)

$J_{\alpha}(\rho$, $Y)$ is also positive from Eq. (A9) in this paper.

Adopting the Luo’s interpretations, by the following properties of Wigner-
Yanase-Dyson information we interpret $I_{\alpha}(\rho,X)$ as quantum
uncertainty of $X$ in $\rho$, $V(\rho,X)-I_{\alpha}(\rho$, $X)$ as the
classical mixing uncertainty, and
$U_{\alpha}(\rho,X)=\sqrt{V^{2}(\rho,X)-[V(\rho,X)-I_{\alpha}(\rho,X)]^{2}}$
as a measure of quantum uncertainty. Lieb studied the properties of Wigner-
Yanase-Dyson information in [4]. Wigner-Yanase-Dyson information satisfies the
following requirements.

(1). Wigner-Yanase-Dyson conjecture about the convexity of
$I_{\alpha}(\rho,X)$ with respect to $\rho$ was proved by Lieb [4].

(2). Wigner-Yanase-Dyson information $I_{\alpha}(\rho,X)$ is additive under
the following sense (See [2] and [4]). Let $\rho_{1}$ and $\rho_{2}$ be two
density operators of two subsystems, and $A_{1}$ (resp. $A_{2}$) be a self-
adjoint operator on $H^{1}$ (resp. $H^{2}$). $I_{\alpha}(\rho,X)$ is additive
if $I_{\alpha}(\rho_{1}\otimes\rho_{2}$, $A_{1}\otimes I_{2}+I_{1}\otimes
A_{2})=I_{\alpha}(\rho_{1}$, $A_{1})+I_{\alpha}(\rho_{2}$, $A_{2})$, where
$I_{1}$ and $I_{2}$ are the identity operators for the first and second
systems, respectively. For the proof see Appendix B.

(3). $J_{\alpha}(\rho$, $Y)$ is also additive under the above sense. For the
proof see Appendix B.

(4). However, Hansen showed that Wigner-Yanase-Dyson information is not
subadditive [11]. For the definition of subadditivity see [4] and [11].

(5). $J_{\alpha}(\rho$, $Y)$ is concave with respect to $\rho$. This is
because $tr(\rho Y_{0}^{2})$ is linear operator with respect to $\rho$ and
$tr(\rho^{\alpha}Y_{0}\rho^{1-\alpha}Y_{0})$ is concave with respect to
$\rho$.

(6). When $\rho$ is pure, $V(\rho,X)=I_{\alpha}(\rho$, $X)$. Thus, Wigner-
Yahase-Dyson information reduces to the variance. That is, the variance
$V(\rho,X)$ does not include the classical mixing uncertainty because of no
mixing. In other words, the variance only includes the quantum uncertainty of
$X$ in $\rho$. The case in which $\alpha=1/2$ was discussed in [7].

The above fact can be argued as follows. When $\rho$ is pure,
$tr(\rho^{\alpha}X_{0}\rho^{1-\alpha}X_{0})=$ $(tr(\rho X_{0}))^{2}=0$. Thus,
$I_{\alpha}(\rho,X)=tr(\rho X_{0}^{2})=V(\rho,X)$.

(7). When $\rho$ is a mixed state, $V(\rho,X)\geq I_{\alpha}(\rho$, $X)$. This
is because $tr(\rho^{\alpha}X\rho^{1-\alpha}X)$ $=$
$tr((\rho^{\alpha/2}X\rho^{(1-\alpha)/2})$
$(\rho^{\alpha/2}X\rho^{(1-\alpha)/2})^{\dagger})\geq 0$. Also, see Eq. (A3)
in this paper. The case in which $\alpha=1/2$ was discussed in [7].

(8). When $\rho$ and $A$ commute, according to the discussion for the skew
information in [6, 8], the quantum uncertainty should vanish and thus, the
variance only includes the classical uncertainty. We can argue that the above
conclusion is also true for Wigner-Yanase-Dyson information. When $\rho$ and
$A$ commute, it is well known that $\rho$ and $A$ have the same orthonormal
eigenvector basis [9]. Hence, $\rho^{\alpha}$ and $A$ also commute. By the
definition in Eq. (5), Wigner-Yanase-Dyson information $I_{\alpha}(\rho,X)$
vanishes.

However, $I_{\alpha}(\rho,X)$ and $J_{\alpha}(\rho$, $Y)$ do not satisfy Eq.
(3). We give the following counter example for Eq. (3).

Let $n=2$, $\alpha=1/4$, and $\rho$ have the eigenvalues $\lambda_{1}=1/4$ and
$\lambda_{2}=3/4$. Since $A$ and $B$ are self-adjoint, then we write
$A=\left(\begin{tabular}[]{ll}$x$&$u+iv$\\\ $u-iv$&$y$\end{tabular}\right)$,
$B=\left(\begin{tabular}[]{ll}$a$&$c+di$\\\ $c-di$&$b$\end{tabular}\right)$.
In this example, $u=4$, $v=2$, $a=b=0$, $c=1$, and $d=-5$. By calculating
$I_{\alpha}(\rho,A)$ in Eq. (LABEL:q-info-2) and $J_{\alpha}(\rho,B)$ in Eq.
(A8),
$I_{\alpha}(\rho,A)J_{\alpha}(\rho,B)=[1-(\lambda_{1}^{\alpha}\lambda_{2}^{1-\alpha}+\lambda_{2}^{\alpha}\lambda_{1}^{1-\alpha})^{2}](u^{2}+v^{2})(c^{2}+d^{2})=99.83$.
By calculating $\mbox{Tr}(\rho[A,B]$ in Eq. (A11),
$\frac{1}{4}|\mbox{Tr}(\rho[A,B]|^{2}=(\lambda_{1}-\lambda_{2})^{2}(cv-
du)^{2}=121$. Hence, it violates Eq. (3). It implies that the bound on the
right side of the inequality in Eq. (3) is too large in this example. We need
to get the appropriate lower bound for Wigner-Yanase-Dyson information, i.e.,
we need to modify the term on RHS of the inequality.

##  3 The general uncertainty relation

We replace $\mbox{Tr}(\rho[X,Y]$ with $l_{\alpha}(\rho$, $X$, $Y)$ which is
defined as follows:

$l_{\alpha}(\rho,X,Y)=\mbox{Tr}(\rho[X,Y])-\mbox{Tr}\rho^{\left|2\alpha-1\right|}[X,Y].$
(7)

When $\alpha=1/2$, $l_{\alpha}(\rho$, $X$, $Y)$ reduces to $\mbox{Tr}(\rho[X$,
$Y])$. In [6], Luo defined $k=\mathrm{i}[\rho^{1/2}$,
$X_{0}]t+\\{\rho^{1/2},Y_{0}\\}$, where $t\in R$ and $\mathrm{i}$ is an
imaginary number. From $\mbox{Tr}(kk^{\dagger})\geq 0$, by expanding
$\mbox{Tr}(kk^{\dagger})$, he derived $\mbox{Tr}(kk^{\dagger})=2(I[\rho$,
$X]t^{2}+\mathrm{i}(tr(\rho[X$, $Y])t+J[\rho$, $Y])\geq 0$. Since the above
inequality is true for any real $t$, Luo obtained the inequality in Eq. (3).
However, unlike his previous case, the form of $I_{\alpha}(\rho,X)$ does not
allow us to employ the trick $k=\mathrm{i}[\rho^{\alpha}$,
$X_{0}]t+\\{\rho^{\alpha},Y_{0}\\}$ nor $k=\mathrm{i}[\rho^{1-\alpha}$,
$X_{0}]t+\\{\rho^{1-\alpha},Y_{0}\\}$ to derive the uncertainty relation from
$\mbox{Tr}(kk^{\dagger})\geq 0$. The proof becomes more involved and one needs
to modify the RHS of the previous uncertainty relation.

In Appendix A, we see that if $A$ and $B$ are self-adjoint observables, then

$I_{\alpha}(\rho,A)J_{\alpha}(\rho,B)\geq\frac{1}{4}||l_{\alpha}(\rho,A,B)||^{2}\text{,
}$ (8)

and

$I_{\alpha}(\rho,B)J_{\alpha}(\rho,A)\geq\frac{1}{4}||l_{\alpha}(\rho,A,B)||^{2}\text{.
}$ (9)

If we denote $U_{\alpha}(\rho,\mathcal{O})$ as
$\sqrt{V^{2}(\rho,\mathcal{O})-[V(\rho,\mathcal{O})-I_{\alpha}(\rho,\mathcal{O})]^{2}}$,
we see that by Eq. (2) and Eq.(5) (and the analogous form for
$J_{\alpha}(\rho,\mathcal{O})$),
$U_{\alpha}(\rho,\mathcal{O})=\sqrt{I_{\alpha}(\rho,\mathcal{O})J_{\alpha}(\rho,\mathcal{O})}$,
where $\mathcal{O}$ is either the operator $A$ or $B$. Thus, we obtain our
main result from Eqs. (8) and (9),

$U_{\alpha}(\rho,A)U_{\alpha}(\rho,B)\geq\frac{1}{4}||l_{\alpha}(\rho,X,Y)||^{2}.$
(10)

For the counter example in Sec. 2, a direct calculation of Eq. (A13) yields
$\frac{1}{4}||l_{\alpha}(\rho,A,B)||^{2}=$ $8.\,\allowbreak 687\,4$.
Therefore, the inequality in Eq. (8) holds in this case.

##  4 Summary

In [6], Luo presented a refined Heisenberg uncertainty relation. In this
paper, we demonstrate some properties of Wigner-Yanase-Dyson information and
provide a counter example to show that Wigner-Yanase-Dyson information does
not in general satisfy Heisenberg uncertainty relation. We have also proposed
a new general uncertainty relation of mixed states based on Wigner-Yanase-
Dyson information. Bell-type inequalities based on the skew information have
been proposed as nonlinear entanglement witnesses [12]. We note here that
similar Bell-type inequalities with the advantage of an additional $\alpha$
parameter for fine adjustments could also be constructed from the uncertainty
principle derived from the Wigner-Yanase-Dyson information.

## Appendix A. Proof of uncertainty relation

By spectral decomposition, there exists an orthonormal basis $\\{x_{1}$,…,
$x_{n}\\}$ consisting of eigenvectors of $\rho$. Let $\lambda_{1}$, …,
$\lambda_{n}$ be the corresponding eigenvalues, where
$\lambda_{1}+...+\lambda_{n}=1$ and $\lambda_{i}\geq 0$. Thus, $\rho$ has a
spectral representation

$\rho=\lambda_{1}|x_{1}\rangle\langle
x_{1}|+....+\lambda_{n}|x_{n}\rangle\langle x_{n}|.$ (A1)

### 1\. Calculating $I_{\alpha}(\rho$, $A)$

By Eq. (A1), $\rho A^{2}=\lambda_{1}|x_{1}\rangle\langle
x_{1}|A^{2}+....+\lambda_{n}|x_{n}\rangle\langle x_{n}|A^{2}$ and

$\displaystyle\mbox{Tr}\rho A^{2}$ $\displaystyle=$
$\displaystyle\lambda_{1}\langle
x_{1}|A^{2}|x_{1}\rangle+....+\lambda_{n}\langle x_{n}|A^{2}|x_{n}\rangle$
(A2) $\displaystyle=$
$\displaystyle\lambda_{1}||A|x_{1}||^{2}+....+\lambda_{n}||A|x_{n}||^{2}.$

Moreover, since $\rho^{\alpha}A=\lambda_{1}^{\alpha}|x_{1}\rangle\langle
x_{1}|A+....+\lambda_{n}^{\alpha}|x_{n}\rangle\langle x_{n}|A$ and
$\rho^{1-\alpha}A=\lambda_{1}^{1-\alpha}|x_{1}\rangle\langle
x_{1}|A+....+\lambda_{n}^{1-\alpha}|x_{n}\rangle\langle x_{n}|A$, we have,
$\rho^{\alpha}A\rho^{1-\alpha}A=\sum_{i,j=1}\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}|x_{i}\rangle\langle
x_{i}|A|x_{j}\rangle\langle x_{j}|A$. Thus

$\displaystyle\mbox{Tr}\rho^{\alpha}A\rho^{1-\alpha}A$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}\langle
x_{i}|A|x_{j}\rangle\langle x_{j}|A|x_{i}\rangle$ (A3) $\displaystyle=$
$\displaystyle\sum_{i,j=1}\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}||\langle
x_{i}|A|x_{j}\rangle||^{2}.$

From Eqs. (5), (A2) and (A3),

$I_{\alpha}(\rho,A)=\sum_{i=1}\lambda_{i}||A|x_{i}||^{2}-\sum_{i,j=1}\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}||\langle
x_{i}|A|x_{j}\rangle||^{2}.$ (A4)

Let $A=\\{A_{ij}\\}$ (resp. $B=\\{B_{ij}\\}$) be the matrix representation of
the operator $A$ (resp. $B$) corresponding to the orthonormal basis
$\\{x_{1}$,…, $x_{n}\\}$. Then $\langle x_{i}|A|x_{j}\rangle=A_{ij}$, and

$\displaystyle I_{\alpha}(\rho,A)$ $\displaystyle=$ $\displaystyle\sum_{i\neq
j}(\lambda_{i}-\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha})\left|\left|A_{ij}\right|\right|^{2}$
$\displaystyle=$
$\displaystyle\sum_{i<j}(\lambda_{i}+\lambda_{j}-\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}-\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})\left|\left|A_{ij}\right|\right|^{2}\text{.}$

### 2\. Calculating $J_{\alpha}(\rho,B)$

Similarly, from Eqs. (6) and (A1), we can obtain

$\displaystyle J_{\alpha}(\rho,B)$ $\displaystyle=$
$\displaystyle\sum_{i=1}\lambda_{i}||B|x_{i}||^{2}+\sum_{i,j=1}\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}||\langle
x_{i}|B|x_{j}\rangle||^{2}$ (A6) $\displaystyle-2(\sum\lambda_{i}\langle
x_{i}|B|x_{i}\rangle)^{2}.$

Let $\langle x_{i}|B|x_{j}\rangle=B_{ij}$. Then, from Eq. (A6),

$\displaystyle J_{\alpha}(\rho,B)$ $\displaystyle=$ $\displaystyle
2\sum_{i=1}\lambda_{i}\left|B_{ii}\right|^{2}-2(\sum_{i=1}\lambda_{i}B_{ii})^{2}$
(A7) $\displaystyle+\sum_{i\neq
j}(\lambda_{i}+\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha})\left|\left|B_{ij}\right|\right|^{2}\text{.}$

By simplifying,

$\displaystyle J_{\alpha}(\rho,B)$ $\displaystyle=$ $\displaystyle
2\sum_{i=1}\lambda_{i}\left|B_{ii}\right|^{2}-2(\sum_{i=1}\lambda_{i}B_{ii})^{2}$
(A8)
$\displaystyle+\sum_{i<j}(\lambda_{i}+\lambda_{j}+\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})\left|\left|B_{ij}\right|\right|^{2}\text{.}$

Since $x^{2}$ is convex,
$(\sum_{i=1}\lambda_{i}B_{ii})^{2}\leq\sum_{i=1}\lambda_{i}\left|B_{ii}\right|^{2}$.
So from Eq. (A8),

$J_{\alpha}(\rho,B)\geq\sum_{i<j}(\lambda_{i}+\lambda_{j}+\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})\left|\left|B_{ij}\right|\right|^{2}\text{.}$
(A9)

### 3\. Calculating $l_{\alpha}(\rho$, $A$, $B)$

First we calculate $\mbox{Tr}(\rho[A$, $B])$. By Eq. (A1),
$\rho[A,B]=\lambda_{1}|x_{1}\rangle\langle x_{1}|[A$,
$B]+....+\lambda_{n}|x_{n}\rangle\langle x_{n}|[A$, $B]$ and
$\mbox{Tr}(\rho[A$, $B])=\lambda_{1}\langle x_{1}|[A$,
$B]|x_{1}\rangle+....+\lambda_{n}\langle x_{n}|[A$, $B]|x_{n}\rangle$. It is
well known that $Re\langle x_{i}|[A$, $B]|x_{i}\rangle=0$ and $\langle
x_{i}|[A$, $B]|x_{i}\rangle=\mathrm{i}(2Im\langle x_{i}|AB|x_{i}\rangle)$,
where $\mathrm{i}$ is an imaginary number. Consequently, $\mbox{Tr}(\rho[A$,
$B])=2\mathrm{i}(\lambda_{1}Im\langle
x_{1}|AB|x_{1}\rangle+....+\lambda_{n}Im\langle x_{n}|AB|x_{n}\rangle)$.
Therefore we obtain

$\displaystyle\mbox{Tr}(\rho[A,B])$ $\displaystyle=$ $\displaystyle
2\mathrm{i}Im(\lambda_{1}\langle x_{1}|AB|x_{1}\rangle+...+\lambda_{n}\langle
x_{n}|AB|x_{n}\rangle)$ (A10) $\displaystyle=$ $\displaystyle
2\mathrm{i}Im\sum_{j\neq i}\lambda_{i}A_{ij}B_{ji}.$

Note that in Eq. (A10)$\ A_{ii}$ and $B_{ii}$ are real because $A$ and $B$ are
self-adjoint. Since $A_{ij}B_{ji}=(A_{ji}B_{ij})^{\ast}$, $\Im\sum_{j\neq
i}\lambda_{i}A_{ij}B_{ji}=Im\sum_{i<j}(\lambda_{i}-\lambda_{j})A_{ij}B_{ji}$.
Thus, by simplifying,

$\mbox{Tr}(\rho[A,B])=2\mathrm{i}Im\sum_{i<j}(\lambda_{i}-\lambda_{j})A_{ij}B_{ji}.$
(A11)

Moreover,

$\mbox{Tr}\rho^{\left|2\alpha-1\right|}[A,B]=2\mathrm{i}Im\sum_{i<j}(\lambda_{i}^{\left|2\alpha-1\right|}-\lambda_{j}^{\left|2\alpha-1\right|})A_{ij}B_{ji}.$
(A12)

Hence, from Eqs. (7), (A11) and (A12),

$l_{\alpha}(\rho,A,B)=2\mathrm{i}\sum_{i<j}(\lambda_{i}-\lambda_{j}-(\lambda_{i}^{\left|2\alpha-1\right|}-\lambda_{j}^{\left|2\alpha-1\right|}))Im(A_{ij}B_{ji}).$
(A13)

### 4\. The proof of the uncertainty relation

From Eqs. (LABEL:q-info-2), (A9) and (A13), for Eq. (8) we need to show

$\displaystyle[\sum_{i<j}(\lambda_{i}+\lambda_{j}-\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}-\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})\left|\left|A_{ij}\right|\right|^{2}][\sum_{i<j}(\lambda_{i}+\lambda_{j}+\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})\left|\left|B_{ij}\right|\right|^{2}]$
(A14) $\displaystyle\geq$
$\displaystyle\\{\sum_{i<j}[\lambda_{i}-\lambda_{j}-(\lambda_{i}^{\left|2\alpha-1\right|}-\lambda_{j}^{\left|2\alpha-1\right|})]Im(A_{ij}B_{ji})\\}^{2}.$

It is easy to know
$[Im(A_{ij}B_{ji})]^{2}\leq\left|\left|A_{ij}\right|\right|^{2}\left|\left|B_{ij}\right|\right|^{2}$.
Note that
$\lambda_{i}+\lambda_{j}-\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}-\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha}=(\lambda_{i}^{\alpha}-\lambda_{j}^{\alpha})(\lambda_{i}^{1-\alpha}-\lambda_{j}^{1-\alpha})\geq
0$. By the Cauchy-Schwartz inequality, the LHS of the inequality in Eq. (A14)
$\geq\\{\sum[(\lambda_{i}+\lambda_{j})^{2}-(\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})^{2}]^{1/2}Im(A_{ij}B_{ji})\\}^{2}$.
Finally, what needs to be shown is

$\displaystyle(\lambda_{i}+\lambda_{j})^{2}-(\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})^{2}$
(A15) $\displaystyle\geq$
$\displaystyle|(\lambda_{i}-\lambda_{j})-(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})|^{2}\text{.}$

It is easy to see that

$\displaystyle(\lambda_{i}+\lambda_{j})^{2}-(\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}+\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})^{2}$
$\displaystyle=$
$\displaystyle(\lambda_{i}-\lambda_{j})^{2}-(\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}-\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})^{2}.$

When $\alpha\geq 1/2$,

$\displaystyle(\lambda_{i}-\lambda_{j})^{2}-(\lambda_{i}^{\alpha}\lambda_{j}^{1-\alpha}-\lambda_{i}^{1-\alpha}\lambda_{j}^{\alpha})^{2}$
$\displaystyle=$
$\displaystyle(\lambda_{i}-\lambda_{j})^{2}-\lambda_{i}^{2(1-\alpha)}\lambda_{j}^{2(1-\alpha)}(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})^{2}$
$\displaystyle\geq$
$\displaystyle(\lambda_{i}-\lambda_{j})^{2}-(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})^{2}$
$\displaystyle=$
$\displaystyle|(\lambda_{i}-\lambda_{j})-(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})|$
$\displaystyle\times\left|(\lambda_{i}-\lambda_{j})+(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})\right|$
$\displaystyle\geq$
$\displaystyle|(\lambda_{i}-\lambda_{j})-(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})|^{2}.$

Note that the last inequality holds because $(\lambda_{i}-\lambda_{j})$ and
$(\lambda_{i}^{2\alpha-1}-\lambda_{j}^{2\alpha-1})$ have the same sign. Also,
when $0<\alpha\leq 1/2$, we can prove the inequality in Eq. (A15) as follows:
Let $\beta=1-\alpha$ with $1/2\leq\beta<1$. Replacing $\alpha$ in Eq. (A15)
with $1-$ $\beta$, we obtain
$(\lambda_{i}+\lambda_{j})^{2}-(\lambda_{i}^{1-\beta}\lambda_{j}^{\beta}+\lambda_{i}^{\beta}\lambda_{j}^{1-\beta})^{2}\geq|(\lambda_{i}-\lambda_{j})-(\lambda_{i}^{2\beta-1}-\lambda_{j}^{2\beta-1})|^{2}$.
This ends the proof.

## Appendix B. Additivity

The quantity $J_{\alpha}(\rho$, $B)$ is additive in the following sense:
$J_{\alpha}(\rho_{1}\otimes\rho_{2}$, $B_{1}\otimes I_{2}+I_{1}\otimes
B_{2})=J_{\alpha}(\rho_{1}$, $B_{1})+J_{\alpha}(\rho_{2}$, $B_{2})$. Using the
notation in [4], the proof proceeds by letting
$\rho_{12}=\rho_{1}\otimes\rho_{2}$ and $L=B_{1}\otimes I_{2}+I_{1}\otimes
B_{2}$. Setting
$\rho_{12}^{\alpha}=\rho_{1}^{\alpha}\otimes\rho_{2}^{\alpha}$, we have

$\rho_{12}^{\alpha}L\rho_{12}^{1-\alpha}L$
$=\rho_{1}^{\alpha}B_{1}\rho_{1}^{1-\alpha}B_{1}\otimes\rho_{2}+\rho_{1}^{\alpha}B_{1}\rho_{1}^{1-\alpha}\otimes\rho_{2}B_{2}$
$+\rho_{1}B_{1}\otimes\rho_{2}^{\alpha}B_{2}\rho_{2}^{1-\alpha}+\rho\otimes\rho_{2}^{\alpha}B_{2}\rho_{2}^{1-\alpha}B_{2}$,
and

$\displaystyle\mbox{Tr}(\rho_{12}^{\alpha}L\rho_{12}^{1-\alpha}L)$ (B1)
$\displaystyle=$
$\displaystyle\mbox{Tr}(\rho_{1}^{\alpha}B_{1}\rho_{1}^{1-\alpha}B_{1})+2\mbox{Tr}(\rho_{1}B_{1})\mbox{Tr}(\rho_{2}B_{2})+\mbox{Tr}(\rho_{2}^{\alpha}B_{2}\rho_{2}^{1-\alpha}B_{2}).$

Similarly,

$\mbox{Tr}(\rho_{12}^{\alpha}L^{2})=\mbox{Tr}(\rho_{1}B_{1}^{2})+2\mbox{Tr}(\rho_{1}B_{1})\mbox{Tr}(\rho_{2}B_{2})+\mbox{Tr}(\rho_{2}B_{2}^{2}).$
(B2)

From the above Eqs. (B1) and (B2), we can derive $I_{\alpha}(\rho$, $B)$ is
additive.

Similarly,

$\mbox{Tr}(\rho_{12}^{\alpha}L)=\mbox{Tr}(\rho_{1}B_{1})+\mbox{Tr}(\rho_{2}B_{2}).$
(B3)

By Eqs. (B1), (B2), and (B3), and the definition of $J_{\alpha}(\rho$, $B)$ in
Eq. (6), we can conclude that $J_{\alpha}(\rho$, $B)$ is additive.

Acknowledgments

The first author wants to thank Prof. Jinwen Chen for his helpful discussion
about the inequality in Eq. (A15) and Mr Qin Zhang for the discussion about
the idea for $\mbox{Tr}(KK^{\prime})\geq 0$. The paper was supported by
NSFC(Grants No.10875061, 60673034. KLC would like to acknowledge financial
support by the National Research Foundation & Ministry of Education,
Singapore, for his visit and collaboration at Tsinghua University.

## References

  * [1] Fisher information has been discussed extensively in literature on statistical estimation theory. A particularly insightful survey of the connection between Fisher information and the Heisenberg uncertainty principle can be found in arXiv: quant-ph//0309184
  * [2] S. Luo, Phys. Rev. Lett., 91, 180403 (2003).
  * [3] E.P. Wigner and M. M. Yanase, Proc. Nat. Acad. Sci. U.S.A. 49, 910-918 (1963).
  * [4] E.H. Lieb, Adv. Math. 11, 267 (1973).
  * [5] E.H. Lieb and M.M. Ruskai, Phys. Rev. Lett. 30, 434 (1973).
  * [6] S. Luo, Phys. Rev. A 72, 042110 (2005).
  * [7] S. Luo, Phys. Rev. A 73, 022324 (2006).
  * [8] S. Luo, Theor. Math. Phys. 143, 681 (2005).
  * [9] M. Hirvensalo, Quantum Computing, Springer-Verlag, Berlin, (2001).
  * [10] M. A. Nielsen and I. L. Chuang, see p.89, Quantum Computation and Quantum Information, Cambridge Univ. Press, Cambridge (2000).
  * [11] F. Hansen, J. Stat. Phys., 126, 643 (2007).
  * [12] Z. Chen, Phys. Rev. A 71, 052302 (2005).


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Thermal melting of discrete time crystals: a dynamical phase transition
induced by thermal fluctuations

Mingxi Yue  Wilczek Quantum Center and Key Laboratory of Artificial Structures
and Quantum Control, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China  Xiaoqin Yang  Wilczek Quantum Center and
Key Laboratory of Artificial Structures and Quantum Control, School of Physics
and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Zi Cai
<EMAIL_ADDRESS>Wilczek Quantum Center and Key Laboratory of Artificial
Structures and Quantum Control, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Shanghai Research Center for Quantum
Sciences, Shanghai 201315, China

###### Abstract

The stability of a discrete time crystal against thermal fluctuations has been
studied numerically by solving a stochastic Landau-Lifshitz-Gilbert equation
of a periodically-driven classical system composed of interacting spins, each
of which couples to a thermal bath. It is shown that in the thermodynamic
limit, even though the long-range temporal crystalline order is stable at low
temperature, it is melting above a critical temperature, at which the system
experiences a non-equilibrium phase transition. The critical behaviors of the
continuous phase transition have been systematically investigated, and it is
shown that despite the genuine non-equilibrium feature of such a periodically
driven system, its critical properties fall into the 3D Ising universality
class with a dynamical exponent ($z=2$) identical to that in the critical
dynamics of kinetic Ising model without driving.

Introduction – Spontaneous symmetry breakings (SSB) and universality classes
are among the most fundamental concepts in modern physics: the former is used
to characterize different phases of matter while the latter is to categorize
the transitions between them. Compared to equilibrium systems, the physics in
out-of-equilibrium systems is much richer but less known in general. The
evolutionary processes give rise to increasing richness of the paradigm of
spontaneous symmetry breaking and universality classes, which could take place
not only in space but also in time, the latter opens up new opportunities to
explore novel states of matter (the time crystal(TC)Wilczek (2012) for
instance) and dynamical universality classes ( the critical dynamical
phenomenaHohenberg and Halperin (1977); Taeuber (2014) or the Kardar-Parisi-
Zhang (KPZ) universality classKardar et al. (1986) for instance) beyond the
scope of equilibrium physics.

As a prototypical example of spontaneous symmetry breaking in time domain, the
time crystal phaseWilczek (2012), in its different formsShapere and Wilczek
(2012); Li et al. (2012); Wilczek (2013); Sacha (2015); Else et al. (2016);
Khemani et al. (2016); Yao et al. (2017); Russomanno et al. (2017); Gong et
al. (2018); Huang et al. (2018); Iemini et al. (2018); Das et al. (2018); Zhu
et al. (2019); Liao et al. (2019); Kozin and Kyriienko (2019); Cai et al.
(2020); Chinzei and Ikeda (2020); Lyu et al. (2020); Choudhury and Liu (2021),
has attracted considerable interests in past decades. In spite of being proven
absent in thermal equilibriumBruno (2013); Watanabe and Oshikawa (2015), such
an intriguing phase has been observed in periodically driven non-equilibrium
settingsChoi et al. (2017); Zhang et al. (2017); Autti et al. (2018); Träger
et al. (2021); Keßler et al. (2021); Mi et al. (2021); Frey and Rachel (2021).
These systems exhibit oscillations with period doubling with respect to that
of the external driving, thus spontaneously break the discrete time
translational symmetry(TTS) from symmetry group $\mathbb{Z}$ to $2\mathbb{Z}$.
As a phase of matter, a profound question is its stability against
perturbations, especially the thermal fluctuations inevitable in almost all
the realistic experimental setups. From a theoretical point of view, this
problem can be considered as a non-equilibrium analogue of thermal melting of
the crystalline order, and is of fundamental interest due to its relevance to
broader questions of the robustness of the temporal order against
fluctuations, as well as the dynamical universality classes in non-equilibrium
matter.

Recently, Yao et al explored this question by studying a classical Hamiltonian
dynamics coupled to a finite-temperature bath, and found an activated discrete
time crystal(DTC) whose crystalline order survives to long, but not infinite
times at low temperature Yao et al. (2020). A question naturally raised is
whether there exists DTC phases with true long-range crystalline order
persisting forever Zhuang et al. (2021)? Is it possible for the thermal
fluctuation to melt such a dynamical crystalline order, just as it does for
conventional crystals? If so, how to characterize such a non-equilibrium phase
transition and what’s the corresponding universality class?

In this paper, we attempt to answer these questions, focusing for simplicity
on a three-dimensional (3D) classical periodically-driven interacting spin
model, which could exhibit period doubling with respect to the driving (DTC
phase) at zero temperature. The thermal fluctuations are introduced by
coupling each spin to a heat bath, which provides both dissipation and noise.
The problem is approached by solving a stochastic equation of motion (EOM),
where the dynamics of each spin is governed by a Landau-Lifshitz-Gilbert (LLG)
equationLakshmanan (2010) augmented by a random thermal forceBrown (1963). It
has been shown that at high temperature, the DTC order parameter decays
exponentially with time, indicating a finite life time of the DTC phase as in
the 1D case (activated DTC). With decreasing temperature, the non-equilibrium
system experiences a dynamical phase transition from an activated DTC to a DTC
phase with true long-range crystalline order. It is interesting to show that
despite the genuine non-equilibrium feature of such periodically driven
systems, its critical properties are characterized by an Ising universality
class with a dynamical exponent ($z=2$) identical with the value in the
critical dynamics of undriven systems.

Figure 1: (Color online) Zero temperature phase diagram of the long-time
behavior of the classical EOM.(7) with $h_{z}=\lambda=J$ and $\mathcal{D}=0$.

Model and method – The time crystal phases in classical many-body systems have
attracted considerable interest recentlyLupo and Weber (2019); Yao et al.
(2020); Gambetta et al. (2019); Khasseh et al. (2019); Hurtado-Gutiérrez et
al. (2020); Ye et al. (2021); Pizzi et al. (2021). We start with a 3D
classical spin model (transverse Ising model) with a time-dependent coupling
strength. The system Hamiltonian reads:

$H_{s}(t)=-V(t)\sum_{\langle
ij\rangle}s_{i}^{x}s_{j}^{x}+\sum_{i}h_{z}s_{i}^{z}$ (1)

where the dynamical variable $\mathbf{s}_{i}$ is a classical vector with fixed
length $|\mathbf{s}_{i}|=1$. The summation $\langle ij\rangle$ is over the
nearest-neighboring (NN) sites of site $i$ in the cubic lattice with length
$L$. The strength of the interaction $V(t)$ periodically oscillates as
$V(t)=V_{0}\cos\omega_{0}t$. $h_{z}$ is the strength of a uniform magnetic
field along z-direction. Throughout this paper we fix $h_{z}=J$. In the
absence of thermal bath, the dynamics of each spin can be described by the
EOM: $\dot{\mathbf{s}}_{i}=\mathbf{h}^{0}_{i}\times\mathbf{s}_{i}$, where the
effective magnetic field $\mathbf{h}^{0}_{i}=[V(t)\bar{s}_{i}^{x},0,h_{z}]$
with $\bar{s}_{i}^{x}=\frac{1}{2}\sum_{\langle j\rangle}s_{j}^{x}$ (the
summation is over the NN sites of site i).

Even though our driven spin system is a genuine non-equilibrium system whose
temperature is ill-defined, usually the degrees freedom of bath are much
larger than those of the system thus the back action of the system to the bath
can be neglected, which allows us to consider a thermal bath with a well-
defined temperature (T) that doesn’t change during the dynamical process. The
effect of the thermal bath can be modeled by methods familiar in the theory of
the Brownian motion and other stochastic processes, where the EOM of each spin
is described by a stochastic LLG equationLakshmanan (2010) as:

$\dot{\mathbf{s}}_{i}=\mathbf{h}_{i}\times\mathbf{s}_{i}-\lambda\mathbf{s}_{i}\times(\mathbf{s}_{i}\times\mathbf{h}_{i})$
(2)

where $\lambda$ is the strength of the friction provided by the local thermal
bath, which is fixed as $\lambda=J$ throughout this paper.
$\mathbf{h}_{i}=\mathbf{h}^{0}_{i}+\bm{\xi}_{i}(t)$ is the effective magnetic
field, where $\bm{\xi}_{i}(t)$ is a three-dimensional(3D) stochastic magnetic
field representing the thermal noise, which can be approximated as white noise
if assuming that the typical time scale of the bath is much shorter than that
of the system. We further assume the local bath around each spin is
independent with each other thus the stochastic variables satisfy:

$\displaystyle\langle\xi_{i}^{\alpha}(t)\rangle_{\bm{\xi}}$ $\displaystyle=$
$\displaystyle 0$ (3)
$\displaystyle\langle\xi_{i}^{\alpha}(t)\xi_{j}^{\beta}(t^{\prime})\rangle_{\bm{\xi}}$
$\displaystyle=$
$\displaystyle\mathcal{D}^{2}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$
(4)

where $\alpha,\beta=x,y,z$, $\mathcal{D}$ is the strength of the noise, and
the average $\langle\rangle_{\bm{\xi}}$ is over all the noise trajectories. If
the bath is in thermal equilibrium with temperature T, the fluctuation-
dissipation theorem (FDT) indicates that the strengths of the friction and
noise satisfy

$\mathcal{D}^{2}=2T\lambda.$ (5)

Numerically, we adopt Stratonovich’s formula of the stochastic differential
Eq.(7), and solve it using the standard Heun methodAment et al. (2016) with
the time step of $\Delta t=10^{-3}$, the convergence of which has been checked
numerically (see the Supplmentary Material(SM)Sup ). The system size in our
simulation ranges from $L=8$ to $L=28$, which enable us to systematically
analyze the finite-size effect. The ensemble average over the noise
trajectories can be performed by directly sampling over $\mathcal{N}$ sets of
noise realizations. ($\mathcal{N}$ ranges from $4\times 10^{3}$ to $10^{5}$
depending on the simulated system size). In our simulation, we calculate the
evolution of the average magnetization along x-direction
$M(t)=\langle\frac{1}{L^{3}}\sum_{i}s_{i}^{x}(t)\rangle_{\bm{\xi}}$ to
characterize various dynamical behaviors. We choose the initial state as the
ground state of Hamiltonian.(1) with $t=0$, since we start from a spatially
uniform initial state, $M(t)$ is proportional to the auto-correlation function
$C(t)=\langle\frac{1}{L^{2}}\sum_{i}s_{i}^{x}(0)s_{i}^{x}(t)\rangle_{\bm{\xi}}$,
which characterizes the memory effect of the initial state information.
However, we find that the long-time behavior doesn’t change even if we start
from a non-uniform random initial state, that different spins will finally
synchronize with each otherSup .

Zero temperature phase diagram: – It is shown that even at zero temperature
($\mathcal{D}=0$), the classical EOM.(7) could exhibit rich long-time
dynamical behaviors, which can be found in the $V_{0}-\omega_{0}$ phase
diagram as shown in Fig.1. For sufficiently large $V_{0}$ and small
$\omega_{0}$, one can find a $DTC$-2 phase, where the ferromagnetic order
parameter $M(t)$ oscillate with a period twice of that of the driving
($\frac{2\pi}{\omega_{0}}$). Besides that, one can also find other DTC phases
whose periods are other integer multiples of $\frac{2\pi}{\omega_{0}}$
($DTC$-3 and $DTC$-4 phases for instance). If both $V_{0}$ and $\omega_{0}$
are large enough, one can find a synchronization between $M(t)$ ad $V(t)$,
both of which oscillates in the same period(synchronized phase). Besides these
periodic oscillations, one can also find other dynamical phases that $M(t)$
could decay to zero (paramagnetic phase), or oscillate chaotically (chaos).

Figure 2: (Color online). (a) The dynamics of the ferromagnetic order
parameter $M(t)$ at temperatures above ($\mathcal{D}=0.29J$) and below
($\mathcal{D}=0.26J$) the critical temperature $T_{c}=\frac{D_{c}^{2}}{2J}$
with system size $L=20$. The inset is the envelope (the peaks) of the
oscillations in Fig.2 (a) in a semi-log plot. (b) Comparison of $M(t)$ in the
evolutions under a single noise trajectory with $\mathcal{D}=0.26J$ and
$\mathcal{D}=0.29J$. A thermal activated $\pi$-phase shift can be found in the
latter. (c) Qualitatively different long-time behaviors of $M(t)$ at
temperatures above, below and at the critical temperature. (d)The dynamics of
$M(t)$ at the critical point for systems with various system size $L$. The
inset is the finite-size scaling of the exponents $\Delta(L)$ of late time
exponential dynamics at $D=D_{c}$. (e) The renormalized relaxation time as a
function of $\mathcal{D}$ for systems with different system size $L$. The
insets are the data collapse for (e) using the critical exponents $\nu=0.64$
and $z=2$. (f) At the critical temperature, the correlation length grows in
time as $l_{c}(t)\sim t^{\frac{1}{z}}$. We fixed $V_{0}=5J$, $\omega_{0}=2\pi
J$ and $h=\lambda=J$. $L=20$ for (a) and (b), and $L=28$ for (c) and (f).
$D=D_{c}=0.2782J$ for (d)-(f).

The fate of DTC phase at various temperatures: – In the following, we will
focus on the $DTC$-2 phase by fixing $V_{0}=5J$, $\omega_{0}=2\pi J$ and
change the temperature. The noises, in general, will suppress the long-range
order. As shown in Fig.2 (a), $M(t)$ exhibits a damped oscillation whose
amplitude decays exponentially with time ($\sim e^{-\Delta t}$ as shown in the
inset) at a relatively high temperature (e.g $\mathcal{D}=0.29J$ corresponding
to $T=0.042J$) . Such an exponential decay can be understood as a consequence
of the “kinks” activated by thermal fluctuations: the noise can induce
tunnelings from one of the degenerate DTC phases to the other by a $\pi$-phase
shiftingYang and Cai (2021), as shown in Fig.2 (b) for the time evolution of
M(t) under a single noise trajectory. Such topological “excitations”, no
matter how rare they are, will inevitably destroy the long-range crystalline
order along the temporal direction via a similar mechanism responsible for the
absence of spatial long-range order in 1D thermal equilibrium systems. At low
temperature (e.g $\mathcal{D}=0.26J$ corresponding to $T=0.038J$), the thermal
fluctuation is not strong enough to induce a $\pi$-phase shifting, but can
only modulate the amplitude of the oscillation, thus we find a persistent
oscillation of $M(t)$ as shown in Fig.2 (a), which indicates a DTC phase with
true long-range order in time domain.

The critical behavior of the non-equilibrium phase transition– As the
temperature is lowered, the exponent of the exponential decay $\Delta$
decreases and finally vanishes when the temperature drops below a critical
temperature $T_{c}=\frac{\mathcal{D}_{c}^{2}}{2J}$ with $D_{c}=0.2782J$, which
indicates an absence of the exponential decay in the thermodynamic limit. For
a sufficiently large system, it is shown in Fig.2 (c) that at the critical
temperature ($T=T_{c}$), the amplitude of $M(t)$ decays algebraically with
time ($\sim t^{-\alpha_{c}}$ with $\alpha_{c}=0.33(1)$), qualitatively
different from either the exponential decay at the high temperature or the
persistent oscillation at low temperature. Such a critical exponent
$\alpha_{c}=0.33$ is different from the value observed recently
($\alpha_{c}=0.5$) in a mean-field analysis of the critical properties of a
dissipative DTC phaseChinzei and Ikeda (2021).

For a finite system, the dynamics of $M(t)$ at the critical temperature can be
separated into three regimes: a transient regime that depends on the initial
state, a long-time exponential regime dominated by the finite size effect, and
an intermediate universal algebraic regime sandwiched between them, as shown
in Fig.2 (d). The finite size effect leads to an exponential decay after
sufficiently long time. However, as the system size increases, a finite size
scaling of $\Delta_{c}(L)$ (the inset of Fig.2 (d)) suggests that
$\Delta_{c}(L)$ vanishes in the thermodynamic limit as $\Delta_{c}(L)\sim
L^{-z}$ with $z=2$. This result indicates a divergence of the relaxation time
with the system size as $\tau_{c}\sim L^{z}$, where z is the dynamical
critical exponent characterizing the dynamical universality classHohenberg and
Halperin (1977). It is known that in a relaxation dynamics of a critical
kinetic Ising model without periodical driving, $z$ could be either $2$Glauber
(1963) or $3$Huse (1986) depending on the conservation law during the
relaxation dynamics. Here, we find the dynamical critical exponent $z$ of our
model agrees well with the value in relaxation dynamics of undriven kinetic 2D
Ising model without conservation law (Glauber modelGlauber (1963)). In the
thermodynamic limit, $\Delta_{c}(L)$ vanishes and the exponential regime gives
ways to the algebraic regime, which persists forever.

To study other critical exponents of this dynamical phase transition, we focus
on the divergence of the relaxation time $\tau$ (the inverse of $\Delta$)
close to the critical point. We first define a renormalized relaxation time
$\tilde{\tau}=\tau/L^{z}$ with $z=2$, and plot its dependence on $\mathcal{D}$
for systems with various system sizes $L$ around the critical point. As shown
in Fig.2 (e), one can find a scaling invariant point at
$\mathcal{D}=\mathcal{D}_{c}$, a signature of a continuous phase transition.
Near the critical point, the relaxation time diverges as
$\tau\sim|T-T_{c}|^{-z\nu}$. To obtain the critical exponent $\nu$, we perform
the data collapses as shown in the insets of Fig.2 (e), from which we can
extract that the critical exponents $\nu=0.64$, which agrees with the value
($\nu=0.642$) in the 3D Ising universality class .

Another important quantity to characterize the critical behavior is the
correlation length $l$, which can be extracted from the equal-time correlation
function $S(\mathbf{r},t)=\langle
s_{1}^{x}(t)s_{\mathbf{r}}^{x}(t)\rangle_{\bm{\xi}}-\langle
s_{1}^{x}(t)\rangle_{\bm{\xi}}\langle s_{\mathbf{r}}^{x}(t)\rangle_{\bm{\xi}}$
asSandvik (2010):

$l_{x}(t)=\frac{1}{|\mathbf{q}_{x}|}\sqrt{\frac{S(\mathbf{0},t)}{S(\mathbf{q}_{x},t)}-1}$
(6)

where
$S(\mathbf{Q},t)=\frac{1}{L^{3}}\sum_{\mathbf{r}}e^{i\mathbf{Q}\cdot\mathbf{r}}S(\mathbf{r},t)$
is the structure factor of the spin-spin correlation.
$\mathbf{q}_{x}=[\frac{2\pi}{L},0,0]$ and $l_{x}(t)$ is the correlation length
along x-direction at time t. As shown in Fig.2 (f), at the critical point, the
correlation length diverges algebraically in time as $l(t)\sim t^{\beta}$ with
$\beta=0.55(3)$, which agrees with the dynamical critical exponent $\beta=1/z$
in the Glauber model.

Discussion – Despite the presence of the periodical driving of our model, we
find that its dynamical critical behaviors are similar with the relaxation
dynamics of the undriven kinetic Ising model near the critical point (the
model A in Hohenberg-Halperin’s alphabetic classification of dynamical
critical phenomenaHohenberg and Halperin (1977)). This result can be
understood by a time-independent description of stroboscopic dynamics under
periodic driving. Such a time-independent Hamiltonian is similar with the
Floquet Hamiltonian in the periodically-driven quantum systems, and can be
derived using an approximation scheme similar with the Magnus expansion, a
series expansion in terms of the driving periodMagnus (1954). Even though
higher order terms (the four-site or multi-site interactions for instance) may
appear in the effective time-independent Hamiltonian, they could be irrelevant
for the long-time dynamics near the critical point, thus the stroboscopic
dynamics of our periodically driven system is similar with the relaxation
dynamics of a time-independent kinetic Ising model near the critical point.

This similarity reminds us of the dynamical phase transition observed in a
kinetic Ising model in a periodically oscillating magnetic fieldKorniss et al.
(2000), which also experiences a non-equilibrium phase transition from a
symmetry-restoring oscillations to symmetry-breaking oscillations for the
time-dependent magnetization, where the dynamical critical behavior is also
found to be the same as those of the undriven casesFujisaka et al. (2001). In
spite of the similarity of this modelKorniss et al. (2000) and our model,
there is an important difference that in our model, instead of external field,
the periodic driving is imposed on the strength of interaction, which give
rise to the consequence of the spontaneous breaking of discrete TTS. While in
model with oscillating fieldKorniss et al. (2000), the period of the
magnetization and external field is always the same, thus the discrete TTS is
not broken, corresponding to “synchronized phase” in our language.

Another important issue has not been discussed so far is the quantum effect:
what will happen if thermal noises are replaced by quantum fluctuations? In
another word, if the Hamiltonian.(1) is a quantum transverse Ising model, can
we observe DTC behavior in such a periodical-driven quantum system? The effect
of quantum fluctuations on systems with continuous TTS breaking has been
studied in the single-particle level. For instance, for a Van der Pol
oscillator with a classical trajectory of limiting circle, quantum fluctuation
induces a phase diffusion in the limiting circleNavarrete-Benlloch et al.
(2017), which recovers the continuous TTS symmetry thus destroys the long-
range TC order. On the contrary, our model breaks the discrete TTS, which is
expected to be more stable compared to the continuous ones. The competition
between the quantum fluctuations and time crystalline order may give rise to
rich dynamical phenomena, thus is worthwhile to be explored in the future
work, even though a numerical simulation of such a non-equilibrium 3D (open)
quantum many-body system is a formidable, if not impossible, task.

Conclusion and outlook – In conclusion, we study the stability of DTC phase in
the presence of thermal bath, and find that despite its genuine non-
equilibrium feature, this DTC phase shares a lot of common properties with
equilibrium ordered phases with spatially $Z_{2}$ symmetry breaking,
especially its finite-size effect and critical behavior. Future developments
will include an analytic explanation of the critical behavior, which calls for
a coarse-grained effective description and a non-equilibrium field theory
analysis of our system in the Keldysh formalismMartin et al. (1973); Kamenev
(2011); Altland and Simons (2010). Recently, Natsheh et al analytically
studied the critical behavior of a periodically driven quantum O(N) model with
DTC phaseNatsheh et al. (2021), a generalization of this work may shed light
on our problem. Another important question involves the generality of our
conclusion: whether similar behavior can be observed in other DTC models (e.g.
the coupled driven-dissipative nonlinear pendulumsYao et al. (2020))? Last but
not the least, whether it is possible to find the melting transitions with the
dynamical universality that go beyond the undriven cases? After all, the
Magnus expansion could generate long-range interactions in the effective
Hamiltonian, especially in the presence of low-frequency driving, and it is
possible that such long-range interaction may qualitatively change the
critical behavior of the systems and give rise to intriguing dynamical
universality classes.

Acknowledgments.—This work is supported by the National Key Research and
Development Program of China (Grant No.2020YFA0309000), Natural Science
Foundation of China (Grant No.12174251), Natural Science Foundation of
Shanghai (Grant No.22ZR142830), Shanghai Municipal Science and Technology
Major Project (Grant No.2019SHZDZX01).

## References

  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Taeuber (2014) U. Taeuber, _CRITICAL DYNAMICS: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior_ ( Cambridge University Press, Cambridge, 2014).
  * Kardar et al. (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 56, 889 (1986).
  * Shapere and Wilczek (2012) A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012).
  * Li et al. (2012) T. Li, Z.-X. Gong, Z.-Q. Yin, H. T. Quan, X. Yin, P. Zhang, L.-M. Duan, and X. Zhang, Phys. Rev. Lett. 109, 163001 (2012).
  * Wilczek (2013) F. Wilczek, Phys. Rev. Lett. 111, 250402 (2013).
  * Sacha (2015) K. Sacha, Phys. Rev. A 91, 033617 (2015).
  * Else et al. (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Khemani et al. (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Yao et al. (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Russomanno et al. (2017) A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Phys. Rev. B 95, 214307 (2017).
  * Gong et al. (2018) Z. Gong, R. Hamazaki, and M. Ueda, Phys. Rev. Lett. 120, 040404 (2018).
  * Huang et al. (2018) B. Huang, Y.-H. Wu, and W. V. Liu, Phys. Rev. Lett. 120, 110603 (2018).
  * Iemini et al. (2018) F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. 121, 035301 (2018).
  * Das et al. (2018) P. Das, S. Pan, S. Ghosh, and P. Pal, Phys. Rev. D 98, 024004 (2018).
  * Zhu et al. (2019) B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, New Journal of Physics 21, 073028 (2019).
  * Liao et al. (2019) L. Liao, J. Smits, P. van der Straten, and H. T. C. Stoof, Phys. Rev. A 99, 013625 (2019).
  * Kozin and Kyriienko (2019) V. K. Kozin and O. Kyriienko, Phys. Rev. Lett. 123, 210602 (2019).
  * Cai et al. (2020) Z. Cai, Y. Huang, and W. V. Liu, Chinese Physics Letters 37, 050503 (2020).
  * Chinzei and Ikeda (2020) K. Chinzei and T. N. Ikeda, Phys. Rev. Lett. 125, 060601 (2020).
  * Lyu et al. (2020) C. Lyu, S. Choudhury, C. Lv, Y. Yan, and Q. Zhou, Phys. Rev. Research 2, 033070 (2020).
  * Choudhury and Liu (2021) S. Choudhury and W. V. Liu, arXiv e-prints arXiv:2109.05318 (2021), eprint 2109.05318.
  * Bruno (2013) P. Bruno, Phys. Rev. Lett. 111, 070402 (2013).
  * Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
  * Choi et al. (2017) S. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, et al., Nature 543, 221 (2017).
  * Zhang et al. (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, et al., Nature 543, 217 (2017).
  * Autti et al. (2018) S. Autti, V. B. Eltsov, and G. E. Volovik, Phys. Rev. Lett. 120, 215301 (2018).
  * Träger et al. (2021) N. Träger, P. Gruszecki, F. Lisiecki, F. Groß, J. Förster, M. Weigand, H. Głowiński, P. Kuświk, J. Dubowik, G. Schütz, et al., Phys. Rev. Lett. 126, 057201 (2021).
  * Keßler et al. (2021) H. Keßler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, Phys. Rev. Lett. 127, 043602 (2021).
  * Mi et al. (2021) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, et al., Nature 601, 531 (2022).
  * Frey and Rachel (2021) P. Frey and S. Rachel, arXiv e-prints arXiv:2105.06632 (2021), eprint 2105.06632.
  * Yao et al. (2020) N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Nat. Phys. 16, 438 (2020).
  * Zhuang et al. (2021) Q. Zhuang, F. Machado, N. Y. Yao, and M. P. Zaletel, arXiv e-prints arXiv:2110.00585 (2021), eprint 2110.00585.
  * Lakshmanan (2010) M. Lakshmanan, Phil. Trans. R. Soc. A. 369, 1280 (2010).
  * Brown (1963) W. F. Brown, Phys. Rev. 130, 1677 (1963).
  * Lupo and Weber (2019) C. Lupo and C. Weber, Phys. Rev. B 100, 195431 (2019).
  * Gambetta et al. (2019) F. M. Gambetta, F. Carollo, A. Lazarides, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. E 100, 060105 (2019).
  * Khasseh et al. (2019) R. Khasseh, R. Fazio, S. Ruffo, and A. Russomanno, Phys. Rev. Lett. 123, 184301 (2019).
  * Hurtado-Gutiérrez et al. (2020) R. Hurtado-Gutiérrez, F. Carollo, C. Pérez-Espigares, and P. I. Hurtado, Phys. Rev. Lett. 125, 160601 (2020).
  * Ye et al. (2021) B. Ye, F. Machado, and N. Y. Yao, Phys. Rev. Lett. 127, 140603 (2021).
  * Pizzi et al. (2021) A. Pizzi, A. Nunnenkamp, and J. Knolle, Phys. Rev. Lett. 127, 140602 (2021).
  * Ament et al. (2016) S. Ament, N. Rangarajan, A. Parthasarathy, and S. Rakheja, arXiv e-prints arXiv:1607.04596 (2016), eprint 1607.04596.
  * (44) See the Supplementary Material for the details of our numerical method, convergence check of our numerical result as well as an analysis of the finite size effect.
  * Yang and Cai (2021) X. Yang and Z. Cai, Phys. Rev. Lett. 126, 020602 (2021).
  * Chinzei and Ikeda (2021) K. Chinzei and T. N. Ikeda, arXiv e-prints arXiv:2110.00591 (2021), eprint 2110.00591.
  * Glauber (1963) R. J. Glauber, J. Math. Phys. 4, 294 (1963).
  * Huse (1986) D. A. Huse, Phys. Rev. B 34, 7845 (1986).
  * Sandvik (2010) A. W. Sandvik, AIP Conference Proceedings 1297, 135 (2010).
  * Magnus (1954) W. Magnus, Communications on Pure and Applied Mathematics 7, 649 (1954).
  * Korniss et al. (2000) G. Korniss, C. J. White, P. A. Rikvold, and M. A. Novotny, Phys. Rev. E 63, 016120 (2000).
  * Fujisaka et al. (2001) H. Fujisaka, H. Tutu, and P. A. Rikvold, Phys. Rev. E 63, 036109 (2001).
  * Navarrete-Benlloch et al. (2017) C. Navarrete-Benlloch, T. Weiss, S. Walter, and G. J. de Valcárcel, Phys. Rev. Lett. 119, 133601 (2017).
  * Martin et al. (1973) P. C. Martin, E. D. Siggia, and H. A. Rose, Phys. Rev. A 8, 423 (1973).
  * Kamenev (2011) A. Kamenev, _Field theory of non-equilibrium systems_ ( Cambridge University Press, Cambridge, 2011).
  * Altland and Simons (2010) A. Altland and B. Simons, _Condensed matter field theory_ (Cambridge University Press, Cambridge, 2010).
  * Natsheh et al. (2021) M. Natsheh, A. Gambassi, and A. Mitra, Phys. Rev. B 103, 224311 (2021).

Supplemental material

In this supplementary material, we first provide some details of the Heun
algorithm used in our simulation, then we numerically check the convergence of
our results with the finite discrete time step $\Delta t$. Finally, we discuss
the dependence of our results on the initial states as well as the finite size
effect in our simulations.

##  I Solving the stochastic Landau-Lifshitz-Gilbert equation: Heun algorithm

In this section, we first derive the Stratonovich form of the stochastic LLG
equation, then formulate the Heun algorithm for Stratonovich Stochastic
Differential equations (SDE). A stochastic LLG equation reads:

$\dot{\mathbf{s}}_{i}=\mathbf{h}_{i}\times\mathbf{s}_{i}-\lambda\mathbf{s}_{i}\times(\mathbf{s}_{i}\times\mathbf{h}_{i})$
(7)

where $\mathbf{s}_{i}$ is a unit vector.
$\mathbf{h}_{i}(t)=\mathbf{h}^{0}_{i}(t)+\bm{h}^{T}_{i}(t)$ is the effective
magnetic field $\mathbf{h}^{0}_{i}(t)=[V(t)\bar{s}_{i}^{x},0,h_{z}]$ with
$\bar{s}_{i}^{x}=\frac{1}{2}\sum_{\langle j\rangle}s_{j}^{x}$ where the
summation is over the nearest neighboring sites of site i. $\bm{h}^{T}_{i}(t)$
is a three-dimensional(3D) stochastic magnetic field representing the thermal
noise satisfying:

$\langle
h_{i}^{T\alpha}(t)h_{j}^{T\beta}(t^{\prime})\rangle_{\bm{\xi}}=\mathcal{D}^{2}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$
(8)

where $\alpha,\beta=x,y,z$ and $\mathcal{D}$ is the strength of the noise
satisfying $\mathcal{D}^{2}=2T\lambda$, and the average
$\langle\rangle_{\bm{\xi}}$ is over all the trajectories of noises.

To simulate the thermal noise numerically, we discretize the time with the
time step $\Delta t$ of the numerical method. Provided that the spin
configuration in the $m$-th time step ($t_{m}=m\Delta t$) is defined as
$\\{\mathbf{s}_{i}^{m}\\}$, in the Heun algorithm, the calculation of
$\\{\mathbf{s}_{i}^{m+1}\\}$ can be divided into two steps. We first
calculate:

$\tilde{\mathbf{s}}_{i}^{m+1}=\mathbf{s}_{i}^{m}+[\mathbf{h}^{m}_{i}\times\mathbf{s}_{i}^{m}-\lambda\mathbf{s}_{i}^{m}\times(\mathbf{s}_{i}^{m}\times\mathbf{h}^{m}_{i})]\Delta
t$ (9)

with $\mathbf{h}_{i}^{m}=\mathbf{h}^{0}_{i,m}+\bm{\tilde{h}}^{T}_{i,m}$, where
$\mathbf{h}^{0}_{i,m}=\mathbf{h}^{0}_{i}(t_{m})$ and
$\bm{\tilde{h}}^{T}_{i,m}$ is a stochastic magnetic field satisfying:

$\tilde{h}^{T\alpha}_{i,m}=\frac{\mathcal{D}}{\sqrt{\Delta
t}}\xi_{i,m}^{\alpha}$ (10)

where $\xi_{i}^{\alpha}$ is a random number satisfying the Gaussian
distribution with $\mathcal{N}(0,1)$:
$\langle\xi_{i}^{\alpha}\rangle_{\bm{\xi}}=0$,
$\langle(\xi_{i}^{\alpha})^{2}\rangle_{\bm{\xi}}=1$.

In the Heun algorithm, $\mathbf{s}_{i}$ at the $m+1$-th time step can be
expressed as:

$\displaystyle\mathbf{s}_{i}^{m+1}$
$\displaystyle=\mathbf{s}_{i}^{m}+\frac{\Delta
t}{2}\\{\mathbf{h}^{m}_{i}\times\mathbf{s}_{i}^{m}-\lambda\mathbf{s}_{i}^{m}\times(\mathbf{s}_{i}^{m}\times\mathbf{h}^{m}_{i})$
(11) $\displaystyle+$
$\displaystyle\tilde{\mathbf{h}}^{m+1}_{i}\times\tilde{\mathbf{s}}_{i}^{m+1}-\lambda\tilde{\mathbf{s}}_{i}^{m+1}\times(\tilde{\mathbf{s}}_{i}^{m+1}\times\tilde{\mathbf{h}}^{m+1}_{i})\\}$

where $\tilde{\mathbf{s}}_{i}^{m+1}$ has been defined in Eq.(9), and
$\tilde{\mathbf{h}}_{i}^{m+1}=\mathbf{h}^{0}_{i,m+1}+\bm{\tilde{h}}^{T}_{i,m}$.

Figure 3: (Color online) Comparison of the M(t) with parameters $L=8$,
$\mathcal{D}=0.2782$, $\mathcal{N}=1000$ and different $\Delta t$,  Figure 4:
The long-time dynamics of the FM order parameter M(t) in the DTC phase with
parameters $L=20$, $D=0.26J$. Figure 5: (Color online) Comparison of the M(t)
starting from different initial states with parameters $L=20$,
$\mathcal{D}=0$.

Figure 6: (Color online) (a) The envelop of M(t) in the DTC phases of small
systems with various system sizes and parameter $D=0.26J<D_{c}$. (b) The decay
rate $\Delta(L)$ as a function of system size L.  Figure 7: Replotting Fig. 2
(d) in the main text using a semi-log plot

##  II  Convergence of numerical results

In our simulation, we choose the time step $\Delta t=10^{-3}J^{-1}$. In
general, for stochastic differential equations, the dependence of the
numerical results on $\Delta t$ is more subtle than the deterministic ones (as
shown in Eq.(8), the strength of the stochastic magnetic fields depend on
$\Delta t$) thus we need to carefully examine the convergence of our result
(especially the power-law decays) with $\Delta t$, and preclude the
possibility that it is an artifact because of the finite $\Delta t$ we choose.
To this end, we choose different $\Delta t=2\times 10^{-3},10^{-3}$ and
$5\times 10^{-4}$, and compare their results. As shown in Fig.3, the results
with $\Delta t=10^{-3}$ and $5\times 10^{-4}$ agree with each other very well
within the statistical error bar from the ensemble average of the trajectories
of noise, which indicates that the $\Delta t$ we choose in our simulation is
sufficiently small that allows us to ignore the finite-$\Delta t$ induced
errors.

The typical time scale of the simulation in the maintext is up to $\sim
100J^{-1}$ (50 DTC periods). One may wonder whether such a time scale is in a
prethermal regime, and longer simulation may give rise to qualitatively
different dynamics. To verify this point, we extended the simulation time up
to $10^{3}J^{-1}$ (500 periods of the DTC) , which is the maximum accessible
time scale in our numerical simulations considering the accumulated errors in
the Heun algorithm. As shown in Fig.4, the DTC order parameter barely decays
for a sufficiently large system ($L=20$). However, as we will show below, for
a small system, there is an exponential decay even at temperature below the
critical temperature, which can be considered as a finite-size effect.

##  III Dependence on the initial state

In our simulation in the main text, we start from an spatially homogeneous
initial state (the ground state of the system Hamiltonian at $t=0$), where the
FM order parameter $M(t)$ is proportional to the auto-correlation function in
time $C(t)$. It is important to check that our results doesn’t crucially
depend on this specific choice of the initial state. To this end, we choose an
inhomogeneous random initial state (paramagnetic state with $M(t=0)=0$): for
each site, we choose the the initial spin as
$\bm{s}_{i}^{0}=[s_{i}^{x},0,s_{i}^{z}]$, where $s_{i}^{x}$ is an random
number different from site by site and uniformly distributed within [-1,1],
the corresponding z-component spin is chosen as
$s_{i}^{z}=\sqrt{1-[s_{i}^{x}]^{2}}$. We compare the $M(t)$ from an uniform
ferromagnetic (FM) and random initial states. As shown in Fig.5, after a
relaxation dynamics, the time evolution of $M(t)$ in these two cases agree
with each other very well within the statistical error, indicating a rapid
loss of the memory of the initial state information, which can be understood
as a consequence of the coupling to the thermal bath.

##  IV Finite size effect in the DTC phases and critical point

All the numerical simulations in the main text are preformed in finite size
systems. In equilibrium physics, it is well known that the spontaneous
symmetry breaking can only occur in thermodynamic limit, while a symmetry
breaking phase in a finite system has a life time, which exponentially
diverges with the system size. Since the DTC phases in our model also
spontaneously breaks the discrete translational symmetry, it is interesting to
ask whether such a non-equilibrium symmetry breaking phase possesses similar
properties.

To this end, we study the dynamics of the FM order parameter for small systems
with fixed $D=0.26J$ below the critical temperature. As shown in Fig.6 (a),
M(t) in small systems with various system size indeed decay exponentially in
time ($|M(t)|\sim e^{-\Delta(L)t})$, which indicates a finite life time
$\tau_{c}\sim\Delta(L)^{-1}$ for the DTC phase. A finite size scaling of
$\Delta(L)$ is plotted in Fig.6 (b), which shows an exponential decay of
$\Delta(L)$ with system size (an exponential divergence of $\tau_{c}$). This
result indicates that despite the genuine non-equilibrium feature of the DTC
phases, it share with some common features with the spontaneous symmetry
breaking phases in equilibrium systems.

Now we focus on the finite-size effect at the critical point. As shown in
Fig.7, for a finite system at the critical point, the long-time dynamics also
exhibit an exponential decay due to the finite-size effect, whose exponent
$\Delta$ depends on the system size. However, different from the DTC phase
below the critical temperature, at the critical point, $\Delta(L)$ decay
algebraically instead of exponentially with $L$ (see the inset of Fig.2 (d) in
the maintext), from which we can extract a dynamical critical exponent $z$
($\Delta(L)\sim L^{-z}$).


The COSINUS Collaboration

# First measurements of remoTES cryogenic calorimeters:
easy-to-fabricate particle detectors for a wide choice of target materials

G. Angloher  Max-Planck-Institut für Physik, 80805 München - Germany  M.R.
Bharadwaj  Max-Planck-Institut für Physik, 80805 München - Germany  I. Dafinei
INFN - Sezione di Roma, 00185 Roma - Italy  N. Di Marco  Gran Sasso Science
Institute, 67100 L’Aquila - Italy  INFN - Laboratori Nazionali del Gran Sasso,
67010 Assergi - Italy  L. Einfalt  Institut für Hochenergiephysik der
Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria
Atominstitut, Technische Universität Wien, 1020 Wien - Austria  F. Ferroni
INFN - Sezione di Roma, 00185 Roma - Italy  Gran Sasso Science Institute,
67100 L’Aquila - Italy  S. Fichtinger  Institut für Hochenergiephysik der
Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria  A.
Filipponi  Dipartimento di Scienze Fisiche e Chimiche, Università degli Studi
dell’Aquila, 67100 L’Aquila - Italy  INFN - Laboratori Nazionali del Gran
Sasso, 67010 Assergi - Italy  T. Frank  Max-Planck-Institut für Physik, 80805
München - Germany  M. Friedl  Institut für Hochenergiephysik der
Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria  A. Fuss
Institut für Hochenergiephysik der Österreichischen Akademie der
Wissenschaften, 1050 Wien - Austria  Atominstitut, Technische Universität
Wien, 1020 Wien - Austria  Z. Ge  SICCAS - Shanghai Institute of Ceramics,
Shanghai - P.R.China 200050  M. Heikinheimo  Helsinki Institute of Physics,
Univ. of Helsinki, 00014 Helsinki - Finland  K. Huitu  Helsinki Institute of
Physics, Univ. of Helsinki, 00014 Helsinki - Finland  M. Kellermann  Max-
Planck-Institut für Physik, 80805 München - Germany  R. Maji  Institut für
Hochenergiephysik der Österreichischen Akademie der Wissenschaften, 1050 Wien
- Austria  Atominstitut, Technische Universität Wien, 1020 Wien - Austria  M.
Mancuso  Max-Planck-Institut für Physik, 80805 München - Germany  L. Pagnanini
INFN - Laboratori Nazionali del Gran Sasso, 67010 Assergi - Italy  Gran Sasso
Science Institute, 67100 L’Aquila - Italy  F. Petricca  Max-Planck-Institut
für Physik, 80805 München - Germany  S. Pirro  INFN - Laboratori Nazionali del
Gran Sasso, 67010 Assergi - Italy  F. Pröbst  Max-Planck-Institut für Physik,
80805 München - Germany  G. Profeta  Dipartimento di Scienze Fisiche e
Chimiche, Università degli Studi dell’Aquila, 67100 L’Aquila - Italy  INFN -
Laboratori Nazionali del Gran Sasso, 67010 Assergi - Italy  A. Puiu  Gran
Sasso Science Institute, 67100 L’Aquila - Italy  INFN - Laboratori Nazionali
del Gran Sasso, 67010 Assergi - Italy  F. Reindl  Corresponding author,
<EMAIL_ADDRESS>Institut für Hochenergiephysik der
Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria
Atominstitut, Technische Universität Wien, 1020 Wien - Austria  K. Schäffner
Corresponding author<EMAIL_ADDRESS>Max-Planck-Institut für
Physik, 80805 München - Germany  J. Schieck  Institut für Hochenergiephysik
der Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria
Atominstitut, Technische Universität Wien, 1020 Wien - Austria  D.
Schmiedmayer  Institut für Hochenergiephysik der Österreichischen Akademie der
Wissenschaften, 1050 Wien - Austria  Atominstitut, Technische Universität
Wien, 1020 Wien - Austria  C. Schwertner  Institut für Hochenergiephysik der
Österreichischen Akademie der Wissenschaften, 1050 Wien - Austria
Atominstitut, Technische Universität Wien, 1020 Wien - Austria  M. Stahlberg
Corresponding author<EMAIL_ADDRESS>Max-Planck-Institut für
Physik, 80805 München - Germany  A. Stendahl  Helsinki Institute of Physics,
Univ. of Helsinki, 00014 Helsinki - Finland  F. Wagner  Institut für
Hochenergiephysik der Österreichischen Akademie der Wissenschaften, 1050 Wien
- Austria  S. Yue  SICCAS - Shanghai Institute of Ceramics, Shanghai -
P.R.China 200050  V. Zema  Corresponding author<EMAIL_ADDRESS>Max-
Planck-Institut für Physik, 80805 München - Germany  Y. Zhu  SICCAS - Shanghai
Institute of Ceramics, Shanghai - P.R.China 200050  A. Bento  Max-Planck-
Institut für Physik, 80805 München - Germany  LIBPhys-UC, Physics Departments,
University of Coimbra, 3004-516 Coimbra - Portugal  L. Canonica  Max-Planck-
Institut für Physik, 80805 München - Germany  A. Garai  Max-Planck-Institut
für Physik, 80805 München - Germany

###### Abstract

Low-temperature calorimeters based on a readout via transition edge sensors
(TESs) and operated below 100 mK are well suited for rare event searches. We
present first experimental results from two detector prototypes using a novel
thermometer design denoted remoTES. This design facilitates the use of TESs in
combination with absorber materials which, due to their physical and chemical
properties, as e.g. hygroscopicity, low hardness and low melting point,
prevent the direct deposition of the TES onto their surface. In a remoTES
detector, the TES is fabricated onto a separate wafer. The absorber crystal is
then equipped with a gold pad that transmits the phonon signal created from an
interaction in the absorber to the thermometer via a gold bonding wire.

With recent prototype detectors operated in an above-ground R&D facility, we
achieve energy resolutions of $\sigma=87.8$ eV for a silicon absorber and
$\sigma=193.5$ eV for an -TeO2 absorber, respectively. RemoTES calorimeters
offer – besides the wider choice of absorber materials – a simpler production
process combined with a higher reproducibility for large detector arrays and
an enhanced radiopurity standard, which is of particular interest for rare
event searches.

##  I Introduction

To date, low-temperature detectors operated in the temperature region of
10-100 mK make up an important detector class in particle physics. The
detection of particles at temperatures close to absolute zero, in so-called
cryogenic calorimeters, is based on the measurement of the tiny temperature
increase $\mathcal{O}$(µK) caused by a particle interaction in an absorber
material (inorganic crystalline insulators, semiconductors, metals) by using
very sensitive temperature sensors (Transition Edge Sensors (TESs), Neutron
Transmutation Doped-thermistors (NTDs), Kinetic Inductance Detectors (KIDs),
or Metallic Magnetic Calorimeters (MMCs)).

A TES consists of a superconductive thin film (mono-layer or bi-layers)
operated in the state between the normal conducting and the superconducting
phase. Thanks to the strong temperature dependence of the film resistance in
the transition, even tiny temperature increases induced by very small energy
depositions in the absorber lead to measurable changes in resistance. For many
decades, TESs have been a main pillar in the scope of cryogenic imaging
spectrometers for earth and space applications; their application in the field
of rare event searches for direct dark matter detection (e.g. CRESST [1],
(Super)CDMS [2, 3], COSINUS [4]) and CE$\nu$NS (e.g. NUCLEUS [5], RICOCHET
[6]) by using massive absorbers $\mathcal{O}$(1-100 g) is a very active field
of research at the frontier of knowledge.

TESs based on superconducting tungsten thin-films have been developed within
the CRESST dark matter search at the Max-Planck Institute for Physics in
Munich over the last 30 years. In the current stage, called CRESST-III, the
200 nm-thick tungsten film is directly fabricated onto the surface of the
absorber (Al2O3, CaWO4, LiAlO2 and Si). The extraordinary high sensitivity
achieved with TESs is one of the cornerstones of CRESST being the leading
cryogenic experiment for low-mass dark matter search [1, 7].

The fabrication of the TES thermometer directly onto the absorber ensures an
excellent transmission of the non-thermal phonons to the W-film and is the
preferred solution, if the absorber material can withstand the conditions
during the involved steps (electron-beam evaporation, sputtering, chemical
etching, photolithography). However, absorber materials that have a low
melting point and/or are hygroscopic cannot undergo these fabrication
processes.

A work-around solution for such delicate materials, e.g. sodium iodide (NaI)
used in the COSINUS experiment [4], is the so-called composite design [8]
where the TES is fabricated onto a separate substrate named ”carrier crystal”
(Al2O3, Si, CaWO4, CdWO4) which in turn is connected to the absorber via an
amorphous interface (glue, oil, or grease). This approach requires the signal
from a particle interaction in the absorber to propagate through the
interface, the carrier crystal, and finally couple to the TES. Signal loss is
expected in this process due to the acoustic mismatch (reflection vs.
transmission) [9] between the different materials. Furthermore, in the case of
e.g. NaI, which is a very good scintillator, a carrier crystal which is not
fully transparent to the scintillation light will reabsorb a part of the
light, and the resulting pulse shape will be a superposition of both
scintillation and thermal signal.

In [10], the authors propose a TES-based detector concept which does not
suffer from these disadvantages, features a simpler production process, and
additionally promises a better detector reproducibility. To our knowledge,
experimental works using this design in combination with TES are absent to
present date, although the experiment RICOCHET is considering this kind of
detector design to read out superconducting absorbers for precision
measurements of coherent elastic neutrino-nucleus scattering (CENS). The AMORE
experiment [11] utilizes 400 nm thin gold films as phonon absorbers on their
massive 48deplCa100MoO4 target crystals for the search of the neutrinoless
double beta decay of 100Mo. In the AMORE detector design the Au-pad is then
connected to a metallic magnetic calorimeter (MMC) via 25 annealed Au bonding
wires [12].

Marrying ease of fabrication with a wide choice of absorber materials, the
concept denoted remoTES in the following, may be a breakthrough for TES-
instrumented calorimeters based on delicate absorber materials like NaI. In
this article, we present first experimental results which demonstrate the
potential of such detectors.

##  II Experimental technique

###  II.1 Detector design

The design of a remoTES detector is schematically depicted in Figure 1. It
completely avoids the carrier crystal by coupling the TES directly to the
target crystal via a gold pad and a gold bonding wire. A significant advantage
of this design is that it eliminates the potential phonon barriers, i.e. the
loss of signal due to reflections caused by the acoustic mismatch between
absorber crystal, interface and carrier crystal. Even more importantly,
phonons directly couple to the electron system of the gold pad via electron-
phonon-coupling [13]. The gold pad introduces two competing effects: on the
one hand, the collection of phonons should be enhanced as the electron-phonon-
coupling in gold is larger than in tungsten [10]. On the other hand, the gold
adds an additional heat capacity potentially decreasing the sensitivity.
However, an overall net positive effect of the remoTES design is expected.

Figure 1: Schematic view of a remoTES detector. The wafer where the TES is
fabricated onto (right) is separated from the absorber crystal (left). This
crystal is equipped with a gold pad which is linked directly to the TES via a
gold bonding wire.

In this manuscript we report the results from two different absorber
materials, silicon (Si, 2.33 g) and tellurium dioxide (-TeO2, 2.27 g), read
out by a remoTES. These materials were selected to gain a fundamental
understanding of the crucial parameters of this new TES design. We selected Si
as it is a highly-studied standard material for low-temperature application
and, thus, useful as a benchmark. We chose TeO2 as a first step on the way to
NaI, mainly for its similar solid-state properties, in particular regarding
the heat capacity and the gap between the acoustic and optical phonon modes
[14, 15]. We already performed a first measurement with a NaI absorber which
is subject of a future publication based on NaI-absorbers only.

On the Si absorber ((20x10x5) mm3) a gold film of circular shape (3 mm
diameter) of 200 nm thickness was deposited by Magnetron-sputtering. The
connection between the Au-pad on the Si absorber and the TES on the wafer was
provided via a 17 µm thick and 2-3 mm long Au-bonding wire. While the Au-wire
was wedge-bonded onto the Au-pad of the TES, it was connected to the Au-pad on
the Si absorber by using a silver-loaded epoxy ($\sim$25 µg). Avoiding wedge-
bonding on the Au-pad of the absorber opens up the possibility to also
introduce very soft and fragile absorber crystals to this field of research,
which would be damaged by solid phase welding processes (e.g. wedge and ball
bonding). Alternatively, the Au-pad on the TeO2 absorber ((20x10x2) mm3)
consists of a 400 nm thick Au-foil glued onto its surface by using a two-
component epoxy resin. Using Au-foils which reveal high quality and increased
thickness (Residual Resistance Ratio (RRR), listed in Table 1), allow to
exclude the absorber material from any kind of thin-film deposition processes
which require high vacuum conditions. The electrical and thermal contact
between Au-pad on the TeO2 and the TES was realized with a wedge-bond on both
sides, respectively. Table 1 summarizes the condensed information on both
prototype detectors.

The TES on the thin Al2O3 wafer ((10x10x0.4) mm3) was identical for both
measurements. It consists of a 100 nm-thick tungsten film of an area of
(220x300) µm2 and two aluminum pads for biasing of the sensor. The area of the
W-film not overlapping with the pads is (70x300) µm2. The weak thermal link to
the heat bath consists of a 1 mm long Au-stripe (20 µm in width, 80 nm in
thickness).

Absorber | Absorber | Au-pad | Au-wire | TES | Energy
---|---|---|---|---|---
material | volume (mm3) | properties | properties |  | resolution (eV)
Si | 20x10x5 | 200nm | 17 µm | W-TES | $87.8\pm 5.6$
|  | sputtered | glued on pad | on Al2O3 |
|  | RRR=3.79 |  |  |
TeO2 | 20x10x2 | 400nm foil | 17 µm | W-TES | $193.5\pm 3.1$
|  | glued | 2 wedge bonds | on Al2O3 |
|  | RRR=15 |  |  |
Table 1: Summary of the two measurements using the remoTES design. RRR =
Residual Resistance Ratio

###  II.2 Experimental setup

The measurements reported here were carried out in a dilution refrigerator of
the CRESST group at the Max-Planck Institute for Physics in Munich. This is an
above-ground wet dilution cryostat for R&D measurements and TES testing. It is
equipped with four SQUID channels (APS company) and a VME-based data
acquisition using a hardware-triggered readout. The total measuring time was
12.54 h for the Si detector, and 37.15 h for TeO2. This corresponds to an
exposure of 1.06 g$\cdot$d and 2.28 g$\cdot$d, respectively. For means of
calibration both prototype detectors were irradiated by an uncollimated 55Fe-
source shining onto the bottom side of the absorbers. Furthermore, the TeO2
detector was equipped with an additional collimated 55Fe-source irradiating
solely the Au-pad glued onto the absorber.

##  III Experimental results and discussion

###  III.1 Data analysis

We observe different event classes in both detectors. Standard events (SEVs
$=$ averaged pulses scaled to one) for these classes are shown in Figure 2:
the most prevalent pulse shape (depicted in green) features a rise-time of
only a few milliseconds, with a decay time of approximately 50 ms for Si and
several hundred ms for TeO2. The high rate of these events is consistent with
hits in the absorber, considering its large mass in comparison with the mass
of the wafer carrying the TES. Hits in the latter (orange in Figure 2) feature
a lower rate, and much faster signals, with decay times of only few ms.

The TeO2 prototype exhibits an additional event class. We identify these
comparatively rare events (cf. Table 2) as hits from the collimated 55Fe-
source irradiating the Au-pad on the absorber (see Appendix A). The rise time
of the pulse is driven by the time needed to transport energy to the TES which
depends on the energy carriers. The transport is fast for electrons and non-
thermal phonons, which are the dominant energy carriers for the events created
by interactions in the Au-pad or the carrier. Energy depositions in the TeO2
absorber, instead, undergo a stronger thermalization and consequently lead to
slower pulses. The decay time of Au-pad events is faster than the one of
absorber events and slower than the one of wafer events. The explanation can
be two-fold: firstly, energy deposition in the Au-pad are partially
transmitted back and forth to and from the absorber, leading to a larger decay
time; secondly, the heat capacity of Au affects the decay time, leading to
even slower pulses. A definitive interpretation of the pulse shape requires
further studies, which are currently ongoing.

The SEVs in Figure 2 were created from a small, cleaned sample of pulses for
each event class, where all pulses of one set have the same shape, and are
scaled to the same pulse height before averaging. Few basic quality cuts are
applied to the complete dataset: events sitting on a decaying baseline from a
previous pulse are removed by placing a cut on the pre-trigger baseline slope.
Additionally, events with a reconstructed pulse onset more than 1 ms away from
the trigger position are discarded. To separate different event classes
reliably, the root mean square (RMS) from a truncated standard event fit (see
e.g. [7]) is used. We fit each pulse with all three SEVs, and only those
events where the respective fit yields the lowest RMS are kept for each event
class dataset. In the fit amplitude versus RMS plane, pulses which deviate
from the main population at low RMS values are removed for each event class.
Table 2 shows the number of surviving events for each measurement and event
class; as stated before, absorber events are the dominant event class in both
measurements, and Au-hits are less frequent than wafer hits for the TeO2
detector.

Absorber | ”Absorber” | ”Wafer” | ”Gold”
---|---|---|---
material | events | events | events
Si | 6072 | 398 | -
TeO2 | 28555 | 711 | 314
Table 2: Number of surviving events after quality cuts for each measurement
and event class (see text for details).

Apart from events that are caused by particle interactions in detector parts,
electrical pulses with known amplitudes are periodically injected into the
heater of each detector. These heater pulses allow monitoring the detector
response over time, and accounting for the effect of saturation in the TES,
which influences the pulse shape at high energies (see e.g. [7]). A truncated
fit with a dedicated SEV is performed for all heater pulses, and the resulting
amplitude is used in combination with the input pulser setting to compensate
non-linearities and time-variations of the reconstructed (fitted) amplitude
for particle-induced events [16]. The amount of heater pulses surviving the
quality cuts also provides an estimate of the survival probability for real
events: For the TeO2 analysis, 3464 heater pulses out of 5944 survive the
baseline slope cut, the pulse onset cut, and a heater pulse fit RMS cut. This
corresponds to a survival probability of 58.3 %. Similarly, the survival
probability for the Si remoTES analysis is 44.6 % (978 out of 2192 pulses).
These values are only approximate - energy-dependent effects are not
considered.

(a)

(b)

Figure 2: Standard events (SEVs) for the two prototype detectors. (a) Run 566:
Events in the Si absorber (green) and the TES wafer (orange). (b) Run 573:
Events in the TeO2 absorber (green), the TES wafer (orange) and the gold pad
on the TeO2 (violet). The pulse onsets for each event class are slightly
different.

###  III.2 Energy calibration and resolution

The Kα peak from the 55Fe source (5.89 keV) prominent in the spectrum of each
detector is used to calibrate the detector response, which is given by the
truncated fit result for absorber events corrected by the heater response
information. We chose the truncated fit over the usual method of optimal
filtering [17] here, since the noise conditions in both measurements changed
over time, and the filter relies on a stable noise power spectrum. The
resulting spectra for absorber events are given in Figure 3.

(a)

(b)

Figure 3: The energy spectra of the two prototype detectors: (a) Si absorber
and (b) TeO2 absorber. The intensity of the 55Fe-source producing X-rays of
5.89 keV (Kα) and 6.49 keV (Kβ) was significantly stronger for (b). The
additional peaks in the Si detector ($\sim$1.8 keV and $\sim$5 keV) may be
attributed to escape lines from Cu and Si.

The baseline resolution is determined for each detector in the following way:
in each measurement, records without a trigger event are taken in regular
intervals. These ”empty baselines” are then cleaned, removing random
coincidences with pulses and artifacts. Subsequently, we superimpose the SEV
for a given detector with all remaining empty baselines. These artificial
events (122 for the Si detector and 1932 for TeO2) are passed through the
analysis chain and reconstruction algorithm, and the detector resolution is
extracted from the reconstructed amplitudes assuming that the noise
fluctuations follow a Gaussian distribution. The resolutions (all in $\sigma$)
are ($18.4\pm 1.2)$ mV for the Si and ($1.22\pm 0.02)$ mV for the TeO2
absorber, corresponding to ($87.8\pm 5.6$) eV and ($193.5\pm 3.1$) eV,
respectively (cf. Table 1). Typically, the achievable threshold of TESs is
five times the baseline resolution.

##  IV Conclusion

We successfully operated two remoTES detectors using a Si and a TeO2 absorber
in above-ground measurements, achieving baseline resolutions of 87.8 eV and
193.5 eV, respectively. We observed two event classes with distinctively
different pulse shapes in each detector, which result from particle
interactions in the absorber and in the wafer of the TES. Additionally, a
third event class was visible in the TeO2 measurement, which featured a very
low rate and long decay time, consistent with hits in the gold pad glued onto
the TeO2 absorber. These results show that remoTES detectors meet the
requirements of rare event searches and tiny energy depositions: the design
yields a low recoil energy threshold, since the phonon propagation from the
absorber to the sensor is not hindered as much as in ”composite-type”
detectors. Additionally, the remoTES prevents the undesirable collection of
scintillation light in the wafer and also provides a handle on pulse shape
discrimination via modification of the Au link between the absorber and the
TES. Absorber materials that required the use of a carrier crystal in the past
can now be operated without such a crystal by simply adding a gold interface.
The separately fabricated wafer-TES is easier to produce than the combined
system of absorber crystal, carrier and TES, and can even be reused with a
different target – this implies an improved reproducibility for remoTES
detectors with respect to their carrier design counterparts, and greatly
facilitates comparisons between different absorber materials. In summary, our
results show that the TES technology can be successfully used with almost any
absorber material. Even for ”standard” absorbers, the use of a remoTES would
avoid exposing the material to fabrication processes, thus preserving the
initial radiopurity of the absorber crystal.

## Acknowledgments

This work was carried out in the frame of the COSINUS collaboration. We would
like to acknowledge Matt Pyle, Enectali Figueroa-Feliciano and Bernard
Sadoulet who had the idea of developing such-kind of TES-based detectors, as
published already in 2015 on arXiv:1503.01200 [10]. The data analysis was
mainly carried out with collaboration-internal software tools, and partially
with the open source Python package CAIT[18]. We want to thank the MPP
mechanical workshop team for in-time support in all needs regarding the
cryogenic facility and detector holder fabrication. We are grateful to the
colleagues from the CRESST MPP group for sharing the infrastructure for
detector production, the wet dilution refrigerator facility for the
measurements and their cryogenic expertise and experience. In particular we
thank Ahmed Abdelhameed for many insightful discussions on thin-film
technology.

##  Appendix A Identification of events in the gold foil

The TeO2-remoTES is calibrated using two collimated 55Fe-sources, one
irradiating the absorber, and one irradiating the gold foil glued onto the
absorber. The X-rays from the 55Fe-source directly absorbed in the gold foil
are visible in Figure 4(a) as a cluster of events in the violet box. As
expected, they are distributed around a larger pulse height ($\sim$ 0.07 V)
with respect to the pulse height of the 55Fe-events in the absorber ($\sim$
0.04 V). As mentioned in Section III.1, the transmission of the signal to the
thermometer is more efficient for the gold foil events. For this reason, the
rise time of the events in the gold is smaller ($\lesssim$ 0.5 ms) than the
one in the absorber ($>0.8$ ms), as also evident in Figure 2(b). Figure 4(b)
shows the RMS of the gold SEV fit to the TeO2-remoTES data set as a function
of the pulse height. Here, the events from the violet box in Figure 4(a) are
marked as violet dots. Their corresponding RMS is clearly smaller than the RMS
of the absorber events.

(a)

(b)

Figure 4: TeO2-remoTES data set: (a) Rise time versus moving average pulse
height distribution. The violet box encloses the events in the gold foil
produced by the collimated 55Fe- source. (b) Fit RMS for the gold SEV as a
function of the pulse height distribution. The events from the violet box in
panel (a) are tagged and depicted in violet.

## References

  * Abdelhameed et al. [2019] A. Abdelhameed et al.  First results from the CRESST-III low-mass dark matter program.  _Phys. Rev. D_ , 100(10):102002, Nov. 2019.  doi: 10.1103/PhysRevD.100.102002.
  * Alkhatib et al. [2021] I. Alkhatib et al.  Light Dark Matter Search with a High-Resolution Athermal Phonon Detector Operated above Ground.  _Phys. Rev. Lett._ , 127(6):061801, Aug. 2021\.  doi: 10.1103/PhysRevLett.127.061801.
  * Agnese et al. [2019] R. Agnese et al.  Search for low-mass dark matter with CDMSlite using a profile likelihood fit.  _Phys. Rev. D_ , 99(6):062001, Mar. 2019.  doi: 10.1103/PhysRevD.99.062001.
  * Angloher et al. [2016] G. Angloher et al.  The COSINUS project: perspectives of a NaI scintillating calorimeter for dark matter search.  _Eur. Phys. J. C_ , 76(8):441, Aug. 2016.  ISSN 1434-6044, 1434-6052.  doi: 10.1140/epjc/s10052-016-4278-3.
  * Angloher et al. [2019] G. Angloher et al.  Exploring $\hbox{{CE}}\nu\hbox{{NS}}$ with NUCLEUS at the Chooz nuclear power plant.  _Eur. Phys. J. C_ , 79(12):1018, Dec. 2019.  ISSN 1434-6052.  doi: 10.1140/epjc/s10052-019-7454-4.
  * Billard et al. [2017] J. Billard et al.  Coherent neutrino scattering with low temperature bolometers at Chooz reactor complex.  _J. Phys. G: Nucl. Part. Phys._ , 44(10):105101, Aug. 2017.  ISSN 0954-3899.  doi: 10.1088/1361-6471/aa83d0.
  * Angloher et al. [2017] G. Angloher et al.  Results on MeV-scale dark matter from a gram-scale cryogenic calorimeter operated above ground.  _Eur. Phys. J. C_ , 77(9):637, Sept. 2017.  ISSN 1434-6044, 1434-6052.  doi: 10.1140/epjc/s10052-017-5223-9.
  * Angloher et al. [2009] G. Angloher et al.  Composite CaWO4 Detectors for the CRESST‐II Experiment.  In _AIP Conference Proceedings_ , volume 1185, pages 651–654. AIP Publishing, Dec. 2009.  doi: 10.1063/1.3292426.
  * Gray [2012] K. E. Gray.  _Nonequilibrium superconductivity, phonons, and Kapitza boundaries_ , volume 65.  Springer Science & Business Media, 2012.
  * Pyle et al. [2015] M. Pyle et al.  Optimized designs for very low temperature massive calorimeters.  _arXiv:1503.01200_ , 2015.  URL https://arxiv.org/abs/1503.01200.
  * Alenkov et al. [2019] V. Alenkov et al.  First Results from the AMoRE-Pilot neutrinoless double beta decay experiment.  _Eur. Phys. J. C_ , 79(9):791, 2019.  doi: 10.1140/epjc/s10052-019-7279-1.
  * Kim et al. [2004] Y. Kim, H. Eguchi, C. Enss, Y. Huang, R. Lanou, H. Maris, A. Mocharnuk-Macchia, G. Seidel, B. Sethumadhavan, and W. Yao.  Measurements and modeling of the thermal properties of a calorimeter having a sapphire absorber.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 520(1):208–211, 2004.  ISSN 0168-9002.  doi: https://doi.org/10.1016/j.nima.2003.11.230.  URL https://www.sciencedirect.com/science/article/pii/S0168900203031334.  Proceedings of the 10th International Workshop on Low Temperature Detectors.
  * Little [1959] W. Little.  The transport of heat between dissimilar solids at low temperatures.  _Canadian Journal of Physics_ , 37(3):334–349, 1959.  URL https://doi.org/10.1139/p59-037.
  * Jain et al. [2013] A. Jain et al.  The materials project: A materials genome approach to accelerating materials innovation.  _APL Materials_ , 1(011002), 2013.  URL http://link.aip.org/link/AMPADS/v1/i1/p011002/s1&Agg=doi.
  * Hinuma et al. [2017] Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, and I. Tanaka.  Band structure diagram paths based on crystallography.  _Computational Materials Science_ , 128:140–184, 2017.  URL https://doi.org/10.1016/j.commatsci.2016.10.015.
  * Angloher et al. [2002] G. Angloher et al.  Limits on WIMP dark matter using sapphire cryogenic detectors.  _Astropart. Phys._ , 18(1):43–55, Aug. 2002.  ISSN 0927-6505.  doi: 10.1016/S0927-6505(02)00111-1.
  * Gatti and Manfredi [1986] E. Gatti and P. F. Manfredi.  Processing the signals from solid-state detectors in elementary-particle physics.  _Riv. Nuovo Cim._ , 9(1):1–146, Jan. 1986.  ISSN 1826-9850.  doi: 10.1007/BF02822156.
  * Wagner et al. [2021] F. Wagner et al.  fewagner/cait: v1.0.0, July 2021.  URL https://doi.org/10.5281/zenodo.5091416.