I-simulate ang Bose-Hubbard dynamics na may mapipiling encoding

Ang gabay na ito ay bumubuo ng Circuit para sa time dynamics ng one-dimensional na Bose-Hubbard model, mula sa pagsulat ng Hamiltonian hanggang sa pagsuri ng mga evolved na site density laban sa eksaktong reference. Ang dahilan kung bakit kawili-wili ang ruta ay isang ideya lang:

Huwag pumili ng boson-to-qubit encoding kapag binubuo mo ang Circuit. Buuin ang Circuit sa bosonic modes, at hayaan ang Transpiler na magpasya kung paano magiging mga qubit ang mga mode na iyon.

Para sa mga boson, halos palagi mong kailangang ipagpaliban ang encoding. Ang isang bosonic mode ay kayang maglaman ng kahit ilang particle, kaya ang Fock space nito ay infinite-dimensional at walang encoding hangga't hindi ka muna pumipili ng truncation $d$ (ang bilang ng mga occupation level na itatago mo). Sa code, ibinibigay ang truncation bilang pinakamalaking occupation na kayang hawakan ng isang mode, max_occupation $= d - 1$, kaya ang isang mode na nagtatago ng mga level na $0, 1, 2$ ay may max_occupation = 2 at $d = 3$.

Dahil dito, ang truncation ay parameter ng encoding, hindi ng physics, at ito ay parameter na gusto mong baguhin: kapag masyadong maliit, mali ang sagot, at kapag masyadong malaki, hindi kailangang lumapad ang Circuit. Kung ipagpapaliban mo ang encoding hanggang sa transpile time, maaari mong i-sweep ang truncation nang hindi binabago ang Circuit, ang Hamiltonian, o ang grouping.

Tingnan din

Ang gabay sa operators para sa operator data structure, ang gabay sa mappers para sa mga encoding at ang recipe sa pagsulat ng custom mapper, at ang gabay sa transpilation para sa mga yugto ng Transpiler na tinutukoy sa buong gabay.

Ang truncation, at ang halaga nito

Kapag tinruncate mo ang isang mode sa $d$ na level, nasisira mo ang canonical commutator sa pinakamataas na level na itinago [1], kaya walang encoding na isomorphism ng mga algebra: ang isang bosonic encoding ay nagmamapa ng mga state, at mga operator na kumikilos sa loob ng truncated subspace. Dalawang bunga ang mahalaga sa lahat ng susunod:

1. Dapat na ang $d$ ang maglimita sa occupation na maaabot ng dynamics, hindi lang ng initial state.

2. Hindi nagko-commute ang mapping sa normal ordering, kaya magpasya nang minsanan kung ang Hamiltonian mo ay nakadefine sa normal-ordered na anyo, at manatiling pare-pareho.

Ang gabay sa mappers ang naghahango ng depekto at nagpapakita ng parehong bunga sa tumatakbong code.

1. Ang Hamiltonian

Ang Bose-Hubbard model sa isang lattice na may mga bond na $\langle i j \rangle$ ay

$$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) \, ,$$

na may $n_i = b^\dagger_i b_i$. Binibilang ng interaksyon ang mga pares ng boson sa isang site, kaya kailangang magkaroon ang site ng dalawang particle nang sabay. Kapag normal-ordered, $n_i(n_i-1) = b^\dagger_i b^\dagger_i b_i b_i$, na siyang anyong gagamitin sa susunod.

Sa edge list lang pumapasok ang geometry, kaya walang anumang partikular sa isang chain:

>>> 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

Dalawang katangian ang dapat munang suriin. Dapat Hermitian ang Hamiltonian para maging unitary ang $e^{-itH}$, at iniingatan nito ang kabuuang bilang ng boson:

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

Ngayon, ang tanong tungkol sa boson: ilang level ang hinihingi ng operator na ito?

>>> hamiltonian.max_occupation()
2

Ang interaction term ay naglalapat ng dalawang creation operator sa isang mode, kaya kailangang magkaroon ang mode ng dalawang boson para gumana ang term. Iyon ay maximum occupation na 2, na siyang numerong tinatanggap ng LevelSpec bilang max_occupation nito. Ang isang mode na nagtatago ng mga occupation na $0, 1, 2$ ay may $d = \texttt{max\_occupation} + 1 = 3$ na level, kaya ang max_occupation >= 2 ay kapareho ng requirement na $d \geq 3$.

