# A note on nonlocal approximations of sign-unrestricted solutions of
conservation laws

Alexander Keimer  Department of Mathematics
Friedrich-Alexander-University Erlangen-Nuremberg (FAU)
Cauerstraße 11
91058 Erlangen<EMAIL_ADDRESS>Lukas Pflug  Department of Mathematics
& Competence Center Scientific Computing
Friedrich-Alexander-University Erlangen-Nuremberg (FAU)
Cauerstraße 11 & Martensstraße 5a
91058 Erlangen<EMAIL_ADDRESS>
###### Abstract.

We study the singular limit problem for nonlocal conservation laws in which
the sign of the initial datum is unrestricted and the velocity of the
conservation law depends on a nonlocal approximation of the absolute value of
the density. We demonstrate that the nonlocal solutions converge to the local
entropy solution when the nonlocal kernel tends to a Dirac distribution, and
thus obtain an approximation result for local unsigned conservation laws,
generalizing the current results on the so-called sign-restricted singular
limit problem. The considered model class covers special cases like a
generalized Burgers’ equation and scalar versions of the Keyfitz–Kranzer
system.

###### Key words and phrases:

nonlocal conservation laws, singular limit problem, sign-unrestricted initial
datum, entropy solutions

###### 1991 Mathematics Subject Classification:

Primary 35L65; Secondary 35L04

##  1\. Introduction

Nonlocal conservation laws have received increased attention in the last
decade as they can model a variety of important physical and human-driven
dynamics [1, 4, 2, 3, 37, 19, 5, 36, 39, 26, 33, 6, 12, 11], and also possess
interesting mathematical properties that deviate from the theory of local
conservation laws [29, 21, 20]. Usually, the word “nonlocal” refers to the
fact that the velocity of the conservation law depends not only on the
solution at a given point in space and time but also on its integral in a
neighborhood of this point.

In this work, we focus in particular on nonlocal conservation laws in which
the nonlocality enters the flux nonlinearly, and the “locality” linearly; that
is, when the flux can be decomposed into

$f(q,\gamma\ast q)=V(\gamma\ast q)\cdot q$

where $\gamma$ is an integral kernel and $\ast$ denotes the spatial
convolution.

A question that has attracted great interest (first raised in [3]) is whether
one can recover the corresponding local entropy solution when the nonlocal
kernel converges to a Dirac distribution. Several publications have addressed
this question [43, 30, 22, 17, 8, 9, 15, 25, 14, 13, 32, 18, 13, 31]. For an
overview, we refer the reader to [18].

So far, however, all results require the initial datum—and thus the
solution—to be sign restricted. In this work, we will consider sign-
unrestricted initial data for cases in which the previously mentioned velocity
function $V$ depends only on the absolute value of the density, i.e.,

$V(x)=V(|x|)\quad\forall x\in I$

where $I$ represents all the solution’s values. Then, we replace the nonlocal
operator in the nonlocal approximation by its normed equivalent, i.e.,

$\gamma_{\eta}\ast q\text{ by
}W[|q_{\eta}|,\gamma_{\eta}]\coloneqq\gamma_{\eta}\ast|q|.$

On ${\Omega_{T}}$, the local and nonlocal conservation laws read as

$\partial_{t}q+\partial_{x}\big{(}V(|q|)q\big{)}=0\qquad\text{and}\qquad\partial_{t}q_{\eta}+\partial_{x}\big{(}V(W[|q_{\eta}|,\gamma_{\eta}])q_{\eta}\big{)}=0.$

Under some additional constraints, that the involved kernels are one-sided and
that the velocity function is monotone (these are classical assumptions in the
considered model class and are not restrictive), we show that the local
entropy solution of the corresponding local equation is recovered in the limit
$\gamma\rightarrow\delta$. Before we turn to the analysis, we highlight some
important special cases of this model class.

#### Keyfitz–Kranzer systems:

The class includes the scalar case of the Keyfitz–Kranzer system [38], i.e.,

$\displaystyle\partial_{t}q+\partial_{x}(q\phi(\|q\|))$ $\displaystyle=0$
(1.1)

with $q:[0,T]\times\mathbb{R}\rightarrow\mathbb{R}^{n}$ for $n\in\mathbb{N}$
and $\phi\in C^{1}(\mathbb{R}_{\geq 0})$, which is of interest in a number of
applications ranging from physics to engineering [42]. This contribution is
thus a first step toward approximating this type of system nonlocally while—in
contrast to viscosity approximations—preserving the hyperbolic nature of the
problem.

#### Generalized Burgers’ equation:

The approximation result is also applicable to the generalized Burgers’
equation [41], i.e.,

$\partial_{t}q+\partial_{x}\big{(}q^{2m+1}\big{)}=0$

with $m\in\mathbb{N}_{\geq 1}$ since in this case, the velocity is
$V\equiv(\cdot)^{2m}\equiv|\cdot|^{2m}$.

##  2\. Basic assumptions and problem setup

We study the nonlocal version of the following scalar conservation law as a
Cauchy problem on $\mathbb{R}$:

$\displaystyle\partial_{t}q+\partial_{x}\big{(}V(|q|)q\big{)}$
$\displaystyle=0,$ $\displaystyle(t,x)\in{\Omega_{T}}$ (2.1) $\displaystyle
q(0,\cdot)$ $\displaystyle\equiv q_{0},$ $\displaystyle x\in\mathbb{R}$

for $q:{\Omega_{T}}\coloneqq(0,T)\times\mathbb{R}\rightarrow\mathbb{R}$ and
$V:\mathbb{R}\rightarrow\mathbb{R}$. For
$q_{\eta}:(0,T)\times\mathbb{R}\rightarrow\mathbb{R}$, its nonlocal
approximation reads as

$\displaystyle\partial_{t}q_{\eta}+\partial_{x}\big{(}V(W[|q_{\eta}|,\gamma_{\eta}])q_{\eta}\big{)}$
$\displaystyle=0,$ $\displaystyle(t,x)\in{\Omega_{T}}$ (2.2) $\displaystyle
W[|q_{\eta}|,\gamma_{\eta}](t,x)$
$\displaystyle=\int^{x}_{-\infty}\\!\\!\\!\\!\\!\gamma_{\eta}(y-x)|q_{\eta}(t,y)|\,\mathrm{d}y,$
$\displaystyle(t,x)\in{\Omega_{T}}$ $\displaystyle q_{\eta}(0,\cdot)$
$\displaystyle\equiv q_{0},$ $\displaystyle x\in\mathbb{R}$

for a kernel $\gamma_{\eta}:\mathbb{R}_{\leq 0}\rightarrow\mathbb{R}_{\geq 0}$
yet to be specified. We will prove that

$q_{\eta}\overset{*}{\rightharpoonup}q^{*}\text{ in
}{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R})),\
W[|q_{\eta}|,\gamma_{\eta}]\rightarrow|q^{*}|\text{ in
}\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}$

where $q^{*}$ is the (entropy) weak solution of the local conservation law in
Eq. 2.1.

The fundamental difference from the nonlocal conservation laws studied so far
is that we take the absolute value in the nonlocal term, so that we
approximate the solution in the ${\mathsf{L}\\!}^{1}$-norm and can thus also
only expect convergence in the nonlinearity to the absolute value. As stated
before, the initial datum in the considered setup does not need to be sign
restricted. All assumptions that will be used in this contribution are
collected in the following.

######  Assumption 2.1.

We make these assumptions for the involved datum:

  * •

$V\in\mathsf{W}^{1,\infty}_{\;\textnormal{loc}}(\mathbb{R}_{>0}):\
V^{\prime}\geqq 0$,

  * •

$q_{0}\in{\mathsf{L}\\!}^{\infty}(\mathbb{R})\cap\mathsf{TV}(\mathbb{R})$,

  * •

$\gamma_{\eta}\in{\mathsf{BV}}(\mathbb{R}_{\leq 0};\mathbb{R}_{\geq 0})$
satisfying either the assumptions on the kernel in [17, Eq. (1.5) and (1.8)]
or [32, Assumption 2].

######  Remark 2.2 (The meaning of the assumptions on the involved kernel).

Both classes of kernels must have $L^{1}$ mass equal to $1$ and be non-
negative. The assumptions also imply the following:

  * •

The kernels satisfying [17, Eq. (1.5) and (1.8)] need to increase in a convex
manner on $\mathbb{R}_{\leq 0}$.

  * •

The conditions on the kernels in [32, Assumption 2] are more involved as these
kernels are scaled in the exponent. This enables the use of non-monotonic and
rather irregular kernels.

We refer the reader to the corresponding references for further details.

To provide a better understanding of why the claimed convergence will hold, we
point out that the considered problem class conserves the
${\mathsf{L}\\!}^{1}$-norm. This becomes evident when multiplying the
conservation law in Eq. 2.1 or Eq. 2.2 by $\textnormal{sgn}(q)$ or
$\textnormal{sgn}(q_{\eta}),$ respectively.

##  3\. The involved equations

###  3.1. The nonlocal equation

As the type of nonlocal dynamics studied here has not been considered in the
literature before, we will state an existence and uniqueness theorem, which
also guarantees a (weak) form of a maximum principle on $|q_{\eta}|$. However,
we first define what we mean by a weak solution.

######  Definition 3.1 (Weak solution to the nonlocal conservation law).

We call
$q_{\eta}\in\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}\cap{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R}))$
a weak solution to the nonlocal conservation law in Eq. 2.2 iff
$\forall\phi\in\mathsf{C}^{1}_{\text{c}}((-42,T)\times\mathbb{R})$ it holds
that

$\displaystyle\iint_{{\Omega_{T}}}q_{\eta}(t,x)\big{(}\phi_{t}(t,x)+\phi_{x}(t,x)V(W[|q_{\eta}|,\gamma_{\eta}](t,x))\big{)}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\qquad+\int_{\mathbb{R}}\phi(0,x)q_{0}(x)\,\mathrm{d}x=0.$

######  Theorem 3.2 (Existence and uniqueness of the nonlocal conservation law
and a maximum principle).

Let Ass. 2.1 hold. Then, there exists a unique weak solution to Eq. 2.2,

$q_{\eta}\in\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}\cap{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}\cap\mathsf{TV}(\mathbb{R}))$

on any time horizon $T\in\mathbb{R}_{>0}$. Additionally, for the assumptions
on the kernel, we have the following:

  * •

In [17, Eq. (1.5) and (1.8)], the solution satisfies a maximum principle

$\|q_{\eta}(t,\cdot)\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})}\leq\|q_{0}\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})},\forall
t\in[0,T].$

  * •

In [32, Assumption 2], a weak maximum principle holds: for any
$(T,\kappa)\in\mathbb{R}_{>0}^{2}$ there exists
$\eta_{T,\kappa}\in\mathbb{R}_{>0}$ such that

$\forall(t,\eta)\in[0,T]\times(0,\eta_{T,\kappa}):\
\|q_{\eta}(t,\cdot)\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})}\leq(1+\kappa)\|q_{0}\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})},\
\forall t\in[0,T].$

Finally, the following approximation result holds. Assuming
$q_{0,\varepsilon}\in
C^{\infty}(\mathbb{R})\cap\mathsf{TV}(\mathbb{R})\cap{\mathsf{L}\\!}^{\infty}(\mathbb{R})$
and denoting with $q_{\eta,\varepsilon}$ the corresponding solution to the
nonlocal equation, there exists $C=C(T,\eta)\in\mathbb{R}_{>0}$ such that

$\|q_{\eta}-q_{\eta,\varepsilon}\|_{\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}(\mathbb{R}))}\leq
C\|q_{0}-q_{0,\varepsilon}\|_{{\mathsf{L}\\!}^{1}(\mathbb{R})}$

and $q_{\eta,\varepsilon}\in\mathsf{W}^{1,\infty}({\Omega_{T}})$, and as such
a strong solution of the nonlocal conservation law in Eq. 2.2.

###### Proof.

