Quantum criticality and state engineering in the simulated anisotropic quantum Rabi model

Promising applications of the anisotropic quantum Rabi model (AQRM) in broad parameter ranges are explored, which is realized with superconducting flux qubits simultaneously driven by two-tone time-dependent magnetic fields. Regarding the quantum phase transitions (QPTs), with assistant of fidelity susceptibility, we extract the scaling functions and the critical exponents, with which the universal scaling of the cumulant ratio is captured with rescaling of the parameters due to the anisotropy. Moreover, a fixed point of the cumulant ratio is predicted at the critical point of the AQRM. In respect to quantum information tasks, the generation of the macroscopic Schr\"{o}dinger cat states and quantum controlled phase gates are investigated in the degenerate case of the AQRM, whose performance is also investigated by numerical calculation with practical parameters. Therefore, our results pave a way to explore distinct features of the AQRM in circuit QED systems for QPTs, quantum simulations and quantum information processings.


I. INTRODUCTION
The quantum Rabi model (QRM) [1][2][3], describing the interaction between a two-level system and a single bosonic mode, H QRM = ωa † a + ω q σ z /2 + gσ x (a † + a), plays an important role in studying the dynamics of quantized lightmatter interaction. Here, a and a † are the annihilation and creation operators of the bosonic mode with frequencyω, while σ z and σ x = σ + + σ − are Pauli operators associated with a qubit with ground state |g , excited state |e , and transition frequencyω q . Typically, due to the parameter accessibility of most experiments, the QRM could be drastically simplified to the easily solvable Jaynes-Cummings model (JCM) [4], through the celebrated rotating-wave approximation (RWA). In this case, Rabi oscillations inside the JCM doublets or collapses and revivals of the system populations are paradigmatic examples of the intuitive physics behind the JCM dynamics [5].
However, recent experimental progresses in solid-statebased quantum systems have allowed the advent of the socalled ultrastrong coupling (USC) regime [6][7][8] and the deep strong coupling (DSC) regime [9,10] of light-matter interactions, where the coupling strength is comparable to (USC) or larger than (DSC) appreciable fractions of the mode frequency. In these regimes, the RWA breaks down and the QRM is again invoked. In addition to the relatively complex quantum dynamics provided by the QRM, it brings about novel quantum phenomena [11][12][13] and challenges in implementing quantum information tasks [14][15][16][17]. Although exciting, natural implementations of the QRM in the USC/DSC regime in * wlyou@suda.edu.cn † jqyou@csrc.ac.cn other platforms remain very challenging since they are confined by fundamental limitations. However, different schemes have been proposed to simulate the QRM using superconducting circuits [18], quantum optical systems [19], trapped ions [20] and cold atoms [21].
In the other aspect, the fascinating promises of the QRM has trigged many studies of the anisotropic quantum Rabi model (AQRM), see e.g., Refs. [22][23][24][25], H AQRM = ωa † a+ 2ω q σ z + g r (σ − a † +σ + a)+ g ar (σ + a † +σ − a), (1) which is a generalization of the QRM affiliated withg r and g ar denoting the different strengths of the rotating terms and the counter-rotating terms, respectively. The AQRM returns to the JCM withg ar = 0 or the original QRM withg r =g ar . The eigenspectra of the AQRM can be exactly mapped to the solutions of a transcendental equation [22] and approximate methods for solving the AQRM have also been widely studied [23,26,27]. Since individual addressing of the two coupling constants are allowed, the AQRM presents a favorable test bed for many valuable theoretical issues. For example, it provides the opportunity to explore the role of the counterrotating terms [23,28,29] and thus bridge the gap between the JCM and the QRM in the dynamics [30]. It was described in Refs. [22,31] and presented in the Appendix that, the rotating and the counter-rotating terms have unbalanced contributions to the dynamics of the AQRM. In addition, the AQRM plays a key role in establishing the universality scaling of the concept of quantum phase transition in few-body systems [25]. On the other hand, the anisotropy of the AQRM gives rise to the occurrence of the level crossings between the eigenvalues of different parity sectors [22,32] and the enhanced squeezing [26,27]. This leads to promising applications, such as Fisher information [32], quantum state engineering [24] and the squeezing enhancement based measurement [33], which is wildly used in gravitational wave detection [34]. The discussions of the anisotropy in the standard AQRM were further extended to the semi-classical case [35], the multi-qubit case [25,36] (namely the anisotropic Dicke model), and the non-Markovian case [29]. As a consequence, the theoretical advancements bring up experimental requests for individual adjustability of the coupling constantsg r andg cr in wide parameter ranges to demonstrate the innovative features of the AQRM. Although there have already been some experimental proposals for the realization of the AQRM in some systems, i.e., quantum well with spin-orbit coupling [32,37], and circuit QED systems [36,38], they are quite limited on the tunability and the achievable parameter ranges. Therefore, the demand to explore new platforms to study the dynamics of the AQRM is put forward.
In this work, we propose an experimentally feasible scheme to simulate the controllable AQRM demostrating the USC and DSC dynamics in a circuit QED system with only strong coupling between the qubit and the resonator. By resorting to the two-tone time-dependent magnetic fluxes threading the qubit, we show through analytical and numerical calculations that our proposal will not only have access to the regimes of USC and DSC but also realizing the controllable AQRM with individually tunable couplingsg r andg cr . While being important for understanding the transition from the JCM to the QRM in investigating the AQRM and discovering exotic quantum scenarios, our proposal will pave the way for the implementation of quantum simulators [39] for rich coupling regimes of light-matter interaction in systems where they are experimentally inaccessible. As we will discuss below, our scheme can be used to generate superpositions of coherent states with different phases, which can be directly used for the hardware-efficient quantum memory and quantum error corrections [40,41]. Moreover, our scheme finds applications in developing the quantum phase gates when generalized to multi-qubit cases.
The paper is organized as follows. We firstly describe in Sec. II the Hamiltonian of our qubit-resonator setup, where the flux qubit is controlled by bichromatic time-dependent magnetic fluxes. The effective AQRM is obtained when the frequency conditions are well-respected. In Sec. III, we present results of numerical calculations that demonstrate the performance of our scheme by studying the properties of the atomic-state population, the ground-state entanglement entropy, and the Wigner distribution function of the resonator's field. In Sec. IV, we discuss the possibility of preparing the Schrödinger cat states and generating the quantum controlled phase gates with our setup. The conclusions are presented in Sec. V.

