Simulasikan dinamika Bose-Hubbard dengan pilihan encoding

Panduan ini membangun circuit untuk dinamika waktu model Bose-Hubbard satu dimensi, mulai dari menuliskan Hamiltonian sampai memeriksa densitas site hasil evolusi terhadap referensi eksak. Yang membuat jalur ini menarik adalah satu ide:

Jangan memilih encoding boson-ke-qubit saat membangun circuit. Bangun circuit di atas mode bosonik, dan biarkan transpiler yang memutuskan bagaimana mode-mode itu menjadi qubit.

Untuk boson, hampir selalu kamu perlu menunda encoding. Sebuah mode bosonik bisa menampung partikel dalam jumlah berapa pun, jadi ruang Fock-nya berdimensi tak hingga dan tidak ada encoding sampai kamu memilih dulu sebuah truncation $d$ (jumlah level okupansi yang kamu pertahankan). Dalam kode, truncation diberikan sebagai okupansi terbesar yang bisa ditampung sebuah mode, max_occupation $= d - 1$, jadi mode yang mempertahankan level $0, 1, 2$ punya max_occupation = 2 dan $d = 3$.

Itu membuat truncation menjadi parameter encoding, bukan parameter fisika, dan parameter ini memang ingin kamu variasikan: terlalu kecil jawabannya salah, terlalu besar circuit-nya melebar tanpa perlu. Kalau kamu menunda encoding sampai waktu transpile, kamu bisa menyapu truncation tanpa mengubah circuit, Hamiltonian, atau pengelompokannya.

Lihat juga

Panduan operator untuk struktur data operator, panduan mapper untuk encoding dan resep menulis mapper kustom, serta panduan transpilasi untuk tahap-tahap transpiler yang dirujuk di seluruh panduan ini.

Truncation, dan biayanya

Ketika kamu memotong sebuah mode menjadi $d$ level, kamu merusak komutator kanonik pada level tertinggi yang dipertahankan [1], sehingga tidak ada encoding yang merupakan isomorfisme aljabar: encoding bosonik memetakan state, dan operator yang bekerja di dalam subruang terpotong. Dua konsekuensi penting untuk semua yang berikut:

1. $d$ harus membatasi okupansi yang bisa dicapai oleh dinamika, bukan hanya okupansi state awal.

2. Pemetaan tidak komut dengan normal ordering, jadi putuskan sekali apakah Hamiltonian-mu didefinisikan dalam bentuk normal-ordered, lalu tetap konsisten.

Panduan mapper menurunkan cacat tersebut dan menunjukkan kedua konsekuensi pada kode yang berjalan.

1. Hamiltonian

Model Bose-Hubbard pada lattice dengan bond $\langle i j \rangle$ adalah

$$H = -t \sum_{\langle i j \rangle} \left(b^\dagger_i b_j + b^\dagger_j b_i\right) + \frac{U}{2} \sum_i n_i (n_i - 1) \, ,$$

dengan $n_i = b^\dagger_i b_i$. Interaksi ini menghitung pasangan boson pada sebuah site, jadi sebuah site harus bisa menampung dua partikel sekaligus. Dalam bentuk normal-ordered, $n_i(n_i-1) = b^\dagger_i b^\dagger_i b_i b_i$, yang merupakan bentuk yang dipakai selanjutnya.

Geometrinya hanya masuk lewat daftar edge, jadi tidak ada yang spesifik untuk sebuah rantai:

>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>>
>>> def bose_hubbard(edges, num_sites, tunneling, interaction):
...     """Build the normal-ordered Bose-Hubbard Hamiltonian on an arbitrary lattice."""
...     terms = {}
...     for i, j in edges:
...         terms[(cre(i), ann(j))] = -tunneling
...         terms[(cre(j), ann(i))] = -tunneling
...     for i in range(num_sites):
...         # U/2 * n_i (n_i - 1), already normal-ordered
...         terms[(cre(i), cre(i), ann(i), ann(i))] = 0.5 * interaction
...     return BosonOperator.from_dict(terms)
>>>
>>> num_sites = 4
>>> edges = [(i, i + 1) for i in range(num_sites - 1)]
>>> hamiltonian = bose_hubbard(edges, num_sites, tunneling=1.0, interaction=2.0)
>>> print(f"Number of terms: {len(hamiltonian)}")
Number of terms: 10