We only sketch the proof. The existence and uniqueness of weak solutions
follow with classical fixed-point arguments laid out, for instance, in [29,
16]. The integral operator in the velocity function $V$ renders the
corresponding characteristics and the related system of initial value problems
in the ODE Lipschitz-continuous so that on a small time horizon, it can be
shown that a unique weak solution exists using classical fixed-point
arguments.

The long-time-horizon existence is then a consequence of the maximum
principle, which is addressed later. To this end, one can derive a related
equation in the absolute value of the original solution. Indeed, multiplying
Eq. 2.2 by $\textnormal{sgn}(q_{\eta}(t,x)),\ (t,x)\in{\Omega_{T}}$, one
obtains, formally,

$\displaystyle\partial_{t}|q_{\eta}(t,x)|+\partial_{x}\big{(}V(W[|q_{\eta}|,\gamma_{\eta}](t,x))|q_{\eta}(t,x)|\big{)}$
$\displaystyle=0,$ $\displaystyle(t,x)\in{\Omega_{T}}$ $\displaystyle|q(0,x)|$
$\displaystyle=|q_{0}(x)|,$ $\displaystyle x\in\mathbb{R}.$

Thus, the quantity $|q_{\eta}|$ satisfies the well-known sign-restricted
nonlocal conservation law studied in, e.g., [29, 27, 17]. In particular, for
the kernel in [17, Eq. (1.5) and (1.8)], the claimed maximum principle for
$|q_{\eta}|$ (see, e.g., [17, Proposition 2.1], [27, Theorem 2.1], [16,
Corollary 2.1], [29, Corollary 4.3, Remark 4.4]) holds.

For the weak form of maximum principle that applies with a kernel as in [32,
Assumption 2], we refer the reader to [32, Theorem 6]. The maximum principle
on $|q_{\eta}|$ carries over to the maximum/minimum principle for $q_{\eta}$.
For the claimed stability and regularity of $q_{\eta,\varepsilon}$ we refer
the reader to [32, Theorem 3] and to [30, Theorem 3.1], concluding the sketch
of the proof. ∎

###  3.2. The local equation

Concerning the theory of local conservation laws, we give in the following the
definition of weak and entropy solutions, and the basic result on existence
and uniqueness.

######  Definition 3.3 (Weak solution for the local equation).

We call a function
$q^{*}\in{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R}))$
a weak solution to Eq. 2.1 for the initial datum
$q_{0}\in{\mathsf{L}\\!}^{\infty}\cap\mathsf{TV}(\mathbb{R})$ iff
$\forall\phi\in\mathsf{C}^{1}_{\text{c}}((-42,T)\times\mathbb{R})$ the
following integral equation holds:

$\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\\!q^{*}(t,x)\big{(}\partial_{t}\phi(t,x)+\partial_{x}\phi(t,x)V(|q^{*}(t,x)|)\big{)}\,\mathrm{d}x\,\mathrm{d}t+\\!\\!\int_{\mathbb{R}}\\!\\!\phi(0,x)q_{0}(x)\,\mathrm{d}x=0.$

As is well known, weak solutions are not unique, which is why one postulates
an entropy condition (which also contains the weak formulation).

######  Definition 3.4 (Entropy solution).

We call the function
$q\in\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R}))\cap{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R}))$
a weak entropy solution to Eq. 2.1 iff it satisfies the following entropy
inequality for all convex $\alpha\in\mathsf{C}^{2}(\mathbb{R})$ and
$\beta\in\mathsf{C}^{1}(\mathbb{R})$ satisfying
$\beta^{\prime}\equiv\alpha^{\prime}\cdot f^{\prime},$ where

$f(x)=xV(x),\ \forall x\in\mathbb{R},$

and for all
$\phi\in\mathsf{C}_{\text{c}}^{1}\big{(}(-42,T)\times\mathbb{R};\mathbb{R}_{\geq
0}\big{)}$:

$\displaystyle\mathcal{E}[\phi,\alpha,q]$
$\displaystyle\coloneqq\iint_{{\Omega_{T}}}\\!\\!\\!\\!\alpha(q(t,x))\phi_{t}(t,x)+\beta(q(t,x))\phi_{x}(t,x)\,\mathrm{d}x\,\mathrm{d}t$
(3.1)
$\displaystyle\quad+\int_{\mathbb{R}}\\!\\!\alpha(q_{0}(x))\phi(0,x)\,\mathrm{d}x\geq
0.$

######  Theorem 3.5 (Existence and uniqueness of local solutions).

Let the (local) conservation law in Eq. 2.1 be given with the assumptions in
Ass. 2.1. Then, there exists on every time horizon $T\in\mathbb{R}_{>0}$ a
unique weak entropy solution in the sense of Def. 3.4.

###### Proof.

This can be found in [7, Theorem 6.3], [24, Theorem 19.1]. For Kružkov
entropies and entropy flux pairs, we refer to [34, Theorem 2, Theorem 5,
Section 5 Item 4], or [28]. ∎

##  4\. The singular limit problem

As we are interested in the convergence of $q_{\eta}$ for $\eta\rightarrow 0$,
a compactness result on the nonlocal term is required, which is stated in the
following lemma.

######  Lemma 4.1 (Compactness of $W_{\eta}[|q_{\eta}|,\gamma_{\eta}]$ in
$\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R}))$).

Given the assumptions of Thm. 3.2, and assuming that the kernel satisfies [17,
Eq. (1.5) and (1.8)] or [32, Assumption 2], there exists for every
$T\in\mathbb{R}_{>0}$ an $\varepsilon\in\mathbb{R}_{>0}$ such that

$\big{\\{}W[|q_{\eta}|,\gamma_{\eta}]\in\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}:\
\eta\in(0,\varepsilon)\big{\\}}\overset{\text{c}}{\hookrightarrow}\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}.$

###### Proof.

This is a consequence of the $\mathsf{TV}$ bounds uniform in
$\eta\in\mathbb{R}_{>0}$ for both classes of kernels (see [17, Theorem 1] and
[32, Theorem 7]) together with the uniform bounds of the solutions in $\eta$,
thanks to the maximum principle Thm. 3.2. Together with an equicontinuity in
time as required by [40, Lemma 1], and which is also a consequence of the
uniform $\mathsf{TV}$ estimate, the named compactness follows (as also
explicitly stated in [32, Theorem 8]). ∎

Equipped with the results from the previous section, we show here that the
solution to the nonlocal equation converges to the local entropy solution when
the nonlocal operator converges to a point evaluation. We start with a short
lemma showing the convergence to a weak local solution not only of the
nonlocal operator but also weak-star convergence of the solution itself.

######  Lemma 4.2 (Convergence to a weak local solution).

Let Ass. 2.1 be given. Then for any $\eta_{k}\in\mathbb{R}_{>0},\
k\in\mathbb{N}_{\geq 1},\lim_{k\rightarrow\infty}\eta_{k}=0$, there exists a
subsequence (still denoted by $\eta_{k}$) such that the nonlocal operator
converges in the following sense:

$\lim_{k\rightarrow\infty}\|W[|q_{\eta}|,\gamma_{\eta}]-|q^{*}|\|_{\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R}))}=0,$

as well as

$q_{\eta}\overset{*}{\rightharpoonup}q^{*}\ \text{ in }\
{\mathsf{L}\\!}^{\infty}({\Omega_{T}}),$

where
$q^{*}\in\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}$
is a weak solution to the local equation.

###### Proof.

Thanks to the uniform boundedness of $q_{\eta}$, as stated in Thm. 3.2, we can
invoke the Banach–Alaoglu–Bourbaki theorem [10, Theorem 3.16] and we obtain
that there exists a subsequence (still denoted by $\eta$) and a
$q^{*}\in{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R}))$
such that

$q_{\eta}\overset{*}{\rightharpoonup}q^{*}\ \text{ in }\
{\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R})).$

Using this subsequence in the nonlocal operator, we can apply the compactness
result in Lem. 4.1 to deduce that there exists another subsequence, denoted by
$\big{(}\eta_{k}\big{)}_{k\in\mathbb{N}}\subset\mathbb{R}_{>0},\
\lim_{k\rightarrow\infty}\eta_{k}=0$ such that

$\lim_{k\rightarrow\infty}\|W[|q_{\eta_{k}}|,\gamma_{\eta_{k}}]-|q^{*}|\|_{\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R}))}=0.$

Next, we need to show that this limit point $q^{*}$ is a weak local solution.
For $\eta\in\mathbb{R}_{>0}$, we have that $q_{\eta_{k}}$ and the
corresponding nonlocal operator satisfy

$\iint_{{\Omega_{T}}}q_{\eta_{k}}\big{(}\phi_{t}+\phi_{x}V(W[|q_{\eta_{k}}|,\gamma_{\eta}])\big{)}\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}}\phi(0,\cdot)q_{0}\,\mathrm{d}x=0.$

However, thanks to the strong $\mathsf{C}({\mathsf{L}\\!}^{1})$ convergence of
the nonlocal operator,
$V\in\mathsf{W}^{1,\infty}_{\;\textnormal{loc}}(\mathbb{R})$ and $q_{\eta}$
weak-star converges in ${\mathsf{L}\\!}^{\infty}({\mathsf{L}\\!}^{\infty})$,
and thanks to the fact that the product of a strongly converging sequence and
a weak-star converging sequence converges weak-star, we obtain in the limit
$k\rightarrow\infty$ that

$\iint_{{\Omega_{T}}}q^{*}\big{(}\phi_{t}+\phi_{x}V(|q^{*}|)\big{)}\,\mathrm{d}x\,\mathrm{d}t+\int_{\mathbb{R}}\phi(0,\cdot)q_{0}\,\mathrm{d}x=0,$

which is the weak formulation for the local conservation law, as stated in
Def. 3.3. Thus, $q^{*}$ is a weak solution to the local conservation law,
concluding the proof. ∎

######  Remark 4.3 (Convergence and weak solutions).

As is well known, the convergence to a weak solution of the local equation
follows with classical compactness results. However, these weak solutions are
in the local case not unique, requiring some well-known and established
entropy conditions to single out the “reasonable” unique solution. This is
what we will discuss in the following.

###  4.1. Entropy admissibility

The following calculations will lead to the crucial result, that if $q_{\eta}$
converges to a limit then it satisfies the entropy admissibility condition.
However, we will present these calculations only for the one-sided exponential
kernel, strictly convex “flux” and provide some remarks about how this can be
generalized to the arbitrary kernels. In Fig. 1, the mentioned convergence is
numerically exemplified for two initial data and velocity functions.

Figure 1. Solution of the nonlocal balance law stated in Eq. 2.2 for
$\eta\in\\{10^{-1},10^{-2},10^{-3}\\}$ (red to blue) for initial datum
$q_{0}\equiv-\tfrac{1}{2}\chi_{(-0.5,0)}+\chi_{(0,0.5)}$ (top) and
$q_{0}\equiv\tfrac{1}{2}\chi_{(-0.5,0)}-\chi_{(0,0.5)}$ (bottom), and velocity
$V\equiv\text{Id}$ (left) and $V\equiv\text{Id}^{2}$ (right). This illustrates
the proven convergence for $\eta\rightarrow 0$ to entropy solutions to the
local counterpart of the nonlocal conservation law, as stated in Eq. 2.1.

######  Theorem 4.4 (Nonlocal solution entropy admissible in the
limit/convergence).

Choose as nonlocal kernel $\gamma_{\eta}\equiv\exp(\cdot/\eta)/\eta$ and
consider the corresponding nonlocal solution
$q_{\eta}\in\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}$
of Eq. 2.2. Assume in addition that for $\varepsilon\in\mathbb{R}_{>0}$, the
flux $x\mapsto xV(x),\
x\in[-\|q_{0}\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})}-\varepsilon,\|q_{0}\|_{{\mathsf{L}\\!}^{\infty}(\mathbb{R})}+\varepsilon]$,
is strictly concave. Then, this solution $q_{\eta}$ is entropy admissible in
the singular limit, i.e., it satisfies $\forall\phi,\alpha$, as in Def. 3.4,

