סימולציה של דינמיקת Fermi-Hubbard דו-ממדית עם קבוצות זרימה

חשוב

המושגים במדריך הזה זמינים כרגע רק ב-API של Python. פונקציונליות מקבילה תהיה זמינה ב-API של C בגרסה עתידית.

ראו גם

קראו קודם את המדריך לקבוצות זרימה חד-ממדיות. המדריך הזה הוא ההמשך הדו-ממדי שלו ולא חוזר על המידע שפותח שם: אופרטורי העברה וקודקוד, מהי קבוצת זרימה, איך groups שורד אל תוך המעגל, ואיך קידוד מותאם אישית מחובר ל-Transpiler דרך MapperFnEvolutionSynthesis.

ראו את מדריך הטרנספילציה לשלבי ה-Transpiler שמוזכרים לאורך המדריך, ואת מדריך המפרים למתכון הכללי לכתיבת מפר מותאם אישית.

המדריך החד-ממדי הסתיים בהזכרה שדרושים קידודים מקומיים כשעובדים בשני ממדים. המדריך הזה עושה את הצעד הזה, עבור מודל Fermi-Hubbard הדו-ממדי על סריג ריבועי, באמצעות קידוד Verstraete-Cirac (VC) [2].

יש שלושה הבדלים לעומת הבעיה החד-ממדית, והם בסיס העבודה הזו:

• קבוצות הזרימה מתרבות. לשרשרת יש שני כיווני חצים; לסריג ריבועי יש ארבעה (מזרח, מערב, צפון, דרום), וכל אחד הוא איחוד של מסלולים מכוונים זרים בקודקודים (אחד לכל שורה או עמודה בסריג) ולא מסלול יחיד.

• הקידוד מוסיף Qubit עזר לכל אתר, ולא אחד לכל המערכת, מה שנותן

  יחס Qubit-למוד של $2$. אופרטורי העברה הופכים לפאולים ממשקל 3 (אופקיים) וממשקל 4 (אנכיים), וכל אחת מארבע קבוצות הזרימה ממופה לצורת פאולי משלה, ולכן כל אחת צריכה Clifford מלכסן משלה.

• ההמילטוניאן המקודד נאמן רק על תת-מרחב: המרחב העצמי המשותף $+1$

  של אילוצי הפלקט שכל קידוד מקומי דו-ממדי נושא [2] [3]. בחד-ממד הקידוד היה איזומטריה פשוטה והאימות היה בדיקת ספקטרום. כאן ההשוואה צריכה להיות מוטלת, כפי ששלב 9 ממחיש.

התגמול הוא באותה צורה כמו בחד-ממד: עומק הדו-Qubit של צעד טרוטר אחד נעשה קבוע בגודל הסריג, ב-$24$ ($16$ לקפיצה, בהתאמה למאמר [1], ועוד $8$ לאינטראקציה), בעוד שטרוטריזציה איבר-איבר של אותו המילטוניאן מקודד גדלה ליניארית. מספר השערים עדיין גדל עם מספר הקשרים, כמובן; מה שנעשה בלתי תלוי בגודל הוא כמה שכבות השערים האלה תופסים. העומקים לאורך המדריך נספרים על רגיסטר מופשט ומחובר במלואו; מיפוי הקידוד על מפת הצימוד המוגבלת של מכשיר הוא בעיה בפני עצמה, שנידונה בסוף.

הסריג וארבעת הכיוונים שלו

השינוי הראשון נראה בגרף האינטראקציה המכוון. האתרים ממוספרים לפי שורות, site(r, c) = r * cols + c, וכל קשר נושא את שני הכיוונים, כפי שהשרשרת נשאה בחד-ממד, אבל עכשיו הקשרים עוברים לאורך שני צירים:

כלל ההחלפה הוא החד-ממדי, ללא שינוי: שני אופרטורי העברה שחולקים אתר מתחלפים כשהחצים זורמים דרכו ואנטי-מתחלפים כשהם מתנגשים. עם זאת, לאתר יש עכשיו עד ארבעה קשרים נכנסים במקום שניים, ולכן החצים דרכו יכולים להתנגש בעוד דרכים, וקבוצת חצים שלעולם לא מתנגשים כבר אינה מסלול יחיד. קריאת החצים האופקיים משמאל לימין נותנת קבוצה אחת כזו, החצים האנכיים מלמעלה למטה נותנים אחרת, והיפוך כל אחת נותן עוד שתיים. אלה קבוצות הזרימה מזרח, מערב, צפון ודרום שמארגנות את כל מה שבא אחר כך.

קידוד Verstraete-Cirac

Verstraete-Cirac (VC) מזווג כל מוד פרמיוני $j$ עם Qubit עזר משלו. עם $N$ אתרי סריג, מוד $j$ יושב על Qubit $j$ וה-Qubit העזר שלו על Qubit $N + j$, כך שהרגיסטר מחזיק $2N$ Qubits. אופרטורי קודקוד נשארים פשוטים ככל האפשר,

$$V_j = Z_j \, ,$$

וזה מה שהופך את האכלוסים לקריאים ישירות בהמשך. אופרטורי העברה תלויים בכיוון הקשת. עבור קשת אופקית $j \to k$ כש-$j$ האתר המערבי, וקשת אנכית כש-$j$ האתר הצפוני,

$$T_{jk}^{\text{horizontal}} = \tfrac{1}{2} X_j X_k Z_{N+j} \, , \qquad T_{jk}^{\text{vertical}} = \tfrac{1}{2} X_j Y_k X_{N+j} Y_{N+k} \, ,$$

“משולש” ממשקל 3 ו“ריבוע” ממשקל 4. זיהוי כל חץ חשוב: כש-$j$ הוא האתר המערבי, $j \to k$ מצביע מזרחה, אבל כש-$j$ הוא האתר הצפוני, הוא מצביע דרומה. הכיוונים ההפוכים נובעים מהזהות $T_{kj} = -V_j V_k T_{jk}$, שהופכת את משולש ה-XXZ ל-YYZ, ואת ריבוע ה-XYXY ל-YXXY. ארבע התוויות השונות האלה הן ארבע קבוצות הזרימה, והמדריך הזה עוקב אחריהן עד הסוף.

הערה

אלה לא בדיוק משוואה (A1) ממקור. [1], שקובעת "עד כדי סיבובי קיוביט בודד" ושונה בשלושה אופנים שכדאי לציין, כי מי שבודק את המאמר ישים לב אליהם:

• הגורם המקדים: משוואה (A1) לא נושאת גורם כזה. ה-$\tfrac12$ כאן הוא מה שגורם ל-$T_{jk}^2 = \tfrac14$, בהתאם לנרמול $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$ ש-TransferVertexOperator משתמש בו. המאמר צריך את אופרטורי ההעברה שלו רק בתוך $\exp(-i\,dt\,H)$, שם גורם קנה מידה נבלע בזווית הסיבוב. כאן הם צריכים להיות אופרטורים שמקיימים זהות.

• הכיוון: משוואה (A1) בוחרת אותיות לפי זוגיות השורה של $j$, שבקשת אופקית זהה לשני הקצוות, ולכן היא לא מבחינה בין $T_{jk}$ ל-$T_{kj}$. המאמר משחזר את הכיוון בנפרד, באמצעות צמידות עם שכבה של $R^Z(\pi/2)$. המדריך הזה במקום זה משלב אותו בקידוד, ומפיק את ההפוך מתוך $T_{kj} = -V_j V_k T_{jk}$, וזה מה שגורם לשני הכיוונים לנחות על פאולים שונים בשלב 3.

• הצורה האנכית: משוואה (A1) מחליפה בין XXXX ל-YYYY לפי זוגיות השורה, ואילו המדריך הזה משתמש ב-XYXY יחיד. אלה שתי דרכים להשיג את אותו הדבר: אי-סימטריה באותיות שגורמת לקשת אנכית להתנגש נכון עם הקשתות האופקיות שפוגשות אותה. אל תקראו את ההבדל כבחירת בסיס חופשית. אם מוותרים על החלפת הזוגיות ומשתמשים ב-XXXX יחיד, האלגברה נשברת לגמרי, ולכן האותיות השונות נבדקות ולא מניחים אותן מראש.

הערה