Ada dua properti yang layak dicek dulu. Hamiltonian harus Hermitian agar $e^{-itH}$ unitary, dan ia menjaga jumlah total boson tetap:

>>> hamiltonian.is_hermitian()
True
>>> hamiltonian.conserves_particle_number()
True

Sekarang pertanyaan bosoniknya: berapa banyak level yang dituntut operator ini?

>>> hamiltonian.max_occupation()
2

Suku interaksi menerapkan dua operator pembentukan pada satu mode, jadi sebuah mode harus bisa menampung dua boson agar suku itu bisa bekerja sama sekali. Itu berarti okupansi maksimum sebesar 2, yaitu angka yang diterima LevelSpec sebagai max_occupation-nya. Sebuah mode yang mempertahankan okupansi $0, 1, 2$ punya $d = \texttt{max\_occupation} + 1 = 3$ level, jadi max_occupation >= 2 adalah syarat yang sama dengan $d \geq 3$.

Ini hanya batas bawah: ia berasal dari satu kali penerapan operator, dan penerapan berulang selama evolusi waktu bisa menjangkau lebih jauh. Batas ini menyingkirkan cutoff yang pasti terlalu kecil, dan langkah 6 menentukan, untuk model dan state ini, seberapa besar cutoff yang diperlukan. mode_occupations() melaporkan batas untuk tiap mode secara terpisah, yang menjadi pembanding untuk cutoff per mode.

2. Kelompokkan sebelum memetakan

Evolusinya di-Trotterize, dan cara Hamiltonian dipartisi menentukan apa yang diserahkan ke sintesis. Dua pertimbangan mengarah ke arah yang sama.

Pertama, pertimbangan fisika: suku hopping individual $b^\dagger_i b_j$ tidak Hermitian. Hanya pasangan $b^\dagger_i b_j + b^\dagger_j b_i$ yang Hermitian. Kalau kamu mengevolusikan dengan satu suku non-Hermitian, operasinya tidak unitary, dan operator terpetakan akan membawa koefisien kompleks yang ditolak oleh PauliEvolutionGate. Karena itu, dua suku hopping dari tiap bond harus tetap bersama.

Kedua, pertimbangan efisiensi: semua suku interaksi on-site bersifat diagonal pada basis okupansi dan saling komut, jadi mereka bisa berbagi satu grup dan disintesis bersama. Hal yang sama berlaku untuk bond yang tidak menyentuh site yang sama. Pada rantai, setiap bond kedua membentuk himpunan seperti itu: bond 0 dan 2 bekerja pada pasangan mode yang saling lepas, jadi suku hopping mereka komut dan bisa berbagi satu grup, sedangkan bond 1 membentuk grupnya sendiri. Kalau kamu mengelompokkan bond genap dan ganjil, kamu hanya butuh tiga grup, sementara satu grup per bond butuh empat. Grup yang lebih sedikit berarti faktor yang lebih sedikit, dan juga circuit yang lebih dangkal: dengan Gray code pada max_occupation = 2, depth dua-qubit-nya 95 berbanding 139. Lattice menentukan bond mana yang bisa berbagi grup, itulah sebabnya partisi dipilih di sini.

Atribut groups membawa keputusan itu ke dalam circuit: gate Evolution di atas operator yang dikelompokkan terurai menjadi satu sub-evolusi per grup. Partisi dipilih di level bosonik, selagi struktur lattice masih terlihat, dan ia bertahan utuh sepanjang pemetaan. Panduan pengelompokan membahas pemeriksaan yang harus lolos oleh sebuah pengelompokan.

