Stable Many-Body Resonances in Open Quantum Systems

Abstract

Periodically driven quantum many-body systems exhibit novel nonequilibrium states, such as prethermalization, discrete time crystals, and many-body localization. Recently, the general mechanism of fractional resonances has been proposed that leads to slowing the many-body dynamics in systems with both $U\left(1\right)$ and parity symmetry. Here, we show that fractional resonance is stable under local noise models. To corroborate our finding, we numerically study the dynamics of a small-scale Bose–Hubbard model that can readily be implemented in existing noisy intermediate-scale quantum (NISQ) devices. Our findings suggest a possible pathway toward a stable nonequilibrium state of matter, with potential applications of quantum memories for quantum information processing.

1. Introduction

Nonequilibrium quantum phases without static analogs have become an active area of research because programmable quantum simulators [1] such as cold atoms [2,3,4,5] and superconducting circuits [6,7,8,9,10,11] allow the preparation of exotic states of matter into a states out of equilibrium. Paradigmatic examples of nonequilibrium states are discrete time crystals [12,13,14,15,16,17,18,19,20,21] and Floquet prethermalization [22,23,24,25,26,27,28,29]. The search for nonequilibrium states of matter is challenging because one needs to move beyond standard quantum statistical mechanics [30,31,32,33,34], and simulations of quantum many-body systems in classical computers [35,36,37]. In a recent contribution by some of the authors [38], it has been pointed out the emergence of a prethermal and localized nonequilibrium phase, termed “fractional resonance”, appearing in a broad class of many-body Hamiltonians exhibiting U(1) and parity symmetry such as the Bose–Hubbard model [39,40], the $XXZ$ spin-1 model [41,42] or the Jaynes–Cummings–Hubbard model [43,44,45]. Moreover, small-scale versions of these models can be experimentally implemented in NISQ devices. Therefore, we believe it is pertinent to investigate for the stability of such nonequilibrium states of matter under the action of loss mechanisms, which will be relevant for applications in quantum memories [46] and quantum metrology [47].

On the other hand, quantum computing [48] and quantum machine learning [49,50,51] are progressing tremendously. The main driving force behind this progress is that quantum computers promise exponential speedup over their classical counterparts in solving specific tasks [52]. When fault-tolerant general-purpose quantum computers become available in the long term, it is expected they will implement adiabatic state preparation and quantum phase estimation as the standard quantum routine to determine the ground-state energy of a sophisticated physical Hamiltonian [53,54,55,56]. However, such schemes are resource-intensive, and thus they are not appropriate for current NISQ hardware [57,58,59,60,61], thereby shifting the central research theme toward low-depth hybrid quantum-classical algorithms, otherwise known as NISQ algorithms [62]. Recently, by utilizing quantum data processing based on the classical shadow, learning certain information about a many-body quantum state from an experiment can be exponentially reduced [63,64,65]. In contrast, in the standard classical paradigm, one must repeat the experiment many times to build statistical certainty about the measurement data at each experimental run. By loading the entire quantum state of a many-body system into a quantum memory, a shift in the paradigm occurs in a new approach. After some quantum processing is applied to the replicated state, measurements are performed. Then the protocol can be repeated. Furthermore, renowned quantum machine learning algorithms center around the Harrow–Hassidim–Lloyd (HHL) algorithm [66] and quantum principal component analysis [67], which require the physical realization of quantum random access memory [68,69]. Here, we recognize a substantial overlap with localized nonequilibrium states of matter, such as prethermal states [22,23,24,25,26,27,28,29] and many-body localization [70,71,72]. Moreover, because localized nonequilibrium states of matter retain information about their initial state, they may be used as quantum memories [73]. These and the above arguments indicate that developing and realizing a quantum memory represents the pinnacle achievement in quantum technologies and applications.

In this article, we provide a detailed analysis of the stability of fractional many-body resonances under noisy environments within the paradigm of NISQ devices. This investigation represents an essential step before seeking potential applications of fractional resonances and their associated prethermal states as a quantum memory. The article is organized as follows. In Section 2, we briefly recap the emergence of many-body resonances, emphasizing the integer and fractional resonance states using a generic Hamiltonian. In Section 3, we consider a one-dimensional lattice of strongly correlated bosonic particles described by the Bose–Hubbard model (BHM) [39,40], where we study the open quantum system dynamics of the BHM, and discuss the stability of integer and fractional resonance states when considering realistic parameters of NISQ devices such as superconducting circuits, for the sake of simulating the prescribed physics in an experiment. Finally, in Section 4, we present our concluding remarks.

