Quantum collisional thermostats

We consider units drawn from a reservoir and colliding, one by one, with the system [18]. Each unit U is a particle of mass m with internal states and moving in one dimension. Its corresponding Hilbert space is ${\mathcal{H}}_{U}={\mathcal{H}}_{U,\mathrm{p}}\otimes {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}$, where ${\mathcal{H}}_{U,\mathrm{p}}$ refers to the spatial states and ${\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}$ is the space of internal states. The units collide with the system S, which only has internal degrees of freedom and whose states are vectors in the Hilbert space ${\mathcal{H}}_{S}$.

The system is a fixed scatterer located in an interval [−L/2, L/2], as sketched in figure 1. The whole setup is described by the following Hamiltonian H_tot, which is an operator acting on ${\mathcal{H}}_{U}\otimes {\mathcal{H}}_{S}$:

$$\begin{equation}{H}_{\text{tot}}=\frac{{\hat{p}}^{2}}{2m}+{\chi }_{L}(\hat{x}){H}_{US}+{H}_{U}+{H}_{S}.\end{equation} \tag{ 1 }$$

$\hat{p}$ and $\hat{x}$ being, respectively, the momentum and position operators in ${\mathcal{H}}_{U,\mathrm{p}}$. χ_L (x) is the indicator function of the scattering region [−L/2, L/2]: χ_L (x) = 1 if x ∈ [−L/2, L/2] and zero otherwise. Outside the scattering region, the free Hamiltonian H₀ = H_U + H_S rules the evolution of the internal degrees of freedom and is the sum of the Hamiltonian of the system H_S and of the internal degrees of freedom of the unit H_U . Within the scattering region, the Hamiltonian affecting the internal degrees of freedom is H = H_US + H₀, which we will call total internal Hamiltonian. The free and the total internal Hamiltonians, H₀ and H respectively, are operators in ${\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$. We will assume that both have a discrete spectrum with eigenstates:

$$\begin{align}\hfill {H}_{0}\vert {s}_{J}\rangle & ={e}_{J}\vert {s}_{J}\rangle \hfill \\ \hfill H\vert {s}_{J}^{\prime }\rangle & ={e}_{J}^{\prime }\vert {s}_{J}^{\prime }\rangle .\hfill \end{align} \tag{ 2 }$$

Notice that {|s_J ⟩} and $\left\{\vert {s}_{J}^{\prime }\rangle \right\}$ are orthonormal basis of the Hilbert space of the internal states of the unit and the system, ${\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$. Moreover, the eigenvectors of H₀ can be written as $\vert {s}_{J}\rangle =\vert {s}_{{j}_{U}}{\rangle }_{U}\otimes \vert {s}_{{j}_{S}}{\rangle }_{S}\in {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$, with

$$\begin{align}\hfill {H}_{U}\vert {s}_{{j}_{U}}{\rangle }_{U}& ={e}_{{j}_{U}}^{(U)}\vert {s}_{{j}_{U}}{\rangle }_{U}\hfill \\ \hfill {H}_{S}\vert {s}_{{j}_{S}}{\rangle }_{S}& ={e}_{{j}_{S}}^{(S)}\vert {s}_{{j}_{S}}{\rangle }_{S}\hfill \end{align} \tag{ 3 }$$

and total energy ${e}_{J}={e}_{{j}_{U}}^{(U)}+{e}_{{j}_{S}}^{(S)}$. Here and in the rest of the paper, we use capital letters J for the quantum numbers labelling the eigenstates of H₀ and H, and lower case letters, j_U , j_S , for the quantum numbers corresponding to H_U and H_S . In this notation, the quantum number J of an eigenstate of H₀ comprises the two quantum numbers J = (j_S , j_U ).

Zoom In Zoom Out Reset image size

In order to observe well-defined collisions, the spatial state of the units must be a wave packet $\vert {\phi }_{{p}_{0},{x}_{0}}\rangle$ centered around position x₀ and momentum p₀, with momentum dispersion σ_p , as the one depicted in figure 1. An example is the Gaussian wave packet, whose wave function in momentum representation reads [18]

$$\begin{equation}\langle p\vert {\phi }_{{p}_{0},{x}_{0}}\rangle ={(2\pi {\sigma }_{p}^{2})}^{-1/4}\enspace \mathrm{exp}\left[-\frac{{(p-{p}_{0})}^{2}}{4{\sigma }_{p}^{2}}-\mathrm{i}\frac{p{x}_{0}}{\hslash }\right].\end{equation} \tag{ 4 }$$

Here |p⟩ denotes the non-normalizable plane wave with momentum p. The Hamiltonian (1) is invariant under spatial reflection, $(\hat{x},\hat{p})\to (-\hat{x},-\hat{p})$, since the indicator function of the interval [−L/2, L/2] is even: χ_L (x) = χ_L (−x). Consequently, a collision with a unit coming from the left with positive velocity, x₀ < 0 and p₀ > 0, is equivalent to the mirror collision with a unit $\vert {\phi }_{-{p}_{0},-{x}_{0}}\rangle$ coming from the right. Hence, we can limit our discussion to units with positive velocity, without loss of generality (the mathematical consequences of this spatial symmetry in the scattering problem are explained in detail in appendix B2).

The internal state of the unit is disentangled from the system and depends on the properties of the reservoir. For instance, if the reservoir is in thermal equilibrium at inverse temperature β, the internal state is the Gibbs state ${\rho }_{U,\mathrm{e}\mathrm{q}}={\text{e}}^{-\beta {H}_{U}}/{Z}_{U}$, where Z_U is the corresponding partition function. We first analyze the case of a unit in a pure eigenstate of H_U , $\vert {s}_{{j}_{U}}{\rangle }_{U}$, and later on we consider thermal mixtures of these eigenstates.

The effect of the collision on the system is given by a CPTP map, which depends on the incident momentum p₀ and the internal state of the unit. However, it is convenient to consider first the effect of the collision on all the internal degrees of freedom, those of the system and of the internal state of the unit, following reference [18]. The scattering map $\mathbb{S}$ relates the internal state before the collision, ρ, and after, ${\rho }^{\prime }=\mathbb{S}\rho$. Expressing the states in the eigenbasis of H₀, ρ_JK = ⟨s_J |ρ|s_K ⟩ and ${\rho }_{JK}^{\prime }=\langle {s}_{J}\vert {\rho }^{\prime }\vert {s}_{K}\rangle$, the scattering map is given by a tensor ${\mathbb{S}}_{{J}^{\prime }{K}^{\prime }}^{JK}$ such that

$$\begin{equation}{\rho }_{{J}^{\prime }{K}^{\prime }}^{\prime }=\sum\limits _{J,K}\enspace {\mathbb{S}}_{{J}^{\prime }{K}^{\prime }}^{JK}{\rho }_{JK}.\end{equation} \tag{ 5 }$$

In reference [18], we have analyzed in detail this scattering map and found that the behavior of the state ρ crucially depends on the momentum dispersion of the packet, σ_p . If the dispersion is small enough, the outgoing wave packets corresponding to different transitions |s_J ⟩ → |s_J'⟩ are either identical or do not overlap. The precise condition for these narrow wave packets, in terms of the transition energies Δ_JJ' = e_J − e_J', reads

$$\begin{equation}{\sigma }_{p}\ll \frac{m\vert {{\Delta}}_{J{J}^{\prime }}-{{\Delta}}_{K{K}^{\prime }}\vert }{2{p}_{0}}\end{equation} \tag{ 6 }$$

for every pair of transitions |s_J ⟩ → |s_J'⟩ and |s_K ⟩ → |s_K'⟩ with Δ_J'J ≠ Δ_KK'. If the incident packet fulfills this condition, we call it narrow wave packet and the scattering map induced by the collision is given by [18]:

$$\begin{align}\hfill {\mathbb{S}}_{{J}^{\prime }{K}^{\prime }}^{JK}& \simeq {t}_{{J}^{\prime }J}({E}_{{p}_{0}}+{e}_{J}){\left[{t}_{{K}^{\prime }K}({E}_{{p}_{0}}+{e}_{K})\right]}^{\ast }\hfill \\ \hfill & \quad +{r}_{{J}^{\prime }J}({E}_{{p}_{0}}+{e}_{J}){\left[{r}_{{K}^{\prime }K}({E}_{{p}_{0}}+{e}_{K})\right]}^{\ast },\hfill \end{align} \tag{ 7 }$$

whenever

$$\begin{equation}{e}_{{J}^{\prime }}-{e}_{J}={e}_{{K}^{\prime }}-{e}_{K}\end{equation} \tag{ 8 }$$

and ${E}_{{p}_{0}}+{e}_{J}\geqslant {e}_{{J}^{\prime }}$, and zero otherwise. Here t_J'J (E) and r_J'J (E) are the transmission and reflection amplitudes that depend on the total energy, kinetic ${E}_{{p}_{0}}\equiv {p}_{0}^{2}/(2m)$ plus internal e_J . They are defined for all J and J' such that e_J , e_J' ⩽ E, which are the so-called open channels in the collision and span the Hilbert subspace

$$\begin{equation}{\mathcal{H}}_{\text{open}}=\mathrm{l}\mathrm{i}\mathrm{n}\left\{\vert {s}_{J}\rangle :{e}_{J}\leqslant E\right\}\subseteq {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}.\end{equation} \tag{ 9 }$$

For Hamiltonians that are invariant under spatial reflection, $(\hat{x},\hat{p})\to (-\hat{x},-\hat{p})$, the transmission and reflection amplitudes are usually arranged into two matrices t(E) and r(E) that form the scattering matrix (see appendix B2):

$$\begin{equation}\tilde{\mathcal{S}}(E)=\left(\begin{matrix}\hfill \mathbf{r}(E)\hfill & \hfill \mathbf{t}(E)\hfill \\ \hfill \mathbf{t}(E)\hfill & \hfill \mathbf{r}(E)\hfill \end{matrix}\right).\end{equation} \tag{ 10 }$$

The two matrices t(E) and r(E) are defined on the subspace of open channels ${\mathcal{H}}_{\text{open}}$. One important property of the scattering matrix is that it is unitary on the subspace ${\mathcal{H}}_{\text{open}}$ for a given total energy E, that is, ${\tilde{\mathcal{S}}}^{{\dagger}}(E)\tilde{\mathcal{S}}(E)=\mathbb{I}$ for all E, implying

$$\begin{align}\hfill \mathbf{r}(E){\mathbf{r}}^{{\dagger}}(E)+\mathbf{t}(E){\mathbf{t}}^{{\dagger}}(E)& =\mathbb{I}\hfill \\ \hfill \mathbf{r}(E){\mathbf{t}}^{{\dagger}}(E)+\mathbf{t}(E){\mathbf{r}}^{{\dagger}}(E)& =0.\hfill \end{align} \tag{ 11 }$$

The scattering map equation (7) determines the effect of a single collision on the system. If we now bombard the system with a stream of units, the evolution will be given by successive applications of the scattering map followed by the free evolution ruled by the Hamiltonian H_S [18]. If we neglect the free evolution, the behavior of the diagonal terms of the density matrix ρ_JJ , which are the populations of the energy levels ${e}_{J}={e}_{{j}_{U}}^{(U)}+{e}_{{j}_{S}}^{(S)}$, is determined by the coefficients ${\mathbb{S}}_{{J}^{\prime }{J}^{\prime }}^{JK}$ of the scattering map. Condition (8), particularized to J' = K', indicates that these coefficients are different from zero only if e_J = e_K . Moreover, since the initial internal state of the unit is an eigenstate of H_U , j_U = k_U ; hence, ${e}_{{j}_{S}}^{(S)}={e}_{{k}_{S}}^{(S)}$. If the Hamiltonian of the system H_S is non degenerate, this implies j_S = k_S and populations evolve independently of the off-diagonal terms of the density matrix

$$\begin{equation}{\rho }_{{J}^{\prime }{J}^{\prime }}^{\prime }=\sum\limits _{J}\enspace {P}_{{J}^{\prime }J}({p}_{0}){\rho }_{JJ}\end{equation} \tag{ 12 }$$

with the following transition probabilities that depend on the momentum p₀ of the incident unit:

$$\begin{equation}{P}_{{J}^{\prime }J}({p}_{0})\equiv {\mathbb{S}}_{{J}^{\prime }{J}^{\prime }}^{JJ}=\vert {t}_{{J}^{\prime }J}({E}_{{p}_{0}}+{e}_{J}){\vert }^{2}+\vert {r}_{{J}^{\prime }J}({E}_{{p}_{0}}+{e}_{J}){\vert }^{2}\end{equation} \tag{ 13 }$$

if ${p}_{0}^{2}\geqslant 2m{{\Delta}}_{{J}^{\prime }J}$ and zero otherwise. On the other hand, the unitarity of the scattering matrix on the subspace ${\mathcal{H}}_{\text{open}}$ of open channels, equation (11), implies that the off-diagonal terms decay [18], since |t_J'J |² + |r_J'J |² ⩽ 1, and that the trace of the density matrix is preserved, ∑_J' P_J'J (p₀) = 1 for all p₀. From now on, we will focus on the effect of narrow wave packets and only discuss the behavior of the populations, assuming the off-diagonal terms of the density matrix rapidly decay due to the collisions.

In this subsection we explore whether the system thermalizes if the units are in equilibrium at inverse temperature β. This implies that the units are in an internal state $\vert {s}_{{j}_{U}}{\rangle }_{U}$ with probability

$$\begin{equation}{p}_{{j}_{U}}=\frac{{\text{e}}^{-\beta {e}_{{j}_{U}}^{(U)}}}{{Z}_{U}}\end{equation} \tag{ 14 }$$

where Z_U is the internal partition function of the unit. The momentum of the units coming from a thermal bath is distributed as [18, 19]:

$$\begin{equation}\mu (p)=\frac{\beta \vert p\vert }{m}\enspace {\text{e}}^{-\beta {p}^{2}/(2m)}\quad p\in [0,\infty ].\end{equation} \tag{ 15 }$$

This is the effusion distribution describing the momentum of particles in equilibrium that cross a given point or hit a fixed scatterer coming from the left (since our scatterer, as described by the total Hamiltonian (1), is symmetric, there is no need to explicitly consider the case of negative incident velocity). We have shown in reference [19] how this effusion distribution arises from the Maxwellian velocity distribution and a uniform density of classical particles, which characterize an ideal gas at equilibrium.

In this case, the populations p(J) ≡ ρ_JJ , obey the following evolution equation

$$\begin{equation}{p}^{\prime }({J}^{\prime })=\sum\limits _{J}p(J)p(J\to {J}^{\prime })\end{equation} \tag{ 16 }$$

with

$$\begin{equation}p(J\to {J}^{\prime })={\int }_{0}^{\infty }\mathrm{d}{p}_{0}\enspace \mu ({p}_{0}){P}_{{J}^{\prime }J}({p}_{0}).\end{equation} \tag{ 17 }$$

The evolution equation for the state of the system

$$\begin{equation}p({j}_{S})=\sum\limits _{{j}_{U}}p({j}_{S},{j}_{U}),\end{equation} \tag{ 18 }$$

with p(j_S , j_U ) ≡ p(J), reads

$$\begin{equation}{p}^{\prime }({j}_{S}^{\prime })=\sum\limits _{{j}_{S}}p({j}_{S})p({j}_{S}\to {j}_{S}^{\prime })\end{equation} \tag{ 19 }$$

with

$$\begin{equation}p({j}_{S}\to {j}_{S}^{\prime })=\sum\limits _{{j}_{U},{j}_{U}^{\prime }}\frac{{\text{e}}^{-\beta {e}_{{j}_{U}}^{(U)}}}{{Z}_{U}}p(J\to {J}^{\prime })\end{equation} \tag{ 20 }$$

where we recall that J denotes the pair of quantum numbers (j_S , j_U ).

A sufficient condition for thermalization is micro-reversibility or invariance of the collision probabilities under time reversal [18]. In a quantum system, the states are transformed under time reversal by means of an anti-unitary operator $\mathsf{T}$ defined on the corresponding Hilbert space. Any anti-unitary operator can be written as $\mathsf{T}=CU$, where C is the conjugation of coordinates in a given basis and U is a unitary operator [22]. The time reversal operator depends on the physical nature of the system. Consider for instance a qubit with Hilbert space $\mathcal{H}={\mathbb{C}}^{2}$. If the qubit is a 1/2 spin, then time-reversal must change the sign of all the components of the spin, i.e. $\mathsf{T}{\sigma }_{\alpha }{\mathsf{T}}^{{\dagger}}=-{\sigma }_{\alpha }$ for α = x, y, z, where σ_α are the Pauli matrices. The anti-unitary operator that fulfills these transformations is $\mathsf{T}=C{\sigma }_{y}$, where C is the conjugation of the coordinates of the qubit in the canonical basis (the eigenbasis of σ_z ) [22]. On the other hand, if the qubit is a two-level atom whose states are superpositions of real wave functions in the position representation, then the time reversal operator is just $\mathsf{T}=C$, since the time-reversal of spinless particles is the conjugation of the wave function in the position representation.

In our case, the total time-reversal operator acting on the Hilbert space ${\mathcal{H}}_{U,\mathrm{p}}\otimes {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$ can be decomposed into three parts $\mathsf{T}={\mathsf{T}}_{U,\mathrm{p}}\otimes {\mathsf{T}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathsf{T}}_{S}$. The operator ${\mathsf{T}}_{U,\mathrm{p}}$ is the conjugation of the spatial wave function of the unit in the position representation ${\mathsf{T}}_{U,\mathrm{p}}\psi (x)={\psi }^{\ast }(x)$, whereas in momentum representation reads ${\mathsf{T}}_{U,\mathrm{p}}\phi (p)={\phi }^{\ast }(-p)$ [22, 23]. The time-reversal operator for the internal degrees of freedom can in principle be any anti-unitary operator ${\mathsf{T}}_{\text{int}}={\mathsf{T}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathsf{T}}_{S}$.

Micro-reversibility occurs when the total Hamiltonian commutes with the time-reversal operator, $[{H}_{\text{tot}},\mathsf{T}]=0$. Since the kinetic part is already invariant under time reversal, the commutation $[H,{\mathsf{T}}_{\text{int}}]=0$ is a sufficient condition for micro-reversibility. For simplicity, we further assume that $[{H}_{0},{\mathsf{T}}_{\text{int}}]=0$ and that the eigenstates of H₀ are time-reversal invariant: ${\mathsf{T}}_{\text{int}}\vert {s}_{J}\rangle =\vert {s}_{J}\rangle$ for all J. In appendix B we show that, if these conditions are fulfilled, then the scattering matrix obeys ${\tilde{\mathcal{S}}}^{\ast }\tilde{\mathcal{S}}=\mathbb{I}$. Combining this expression with the unitarity of $\tilde{\mathcal{S}}$, we conclude that the matrices t and r are symmetric for a given energy E:

$$\begin{align}\hfill \langle {s}_{{J}^{\prime }}\vert \mathbf{t}(E)\vert {s}_{J}\rangle & =\langle {s}_{J}\vert \mathbf{t}(E)\vert {s}_{{J}^{\prime }}\rangle \hfill \\ \hfill \langle {s}_{{J}^{\prime }}\vert \mathbf{r}(E)\vert {s}_{J}\rangle & =\langle {s}_{J}\vert \mathbf{r}(E)\vert {s}_{{J}^{\prime }}\rangle .\hfill \end{align} \tag{ 21 }$$

If we now apply this symmetry to the transition probabilities given by equation (13), we obtain

$$\begin{equation}{P}_{{J}^{\prime }J}({p}_{0})={P}_{J{J}^{\prime }}\left(\sqrt{{p}_{0}^{2}-2m{{\Delta}}_{{J}^{\prime }J}}\right)\end{equation} \tag{ 22 }$$

for all p₀ satisfying ${p}_{0}^{2}\geqslant 2m{{\Delta}}_{{J}^{\prime }J}$.

Let us prove now that micro-reversibility, as expressed by equation (22) for the transition probabilities, is a sufficient condition for thermalization. We first focus on transitions of internal states including the unit, that is, from |s_J ⟩ to |s_J'⟩. If e_J' ⩾ e_J , then Δ_J'J ⩾ 0 and the transition probability reads

$$\begin{equation}p(J\to {J}^{\prime })={\int }_{\sqrt{2m{{\Delta}}_{{J}^{\prime }J}}}^{\infty }\mathrm{d}{p}_{0}\enspace \frac{\beta {p}_{0}}{m}\enspace {\text{e}}^{-\beta {p}_{0}^{2}/(2m)}{P}_{{J}^{\prime }J}({p}_{0}).\end{equation} \tag{ 23 }$$

Here, the lower limit in the integral is due to the fact that P_J'J (p₀) is zero for ${p}_{0}^{2}\leqslant 2m{{\Delta}}_{{J}^{\prime }J}$. If we change the integration variable to ${p}_{0}^{\prime }=\sqrt{{p}_{0}^{2}-2m{{\Delta}}_{{J}^{\prime }J}}{\Rightarrow}\mathrm{d}{p}_{0}^{\prime }=\vert {p}_{0}\vert \mathrm{d}{p}_{0}/\vert {p}_{0}^{\prime }\vert$, we obtain

$$\begin{equation}p(J\to {J}^{\prime })={\int }_{0}^{\infty }\mathrm{d}{p}_{0}^{\prime }\enspace \frac{\beta {p}_{0}^{\prime }}{m}\enspace {\text{e}}^{-\beta ({p}_{0}^{\prime 2}/(2m)+{{\Delta}}_{{J}^{\prime }J})}\enspace {P}_{{J}^{\prime }J}\left(\sqrt{{p}_{0}^{\prime 2}+2m{{\Delta}}_{{J}^{\prime }J}}\right).\end{equation} \tag{ 24 }$$

Finally, applying the micro-reversibility condition (22),

$$\begin{align}\hfill p(J\to {J}^{\prime })& ={\text{e}}^{-\beta {{\Delta}}_{{J}^{\prime }J}}{\int }_{0}^{\infty }\mathrm{d}{p}_{0}^{\prime }\enspace \frac{\beta {p}_{0}^{\prime }}{m}\enspace {\text{e}}^{-\beta {p}_{0}^{\prime 2}/(2m)}{P}_{J{J}^{\prime }}({p}_{0}^{\prime })\hfill \\ \hfill & ={\text{e}}^{-\beta {{\Delta}}_{{J}^{\prime }J}}p({J}^{\prime }\to J).\hfill \end{align} \tag{ 25 }$$

We can proceed in an analogous way for the case e_J' ⩽ e_J . The final result is the local detailed balance condition

$$\begin{equation}\frac{p(J\to {J}^{\prime })}{p({J}^{\prime }\to J)}={\text{e}}^{-\beta ({e}_{{J}^{\prime }}-{e}_{J})}\quad \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\enspace J,{J}^{\prime }.\end{equation} \tag{ 26 }$$

We now explicitly consider the internal states of the unit. Recall that the subindex J in the previous sections comprises two quantum numbers J = (j_S , j_U ). If the internal states of the unit are in thermal equilibrium at inverse temperature β, then the transition probabilities between the states of the system are given by (20). The detailed balance condition (26) can be written as

$$\begin{equation}p(J\to {J}^{\prime })={\text{e}}^{-\beta \left[{e}_{{j}_{U}^{\prime }}^{(U)}+{e}_{{j}_{S}^{\prime }}^{(S)}-{e}_{{j}_{U}}^{(U)}-{e}_{{j}_{S}}^{(S)}\right]}p({J}^{\prime }\to J).\end{equation} \tag{ 27 }$$

Inserting (27) into (20), one gets

$$\begin{equation}p({j}_{S}\to {j}_{S}^{\prime })={\text{e}}^{-\beta \left[{e}_{{j}_{S}^{\prime }}^{(S)}-{e}_{{j}_{S}}^{(S)}\right]}p({j}_{S}^{\prime }\to {j}_{S})\end{equation} \tag{ 28 }$$

which is the detailed balance condition for the populations of the states of the system and ensures thermalization.

The standard procedure to obtain the scattering matrix $\tilde{S}(E)$ for a given energy E consists in solving the time-independent Schrödinger equation

$$\begin{equation}{H}_{\text{tot}}\vert \psi \rangle =\left[\frac{{\hat{p}}^{2}}{2m}+{H}_{0}+{\chi }_{L}(\hat{x}){H}_{US}\right]\vert \psi \rangle =E\vert \psi \rangle \end{equation} \tag{ 29 }$$

for quantum states |ψ⟩ that behave as plane waves outside the scattering region [−L/2, L/2]. These solutions are called scattering states and are not proper quantum states since they are not normalizable. They can be written in terms of the eigenvectors of H₀ and H: $\vert {s}_{J}\rangle ,\vert {s}_{J}^{\prime }\rangle \in {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$, respectively (see equation (2)). In position representation, the scattering states, $\langle x\vert \psi \rangle \in {\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$, read

$$\begin{equation}\langle x\vert \psi \rangle =\begin{cases}\sum\limits _{J}\left({\alpha }_{J}\enspace {\mathrm{e}}^{\mathrm{i}{k}_{J}x}+{\beta }_{J}\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{J}x}\right)\vert {s}_{J}\rangle \quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace x< -L/2\hfill \\ \sum\limits _{J}\left({\alpha }_{J}^{\prime }\enspace {\mathrm{e}}^{\mathrm{i}{k}_{J}^{\prime }x}+{\beta }_{J}^{\prime }\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{J}^{\prime }x}\right)\vert {s}_{J}^{\prime }\rangle \quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace -L/2< x< L/2\hfill \\ \sum\limits _{J}\left({\alpha }_{J}^{{\prime\prime}}\enspace {\mathrm{e}}^{\mathrm{i}{k}_{J}x}+{\beta }_{J}^{{\prime\prime}}\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{J}x}\right)\vert {s}_{J}\rangle \quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace L/2< x.\hfill \end{cases}\end{equation} \tag{ 30 }$$

Inserting this wave function into the Schrödinger equation (29), one obtains the following energy conservation condition for the wave vectors k_J and ${k}_{J}^{\prime }$:

$$\begin{equation}\frac{{k}_{J}^{2}}{2m}+{e}_{J}=\frac{{k}_{J}^{\prime 2}}{2m}+{e}_{J}^{\prime }=E\quad \text{for}\enspace \text{all}\enspace J.\end{equation} \tag{ 31 }$$

This condition fixes the value of the wave vectors k_J and ${k}_{J}^{\prime }$, which can be real or imaginary depending on the energy E. An imaginary wave vector k_J implies an exponential decay outside the scattering region, which does not describe a scattering event. This is why the scattering matrix is defined only for states |s_J ⟩ with real k_J , which span the subspace of open channels ${\mathcal{H}}_{\text{open}}$ for a given energy E, introduced in equation (9). On the other hand, ${k}_{J}^{\prime }$ can be real or imaginary, the latter case corresponding to channels where the transmission is due to quantum tunneling.

The amplitudes of the scattering state (30) in the different segments of the real line, ${\alpha }_{J},{\beta }_{J},{\alpha }_{J}^{\prime },{\beta }_{J}^{\prime },{\alpha }_{J}^{{\prime\prime}}$, and ${\beta }_{J}^{{\prime\prime}}$, are determined by imposing the continuity and differentiability of the wave function at x = −L/2 and x = L/2. It is convenient to write the amplitudes as coordinates of vectors in the internal Hilbert space ${\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$:

$$\begin{equation}\begin{aligned}\hfill \vert a\rangle =\sum\limits _{J}{\alpha }_{J}\vert {s}_{J}\rangle ;\qquad \vert {a}^{\prime }\rangle =\sum\limits _{J}{\alpha }_{J}^{\prime }\vert {s}_{J}^{\prime }\rangle ;\qquad & \vert {a}^{{\prime\prime}}\rangle =\sum\limits _{J}{\alpha }_{J}^{{\prime\prime}}\vert {s}_{J}\rangle \hfill \\ \hfill \vert b\rangle =\sum\limits _{J}{\beta }_{J}\vert {s}_{J}\rangle ;\qquad \vert {b}^{\prime }\rangle =\sum\limits _{J}{\beta }_{J}^{\prime }\vert {s}_{J}^{\prime }\rangle ;\qquad & \vert {b}^{{\prime\prime}}\rangle =\sum\limits _{J}{\beta }_{J}^{{\prime\prime}}\vert {s}_{J}\rangle .\hfill \end{aligned}\end{equation} \tag{ 32 }$$

We also introduce two operators, acting on the internal Hilbert space ${\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}$, which will play an important role in the rest of the paper:

$$\begin{equation}{\mathbb{K}}_{0}(E)\equiv \sqrt{2m(E-{H}_{0})}\mathbb{K}(E)\equiv \sqrt{2m(E-H)}.\end{equation} \tag{ 33 }$$

We call them WVOs, since their eigenvalues are the wave vectors corresponding to a given energy E: ${\mathbb{K}}_{0}(E)\vert {s}_{J}\rangle ={k}_{J}\vert {s}_{J}\rangle$ and $\mathbb{K}(E)\vert {s}_{J}^{\prime }\rangle ={k}_{J}^{\prime }\vert {s}_{J}^{\prime }\rangle$. Notice that they are self-adjoint only for sufficiently high energy. In particular, ${\mathbb{K}}_{0}$ is self-adjoint when restricted to ${\mathcal{H}}_{\text{open}}$.

The boundary conditions allow us to eliminate the intermediate amplitudes |a'⟩ and |b'⟩ and find a relationship between the rest. The relationship can be written as

$$\begin{equation}\left(\begin{matrix}\hfill \vert {a}^{{\prime\prime}}\rangle \hfill \\ \hfill \vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)=\mathcal{M}\left(\begin{matrix}\hfill \vert a\rangle \hfill \\ \hfill \vert b\rangle \hfill \end{matrix}\right).\end{equation} \tag{ 34 }$$

$\mathcal{M}$ is called the transfer matrix and connects the amplitudes of the plane waves at the right and at the left sides of the scatterer (see figure 2). Notice that it is a matrix defined in the Hilbert space $[{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]\oplus [{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]$. In appendix A, we obtain the following closed expression for the transfer matrix from the boundary conditions:

$$\begin{equation}\mathcal{M}={\mathbb{M}}^{-1}(L/2,{\mathbb{K}}_{0})\mathbb{M}(L/2,\mathbb{K}){\mathbb{M}}^{-1}(-L/2,\mathbb{K})\mathbb{M}(-L/2,{\mathbb{K}}_{0}),\end{equation} \tag{ 35 }$$

where we have introduced the matrix $\mathbb{M}(x,\mathbb{K})$ acting on the Hilbert space $[{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]\oplus [{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]$ and depending on a position x and an operator $\mathbb{K}$:

$$\begin{equation}\mathbb{M}(x,\mathbb{K})\equiv \left(\begin{matrix}\hfill {\text{e}}^{\text{i}\mathbb{K}x}\hfill & \hfill {\text{e}}^{-\text{i}\mathbb{K}x}\hfill \\ \hfill \mathbb{K}\enspace {\text{e}}^{\text{i}\mathbb{K}x}\hfill & \hfill -\mathbb{K}\enspace {\text{e}}^{-\text{i}\mathbb{K}x}\hfill \end{matrix}\right).\end{equation} \tag{ 36 }$$

Zoom In Zoom Out Reset image size

An alternative way of relating the amplitudes of the plane waves is the matrix $\mathcal{S}$ connecting the incoming and outgoing amplitudes (see figure 2):

$$\begin{equation}\left(\begin{matrix}\hfill \vert b\rangle \hfill \\ \hfill \vert {a}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)=\mathcal{S}\left(\begin{matrix}\hfill \vert a\rangle \hfill \\ \hfill \vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\mathcal{S}}_{11}\vert a\rangle +{\mathcal{S}}_{12}\vert {b}^{{\prime\prime}}\rangle \hfill \\ \hfill {\mathcal{S}}_{21}\vert a\rangle +{\mathcal{S}}_{22}\vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right).\end{equation} \tag{ 37 }$$

A direct comparison between (34) and (37) yields [21]:

$$\begin{align}\hfill & {\mathcal{S}}_{11}=-{\mathcal{M}}_{22}^{-1}{\mathcal{M}}_{21}\qquad {\mathcal{S}}_{12}={\mathcal{M}}_{22}^{-1}\hfill \\ \hfill & {\mathcal{S}}_{21}={\mathcal{M}}_{11}-{\mathcal{M}}_{12}{\mathcal{M}}_{22}^{-1}{\mathcal{M}}_{21}\qquad {\mathcal{S}}_{22}={\mathcal{M}}_{12}{\mathcal{M}}_{22}^{-1}.\hfill \end{align} \tag{ 38 }$$

The matrix $\mathcal{S}$ is not exactly the scattering matrix $\tilde{\mathcal{S}}(E)$, defined in equation (10), for two reasons. First, $\mathcal{S}$ acts on the whole Hilbert space $[{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]\oplus [{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]$, whereas the scattering matrix $\tilde{\mathcal{S}}$ introduced in equation (10) is restricted to the space open channels, ${\mathcal{H}}_{\text{open}}$, spanned by the eigenstates with real wave vectors k_J . Second, the entries of $\tilde{\mathcal{S}}$ are the transmission and reflection amplitudes, which are given respectively by the amplitudes of the transmitted and reflected waves multiplied by the ratio of outgoing to incoming momenta [18, 23]. According to equation (10), the scattering matrix can be written as the following operator acting on ${\mathcal{H}}_{\text{open}}\oplus {\mathcal{H}}_{\text{open}}$:

$$\begin{equation}\tilde{\mathcal{S}}\equiv \left(\begin{matrix}\hfill {\mathbb{K}}_{0}^{1/2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathbb{K}}_{0}^{1/2}\hfill \end{matrix}\right)\mathcal{S}\enspace \left(\begin{matrix}\hfill {\mathbb{K}}_{0}^{-1/2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathbb{K}}_{0}^{-1/2}\hfill \end{matrix}\right)={\mathbb{K}}_{0}^{1/2}\mathcal{S}{\mathbb{K}}_{0}^{-1/2}\end{equation} \tag{ 39 }$$

where all the operators ${\mathbb{K}}_{0}$ and ${\mathcal{S}}_{ij}$ are restricted to ${\mathcal{H}}_{\text{open}}$.

The matrix $\mathcal{M}$ can be calculated exactly from equation (35). However, here we introduce an approximation that preserves the symmetries and the unitarity of the scattering matrix and therefore provides a simple implementation of a thermal reservoir. The approximation is valid for incident particles with large kinetic energy. More precisely, if $E\gg {e}_{J},{e}_{J}^{\prime }$, all the wave vectors are approximately equal, ${k}_{J}\simeq {k}_{J}^{\prime }\simeq \sqrt{2mE}$, and we can approximate ${\mathbb{K}}^{-1}{\mathbb{K}}_{0}\simeq \mathbb{I}$ yielding (see appendix C for a detailed calculation)

$$\begin{equation}\mathcal{M}\simeq \left(\begin{matrix}\hfill {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\enspace {\text{e}}^{\text{i}\mathbb{K}L}\enspace {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\text{e}}^{\text{i}{\mathbb{K}}_{0}L/2}\enspace {\text{e}}^{-\text{i}\mathbb{K}L}\enspace {\text{e}}^{\text{i}{\mathbb{K}}_{0}L/2}\hfill \end{matrix}\right).\end{equation} \tag{ 40 }$$

Using equation (38) and the definition of the scattering matrix (39), we find that the scattering matrix $\tilde{\mathcal{S}}$ and the matrix $\mathcal{S}$ in this approximation are

$$\begin{equation}\tilde{\mathcal{S}}\simeq \mathcal{S}\simeq \left(\begin{matrix}\hfill 0\hfill & \hfill {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\enspace {\text{e}}^{\text{i}\mathbb{K}L}\enspace {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\hfill \\ \hfill {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\enspace {\text{e}}^{\text{i}\mathbb{K}L}\enspace {\text{e}}^{-\text{i}{\mathbb{K}}_{0}L/2}\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 41 }$$

We see that, in this approximation, the reflection amplitudes vanish and the transition amplitudes read

$$\begin{equation}\langle {s}_{{J}^{\prime }}\vert \mathbf{t}\vert {s}_{J}\rangle \simeq {\text{e}}^{-\text{i}({k}_{J}+{k}_{{J}^{\prime }})L/2}\langle {s}_{{J}^{\prime }}\vert {\text{e}}^{\text{i}\mathbb{K}L}\vert {s}_{J}\rangle .\end{equation} \tag{ 42 }$$

This matrix $\mathcal{S}$ is unitary and symmetric for a given energy $E={k}_{J}^{2}/(2m)+{e}_{J}={k}_{{J}^{\prime }}^{2}/(2m)+{e}_{{J}^{\prime }}$ if $\mathbb{K}$ is self-adjoint, that is, if all ${k}_{J}^{\prime }$ are real. This occurs if the total energy E is larger than the maximum eigenvalue of H. Hence, for sufficiently high incident kinetic energy, this approximation fulfills all the symmetries of the original collision problem. To see that the matrix is symmetric, notice that $[H,{\mathsf{T}}_{\text{int}}]=0$ implies ${\text{e}}^{\text{i}\mathbb{K}L}{\mathsf{T}}_{\text{int}}={\mathsf{T}}_{\text{int}}\enspace {\text{e}}^{-\text{i}\mathbb{K}L}$ in the subspace where $\mathbb{K}$ is self-adjoint. Therefore, for energies E larger than the maximum eigenvalue of H, we have

$$\begin{align}\hfill \langle {s}_{{J}^{\prime }}\vert {\text{e}}^{\text{i}\mathbb{K}L}\vert {s}_{J}\rangle & =({\mathsf{T}}_{\text{int}}\vert {s}_{{J}^{\prime }}\rangle ,{\text{e}}^{\text{i}\mathbb{K}L}{\mathsf{T}}_{\text{int}}\vert {s}_{J}\rangle )=({\mathsf{T}}_{\text{int}}\vert {s}_{{J}^{\prime }}\rangle ,{\mathsf{T}}_{\text{int}}\enspace {\text{e}}^{-\text{i}\mathbb{K}L}\vert {s}_{J}\rangle )\hfill \\ \hfill & ={(\vert {s}_{{J}^{\prime }}\rangle ,{\text{e}}^{-\text{i}\mathbb{K}L}\vert {s}_{J}\rangle )}^{\ast }={\langle {s}_{{J}^{\prime }}\vert {\text{e}}^{-\text{i}\mathbb{K}L}\vert {s}_{J}\rangle }^{\ast }\hfill \\ \hfill & =\langle {s}_{J}\vert {\text{e}}^{\text{i}\mathbb{K}L}\vert {s}_{{J}^{\prime }}\rangle .\hfill \end{align} \tag{ 43 }$$

Here (⋅, ⋅) is the scalar product in the Hilbert space and we have used the time-reversal invariance of the eigenstates of H₀, ${\mathsf{T}}_{\text{int}}\vert {s}_{J}\rangle =\vert {s}_{J}\rangle$ for all J, and that any anti-unitary operator verifies $(\mathsf{T}\vert a\rangle ,\mathsf{T}\vert b\rangle )={(\vert a\rangle ,\vert b\rangle )}^{\ast }$. Notice that the same symmetry holds if we replace $\mathbb{K}L$ by any real function of H.

Equation (42) is a valid approximation for high incident kinetic energy. To complete our first thermostat model, we need an expression for the transmission and reflection amplitudes at low velocities. In order to preserve micro-reversibility and the unitarity of the scattering matrix, we adopt the simplest assumption for low energies, namely, that the incident unit is reflected without affecting the state of the system. This is also justified by the fact that in a large scatterer the transmission amplitudes corresponding to tunneling vanish. However, from the point of view of the system, it does not matter whether the unit is reflected or transmitted, as long as it does not affect the system. Then, for simplicity, we define our model of a collisional thermostat as given by vanishing reflection amplitudes, r_J'J (E) = 0 for all E, and the following transmission amplitudes:

$$\begin{equation}{t}_{{J}^{\prime }J}(E)=\begin{cases}{\mathrm{e}}^{-\mathrm{i}L({k}_{J}+{k}_{{J}^{\prime }})/2}\langle {s}_{{J}^{\prime }}\vert {\mathrm{e}}^{\mathrm{i}L\mathbb{K}\left(E\right)}\vert {s}_{J}\rangle \quad \hfill & \mathrm{i}\mathrm{f}\enspace E > {e}_{\mathrm{m}\mathrm{a}\mathrm{x}}\hfill \\ {\delta }_{{J}^{\prime }J}\quad \hfill & \mathrm{i}\mathrm{f}\enspace E\leqslant {e}_{\mathrm{m}\mathrm{a}\mathrm{x}}\hfill \end{cases}\end{equation} \tag{ 44 }$$

where $E={p}_{0}^{2}/(2m)+{e}_{J}$ and e_max is the maximum of the eigenvalues of H and H₀. With this choice

$$\begin{equation}\sum\limits _{{J}^{\prime }}\left[\vert {t}_{{J}^{\prime }J}(E){\vert }^{2}+\vert {r}_{{J}^{\prime }J}(E){\vert }^{2}\right]=1\end{equation} \tag{ 45 }$$

for all J and E, ensuring the conservation of the trace of the density matrix Tr(ρ') = Tr(ρ).

The transmission amplitudes defined by equation (44) obey condition (21), as shown in the previous section 3.3. We conclude that our model, based on the WVO $\mathbb{K}$, induces the thermalization of the system. Consequently, it constitutes a simple model of a thermostat. Furthermore, it is also a good approximation of a system colliding with units that escape from a thermal reservoir, specially for large scatterers. In section 5, we check the validity of this approximation in explicit examples.

We now present a second model that also induces thermalization and is more directly related to the repeated interaction schemes considered in the literature [11, 13], where the interaction H_US is switched on for a time interval. This can be done if the incident momentum is large and we can further expand the operator $\mathbb{K}(E)=\sqrt{2m(E-H)}$ as

$$\begin{equation}\mathbb{K}(E)\simeq \sqrt{2mE}-\sqrt{\frac{m}{2E}}\enspace H.\end{equation} \tag{ 46 }$$

An analogous expansion of the wave vectors outside the scattering region yields ${k}_{J}\simeq \sqrt{2mE}-{e}_{J}\sqrt{m/(2E)}$. Inserting these expressions in the transmission amplitudes given by equation (42), we get

$$\begin{equation}{t}_{{J}^{\prime }J}(E)\simeq {\text{e}}^{\text{i}\tau (E)[{e}_{{J}^{\prime }}+{e}_{J}]/2}\enspace \langle {s}_{{J}^{\prime }}\vert {\text{e}}^{-\text{i}\tau (E)H}\vert {s}_{J}\rangle .\end{equation} \tag{ 47 }$$

Here, we have introduced the time

$$\begin{equation}\tau (E)\equiv \frac{L}{{v}_{E}}=\frac{L}{\sqrt{2E/m}},\end{equation} \tag{ 48 }$$

which is the time a classical particle with velocity ${v}_{E}\equiv \sqrt{2E/m}=\sqrt{{({p}_{0}/m)}^{2}+2{e}_{J}/m}$ takes to cross the scattering region of length L. Except for a phase, equation (47) is equivalent to the evolution of the state |s_J ⟩ under the total internal Hamiltonian H = H₀ + H_US during a time τ(E), which is, approximately, the interaction time between the wave packet and the scatterer. We thus recover for the transmission amplitudes the usual picture that ignores the translational part of the unit and considers that the coupling is switched on for a given time τ(E) [11, 13]. Notice however that this velocity does not exactly coincide with the velocity of the wave packet v_packet ≡ p₀/m. In fact, the scattering matrix resulting from setting the interaction time equal to τ_packet ≡ L/v_packet = Lm/p₀ does not obey micro-reversibility and does not thermalize the system, as we show in section 5 for an specific example.

Two important remarks must be made. First, expression (47) can be evaluated for any pair of states, J and J', and any positive energy E. However, the transfer matrix is only defined for open channels, obeying E ⩾ e_J , e_J'. To preserve the unitarity of the scattering matrix we have to restrict the use of (47) to energies E larger than the eigenvalues of H₀. To be consistent with the wave-vector-operator (WVO) model, we adopt here a more conservative strategy, restricting expression (47) to energies higher than e_max, the maximum eigenvalue of both H₀ and H. We also assume that e_max is positive, to avoid a negative total energy E, which would yield an imaginary velocity v_E . Summarizing, our model consists of null reflection amplitudes, r_J'J (E) = 0 for all E, and transmission amplitudes given by

$$\begin{equation}{t}_{{J}^{\prime }J}(E)=\begin{cases}{\text{e}}^{\text{i}\tau (E)[{e}_{{J}^{\prime }}+{e}_{J}]/2}\enspace \langle {s}_{{J}^{\prime }}\vert {\text{e}}^{-\text{i}\tau (E)H}\vert {s}_{J}\rangle \quad \hfill & \mathrm{i}\mathrm{f}\enspace E > {e}_{\mathrm{m}\mathrm{a}\mathrm{x}}\hfill \\ {\delta }_{{J}^{\prime }J}\quad \hfill & \mathrm{i}\mathrm{f}\enspace E\leqslant {e}_{\mathrm{m}\mathrm{a}\mathrm{x}}.\hfill \end{cases}\end{equation} \tag{ 49 }$$

With this definition, our random-interaction-time (RIT) model, like the WVO model, preserves all the properties of the exact scattering matrix and consequently induces the thermalization of the system. As in the previous model, we have set to zero the reflection amplitudes for E ⩽ e_max, although for low energies the unit is most likely reflected. The choice in (49) is equivalent to assume that the system is not affected if E ⩽ e_max, independently of whether the unit is reflected or transmitted.

The second remark is to notice that the interaction time τ(E) depends on the zero of the total energy, that is, if we add a constant E₀ to the total Hamiltonian H_tot, τ(E) changes. Then the transition amplitudes will depend as well on the zero of energy. The reason of this dependency is the Taylor expansion around H = 0 in equation (46). Shifting the internal energies an amount E₀ is equivalent to expanding the square root around H = −E₀ in equation (46). Hence, to minimize the error in the expansion we have to choose the zero of energy in such a way that H is small. There are several criteria to define the 'smallness' of an operator, based on different matrix norms. For the example in section 5, we minimize the spectral norm of H, which is the square root of the largest eigenvalue of H^† H = H². Notice however that all the models obtained by an energy shift with e_max > 0 are effective thermostats when the scatterer is bombarded by equilibrium units, since equation (49) is unitary and fulfills micro-reversibility.

The model given by equation (49) and the corresponding scattering map in equation (7) are similar to the ones previously considered in the literature [11, 13], except for the randomization of the interaction time τ(E), which depends on the initial state |s_J ⟩, and for the removal of coherences due to tracing out the outgoing narrow packets [18].

To further explore the differences and similarities between our thermostats and a repeated-interaction reservoir, it is convenient to use the Kraus representation of the scattering map given by equation (5), together with (7) and condition (8). If we neglect the reflecting amplitudes, a representation of this map is given by the following Kraus operators:

$$\begin{equation}{M}_{l}=\sum\limits _{J,{J}^{\prime }}{t}_{{J}^{\prime }J}({E}_{{p}_{0}}+{e}_{J})\enspace {\delta }_{{{\Delta}}_{{J}^{\prime }J},{{\Delta}}_{l}}\enspace \vert {s}_{{J}^{\prime }}\rangle \langle {s}_{J}\vert \end{equation} \tag{ 50 }$$

where δ is a Kronecker delta and Δ_l runs over all possible Bohr frequencies of the free internal Hamiltonian H₀. Indeed, the map

$$\begin{equation}{\rho }^{\prime }=\sum\limits _{l}{M}_{l}\rho {M}_{l}^{{\dagger}}\end{equation} \tag{ 51 }$$

in the eigenbasis of H₀ is given by the tensor

$$\begin{equation}{\mathbb{S}}_{{J}^{\prime }{K}^{\prime }}^{JK}=\sum\limits _{l}\langle {s}_{{J}^{\prime }}\vert {M}_{l}\vert {s}_{J}\rangle \langle {s}_{K}\vert {M}_{l}^{{\dagger}}\vert {s}_{{K}^{\prime }}\rangle ,\end{equation} \tag{ 52 }$$

which coincides with the one given by equations (5), (7) and (8), if the reflection amplitudes are neglected. If we now use the approximation (49), the Kraus operators in the eigenbasis of H₀ read

$$\begin{equation}\langle {s}_{{J}^{\prime }}\vert {M}_{l}\vert {s}_{J}\rangle ={\text{e}}^{-\text{i}\tau (E)({e}_{J}+{e}_{{J}^{\prime }})/2}\langle {s}_{{J}^{\prime }}\vert {\text{e}}^{-\text{i}\tau (E)H}\vert {s}_{J}\rangle {\delta }_{{{\Delta}}_{{J}^{\prime }J},{{\Delta}}_{l}}\end{equation} \tag{ 53 }$$

with $E={p}_{0}^{2}/(2m)+{e}_{J}$.

On the other hand, the Kraus representation of the unitary evolution in a repeated-interaction scheme consists of a unique unitary operator M given by

$$\begin{equation}\langle {s}_{{J}^{\prime }}\vert M\vert {s}_{J}\rangle =\langle {s}_{{J}^{\prime }}\vert {\text{e}}^{-\text{i}{\tau }_{\text{int}}H}\vert {s}_{J}\rangle \end{equation} \tag{ 54 }$$

where τ_int is the interaction time. Comparing (53) and (54), we see three main differences: first, the Kronecker delta kills all coherences between jumps with different Bohr frequencies. Recall that the map acting on a pure state can be seen as the application of a randomly chosen operator M_l [24]. The Kronecker delta only allows for superpositions with the same energy jump Δ_l . This is a consequence of using narrow packets, which is a necessary condition for thermalization, as proved in reference [18]. Remarkably, this condition has also been shown to be necessary to derive a fluctuation theorem for quantum maps (see equation (12) in reference [24]), and is equivalent to imposing that the energy exchange with the reservoir, i.e. the heat, is well defined for each possible transformation of a pure state given by the Kraus operators. This implies that heat is well defined for any quantum stochastic trajectory [24]. Second, the interaction time in the RIT model, τ(E), depends on the energy of the initial state e_J . Third, there is an extra phase that appears in the solution of the scattering problem, although it does not play a role in thermalization.

We first check whether the two models presented in the previous section are able to reproduce the transition probabilities P_J'J (p₀). We calculate the transfer matrix given by equation (35) and compare the exact transition probabilities in equation (13) for a given incident momentum p₀ with the ones obtained from the WVO model, equation (44), and the random interaction model (RIT), equation (49). The comparison is shown in figure 3 as a function of the kinetic energy p₀/(2m), for the transition |00⟩_US → |11⟩_US . As expected, the two models reproduce with good accuracy the exact transition probabilities for high kinetic energy. It is remarkable that the wave vector operator model is a very good approximation of the scattering problem even for low kinetic energy, as shown in the inset, whereas the RIT model fails in this regime.

Zoom In Zoom Out Reset image size

We now bombard the qubit with narrow wave packets with random momentum, according to the effusion distribution at temperature T, and a random internal state, according to the Boltzmann distribution at the same temperature. We simulate 500 quantum trajectories where the system jumps between pure eigenstates of the Hamiltonian H_S , and calculate the steady population of the two levels of the qubit. Each trajectory is 10 000 collision long, providing sufficient statistics to neglect the uncertainty.

As expected, the system thermalizes not only for the exact solution of the scattering problem, given by equation (35), but also for the two models, equation (44) and (49). In figure 4, we plot the stationary population of the ground state in the three cases and in the thermal state.

Zoom In Zoom Out Reset image size

To stress the importance of micro-reversibility for thermalization, we also plot in the figure with dark blue circles the population when the interaction time is chosen as ${\tau }_{\text{packet}}\equiv \tau ({p}_{0}^{2}/(2m))=L/{v}_{\text{packet}}$, where v_packet = p₀/m is the velocity of the incoming wave packet. In this case, the system does not reach the temperature of the reservoir and can even exhibit population inversion (see appendix D for a detailed discussion of the model and reference [25] for a general discussion of the phenomenon within the repeated interaction framework). From the point of view of the dynamics of the system, it is striking that the replacement of ${v}_{E}=\sqrt{2E/m}$ by v_packet = p₀/m in the calculation of the interaction time has such significant consequences. Notice however that the interaction time in the RIT model with time given by equation (48) depends both on the energy of the system ${e}_{{j}_{S}}^{(S)}$ and of the internal state of the unit ${e}_{{j}_{U}}^{(U)}$.

Finally, we check the two models in a non-equilibrium scenario where the internal states of the unit are in equilibrium at temperature T_int, different from the temperature T_kin of the effusion distribution, i.e. the state of the unit is given by (14) and (15) but with temperatures T_int and T_kin, respectively. In this situation, the system is exchanging heat with two different thermal baths and reaches a non-equilibrium steady state where a heat Q is transferred from the hot to the cold bath in each collision. In our case, the two baths are the internal and the kinetic degrees of freedom of the unit. The irreversible heat transfer induces an entropy production per collision ΔS = Q|1/T_int − 1/T_kin|, which is shown in figure 5 as a function of the kinetic temperature T_kin, for a fixed internal temperature T_int = 20 and for the two different models and the exact solution of the scattering matrix. In this simulation, the heat Q is calculated as minus the change of the internal energy of the unit in each collision, which is then averaged over quantum trajectories in the steady state. We see in the figure that, for this range of temperatures, the two thermostats are accurate approximations of the exact solution of the scattering problem, even far from equilibrium.

Zoom In Zoom Out Reset image size

If we write the scattering states in the following from

$$\begin{equation}\langle x\vert \psi \rangle =\sum\limits _{J}{\psi }_{J}(x)\vert {s}_{J}\rangle \end{equation} \tag{ B1 }$$

and introduce this expression in the Schödinger equation (29), we get

$$\begin{equation}-\frac{1}{2m}\frac{{\partial }^{2}{\psi }_{J}(x)}{\partial {x}^{2}}+({e}_{J}-E){\psi }_{J}(x)+{\chi }_{L}(x)\sum\limits _{K}{\psi }_{K}(x)\langle {s}_{J}\vert {H}_{US}\vert {s}_{K}\rangle =0.\end{equation} \tag{ B2 }$$

The following generalization of the Wronskian

$$\begin{equation}W(x)\equiv \sum\limits _{J}\left[\frac{\partial {\psi }_{J}(x)}{\partial x}{\psi }_{J}^{\ast }(x)-{\psi }_{J}(x)\frac{\partial {\psi }_{J}^{\ast }(x)}{\partial x}\right]\end{equation} \tag{ B3 }$$

can be interpreted as a total current of particles and is independent of x. To prove it, we use the Schrödinger equation to compute the derivative

$$\begin{align}\hfill \frac{\mathrm{d}W(x)}{\mathrm{d}x}& =\sum\limits _{J}\left[\frac{{\partial }^{2}{\psi }_{J}(x)}{\partial {x}^{2}}{\psi }_{J}^{\ast }(x)-{\psi }_{J}(x)\frac{{\partial }^{2}{\psi }_{J}^{\ast }(x)}{\partial {x}^{2}}\right]\hfill \\ \hfill & =2m{\chi }_{L}(x)\sum\limits _{J,K}\left[{\psi }_{K}(x){\psi }_{J}^{\ast }(x)\langle {s}_{J}\vert {H}_{US}\vert {s}_{K}\rangle -{\psi }_{J}(x){\psi }_{K}^{\ast }(x){\langle {s}_{J}\vert {H}_{US}\vert {s}_{K}\rangle }^{\ast }\right]=0.\hfill \end{align} \tag{ B4 }$$

To obtain the last equality, we have taken into account that H_US is self-adjoint and, consequently, the two terms in the sum are equal under a permutation of the indexes.

For a given total energy E, the solution (30) at x → −∞ corresponds to ${\psi }_{J}(x)={\alpha }_{J}\enspace {\mathrm{e}}^{\mathrm{i}{k}_{J}x}+{\beta }_{J}\enspace {\mathrm{e}}^{-\mathrm{i}{k}_{J}x}$ if k_J is real (e_J ⩽ E) and ψ_J (x) ≃ 0 if k_J is imaginary (e_J ⩾ E). Then, the Wronskian reads

$$\begin{equation}W(x)\simeq \sum\limits _{J:{e}_{J}\leqslant E}2\mathrm{i}{k}_{J}\left[\vert {\alpha }_{J}{\vert }^{2}-\vert {\beta }_{J}{\vert }^{2}\right].\end{equation} \tag{ B5 }$$

Similarly, for x → ∞:

$$\begin{equation}W(x)\simeq \sum\limits _{J:{e}_{J}\leqslant E}2\mathrm{i}{k}_{J}\left[\vert {\alpha }_{J}^{{\prime\prime}}{\vert }^{2}-\vert {\beta }_{J}^{{\prime\prime}}{\vert }^{2}\right].\end{equation} \tag{ B6 }$$

Therefore

$$\begin{equation}\sum\limits _{J:{e}_{J}\leqslant E}{k}_{J}\left[\vert {\alpha }_{J}{\vert }^{2}-\vert {\beta }_{J}{\vert }^{2}\right]=\sum\limits _{J:{e}_{J}\leqslant E}{k}_{J}\left[\vert {\alpha }_{J}^{{\prime\prime}}{\vert }^{2}-\vert {\beta }_{J}^{{\prime\prime}}{\vert }^{2}\right].\end{equation} \tag{ B7 }$$

Notice that the sums in the previous expressions run only over the states |s_J ⟩ with k_J real, i.e. the incoming and outgoing plane waves.

The conservation of the total probability current (B7) imposes some constraints on the matrices $\mathcal{M}$ and $\mathcal{S}$. Let ${\mathbb{P}}_{\text{open}}$ be the projector onto ${\mathcal{H}}_{\text{open}}$, i.e. onto the eigenstates |s_J ⟩ with k_J real:

$$\begin{equation}{\mathbb{P}}_{\text{open}}=\sum\limits _{J:{e}_{J}\leqslant E}\vert {s}_{J}\rangle \langle {s}_{J}\vert \end{equation} \tag{ B8 }$$

and let us define the operator acting on $[{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]\oplus [{\mathcal{H}}_{U,\mathrm{i}\mathrm{n}\mathrm{t}}\otimes {\mathcal{H}}_{S}]$

$$\begin{equation}\mathcal{P}=\left(\begin{matrix}\hfill {\mathbb{P}}_{\text{open}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathbb{P}}_{\text{open}}\hfill \end{matrix}\right),\end{equation} \tag{ B9 }$$

which verifies ${\mathcal{P}}^{2}=\mathcal{P}$. Condition (B7) can be written as

$$\begin{equation}\left(\langle a\vert \enspace \langle b\vert \right)\left(\begin{matrix}\hfill {\mathbb{K}}_{0}{\mathbb{P}}_{\text{open}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}{\mathbb{P}}_{\text{open}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \vert a\rangle \hfill \\ \hfill \vert b\rangle \hfill \end{matrix}\right)=\left(\langle {a}^{{\prime\prime}}\vert \enspace \langle {b}^{{\prime\prime}}\vert \right)\left(\begin{matrix}\hfill {\mathbb{K}}_{0}{\mathbb{P}}_{\text{open}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}{\mathbb{P}}_{\text{open}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \vert {a}^{{\prime\prime}}\rangle \hfill \\ \hfill \vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right).\end{equation} \tag{ B10 }$$

Applying the relationship (34) between the amplitudes of the waves at the right and left sides of the scatterer via the transfer matrix, we obtain

$$\begin{equation}\left(\begin{matrix}\hfill {\mathbb{K}}_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}\hfill \end{matrix}\right)\mathcal{P}={\mathcal{M}}^{{\dagger}}\left(\begin{matrix}\hfill {\mathbb{K}}_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}\hfill \end{matrix}\right)\mathcal{P}\mathcal{M}\end{equation} \tag{ B11 }$$

and, multiplying by $\mathcal{P}$ from right, we get

$$\begin{equation}{\mathcal{M}}^{{\dagger}}\left(\begin{matrix}\hfill {\mathbb{K}}_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}\hfill \end{matrix}\right)\mathcal{P}\mathcal{M}={\mathcal{M}}^{{\dagger}}\left(\begin{matrix}\hfill {\mathbb{K}}_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\mathbb{K}}_{0}\hfill \end{matrix}\right)\mathcal{P}\mathcal{M}\mathcal{P}.\end{equation} \tag{ B12 }$$

Since ${\mathcal{M}}^{{\dagger}}$ and ${\mathbb{K}}_{0}$ are both invertible in their respective Hilbert spaces (we assume that k_J ≠ 0 for all J), we conclude that $\mathcal{P}\mathcal{M}\mathcal{P}=\mathcal{P}\mathcal{M}$ or $\mathcal{P}\mathcal{M}(\mathbb{I}-\mathcal{P})=0$. This relationship indicates that the amplitudes of the real exponentials (k_J imaginary) do not affect the amplitudes of the plane waves (k_J real) and that we can restrict ourselves to ${\mathcal{H}}_{\text{open}}$. Notice however that $\mathcal{P}$ and $\mathcal{M}$ do not necessarily commute, i.e. ${\mathcal{H}}_{\text{open}}$ is not in general invariant under the transfer matrix $\mathcal{M}$. However, the action of $\mathcal{M}$ on vectors in ${\mathcal{H}}_{\text{open}}$ is entirely determined by its restriction to this subspace $\mathcal{P}\mathcal{M}\mathcal{P}$. In particular any power n of $\mathcal{M}$ verifies $\mathcal{P}{\mathcal{M}}^{n}\mathcal{P}$ = ${(\mathcal{P}\mathcal{M}\mathcal{P})}^{n}$ and the inverse of $\mathcal{M}$ in ${\mathcal{H}}_{\text{open}}$ is $\mathcal{P}{\mathcal{M}}^{-1}\mathcal{P}$, that is $[\mathcal{P}\mathcal{M}\mathcal{P}][\mathcal{P}{\mathcal{M}}^{-1}\mathcal{P}]=[\mathcal{P}{\mathcal{M}}^{-1}\mathcal{P}][\mathcal{P}\mathcal{M}\mathcal{P}]=\mathcal{P}$. The same arguments apply to the matrix $\mathcal{S}$, which obeys $\mathcal{P}\mathcal{S}(\mathbb{I}-\mathcal{P})=0$. Hence, from now on, we can neglect the eigenstates |s_J ⟩ with imaginary k_J and explore the properties of the matrices $\mathcal{M}$ and $\mathcal{S}$ restricted to ${\mathcal{H}}_{\text{open}}$. Nevertheless, we will keep the same notation, for simplicity. Notice also that ${k}_{J}^{\prime }$, the wave vectors within the scattering region [−L/2, L/2], can be imaginary, indicating that the corresponding channel is associated with tunneling.

A second important consequence of the conservation of probability current is the unitarity of the scattering matrix. Condition (B7) can also be written as

$$\begin{equation}\left(\langle a\vert \enspace \langle {b}^{{\prime\prime}}\vert \right){\mathbb{K}}_{0}\left(\begin{matrix}\hfill \vert a\rangle \hfill \\ \hfill \vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)=\left(\langle b\vert \enspace \langle {a}^{{\prime\prime}}\vert \right){\mathbb{K}}_{0}\left(\begin{matrix}\hfill \vert b\rangle \hfill \\ \hfill \vert {a}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)=\left(\langle a\vert \enspace \langle {b}^{{\prime\prime}}\vert \right){\mathcal{S}}^{{\dagger}}{\mathbb{K}}_{0}\mathcal{S}\left(\begin{matrix}\hfill \vert a\rangle \hfill \\ \hfill \vert {b}^{{\prime\prime}}\rangle \hfill \end{matrix}\right)\end{equation} \tag{ B13 }$$

where all vectors and operators are restricted to ${\mathcal{H}}_{\text{open}}$. Hence, ${\mathbb{K}}_{0}={\mathcal{S}}^{{\dagger}}{\mathbb{K}}_{0}\mathcal{S}$ in this subspace and we finally obtain

$$\begin{equation}{\tilde{\mathcal{S}}}^{{\dagger}}\tilde{\mathcal{S}}={\mathbb{K}}_{0}^{-1/2}{\mathcal{S}}^{{\dagger}}{\mathbb{K}}_{0}^{1/2}{\mathbb{K}}_{0}^{1/2}\mathcal{S}\enspace {\mathbb{K}}_{0}^{-1/2}=\mathbb{I}\end{equation} \tag{ B14 }$$

that is, the scattering matrix $\tilde{\mathcal{S}}$ is unitary.

The collision problem that we consider in this paper is invariant under spatial inversion (x, p) → (−x, −p), which is equivalent to the following transformation of vectors (see figure 2):

$$\begin{equation}\vert a\rangle {\leftrightarrow}\vert {b}^{{\prime\prime}}\rangle \qquad \vert b\rangle {\leftrightarrow}\vert {a}^{{\prime\prime}}\rangle .\end{equation} \tag{ B15 }$$

This transformation converts equation (37) into

$$\begin{equation}\left(\begin{matrix}\hfill \vert {a}^{{\prime\prime}}\rangle \hfill \\ \hfill \vert b\rangle \hfill \end{matrix}\right)=\mathcal{S}\left(\begin{matrix}\hfill \vert {b}^{{\prime\prime}}\rangle \hfill \\ \hfill \vert a\rangle \hfill \end{matrix}\right)\end{equation} \tag{ B16 }$$

and comparing this expression with equation (37), we get

$$\begin{equation}\left(\begin{matrix}\hfill 0\hfill & \hfill \mathbb{I}\hfill \\ \hfill \mathbb{I}\hfill & \hfill 0\hfill \end{matrix}\right)\mathcal{S}\left(\begin{matrix}\hfill 0\hfill & \hfill \mathbb{I}\hfill \\ \hfill \mathbb{I}\hfill & \hfill 0\hfill \end{matrix}\right)=\mathcal{S}\end{equation} \tag{ B17 }$$

which yields ${\mathcal{S}}_{11}={\mathcal{S}}_{22}$ and ${\mathcal{S}}_{12}={\mathcal{S}}_{21}$. The same symmetry applies to the scattering matrix $\tilde{\mathcal{S}}$. This symmetry allows us to write the scattering matrix as in equation (10):

$$\begin{equation}\tilde{\mathcal{S}}=\left(\begin{matrix}\hfill \mathbf{r}\hfill & \hfill \mathbf{t}\hfill \\ \hfill \mathbf{t}\hfill & \hfill \mathbf{r}\hfill \end{matrix}\right)\end{equation} \tag{ B18 }$$

where r and t are matrices whose elements are the reflection and transmission amplitudes, respectively. The unitarity of $\tilde{\mathcal{S}}$ derived in equation (B14), can be written now as

$$\begin{equation}\mathbf{r}{\mathbf{r}}^{{\dagger}}+\mathbf{t}{\mathbf{t}}^{{\dagger}}=\mathbb{I}\mathbf{r}{\mathbf{t}}^{{\dagger}}+\mathbf{t}{\mathbf{r}}^{{\dagger}}=0\end{equation} \tag{ B19 }$$

which is equation (B14) in the main text.

The symmetry under time reversal implies that there is an anti-unitary operator ${\mathsf{T}}_{\text{int}}$ in the Hilbert space of internal states that commutes with H₀ and H_US . As discussed in the main text, the total time-reversal operator is $\mathsf{T}={\mathsf{T}}_{U,\mathrm{p}}\otimes {\mathsf{T}}_{\text{int}}$ where ${\mathsf{T}}_{U,\mathrm{p}}$ is the conjugation of the wave function in the position representation. Hence, if |ψ⟩ is given by (B1), then

$$\begin{equation}\mathsf{T}\vert \psi \rangle =\sum\limits _{J}{\psi }_{J}^{\ast }(x){\mathsf{T}}_{\text{int}}\vert {s}_{J}\rangle .\end{equation} \tag{ B20 }$$

For simplicity, we assume that the eigenstates of H₀ are invariant under time reversal, i.e. ${\mathsf{T}}_{\text{int}}\vert {s}_{J}\rangle =\vert {s}_{J}\rangle$. In this case, $[{H}_{US},{\mathsf{T}}_{\text{int}}]=0$ implies that ⟨s_J |H_US |s_K ⟩ is real, that is, the matrix of the interaction Hamiltonian H_US in the eigenbasis of H₀ is real and symmetric. To prove this property, take into account that an anti-unitary operator verifies $(\mathsf{T}\vert a\rangle ,\mathsf{T}\vert b\rangle )={(\vert a\rangle ,\vert b\rangle )}^{\ast }$, where (⋅, ⋅) is the scalar product in the Hilbert space. Hence, we can take the complex conjugate of the Schrödinger equation (B2) and obtain the following transformation under time reversal for the vectors restricted to ${\mathcal{H}}_{\text{open}}$:

$$\begin{equation}\vert a\rangle {\leftrightarrow}\vert {b}^{\ast }\rangle \qquad \vert {a}^{\prime \prime }\rangle {\leftrightarrow}\vert {b}^{\prime \prime \ast }\rangle \end{equation} \tag{ B21 }$$

where $\vert {b}^{\ast }\rangle =\mathsf{T}\vert b\rangle ={\sum }_{J}{\beta }_{J}^{\ast }\vert {s}_{J}\rangle$. This symmetry implies

$$\begin{equation}\left(\begin{matrix}\hfill \vert {a}^{\ast }\rangle \hfill \\ \hfill \vert {b}^{\prime \prime \ast }\rangle \hfill \end{matrix}\right)=\mathcal{S}\left(\begin{matrix}\hfill \vert {b}^{\ast }\rangle \hfill \\ \hfill \vert {a}^{\prime \prime \ast }\rangle \hfill \end{matrix}\right).\end{equation} \tag{ B22 }$$

Using (37), we obtain ${\mathcal{S}}^{\ast }\mathcal{S}={\mathbb{K}}_{0}^{-1/2}{\tilde{\mathcal{S}}}^{\ast }\tilde{\mathcal{S}}\enspace {\mathbb{K}}_{0}^{1/2}=\mathbb{I}$, implying ${\tilde{\mathcal{S}}}^{\ast }\tilde{\mathcal{S}}=\mathbb{I}$, and

$$\begin{align}\hfill \mathbf{r}{\mathbf{r}}^{\ast }+\mathbf{t}{\mathbf{t}}^{\ast }=& \mathbb{I}\hfill \\ \hfill \mathbf{r}{\mathbf{t}}^{\ast }+\mathbf{t}{\mathbf{r}}^{\ast }=& 0.\hfill \end{align} \tag{ B23 }$$

Combining this symmetry with the unitarity of the scattering matrix, equation (B19), we conclude t^† = t* and r^† = r*, i.e. the matrices t and r are symmetric.