Why Do Thermal Operations Form a Semigroup? Given a fixed "background temperature" $T$ and an $n$-dimensional system, the set of thermal operations is usually defined to be
$$
\bigcup_{m\in\mathbb N}\Big\{ \operatorname{tr}_B\Big(U\Big((\cdot)\otimes \frac{e^{-H_B/kT}}{\operatorname{tr}(e^{-H_B/kT})}\Big)U^*\Big) :\ \substack{H_B\in\mathbb C^{m\times m}\text{ hermitian}\\ U\in\mathbb C^{nm\times nm}\text{ unitary, }[U,H_S\otimes\textbf{1}_B+\textbf{1}_S\otimes H_B]=0 }\Big\}
$$
cf. Sec. 1.3 in Lostaglio, Rep. Prog. Phys. 82 114001 (2019) (arXiv version). In other words---as Lostaglio writes on the bottom of page 6---thermal operations are compositions of "preparing thermal states at a fixed temperature $T$ and with arbitrary Hamiltonian $H_B$, performing energy-preserving unitaries, and tracing out subsystems" where the dimension of the bath is arbitrary (but finite).
In the same review article Lostaglio recalls that a resource theory is defined by---among other things--a subset $\mathcal A$ of all quantum channels ("allowed operations") which is closed under composition and includes the identity; in other words $\mathcal A$ has to be a semigroup with identity which does make sense from an operational point of view. Thus the question I had for a while was the following:

Is the set of thermal operations---as defined above---closed under composition? And if so, why?

On the one hand, this whole construction is known under the name "resource theory of thermal operations" so the implied answer is yes, yet I have seen neither this explicit statement nor any computation justifying this result anywhere until now. It turns out that the computation is rather elementary---although non-trivial---so I would like to answer my own question by providing a short proof below.
 A: Consider two elements of the set of thermal operations, with unitaries $U$ and $U'$ and Hamiltonians $H_B$ and $H_{B'}$, respectively. Following your convention, I will write $\gamma_\bullet=\frac{e^{-H_\bullet/kT}}{\operatorname{tr}(e^{-H_\bullet/kT})}$.
