Universal scaling and critical exponents of the anisotropic quantum Rabi model

TL;DR

This paper analyzes the phase transitions and critical behavior of the anisotropic quantum Rabi model, revealing universal properties and novel features across different coupling regimes relevant to light-matter systems.

Contribution

It provides a comprehensive phase diagram, scaling functions, and critical exponents for the anisotropic quantum Rabi model, demonstrating universality across anisotropy levels.

Findings

Universality class remains the same at finite anisotropy and isotropic limit.

Identification of superradiance-induced effective mass freezing.

Discontinuous scaling functions observed in the Jaynes-Cummings limit.

Abstract

We investigate first- and second-order quantum phase transitions of the anisotropic quantum Rabi model, in which the rotating- and counter-rotating terms are allowed to have different coupling strength. The model interpolates between two known limits with distinct universal properties. Through a combination of analytic and numerical approaches we extract the phase diagram, scaling functions, and critical exponents, which allows us to establish that the universality class at finite? anisotropy is the same as the isotropic limit. We also reveal other interesting features, including a superradiance-induced freezing of the effective mass and discontinuous scaling functions in the Jaynes-Cummings limit. Our findings are relevant in a variety of systems able to realize strong coupling between light and matter, such as circuit QED setups where a finite anisotropy appears quite naturally.

Tables

Table 1. Table 1: Critical point g c ~ ~ subscript 𝑔 𝑐 \tilde{g_{c}} , critical exponent ratio β / ν 𝛽 𝜈 \beta/\nu , and critical exponent 1 / ν 1 𝜈 1/\nu at different values of the anisotropy parameter λ 𝜆 \lambda . Both numerical and analytical results are presented. In the numerical calculations, the error for g ~ c subscript ~ 𝑔 𝑐 \tilde{g}_{c} is less than 0.0005 0.0005 0.0005 and the error for β / ν 𝛽 𝜈 \beta/\nu ( 1 / ν 1 𝜈 1/\nu ) is less than 0.02 0.02 0.02 .

$\lambda$	$0.1$	$0.2$	$0.3$	$0.4$	$0.5$	$0.6$	$0.7$	$0.8$	$0.9$	$1.0$	$1.5$	$2.0$	$5.0$
${\tilde{g}}_c$ (Analytical)	$20/11$	$5/3$	$20/13$	$10/7$	$3/2$	$4/3$	$20/17$	$10/9$	$20/19$	$1$	$4/5$	$2/3$	$1/3$
${\tilde{g}}_c$ (Numerical)	$1.18182$	$1.6667$	$1.5384$	$1.4286$	$1.4999$	$1.3334$	$1.1765$	$1.1111$	$1.0526$	$1.0002$	$0.7999$	$0.6666$	$0.3334$
$\beta /\nu$ (Analytical)	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$
$\beta /\nu$ (Numerical)	$0.34$	$0.33$	$0.34$	$0.32$	$0.33$	$0.33$	$0.35$	$0.33$	$0.32$	$0.33$	$0.34$	$0.33$	$0.33$
$1/\nu$ (Analytical)	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$	$1/3$
$1/\nu$ (Numerical)	$0.33$	$0.32$	$0.33$	$0.32$	$0.34$	$0.33$	$0.34$	$0.33$	$0.33$	$0.33$	$0.32$	$0.33$	$0.33$

Equations

$$H=\omega a^{†}a+\frac{\Omega }{2}{\sigma }_x+g[({\sigma }_{+}a+{\sigma }_{-}{a}^{†})+\lambda ({\sigma }_{+}{a}^{†}+{\sigma }_{-}a)],$$

$$H=\frac{p^2+x^2}{2\eta }+\frac{{\sigma }_x}{2}+\tilde{g}\left[\frac{1+\lambda }{\sqrt{8\eta }}{\sigma }_{z}x+\frac{1-\lambda }{\sqrt{8\eta }}{\sigma }_{y}p\right],$$

$$H_{eff}\simeq \frac{p^2+x^2}{2\eta }-{\tilde{g}}^2\frac{(1+\lambda {)}^2{x}^2+(1-\lambda {)}^2{p}^2}{8\eta }+\dots ,$$

$$\epsilon =\sqrt{(1-{\xi }^2)(1-{{\xi }^{\prime }}^2)},$$

$$\tilde{H}=\frac{{\tilde{p}}^2}{2{\eta }^2}+\frac{1}{2}{\tilde{x}}^2+\frac{1}{2}{\sigma }_x+\frac{\xi }{\sqrt{2}}{\sigma }_z\tilde{x}+\frac{{\xi }^{\prime }}{\sqrt{2\eta }}{\sigma }_y\tilde{p}.$$

$${\tilde{E}}_{\pm }=\frac{1}{2}\left({\tilde{x}}^2\pm \sqrt{1+2{\xi }^2{\tilde{x}}^2}\right),$$

$${\tilde{x}}_{0,\pm }=\pm \sqrt{\frac{{\xi }^2-{\xi }^{-2}}{2}}\theta (\xi -1),$$

$$H_{eff}\simeq \frac{p^2+x^2}{2\eta }-\sqrt{\frac{1}{4}+\frac{{\xi }^2{x}^2+{{\xi }^{\prime }}^2{p}^2-\xi {\xi }^{\prime }}{2\eta }}.$$

$${\tilde{H}}_{eff}\simeq -\frac{1}{2}+\frac{\xi {\xi }^{\prime }}{2\eta }+\frac{2\lambda }{(1+\lambda {)}^2}\frac{{\tilde{p}}^2}{{\eta }^2}+\frac{1-{\xi }^2}{2}{\tilde{x}}^2+\frac{{\xi }^4}{4\eta }{\tilde{x}}^4.$$

$$M_{\lambda }={\eta }^2\frac{(1+\lambda {)}^2}{4\lambda },$$

$$\left(-\frac{1}{2}\frac{{\partial }^2}{\partial u^2}-v{u}^2+\frac{u^4}{4}\right){\varphi }_0(u,v)=E_0(v){\varphi }_0(u,v),$$

$$E_G(\lambda )=-\frac{1}{2}+\frac{1}{2\eta }\frac{1-\lambda }{1+\lambda }+M_{\lambda }^{-2/3}{E}_0(t{M}_{\lambda }^{1/3})+\dots .$$

$$⟨x^{2n}⟩={\left(\eta M_{\lambda }^{-1/3}\right)}^n{X}_n\left(t{M}_{\lambda }^{1/3}\right),$$

