\section{Supplmental material}

\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{equation}{0}
\setcounter{figure}{0}

\subsection{Numerical method to solve the Schr\"odinger equation}
The state vector at time $t$ is given by 
$|\Psi(t)\rangle = \hat{U}(t) |\Psi(0)\rangle$ with
the unitary time evolution operator 
\begin{equation}
\hat{U}(t) \equiv e^{-i \hat{H}t} .
\end{equation}
When the whole set of eigenvectors $\{|a\rangle\}$ of 
$\hat{H}$ is available, 
it is given by $|\Psi(t)\rangle = \sum_a c_a e^{-i E_a t} 
| a \rangle$ with $c_a = \langle a|\Psi(0)\rangle$. 
This is the exact method.

\begin{figure}[h]
\includegraphics*[width=0.5\columnwidth]{sfig1.eps}
\caption{Illustration of the XXZ Hamiltonian with the nearest and next
nearest interactions under the open boundary condition.}
\label{sfig1}
\end{figure}

When the dimension of the Hilbert space is too large, the exact method is
impractical and one needs to rely on an approximate method. 
We will explain an efficient method to simulate the time evolution.
The XXZ Hamiltonian of $N$ spins under the open boundary condition 
is given by~(see Fig.~\ref{sfig1})
\begin{equation}
\hat{H} = \sum_{l=1}^{N-1} \hat{h}_{l,l+1} + \sum_{l=1}^{N-2} \hat{h}_{l,l+2} 
\end{equation}
with the short-hand notation
\begin{equation}
\begin{aligned}
\hat{h}_{l,m} \equiv h(\hat{\sigma}_l,\hat{\sigma}_m) &= -\frac{J}{2} \left( \hat{\sigma}_l^x\hat{\sigma}_m^x + 
                                      \hat{\sigma}_l^y\hat{\sigma}_m^y + 
                        \Delta \hat{\sigma}_l^z \hat{\sigma}_m^z \right) \\
  & = -J \left( \hat{\sigma}_l^+ \hat{\sigma}_m^- + 
                \hat{\sigma}_l^- \hat{\sigma}_m^+ 
      + \frac{\Delta}{2} \hat{\sigma}_l^z \hat{\sigma}_m^z \right) ,
\end{aligned}
\end{equation}
where $\hat\sigma_l^{\pm} = (\hat\sigma_l^x \pm i \hat\sigma_l^y)/2$.

Note that the Hamiltonian is slightly different from that considered in the
paper. 
Generalizations to the Hamiltonian with bond-dependent $J$ and $\Delta$ 
are straightforward.

The time evolution operator satisfies $\hat{U}(\epsilon M) =
\hat{U}(\epsilon)^{M}$. Thus, the state vector at $t = \epsilon M$
is obtained by multiplying $\hat{U}(\epsilon)$ to the initial 
state vector $M$ times repeatedly. 
For infinitesimal $\epsilon$, $\hat{U}(\epsilon)$ can be approximated
efficiently.  First, we decompose the Hamiltonian as the sum 
$\hat{H} = \hat{H}_A + \hat{H}_B + \hat{H}_C + \hat{H}_D$ with
\begin{equation}\label{decomp}
\begin{split}
\hat{H}_A & =  \sum_{l} \hat{h}_{2l-1,2l} \\
\hat{H}_B & =  \sum_{l} \hat{h}_{2l,2l+1} \\
\hat{H}_C & =  \sum_{l} ( \hat{h}_{4l-3,4l-1} + \hat{h}_{4l-2,4l})\\
\hat{H}_D & =  \sum_{l} ( \hat{h}_{4l-1,4l+1} + \hat{h}_{4l,4l+2} ) .
\end{split}
\end{equation}
In Fig.~\ref{sfig1}, the interaction bonds contained in 
$\hat{H}_A$, $\hat{H}_B$, $\hat{H}_C$, and $\hat{H}_D$
are represented by the solid, dashed, dotted, and dash-dotted lines,
respectively.
Then, using the generalized Lie-Suzuki-Trotter formula of the second order
$e^{\epsilon(\hat{X}+\hat{Y})} = e^{\frac{\epsilon}{2}\hat{X}} e^{\epsilon
\hat{Y}} e^{\frac{\epsilon}{2}\hat{X}} +
O(\epsilon^3)$~[30], we obtain
\begin{equation}
\hat{U}(\epsilon) =  e^{-i\frac{\epsilon}{2}\hat{H}_A}
                     e^{-i\frac{\epsilon}{2}\hat{H}_B}
                     e^{-i\frac{\epsilon}{2}\hat{H}_C}
                     e^{-{i\epsilon}\hat{H}_D}
                     e^{-i\frac{\epsilon}{2}\hat{H}_C}
                     e^{-i\frac{\epsilon}{2}\hat{H}_B}
                     e^{-i\frac{\epsilon}{2}\hat{H}_A} + O(\epsilon^3)
\end{equation}

The important feature of the decomposition \eqref{decomp} is that all
$\hat{h}_{l,m}$'s contained in a partial Hamiltonian commute with each
others. It yields
\begin{equation}
e^{-i\gamma \hat{H}_A} = \prod_{l} e^{-i \gamma \hat{h}_{2l-1,2l}}
\end{equation}
and similar relations for the others. Consequently, $\hat{U}(\epsilon)$ can
be approximated as the ordered product of $\hat{U}_{l,m}(\gamma) = e^{-i \gamma
\hat{h}_{l,m}}$ with $\gamma=\epsilon$ or $\epsilon/2$ with an error $O(\epsilon^3)$.
If one chooses
$\{ |-\rangle_l\otimes |-\rangle_m, 
|-\rangle_l\otimes |+\rangle_m, 
|+\rangle_l\otimes |-\rangle_m,
|+\rangle_l\otimes |+\rangle_m \}$ as the basis states for the spins at
sites $l$ and $m$, $\hat{U}_{l,m}(\gamma) = \hat{u}_{l,m}(\gamma) \otimes
\mathsf{I}_{N-2}$ with the $4\times 4$ unitary matrix
\begin{equation}
\hat{u}_{l,m}(\gamma) = 
\begin{pmatrix}
e^{i \phi } & 0 & 0 & 0 \\
0 & e^{-i\phi} \cos\gamma J& - i e^{-i\phi} \sin \gamma J & 0 \\
0 & -i e^{-i\phi} \sin\gamma J&  e^{-i\phi} \cos \gamma J & 0  \\
0 & 0 & 0 & e^{i \phi} 
\end{pmatrix}
\end{equation}
with $\phi = \gamma J \Delta /2$.
The overall error of this method is of the order of $O(\epsilon^2)$.

Computationally, the infinitesimal time evolution is achieved by the
multiplication of a state vector with the $4\times 4$ matrices $O(N)$ times.
All the operators $\hat{U}_{l,m}$'s conserve the total number of up spins. 
Thus, one can work within the subspace consistent with the conservation,
which make numerical calculations efficient.
In the numerical calculations, we used the exact time evolution when $N\leq
16$, and the decomposition method for $N>16$ with $\epsilon = 0.01$. We
examined whether $\epsilon = 0.01$ is small enough. As can be seen in
Fig.~\ref{sfig2}, the data with $\epsilon=0.1, 0.05$, and $0.01$ are 
indistinguishable in the plot. More quantitatively, the values of 
$Q_{\rm L}$ at time $t=10^3$ are $6.004072339$, $6.003279566$, and 
$6.003025718$ for $\epsilon=$ $0.1$, $0.05$, and $0.01$, respectively. 
They are well fitted to the function $Q_{\rm L} = Q_0 + a \epsilon^{1.99926}$ with
$Q_0 = 6.003015128$. The relative error of the value at $\epsilon=0.01$ is
$\simeq 1.7\times 10^{-6}$. It proves that $\epsilon=0.01$ is small enough.

\begin{figure}
\includegraphics*[width=0.5\columnwidth]{sfig2.eps}
\caption{Numerical data for the energy density and the particle density for
the subsystem L using the same condition as Fig.~3 of the main text. The
data obtained by using $\epsilon=0.1, 0.05,$ and $0.01$ cannot be
distinguished.}
\label{sfig2}
\end{figure}

\subsection{Singular value decomposition for reduced density matrix}
Suppose that one separates the $N$ spins into two sets $A$ of $N_A$ spins 
and $B$ of $N_B = N-N_A$ spins. We assume that $N_A \leq N_B$ without loss
of generality.
Then, the state vector can be written as
\begin{equation}
|\Psi\rangle = \sum_{a,b} \psi_{a,b} |a\rangle\otimes |b\rangle ,
\end{equation}
where $\{|a\rangle\}$ and $\{|b\rangle\}$ span the Hilbert spaces of
dimension $d_A=2^{N_A}$ and $d_B = 2^{N_B}$ for spins in $A$ and $B$, 
respectively. The reduced density operator for the spins in $A$ is
given by 
\begin{equation}
\hat{\rho}_A = {\rm Tr}_B |\Psi\rangle\langle \Psi| = \sum_{b}
\psi_{a,b}\psi_{a',b}^* |a\rangle \langle a'| 
\end{equation}
and represented by the $({d_A}\times {d_A})$ reduced density matrix 
is given by
\begin{equation}
(\rho_A )_{a,a'} = \sum_b \psi_{a,b} \psi_{a',b}^* .
\end{equation}
The reduced density matrix is necessary when one calculates the local 
observables and the probability distribution for the spins in $A$. 
The eigenvalues are necessary
for the entanglement entropy. The singular value decomposition~(SVD) is a
convenient method to handle the reduced density matrix.

The probability amplitudes $\psi_{a,b}$ can be regarded as the 
elements of a $({d_A}\times {d_B})$ positive semidefinite complex matrix 
$M$, $M_{ab} =
\psi_{a,b}$. According to the SVD, such a matrix $M$ can be written as
\begin{equation}
M = V_1 \Sigma V_2^\dagger , 
\end{equation}
where $V_1$ is a $({d_A}\times {d_A})$ unitary matrix, $\Sigma$ is a $({d_A}
\times {d_B})$ rectangular diagonal matrix with non-negative real diagonal
elements, and $V_2$ is a $({d_B}\times {d_B})$ unitary matrix. 
The SVD can be done with a numerical analysis software. 

Once the SVD is done, the reduced density matrix is obtained from the
relation
\begin{equation}
\rho_A = M M^\dagger = V_1 \Sigma V_2^\dagger V_2 \Sigma^\dagger V_1^\dagger
       = V_1 \Sigma \Sigma^\dagger V_1^\dagger .
\end{equation}
That is, the reduced density matrix is given by the similarity transformation 
of the diagonal $(d_A\times d_A)$ matrix
\begin{equation}
\Sigma \Sigma^\dagger = {\rm
diag}(|\lambda_1|^2,|\lambda_2|^2,\cdots,|\lambda_{d_A}|^2) .
\end{equation}
The SVD allows one to calculate the entanglement entropy easily. 
It is given by
\begin{equation}
S_{ent} = -\sum_{n=1}^{d_A} |\lambda_n|^2 \ln |\lambda_n|^2 .
\end{equation}

\subsection{Temperature-energy relation in the $Q=1$ sector}
Consider the XXZ Hamiltonian of $N$ spins under the periodic boundary
condition. In the $Q=0$ sector where all spins are down~(or all sites are
empty), there is a single eigenstate with the energy 
\begin{equation}
E_{0,1} = -\frac{\Delta}{2} N .
\end{equation}
The Hilbert space of the $Q=1$ sector is spanned by the states $\{|1\rangle,
\cdots, |x\rangle,\cdots,|N\rangle\}$ where $|x\rangle$ represents the state
where the particle is at site $x$. The plane
wave $|k_n\rangle = \frac{1}{\sqrt{N}} \sum_x e^{ik_n x} |x\rangle$
is the eigenstate and the energy eigenvalue is given by
\begin{equation}
E_{1,k_n} = -\frac{\Delta}{2}(N-4) - \frac{2}{1+\lambda} \left( \cos k_n +
\lambda \cos 2k_n \right) = -\frac{\Delta}{2}(N-4) +
\varepsilon(k_n)
\end{equation}
with $k_n = 2\pi n / N$~($n=0,\cdots,N-1$). The partition function in the
$Q=1$ sector is given by
\begin{equation}
Z_1(\beta) = 2^{\beta \Delta (N-4)/2}\sum_{n=0}^{N-1}  e^{-\beta
\varepsilon(k_n)} 
\end{equation}
with the mean energy
\begin{equation}
\langle E\rangle = - \frac{\partial}{\partial \beta} \ln Z_1(\beta)
= -\frac{\Delta}{2}(N-4) + \frac{\sum_n \varepsilon(k_n) e^{-\beta
\varepsilon(k_n)}}{\sum_n e^{-\beta \varepsilon(k_n)}} .
\end{equation}

We estimate the location of the inversion temperature $T_{\rm I}$ 
by requiring that
the mean energy in the $Q=1$ sector should be equal to the energy in the $Q=0$
section. It yields that
\begin{equation}
-2 \Delta =  \frac{\sum_n \varepsilon(k_n) e^{-\beta_{\rm I}
\varepsilon(k_n)}}{\sum_n e^{-\beta_{\rm I} \varepsilon(k_n)}} 
\end{equation}
with $\beta_{\rm I} = 1/T_{\rm I}$. In the $N\to\infty$ limit, the summation is replaced
by the integration to yield the expression
\begin{equation}
-2\Delta = \frac{ \int_{-\pi}^\pi dk \ \varepsilon(k) e^{-\beta_{\rm I}
\varepsilon(k)}}{ \int_{-\pi}^\pi dk \ e^{-\beta_{\rm I} \varepsilon(k)}}  \equiv
\langle \varepsilon(k)\rangle_{\beta_{\rm I}}.
\end{equation}

Note that $\langle \varepsilon(k)\rangle_{\beta_{\rm I}}$ increases as
$\beta_{\rm I}
\to 0^+$ or $T_{\rm I} \to +\infty$. For small $\beta_{\rm I}$, we can approximate
\begin{equation}
\langle \varepsilon(k)\rangle_{\beta_{\rm I}} \simeq -\beta_{\rm I} \int \frac{dk}{2\pi}
\varepsilon(k)^2  = -2 \beta_{\rm I} \frac{1+\lambda^2}{(1+\lambda)^2} . 
\end{equation}
Thus, the inversion temperature has the limit form
\begin{equation}
T_{\rm I} \simeq \frac{1+\lambda^2}{(1+\lambda)^2} \Delta^{-1}
\end{equation}
as $\Delta \to 0^{+}$.
In the opposite limit where $\beta_{\rm I} \to \infty$~($T_{\rm I}\to 0^+$), $\langle
\varepsilon(k)\rangle_{\beta_{\rm I}}$ is dominated by the minimum value of
$\varepsilon(k)$. For positive $\lambda$, $\varepsilon(k)$ takes the minimum
value $-2$ at $k=0$. It yields that $T_{\rm I} \simeq (1-\Delta)/4$.

\subsection{Joule expansion}
In the main text, we present the Joule expansion of a gas into the
vacuum. We have performed the numerical simulation of the Joule expansion
in a mixture of a dense gas and a dilute gas.
The subsystems ${\rm L}$ and ${\rm R}$ is prepared in the energy eigenstates
$|Q_{\rm L}^0 =
q_{\rm L} N, n_{\rm L}\rangle$ and $|Q_{\rm R}^0 = q_{\rm R} N, n_{\rm R}\rangle $, respectively, 
with $q_{\rm L} = q = 1/4$ and $q_{\rm R} = 1-q = 3/4$ with the same quantum
number $n_{\rm L}
= n_{\rm R} = n_0$. 
The two states are tied to each other by the particle-hole symmetry.
We select those two states in order to keep them in the same initial
temperature.

When the total system reaches a stationary state, we measure the probability
distribution $P_{\rm L}(Q,n)$ that the subsystem ${\rm L}$ is in the state
$|Q,n\rangle$. Then, we fit the probability distribution to the function
$P_{\rm L}(Q,n) \propto \exp[ -(E_{Q,n}-\overline{E_{\rm L}})/T_f - a_{11}
(E_{Q,n}-\overline{E_{\rm L}})^2]$ at the most probable sector $Q=N/2$ to find the
final temperature $T_f$. The numerical result is presented in
Fig.~\ref{sfig3}. We also observe that the gas heats up (cools down) when
the initial temperature is lower (higher) than a threshold temperature.
The finite size effect is also seen in the plot.

\begin{figure}
\includegraphics*[width=0.5\columnwidth]{sfig3.eps}
\caption{Heating/cooling curve when the dense gas of density $\rho_{\rm R}=3/4$ is
mixed with the dilute gas of density $\rho_{\rm L} = 1/4$. The filled symbols
represent the numerical results obtained by fitting the probability 
distribution. The open symbols represent the expected result from 
the isotherm curve.}
\label{sfig3}
\end{figure}

\end{document}