האנצילה מחוברת לאתר, מה שנותן יחס של $2$. קידוד Derby-Klassen הוא זה שמציב קיוביט עזר אחד לכל פלקט (אי-זוגי), ליחס קטן יותר של $3/2$. כאן משתמשים ב-VC כי אופרטורי ההעברה שלו הם אלה שמקור [1] נותן את הקליפורדים המלכסנים שלהם בצורה סגורה.

1. ההמילטוניאן על סריג ריבועי

הבנייה זהה לזו של 1D, אבל עם רשימת קשתות דו-ממדית. האתרים שומרים על המספור לפי שורות שבו השתמשנו באיור שלמעלה, site(r, c) = r * cols + c, כך שקשרים אופקיים מחברים אינדקסים עוקבים, וקשרים אנכיים מחברים אינדקסים שנבדלים ב-cols:

>>> from collections import defaultdict
>>> from qiskit_fermions.operators import TransferVertexOperator
>>>
>>> def fermi_hubbard_2d(rows, cols, tunneling, interaction):
...     site = lambda r, c: r * cols + c
...     edges = [(site(r, c), site(r, c + 1)) for r in range(rows) for c in range(cols - 1)]
...     edges += [(site(r, c), site(r + 1, c)) for r in range(rows - 1) for c in range(cols)]
...     data: defaultdict[tuple[tuple[int, int], ...], complex] = defaultdict(complex)
...     for j, k in edges:
...         # hopping: -t (a^dag_j a_k + a^dag_k a_j) = t (T_jk + T_kj)
...         data[((j, k),)] += tunneling
...         data[((k, j),)] += tunneling
...         # interaction: U n_j n_k with n_j = (1 - V_j) / 2
...         data[()] += interaction / 4
...         data[((j, j),)] -= interaction / 4
...         data[((k, k),)] -= interaction / 4
...         data[((j, j), (k, k))] += interaction / 4
...     return TransferVertexOperator.from_dict(data)
>>>
>>> rows, cols = 3, 3
>>> num_sites = rows * cols
>>> hamiltonian = fermi_hubbard_2d(rows, cols, tunneling=1.0, interaction=2.0)
>>> print(f"Number of terms: {len(hamiltonian)}")
Number of terms: 46

2. ארבע קבוצות זרימה במקום שתיים

כל איבר קפיצה מסווג לפי הכיוון שאליו החץ שלו מצביע על הסריג. השוואת קואורדינטות השורה/עמודה divmod של שני הקצוות מפרידה בין קשתות אופקיות לאנכיות, והסימן של ההפרש קובע את הכיוון. כל מה שאלכסוני (הזהות, האיברים $V_j$ ו-$V_j V_k$) נכנס לקבוצה חמישית:

הערה

מוסכמת הכיוון. השמות נקראים פשוט מהאיור למעלה: חץ הוא מזרח או מערב כשהוא מצביע ימינה או שמאלה לאורך שורה, וצפון או דרום כשהוא מצביע למעלה או למטה לאורך עמודה. מכיוון שהמספור רץ לפי שורות, מזרח ודרום הם הכיוונים שבהם האינדקס גדל. שום דבר בהמשך לא תלוי בשמות (הם רק מסמנים קבוצות), אבל קל לבלבל ביניהם, ולכן כדאי לקבע אותם מול התמונה. מקור [1] קורא לקבוצות הזרימה שלו על סריג ריבועי באותו אופן.

>>> EAST, WEST, NORTH, SOUTH, DIAGONAL = range(5)
>>>
>>> def flow_set_groups(operator, cols):
...     groups = []
...     for terms, _ in operator.iter_terms():
...         if len(terms) == 1 and terms[0][0] != terms[0][1]:
...             j, k = terms[0]
...             row_j, col_j = divmod(j, cols)
...             row_k, col_k = divmod(k, cols)
...             if row_j == row_k:
...                 groups.append(EAST if col_k > col_j else WEST)
...             else:
...                 groups.append(SOUTH if row_k > row_j else NORTH)
...         else:
...             groups.append(DIAGONAL)
...     return groups
>>>
>>> hamiltonian.groups = flow_set_groups(hamiltonian, cols)
>>>
>>> names = ["east", "west", "north", "south", "diagonal"]
>>> for index, name in enumerate(names):
...     count = sum(1 for group in hamiltonian.groups if group == index)
...     print(f"{name:9s} {count:2d}")
east       6
west       6
north      6
south      6
diagonal  22

בתוך קבוצת זרימה אחת כל חץ או פוגש חץ אחר ראש-אל-זנב או לא פוגש אותו בכלל, כמו ב-1D, ולכן קבוצת זרימה עדיין מתחלפת איבר-איבר וההתפתחות שלה לא נושאת שגיאת Trotter. החידוש הוא שקבוצת זרימה כעת לא קשירה: קבוצת המזרח של סריג rows x cols מורכבת מ-rows מסלולים עצמאיים משמאל לימין. הניתוק הזה הוא מתנה ולא סיבוך, כפי שיוצג בשלב 6.

3. קידוד האופרטור

הקידוד הוא פונקציה אחת שממפה אופרטור העברה מוכלל יחיד למחרוזת פאולי, שמועברת ל-map_transfer_vertex_generators(). בהשוואה לגרסת ה-1D, בגוף יש יותר מקרים: אלכסוני, אופקי, אנכי, והכיוונים ההפוכים דרך $T_{kj} = -V_j V_k T_{jk}$.

יש בעיה אחת. map_transfer_vertex_generators() קוראת חזרה עם מחולל ותו לא, בעוד שקביעת מחרוזת הפאולי של קשר דורשת גם את cols כדי לשחזר קואורדינטות שורה ועמודה מאינדקס מוד. במקום לקנן את הפונקציות כדי לסגור עליו, functools.partial() קושר אותו, וכך שתי הפונקציות נשארות ברמה העליונה וניתנות לקריאה עצמאית.

>>> from functools import partial
>>>
>>> from qiskit.quantum_info import SparseObservable
>>> from qiskit_fermions.mappers import map_transfer_vertex_generators
>>>
>>> def ancilla_qubit(mode, num_sites):
...     """The auxiliary qubit paired with fermionic mode `mode`."""
...     return num_sites + mode
>>>
>>> def vc_action(generator, width, cols):
...     """Map one generalized transfer operator to its Verstraete-Cirac Pauli."""
...     num_sites = width // 2
...     ancilla = partial(ancilla_qubit, num_sites=num_sites)
...     j, k = generator
...     if j == k:  # vertex operator V_j = Z_j
...         return SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     row_j, col_j = divmod(j, cols)
...     row_k, col_k = divmod(k, cols)
...     if row_j == row_k:  # horizontal: weight-3 triangle
...         left, right = (j, k) if col_j < col_k else (k, j)
...         operator = SparseObservable.from_sparse_list(
...             [("XXZ", [left, right, ancilla(left)], 0.5)], num_qubits=width
...         )
...         forward = (j, k) == (left, right)
...     else:  # vertical: weight-4 square
...         up, down = (j, k) if row_j < row_k else (k, j)
...         operator = SparseObservable.from_sparse_list(
...             [("XYXY", [up, down, ancilla(up), ancilla(down)], 0.5)], num_qubits=width
...         )
...         forward = (j, k) == (up, down)
...     if forward:
...         return operator
...     # reverse orientation: T_kj = -V_j V_k T_jk
...     v_j = SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     v_k = SparseObservable.from_sparse_list([("Z", [k], 1.0)], num_qubits=width)
...     return (-1.0 * v_j.compose(v_k).compose(operator)).simplify()
>>>
>>> def vc_encoding(operator, width, cols):
...     """Map a TransferVertexOperator to a qubit observable, Verstraete-Cirac style."""
...     return map_transfer_vertex_generators(
...         operator,
...         partial(vc_action, width=width, cols=cols),
...         identity=lambda: SparseObservable.identity(width),
...         compose=SparseObservable.compose,
...     ).simplify()
>>>
>>> # bind the lattice width once; the transpiler only ever passes (operator, width)
>>> encoding = partial(vc_encoding, cols=cols)
>>> num_qubits = 2 * num_sites
>>> num_qubits
18

את הערך של קיבוץ לפני המיפוי אפשר לראות מיד: כל קבוצת זרימה מקודדת לצורת פאולי אחת, אחידה בכל הקבוצה.