2. The Model

Let us consider a one-dimensional lattice with open boundary conditions described by the generic lattice Hamiltonian. We have

$$\widehat{H}=\hslash \sum _{j=1}^L\left(\omega {\widehat{O}}_j+\frac{U}{2}{\widehat{O}}_j^2\right)-\hslash J_{0}cos\left(\mathsf{\Omega }t\right)\sum _{j=1}^{L-1}\left({\widehat{A}}_j^{†}{\widehat{A}}_{j+1}+{\widehat{A}}_{j+1}^{†}{\widehat{A}}_j\right),$$

(1)

which is composed of a local energy term ${\widehat{H}}_0\equiv \hslash {\sum }_{j=1}^L\left(\omega {\widehat{O}}_j+U{\widehat{O}}_j^2/2\right)$ with local operators ${\widehat{O}}_j$, and a time-dependent hopping interaction ${\widehat{H}}_I\left(t\right)=\hslash J_{0}cos\left(\mathsf{\Omega }t\right){\sum }_{j=1}^{L-1}\left({\widehat{A}}_j^{†}{\widehat{A}}_{j+1}+{\widehat{A}}_{j+1}^{†}{\widehat{A}}_j\right)$, allowing for the exchange of particles/excitations between the nearest neighboring sites via the local ladder operators ${\widehat{A}}_j$ and ${\widehat{A}}_j^{†}$. These operators satisfy the commutation relations $[{\widehat{O}}_i,{\widehat{A}}_j^{†}]={\delta }_{ij}{\widehat{A}}_j^{†}$, $[{\widehat{O}}_i,{\widehat{A}}_j]=-{\delta }_{ij}{\widehat{A}}_j$. Here, we restrict ourselves to studying bosonic particles. However, the treatment is also applicable to spin-d systems. If the many-body system is isolated from an external environment, the term ${\sum }_{j=1}^N\omega {\widehat{O}}_j$ is a constant of motion because the Hamiltonian exhibits U(1) symmetry. The Hermitian operator ${\widehat{O}}_j$ satisfies the eigenvalue equation ${\widehat{O}}_j|m_j\rangle =m_j|m_j\rangle$, where $m_j$ is a quantum number that labels the local quantum state of the jth lattice site. For instance, $m_j$ may represent the occupation number of a bosonic system or the spin component along the z direction of spin-d systems. In this work, we focus on finite-size lattices with reflection symmetry characterized by the parity operator $\widehat{\mathcal{P}}$, which satisfies $\widehat{\mathcal{P}}|m_1,\dots ,m_j,\dots ,m_L\rangle =|m_L,\dots ,m_j,\dots ,m_1\rangle$, where $|m_1,\dots ,m_j,\dots ,m_L\rangle \equiv {⨂}_{l=1}^L|m_l\rangle$.

To understand the processes that may occur due to the hopping of particles/excitations, let us move to a rotating frame with respect to the free Hamiltonian ${\widehat{H}}_0$. We obtain

$$\begin{array}{cc}\hfill {\tilde{H}}_I\left(t\right)=& -\hslash J_{0}cos\left(\mathsf{\Omega }t\right){\displaystyle \sum _{j=1}^{L-1}}[e^{iUt\left({\widehat{O}}_{j+1}-{\widehat{O}}_j-1\right)}{\widehat{A}}_j^{†}{\widehat{A}}_{j+1}\hfill \\ \hfill & +{\widehat{A}}_{j+1}^{†}{\widehat{A}}_j{e}^{-iUt\left({\widehat{O}}_{j+1}-{\widehat{O}}_j-1\right)}],\hfill \end{array}$$

(2)