>>> def group_of(actions, edges):
...     """Assign each term to the even bonds, the odd bonds, or the on-site group."""
...     creations = sorted(mode for action, mode in actions if action)
...     annihilations = sorted(mode for action, mode in actions if not action)
...     if creations == annihilations:
...         return 2  # diagonal: the on-site interaction group
...     bond = tuple(sorted(set(creations) | set(annihilations)))
...     return edges.index(bond) % 2  # bonds 0 and 2 share a group, bond 1 has its own
>>>
>>> hamiltonian.groups = [group_of(actions, edges) for actions, _ in hamiltonian.iter_terms()]
>>> hamiltonian.num_groups()
3

Setiap grup harus Hermitian juga, dan itulah yang didapat dari pemasangan sebelumnya. groups_are_hermitian() memeriksa semua grup dalam satu panggilan:

>>> from qiskit_bosons.operators.terms import groups_are_hermitian
>>> groups_are_hermitian(hamiltonian)
[True, True, True]

Untuk melihat isi tiap grup, pisahkan grup-grupnya:

>>> for index, group in enumerate(hamiltonian.split_out_groups()):
...     print(f"group {index}: {len(group)} terms, hermitian={group.is_hermitian()}")
group 0: 4 terms, hermitian=True
group 1: 2 terms, hermitian=True
group 2: 4 terms, hermitian=True

Pengelompokan genap-dan-ganjil adalah sebuah pilihan, dan alternatifnya mudah diukur. Kalau kamu memberi setiap bond grupnya sendiri, kamu butuh empat grup, dan circuit yang sama kemudian memakan depth dua-qubit lebih banyak. Depth dihitung setelah circuit diturunkan dengan Gray code pada max_occupation = 2, yang dijelaskan di langkah 4:

>>> from qiskit.transpiler import PassManager
>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import Evolution
>>> from qiskit_bosons.mappers import LevelSpec
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> def two_qubit_depth(operator):
...     """The two-qubit depth of one evolution under ``operator``, lowered with the Gray code."""
...     probe = BosonicCircuit(num_sites)
...     probe.append(Evolution(num_sites, operator, time=0.4), probe.modes)
...     pass_manager = generate_preset_occupation_pass_manager(LevelSpec(2, "gray"))
...     pass_manager.qubit = PassManager()
...     lowered = pass_manager.run(probe.decompose())
...     return lowered.decompose(reps=6).depth(lambda instruction: len(instruction.qubits) == 2)
>>>
>>> per_bond = bose_hubbard(edges, num_sites, tunneling=1.0, interaction=2.0)
>>> per_bond.groups = [
...     len(edges) if sorted(m for a, m in actions if a) == sorted(m for a, m in actions if not a)
...     else edges.index(tuple(sorted({m for _, m in actions})))
...     for actions, _ in per_bond.iter_terms()
... ]
>>> per_bond.num_groups(), two_qubit_depth(per_bond)
(4, 139)
>>> hamiltonian.num_groups(), two_qubit_depth(hamiltonian)
(3, 95)

3. Circuit bosonik

Perhatikan apa yang tidak ada di circuit ini. Tidak ada truncation, tidak ada encoding, dan tidak ada jumlah qubit. Circuit ini bekerja pada empat mode bosonik:

>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import Evolution, InitializeModes
>>>
>>> total_time = 0.4
>>> occupation = [1, 1, 0, 0]  # one boson on each of the first two sites
>>>
>>> circuit = BosonicCircuit(num_sites)
>>> circuit.append(InitializeModes(occupation), circuit.modes)
>>> circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
>>> dict(circuit.count_ops())
{'InitializeModes': 1, 'Evolution': 1}

Kabel-kabelnya adalah mode bosonik, bukan qubit. Perhatikan bahwa InitializeModes menerima bilangan bulat, bukan flag: sebuah mode bisa menampung boson dalam jumlah berapa pun, jadi state Fock adalah daftar bilangan okupansi.

Satu kali decompose() mengubah pengelompokan dari langkah 2 menjadi satu gate per grup:

>>> decomposed = circuit.decompose()
>>> dict(decomposed.count_ops())
{'Evolution': 3, 'InitializeModes': 1}

Grup pertama memuat bond pada mode 0 dan 1 serta pada mode 2 dan 3, jadi gate-nya mencakup keempat mode, begitu juga grup on-site. Bond tengah punya gate sendiri pada mode 1 dan 2.

