# IBM 2022 春季量子挑戰賽回顧 part2
文/林昱誠
### 第三題 Many-body quantum dynamics
這題我們要探討晶格無序（lattice disorder）和粒子-粒子交互作用。一個一開始處於非平衡態的封閉多體量子系統最終會達到平衡狀態，這現象我們稱之為 thermalization。這是統計力學的結果，就像是在沒有受外力影響情況下，桌上的熱咖啡最終會冷卻至跟周圍溫度相同。

然而，lattice disorder 的存在會阻止系統演化盡 egodic thermalized 狀態，這種 disorder 和粒子間的交互作用會導致 many-body localization（MBL），讓系統始終保持原本初始的狀態。欲了解更多有關 MBL，可以參看[此處](http://arxiv.org/abs1804.11065。)

#### **Box 1. Many-body localization 簡介**
如下圖，我們有一排的空格，在奇數格放上粒子，偶數格則空，即 $ψ_i=|101010..⟩$，定義粒子數算符為 $〈n ̂(x)〉=1,0,1,0…$ \（x 為格子座標），經過一段時間後粒子會在格子之間來回穿梭，我們可以預見最終會達成兩種平衡，一個如圖右上角，粒子均勻分配，$〈n ̂(x)〉=0.5,0.5,0.5,0.5…$，系統失去原本初態的所有訊息，這現象稱之為 thermalization

一般情況下，thermalization 後的狀態會滿足 Eigenstate thermalization hypothesis，因此這狀態也被稱為 ETH phase；第二種情況如圖右下角，粒子分佈呈現如 $〈n ̂(x)〉=0.9,0.1,0.9,0.1…$ （數字為舉例），跟一開始的狀態十分接近，仍保留原本初態的訊息，這現象我們稱之為 many-body localization (MBL)，物理學家也把 MBL看成一種新的物相，稱之為 MBL phase。


回顧前面兩題提到的 tight-binding Hamiltonian：

$$\frac{H_{tb}}{ℏ}=\sum_{i=0}^{3}(X_iX_{i+1}+Y_iY_{i+1})+\sum_{i}\epsilon_i Z_i$$

如同第二題，我們可以透過設定 $ϵ_i=W cos⁡(2πβ_i)$ 來模仿系統中的 disorder，其中 W 是 disorder 強度，β 是準晶體（介於晶體和非晶體之間的固體）週期因子。

阻礙 thermalization 的原因其中一個是出現不平衡現象。在 thermalized 系統中，我們預期系統達到穩定狀態時，晶格中每一個位置都被平均數量個粒子佔據（見 Box 1）。然而，如果給系統添加 disorder，會產生與 ETH phase 稍微有偏差的結果。一個系統的不平衡程度可以用下列公式做量化：

$$I=〈\frac{(N_e-N_o)}{(N_e+N_o )}〉$$

其中 $N_e$ 和 $N_o$ 代表系統基數位與偶數位的粒子數。對於 thermalized system（即 ETH phase），$I=0$，而 disorder 會對這個值產生偏差（$I≠0$）。
<br/>
:::danger
Challenge question 3a
完成可以產生該量子態的不平衡值的函數
:::

def get_imbalance(state):
    ###EDIT CODE BELOW
    ### HINT: MAKE SURE TO SKIP CALCULATING IMBALANCE OF THE |00...0> STATE
    imbalance_val=0
    
    ###DO NOT EDIT BELOW
    states_number = len(state.probabilities_dict())
    keys = []
    values = []
    states_imbalance_val = []
    
    for key, value in state.probabilities_dict().items():
        keys.append(key)
        values.append(value)
    
    for i in range(len(keys)):
        even = 0
        odd = 0
        for j in range(len(keys[i])):
            if j%2: #odd
                odd = odd + int(keys[i][j])
            else:
                even = even + int(keys[i][j])
            
        if (odd+even)== 0:
            state_imbalance = (odd-even)
            states_imbalance_val.append(state_imbalance)
        else:
            state_imbalance = (odd-even)/(odd+even)
            states_imbalance_val.append(state_imbalance)
            
    for k in range(len(states_imbalance_val)):
        imbalance_val = imbalance_val + states_imbalance_val[k] * values[k]
            
    return -imbalance_val

接著，我們要探討量子資訊和糾纏 entropy 如何在系統中增長。在糾纏系統中，忽略會使子系統進入不同數量態之混合態的其他系統，我們可以用 von Neuman entropy 來深入了解子系統 A 和其他系統之間的糾纏程度：

$$S_vn (ρ_A )=-tr(ρ_A \ln⁡ρ_A )$$

在這裡，$ρ_A$ 是子系統 A 的 reduced density matrix。如果子系統 A 跟其他系統是 fully entangled，則 $S_vn (ρ_A )=\ln⁡2$，如果子系統跟其他系統是完全隔開（即純態），則 $S_vn (ρ_A )=0$。讓我們看一個範例：

bell_state = qi.Statevector(np.array([0,1,1,0])/np.sqrt(2))

rho_0 = qi.partial_trace(bell_state,[1]) # We trace out qubit 1
rho_1 = qi.partial_trace(bell_state,[0]) # We trace out qubit 0

print('QB0 vn entropy: ', qi.entropy(rho_0, base=np.exp(1)))
print('QB1 vn entropy: ', qi.entropy(rho_1, base=np.exp(1)))

> QB0 vn entropy: 0.6931471805599454
> QB1 vn entropy: 0.6931471805599454

探討 12 個 qubits 之直鏈的多體量子動力學

t = Parameter(‘t’) 

num_qubits=12
deltas=[Parameter('delta_{:d}'.format(idx)) for idx in range(num_qubits)]
disorder_trot_step=generate_disordered_tb_instruction(t, deltas, num_qubits)

# Here we define the disorder pattern

beta=(np.sqrt(5)-1)/2 # DO NOT EDIT
AA_pattern=np.cos(2*np.pi*beta*np.arange(num_qubits)) # DO NOT EDIT

<br/>

:::danger
Challenge question 3b
將 qubits 0, 4, 和 8 設定為 |1⟩ 設定系統初始狀態
:::

delta_t=0.1
time_steps=np.arange(0,21,2)

circuits={}
Ws=[1,4,10]

for W in Ws:
    disorders=W*AA_pattern
    
    circuits[W]=[]

    for n_steps in time_steps:

        qr = QuantumRegister(num_qubits)
        qc = QuantumCircuit(qr)

        ###EDIT CODE BELOW
        qc.x(0)
        qc.x(4)
        qc.x(8)
        
        ###DO NOT EDIT BELOW     

        for _ in range(n_steps):
            qc.append(disorder_trot_step, [i for i in range(num_qubits)])
        
        if n_steps!=0:
            qc = qc.bind_parameters({t: delta_t})
            qc = qc.bind_parameters({deltas[idx]: disorders[idx] for idx in range(num_qubits)})

        circuits[W].append(qc)

<br/>

:::danger
Challenge question 3c
計算 qubit 0 於不同 disorder strengths與不同演化時間點的 von Neumann entropy
:::

<br/>

:::danger
Challenge question 3d
計算不同 disorder strengths 與不同演化時間點，晶格的不平衡值
:::

state_vector_imbalances={}
vn_entropies={}

for W in tqdm(Ws):
    probability_densities[W]=[]
    state_vector_imbalances[W]=[]
    vn_entropies[W]=[]
    
    for circ in circuits[W]:

        transpiled_circ=transpile(circ, backend_sim, optimization_level=3)

        job_sim = backend_sim.run(transpiled_circ)

        # Grab the results from the job.
        result_sim = job_sim.result()
        outputstate = result_sim.get_statevector(transpiled_circ, decimals=6)
        ps=[]
        for idx in range(num_qubits):
            ps.append(np.abs(qi.partial_trace(outputstate,[i for i in range(num_qubits) if i!=idx]))[1,1]**2)
        
        entropy=0
        ### EDIT CODE BELOW (extract the density matrix of qubit 0 by tracing out all other qubits)
        rho_3 = qi.partial_trace(outputstate,[i for i in range(1,num_qubits)])
    
        entropy = qi.entropy(rho_3, base=np.exp(1))
        ###DO NOT EDIT BELOW 
        
        imbalance=0
        ### EDIT CODE BELOW
        keys = []
        values = []
        states_imbalance_val = []
    
        for key, value in outputstate.probabilities_dict().items():
            keys.append(key)
            values.append(value)
    
        for i in range(len(keys)):
            even = 0
            odd = 0
            for j in range(len(keys[i])):
                if j%2: #odd
                    odd = odd + int(keys[i][j])
                else:
                    even = even + int(keys[i][j])
            
            if (odd+even)== 0:
                state_imbalance = (odd-even)
                states_imbalance_val.append(state_imbalance)
            else:
                state_imbalance = (odd-even)/(odd+even)
                states_imbalance_val.append(state_imbalance)
            
        for k in range(len(states_imbalance_val)):
            imbalance = imbalance + states_imbalance_val[k] * values[k]
        imbalance = -imbalance
        
        ###DO NOT EDIT BELOW
        
        vn_entropies[W].append(entropy)
        probability_densities[W].append(ps)
        state_vector_imbalances[W].append(imbalance)  

fig, axs = plt.subplots(1,3,figsize=(15,5), facecolor='white', sharey=True)

for i,W in enumerate(Ws):
    ax=axs[i]
    ax.pcolormesh(np.arange(0,num_qubits,1), time_steps*delta_t ,probability_densities[W])
    ax.set_xlabel('Qubit index')
    ax.set_xticks(np.arange(1,num_qubits+1,1))

axs[0].set_ylabel('Time (1/J)')

plt.show()

- **隨著時間，每一個 qubit 的機率密度**


- **隨著時間，不同 W 值的 entropy 值**


- **不同 disorder strength 之系統的不平衡值隨著時間的進展**