from which we can identify two characteristic time scales: one is provided by the driving frequency $\mathsf{\Omega }$, and another is the onsite interaction U, which leads to a local anharmonic spectrum. The Hamiltonian (2) is not strictly periodic either in $\mathsf{\Omega }$ or U. However, we will consider integer and fractional driving frequencies in the unit of U defined as $\mathsf{\Omega }=U$ and $\mathsf{\Omega }=U/2$, respectively. In both cases, it can be shown that ${\tilde{H}}_I(t+T)={\tilde{H}}_I\left(t\right)$ is periodic with period $T=2\pi /\mathsf{\Omega }$ [38]. Moreover, we consider the strongly interacting regime, where $U\gg J_0$ [3], where $J_0$ is the bare (static) hopping rate. This regime allows us to describe the particle/excitation hopping processes within the semi-classical picture as discussed in Refs. [29,38].

Let us briefly describe the aforementioned semi-classical picture to identify many-body resonances [74] in the closed system scenario. To find many-body resonances, we focus on how specific quantum states, referred to as configurations, are connected via the hopping term. Let us suppose the many-body system is initialized in the configuration $|m_1,\dots ,m_j,m_k,m_l,\dots ,m_L\rangle$, a hopping event from the site j to k will lead the system to occupy the state $|m_1,\dots ,m_j-1,m_k+1,m_l,\dots ,m_L\rangle$. Now, let us compute the local energy of those states using the unperturbed Hamiltonian ${\widehat{H}}_0$. The result reads

$${\widehat{H}}_0|m_1,\dots ,m_j,m_k,m_l,\dots ,m_L\rangle =\hslash \sum _{j=1}^L\left(\omega m_j+\frac{U}{2}{m}_j^2\right)|m_1,\dots ,m_j,m_k,m_l,\dots ,m_L\rangle$$

(3a)

$${\widehat{H}}_0|m_1,\dots ,m_j-1,m_k+1,m_l,\dots ,m_L\rangle =\hslash \left[\sum _{j=1}^L\left(\omega m_j+\frac{U}{2}{m}_j^2\right)+U(m_k-m_j+1)\right]|m_1,\dots ,m_j-1,m_k+1,m_l,\dots ,m_L\rangle .$$

(3b)

The energy difference between these two configurations is $\Delta E=\hslash U(m_k-m_j+1)$. Therefore, to evolve from the state $|m_1,\dots ,m_j,m_k,m_l,\dots ,m_L\rangle$ to $|m_1,\dots ,m_j-1,m_k+1,m_l,\dots ,m_L\rangle$, the driving frequency should match the condition $m\mathsf{\Omega }=U(m_k-m_j+1)$. Analogously, it can be shown that a hopping event from $k\to j$ leads to the condition $m\mathsf{\Omega }=U(m_j-m_k+1)$, so we can write generically $m\mathsf{\Omega }=U[\pm (m_k-m_j)+1]$, with $m\in \mathbb{Z}$, which defines the integer many-body resonance [29,38]. In this case, the emerging many-body dynamics is ruled by nearest-neighbor interactions where the time scale for spreading particles/excitations is $J_0^{-1}$. In our recent work [38], the emergence of fractional resonances in many-body quantum systems has been proposed where second-order processes rule the nonequilibrium dynamics.

In this article, we have rigorously demonstrated the emergence of fractional resonance condition by using the Magnus expansion [75,76,77]. Here, we will follow the semi-classical picture described above to find the fractional resonance condition. If the initial configuration is $|m_1,\dots ,m_j,m_k,m_l,\dots ,m_L\rangle$, we want to connect this state with $|m_1,\dots ,m_j-1,m_k,m_l+1,\dots ,m_L\rangle$ via two hopping events. Notice that the latter configuration has local energy $\hslash {\sum }_{j=1}^L\left(\omega m_j+\frac{U}{2}{m}_j^2\right)+\hslash U(m_l-m_j+1)$. The energy difference between this configuration and the initial one is $\Delta E=\hslash U(m_l-m_j+1)$. Therefore, to connect both configurations, two identical hopping events are necessary, each kicking the system with energy $\hslash U(m_l-m_j+1)/2$. The same analysis can be done for a hopping from j→l leading to the condition $\mathsf{\Omega }=U(m_j-m_l+1)/2$, so we can write generically $m\mathsf{\Omega }=U[\pm (m_j-m_l)+1]/2$, where $m\in \mathbb{Z}$.