Lewati dekomposisi itu dan pengelompokan akan diabaikan begitu saja. Metode yang melakukan pemecahan adalah synthesis dari gate, dengan BosonicLieTrotter sebagai default-nya; Sintesis gate Evolution di ruang bosonik menunjukkan alternatifnya.

4. Pilih encoding saat waktu transpile

Inilah hasilnya. Circuit yang sama diturunkan empat kali, dengan empat nilai LevelEncoding berbeda, hanya dengan mengubah satu argumen. Sebuah LevelSpec memasangkan encoding dengan truncation setiap mode, dan preset berikut menurunkan layout, persiapan state, dan pemetaan operator dari satu spec itu (lihat Menjelaskan truncation dengan LevelSpec). String seperti "gray" diterima di mana pun sebuah LevelEncoding diterima.

>>> from qiskit.transpiler import PassManager
>>> from qiskit_bosons.mappers import LevelEncoding, LevelSpec
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> max_occupation = 2  # the largest occupation per mode, so d = max_occupation + 1 = 3 levels
>>>
>>> encodings = (
...     LevelEncoding.UNARY_PROJECTED,
...     LevelEncoding.UNARY_UNPROJECTED,
...     LevelEncoding.STANDARD_BINARY,
...     LevelEncoding.GRAY,
... )
>>> for encoding in encodings:
...     pass_manager = generate_preset_occupation_pass_manager(
...         LevelSpec(max_occupation, encoding)
...     )
...     pass_manager.qubit = PassManager()  # keep the circuit as the pipeline synthesized it
...     qubit_circuit = pass_manager.run(decomposed)
...     basis = qubit_circuit.decompose(reps=6)
...     two_qubit = basis.depth(lambda instruction: len(instruction.qubits) == 2)
...     print(f"{encoding!s:18s} qubits={qubit_circuit.num_qubits:2d}  2q-depth={two_qubit:4d}")
unary_projected    qubits=12  2q-depth= 645
unary_unprojected  qubits=12  2q-depth= 383
standard_binary    qubits= 8  2q-depth= 218
gray               qubits= 8  2q-depth=  95

Angka-angka ini bertentangan dengan tebakan yang wajar. Encoding unary memakai $d$ qubit per mode, dibandingkan $\lceil \log_2 d \rceil$ untuk kode kompak, dan pada $d = 3$ kode kompak menang di kedua sisi sekaligus: qubit lebih sedikit dan gate entangling jauh lebih sedikit. Kedua varian unary berbeda pada operator yang mereka hasilkan, bukan pada layout-nya: LevelEncoding.UNARY_PROJECTED membawa proyektor pada setiap qubit selang-seling dalam sebuah mode, yang membuatnya paling mahal untuk dievolusikan di antara keempatnya, dan LevelEncoding.UNARY_UNPROJECTED menghilangkan proyektor dan sepakat dengannya pada subruang one-hot. Panduan mapper menampilkan kedua operator berdampingan, dan menjelaskan kenapa Gray code mengalahkan standard binary (SB) pada suku hopping.

Sawaya et al. melaporkan urutan yang sama [2]:

meskipun kode unary cenderung membutuhkan lebih banyak qubit tetapi lebih sedikit operasi, sering kali kode SB atau Gray lebih efisien, baik dari sisi jumlah qubit maupun jumlah operasi

Kesimpulan mereka adalah bahwa kode terbaik bergantung pada Hamiltonian dan pada $d$, dan itu argumen terkuat untuk tidak mengunci pilihan di dalam circuit.

Catatan

Di bawah level encoding, sebuah mode mencakup beberapa qubit, jadi register qubit lebih lebar daripada register mode. Preset memakai OccupationB2QLayout, yang menurunkan lebar register dari LevelSpec. Panduan transpilasi menelusuri tahap-tahap yang dijalankan preset.

5. Verifikasi terhadap dinamika eksak

Tidak ada langkah sebelumnya yang membuktikan bahwa circuit-nya benar. Untuk itu, bandingkan densitas site $\langle n_m \rangle$ setelah evolusi dengan referensi eksak yang dihitung pada Hamiltonian terenkode.