II. THE QUBIT-RESONATOR CIRCUIT
For simplicity, but here without loss of generality, we use three-junction flux qubits (e.g., [42][43][44]) in our scheme. As shown in Fig. 1(a), a flux qubit is coupled to a LC circuit with an inductance L and a capacitance C. The mutual inductance between the flux qubit and the LC circuit is M. The applied magnetic flux Φ through the flux qubit loop in Fig. 1(a), which controls the qubit-resonator couplings, is assumed to include a static magnetic flux Φ e 0 , and also two time-dependent magnetic fields (TDMFs), Φ e j (t) = A j cos(ω j t + ϕ j ). Here j = r, d label the two TDMFs, individually. Considering one threejunction flux qubit, the qubit's Hamiltonian H q reads where we have assumed each junction in the flux qubit has a capacitance C Ji , phase drop φ i = 2πΦ i /Φ 0 , Josephson energy E Ji , and critical current I 0i = 2πE Ji /Φ 0 . Here, Φ 0 = h/2e is the magnetic flux quantum. With current-phase relation, the super-current for each junction reads I i = I 0i sin φ i .
And thus the persistent current in the qubit loop is [44,45], where C q is the total capacitance of the flux qubit, C −1 Ji , with the convention C J3 = ηC J1 = ηC J2 and η being the relative size of the Josephson junction. Taking into account the TDMFs, the flux quantization around the qubit's loop imposes a constraint on the phase drop across the three junctions [42,44,45], In order to define an effective qubit within the junction architecture, we diagonalize the Hamiltonian containing only the junctions Eq.(2) in absence of the TDMFs, i.e., Φ e j (t) = 0. The two lowest eigenstates are labeled as the eigenstates of σ z , i.e., |g and |e , and the spanned two-dimensional subspace describes the effective qubit.
In the other aspect, the Hamiltonian of the total system is written as H = H q + H c + I MI q , where H c = Q 2 /2C +Φ 2 /2L is the Hamiltonian of the LC circuit, with Q being the capacitor's charge, Φ = IL being the magnetic flux through the LC circuit loop and I is the inductor's current. The Hamiltonian of the LC circuit H c can be simply quantized by introducing the annihilation and creation operators Φ = √ /(2Cω)(a + a † ) and  {|g , |e }, we obtaiñ The first two terms in Eq.(3) denote the free Hamiltonians of both the qubit and the LC circuit, where ω q is the transition frequency of the effective qubit. The third term in Eq. (3) represents the qubit-resonator interaction with the coupling strength being HereĨ and f e ≡ Φ e 0 /Φ 0 is the reduced dc bias magnetic flux. The fourth term in Eq.(3) plays the role of a driving Hamiltonian representing the interaction between the qubit and its TDMFs with the respective driving strength being The fifth term of Eq. (3) is the controllable nonlinear interaction among the qubit, the resonator, and the TDMFs, with the respective coupling strength being As noticed above, the TDMFs Φ e j (t) equal to zero when calculating the coupling strengths g, Ω j , and Λ j . It is worthy noting that in the above derivations, we keep the time-dependent amplitudes small such that the reduced time-dependent magnetic fluxes satisfy | f j (t)| 10 −3 . This leads to: (1) the approximation of sin[2π f j (t)] ∼ 2π f j (t) and cos[2πΦ e j (t)/Φ 0 ] ∼ 1; (2) the ignorance of the interaction terms controlled by two simultaneously applied TDMFs in the form of ∼ f 1 (t) f 2 (t). As a result, when expanding the potential energy in Hamiltonian Eq.(2) and the qubit's loop current I q , we only need to keep the first order of the small reduced flux Φ e j (t)/Φ 0 . This weak-amplitude approximation of the control fields also prevents the qubit from deviating much from the dc bias point that is set to be the optimal coherence point in our case. Moreover, this weak-amplitude approximation helps to reduce unwanted possible excitations to the higher-energy states outside the computational state space when we make the two-level (qubit) approximation.
When the conditions g |ω q ± ω| and Ω j / |ω q ± ω j | are satisfied, Eq.(3) can be written as where we have performed the RWAs and neglected all terms that are fastly oscillating in the interaction picture with respect to the system's free Hamiltonian ω q σ z /2 + ωa † a. We next consider the case where the TDMFs are inducing the respective first-order red (r) and blue (b) sideband transitions with detunings δ r and δ b onto the qubit-resonator system, e.g., ω r = ω q − ω − δ r and ω b = ω q + ω − δ b . In such a scenario, when the rest of the frequency detunings are large compared to the coupling parameters, i.e., |ω − ω r ± ω q | Λ r and |ω − ω q ± ω b | Λ b , one may neglect the rest of the fastoscillating terms. These approximations lead to a simplified time-dependent Hamiltoniañ It is worth noting that Eq.(8) corresponds to the interaction picture of the generalized AQRM with respect to the uncoupled [20], such that with the effective parameters being Here the qubit's and the resonator's frequencies are represented by the sum and the difference of the two detunings, respectively. The tunability of these parameters permits the study of all coupling regimes of the AQRM via suitable choices of the amplitudes and the detunings of the TDMFs. It is noteworthy that the complex coupling strengthsg r and g ar can be realized by choosing the phases ϕ r and ϕ b of the TDMFs. For example, Eq.(9) reduces to the standard AQRM in Eq.(1) when ϕ r = ϕ b = 0. To go beyond the USC, we only need the condition that Λ j |δ b − δ r |.