Let us consider the schematic representation of Figure 1 to illustrate integer and fractional resonance conditions better. Here, we display a schematic representation of many-body resonances in the driven lattice model. In Figure 1a, we show two configurations with all particles (green circles) occupying a lattice site and the leftmost site occupied by two particles. The energy difference between these configurations is $\Delta E=\hslash U$. Therefore, to evolve from the lower to the upper configuration, one needs to activate a hopping event with frequency $\mathsf{\Omega }=U$ to increase the system energy. In Figure 1b, we show configurations that involve two hopping events where the leftmost particle moves to the rightmost lattice site. The energy difference between these configurations is $\Delta E=\hslash U$, so each hopping event increases the system energy by a $\hslash U/2$.

Along this line of thinking, we will prove that the emergent many-body fractional resonance is a stable phenomenon under loss mechanisms inherent to NISQ devices. In particular, we will numerically study the nonequilibrium dynamics of the Bose–Hubbard model [39,40,78], which allows us to describe implementations of strongly interacting lattice systems in superconducting circuits [6,7,8,10,79,80,81].

The Bose–Hubbard Model

The Bose–Hubbard model (BHM) describes strongly interacting bosonic systems on a lattice, where we recognize operators ${\widehat{O}}_j={\widehat{n}}_j={\widehat{a}}_j^{†}{\widehat{a}}_j$, ${\widehat{a}}_j^{†}={\widehat{A}}_j^{†}$, and ${\widehat{a}}_j={\widehat{A}}_j$. Replacing these operators in the generic Hamiltonian (1), we can generate the Bose–Hubbard Hamiltonian

$${\widehat{H}}_{\mathrm{BH}}\left(t\right)=\hslash \sum _{j=1}^L\left(\omega {\widehat{a}}_j^{†}{\widehat{a}}_j+\frac{U}{2}{\widehat{a}}_j^{†}{\widehat{a}}_j^{†}{\widehat{a}}_j{\widehat{a}}_j\right)-\hslash J\left(t\right)\sum _{j=1}^{L-1}({\widehat{a}}_j^{†}{\widehat{a}}_{j+1}+\mathrm{H}.\mathrm{c}.),$$

(4)

where ${\widehat{a}}_j$(${\widehat{a}}_j^{†}$) is the annihilation (creation) bosonic operator at site jth, $\omega$ is the single site frequency, U the on-site interaction, $J\left(t\right)=J_{0}cos\left(\mathsf{\Omega }t\right)$ the modulated hopping strength between neighboring sites, whereas $J_0$ is the bare hopping rate, and $\mathsf{\Omega }$ is the driving frequency. We stress that modulating the hopping rate can be achieved in superconducting circuits with transmons [11]. The symmetries of the Hamiltonian (4) play a crucial role in describing the many-body dynamics. In particular, the BHM exhibits $U\left(1\right)$ symmetry characterized by the conservation of the total number of particles/excitations $[{\widehat{H}}_{\mathrm{BH}},\widehat{N}]=0$ with $\widehat{N}={\sum }_{j=1}^L{\widehat{a}}_j^{†}{\widehat{a}}_j$. Moreover, because we consider a lattice with open boundary conditions, the model preserves the parity $[{\widehat{H}}_{\mathrm{BH}},\widehat{P}]=0$, where $\widehat{P}|n_1,n_2,\dots ,n_L\rangle =|n_L,\dots ,n_2,n_1\rangle$, where $n_j$ stands for the number of particles/excitations at the jth lattice site. Hereafter, we will consider an initial state with unit filling, that is, $|{\psi }_0\rangle ={⨂}_{j=1}^L|1_j\rangle$. In general, the total number of states that may be involved in the dynamics correspond to all possible configurations of N particles distributed in L lattice sites $D_{N,L}=(N+L-1)!/N!(L-1)!$.

3. Open Quantum System Dynamics

Practical implementations of strongly correlated bosonic systems in superconducting circuits always involve interacting with some electromagnetic environment which leads to noisy dynamics. NISQ devices made of superconducting circuits have proven successful in stabilizing nonequilibrium many-body states [6,7,8,9,10,11,28], whose features are well captured by the BHM [39,40]. In these experiments, the open system dynamics is well described by the Lindblad master equation. However, because superconducting circuits such as transmons are anharmonic oscillators, we must use the microscopic master equation to fully capture the colored nature of the bath [82]. Considering a single anharmonic oscillator located at the jth lattice site, and described by the local Hamiltonian ${\widehat{H}}_S^{\left(j\right)}=\hslash {\sum }_n{\omega }_n^{\left(j\right)}|n_j\rangle \langle n_j|$, the master equation in the Lindblad form reads