Lower bound lang ito: nagmumula ito sa isang paglalapat ng operator, at ang paulit-ulit na paglalapat sa buong time evolution ay maaaring umabot nang mas malayo. Inaalis nito ang mga cutoff na tiyak na masyadong maliit, at ang hakbang 6 ang magsasabi, para sa model at state na ito, kung gaano kalaki dapat ang cutoff. Iniuulat ng mode_occupations() ang bound para sa bawat mode nang hiwalay, na siyang ikinukumpara sa cutoff kada mode.

2. Mag-group bago mag-map

Trotterized ang evolution, at ang paraan ng paghahati ng Hamiltonian ang nagpapasya kung ano ang ibibigay sa synthesis. Dalawang konsiderasyon ang nagtutulak sa iisang direksyon.

Una, isang pisikal: ang indibidwal na hopping term na $b^\dagger_i b_j$ ay hindi Hermitian. Ang pares lang na $b^\dagger_i b_j + b^\dagger_j b_i$ ang Hermitian. Kung ie-evolve mo sa ilalim ng iisang non-Hermitian na term, hindi unitary ang operasyon, at ang na-map na operator ay magkakaroon ng mga complex coefficient na tinatanggihan ng PauliEvolutionGate. Kaya dapat manatiling magkasama ang dalawang hopping term ng bawat bond.

Pangalawa, isang pang-episyensiya: ang lahat ng on-site interaction term ay diagonal sa occupation basis at nagko-commute sa isa't isa, kaya maaari silang magbahagi ng iisang group at i-synthesize nang magkakasama. Ganoon din para sa mga bond na hindi humahawak sa iisang site. Sa isang chain, bawat ikalawang bond ay ganoong set: ang mga bond 0 at 2 ay kumikilos sa magkahiwalay na pares ng mode, kaya nagko-commute ang mga hopping term nila at maaaring magbahagi ng group, at ang bond 1 ay bumubuo ng sarili nitong group. Kung igagrupo mo ang even at odd na bond, tatlong group ang kailangan mo, samantalang apat ang kailangan kung isang group kada bond. Mas kaunting group ay mas kaunting factor, at mas mababaw din na Circuit: sa Gray code na may max_occupation = 2, ang two-qubit depth ay 95 laban sa 139. Ang lattice ang nagpapasya kung aling mga bond ang maaaring magbahagi ng group, kaya dito pinipili ang partition.

Dinadala ng attribute na groups ang desisyong iyon sa Circuit: ang isang Evolution gate sa ibabaw ng grouped na operator ay nade-decompose sa isang sub-evolution kada group. Pinipili ang partition sa bosonic na antas, habang nakikita pa ang estruktura ng lattice, at nananatili itong buo sa buong mapping. Saklaw ng gabay sa grouping ang mga pagsusuring dapat pasado sa isang grouping.

>>> 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

Dapat Hermitian din ang bawat group, na siyang naibibigay ng naunang pagpapares. Sinusuri ng groups_are_hermitian() ang lahat ng group sa isang tawag:

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

Para makita kung ano ang laman ng bawat group, paghiwalayin ang mga group:

>>> 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

Pagpili ang even-and-odd na grouping, at madaling sukatin ang alternatibo. Kung bibigyan mo ang bawat bond ng sarili nitong group, apat na group ang kailangan mo, at mas malaking two-qubit depth ang magagastos ng parehong Circuit. Binibilang ang depth pagkatapos ibaba ang Circuit gamit ang Gray code na may max_occupation = 2, na ipinapaliwanag ng hakbang 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. Ang bosonic Circuit

Pansinin kung ano ang wala sa Circuit. Walang truncation, walang encoding, at walang bilang ng qubit. Kumikilos ito sa apat na bosonic mode:

>>> 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}

Ang mga wire ay mga bosonic mode, hindi mga qubit. Tandaan na ang InitializeModes ay tumatanggap ng mga integer, hindi flag: ang isang mode ay maaaring maglaman ng kahit ilang boson, kaya ang Fock state ay listahan ng mga occupation number.

Ginagawa ng iisang decompose() na isang gate kada group ang grouping mula sa hakbang 2:

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

Ang unang group ay naglalaman ng mga bond sa mode 0 at 1 at sa mode 2 at 3, kaya sumasaklaw ang gate nito sa lahat ng apat na mode, at ganoon din ang on-site group. Ang gitnang bond ay may sariling gate sa mode 1 at 2.