$\lim_{\eta\rightarrow 0}\mathcal{E}[\phi,\alpha,q_{\eta}]\geq 0.$

In particular, the nonlocal $q_{\eta}$ converges weak-star in
${\mathsf{L}\\!}^{\infty}((0,T);{\mathsf{L}\\!}^{\infty}(\mathbb{R}))$ to the
local entropy solution $q^{*}$, and the nonlocal operator
$W[|q_{\eta}|,\gamma_{\eta}]$ converges strongly in
$\mathsf{C}\big{(}[0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R})\big{)}$
to $|q^{*}|$.

###### Proof.

We first remark that in the case of the considered exponential kernel, we
obtain the following identity $\forall(t,x)\in{\Omega_{T}}$:

$\partial_{x}W_{\eta}(t,x)=\tfrac{1}{\eta}\big{(}|q_{\eta}(t,x)|-W_{\eta}(t,x)\big{)}.$
(4.1)

Here we use the abbreviation $W_{\eta}\coloneqq
W_{\eta}[|q_{\eta}|,\tfrac{1}{\eta}\exp(\tfrac{\cdot}{\eta})]$, which will be
useful later.

Assuming sufficient regularity on the nonlocal equation, which is possible
thanks to the stability in ${\mathsf{L}\\!}^{1}$ (see again Thm. 3.2), we can
perform the following manipulations when plugging $q_{\eta}$ into the entropy
admissibility condition:

$\displaystyle\mathcal{E}[\phi,\alpha,q_{\eta}]$
$\displaystyle=-\iint_{{\Omega_{T}}}\Big{(}\alpha^{\prime}(q_{\eta})\partial_{t}q_{\eta}+\beta^{\prime}(q_{\eta})\partial_{x}q_{\eta}\Big{)}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=-\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\alpha^{\prime}(q_{\eta})\Big{(}\partial_{t}q_{\eta}+f^{\prime}(q_{\eta})\partial_{x}q_{\eta}\Big{)}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\alpha^{\prime}(q_{\eta})\Big{(}\partial_{x}\big{(}V(W_{\eta})q_{\eta}\big{)}-V(|q_{\eta}|)\partial_{x}q_{\eta}-V^{\prime}(|q_{\eta}|)\partial_{x}|q_{\eta}|q_{\eta}\Big{)}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=\iint_{{\Omega_{T}}}\alpha^{\prime\prime}(q_{\eta})\partial_{x}q_{\eta}\big{(}V(|q_{\eta}|)-V(W_{\eta})\big{)}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t$
(4.2)
$\displaystyle\quad+\iint_{{\Omega_{T}}}\alpha^{\prime}(q_{\eta})\big{(}V(|q_{\eta}|)-V(W_{\eta})\big{)}q_{\eta}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t.$

To show a lower bound for $\mathcal{E},$ we first focus on the first integral
in the latter expression. Take $\alpha^{\prime\prime}=1$ (this is not a
restriction for the entropy admissibility as we can take advantage of the
results in [23, 35], which guarantee the uniqueness of solutions when only
using one entropy flux pair) to obtain