$$⟨p^{2n}⟩={\left(\eta M_{\lambda }^{-1/3}\right)}^{-n}{P}_n\left(t{M}_{\lambda }^{1/3}\right),$$

$$\frac{(2n-1)!!}{2^n}\left[1+\sum _{q=0}^{\infty }\theta (\eta t-q)(D(n,q+1)-D(n,q))\right],$$

$$H_{eff}$$

$$P=|-⟩⟨-|,Q=|+⟩⟨+|,$$

$$H_0=\frac{1}{2}{\sigma }_x+\frac{p^2+x^2}{2\eta },V_{od}=\tilde{g}\left(\frac{1+\lambda }{\sqrt{8\eta }}x{\sigma }_z+\frac{1-\lambda }{\sqrt{8\eta }}p{\sigma }_y\right).$$

$$H_{eff}^{(2)}=\frac{1}{2}P[S_1,V_{od}]P=\frac{1}{2}P[L_0^{-1}{V}_{od},V_{od}]P,$$

$$S_1$$

$$=-i\tilde{g}\underset{s\to 0}{\lim }{\int }_0^{\infty }{e}^{-st}\left[\frac{1+\lambda }{\sqrt{8\eta }}x(t){\sigma }_z(t)+\frac{1-\lambda }{\sqrt{8\eta }}p(t){\sigma }_y(t)\right]dt$$

$$\begin{array}[]{l}x(t)=x\cos(t/\eta)+p\sin(t/\eta),\vspace{0.1cm}\\ p(t)=p\cos(t/\eta)-x\sin(t/\eta),\end{array}\qquad\quad\begin{array}[]{l}\sigma_{z}(t)=\sigma_{z}\cos t+\sigma_{y}\sin t,\vspace{0.1cm}\\ \sigma_{y}(t)=\sigma_{y}\cos t-\sigma_{z}\sin t.\end{array}$$

$$H_{eff}^{(4)}=-\frac{1}{6}P\left[L_0^{-1}\left[S_1,\left[S_1,\left[S_1,H_0\right]\right]\right],V_{od}\right]P-\frac{1}{24}P\left[S_1,\left[S_1,\left[S_1,V_{od}\right]\right]\right]P.$$

$$S_1\simeq -i\tilde{g}\left(\frac{1+\lambda }{\sqrt{8\eta }}x{\sigma }_y-\frac{1-\lambda }{\sqrt{8\eta }}p{\sigma }_z\right),$$

$$H_{eff}\simeq -\frac{1}{2}+\frac{p^2+x^2}{2\eta }-V^2+V^4+\dots$$

$$V^2=⟨-|V_{\mathrm{od}}^2|-⟩=\frac{{\tilde{g}}^2}{8\eta }\left[(1+\lambda {)}^2{x}^2+(1-\lambda {)}^2{p}^2-(1-{\lambda }^2)\right].$$

$$⟨m|H_{eff}^{(2)}|m^{\prime }⟩=\frac{1}{2}\sum _l⟨m|V_{od}|l⟩⟨l|V_{od}|m^{\prime }⟩\left[\frac{1}{E_{ml}}+\frac{1}{E_{m^{\prime }l}}\right],$$

$$⟨m|H_{eff}^{(2)}|m^{\prime }⟩\simeq -\sum _l⟨m|V_{od}|l⟩⟨l|V_{od}|m^{\prime }⟩=-⟨m|V_{od}^2|m^{\prime },⟩$$

$$⟨m|H_{eff}^{(4)}|m^{\prime }⟩=-\frac{1}{24}\sum _{l,l^{\prime },m^{\prime \prime }}$$

$$+\frac{4}{E_{m{l}^{\prime }}{E}_{m^{\prime \prime }l}{E}_{ml}}+\frac{4}{E_{m{l}^{\prime }}{E}_{m^{\prime \prime }l}{E}_{m^{\prime \prime }{l}^{\prime }}}+\frac{4}{E_{m^{\prime }l}{E}_{m^{\prime \prime }{l}^{\prime }}{E}_{m^{\prime }{l}^{\prime }}}+\frac{4}{E_{m^{\prime }l}{E}_{m^{\prime \prime }{l}^{\prime }}{E}_{m^{\prime \prime }l}}$$

$$-\frac{1}{E_{m^{\prime \prime }l}{E}_{m^{\prime \prime }{l}^{\prime }}{E}_{ml}}-\frac{1}{E_{m^{\prime \prime }l}{E}_{m^{\prime \prime }{l}^{\prime }}{E}_{m^{\prime }{l}^{\prime }}}-\frac{3}{E_{ml}{E}_{m^{\prime }{l}^{\prime }}{E}_{m^{\prime \prime }l}}-\frac{3}{E_{ml}{E}_{m^{\prime }{l}^{\prime }}{E}_{m^{\prime \prime }{l}^{\prime }}}).$$

$$⟨m|H_{eff}^{(2n)}|m^{\prime }⟩\simeq {\alpha }_n\sum ⟨m|V_{od}|l^{(n)}⟩⟨l^{(n)}|V_{od}|m^{(n)}⟩\dots ⟨m^{\prime \prime }|V_{od}|l^{\prime }⟩⟨l^{\prime }|V_{od}|m^{\prime }⟩={\alpha }_n⟨m|V_{\mathrm{od}}^{2n}|m^{\prime }⟩.$$

Full text

Introduction.– While critical phenomena are traditionally associated with collective behavior in the thermodynamic limit, quantum phase transitions in systems with few degrees of freedom were recently brought to prominence Hwang et al. (2015); Hwang and Plenio (2016). As it turns out, the topic is of great relevance for ongoing efforts on enhancing and engineering light-matter interactions. By achieving the strong Wallraff et al. (2004); Devoret et al. (2007), ultrastrong Bourassa et al. (2009); Niemczyk et al. (2010); Peropadre et al. (2010); Ridolfo et al. (2012); Todorov and Sirtori (2014); Baust et al. (2016); Forn-Diaz et al. (2016), and even deep strong coupling regime Casanova et al. (2010); De Liberato (2014); Fornd az et al. (2017); Yoshihara et al. (2017); Chen et al. (2016), atomic and solid-state resonances are able to induce profound modifications of the photon fields they interact with. The quantum Rabi model (QRM), describing a two-level system coupled to a single electromagnetic mode, represents the simplest realization of such light-matter interactions Rabi (1936, 1937). Thus, it has served as a paradigmatic example to explore this kind of strong-coupling phenomena, and has received renewed attention in recent years Forn-Díaz et al. (2010); Braak (2011); Solano (2011); Bakemeier et al. (2012); Chen et al. (2012); Ashhab (2013); Wolf et al. (2013); Ying et al. (2015); Hwang et al. (2015); Braak et al. (2016).

