\section{Introduction}
The exploration of
quantum phenomena
is an outstanding challenge and has been one of the most active   arenas in condensed matter physics~\cite{sachdev2008quantum}.
Novel forms of phases
are enriched by
ongoing discoveries in transition metal compounds~\cite{Fe,MnCO,zhang2021strain,nature2021},
such as spin-orbit-entangled electronic phases~\cite{muhlbauer09skyrmion,yu2010real}. In particular, the quantum spin liquid (QSL) has recently emerged as a new paradigm in correlated electron physics,
as it holds promise for the potential application of quantum computing and quantum information.
One avenue towards QSL is the focus on highly frustrated materials, exemplified by either the geometrical frustration or the exchange frustration~\cite{Kiatev2014,kiatev2014b,HK1,HK2}.
The triangular, kagome, and pyrochlore structures are categorized into the first type,
while the actively sought-after Kitaev magnets belong to the second type. In this context,
much attention has been recently devoted towards $4d$ and $5d$ transition metal compounds,
due to the interplay between spin-orbit couplings and electronic correlations~\cite{material2018,material2020,Ru2019}.
The playground to search for the QSL was recently extended to $3d$ transition metal compounds~\cite{xu20possible}.
Despite these endeavors, the Kitaev QSL state has not been conclusively identified among most
candidate materials.
Certain long-range order (LRO) is unexpectedly found in a variety of
Mott insulators, indicating the existence of more conventional types of exchanges beyond a dominant bond-directional Kitaev interaction in these nonideal materials.

Since engineering robust QSL states remains challenging in spin-orbit-entangled candidate materials, an
alternative route for exploring novel phases of matter and forms of entanglement is the experimental implementation
with the help of other quantum systems in the laboratory. After Feynman's proposal in 1982~\cite{feynman1982computer}, the field of quantum simulation has been developing rapidly for decades and nowadays enables the investigation of quantum systems in a programmable fashion. Especially, recent
advances of the
laser technology and the laser
manipulation of atomic gases have made it feasible to implement a wide class of analog quantum simulations in atomic-molecular-optical (AMO) laboratories. Quantum simulators have been realized on
a few platforms,
e.g., ultracold atoms~\cite{trotzky2008time}, polar molecules~\cite{yanbo13nature},
trapped ions~\cite{roos14nature,monroe14nature},
photonic systems\,\cite{hung2016quantum,Douglas15quantum,Gonz15quantum},
and Rydberg atom arrays~\cite{ABrowaeys20np,lukin21nature256}, etc.
These systems can be finely tuned in a sufficiently precise way and observed in real time. The effective many-body Hamiltonian can be incorporated from recent developments of simulating quantum  magnetism and related quantum dynamics using atoms interacting with the same quantum modes, wherein the quantum channel can be the guided modes in the photonic crystal waveguides~\cite{hung2016quantum}, the photon of cavity modes~\cite{monika20heisenberg,chiocchetta2021cavity},
and the Rydberg dressing states~\cite{zoller15designing}.
With present architectures of quantum simulators, a generic  Hamiltonian consisting of flexible coupling graphs
can be freely realized, offering the opportunity to implement, simulate, and experimentally test fundamental paradigmatic model Hamiltonians.


This paper is organized as follows. In Sec.~\ref{sec:simulation}, we discuss the possible engineering of the
\emph{XY}-Gamma model using the AMO system facilitated by
coherent photon-mediated Raman transitions.
In our scenario, the independent control of  \emph{XX}-,\emph{ YY}-,\emph{ XY}- and \emph{YX}- terms can be achieved by a double $\Lambda$ scheme in neutral atoms.
Section~\ref{model} is devoted to exact solutions of the one-dimensional (1D) short-range interacting Hamiltonian using the Jordan-Wigner transformation. Two-site correlations and dimer correlations are then analyzed.
Next, the quantum steered coherence is investigated in Sec.~\ref{SQC}. A discussion and summary  follow in Sec.~\ref{summary}.


\section{Effective model with Photon-mediated atom-atom Interactions}\label{sec:simulation}
Recently, using photons to mediate controllable atom-atom interactions has been a well-established paradigm in AMO quantum simulation~\cite{zoller15designing,Mivehvar19toolbox,Mivehvar21review,chiocchetta2021cavity}.
The constituents of the engineered quantum spin models are not restricted to only Heisenberg-type
interactions~\cite{lev11frustration,Mivehvar17disorder,Esslinger18formation}.
The complex photon-mediated interaction graphs, including both their amplitudes and interaction ranges,
can be arbitrarily programed by state-of-the-art techniques involving more laser beams, posing envisioned synthetic quantum matter.
For instance, it is well known that the LRO is prohibited by the Mermin-Wagner theorem for short-range interacting models with a continuous symmetry in one spatial dimension~\cite{ren22longrange},
while the counterpart on a two-dimensional (2D) bipartite lattice is generally expected to host the long-range N\'{e}el order for any spin magnitude $S$,
although a rigorous proof of the existence of LRO in a 2D $S=1/2$ Heisenberg model is still lacking~\cite{you2009long,tang1989long,lin92long}.
The long-range \emph{XY} order is induced in an \emph{XXZ} chain by single-mode-cavity-mediated infinite-range interactions~\cite{li21long}. Instead, QSLs are stabilized in a 2D isotropic Heisenberg model by power-law decaying interactions in  multimode cavities~\cite{chiocchetta2021cavity}.