Kung lalaktawan ang decomposition na iyon, hindi papansinin ang grouping. Ang synthesis ng gate ang paraang gumagawa ng paghahati, na ang BosonicLieTrotter ang default; ipinapakita ng Synthesize an Evolution gate in bosonic space ang mga alternatibo.

4. Piliin ang encoding sa transpile time

Narito ang bentahe. Ang parehong Circuit ay ibinababa nang apat na beses, sa ilalim ng apat na magkakaibang value ng LevelEncoding, sa pamamagitan ng pagbabago ng isang argument. Ipinapares ng LevelSpec ang encoding sa truncation ng bawat mode, at hinahango ng sumusunod na preset ang layout, ang state preparation, at ang operator mapping mula sa isang spec na iyon (tingnan ang Describe the truncation with a LevelSpec). Tinatanggap ang string gaya ng "gray" saanman tinatanggap ang LevelEncoding.

>>> 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

Sumasalungat ang mga numerong ito sa natural na hula. Gumagamit ang mga unary encoding ng $d$ qubit kada mode, kumpara sa $\lceil \log_2 d \rceil$ para sa mga compact code, at sa $d = 3$ nananalo ang mga compact code sa dalawang aspeto nang sabay: mas kaunting qubit at mas kaunting entangling gate. Magkaiba ang dalawang unary variant sa operator na ibinibigay nila, hindi sa layout nila: ang LevelEncoding.UNARY_PROJECTED ay may projector sa bawat ikalawang qubit ng isang mode, kaya ito ang pinakamahal i-evolve sa apat, at ang LevelEncoding.UNARY_UNPROJECTED ay walang projector at sumasang-ayon dito sa one-hot subspace. Ipinapakita ng gabay sa mappers ang dalawang operator nang magkatabi, at ipinapaliwanag kung bakit mas mahusay ang Gray code kaysa sa standard binary (SB) sa mga hopping term.

Iniulat ni Sawaya et al. ang parehong pagkakasunod-sunod [2]:

bagaman ang unary code ay may tendensiyang mangailangan ng mas maraming qubit pero mas kaunting operasyon, madalas na ang SB o Gray code ay mas episyente pareho sa bilang ng qubit at sa bilang ng operasyon

Ang konklusyon nila ay nakadepende ang pinakamahusay na code sa Hamiltonian at sa $d$, na siyang pinakamalakas na argumento para hindi ayusin ang pagpili sa Circuit.

Tandaan

Sa ilalim ng level encoding, sumasaklaw ang isang mode sa ilang qubit, kaya mas malapad ang qubit register kaysa sa mode register. Gumagamit ang preset ng OccupationB2QLayout, na hinahango ang lapad ng register mula sa LevelSpec. Dinadaanan ng gabay sa transpilation ang mga yugtong pinapatakbo ng preset.

5. I-verify laban sa eksaktong dynamics

Wala sa mga naunang hakbang ang nagpapatunay na tama ang Circuit. Para doon, ihambing ang mga site density na $\langle n_m \rangle$ pagkatapos ng evolution laban sa eksaktong reference na kinalkula sa na-encode na Hamiltonian.

Dapat kalkulahin ang reference sa parehong truncated space na ginagamit ng Circuit, dahil ang paghahambing sa untruncated na kalkulasyon ay maghahalo sa Trotter error at truncation error, na magkaibang bagay at sulit paghiwalayin (tingnan ang Truncate the Fock space).

>>> 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

Taglay ang mga iyon, ang eksaktong evolution sa truncated space:

>>> 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

Kumakalat ang mga boson mula sa okupadong dulo ng chain patungo sa walang laman, at napapanatili ang dalawang boson:

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

Ngayon, ang Circuit. Pansinin na ang naunang reference ay binuo mula sa initial nang direkta, kaya ang Circuit na ikinukumpara dito ay dapat naglalaman lang ng evolution, dahil ang gate na InitializeModes ay muling maghahanda ng state na nakahanda na.

>>> 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

Nagiging kalahati ang error kapag dinoble ang bilang ng step, na siyang pirma ng $O(t^2/\text{steps})$ ng isang first-order product formula. Iyan ang ebidensya na Trotter error ang pagkakaiba at tapat ang encoding.

Mahalaga