III. NUMERICAL ANALYSIS
A. The probability and the entanglement entropy To study the feasibility of our proposal, we have performed numerical calculations with realistic parameters in circuit QED systems [44][45][46]. As an illustration, we make comparisons between the original Hamiltonian in Eq.(3) and the effective Hamiltonian in Eq.(9) for the ground-state probability P g (t) = | Ψ G |g | 2 and the ground-state entanglement entropy S G = − Tr ρ q G log 2 (ρ q G ) , where |Ψ G is the ground state of the total system, ρ q G = Tr f ρ G is the reduced density matrix of the qubit's subsystem by tracing out the field's degree of freedom. The ground-state probability P g (t) indicates the 005ω. This leads to simulated effective parameters ofω q =ω =g r = 1,g cr = 2 for (a, b); andω q =ω =g cr = 1,g r = 2 for (c, d). For the simulation, the system is initially prepared in the ground state of the whole system |g, 0 , and the rest of the parameters for the simulation are chosen as ω/2π = 3GHz, ω q = 3ω, ω r = 2ω, ω b = 3.995ω, g = 0.01ω, and φ r = φ b = 0.
atomic-excitation probability in the ground state |Ψ G of the total system, and the entanglement entropy S G measures the entanglement between the qubit and the resonator in |Ψ G . For the simulations, we assume that the original undriven qubit-resonator system is prepared in its ground state |g, 0 , which works only in the strong coupling regime. Then at time t = 0, we switch on the external TDMFs, which act as external drivings onto the qubits. Taking the close-system assumption, the whole system evolves according to the unitary operator, which is computed by integrating the time-dependent Hamiltonian equation Eq.(3). We consider two sets of parameters in Fig. 2: (a, b) Ω r = Λ r = 2π × 15 MHz and MHz, which result in the simulated values:ω q =ω =g r = 1,g cr = 2 for Fig. 2(a, b); and ω q =ω =g cr = 1,g r = 2 for Fig. 2(c, d), respectively. Clearly, the results from the original Hamiltonian (red solid lines) are completely consistent with the ones from the effective Hamiltonian (blue dashed lines with circles) for the two sets of parameters. It is also obvious that such strong driving amplitudes of about tens of megahertz are sufficient enough to simulate the dynamics of a strongly-coupled system reaching and even beyond the USC regime. As seen from Fig. 2(b,  d), the ground state |Ψ G can be highly entangled between the qubit and the resonator for a long period of evolution. For the simulation, the original qubit-resonator coupling strength is chosen as g = 2π × 30 MHz, which ensures that the original system is in the strong coupling regime, and the rest of the parameters are set to ω = 2π × 3 GHz, ω q = 2π × 9 GHz, ω r = 2π × 6 GHz, ω b = 2π × 11.985 GHz, and φ r = φ b = 0.
A very special and interesting aspect of our scheme is that it can be used to simulate the dynamics of the controllable AQRM with a degenerate qubit and a degenerate resonator, i.e.,ω =ω q = 0, which can obtained by tuning the frequencies of the TDMFs such that δ r = δ b = 0. Similar to the above discussions, numerical calculations are performed for the atomic-ground-state probability P g (t) and the ground-state entanglement entropy S G for three sets of parameters, which effectively gives the simulated parameters ofω q =ω = 0 with 005ω. This lead to simulated effective parameters ofω q =ω = 0 withg r = 1,g cr = 0.2 for (a, b);g r =g cr = 1 for (c, d); andg r = 0.2, g cr = 1 for (e, f). For this simulation, the system is also initially prepared in the ground state of the whole system |g, 0 , and the rest of the parameters for the simulation are chosen as ω/2π = 3 GHz, ω q = 3ω, ω r = 2ω, ω b = 4ω, g = 0.001ω, and φ r = φ b = 0. g r = 1,g cr = 0.2 for Fig. 3(a, b);g r =g cr = 1 for Fig. 3(d,  d); andg r = 0.2,g cr = 1 for Fig. 3(e, f). In Fig. 3, the relative ratio ofg cr /g r is increased gradually in top-down order for each column and the corresponding curves for the groundstate probabilities show that periodic and complete atomic population transfer occurs wheng cr /g r = 1 andg cr /g r = 5, as displayed in Fig. 3(c) and Fig. 3(e), respectively. By carefully choosing the parameters to well respect the required conditions, the curved lines for the original Hamiltonian in Eq.(3) (red solid line) reproduce the ones calculated for the effective Hamiltonian in Eq.(9) (blue dashed lines with circles) with high accuracy. The numerical agreements shown in both Fig. 2 and Fig. 3 prove that our scheme has excellent performance in simulating static properties and the dynamics of the AQRM in both the USC and the DSC regimes.