Surprisingly, an analytic solution of the QRM was only found recently, and has also motivated proposing a novel operational criterion of integrability Braak (2011). More directly related to the present study are several recent analyses on the dependence of QRM ground-state properties on the coupling strength Ashhab (2013); Ying et al. (2015); Hwang et al. (2015). An ansatz for the ground state based on the polaron and antipolaron concept was introduced in Ref. Ying et al. (2015), which has presented a phase diagram where the quadpolaron dominates in the weak coupling regime and the bipolaron dominates at strong coupling. The crossover between these two types of ground state becomes sharper by reducing the bosonic frequency Ashhab (2013); Ying et al. (2015), namely, in the classical oscillator limit Bakemeier et al. (2012). It was proved that this behavior indeed reflects the existence of a true quantum phase transition (QPT), and the static and dynamical properties of the critical point were studied in detail Hwang et al. (2015). Furthermore, it was later found that also the Jaynes-Cummings (JC) model Jaynes and Cummings (1963) exhibits a second-order QPT in the same limit Hwang and Plenio (2016).

These findings have motivated us to study the physics of the QPT in the anisotropic QRM Xie et al. (2014), which includes the two known limits as special cases. The model reads

[TABLE]

where $\hbar=1$, $a^{\dagger}$($a$) is the creation (annihilation) operator of the bosonic field oscillator, and $\boldsymbol{\sigma}$ are the Pauli matrices with $\sigma_{\pm}=\frac{1}{2}(\sigma_{z}\mp i\sigma_{y})$. As known, the anisotropic QRM is relevant in a variety of systems including quantum well with spin-orbit coupling Schliemann et al. (2002); Wang et al. (2016) and circuit QED, where strong interactions were realized Wallraff et al. (2004); Devoret et al. (2007); Niemczyk et al. (2010); Forn-Díaz et al. (2010); Forn-Diaz et al. (2016); Fornd az et al. (2017); Yoshihara et al. (2017); Chen et al. (2016). In this case, the asymmetry between rotating- and counter-rotating terms is often the typical scenario. For example, in the proposal of Ref. Baksic and Ciuti (2014) the degree of anisotropy is simply given by the relative strength of the inductive and capacitive couplings to the cavity. Also in a setup of two coupled SQUIDs, which was not specifically designed to implement Eq. (1), measurements of the Bloch-Siegert shift are in good agreement with an intermediate value $\lambda\simeq 0.5$ Forn-Díaz et al. (2010); Xie et al. (2014). The anisotropic QRM could also be implemented with trapped-ions Puebla et al. (2016) and is equivalent to a Hamiltonian of spin-orbit coupled electrons in semiconductors Xie et al. (2014) (possibly emulated by fermionic gases).

From the theoretical perspective, Eq. (1) allows us to address systematically the role of the counter-rotating terms in the QPT, since the parameter $\lambda$ can interpolate between the QRM ($\lambda=1$) and JC model ($\lambda=0$). The intermediate case is of special interest because the two known limits have different critical exponents and types of broken symmetry phase Hwang et al. (2015); Hwang and Plenio (2016). Thus, the question about intermediate values of $\lambda$ arises very naturally. In the following, we develop a unified treatment which allows us to identify the appropriate scaling parameters and establish the complete phase diagram in the presence of anisotropy. We find that the QPT at fixed $\lambda$ exhibits a remarkable degree of universality: it occurs between a normal and superradiant (x- or p-type) phase, with critical exponents and scaling functions which are independent of $\lambda\neq 0$. Therefore, the generic realization of the QPT belongs to a well-defined universality class, which is the same of $\lambda=1$ limit. On the other hand, by extending previous work on the JC model Hwang and Plenio (2016), we expose the singular character of the QPT at $\lambda=0$.

Phase diagram.– To understand the main features of the phase diagram, we first rewrite Eq. (1) by using dimensionless coordinate and momentum operators, $x=\frac{1}{\sqrt{2}}(a^{\dagger}+a)$ and $p=\frac{i}{\sqrt{2}}(a^{\dagger}-a)$. Within a constant term:

[TABLE]

where $\eta=\Omega/\omega$ and $\tilde{g}=g/g_{c,0}$, with $g_{c,0}=\sqrt{\omega\Omega}/2$ Ashhab (2013). To simplify the notation, we set $\Omega=1$ in Eq. (2) and hereafter. First of all, we note that the second term becomes dominant in the $\eta\to\infty$ limit of interest, thus the relevant low-energy states have $\langle\sigma_{x}\rangle\simeq-1$. Within this subspace, the ground-state is determined by the competition between the first term (a conventional oscillator) and the last term (the coupling between the bosonic mode and the two-level system). The coupling term has a larger prefactor, proportional to $\eta^{-1/2}$, but is off-diagonal in $\sigma_{x}$. Treating it with second-order perturbation theory gives:

[TABLE]

which shows that the oscillator term is dominant in the weak coupling regime, corresponding to the normal phase. On the other hand, at sufficiently large $\tilde{g}$ the coupling term will dominate and the Hamiltonian in Eq. (3) becomes unbounded. To this order of approximation, $H_{eff}$ implies divergent values of $\langle x^{2}\rangle,\langle p^{2}\rangle$, which in turn signals the onset of a superradiant phase. By introducing $\xi=\tilde{g}(1+\lambda)/2$ and $\xi^{\prime}=\tilde{g}(1-\lambda)/2$, the gap $\varepsilon/\eta$ of the normal phase is simply given by:

[TABLE]

which becomes zero at $\tilde{g}_{c}=\frac{2}{1+|\lambda|}$. This phase boundary is plotted in Fig. 1(a) as a solid (blue) line and recovers the known results for the isotropic QRM Ashhab (2013); Hwang et al. (2015) and JC model Hwang and Plenio (2016) at $\lambda=1$ and $\lambda=0$, respectively. The onset of instability in Eq. (3) is due to the $x^{2}$ ($p^{2}$) terms when $\lambda>0$ ($\lambda<0$), suggesting that the superradiant phase should be divided into two regimes according to the sign of the anisotropy parameter $\lambda$. This observation is reflected by the vertical (red) phase boundary of Fig. 1(a). All these conclusions are further confirmed by the numerical analysis (see Fig. 1) and will be more rigorously justified and extended in the rest of the paper.

