## 2 Model
### Nearest-Neighbor Coupling Interactions in Superconductivity  
The coupling of qubits can be achieved through various types of interactions. For superconducting qubits, coupling can be achieved through interactions such as inductance, capacitance, or Josephson junctions. The interaction models can be classified into four categories: XY, Heisenberg, transverse Ising, and longitudinal Ising. These four types of Hamiltonians can be written together as:  
\[H=H_{0}+H_{I} \tag{3}\]  
\(H_{0}\) is the single qubit Hamiltonian,  
\[H_{0}=\sum\
olimits_{l=1}^{N}\frac{1}{2}\omega_{l}\sigma_{l}^{z} \tag{4}\]  
and \(H_{I}\) denotes the interaction Hamiltonian,  
\[H_{I}=\sum\
olimits_{l=1}^{N-1}\left(j_{l}^{x}\sigma_{l}^{x}\sigma_{l+1}^{x}+ J_{l}^{y}\sigma_{l}^{y}\sigma_{l+1}^{y}+j_{l}^{z}\sigma_{l}^{z}\sigma_{l+1}^{z}\right) \tag{5}\]  
Here, \(\sigma_{l}^{x},\sigma_{l}^{y},\sigma_{l}^{z}\) are the Pauli matrices in the basic of \(\sigma_{l}^{z}\) and \(\mathrm{l}\) denotes the \(l^{th}\) qubit. The coupling strength between the \(l^{th}\) and the \(l+1^{th}\) qubits is denoted by \(J_{l,l+1}\) (\(\mathrm{l}=1,\ldots,N-1\)). The parameters of the previous two equations are chosen to obtain the Hamiltonian of four different models and \(J_{l}\) is fixed accordingly. The different types of interaction hamiltonians are shown in Table 1. It can be seen that the pairing of models in Equation (2) does not follow the form of the Hamiltonians in Equation (3). Thus, it is imperative to design algorithms to simulate these pairing Hamiltonians using the four types of interaction Hamiltonians. For the operators \(\sigma_{m}^{x},\sigma_{m}^{y},\sigma_{m}^{z}\), care should be taken when simulating the Suzuki-Trotter decomposition, as the operators may not commute. It should be noted that the tunability of parameters \(\omega_{l}(l=1,\ldots,N)\) and \(j_{l}^{z}(i=x,y,z,l=1,\ldots,N-1)\) affect the efficiency of the algorithms. Figure 1: Pictorial representation of the BCS pairing theory. Given the initial state at time \(t=0\), it can be solved to obtain the state at any later instant [53]. If \(\hat{H}\) is independent of time, then  
\[|\psi(x,t)>\ =e^{-i\hat{H}t}|\psi(0)> \tag{6}\]  
where \(U(t)=e^{-i\hat{H}t}\) is the time evolution operator. \(U(t)\) is a unitary operator, i.e., \(U^{t}U=I\). For our case, the simulation of four different Hamiltonian was performed to find the corresponding unitary operations to build the quantum circuit. In general, to find the matrices for the time evolution operator, what the terms \(\hat{\sigma}^{x}\), \(\hat{\sigma}^{y}\), \(\hat{\sigma}^{z}\), \(\hat{\sigma}^{x}\hat{\sigma}^{x}\), \(\hat{\sigma}^{y}\hat{\sigma}^{y}\), and \(\hat{\sigma}^{z}\hat{\sigma}^{z}\) look like needs to be determined. Then, the terms can be separated for \(x,y,z\) as \(\hat{\sigma}^{x}\hat{\sigma}^{x}\), \(\hat{\sigma}^{y}\hat{\sigma}^{y}\), and \(\hat{\sigma}^{z}\hat{\sigma}^{z}\) to commute with each other. The Hamiltonian of our choice is  
\[H=\sum_{l=1}^{N}\frac{1}{2}\omega_{l}\sigma_{l}^{z}+\sum_{l=1}^{N-1}\left(j_ {l}^{x}\sigma_{l}^{x}\sigma_{l+1}^{x}+j_{l}^{y}\sigma_{l}^{y}\sigma_{l+1}^{y} +j_{l}^{z}\sigma_{l}^{z}\sigma_{l+1}^{z}\right) \tag{7}\]  
Figure 4: (**a**) Theoretical density matrix for the Heisenberg model, (**b**) experimental density matrix for 2 iterations, (**c**) experimental density matrix for 10 iterations, and (**d**) fidelity as a function of iteration. The corresponding time evolution operator is given by  
\[U(t)=e^{-i(\sum_{i=1}^{N}\frac{1}{2}\omega_{1}\sigma_{i}^{2})t}e^{-i\sum_{i=1}^{N -1}\left(l_{l}^{x}\sigma_{i}^{x}\sigma_{i+1}^{x}+\vec{l}_{l}^{y}\sigma_{i}^{y} \sigma_{i+1}^{y}+\vec{l}_{l}^{z}\sigma_{i}^{z}\sigma_{i+1}^{z}\right)t} \tag{8}\]  
Here, \(\vec{\sigma^{x}}= \begin{bmatrix}0&1\\ 1&0\end{bmatrix} \), \(\hat{\sigma^{y}}= \begin{bmatrix}0&-i\\ i&0\end{bmatrix} \), and \(\hat{\sigma^{z}}= \begin{bmatrix}0&-i\\ i&0\end{bmatrix} \). The first part of the operator can be written as per the value of Pauli matrix \(\hat{\sigma^{z}}\) as shown.