# Entropy and Irreversibility in the Quantum Realm

## Model
The author accompanies the symbolic discussions of the topic with vivid simulations depicting the physical effects of interest. The simulations all concern a system of $8$ electrons in which only their spins are taken into account, and the simulations are performed according to the scheme below. I precede the illustration of the scheme with the rationale for adopting such a setting.

### Rationale
1. An electron only concerned of its spin is a two-level system which is the simplest to describe.
2. A composite system of $8$ such two-level subsystems has $2^8 = 256$ eigenstates which serves narrative purposes well while keeping computations manageable.
3. Though any composite system of two-level subsystems would suffice, the author favors a system of electrons for concreteness.

### Scheme
Let $n = 8$, $N = 2^n = 256$, and denote the up-spin and down-spin of the $j^{\text{ th}}$ electron by $\left|1\right>_j$ and $\left|0\right>_j$ respectively where $j = 0,\dots,n-1$. Let $\{\left|k\right>\}_{k=0}^{N-1}$ be the ordered basis for the composite system such that $\displaystyle
\left|k\right> = \bigotimes_{j=0}^{n-1}\left|k_j\right>_j$ where $\displaystyle
k = \sum_{j=0}^{n-1} k_j\cdot 2^{n-j-1}$ and $k_j\in\{0,1\}$.
E.g. $\
\left|31\right>=
\left|0\right>_0
\left|0\right>_1
\left|0\right>_2
\left|1\right>_3
\left|1\right>_4
\left|1\right>_5
\left|1\right>_6
\left|1\right>_7$. Then the state vector of the system takes the form $\displaystyle
\left|\psi\right> = \sum_{k=0}^{N-1} \psi_k \left|k\right>$. Introduce the $N\times N$ unitary matrix $U=U_{N-1}\cdots U_1U_0$ where
$$
SU_kS^* = U_{k+1},\,
S = \begin{bmatrix} 0 & 1 \\I_{N-1} & 0 \end{bmatrix},\,
U_0 = \begin{bmatrix} Q & 0 \\ 0 & I_{N-2} \end{bmatrix},\,
Q = \begin{bmatrix}
\sqrt{1 - |\beta|^2} & \beta \\
-\beta^* & \sqrt{1 - |\beta|^2} \\
\end{bmatrix}
$$
$U_0$ is unitary due to $Q$ being unitary. $S$ is unitary being a permutation matrix. By induction, each $U_k$ is unitary. Thus, $U$ is unitary.
1. $\beta$ is taken to be $0.15$ for reasons unclear.
2. Initially, prepare the system in the state $\left|\psi\right> \leftarrow \left|0\right>$, i.e. all spins down.
4. Evolve according to $\left|\psi\right> \leftarrow U\left|\psi\right>$ for 199 times.
6. Do something, such as applying a measurement, or nothing.
7. Evolve according to $\left|\psi\right> \leftarrow U^*\left|\psi\right>$ for 199 times.
9. Plot the statistics of interest for all iterations. The author presented these.
#### Expected Spin of the $j^{\text{ th}}$ Electron
Given the state vector of the system being $\left|\psi\right>$ , the author defines it by
$$
s_j = \frac{1}{2} + \frac{1}{2}
\left<\psi\right|
\left(
I_{2^{j-1}}
\otimes\begin{bmatrix}
-1 & 0 \\
0 & 1 \\
\end{bmatrix}\otimes
I_{2^{n-j-1}}
\right)
\left|\psi\right>
$$
Note that the matrix in the middle differs from the conventional Pauli z matrix by a negation, as the author assigns the spin-down eigenvector to the first basis vector. Furthermore, the author applies extra scaling and translation by half.
#### Expected Total Magnetization
The author defines it by $\displaystyle
\mu = \displaystyle\sum_{j=0}^{n-1} s_j$ .
#### Entropy
The definition given by the author follows that of Von Neumann. Omitting Boltzmann's constant, $\displaystyle S_{\Gamma}
= -\sum_{\gamma\in\Gamma}
\operatorname{Tr}(\rho P_{\gamma})
\log\left(\frac
{\operatorname{Tr}(\rho P_{\gamma})}
{\operatorname{Tr}(P_{\gamma})}
\right)$ where $\Gamma$ is a partition over $[N]$ and $\displaystyle
P_{\gamma} = \sum_{k\in\gamma}\left|k\right>\left<k\right|$. If $\Gamma$ is the finest partition over $[N]$, the author calls it the microscopic entropy; otherwise, it's called a macroscopic entropy.

