1. Introduction
Phase transition (PT) is a fundamental issue in statistical and condensed matter physics. Recent advances in solid-state quantum devices and quantum simulation have pushed light–matter interaction systems into the ultrastrong and deep-strong coupling regimes, leading to novel collective phenomena [1, 2]. Unlike correlated solid-state materials, which are sensitive to impurities, environments, and temperature effects that complicate parameter tuning, light–matter interaction systems (e.g., superconducting circuit quantum electrodynamics (QED) [3], ion traps [4, 5], and cavity-embedded Bose–Einstein condensates [6]) offer advantages such as high tunability and scalability, making them ideal for studying quantum phase transitions (QPTs) and their critical properties. Theoretically, the Dicke model (DM) [7–14] and the spin-boson model [15–17, 18] in light-matter interacting systems can exhibit QPTs at strong coupling between the two-level systems (qubits) and the cavity or bosonic baths.
In recent decades, it has been proposed that the quantum Rabi model (QRM), consisting of a single qubit and a single-mode cavity, can undergo a QPT at an infinite ratio of qubit to cavity frequencies [19], sparking a surge of studies on the so-called finite-component QPT [20–28]. It is widely accepted that both the DM and the QRM undergo a single QPT from the normal phase to the superradiant phase (SRP), exhibiting the same critical behavior [20].
Unlike the conventional cavity QED system, artificial qubits in modern solid-state devices typically feature both the energy splitting and bias (ϵ) between two qubit states [1, 3], making the asymmetric QRM and Dicke models prevalent. The Z2 symmetry (parity) is typically broken by finite biases [29], which prohibits the QPTs in these models. Fortunately, if staggered biases—where half of the qubits have a positive bias and the other half a negative bias (−ϵ) with identical values—are applied, the Z2 symmetry of the entire system is restored, making the phase transition (PT) possible [18, 30, 31]. Advances in experimental techniques and quantum information science have enabled precise control over the finite bias energies of qubits, leading to the emergence of rich quantum phases and multiple criticalities. Especially in the quantum simulation of a single ytterbium ion trap, the QPT of the QRM has been convincingly observed by measuring the population of the spin-up state and the average phonon number of the ions [5]. Designing extended yet simplified light–matter interaction systems, such as the generalized QRM, to realize complex quantum phase diagrams similar to those in correlated solid-state materials, will significantly enhance our understanding of QPTs.
In this paper, we investigate the superradiant phase transitions (SRPTs) of the two-qubit QRM with staggered qubit biases (±ϵ). In the limit of an infinite qubit-to-cavity frequency ratio, equivalent to the thermodynamic limit, we analytically derive the mean-field Hamiltonian and the order-parameter-dependent energy density functional, which provides the ground energy and the ground-state order parameter. By combining the analysis based on the mean-field energy density functional with numerical results at large effective sizes, we derive and confirm the rich SRPTs, including both second- and first-order SRPTs and a tricritical point (TCP). Specifically, we derive the analytical expression for all PT points, including the nonperturbative one of the first-order SRPT. Analytical methods are crucial for understanding the intrinsic physics of such systems. Finally, we perform finite-size scaling analysis and surprisingly find different critical exponents at the TCP, including both the correlation-length and order-parameter exponents, in contrast to those of the original second-order SRPTs. This implies that the TCP of the two-qubit QRM with staggered qubit biases belongs to a distinct universality class compared to the second-order SRPTs in the original QRM. This work provides reliable theoretical frameworks for designing new experimental platforms to explore QPTs.
The paper is organized as follows: section 2 provides a brief introduction to the two-qubit QRM with staggered qubit biases. In section 3 , we derive the analytical expressions for all critical points, including that of the nonperturbative first-order QPT. A detailed phase diagram is constructed in the large limit of the frequency ratio. These analytical findings are numerically verified in section 4 , where the critical behavior at the triple point is distinguished through finite-size scaling analysis. A brief summary is presented in section 5 .
2. Model and its symmetry
The two-qubit Rabi model with staggered biases describes a single-mode cavity field interacting with two identical two-level systems (qubits) with opposite biases (±ϵ), and its Hamiltonian is given by:
$\begin{eqnarray}\begin{array}{rcl}H & = & \omega {a}^{\dagger }a+\frac{{\rm{\Delta }}}{2}\left({\sigma }_{z}^{1}+{\sigma }_{z}^{2}\right)+\frac{\varepsilon }{2}\left({\sigma }_{x}^{1}-{\sigma }_{x}^{2}\right)\\ & & +g(a+{a}^{\dagger })\left({\sigma }_{x}^{1}+{\sigma }_{x}^{2}\right),\end{array}\end{eqnarray}$
where a and a† are the annihilation and creation operators of cavity field, ${\sigma }_{x}^{n}$ and ${\sigma }_{z}^{n}$ are the Pauli matrices of the corresponding qubit labeled by n = 1, 2, ω denotes the frequency of cavity field, Δ is the common energy splitting of the qubits, ϵ is the strength of the qubit bias and g is the common qubit-photon coupling strength.Both the QRM and DM are well-known for their robustness against parity transformation $\left(a,{a}^{\dagger },{\sigma }_{x},{\sigma }_{y}\right)\to \left(-a,-{a}^{\dagger },-{\sigma }_{x},-{\sigma }_{y}\right)$. The parity symmetry divides the Hilbert space into two independent subspaces based on the even and odd total excitation numbers, with $\hat{N}={a}^{\dagger }a+\left({\sigma }_{z}+1\right)/2$. Moreover, this symmetry enables the SRPT of the DM in the limit of infinite number of qubits [7–10] and the SRPT of the QRM in the limit η = Δ/ω → ∞ [19, 20], which has been studied for decades. However, the qubit bias mixes the even and odd subspaces of the total excitation number, resulting in the vanishing of the parity symmetry, represented by the operator ${{\rm{\Pi }}}_{\varepsilon =0}={{\rm{e}}}^{{\rm{i}}\pi \hat{N}}$. Despite this, staggered biases can recover the ${{\mathbb{Z}}}_{2}$ symmetry if both qubits and their couplings to the cavity field are tuned identically. The symmetry operator of the two-qubit QRM with staggered biases is given by 1 ), remains invariant under the above transformation, and the ${{\mathbb{Z}}}_{2}$ symmetry is restored. This paves the way for the SRPT, which is why the novel QPTs induced by the staggered biases have been explored in both the Dicke model [30–32] and two-qubit Spin-Boson model [18].
$\begin{eqnarray}{\rm{\Pi }}=\left[\frac{{\sigma }_{z}^{1}{\sigma }_{z}^{2}+1}{2}-\left({\sigma }_{+}^{1}{\sigma }_{-}^{2}+{\sigma }_{-}^{1}{\sigma }_{+}^{2}\right)\right]{{\rm{e}}}^{{\rm{i}}\pi {a}^{\dagger }a},\end{eqnarray}$
where the raising and lowering operators of the Pauli matrices are ${\sigma }_{\pm }^{1(2)}=\left[{\sigma }_{x}^{1(2)}\pm {\rm{i}}{\sigma }_{y}^{1(2)}\right]/2$. It is clear that the transformation associated with the above symmetry leaves the photon number operator a†a and the collective spin operator ${\sigma }_{z}^{1}+{\sigma }_{z}^{2}$ unchanged. Furthermore, it performs a parity transformation along with qubit exchange, i.e., $\left(a,{a}^{\dagger },{\sigma }_{x}^{1},{\sigma }_{x}^{2}\right)\to \left(-a,-{a}^{\dagger },-{\sigma }_{x}^{2},-{\sigma }_{x}^{1}\right)$ [31]. It is clear that the present model, as denoted in equation (3. Analytical results
3.1. Mean-field Hamiltonian
We employ mean-field theory to explore the ground state phase diagram of this model. To derive the mean-field energy density functional, we first rescale the Hamiltonian by the effective system size η = Δ/ω to make it dimensionless, i.e.,
$\begin{eqnarray}\begin{array}{rcl}\tilde{H} & = & 2H/{\rm{\Delta }}=\frac{2}{\eta }{a}^{\dagger }a+\left({\sigma }_{z}^{1}+{\sigma }_{z}^{2}\right)+\tilde{\varepsilon }\left({\sigma }_{x}^{1}-{\sigma }_{x}^{2}\right)\\ & & +\frac{\tilde{g}}{\sqrt{\eta }}(a+{a}^{\dagger })\left({\sigma }_{x}^{1}+{\sigma }_{x}^{2}\right),\end{array}\end{eqnarray}$
with the rescaled qubit-photon coupling $\tilde{g}=2g/\sqrt{{\rm{\Delta }}\omega }$ and the rescaled bias strength $\tilde{\varepsilon }=\varepsilon /{\rm{\Delta }}$. In terms of the field quadratures, the position operator $x=(a+{a}^{\dagger })/\sqrt{2}$ and momentum operator ${p}_{x}={\rm{i}}({a}^{\dagger }-a)/\sqrt{2}$, the rescaled Hamiltonian is rewritten as $\begin{eqnarray}\begin{array}{rcl}\tilde{H} & = & \frac{1}{\eta }({x}^{2}+{p}_{x}^{2})+\left({\sigma }_{z}^{1}+{\sigma }_{z}^{2}\right)+\tilde{\varepsilon }\left({\sigma }_{x}^{1}-{\sigma }_{x}^{2}\right)\\ & & +\frac{\sqrt{2}\tilde{g}}{\sqrt{\eta }}x\left({\sigma }_{x}^{1}+{\sigma }_{x}^{2}\right).\end{array}\end{eqnarray}$
Inspired by previous works where the SRPT induced by the xσx-type coupling, also known as the isotropic coupling, is characterized by macroscopically excited ⟨x2⟩ [20, 22], we introduce a new set of field quadratures: the position operator $y=x/\sqrt{\eta }$ and the corresponding momentum operator ${p}_{y}=\sqrt{\eta }{p}_{x}$. The rescaled Hamiltonian is then expressed in terms of the new field quadratures as
$\begin{eqnarray}\begin{array}{rcl}\tilde{H} & = & {y}^{2}+\frac{1}{{\eta }^{2}}{p}_{y}^{2}+\left({\sigma }_{z}^{1}+{\sigma }_{z}^{2}\right)+\tilde{\varepsilon }\left({\sigma }_{x}^{1}-{\sigma }_{x}^{2}\right)\\ & & +\sqrt{2}\tilde{g}y\left({\sigma }_{x}^{1}+{\sigma }_{x}^{2}\right),\end{array}\end{eqnarray}$
where the new momentum term is proportional to 1/η2 and can thus be neglected. It should be noted that the omission of the new momentum term is a type of mean-field approximation. Using this approximation, the effective Hamiltonian is given by $\begin{eqnarray}{\tilde{H}}_{\,\rm{eff}\,}={y}^{2}+\left({\sigma }_{z}^{1}+{\sigma }_{z}^{2}\right)+\tilde{\varepsilon }\left({\sigma }_{x}^{1}-{\sigma }_{x}^{2}\right)-\sqrt{2}\tilde{g}y\left({\sigma }_{x}^{1}+{\sigma }_{x}^{2}\right),\end{eqnarray}$
where the new position operator y corresponds to a classical coordinate and can be further regarded as the order parameter describing the coarse-grained state space in Landau's theory of phase transitions.3.2. Mean-field phase diagram including a tricritical point
Based on the mean-field Hamiltonian (6 ), four branches of the energy spectrum can be obtained as
$\begin{eqnarray}{E}_{1,2}={y}^{2}-\sqrt{2}\sqrt{{f}_{1}\pm \sqrt{{f}_{1}^{2}-{f}_{2}}},\end{eqnarray}$
$\begin{eqnarray}{E}_{3,4}={y}^{2}+\sqrt{2}\sqrt{{f}_{1}\pm \sqrt{{f}_{1}^{2}-{f}_{2}}},\end{eqnarray}$
with ${f}_{1}={\tilde{\varepsilon }}^{2}+2{\tilde{g}}^{2}{y}^{2}+1$ and ${f}_{2}=8{\tilde{\varepsilon }}^{2}{\tilde{g}}^{2}{y}^{2}$. Focusing on the ground energy from the first branch, the mean-field energy density functional E(y) is denoted as $\begin{eqnarray}E(y)={y}^{2}-\sqrt{2}\sqrt{{f}_{1}+\sqrt{{f}_{1}^{2}-{f}_{2}}}.\end{eqnarray}$
Next, minimizing E(y) with respect to y yields the ground energy, and the location of the minimum provides the order parameter at the ground state. It is clear that E(y) is an even function of y, and the series expansion, retaining terms up to the sixth-order perturbation in y, is then
$\begin{eqnarray}E(y)\approx -2\sqrt{{\tilde{\varepsilon }}^{2}+1}+{c}_{2}{y}^{2}+{c}_{4}{y}^{4}+{c}_{6}{y}^{6},\end{eqnarray}$
where the expansion coefficients are given by $\begin{eqnarray}{c}_{2}=-2\left[\frac{{\tilde{g}}^{2}}{{(1+{\tilde{\varepsilon }}^{2})}^{2}}-\frac{1}{2\sqrt{1+{\tilde{\varepsilon }}^{2}}}\right]\sqrt{1+{\tilde{\varepsilon }}^{2}},\end{eqnarray}$
$\begin{eqnarray}{c}_{4}=-\frac{2{\tilde{g}}^{4}}{{(1+{\tilde{\varepsilon }}^{2})}^{4}}\left({\tilde{\varepsilon }}^{2}-\frac{1}{2}\right)\sqrt{1+{\tilde{\varepsilon }}^{2}},\end{eqnarray}$
$\begin{eqnarray}{c}_{6}=-\frac{2{\tilde{g}}^{6}}{{(1+{\tilde{\varepsilon }}^{2})}^{6}}\left[\frac{1}{2}-6{\tilde{\varepsilon }}^{2}+4{\tilde{\varepsilon }}^{4}\right]\sqrt{1+{\tilde{\varepsilon }}^{2}}.\end{eqnarray}$
Using Landau's theory of phase transitions, the condition for a second-order QPT is c2 = 0, c4 > 0. The critical coupling ${\tilde{g}}_{c}$ is thus given by the following expression: 9 ) may have three local minima at y = 0 and y±. Across the PT point, E(y = 0) shifts from being the smallest minimum to a slightly larger one (E(y = 0) > E(y±)). Despite the minor change in Ey, it results in the dramatic jump in the the minimum point, implying that the ground-state order parameter jumps from zero to y±. As shown in the right side of figure 1 for $\tilde{\varepsilon }\gt 1/2$, the cavity field excitation ⟨x2⟩ jumps from nearly zero to a large value of the same order as the effective size, indicating a first-order SRPT. We will derive the first-order PT point nonperturbatively below.
$\begin{eqnarray}{\tilde{g}}_{c}=\frac{\sqrt{2}}{2}{({\tilde{\varepsilon }}^{2}+1)}^{3/4},\tilde{\varepsilon }\lt \frac{1}{2}.\end{eqnarray}$
As shown on the left side of figure 1 with $\tilde{\varepsilon }\leqslant 1/2$, when $\tilde{g}$ exceeds the critical coupling ${\tilde{g}}_{c}$ (denoted by the red line), the ground-state photon field is gradually excited, with ${\langle {x}^{2}\rangle }_{0}$ of the same order as the effective size η = 1000, implying the presence of the second-order SRPT [19, 20]. Moreover, the condition for the TCP is c2 = c4 = 0, c6 > 0, and further the TCP $({\tilde{g}}_{tc}$, ${\tilde{\varepsilon }}_{tc})$ is obtained as $\begin{eqnarray}{\tilde{g}}_{tc}=\frac{{5}^{\frac{3}{4}}}{4},{\tilde{\varepsilon }}_{tc}=\frac{1}{2}.\end{eqnarray}$
It is well-known that, along the phase boundary, the second-order PT changes to a first-order PT upon crossing the TCP. Specifically, if c4 becomes negative and c6 > 0 for the stability of the mean-field functional, equation (
Figure 1. Numerical calculation of ${\langle {x}^{2}\rangle }_{0}$ for the ground state in the plane defined by the rescaled qubit-photon coupling strength, $\tilde{g}$, and the rescaled staggered biases, $\tilde{\varepsilon }$, for a large-size system with ω = 1 and Δ = 1000 using exact diagonalization. The red line indicates the QPT point of the second-order SRPT, as given by equation ( |
3.3. Nonperturbative first-order phase transition point
Proceeding to the first-order SRPT, considering the jump of the ground-state order parameter across the PT point, the approximation of the weakly excited cavity field, and further the perturbative series expansion (equation (9 )), may become invalid. We thus turn to derive the first-order phase boundary based on the properties of the solutions to the inequality E(y) ≤ E(y = 0) using the complete expression, equation (8 ).
The inequality $E(y)\leqslant E\left(y=0\right)$ admits nontrivial solutions only at and after the QPTs. In particular, the nontrivial solutions of the inequality $E(y)\leqslant E\left(y=0\right)$ distinguish between second-order and first-order QPTs. Specifically, just after the second-order PT, where E(y = 0) transitions from a minimum (before the second-order PT) to a local maximum, a neighborhood around y = 0 must satisfy the inequality $E(y)\leqslant E\left(y=0\right)$. In contrast, just after the first-order QPT, where E(y = 0) remains a local minimum, such a neighborhood around y = 0 satisfying the inequality $E(y)\leqslant E\left(y=0\right)$ must vanish. The distinguishing properties of the solutions can determine the first-order SRPT point. We now proceed to analyze the solution of $E(y)\leqslant E\left(y=0\right)$ below.
The inequality $E(y)\leqslant E\left(y=0\right)$ can be equivalently written as
$\begin{eqnarray}{\left[{\tilde{y}}^{2}-E\left(y=0\right)\right]}^{2}-2{f}_{1}\leqslant 2\sqrt{{f}_{1}^{2}-{f}_{2}}.\end{eqnarray}$
The solutions can be divided into two scenarios based on the sign of the left-hand side of the inequality. In the first scenario, the left-hand side of the inequality is negative, which is equivalent to $\begin{eqnarray}\tilde{g}\gt {\tilde{g}}_{-}=\sqrt{\left(1+\sqrt{2}/2\right)}\sqrt[4]{1+{\tilde{\varepsilon }}^{2}}.\end{eqnarray}$
In this scenario, the coupling is so strong that the system has entered the SRP deeply, making it independent of the QPT, as will be clarified once the QPT point is derived below. In the other scenario, the left-hand side of the inequality is nonnegative but smaller than the right-hand side when $\tilde{g}\leqslant {\tilde{g}}_{-}$, leading to $\begin{eqnarray}{y}^{2}{f}_{3}={y}^{2}\left({y}^{6}+{d}_{4}{y}^{4}+{d}_{2}{y}^{2}+{d}_{0}\right)\leqslant 0,\end{eqnarray}$
where the coefficients are defined as $\begin{eqnarray}{d}_{0}=16\left[{\left(1+{\tilde{\varepsilon }}^{2}\right)}^{3/2}-2{\tilde{g}}^{2}\right],\end{eqnarray}$
$\begin{eqnarray}{d}_{2}=4\sqrt{1+{\tilde{\varepsilon }}^{2}}\left(5\sqrt{1+{\tilde{\varepsilon }}^{2}}-8{\tilde{g}}^{2}\right),\end{eqnarray}$
$\begin{eqnarray}{d}_{4}=8\left(\sqrt{1+{\tilde{\varepsilon }}^{2}}-{\tilde{g}}^{2}\right).\end{eqnarray}$
Clearly, the factor y2 yields the trivial solution, while the nontrivial solution is given by the inequality f3(y2) ≤ 0. When the system parameter satisfies d0 ≤ 0, the real solution of f3(y2) ≤ 0 is a neighborhood around y = 0. Therefore, the second-order QPT point is derived as d0 = 0, which is identical to equation ( 11 ) based on the Landau's theory of phase transitions. On the other hand, when the system parameter gives d0 > 0, the real solution of f3(y2) ≤ 0 requires that the function f3(y2) has a zero or negative local minimum at a positive value ${y}_{\pm }^{2}=\left(-{d}_{4}+\sqrt{{d}_{4}^{2}-3{d}_{2}}\right)/3$. This scenario clearly corresponds to the first-order QPT, and the first-order QPT point ${\tilde{g}}_{c}^{1}$ is then given by a vanishing local minimum, i.e.,
$\begin{eqnarray}{\left.{f}_{3}({y}^{2}={y}_{\pm }^{2})\right|}_{\tilde{g}={\tilde{g}}_{c}^{1}}=0.\end{eqnarray}$
The equation may be solved numerically, and the first-order QPT point given by its solution is shown in figure 1 as the magenta line. This is consistent with the numerical result for the system with an effective size η = 1000, obtained by exact diagonalization, as shown in figure 1.4. Numerical verification via finite-frequency scaling analysis
In this section, we explore the critical behavior of the second-order SRPT and the TCP through finite-size scaling analysis. In the critical regime of a continuous phase transition, the finite-size scaling function of the physical observable O can be formulated as
$\begin{eqnarray}O\left(\eta ,\tilde{g}\right)={\left|1-\tilde{g}/{\tilde{g}}_{c}\right|}^{{\beta }_{O}}f\left(\left|1-\tilde{g}/{\tilde{g}}_{c}\right|{\eta }^{\nu }\right),\end{eqnarray}$
where η is the effective size of the system, βO is the critical exponent of the observable, and ν is the observable-independent correlation-length exponent.Based on the finite-size scaling function described above, we present the scaling behaviors of the order parameter ${n}_{0}=\langle {\psi }_{0}\left|{a}^{\dagger }a\right|{\psi }_{0}\rangle /\eta $ in figures 2(a) and (c), and the energy gap ε = E1 − E0 in figures 2(b) and (d), where En and ψn represent the eigenvalue and eigen-wavefunction of the system. Figures 2(a) and (b) show the critical behaviors of the second-order SRPT, where the curves for significantly varying effective sizes η collapse onto a single curve. The correlation-length exponent ν = 2/3, and the exponents ${\beta }_{{n}_{0}}=1$ (for the order parameter ) and βε = 1/2 (for the energy gap) match those of the QRM and the DM [7–10, 19, 20]. Moreover, the critical behaviors at the TCP with $\tilde{\varepsilon }=1/2$ are shown in figures 2(c) and (d). Although the critical exponent of the energy gap remains unchanged, the correlation-length (order-parameter) exponents change to νT = 1 and ${\beta }_{{n}_{0}}^{{\rm{T}}}=1/2$, which differ from those of the original SRPTs in the QRM and the DM, as well as from those of the TCP in their anisotropic extensions [21, 33]. Additionally, the insets display the power-law behaviors of the order parameter and the energy gap at the PT point, i.e., ${O}_{c}\propto {\eta }^{-{\gamma }_{O}}$. The fitted values of the exponent γO satisfy the scaling law γO = βOν quite well, further confirming the exponents extracted above.
Figure 2. Finite size scaling of the ground-state order parameter ${n}_{0}=\langle {\psi }_{0}\left|{a}^{\dagger }a\right|{\psi }_{0}\rangle /\eta $ (left panel) and the energy gap ε = E1 − E0 (right panel) for the second-order SRPT with $\tilde{\varepsilon }=0.3$ (upper panel) and at the TCP with $\tilde{\varepsilon }=0.5$ (lower panel). Additionally, the insets show the size η dependence of the order parameter nc (left panel) and the energy gap εc (right panel) at the corresponding critical point, with the slope of the linear fit presented in a log-log scale. |
Aside from the order parameter and the energy gap, the ground-state fidelity and its susceptibility can also detect the QPT, even without the detailed knowledge of the QPT. The ground-state susceptibility is defined as [34, 35]
$\begin{eqnarray}{\chi }_{F}(g)=\mathop{\mathrm{lim}}\limits_{\delta g\to 0}\frac{-2\,{\rm{ln}}\,F(g,\delta g)}{\delta {g}^{2}},\end{eqnarray}$
where the ground-state fidelity is the overlap of two ground states at close parameters, i.e., $F(g,\delta g)=\left|\langle {\psi }_{0}\left(g\right)| {\psi }_{0}\left(g+\delta g\right)\rangle \right|$. A continuous QPT can be characterized by the fidelity susceptibility, which exhibits a dramatic increase around the QPT and reaches a maximum shortly after the critical coupling [34–36]. As shown in figures 3(b) and (d), the maximum of the fidelity susceptibility ${\chi }_{F}^{\max }$ and its location deviation ${\tilde{g}}_{m}-{\tilde{g}}_{c}$ from the critical coupling both satisfy a nice power law at large effective size η, i.e., ${\chi }_{F}^{\max }\propto {\eta }^{\mu }$ with μ = 1/3 for the second-order SRPT and μT = 1 for the TCP, and $\left({\tilde{g}}_{m}-{\tilde{g}}_{c}\right)\propto {\eta }^{-\zeta }$ with ζ = 2/3 for the second-order SRPT and ζT = 1 for the TCP. Moreover, the critical behaviors of the fidelity susceptibility can also offer the correlation-length critical exponent, which should be identical to the value given above. Its finite-size scaling function may be denoted as [11–14, 37] $\begin{eqnarray}\frac{{\chi }_{F}^{\max }-{\chi }_{F}(g)}{{\chi }_{F}(g)}=f\left[\left(1-\tilde{g}/{\tilde{g}}_{m}\right){\eta }^{\nu }\right],\end{eqnarray}$
where ${\tilde{g}}_{m}$ represents the location of the maximum of the fidelity susceptibility. Figures 3(a) and (b) illustrate the critical behaviors of the second-order SRPT and at the TCP, where the curves for significantly varying effective sizes η collapse well onto a single curve. The corresponding correlation-length exponents are ν = 2/3 for the second-order SRPT and νT = 1 at the TCP. Notably, these exponents match those derived earlier. Based on this numerical analysis, it can be concluded that the TCP introduces a new criticality that differs from that of the ordinary SRPT.
Figure 3. Finite-size scaling of the ground-state fidelity susceptibility for the second-order SRPT (a) and at the TCP (c). Additionally, the power-law behaviors of the maximum fidelity susceptibility (left y-axis in blue) and its deviation from the critical coupling (right y-axis in orange) are shown in a log-log scale with a dual y-axis plot for the second-order SRPT (b) and at the TCP (d). |
5. Summary
This work explores the rich criticality in the two-qubit QRM with staggered qubit biases. Although the qubit biases typically make the model asymmetric, using identical positive and negative biases for the two qubits restores the parity symmetry of the generalized QRM, enabling symmetry breaking with strong qubit-cavity coupling. The mean-field theory can be applied to this model, and a rescaled Hamiltonian is expressed in terms of the new field quadratures. Using Landau's theory of phase transitions, we derive a TCP that separates the first- and second-order QPTs. The position of the TCP in the rescaled qubit–cavity coupling (${\tilde{g}}_{tc}$) and qubit bias (${\tilde{\varepsilon }}_{tc}$) plane is derived analytically : ${\tilde{g}}_{tc}=\frac{{5}^{\frac{3}{4}}}{4},{\tilde{\varepsilon }}_{tc}=\frac{1}{2}$. All phase boundaries in the phase diagram are analytically derived based on mean-field theory, including the nonperturbative first-order SRPT. The analytical phase diagram is numerically confirmed by finite-size scaling analysis of the macroscopic excitation of the photon field. The extracted critical exponents for the correlation length and order parameter at the TCP are νT = 1 and ${\beta }_{{n}_{0}}^{{\rm{T}}}=1/2$, respectively, different from ν = 2/3 and ${\beta }_{{n}_{0}}=1$ in the ordinary second-order SRPT. This implies that the TCP belongs to a new universality class. Further evidence of the rich phase diagram is provided by an analysis of the ground-state fidelity susceptibility, which is essentially independent of the order parameter. The critical coupling strength of the SRPTs and the correlation-length exponent for both the second-order SRPT and TCP are recovered independently. Since qubit bias can be easily manipulated in solid-state quantum devices, such as circuit QED [38, 39], this work provides a reliable theoretical framework for designing simple experimental platforms to explore rich QPTs.