<jupyter_start><jupyter_text>Entanglement in the Stern-Gerlach ExperimentIn this problem we want to consider Stern-Gerlach experiment with a more realistic approach.Assume the electrons are shot towards the apparatus. The hamiltonian is as follows:$H = \frac{P_z^2}{2m} - \mu \lambda z \sigma_z$where the second term is from the interaction of the spin of the electron with the linear magnetic field.The state of the system is a member of the tensor product space of spin and z-space. Let's make a concretez-space of 501 sites with $a = 0.01$ where he middle site is for $z=0$. Then, the state is a vector with 1002coefficients.For making the evolution matrix let's take the time steps to be $\epsilon = 0.01$ and $\hbar = 1$. Also let's assume $\mu \lambda = 1, m = 1$.<jupyter_code># constants
mu_lambda = 0.05
n = 501
a = 0.01
m = 1
epsilon = 0.01

# the z-space operator
z_op = np.asmatrix(np.kron(np.diag(np.linspace(-2.5, 2.5, 501)), np.eye(2)))

# spin operator
z_pauli = np.kron(np.matrix([[1, 0], [0, -1]]), np.eye(n))
z_pauli
## making the momentum operator
# making shift_op
shift_op = np.eye(n, k=1)
shift_op[n-1, 0] = 1
shift_op = np.asmatrix(shift_op)

# diagonalizing the shift_op
vals, vecs = np.linalg.eig(shift_op)

# making P
P = vecs * np.diag(np.angle(vals)) * vecs.H / a
P = np.kron(P, np.eye(2))
P.shape
H = P ** 2 / (2 * m) - mu_lambda * z_op * z_pauli
H.shape<jupyter_output><empty_output><jupyter_text>Now we have the hamiltonian, let's generate the evolution operator $U = \exp(-i H \epsilon$, taking $\hbar=1$.<jupyter_code>evolution = expm(-1j * H * epsilon)
evolution = np.asmatrix(evolution)
evolution<jupyter_output><empty_output><jupyter_text>Section ALet's take the initial state to be at $z = 0$ and spin up. Let's find the evolution to $t=10$<jupyter_code># making the z=0 state in z-space
z_0 = np.zeros((501, 1))
z_0[250, 0] = 1
z_0 = np.asmatrix(z_0)

# making the spin up state in spin space
spin_up = np.array([[1], [0]])

# making the initial state of the system.
init_state = np.kron(z_0, spin_up)

init_state.shape
# now to evolve the system to t = 10 = 1000 * epsilon
state = init_state

for _ in range(1000):
    state = evolution * state
    
state<jupyter_output><empty_output><jupyter_text>Section BNow for state $z=0$ and spin down. We follow the same procedure:<jupyter_code># making spin down in the spin space
spin_down = np.matrix([[0], [1]])


# making initial state
init_state = np.kron(z_0, spin_down)
init_state.shape
# now we evolve the system
state = init_state

for _ in range(1000):
    state = evolution * state

state<jupyter_output><empty_output><jupyter_text>Section CNow we try $z = 0$ and spin up in the x axis. In the spin space, this state is given as:$|x, +> = \frac{1}{\sqrt{2}} (1, 1)^t$<jupyter_code># making spin +x
spin_x_up = np.matrix([[1], [1]]) / (2 ** 0.5)

# making the initial state
init_state = np.kron(z_0, spin_x_up)
init_state.shape
# evolve the state to t = 10
state = init_state

for _ in range(1000):
    state = evolution * state
    
state
plt.plot(np.linspace(0, 1001, 1002), np.absolute(np.array(state)[:, 0]) ** 2)
plt.show()<jupyter_output><empty_output><jupyter_text>Now let's check to see if this final state is entangled.In order to do this, we need to generate the reduced density operator of the state and take the trace of $\rho^2$.If the trace is less than 1, then the state is entangled.$\hat{\rho} = |\psi><\psi|, tr(\hat{\rho}^2) \leq 1$Let's take a look at our state right now. The first sub-space is for z-space which has a dimension of 501 and the second is the spin sub-space which has a dimension of 2.There's a technique with which we can generate the reduced density matrix; and that is the `reshape` method of the `matrix` and `array` classes in `numpy`. We need to reshape the matrix from $1002 \times 1$ to $2 \times 501$ and take its transpose. according to the `reshape` algorithm by `numpy`, if we multiply this ket with its bra, we get the reduced density matrix.$|\psi'> = `state`, \rho' = |\psi><\psi| = `state` * `state.H`$<jupyter_code># generating the reduced density operator
altered_state = state.reshape((2, 501))    # reshaping the state

reduced_rho = altered_state.T * altered_state.conjugate()

reduced_rho<jupyter_output><empty_output><jupyter_text>Now we need to take the trace of the square of the reduced density operator. If this parameter is significantlyless than 1, then the state is entangled.<jupyter_code>param = np.trace(reduced_rho ** 2)
param<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11775121735281098, 'fraction_numerical': 0.03142983621071271, 'mean_word_length': 4.457729468599034, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '30228', 'license': 'MIT', 'n_tokens_mistral': 1684, 'n_tokens_neox': 1552, 'n_words': 578, 'path': 'quantum_numerical_analysis/module4/stern_gerlach_entanglement.ipynb', 'repo_name': 'ali-manhani/fun-simulations', 'size': '4518'}