### Simple Simulation
<!-- [Click here for the plot in action.](https://www.csie.ntu.edu.tw/~b05902035/simulation.html) -->
[Click here for the plot in action.](https://observablehq.com/d/01b2cf913e5863fc)

Here is an implementation in javascript which runs in modern browsers.

function simulate(n, beta, halfIterations) {
    const N = 1 << n;

Let $\alpha = \sqrt{1-|\beta|^2}$ and $k\in[N]$. Then,
$\displaystyle
U_k \left|\psi\right>
= \sum_{k'\in[N]}\psi_{k'}U_k\left|k'\right>
= \left(
\alpha\,\psi_k+\beta\,\psi_{k+1}
\right) \left|k\right>
+\left(
\alpha\,\psi_{k+1}-\beta^*\psi_k
\right) \left|k+1\right>
+\sum_{\substack{k'\in[N]} \\ k'\neq k,k+1}\psi_{k'}\left|k'\right>$

    const U = (beta) => {
        const alpha = Math.sqrt(1 - beta * beta);
        return (k, psi) => {
            const k_ = (k + 1) % N;
            const x = psi[k], y = psi[k_];
            psi[k ] = alpha * x + beta * y;
            psi[k_] = alpha * y - beta * x;
        }
    };

Thus, $U\left|\psi\right> = U_{N-1}\cdots U_1U_0\left|\psi\right>$ yields

    const evolve = (psi) => {
        const u = U(beta);
        psi.forEach((_, k) => u(k, psi));
    };

Similary, $\displaystyle
U^*\left|\psi\right> = U_0^*U_1^*\cdots U_{N-1}^*\left|\psi\right>$ and $\begin{bmatrix}\alpha & \beta \\ -\beta^* & \alpha\end{bmatrix}^*
= \begin{bmatrix}\alpha & -\beta^* \\ \beta & \alpha\end{bmatrix}$ give

    const reverse = (psi) => {
        const u = U(-beta);
        psi.forEach((_, k) => u(N - k - 1, psi));
    };

$$
\begin{align}
& \left<\psi\right|
\left(
I_{2^{j-1}}
\otimes\begin{bmatrix}
-1 & 0 \\
0 & 1 \\
\end{bmatrix}\otimes
I_{2^{n-j-1}}
\right)
\left|\psi\right>
\\=\;&
\left(
\sum_{k_1,k_2,k}
\psi_{k_1,k,k_2}^*
\left<k_1\right|\left<k\right|\left<k_2\right|
\right)
\left(
I_{2^{j-1}}
\otimes\begin{bmatrix}-1 & 0 \\0 & 1 \end{bmatrix}\otimes
I_{2^{n-j-1}}
\right)
\left(
\sum_{k_1',k_2',k'}
\psi_{k_1',k',k_2'}
\left|k_1'\right>\left|k'\right>\left|k_2'\right>
\right)
\\=\;&
\sum_{k_1,k_1'}\sum_{k_2,k_2'}\sum_{k,k'}
\psi_{k_1,k,k_2}^* \psi_{k_1',k',k_2'}
\left<k_1\right|\left.k_1'\right>
\left<k\right|\begin{bmatrix}-1 & 0 \\0 & 1 \end{bmatrix}\left|k'\right>
\left<k_2\right|\left.k_2'\right>
\\=\;&
\sum_{k_1,k_2}\sum_{k,k'}
\psi_{k_1,k,k_2}^* \psi_{k_1,k',k_2}
(-1)^{k+1}\left<k\right|\left.k'\right>
\\=\;&
\sum_{k_1,k_2} \left(
\psi_{k_1,1,k_2}^* \psi_{k_1,1,k_2} - \psi_{k_1,0,k_2}^* \psi_{k_1,0,k_2}
\right)
\\=\;&
\sum_{k_j=1}|\psi_k|^2 - \sum_{k_j=0}|\psi_k|^2
\end{align}
$$

    const spin = (psi) => Array.from({length: n}, (_, j) => {
        const mask = 1 << (n - j - 1);
        return psi.reduce((s, p, k) => s + (k & mask ? p * p : -p * p), 0.0);
    });

This is the actual simulation.

    const psi = Array(N).fill(0.0).fill(1.0, 0, 1);
    return [spin(psi),
        ...Array.from({length: halfIterations}, (_) => (evolve(psi), spin(psi))),
        ...Array.from({length: halfIterations}, (_) => (reverse(psi), spin(psi))),
    ];
}

The plot is powered by `Plotly,js`.

function plot(z)
{
    const addPlot = (name, data, layout) => {
        document.body.insertAdjacentHTML('beforeend', `<div id="${name}""></div>`);
        Plotly.newPlot(name, [data,], layout);
    };

    addPlot('spin', {
        z: z,
        type: 'surface',
        contours: {x: {show: true}},
    }, {
        title: 'Expected Spin of Each Electron',
        width: 900,
        height: 600,
        scene: {
            aspectmode: 'manual',
            aspectratio: {x:2, y:3, z:1},
            camera: {
                up: {x:0, y:0, z:1},
                center: {x:0, y:0, z:-0.8},
                eye: {x:5/2.2, y:-4/2.2, z:3/2.2},
            },
        },
    });

    addPlot('magnetization', {
        y: z.map((eight) => eight.reduce((s, v) => s + v, 0.0)),
        type: 'scatter',
    }, {
        title: 'Expected Total Magnetization',
        width: 900,
    });
}

Putting everything together.

plot(simulate(8, 0.15, 199));

The simulation above serves as the basis for comparison for more interesting simulations where particular measurements are made at the $200^{\text{ th}}$ iteration. To better study the effect of measurements, the author maintains an ensemble of identically prepared systems instead of a specific system.
* The distribution of states of an ensemble of identically prepared systems after a measurment is exactly the propability distribution of the outcome of a measurement on a single system. If there is a way of maintaining such an ensemble, one can obtain much more intriguing data, such as the averaged expected total magnetization of said ensemble, efficiently compared to the uninteresting collapse of a specific single system.
* The author suggests the density operator approach.
* In general, the density operator will be in a mixed state after a measurement, i.e. $\
\nexists\left|\psi\right>\in
\operatorname{span}\{\left|k\right>\}_{k=0}^{N-1}$ s.t. $\rho=\left|\psi\right>\left<\psi\right|$. Thus, my previous approach, maintaining a state vector, is no longer feasible. Instead, the whole density operator must be maintained.

### Density Operator Approach
Let $n = 8$, $N = 2^n = 256$, and denote the up-spin and down-spin of the $j^{\text{ th}}$ electron by $\left|1\right>_j$ and $\left|0\right>_j$ respectively where $j\in[n]$. Let $\{\left|k\right>\}_{k\in[N]}$ be the ordered basis for the composite system such that $\displaystyle \left|k\right> = \bigotimes_{j=0}^{n-1}\left|k_j\right>_j$ where $\displaystyle k = \sum_{j=0}^{n-1} k_j\cdot 2^{n-j-1}$ and $k_j\in\{0,1\}$. Then the state takes the form $\displaystyle \rho = \sum_{i,j=0}^{N-1}\rho_{ij}\left|i\right>\left<j\right|$. Introduce $\displaystyle U=\prod_{k=0}^{N-1}U_{N-k-1}$ s.t. $SU_kS^* = U_{k+1}$ where
$$
S = \begin{bmatrix} 0 & 1 \\I_{N-1} & 0 \end{bmatrix},\,
U_0 = \begin{bmatrix} Q & 0 \\ 0 & I_{N-2} \end{bmatrix},\,
Q = \begin{bmatrix}
\sqrt{1-|\beta|^2} & \beta \\
-\beta^* & \sqrt{1-|\beta|^2}
\end{bmatrix}
$$
$\prod_{k=0}^{m-1}U_k = S^m(S^*U_0)^m \implies \prod_{k=0}^{m}U_k = S^mU_0(S^*)^mS^m(S^*U_0)^m = S^{m+1}(S^*U_0)^{m+1}$ and $U_1U_0=SU_0S^*U_0=S^2(S^*U_0)^2$ yields $U = S^N(S^*U_0)^N = (S^*U_0)^N$. This form is more concise but more computationally expensive. In this case,
$$
U = \begin{bmatrix}
-\beta^* & \sqrt{1-|\beta|^2} & 0 \\
0 & 0 & I_{N-2} \\
\sqrt{1-|\beta|^2} & \beta & 0 \\
\end{bmatrix}^N
$$
1. Initially, prepare the system where all spins are downward, i.e. $\rho \leftarrow \left|0\right>\left<0\right|$.
2. Evolve according to $\rho \leftarrow U\rho U^*$ for 199 times. Since the system are in a pure states, calculations are identical that of the previous section.
4. Apply a measurement $M$, i.e. $\displaystyle \rho \leftarrow \sum_m \left<m\right|\rho\left|m\right> \left|m\right>\left<m\right|$. Before the measurement, the system is in a pure state $\rho = \left|\psi\right>\left<\psi\right|$ where $\displaystyle \left|\psi\right> = \sum_{k=0}^{N-1}\psi_k\left|k\right>$. $\displaystyle
Z^{\otimes n}\left|k\right>
= \bigotimes_{j=0}^{n-1} Z\left|k_j\right>
= \bigotimes_{j=0}^{n-1} (-1)^{k_j + 1}\left|k_j\right>
= (-1)^{\sum_{j=0}^{n-1}k_j + n}\left|k\right>$, implies the basis vectors are the eigenvector of the complete measurement of spins, and the state will become
$$
\rho
\leftarrow \sum_{k=0}^{N-1}
\left<k\right|\left.\psi\right>\left<\psi\right|\left.k\right>
\left|k\right>\left<k\right|
= \sum_{k=0}^{N-1}
|\psi_k|^2 \left|k\right>\left<k\right|
$$
6. Evolve according to $\rho \leftarrow U\rho U^*$ for 199 times.
$$
U_j\rho U_j^*
= \begin{bmatrix}
I & 0 & 0 \\
0 & Q & 0 \\
0 & 0 & I \\
\end{bmatrix}
\begin{bmatrix}
R_{11} & R_{12} & R_{13} \\
R_{21} & R_{22} & R_{23} \\
R_{31} & R_{32} & R_{33} \\
\end{bmatrix}
\begin{bmatrix}
I & 0 & 0 \\
0 & Q^* & 0 \\
0 & 0 & I \\
\end{bmatrix}
= \begin{bmatrix}
R_{11} & R_{12}Q^* & R_{13} \\
QR_{21} & QR_{22}Q^* & QR_{23} \\
R_{31} & R_{32}Q^* & R_{33} \\
\end{bmatrix}
$$
Only, $8N$ entries requires update.
7. Plot the statistics of interest for all iterations. The author presented these.