$$\begin{array}{cc}\hfill \frac{d{\widehat{\rho }}_j}{dt}=& -\frac{i}{\hslash }[{\widehat{H}}_S^{\left(j\right)},{\widehat{\rho }}_j]\hfill \\ \hfill +& {\displaystyle \sum _{m,n>m}}{\kappa }_{nm}^{\left(j\right)}\mathbb{L}\left[\right|m_j\rangle \langle n_j\left|\right]{\widehat{\rho }}_j+{\displaystyle \sum _n}{\gamma }_{nn}^{\left(j\right)}\mathbb{L}\left[\right|n_j\rangle \langle n_j\left|\right]{\widehat{\rho }}_j,\hfill \end{array}$$

(5)

where $\mathbb{L}\left[\widehat{\mathcal{O}}\right]\widehat{\rho }=\widehat{\mathcal{O}}\widehat{\rho }\widehat{{\mathcal{O}}^{†}}-\frac{1}{2}({\widehat{\mathcal{O}}}^{†}\widehat{\mathcal{O}}\rho +\rho {\widehat{\mathcal{O}}}^{†}\widehat{\mathcal{O}})$. Here, ${\kappa }_{nm}^{\left(j\right)}$ and ${\gamma }_{nn}^{\left(j\right)}$ define decay and dephasing rates, respectively. In the decay rates, the subscripts $n,m$ refer to the decay from the state $|n_j\rangle$ to $|m_j\rangle$, whereas the dephasing rates refer to the dephasing for the superpositions of states $|m_j\rangle$ and $|n_j\rangle$, within each lattice site j. In this work, we consider finite lattices with $L=3$ and $L=4$ sites with up to $n_{\max }=3$ and $n_{\max }=4$ particles per site with local Hilbert space dimension $\dim \left({\mathcal{H}}_{\ell }\right)=4$ and $\dim \left({\mathcal{H}}_{\ell }\right)=5$, respectively. In these cases, a multilevel approach with local states $\left\{\right|0_j\rangle ,|1_j\rangle ,|2_j\rangle ,|3_j\rangle \}$ and $\left\{\right|0_j\rangle ,|1_j\rangle ,|2_j\rangle ,|3_j\rangle ,|4_j\rangle \}$ must be included in the Lindblad master Equation (5).

By using the multi-level approach for a single anharmonic oscillator, in this work, we consider the Bose–Hu+ter equation

$$\begin{array}{cc}\hfill \frac{d\widehat{\rho }}{dt}& =-\frac{i}{\hslash }[{\widehat{H}}_{\mathrm{BH}}\left(t\right),\widehat{\rho }]+{\displaystyle \sum _{j=1}^L\sum _{m=0,n>m}^{n_{\max }}}{\kappa }_{nm}^{\left(j\right)}\mathbb{L}\left[\right|m_j\rangle \langle n_j\left|\right]\widehat{\rho }\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^L\sum _{n=0}^{n_{\max }}}{\gamma }_{nn}^{\left(j\right)}\mathbb{L}\left[\right|n_j\rangle \langle n_j\left|\right]\widehat{\rho },\hfill \end{array}$$

(6)

where $\widehat{\rho }$ represents the density matrix of the strongly interacting lattice.

3.1. Three-Site Bose–Hubbard Lattice

