Multi-qubit Quantum Rabi Model and Multi-partite Entangled States in a Circuit QED System

TL;DR

This paper proposes a method to simulate the multi-qubit quantum Rabi model in circuit QED systems using external driving fields, enabling the generation of entangled states and exploring strong coupling regimes.

Contribution

A theoretical scheme utilizing external drives to realize the multi-qubit quantum Rabi model and generate entangled states in circuit QED systems.

Findings

Effective Hamiltonian derived for multi-qubit Rabi model

Ability to generate high-fidelity entangled states

Simulation shows robustness across parameter ranges

Abstract

Multi-qubit quantum Rabi model, which is a fundamental model describing light-matter interaction, plays an important role in various physical systems. In this paper, we propose a theoretical method to simulate multi-qubit quantum Rabi model in a circuit quantum electrodynamics system. By means of external transversal and longitudinal driving fields, an effective Hamiltonian describing the multi-qubit quantum Rabi model is derived. The effective frequency of the resonator and the effective splitting of the qubits depend on the external driving fields. By adjusting the frequencies and the amplitudes of the driving fields, the stronger coupling regimes could be reached. The numerical simulation shows that our proposal works well in a wide range of parameter space. Moreover, our scheme can be utilized to generate two-qubit gate, Schr\"odinger states, and multi-qubit GHZ states. The maximum displacement of the Schr\"odinger cat states can be enhanced by increasing the number of the qubits and the relative coupling strength. It should be mention that we can obtain high fidelity Schr\"odinger cat states and multi-qubit GHZ states even the system suffering dissipation. The presented proposal may open a way to study the stronger coupling regimes whose coupling strength is far away from ultrastrong coupling regimes.

Tables

Table 1. Table 1: The realistic parameters in the numerical simulation are listed in the following table [ 67 , 68 ] .

$\epsilon /2\pi$	${\omega }_r/2\pi$	$g/2\pi$	${\Omega }_x/2\pi$	${\omega }_x/2\pi$	$\gamma /2\pi$	$\kappa /2\pi$
10 GHz	10 GHz	20 MHz	2 GHz	9.98 GHz	0.05 MHz	0.012 MHz

Table 2. Table 2: The fidelity of the GHZ states at time t = 2 ​ π / ω ~ r 𝑡 2 𝜋 subscript ~ 𝜔 𝑟 t=2\pi/\tilde{\omega}_{r} is list in the following table.

	$N=2$	$N=3$	$N=4$	$N=5$	$N=6$
$F_{\mathrm{full}}$	0.9971	0.9924	0.9886	0.9816	0.9761
$F_{\mathrm{diss}}$	0.9806	0.9644	0.9497	0.9321	0.9162

Equations

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\mathrm{int}}+{\hat{H}}_{\mathrm{d}},$$

$${\hat{H}}_0$$

$${\hat{H}}_{\mathrm{int}}$$

$${\hat{H}}_{\mathrm{d}}$$

$$U(t)=\exp \left(-i\frac{{\Omega }_x}{2}{\hat{J}}_{x}t\right)\exp \left(-i{\omega }_x{\hat{J}}_{z}t-i{\omega }_x{\hat{a}}^{†}\hat{a}t\right)$$

$${\Omega }_x=2{\omega }_z,{\omega }_x\gg {\Omega }_x\gg g,{\omega }_z\gg {\Omega }_z,{\omega }_z\gg \epsilon -{\omega }_x,$$

$${\hat{H}}_{\mathrm{eff}}={\tilde{\omega }}_r{\hat{a}}^{†}\hat{a}+\tilde{\epsilon }{\hat{J}}_z+\tilde{g}(\hat{a}+{\hat{a}}^{†}){\hat{J}}_x,$$

$$\frac{\tilde{g}}{{\tilde{\omega }}_r}=\frac{g}{2({\omega }_r-{\omega }_x)},$$

$${\hat{H}}_{\mathrm{eff}}^I(t)=\frac{\tilde{g}}{2}\left({\hat{J}}_{+}{e}^{i\tilde{\epsilon }t}+{\hat{J}}_{-}{e}^{-i\tilde{\epsilon }t}\right)\left(\hat{a}{e}^{-i{\tilde{\omega }}_{r}t}+{\hat{a}}^{†}{e}^{i{\tilde{\omega }}_{r}t}\right).$$

$${\hat{H}}_{\mathrm{eff}}^I(t)=\tilde{g}\left(\hat{a}{e}^{-i{\tilde{\omega }}_{r}t}+{\hat{a}}^{†}{e}^{i{\tilde{\omega }}_{r}t}\right){\hat{J}}_x.$$

$$\mathcal{U}(t)=\exp \left({\Omega }_1(t)+{\Omega }_2(t)\right),$$

$${\Omega }_1(t)$$

$${\Omega }_2(t)$$

$$\mathcal{U}(t)=D\left(\beta (t){\hat{J}}_x\right)\exp \left(i\varphi (t){\hat{J}}_x^2\right),$$

