Solvable case of \(p=0\), generic initial coherent states

It is both instructive and interesting to first examine the evolution of the model in Eq. (1) without the periodic kicks (\(p=0\)). This is the well-known spin-glass model, the Sherrington-Kirkpatrick model82,83,84, with a nonzero average interaction strength. The thermodynamics of this model has long been of great interest. One can exactly solve the time evolution of disorder averaged \(J^2\) denoted as \(\langle J^2(t) \rangle _{w} = \langle J_x^{2} (t)+ J_y^{2}(t) + J_z^{2}(t) \rangle _w\), where \(J^2(t)=\langle \psi (t)|J^2|\psi (t) \rangle\). Here, \(|\psi (t)\rangle\) is the state at time t starting from an initial state \(|\theta , \phi \rangle\) in the PSS. Since \(J_x^2\) commutes with the \(p=0\) Hamiltonian,

$$\begin{aligned} \langle J_x^{2} (t) \rangle _w= \frac{N}{4} + \frac{N(N-1)}{4}\cos ^2{\phi }\sin ^2{\theta } \end{aligned}$$

(9)

The calculation of \(\langle J_y^{2}(t)+J_z^{2} (t)\rangle _w\) is more involved. However, a detailed analysis (see Supplementary materials for details) leads to:

$$\begin{aligned} \langle J_y^{2} (t)+ J_z^{2} (t) \rangle _w = \frac{2N}{4}+\frac{N(N-1)}{4}(1-\cos ^2 \phi \sin ^2 \theta ) \left[ \int _{-\infty }^{\infty } \frac{1}{\sqrt{2 \pi w^{2}}}\cos \left( \frac{kt}{N}\epsilon _{ij}\right) \exp \left( -\frac{\epsilon _{ij}^2}{2w^2}\right) d\epsilon _{ij} \right] ^{2(N-2)} \end{aligned}$$

$$\begin{aligned} = \frac{2N}{4}+\frac{N(N-1)}{4} (1-\cos ^2 \phi \sin ^2 \theta ) \exp \left[ -w^2 k^2 t^2 (N-2)/N^2\right] , \end{aligned}$$

(10)

and hence finally

$$\begin{aligned} \langle J^2 (t) \rangle _w = \frac{3N}{4} + \frac{N(N-1)}{4} \left[ \cos ^2 \phi \sin ^2 \theta \right. + \left. (1-\cos ^2 \phi \sin ^2 \theta ) \exp \left[ -w^2 k^2 t^2 (N-2)/N^2\right] \right] . \end{aligned}$$

(11)

Clearly, for a generic initial coherent state, \(\langle J^2 (t) \rangle _w\) decreases exponentially with time and saturates to a disorder independent, but initial state dependent value, which is \(O(N^2)\). The irreversible nature of the dynamics arises from the disorder averaging even in finite systems. When \(\theta =\pi /2\) and \(\phi =0\) or \(\pi\), there is no decay and \(\langle J^2 (t) \rangle _w=\langle J^2 (0) \rangle =j(j+1)\), where \(j=N/2\). These special initial states have no evolution except for a phase as they are eigenstates of the Hamiltonian. Interestingly, when \(\phi =\pm \pi /2\) (any value of \(\theta\)), \(\langle J^2(t) \rangle _w\) saturates to the RMT value given by 3N/4 for any non-zero value of disorder. These initial states are in the \(y-z\) plane of the Bloch sphere, and are maximally coherent (in the sense of quantum information theory, as discussed below) in the preferred x basis as this is the direction of the interaction. These are the only initial states for which the saturation value of \(J^2\) scales linearly as N, rather than quadratically, and are unstable in the sense that any deviation in the \(\phi\) value changes this from the RMT scaling. We will find that in the presence of the transverse field kicks when \(p \ne 0\), these behaviours are interestingly modified, while some persist. Below, we present our results in the presence of kicks, first for a generic coherent state and then for a non-generic, initial coherent states.

Fig. 3figure 3