Classical oscillator limit.– For $\lambda>0$, the $\sigma_{z}x$ coupling plays a dominant role and, to proceed, it is helpful to rescale the coordinate by $\tilde{x}=x/\sqrt{\eta}$. Correspondingly, $\tilde{p}=-i\partial/\partial\tilde{x}$. Thus, $H$ becomes

[TABLE]

By taking the classical oscillator limit Bakemeier et al. (2012), i.e., $\eta\rightarrow\infty$ (but keeping $\tilde{g}$ finite), one can drop the first and last terms of Eq. (5). The remaining Hamiltonian is readily diagonalized and has eigenvalues

[TABLE]

where $\tilde{x}$ is now a classical coordinate. The lower branch $\tilde{E}_{-}$ has the standard behaviour of the Landau potential across a continuous phase transition, where $\xi$ plays a role analogous to the inverse temperature: $\tilde{E}_{-}$ has one minimum around $\tilde{x}=0$ when $\xi<1$ and two minima at finite $\tilde{x}$ when $\xi>1$. The order parameter is:

[TABLE]

where $\theta(x)$ is the step function and the critical value $\xi_{c}=1$ is in agreement with the $\lambda>0$ side of Fig. 1. The ground state energy $\tilde{E}_{gs}=-\frac{1}{2}-\frac{1}{4}\left(\xi^{2}+\xi^{-2}-2\right)\theta(\xi-1)$ is easily obtained from Eqs (5) and (7) and indicates a second-order QPT at $\xi_{c}$, since the first-order derivative $\partial\tilde{E}_{gs}/\partial\xi$ is continuous but the second-order derivative $\partial^{2}\tilde{E}_{gs}\partial\xi^{2}$ is not [see Fig. 1(c)]. We also emphasize that, although the anisotropy parameter $\lambda$ obviously influences the critical coupling strength $\tilde{g}_{c}$ as well as $\tilde{x}_{0,\pm}$ and $\tilde{E}_{gs}$, the functional dependence of these physical quantities becomes universal – in the sense of being independent of $\lambda$ – once the problem is formulated in terms of the rescaled coupling $\xi$.

The other physical quantities have also been calculated Sup . In particular, $\langle\tilde{x}^{2}\rangle=\tilde{x}_{0,\pm}^{2}=(\xi^{2}-\xi^{-2})/2$ and the correlation function $\langle\sigma_{z}\tilde{x}\rangle=-(\xi-\xi^{-3})/\sqrt{2}$ have finite values as the system enters into the superradiant phase. On the other hand, $\langle\tilde{p}^{2}\rangle/\eta^{2}=\frac{1+\lambda}{4\eta\sqrt{\lambda}}\sqrt{1-\xi^{-4}}$ and $\langle\tilde{p}\sigma_{y}\rangle/\eta=\frac{\xi^{-3}}{\sqrt{2}\eta}\left(1-\frac{1-\lambda}{2\sqrt{\lambda}}\sqrt{\xi^{4}-1}\right)$ tend to zero as $\eta\rightarrow\infty$, which is compatible with dropping terms involving the momentum in the classical oscillator limit.

The case $\lambda<0$ can be treated in a similar way, or by using the fact that Eq. (2) remains unchanged under the following mapping $\mathscr{F}=\mathscr{F}_{1}\otimes\mathscr{F}_{2}$, which contains a $Z_{2}$ symmetric mapping $\mathscr{F}_{1}$ : $\{H(\lambda)\rightarrow H(-\lambda)\}$, and a unitary transformation $\mathscr{F}_{2}$ : $\{H\rightarrow U^{{\dagger}}HU\}$, where $U=e^{-i\frac{\pi}{2}a^{{\dagger}}a}\otimes e^{-i\frac{\pi}{4}\sigma_{x}}$. The Hamiltonian in Eq. (2) is a fixed point of the functional $\mathscr{F}$. Therefore, the properties we have discussed so far are readily translated to the $\lambda<0$ side of the phase diagram. In particular, the critical value $\xi^{\prime}=1$ gives the phase boundary of Fig. 1 for $\lambda<0$ (blue solid line). Furthermore, as made clear by the exact mapping, the superradiant phase at $\lambda<0$ is caused by the momentum $p$ rather than the position $x$. Thus, one gets two types of superradiant phases (x-type and p-type) and a sharp jump of both $\langle x^{2}\rangle$ and $\langle p^{2}\rangle$ at $\lambda=0$. One can also check that the first derivative of $\tilde{E}_{gs}$ is discontinuous at $\lambda=0$, since $\partial\tilde{E}_{gs}/\partial\lambda|_{\lambda\rightarrow 0^{\pm}}=\mp\left(\xi-1/\xi^{3}\right)\tilde{g}/4$ [see Fig. 1(d)]. This indicates the existence in the superradiant regime of a first-order QPT dependent on the sign of $\lambda$. We stress that this new transition is only revealed by considering the complete phase diagram of the QRM, extended to include anisotropy.

Effective Hamilltonians.– We would like next to address the critical scaling at the second-order QPT. It begins with our analytic approach which leads to the following effective Hamiltonian Sup :

[TABLE]

The above expression was obtained by performing an exact resummation of leading perturbative terms Sup . As we will see shortly, $H_{eff}$ is fully consistent with our previous discussion and, in fact, represents a generalization of Eq. (3) and the classical oscillator limit.

Firstly, we note that the order parameter in the superradiant phase ($\langle x^{2}\rangle$ or $\langle p^{2}\rangle$) is $\propto\eta$, thus an expansion of the square root in Eq. (8) becomes justified. Focusing on the x-type phase, we can neglect to first approximation $p^{2}$ and $\xi\xi^{\prime}$, to recover the classical potential $\tilde{E}_{-}$ of Eq. (6). If higher orders in $\eta^{-1}$ are considered we obtain the following approximation (in rescaled coordinates), which is suitable to characterize the scaling properties at $\xi\simeq 1$:

[TABLE]