B. The statistics of the fields
What coming along with the atomic population transfers are the collapses and revivals of the photon wave packets and the variation of the photon statistics. In the following, by employing the Wigner quasi-probability distribution function (WF), we show some interesting features of the field statistical properties of the degenerate AQRM withω =ω q = 0.
The non-classicality of a bosonic field can be signaled by the WF, which is defined as where ρ f = Tr q [ρ] is the reduced field-density matrix after the qubit is traced out, and D(α) = exp [αa † − α * a] is the coherent displacement operator with amplitude α. In Fig. 4, we plot the WF of the AQRM at different time intervals for four sets of parameters with the initial state |g, 0 . The top row of Fig. 4(a-d) depicts the evolution of the WF of the field generated whenω q =ω =g r = 0,g cr = 1, which corresponds to the population transfer between the states of |g, 0 and |e, 1 , and the WF of the single photon Fock state is shown in Fig. 4(c) at timeg cr t/2π = 0.25. The third row of Fig. 4(i-l) shows the evolution of the WF of the field generated whenω q =ω, g r = 1, andg cr = 1, which describes a mixture of two coherent states |±α with time-dependent displacement amplitude of α(t) = −iḡ r t [47]. The amplitudes of the coherent states ideally increase linearly and practically, they will be prevented from diverging into instability by the damping of the oscillator and the finite duration of the evolution. It is noted that the small distortion of the WF from the ones of the ideal coherent state is due to a small deviation of our scheme from the effective ones for longer evolution time. The second row and the bottom row of Fig. 4 display the field properties with unbalanced and nonzero rotating and counter-rotating coupling terms in the degenerate AQRM, i.e.,ω q =ω,g r = 0.5 and g cr = 1 for Fig. 4(e-h) andω q =ω,g r = 1, andg cr = 0.5 for Fig. 4(m-p). In these two cases, both the rotating and counterrotating terms contribute to the dynamics of the system, but unbalanced. An intuitive picture to understand these figures could be the following. Started from |g, 0 , the photons spreadg cr t/2π = 0.13  ) with |g, 0 being the initial state of the system. We have considered four cases: (a-d) Ω r = Λ r = 0 and Ω b = Λ b = 0.005ω, which corresponds to simulated effective parameters ofω q =ω =g r = 0,g cr = 1; (e-h) Ω r = Λ r = 0.0025ω and Ω b = Λ b = 0.005ω, which corresponds to simulated effective parameters ofω q =ω = 0,g r = 0.5,g cr = 1; (i-l) Ω r = Λ r = 0.005ω and Ω b = Λ b = 0.005ω, which corresponds to simulated effective parameters ofω q =ω = 0,g r =g cr = 1; (m-p) Ω r = Λ r = 0.005ω and Ω b = Λ b = 0.0025ω, which corresponds to simulated effective parameters ofω q =ω = 0,g r = 1,g cr = 0.5. The rest of the parameters for the simulation are chosen as ω/2π = 3GHz, ω q = 3ω, ω r = 2ω, ω b = 4ω, g = 0.001ω, and φ r = φ b = 0.
independently along the even parity chain, and thus produce a qubit-resonator entangled state. Such entangled state has the properties of the displaced squeezed states, whose squeezing parameters are functions of the relative ratiog cr /g r .

