Transverse Field Ising Model
Index
- Transverse Field Ising Model
The goal of this tutorial will be to implement Transverse Field Ising Models (TFIM) using QURI-SDK. We will approximate the ground state of the TFIM Hamiltonian using VQE and further improving it using SPE, QSCI. Finally we will calculate the magnetization and magnetic susceptibility. The magnetic susceptibility diverges in a way that is consistent with the expected critical behaviour in the thermodynamic limit near the critical point.
Introduction
The Ising Model is a theoretical model to explain ferromagnetism in solids. It assumes only nearest neighbor interaction between lattice spins. For details, see Ising model.
Over time, Ising models have been applied to various other fields such as Machine Learning (Hopfield networks), Biology (Protein folding), etc. The Transverse Field Ising Model (TFIM), which is the main part of this notebook, is a quantum version of the Ising model, where the correlation between spins is affected by a transverse magnetic field. This model is used to study quantum phase transitions and has applications in Quantum Computing and Condensed Matter Physics.
From a Computer Science perspective, NP-hard problems, including all of Karp's 21 NP-complete problems can be mapped to the Ising model with a polynomial overhead. Thus, any NP-hard problem can be reduced to finding the ground state of the Ising model. For further reference, look at Ising formulations of NP-Complete problems.
The Hamiltonian of the Transverse Field Ising Model is given by:
- : The Pauli matrix (for spin measurement)
- : The Pauli matrix (for transverse field)
- Ferromagnetic coupling (): Favors alignment of neighboring spins along X.
- Transverse field (): Induces quantum fluctuations, polarizing spins along J.
- The boundary condition is periodic, meaning the last spin interacts with the first spin.
Note: In most references, the Hamiltonian will have Pauli matrices in the opposite order, i.e., . This is due to the convention of defining the spin along Z axis in Condensed Matter Physics. The reason we switch Pauli X and Z is due to the Z2 ansatz we will be using for the VQE. Everything else remains the same, if you are doing any calculation wrt to the standard TFIM Hamiltonian, just switch the Pauli matrices.
Ground States
h >> J (Spin Polarised Phase):
Spin-polarized states refers to states where electron spins in a material are not randomly oriented but rather exhibit a net alignment due to an applied field, leading to a non-zero magnetization. When is much larger than , the spins tend to align along the Z-axis (spin-polarized states). One of the fully spin-polarised states is:
The energy eigenvalue is and magnetization = , where is the number of spins.
J >> h (Ferromagnetic Phase):
When is much larger than , the spins tend to be aligned along the X-axis, forming an equal superposition of and . The two ground states are:
Formally, any superposition of these two states is also a ground state. But in the thermodynamic limit , such symmetric superpositions are not stable. To tunnel between the two states requires flipping spins, so the tunneling amplitude goes to zero exponentially with N. As a result, the system chooses one of the two magnetization directions even at zero temperature. This is known as Spontaneous Symmetry breaking.
Symmetry
The symmetry corresponds to a global spin-flip operation: which satisfies , forming the cyclic group .
Interaction term :
Transverse field term :
The Hamiltonian is invariant under , confirming symmetry.
Setting up the QubitHamiltonian in QURI
import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence
from quri_parts.core.operator import Operator, pauli_label, get_sparse_matrix
def construct_TFIM_hamiltonian(J: float, h: float, n_qubits: int) -> Operator:
hamiltonian = Operator()
# Add Ising interaction terms (-J Σ X_i X_{i+1}) (Assuming periodic boundary conditions)
for i in range(n_qubits):
pauli_index = pauli_label('X'+str(i)+' X'+str((i+1)%n_qubits))
hamiltonian.add_term(pauli_index, -J)
# Add transverse field terms (-h Σ Z_i)
for i in range(n_qubits):
pauli_index = pauli_label('Z'+str(i))
hamiltonian.add_term(pauli_index, -h)
return hamiltonian
# Parameters
n_qubits = 12 # Number of spins/qubits
J = 1.0 # Interaction strength
h = 1.0 # Transverse field strength
hamiltonian = construct_TFIM_hamiltonian(J, h, n_qubits)
Since n_qubits is small, we can find the smallest eigenvalue (ground state energy) of the Hamiltonian matrix via exact diagonalization. This is only performed for benchmarking. For relistic physical systems involving large n_qubits, this computation is not feasible and we have to rely on Quantum Computing methods to find the approximate ground state.
## Ground eigenvalue and eigenvector computation using numpy (Will be later used for benchmarking)
vals, vecs = np.linalg.eigh(get_sparse_matrix(hamiltonian).toarray())
EXACT_GS_ENERGY = np.min(vals)
EXACT_GAP = vals[np.argsort(vals)][:2] @ [-1, 1]
print("E_0:", EXACT_GS_ENERGY)
print("Delta (E_1 - E_0)):", EXACT_GAP)
E_0: -15.322595151080774
Delta (E_1 - E_0)): 0.13108692563047342
VQE Optimization
In case you're new to QURI SDK's implementation of Variational algorithms, you may refer to this tutorial Variational Algorithms.
Cost Function and Gradient Estimators
We will first define closures (higher-order functions) that takes the Hamiltonian and parametric quantum state as the parameter and returns cost function and gradient estimators for the VQE optimization. The cost estimator
computes the expectation value of the Hamiltonian with respect to the parametric state, while the gradient estimator
computes the gradient of the expectation value with respect to the parameters of the parametric state.
from quri_parts.core.estimator.gradient import create_parameter_shift_gradient_estimator
from quri_parts.qulacs.estimator import create_qulacs_vector_concurrent_parametric_estimator, create_qulacs_vector_parametric_estimator
from quri_parts.core.state import ParametricCircuitQuantumState
from typing import Callable
import numpy.typing as npt
estimator = create_qulacs_vector_parametric_estimator()
def cost_fn_estimator(hamiltonian : Operator, parametric_state: ParametricCircuitQuantumState) -> Callable[[Sequence[float]], Sequence[float]]:
return lambda param_values : estimator(hamiltonian, parametric_state, param_values).value.real
concurrent_parametric_estimator = create_qulacs_vector_concurrent_parametric_estimator()
gradient_estimator = create_parameter_shift_gradient_estimator(concurrent_parametric_estimator)
def grad_fn_estimator(hamiltonian : Operator, parametric_state: ParametricCircuitQuantumState) -> Callable[[Sequence[float]], npt.NDArray[np.float64]]:
return lambda param_values: np.real(gradient_estimator(hamiltonian, parametric_state, param_values).values)
Setting up VQE Optimization loop
Next, we will set up the VQE optimization loop essentially updating the parameters of the parametric state using the Gradient descent method and a momentum based optimizer (e.g. Adam). QURI Parts provides a convenient OptimizerStatus
class that will handle the optimizer convergence for us.
The below vqe
function takes the parametric circuit (Ansatz), hamiltonian, initial parameters, optimizer and returns the optimzied ground state, optimized params , log of ground state energy, and the number of iterations taken to converge.
from quri_parts.algo.optimizer import Adam, Optimizer, OptimizerStatus
from quri_parts.circuit import ImmutableLinearMappedParametricQuantumCircuit
from quri_parts.core.state import quantum_state, apply_circuit
def vqe(
parametric_circuit : ImmutableLinearMappedParametricQuantumCircuit,
hamiltonian : Operator,
init_params: Sequence[float] = None,
optimizer : Optimizer = Adam()
):
# Setting up the parametric state
cb_state = quantum_state(n_qubits, bits=0)
parametric_state = apply_circuit(parametric_circuit, cb_state)
# Defining the cost and gradient estimators based on the hamiltonian and parametric state
cost_fn = cost_fn_estimator(hamiltonian, parametric_state)
grad_fn = grad_fn_estimator(hamiltonian, parametric_state)
# Intializing to random parameters if params not provided
if init_params is None:
init_params = np.random.random(parametric_circuit.parameter_count)
opt_state = optimizer.get_init_state(init_params)
energy = []
# Running the VQE Optimization loop
while True:
opt_state = optimizer.step(opt_state, cost_fn, grad_fn)
energy.append(opt_state.cost)
if opt_state.status == OptimizerStatus.FAILED:
print("Optimizer failed")
break
if opt_state.status == OptimizerStatus.CONVERGED:
print("Optimizer converged")
break
bound_state = parametric_state.bind_parameters(opt_state.params)
return bound_state, opt_state.params, energy, opt_state.niter