Let us consider a three-site Bose–Hubbard lattice. We will present results about the stability of the emerging many-body dynamics when considering integer (${\mathsf{\Omega }}_1=U$) and fractional (${\mathsf{\Omega }}_2=U/2$) resonance conditions in an open system scenario. Our study involves the numerical solution of the Lindblad master Equation (6) by using the fourth-order Runge–Kutta algorithm. We consider a product state with unit filling $|{\psi }_0\rangle =|111\rangle$ as an initial condition. In a closed system scenario, the system will only populate states within the positive parity subspace $|{\psi }_0\rangle$, $|{\psi }_1\rangle =\frac{1}{\sqrt{2}}\left(\right|120+021\rangle )$, $|{\psi }_2\rangle =\frac{1}{\sqrt{2}}\left(\right|102+201\rangle )$, $|{\psi }_3\rangle =\frac{1}{\sqrt{2}}\left(\right|210+012\rangle )$, $|{\psi }_4\rangle =|030\rangle$, $|{\psi }_5\rangle =\frac{1}{\sqrt{2}}\left(\right|300+003\rangle )$ [38]. In an open quantum system scenario, we expect that U(1) symmetry will no longer be preserved. Our model will consider identical loss mechanisms for each lattice site, implying parity symmetry is still held.

In Figure 2 and Figure 3, we plot the populations of states $|{\psi }_i\rangle$ with $i=0,1,2,3,4,5$ for the integer resonance condition ${\mathsf{\Omega }}_1=U$ and fractional resonance condition ${\mathsf{\Omega }}_2=U/2$, respectively. We identify populations as $P_i\left(t\right)={\left|\langle {\psi }_i|\psi \left(t\right)\rangle \right|}^2$. In both figures, $P_i$ are populations numerically computed considering a closed system scenario governed by the Hamiltonian (4), and $P_i^D$ are populations numerically obtained considering decay and dephasing mechanisms acting upon each lattice site via the master Equation (6). We stress that the populations of states $|{\psi }_1\rangle$ and $|{\psi }_2\rangle$ are the same, so we only show $P_1\left(t\right)$.

In our numerical simulations, we use realistic values for decay (${\kappa }_{nm}$) and dephasing (${\gamma }_{nn}$) rates of superconducting circuit experiments [83]. In addition, single site frequency $\omega$ and hopping rate $J_0$ are taken from Ref. [10], and $U=40J_0$ (see the caption of Figure 2 and Figure 3 for details). It is worth mentioning that despite the long time needed for the population of state $|{\psi }_3\rangle$ to occur with the highest probability in the fractional resonance case, see Figure 3b, the system is robust under loss mechanisms. By using realistic parameters, we estimate driving periods of about $T_1\approx 2$ ns and $T_2\approx 4$ ns, for integer and fractional resonance conditions, respectively. Although relaxation and dephasing affect the populations within the positive parity subspace, strong oscillations of probabilities still survive within our simulating time. This is a signature of the stability of many-body resonances. Strong oscillations still survive in local observables, such as the average occupation number per site. Figure 4 shows the average occupation number of each lattice site $N_j\left(t\right)$ for the integer resonance ${\mathsf{\Omega }}_1=U$. The upper panel shows the dynamic evolution of the closed system, and the lower panel shows the dynamics in the open system scenario. As expected, $N_1\left(t\right)=N_3\left(t\right)$ due to the reflection symmetry of the lattice with identical parameters for loss mechanisms. Within the simulating time $t=8\mu$s ($t=4000T_1)$, we still see strong oscillations of the occupation numbers and Poincare recurrences. Because we consider zero-temperature baths for relaxation, we expect $N_j\left(t\right)\to 0$ in the long-time dynamics. We stress that due to the memory consumption of the Runge–Kutta algorithm, we cannot simulate longer times by using the realistic parameters of Ref. [83].

Figure 5 shows the average occupation number of each lattice site $N_j\left(t\right)$ for the fractional resonance ${\mathsf{\Omega }}_2=U/2$. The upper panel shows the dynamic evolution of the closed system, and the lower panel shows the dynamics in the open system scenario. Although the average occupation numbers start decaying, within the simulating time $t=4\mu$s ($t=1000T_2$), strong oscillations around the initial value $N_j\left(0\right)=1$ are still present. As in the integer resonance, we expect $N_j\left(t\right)\to 0$ in the long-time dynamics because zero-temperature baths act upon individual bosonic particles. Notice that $N_1\left(t\right)$ and $N_3\left(t\right)$ decay faster in comparison to $N_2\left(t\right)$. This happens because second-order processes dominate the dynamics, and lattice edges can be populated with more than one particle, which results in faster relaxation processes.