Here, the quartic potential is simply the small-x expansion of $\tilde{E}_{-}$ [see Eq. (6)] and is valid for $\langle\tilde{x}^{2}\rangle\ll\eta|1-\xi|$ (i.e., $\xi\to 1$). The kinetic term was also derived under the asumption $\xi=1$. Since $\langle x^{2}\rangle$, $\langle p^{2}\rangle$ are small in the critical regime, low-order perturbation theory becomes accurate and Eq. (9) can be confirmed by a direct calculation. To do that, we have chosen $H_{0}=\sigma_{x}/2+(p^{2}+x^{2})/(2\eta)$ as unperturbed Hamiltonian and $V_{od}=(\xi x\sigma_{z}+\xi^{\prime}p\sigma_{y})/\sqrt{2\eta}$ as off-diagonal perturbation and computed the Schrieffer-Wolff transformation up to fourth-order Sup .

Equation (9) also defines the effective mass:

[TABLE]

which plays a central role for the critical scaling and the stability of the superradiant phase. To extend Eq. (10) beyond $\xi\simeq 1$, one can derive from Eq. (8) a kinetic term of the form $\frac{\tilde{p}^{2}}{2\eta^{2}}\left(1-{\xi^{\prime}}^{2}/\sqrt{1+2\xi^{2}\tilde{x}_{0,\pm}^{2}}\right)$, which is in agreement with Eq. (3) if $\xi<1$ (i.e., $\tilde{x}_{0,\pm}=0$). This expression implies that fluctuations of $p$ are promoted by the interaction in the normal regime, where a larger $\tilde{g}$ enhances the effective mass. In the superradiant phase, however, the second-order terms of Eq. (3) give an incorrect result and, in particular, naively predict a phase transition at $\xi^{\prime}=1$. This instability is prevented by the finite value of $\tilde{x}_{0,\pm}$, which takes into account non-perturbative effects beyond second-order. Using Eq. (7) leads to the remarkable result of a constant effective mass in the superradiant phase: Although the mass enhancement is actually due to the interaction, it becomes independent on the coupling strength for $\xi>1$ and only reflects the interaction anisotropy. In other words, the formation of the x-type superradiant phase freezes the effective mass renormalization at the $\xi=1$ value.

As a final remark on Eq. (10) we note that in the isotropic QRM there is no renormalization effect ($M_{\lambda=1}=\eta^{2}$). Thus, the peculiar interplay of superradiance and effective mass renormalization between the two quadratures is specifically related to the intermediate values of $\lambda$ considered here.

Finite-$\eta$ scaling.– We are now ready to discuss the scaling properties, based on Eq. (9) and the definition of $M_{\lambda}$. To reveal the universal properties of the phase transition, it is appropriate to introduce the scaling variables $u=\tilde{x}M_{\lambda}^{1/6}$ and $v=tM_{\lambda}^{1/3}$ (where $t=\xi-1$). Although the ground state wavefunction $\phi_{0}$ depends in general on four variables ($t,~{}x,~{}\eta,$ and $\lambda$), it is described by a simple equation in terms of $u,v$:

[TABLE]

where $E_{0}(v)$ gives the ground-state energy:

[TABLE]

The scaling law obeyed by a certain observable is easily derived from the scaling form $\phi_{0}(\tilde{x}M_{\lambda}^{1/6},tM_{\lambda}^{1/3})$ of the ground state. In the important case of $x^{2n},p^{2n}$ we obtain:

[TABLE]

where the universal functions $X_{n}(v),P_{n}(v)$ are given by expectation values over $\phi_{0}(u,v)$ and can be readily evaluated Sup . We have confirmed the validity of our treatment by direct numerical solution of $H$ [see Eq. (1)]. Figure 2 shows that the numerical values of $\langle x^{2}\rangle$ and $\langle p^{2}\rangle$ at large $\eta$ and different values of $\lambda$ all collapse into a single curve when appropriately scaled. The two numerical scaling functions agree with the $X_{1}(v),P_{1}(v)$ obtained from Eq. (11). We thus conclude that the presence of anisotropy does not modify either the critical exponents or the scaling behavior and identify the whole second-order phase transition line as belonging to the same universality class Larson and Irish (2016). On the other hand, the JC model ($\lambda=0$) is a special case which will be discussed next.

The JC critical line.– Due to the presence of the first-order transition at $\lambda=0$, the scaling functions have an abrupt change on that critical line. Since the exact mapping between opposite values of $\lambda$ interchanges the roles of $x$ and $p$, the right-hand sides of Eqs. (13) and (14) must be switched when $\lambda<0$ (together with the changes $M_{\lambda}\to M_{-\lambda}$ and $t\to\xi^{\prime}-1$). At $\lambda=0$ all the quadratures are equivalent since $\xi=\xi^{\prime}$ and Eq. (8) becomes a function of $p^{2}+x^{2}$. Therefore, we recover the Mexican-hat potential in phase space discussed in Ref. Hwang and Plenio (2016), together with the related spontaneous breaking of $U(1)$ symmetry and gapless Goldstone mode.

The anomalous nature of the JC model in the phase diagram is reflected by other critical properties. For example, at $\lambda=0$ the critical exponent $\alpha$ for the excitation gap is different from the rest of the phase diagram ($\alpha=1$ instead of $1/2$). The different value of $\alpha$ in the QRM and JC model was noted already in Ref. Hwang and Plenio (2016), but our analysis is much more general as it extends the $\alpha=1/2$ exponent to all values $\lambda\neq 0$. Furthermore, while our Eq. (7) is universally applicable for $\lambda>0$, we find that the corresponding result at $\lambda=0$ is half of Eq. (7). These two properties are illustrated by Fig. 3(a), where the energy gap ($\xi<1$) and the order parameter ($\xi>1$) are plotted.

Finally, we have also considered the scaling functions of the JC model. The expectation values of $\langle x^{2n}\rangle=\langle p^{2n}\rangle$ in the $\eta\to\infty$ limit are given by:

[TABLE]

where $D(n,q)$ are the Delannoy numbers Sup . The case $n=1$ is shown in Fig. 3(b) and should be compared to Fig. 2. We see that $\eta$ plays the role of $M_{\lambda}^{1/3}$ (which is not defined at $\lambda=0$) and that the functional dependence is completely different from the rest of the phase diagram. In particular, the scaling function has a discontinuous nature. This reflects the fact that the QPT of the JC model is given by a succession of level crossings Hwang and Plenio (2016); Sup (as seen, the spacing is $\Delta\xi\simeq 1/\eta$). Such scaling functions highlight again the singular nature of the $\lambda=0$ line within the phase diagram.