IV. QUANTUM INFORMATIONS APPLICATIONS OF OUR SCHEME
Without lose of generality, our scheme of simulating the controllable AQRM can be generalized to the multi-qubit case, where multiple flux qubits are coupled to the LC circuit as shown in Fig. 1(b). When the corresponding conditions for each qubit to realize the effective anisotricpic Rabi model as in Eq.(9) are well satisfied, and the parameters for the lth qubit (l = 1, 2, ...N) are chosen such that, δ l b = −δ l r = δ, ϕ l b = −ϕ l r = ϕ, and Λ l r = Λ l b = Λ, the effective multi-qubit Hamiltonian can be written as where J x = N l=1 σ l x are the collective qubit operators and we have definedḡ ≡ Λ/2. The evolution operator U δ,ϕ can then be found as whereΦ(t) = (ḡ/δ) 2 (δt − sin δt) and D(χ) = exp [χ(t)a † − χ * (t)a] is the displacement operator with χ(t) = (ḡ/δ) e iϕ (1 − e iδt ) being the displacement amplitude of the oscillator.

A. The generation of Schrödinger cat states
It has been proved that, the Schrödinger cat states have promising applications in hardware-efficient quantum memory and quantum error corrections [40,41]. In this section, we show the performance of our scheme in generating this class of non-classical states with both theoretical and numerical approaches. In the single-qubit case, J x in Eqs. (13,14) are replaced by J x = σ x , and when the initial state of the whole system is prepared in the ground state as |Ψ 1 (0) = |g, 0 , we obtain the final state at time t as where |± = (|e ± |g ) / √ 2, and |±α(t) being the coherent states of the harmonic oscillator, which are of the same amplitude but opposite phase in the phase space. α(t) = (1 − e iδt )e iϕ /δ is the coherent-state amplitude for the singlequbit case. It is worthy noting that from Eq.(15) that, since σ 2 x = I, the first term in Eq. (14) behaves only as a global phase factor in Eq. (15). Obviously, depending on the states of the flux qubit |± , the coherent states undergo different displacements |±α(t) , respectively. In the bases of the {|e , |g }, the state in Eq. (15) can be rewritten as where |C ± α ≡ N(|α ± |−α ) with N = 1/ 2(1 ± e −2|α| 2 ) ≈ 1/ √ 2 are the so-called even (|C + α ) and odd (|C − α ) Schrödinger cat states. By choosing the phase ϕ, the four types of the quasi-orthogonal states |±α and |±iα [40], i.e., | α|iα | 2 << 1 (note that for α = 2, | α|iα | 2 < 10 −3 ), can be generated by measuring the qubit in the |± bases. By performing projective measurements in the qubit {|e , |g } bases, the oscillator will collapse into the Schrödinger cat states with probability of (1 ± e −2|α| 2 )/2, respectively. In Fig. 5, we show that the Wigner distribution of the even cat state |C + α , which is generated from a projective measurement of the system's state onto the qubit state |g . One can easily find that Fig. 5 is different from Fig. 4(i), in the sense that, Fig. 5 displays the quantum coherence between the two coherent states |±α with opposite phases, while Fig. 4(i) only indicates a mixture of the coherent states |±α .
Indicated from Eq.(16) that, the maximum displacement amplitude is |α| max = 2ḡ/δ, and it can be obtained at the times t = (2m + 1)π/δ for natural number m. By choosing a small value for δ and a large effective coupling strengthḡ, we can create macroscopically distinct Schrödinger cat states of considerable size of |α| > 1.
In another aspect, the displacement amplitude of the Schrödinger cat states can be further enhanced with even number of flux qubits by exploring the multi-qubit case and preparing the system in the state of |Ψ N (0) = (|+, +, ..., + − |−, −, ..., − ) |0 / √ 2. In this case, the state after evolution is given by where the coherent-state amplitude is enhanced by a factor N, and the first term in Eq. (14) remains as a global phase factor. However, collective measurements on the flux qubit in the bases of (|+, +, ..., + ± |−, −, ..., − )/ √ 2 are required to obtain the Schrödinger cat states with an enhanced amplitude.
B. The two-qubit controlled quantum phase gate generation As seen from the evolution operator Eq. (14), the Hamiltonian in Eq.(13) introduces qubit-qubit interaction between any pair of qubits. And thus our circuit can be used to generate quantum gates and produce highly-entangled states between qubits. Let the system evolve for a time period of T = 2π/δ, we obtain χ(T ) = 0 and up to an overall phase factor, the evolution operator can be recast as withθ = 2Φ(T ) = 4πḡ 2 /δ 2 . In the following, we show that the generation of a two-qubit quantum controlled-NOT gate is straightforward from Eq. (18). In the two-qubit bases of {|ee , |eg , |ge , |gg }, the evolution operator can be expressed as which represents the non-trivial two-qubit gates whenθ mπ (m = 0, 1, 2, · · · ). Specifically, whenθ = π/4 (i.e.,ḡ = δ/4), U(T ) is locally equivalent to the controlled-NOT (CNOT) gate.