When considering noisy dynamics, the U(1) symmetry is broken, and the system could populate states outside the unit-filling subspace. In this case, the wave function may be written as a linear combination of all possible configurations $\left|\psi \right(t)\rangle ={\sum }_{l=1}^{\mathcal{M}}{c}_l\left(t\right)|l\rangle$, with $\mathcal{M}={\sum }_{N=0}^L{D}_{N,L}$ for a fixed number of sites L. Therefore, it is necessary to compute the probability of all accessible configurations the system may visit along with its dynamical evolution. In the three-site Bose–Hubbard lattice, there are $\mathcal{M}={\sum }_{N=0}^3{D}_{N,3}=64$ accessible configurations. Figure 6 shows the distribution of populations of each configuration $|c_l{\left(t\right)|}^2$ for the integer case ${\mathsf{\Omega }}_1=U$. The upper panel shows the system dynamics considering a closed system, whereas the lower panel is an open system scenario. In the long-time dynamics, a fraction of configurations outside the unit filling subspace start to be populated; however, signatures of the stability of the integer resonance are visible in the average occupation number shown in Figure 4. In analogy, Figure 7 shows the distribution of populations of each configuration $|c_l{\left(t\right)|}^2$ for the fractional case ${\mathsf{\Omega }}_2=U/2$. The upper panel shows the system dynamics considering a closed system, whereas the lower panel is an open system scenario.

3.2. Four-Site Bose-Hubbard Lattice

Let us consider the four-site Bose–Hubbard lattice initialized in the state $|{\psi }_0\rangle =|1111\rangle$. In this case, there are $\mathcal{M}={\sum }_{N=0}^4{D}_{N,4}=625$ accessible configurations. As in the three-site Bose–Hubbard lattice, in Figure 8, we plot the distribution of populations of each configuration $|c_l{\left(t\right)|}^2$ for the fractional case ${\mathsf{\Omega }}_2=U/2$. The upper panel shows the system dynamics considering a closed system, whereas the lower panel is an open system scenario. In both cases, we see that the system populates the same configurations. Therefore, we conclude that, for finite lattice sites, and for short times up to $t=20T_2$ (t = 80 $\mathsf{\mu }$s), the slowing-down characteristic of the fractional resonance is a stable phenomenon under noisy dynamics. We emphasize that due to resource consumption of the Runge–Kutta algorithm and our limited computational resources, and our current computational resources, we can only simulate up to $t=20T_2$ for a lattice of $L=4$ sites. In the short term, we expect to extend our study to larger lattice sizes $L>10$ using the adaptive time-dependent density matrix renormalization group approach to confirm further the findings presented here. Notice that there is strong evidence about the stability of many-body resonances when increasing the lattice size, at least for the closed system dynamics, see Ref. [38] and Appendix A.

3.3. The Linear Entropy

Another way to characterize the stability of the system under noisy mechanisms is by using the linear entropy defined as $S\left(\widehat{\rho }\right)=1-\mathrm{tr}({\widehat{\rho }}^2)$, where $\widehat{\rho }$ is the system density matrix. In a closed system scenario, the linear entropy $S\left(\widehat{\rho }\right)=0$ at all times because $\widehat{\rho }$ is a pure state. However, a realistic situation will necessarily imply the density matrix to be a statistical mixture. In the previous section, we demonstrated that the fractional many-body resonance and its characteristic slowing down of the dynamics [38] is a stable phenomenon under loss mechanisms. The latter is also reflected in linear entropy as shown in Figure 9, where the upper (lower) panel shows the linear entropy as a function of time for a lattice of $L=3$ ($L=4$) sites. In both cases, $S\left(\widehat{\rho }\right)$ is is less than 1 within the simulating time.

4. Conclusions

Considering small lattices of strongly interacting bosonic particles, we have provided robust evidence about the stability of many-body resonances and their characteristic slowing down under realistic parameters of NISQ devices implemented in superconducting circuits. Our investigation includes decay and dephasing mechanisms acting locally on each bosonic particle using the Lindblad master equation, which has been proven helpful in describing state-of-the-art superconducting circuit experiments. In the short term, we expect to extend our study to larger lattice sizes $L>10$ using the adaptive time-dependent density matrix renormalization group approach to confirm further the findings presented here. Our conclusions are essential before seeking potential fractional resonance applications and their associated prethermal states as a quantum memory device.