Conclusion.– We have characterized the QPTs of the QRM as a function of coupling strength and anisotropy and established the universal character of the second-order phase transition for $\lambda\neq 0$. A first-order critical line at $\lambda=0$ separates two types of superradiant phases. Besides universality, we have found other interesting features such as the freezing of the effective mass (induced by the broken-symmetry in $x$) and the discontinuous scaling functions of the JC model. Our results emphasize the critical role played by counter-rotating terms, whose current experimental relevance is due to the rapid progress in enhancing light-matter interactions Wallraff et al. (2004); Devoret et al. (2007); Bourassa et al. (2009); Niemczyk et al. (2010); Peropadre et al. (2010); Ridolfo et al. (2012); Todorov and Sirtori (2014); Baust et al. (2016); Forn-Diaz et al. (2016); Casanova et al. (2010); De Liberato (2014); Fornd az et al. (2017); Yoshihara et al. (2017); Chen et al. (2016); Gröblacher et al. (2009); Günter et al. (2009); Niemczyk et al. (2010). In particular, the exposed singularity of the JC limit implies that even tiny counter-rotating terms lead to dramatic changes of the scaling behavior. Since a finite anisotropy is quite natural Forn-Díaz et al. (2010); Xie et al. (2014), the universal scaling features discussed here may be tested by practical implementations. A study of the universality from the point of view of critical dynamics would be another interesting extension of this work. Finally, we note that the superradiant phase leads to a strongly squeezed ground state of light, which has potential value for metrology and enhanced sensing applications.

Acknowledgment.– This work is supported by NSFC (Grants No. 11604009, No. 11574025, No. 11325417, No. 11674139, No. 1121403), NSAF (Grant No. U1530401). SC acknowledges support from the National Key Research and Development Program of China (Grant No. 2016YFA0301200). We also acknowledge the computational support in Beijing Computational Science Research Center (CSRC).

Supplemental Material for “Universal scaling and critical exponents of the anisotropic quantum Rabi model”

Appendix A Effective Hamiltonians

We discuss here a perturbative treatment of the anisotropic QRM [see Eq.(2) in the main text]. We first consider the Sch rieffer-Wolff (SW) transformation, which provides a general method to derive low-energy effective Hamiltonians. The lower orders are relatively straightforward to evaluate and results up to fourth-order are given below, together with some details of the derivation. To extend such low-order SW Hamiltonian, an exact resummation of the leading contributions can be performed, which is described later in this Section.

A.1 SW effective Hamiltonian

For convenience, we give immediately the final result of the SW transformation (with $\Omega=1$):

[TABLE]

where the first line is the unperturbed Hamiltonian in the $\sigma_{x}=-1$ subspace. The second line is the full second-order result and the third line is the leading contribution from the fourth-order term (odd terms are all identically zero). It is easily checked that the first two lines recover Eq. (3) of the main text, i.e., omitting unnecessary constants and considering the limit of large $\eta$.

In discussing the derivation of Eq. (16), we refer to the brief summary of the SW transformation given in Appendix A of Ref. SC . In particular, we adopt the same notation. To apply the SW transformation, we define a low-energy subspace by introducing the following projectors:

[TABLE]

where $|\pm\rangle$ are eigenstates of $\sigma_{x}$. We partition the anisotropic QRM into an unperturbed Hamiltonian and off-diagonal perturbation as follows:

[TABLE]

In general, the SW transformation is formulated in the presence of a diagonal perturabtion $V_{d}$ which, however, is zero in our case. This simplifies things and, in particular, the first-order correction $H_{eff}^{(1)}=PV_{d}P=0$ is absent.

The second-order term is given as follows, in terms of the $L_{0}$ superoperator (defined by $L_{0}A=[H_{0},A]$):

[TABLE]

and can be evaluated by a straightforward calculation. The following quantity is useful:

[TABLE]

where

[TABLE]

After performing elementary time integrations, Eqs. (19) and (A.1) yield the second line of Eq. (16).

Since the third-order term $H_{eff}^{(3)}=\frac{1}{2}P[L_{0}^{-1}[S_{1},V_{d}],V_{od}]P$ is trivially zero, we focus on the fourth-order terms:

[TABLE]

which can be computed in a similar way to the second-order term. However, the full result is more cumbersome and is not given here. Instead, we concentrate on an approximation which useful in the limit of large $\eta$ and can be extended to higher orders (see the next Section). As seen in Eq. (21), the time dependence of $x(t),p(t)$ is much slower than $\sigma_{y,z}(t)$. Thus, in evaluating $S_{1}$ to leading order, we can substitute the $t=0$ value $x(t),p(t)\to x,p$ and obtain:

[TABLE]

which allows us to directly compute the second term of Eq. (22). In the first term, there is an additional time integral induced by $L_{0}^{-1}$. Performing the same approximation on $x(t),p(t)$ discussed above, we obtain the last line of Eq. (16) as final result.

A.2 Infinite resummation of the leading-order terms

By considering Eq. (16), we notice that we can approximate the second-order term at large $\eta$ and obtain a simple form for $H_{eff}$:

[TABLE]

where:

[TABLE]

The origin of the series in power of $V^{2}$, appearing in Eq. (24), is not difficult to understand if we consider the matrix elements of the perturbative corrections. For the second-order term:

[TABLE]

where $m,m^{\prime},l$ label eigenstates of $H_{0}$ with energies $E_{l},E_{m^{\prime}},E_{l}$. For the energy denominators, we used the notation $E_{ml}=E_{m}-E_{l}$. Since $m,m^{\prime}\in P$ and $l\in Q$, the two energy denominators are both equal to $-1+O(1/\eta)$ and, neglecting small corrections, we have:

[TABLE]

implying that $H_{eff}^{(2)}\simeq-V_{od}^{2}$ in the low-energy subspace. The effective Hamiltonian entering Eq. (24) is \langle-|H_{eff}^{(2)}|-\rangle\, which leads to the $-V^{2}$ correction in Eq. (24).

Similarly, we can consider the formula for the the fourth-order matrix elements Winkler_book :

[TABLE]

This again can be simplified by approximating all the energy denominators as $-1$, to yield the last term of Eq. (24). It should be clear that, neglecting the $O(1/\eta)$ corrections of the denominators, the generic form of the $n$-th order is:

[TABLE]

where the coefficient $\alpha_{n}$ depends on the detailed form of the perturbation theory formula at order $n$. As a consequence, we have:

[TABLE]