V. CONCLUSIONS
In conclusion, we have proposed an experimentally realizable scheme to simulate the dynamical properties of the anisotropic quantum Rabi model in both the ultrastrong and the deep strong coupling regimes with superconducting circuits, which is composed of flux qubits strongly-coupled to microwave resonators. This is achieved by applying twotone time-dependent magnetic fields into the loop of the flux qubits. By carefully choosing the parameters of the timedependent magnetic fields, our model provides excellent controllability of all the parameters in the simulated anisotropic quantum Rabi model, where the effective qubit's and resonator's frequencies can be tuned by the frequency of the time-dependent magnetic fields matching or mismatching to the detuning (or sum) of the qubit-resonator frequencies; the rotating and counter-rotating coupling terms can be individually tuned by the amplitudes of the time-dependent magnetic fields. Along with theoretical arguments, we study the performance of our setup by numerically comparing the atomicstate probability and the ground-state entanglement entropy of the effective Hamiltonian with those of the original Hamiltonian. The Wigner distributions are also investigated to demonstrate the nonclassial properties of the resonator's field in the anisotropic quantum Rabi model. Therefore, our proposed quantum simulation of the anisotropic quantum Rabi model in broad parameter ranges of the light-matter interaction may become a building block in simulations of physics inaccessible in standard quantum optics. Last but not least, as for the applications for quantum information processing, we show that our scheme can be easily scaled up to generate macroscopic Schrödinger cat states and quantum controlled phase gates. As presented in Refs. [22,31] and shown in Fig. A1, the rotating and the counter-rotating coupling terms have unequal contributions to the dynamics of the AQRM. For further illustration, we show the properties of the AQRM's ground-state |0 on the rotating coupling coefficientg r and the counter-rotating coupling coefficientg ar in Fig. A1, i.e., the ground-state energỹ E˜0, the ground-state entanglement entropyS˜0, the ground-state photon-populationÑ r 0 = a † a 0, and the ground-state atomic-populationP ẽ 0 = σ + σ − 0. It is clear that, all of the properties we investigated are more sensitive to the counter-rotating coupling coefficientg ar than the rotating coupling coefficientg r . The entanglement entropy in Fig. A1(b) indicates that the ground state |0 is highly entangled when the AQRM works in the DSC regime with bothg r /ω 1 andg ar /ω 1.

VI. ACKNOWLEDGEMENTS
On the other hand, the QRM has substantial differences as compared to the JCM, since the QRM respects a discrete Z 2 symmetry instead of the continuous U(1) symmetry in the JCM. To be more precise, the parity P = e iπ(a † a+σ + σ − ) is conserved in the QRM while the total excitation number N = a † a + σ + σ − is maintained in the JCM. In addition, the competition between the rotating and counter-rotating interaction terms may give rise to the occurrence of the level crossings between the eigenvalues of different parity sectors [22]. This further leads to a sharp discontinuity of the ground-state entropy at the level-crossing point. This phenomenon does not occur in the isotropic QRM. , the photon-population (c), and the atomic-population (d) in the ground state of the AQRM as a function of the normalized rotating (g r /ω) and counter-rotating coefficients (g ar /ω). The parameters used for the numerical analysis areω q = 0.8ω.