Evolution of \(\left\langle J^{2}(n) \right\rangle _w\) with time n for the disordered system with kicks \((p=\pi /2)\) and without kicks \((p=0)\) for two different values of disorder strength \(w=1.5,~2.6\), initial state \({|{\theta ,\phi }\rangle }={|{2.25,1.1}\rangle }\), \(k=1\), and \(N=12\). The points correspond to numerical data and the dotted lines correspond to the analytical expression given in Eq. (11) for the unkicked Hamiltonian. The black horizontal lines is \(\langle J^2 _{RMT} \rangle =3N/4=9\) which is well below the saturation value. For \(p=0\), the saturation value clearly is independent of the strength of disorder, but depends on the initial state, as also predicted by Eq. (11). On the other hand, when \(p=\pi /2\), saturation value depends upon w and initial state.

Before doing that, it is noted that the kicked top unitary as defined in Eq. (2), in the absence of disorder, will be strictly periodic in k, with the period \(4\pi N\). Thus, the limits \(k=\infty\) and \(p=0\) are not equivalent, and the analytical scaling found for the latter case does not continue to hold for large k. In fact, we will restrict values of k to be much smaller than this periodicity. Also comparing the limiting cases for p and k is complicated by the Dirac delta kicks. In the usual kicked top without disorder, the relevant limit is \(N \rightarrow \infty\) with a fixed k and p.

Fig. 4figure 4

Left: \(\left\langle J^2(n) \right\rangle _w\) as a function of number of kicks n for different disorder strengths w with \(k=1\), \(p=\pi /2\), \(N=12\) qubits, averaged over 100 disorder realizations. Right: shows the von Neumann entanglement entropy \(\langle S_{N/2} \rangle _w\) as a function of n for the same parameters. The two vertical lines correspond to the two Heisenberg time scales discussed in the text, \(t_{PSS}^H\) (blue) and \(t_{FHS}^H\) (green). The black dashed horizontal line in both figures correspond to the RMT values in full Hilbert space. The initial state considered is \({|{2.25,1.1}\rangle }\).

Case of \(p=\pi /2\), generic initial states

We present the results of numerical calculations which are performed with an initial state \({|{\psi _0}\rangle }\) in the PSS as given in Eq. (5) with \(\theta =2.25,~\phi =1.1\) (\(\blacktriangle\) in Fig. 2) and is evolved to \({|{\psi _n}\rangle }=U_w^n {|{\psi _0}\rangle }\) using the Floquet operator \(U_{w}\) given in Eq. (2). This initial state is “generic” in the sense that it has no special dynamical significance in the permutation symmetric classical limit. It is still a coherent state and not generic in the sense of Haar measure. We use the fast Walsh Hadamard transform to reach system sizes of \(N=16\), and in one case to \(N=18\). Figure 3 compares the dynamics of the system, both, in the absence \((p=0)\) and in the presence of kicks \((p=\pi /2)\) for two different values of disorder strength w. In the absence of kicks, and for the generic initial coherent state with \(\cos \phi \ne 0\), the system reaches its saturation value when time \(t>t^* \sim \sqrt{N}/wk\) in the large N limit (see Eq. 11). The relevant time scale for saturation in the presence of both disorder and kicks is much more and is shown in Fig. 4 using both the average total angular momentum \(\langle J^2(n) \rangle _w\) and the entanglement entropy \(\langle S_{N/2} (n) \rangle _w\) for \(k=1\).

In contrast to the unkicked model, there is a disorder dependent steady state value eventually reached for all values of the disorder; however we find an interesting feature consisting of two regimes, first is till time \(n \sim t^H_{PSS}= N+1\) corresponding to the Heisenberg time relevant to the PSS, when there is a quasi steady state value. The other is at \(n \sim t^H_{FHS}= 2^N\), the Heisenberg time relevant to the FHS where \(\langle J^2(n) \rangle _w\) and \(\left\langle S_{N/2}(n) \right\rangle _w\) again starts saturating, which implies that the system has “realized” the presence of the full Hilbert space. Notice that with large disorder strength such as \(w=5\), when the full Hilbert space is being accessed, the time scale at \(t^H_{FHS}\) has become irrelevant for this quantity. Instead the Heisenberg time for the PSS provides a time scale for thermalization, when \(\langle J^2(n) \rangle _w\) and \(\left\langle S_{N/2}(n) \right\rangle _w\) saturates close to the RMT value corresponding to FHS. This is to be expected as the Heisenberg time for the PSS is also the Ehrenfest or log-time for the dynamics in the FHS. For intermediate values of the disorder strength w, the nonequilibrium state’s \(J^2\) as well as the entanglement saturates to non-thermal (non-RMT) values, signifying some kind of localization. It would be interesting to compare these with the well-studied case of MBL15,39,40.