Ang mga pag-uulit ay nagmumula sa synthesis ng gate, na nagta-Trotterize sa bosonic space kung saan umiiral pa ang grouping: bawat step ay naglalapat ng bawat group nang isang beses, at sunud-sunod ang mga step. Iba iyon sa pagbibigay ng LieTrotter na may reps > 1 sa qubit-level synthesis plugin. Doon, ang bawat group ay sarili nang PauliEvolutionGate, kaya hinahati lang ng formula ang mga term sa loob ng isang group at hindi kailanman pinagsasalitan ang mga group. Hindi nagko-commute ang mga group sa isa't isa, kaya nag-iiwan ito ng iisang outer Trotter step na ang error ay hindi lumiliit gaano man kalaki ang reps. Ikinukumpara ng Synthesize an Evolution gate in bosonic space ang dalawa at saklaw ang mas mataas na order.

Babala

Mahalaga ang naunang tawag sa decompose(). Ginagamit ng evolve() ang eksaktong definition ng isang PauliEvolutionGate at llumalampas sa product formula. Kung magve-verify ka laban sa Circuit na hindi decomposed, makikita mo ang halos perpektong pagsang-ayon kahit sa sadyang maling synthesis. Pinipilit din ng un-decomposed na Circuit na bumuo ng siksik na $2^{n} \times 2^{n}$ na operator, samantalang ang decomposed na Circuit ay naglalapat ng mga gate nang isa-isa. Maaaring iba ang pagtrato ng ibang simulator sa gate, kaya nakadepende sa simulator ang pagsusuring ito, samantalang laging kailangan ng hardware ang decomposed na anyo.

Panghuli, ang sentral na pagsusuri: hindi dapat nakadepende ang physics sa kung aling encoding ang pinili.

>>> 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

Magkakaiba ang bilang ng qubit, ang mga Pauli structure, at ang mga depth ng Circuit, at nagkakasundo ang physics hanggang sampung decimal place. Ang encoding ay pagpili sa implementasyon, kaya ito ay nabibilang sa Transpiler. Nagkakasundo rito ang dalawang unary variant dahil nananatili ang na-simulate na state sa one-hot subspace; sa labas nito ay magkakaiba sila.

6. Ang truncation ay isang convergence parameter

Ang huling hakbang ay partikular sa mga boson. Sa lahat ng nauna, nakatakda ang max_occupation = 2, ibig sabihin $d = 3$. Gayunpaman, hindi malayang pagpili ang truncation: ito ay isang convergence parameter, at mapagkakatiwalaan lang ang sagot kapag tumigil na itong gumalaw (tingnan ang Truncate the Fock space).

Walang truncation ang Circuit, kaya hindi nagbabago ang Circuit mismo kapag nagbago ang cutoff: ang transpilation lang ang inuulit. Ang nagbabago ay ang laki ng lumalabas, at ang sweep ay mura lang dito dahil napakaliit ng sistema. Ang sumusunod na reference ay isang eksaktong classical simulation ng truncated dynamics, na posible para sa apat na site at lumalaki nang exponential kasabay ng bilang ng qubit. Para sa problemang lampas sa abot nito, walang eksaktong reference na mapaghahambingan, at kailangang husgahan ang convergence sa ibang paraan.

>>> 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

Sa max_occupation = 1 ($d = 2$), isang boson lang ang kayang hawakan ng bawat site, kaya hindi talaga makakilos ang interaction term, na nangangailangan ng dalawa sa isang site, at mali ng humigit-kumulang $7 \times 10^{-4}$ ang mga density. Kapag idinagdag ang ikatlong level, nakukuha ng resulta ang interaksyon, at mula sa max_occupation = 3 pataas ay stationary ang resulta hanggang machine precision: converged na.

Ang dahilan kung bakit tumitigil ito sa max_occupation = 2 ay ang bilang ng particle. Iniingatan ng Hamiltonian ang bilang ng boson, at dalawa ang laman ng state na ito, kaya walang site na maaaring humawak ng higit sa dalawa, anuman ang evolution time at maging repulsive man o attractive ang interaksyon. Kung ie-evolve mo hanggang $t = 6$ sa halip na $t = 0.4$, o babaligtarin ang sign ng $U$, nananatiling hindi nagbabago ang mga density sa max_occupation = 2 at mas mataas hanggang $10^{-13}$.

