Simulating the theory on a quantum computer
In [1] and [2], algorithms for simulating lattice scalar field lattice on a quantum computer was introduced. This can serve as a first example of simulating general lattice field theory on quantum computers.
The Hamiltonian under consideration is given by:
Here and satisfy the commutation relation:
To simulate it on a quantum computer, we need to discretize it so that the lattice Hamiltonian takes the form:
The lattice field operators satisfy:
Rescaling with and makes the field and the conjugate momentum dimensionless.
def get_qubit_label(in_site_label: int, site_label: int, n_phi_qubit: int) -> int:
"""Converts a local qubit index to the full system's qubit index.
in_site_label: The qubit label on the local site.
site_label: The label of the site.
n_phi_qubit: Number of qubits assigned to each site.
"""
return n_phi_qubit * site_label + in_site_label
In QURI Algo, we can define a subclass of Problem
containing the parameters that characterize the system. So, we define DiscreteScalarField1D
that represents the discrete scalar field Hamiltonian in 1 spatial dimension.
from quri_algo.problem import HamiltonianInput
from dataclasses import dataclass, field
import numpy as np
@dataclass
class DiscreteScalarField1D(HamiltonianInput):
"""Represents the Hamiltonian of the 1D scalar field:
Hamiltonian given by:
H = 1/2 Π_j^2 + 1/2 mb^2 Φ_j^2 + 1/2 (Φ_j - Φ_{j+1})^2 + λ/4! Φ_j^4 + J_j Φ.
Note:
1D here means 1 spatial dimension.
Args:
n_state_qubit
n_discretize: Number of points discretizing the field.
n_phi_qubit: Number of qubits per site.
mb: Boson mass.
lam: coupling costant of the phi^4 term.
external_field: External field strength J.
"""
n_state_qubit: int = field(init=False)
n_discretize: int
n_phi_qubit: int
mb: float
lam: float
external_field: float = 0.0
def __post_init__(self) -> None:
self.n_state_qubit = self.n_discretize * self.n_phi_qubit
@property
def n_phi_dimension(self) -> int:
return 2**self.n_phi_qubit
@property
def delta_phi(self) -> float:
return np.sqrt(2 * np.pi * self.mb / self.n_phi_dimension)
Discrete scalar field
The discrete scalar field is designed to satisfy the quantization condition:
where with being the number of qubits assigned to site . The site index will be suppressed from now and we adopt the conventions:
The discrete field operator on the site can be expressed as:
which can be implemented as a QURI Parts Operator
.
from quri_parts.core.operator import Operator, pauli_label
def get_scalar_field_operator(
site_label: int, n_phi_qubit: int, mass: float
) -> Operator:
phi = Operator({})
delta_phi = np.sqrt(2 * np.pi * mass / 2**n_phi_qubit)
for q in range(n_phi_qubit):
coeff = - delta_phi * 2**(q) / 2
l = get_qubit_label(q, site_label, n_phi_qubit)
phi.add_term(pauli_label(f"Z {l}"), coeff)
return phi
Example: Check that the operator satisfies the field operator quantization condition.
from quri_parts.core.state import quantum_state
from quri_parts.qulacs.estimator import create_qulacs_general_vector_estimator
site_label = 0
n_qubits = 4
mb = 1
estimator = create_qulacs_general_vector_estimator()
field_op = get_scalar_field_operator(site_label, n_qubits, mb)
b = 0b0101
estimator(field_op, quantum_state(n_qubits, bits=b)), (b - (2**n_qubits - 1)/2) * np.sqrt(2 * np.pi * mb / 2**n_qubits)
(_Estimate(value=(-1.5666426716443749+0j), error=0.0), -1.566642671644375)
Discrete conjugate momentum
The discrete conjugate momentum is defined as
The terms is defined as ( is suppressed.)
The eigenstate of the conjugate momentum operator is defined as:
where is the eigenstate of the field operator so that
The transformation is harder to be represented as a qubit operator. Instead, it can be represented as
with and QFT corresponds to the quantum Fourier transform:
Quantum Fourier Transform
The quantum Fourier transform can be performed by the following circuit
where
and the controlled U gate can be decomposed into CNOTs and U gates as:
Here denotes the controlled qubit and denotes the target qubit. This can be implemented in QURI Parts as:
from quri_parts.circuit import QuantumCircuit, NonParametricQuantumCircuit, ImmutableBoundParametricQuantumCircuit
import numpy as np
def add_controlled_U1_gate(
circuit: QuantumCircuit, control: int, target: int, angle: float
) -> None:
circuit.add_U1_gate(control, angle/2)
circuit.add_CNOT_gate(control, target)
circuit.add_U1_gate(target, -angle/2)
circuit.add_CNOT_gate(control, target)
circuit.add_U1_gate(target, angle/2)
Then the QFT circuit can be implemented with:
def create_qft_gate(qubit_count: int) -> ImmutableBoundParametricQuantumCircuit:
circuit = QuantumCircuit(qubit_count)
for i in range(qubit_count//2):
circuit.add_SWAP_gate(i, qubit_count-i-1)
for target in range(qubit_count):
circuit.add_H_gate(target)
for l, control in enumerate(range(target+1, qubit_count)):
angle = 2 * np.pi/2**(l+2)
add_controlled_U1_gate(circuit, control, target, angle)
return circuit.freeze()
Example: execute the quantum Fourier transform
from quri_parts.qulacs.simulator import evaluate_state_to_vector
n_qubits = 4
n = 2**n_qubits
QFT = create_qft_gate(n_qubits)
for b in range(n):
qft_state = quantum_state(n_qubits, circuit=QFT, bits=b)
qft_state_vector = evaluate_state_to_vector(qft_state).vector
expected_qft_vector = np.array([np.exp(2j * np.pi / n * i * b) for i in range(n)]) / np.sqrt(n)
assert np.allclose(qft_state_vector, expected_qft_vector)
The linear rotation term
Next, implement the linear rotation terms: .
def create_linear_rotation(qubit_count: int) -> ImmutableBoundParametricQuantumCircuit:
circuit = QuantumCircuit(qubit_count)
n = 2**qubit_count
delta = (n - 1) * np.pi / n
for q in range(qubit_count):
angle = - 2**q * delta
circuit.add_RZ_gate(q, angle)
return circuit.freeze()
Example: the linear rotation term
import itertools
n_qubits = 4
n = 2**n_qubits
linear_rotation = create_linear_rotation(n_qubits)
delta = (n - 1) / n * np.pi
for b1, b2 in itertools.product(range(n), repeat=2):
linear_rotation_state_1 = quantum_state(n_qubits, circuit=linear_rotation, bits=b1)
linear_rotation_vector_1 = evaluate_state_to_vector(linear_rotation_state_1).vector
linear_rotation_state_2 = quantum_state(n_qubits, circuit=linear_rotation, bits=b2)
linear_rotation_vector_2 = evaluate_state_to_vector(linear_rotation_state_2).vector
n1 = linear_rotation_vector_1[linear_rotation_vector_1.nonzero()]
n2 = linear_rotation_vector_2[linear_rotation_vector_2.nonzero()]
assert np.isclose(
n1/n2, np.exp(-1j * delta * (b1 - b2))
)
The operator
Finally, we can combine the functions in the last 2 subsections and finally implement the full circuit for the operation.
def create_F_gate(qubit_count: int) -> NonParametricQuantumCircuit:
circuit = QuantumCircuit(qubit_count)
qft = create_qft_gate(qubit_count)
linear_rotation = create_linear_rotation(qubit_count)
circuit.extend(linear_rotation)
circuit.extend(qft)
circuit.extend(linear_rotation)
return circuit
Example: full operator
from quri_parts.qulacs.simulator import evaluate_state_to_vector
n_qubits = 4
N = 2**n_qubits
const = (N - 1)/2
F = create_F_gate(n_qubits)
b = 0b1001
conj_eig_state = quantum_state(n_qubits, circuit=F, bits=b)
conj_eig_state_vector = evaluate_state_to_vector(conj_eig_state).vector
expected = np.zeros(2**n_qubits, dtype=np.complex128)
for i in range(N):
expected[i] = np.exp(
2j * np.pi / N * (b - const) * (i - const)
) / np.sqrt(N)
expected / conj_eig_state_vector # Only off by overall phase
array([-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j,
-0.99518473+0.09801714j, -0.99518473+0.09801714j])
The Exact Hamiltonian
Here we build the functions for generating the matrix representations of , , and the hamiltonian and the local hamiltonian . These will be used extensively in later sections.
import math
import functools
import numpy.typing as npt
from quri_parts.core.operator import get_sparse_matrix
def get_full_matrix(
n_sites: int, phi_qubit: int, on_site_matrix: npt.NDArray[np.complex128], site_idx: int
) -> npt.NDArray[np.complex128]:
"""Expand an operator on a local site to a matrix of the whole system.
"""
full = [np.eye(2 ** phi_qubit, dtype=np.complex128) for _ in range(n_sites)]
full[n_sites - 1 - site_idx] = on_site_matrix
return functools.reduce(lambda a, b: np.kron(a, b), full)
def get_F_matrix(n_qubits: int) -> npt.NDArray[np.complex128]:
"""The matrix representation of the F operator on a single site.
"""
dim = 2**n_qubits
f_matrix = np.zeros((dim, dim), dtype=np.complex128)
for a, b in itertools.product(range(dim), repeat=2):
f_matrix[a, b] = np.exp(2j * np.pi/dim * (a - (dim - 1)/2) * (b - (dim - 1)/2))
return f_matrix/np.sqrt(dim)
def get_phi_matrix(n_qubits: int, mass: float) -> npt.NDArray[np.complex128]:
"""The matrix representation of the field operator on a single site.
"""
op = get_scalar_field_operator(0, n_qubits, mass)
return get_sparse_matrix(op, n_qubits).toarray()
def get_hamiltonian_matrix(system: DiscreteScalarField1D) -> npt.NDArray[np.complex128]:
"""Matrix representation of full Hamiltonian.
"""
n_qubits = system.n_phi_qubit
n_site = system.n_discretize
mb = system.mb
J = system.external_field
lam = system.lam
F_matrix = get_F_matrix(n_qubits)
phi = get_phi_matrix(n_qubits, mb)
full_dim = 2**system.n_state_qubit
hamiltonian = np.zeros((full_dim, full_dim), dtype=np.complex128)
_phi_cache = {
site: get_full_matrix(n_site, n_qubits, phi, site) for site in range(n_site)
}
for site in range(n_site):
full_phi = _phi_cache[site]
full_F = get_full_matrix(n_site, n_qubits, F_matrix, site)
# conjugate momentum: mb^2 F Φ_j^2 F^{\dagger}
hamiltonian += 0.5 * mb**2 * (full_F @ full_phi @ full_phi @ full_F.conj().T)
# J Φ_j
hamiltonian += J * full_phi
# 1/2 mb^2 Φ_j^2
hamiltonian += 0.5 * mb**2 * (full_phi @ full_phi)
# lam/4! * Φ_j^4
hamiltonian += (lam/math.factorial(4)) * (full_phi @ full_phi @ full_phi @ full_phi)
# kin: 1/2 * (Φ_j - Φ_{j+1})^2
if site < n_site - 1:
kin_root = _phi_cache[site] - _phi_cache[site + 1]
hamiltonian += 0.5 * kin_root @ kin_root
return hamiltonian
def get_local_hamiltonian_matrix(system: DiscreteScalarField1D) -> npt.NDArray[np.complex128]:
"""Matrix representation of local Hamiltonian (on the zeroth site).
"""
n_qubits = system.n_phi_qubit
mb = system.mb
J = system.external_field
lam = system.lam
F_matrix = get_F_matrix(n_qubits)
phi = get_phi_matrix(n_qubits, mb)
dim = 2**n_qubits
hamiltonian = np.zeros((dim, dim), dtype=np.complex128)
# conjugate momentum: mb^2 F Φ_j^2 F^{\dagger}
hamiltonian += 0.5 * mb**2 * (F_matrix @ phi @ phi @ F_matrix.conj().T)
# J Φ_j
hamiltonian += J * phi
# 1/2 mb^2 Φ_j^2
hamiltonian += 0.5 * mb**2 * (phi @ phi)
# lam/4! * Φ_j^4
hamiltonian += (lam/math.factorial(4)) * (phi @ phi @ phi @ phi)
return hamiltonian
Local state preparation
As a demonstration, we prepare the ground state of the local Hamiltonian on site :
with VQE. We assume the ansatz to be the Hardware Efficient Ansatz . The full wave function can then be prepared by:
import math
from quri_parts.algo.ansatz import HardwareEfficient
from quri_parts.algo.optimizer import CostFunction, Params
from quri_parts.qulacs.estimator import create_qulacs_general_vector_estimator
from quri_parts.core.state import ParametricCircuitQuantumState
from quri_parts.circuit import inverse_circuit
estimator = create_qulacs_general_vector_estimator()
def get_cost_function(system: DiscreteScalarField1D, n_layers: int) -> tuple[CostFunction, ParametricCircuitQuantumState]:
n_local_qubit = system.n_phi_qubit
F = create_F_gate(n_local_qubit)
ansatz = HardwareEfficient(n_local_qubit, n_layers)
ansatz_state = ParametricCircuitQuantumState(n_local_qubit, circuit=ansatz)
scalar_field_op = get_scalar_field_operator(0, n_local_qubit, system.mb)
def cost_function(param: Params) -> float:
bound_state = ansatz_state.bind_parameters(param)
transformed_state = bound_state.with_gates_applied(inverse_circuit(F))
mb = system.mb
lam = system.lam
J = system.external_field
# Π^2 = m^2 F Φ^2 F†
# <Ψ| Π^2 |Ψ> = m^2 <Ψ|F Φ^2 F†|Ψ>
pi_sq = estimator(0.5 * mb**2 * scalar_field_op * scalar_field_op, transformed_state)
# 1/2 m^2 Φ^2 + lambda/4! Φ^4 + J Φ
potential_op = 0.5 * mb**2 * scalar_field_op * scalar_field_op
potential_op += (lam / math.factorial(4)) * scalar_field_op * scalar_field_op * scalar_field_op * scalar_field_op
potential_op += J * scalar_field_op
non_dyn = estimator(potential_op, bound_state)
cost = pi_sq.value + non_dyn.value
return cost.real
return cost_function, ansatz_state
Running the optimization loop
from quri_parts.algo.optimizer import NFT, OptimizerState, Optimizer, OptimizerStatus
N_DISCRETIZE = 2
N_LOCAL_QUBIT = 4
MASS = 1.0
LAM = 1.0
system = DiscreteScalarField1D(N_DISCRETIZE, N_LOCAL_QUBIT, MASS, LAM,)
n_layers = 2
def vqe(cost: CostFunction, optimizer: Optimizer, init_param: Params) -> OptimizerState:
it = 0
opt_state = optimizer.get_init_state(init_param)
while opt_state.status != OptimizerStatus.CONVERGED:
opt_state = optimizer.step(opt_state, cost)
it += 1
print(f"{it}-th iteration", opt_state.cost)
return opt_state
cost_fn, ansatz_state = get_cost_function(system, n_layers)
init_param = np.random.random(ansatz_state.parametric_circuit.parameter_count)
opt_result = vqe(cost_fn, NFT(), init_param)
1-th iteration 2.105098883857522
2-th iteration 1.868918383240643
3-th iteration 1.785920126439606
4-th iteration 1.6897866940576605
5-th iteration 1.6180315895934418
``````output
6-th iteration 1.5492399998468112
7-th iteration 1.4792533028063701
8-th iteration 1.4032048719674437
9-th iteration 1.327220025127388
``````output
10-th iteration 1.2584036727909962
11-th iteration 1.1995078205325649
12-th iteration 1.14960996221769
13-th iteration 1.1071663286132931
14-th iteration 1.0714078889266458
``````output
15-th iteration 1.041975706094989
16-th iteration 1.0182410870036964
17-th iteration 0.9992336332605014
18-th iteration 0.9837289632463216
19-th iteration 0.9705962932444469
``````output
20-th iteration 0.9592350240621574
21-th iteration 0.949260497671048
22-th iteration 0.9405438744840359
23-th iteration 0.9329461466741175
24-th iteration 0.9261167237341027
``````output
25-th iteration 0.9197133546716935
26-th iteration 0.9135813783660476
27-th iteration 0.9076701134536004
28-th iteration 0.9019558081885084
``````output
29-th iteration 0.8964122536798823
30-th iteration 0.8909928673929839
31-th iteration 0.8856287064106325
32-th iteration 0.8802368468684536
``````output
33-th iteration 0.8747290552633421
34-th iteration 0.8690167958751889
35-th iteration 0.8630132220167882
36-th iteration 0.8566340942549339
37-th iteration 0.8497992634185597
``````output
38-th iteration 0.8424356902998342
39-th iteration 0.8344823655179884
40-th iteration 0.8258969246285873
41-th iteration 0.8166631164000205
42-th iteration 0.8067975911651573
``````output
43-th iteration 0.7963539610155999
44-th iteration 0.7854221209176897
45-th iteration 0.774121682533871
46-th iteration 0.7625899205251458
47-th iteration 0.7509662349539434
``````output
48-th iteration 0.7393760282821533
49-th iteration 0.7279167372504098
50-th iteration 0.7166478514929697
51-th iteration 0.7055857045192517
52-th iteration 0.6947030446911531
``````output
53-th iteration 0.6839329274698245
54-th iteration 0.6731762491847584
55-th iteration 0.6623123465343141
56-th iteration 0.6512126996551968
57-th iteration 0.6397589328206792
``````output
58-th iteration 0.6278676017850735
59-th iteration 0.6155245388632262
60-th iteration 0.6028285686224029
61-th iteration 0.5900353028078362
62-th iteration 0.5775761635544661
``````output
63-th iteration 0.5660156701815144
64-th iteration 0.5559251708335791
65-th iteration 0.5477086124427533
66-th iteration 0.5414763356570298
67-th iteration 0.5370473280814598
``````output
68-th iteration 0.5340641385062819
69-th iteration 0.5321317929752634
70-th iteration 0.5309096116084417
71-th iteration 0.5301435110328949
72-th iteration 0.5296606998480333
``````output
73-th iteration 0.5293504729442999
74-th iteration 0.529144568147051
75-th iteration 0.52900188510774
76-th iteration 0.5288980094033495
77-th iteration 0.528818518189742
``````output
78-th iteration 0.528754878990803
79-th iteration 0.5287020019393596
80-th iteration 0.5286568023683746
81-th iteration 0.5286173648841361
82-th iteration 0.5285824597538293
``````output
83-th iteration 0.5285512636178
84-th iteration 0.5285231980182197
85-th iteration 0.5284978356683727
86-th iteration 0.528474845644187
87-th iteration 0.5284539609658443
``````output
88-th iteration 0.5284349590977298
89-th iteration 0.5284176499407014
90-th iteration 0.5284018682042708
91-th iteration 0.5283874683695011
92-th iteration 0.5283743212104481
``````output
93-th iteration 0.5283623112750935
94-th iteration 0.5283513349755249
95-th iteration 0.5283412990798031
96-th iteration 0.5283321194807424
97-th iteration 0.5283237201644362
``````output
98-th iteration 0.5283160323295841
99-th iteration 0.5283089936252741
100-th iteration 0.5283025474847542
101-th iteration 0.5282966425392321
102-th iteration 0.5282912320993431
``````output
103-th iteration 0.5282862736948211
We can compare the VQE state against the exact ground state of the local Hamiltonian
from quri_algo.core.estimator.hadamard_test import shift_state_circuit
from quri_parts.qulacs.simulator import evaluate_state_to_vector
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import eigsh
local_state = ansatz_state.bind_parameters(opt_result.params)
local_state_vector = evaluate_state_to_vector(local_state).vector
local_hamiltonian_matrix = get_local_hamiltonian_matrix(system)
local_field_operator = get_scalar_field_operator(0, system.n_phi_qubit, system.mb)
local_field_operator_matrix = get_sparse_matrix(local_field_operator, system.n_phi_qubit)
sparse_local = csc_matrix(local_hamiltonian_matrix)
sparse_local.eliminate_zeros()
_, exact_local_gs_vec = eigsh(sparse_local, k=1, which="SA")
print(
"Overlap with exact local Hamiltonian GS vector:",
np.abs(exact_local_gs_vec.flatten().conj() @ local_state_vector)**2
)
print("------------------------------------------------------------")
print(
"Exact GS energy:",
(exact_local_gs_vec.conj().T @ local_hamiltonian_matrix @ exact_local_gs_vec)[0, 0].real
)
print("VQE energy:", cost_fn(opt_result.params))
print("------------------------------------------------------------")
print("Exact state <Φ0>:", np.round(exact_local_gs_vec.conj().T @ local_field_operator_matrix @ exact_local_gs_vec, 12)[0, 0].real)
print("VQE state <Φ0>:", estimator(local_field_operator, local_state).value.real)
Overlap with exact local Hamiltonian GS vector: 0.9999129953128841
------------------------------------------------------------
Exact GS energy: 0.5277361259588181
VQE energy: 0.5282862736948237
------------------------------------------------------------
Exact state <Φ0>: -0.0
VQE state <Φ0>: 0.004296736040465901
Now, we build a state for the full system, i.e.
init_state_circuit = QuantumCircuit(system.n_state_qubit)
for site in range(system.n_discretize):
shift = site * system.n_phi_qubit
local_circuit = shift_state_circuit(local_state.circuit, shift)
init_state_circuit.extend(local_circuit)
INIT_STATE = quantum_state(system.n_state_qubit, circuit=init_state_circuit)
We may examine some properties of the prepared state
from quri_parts.qulacs.simulator import evaluate_state_to_vector
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import eigsh
import itertools
import functools
init_state_vector = evaluate_state_to_vector(INIT_STATE).vector
full_hamiltonian = get_hamiltonian_matrix(system)
sparse_full = csc_matrix(full_hamiltonian)
sparse_full.eliminate_zeros()
exact_gs_energy, exact_gs_vec = eigsh(sparse_full, k=1, which="SA")
field_operator = functools.reduce(
lambda a, b: a + b,
[get_scalar_field_operator(n, system.n_phi_qubit, system.mb) for n in range(system.n_discretize)]
)
field_operator_matrix = get_sparse_matrix(field_operator).toarray()
print("Overlap with exact GS:", np.abs(exact_gs_vec.flatten().conj() @ init_state_vector)**2)
print("----------------------------------------------------------------")
print("Exact GS energy:", exact_gs_energy[0])
print("Prepared state <H>:", (init_state_vector.conj() @ full_hamiltonian @ init_state_vector).real)
print("----------------------------------------------------------------")
print("Exact <Φ>:", (exact_gs_vec.flatten().conj() @ field_operator_matrix @ exact_gs_vec).real[0])
print("Prepared <Φ>:", (init_state_vector.conj() @ field_operator_matrix @ init_state_vector).real)
Overlap with exact GS: 0.9720252441279884
----------------------------------------------------------------
Exact GS energy: 1.4021160736167078
Prepared state <H>: 1.5068930934213345
----------------------------------------------------------------
Exact <Φ>: 1.4051260155412137e-16
Prepared <Φ>: 0.008593472080931645
There's a nonzero for the prepared state, so there should be non-trivial time evolution. The ground state with a non-trivial VEV needs to be prepared with adiabatic time evolution, which is not covered in this notebook, but is introduced in section IV of [1].
Time Evolution
Here, we evolve the prepared state from the last section.
Trotter time evolution
It is generally hard to build the exact time evolution due to the non-trivial commutator . Usually, Trotterization is used to approximate the exact time evolution on a quantum computer:
In QURI Algo, circuit implementation of time evolution operators are represented by subclasses of TimeEvolutionCircuitFactory
, where we build one for Trotterization implementation of the time evolution operator.
import math
from quri_algo.circuit.time_evolution.interface import TimeEvolutionCircuitFactory
from quri_algo.circuit.utils.transpile import apply_transpiler
from quri_parts.circuit import NonParametricQuantumCircuit, LinearMappedUnboundParametricQuantumCircuit
from quri_parts.core.circuit import add_parametric_commuting_paulis_exp_gate
@dataclass
class DiscreteScalarFieldTrotterTimeEvoFactory(TimeEvolutionCircuitFactory[DiscreteScalarField1D]):
n_trotter: int
def __post_init__(self) -> None:
super().__post_init__()
self._evo_circuit = LinearMappedUnboundParametricQuantumCircuit(
self.encoded_problem.n_state_qubit
)
self._t = self._evo_circuit.add_parameter("t")
self._construct_evolution_circuit()
def field_operator(self, site_label: int) -> Operator:
"""The field operator on a given site.
"""
return get_scalar_field_operator(
site_label, self.encoded_problem.n_phi_qubit, self.encoded_problem.mb
)
def _construct_evolution_circuit(self) -> None:
for _ in range(self.n_trotter):
self._add_conj_momentum_square_evolution()
self._add_ext_field_evolution()
self._add_kinetic_evolution()
self._add_phi_square_evolution()
self._add_phi4_evolution()
def _add_ext_field_evolution(self) -> NonParametricQuantumCircuit:
"""Add On site Φ_j^2 evolution term
"""
J = self.encoded_problem.external_field
n_site = self.encoded_problem.n_discretize
for site in range(n_site):
field_operator = self.field_operator(site)
operator = J * field_operator
add_parametric_commuting_paulis_exp_gate(
self._evo_circuit, {self._t: -1.0/self.n_trotter}, operator,
)
def _add_phi_square_evolution(self) -> NonParametricQuantumCircuit:
"""Add On site Φ_j^2 evolution term
"""
mb = self.encoded_problem.mb
n_site = self.encoded_problem.n_discretize
for site in range(n_site):
field_operator = self.field_operator(site)
operator = 0.5 * mb**2 * field_operator * field_operator
add_parametric_commuting_paulis_exp_gate(
self._evo_circuit, {self._t: -1.0/self.n_trotter}, operator,
)
def _add_conj_momentum_square_evolution(self) -> None:
"""On site Π_j^2 evolution term.
"""
mb = self.encoded_problem.mb
f_gate = QuantumCircuit(
self.encoded_problem.n_state_qubit, gates=create_F_gate(n_qubits).gates
)
n_site = self.encoded_problem.n_discretize
n_site_qubit = self.encoded_problem.n_phi_qubit
for site in range(n_site):
site_f = shift_state_circuit(f_gate, shift=site * n_site_qubit)
inverse_f = inverse_circuit(site_f)
field_operator = self.field_operator(site)
operator = 0.5 * mb**2 * field_operator * field_operator
self._evo_circuit.extend(inverse_f.gates)
add_parametric_commuting_paulis_exp_gate(
self._evo_circuit, {self._t: -1.0/self.n_trotter}, operator,
)
self._evo_circuit.extend(site_f.gates)
def _add_phi4_evolution(self) -> None:
"""Add Φ_j^4 evolution term to circuit.
"""
lam = self.encoded_problem.lam
n_site = self.encoded_problem.n_discretize
for site in range(n_site):
field_operator = self.field_operator(site)
operator = (lam / math.factorial(4)) * field_operator * field_operator * field_operator * field_operator
add_parametric_commuting_paulis_exp_gate(
self._evo_circuit, {self._t: -1.0/self.n_trotter}, operator,
)
def _add_kinetic_evolution(self) -> None:
"""Add (Φ_j - Φ_{j+1})^2 evolution term to circuit.
"""
n_site = self.encoded_problem.n_discretize
for site in range(n_site - 1):
field_operator = self.field_operator(site)
next_field_operator = self.field_operator(site + 1)
diff = next_field_operator - field_operator
operator = 0.5 * diff * diff
add_parametric_commuting_paulis_exp_gate(
self._evo_circuit, {self._t: -1.0/self.n_trotter}, operator,
)
@apply_transpiler
def __call__(self, evolution_time: float) -> NonParametricQuantumCircuit:
return self._evo_circuit.bind_parameters([evolution_time])
Exact time evolution
We can also build one exact time evolution operator using one unitary matrix gate
import functools
import numpy.typing as npt
from scipy.linalg import expm
@dataclass
class DiscreteScalarFieldExactTimeEvoFactory(TimeEvolutionCircuitFactory[DiscreteScalarField1D]):
def __post_init__(self) -> None:
super().__post_init__()
@functools.cached_property
def hamiltonian_matrix(self) -> npt.NDArray[np.complex128]:
return get_hamiltonian_matrix(self.encoded_problem)
def get_time_evolution_matrix(self, evolution_time: float) -> npt.NDArray[np.complex128]:
return expm(-1j * evolution_time * self.hamiltonian_matrix)
@apply_transpiler
def __call__(self, evolution_time: float) -> NonParametricQuantumCircuit:
n_qubits = self.encoded_problem.n_state_qubit
evo = self.get_time_evolution_matrix(evolution_time)
circuit = QuantumCircuit(self.encoded_problem.n_state_qubit)
circuit.add_UnitaryMatrix_gate(list(range(n_qubits)), evo)
return circuit