Fig. 5figure 5

Left: Change in scaling of \(\overline{\langle J^2 \rangle }_w\) from \(N^{2}/4+N/2\) (red dashed line) to 3N/4 (black dotted line) as disorder strength w is increased with k set to 1 and \(p=\pi /2\). Right: Same for \(\overline{\langle S_{N/2} \rangle }_w\) as a function of N at \(k=1\) where the red dashed line corresponds to fitted line \(0.41+0.54 \log _2 (N/2+1)\) and the black dotted line is the RMT expression given by \(N/2- 1/2 \ln 2\). The initial state is \({|{2.25,1.1}\rangle }\).

To quantify the saturation value as a function of the disorder and system size, we now study the behavior of long time averaged values (averaged from \(n=10^5\) to \(n=3 \times 10^5\) for 100 disorder realizations) of \(J^2\) and \(S_{N/2}\), denoted as \(\overline{\langle J^2\rangle }_w\) and \(\overline{\langle S_{N/2} \rangle }_w\), respectively. From the discussions above, we expect that as \(w \rightarrow 0\) or small when the dynamics is predominantly within the permutation symmetric subspace, the scaling of long time limit of various quantities will follow the PSS scaling, i.e.,  \(\overline{\langle J^2\rangle }_w\sim N^2\) and \(\overline{\langle S_{N/2} \rangle }_w \sim \log _{2}(N/2+1)\). In the other extreme limit, when \(w \rightarrow \infty\) and the dynamics is described by RMT in the full \(2^N-\)dimensions, \(\overline{\langle J^2\rangle }_w=3N/4\) and \(\overline{\langle S_{N/2} \rangle }_w\sim N/2\). Figure 5 highlights the above discussed change in scaling when w is increased. Note that for the solvable case in the absence of kicks, while Eq. (11), indicates an earlier saturation time scale \(\sim \sqrt{N}\) for generic initial coherent states, the saturation value itself is still similar to the PSS case with \(\langle J^2 \rangle \sim N^2\).

It is tempting to associate this shift in the scaling of \(\overline{\langle J^2 \rangle }_w\) and \(\overline{\langle S_{N/2} \rangle }_w\) with a second order phase transition. Therefore, we now investigate the possible existence of a critical \(w_{c}\) below which the dynamics is predominantly within the PSS and above which the dynamics captures the properties of the FHS described by RMT, and obtain relevant critical exponents of the associated phase transition, if any. We propose a finite size scaling of the form38,85,86,87,88,89

$$\begin{aligned} {\overline{\langle J^2 \rangle }_w}=N^{\zeta /\nu }F((w-w_c)N^{1/\nu }), \end{aligned}$$

(12)

where \(\nu\) is the correlation length exponent. Figure 6a shows the crossing of the data around \(w_c=2.11\) when \(k=1\) whereas the collapse of the data with \(\nu =0.50\) and \(\zeta =0.57\) is shown in Fig. 6b. In order to confirm this phase transition further, we also study the variance of \(J^2\), var\((J^2)\)=\(\langle (\overline{J^2}-{\langle \overline{J^2} \rangle })^2 \rangle\) where \(\overline{J^2}\) is the time averaged value for each disorder. This quantity also shows a peak around \(w_c=2.0\) consistent with finite size scaling results, see Fig. 6c90,91,92. Next, we repeated the calculations for other values of k, the results of which are shown in Fig. 7.

Table 1 lists the critical disorder strength \(w_c\) along with the critical exponents \(\nu\) and \(\zeta\) for different values of the interaction strength or chaos parameter k. Intuitively, the chaotic dynamics at larger k values will require smaller disorder to take the dynamics to the full Hilbert space, which is also what we observe. Interestingly, we also find a weak parameter dependence on the critical exponent \(\nu\) for this phase transition. The exponent \(\nu\) lies within a range of 0.6 to 0.3 for the parameters studied. Such a variation seems to be a characteristic of a finite size disordered system such as the well studied many-body-localization40,93. We also find that in most cases, \(\nu \sim \zeta\). The quality of the data collapse for different values of k confirm the presence of a second order phase transition. At the same time, we do observe that identifying correct \(w_c\) values for larger k is difficult since \(w_c\) then is very small so that a single crossing point for all system sizes is difficult to obtain, and hence we restrict our studies for k values upto 2. Numerics show a decrease in the value of the critical strength as k is increased, indicating zero critical strength for large value of k. On the other hand when \(k \rightarrow 0\), the dynamics can get highly non-ergodic and these parameters need more careful examinations.