>>> for index, name in enumerate(names):
...     observable = encoding(hamiltonian.split_out_groups()[index], num_qubits)
...     shapes = sorted({label for label, _, _ in observable.to_sparse_list()})
...     print(f"{name:9s} {shapes}")
east      ['XXZ']
west      ['YYZ']
north     ['YXXY']
south     ['XYXY']
diagonal  ['', 'Z', 'ZZ']

ארבע צורות הקפיצה שונות זו מזו, ולכן כל קבוצת זרימה מקבלת קליפורד משלה ב-שלב 5. הכיוונים ההפוכים הם פאולים שונים, לא אותם פאולים בשם אחר. לכן סינתזה של קבוצות זרימה ב-2D דורשת ארבעה מעגלים, בעוד שב-1D היא דורשת אחד.

4. טרנספילציה עם הקידוד המותאם

השלב הזה הוא כמו ה-1D עם אוגר קיוביטים רחב יותר: CustomF2QLayout ממפה את אוגר $N$ המודים על $2N$ קיוביטים, MapperFnEvolutionSynthesis נושא את הקידוד, והארגומנט product_formula שלו בוחר איך ה-PauliEvolutionGate המתקבל מפורק. השוו עם שלב 6 במדריך ה-1D, שבונה את אותו צינור עבור אנצילה יחידה.

>>> from qiskit.circuit import QuantumRegister
>>> from qiskit.passmanager import MultiStagePassManager
>>> from qiskit.synthesis import LieTrotter
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution
>>> from qiskit_fermions.circuit.library.synthesis import FermionicLieTrotter
>>> from qiskit_fermions.transpiler import FermionicCircuitToDAG, QuantumDAGToCircuit
>>> from qiskit_fermions.transpiler.passes import (
...     CustomF2QLayout,
...     F2QSynthesis,
...     MapperFnEvolutionSynthesis,
... )
>>>
>>> total_time = 1.0
>>>
>>> def flow_set_pass_manager(mode_register, encoding, product_formula=None):
...     """Build the pipeline mapping `mode_register` with the VC `encoding`."""
...     synthesis = F2QSynthesis()
...     synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(
...         encoding, product_formula=product_formula
...     )
...     # this encoding spends one ancilla per fermionic mode
...     qubit_register = QuantumRegister(2 * mode_register.size, "q")
...     return MultiStagePassManager(
...         input=FermionicCircuitToDAG(),
...         layout=CustomF2QLayout({mode_register: qubit_register}),
...         synthesis=synthesis,
...         output=QuantumDAGToCircuit(),
...     )
>>>
>>> circuit = FermionicCircuit(num_sites)
>>> circuit.append(
...     Evolution(
...         num_sites, hamiltonian, time=total_time, synthesis=FermionicLieTrotter()
...     ),
...     circuit.modes,
... )
>>>
>>> # as in 1D, one decompose() splits the Evolution group by group (five gates here,
>>> # the four hopping flow sets in group order, then the diagonal interaction)
>>> flow_set_circuit = circuit.decompose()
>>> flow_set_circuit.count_ops()["Evolution"]
5

חמש שערים, שנושאים את 6, 6, 6, 6 ו-22 האיברים שנספרו בשלב 2. כל אחד ממופה ומסונתז בנפרד, כך שהחלוקה שנקבעה ברמה הפרמיונית שורדת עד לתוך המעגל:

>>> pass_manager = flow_set_pass_manager(flow_set_circuit.register, encoding, LieTrotter())
>>> qubit_circuit = pass_manager.run(flow_set_circuit)
>>> qubit_circuit.num_qubits
18
>>>
>>> def two_qubit_stats(circuit):
...     decomposed = circuit.decompose(reps=6)
...     two_qubit = lambda instruction: len(instruction.qubits) == 2
...     return decomposed.count_ops().get("cx", 0), decomposed.depth(two_qubit)
>>>
>>> two_qubit_stats(qubit_circuit)
(144, 47)

זה קו הבסיס: הקידוד והקיבוץ שניהם במקומם, אבל LieTrotter מסנתז כל קבוצת זרימה איבר-איבר, והעומק גדל עם הסריג.

5. קליפורד בעומק 2 לכל קבוצת זרימה

המנגנון זהה לזה של 1D: צמידה של קבוצת זרימה מתחלפת שלמה בקליפורד שממפה כל אחד מאיבריה לפאולי נפרד במשקל 1, כך שההתפתחות של הקבוצה כולה עולה רק סיבובי קיוביט בודד. ב-1D, brickwork אחד של שערי CZ עשה את זה.

ב-2D, כל קבוצת זרימה צריכה מעגל משלה, לפי סעיף IV C ואיור 7 של מקור [1]. כל הארבעה פועלים על קו אחד של הסריג בכל פעם (שורה עבור מזרח/מערב, עמודה עבור צפון/דרום) יחד עם האנצילות של אותו קו, ולכולם אותה צורה בת שלושה חלקים:

• שכבה של סיבובי קיוביט בודד, שמציבה כל קיוביט בבסיס שקבוצת הזרימה שלו צריכה.

• שכבת CX ראשונה, תמיד מאתר לאנצילה באותו אתר.

• שכבת CX שנייה, שמצמדת שכנים לאורך הקו.

שלבים 2 ו-3 הם שכבות שני הקיוביטים היחידות, ולכן עומק שני הקיוביטים הוא $2$, בלי קשר לאורך הקו. הן כתובות למטה בפונקציה אחת לכל קבוצת זרימה. הסיבובים הם מה שהם כי הם נמצאו, על ידי צמידת האיברים של כל קבוצת זרימה ודרישה לתמונות נפרדות במשקל 1, ולכן התייחסו אליהם כנתונים והסתמכו על הבדיקה שבאה אחריהם.

שתי הקבוצות האופקיות הן הפשוטות ביותר. האיברים שלהן במשקל 3, XXZ ו-YYZ, ונבדלים רק ב-$S$ הנוסף שהופך $X$ ל-$Y$ באתרים:

>>> from qiskit.circuit import QuantumCircuit
>>>
>>> def east_clifford(circuit, line, ancilla):
...     """Diagonalize the east flow set (XXZ) on one lattice row."""
...     for j in line:
...         circuit.h(ancilla(j))
...     for j in line[:-1]:  # layer 1: each site with its own ancilla
...         circuit.cx(ancilla(j), j)
...     for j, k in zip(line, line[1:]):  # layer 2: ancilla with the next site along
...         circuit.cx(k, ancilla(j))
>>>
>>> def west_clifford(circuit, line, ancilla):
...     """Diagonalize the west flow set (YYZ) on one lattice row."""
...     for j in line:
...         circuit.s(j)  # the only difference from east: X -> Y on the sites
...         circuit.h(ancilla(j))
...     for j in line[:-1]:
...         circuit.cx(ancilla(j), j)
...     for j, k in zip(line, line[1:]):
...         circuit.cx(k, ancilla(j))

שתי הקבוצות האנכיות הן במשקל 4, XYXY ו-YXXY, ובהתאם עמוסות יותר: הסיבובים מתחלפים לפי המיקום לאורך הקו, וגם שכבת ה-CX השנייה מתחלפת, ומצמדת אתרים במקומות זוגיים ואנצילות באי-זוגיים (או להפך):

>>> def south_clifford(circuit, line, ancilla):
...     """Diagonalize the south flow set (XYXY) on one lattice column."""
...     for position, j in enumerate(line):
...         if position % 2 == 0:
...             circuit.sdg(ancilla(j))
...             circuit.h(ancilla(j))
...         else:
...             circuit.h(j)
...             circuit.s(j)
...             circuit.h(ancilla(j))
...             circuit.s(ancilla(j))
...     for j in line:  # layer 1
...         circuit.cx(ancilla(j), j)
...     for position, (j, k) in enumerate(zip(line, line[1:])):  # layer 2
...         if position % 2 == 0:
...             circuit.cx(ancilla(j), ancilla(k))  # ancilla to ancilla
...         else:
...             circuit.cx(j, k)  # site to site
>>>
>>> def north_clifford(circuit, line, ancilla):
...     """Diagonalize the north flow set (YXXY) on one lattice column."""
...     for position, j in enumerate(line):
...         if position % 2 == 0:
...             circuit.h(ancilla(j))
...             circuit.s(ancilla(j))
...         else:
...             circuit.h(j)
...             circuit.s(j)
...             circuit.sdg(ancilla(j))
...             circuit.h(ancilla(j))
...     for j in line:  # layer 1
...         circuit.cx(ancilla(j), j)
...     for position, (j, k) in enumerate(zip(line, line[1:])):  # layer 2
...         if position % 2 == 0:
...             circuit.cx(k, j)  # site to site, and the other way round from south
...         else:
...             circuit.cx(ancilla(j), ancilla(k))