To find the values of $\alpha_{n}$, the simplest way is to consider the auxiliary problem with $H_{0}=\frac{1}{2}\sigma_{z}$ and $V_{od}=\epsilon\sigma_{x}$. The advantage of this Hamiltonian is that, in applying the perturbation theory formulas, the denominators are exactly equal to $-1$. The exact energy of the ground state is $-\sqrt{1/4+\epsilon^{2}}$ which can be expanded in Taylor series and allows to extract $\alpha_{n}=(-2)^{n-1}(2n-3)!!/n!$. This expression is in agreement with the first few values $\alpha_{0}=1/2$, $\alpha_{1}=1$ and $\alpha_{2}=-1$, known from Eq. (24).

Actually, the ground state energy of $\frac{1}{2}\sigma_{z}+\ \epsilon\sigma_{x}$ also indicates that the closed form of the series in Eq. (29) is simply $\sqrt{1/4+V^{2}}$. Thus, we obtain the final result of this section:

[TABLE]

Appendix B Expectation values in the superradiant phase

In this section we compute several relevant expectation values in the $\eta\to\infty$ limit. In particular, we consider $\langle x^{2}\rangle$, $\langle p^{2}\rangle$, $\langle x\sigma_{z}\rangle$, and $\langle p\sigma_{y}\rangle$, where the average is over the ground state of the anisotropic QRM. The behavior of these quantities at large $\eta$ justifies our discussion in the main text in terms of a “Classical oscillator limit”.

Focusing on the x-type superradiant phase, we can make use of the effective Hamiltonian Eq. (30) to find the ground state. Following the discussion after Eq. (10) of the main text, we have that $H_{eff}$ can be written as follows in the superradiant phase:

[TABLE]

where the effective mass corresponds to Eq. (10) of the main text, and the potential corresponds to $\tilde{E}_{-}$ [see Eq. (6) of the main text]. In the superradiant phase, the potential has two minima which become equivalent when $\eta\to\infty$. For simplicity we focus on the one at $x_{0}\simeq\sqrt{(\xi^{2}+\xi^{-2})\eta/2}$. Performing an harmonic approximation (valid for $\eta\to\infty$) we have:

[TABLE]

which allows us to immediately write the ground state $|\psi_{0}\rangle$ for $H_{eff}$. In fact, remembering that $H_{eff}$ is valid in the $\sigma_{x}=-1$ subspace, the ground state is actually $|\Phi_{0}\rangle=|\psi_{0}\rangle|-\rangle$.

Given an operateor $A$, the expectation value $\langle\Phi_{0}|A|\Phi_{0}\rangle$ is easily computed. However, one should pay attention that $H_{eff}$ is obtained from the original $H$ after a unitary rotation $e^{S}$, thus the ground-state expectation value for the Rabi model is actually given by:

[TABLE]

We then must consider first the derivation of suitable expressions for $S$ and $e^{S}$.

As it is standard in the Schrieffer-Wolff transformation, we expand $S$ in a series $S=\sum_{n}S_{n}$, where explicit expressions for the low-order $S_{n}$ are readily available (see, e.g., Ref. Winkler_book ). Due to the structure of $S_{n}$, one can use an argument similar to the one we have developed for $H_{eff}$. Within the same approximation of Eq. (30), it easy to see that:

[TABLE]

while $S_{2n}=0$. To determine the coefficients $\beta_{n}$, we can consider again the auxiliary Hamiltonian $\frac{1}{2}\sigma_{z}+\ \epsilon\sigma_{x}$, for which $\langle\uparrow|S|\downarrow\rangle=1/2\arctan(2\epsilon)$. This leads to: $\langle l|S|m\rangle\simeq 1/2\langle l|\arctan(2V_{od})|m\rangle=1/2\langle l|\sigma_{x}\arctan(2V_{od})|m\rangle$. The last step, i.e., inserting the $\sigma_{x}$ operator, is justified because $|l\rangle$ is a $+1$ eigenstate of $\sigma_{x}$. This step allows us to write the matrix elements in terms of an anti-hermitian operator, as appropriate for $S$. Based on this we conclude:

[TABLE]

which is consistent with $S_{1}$ given in Eq. (23).

To proceed further, we write $S$ as a leading-order term plus a small $O(\eta^{-1/2})$ correction:

[TABLE]

This expression allows us to approximate the unitary transformation to the same order. We obtain:

[TABLE]

which can be directly used to compute the leading-order of the various expectation values.

In fact, in several cases the approximation $e^{S}\simeq e^{-i\frac{\sigma_{y}}{2}\arctan\sqrt{\xi^{4}-1}}$ is sufficient. This includes the evaluation of $\langle x^{2}\rangle$ and $\langle p^{2}\rangle$. Since the leading order of $e^{S}$ commutes with $x^{2},~{}p^{2}$, we have:

[TABLE]

where we have made use of Eq. (32). These results are in agreement with those cited in the main text. Using $e^{S}\simeq e^{-i\frac{\sigma_{y}}{2}\arctan\sqrt{\xi^{4}-1}}$ we can also compute $\langle x\sigma_{z}\rangle$:

[TABLE]

Finally, we consider $\langle p\sigma_{y}\rangle$, for which taking $e^{S}\simeq e^{-i\frac{\sigma_{y}}{2}\arctan\sqrt{\xi^{4}-1}}$ gives a vanishing result. In this case, we have to use the full expression given in Eq. (37). A straightforward calculation gives:

[TABLE]

which is the result cited in the main text. We have also confirmed that these analytical expressions reproduce accurately the expectation values obtained by direct numerics at large $\eta$.

Appendix C Analytical consideration on the finite-$\eta$ scaling behavior

We consider here the scaling form of the expectation values $\langle x^{2n}\rangle$ and $\langle p^{2n}\rangle$, at $\lambda\neq 0$ as well as $\lambda=0$ (the scaling behavior of the JC model was not discussed in Ref. Hwang and Plenio (2016)). For $\lambda>0$ we can immediately generalize the argument of the previous section to write:

[TABLE]

where $|\phi_{0}\rangle$ is the ground state of Eq. (31) at $\xi\simeq 1$. As discussed in the main text, the wavefinction has the scaling form $\phi_{0}(xM_{\lambda}^{1/6}/\sqrt{\eta},tM_{\lambda}^{1/3})$ and after a simple change of variable the expectation values are written as:

[TABLE]