Table 1 Critical disorder strength \(w_c\) and exponents \(\nu\), \(\zeta\) for different chaos parameter, k. We have used \(\overline{\langle J^2 \rangle }_w\) with \(N=12,14\) and 16 to identify the critical points with time averaging ranging from \(n=10^5\) to 3 \(\times 10^5\) for 100 different realizations. The initial state is fixed to \(|\theta , \phi \rangle = |2.25, 1.1\rangle\).Fig. 6figure 6

(a) \(\overline{\langle J^{2} \rangle }_w/N^{\zeta /\nu }\) (where \(\zeta\) and \(\nu\) are the critical exponents) as a function of disorder strength w at \(k=1, p=\pi /2\) and initial state \({|{2.25,1.1}\rangle }\) for different system sizes. These curves cross each other at \(w=w_c \approx 2.11\). (b) The collapse of the data with \(\zeta =0.57\) and \(\nu =0.52\) consistent with the scaling proposed in Eq. (12). (c) var\(({J^{2}})\) with respect to disorder strength w for \(k=1\) that shows a peak around 2.0.

Fig. 7figure 7

(Left): The main figure shows the data of \(\overline{\langle J^{2} \rangle }_w/N^{\zeta /\nu }\) with respect to w corresponding to different system sizes for \(k=0.5, p=\pi /2\) and initial state \({|{\theta ,\phi }\rangle }={|{2.25,1.1}\rangle }\) crossing each other at \(w=w_c \sim 3.16\) whereas the inset shows the collapse of the data with \(\nu =0.58\) and \(\zeta =0.62\). (Right): Same as above but for \(k=2\) with \(w=w_c \sim 0.47\), \(\nu =0.37\) and \(\zeta =0.40\).

We have further confirmed this phase transition using entanglement entropy \(\overline{\langle S_{N/2} \rangle }_w\) with a similar finite size scaling as given in Eq. (12) and a different exponent \(\zeta ^{‘}\). The analysis gives \(\nu =0.60\), \(\zeta ^{‘}=0.68\) and \(w_c=1.93\) for \(k=1\), as shown in Fig. 8. We do acknowledge the fact that these numerics are limited by finite size systems and may not fully capture the behavior in the thermodynamic limit. This may also be the reason for slightly different values of \(w_c\) obtained using \(\overline{\langle S_{N/2} \rangle }_w\) and \(\overline{\langle J^2 \rangle }_w\).

Fig. 8figure 8

(a) Plot of \(\overline{\langle S_{N/2} \rangle }_w/N^{\zeta ‘/\nu }\) vs w for \(k=1, p=\pi /2\) and initial state \({|{2.25,1.1}\rangle }\) corresponding to different system sizes cross each other at \(w=w_c\sim 1.93 \pm 0.06\). Inset shows the collapse of the data with \(\nu =0.60\pm 0.06\) and \(\zeta ‘=0.68 \pm 0.07\).

Case of \(p=\pi /2\), special initial states

The above calculations were performed for a generic initial state specified by a \((\theta , \phi )\) which has no special significance. It is interesting to study the effect of changing this initial state on the nonequilibrium dynamics and on the critical behavior. For example, in the solvable case \(p=0\), we found that for all initial states in the \(y-z\) plane of the Bloch sphere, \(\langle J^2(n) \rangle _w\) eventually saturates to the RMT value, independent of the strength of the disorder. Figure 9 shows the scaling of the saturation value \(\langle J^2 \rangle\) for four different initial coherent states (see Fig. 2 for their location in phase space), including \((\pi /2,\pi /2)\), \((\pi /2,0)\), and two others \((\theta ,\phi )=(\pi /2, 1.9)\), \((\theta ,\phi )= (\pi /2, -0.56)\), which are somewhat close to \((\pi /2,\pi /2)\) and \((\pi /2,0)\), respectively. We observe a shift in \(w_c\) for certain initial states as compared to the one reported in Table 1.

Fig. 9figure 9