Ang nagtutulak pataas ng requirement ay ang bilang ng boson. Ang tatlong boson sa parehong chain ay nagbabago ng humigit-kumulang $2 \times 10^{-3}$ sa pagitan ng max_occupation = 2 at 3 pagkatapos ng parehong maikling evolution, at ng humigit-kumulang $0.3$ pagkatapos ng $t = 6$. Pahayag iyan tungkol sa model na ito at sa isang conserved na dami, at bound lang ito para sa mga operator na nagpapanatili ng bilang ng particle; hindi ito pangkalahatang tuntunin, kaya kailangang i-sweep ang cutoff, hindi hulaan.

May presyo rin ang sweep sa panig ng Circuit, na hindi ipinapakita ng eksaktong reference. Ang parehong bosonic Circuit, na ibinaba gamit ang Gray code sa bawat 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

Bawat dagdag na level ay may gastos sa depth, at tumatalon ang bilang ng qubit kapag lumampas ang $d$ sa power ng dalawa. Kapag nag-converge ka sa $d$, ikinakalakal mo ang accuracy laban sa laki ng Circuit, kaya higit pa sa numerical check ang sweep.

Pansinin din na ang pagiging halos tama ng $d = 2$ ay nagkataon lang dahil sa observable: halos hindi gumagalaw ang mga site density, habang qualitatively mali ang pinagbabatayang state (wala itong doubly-occupied na component). Kung nag-converge ang isang observable, hindi ito nagpapatunay na nag-converge na ang iba.

Ano ang dapat tandaan

• Ang encoding ay usapin ng transpiler. Walang binabanggit na truncation, encoding, o bilang ng qubit ang bosonic circuit sa step 3. Apat na encoding ang nag-lower nito sa apat na magkakaibang circuit, at pare-pareho ang physics na nai-reproduce ng lahat. Walang kailangang baguhin sa anumang nasa itaas ng transpiler.

• Para sa mga boson, hindi lang ito basta maginhawa. Ang truncation ay isang convergence parameter na kailangang i-sweep, at nasa loob ito ng encoding. Ang circuit na binuo laban sa isang fixed na encoding ay kailangang buuin muli para sa bawat halaga ng $d$.

• Mag-group bago mag-mapping, at mag-group para sa Hermiticity. Ang isang hopping term ay hindi Hermitian at hindi maaaring i-exponentiate; ang bond pair sa kamay ay Hermitian kaya maaari itong i-exponentiate. Piliin ang partition sa bosonic level, kung saan nakikita pa ang lattice. Kailangan ito ayon sa physics, at ang groups attribute ay nariyan para dalhin pasulong ang desisyong iyon.

• Huwag ipagpalagay na ang low-weight encoding ang mura. Sa $d = 3$, mas mahusay ang mga compact code kaysa sa parehong unary variant pagdating sa qubit at gate nang sabay, at ang projected variant ang pinakamahal sa lahat. Sukatin mo; nakadepende sa $d$ ang crossover.

• Sinisira ng truncation ang algebra, sa paraang masusukat. Sa pinakamataas na retained level lang, ng $d$. Kaya rin hindi nagko-commute ang mapping sa normal ordering.

Mga susunod na hakbang

May dalawang direksyon para palawakin ang halimbawang ito. Kung papalitan mo ang edge list ng chain ng mga bond ng isang square lattice, wala ka nang ibang babaguhin, dahil geometry-agnostic na ang bose_hubbard function na dinefine kanina. Ang sarili mong custom na encoding (halimbawa, isang block-unary code, na nag-iinterpolate sa pagitan ng unary at compact) ay maisasaksak din sa parehong MapperFnEvolutionSynthesis seam na ginagamit ng preset.

• Ipinapaliwanag ng Mga prinsipyo ng disenyo ng mga representasyon ng operator ang data structure ng operator sa likod ng Hamiltonian.

• Sinasaklaw ng I-map ang mga representasyon ng operator ang mga encoding at ang pagsulat ng sarili mong mapper.

• Ipinapakita ng I-transpile ang mga bosonic circuit ang mga stage na pinapatakbo ng preset.

• Ipinapakita ng Mag-synthesize ng Evolution gate sa bosonic space kung paano bawasan ang Trotter error sa bosonic space.

Mahalaga

Ang bawat iniulat na depth ay nasa abstract at fully connected na register. Ang pag-map sa tunay na hardware connectivity ay hiwalay na stage, at ang mga per-mode na bloke ng qubit ng isang occupation encoding ay nagbibigay sa router ng medyo ibang estruktura kaysa sa isang qubit bawat 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.