$$\dot{\rho }=-i\left[\hat{H}(t),\rho \right]+\gamma \sum _{k=1}^N\mathcal{L}[{\hat{\sigma }}_k^{-}]\rho (t)+\kappa \mathcal{L}[\hat{a}]\rho (t),$$

$$e_p(\mathcal{U})={\left(\frac{d}{d+1}\right)}^2\left[E(\mathcal{U})+E(\mathcal{U}{S}_{12})-E(S_{12})\right],$$

$$\mathrm{CNOT}=(u_1\otimes u_2)\mathcal{U}(u_3\otimes u_4),$$

$$\begin{array}[]{llll}u_{1}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&-1\\ -1&-1\\ \end{array}\right),&u_{2}=\left(\begin{array}[]{cc}1&0\\ 0&1\\ \end{array}\right),&u_{3}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&i\\ -1&i\\ \end{array}\right),&u_{4}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}1&i\\ i&1\\ \end{array}\right).\end{array}$$

$$\begin{array}{r}{\mathrm{F}}_{\mathrm{pro}}(\mathcal{E},\mathcal{U})=\frac{1}{d^3}\sum _{j=1}^{d^2}\mathrm{Tr}\left[\mathcal{U}{W}_j^{†}{\mathcal{U}}^{†}\mathcal{E}\left(W_j\right)\right]\end{array}.$$

$$\begin{array}{r}|\psi (t)⟩=\frac{1}{\sqrt{2}}\exp \left(i{N}^2\varphi (t)/4\right)\left({|\frac{N}{2},\frac{N}{2}⟩}_x\otimes {|\frac{N}{2}\beta (t)⟩}_r+{|\frac{N}{2},-\frac{N}{2}⟩}_x\otimes {|-\frac{N}{2}\beta (t)⟩}_r\right),\end{array}$$

$$|\psi (t)⟩$$

$$\mathcal{U}(T)=\exp \left(i\varphi (T){\hat{J}}_x^2\right).$$

$${|\frac{N}{2},-\frac{N}{2}⟩}_z$$

$${|\frac{N}{2},\frac{N}{2}⟩}_z$$

$$|\psi (T)⟩=\mathcal{U}(T){|\frac{N}{2},-\frac{N}{2}⟩}_z=\sum _{M=-N/2}^{N/2}{C}_M{e}^{i\varphi (T)M^2}{|\frac{N}{2},M⟩}_x.$$

$$|\psi (T)⟩=\mathcal{U}(T){|\frac{N}{2},-\frac{N}{2}⟩}_z=\sum _{M=-N/2}^{N/2}{C}_M{e}^{i\frac{M^2}{2}\pi }{|\frac{N}{2},M⟩}_x.$$

$${|\psi (T)⟩}_e=\sum _{M=-N/2}^{N/2}\frac{C_M}{\sqrt{2}}\left(e^{i\pi /4}+(-1{)}^{-M}{e}^{-i\pi /4}\right){|\frac{N}{2},M⟩}_x.$$

$${|\psi (T)⟩}_e$$

$${|\psi (T)⟩}_o$$

$${|{\psi }^{\prime }(T)⟩}_o$$

$${\hat{H}}^{\prime }(t)$$

$${\hat{H}}^{\prime }(t)\approx ({\omega }_r-{\omega }_x){\hat{a}}^{†}\hat{a}+[\epsilon -{\omega }_x+{\Omega }_z\cos ({\omega }_{z}t)]{\hat{J}}_z+\frac{{\Omega }_x}{2}{\hat{J}}_x+\frac{g}{2}({\hat{a}}^{†}{\hat{J}}_{-}+\hat{a}{\hat{J}}_{+}),$$

Full text

Introduction