$\displaystyle\eqref{eq:42}$
$\displaystyle=\iint_{{\Omega_{T}}}\partial_{x}q_{\eta}\big{(}V(|q_{\eta}|)-V(W_{\eta})\big{)}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=\iint_{{\Omega_{T}}}V(|q_{\eta}|)|q_{\eta}|\partial_{x}|q_{\eta}|\phi\,\mathrm{d}x\,\mathrm{d}t-\iint_{{\Omega_{T}}}V(W_{\eta})q_{\eta}\partial_{x}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t.$
With $G(y)\coloneqq\int_{0}^{y}sV(s)\,\mathrm{d}s$ and by integration by parts
in the first integral, we continue: $\displaystyle\eqref{eq:42}$
$\displaystyle=-\iint_{{\Omega_{T}}}G(|q_{\eta}|)\phi_{x}\,\mathrm{d}x\,\mathrm{d}t-\iint_{{\Omega_{T}}}V(W_{\eta})q_{\eta}\partial_{x}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=-\iint_{{\Omega_{T}}}\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t-\iint_{{\Omega_{T}}}G(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad-\iint_{{\Omega_{T}}}V(W_{\eta})q_{\eta}\partial_{x}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t.$
Another integration by parts in the second term yields
$\displaystyle\eqref{eq:42}$
$\displaystyle=-\iint_{{\Omega_{T}}}\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\iint_{{\Omega_{T}}}V(W_{\eta})W_{\eta}\partial_{x}W_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t-\iint_{{\Omega_{T}}}V(W_{\eta})q_{\eta}\partial_{x}q_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t,$
and integration by parts in the second and third expressions leads to
$\displaystyle\eqref{eq:42}$
$\displaystyle=-\iint_{{\Omega_{T}}}\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad-\tfrac{1}{2}\iint_{{\Omega_{T}}}V(W_{\eta})W_{\eta}^{2}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t-\tfrac{1}{2}\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})W_{\eta}^{2}\partial_{x}W_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\tfrac{1}{2}\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})\partial_{x}W_{\eta}q_{\eta}^{2}\phi\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\iint_{{\Omega_{T}}}V(W_{\eta})q_{\eta}^{2}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=-\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\\!\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}q_{\eta}^{2}-W_{\eta}^{2}\big{)}V(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad-\tfrac{1}{2}\iint_{{\Omega_{T}}}W_{\eta}^{2}V^{\prime}(W_{\eta})\partial_{x}W_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})\partial_{x}W_{\eta}q_{\eta}^{2}\phi\,\mathrm{d}x\,\mathrm{d}t.$
Using the identity for the nonlocal operator in Eq.
4.1$|q_{\eta}|=W_{\eta}+\eta\partial_{x}W_{\eta}$, which implies
$q_{\eta}^{2}=W_{\eta}^{2}+2\eta
W_{\eta}\partial_{x}W_{\eta}+\eta^{2}\big{(}\partial_{x}W_{\eta}\big{)}^{2}$,
in the last term yields $\displaystyle\eqref{eq:42}$
$\displaystyle=-\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\Big{(}G(|q_{\eta}|)-G(W_{\eta})\Big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\\!\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}q_{\eta}^{2}-W_{\eta}^{2}\big{)}V(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad-\tfrac{1}{2}\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})\partial_{x}W_{\eta}W_{\eta}^{2}\phi\,\mathrm{d}x+\tfrac{1}{2}\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})\partial_{x}W_{\eta}W_{\eta}^{2}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\eta\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!V^{\prime}(W_{\eta})\big{(}\partial_{x}W_{\eta}\big{)}^{2}W_{\eta}\phi\,\mathrm{d}x\,\mathrm{d}t+\tfrac{\eta^{2}}{2}\\!\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!V^{\prime}(W_{\eta})\big{(}\partial_{x}W_{\eta}\big{)}^{3}\phi\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle=-\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\\!\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}q_{\eta}^{2}-W_{\eta}^{2}\big{)}V(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\eta\iint_{{\Omega_{T}}}V^{\prime}(W_{\eta})\big{(}\partial_{x}W_{\eta}\big{)}^{2}\big{(}W+\tfrac{\eta}{2}W_{x}\big{)}\phi\,\mathrm{d}x\,\mathrm{d}t.$
Again, thanks to the identity in Eq. 4.1, as well as
$|q_{\eta}|=W_{\eta}+\eta\partial_{x}W_{\eta}$ and $V^{\prime}\geqq 0$, we can
estimate the last term by zero from below and obtain
$\displaystyle\eqref{eq:42}$
$\displaystyle\geq-\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t+\tfrac{1}{2}\\!\\!\iint_{{\Omega_{T}}}\\!\\!\\!\\!\\!\big{(}q_{\eta}^{2}-W_{\eta}^{2}\big{)}V(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t.$

Thus, in total, we obtain

$\displaystyle\mathcal{E}[\phi,\alpha,q_{\eta}]$
$\displaystyle\geq-\iint_{{\Omega_{T}}}\big{(}G(|q_{\eta}|)-G(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\tfrac{1}{2}\iint_{{\Omega_{T}}}\big{(}q_{\eta}-W_{\eta}\big{)}\big{(}q_{\eta}+W_{\eta}\big{)}V(W_{\eta})\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\quad+\iint_{{\Omega_{T}}}q_{\eta}^{2}\big{(}V(|q_{\eta}|)-V(W_{\eta})\big{)}\phi_{x}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\overset{\eta\rightarrow 0}{=}0,$

which holds when recalling that, by definition, we have the following:

  * •

$V,G\in\mathsf{W}^{1,\infty}_{\;\textnormal{loc}}(\mathbb{R})$ according to
Ass. 2.1 and $G(\cdot)\coloneqq\int_{0}^{\cdot}xV(x)\,\mathrm{d}x$,

  * •

$\||q_{\eta}|-W_{\eta}\|_{\mathsf{C}([0,T];{\mathsf{L}\\!}^{1}_{\;\textnormal{loc}}(\mathbb{R}))}\overset{\eta\rightarrow
0}{\longrightarrow}0$ according to Lem. 4.2,

  * •

$q_{\eta}$ and $W_{\eta}$ are uniformly bounded in $\eta$ thanks to Thm. 3.2.

Up to now, this convergence only holds for a suitable subsequence. But as the
local entropy solution is unique (see Thm. 3.5) and we can show for all
subsequences that there exists a subsequence converging to the local entropy
solution, i.e., to the same limit, we arrive at the convergence to the entropy
solution for any sequence. This concludes the proof. ∎

######  Remark 4.5 (Generalization).

Similarly to the proof of the latter theorem, one can show, as in [32] and
[17], that the convergence to the entropy solution in Thm. 4.4 holds for all
kernels defined in Ass. 2.1 and that the strict concavity assumption can be
removed as well, so that in the end, we obtain the general convergence result
for the problem class in Eq. 2.2.

### Acknowledgment

L. Pflug acknowledges funding by the Deutsche Forschungsgemeinschaft (German
Research Foundation)—Project-ID 416229255—SFB 1411.

## References

  * [1] D. Amadori and W. Shen.  An integro-differential conservation law arising in a model of granular flow.  Journal of Hyperbolic Differential Equations, 9(1):105–131, 2012.
  * [2] P. Amorim.  On a nonlocal hyperbolic conservation law arising from a gradient constraint problem.  Bulletin of the Brazilian Mathematical Society, New Series, 43(4):599–614, 2012.
  * [3] P. Amorim, R.M. Colombo, and A. Teixeira.  On the numerical integration of scalar nonlocal conservation laws.  ESAIM: Mathematical Modelling and Numerical Analysis, 49(1):19–37, 2015.
  * [4] D. Armbruster, P. Degond, and C. Ringhofer.  A model for the dynamics of large queuing networks and supply chains.  SIAM Journal on Applied Mathematics, 66(3):896–920, 2006.
  * [5] A. Bayen, J. Friedrich, A. Keimer, L. Pflug, and T. Veeravalli.  Modeling multilane traffic with moving obstacles by nonlocal balance laws.  SIAM Journal on Applied Dynamical Systems, 21(2):1495–1538, June 2022.
  * [6] S. Blandin and P. Goatin.  Well-posedness of a conservation law with non-local flux arising in traffic flow modeling.  Numerische Mathematik, 132(2):217–241, 2016.
  * [7] A. Bressan.  Hyperbolic Systems of Conservation Laws.  Oxford University Press, Oxford, 2000.
  * [8] A. Bressan and W. Shen.  On traffic flow with nonlocal flux: a relaxation representation.  Archive for Rational Mechanics and Analysis, 237, 2020.
  * [9] A. Bressan and W. Shen.  Entropy admissibility of the limit solution for a nonlocal model of traffic flow.  Communications in Mathematical Sciences, 19(5):1447–1450, 2021.
  * [10] H. Brezis.  Functional analysis, Sobolev spaces and partial differential equations.  Universitext. Springer, New York, 2011.
  * [11] F. A. Chiarello, P. Goatin, and E. Rossi.  Stability estimates for non-local scalar conservation laws.  Nonlinear Analysis: Real World Applications, 45:668–687, 2019.
  * [12] F.A. Chiarello and P. Goatin.  Global entropy weak solutions for general non-local traffic flow models with anisotropic kernel.  ESAIM: Math. Modelling and Numerical Analysis, 52(1):163–180, 2018.
  * [13] F.A. Chiarello and A. Keimer.  On the singular limit problem in nonlocal balance laws: Applications to nonlocal lane-changing traffic flow models, 2023.
  * [14] G. M. Coclite, M. Colombo, G. Crippa, N. De Nitti, A. Keimer, E. Marconi, L. Pflug, and L.V. Spinolo.  Oleĭnik-type estimates for nonlocal conservation laws and applications to the nonlocal-to-local limit, 2023.
  * [15] G. M. Coclite, J.-M. Coron, N. De Nitti, A. Keimer, and L. Pflug.  A general result on the approximation of local conservation laws by nonlocal conservation laws: The singular limit problem for exponential kernels.  Annales de l'Institut Henri Poincaré C, Analyse non linéaire, November 2022.
  * [16] G. M. Coclite, N. De Nitti, A. Keimer, and L. Pflug.  On existence and uniqueness of weak solutions to nonlocal conservation laws with bv kernels.  Zeitschrift für angewandte Mathematik und Physik, 73(6):241, 2022.
  * [17] M. Colombo, G. Crippa, E. Marconi, and L. V. Spinolo.  Nonlocal traffic models with general kernels: Singular limit, entropy admissibility, and convergence rate.  Archive for Rational Mechanics and Analysis, 247(2), February 2023.
  * [18] M. Colombo, G. Crippa, E. Marconi, and L.V. Spinolo.  An overview on the local limit of non-local conservation laws, and a new proof of a compactness estimate, 2023.
  * [19] R.M. Colombo, M. Garavello, and M. Lécureux-Mercier.  A class of nonlocal models for pedestrian traffic.  Mathematical Models and Methods in Applied Sciences, 22(04):1150023, 2012.
  * [20] J.-M. Coron, M. Kawski, and Z. Wang.  Analysis of a conservation law modeling a highly re-entrant manufacturing system.  Discrete and Continuous Dynamical Systems - Series B, 14(4):1337–1359, 2010.
  * [21] G. Crippa and M. Lécureux-Mercier.  Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow.  Nonlinear Differential Equations and Applications NoDEA, 20(3):523–537, 2013.
  * [22] G. Crippa, E. Marconi, L.V. Spinolo, and M. Colombo.  Local limit of nonlocal traffic models: Convergence results and total variation blow-up.  Annales de l'Institut Henri Poincaré C, Analyse non linéaire, 38(5):1653–1666, October 2021.
  * [23] C. De Lellis, F. Otto, and M. Westdickenberg.  Minimal entropy conditions for Burgers equation.  Quarterly of Applied Mathematics, 62(4):687–700, 2004.
  * [24] R. Eymard, T. Gallouët, and R. Herbin.  Finite volume methods.  Handbook of Numerical Analysis, 7:713–1018, 2000.
  * [25] J. Friedrich, S. Göttlich, A. Keimer, and L. Pflug.  Conservation laws with nonlocal velocity – the singular limit problem, 2022.
  * [26] J. Friedrich, O. Kolb, and S. Göttlich.  A Godunov type scheme for a class of LWR traffic flow models with non-local flux.  Networks & Heterogeneous Media, 13:531, 2018.
  * [27] P. Goatin and S. Scialanga.  Well-posedness and finite volume approximations of the LWR traffic flow model with non-local velocity.  Networks and Hetereogeneous Media, 11(1):107–121, 2016.
  * [28] E. Godlewski and P.A. Raviart.  Hyperbolic Systems of Conservation Laws.  Mathématiques & applications. Ellipses, Paris, 1991.
  * [29] A. Keimer and L. Pflug.  Existence, uniqueness and regularity results on nonlocal balance laws.  Journal of Differential Equations, 263:4023–4069, 2017.
  * [30] A. Keimer and L. Pflug.  On approximation of local conservation laws by nonlocal conservation laws.  Journal of Mathematical Analysis and Applications, 475(2):1927–1955, 2019.
  * [31] A. Keimer and L. Pflug.  On the singular limit problem for a discontinuous nonlocal conservation law.  Nonlinear Analysis, 237:113381, December 2023.
  * [32] A. Keimer and L. Pflug.  On the singular limit problem for nonlocal conservation laws: A general approximation result for kernels with fixed support, 2023.
  * [33] P. E. Kloeden and T. Lorenz.  Nonlocal multi-scale traffic flow models: analysis beyond vector spaces.  Bulletin of Mathematical Sciences, 6(3):453–514, 2016.
  * [34] S. N. Kružkov.  First order quasilinear equations in several independent variables.  Mathematics of the USSR-Sbornik, 10(2):217, 1970.
  * [35] E. Y. Panov.  Uniqueness of the solution of the cauchy problem for a first order quasilinear equation with one admissible strictly convex entropy.  Mathematical Notes, 55(5):517–525, 1994.
  * [36] L. Pflug, T. Schikarski, A. Keimer, W. Peukert, and M. Stingl.  eMoM: Exact method of moments—nucleation and size dependent growth of nanoparticles.  Computers & Chemical Engineering, 136:106775, 2020.
  * [37] B. Piccoli, N.P. Duteil, and E. Trélat.  Sparse control of Hegselmann–Krause models: Black hole and declustering.  SIAM Journal on Control and Optimization, 57(4):2628–2659, 2019.
  * [38] N. H. Risebro and F. Weber.  A note on front tracking for the Keyfitz–Kranzer system.  Journal of Mathematical Analysis and Applications, 407(2):190–199, 2013.
  * [39] E. Rossi, J. Weißen, P. Goatin, and S. Göttlich.  Well-posedness of a non-local model for material flow on conveyor belts.  ESAIM: Mathematical Modelling and Numerical Analysis, 54(2):679–704, 2020.
  * [40] J. Simon.  Compact sets in the space $L^{p}(0,T;B)$.  Annali di Matematica Pura ed Applicata (4), 146:65–96, 1987.
  * [41] A. S. Tersenov.  On the generalized Burgers equation.  Nonlinear Differential Equations and Applications NoDEA, 17(4):437–452, February 2010.
  * [42] Aslak Tveito and Ragnar Winther.  Existence, uniqueness, and continuous dependence for a system of hyperbolic conservation laws modeling polymer flooding.  SIAM journal on mathematical analysis, 22(4):905–933, 1991.
  * [43] K. Zumbrun.  On a nonlocal dispersive equation modeling particle suspensions.  Quarterly of Applied Mathematics, 57(3):573–600, 1999.


# Quantum fluctuations and unusual critical exponents in a quantum Rabi
Triangle

Xiao Qin  Department of Physics, Chongqing Key Laboratory for strongly coupled
Physics, Chongqing University, Chongqing 401330, China  Yu-Yu Zhang
<EMAIL_ADDRESS>Department of Physics, Chongqing Key Laboratory for strongly
coupled Physics, Chongqing University, Chongqing 401330, China

###### Abstract

Quantum fluctuations of a quantum Rabi triangle are studied using an
analytical approach beyond the mean-field theory. By applying an artificial
magnetic field among three cavities, time-reversal symmetry breaking is
manifested through a directional transfer dynamics of photons. In contrast to
previous studies, we focus on the scaling exponents of the fluctuations of the
local photon number and the position variance near the critical point. By
accurate calculation using Bogoliubov transformation we show that two scaling
laws emerge respectively for the frustrated cavity and the remaining cavities,
which are associated with the geometric frustrations. Especially, for the
frustrated cavity, the scaling exponent in the chiral superradiant phase is
different from that in the frustrated antiferromagnetic superradiant phase
without an artificial magnetic field. The unusual scaling exponents predict
distinct universality classes from the single-cavity Rabi universality. We
suggest that the accurate critical exponents in few-body system is useful for
identifying exotic quantum phase transition in light-matter coupling system.

##  I Introduction

The interaction between light and matter has brought forth a new class of
quantum many-body system in understanding of strongly correlated systems and
quantum phase transition [1, 2, 3, 4]. The effect of quantum fluctuations
driving the quantum phase transitions is especially pronounced in
characterizing singularity and universality classes by universal scaling laws
[5, 6, 7, 8]. A superradiant phase transition (SPT) is a remarkable phenomenon
in the Dicke mode [9, 10, 11, 12, 13], which describes an ensemble of two-
level atoms interacting with a quantized single-mode cavity. Such SPT has been
found in few-body systems such as the quantum Rabi model in the limit where
the atomic transition frequency in the unit of the cavity frequency tends to
infinity [14, 15, 16, 17, 18, 19], exhibting the same critical exponents as
the Dicke model [14, 12, 13]. A few-body system of light-atom interactions
sheds new light in quantum simulation of quantum phase transitions due to its
high control and tunability.

To explore intriguing quantum phases, synthetic magnetic fields are applied in
atom-cavity coupling systems and photon lattices [20], chiral edge currents in
atoms [21, 22], chiral ground-state currents of interacting photons [23], and
fractional quantum Hall physics in the Jaynes-Cummings Hubbard lattice [24,
25]. In the presence of an artificial magnetic field, unusual superradiant
phases have been found in the generalized Rabi and Dicke systems, including a
chiral phase in a quantum Rabi triangle [26], chiral superradiant phases
analogies to quantum magnetism in a quantum Rabi ring [27] and anomalous
superradiant phases in a Dicke lattice model [28]. Such exotic SPTs exhibit
unusual scaling behaviors by comparing to the conventional SPT in the Dicke
model. Besides the artifical magnetic field, geometric frustration has been
proved to induce counterintuitive critical phenomena. Different critical
exponents of the excitation energy have been found on the two sides of phase
transitions in the Rabi triangle [27] and in the Dicke trimer [29, 28], which
are associated the frustrated geometry. However, there appears to be no clear
scaling exponents of quantum fluctuations near the critical value. Since there
is a challenge of an accurate solution beyond a mean-field approach for the
Rabi and Dicke lattice including interactions between cavities.

In this paper, we perform an analytical solution beyond the mean-field
approximation for the quantum Rabi triangle with an artificial magnetic filed.
In contrast to previous studies, we present an analytical expression of
quantum fluctuations and obtain accurate scaling exponents near the critical
point. In the chiral superradiant phase (CSP), we find the fluctuation of the
mean photons in the frustrated cavity diverges with an anomalous critical
exponent, which is distinct from the conventional exponent of the remaining
cavities. In the absence of an artificial magnetic field, the scaling exponent
of the frustrated cavity in the frustrated antiferromagnetic superradiant
phase (AFSP) is different from that in the CSP. Thus, the scaling behavior of
quantum fluctuations in the CSP and the frustrated AFSP fall into two classes,
one for the frustrated cavity and the other for the remaining cavities.
Moreover, the scaling behavior below and above the transition point is
different due to the frustrated local photons. In contrast, the scaling
exponents are equal to each other on the two sides of the transition of the
ferromagnetic superradiant phase (FSP), which are the same as that in the
Dicke model.

The rest of this paper is organized as following: In Sec. II, we introduce the
Hamiltonian of the quantum Rabi triangle and the dynamics of photons with an
artificial magnetic field. In Sec. III, the quantum fluctuations of the mean
photons and the position variance are given analytically in the normal phase
below the critical value. Sec. IV gives the analytical solutions for different
superradiant phases and the expressions for the quantum fluctuations. Sec. V
gives the scaling exponents of the quantum fluctuations and excitation
energies beyond the mean-field approximation. The conclusion is given in Sec.
VI.

##  II Hamiltonian

A Rabi triangle model consists of $N=3$ coupled cavities, and the Hamiltonian
is

$H_{\mathtt{QRT}}=\sum_{n}^{3}H_{R,n}+\sum_{n}^{3}J(e^{i\mathcal{\theta}}a_{n}^{\dagger}a_{n+1}+e^{-i\mathcal{\theta}}a_{n}a_{n+1}^{\dagger}),$
(1)

where the quantum Rabi Hamiltonian of each cavity coupled to a two-level atom
is desribed as $H_{R,n}=\omega
a_{n}^{\dagger}a_{n}+g\left(a_{n}^{\dagger}+a_{n}\right)\sigma_{n}^{x}+\frac{\Delta}{2}\sigma_{n}^{x}$.
The second term of the Hamiltonian $H_{\mathtt{QRT}}$ describes the photons
hopping between the nearest-neighbor cavities with a phase $\theta$, which can
be realized on cavities by dynamical modulating the cavity-atom coupling and
photon hopping strengths [26]. An artificial vector potential $A(r)$ leads the
photon hopping terms between nearby cavities $n$ and $m$ to become complex
with the phase given by $\theta=\int_{r_{n}}^{r_{m}}$ $A(r)dr$. The time-
reversal symmetry (TRS) is artificially broken when $\theta\neq
m\pi\,(m\in\mathbb{Z})$. The effective magnetic flux is $3\theta$ in three
cavities, which form a closed lopp.

To explore the effects of the artificial magnetic field we study the dynamics
of photon flow in our system. At $t=0$, a photon is prepared in the $1$st
cavity and the two-level atom in each cavity is in the down state, giving the
initial state
$|\varphi(0)\rangle=|100\rangle|\downarrow\downarrow\downarrow\rangle$. The
time-evolving wave function is
$|\varphi(t)\rangle=e^{-iH_{\mathtt{QRT}}t}|\phi_{0}\rangle|\varphi$. The mean
photon in the $i$-th cavity is given by
$N_{i}=\langle\varphi(t)|a_{i}^{\dagger}a_{i}|\varphi(t)\rangle$ .

Fig.1 shows the mean photon number evolution in each cavity dependent on
$\theta$. The photon number is initially occupied in the first cavity with
$N_{1}=1$. For $\theta=0$ in Fig.1(a), the photon transfer symmetrically from
cavity $1$ to cavity $2$ and cavity $3$ simultaneously, then back to cavity
$1$. There is no preferred circulation direction. Since TRS is preserved for
the trivial case $\theta=0$, the states at time $t$ and $T-t$ satisfy
$\varphi(t)=\varphi(T-t)$ with a time-evolving period $T$. A completely
different dynamics is observed for $\theta=\pi/2$ in Fig.1(b). The photon is
flowing unidirectionally, first from cavity $1$, to cavity $2$, to cavity $3$,
eventually back to cavity $1$. Such chiral current direction is a signature of
the breaking of TRS. Because the evolution of the state from $t=T$ backwards
is different as going forwards from $t=0$. Choosing $\theta=-\pi/2$ leads to
the opposite direction of the chiral photon flow in Fig.1(c). It demonstrate
that the artificial flux $\theta$ leads to the breaking of TRS, which behaves
similarly to a magnetic flux.

Figure 1: Mean photon number in each cavity $N_{1}$ (black solid line),
$N_{2}$ (red dashed line), $N_{3}$ (blue dotted line) as a function of time
$t$ for three values of (a)$\theta=0$, (b)$\pi/2$ and (c)$-\pi/2$ for the
scaled coupling strength $g_{1}=0.1$. In this paper we use $\Delta/\omega=50$
and $J/\omega=0.05$ by choosing $\omega=1$ as the unit for frequency.

Besides the TRS breaking induced by the artificial magnetic flux, rich
superradiant phase transitions have been explored dependent on $\theta$ in the
quantum Rabi triangle [26, 27]. The phase transitions occur in the limit where
$\Delta$ is much larger than frequency scales in the system
$\Delta/\omega\rightarrow\infty$. Different from previous studies, we extract
unusual critical exponents to classify the universality classes of the
nontrivial phase transition by involving the geometric frustrations and the
effects of the magnetic flux. In particular, quantum fluctuations are crucial
to capture the singularity and scaling exponents near the critical point using
an accurate solution beyond a mean-field approximation.

##  III Quantum fluctuations in the normal phase

We perform a Schrieffer-Woff transformation with an unitary operator
$S_{n}=\exp[-i\sigma_{n}^{y}g_{1}\sqrt{\omega/\Delta}\left(a_{n}^{\dagger}+a_{n}\right)]$.
The lower-energy Hamiltonian is obtained by projecting $H_{\mathtt{iCP}}$ to
the spin subspace $|\downarrow\rangle$, giving

$\displaystyle H_{\mathtt{NP}}^{\downarrow}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{3}(\omega-2\omega
g_{1}^{2})a_{n}^{\dagger}a_{n}-\omega g_{1}^{2}(a_{n}^{2}+a_{n}^{\dagger 2})$
(2)
$\displaystyle+J\sum_{<nn^{\prime}>}^{3}(e^{i\theta}a_{n}^{\dagger}a_{n^{\prime}}+h.c.)+E_{0},$

where the energy constant is $E_{0}=3[-\Delta/2-\omega
g_{1}^{2}+(\omega+J)\omega^{2}g_{1}^{2}/\Delta]$.

By introducing the Fourier transformation
$a_{n}^{\dagger}=\sum_{q}e^{inq}a_{q}^{\dagger}/\sqrt{N}$ with the quasi-
momentum $q=0,\pm 2\pi/3$, the transformed Hamiltonian is expressed as
$H_{\mathtt{iCP}}^{\downarrow}=\sum_{q}\omega_{q}a_{q}^{\dagger}a_{q}-\omega
g_{1}^{2}(a_{q}a_{-q}+a_{q}^{\dagger}a_{-q}^{\dagger})+E_{0}$, where
$\omega_{q}=\omega-2\omega g_{1}^{2}+2J\cos(\theta-q)$. By performing a
unitary transformation
$S_{q}=\exp[\lambda_{q}(a_{q}^{\dagger}a_{-q}^{\dagger}-a_{q}a_{-q})]$ with a
variational squeezing parameter
$\lambda_{q}=-\frac{1}{8}\ln\frac{\omega_{q}+\omega_{-q}-4\omega
g_{1}^{2}}{\omega_{q}+\omega_{-q}+4\omega g_{1}^{2}}$, we derive the
Hamiltonian in diagonalized form as
$H_{\mathtt{iCP}}^{\downarrow}=\sum_{q}\varepsilon_{q}a_{q}^{\dagger}a_{q}+E_{g}$
with the ground-state energy
$E_{g}=E_{0}+\frac{1}{2}\sum_{q}(\varepsilon_{q}-\omega_{q})$. The excitation
spectrum is obtained as

$\varepsilon_{q}^{\mathtt{NP}}=\frac{1}{2}(\sqrt{\Omega_{+,q}^{2}-16\omega^{2}g_{1}^{4}}+\Omega_{-,q})$
(3)

with the dispersion $\Omega_{\pm,q}=\omega_{q}\pm\omega_{-q}$. The lowest
excitation energy is associated with the qusi-momentum $q$ dependent on
$\theta$. For $\theta=0$, the excitation energy $\varepsilon_{q=\pm
2/3}^{\mathtt{NP}}$ with $q=\pm 2\pi/3$ becomes the lowest one, and decreases
to zero for $g_{1}$ approaching $g_{c}$ from below in Fig. 2(a). As $\theta$
increases, the lowest energy is determined by $q=-2\pi/3$ in Fig. 2(b). Then
it changes to $\varepsilon_{q=0}^{\mathtt{NP}}$ with $q=0$ in Fig. 2(d). The
vanishing of $\varepsilon_{q}^{\mathtt{NP}}=0$ dependent on $\theta$ and
$q$-momentum gives the critical scaled coupling strength
$g_{1c}(q)=\frac{1}{2}\sqrt{\frac{1+4J/\omega\cos\theta\cos
q+4J^{2}/\omega^{2}\cos(\theta+q)\cos(\theta-q)}{1+2J/\omega\cos\theta\cos
q}}$, which gives the phase boundary. Especially for $\theta=0$, the critical
coupling strength is $g_{c}(\pm 2\pi/3)$ originating from the lowest
excitation energy with $q=\pm 2\pi/3$. By increasing the magnetic flux, the
critical flux is given by $\theta_{c}=\pm
cos^{-1}[-2J/(\sqrt{8J^{2}+\omega^{2}}+\omega)]$ by solving $g_{c}(\pm
2\pi/3)=g_{c}(0)$, which classifies phase boundary by tuning $\theta$.

The ground-state in the NP is

$|\varphi_{np}\rangle=\prod_{q}e^{\lambda_{q}(a_{q}^{\dagger}a_{-q}^{\dagger}-a_{q}a_{-q})}|0\rangle_{q}|\downarrow\rangle,$
(4)

where $|\downarrow\rangle$ is the lowest state of the atom. The average photon
number in the ground state can be obtained as $\langle
a_{q}^{\dagger}a_{q}\rangle_{np}=\langle\varphi_{np}|a_{q}^{\dagger}a_{q}|\varphi_{np}\rangle=(\cosh
4\lambda_{q}-1)/2$, which desribes the quantum fluctuation of photons for
$g_{1}$ approaching the critical value $g_{1c}$ from below. The local mean
photon number in the $n$th cavity is obtained by

$\langle a_{n}^{\dagger}a_{n}\rangle_{np}=\frac{1}{N}\sum_{q}\langle
a_{q}^{\dagger}a_{q}\rangle_{np}=\frac{1}{2N}\sum_{q}[\frac{\Omega_{+,q}}{2\varepsilon_{q}-\Omega_{-,q}}-1].$
(5)

The variance of $x_{q}=a_{q}+a_{q}^{\dagger}$ and
$p_{q}=i(a_{q}^{\dagger}-a_{q})$ are derived as $(\Delta
x_{q})^{2}=\langle\varphi_{np}|x_{q}^{2}|\varphi_{np}\rangle-\langle\varphi_{np}|x_{q}|\varphi_{np}\rangle^{2}=e^{4\lambda_{q}}$,
and $(\Delta p_{q})^{2}=\langle p_{q}^{2}\rangle-\langle
p_{q}\rangle^{2}=e^{-4\lambda_{q}}$. The local variance of quadratures in the
$n$th cavity is

$\displaystyle(\Delta x_{n})^{2}=\frac{1}{N}\sum_{q}(\Delta
x_{q})^{2}=\frac{\Omega_{+,q}+4\omega
g_{1}^{2}}{2\varepsilon_{q}-\Omega_{-,q}},$ (6)

and $(\Delta p_{n})^{2}=\frac{1}{N}\sum_{q}e^{-4\lambda_{q}}$.

It indicates that the singularity of the mean photon number $\langle
a_{n}^{\dagger}a_{n}\rangle_{np}$, the variance of $(\Delta x_{n})^{2}$ and
$(\Delta p_{n})^{2}$ is determined by the excitation energy $\varepsilon_{q}$
and $\Omega_{-,q}$ in the denominator. Obviously, $\varepsilon_{q}$ becomes
zero at the critical value $g_{1c}$. It is interesting to understand the
divergence from $\Omega_{-,q}=4J\texttt{sin}\theta\texttt{sin}q$ dependent on
$q$ and $\theta$. In the absence of the magnetic field with $\theta=0$ or the
momentum $q=0$, one has $\Omega_{-,q}=0$. It leads to the singularity of
quantum fluctuations of $\langle a_{n}^{\dagger}a_{n}\rangle_{np}$ at the
critical point in Fig. 3(a) (d) as well as $(\Delta x_{n})^{2}$ in Fig.
4(a)(d). In contrast, for $\theta\neq 0$ and $q=\pm 2\pi/3$, it gives nonzero
value of $\Omega_{-,q}$, which results in a finite value of $\langle
a_{n}^{\dagger}a_{n}\rangle_{\texttt{NP}}$ and $(\Delta
x_{n})^{2}_{\texttt{NP}}$, respectively. There appear non-divergent
fluctuations in Fig. 3(b) and Fig. 4(b) for $g_{1}$ approaching $g_{1c}$ from
below.

Figure 2: The lowest and second excitation energy $\varepsilon_{1(2)}$ below
and above the critical coupling strength $g_{1c}$ for the NP-AFSP transition
($\theta=0$)(a), the NP-CSR transition ($\theta=0.1<\theta_{c}$)(b), the
triple point ($\theta=\theta_{c}$)(c) and the NP-FSP transition
($\theta=1.7>\theta_{c}$) (d). The parameters are $\Delta/\omega=100$ and
$J/\omega=0.05$ by setting $\omega=1$ as the unit for the frequency. The inset
(b) shows the second excitation energy $\varepsilon_{2}$ for the NP-CSP
transition, which exhibits an obvious energy gap in comparison to
$\varepsilon_{1}$.

##  IV Quantum fluctuations in the superadiant phase

It is interesting to explore the quantum fluctuations and scaling exponents
when $g_{1}$ approaches $g_{1c}$ from above. As the atom-cavity coupling
increases $g_{1}>g_{1c}$, rich superradiant phases emerge by adjusting the
flux $\theta$ [27, 26]. Since the photon population in each cavity becomes
macroscopic. The bosonic operator $a_{n}^{\dagger}$ $\left(a_{n}\right)$ is
expected to be shifted as
$\tilde{a}_{n}=D^{\dagger}(\alpha_{n})a_{n}D(\alpha_{n})=a_{n}+\alpha_{n}$
with the complex displacement $\alpha_{n}=A_{n}+iB_{n}$. It is different from
the mean field approximation, for which the bosonic operator $a_{n}$ is
replaced by its mean value $\langle a_{n}\rangle=\alpha_{n}$. The Hamiltonian
in the superradiant phases becomes

$\displaystyle H_{\text{SR}}^{\downarrow}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{3}\omega\tilde{a}_{n}^{\dagger}a_{n}-\frac{\Delta_{n}}{2}\tau_{n}^{z}+\lambda_{n}\left(\tilde{a}_{n}^{\dagger}+\tilde{a}_{n}\right)\tau_{n}^{x}$
(7)
$\displaystyle+Ja_{n}^{\dagger}(e^{i\theta}\tilde{a}_{n+1}+e^{-i\theta}\tilde{a}_{n-1})+V_{\text{off}}+E_{0},$

where the transformed Pauli matrix is
$\tau_{n}^{z}=\Delta/\Delta_{n}\sigma_{n}^{z}+4gA_{n}/\Delta_{n}\sigma_{n}^{x}$
with the renormalized parameter
$\Delta_{n}=\sqrt{\Delta^{2}+16g^{2}A_{n}^{2}}$, and the effective coupling
strength is $\lambda_{n}=g\Delta/\Delta_{n}$. The off-diagonal term is
expressed as
$V_{\text{off}}=\sum_{n}\omega(\alpha_{n}a_{i}^{\dagger}+\alpha_{n}^{\ast}a_{n})+g\left(a_{n}^{\dagger}+a_{n}\right)\sin(2\gamma_{n})\sigma_{n}^{z}+J[a_{n}^{\dagger}(e^{i\theta}\alpha_{n+1}+e^{-i\theta}\alpha_{n-1})+h.c.]$.
By eliminating $V_{\text{off}}$ term, we obtain equations explicitly

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\omega
A_{n}-g\sin(2\gamma_{n})+J[(A_{n+1}+A_{n-1})\cos\theta$ (8) $\displaystyle+$
$\displaystyle(B_{n-1}-B_{n+1})\sin\theta]=0,$

and

$0=\omega B_{n}+J[\cos\theta(B_{n+1}+B_{n-1})+\sin\theta(A_{n+1}-A_{n-1})].$
(9)

$\alpha_{n}$ can be accurately obtained by solving above equations. The
lowest-energy Hamiltonian is obtained by projecting to the spin subspace
$|\downarrow\rangle$, giving

$\displaystyle H_{\text{eff}}^{\downarrow}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{3}\omega
a_{n}^{\dagger}a_{n}-\frac{\lambda_{n}^{2}}{\Delta_{n}}\left(a_{n}^{\dagger}+a_{n}\right)^{2}$
(10)
$\displaystyle+Ja_{n}^{\dagger}(e^{i\theta}a_{n+1}+e^{-i\theta}a_{n-1})+E_{g}.$

where the ground-state energy is
$E_{g}=\sum_{n}\omega\alpha_{n}^{\ast}\alpha_{n}+J\sum_{n}\alpha_{n}^{\ast}(e^{i\theta}\alpha_{n+1}+e^{-i\theta}\alpha_{n-1})-\Delta_{n}/2$.

Superradiant phases can be characterized by the order parameter $\alpha_{n}$.
Firstly, for $\theta=0$ and $g_{1}>g_{1c}^{\text{AFSP}}(\pm 2\pi/3)$, the
system is in the frustrated AFSP. $\alpha_{n}$ is real with opposite signs for
the neighboring cavities. It corresponds to an AFSP with antiferromagnetic
order. Moreover, geometric frustration in three cavities yields site-dependent
$\alpha_{n}$, so-called frustrated AFSP. The ground state breaks the $C_{3}$
symmetry. Secondly, for $0<\theta<\theta_{c}$, the system enters the CSP
regime above the critical coupling strength $g_{c}^{\text{CSP}}(-2\pi/3)$
[26]. $\alpha_{n}$ is complex and $n$ dependent, which is different from that
in the AFSP. For $\theta_{c}<\theta\leq\pi$ the displacement in each cavity
$\alpha_{n}$ is real and is independent of $n$, which gives
$\alpha_{n}=\pm\sqrt{16g^{4}/(\omega+2Jcos\theta)^{2}-\Delta^{2}}/(4g)$. It
forms a ferromagnetic order with the same displacement of the neighboring
cavities for $g_{1}>g_{1c}^{\text{FSP}}(q=0)$. It corresponds to the FSP.

The Hamiltonian in Eq. (10) is bilinear in the creation and annihilation
operators $a_{i}^{\dagger}$ and $a_{i}$, which captures the quantum
fluctuations especially near the critical coupling strength $g_{c}$. It can be
diagonalized by the bosonic Bogoliubov transformation. By using the denotation
$\alpha=\\{a_{1},a_{2},a_{3},a_{1}^{\dagger},a_{2}^{\dagger},a_{3}^{\dagger}\\}$,
the Hamiltonian in Eq. (10) can be written in matrix form as
$H_{\text{eff}}^{\downarrow}=\alpha
M\alpha^{\dagger}-3(\omega-\lambda_{n}^{2}/\Delta_{n})/2$, where a transformed
matrix $M$ is

$M=\left(\begin{array}[]{cccccc}\omega/2-\lambda_{1}^{2}/\Delta_{1}&Je^{-i\theta}/2&Je^{i\theta}/2&-\lambda_{1}^{2}/\Delta_{1}&0&0\\\
Je^{i\theta}/2&\omega/2-\lambda_{2}^{2}/\Delta_{2}&Je^{-i\theta}/2&0&-\lambda_{2}^{2}/\Delta_{2}&0\\\
Je^{-i\theta}/2&Je^{i\theta}/2&\omega/2-\lambda_{3}^{2}/\Delta_{3}&0&0&-\lambda_{3}^{2}/\Delta_{3}\\\
-\lambda_{1}^{2}/\Delta_{1}&0&0&\omega/2-\lambda_{1}^{2}/\Delta_{1}&Je^{i\theta}/2&Je^{-i\theta}/2\\\
0&-\lambda_{2}^{2}/\Delta_{2}&0&Je^{-i\theta}/2&\omega/2-\lambda_{2}^{2}/\Delta_{2}&Je^{i\theta}/2\\\
0&0&-\lambda_{3}^{2}/\Delta_{3}&Je^{i\theta}/2&Je^{-i\theta}/2&\omega/2-\lambda_{3}^{2}/\Delta_{3}\end{array}\right).$
(11)

We perform a Bogoliubov’s transformation to give bosonic operators
$\beta=\\{b_{1}^{\dagger},b_{2}^{\dagger},b_{3}^{\dagger},b_{1},b_{2},b_{3}\\}$
as a linear combination of
$\alpha=\\{a_{1}^{\dagger},a_{2}^{\dagger},a_{3}^{\dagger},a_{1},a_{2},a_{3}\\}$,
which satisfies $\alpha^{\dagger}=T\beta^{\dagger}$ with a paraunitary matrix
$T$ rather than a unitary matrix. To ensure the bosonic commutation relation,
the paraunitary matrix $T$ satisfies the relations of $T^{\dagger}\Lambda
T=T\Lambda T^{\dagger}=\Lambda$, where $\Lambda=$
$\left(\begin{array}[]{cc}I_{3\times 3}&0\\\ 0&-I_{3\times
3}\end{array}\right)$ with $I_{3\times 3}$ being the identity matrix of order
$3$. Substituting for $\alpha$ and $\alpha^{\dagger}$ in terms of
$\beta^{\dagger}$ and $\beta$, one obtains the diagonalized form as

$H_{\text{eff}}^{\downarrow}=\beta
T^{\dagger}MT\beta^{\dagger}=2\sum_{k=1}^{3}\varepsilon_{k}b_{k}^{\dagger}b_{k}+\frac{\varepsilon_{k}-\omega}{2},$
(12)

where $\varepsilon_{k}$ is the corresponding eigenvalue for each excitation
mode. The eigenvalues $\pm\varepsilon_{n}$ are obtained by diagonalizing the
matrix $\Lambda M$ as $T^{-1}\Lambda MT=\Lambda\varepsilon$. The corresponding
eigenvector gives the $k$-th column vector of the paraunitary matrix as
$T_{k}=[T_{1k},T_{2k},\cdots,T_{6k}]^{T}$.

The ground state in the superradiant phases is the vacuum state

$|\varphi_{\texttt{SR}}\rangle=|0\rangle_{b_{1}}|0\rangle_{b_{2}}|0\rangle_{b_{3}},$
(13)

where $|0\rangle_{b_{n}}$ satisfies $b_{n}|0\rangle_{b_{n}}=0$. Note that the
operator $b_{n}$ corresponds to the original operator $a_{n}$ with the
displacement transformation $D(\alpha)$ and the transformation
$\alpha^{\dagger}=T\beta^{\dagger}$. One obtains
$a_{n}=\sum_{i=1}^{3}T_{n,i}b_{i}+T_{n,i+3}b_{i}^{\dagger}+\alpha_{n}$. The
local photon number in the $n$th cavity of the ground state
$|\varphi_{\texttt{SR}}\rangle$ is derived explicitly as

$\displaystyle\langle
a_{n}^{\dagger}a_{n}\rangle=|T_{n4}|^{2}+|T_{n5}|^{2}+|T_{n6}|^{2}+|\alpha_{n}|^{2}.$
(14)

The variance of $x_{n}$ and $p_{n}$ are

$(\Delta
x_{n})^{2}=|T_{n1}+T_{n4}^{\ast}|^{2}+|T_{n2}+T_{n5}^{\ast}|^{2}+|T_{n3}+T_{n6}^{\ast}|^{2},$
(15)

and $(\Delta
p_{n})^{2}=|T_{n4}^{\ast}-T_{n1}|^{2}+|T_{n5}^{\ast}-T_{n2}|^{2}+|T_{n6}^{\ast}-T_{n3}|^{2}$,
which are derived in detail in Appendix. Thus the quantum fluctuations of
$(\Delta x_{n})^{2}$ and $\langle a_{n}^{\dagger}a_{n}\rangle$ are obtained
beyond the mean-field approaximation. We calculate the scaling behaviors of
the critical fluctuations in the vicinity of the critical value $g_{c}$ in the
following.

##  V Unusual critical exponents

Our major interests are the scaling exponents of the quantum fluctuations as
$g_{1}$ approaches $g_{c}$ for the second-order phase transitions from the NP
to different superradiant phases. Generally, the energy gap which is measured
by the lowest excitation energy $\varepsilon_{1}$ vanishes as $g_{1}$
approaches $g_{1c}$

$\varepsilon_{1}\propto|g_{1}-g_{1c}|^{\gamma}.$ (16)

Here the critical exponent $\gamma=z\nu$ is usually universal [6], which is
independent of most of the parameters of the Hamiltonian. In addition to the
vanishing energy gap, a length scale which measures the correlations at the
longest distances becomes diverge at the critical value. The variance of
position quadrature plays an analogous role of the diverging length scale,
which diverges as

$\Delta x\propto|g_{1}-g_{1c}|^{-\nu}$ (17)

with a critical exponent $\nu$. It leads to the dynamical critical exponent
$z=\gamma/\nu$. Meanwhile, the fluctuations of the local mean photons near the
critical value diverges as

$\langle a_{n}^{\dagger}a_{n}\rangle\propto|g_{1}-g_{1c}|^{-\beta},$ (18)

where $\beta$ is a critical exponent. Using the Bogoliubov’s diagonalization
method, we calculate the scaling exponents of the lowest excitation energy
$\varepsilon_{1}$, the fluctuation of position variance $(\Delta x)^{2}$ in
Eq.(15) and photons $\langle a_{n}^{\dagger}a_{n}\rangle$ in Eq.(14).

Figure 3: Quantum fluctuations of the local mean photons for the $n$th cavity
$a^{\dagger}_{n}a_{n}$ as a function of $g_{1}/g_{1c}$ for the NP-AFSP
transition ($\theta=0$)(a), the NP-CSP transition ($\theta=0.1<\theta_{c}$)
(b), the triple point ($\theta=\theta_{c}$) (c), and the NP-FSP transition
($\theta=1.7>\theta_{c}$) (d). The inset (a) shows the scaling behavior of
$a^{\dagger}_{2(3)}a_{2(3)}$. The right side of the inset (b) shows the
behavior of $a^{\dagger}_{1}a_{1}$, and the left side shows the scaling
exponents. The inset (c) shows the scaling results by the solution in the CSP,
which is same as one obtained by the solution in the FSP. The inset (d) shows
the exponent .

For the second-order phase transition between the NP and the frustrated AFSP
with $\theta=0$, Fig. 2(a) shows $\varepsilon_{1}$ closing the gap at $g_{1c}$
with the exponent $1/2$ for $g_{1}<g_{1c}$. In contrast, two excitation
energies $\varepsilon_{1}$ and $\varepsilon_{2}$ vanish at the critical point
for $g_{1}>g_{1c}$

$\displaystyle\varepsilon_{1,\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1},\varepsilon_{2,\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1/2}.$
(19)

It exhibits two different exponents $\gamma_{+}$ ($\gamma_{-}$) at the two
sides of phase transition for $g_{1}>g_{1c}$ ($g<g_{1c}$), which is
distinguished from the conventional second-order phase transitions. The
universal exponent $\gamma_{-}=1/2$ is the same as a mean-field exponent of
the Rabi model [14]. While two different exponents $\gamma_{+}=1/2$ and
$\gamma_{+}=1$ in Eq.(19) appear for $g_{1}>g_{1c}$ , for which the latter is
the unconventional scaling exponent beyond the mean-field one. The distinct
scaling behavior are consistent with results in the Dicke trimer [29].
Meanwhile, Fig. 3(a) shows the fluctuations of the mean photons in each cavity
diverging with the same exponent $\beta_{-}=1/2$ for $g_{1}<g_{1c}$. However,
for $g_{1}>g_{1c}$, the photon number diverges locally dependent on $n$, which
is associated with the frustrated geometry. $\langle
a_{n}^{\dagger}a_{n}\rangle$ for the $n$th cavity exhibits different scaling
law as

$\displaystyle\langle
a_{1}^{\dagger}a_{1}\rangle_{\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1},$
$\displaystyle\langle
a_{2(3)}^{\dagger}a_{2(3)}\rangle_{\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1/2}.$
(20)

The scaling of the frustrated cavity ($n=1$) yields an unusual exponent
$\beta=1$ , while the remaining cavities ($n=2,3$) diverges with the same
exponent $\beta=1/2$ as that in the Rabi model, respectively. The similar
$n$-dependent fluctuation of $(\Delta x_{n})^{2}$ is shown in Fig. 4(a). The
scaling exponent of $\Delta x_{n}$ of each cavity below the critical value
$g_{1c}$ is the same $\nu_{-}=1/4$ for $g_{1}<g_{1c}$. By contrast, for
$g_{1}>g_{1c}$ the fluctuation diverges respectively

$\displaystyle(\Delta x_{1})^{2}_{\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1},$
$\displaystyle(\Delta
x_{2(3)})^{2}_{\texttt{AFSP}}\propto|g_{1}-g_{1c}|^{-1/2}.$ (21)

The scaling exponents are extracted $\nu_{+}=1/2$ for the frustrated cavity
and $\nu_{+}=1/4$ for the remaining cavities.

Figure 4: Quantum fluctuation of the variance of position for the $n$th cavity
$(\Delta x_{n})^{2}$ as a function of $g_{1}/g_{1c}$ for the NP-AFSP
transition ($\theta=0$)(a), the NP-CSP transition ($\theta=0.1<\theta_{c}$)
(b), the triple point ($\theta=\theta_{c}$) (c), and the NP-FSP transition
($\theta=1.7>\theta_{c}$) (d). The inset (a) shows the scaling behavior of
$(\Delta x_{2(3)})^{2}$. The inset (b) shows the fluctuations of the first
frustrated cavity $(\Delta x_{1})^{2}$ above $g_{1c}$. The inset (c) shows the
scaling results by the solution in the CSP, which is same as one obtained by
the solution in the FSP.

For the NP-CSP transition with $0<\theta<\theta_{c}$, the unusual mode
$\varepsilon_{1}$ closes the gap for $g_{1}$ approaching $g_{1c}$ from below
or above in Fig. 2(b). It behaves as

$\displaystyle\varepsilon_{1,\texttt{NP}}\propto|g_{1}-g_{1c}|^{-1}.\varepsilon_{1,\texttt{CSP}}\propto|g_{1}-g_{1c}|^{-3/2}.$
(22)

It gives two unconventional exponents $\gamma_{-}=1$ and $\gamma_{+}=3/2$
below and above $g_{1c}$, which are different from results in the NP-AFSR
transition. In contrast to the divergence at $g_{1c}$, one observes a finite
value of $\langle a_{n}^{\dagger}a_{n}\rangle_{\texttt{NP}}$ and $(\Delta
x_{n})_{\texttt{NP}}$ when $g_{1}$ approaches $g_{1c}$ from below in Fig. 3
(b) and Fig. 4(b). The results are consistent with the analysis from Eq.(5)
and (6), for which the denominator is nonzero at the critical point for $q=\pm
2\pi/3$. For $g_{1}>g_{1c}$, Fig. 3(b) shows $n$-dependent scaling behaviors
of the local photon number

$\displaystyle\langle
a_{1}^{\dagger}a_{1}\rangle_{\texttt{CSP}}\propto|g_{1}-g_{1c}|^{-1/3},$
$\displaystyle\langle
a_{2(3)}^{\dagger}a_{2(3)}\rangle_{\texttt{CSP}}\propto|g_{1}-g_{1c}|^{-1/2}.$
(23)

$\langle a_{1}^{\dagger}a_{1}\rangle$ of the frustrated cavity diverges with
the exponents $1/3$, which is different from that in the AFSP. The anomalous
exponent $\nu_{+}=1/3$ of the frustrated cavity is more accurate by comparing
to that in the Dicke lattice [28]. Meanwhile, Fig. 4(b) shows the fluctuation
of the position variance of the $n$th cavity, which diverges as

$\displaystyle(\Delta x_{1})^{2}_{\texttt{CSP}}\propto|g_{1}-g_{1c}|^{-1/3},$
$\displaystyle(\Delta
x_{2(3)})^{2}_{\texttt{CSP}}\propto|g_{1}-g_{1c}|^{-1/2}.$ (24)

The unusual exponent $\nu_{+}=1/6$ of the frustrated cavity is also different
from $\nu_{+}=1/2$ in the AFSP transition. Since the CSR phase transition is
associated with both of the geometry frustration and the artificial magnetic
field. It exhibits distinct exponents of the frustrated cavity by comparing to
the frustrated AFSP transition with $\theta=0$. They belong to different
universality classes of phase transitions, which are distinguished from the
conventional superradiant phase transition in the Dicke and Rabi model.

At the triple point $\theta_{c}$, the CSR phase, FSP and NP coexist. There
appear two excitation modes closing the gaps at $g_{1c}$. Fig. 2(c) show two
different power laws behaviors

$\displaystyle\varepsilon_{1,\texttt{TP}}\propto|g_{1}-g_{1c}|^{-1},\varepsilon_{2,\texttt{TP}}\propto|g_{1}-g_{1c}|^{-1/2}.$
(25)

We obtain exponents $\gamma_{\pm}=1$ and $\gamma_{\pm}=1/2$. It is a signature
of the coexistence of both CSP and FSP. And scaling exponents are the same at
the two sides of the TP. We calculate the the fluctuation of $\langle
a_{n}^{\dagger}a_{n}\rangle$ by the analytical solutions in the FSP phase in
Fig. 3(c), which behaves the same as the solution in the CSR in the inset.
Fig. 4(c) show the scaling laws of $(\Delta x_{n})^{2}_{\texttt{TP}}$. We
obtain the scaling behaviors near the critical point

$\displaystyle\langle
a_{n}^{\dagger}a_{n}\rangle_{\texttt{TP}}\propto|g_{1}-g_{1c}|^{-1/2},(\Delta
x_{n})^{2}_{\texttt{TP}}\propto|g_{1}-g_{1c}|^{-1/2}.$ (26)

Thus, we obtain the scaling exponents $\beta=1/2$ and $\nu=1/4$ for the
fluctuations of photon number and the position variance at the triple point,
respectively.

For the NP-FSP transition, the scaling laws holds with the same value of the
exponent $\gamma$ both for $g_{1}>g_{1c}$ and $g<g_{1c}$. The exponent of the
excitation energy at the two sides of the phase transition are the same
$\gamma_{+}=\gamma_{-}=1/2$ in Fig. 2(d). And the mean photons of the $n$th
cavity diverges with the same exponent $\beta_{\pm}=1/2$ in Fig. 3(d). Fig.
3(d) show the exponent $\nu_{\pm}=1/4$ for the fluctuation of the position
variance.

Various critical exponents for $g_{1}$ approaching $g_{1c}$ from below and
above are listed in Table I for different superradiant phase transitions.
Scaling exponents below and above $g_{1c}$ are different for the NP-AFSP and
NP-CSP transition, while are the same for the NP-FSP transition. For the
special frustrated cavity, the unusual critical exponent of the photon
fluctuations are $\beta_{+}=1$ and $\beta_{+}=1/3$ for the NP-AFSP and NP-CSP
transition, respectively. Meanwhile, the anomalous exponents of the
fluctuations of the position variance are $\nu_{+}=1/2$ and $\nu_{+}=1/6$.
They are distinguished from the exponents of the remaining cavities with
$\beta_{+}=1/2$ and $\nu_{+}=1/4$, which are the same exponents as that in the
conventional quantum Rabi model. It demonstrates that the artificial magnetic
field and geometric fluctuations yield unusual exponents, which predict
nontrivial universality classes of superradiant phase transitions.

Table 1:  Scaling exponents $\gamma_{\pm}$ of the excitation energy $\varepsilon_{n}$, $\nu_{\pm}$ and $\beta_{\pm}$ for quantum fluctuations of the $\Delta x_{n}$ and $\langle a_{n}^{\dagger}a_{n}\rangle$ on two sides of phase transitions of the NP-AFSP($\theta=0$), NP-CSP ($\theta<\theta_{c}$), and NP-FSP ($\theta>\theta_{c}$) as well as the triple point (TP) ($\theta=\theta_{c}$). |  $\gamma_{-}$ |  $\gamma_{+}$ |  $\nu_{-}$ |  $\nu_{+}$ |  $\beta_{-}$ |  $\beta_{+}$
---|---|---|---|---|---|---
NP-AFSP |  $1/2$ |  $1,1/2$ |  $1/4$ |  $1/2,1/4$ |  $1/2$ |  $1,1/2$
NP-CSP |  $1$ |  $3/2$ |  $/$ |  $1/6,1/4$ |  $/$ |  $1/3,1/2$
TP |  $1,1/2$ |  $1,1/2$ |  $1/4$ |  $1/4$ |  $1/2$ |  $1/2$
NP-FSP |  $1/2$ |  $1/2$ |  $1/4$ |  $1/4$ |  $1/2$ |  $1/2$

##  VI  Conclusion

we have demonstrated quantum fluctuations and scaling behaviors by using an
analytical approach and a Bogoliubov transformation in the Rabi triangle with
an artificial magnetic field. The time-reversal symmetry breaking is observed
from the photon population dynamics, exhibiting the effects of an artificial
magnetic field. Due to the frustrations of the triangular geometry, there
emerges two scaling laws, one for the frustrated cavity with an unusual
exponent and the other for the remaining cavities with the same exponent as
the conventional Dicke model. For the fluctuations of the local photon number
and the position variance, the scaling exponents of the frustrated cavity are
$\beta_{+}=1/3$ and $\nu_{+}=1/6$ for the CSP transition, which are distinct
from $\beta_{+}=1$ and $\nu_{+}=1/2$ for the frustrated AFSP transition
without an artificial magnetic flux. Moreover, the scaling exponents are
different on two sides of the phase transition from the NP to the CSP and
AFSP. The unconventional critical exponents predict different universality
classes beyond the conventional single-cavity Rabi model. Therefore, our work
paves a way for exploring unconventional phase transitions in the few-body
light-matter interacting system.

###### Acknowledgements.

The authors thank Qing-Hu Chen and Xiang-You Chen for useful discussions. This
work was supported by NSFC under Grant No.12075040 and No. 12347101.

##  Appendix A Derivations of quantum fluctuations

In the superradiant phase, the operator $a_{n}$ is shifted as
$a_{n}+\alpha_{n}$. Using the Bogoliubov’s diagonalization method, the ground
state of the photon part is obtained as
$|\varphi_{\texttt{SR}}\rangle=|0\rangle_{b_{1}}|0\rangle_{b_{2}}|0\rangle_{b_{3}}$
with $b_{n}|0\rangle_{b_{n}}=0$. With the transformation $D(\alpha)$ and the
transformation $\alpha^{\dagger}=T\beta^{\dagger}$, the operator $a$ becomes
$a_{n}=\sum_{i=1}^{3}T_{n,i}b_{i}+T_{n,i+3}b_{i}^{\dagger}+\alpha_{n}$. The
expected value of $x_{n}=(a_{n}+a_{n}^{\dagger})$ is given by

$\displaystyle\langle x_{n}\rangle$ $\displaystyle=$ ${}_{b_{n}}\langle
0|(a_{n}+a_{n}^{\dagger})|0\rangle_{b_{n}}=\alpha_{n}+\alpha_{n}^{*}$ (27)
$\displaystyle+$ ${}_{b_{n}}\langle
0|\sum_{i=1}^{3}T_{n,i}b_{i}+T_{n,i+3}b_{i}^{\dagger}+T^{*}_{n,i}b_{i}^{\dagger}+T_{n,i+3}^{*}b_{i}|0\rangle_{b_{n}}$
$\displaystyle=$ $\displaystyle\alpha_{n}^{\ast}+\alpha_{n}.$

Then we derive the mean value of $x_{n}^{2}$ as

$\displaystyle\langle x_{n}^{2}\rangle$ $\displaystyle=$ ${}_{b_{n}}\langle
0|(a_{n}+a_{n}^{\dagger}+\alpha_{n}^{\ast}+\alpha_{n})^{2}|0\rangle_{b_{n}}$
$\displaystyle=$
$\displaystyle\langle(a_{n}+a_{n}^{\dagger})^{2}+(\alpha_{n}^{\ast}+\alpha_{n})(a_{n}+a_{n}^{\dagger})\rangle+(\alpha_{n}^{\ast}+\alpha_{n})^{2}$
$\displaystyle=$
$\displaystyle\langle\varphi_{\texttt{SR}}|[(T_{n1}+T_{n4}^{\ast})b_{1}+(T_{n2}+T_{n5}^{\ast})b_{2}$
$\displaystyle+(T_{n3}+T_{n6}^{\ast})b_{3}+h.c.]^{2}|\varphi_{\texttt{SR}}\rangle+(\alpha_{n}^{\ast}+\alpha_{n})^{2}$
$\displaystyle=$
$\displaystyle|T_{n1}+T_{n4}^{\ast}|^{2}+|T_{n2}+T_{n5}^{\ast}|^{2}$
$\displaystyle+|T_{n3}+T_{n6}^{\ast}|^{2}+(\alpha_{n}^{\ast}+\alpha_{n})^{2}.$

The expected value of the variance of $\Delta x_{n}$ is obtained as

$\displaystyle(\Delta x_{n})^{2}$ $\displaystyle=\langle
x_{n}^{2}\rangle-\langle x_{n}\rangle^{2}$
$\displaystyle=|T_{n1}+T_{n4}^{\ast}|^{2}+|T_{n2}+T_{n5}^{\ast}|^{2}+|T_{n3}+T_{n6}^{\ast}|^{2}.$

Additionally, the expected value of the momentum quadrature
$p_{n}=i(a_{n}^{\dagger}-a_{n})$ is

$\langle p_{n}\rangle=_{b_{n}}\langle
0|i(a_{n}^{\dagger}-a_{n}+\alpha_{n}^{\ast}-\alpha_{n})|0\rangle_{b_{n}}=i(\alpha_{n}^{\ast}-\alpha_{n}).$
(29)

The expected value of $p_{n}^{2}$ is derived as

$\displaystyle\langle p_{n}^{2}\rangle$ $\displaystyle=$ ${}_{b_{n}}\langle
0|-(a_{n}^{\dagger}-a_{n}+\alpha_{n}^{\ast}-\alpha_{n})^{2}|0\rangle_{b_{n}}$
(30) $\displaystyle=$
$\displaystyle-\langle(a_{n}^{\dagger}-a_{n})^{2}+(\alpha_{n}^{\ast}-\alpha_{n})(a_{n}^{\dagger}-a_{n})\rangle-(\alpha_{n}^{\ast}-\alpha_{n})^{2}$
$\displaystyle=$
$\displaystyle-\langle|[(T_{n4}^{\ast}-T_{n1})b_{1}+(T_{n5}^{\ast}-T_{n2})b_{2}+(T_{n6}^{\ast}-T_{n3})b_{3}$
$\displaystyle-(T_{n4}-T_{n1}^{\ast})b_{1}^{\dagger}-(T_{n5}-T_{n2}^{\ast})b_{2}^{\dagger}-(T_{n6}-T_{n3}^{\ast})b_{3}^{\dagger}]^{2}\rangle$
$\displaystyle-(\alpha_{n}^{\ast}-\alpha_{n})^{2}$ $\displaystyle=$
$\displaystyle|T_{n4}^{\ast}-T_{n1}|^{2}+|T_{n5}^{\ast}-T_{n2}|^{2}+|T_{n6}^{\ast}-T_{n3}|^{2}$
$\displaystyle-(\alpha_{n}^{\ast}-\alpha_{n})^{2}.$

Then one obtains the variance of $\Delta p_{n}$

$\displaystyle(\Delta p_{n})^{2}$ $\displaystyle=\langle
p_{n}^{2}\rangle-\langle p_{n}\rangle^{2}$
$\displaystyle=|T_{n4}^{\ast}-T_{n1}|^{2}+|T_{n5}^{\ast}-T_{n2}|^{2}+|T_{n6}^{\ast}-T_{n3}|^{2}.$

## References

  * Hartmann  _et al._ [2006] M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, Nature Physics 2, 849 (2006).
  * Rossini and Fazio [2007] D. Rossini and R. Fazio, Phys. Rev. Lett. 99, 186401 (2007).
  * Carusotto  _et al._ [2009] I. Carusotto, D. Gerace, H. E. Tureci, S. De Liberato, C. Ciuti, and A. Imamoǧlu, Phys. Rev. Lett. 103, 033601 (2009).
  * Schiró  _et al._ [2012] M. Schiró, M. Bordyuh, B. Öztop, and H. E. Türeci, Phys. Rev. Lett. 109, 053601 (2012).
  * Vojta [2003] M. Vojta, Rep. Prog. Phys. 66, 2069 (2003).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, Cambridge, 2011).
  * Cardy [1996] J. Cardy, (1996).
  * Sondhi  _et al._ [1997] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar, Rev. Mod. Phys. 69, 315 (1997).
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Emary and Brandes [2003] C. Emary and T. Brandes, Phys. Rev. E 67, 066203 (2003).
  * Baumann  _et al._ [2011] K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger, Phys. Rev. Lett. 107, 140402 (2011).
  * Chen  _et al._ [2008] Q.-H. Chen, Y.-Y. Zhang, T. Liu, and K.-L. Wang, Phys. Rev. A 78, 051801 (2008).
  * Liu  _et al._ [2009] T. Liu, Y.-Y. Zhang, Q.-H. Chen, and K.-L. Wang, Phys. Rev. A 80, 023810 (2009).
  * Hwang  _et al._ [2015] M.-J. Hwang, R. Puebla, and M. B. Plenio, Phys. Rev. Lett. 115, 180404 (2015).
  * Ashhab [2013] S. Ashhab, Phys. Rev. A 87, 013826 (2013).
  * Liu  _et al._ [2017] M. Liu, S. Chesi, Z.-J. Ying, X. Chen, H.-G. Luo, and H.-Q. Lin, Phys. Rev. Lett. 119, 220601 (2017).
  * Chen  _et al._ [2020] X.-Y. Chen, Y.-Y. Zhang, L. Fu, and H. Zheng, Phys. Rev. A 101, 033827 (2020).
  * Chen  _et al._ [2021] X. Chen, Z. Wu, M. Jiang, X.-Y. Lü, X. Peng, and J. Du, Nat. Commun. 12, 1 (2021).
  * Cai and et al. [2021] M. L. Cai and et al., Nat. Commun. 12, 1126 (2021).
  * Dalibard  _et al._ [2011] J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Cai  _et al._ [2019] H. Cai, J. Liu, J. Wu, Y. He, S.-Y. Zhu, J.-X. Zhang, and D.-W. Wang, Phys. Rev. Lett. 122, 023601 (2019).
  * Li  _et al._ [2020] Y. Li, H. Cai, D.-w. Wang, L. Li, J. Yuan, and W. Li, Phys. Rev. Lett. 124, 140401 (2020).
  * Roushan  _et al._ [2017] P. Roushan, C. Neill, A. Megrant, Y. Chen, R. Babbush, R. Barends, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, _et al._ , Nat. Phys. 13, 146 (2017).
  * Hayward  _et al._ [2012] A. L. Hayward, A. M. Martin, and A. D. Greentree, Phys. Rev. Lett. 108, 223602 (2012).
  * Hayward and Martin [2016] A. L. Hayward and A. M. Martin, Phys. Rev. A 93, 023828 (2016).
  * Zhang  _et al._ [2021] Y.-Y. Zhang, Z.-X. Hu, L. Fu, H.-G. Luo, H. Pu, and X.-F. Zhang, Phys. Rev. Lett. 127, 063602 (2021).
  * Fallas Padilla  _et al._ [2022] D. Fallas Padilla, H. Pu, G.-J. Cheng, and Y.-Y. Zhang, Phys. Rev. Lett. 129, 183602 (2022).
  * Zhao and Hwang [2023] J. Zhao and M.-J. Hwang, Phys. Rev. Res. 5, L042016 (2023).
  * Zhao and Hwang [2022] J. Zhao and M.-J. Hwang, Phys. Rev. Lett. 128, 163601 (2022).