כשמציירים אותם על קו יחיד של שלושה אתרים (קודם האתרים ואז האנצילות שלהם), כל המנגנון נכנס בתמונה אחת. צורת המשולש היא הברורה מבין השתיים: שינוי בסיס על האנצילות, ואז שתי שכבות CX, הראשונה מזווגת כל אתר עם האנצילה שלו, והשנייה מגיעה לאתר הבא לאורך הקו.

>>> line = [0, 1, 2]  # three sites, with ancillas on qubits 3, 4, 5
>>> triangle = QuantumCircuit(6)
>>> east_clifford(triangle, line, lambda mode: 3 + mode)
>>> triangle.draw("mpl")
<Figure size ... with 1 Axes>

צורת הריבוע צריכה יותר עבודת קיוביט בודד (הסיבובים מתחלפים לאורך הקו) ושכבת ה-CX השנייה שלה מתחלפת בין צימוד אתרים לצימוד אנצילות:

>>> square = QuantumCircuit(6)
>>> south_clifford(square, line, lambda mode: 3 + mode)
>>> square.draw("mpl")
<Figure size ... with 1 Axes>

שתיהן בעומק של שתי שכבות CX, והוספת אתרים לקו מרחיבה אותן בלי להעמיק אותן. זו כל טענת העומק הקבוע, שנראית לעין. הקליפורדים של מערב וצפון נבדלים רק בשכבות הקיוביט הבודד שלהם, ולכן הם לא מצוירים בנפרד.

מה שהסינתזה עדיין צריכה הוא לקבוע לאן כל איבר מגיע. קו של $L$ אתרים נושא $L - 1$ קשרים, וכל אחד תורם איבר אחד, ולכן הפונקציה העוזרת למטה מחזירה רשומה אחת לכל קשר: הקיוביט שעליו האיבר נחת, ופאולי הקיוביט הבודד שהוא הפך להיות. חשוב: אלה נקראים מהצמידה ולא מנוחשים מראש. בניגוד ל-1D, שבו brickwork אחד של CZ ממפה $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ בצורה סגורה, כאן אין נוסחה כזו, והאנצילות של קו עמודה אינן רציפות. הדבר היחיד שנבדל בין קבוצות הזרימה הוא מחרוזת הפאולי שנכנסת, ששלב 3 כבר הדפיס:

>>> from qiskit.quantum_info import Clifford, SparsePauliOp
>>>
>>> def clifford_supports(clifford, pauli, line, ancilla, num_qubits):
...     """Conjugate a flow set's terms to find which qubit each one lands on.
...
...     `pauli` is the flow set's Pauli string, as printed in step 3. Its letters apply
...     to the i-th term's qubits in the order below: the two sites of the bond, then
...     their ancillas.
...
...     This is the one check that needs SparsePauliOp: Clifford conjugation wants a
...     `Pauli` and its symplectic bits, and SparseObservable exposes neither -- its own
...     `evolve` conjugates by a Pauli, not by a Clifford.
...     """
...     supports = []
...     for i in range(len(line) - 1):
...         qubits = [line[i], line[i + 1], ancilla(line[i]), ancilla(line[i + 1])]
...         term = SparsePauliOp.from_sparse_list(
...             [(pauli, qubits[: len(pauli)], 1.0)], num_qubits=num_qubits
...         )
...         image = term.paulis[0].evolve(clifford, frame="s")
...         # a weight-1 image has exactly one qubit with an X and/or Z bit set
...         (support,) = np.nonzero(image.x | image.z)
...         assert len(support) == 1, f"term {i} is not weight-1"
...         supports.append((int(support[0]), image[support].to_label()))
...     return supports

התמונות חייבות לצאת נפרדות ובמשקל 1, כך שאפשר לסובב כל איבר באופן עצמאי על הקיוביט שלו. עומק שני הקיוביטים של $2$ כבר נראה בארבע הבניות שלמעלה, שכל אחת מהן פולטת שתי שכבות CX ללא קשר לאורך הקו. למטה שתי התכונות נבדקות על פני כל קבוצת זרימה, כל קו וכמה גדלי סריג. שלב 6 קורא אחר כך ל-clifford_supports, כך שזה חלק עובד מהסינתזה ולא רק הדגמה.

6. סינתזה מודעת לקבוצות זרימה

הסינתזה מקבלת יתרון פתיחה שימושי. מכיוון ש-Evolution מפצלת את עצמה קבוצה אחר קבוצה, synthesize נקראת פעם אחת לכל קבוצת זרימה (חמש פעמים כאן, ארבע קבוצות קפיצה ועוד הקבוצה האלכסונית), כך שלכל איבר קפיצה שהיא רואה יש אותה צורת פאולי, וקליפורד אחד מלכסן את כולם. היא רק צריכה לקרוא את הצורה הזו מהאיבר הראשון ולמיין את השאר לפי הקו שעליו הם יושבים. ואז מגיעה התובנה שהופכת שני ממדים לזולים יותר מכפי שאפשר לחשוש:

חשוב

הרכיבים של קבוצת זרימה יחידה זרים בקודקודים: השורות של קבוצת המערב לא חולקות קיוביטים. לכן אפשר לפלוט את כל הקליפורדים של קבוצת זרימה יחד והם חולקים את אותן שתי שכבות שני קיוביטים, במקום לעמוד בתור שורה אחרי שורה. קבוצת זרימה עולה אפוא מספר קבוע של שכבות (שתיים עבור הקליפורד ושתיים עבור ההופכי שלו) לא משנה כמה שורות או עמודות יש לסריג, וזו הסיבה לכך שהעומק הכולל יוצא קבוע.