The quantum Rabi model (QRM) [1, 2, 3], which is a fundamental model to describe the interactions between light and matter, occupies a crucial position in various physical systems, such as quantum optics [4], solid state system [5], molecular system [6], and so on. When the ratio between the coupling strength ($g$) and the mode frequency ($\omega$) satisfies $g/\omega\ll 1$, the rotating-wave approximation (RWA) is suitable, and the counter-rotating term (CRT) can be ignored. In this case, the QRM is reduced to the Jaynes-Cummings (JC) model [7], which has been applied to explain many physical phenomena, such as the revivals of the atomic population inversion after its collapse [8, 9], vacuum Rabi splitting [10, 11], and so on. Recently, new coupling regimes, such as ultra-strong coupling (USC) and deep-strong coupling (DSC) regimes, have been reached in some circuit QED systems [12, 13, 14, 15, 16, 17]. In this case, the CRT cannot be neglected. Consequently, many interesting effects induced by CRT appear in these regimes [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. The implementations of QRM in USC and DSC regimes have also motivated new applications to the quantum information processing [30, 31, 32, 33, 34, 35]. It should be mentioned that, though great progress has been achieved, it is also challenging to implement QRM in USC, and DSC regimes experimentally.

Of particular interest is how to simulate the QRM in USC, and even DSC regimes when the system is far from USC regime. Motivated by this consideration, some quantum simulation approaches have been proposed in various physical systems, such as superconducting circuits [36, 37, 38, 39, 40], quantum optical [41], trapped ions [42, 43], cold atoms [44, 45, 46], and so on. These quantum simulation proposals provide us with experimental feasible methods to implement QRM in USC and DSC regimes. Very recently, the quantum simulations of USC and DSC regimes extend to the multi-qubit case with trapped ions and anisotropic quantum Rabi model with superconducting circuits. The simulations of the generalized models provide us with platforms to study concerning physical issues, such as quantum critical phenomena, multi-partite entanglement, and so on.

On the other hand, the qubit-dependent displacement interaction describes a quantum resonator conditionally displaced according to qubit(s)’ states. Such interaction plays an important role in understanding the fundamentals of quantum physics [47, 48, 49, 50, 51, 52, 53]. Based on such type interaction, the superposition of the coherent states can be prepared in various of systems [48, 49, 50, 51, 52, 53]. The qubit-dependent displacement interaction also has been used to the quantum information processing [54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65], such as generation of unconventional phase gate [60] and multipartite entangled states [61, 62, 63]. Following the theoretical and experimental study of the QRM[38], we focus on the simulation of multi-qubit QRM, and we will study its applications to generation of two-qubit quantum gate, Schrödinger cat states and multi-qubit GHZ states.

In this paper, we propose an alternative scheme to simulate multi-qubit QRM in USC regime, and even DSC regime with a circuit QED setup. The system consists of multiple flux qubits, which strongly coupled to a resonator. To obtain the tunable multi-qubit QRM, we apply transversal and longitudinal external driving fields on the qubits. We show the stronger coupling regimes can be reached by tuning the driving amplitudes and frequencies. Additionally, we study some applications of simulated Hamiltonian on two-qubit quantum gate, superposition coherent states and multi-qubit entangled states. The results show that the non-trivial two-qubit gate is equivalent to the controlled-NOT (CNOT) gate. Based on the multi-qubit conditional interaction Hamiltonian, the Schrödinger cat states and multi-qubit GHZ states can be generated. The maximum displacement of the Schrödinger cat states depends on the number of qubits and the relative coupling strength, which indicates the maximum displacement can be enhanced by increasing number of the qubits and the relative coupling strength.

The derivation of the effective Hamiltonian

In this section, we first derive a effective QRM, in which the relative coupling strength can be adjusted by tuning the frequency of the external driving fields. We also show the fidelity of the simulated Hamiltonian. We consider $N$ qubits strongly coupled to a single-mode harmonic oscillator. The qubits are driven by the longitudinal and transversal external driving fields. Such model can be realized in a variety of physical systems. Here we adopt a circuit QED setup to demonstrate our proposal. We consider $N$ flux qubits are coupled to a transmission line resonator, which can be modeled as a single mode harmonic oscillator. Assuming the qubits are tuned to the degeneracy point, then the Hamiltonian in this case reads (here and after, we set $\hbar=1$)

[TABLE]

where

[TABLE]

Here the operator $\hat{a}$ ($\hat{a}^{{\dagger}}$) is the annihilation (creation) operator of the bosonic field with frequency $\omega_{r}$. The qubits are described by Pauli matrices $\hat{\sigma}^{\alpha}_{k}$ ($\alpha=x,y,z$), which denotes $\alpha$ component of the $k$-th Pauli matrix. For simplicity, we consider all the qubits possess the same energy splitting $\varepsilon$ (i.e., $\varepsilon_{k}=\varepsilon$), and the qubits couple to the bosonic field with unified coupling strength $g$ (i.e. $g_{k}=g$). $\hat{H}_{\rm int}$ shows the interaction between the resonator and the qubits. All the qubits are driven by two classical fields with the frequencies $\omega_{z}$ and $\omega_{x}$, and the corresponding amplitudes are denoted by $\Omega_{z}$ and $\Omega_{x}$. In this case, we introduce the collective operators $\hat{J}_{\alpha}=\frac{1}{2}\sum_{k=1}^{N}\hat{\sigma}_{\alpha}$ to simplify the Hamiltonian (1) as $\hat{H}_{0}=\omega_{r}\hat{a}^{{\dagger}}\hat{a}+\varepsilon\hat{J}_{z}$, $\hat{H}_{\rm int}=g(\hat{a}^{{\dagger}}+\hat{a})\hat{J}_{x}$, and $\hat{H}_{\rm d}=\Omega_{z}\cos(\omega_{z}t)\hat{J}_{z}+\Omega_{x}\cos(\omega_{x}t)\hat{J}_{x}$. Choosing the rotating framework defined by

[TABLE]

and considering the following conditions

[TABLE]

we can neglect the fast oscillating terms and obtain the following time-independent effective Hamiltonian (the detailed derivation is shown in the Methods section)

[TABLE]

where $\tilde{\omega}_{r}=\omega_{r}-\omega_{x}$, $\tilde{\varepsilon}=\Omega_{z}/2$, and $\tilde{g}=g/2$ are the effective energy splitting of the qubits, effective frequency of the resonator, and effective coupling strength, respectively. Such effective Hamiltonian describes a multi-qubit generalization of quantum Rabi model (i.e., Dicke model), in which the frequency of the resonator and the energy splitting of qubits can be adjusted by tuning the frequencies and amplitudes of external driving fields. The relative coupling strength reads

[TABLE]

The relative coupling strength can be adjusted by tuning the frequency of the transversal driving fields. Thus we can obtain the multi-qubit QRM in different coupling regimes.

In order to assess the validity of the effective Hamiltonian. We compare the time-dependent evolution states governed by the full Hamiltonian (1) and the effective Hamiltonian (5). Let $\left|\psi(0)\right\rangle=\left|gg\cdots g\right\rangle\otimes\left|0\right\rangle_{r}$ be the initial state and the evolution states governed by the Hamiltonian (1) and (5) are denoted by $\left|\psi(t)\right\rangle$ and $\left|\tilde{\psi}(t)\right\rangle_{\rm ideal}$, respectively. We denote the evolution state governed by Hamiltonian (1) in the rotating framework defined by $U(t)$ with $\left|\tilde{\psi}(t)\right\rangle=U^{{\dagger}}(t)\left|\psi(t)\right\rangle$. The fidelity of the evolution states $\left|\tilde{\psi}(t)\right\rangle$ and $\left|\tilde{\psi}(t)\right\rangle_{\rm ideal}$ reads ${\rm F}(t)=\left|\left\langle\tilde{\psi}(t)|\tilde{\psi}(t)\right\rangle_{\rm ideal}\right|^{2}$. Considering the approximate conditions, we choose the following parameters: $\varepsilon=\omega_{r}$, $\Omega_{z}=0.004\omega_{r}$, $\Omega_{x}=2\omega_{z}=0.2\omega_{r}$, $g=0.002\omega_{r}$ and $\omega_{x}=\left\{0.996,0.998,0.999,0.9995\right\}\omega_{r}$. Under such parameters, the relative coupling strength are $\tilde{g}/\tilde{\omega}_{r}=\left\{0.25,0.5,1,2\right\}$ and the system is driven to stronger coupling regimes. In Fig. 1, we plot the fidelity of evolution states for $N=2$ (black solid line), $N=3$ (blue dash-dotted line), $N=4$ (red dashed line), $N=5$ (green dotted line) and $N=6$ (cyan solid line). The Fig. 1(a)-1(d) show the fidelity when the relative coupling strength $\tilde{g}/\tilde{\omega}_{r}=\left\{0.25,0.5,1,2\right\}$, respectively. The results show that the effective Hamiltonian is validity when the number of the qubits and the relative coupling strength are not very large.

The applications of the effective Hamiltonian

In this section, we will illustrate some applications to the simulated multi-qubit QRM on quantum information processing. Such as the generation of quantum gate, the Schrödinger cat states, and multi-qubit GHZ states. Moving to the rotating frame associated with $U^{\prime}_{3}=\exp\left(-i\tilde{\varepsilon}\hat{J}_{z}t-i\tilde{\omega}_{r}\hat{a}^{{\dagger}}\hat{a}t\right)$, the effective Hamiltonian is recast as following form

[TABLE]

If we consider all the qubits have zero effective energy splitting (i.e., $\tilde{\varepsilon}=0$), the Eq. (7) can be reduced to the following form

[TABLE]

This is a periodic Hamiltonian with period $T=2\pi/|\tilde{\omega}_{r}|$. The evolution operator for Hamiltonian (8) can be obtained by means of the Magnus expansion [66]

[TABLE]

where

[TABLE]

Considering the commutator $[\Omega_{1}(t),\Omega_{2}(t)]=0$, the evolution operator can be recast as follows

[TABLE]

where the displacement operator is given by $D(\beta)=\exp\left(\beta\hat{a}^{{\dagger}}-\beta^{*}\hat{a}\right)$. The parameters $\beta(t)$ and $\phi(t)$ are defined as $\beta(t)=(\tilde{g}/\tilde{\omega}_{r})(1-e^{i\tilde{\omega}_{r}t})$ and $\phi(t)=(\tilde{g}/\tilde{\omega_{r}})^{2}(\tilde{\omega}_{r}t-\sin\tilde{\omega}_{r}t)$. For the following convenience, we introduce the collective states, which is the eigenstates of the collective operators $\left\{\hat{J}^{2},\hat{J}_{\alpha}\right\}$. Let the collective states $\left\{\left|j,j_{\alpha}\right\rangle_{\alpha};\leavevmode\nobreak\ j_{\alpha}=-j,-j+1,\cdots j;\leavevmode\nobreak\ \alpha=x,y,z\right\}$ be the eigenstates of operator set $\{\hat{J}^{2},\hat{J}_{\alpha}\}$, and they satisfy the following equations: $\hat{J}^{2}\left|j,j_{\alpha}\right\rangle_{\alpha}=j(j+1)\left|j,j_{\alpha}\right\rangle_{\alpha}$, $\hat{J}_{\alpha}\left|j,j_{\alpha}\right\rangle_{\alpha}=j_{\alpha}\left|j,j_{\alpha}\right\rangle_{\alpha}$.

In the following, we will use the evolution operator given in Eq. (11) to generate two-qubit quantum gate, Schrödinger cat state, and $N-$qubit GHZ states. To describe the dynamics of the system under dissipation, we utilize the following master equation

[TABLE]

where $\rho(t)$ is the time-dependent density matrix. The time-dependent density matrix in the rotating framework can be obtained by $\tilde{\rho}=U(t)\rho U^{{\dagger}}(t)$ and its dynamics is governed by the Hamiltonian $\hat{\tilde{H}}(t)=U^{{\dagger}}H(t)U(t)-iU^{{\dagger}}(t)\partial_{t}U(t)$, which is full Hamiltonian in the rotating framework. The qubits decay rate and resonator loss rate are denoted by $\gamma$ and $\kappa$, respectively. $\mathcal{L}[\hat{A}]\rho=\frac{1}{2}\left(2\hat{A}\rho\hat{A}^{{\dagger}}-\hat{A}^{{\dagger}}\hat{A}\rho-\rho\hat{A}^{{\dagger}}\hat{A}\right)$ is the Lindblad superoperator describing the losses of the system. In the following numerical simulation, we adopt the following realistic parameters [67, 68]: $\varepsilon=\omega_{r}=2\pi\times 10$ GHz, $\Omega_{x}=2\omega_{z}=2\pi\times 2$ GHz, $g=2\pi\times 20$ MHz and $\omega_{x}=2\pi\times 9.98$ GHz. The decay rate of the qibit and resonator loss rate are taken as $\gamma=2\pi\times 0.05$ MHz and $\kappa=2\pi\times 0.012$ MHz. We switch off the longitudinal driving fields (i.e., $\Omega_{z}=0$). The parameters are list in Table (1). Under such parameters, the relative coupling strength is $\tilde{g}/\tilde{\omega}_{r}=0.5$ and the effective energy splitting is $\tilde{\varepsilon}=0$.

The realization of the quantum gate.

To obtain the two-qubit quantum gate, we consider $N=2$ and evolution time $t=T=2\pi/\tilde{\omega}_{r}$. In this case, $\hat{J}_{x}=(\hat{\sigma}_{1}^{x}+\hat{\sigma}_{2}^{x})/2$ and the evolution operator (11) reduces to $\mathcal{U}(T)=\cos(\phi/2)\mathcal{I}+\sin(\phi/2)\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{x}$ with $\phi=2\pi(\tilde{g}/\tilde{\omega_{r}})^{2}$, where $\mathcal{I}$ is the identity operator for two-qubit system. Here, we have omitted a global phase. Obviously, such quantum gate is capable to generate entanglement when $\phi\neq m\pi$ ($m$ is an integer). To describe the entanglement generation capacity of the unitary operator, we utilize the entangling power given by Zanardi et.al. [69, 70, 71, 72]. The entangling power defined on $d\times d$ system can be expressed in terms of the linear entropy of operators $\mathcal{U}$, $\mathcal{U}S_{12}$, and $S_{12}$ as follows

[TABLE]

where $S_{12}=\sum_{i,j}|ij\rangle\langle ji|$ is the swapping operator acting on the tensor product space and the linear entropy of the $\mathcal{U}$ is given by $E(\mathcal{U})=1-\frac{1}{d^{4}}Tr[\mathcal{U}^{R}(\mathcal{U}^{R})^{{\dagger}}\mathcal{U}^{R}(\mathcal{U}^{R})^{{\dagger}}]$. The rearrangement of $\mathcal{U}$ is defined as $\mathcal{U}^{R}_{ij,kl}=\mathcal{U}_{ik,jl}$ [72]. The entangling power of the quantum gate can be obtained as $e_{p}(\mathcal{U})=\frac{2}{9}\sin^{2}\phi$. When $\phi=\pi/2$ (i.e. $\tilde{g}/\tilde{\omega_{r}}=0.5$), we obtain a quantum gate with maximum quantum entangling power. Such non-trivial quantum gate is local equivalent to the CNOT gate [73, 74]. We can check the following local equivalence relation

[TABLE]

where the local unitary operators are as follows

[TABLE]

In order to assess the performance of our proposal to generate CNOT equivalent gate against sources of error, we adopt the process fidelity ${\rm F}_{\rm pro}$, which measures the difference between ideal and real quantum processes. For an ideal unitary process $\mathcal{U}$ and its real process $\mathcal{E}(\mathcal{U})$, the process fidelity reads

[TABLE]

For two-qubit system, $d=4$ and $W_{j}$ is the operator basis acting on the 4-dimensional Hilbert space. The operator basis can be represented with the Pauli matrices (i.e., $W_{j}\in\left\{\mathcal{I},\hat{\sigma}_{1}^{x},\cdots,\hat{\sigma}_{1}^{z}\hat{\sigma}_{2}^{z}\right\}$). If we adopt the full Hamiltonian without dissipation (i.e., Eq. (1)) under the parameters listed in Table (1), the process fidelity can reach $99.57\%$. If we adopt the full Hamiltonian with dissipation (i.e., Eq. (12)), the process fidelity of the quantum gate is $96.32\%$. The higher performance of the quantum gate needs to resort to adopt superconducting qubit with lower decay rate.

The generation of Schrödinger cat states.

The conditional interaction is also crucial in creating superposed coherent states and hence exploring the superposition rule. We then investigate the creation of the superposed coherent states with the derived effective Hamiltonian. We choose $|\psi(0)\rangle=\frac{1}{\sqrt{2}}\left(\left|\frac{N}{2},\frac{N}{2}\right\rangle_{x}+\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{x}\right)\otimes\left|0\right\rangle_{r}$ as initial state. Acting the evolution operator in Eq. (11) on the initial state, we obtain

[TABLE]

where the coherent states $|\pm\frac{N}{2}\beta(t)\rangle_{r}=\hat{D}\left[\pm\frac{N}{2}\beta(t)\right]|0\rangle_{r}$ with the coherent state amplitude $\pm\frac{N}{2}\beta(t)=\pm\frac{N}{2}(\tilde{g}/\tilde{\omega}_{r})(1-e^{i\tilde{\omega}_{r}t})$. Obviously, the spin states $\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{x}$ and $\left|\frac{N}{2},\frac{N}{2}\right\rangle_{x}$ undergo different dynamics, which depends on the spin states. Let us introduce the states $\left|\pm\right\rangle=\frac{1}{\sqrt{2}}\left(\left|\frac{N}{2},\frac{N}{2}\right\rangle_{x}\pm\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{x}\right)$, then the evolution state can be rewritten as

[TABLE]

where $\mathcal{N}_{\pm}=\left[2(1\pm\exp(-\frac{1}{2}N^{2}|\beta(t)|^{2}))\right]^{-1/2}$ and $\left|{\rm Cat}_{\pm}\right\rangle=\mathcal{N}_{\pm}\left(\left|\frac{N}{2}\beta(t)\right\rangle_{r}\pm\left|-\frac{N}{2}\beta(t)\right\rangle_{r}\right)$. The superposition coherent states $\left|{\rm Cat}_{\pm}\right\rangle$ are the so-called even and odd Schrödinger cat states [49, 50, 51, 52, 53]. After measurement is performed on the states $\left|+\right\rangle$ and $\left|-\right\rangle$, the final state in Eq. (18) collapses to the states $\left|{\rm Cat}_{+}(t)\right\rangle$ or $\left|{\rm Cat}_{-}(t)\right\rangle$. The probability of obtaining even and odd cat states are $\frac{1}{2}\left[1+\exp{(-\frac{1}{2}N^{2}|\beta(t)|^{2})}\right]$ and $\frac{1}{2}\left[1-\exp{(-\frac{1}{2}N^{2}|\beta(t)|^{2})}\right]$, respectively. The magnitude of the displacement $\frac{N}{2}|\beta(t)|$ changes depending on the evolution time. When $t_{0}=\pi/\tilde{\omega}_{r}$, the displacement reaches its maximum value $N\tilde{g}/\tilde{\omega}_{r}$, which indicates the maximum displacement can be enhanced by increasing number of the qubits $N$ and the relative coupling strength $\tilde{g}/\tilde{\omega}_{r}$.

In order to study the Schrödinger states generation when the system subjects to dissipation, we compare the evolution states under effective Hamiltonian with quantum states governed by full Hamiltonian with and without dissipation. Let $\left|\psi(t_{0})\right\rangle$ be target state. We denote time-dependent density matrix governed by the effective Hamiltonian, full Hamiltonian without dissipation and master equation with $\tilde{\rho}_{\rm ideal}(t)$, $\tilde{\rho}_{\rm full}(t)$ and $\tilde{\rho}_{\rm diss}(t)$, respectively. We compare expected state $\left|\psi(t_{0})\right\rangle$ with evolution states by using the fidelities ${\rm F}_{\rm ideal}=\left\langle\psi(t_{0})\right|\tilde{\rho}_{\rm ideal}(t)\left|\psi(t_{0})\right\rangle$, ${\rm F}_{\rm full}=\left\langle\psi(t_{0})\right|\tilde{\rho}_{\rm full}(t)\left|\psi(t_{0})\right\rangle$ and ${\rm F}_{\rm diss}=\left\langle\psi(t_{0})\right|\tilde{\rho}_{\rm diss}(t)\left|\psi(t_{0})\right\rangle$. The Fig. 2 shows the numerical results for ${\rm F}_{\rm ideal}$ (black dotted line), ${\rm F}_{\rm full}$ (red dash-dotted line) and ${\rm F}_{\rm diss}$ (blue solid line). The results show that when evolution time $t=\pi/\tilde{\omega}_{r}$, the target state is reached. Even when the system subjects to dissipation, we also can obtain Schrödinger cat states when the number of the qubits is not very large.

The generation of multi-qubit GHZ states.

The derived effective Hamiltonian in Eq. (5) also can be used to generate the multi-qubit GHZ states [62, 63]. Let $\tilde{\varepsilon}=0$ and the evolution time $t=T=2\pi/|\tilde{\omega}_{r}|$, we get $\beta(T)=0$ and $\phi(T)=2\pi\left(\tilde{g}/\tilde{\omega}_{r}\right)^{2}$. Then the evolution operator (11) reduces to

[TABLE]

The multi-qubit states $\left|gg\cdots g\right\rangle$ and $\left|ee\cdots e\right\rangle$ can be recast in terms of the collective states as $\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{z}$ and $\left|\frac{N}{2},\frac{N}{2}\right\rangle_{z}$, respectively. The collective states $\left|\frac{N}{2},\pm\frac{N}{2}\right\rangle_{z}$ can be expressed in terms of eigenstates of the $\{\hat{J}^{2},\hat{J}_{x}\}$ as follows

[TABLE]

Let $\left|\psi(0)\right\rangle=\left|gg\cdots g\right\rangle\equiv\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{z}$ be the initial state. Acting the unitary operator (19) on the initial state, we obtain

[TABLE]

If we set $\phi(T)=\pi/2$ (i.e., $\tilde{g}/\tilde{\omega}_{r}=1/2$), the above final state reads

[TABLE]

In the following, we proof the above final state is local equivalent to the $N-$qubit GHZ state. When $N$ is an even integer, $M$ are integers ranging from $-\frac{N}{2}$ to $\frac{N}{2}$. We also check that $e^{iM^{2}\pi/2}$ is equal to $1$ for even $M$ and $i$ for odd $M$. Then the final state in this case reads

[TABLE]

Such state can be expressed as superposition of the collective states $\left|\frac{N}{2},\frac{N}{2}\right\rangle_{z}$ and $\left|\frac{N}{2},-\frac{N}{2}\right\rangle_{z}$ as follows

[TABLE]

If $N$ is an odd integer, $M$ are half integers. We can introduce an integer $M^{\prime}$ with $M^{\prime}=M-\frac{1}{2}$. Then the final state (22) are

[TABLE]

where $\left|\psi^{\prime}(T)\right\rangle_{o}=\sum_{M^{\prime}=-(N+1)/2}^{(N-1)/2}C_{M}e^{i\frac{M^{\prime 2}}{2}\pi}\left|\frac{N}{2},M\right\rangle_{x}$. The state $\left|\psi^{\prime}(T)\right\rangle_{o}$ is local equivalent to the state $\left|\psi(T)\right\rangle_{o}$. Considering even or odd integer $M^{\prime}$, the state $\left|\psi^{\prime}(T)\right\rangle_{o}$ can be rewritten as

[TABLE]

Based on the Eqs. (24) and (26), the final state is equivalent to the GHZ state for even or odd integer $N$. The above results apply to an ideal situation, namely, dissipation-free environment. To assess the experimental feasibility of our proposal, we compare multi-qubit GHZ states $\left|{\rm GHZ}\right\rangle^{(N)}_{o/e}$ (we denote $\left|{\rm GHZ}\right\rangle$ for simplicity) with evolution states governed by the effective Hamiltonian (i.e., Eq. (5)), the full Hamiltonian without dissipation (i.e., Eq. (1)), the full Hamiltonian with dissipation (i.e., Eq. (12)). We denote the evolution density matrices governed by Eq. (5), Eq. (1) and Eq. (12) with $\tilde{\rho}_{\rm ideal}(t)$, $\tilde{\rho}_{\rm full}(t)$ and $\tilde{\rho}_{\rm diss}(t)$, respectively. The fidelity between multi-qubit GHZ states $\left|{\rm GHZ}\right\rangle$ and evolution states are denoted by ${\rm F}_{\rm ideal}=\left\langle\rm GHZ\right|\tilde{\rho}_{\rm ideal}(t)\left|{\rm GHZ}\right\rangle$, ${\rm F}_{\rm full}=\left\langle\rm GHZ\right|\tilde{\rho}_{\rm full}(t)\left|{\rm GHZ}\right\rangle$ and ${\rm F}_{\rm diss}=\left\langle\rm GHZ\right|\tilde{\rho}_{\rm diss}(t)\left|{\rm GHZ}\right\rangle$. The Fig. 3 shows the numerical results for ${\rm F}_{\rm ideal}$ (black dotted line), ${\rm F}_{\rm full}$ (red dash-dotted line) and ${\rm F}_{\rm diss}$ (blue solid line). The Fig. 3 (a)-3(e) are the numerical results for $N=2,3,4,5,6$, respectively. The fidelity for ${\rm F}_{\rm full}$ and ${\rm F}_{\rm diss}$ at time $t=2\pi/\tilde{\omega}_{r}$ are shown in Table (2). The results show that we can obtain high fidelity multi-qubit GHZ state even the system subjecting to dissipation.

Discussion

In summery, we have proposed a scheme to simulate the multi-qubit quantum Rabi model in circuit QED system. The effective Hamiltonian for multi-qubit quantum Rabi model can be derived. Based on unitary dynamics, the fidelity of effective Hamiltonian is discussed in detail. The results show that the system can reach stronger coupling regimes by adjusting the external driving amplitudes and frequencies. With this tunable effective Hamiltonian, the qubit-dependent displacement interaction Hamiltonian can be obtained by tuning the driving parameters. Based on such Hamiltonian, we also discuss the applications to constructing nontrivial quantum gate, the Schrödinger cat states and multi-qubit GHZ states. With the effective Hamiltonian, we can generate the quantum gate with the maximum two-qubit entangling power. The local equivalence between the achieved quantum gate and the CNOT gate has been discussed in detail. The numerical calculation shows that the process fidelity of the quantum gate reaches $96.32\%$ under the chosen parameters. The Schrödinger cat states can be generated with the effective Hamiltonian, and the magnitude of the displacement can be enhanced by increasing the number of the qubits and relative coupling strength. In the case of multiple quantum qubits, we generate high fidelity multi-qubit GHZ states for even and odd $N$. We show that the high fidelity Schrödinger cat state and multi-qubit GHZ state can be obtained even the system subjecting to dissipation.

The presented proposal may open a way to study the stronger coupling regimes whose coupling strength is far away from ultrastrong coupling regimes. We should note that the effective Hamiltonian is not validity when the number of the qubits and the relative coupling strength are very large. Even so, our scheme may also provide potential applications to the quantum computation and quantum state engineering.

Methods

In this part, we will show how to obtain the effective Hamiltonian in Eq. (5). We choose the rotating framework related to the time-dependent unitary transformation $U_{1}(t)=e^{-i\omega_{x}\hat{J}_{z}t}e^{-i\omega_{x}\hat{a}^{{\dagger}}\hat{a}t}$. The transformed Hamiltonian reads

[TABLE]

where $\hat{J}_{\pm}=\hat{J}_{x}\pm i\hat{J}_{y}$ and ${\rm H.c.}$ denotes the Hermitian conjugate. Considering the condition $\omega_{x}\gg\Omega_{x}$ and $\omega_{x}\gg g$, the fast oscillating terms can be ignored by performing RWA. Then the simplified Hamiltonian reads

[TABLE]

Moving to the rotation framework with respect to $U_{2}(t)=e^{-i\hat{H}^{\prime}_{0}t}$ with $\hat{H}^{\prime}_{0}=\frac{\Omega_{x}}{2}\hat{J}_{x}$, we obtain the following transformed Hamiltonian

[TABLE]

Let $\Omega_{x}=2\omega_{z}$. The Hamiltonian in Eq. (29) can be rewritten as

[TABLE]

When the parameters satisfy the conditions: $\omega_{z}\gg\Omega_{z}$, $\omega_{z}\gg\varepsilon-\omega_{x}$, $\omega_{z}\gg g$, we can neglect the fast oscillating terms and recast above Hamiltonian as

[TABLE]

where $\tilde{\omega}_{r}=\omega_{r}-\omega_{x}$, $\tilde{\varepsilon}=\Omega_{z}/2$, and $\tilde{g}=g/2$. Thus we obtain the effective Hamiltonian shown in Eq. (5), which is effective under the following conditions: $\Omega_{x}=2\omega_{z}$, $\omega_{x}\gg\Omega_{x}\gg g$, $\omega_{z}\gg\Omega_{z}$, $\omega_{z}\gg\varepsilon-\omega_{x}$. Such Hamiltonian is multi-qubit extension of the QRM with tunable parameters (i.e., the tunable Dicke model). The simulated coupling ratio is $\tilde{g}/\tilde{\omega}_{r}=g/[2(\omega_{r}-\omega_{x})]$, which is also turnable by adjusting the frequency of the transverse driving.

Acknowledgements

The work is supported by the NSF of China (Grant No. 11575042).

Author Contributions

G.W., C.W., and K.X. initiated the idea. C.S., J.L. and R.X. developed the model and performed the calculations. J.L. and R.X. provided numerical results. All authors developed the scheme and wrote the main manuscript text.

Additional Information

Competing Interests: The authors declare no competing interests.