\begin{figure
  \centering
  \includegraphics[width=0.56\columnwidth]{Fig1.pdf}
  \caption{
Schematic setups. Atoms are trapped in 1D optical lattice with four-level energy diagram (inset).
Two pump lasers with Rabi frequencies $\Omega_1$ (blue) and $\Omega_2$ (green) induce $\sigma$-transitions between atomic ground and excited states $|g\rangle\leftrightarrow|1\rangle$ and
$|s\rangle\leftrightarrow|2\rangle$, respectively.
A
cavity mode (red)
induces $\pi$-transitions between $|s\rangle \leftrightarrow |1\rangle$ and  $|g\rangle\leftrightarrow |2\rangle$.}
  \label{fig:cavity_scheme}
\end{figure}

In the following,
we consider $N$ atoms that are trapped tightly in a 1D optical lattice,
as depicted in Fig.~\ref{fig:cavity_scheme}.
A double $\Lambda$ scheme of an atomic-level diagram is assumed,
where two
internal atomic states $\{|s\rangle, |g\rangle\}$ in the ground-state manifold
represent the pseudospin-$1/2$ states with energy $\{0, \omega_g\}$ and two auxiliary excited states $\{|1\rangle, |2\rangle\}$ exist with energy
$\{\omega_1, \omega_2\}$, respectively.
A photon mode at frequency $\omega_\mathbf{k}$ with the field operator $\hat{a}_\mathbf{k}$
induces $\pi$-transitions between atomic ground and excited states $|s\rangle\leftrightarrow|1\rangle$ and $|g\rangle\leftrightarrow|2\rangle$,
where $\mathbf{k}$ denotes the index of bosonic modes.
$G_\mathbf{k}(\bf{r})$ is the corresponding spatially dependent coupling strength.
Two pump lasers, denoted as $(L1, L2)$ with Rabi strengths $(\Omega_1, \Omega_2)$
and frequencies $(\omega_{L1}, \omega_{L2})$,
are implemented to induce $\sigma$ transitions between atomic ground and excited states, i.e., $|g\rangle\leftrightarrow|1\rangle$ and
$|s\rangle\leftrightarrow|2\rangle$, respectively.

Under rotating wave approximation, the atom-light hybrid system is described by the Hamiltonian $\hat{H} =\hat{H}_0+\hat{H}_{AL}$,
where $\hat{H}_0$ is the free Hamiltonian consisting of photon fields $\hat{a}_\mathbf{k}$ and atomic Zeeman energy levels $\hat{\sigma}_j^{aa}$,
and $\hat{H}_{AL}$ is the atom-light interaction Hamiltonian ($\hbar=1$ throughout):
\begin{eqnarray}
  \hat{H}_{0}&=&\sum_{\mathbf{k}} \omega_{\mathbf{k}}\hat{a}^\dagger_\mathbf{k}\hat{a}_\mathbf{k}
  +\sum_{j}\sum_{\alpha}\omega_\alpha\hat{\sigma}^{\alpha\alpha}_j,\\
  \hat{H}_{\rm AL}&=&\sum_j \left[{\Omega_1(\mathbf{r}_j)}e^{-i\omega_{L1}t}\hat{\sigma}_j^{1g}
  +{\Omega_2(\mathbf{r}_j)}e^{-i\omega_{L2}t}\hat{\sigma}_j^{2s}+{\rm H.c.}\right]\nonumber\\
  &&+\sum_j\sum_{\mathbf{k}}\left[G_{\mathbf{k}}(\mathbf{r}_j)\hat{a}_{\mathbf{k}}(\hat{\sigma}_j^{1s}+\hat{\sigma}_j^{2g})+{\rm H.c.}\right].
\end{eqnarray}
Here the atomic transition operators are defined as $\hat{\sigma}_j^{ab}\equiv|a\rangle_j\langle b|$ with four atomic energy levels  $a,b=\{g,s,1,2\}$ and $j$ labels the site.
Two-atom interaction is synthesized by two Raman transitions, where the photon field $\hat{a}_\mathbf{k}$
provides the two-body correlation between the two atoms, as depicted in Fig.\,\ref{fig:cavity_scheme}.



Working in a rotating frame defined by
$\hat{U}$=$\exp\{i[\sum_{j}(\omega_{-}\hat{\sigma}_j^{gg}
    +\omega_{+}\hat{\sigma}^{11}_j+\omega_{L2}\hat{\sigma}^{22}_j)
   +\sum_\mathbf{k}\omega_{+}\hat{a}_\mathbf{k}^\dagger\hat{a}_\mathbf{k}]t\}$ with $\omega_{\pm}$=$(\omega_{L2}\pm\omega_{L1})/2$, the transformed Hamiltonian
$\tilde{H}=U \hat{H} U^\dagger+i(\partial_t U) U^\dagger$ reads
\begin{eqnarray}
\tilde{H}&=&-\sum_{\mathbf{k}} \Delta_\mathbf{k} \hat{a}^\dagger_\mathbf{k}\hat{a}_\mathbf{k}
+\sum_j \left[\delta \hat{\sigma}^{gg}_j  -\Delta_1\hat{\sigma}^{11}_j-\Delta_2\hat{\sigma}^{22}_j\right]\nonumber\\
&&+\sum_j\left[{\Omega_1(\mathbf{r}_j)}\hat{\sigma}_j^{1g}
+{\Omega_2(\mathbf{r}_j)}\hat{\sigma}_j^{2s}+{\rm H.c.}\right]\nonumber\\
&&+\sum_j\sum_{\mathbf{k}}
\left[G_{\mathbf{k}}(\mathbf{r}_j)\hat{a}_{\mathbf{k}}(\hat{\sigma}_j^{1s}+\hat{\sigma}_j^{2g})+{\rm H.c.}\right]\label{eq:Hrot},
\end{eqnarray}
where
$\Delta_{\mathbf{k}}\equiv(\omega_{L1}+\omega_{L2})/2-\omega_{\mathbf{k}}$,
$\Delta_1\equiv (\omega_{L1}+\omega_{L2})/2-\omega_1$,
$\Delta_2\equiv \omega_{L2}-\omega_2$ and
$\delta=\omega_g-(\omega_{L2}-\omega_{L1})/2$.
Supposing the frequencies of pump lasers and bosonic modes are all far detuned from the atomic transitions, i.e., $\Delta_1$ and $\Delta_2$ are much larger
than the Rabi coupling coefficients
$|\Omega_1|, |\Omega_2|$, and $|G_\mathbf{k}|$,
we can safely eliminate the excited states $|1\rangle$ and $|2\rangle$ to obtain the effective
Hamiltonian in the ground-state manifold (see Appendix \ref{DerivationOfEffectiveHamiltonian}),
\begin{eqnarray}
  \hat{H}_{gs}&=&\sum_j\left[(\delta+V_g(\mathbf{r}_j))\hat{\sigma}_j^{gg}+
  V_s(\mathbf{r}_j)\hat{\sigma}_j^{ss}\right]\nonumber\\
&&  +\sum_j\left[\hat{\sigma}_j^{gs}\hat{\Omega}(\mathbf{r}_j)
  +\hat{\Omega}^\dagger(\mathbf{r}_j)\hat{\sigma}_j^{sg}\right],\label{eq:Hgs_a}
\end{eqnarray}
where $\hat{\Omega}(\mathbf{r})= \Omega_1^\ast \sum_\mathbf{k}G_\mathbf{k}(\mathbf{r})\hat{a}_{\mathbf{k}}/\Delta_1
+ \sum_\mathbf{k}G^\ast_\mathbf{k}(\mathbf{r})\hat{a}^\dagger_{\mathbf{k}}/\Omega_2 $  and the light shifts
for ground states $V_g(\mathbf{r})={|\Omega_1|^2}/{\Delta_1}+{\sum_{kk^\prime}G^\ast_\mathbf{k}G_{\mathbf{k}^\prime}\hat{a}^\dagger_{\mathbf{k}}\hat{a}_{\mathbf{k}^\prime}}/{\Delta_2}$,
$V_s(\mathbf{r})={|\Omega_2|^2}/{\Delta_2}+{\sum_{kk^\prime}G^\ast_\mathbf{k}G_{\mathbf{k}^\prime}\hat{a}^\dagger_{\mathbf{k}}\hat{a}_{\mathbf{k}^\prime}}/{\Delta_2}$.
In the adiabatic limit due to large detuing or large dissipation $\kappa$
of modes $\hat{a}_\mathbf{k}$,
the photon field can be approximated by its steady-state value,
\begin{eqnarray}\hat{a}_{\mathbf{k}}^{\rm ss}=
\sum_j\, G^\ast_{\mathbf{k}}(\mathbf{r}_j)
\left(\frac{\Omega_1}{\Delta_1}\hat{\sigma}^-_j+
\frac{\Omega_2}{\Delta_2}\hat{\sigma}^+_j\right)/{\tilde{\Delta}_\mathbf{k}}, \end{eqnarray}
with $\tilde{\Delta}_\mathbf{k}=\Delta_\mathbf{k}+i\kappa-
\sum_j \,|G_\mathbf{k}(\mathbf{r}_j)|^2\left(\hat{\sigma}_j^{ss}/{\Delta_1}
+\hat{\sigma}_j^{gg}/{\Delta_2}\right)$.
Note that here we have included the
cavity dissipation $\kappa$ phenomenologically.
By adiabatically eliminating the photon degree of freedom,
in terms of Pauli operators,  $\sigma^x_j=\hat{\sigma}_j^{sg}+\hat{\sigma}_j^{gs}, \sigma_j^{y}=i(\hat{\sigma}^{sg}_j-\hat{\sigma}^{gs}_j), \sigma^z_j=\hat{\sigma}_j^{gg}-\hat{\sigma}_j^{ss}$,
we obtain the effective spin
Hamiltonian in a compact form,
\begin{eqnarray}
  \hat{H}_{\rm int}\!
  &=&\sum_{ij}\!(J_{ij}^x \sigma_i^x \sigma_j^x
  +J_{ij}^y\sigma_i^y \sigma_j^y)\!+\!J_{ij}^{\rm DM}(\sigma_i^x \sigma_j^y
  -\sigma_i^y\sigma_j^x)\nonumber\\
  &+&\sum_{ij} J_{ij}^{\rm SO}\left(\sigma_i^y \sigma_j^x
  + \sigma_i^x \sigma_j^y\right) + \sum_{j} h_j^z {\sigma}_j^z.\label{eq:Hamiltonian0}
\end{eqnarray}
Here $h_j^z=\delta/2+ V_g(\mathbf{r}_j)-V_s(\mathbf{r}_j)$, $J_{ij}^x=2$({Re}${[\Lambda_0]}$+{Re}${[\Lambda_1]}$),
$J_{ij}^y=$2({Re}${[\Lambda_0]}$$-${Re}${[\Lambda_1]}$),
$J_{ij}^{\rm DM}=$2{Im}$[\Lambda_0]$, and $J_{ij}^{\rm SO}=$-2{Im}$[\Lambda_1]$, where {Re} ({Im}) indicates the real (imaginary) part of a complex variable,  $\Lambda_{0,1}\equiv\Lambda_{0,1}(\mathbf{r}_i,\mathbf{r}_j)$ with
$\Lambda_0(\mathbf{r},\mathbf{r}^\prime)=\sum_\mathbf{k}
[\frac{\Omega_1^\ast(\mathbf{r})\Omega_1(\mathbf{r}^\prime)}{\Delta_1^2\tilde{\Delta}_\mathbf{k}}G_\mathbf{k}(\mathbf{r})G^\ast_{\mathbf{k}}(\mathbf{r}^\prime)
+\frac{\Omega_2(\mathbf{r})\Omega_2^\ast(\mathbf{r}^\prime)}{\Delta^2_2\tilde{\Delta}^\ast_\mathbf{k}}G_\mathbf{k}^\ast(\mathbf{r})G_{\mathbf{k}}(\mathbf{r}^\prime)]$
and $\Lambda_1(\mathbf{r},\mathbf{r}^\prime)=\sum_\mathbf{k}[
\frac{\Omega_1^\ast(\mathbf{r})\Omega_2(\mathbf{r}^\prime)}{\Delta_1\Delta_2\tilde{\Delta}_\mathbf{k}}G_\mathbf{k}(\mathbf{r})G^\ast_{\mathbf{k}}(\mathbf{r}^\prime)
+\frac{\Omega_2(\mathbf{r})\Omega_1^\ast(\mathbf{r}^\prime)}{\Delta_1\Delta_2\tilde{\Delta}^\ast_\mathbf{k}}G_\mathbf{k}^\ast(\mathbf{r})G_{\mathbf{k}}(\mathbf{r}^\prime)]$.
The detailed definitions of coefficients in Eq.\,(\ref{eq:Hamiltonian0}) are derived in Appendix\,\ref{DerivationOfEffectiveHamiltonian}.

The first bracketed term of Eq.~\eqref{eq:Hamiltonian0} corresponds to the conventional \emph{XY}-type interactions, and the term in the second brackets denotes the $z$ component of
Dzyaloshinskii-Moriya interactions (DMIs).
The cross-couplings in the third brackets between $x$ and $y$  spin components are referred to as symmetric off-diagonal $\Gamma$ interactions.
Exotic forms like DMIs
originated from spin-orbit couplings~\cite{dzyaloshinsky1958thermodynamic,DM1959,XYDM1960,DM1960b} and were firstly devised to account for the
weak ferromagnetism in antiferromagnetic crystals~\cite{Fe,MnCO,zhang2021strain,nature2021},
favoring chiral states such as spin spirals and skyrmions~\cite{muhlbauer09skyrmion,yu2010real}. Concurrently,
the importance of pervasive off-diagonal $\Gamma$ interactions can be traced back to the study of the
Kitaev-Heisenberg model~\cite{Kiatev2014,kiatev2014b,HK1,HK2}.
Further research suggests that the symmetric off-diagonal $\Gamma$ interactions should also be taken into account to explain the possible QSLs observed in experiments~\cite{Luo2021unveiling,Luo2021unusual,you2020Lifshitz,Alberto2020phase,Erik2021entanglement,Yang2020comprehensive}





\section{Exact solution and correlations}\label{model}
Note that the cavity-mediated spin-spin interactions in Eq.~(\ref{eq:Hamiltonian0}) have infinite range
if only a single cavity mode is involved. The finite-range interactions are achieved by using a multimode
cavity. The photon modes can be
the guide modes in photonic crystal waveguides with quasimomentum $\mathbf{k}$\,\cite{,Douglas15quantum,Gonz15quantum,hung2016quantum},
or the near-degenerate transverse cavity  modes\,\cite{Lev18tunable,Guo21nature}.
In particular, the multifrequency driving also enables finite-range interactions between intracavity atoms
in a single-mode cavity, which has been recently realized experimentally\,\cite{periwal2021programmable}.
Therefore,
an effective Hamiltonian with tunable interaction strength and interaction range can be constructed by using multimodes $\{\hat{a}_\mathbf{k}\}$.
In this case,
the interference of cavity modes may render the beyond-nearest-neighbor
couplings to be negligibly small.

In the following, we concentrate on the spin models
for an ensemble of spin-1/2 interacting particles
on a 1D lattice with nearest‐neighbor interactions only.   The spin Hamiltonian can be rewritten as
\begin{eqnarray}
\hat{H}&=&\sum_{j=1}^{N}J\left(\frac{1+\gamma}{2}\sigma_j^x\sigma_{j+1}^x+\frac{1-\gamma}{2}\sigma_j^y\sigma_{j+1}^y\right)\nonumber\\
&&+\sum_{j=1}^{N} \left[\Gamma\left(\sigma_j^x\sigma_{j+1}^y+\alpha\sigma_j^y\sigma_{j+1}^x\right)+h\sigma_j^z\right],
\label{eq:Hamiltonian1}
\end{eqnarray}
where the antiferromagnetic coupling  $J\equiv J_{j,j+1}^x+J_{j,j+1}^y$
between the nearest-neighbor atoms is set up as an energy unit for simplicity unless otherwise stated, i.e., $J=1$, $\gamma\equiv(J_{j,j+1}^x-J_{j,j+1}^y)/J$
serves as the anisotropy parameter, $\Gamma \equiv J_{j,j+1}^{\rm DM}+J_{j,j+1}^{\rm SO}$ characterizes the amplitude of off-diagonal exchange interactions,
$\alpha\equiv (J_{j,j+1}^{\rm SO}-J_{j,j+1}^{\rm DM})/ \Gamma$ denotes the relative coefficient of off-diagonal exchange couplings, and $h$ represents the strength of the uniform transverse field.
The $\Gamma$ term reduces to the DMI for $\alpha=-1$ and the symmetric off-diagonal exchange interaction for $\alpha=1$.
In what follows, we impose periodic boundary conditions (PBCs) with $\sigma_{N+1}^a\equiv\sigma_{1}^a$ ($a=x,y,z$).


The motivations of exploring the quantum criticality in the \emph{XY}-Gamma model are two fold. On one hand, the low-dimensional quantum magnets have been particularly of concern owing to their evident quantum aspects and substantial corrections to classical counterparts. The quantum criticalities have been explored in a few magnetic materials. The notable examples range from the spin-1/2 Ising ferromagnet LiHoF$_4$~\cite{bitko1996quantum},
SrCo$_2$V$_2$O$_8$~\cite{cui2019quantum},
Cs$_2$CoCI$_4$~\cite{kenzelmann2002order}, and CoNb$_2$O$_6$~\cite{coldea2010quantum} to  BaCo$_2$V$_2$O$_8$~\cite{PhysRevLett.123.027204}
as well as the spin-1 ferromagnetic Heisenberg chain NiNb$_2$O$_6$~\cite{chauhan20tunable}. To date, quantum phase transitions (QPTs) of analog models have been studied in different contexts, such as the \emph{XY} model with DMIs~\cite{XYDMYTC}.
The simultaneous appearance of off-diagonal exchange $\Gamma$ interactions
and \emph{XY}-type interaction in the presence of external fields, especially counteracting the disordered state,
is less systematically clear.
One the other hand, the merit of
Eq.~(\ref{eq:Hamiltonian0}) resides in its exact solvability. The analytical results render the possibility to calculate accurately the experimental measurable quantities, in particular various dynamic ones, and thus serve as a benchmark for more sophisticate models. As we shall demonstrate, the extracted critical exponents
can be relevant to the experimental measurement of thermodynamic quantities and information measures.



\begin{figure}[ht]
\centering
\includegraphics[width=0.9\columnwidth]{Fig2.pdf}
\caption{The contour map of excitation gap with respect to $\alpha$ and $h$.
The white dashed lines correspond to the critical lines separating
 phases I (AFM), II (PM), and III (Spiral) of the \emph{XY}-Gamma model with $J=1.00$, $\Gamma=0.60, \gamma=0.60$.}
\label{phasediagram}
\end{figure}

The Hamiltonian in Eq.~(\ref{eq:Hamiltonian1}) is analytically solved by means of Jordan-Wigner, Fourier, and Bogoliubov transformations.
The detailed diagonalization procedure is shown in  Appendix \ref{energyspect1}.
Ultimately, the Hamiltonian can be brought into a diagonal form of a 
less fermion in the momentum space,
\begin{eqnarray}
\hat{H}&=& \sum_{k}{\varepsilon}_{k}({b}_{k}^{+}{b}_{k}-\frac{1}{2})
\end{eqnarray}
where the energy spectrum of fermionic quasiparticles is
\begin{eqnarray}
{\varepsilon}_{k}&=&
2\sqrt{[\Gamma^2(\alpha+1)^2+\gamma^2]\sin^2{k}+(\cos{k}-h)^2} \nonumber\\ &&-2\Gamma(1-\alpha)\sin{k} 
\label{varepsilon}
\end{eqnarray}
With the excitation energy at hand, the energy gap $\Delta$=$\min_k {\varepsilon}_{k}$ can be
determined. As shown in Fig. \ref{phasediagram}, the ground-state phase diagram consists of three phases by varying $\alpha$ and $h$.  The horizontal segment $h_{c,1}=1$ for $\alpha > - {\gamma^2}/{4\Gamma^2}$ separates the gapped phase I and phase II, while the gapless phase III is encompassed by the separatrices
 $\alpha_{c,1}=- \gamma^2/(4\Gamma^2)$ for $h\le 1$
and $\alpha_{c,2}=(1-h^2-\gamma^2)/(4\Gamma^2)$ for $h>1$, or equivalently, $h_{c,2}=\sqrt{1-\gamma^2-4\Gamma^2\alpha}$  for $\alpha<- \gamma^2/(4\Gamma^2)$.

Typically, one can define spatial and temporal characteristic lengths that have diverging behavior as the control parameter $\lambda$ approaches the threshold value $\lambda_c$. This diverging property of characteristic lengths with the correlation-length critical exponent $\nu$  and the dynamical critical exponent $z$ enables one to define universality classes.   The critical behavior
is determined by those low-energy states near
the critical mode.
The dynamical exponent $z$ relates the scaling of energy
to length scales,
which can be retrieved by the shape of the spectra near the gap closing mode $\Delta \sim (k-k_c)^{z}$. As $\lambda$ approaches $\lambda_{c}$, the gap vanishes as $\Delta\sim(\lambda-\lambda_{c})^{\nu z}$.  To this end, we expand
Eq.~(\ref{varepsilon}) at the critical line $h_{c,1}$ around the gap closing momentum $k_c=0$,
\begin{eqnarray}
   {\varepsilon}_{k} \simeq  2\sqrt{\gamma^2+\Gamma^2(1+\alpha)^2} \vert k\vert-2\Gamma(1-\alpha) k.
\end{eqnarray}
The relativistic spectra at the critical line $h_{c,1}$ imply that $z$ = 1.
The gap near $h_{c,1}$ is approximated as
\begin{equation}
\label{Delta}
\Delta \simeq
\,2\vert h-h_{c,1}\vert,
\end{equation}
and one then finds $\nu z=1$. In this case, the quantum critical point (QCP) between phase I and phase II belongs to 2D Ising universality class characterized by $\nu=1$, $z=1$. On the verge between gapped phase II and  gapless phase III, one can find the spectra vanish at an incommensurate momentum $k_c=\arccos(h_{c,2}^{-1})$,
 \begin{eqnarray}
\varepsilon_k \simeq\left[\frac{\Gamma(1-\alpha)\cos^2 k_c}{2\sin k_c}+\frac{(h_{c,2}^{-1}-h_{c,2})^2}{2\Gamma(1-\alpha)\sin k_c}\right](k-k_c)^2.\quad \quad\label{quadratic}
\end{eqnarray}
The above quadratic dispersion indicates the dynamical exponent $z = 2$.
While expanding the gap around the QCP from the upper threshold
one obtains the excitation as
\begin{equation}
\label{Delta2}
\Delta \simeq \frac{2(h_{c,2}-h_{c,2}^{-1})}{\Gamma(1-\alpha)\sin k_c} (h-h_{c,2}).
\end{equation}
The critical exponents $\nu=1/2$ and $z=2$ annotate that the QPT is in the so-called Lifshitz universality class~\cite{Rufo2019Multicritical}, which corresponds to the universality class of quantum criticality of free fermions.  %
In the case of  the I-III transition $\alpha_{c,1}$, the spectra are found to be quadratic in $k$ around the gap closing mode $k_c=\arccos h$ as
  \begin{eqnarray}
\varepsilon_k \simeq \frac{\sqrt{1-h^2}}{\Gamma(1-\alpha)} (k-k_c)^2, \label{quadratic2}
\end{eqnarray}
which yields $z$=2. Similarly, the gap around the critical point $\alpha_{c,1}$ from above obeys a power-law relation as
\begin{equation}
\label{Delta3}
\Delta \simeq \frac{4 \Gamma\sin k_c}{1-\alpha_{c,1}} (\alpha-\alpha_{c,1}).
\end{equation}
The scaling form
in Eq.~(\ref{Delta3}) reveals that the transition also belongs to the Lifshitz universality class with $z$ = 2 and $\nu$ = 1/2.

We also calculate the second derivative of the ground-state
energy density $e_0=-\sum_k \vert {\varepsilon}_{k} \vert/(2N)$ in
Fig.~\ref{energyderivative}, which showcases extreme values around critical points.
With increase of the system sizes, the peaks of $-\partial^2e_0/\partial h^2$ for $\alpha=0.50$ become
more pronounced. To be concrete, a logarithmic singularity across the QPT between phase I and phase II is identified as
\begin{eqnarray}
\left(-\frac{\partial^2{e_0}}{\partial{h^2}}\right)_{\rm max}&=&a_{\rm E} \ln N+c_1.
\end{eqnarray}
Meanwhile, in the vicinity of the critical point in the thermodynamic limit, one finds
 \begin{eqnarray}
\left(-\frac{\partial^2{e_0}}{\partial{h^2}}\right)&=& b_{\rm E}\ln |h-h_c|+c_2.
\end{eqnarray}
The numerical fittings in Fig.~\ref{energyderivative}(a) yield
$a_{\rm E}=0.2871\pm0.0058$, $c_1=0.1878$, $b_{\rm E}=-0.2887\pm0.0009$, and $c_2=0.1162$.
According to the logarithmic scaling ansatz~\cite{scalling2006},
the ratio $|a_E/b_E|$ equals the correlation-length exponent $\nu$ $\simeq 1$,
confirming that the QPT from phase I to  phase II coincides with a second-order transition.
The retrieved specific heat exponent $\alpha$ = $2-(d + z)\nu$=0 validates the scaling relation for the logarithmic scaling in $d=1$ dimension.
In contrast, one observes in Fig.~\ref{energyderivative}(b)
that $-\partial^2e_0/\partial \alpha^2$ exhibits a size-independent discontinuity at the critical points for $h=0.50$ and $h=1.17$,
which is a common feature of the transition between the gapless phase and the gapped phase~\cite{GammaLZA}.

\begin{figure}[ht]
  \includegraphics[width=0.85\columnwidth]{Fig3a.pdf}
  \includegraphics[width=0.85\columnwidth]{Fig3b.pdf}
  \caption{(a)  ${-\partial^2e/\partial h^2}$ as a function of $h$ with $\alpha=0.50$  for various system sizes.
  The left inset shows the scaling behavior between the maximum value of ${-\partial^2e/\partial h^2}$
  and the system size $N$. The right inset shows $-{\partial^2e/\partial h^2}$ in the vicinity of the critical point.  The symbols denote the numerical results, and the solid lines correspond to the linear fittings.
  (b) ${-\partial^2e/\partial \alpha^2}$ vs $\alpha$ with $h=0.50$ and $h=1.17$.
  Other parameters: $J=1.00$, $\gamma=0.60$, $\Gamma=0.60$.}
  \label{energyderivative}
\end{figure}

The nature of the
ground state can be gained from the two-qubit correlation functions \begin{eqnarray}
G_{i,j}^{a,b}=\langle\sigma_i^a\sigma_j^b\rangle-\langle\sigma_i^a \rangle \langle\sigma_j^b\rangle
\end{eqnarray}
with $a$, $b$=$x$, $y$, $z$.
In fact, $G_{i,j}^{a,b}$ can be abbreviated as $G_{r}^{a,b}$ with $r=i-j$ due to the translational invariance of
Eq.~(\ref{eq:Hamiltonian1}). A generic correlation ${\langle { \sigma}_{i}^{a}{ \sigma}_{j}^{b}\rangle}$
can be expressed as a Pfaffian form in terms of Wick's theorem~\cite{Pfaffian}, which is the determinant of the $2r\times 2r$ dimensional antisymmetric matrix.
The detailed calculation is exhibited in Appendix \ref{correlation function1}. One observes that the nearest-neighbor correlation functions display kinks across QCPs between the gapless phase and gapped phases. With increasing $\alpha$ for $h=0.50$ in Fig.~\ref{G1}(a), the dominant nearest-neighbor correlation changes from a positive value of $G_1^{yx}$ to a negative value of $G_{1}^{xx}$, implying a QPT from the gapless spiral phase to the gapped antiferromagnetic (AFM) phase.  Instead, the ruling correlations $G_{1}^{zz}$ in Fig.~\ref{G1}(b) for $h=1.17$ suggest that phase II belongs to the paramagnetic (PM) phase.  Similar trends of correlations for $\alpha=0.50$ are displayed in Fig.~\ref{Gxx}. The dominating nearest-neighbor correlations change from a negative value of $G_{1}^{xx}$ to a positive value of $G_{1}^{zz}$ across $h_{c,1}$. Therefore, the first-order derivative of $G_{1}^{xx}$ presents a pronounced peak at $h_{c,1}$ in Fig.~\ref{Gxx}(b).  We can further find the first-order derivative of $G_{1}^{xx}$ also follows a logarithmic divergence across the second-order QPT as
 \begin{eqnarray}
&&\left(\frac{\partial G^{xx}}{\partial{h}}\right)_{\rm max}= a_{\rm G} \ln{N}+c_3, \\
&&\left(\frac{\partial G^{xx}}{\partial{h}}\right)= b_{\rm G}\ln{|h-h_c|}+c_4,
\end{eqnarray}
where $a_{\rm G}=0.2789\pm0.0045$, $b_{\rm G}=-0.2811\pm0.0001$, $c_3=0.1758, $ and $c_4=0.0265$.
 In this case, one can speculate that $\nu\approx|a_{\rm G}/b_{\rm G}|$=$0.9922\pm0.0209$, which is quite close to the value retrieved from the second derivative of the ground-state energy density.


\begin{figure}[ht]
\includegraphics[width=0.85\columnwidth]{Fig4a.pdf}
 \includegraphics[width=0.85\columnwidth]{Fig4b.pdf}
\caption{ The nearest-neighbor correlation functions
with respect to $\alpha$ for (a) $h=0.50$ and  (b) $h=1.17$.  Other parameters: $N$=2000, $J=1.00$, $\Gamma=0.60$, $\gamma=0.60$.}
\label{G1}
\end{figure}


\begin{figure}[ht]
\includegraphics[width=0.9\columnwidth]{Fig5a.pdf}
 \includegraphics[width=0.9\columnwidth]{Fig5b.pdf}
\caption{(a) The correlation functions $G_{1}^{xx}$, $G_{1}^{zz}$, $G_{1}^{xy}$, $G_{1}^{yx}$ as a function of $h$ for $\alpha=0.50$. The inset shows $\vert G_r^{zz}\vert$ vs $r$ for $h=0.50$ and $1.50$ with $N=2000$.
 (b) The first-order derivative of $G_{1}^{xx}$ with $N=200$, $400$, $600$, $800$, $1000$. The left inset shows the scaling behavior between the maximum value of ${\partial G_{1}^{xx}/\partial h}$ and the system size $N$. The right inset shows $\partial G_{1}^{xx}/\partial h$ in the vicinity of the critical point. The symbols denote the numerical results, and the solid lines correspond to the linear fittings. Other parameters: $J=1.00$, $\gamma=0.60$, $\Gamma=0.60$.
}
\label{Gxx}
\end{figure}
\begin{figure}[ht]
  \centering
  \includegraphics[width=0.92\columnwidth]{Fig6a.pdf}
  \includegraphics[width=0.92\columnwidth]{Fig6b.pdf}
  \caption{(a) The correlation function $|G_{i,i+r}^{xy}|-|G_{i,i+r}^{yx}|$ with respect to $h$ for different $r$. The inset shows $|G_{i,i+r}^{xy}|-|G_{i,i+r}^{yx}|$ vs $r$ for $h=0.50$. (b) The dimer correlation $\langle \kappa_{i} \kappa_{i+r}\rangle-\langle \kappa_{i} \rangle \langle \kappa_{i+r}\rangle $ with respect to $h$ for different $r$. The inset displays the dimer correlation vs $r$ for $h=0.50$. Other parameters: $N$=2000, $J=1.00$, $\alpha=-0.50$,$\gamma=0.60$, $\Gamma=0.60$.}
\label{Gr}
\end{figure}

It is well known that the AFM phase hosts N\'{e}el LRO, while the LRO is absent in the PM phase, as is unraveled in the inset of Fig. \ref{Gxx}(a). One can further notice that the amplitudes of $G_{1}^{xy}$ and $G_{1}^{yx}$ coincide in the gapped phases, while become unbalanced in the gapless phase. Hence, such a feature suggests that  $\vert G_{1}^{xy}\vert-\vert G_{1}^{yx}\vert$ is a well-defined order parameter to identify the spiral phase for the \emph{XY}-Gamma model. A natural question arises whether there is LRO in the gapless spiral phase. To probe this question, we numerically calculate the vector-chiral correlations $|G_{r}^{xy}|-|G_{r}^{yx}|$  for different distances $r$ in Fig.~\ref{Gr}(a), in which the absolute value is taken in order to remove the indeterminate sign of the numerical Pfaffian calculation. One can find that
$|G_{r}^{xy}|-|G_{r}^{yx}|$ is always zero in the gapped phase as a consequence of $G_r^{xy}=G_r^{yx}$,  while it remains finite in the gapless spiral phase. Upon increasing the distance $r$, the correlation presents an oscillating decline as $r^{-1/2}$ shown in Fig.~\ref{Gr}(a)~\cite{XYDMYTC},
suggesting the existence of a quasi-long-range order of an incommensurate spiral order. Next, to delve more deeply into the spiral order, we consider four-qubit correlations as exemplified
by the dimer correlation
\begin{eqnarray}
D_{j,j+r}=\langle\kappa_j \kappa_{j+r} \rangle- \langle\kappa_j\rangle \langle \kappa_{j+r} \rangle,
\end{eqnarray}
where the $z$-component vector chiral order parameter is defined by~\cite{mcculloch08vector,ueda2014vector}
\begin{eqnarray}
\kappa_j=\left(\vec{\sigma}_j \times \vec{\sigma}_{j+1}\right)^z.
\end{eqnarray}
We thus calculate the dimer correlation $D_{j,j+r}$ as a function of $h$ for $\alpha=0.50$ in Fig.~\ref{Gr}(b).
One observes that the dimer correlations oscillate in the spiral phase and tend to decay with increasing the distance $r$ between the dimers.
We further find that dimer correlations persist only for a few sites, and thus the four-qubit correlations decay more rapidly than the two-qubit counterparts.




\section{STEERED QUANTUM COHERENCE}\label{SQC}
\begin{figure}[h]
  \centering
  \includegraphics[width=0.9\columnwidth]{Fig7a.pdf}\\
  \includegraphics[width=0.9\columnwidth]{Fig7b.pdf}
  \caption{(a) The SQC $C_{re}^{st}$ vs the transverse field $h$. The dash-dotted, dotted, solid and dashed lines correspond to $\alpha=-0.80$, $-0.50$, $0.50$, and $0.80$, respectively.(b) $C_{re}^{st}$ vs $h$ for $\alpha=-0.50$. The dotted, solid  and dashed lines correspond to different distances between the two qubits $r=1$, $2$, and $3$. Inset shows the first derivative of $C_{re}^{st}$ with respect to $h$.  Other parameters: $N$=2000, $J=1.00$, $\gamma=0.60, \Gamma=0.60$.}
  \label{SQC1}
\end{figure}

\begin{figure}[h]
\includegraphics[width=0.9\columnwidth]{Fig8a.pdf}
 \includegraphics[width=0.9\columnwidth]{Fig8b.pdf}
\caption{(a)$C_{re}^{st}$ vs $h$ for $\alpha=0.50$.  The left inset shows the scaling behavior between the maximum value of ${\partial C_{re}^{st}/\partial h}$ and the system size $N$, and the right inset shows $\partial C_{re}^{st}/\partial h$ in the vicinity of the critical point. (b) The first derivative of $C_{re}^{st}$ with respect to $h$ for $\alpha=0.50$ with $r=1$, $2$, and $r=3$. Other parameters: $N$=2000, $J=1.00$, $\gamma=0.60, \Gamma=0.60$.}
\label{SQCr}
\end{figure}

In recent years, a few approaches inherited from quantum information have been employed to characterize the QPTs, such as quantum entanglement~\cite{vidal2003entanglement,osterloh2002scaling,gu2003entanglement},
quantum discord~\cite{ollivier2001quantum,modi2012classical},
quantum coherence~\cite{chen2016coherence}, and
fidelity susceptibility~\cite{you2007fidelity}.
These information measures have played the role of universal order parameters
and deepened understanding of the
quantumness of correlations at quantum criticality.
Quantum coherence, a landmark manifestation of quantum superposition, has been widely
recognized as a common necessary condition for both entanglement and other types of quantum correlations.
Given that, the precise control over the photon-mediated  interactions between atoms means that  any pair of constituent particles could be ideally coaxed into any desired quantum-mechanical superposition state. In particular, the neutral atoms couple weakly to
the environment, allowing relatively long coherence times. Hence, such an array of atoms can function as a versatile model
for the study of quantum coherence, which may ultimately help
alleviate
the adverse effects of decoherence in quantum computation and quantum information processing~\cite{Bloch2008quantum}.
Based on a rigorous framework to quantify coherence~\cite{Baumgratz14quantify},
a few measures have been proposed, including the relative entropy,  $l_1$ norm of coherence,
Wigner-Yanase-Dyson skew information, and Jensen-Shannon divergence.  These coherence-based indicators have been applied to identify QPTs in many-body systems~\cite{you17emergent,you18quantum,yi19criticality}.
However, their feasibility strongly depends on a careful choice of the specific basis in advance, so these conventional quantum coherence measurements may extract useless information if the reference basis is inappropriate.
To overcome the irrationality, the steered quantum coherence (SQC) was proposed recently~\cite{M.L.Hu2018,M.L.Hu2020,M.L.Hu2021}. The SQC is based on the mutually unbiased bases and shows its figure of merit in characterizing quantum criticality.

For a bipartite state $\rho_{AB}$ shared by  Alice and Bob, the SQC was defined by Alice's local measurements and classical communication between Alice and Bob.
To be explicit, Alice carries out one of some preagreed measurements ${\sigma^\mu}$ $({\mu=x,y,z})$ on qubit $A$ and informs Bob of  the chosen observable $\sigma^\mu$
and the outcome $a\in\{0,1\}$.  Bob's system then collapses to the ensemble states \{${p_{\mu,a},\rho_{B|\Pi_\mu^a}}$\}
with $p_{\mu,a}={\rm tr}(\Pi_\mu^a\rho_{AB})$ being the probability of Alice's outcome $a$,
and $\rho_{B|\Pi_\mu^a}={\rm tr}_A(\Pi_\mu^a\rho_{AB})/p_{\mu,a}$ being Bob's conditional state.
Bob can measure the coherence of the ensemble \{${p_{\mu,a},\rho_{B|\Pi_\mu^a}}$\} with respect to
the eigenbasis of either one of the remaining two Pauli operators $\sigma^\nu$ ($ \nu\neq\mu$).
After  all Alice's  possible measurements
\{${\Pi_\mu^a}$\}$_{\mu=x,y,z}$ with equal probability,
the SQC of qubit $B$ can be defined as the following averaged quantum coherence:
\begin{eqnarray}
  C_{re}^{st}(\rho_{AB})\!&=&\!\frac{1}{2}\Sigma_{{\mu\neq\nu,a}}  p_{\mu,a}C_{re}^{\sigma^\nu}(\rho_{B|\Pi_\mu^a}),  \quad
\end{eqnarray}
where
\begin{eqnarray}
  C_{re}^{\sigma_\nu}(\rho) &=& S(\rho_d)-S(\rho).
\label{Cre1}
\end{eqnarray}
Here the relative entropy of coherence is used due to its clear physical meaning~\cite{you17emergent}, where
 $S(\rho)=-{\rm tr}(\rho \log_2\rho)$ stands for the von Neumann entropy of $\rho$ and $\rho_d$ is obtained from $\rho$ by removing all its off-diagonal entries.


Regarding $\mathbb{Z}_2$ and translation symmetries of Eq. (\ref{eq:Hamiltonian1}), in the bases spanned by the two-qubit product
states of eigenstate of $\sigma^z$, i.e., $\{\vert 0\rangle_i\otimes\vert 0\rangle_j$,
$\vert 0\rangle_i\otimes\vert 1\rangle_j$,
$\vert 1\rangle_i\otimes\vert 0\rangle_j$,
$\vert 1\rangle_i\otimes\vert 1\rangle_j\}$, where $\vert 0\rangle$ ($\vert 1\rangle$) denotes a spin-up (-down) state, the reduced density matrix $\rho_{ij}$ of two qubits $i$ and $j$ can be cast into an \emph{X}-state form,
\begin{equation}
\rho_{ij}=\left(\begin{array}{cccc}
  {u}^+&0&0&{z}_1\\
   0&{\omega}^+&{z}_2&0\\
    0&{z}_2^*&{\omega}^-&0\\
     {z}_1^*&0&0&{u}^-
  \end{array}\right),
  \label{jz}
\end{equation}
with
\begin{eqnarray}
&&{u}^{\pm}=\frac{1}{4}(1\pm2\langle\sigma_{i}^z\rangle+\langle\sigma_{i}^z\sigma_{j}^z\rangle), \\
&&{z}_1=\frac{1}{4}(\langle\sigma_{i}^x\sigma_{j}^x\rangle-\langle\sigma_{i}^y\sigma_{j}^y\rangle-i\langle\sigma_{i}^x\sigma_{j}^y\rangle-i\langle\sigma_{i}^y\sigma_{j}^x\rangle), \\
&&{z}_2=\frac{1}{4}(\langle\sigma_{i}^x\sigma_{j}^x\rangle+\langle\sigma_{i}^y\sigma_{j}^y\rangle+i\langle\sigma_{i}^x\sigma_{j}^y\rangle-i\langle\sigma_{i}^y\sigma_{j}^x\rangle), \\
&&{\omega}^{\pm}=\frac{1}{4}(1
-\langle\sigma_{i}^z\sigma_{j}^z\rangle).
\label{ys}
\end{eqnarray}

The SQC of two-qubit states as a function of $h$ for different $\alpha$ and  $r$ is plotted in Fig.~\ref{SQC1}.  $C_{re}^{st}$ shows a monotonic increase with respect to $h$, which
tends towards the maximum value 2.00, in contrast to a monotonically decreasing behavior of the relative entropy~\cite{XYDMYTC}. With increasing $h$, the SQC of two adjacent spins shows a smooth transition from the AFM phase to the PM phase, while $C_{re}^{st}$ exhibits a salient point across the transition from the spiral phase to the PM phase [cf. Fig.~\ref{SQC1}(a)].
For two qubits farther than the nearest neighbor, $C_{re}^{st}$ decreases with increasing $r$, but the positions of salient points are unchanged.
As is shown in Fig.~\ref{SQC1}(b), the nonanalyticity of the ground state at QCPs can be pinpointed by the discontinuity of the first-order derivative of the SQC. One can find that the coherence susceptibility, i.e., $\chi_{re}^{st}$ $\equiv$ $\partial C_{re}^{st}/\partial h$~\cite{chen2016coherence},  almost superposes onto each other around the QCPs for different $r$. Similarly, $\chi_{re}^{st}$ presents a pronounced peak at $h_{c,1}=1$ for $\alpha=0.50$ in Fig.~\ref{SQCr}(a) and the peaks become sharper and sharper as the system size increases, and it is expected to diverge in the thermodynamic limit. The singularity of $\chi_{re}^{st}$ is manifested in
the logarithmic scaling as
 \begin{eqnarray}
&&\left(\chi_{re}^{st}\right)_{\rm max}= a_{\rm C} \ln{N}+c_5, \label{Cst5}\\
&& \chi_{re}^{st}= b_{\rm C}\ln{|h-h_c|}+c_6.\label{Cst6} \quad
\end{eqnarray}
One further finds in Fig.~\ref{SQCr}(b) that the coherence susceptibilities $\chi_{re}^{st}$ for different $r$ also obey the logarithmic scaling.
The fitting parameters are listed in Table \ref{Fittingparameter}, and the extracted values of $\nu$ agree well with each other, although the deterioration of the precision with increasing $r$ can be easily noticed.

\begin{table}[bh!]
\caption{Fitting parameters $\{$$a_{\rm C}$, $b_{\rm C}$$\}$ of the slopes in logarithmic scaling of coherence susceptibility through relations Eqs.~(\ref{Cst5}) and (\ref{Cst6}) with $\nu=|a_{\rm C}/b_{\rm C}|$. Other parameters: $J=1.00$, $\gamma=0.60, \Gamma=0.60$, $\alpha= 0.50$.}
\label{Fittingparameter}
\begin{ruledtabular}
\begin{tabular}{ c c c c }
	$r$	&$a_{\rm C}$&$b_{\rm C}$&$\nu$ \\
\hline\\
1 &0.5084$\pm$0.0087&-0.5134$\pm$0.0091&0.9903$\pm$0.0355\\
2 &0.2065$\pm$0.0010&-0.2125$\pm$0.0049&0.9718$\pm$0.0272 \\
3 &0.0889$\pm$0.0027&-0.0933$\pm$0.0052& 0.9528$\pm$0.0777\\
\end{tabular}
\end{ruledtabular}
\end{table}



\section{Conclusion and Discussion}\label{summary}
In this work, we show that an elaborate scheme in atom-cavity systems can engineer an
effective Hamiltonian composed of various types of couplings, including \emph{XY}, Dzyaloshinskii-Moriya, and symmetric off-diagonal $\Gamma$ interactions.
We explore the quantum criticality in the so-called \emph{XY}-Gamma model, in which only nearest-neighbor interactions
between the particles are allowed. The
intricate interplay of diverse controlled exchange interactions between atoms in the presence of external fields
counteracts a rich variety of quantum phases at equilibrium.
The Hamiltonian can be rigorously solved through Jordan-Wigner and Bogoliubov transformation.
The exact solutions endow us with precise knowledge of ground-state properties.
For generic values of the parameters,
the phase diagram consists of the antiferromagnetic (AFM) phase, the paramagnetic (PM) phase,  and the gapless spiral phase. The second derivative of the ground-state energy diverges logarithmically across the AFM-PM transition, while it displays
a discontinuity at the critical point between the gapless spiral phase and gapped phases. Similar characteristics can be confirmed by the nearest-neighbor correlations and steered quantum coherence (SQC).
Moreover, we show that the gapless phase is characterized by a quasi-long-range order of an incommensurate spin spiral.
In a sense, the vector-chiral operator for two qubits is proven to act as a suitable order parameter for discerning the Tomonaga-Luttinger liquid. In contrast, the dimer order correlations vanish rapidly with the distance between the dimers.
The findings reveal that the incommensurate phase transitions away from the Tomonaga-Luttinger-liquid (TLL) phase all belong to the second-order transitions.
As a hallmark of critical phenomena in a continuous quantum phase transitions, critical points with the same set of critical exponents are categorized into a universality class, among which the dynamical exponent $z$ and the correlation-length exponent $\nu$ are most crucial.
For the quantum phase transitions between the AFM phase and the PM phase, one finds the correlation-length exponent can be obtained from several scaling forms including the gap and the correlation function as well as the SQC. The critical exponents $\nu = 1$ and $z = 1$ clearly indicate the transition  belong to two-dimensional Ising universality class. As regards the transitions from the TLL phase driven by either the off-diagonal exchange coupling ratio $\alpha$ or the magnetic field $h$, we obtain explicit forms for the energy gap near the critical points yielding
 $z$=2 and  $\nu$=1/2, signaling the critical points on the boundaries of the gapless phase are in the Lifshitz universality class. Thus, an experimental measurement of the correlation-length exponent $\nu$ and the dynamical exponent $z$ becomes tractable. The critical exponents $z$ and $\nu$ can be extracted from specific heat exponent $\alpha$ according to the hyperscaling relation $\nu+\nu z$=$2-\alpha$ or through
the Kibble–Zurek exponent $\mu=\nu/(1+\nu z)$~\cite{PhysRevB.99.094203,chepiga2021kibble}.  The finite-size effect theory developed for the Tomonaga-Luttinger liquid and the associated Lifshitz universality class speaks to the experimental verification for a finite number of atoms.

To conclude, the reported results may serve to test other approximate techniques used to study more realistic models.  This provides an interesting platform to understand the validity of characterizing tools in identifying unconventional transitions. The emergent phenomena in quantum many-body systems [Eq.~(\ref{eq:Hamiltonian0})] with cavity-induced long-range interactions await  further study. In particular, it becomes possible to realize nonequilibrium many-body phenomena in a controlled way~\cite{PhysRevLett.123.027204}, which are inaccessible for conventional solid-state materials. From both an experimental and a theoretical point of view, quantum simulations in the AMO systems offer outstanding possibilities for measuring  quantum coherence encoded in the many-body systems. Thus, the nature of the ground state and the coherence dynamics of the many-body systems can be unveiled, providing a hallmark of the TLL spin dynamics in the one-dimensional AFM chain. However, controlling atoms and coherence distillations with single-site resolution in an optical lattice remains a huge challenge.


\section*{ACKNOWLEDGMENTS}
The authors appreciate insightful discussions with Z.-A. Liu, X. Li, and T. Lv. W.-L. Y. is supported by the National Natural Science Foundation of China (NSFC) under Grant No.~12174194,
the startup fund (Grant No.~1008-YAH20006) of Nanjing University of Aeronautics and Astronautics (NUAA) , Top-notch Academic Programs Project of Jiangsu Higher Education Institutions (TAPP), and stable supports for basic institute research (Grant No.~190101). M. X. acknowledges support by the startup fund (Grant No.~1008-YAT21004) of NUAA
and the Open Research Fund Program of the State Key Laboratory of Low-Dimensional Quantum Physics (Grant No.~KF202111).