Finite size scaling for \(k=1, p=\pi /2\) and different initial conditions \(|\theta , \phi \rangle\): (a) \(|\pi /2, \pi /2\rangle\) shows \(w_c=1.2\), \(\nu =0.55\) and \(\zeta =0.58\), (b) \(|\pi /2,0\rangle\) with \(w_c=2.1\), \(\nu =0.47\) and \(\zeta =0.66\), (c) \(|\pi /2,1.9\rangle\), with \(w_c=1.10\pm 0.09\), \(\nu =0.54 \pm 0.01\) and \(\zeta =0.64 \pm 0.06\). This initial condition is closer to \((\pi /2,\pi /2)\), and hence has a smaller \(w_c\), (d) \(|\pi /2,-0.56\rangle\) with \(w_c=2.31\pm 0.05\), \(\nu =0.56\pm 0.03\) and \(\zeta =0.61 \pm 0.02\) is similar to Fig. 6, since both the initial states are closer to \((\pi /2,0)\).

To understand this initial state dependence behavior better, at least qualitatively, let us start from an initial state given by \((\theta , \phi )=(\pi /2,\pi /2)\), which is an eigenstate of \(\sigma _y\), so that the many body state from Eq. (5) is fully delocalized along the x-direction, in which the Hamiltonian has both the interaction and the disorder. In other words the quantum coherence, as measured by the off-diagonal matrix elements of the initial state is the largest possible in the \(\sigma ^x\) basis. Several recent studies have shown the importance of initial state coherence as a resource for the process of thermalization94,95,96. Thus these initial states with large coherence tend to delocalize preferentially. In the absence of disorder, the classical limit corresponding to these initial states, for \(k=1\), are stable fixed points, and it is interesting that the quantum effects of disorder tend to preferentially delocalize these states. This is consistent with a lower value of the critical strength \(w_c\) of the disorder. However, the crossing points in Fig. 9 are drifting most significantly in this case towards smaller values of \(w_c\) as N increases, and the scaling does not seem to work well. But with the system sizes we are able to go to, \(w_c\) is not larger than \(\approx 1.2\). Infact, our heuristic argument given in the next section discusses the peculiarity of this particular initial condition which may be the reason for not getting a good value of critical strength of disorder.

On the other hand when \((\theta ,\phi )=(\pi /2,0)\), the initial state is an eigenstate of \(\sigma _x\), with no off-diagonal matrix elements in \(\sigma ^x\) basis, or zero quantum coherence. In this case, stronger disorder along the x interaction direction is required to equilibriate in the full Hilbert space. The scaling behavior seems to be more convincing as well, and the critical disorder strength \(w_c\) is not larger than \(\approx 2.1\). Classically this corresponds still to a stable point, however, there is a small resonance island with nearby unstable points. Thus, the coherence of the initial states seems to play a crucial role when the dynamics is regular. We have numerically verified this argument using other values of \((\theta ,\phi )\) as well. The initial states with \((\theta ,\phi )\) closer to the \((\pi /2,0)\) has a larger \(w_c\) as compared to the ones closer to \((\pi /2, \pi /2)\), but with similar critical exponents, as shown in Fig. 9c,d. Both these figures show relatively good crossing and collapse for \(N=12, 14, 16\).

In summary, we performed a detailed analysis for \(k=1\) with different initial conditions and obtain the following critical disorder strength and the exponents by the scaling of \(J^2\) and \(S_{N/2}\) as shown in Table 2:

Table 2 Critical disorder strength \(w_c\) and exponents \(\nu\), \(\zeta\) for different initial conditions (\(\theta ,\phi\)) for the quantities, square of total angular momentum operator \(\left\langle J^2 \right\rangle\) and von Neumann entropy \(S_{N/2}\) at \(k=1\), \(p=\pi /2\).

Clearly, the values of exponent \(\nu\) are close to each other within the error bar, even for different initial conditions considered. The critical disorder strengths are also similar. The reason for smaller value of \(w_c\) for one of the initial conditions has been explained before.

Discussions