>>> from qiskit.synthesis.evolution import EvolutionSynthesis
>>>
>>> # The interaction couples nearest neighbors only, so every interior site carries four
>>> # Rzz gates -- one per bond -- and can only do one at a time: four layers is the best
>>> # any schedule can achieve. Horizontal bonds split by column parity, vertical ones by
>>> # row parity, which achieves it.
>>> def rzz_layers(circuit, angles, cols):
...     """Emit the diagonal Rzz gates in four layers of disjoint qubit pairs."""
...     horizontal = [(j, k) for j, k in angles if k - j == 1]
...     vertical = [(j, k) for j, k in angles if k - j == cols]
...     schedule = [
...         [(j, k) for j, k in horizontal if (j % cols) % 2 == 0],
...         [(j, k) for j, k in horizontal if (j % cols) % 2 == 1],
...         [(j, k) for j, k in vertical if (j // cols) % 2 == 0],
...         [(j, k) for j, k in vertical if (j // cols) % 2 == 1],
...     ]
...     # emit layer by layer, since Qiskit's depth follows the order gates were added
...     for layer in schedule:
...         for j, k in layer:
...             circuit.rzz(angles[j, k], j, k)
>>>
>>> class FlowSetSynthesis(EvolutionSynthesis):
...     """Synthesize one 2D line flow set at a time, via its depth-2 Clifford."""
...
...     preserve_order = True
...
...     def __init__(self, rows, cols):
...         self.rows, self.cols = rows, cols
...
...     def rows_of_lattice(self):
...         return [[r * self.cols + c for c in range(self.cols)] for r in range(self.rows)]
...
...     def columns_of_lattice(self):
...         return [[r * self.cols + c for r in range(self.rows)] for c in range(self.cols)]
...
...     def synthesize(self, evolution):
...         num_qubits = 2 * self.rows * self.cols
...         ancilla = partial(ancilla_qubit, num_sites=self.rows * self.cols)
...         circuit = QuantumCircuit(evolution.operator.num_qubits)
...         terms = [
...             (label, indices, 2 * evolution.time * coeff.real)
...             for label, indices, coeff in evolution.operator.to_sparse_list()
...         ]
...
...         # `Evolution.decompose()` hands us one gate per flow set, so every hopping term
...         # here shares a single Pauli shape and one Clifford diagonalizes all of them
...         shapes = tuple({label for label, _, _ in terms} - {"", "Z", "ZZ"})
...         match shapes:
...             case ("XXZ",):
...                 clifford_function, lines = east_clifford, self.rows_of_lattice()
...                 along_rows = True
...             case ("YYZ",):
...                 clifford_function, lines = west_clifford, self.rows_of_lattice()
...                 along_rows = True
...             case ("XYXY",):
...                 clifford_function, lines = south_clifford, self.columns_of_lattice()
...                 along_rows = False
...             case ("YXXY",):
...                 clifford_function, lines = north_clifford, self.columns_of_lattice()
...                 along_rows = False
...             case ():  # the diagonal group: no conjugation needed
...                 self.diagonal(circuit, terms)
...                 return circuit
...             case _:
...                 raise ValueError(f"expected one known flow set per gate, got {shapes}")
...
...         # a term's line is the row it sits in and its slot the column, or vice versa
...         rotations: dict[int, list] = defaultdict(list)
...         for _, indices, angle in terms:
...             row, col = divmod(indices[0], self.cols)
...             line, slot = (row, col) if along_rows else (col, row)
...             rotations[line].append((slot, angle))
...
...         active = sorted(rotations)
...         # the lines of one flow set are vertex-disjoint, so all of their Cliffords
...         # are emitted together and share the same two two-qubit layers
...         conjugation = QuantumCircuit(num_qubits)
...         for line in active:
...             clifford_function(conjugation, lines[line], ancilla)
...         circuit.compose(conjugation, inplace=True)
...         for line in active:
...             supports = clifford_supports(
...                 Clifford(conjugation), shapes[0], lines[line], ancilla, num_qubits
...             )
...             for slot, angle in rotations[line]:
...                 match supports[slot]:
...                     case qubit, "X":
...                         circuit.rx(angle, qubit)
...                     case qubit, "Y":
...                         circuit.ry(angle, qubit)
...                     case qubit, "Z":
...                         circuit.rz(angle, qubit)
...         circuit.compose(conjugation.inverse(), inplace=True)
...         return circuit
...
...     def diagonal(self, circuit, terms):
...         """Synthesize the diagonal group: a phase, single-qubit Z and paired Rzz."""
...         rzz = {}
...         for label, indices, angle in terms:
...             match label, indices:
...                 case "", _:  # the identity contributes only a global phase
...                     circuit.global_phase -= angle / 2
...                 case "Z", [qubit]:  # a lone vertex operator
...                     circuit.rz(angle, qubit)
...                 case "ZZ", [j, k]:  # interaction, see rzz_layers above
...                     rzz[min(j, k), max(j, k)] = angle
...                 case _:
...                     raise ValueError(f"unexpected Pauli term: {label} on {indices}")
...         rzz_layers(circuit, rzz, self.cols)

כמו ב-1D, זה נכתב במכוון עבור הקידוד הזה: הוא מזהה את ארבע צורות הפאולי וזורק שגיאה על כל דבר אחר במקום ליפול חזרה בשקט. הוא גם מניח את הפיצול של קבוצת זרימה אחת לכל שער מלמעלה, ואומר זאת. יותר מצורת קפיצה אחת בשער יחיד אומר שהקיבוץ לא הוחל, וזו שגיאה ולא משהו שצריך לעקוף.

הערה

כמו מקבילתה ב-1D, הסינתזה הזו לא מקבלת מספר צעדי Trotter. צעד Trotter צריך לעבור בין קבוצות הזרימה, אבל transpile_trotter_step למטה מפרק את המעגל הפרמיוני לשער Evolution אחד לכל קבוצה, ולכן השיטה הזו רואה רק קבוצה אחת בכל פעם ויכולה לחזור רק בתוכה, מה שלא משנה דבר: חזרה על אופרטור אחד רק מפצלת את ההתפתחות שלו לחלקים שווים.

לכן Trotterization מרובה צעדים צריך נוסחת מכפלה שמכבדת את ה-groups של האופרטור בתוך שער יחיד, לפני הפירוק. FermionicLieTrotter עושה את זה, ו-שלב 10 משתמש בו כדי למדוד את ההתכנסות.

שימו לב לשתי המחלקות בעלות השמות הדומים שמעורבות כאן. LieTrotter של Qiskit, ששימשה כקו הבסיס למעלה, פועלת על ההמילטוניאן הממופה, שבו הקיבוץ כבר אבד. FermionicLieTrotter היא המקבילה הפרמיונית שלה, והיא זו שעדיין רואה את קבוצות הזרימה.

7. התוצאה: עומק קבוע

כל מה שעד כה רץ על סריג $3 \times 3$ יחיד שנבנה בשלב 1. השוואת שתי דרכי הסינתזה על פני סריגים הולכים וגדלים דורשת לחזור על כל הרצף הזה (אופרטור, קיבוץ, קידוד, מעגל, מנהל מעברים) עבור כל גודל, ולכן כדאי לעטוף אותו פעם אחת:

>>> def transpile_trotter_step(rows, cols, product_formula, steps=1):
...     """Transpile `steps` Trotter steps on a `rows` x `cols` lattice."""
...     num_sites = rows * cols
...     operator = fermi_hubbard_2d(rows, cols, 1.0, 2.0)
...     operator.groups = flow_set_groups(operator, cols)
...     encoding = partial(vc_encoding, cols=cols)
...
...     circuit = FermionicCircuit(num_sites)
...     circuit.append(
...         Evolution(
...             num_sites,
...             operator,
...             time=total_time,
...             synthesis=FermionicLieTrotter(reps=steps),
...         ),
...         circuit.modes,
...     )
...     circuit = circuit.decompose()  # split into one Evolution per flow set
...
...     pass_manager = flow_set_pass_manager(circuit.register, encoding, product_formula)
...     return pass_manager.run(circuit)

איסוף שני המדדים על פני טווח של גדלי סריג והצגתם מול קצה הסריג מבהיר את ההבדל:

>>> sizes = [3, 4, 5, 6, 7, 8]
>>> lie_trotter = [
...     two_qubit_stats(transpile_trotter_step(n, n, LieTrotter())) for n in sizes
... ]
>>> flow = [
...     two_qubit_stats(transpile_trotter_step(n, n, FlowSetSynthesis(n, n))) for n in sizes
... ]
>>>
>>> figure, axes = plt.subplots(1, 2, figsize=(9, 3.5), layout="constrained")
>>> for panel, index, title in zip(axes, (1, 0), ("two-qubit depth", "CX count")):
...     _ = panel.plot(sizes, [stat[index] for stat in lie_trotter], "o-", label="qubit-level")
...     _ = panel.plot(sizes, [stat[index] for stat in flow], "s-", label="flow set")
...     _ = panel.set(xlabel="lattice edge $L$ (an $L \\times L$ lattice)", title=title)
...     _ = panel.legend()
>>> figure
<Figure size ... with 2 Axes>

עומק שני הקיוביטים הוא 24 וקבוע מ-$3 \times 3$ עד $8 \times 8$ (כלומר מ-18 קיוביטים עד 128), בעוד ש-LieTrotter שמוחל על ההמילטוניאן המקודד באופן זהה גדל ליניארית עם קצה הסריג, מ-47 עד 92. מספר השערים גדל בשתיהן, כי יש $O(N)$ איברים להחיל בכל מקרה. מה שקבוצות הזרימה מסירות הוא הגידול במספר השכבות שהשערים האלה תופסים.

התקציב מתחלק ל-16 לקפיצה ושמונה לאינטראקציה:

• כל אחת מארבע קבוצות הזרימה עולה $2 + 2$ שכבות שני קיוביטים: הקליפורד שלה בכניסה וההופכי שלו ביציאה, שניהם בעומק שניים מעצם הבנייה, עם סיבובי קיוביט בודד בלבד ביניהם. זה $4 \times 4 = 16$, והוא הנתון שמקור [1] מצטט, מכיוון שהוא סופר רק את איברי הקפיצה. הצבת $U = 0$ בקוד שלמעלה משחזרת אותו.

• הקבוצה האלכסונית מוסיפה את שמונה הנותרים. שערי ה-Rzz שלה נארזים בארבע שכבות, וכל Rzz מתפרק לשני שערי CX סביב סיבוב קיוביט בודד, ולכן $4 \times 2 = 8$. ארבע שכבות זה אופטימלי: על סריג ריבועי האינטראקציה מגיעה רק לשכנים הקרובים, ולכן בכל אתר פנימי יש ארבעה שערי Rzz (אחד לכל קשר) והוא יכול להשתתף רק באחד מהם בכל פעם.

אף אחד משני הסכומים לא מתייחס לגודל הסריג, וזו כל הנקודה.

8. אימות המעגל

סינתזת קבוצות הזרימה ו-LieTrotter מתארות את אותו צעד Trotter של אותו המילטוניאן מקודד, ולכן היוניטריים שלהן חייבים להסכים עד כדי פאזה גלובלית. וקטורי מצב אקראיים מגלים כל אי-התאמה בהרבה פחות מאמץ מבניית היוניטרי המלא:

הערה

שני המעגלים מפורקים לשערים של קיוביט אחד ושני קיוביטים לפני שהם מוחלים על מצב. אם משאירים אותם לא מפורקים, כל הוראת Evolution ברמה גבוהה מאלצת את Statevector לבנות את האופרטור הצפוף המלא שלה בגודל $2^{2N} \times 2^{2N}$. פירוק מראש מאפשר לו להחיל שערים אחד אחד במקום זאת, וזה מהיר כאן פי כמה מאות.

>>> import numpy as np
>>> from qiskit.quantum_info import Statevector
>>>
>>> rng = np.random.default_rng(1234)
>>>
>>> def agrees_up_to_one_phase(rows, cols, samples=3):
...     """Compare both syntheses on random states.
...
...     Returns whether they agree up to a phase, and whether it is the *same* phase
...     for every input state.
...     """
...     lie_trotter = transpile_trotter_step(
...         rows, cols, LieTrotter()
...     ).decompose(reps=6)
...     flow = transpile_trotter_step(
...         rows, cols, FlowSetSynthesis(rows, cols)
...     ).decompose(reps=6)
...     dimension = 2**lie_trotter.num_qubits
...     residuals, phases = [], []
...     for _ in range(samples):
...         vector = rng.normal(size=dimension) + 1j * rng.normal(size=dimension)
...         state = Statevector(vector / np.linalg.norm(vector))
...         a = state.evolve(lie_trotter).data
...         b = state.evolve(flow).data
...         overlap = np.vdot(a, b)
...         phases.append(overlap / abs(overlap))
...         residuals.append(np.linalg.norm(b - a * overlap / abs(overlap)))
...     # the SAME global phase must work for every input state
...     spread = max(abs(phase - phases[0]) for phase in phases)
...     return bool(max(residuals) < 1e-9), bool(spread < 1e-9)
>>>
>>> agrees_up_to_one_phase(rows, cols)
(True, True)

הערה

שני חצאי הבדיקה הזו חשובים. חפיפה בעלת מודולוס אחד על מצב יחיד היא תנאי הכרחי אבל לא מספיק: סינתזה שטעתה בפאזה התלויה במצב עדיין הייתה מראה זאת. הדרישה לאותה פאזה על פני כמה מצבים אקראיים היא מה שסוגר את הפער.

9. אימות הפיזיקה על תת-המרחב היציבילי

זו הקרחת הדו-ממדית. ב-1D הקידוד היה איזומטריה, והשוואת ספקטרום מול ייחוס Jordan-Wigner הספיקה. VC אינו איזומטריה: הוא משתמש ב-$2N$ קיוביטים עבור $N$ מודים, כך שמרחב בעל $2^{2N}$ ממדים מכיל מרחב פיזיקלי בעל $2^N$ ממדים. זה כללי לקידודים מקומיים בשני ממדים (הקיוביטים העודפים מוגבלים על ידי סימטריות מתחלפות מקומיות, אחת לכל פלקט [2] [3]) וזה אומר שההמילטוניאן המקודד מייצג את הפרמיוני בנאמנות רק על המרחב העצמי המשותף של $+1$ של האילוצים האלה.

אזהרה

המילה stabilizer (מייצב) משרתת שתי מטרות בתחום הזה, ושתי המשמעויות אינן קשורות. מקור [1] משתמש בה עבור אופרטורי ההעברה המתחלפים של קבוצת זרימה, שיוצרים חבורת מייצבים שהקליפורדים של שלב 5 מלכסנים (התקן קומפילציה, והמשמעות בכל מקום למעלה). הסעיף הזה משתמש בה במובן של תיקון שגיאות: אופרטורים שהמרחב העצמי של $+1$ שלהם הוא המרחב הפיזיקלי. מקור [1] מסמן במפורש את ההתנגשות הזו ומניח בצד את המשמעות השנייה. אילוצי הפלקט שלמטה הם תכונה של קידוד VC [2], לא תוצאה של מסגרת קבוצות הזרימה.

המחוללים בנויים כאן ולא מצוטטים: אחד לכל פאה ריבועית, וכל אחד הוא פשוט לולאת אופרטורי הקשת המקודדים סביב אותה פאה.

אופרטורי קשת הם האובייקט הטבעי לכך. לולאה סגורה שלהם היא סימטריה של כל המילטוניאן ששומר על מספר החלקיקים, בעוד שלולאת אופרטורי ההעברה המתאימה אינה כזו: ההפרש ביניהם הוא אופרטור קודקוד על אתר המקור של כל קשר,

$$T_{jk} = \tfrac{i}{2} \, V_j E_{jk} \qquad \Longleftrightarrow \qquad E_{jk} = -2i \, V_j T_{jk} \, ,$$

סביב לולאה סגורה כל אתר הוא המקור של קשר אחד בדיוק, ולכן הלולאה אוספת $V_j$ אחד לכל אתר, שהוא הזוגיות המקומית של ארבעת אתרי הפלקט. הזוגיות המקומית היא מה שעושה את ההבדל: ההמילטוניאן מקפיץ חלקיקים מעל גבול הפלקט, ולכן הוא שומר רק על הזוגיות הגלובלית. לולאת ההעברה החשופה לכן לא מתחלפת איתו, בעוד שלולאת הקשת (שנושאת את גורם הזוגיות) מתחלפת.

מכיוון שהקידוד מוגדר על אופרטורי העברה, מיפוי אופרטור קשת צריך את היחס שלמעלה. זו פונקציית callback בת חמש שורות, והיא משתלבת ב-map_edge_vertex_generators() בדיוק כמו ש-vc_action משתלבת ב-map_transfer_vertex_generators():

>>> from qiskit_fermions.mappers import map_edge_vertex_generators
>>> from qiskit_fermions.operators import EdgeVertexOperator
>>>
>>> def vc_edge_action(generator, width, cols):
...     """Map one generalized *edge* operator to its Verstraete-Cirac Pauli."""
...     j, k = generator
...     if j == k:  # a vertex operator is the same object in both algebras
...         return SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     # The library decomposes T_jk = (i/2) V_j E_jk, so E_jk = -2i V_j T_jk. Reuse the
...     # transfer encoding for T_jk rather than repeating its four cases here. Note that
...     # `compose` left-multiplies its argument, so this composes to V_j T_jk; the other
...     # order carries the opposite sign, since a vertex and a transfer operator sharing a
...     # mode anticommute.
...     vertex = SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     return (-2.0j * vc_action(generator, width, cols).compose(vertex)).simplify()
>>>
>>> def vc_edge_encoding(operator, width, cols):
...     """Map an EdgeVertexOperator to a qubit observable, Verstraete-Cirac style."""
...     return map_edge_vertex_generators(
...         operator,
...         partial(vc_edge_action, width=width, cols=cols),
...         identity=lambda: SparseObservable.identity(width),
...         compose=SparseObservable.compose,
...     ).simplify()
>>>
>>> def plaquette_loops(rows, cols):
...     """The directed edge loop around each square face."""
...     at = lambda r, c: r * cols + c
...     return [
...         [(at(r, c), at(r, c + 1)), (at(r, c + 1), at(r + 1, c + 1)),
...          (at(r + 1, c + 1), at(r + 1, c)), (at(r + 1, c), at(r, c))]
...         for r in range(rows - 1)
...         for c in range(cols - 1)
...     ]
>>>
>>> def stabilizers(rows, cols):
...     """One plaquette stabilizer per square face of the lattice."""
...     num_qubits = 2 * rows * cols
...     out = []
...     for loop in plaquette_loops(rows, cols):
...         loop_operator = EdgeVertexOperator.from_dict({tuple(loop): 1.0})
...         encoded_loop = vc_edge_encoding(loop_operator, num_qubits, cols)
...         out.append(encoded_loop / complex(encoded_loop.coeffs[0]))
...     return out
>>>
>>> plaquette_stabilizers = stabilizers(rows, cols)
>>> for stabilizer in plaquette_stabilizers:
...     print(stabilizer.to_sparse_list())
[('ZZYXXY', [1, 4, 9, 10, 12, 13], (1+0j))]
[('ZZYXXY', [2, 5, 10, 11, 13, 14], (1+0j))]
[('ZZYXXY', [4, 7, 12, 13, 15, 16], (1+0j))]
[('ZZYXXY', [5, 8, 13, 14, 16, 17], (1+0j))]

לכל מחולל יש משקל $6$, פרוש על שני קיוביטים פיזיקליים של הפלקט שלו וארבע אנצילות. הם אינוולוטוריים ומתחלפים עם ההמילטוניאן המקודד, ולכן המרחב הפיזיקלי נשמר על ידי הדינמיקה. שתי התכונות הן הצהרות על אופרטורים, ולכן הן נבדקות סימבולית על האלגברה של SparseObservable, בלי צורך במטריצות $2^{18} \times 2^{18}$:

>>> encoded = encoding(hamiltonian, num_qubits)
>>> identity = SparseObservable.identity(num_qubits)
>>> for stabilizer in plaquette_stabilizers:
...     square = stabilizer.compose(stabilizer).simplify()
...     print(square == identity, stabilizer.commutes(encoded))
True True
True True
True True
True True

עם המייצבים ביד, ההשוואה לא צריכה התאמת בסיסים, כי $V_j = Z_j$ אומר שהאכלוסים חיים על אותם קיוביטים פיזיקליים בשתי התמונות. אז: כותבים מצב מכפלה ישירות על הקיוביטים הפיזיקליים עם האנצילות ב-$|0\rangle$, מטילים אותו על תת-מרחב המייצבים, מפתחים עם ההמילטוניאן המקודד, וקוראים את צפיפויות האתרים $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$ משני הצדדים.

כל שלב כאן נשאר ברמת וקטור. ההטלה מוחלת גורם $(1 + S)/2$ אחד בכל פעם במקום להיכפל החוצה, וההתפתחות בזמן משתמשת ב-expm_multiply(), שמחשבת את $e^{-iHt} |\psi\rangle$ בלי לבנות אי פעם את $e^{-iHt}$. מכיוון שההמילטוניאן המקודד הוא אופרטור דליל בעל $2^{2N}$ ממדים, ההבחנה הזו היא מה ששומר על הסעיף הזה זול.

>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_fermions.mappers.library import transfer_vertex_jordan_wigner
>>>
>>> # an exact evolution time, unrelated to the Trotter step above
>>> evolution_time = 0.6
>>> reference = transfer_vertex_jordan_wigner(hamiltonian, num_sites).simplify()
>>>
>>> def to_sparse_matrix(observable):
...     # SparseObservable carries no dense- or sparse-matrix method, so route through
...     # SparsePauliOp. This is only ever needed for the numerical checks below --
...     # the encoding itself never leaves SparseObservable.
...     return SparsePauliOp.from_sparse_observable(observable).to_matrix(sparse=True)
>>>
>>> def basis_state(occupation, num_qubits):
...     """The computational basis state with `occupation` on the first qubits.
...
...     Qubit `q` carries bit `2**q`, and the ancillas are the *high* qubits `N..2N-1`,
...     so the same integer indexes the same site occupations on either register and
...     leaves every ancilla in |0>.
...     """
...     index = int("".join(str(bit) for bit in reversed(occupation)), 2)
...     state = np.zeros(2**num_qubits, dtype=complex)
...     state[index] = 1.0
...     return state
>>>
>>> def project_onto_codespace(state, stabilizers):
...     """Apply the projector one (1 + S) factor at a time, then renormalize.
...
...     The factors commute, so their product is itself a projector and the order does
...     not matter. Multiplying it out would build a dense 2^2N x 2^2N operator; applied
...     factor by factor this stays a handful of sparse matrix-vector products.
...
...     Each factor is really (1 + S)/2, but every one contributes the same scalar, so the
...     halvings only rescale the result and the final normalization absorbs them.
...     """
...     for stabilizer in stabilizers:
...         state += to_sparse_matrix(stabilizer) @ state
...     return state / np.linalg.norm(state)
>>>
>>> occupation = [1, 0, 1, 0, 1, 0, 0, 1, 0]  # site occupations of the initial state
>>> initial_reference = basis_state(occupation, num_sites)
>>> initial_encoded = project_onto_codespace(
...     basis_state(occupation, num_qubits), plaquette_stabilizers
... )
>>>
>>> evolved_encoded = expm_multiply(
...     -1j * evolution_time * to_sparse_matrix(encoded), initial_encoded
... )
>>> evolved_reference = expm_multiply(
...     -1j * evolution_time * to_sparse_matrix(reference), initial_reference
... )
>>>
>>> def density(state, qubit, width):
...     """Site density <n_j> = (1 - <Z_j>) / 2 of a statevector."""
...     z = to_sparse_matrix(
...         SparseObservable.from_sparse_list([("Z", [qubit], 1.0)], num_qubits=width)
...     )
...     return (1 - (state.conj() @ (z @ state)).real) / 2
>>>
>>> encoded_densities = [density(evolved_encoded, j, num_qubits) for j in range(num_sites)]
>>> reference_densities = [density(evolved_reference, j, num_sites) for j in range(num_sites)]
>>>
>>> for j, density_value in enumerate(encoded_densities):
...     print(f"site {j}: {density_value:.10f}")
site 0: 0.6750993356
site 1: 0.3689791258
site 2: 0.6750993356
site 3: 0.3818130396
site 4: 0.4868585996
site 5: 0.3818130396
site 6: 0.2693544696
site 7: 0.4916285849
site 8: 0.2693544696
>>> np.allclose(encoded_densities, reference_densities, atol=1e-12)
True

אלה הצפיפויות המקודדות. הן תואמות את ייחוס Jordan-Wigner עד דיוק מכונה, על סריג שהקידוד שלו נאמן רק אחרי הטלה.

10. מדידת ההתכנסות של Trotter

שלב 8 קבע ששתי הסינתזות מסכימות זו עם זו, וזו אמירה על שקילות ולא על דיוק: שתיהן מתארות את אותו צעד Trotter יחיד, ולכן שתיהן נושאות את אותה שגיאת Trotter. מדידת השגיאה הזו צריכה ייחוס מדויק, ושלב 9 בנה כזה זה עתה.

מספר הצעדים מגיע מ-FermionicLieTrotter, המקבילה הפרמיונית של הנוסחה ברמת הקיוביט ששימשה קודם כקו בסיס. כשהוא מוחל על כל האופרטור המקובץ לפני הפירוק, החזרה שלו עוברת בין ארבע קבוצות הזרימה והקבוצה האלכסונית ולא בתוך אחת מהן:

>>> initial_state = project_onto_codespace(
...     basis_state(occupation, num_qubits), plaquette_stabilizers
... )
>>> exact = expm_multiply(
...     -1j * total_time * to_sparse_matrix(encoded), initial_state
... )
>>>
>>> def trotter_fidelity(steps):
...     # `steps` goes to the fermionic formula inside transpile_trotter_step, so the
...     # repetition cycles between the flow sets rather than within any one of them
...     circuit = transpile_trotter_step(
...         rows, cols, FlowSetSynthesis(rows, cols), steps=steps
...     )
...     evolved = Statevector(initial_state).evolve(circuit.decompose(reps=6))
...     return abs(np.vdot(evolved.data, exact)) ** 2
>>>
>>> for steps in (1, 2, 4, 8):
...     print(f"{steps} step(s): fidelity {trotter_fidelity(steps):.6f}")
1 step(s): fidelity 0.053820
2 step(s): fidelity 0.704014
4 step(s): fidelity 0.932937
8 step(s): fidelity 0.982687

הנאמנות עולה מונוטונית לקראת אחד ככל שהצעדים מתקצרים, וזו ההתנהגות שנוסחת מכפלה חייבת להראות. קיצור הצעדים הוא רק אחד משני הכפתורים. FermionicSuzukiTrotter מחליף ביניהם בסדר גבוה יותר, כפי ש-סינתזה של שער Evolution במרחב פרמיוני מודד.

צעד יחיד של ההמילטוניאן הזה ב-$t = 1$ הוא קירוב גרוע (ארבע קבוצות הקפיצה והאינטראקציה לא מתחלפות), ולכן השורה הראשונה רחוקה מאחד, וזו הנקודה: ההתכנסות היא מה שמאשר את הבנייה, לא שורה בודדת.

הערה

כל גורם בלוח הזמנים פורש כל מוד של הסריג הזה, ולכן צמצום גורם לתמיכה שלו אינו עושה דבר כאן, והאינדקסים המוחלטים של הקיוביטים ש-FlowSetSynthesis מתאימה להם נשארים ללא שינוי. סינתזה שנכתבה עבור אופרטור שהקבוצות שלו נוגעות רק בחלק מהאוגר הייתה צריכה להתחשב בשינוי התיוג הזה.

מה כדאי לקחת מכאן

• שני ממדים משנים את הבִּיּוּק, לא את המתכון. כל מרכיב ממדריך ה-1D עובר ללא שינוי: groups לקיבוץ, map_transfer_vertex_generators() לקידוד, CustomF2QLayout לאוגר הלא תואם, ו-MapperFnEvolutionSynthesis כדי לחבר אותו ל-Transpiler. מה שגדל הוא מספר המקרים בתוך כל אחד מהם.

• זרות היא מה שהופך את ה-2D לעומק קבוע. קבוצת זרימה מתפרקת לקווים זרים בקודקודים, ולכן הקליפורדים שלה נפלטים במקביל וחולקים שתי שכבות שני קיוביטים לא משנה כמה גדול הסריג. ארבע קבוצות זרימה וקבוצה אלכסונית נותנות אז תקציב קבוע: כאן $16$ עבור ארבע קבוצות הזרימה ועוד $8$ עבור קבוצת האינטראקציה האלכסונית.

• לכל כיוון יש צורת פאולי משלו. $T_{kj} = -V_j V_k T_{jk}$ הופך XXZ ל-YYZ ו-XYXY ל-YXXY, ולכן ארבע קבוצות הזרימה צריכות ארבעה קליפורדים. להניח שקליפורד אחד מכסה את שני הכיוונים של ציר הוא הדרך הקלה ביותר לטעות בסינתזת קבוצות זרימה ב-2D.

• קראו תמיכות מהצמידה. ב-1D, brickwork אחד של CZ נותן את התמונות בצורה סגורה. ב-2D אין נוסחה כזו, ולקווי עמודה יש אנצילות לא רציפות, ולכן ניחוש שנגזר ביד לגבי המקום שבו כל איבר נוחת הוא שביר. צמידה ובדיקת התוצאה עולות בנייה אחת של Clifford ואי אפשר לטעות בהן. זה גם אומר שקליפורד תקף באותה מידה שמנחית את האיברים במקום אחר ממשיך לעבוד.

• תנו לקיבוץ למיין. מכיוון ש-Evolution מפצלת את עצמה לפי קבוצה, הסינתזה מקבלת קבוצת זרימה אחת בכל פעם ולא צריכה לגזור מחדש את החלוקה: היא קוראת צורת פאולי אחת מהאיבר הראשון ובוחרת את הקליפורד המתאים. ההנחה הזו, וזריקת שגיאה כשהיא לא מתקיימת, היא מה ששומר על synthesize קצרה.

• מספר הצעדים יושב מעל הפיצול. הסינתזה כאן רואה קבוצת זרימה אחת לכל שער, ולכן היא לא יכולה לחלק את ההתפתחות בין הקבוצות. FermionicLieTrotter עושה זאת במרחב פרמיוני, לפני הפירוק, וזה מה שמאפשר בכלל בדיקת התכנסות.

• אמתו על תת-המרחב. קידוד של $N$ מודים ב-$2N$ קיוביטים נאמן רק במקום שהמייצבים אומרים שהוא נאמן. מכיוון ש-VC שומר על $V_j = Z_j$, הבדיקה האמינה והזולה ביותר היא להטיל מצב מכפלה ולהשוות צפיפויות: בלי intertwiner ובלי התאמת ספקטרום, שאינה אמינה כשההמילטוניאן המוגבל מנוון.

הערה

ההמילטוניאן המקודד מתחלף עם המייצבים, ולכן צריך להחיל את ההטלה רק על המצב ההתחלתי. ההתפתחות שומרת על תת-המרחב מעצמה. על חומרה אמיתית זה גם מה שהופך את הקידוד למשאב שימושי לזיהוי שגיאות: אפשר למדוד את המייצבים.

הבא

שני כיוונים מרחיבים את מה שלמעלה, ושניהם שינויי מרכיב ולא שינויי שיטה. הראשון הוא הקידוד: Verstraete-Cirac נבחר כאן כי מקור [1] נותן את הקליפורדים המלכסנים שלו בצורה סגורה, אבל Derby-Klassen מוציא אנצילה אחת לכל פלקט אי-זוגי במקום אחת לכל אתר, ליחס קיוביטים-למודים של $3/2$ במקום $2$. לאופרטורי ההעברה שלו יש צורות פאולי שונות, ולכן העברה אליו מסתכמת בכתיבה מחדש של פונקציית הקידוד של שלב 3 ובמציאת הקליפורד לכל צורה חדשה. הקיבוץ, הפריסה ושלד הסינתזה עוברים ללא נגיעה.

השני הוא הסריג עצמו. כל מה שכאן מניח גבולות פתוחים ומין פרמיוני יחיד. גבולות מחזוריים סוגרים כל קו למעגל: החצים עדיין כולם מתחלפים, כי כל אחד פוגש את הבא ראש-אל-זנב, אבל שתי צורות הקליפורד של שלב 5 בנויות לקו עם קצוות. מעגל של $L$ אתרים נושא $L$ קשרים לכל כיוון בעוד שקו נושא $L - 1$, ולכן לקשר ההיקפי אין מקום בבנייה והוא צריך טיפול משלו. מודל עם ספין קל יותר: הוא מכפיל את האוגר ומוסיף איברי אינטראקציה באתר לקבוצה האלכסונית, בלי לגעת בקבוצות הזרימה. אף אחד מהם לא פוגע בטיעון העומק הקבוע, אבל שניהם משנים את הבִּיּוּק, שכפי שהמסקנה הראשונה אומרת, הוא החלק ששני ממדים מגדילים.

חשוב

קישוריות החומרה לא מטופלת כאן. כל מעגל למעלה בנוי על אוגר מופשט שבו כל זוג קיוביטים יכול לקיים אינטראקציה, והעומקים המדווחים נספרים במסגרת הזו. להתקנים אמיתיים יש מפת צימוד מוגבלת, ומיפוי הקידוד הזה עליה הוא בעיה משמעותית בפני עצמה: אוגר ה-VC הוא גרף של $2N$ קיוביטים ששעריו רצים בין אתר לאנצילה שלו, בין אתרים שכנים ובין אנצילות שכנות, ולכן פריסה צריכה לשמור על שלושת סוגי הזוגות סמוכים בו-זמנית. אם היא לא יכולה, ה-Transpiler מכניס שערי SWAP, והשכבות שהם מוסיפים הן מה שבניית העומק הקבוע קנתה בחזרה. מקור [1] מצטט גם הוא את נתון העומק שלו תחת גרף קישוריות ריבועי מונח ועם CX כמסבך היחיד המקורי. מציאת פריסות ואסטרטגיות ניתוב טובות לקידודים מקומיים על טופולוגיית התקן נתונה היא מעבר להיקף המדריך הזה.

[1] (1,2,3,4,5,6,7,8,9,10)

A. Gandon, S. Piccinelli, M. Rossmannek, F. Tacchino, A. Baiardi, J. Nys, and I. Tavernelli, Stabilizer-based quantum simulation of fermion dynamics with local qubit encodings (2026), arXiv:2512.11418v2.

[2] (1,2,3,4)

F. Verstraete and J. I. Cirac, Mapping local Hamiltonians of fermions to local Hamiltonians of spins, J. Stat. Mech. 2005, P09012 (2005).

[3] (1,2)

Y.-A. Chen and Y. Xu, Equivalence between fermion-to-qubit mappings in two spatial dimensions, PRX Quantum 4, 010326 (2023).