Referensi harus dihitung di ruang terpotong yang sama dengan tempat circuit bekerja, karena membandingkan dengan perhitungan tanpa truncation akan mencampur error Trotter dengan error truncation, yang merupakan dua hal berbeda dan layak dipisahkan (lihat Memotong ruang Fock).

>>> import numpy as np
>>> from qiskit.quantum_info import SparsePauliOp, Statevector
>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>>
>>> def sparse_matrix(operator, encoding):
...     """Map a bosonic operator and hand back a sparse qubit matrix."""
...     spec = LevelSpec(max_occupation, encoding).resolve(num_sites)
...     observable = occupation_encoding(operator, spec)
...     return SparsePauliOp.from_sparse_observable(observable).to_matrix(sparse=True)
>>>
>>> def fock_index(occupation, encoding):
...     """The computational basis index encoding a product Fock state."""
...     spec = LevelSpec(max_occupation, encoding).resolve(num_sites)
...     return sum(
...         encoding.level_to_bits(n) << offset for n, offset in zip(occupation, spec.offsets)
...     )
>>>
>>> def densities(state, encoding):
...     """The per-site boson number <n_m> of an encoded state vector."""
...     out = []
...     for mode in range(num_sites):
...         number = BosonOperator.from_dict({(cre(mode), ann(mode)): 1.0})
...         matrix = sparse_matrix(number, encoding)
...         out.append(float((state.conj() @ (matrix @ state)).real))
...     return out

Dengan semua itu di tangan, berikut evolusi eksak di ruang terpotong:

>>> encoding = LevelEncoding.GRAY
>>> num_qubits = LevelSpec(max_occupation, encoding).resolve(num_sites).total_width
>>>
>>> initial = np.zeros(2**num_qubits, dtype=complex)
>>> initial[fock_index(occupation, encoding)] = 1.0
>>>
>>> evolved_exact = expm_multiply(
...     -1j * total_time * sparse_matrix(hamiltonian, encoding), initial
... )
>>> reference = densities(evolved_exact, encoding)
>>> for site, value in enumerate(reference):
...     print(f"site {site}: {value:.6f}")
site 0: 0.993834
site 1: 0.858853
site 2: 0.141309
site 3: 0.006004

Boson menyebar dari ujung rantai yang terisi menuju ujung yang kosong, dan kedua boson tetap terjaga:

>>> bool(np.isclose(sum(reference), 2.0))
True

Sekarang circuit-nya. Perhatikan bahwa referensi sebelumnya dibangun langsung dari initial, jadi circuit yang dibandingkan dengannya harus hanya berisi evolusi, karena gate InitializeModes akan menyiapkan ulang state yang sudah disiapkan.

>>> from qiskit_bosons.circuit.library.synthesis import BosonicLieTrotter
>>>
>>> def trotter_densities(steps, encoding):
...     """Site densities after `steps` first-order Trotter steps of the grouped evolution."""
...     circuit = BosonicCircuit(num_sites)
...     circuit.append(
...         Evolution(
...             num_sites, hamiltonian, time=total_time, synthesis=BosonicLieTrotter(reps=steps)
...         ),
...         circuit.modes,
...     )
...     pass_manager = generate_preset_occupation_pass_manager(
...         LevelSpec(max_occupation, encoding)
...     )
...     qubit_circuit = pass_manager.run(circuit.decompose())
...     # decompose BEFORE evolving: see the warning that follows
...     basis = qubit_circuit.decompose(reps=6)
...     state = np.asarray(Statevector(initial).evolve(basis).data)
...     return densities(state, encoding)
>>>
>>> for steps in (1, 2, 4, 8, 16):
...     error = max(
...         abs(a - b) for a, b in zip(trotter_densities(steps, encoding), reference)
...     )
...     print(f"steps={steps:2d}  max density error = {error:.4f}")
steps= 1  max density error = 0.0843
steps= 2  max density error = 0.0389
steps= 4  max density error = 0.0183
steps= 8  max density error = 0.0088
steps=16  max density error = 0.0043