Then, the concatenation of both thermal operators is the map
$$
\rho\to \mathrm{tr}\big[U_{SB}U'_{SB'}(\rho_S\otimes\gamma_B\otimes\gamma_{B'})(U_{SB}U'_{SB'})^\dagger\big]\ .
$$
That is, the map is of the form of a thermal operation
$$
\rho\to \mathrm{tr}\big[V(\rho_S\otimes\gamma_{BB'}V^\dagger\big]\ ,
$$
where we define $\gamma_{BB'}:=\gamma_{B}\otimes\gamma_{B'}$ and $V:=U_{SB}U'_{SB'}$. It remains to show that $V$ and $\gamma_{BB'}$ have the necessary properties:

*

*$\gamma_{BB'}%=\gamma_{B}\otimes\gamma_{B'} = \frac{e^{-H_B/kT}}{\operatorname{tr}(e^{-H_B/kT})} \otimes \frac{e^{-H_{B'}/kT}}{\operatorname{tr}(e^{-H_{B'}/kT})}=\frac{e^{-H_{BB'}/kT}}{\operatorname{tr}(e^{-H_{BB'}/kT})}$ is the thermal state of $H_{BB'}=H_B\otimes I_{B'}+I_B\otimes H_{B'}$.

*$V=U_{SB}U'_{SB'}$ and $H_S+H_{BB'}=H_S + H_{B}+H_{B'}$
(I omit the identities) commute: $U'_{SB'}$ commutes with $H_S+H_{B'}$ and with $H_B$ (the latter trivially), and vice versa.

A: All we have to do is take two thermal operations and show that their composition again has the form of a thermal operation. Thus given two bath Hamiltonians $H_B\in\mathbb C^{m\times m}$, $H_{B'}\in\mathbb C^{m'\times m'}$ and unitaries $U\in\mathbb C^{nm\times nm}$, $V\in\mathbb C^{nm'\times nm'}$ which satisfy
$$
[U,H_S\otimes\textbf{1}_B+\textbf{1}_S\otimes H_B]=0\quad\text{ and }\quad [V,\textbf{1}_{B'}\otimes H_S+H_{B'}\otimes\textbf{1}_S]=0\tag{1}
$$
we consider the map
$$
\rho\mapsto \operatorname{tr}_{B'}(V(\gamma_{B'}\otimes \boxed{\operatorname{tr}_B(U(\rho\otimes\gamma_B)U^*)})V^*)\tag{2}
$$
Here we used the isomorphism $S\otimes B\otimes B'\simeq B'\otimes S\otimes B$ to shift the second bath to the first tensor factor; this will simplify the computation but of course does not affect the result itself. Also $\gamma_B$ ($\gamma_{B'}$) is short for the thermal state on $B$ ($B'$) with respect to the Hamiltonian $H_B$ ($H_{B'}$) and the fixed temperature $T$.
A common trick to show that two states (or two channels) coïncide is to show that all their expectation values coïncide; this is especially useful here as the partial trace is characterized by the duality relation $\operatorname{tr}(X\operatorname{tr}_B(A))=\operatorname{tr}((X\otimes\textbf{1}_B)A)$. Thus let arbitrary $X\in\mathbb C^{n\times n}$ be given. We compute
$$
\operatorname{tr}\big(X\operatorname{tr}_{B'}(V(\gamma_{B'}\otimes \operatorname{tr}_B(U(\rho\otimes\gamma_B)U^*))V^*)\big)=\operatorname{tr}\big(V^*(\textbf{1}_{B'}\otimes X)V\big(\gamma_{B'}\otimes \operatorname{tr}_B(U(\rho\otimes\gamma_B)U^*)\big)\big)\,.
$$
The non-trivial (yet somehow intuitive) result we now need is that the partial trace over $B$ "commutes with" attaching the bath $B'$ if the $B'$-part of the system is "in a product state". More formally we use that for all $A\in\mathcal B^1(\mathcal H_1)$ and all $C\in\mathcal B^1(\mathcal H_2\otimes\mathcal H_3)$ one has $\operatorname{tr}_3 (A\otimes C)=A\otimes\operatorname{tr}_3(C)$ (where $\operatorname{tr}_3$ is the partial trace over $\mathcal H_3$), cf. also the footnote at the bottom. This yields
\begin{align*}
\operatorname{tr}\big(X\operatorname{tr}_{B'}(...)\big)&=\operatorname{tr}\big(V^*(\textbf{1}_{B'}\otimes X)V \operatorname{tr}_B(\gamma_{B'}\otimes U(\rho\otimes\gamma_B)U^*))\big)\\
&=\operatorname{tr}\big(\big((V^*(\textbf{1}_{B'}\otimes X)V)\otimes\mathbf{1}_B\big) (\gamma_{B'}\otimes U(\rho\otimes\gamma_B)U^*))\big)\\
&=\operatorname{tr}\big((V^*\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes X\otimes \textbf{1}_B)(V\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes U) (\gamma_{B'}\otimes \rho\otimes\gamma_B)(\textbf{1}_{B'}\otimes U^*)\big)\\
&=\operatorname{tr}\Big(X\operatorname{tr}_{B,B'}\big((V\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes U) (\gamma_{B'}\otimes \rho\otimes\gamma_B)(\textbf{1}_{B'}\otimes U^*)(V^*\otimes \textbf{1}_B)\big)\Big)
\end{align*}
Finally, because $X$ was chosen arbitrarily, we see that the channel from eq.(2) is the same as
$$
\rho\mapsto\operatorname{tr}_{B,B'}\big((V\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes U) (\gamma_{B'}\otimes \rho\otimes\gamma_B)(\textbf{1}_{B'}\otimes U^*)(V^*\otimes \textbf{1}_B)\big)\,.
$$
Most importantly this is a thermal operation generated by the "bath Hamiltonian" $H_B\otimes\textbf{1}_{B'}+\textbf{1}_B\otimes H_{B'}\in\mathbb C^{mm'\times mm'}$ (so the corresponding thermal state is $\gamma_B\otimes\gamma_{B'}$; direct computation) and the unitary $(V\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes U)\in\mathbb C^{nmm'\times nmm'}$. Finally, a straightforward computation shows that the "large" unitary is energy-conserving:
\begin{align*}
(V\otimes \textbf{1}_B)&(\textbf{1}_{B'}\otimes U)(\mathbf{1}_{B'}\otimes H_S\otimes\mathbf{1}_B+\textbf{1}_{B'}\otimes\mathbf{1}_S\otimes H_B+H_{B'}\otimes\textbf{1}_S\otimes\textbf{1}_B)\\
&=(V\otimes \textbf{1}_B)\big( \mathbf{1}_{B'}\otimes(U(H_S\otimes\mathbf{1}_{B}+\mathbf{1}_{S}\otimes H_B))+H_{B'}\otimes U\big)\\
&\overset{(1)}=(V\otimes \textbf{1}_B)\big( \mathbf{1}_{B'}\otimes((H_S\otimes\mathbf{1}_{B}+\mathbf{1}_{S}\otimes H_B)U)+(H_{B'}\otimes\textbf{1}_S\otimes\textbf{1}_B)(\textbf{1}_{B'}\otimes U)\big)\\
&= (V\otimes \textbf{1}_B)(\mathbf{1}_{B'}\otimes H_S\otimes\mathbf{1}_B+\textbf{1}_{B'}\otimes\mathbf{1}_S\otimes H_B+H_{B'}\otimes\textbf{1}_S\otimes\textbf{1}_B)(\textbf{1}_{B'}\otimes U) \\
=\ldots&= (\mathbf{1}_{B'}\otimes H_S\otimes\mathbf{1}_B+\textbf{1}_{B'}\otimes\mathbf{1}_S\otimes H_B+H_{B'}\otimes\textbf{1}_S\otimes\textbf{1}_B)(V\otimes \textbf{1}_B)(\textbf{1}_{B'}\otimes U)
\end{align*}
This shows that the thermal operations indeed form a semigroup which concludes the proof. $\square$

 Footnote: For anyone who wants to see why this partial trace identity holds: Recall that the partial trace over a space with orthonormal basis $\langle g_n\rangle_{n\in\mathcal N}$ satisfies $\langle \psi|\operatorname{tr}_B(A)|\phi\rangle=\sum_{n\in\mathcal N}\langle \psi\otimes g_n|A|\phi\otimes g_n\rangle$. This lets us compute
\begin{align*}
\langle x\otimes y, \operatorname{tr}_3 (A\otimes C) (z\otimes w)\rangle&= \sum_{n\in\mathcal N}\langle x\otimes y\otimes g_n,(A\otimes C)(z\otimes w\otimes g_n)\rangle\\
&=\sum_{n\in\mathcal N}\langle x,Az\rangle\langle y\otimes g_n,C(w\otimes g_n)\rangle=\langle x,Az\rangle\langle y,\operatorname{tr}_3(C)w\rangle =\langle x\otimes y,( A\otimes\operatorname{tr}_3(C)) (z\otimes w)\rangle
\end{align*}
for all $x,z\in\mathcal H_1$, $y,w\in\mathcal H_2$. By linearity this extends to all expectation values and thus the two operators coïncide.