where the expressions in the square parentheses are the explicit form for $X_{n}(tM_{\lambda}^{1/3})$, $P_{n}(tM_{\lambda}^{1/3})$ of the main text. For $\lambda<0$, using the exact mapping, we readily get the following scaling behavior (with $t^{\prime}=\xi^{\prime}-1$):

[TABLE]

Next, we consider the singular behavior at $\lambda=0$, where the well-known solution of the JC model allows for a relatively straighforward treatment. There are two eigenstates of the form $\alpha_{q}|q\rangle|-\rangle+\beta_{q}|q-1\rangle|+\rangle$ ($q=1,2,\ldots$) and the one with lower energy gives:

[TABLE]

At $\xi<1$ the ground state is $|0\rangle|-\rangle$ while at $\xi>1$ the ground state is given by the $q_{0}$ which minimizes Eq. (46).

As a side remark, we note that that Eq (46) has the same form of $\tilde{E}_{-}$ if we approximate $1-1/\eta\simeq 1$ [see Eq. (6) of the main text]. The only difference is that if $\tilde{x}^{2}$ is substituted here by $2q/\eta$. Therefore, by taking the $\eta\to\infty$ limit, $q_{0}/\eta$ is given by one half of Eq. (7) of the main text. It is also easy to check that $\beta_{q}\to 0$ when $\eta\to\infty$, thus the ground state of the JC model approaches $|q_{0}\rangle|-\rangle$. One immediate consequence, illustrated in Fig. 3(a) of the main text, is that $\langle\tilde{x}^{2}\rangle\simeq\langle q_{0}|\tilde{x}^{2}|q_{0}\rangle\simeq q_{0}/\eta$ is exactly half of the value which one would obtain for $\lambda>0$.

Returning to the derivation of the scaling function, we should determine the values $\xi_{q}$ at which the ground state changes from $|q\rangle$ to $|q+1\rangle$. These are easily found from $E_{JC}(q)=E_{JC}(q+1)$ which, to leading order in $\eta$, gives:

[TABLE]

Thus the $\xi_{q}$ are equally spaced on the $\xi$ axis, as illustrated in Fig. 3(b) of the main text. In each interval $[\xi_{q-1},\xi_{q}]$, the expectation values $\langle q|x^{2n}|q\rangle=\langle q|p^{2n}|q\rangle$ lead to Eq. (15) of the main text.

Appendix D Numerical scaling analysis

In the following we describe a numerical analysis of the scaling behavior. Having noticed the second-order nature of the phase transition at $\tilde{g}_{c}$, it is natural to speculate the following scaling law for a physical quantity $Q$ in the critical region:

[TABLE]

where $M=\eta^{2}$ is large enough but finite, $t=(\tilde{g}-\tilde{g_{c}})/\tilde{g_{c}}$ is the reduced coupling, $\beta_{Q}$ is the critical exponent for $Q$, $\nu$ is the critical exponent for $M$, and $\tilde{Q}_{\lambda}$ is the scaling function of $Q$ at a given $\lambda$. This scaling form is in the same spirit of the finite size scaling in classical thermodynamics phase transitions. Notice also that, differently from the main text, we take here a blind numerical approach, i.e., we neglect the $\lambda$ renormalization of $M$ revealed by our analytical study. It is also convenient to take the logarithm on both side of Eq. (48), to obtain the following $\log$-$\log$ relation between $\ln Q$ and $\ln M$:

[TABLE]

At critical point $t=0$, there is linear relation between $\ln Q$ and $\ln M$:

[TABLE]

where the slope of the linear dependence is the critical exponent ratio $\beta_{Q}/\nu$. Based on Eq. (50), we can determine the critical point as well as the critical exponent ratio $\beta_{Q}/\nu$.

We illustrate this procedure by considering $Q=\langle\tilde{x}^{2}\rangle$ at $\lambda=0.8$ as an example, for which we define:

[TABLE]

The left panel of Fig. 4 shows the $\log$-$\log$ relation with coupling strengths around the critical value $\tilde{g}_{c}=1.1111$. At the critical coupling $\tilde{g_{c}}$, we can see that the curve of $\ln\langle\tilde{x}^{2}\rangle$ as a function of $\ln M$ becomes a straight line, confirming the linear relation of Eq. (50). For small deviations from the critical coupling (i.e., $\tilde{g}=1.1108$ and $\tilde{g}=1.1112$), the linear character of $\ln\langle\tilde{x}^{2}\rangle$ breaks down. As seen in the first panel of Fig. 4, the curves at both sides of the critical point have a finite and opposite curvature. Thus, we can estimate the critical coupling as $\tilde{g}_{c}=1.1111\pm 0.0003$. At the same time, the slope of the straight line at $\tilde{g}=\tilde{g}_{c}$ gives the critical exponent ratio $\beta/\nu=0.33$. Furthermore, to confirm the scaling relation of the type Eq. (48) and determine the critical exponent $\nu$, we can set $v=tM^{1/\nu}$ as the variable of the horizontal axis and $M^{2\beta/\nu}\langle\tilde{x}^{2}\rangle$ as the variable of the vertical axis. The middle panel of Fig. 4 shows that, if the value of $\nu$ is properly chosen, curves with different scales of $M$ collapse into a single curve which corresponds to the scaling function of $\tilde{Q}_{\lambda}$. Doing this, we obtain $1/\nu=0.33\pm 0.02$.

We have performed this numerical analysis for several values of $\lambda$ and collected the results in Table 1, which also shows a comparison to the critical couplings and exponents obtained by our analytical method. We see that the two methods are in good agreement. The advantage of the numeirical approach is that it can be always applied, even when an analytical treatment might be difficult to achieve (i.e., for other more complicated models). However, the analytical approach allows here to reach several important conclusions regarding the universlity of the phase transition. First, the analytical results shows that the critical exponents are exactly identical for different values of $\lambda$ (while the numerical method can only establish the equality within errors). Second, the renormalization of $M$ into $M_{\lambda}$ is not obvious from the direct numerical scaling. To make this point explicit we show in the right panel of Fig. 4 that, within this traditional framework, the scaling functions at different $\lambda$ are actually different. Only by substituting $M$ with the renormalized value $M_{\lambda}$ the scaling functions at different $\lambda$ become identical, as shown in Fig. 2(a) of the main text. To establish that two critical points belong to the same universality class, it is widely acknowledged that one should prove the identity of both critical exponents and scaling functions. Thus, the analytical results are instrumental here to conclude that models with different values of $\lambda$ belong to the same universality class.