Error-nya turun setengah ketika jumlah langkah digandakan, yang merupakan ciri khas $O(t^2/\text{steps})$ dari product formula orde pertama. Itu bukti bahwa selisihnya adalah error Trotter dan encoding-nya setia.

Penting

Pengulangan berasal dari synthesis gate, yang melakukan Trotterize di ruang bosonik tempat pengelompokan masih ada: setiap langkah menerapkan setiap grup sekali, dan langkah-langkahnya berurutan. Itu berbeda dengan memberikan LieTrotter dengan reps > 1 ke plugin sintesis tingkat qubit. Di sana, setiap grup sudah menjadi PauliEvolutionGate sendiri, sehingga formula hanya membagi suku-suku di dalam satu grup dan tidak pernah menyisipkan grup-grup secara bergantian. Grup-grup itu tidak saling komut, jadi yang tersisa adalah satu langkah Trotter luar yang error-nya tidak mengecil sebesar apa pun reps. Sintesis gate Evolution di ruang bosonik membandingkan keduanya dan membahas orde yang lebih tinggi.

Peringatan

Pemanggilan decompose() sebelumnya itu penting. evolve() memakai definisi eksak dari PauliEvolutionGate dan melewati product formula sepenuhnya. Kalau kamu memverifikasi terhadap circuit yang belum didekomposisi, kamu akan melihat kecocokan nyaris sempurna bahkan untuk sintesis yang sengaja dibuat salah. Circuit yang belum didekomposisi juga memaksa pembangunan operator padat $2^{n} \times 2^{n}$, sedangkan circuit yang didekomposisi menerapkan gate satu per satu. Simulator lain bisa memperlakukan gate secara berbeda, jadi pemeriksaan ini bergantung pada simulator, sedangkan hardware selalu butuh bentuk yang didekomposisi.

Terakhir, pemeriksaan utama: fisikanya tidak boleh bergantung pada encoding yang dipilih.

>>> def exact_densities(encoding):
...     """Exact truncated-space densities, computed under a given encoding."""
...     width = LevelSpec(max_occupation, encoding).resolve(num_sites).total_width
...     start = np.zeros(2**width, dtype=complex)
...     start[fock_index(occupation, encoding)] = 1.0
...     evolved = expm_multiply(
...         -1j * total_time * sparse_matrix(hamiltonian, encoding), start
...     )
...     return densities(evolved, encoding)
>>>
>>> all(
...     bool(np.allclose(exact_densities(name), reference, atol=1e-10))
...     for name in encodings
... )
True

Jumlah qubit, struktur Pauli, dan depth circuit semuanya berbeda, dan fisikanya cocok sampai sepuluh angka desimal. Encoding adalah pilihan implementasi, itulah kenapa ia menjadi urusan transpiler. Kedua varian unary cocok di sini karena state yang disimulasikan tetap berada di subruang one-hot; di luar itu keduanya akan berbeda.

6. Truncation adalah parameter konvergensi

Langkah terakhir ini spesifik untuk boson. Sejauh ini semuanya memakai max_occupation = 2 yang tetap, yaitu $d = 3$. Namun, truncation bukan pilihan bebas: ia adalah parameter konvergensi, dan sebuah jawaban baru bisa dipercaya setelah ia berhenti bergerak (lihat Memotong ruang Fock).

Circuit tidak membawa truncation, jadi circuit itu sendiri tidak berubah ketika cutoff berubah: hanya transpilasinya yang diulang. Yang berubah adalah ukuran hasilnya, dan sapuan seperti ini hanya semurah di sini karena sistemnya sangat kecil. Referensi berikut adalah simulasi klasik eksak dari dinamika terpotong, yang mungkin dilakukan untuk empat site dan tumbuh secara eksponensial terhadap jumlah qubit. Untuk masalah di luar jangkauan itu, tidak ada referensi eksak untuk dibandingkan, dan konvergensinya harus dinilai dengan cara lain.