Given the small system sizes that we rely on in the finite size scaling analysis, a heuristic understanding may help. First, we notice that all the initial states that we take are in the permutation symmetric subspace and the expectation value of \(J^2\) scales as \(N^2\), which remains exactly true at all times in the absence of disorder. When there is no rotation (\(p=0\)), we have seen that the disorder decreases the \(J^2\) expectation value with time, but for almost all initial conditions the scaling remains \(N^2\) in the large time limit. The exceptions are initial states lying on the great circle, \(\phi =\pi /2\). Turning on the rotation about the y axis in presence of small disorder, a generic permutation symmetric initial state evolves essentially approximately to other symmetric states, and the long-time expectation value of \(J^2\) continues to scale as \(N^2\). This is due to the fact that the rotation is itself the same for all qubits involved, in other words it commutes with \(J^2\).

To make this more quantitative, notice that the Floquet operator in Eq. (2) can be written as \(U_w=U_{\epsilon } U_0\), where \(U_0\) is the kicked top without disorder and hence \([U_0,J^2]=0\). The disordering part is \(U_{\epsilon }\) which corresponds to a torsion about the \(x-\) axis that affects individual qubits differently. If \(w \ll 1\), \(U_w^n \approx U_{\epsilon }^n U_0^n\). Hence, if \({|{\psi _0}\rangle }\) is a permutation symmetric initial state, \(U_0^n {|{\psi _0}\rangle }\) remains in this subspace. From the analytic analysis with \(p=0\) we have seen that such states which will not generally be in the exceptional great circle or even be a simple coherent state, when acted upon by \(U_{\epsilon }^n\) will continue to have an expectation value of \(J^2\) scaling as \(N^2\). Thus, it is plausible that for small disorder values, the scaling of \(J^2\) in the large time limit remains \(N^2\) even when \(p=\pi /2\), pointing towards a non-zero \(w_c\). Results not shown here indicate that the time-evolving state’s overlap with the PSS is itself an order parameter that shows transitions similar to \(J^2\) and is consistent with the picture that before the transition, the state is largely in the PSS.

However, there is a qualitative difference if we start with a \(J_y\) eigenstate such that every qubit is in the coherent state \(\theta =\pi /2, \phi =\pi /2\). Here, the rotations about the y axis have no impact on this state and the rotations about the \(x-\)axis keeps the states on the \(y-z\) plane. We have shown that in the presence of disorder and \(p=0\) these are exceptional and that the \(J^2\) expectation value scales as N, even in the small w case. Hence we are obtaining values of \(w_c\) that seem to be flowing to a small value that may or may not be 0, see Fig. 9a. On the other hand the initial condition with all qubits starting from \(\theta =0\), or the collective \(J_z\) eigenstate has a larger value of \(w_c\), although it is also on the exceptional great circle \(\phi =\pi /2\). However, this initial state is strongly influenced by the rotation about the \(y-\) axis and is immediately removed from this circle, and hence it behaves differently from the \(J_y\) eigenstate.

Returning to the generic initial states, for w large and O(1), the time-evolving states quickly leave the PSS and are generic superpositions of the basis states in any basis, the above approximation breaks down and numerical results indicate the formation of random states. We have also confirmed the generation of random states by studying the spectral level spacing statistics which shows Wigner Dyson distribution at large disorder97. Indeed in this case the dynamics is similar to that of random circuits (for a recent review see98) that are known to generate psuedorandom quantum states that mimic true random states up to some order. Studying systems as in this paper and related transitions from the point of view of such circuits for which there is a considerable current interest is promising.

It may also be useful to compare with the large body of work around many-body-localization, wherein strong disorder drives eigenstates into localized phases15,38,39,40. In contrast to MBL, here we observe that the disorder breaks the permutation symmetry and takes the system to a chaotic dynamics in a larger Hilbert space dimension that leads to thermalisation of the system as opposed to localization in MBL. Although the canonical model of MBL, the Heisenberg 1D chain with a disordered onsite magnetic field39 is a time-independent system, MBL has been noted in Floquet systems as well99. In fact the finite size scaling effects that plague such studies and make the phenomenon itself suspect seems to be milder for Floquet systems100. The system studied here differs crucially with all-to-all interactions, scaling of the disorder with the number of particles, and the existence of a classical limit in the permutation symmetric, disorder free case. If there are further phase transitions in the present work, perhaps to a localized phase in the full Hilbert space remains to be seen. Despite the limitations of finite size systems, both MBL and the transitions observed in this work have clear relevance to the many-body systems being studied experimentally in this so-called NISQ (Noisy Intermediate-Scale Quantum) era.