>>> def densities_at_cutoff(max_occ, encoding=LevelEncoding.GRAY):
...     """Exact truncated-space densities at a given cutoff."""
...     spec = LevelSpec(max_occ, encoding).resolve(num_sites)
...
...     start = np.zeros(2**spec.total_width, dtype=complex)
...     start[
...         sum(encoding.level_to_bits(n) << offset for n, offset in zip(occupation, spec.offsets))
...     ] = 1.0
...
...     def matrix(operator):
...         observable = occupation_encoding(operator, spec)
...         return SparsePauliOp.from_sparse_observable(observable).to_matrix(sparse=True)
...
...     evolved = expm_multiply(-1j * total_time * matrix(hamiltonian), start)
...     out = []
...     for mode in range(num_sites):
...         number = BosonOperator.from_dict({(cre(mode), ann(mode)): 1.0})
...         out.append(float((evolved.conj() @ (matrix(number) @ evolved)).real))
...     return out
>>>
>>> previous = None
>>> for max_occ in (1, 2, 3, 4):
...     current = densities_at_cutoff(max_occ)
...     # Report the drift bucketed rather than to full precision: after it reaches machine
...     # epsilon, the exact digits are arithmetic noise and vary between BLAS builds.
...     if previous is None:
...         drift = "--"
...     else:
...         gap = max(abs(a - b) for a, b in zip(current, previous))
...         drift = f"{gap:.0e}" if gap > 1e-12 else "< 1e-12"
...     print(f"d={max_occ + 1}  "
...           + " ".join(f"{x:.6f}" for x in current)
...           + f"   drift={drift}")
...     previous = current
d=2  0.993986 0.858173 0.141827 0.006014   drift=--
d=3  0.993834 0.858853 0.141309 0.006004   drift=7e-04
d=4  0.993834 0.858853 0.141309 0.006004   drift=< 1e-12
d=5  0.993834 0.858853 0.141309 0.006004   drift=< 1e-12

Pada max_occupation = 1 ($d = 2$) setiap site menampung paling banyak satu boson, sehingga suku interaksi, yang butuh dua boson pada satu site, sama sekali tidak bisa bekerja, dan densitasnya meleset sekitar $7 \times 10^{-4}$. Ketika kamu menambahkan level ketiga, hasilnya menangkap interaksi, dan mulai dari max_occupation = 3 hasilnya stasioner sampai presisi mesin: sudah konvergen.

Alasan ia berhenti di max_occupation = 2 adalah jumlah partikel. Hamiltonian menjaga jumlah boson, dan state ini memuat dua, jadi tidak ada site yang pernah bisa menampung lebih dari dua, berapa pun waktu evolusinya dan baik interaksinya tolak-menolak maupun tarik-menarik. Kalau kamu mengevolusikan sampai $t = 6$ alih-alih $t = 0.4$, atau membalik tanda $U$, densitas pada max_occupation = 2 dan lebih tinggi tetap tidak berubah sampai $10^{-13}$.

Yang justru menaikkan kebutuhan adalah jumlah boson. Tiga boson pada rantai yang sama berubah sekitar $2 \times 10^{-3}$ antara max_occupation = 2 dan 3 setelah evolusi pendek yang sama, dan sekitar $0.3$ setelah $t = 6$. Itu pernyataan tentang model ini dan satu besaran terjaga, dan ia hanya batas untuk operator yang menjaga jumlah partikel; itu bukan aturan umum, jadi cutoff harus disapu, bukan ditebak.

Sapuan ini punya harga di sisi circuit juga, yang tidak ditunjukkan oleh referensi eksak. Circuit bosonik yang sama, diturunkan dengan Gray code pada tiap cutoff:

>>> for max_occ in (1, 2, 3, 4):
...     pass_manager = generate_preset_occupation_pass_manager(
...         LevelSpec(max_occ, LevelEncoding.GRAY)
...     )
...     pass_manager.qubit = PassManager()
...     lowered = pass_manager.run(decomposed)
...     depth = lowered.decompose(reps=6).depth(lambda instruction: len(instruction.qubits) == 2)
...     print(f"d={max_occ + 1}  qubits={lowered.num_qubits:2d}  2q-depth={depth:4d}")
d=2  qubits= 4  2q-depth=   8
d=3  qubits= 8  2q-depth=  95
d=4  qubits= 8  2q-depth= 217
d=5  qubits=12  2q-depth=1176

Setiap level tambahan menambah depth, dan jumlah qubit melonjak setiap kali $d$ melewati pangkat dua. Ketika kamu mengonvergenkan $d$, kamu jadi menukar akurasi dengan ukuran circuit, jadi sapuan ini lebih dari sekadar pemeriksaan numerik.

Perhatikan juga bahwa $d = 2$ yang hampir benar hanyalah kebetulan dari observabelnya: densitas site nyaris tidak bergerak, sementara state yang mendasarinya salah secara kualitatif (tidak punya komponen yang terisi ganda). Kalau satu observabel konvergen, itu tidak menunjukkan bahwa observabel lain sudah konvergen.

Yang perlu dibawa pulang

• Encoding adalah urusan transpiler. Circuit bosonik di langkah 3 tidak menyebut truncation, encoding, ataupun jumlah qubit. Empat encoding lalu menurunkannya menjadi empat circuit berbeda, dan semuanya menghasilkan fisika yang sama. Tidak ada yang perlu diubah di bagian hulu transpiler.

• Untuk boson, ini bukan sekadar kenyamanan. Truncation adalah parameter konvergensi yang harus disapu (sweep), dan letaknya di dalam encoding. Circuit yang dibangun terhadap encoding tetap harus dibangun ulang untuk setiap nilai $d$.

• Kelompokkan sebelum mapping, dan kelompokkan demi Hermitisitas. Satu suku hopping tidak Hermitian dan tidak bisa dieksponensialkan; pasangan bond di sisi lain bersifat Hermitian sehingga bisa dieksponensialkan. Tentukan partisinya di level bosonik, saat kisi (lattice) masih terlihat. Itu perlu secara fisika, dan atribut groups ada untuk membawa keputusan itu ke tahap berikutnya.

• Jangan berasumsi encoding berbobot rendah itu yang murah. Pada $d = 3$ kode kompak mengalahkan kedua varian unary sekaligus dalam jumlah qubit dan gate, dan varian terproyeksi adalah yang paling mahal dari semuanya. Ukur sendiri; titik silangnya bergantung pada $d$.

• Truncation merusak aljabar, dengan cara yang bisa dikuantifikasi. Hanya pada level teratas yang dipertahankan, sebesar $d$. Itu juga alasan mapping tidak komut dengan normal ordering.

Langkah berikutnya

Dua arah bisa memperluas contoh ini. Kalau kamu mengganti daftar edge rantai dengan bond kisi persegi, tidak ada perubahan lain yang dibutuhkan, karena fungsi bose_hubbard yang didefinisikan sebelumnya sudah agnostik terhadap geometri. Encoding buatanmu sendiri (misalnya kode block-unary, yang menginterpolasi antara unary dan kompak) juga bisa dipasang lewat seam MapperFnEvolutionSynthesis yang sama dengan yang dipakai preset.

• Prinsip desain representasi operator menjelaskan struktur data operator di balik Hamiltonian.

• Memetakan representasi operator membahas encoding dan cara menulis mapper sendiri.

• Transpile circuit bosonik menelusuri tahap-tahap yang dijalankan preset.

• Mensintesis Evolution gate di ruang bosonik menunjukkan cara mengurangi error Trotter di ruang bosonik.

Penting

Semua depth yang dilaporkan berasal dari register abstrak yang terhubung penuh. Mapping ke konektivitas hardware nyata adalah tahap terpisah, dan blok qubit per mode pada occupation encoding memberi router struktur yang cukup berbeda dibandingkan satu qubit per mode.

[1]

R. Somma, G. Ortiz, E. Knill, and J. Gubernatis, Quantum Simulations of Physics Problems (2003), arXiv:quant-ph/0304063, Sec. 4, Eq. 22.

[2]

N. P. D. Sawaya, T. Menke, T. H. Kyaw, S. Johri, A. Aspuru-Guzik, and G. G. Guerreschi, Resource-efficient digital quantum simulation of d-level systems for photonic, vibrational, and spin-s Hamiltonians, npj Quantum Information 6, 49 (2020), arXiv:1909.12847.