# A comparison between higher-order nonclassicalities of SUP-engineered
coherent and thermal states

Deepak<EMAIL_ADDRESS>Arpita Chatterjee<EMAIL_ADDRESS>Department of Mathematics, J. C. Bose University of Science and Technology,
YMCA, Faridabad 121006, India

###### Abstract

We consider an experimentally obtainable SUP operator, defined by using a
generalized superposition of products of field annihilation ($a$) and creation
($a^{\dagger}$) operators of the type, $A=saa^{\dagger}+t{a^{\dagger}}a$ with
$s^{2}+t^{2}=1$. We apply this SUP operator on coherent and thermal quantum
states, the states thus produced are referred as SUP-operated coherent state
(SOCS) and SUP-operated thermal state (SOTS), respectively. In the present
work, we report a comparative study between the higher-order nonclassical
properties of SOCS and SOTS. The comparison is performed by using a set of
nonclassicality witnesses (e.g., higher-order antiubunching, higher-order sub-
Poissonian photon statistics, higher-order squeezing, Agarwal-Tara parameter,
Klyshko’s condition). The existence of higher-order nonclassicalities in SOCS
and SOTS have been investigated for the first time.  In view of possible
experimental verification of the proposed scheme, we present exact
calculations to reveal the effect of non-unit quantum efficiency of quantum
detector on higher-order nonclassicalities.

###### keywords:

Higher-order nonclassicality, superposition operator, superposed coherent
state, superposed thermal state, non-unit quantum efficiency

##  1 Introduction

It is a well-known fact that the superposition principle lies at the heart of
quantum mechanics [1]. Being a linear theory, quantum mechanics allows the
existence of various types of superposition states. However, linearity is not
a unique property of quantum mechanics only. Even in classical optics,
interference is being observed which is a clear evidence of superposition
principle via first order correlation [2]. But first order correlation or
first order interference are not strong enough to manifest some unique
(nonclassical) properties of quantum mechanics which have no classical
analogue. To realize the actual beauty, power and mystery of quantum
mechanics, we need to study the second or higher-order correlations.

A quantum state is identified as nonclassical if the corresponding Glauber-
Sudarshan $P$-function has negative values [3, 4]. But unfortunately, except a
single proposal for the measurement of the $P$-function in a special case,
there is no experimental method for determining the $P$-function [5].
Consequently, several operational measures to witness the nonclassicality have
been developed ([6] and references therein). Any finite set of these witnesses
of nonclassicality is only sufficient [7]. So, failure of any of these
criteria does not imply classicality of the state, but satisfying a condition
signifies nonclassicality. Almost all of these witnesses can be expressed as
the moments of annihilation and creation operators, and a nonclassical
property detected by a second-order correlation (i.e., a moment-based
criterion involving terms up to fourth order in annihilation and creation
operators) is called a lower-order nonclassical property. Naturally, higher-
order nonclassicality refers to the nonclassical features related to higher-
order correlations. Most frequently used higher-order nonclassical features
are higher-order antibunching (HOA) [8], higher-order sub-Poissonian photon
statistics (HOSPS) [9, 10], higher-order squeezing (HOS) of Hillery type [11]
and Hong-Mandel type [12]. All of these nonclassical features have lower-order
counterparts which have been studied more rigorously in [13, 14]. In addition
to these, Agarwal-Tara [15, 16] and Klyshko’s [17] conditions are two measures
of utmost interest. Recently, many investigations have been reported in the
direction of higher-order nonclassicalities because of a large number of
successful experimental characterization of such states [18, 19]. Also the
fact that higher-order criteria can identify weak nonclassicalities, not
detected by their lower-order counterparts, has turned this area a promising
field of research [20, 21].

The generation and manipulation of nonclassical states of the electromagnetic
field, at the single-photon level, have earned the attention of quantum optics
and quantum information community [22, 23], in the context of physical
realization of quantum tasks, protocols and communications [24, 25, 26, 27]
using continuous quantum variables. Different methods for preparing
nonclassical states have been suggested, based on operations of photon
addition ($a^{\dagger}$) [28] and subtraction ($a$) [29] on a classical field.
Agarwal and Tara [30] first proposed an $m$-photon-added scheme to create a
nonclassical state from any classical state. Zavatta et. al. [31] demonstrated
a single photon-added coherent state by homodyne tomography technology. In
another real scheme, Zavatta et. al. [32] used single-photon interference for
implementing superpositions of distinct quantum operations. Hu et. al. [33]
investigated the nonclassical properties of the field states generated by
subtracting any number of photons from a squeezed thermal state. Lee and Nha
[34] studied the action of an elementary coherent superposition of $a$ and
$a^{\dagger}$ on continuous variable systems. A SUP operator
$saa^{\dagger}+t{a^{\dagger}}a$ is proposed in [35], where $s$ and $t$ are
scalars with $s^{2}+t^{2}=1$. This experimentally feasible operation
introduced nonclassicality in the classical coherent state and non-Gaussianity
in the classical thermal state. The nonclassicality in the SUP-operated
coherent (SOCS) and thermal (SOTS) states can be quantified and analysed using
quasiprobability distributions in the phase space. However, to the best of our
knowledge, no effort has been made so far to investigate the higher-order
nonclassical properties of SOCS and SOTS. In this paper, we concentrate on
higher-order nonclassicalities of SUP-operated coherent and thermal states.

Apart from what stated above, the work is also motivated by the success of
detecting higher-order nonclassicality in laboratory using various
technologies [19, 20]. Recently, HOA, HOSPS and HOA have been reported in
optomechanical and optomechanical-like systems [36], finite dimensional
coherent state [37], most general form of the hyper Raman process [38], photon
added and subtracted squeezed coherent states [39] etc. In addition,
applications of nonclassical states [40] and effect of non-Gaussianity
inducing operations (e.g., photon addition, photon subtraction, photon
addition followed by subtraction and vice-versa) on the nonclassical
properties [39] have been studied in detail in the recent past. Keeping the
above in mind, we investigate the possibilities of observing HOA, HOSPS, HOS
of Hong-Mandel type etc. in the SUP-operated coherent and thermal states.  We
also discuss how imperfection in single-photon detector [41] influences the
characterization of higher-order nonclassicalities of SOCS and SOTS.

The structure of the paper is as follows: in Sec. 2, we have introduced SOCS
and SOTS and corresponding normalization constants. In Sec. 3, we have checked
the existence of various higher-order nonclassicalities in the superposition
states described in the previous section. Finally, the article is concluded in
Sec. 5.

##  2 Quantum states under consideration

In this section, we focus on SUP-operated coherent and thermal states. Being
an eigenstate of the annihilation operator $a$, the coherent state
corresponding to the unique eigenvalue $\alpha$, can be simply written as
$a\ket{\alpha}=\alpha\ket{\alpha}$. On the other hand, the Fock state
representation of a single-mode thermal field at an absolute temperature $T$
and with frequency $\omega$ is
$\rho_{n}=\sum_{n}a_{n}\ket{n}\bra{n}=\frac{1}{1+\bar{n}}\sum_{n}\Big{(}\frac{\bar{n}}{1+\bar{n}}\Big{)}^{n}\ket{n}\bra{n}$,
where $\bar{n}=[e^{h\omega/kT}-1]^{-1}$ is the mean photon number, $k$ being
the Boltzmann constant. Given a coherent (thermal) state as an input field,
the SUP-operated coherent (thermal) state can be described as [34]

$\displaystyle\left.\begin{array}[]{lcl}\ket{\psi_{1}}&=&N_{1}^{-1/2}A\ket{\alpha},\\\
\\\ \rho_{\psi_{2}}&=&N_{2}^{-1}A^{\dagger}\rho_{n}A,\end{array}\right\\}$ (4)

where $A=saa^{\dagger}+t{a^{\dagger}}a$, $s^{2}+t^{2}=1$ is the SUP operator,
$N_{1}=\Big{[}s^{2}+(1+2s^{2}+4st)|\alpha|^{2}+(1+2st)|\alpha|^{4}\Big{]}$ and
$N_{2}=\Big{[}s^{2}(1+\bar{n})(1+2\bar{n})+4st\bar{n}(1+\bar{n})+t^{2}\bar{n}(1+2\bar{n})\Big{]}$
are the normalization constants for SOCS and SOTS, respectively [35].

##  3 Nonclassical features of the superposed states

Any quantum state $\ket{\psi}$ is referred as nonclassical if the
corresponding Glauber-Sudarshan $P$-function is not a classical probability
distribution [3, 4]. That means the state having the negative $P$-function is
not of classical nature and can be considered as a nonclassical one. As $P$
function is very rarely detectable in laboratories [5], a number of
experimentally measurable nonclassicality witnesses like Mandel’s $Q_{M}$
parameter, higher-order antibunching (HOA), higher-order sub-Poissonian photon
statistics (HOSPS), higher-order squeezing (HOS) of Hong-Mandel type, Agarwal-
Tara and Klyshko’s conditions have been developed. These operational criterion
for identifying nonclassicality can be expressed as the moments of
annihilation and creation operators [37, 39]. To utilize these witnesses for
detecting the nonclassicalities of SOCS and SOTS, it is better to compute an
analytic expression for the most general moment $\langle{a^{\dagger
m}a^{n}}\rangle$, where $m$ and $n$ are the non-negative integers. This is the
most general moment in the sense that any other moment can be derived as a
special case of it. Using the algebraic properties of the boson operators [42]

$\displaystyle[a,a^{\dagger l}]=\frac{\partial a^{\dagger l}}{\partial
a^{\dagger}}$

and

$\displaystyle[a^{\dagger},a^{l}]=-\frac{\partial a^{l}}{\partial a}$

where $l$ is an integer, we obtain

$\displaystyle aa^{\dagger p}a^{q}a^{\dagger}$ $\displaystyle=$ $\displaystyle
aa^{\dagger p}(qa^{q-1}+a^{\dagger}a^{q})$ $\displaystyle=$ $\displaystyle
qaa^{\dagger p}a^{q-1}+aa^{\dagger{p+1}}a^{q}$ $\displaystyle=$ $\displaystyle
q(pa^{\dagger{p-1}}+a^{\dagger p}a)a^{q-1}+\big{\\{}(p+1)a^{\dagger
p}+a^{\dagger{p+1}}\big{\\}}a^{q}$ $\displaystyle=$ $\displaystyle
pqa^{\dagger{p-1}}a^{q-1}+(p+q+1)a^{\dagger p}a^{q}+a^{\dagger{p+1}}a^{q+1}$

Thus the expectation of $a^{\dagger m}a^{n}$ with respect to SOCS is given by

$\displaystyle{\langle{a^{\dagger m}a^{n}}\rangle}_{1}$ $\displaystyle=$
$\displaystyle N_{1}^{-1}\bra{\alpha}Aa^{\dagger m}a^{n}A\ket{\alpha}$
$\displaystyle=$ $\displaystyle
N_{1}^{-1}\alpha^{*m}\alpha^{n}\big{[}(s+t)^{2}\\{mn+(m+n+1)|\alpha|^{2}+|\alpha|^{4}\\}$
$\displaystyle+s(s+t)(m+n+2|\alpha|^{2})+s^{2}\big{]},$

while with respect to SOTS, it is

$\displaystyle{\langle{a^{\dagger m}a^{n}}\rangle}_{2}$ $\displaystyle=$
$\displaystyle N_{2}^{-1}\,\mathrm{tr}(a^{\dagger m}a^{n}A\rho_{r}A)$ (7)
$\displaystyle=$
$\displaystyle\delta_{mn}N_{2}^{-1}\sum_{r}\frac{a_{r}~{}r!}{(r-n)!}(s+(s+t)r)^{2},$

where $\delta_{mn}$ is the usual delta function.

###  3.1 Higher-order photon statistics

The Mandel’s parameter $Q_{M}$ [43] illustrates the nonclassicality of a
quantum state through its photon number distribution. The introductory
definition of $Q_{M}$ can be generalized to an arbitrary order $l$ as [44]

$\displaystyle Q_{M}^{(l)}$ $\displaystyle=$
$\displaystyle\frac{\langle{(\Delta{\mathcal{N}})^{l}}\rangle}{\langle{a^{\dagger}a}\rangle}-1,$
(8)

where $\Delta{\mathcal{N}}\,=\,a^{\dagger}a-\langle{a^{\dagger}a}\rangle$ is
the dispersion in the number operator $\mathcal{N}=a^{\dagger}a$. Using the
identity [44]

$\displaystyle\langle{(\Delta{\mathcal{N}})^{l}}\rangle=\sum_{k=0}^{l}{l\choose
k}(-1)^{k}\langle(a^{\dagger}a)^{l-k}\rangle{\langle a^{\dagger}a\rangle}^{k}$

and [45]

$(a^{\dagger}a)^{r}=\sum_{n=0}^{r}S_{r}^{(n)}a^{\dagger n}a^{n},$

where $S_{r}^{(n)}$ is the Stirling number of second kind [46]

$S_{r}^{(n)}=\frac{1}{n!}\sum_{j=0}^{n}(-1)^{n-j}{n\choose j}j^{r},$

the higher-order Mandel’s parameter $Q_{M}^{(l)}$ can be evaluated explicitly
up to order $l$. The negativity of $Q_{M}^{(2)}$ signifies the negativity of
the conventional Mandel’s $Q_{M}$. All expectations in (8) have been
calculated with help of (3) and (7).

Figure 1: (Color online) Comparison of $Q_{M}^{(l)}$ with respect to the
variable $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}}$ and $\bar{n}$ for thermal state $\rho_{\psi_{2}}$ and with
(a)-(c) $s=0.2$ and $l=2,5,7$, (d)-(f) $s=0.5$ and $l=2,5,7$, (g)-(i) $s=0.8$
and $l=2,5,7$, respectively. Figure 2: (Color online) Comparison of
$Q_{M}^{(l)}$ as a function of $s$ and with (a) $|\alpha|=\bar{n}=3$ and
$l=2$, (b) $|\alpha|=\bar{n}=3.37$ and $l=5$, (c) $|\alpha|=\bar{n}=3.69$ and
$l=7$, respectively.

The negative values of $Q_{M}^{(l)}$ parameter is a sufficient condition [16]
for the negativity of the $P$ function and hence it gives a witness for
nonclassicality. The photon number distribution is classified as sub-
poissionian, poissionian or superpoissionian if the Mandel’s parameter
$Q_{M}^{(l)}$ for any $l>1$, is less than, equal to or greater than zero,
respectively. In Fig. 1, a comparison between Mandel’s criteria of a SUP-
operated coherent state $\ket{\psi_{1}}$ with respect to the state parameter
$\alpha$ and a SUP-operated thermal state $\rho_{\psi_{2}}$ with respect to
the state parameter $\bar{n}$ is shown for different values of $s$. It is
clear that the nonclassicality, not detected by the lower-order parameter
$Q_{M}^{(2)}$, can be identified by the higher-order criteria (when $l=5$ or
$7$). Thus the presence of higher-order nonclassicality while its lower-order
counterpart is absent proves the relevance of the present study. It is also
notified that SOTS is more nonclassical than SOCS in a certain range of
$\alpha$ or $\bar{n}$ and after crossing it, SOCS displays more
nonclassicality than SOTS. The upper limit of this certain range increases
from $\alpha$ or $\bar{n}\approx 3.5$ to 3.8 as $l$ increases from $5$ to $7$,
keeping $s$ fixed. In Fig. 2, $Q_{M}^{(l)}$ is plotted as a function of $s$
and for different values of $\alpha$, $\bar{n}$ and $l$. Here SOCS shows mere
changes in nonclassicality with respect to $s$ while a notable change is
observed in case of SUP-operated thermal state.

###  3.2 Higher-order antibunching

HOA is another criteria which can be used to check the nonclassicality of any
quantum state. Lee introduced HOA by using the concept of majorization [47].
Later Pathak and Gracia [8] modified this to give a clear physical meaning and
to simplify the expression. The $(l-1)$-th order antibunching of a quantum
state is given as following

${\color[rgb]{0,0,0}d_{h}^{(l-1)}}=\langle a^{\dagger l}a^{l}\rangle-{\langle
a^{\dagger}a\rangle}^{l}\,\,<\,\,0$ (9)

Since the negativity of $d_{h}^{(l-1)}$ indicates that the probability of
photons coming bunched is less compared to that of coming independently,
therefore the nonclassicality feature (9) typifies how suitable a state
$\ket{\psi}$ is as a single photon resource. In this subsection, we have
calculated the values of $d_{h}^{(l-1)}$ for SOCS and SOTC using (3) and (7)
and plotted them. The signature of lower-order antibunching can be obtained as
a special case of (9) for $l=2$ and for $l\geq 3$, the negative values of
$d_{h}^{(l-1)}$ correspond to the higher-order antibunching of order $(l-1)$.

Figure 3: (Color online) Comparison of $d_{h}^{(l-1)}$ with respect to the
variable $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}}$ and $\bar{n}$ for thermal state $\rho_{\psi_{2}}$ and with
(a)-(c) $s=0.2$ and $l=2,3,4$, (d)-(f) $s=0.5$ and $l=2,3,4$, (g)-(i) $s=0.8$
and $l=2,3,4$, respectively. Figure 4: (Color online) Comparison of
$d_{h}^{(l-1)}$ as a function of $s$, for coherent state $\ket{\psi_{1}}$ and
thermal state $\rho_{\psi_{2}}$ with (a) $|\alpha|=\bar{n}=0.1$, $l=2$, (b)
$|\alpha|=\bar{n}=0.2$, $l=3$ and (c) $|\alpha|=\bar{n}=0.3$, $l=4$.

In Fig. 3, the variation in antibunching of SOCS (SOTS) is shown with respect
to the state variable $\alpha$ ($\bar{n}$). Here lower as well as higher-order
antibunching criteria demonstrate the nonclassicality for both the states. It
is observed that the depth of nonclassicality witness increases with the order
of nonclassicality as depicted in Figures 3 (a)-(c). This fact is consistent
with the earlier observations ([38, 39] and references therein) that higher-
order nonclassicality criteria is useful in detecting weaker nonclassicality.
Specifically, the depth of nonclassicality present in SOCS and SOTS, lowers as
$s$ increases from 0.2 to 0.8 keeping $l$ fixed [see Figures 3 (b), (e), (h)
and 3(c), (f), (i)].

###  3.3 Higher-order sub-Poissonian photon statistics

HOSPS is one of the important criteria of nonclassicality and sometimes gives
better result than HOA. It is further investigated that HOA and HOSPS do not
depend on each other [38] and both of them can exist irrespective of whether
their lower-order counterparts exist or not [37].

The generalized criteria for observing the $(l-1)$-th order sub-Poissonian
photon statistics (for which
$\langle(\Delta\mathcal{N})^{l}\rangle<\langle(\Delta\mathcal{N})^{l}\rangle_{\ket{\mathrm{Poissonian}}}$
is given by [10]

$\mathcal{D}_{h}^{(l-1)}=\sum_{e=0}^{l}\sum_{f=1}^{e}S_{2}(e,f)^{l}C_{e}(-1)^{e}d_{h}^{(f-1)}{\langle
a^{\dagger}a\rangle}^{l-e}\,\,<\,\,0$ (10)

where $S_{2}(e,f)=\sum_{r=0}^{f}{{}^{f}C_{r}}(-1)^{r}r^{e}$ is the Stirling
number of second kind, ${}^{l}C_{e}$ is the usual binomial coefficient. The
analytic expressions of HOSPS for the SOCS and SOTS can be obtained by
substituting (3) and (7) in (10).

Figure 5: (Color online) Comparison of $D_{h}^{(l-1)}$ with respect to the
variable $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}}$ and $\bar{n}$ for thermal state $\rho_{\psi_{2}}$ and with
(a)-(c) $s=0.2$ and $l=2,3,4$, (d)-(f) $s=0.5$ and $l=2,3,4$, (g)-(i) $s=0.8$
and $l=2,3,4$, respectively. Figure 6: (Color online) Comparison of
$D_{h}^{(l-1)}$ as a function of $s$, for coherent state $\ket{\psi_{1}}$ and
thermal state $\rho_{\psi_{2}}$ with $|\alpha|=\bar{n}=0.2$ and (a) $l=2$, (b)
$l=3$ and (c) $l=4$.

$\mathcal{D}_{h}^{(l-1)}$ is plotted in Fig. 5 with respect to the state
variable $\alpha$ ($\bar{n}$) for SUP-operated coherent (thermal) state and
for different values of $l$ and $s$. HOSPS clearly detects nonclassicality for
both the states under consideration. The lower-order sub-Poissonian photon
statistics is merely present in SOTS unlike lower-order antibunching of SOTS
which is detected for all values of $s$ considered here. It is also noticed
that HOSPS behaves in a similar manner as HOA does while changing $l$ and $s$.

###  3.4 Higher-order squeezing

The criterion of higher-order squeezing is given by the pioneering work of
Hong and Mandel [12]. According to them, the higher-order squeezing of order
$l$ ($>2$) is obtained while the $l$-th order moment of a field quadrature
operator is less than the corresponding coherent state value. Hong-Mandel’s
criteria for higher-order squeezing can be described by the following
inequality

$S^{(l)}=\frac{\langle(\Delta
X)^{l}\rangle-{\left(\frac{1}{2}\right)}_{\left(\frac{l}{2}\right)}}{{\left(\frac{1}{2}\right)}_{\left(\frac{l}{2}\right)}}\,\,<\,\,0,$
(11)

where $(x)_{l}$ is the conventional Pochhammer symbol and the quadrature
variable is defined as $X=\frac{1}{\sqrt{2}}(a+a^{\dagger})$. The inequality
in (11) can also be rewritten as

$\langle(\Delta
X)^{l}\rangle\,\,<\,\,{\left(\frac{1}{2}\right)}_{\left(\frac{l}{2}\right)}=\frac{1}{2^{\frac{l}{2}}}(l-1)!!,$
(12)

with

$\displaystyle\langle(\Delta
X)^{l}\rangle=\sum_{r=0}^{l}\sum_{i=0}^{\frac{r}{2}}\sum_{k=0}^{r-2i}(-1)^{r}\frac{1}{2^{\frac{1}{2}}}(2i-1)!^{2i}C_{k}^{l}C_{r}^{r}C_{2i}\langle
a^{\dagger}+a\rangle^{l-r}\langle a^{\dagger k}a^{r-2i-k}\rangle,$ (13)

where $l$ is an even number and

$\displaystyle n!!=\left\\{\begin{array}[]{lll}&n(n-2)(n-4)\ldots 4.2&\mbox{if
$n$ is even},\\\ \\\ &n(n-2)(n-4)\ldots 3.1&\mbox{if $n$ is
odd},\end{array}\right.$

The analytic expression for the Hong-Mandel type HOS can be obtained by using
(4) in (11)-(13).

Figure 7: (Color online) Comparison of $S^{(l)}$ with respect to the variable
$\gamma$ which is equal to $|\alpha|$ for coherent state $\ket{\psi_{1}}$ and
$\bar{n}$ for thermal state $\rho_{\psi_{2}}$ and with (a)-(c) $s=0.2$ and
$l=2,4,6$, (d)-(f) $s=0.5$ and $l=2,4,6$, (g)-(i) $s=0.8$ and $l=2,4,6$,
respectively. Figure 8: (Color online) Comparison of $S^{(l)}$ as a function
of $s$, for coherent state $\ket{\psi_{1}}$ and thermal state
$\rho_{\psi_{2}}$ with $|\alpha|=\bar{n}=0.01$ and (a) $l=2$, (b) $l=4$ and
(c) $l=6$, respectively.

Fig. 7 illustrates the existence of Hong-Mandel type HOS in the superposed
states, assuming arbitrary values of the state variables $\alpha$ or $\bar{n}$
and for different orders of squeezing ($l=2,\,4,\,6$). Here the lower-order
squeezing criteria fails to detect the nonclassicality of SUP-operated states
[see Figures 7 (a), (d), (g)] while its higher-order counterpart hardly
witnesses the nonclassicality of SOTS. Thus HOS is not much advantageous for
identifying the nonclassicalities of the considered states.

###  3.5 $Q$ function

A direct phase space description of a quantum mechanical system is not
possible due to the uncertainty principle. This fact leads to the construction
of quasiprobability distributions which are very useful in quantum mechanics
as they provide a quantum classical correspondence and facilitate the
calculation of quantum mechanical averages in close analogy to classical phase
space averages [39]. One such quasiprobability distribution is $Q$ function,
and zeros of this function are a signature of nonclassicality [48]. The
Husimi-$Q$ function is defined as

$Q=\frac{1}{\pi}\bra{\beta}\rho\ket{\beta},$

where $\ket{\beta}$ is the usual coherent state. A simple calculation finds
the Husimi-$Q$ for SOCS and SOTS as following:

$\displaystyle\begin{array}[]{rcl}Q_{1}&=&\frac{1}{\pi}N_{1}^{-1}|\langle{\beta|\psi_{1}}\rangle|^{2}\\\
\\\
&=&\frac{1}{\pi}N_{1}^{-1}\left|s+(s+t)\alpha\beta^{*}\right|^{2}e^{-|\alpha|^{2}-|\beta|^{2}+\alpha\beta^{*}+\alpha^{*}\beta},\end{array}$
(18)
$\displaystyle\begin{array}[]{rcl}Q_{2}&=&\frac{1}{\pi}N_{2}^{-1}\langle{\beta|\rho_{\psi_{2}}|\beta}\rangle\\\
\\\
&=&\frac{1}{\pi(1+\bar{n})}N_{2}^{-1}e^{-\frac{|\beta|^{2}}{1+\bar{n}}}\left[\left\\{s+(s+t)\left(\frac{\bar{n}|\beta|^{2}}{1+\bar{n}}\right)\right\\}^{2}+(s+t)^{2}\left(\frac{\bar{n}|\beta|^{2}}{1+\bar{n}}\right)\right]\end{array}$
(22)

The zeros of Husmi $Q$ function for the state $\ket{\psi_{1}}$ can be
calculated by using the condition $s+(s+t)\alpha\beta^{*}=0$ which gives
$s=\frac{\alpha\beta^{*}}{\sqrt{1+2\alpha\beta^{*}+2(\alpha\beta^{*})^{2}}}$.

Figure 9: (Color online) Comparison of Husimi-$Q$ as a function of $\beta$,
for different values of state variables $|\alpha|$ for coherent state
$\ket{\psi_{1}}$ and $\bar{n}$ for thermal state $\rho_{\psi_{2}}$ with (a)
$|\alpha|=\bar{n}=1$ and $s=0.2$, (b) $|\alpha|=\bar{n}=0.1$ and $s=0.5$, (c)
$|\alpha|=\bar{n}=0.01$ and $s=0.8$.

Fig. 9 shows that, in general we have failed to observe the nonclassical
features reflected beyond moments based nonclassicality criteria through a
quasiprobability distribution, i.e., zeros of the $Q$ function. Particularly
when $s=0.2$ and state variable $\alpha=1$ ($\bar{n}=1$), the $Q$ function has
zero for SOCS (SOTS). Thus Husimi-$Q$ depicts the nonclassical nature of both
the states in some specific region.

###  3.6 Agarwal-Tara criterion

Agarwal and Tara proposed a moment based criteria for witnessing the
nonclassical character of a given quantum state [15]. They defined $A_{3}$
which consists of the moments of the number distribution $\mu_{j}$ and the
normal ordered moments $m_{j}$. The analytic expression of $A_{3}$ in terms of
these higher ordered moments is [16]

$A_{3}=\frac{{\mbox{det}\,\,m}^{(3)}}{{\mbox{det}\,\,\mu}^{(3)}-{\mbox{det}\,\,m}^{(3)}}\,\,<\,\,0,$
(23)

where

$\displaystyle m^{(3)}=\begin{bmatrix}1&m_{1}&m_{2}\\\ m_{1}&m_{2}&m_{3}\\\
m_{2}&m_{3}&m_{4}\end{bmatrix},$

and

$\displaystyle\mu^{(3)}=\begin{bmatrix}1&\mu_{1}&\mu_{2}\\\
\mu_{1}&\mu_{2}&\mu_{3}\\\ \mu_{2}&\mu_{3}&\mu_{4}\end{bmatrix}.$

The matrix elements are defined by $m_{j}=\langle{a^{\dagger j}a^{j}}\rangle$
and $\mu_{j}=\langle{(a^{\dagger}a)^{j}}\rangle$.

Figure 10: (Color online) Variation of Agarwal-Tara parameter with respect to
the state parameter $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}}$ and $\bar{n}$ for thermal state $\rho_{\psi_{2}}$ with (a)
$s=0.01$ (b) $s=0.1$ (c) $s=0.2$, respectively. Figure 11: (Color online)
Variation of Agarwal-Tara parameter with respect to $s$ for coherent state
$\ket{\psi_{1}}$ and thermal state $\rho_{\psi_{2}}$ where (a)
$|\alpha|=\bar{n}=0.1$ (b) $|\alpha|=\bar{n}=0.5$ and (c)
$|\alpha|=\bar{n}=0.7$, respectively.

For a classical-like coherent state, the value of $A_{3}$ is 0 while for the
most nonclassical Fock state, it is -1. Thus for any nonclassical state,
$A_{3}$ is negative and bounded below by -1 (if $A_{3}=-1$, the corresponding
state is called maximally nonclassical). To investigate the nonclassicality of
the considered state via the Agarwal-Tara parameter, we have plotted $A_{3}$
[see Fig. 10] as a function of $\alpha$ ($\bar{n}$) for SOCS (SOTS). Here
$A_{3}$ lies between 0 and -1 for all values of $s$. The zero value of $A_{3}$
of SOCS for the displacement parameter $\alpha=0$ is consistent with the fact
that vacuum state is a classical state having non-negative $P$-function. For
some specific mean values, Agarwal-Tara criterion for SOCS shows a singularity
[cf. Figures 10-11], which indicates that the probability of photons coming
bunched is equal to that of coming independently. Nonclassicality reflected
through $A_{3}$ parameter decreases (increases) while increasing the SUP-
operator parameter $s$, and the mean photon number $\alpha$ ($\bar{n}$) for
SOCS (SOTS).

###  3.7 Klyshko’s criterion

One of the witness to check the nonclassicality of any quantum state was given
by Klyshko [17], which is based only on three consecutive photon-number
probabilities. If $p_{m}=\langle{{m}|\rho|{m}}\rangle$ is the photon-number
probability of a state with density matrix $\rho$, then the Klyshko’s
inequality can be written as

$B(m)=(m+2)p_{m}p_{m+2}-(m+1)(p_{m+1})^{2}\,\,<\,\,0$ (24)

Using (4) and (24), the functions $B(m)$ for both the states are obtained as
follows:

$\displaystyle\begin{array}[]{lcl}B(m)_{1}&=&N_{1}^{-2}|\alpha|^{4m+4}\left[(m+2)\left\\{s+(s+t)m\right\\}^{2}\left\\{s+(s+t)(m+2)\right\\}^{2}\right.\\\
&&\left.-(m+1)\left\\{s+(s+t)(m+1)\right\\}^{4}\right]\end{array}$ (27)
$\displaystyle\begin{array}[]{lcl}B(m)_{2}&=&N_{2}^{-2}\left(\frac{1}{1+\bar{n}}\right)^{2}\left(\frac{\bar{n}}{1+\bar{n}}\right)^{2m+2}\left[(m+2)\left\\{s+(s+t)m\right\\}^{2}\right.\\\
&&\left.\times\left\\{s+(s+t)(m+2)\right\\}^{2}-(m+1)\left\\{s+(s+t)(m+1)\right\\}^{4}\right]\end{array}$
(30)

where $B(m)_{1}$, $B(m)_{2}$ are corresponding to SOCS and SOTS respectively.

Figure 12: (Color online) Illustration of Klyshko’s criterion $B(m)$ as a
function of $m$ for coherent state $\ket{\psi_{1}}$ and thermal state
$\rho_{\psi_{2}}$ with (a) $|\alpha|=\bar{n}=2$, $s=0.2$, (b)
$|\alpha|=\bar{n}=1.5$, $s=0.5$, (c) $|\alpha|=\bar{n}=1$, $s=0.8$.

In this criterion, very less information is needed with respect to the other
conditions as probability of only three successive photon numbers is
sufficient to investigate this nonclassical property. The negative values of
$B(m)$ indicate the nonclassicality of the corresponding state. It is clear
from Fig. 12 that Klyshko’s condition signifies nonclassicality only for
$m=0,\,1$ and as $s$ increases from 0.2 to 0.8, $B(m)$ fails to detect the
nonclassicality.

Figure 13: (Color online) Illustration of all criteria values ($C_{v}$) as a
function of (a) $|\alpha|$ for coherent state $\ket{\psi_{1}}$ and (b)
$\bar{n}$ for thermal state $\rho_{\psi_{2}}$ with arbitrary parametric values
$s=0.5$, $l=4$ and for Klyshko’s criterion $m=4$.

Fig. 13 describes a comparative plot for all the seven criteria studied so
far, as a function of $\alpha$ for coherent state $\ket{\psi_{1}}$ and
$\bar{n}$ for thermal state $\rho_{\psi_{2}}$ with arbitrary parametric values
$s=0.5$ and $l=4$. Almost all of the witnesses signify the nonclassical nature
of SUP-operated coherent state and SUP-operated thermal state except $S^{(l)}$
for SOTS.

##  4 Effect of non-unit quantum efficiency

Light comes in the form of photons. The quantum efficiency states that for
every photon coming in there is a probability that the photon will generate an
electron which is measured by the detector. An ideal single-photon detector is
a measurement device that perfectly discriminates the photon-number states,
i.e. a box that produces $n$ clicks when the number state $\ket{n}$ comes in
[49]. The vacuum state $\ket{0}$ can be absolutely extricated from all the
other Fock states $\ket{n_{n\geq 1}}$ by a unit-efficiency and noise-free
detector [50]. To account for the non-unit efficiency $\eta$, the mode $a$ is
sent through a beam-splitter with transmission probability $\eta$. The output
modes are labelled as $a$ and $l$ and the corresponding operator is given by
$C_{\eta}=e^{{\sqrt{1-\eta}/\eta}al^{\dagger}}{\eta}^{\frac{1}{2}a^{\dagger}a}\ket{0_{l}}$.
Thus an inefficient detector is well-described by the operator [41]

$\displaystyle D(\eta)=(1-\eta)^{a^{\dagger}a}$ (31)

Applying the operator $D(\eta)$ on SOCS and SOTS, we obtain $D(\eta)$-operated
SOCS (DSOCS) and $D(\eta)$-operated SOTS (DSOTS) as follows

$\displaystyle\left.\begin{array}[]{lcl}\ket{\psi_{1}^{\eta}}&=&(N_{1}^{\eta})^{-1/2}D(\eta)A\ket{\alpha},\\\
\\\
\rho_{\psi_{2}^{\eta}}&=&(N_{2}^{\eta})^{-1}D^{\dagger}(\eta)A^{\dagger}\rho_{n}AD(\eta),\end{array}\right\\}$
(35)

where

$\displaystyle N_{1}^{\eta}$ $\displaystyle=$
$\displaystyle\left[\left\\{s+(s+t)|\alpha|^{2}\right\\}^{2}+(s+t)^{2}|\alpha|^{2}\right]+\left[\left\\{s+(s+t)(1-\eta)|\alpha|^{2}\right\\}^{2}+(s+t)^{2}(1-\eta)|\alpha|^{2}\right]e^{-\eta|\alpha|^{2}}$
$\displaystyle+\left[\left\\{s+(s+t)(1-\eta)^{2}|\alpha|^{2}\right\\}^{2}+(s+t)^{2}(1-\eta)^{2}|\alpha|^{2}\right]e^{(\eta^{2}-2\eta)|\alpha|^{2}},$

and

$\displaystyle N_{2}^{\eta}$ $\displaystyle=$
$\displaystyle\frac{1}{(1+\bar{n})}\left[\left[\left\\{s+(s+t)\mu^{2}\right\\}^{2}+(s+t)^{2}\mu^{2}\right]e^{\mu}-2\left[\left\\{s+(s+t)\theta^{2}\right\\}^{2}+(s+t)^{2}\theta^{2}\right]e^{\theta}\right.$
$\displaystyle\left.+\left[\left\\{s+(s+t)\zeta^{2}\right\\}^{2}+(s+t)^{2}\zeta^{2}\right]e^{\zeta}\right],$

are the normalization constants for DSOCS and DSOTS, respectively and
$\mu=\frac{\bar{n}}{1+\bar{n}}$, $\theta=\mu(1-\eta)$,
$\zeta=\mu(1-\eta)^{2}$. The expectation of $a^{\dagger m}a^{n}$ with respect
to DSOCS is given by

$\displaystyle{\langle{a^{\dagger m}a^{n}}\rangle}_{1}^{\eta}$
$\displaystyle=$ $\displaystyle\bra{\psi_{1}^{\eta}}a^{\dagger
m}a^{n}\ket{\psi_{1}^{\eta}}$ (36) $\displaystyle=$
$\displaystyle(N_{1}^{\eta})^{-1}\alpha^{*m}\alpha^{n}\left[\left\\{\big{(}s+(s+t)|\alpha|^{2}\big{)}^{2}+(s+t)^{2}|\alpha|^{2}\right\\}+\left\\{(1-\eta)^{m}+(1-\eta)^{n}\right\\}\right.$
$\displaystyle\left.\times\left\\{\big{(}s+(s+t)(1-\eta)|\alpha|^{2}\big{)}^{2}+(s+t)^{2}(1-\eta)|\alpha|^{2}\right\\}e^{-\eta|\alpha|^{2}}+\big{(}1-\eta\big{)}^{m+n}\right.$
$\displaystyle\left.\times\left\\{\big{(}s+(s+t)(1-\eta)^{2}|\alpha|^{2}\big{)}^{2}+(s+t)^{2}(1-\eta)^{2}|\alpha|^{2}\right\\}e^{(\eta^{2}-2\eta)|\alpha|^{2}}\right],$

and with respect to DSOTS, it is

$\displaystyle{\langle{a^{\dagger m}a^{n}}\rangle}_{2}^{\eta}$
$\displaystyle=$ $\displaystyle(N_{2}^{\eta})^{-1}\,\mathrm{tr}(a^{\dagger
m}a^{n}D^{\dagger}(\eta)A^{\dagger}\rho_{r}AD(\eta))$ (37) $\displaystyle=$
$\displaystyle\delta_{mn}(N_{2}^{\eta})^{-1}\sum_{r}\frac{a_{r}~{}r!}{(r-n)!}(s+(s+t)r)^{2}\left\\{1-(1-\eta)^{r}\right\\}^{2},$

where $\delta_{mn}$ is the usual delta function. Using (36) and (37), the
analytic expressions of several higher-order nonclassicality criterion for
DSOCS and DSOTS are evaluated and plotted as follows:

Figure 14: (Color online) Comparison of $Q_{M}^{(l)}$ with respect to the
variable $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}^{\eta}}$ and $\bar{n}$ for thermal state
$\rho_{\psi_{2}^{\eta}}$ and with $s=0.2$ (a) $l=2$, $\eta=0.25$, (b) $l=5$,
$\eta=0.5$, (c) $l=7$, $\eta=0.75$, respectively. Figure 15: (Color online)
Comparison of $d_{h}^{(l-1)}$ with respect to the variable $\gamma$ which is
equal to $|\alpha|$ for coherent state $\ket{\psi_{1}^{\eta}}$ and $\bar{n}$
for thermal state $\rho_{\psi_{2}^{\eta}}$ and with $s=0.2$ (a) $l=2$,
$\eta=0.75$, (b) $l=3$, $\eta=0.5$,(c) $l=4$, $\eta=0.25$, respectively.
Figure 16: (Color online) Comparison of $D_{h}^{(l-1)}$ with respect to the
variable $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}^{\eta}}$ and $\bar{n}$ for thermal state
$\rho_{\psi_{2}^{\eta}}$ and with $s=0.2$ and (a) $l=2$, $\eta=0.3$, (b)
$l=3$, $\eta=0.6$ and (c) $l=4$, $\eta=0.9$, respectively. Figure 17: (Color
online) Comparison of $S^{(l)}$ with respect to the variable $\gamma$ which is
equal to $|\alpha|$ for coherent state $\ket{\psi_{1}^{\eta}}$ and $\bar{n}$
for thermal state $\rho_{\psi_{2}^{\eta}}$ and with $s=0.2$ and (a) $l=2$,
$\eta=0.4$, (b) $l=4$, $\eta=0.6$ and (c) $l=6$, $\eta=0.8$, respectively.
Figure 18: (Color online) Comparison of Husimi-$Q$ as a function of $\beta$,
for different values of state variables $|\alpha|$ for coherent state
$\ket{\psi_{1}^{\eta}}$ and $\bar{n}$ for thermal state
$\rho_{\psi_{2}^{\eta}}$ with (a) $|\alpha|=\bar{n}=1$, $s=0.2$ and
$\eta=0.1$, (b) $|\alpha|=\bar{n}=0.1$, $s=0.5$ and $\eta=0.4$, (c)
$|\alpha|=\bar{n}=0.01$, $s=0.8$ and $\eta=0.7$, respectively. Figure 19:
(Color online) Variation of Agarwal-Tara parameter with respect to the state
parameter $\gamma$ which is equal to $|\alpha|$ for coherent state
$\ket{\psi_{1}^{\eta}}$ and $\bar{n}$ for thermal state
$\rho_{\psi_{2}^{\eta}}$ with (a) $s=0.01$, $\eta=0.8$, (b) $s=0.1$,
$\eta=0.5$, (c) $s=0.2$, $\eta=0.2$, respectively. Figure 20: (Color online)
Illustration of Klyshko’s criterion $B(m)$ as a function of $m$ for coherent
state $\ket{\psi_{1}^{\eta}}$ and thermal state $\rho_{\psi_{2}^{\eta}}$ with
(a) $|\alpha|=\bar{n}=2$, $s=0.2$, $\eta=0.4$, (b) $|\alpha|=\bar{n}=1.5$,
$s=0.5$, $\eta=0.6$, (c) $|\alpha|=\bar{n}=1$, $s=0.8$, $\eta=0.8$,
respectively.

Here, one can visualize that due to the imperfections of real detector, the
depth of nonclassicality of the superposed states increases in case of higher-
order photon statistics. The same is true for HOA, HOSPS and Husimi-$Q$. We
may note that HOS is not observed for the SUP-operated thermal state while
non-unit quantum efficiency is taken into account. Additionally, the negative
peaks in the values of $B(m)$ disappear. That means, the typical realistic
detector can control higher-order nonclassicalities of the states studied
here.

##  5 Conclusion

In the present work, we have introduced a quantum state by applying a
combination of product of two operators $a$ and $a^{\dagger}$ to a coherent
(thermal) state $\ket{\psi_{1}}$ ($\rho_{\psi_{2}}$) with $s^{2}+t^{2}=1$. The
scalars $s$ and $t$ act as control parameters for manipulation of the
nonclassical character of the output state. We have focused on the higher-
order nonclassical features of the states. A set of different measurement
techniques is used to check the existence of higher-order nonclassicality in
the superposed states. It is found that higher-order measurements (except HOS
and Klyshko) perform better to detect the nonclassicality of SOCS and SOTS. In
brief, various types of higher-order nonclassical features of SUP-operated
coherent as well as thermal states are investigated and their variations with
experimentally controllable parameters are studied. We have shown that the
depth of higher-order nonclassicality witnesses increase with the order of
nonclassicality which establishes the fact that higher-order nonclassicality
criteria has an advantage in detecting weaker nonclassicality. The effect of
non-unit quantum efficiency of photon detector is also discussed which ensures
that the proposed work can be realized experimentally in the near future. We
can wind up this article with a hope that the allowed control over the depth
of nonclassicality witnesses will be put to further applications in the
context of quantum information processing.

ACKNOWLEDGEMENT

Deepak’s work is supported by the Council of Scientific and Industrial
Research (CSIR), Govt. of India (Award no. 09/1256(0006)/2019-EMR-1).

## References

  * [1] V. Bužek and P. L. Knight, Prog. Opt. 34, 1 (1995).
  * [2] A. W. Harrow and A. Montanaro, Nature 549(7671), 203 (2017).
  * [3] R. J. Glauber, Phys. Rev. 131(6), 2766 (1963).
  * [4] E. C. G. Sudarshan, Phys. Rev. Lett. 10(7), 277 (1963).
  * [5] T. Kiesel, W. Vogel, V. Parigi, A. Zavatta and M. Bellini, Phys. Rev. A 78(2), 021804 (2008).
  * [6] A. Miranowicz, M. Bartkowiak, X. Wang, Y. X. Liu and F. Nori, Phys. Rev. A 82(1), 013824 (2010).
  * [7] Th. Richter and W. Vogel, Phys. Rev. Lett. 89(28), 283601 (2002).
  * [8] A. Pathak and M. E. Garcia, Appl. Phys. B 84(3), 479 (2006).
  * [9] H. Prakash and D. K. Mishra, J. Phys. B: At. Mol. Opt. Phys. 39(9), 2291 (2006).
  * [10] A. Verma and A. Pathak, Phys. Lett. A 374(8), 1009 (2010).
  * [11] M. Hillery, Phys. Rev. A 36(8), 3796 (1987).
  * [12] C. K. Hong and L. Mandel, Phys. Rev. Lett. 54(4), 323 (1985).
  * [13] A. Pathak and A. Ghatak, J. Electro. Waves and Appl. 32(2), 229 (2018).
  * [14] V. V. Dodonov, J. Opt. B: Quant. and Semiclass. Opt. 4(1), R1 (2002).
  * [15] G. S. Agarwal and K. Tara, Phys. Rev. A 46(1), 485 (1992).
  * [16] P. Malpani, K. Thapliyal, N. Alam, A. Pathak, V. Narayanan and S. Banerjee, Opt. Commun. 459, 124964 (2020).
  * [17] D. N. Klyshko, Phys. Lett. A 213(1-2), 7 (1996).
  * [18] M. Avenhaus, K. Laiho, M. V. Chekhova and Ch. Silberhorn, Phys. Rev. Lett. 104(6), 063602 (2010).
  * [19] A. Allevi, S. Olivares and M. Bondani, Intl. J. Quant. Inform. 10(08), 1241003 (2012).
  * [20] A. Allevi, S. Olivares and M. Bondani, Phys. Rev. A 85(6), 063835 (2012).
  * [21] J. Peřina Jr., V. Michálek and O. Haderka, Phys. Rev. A 96(3), 033852 (2017).
  * [22] P. Kok, Nature Photon. 4, 504 (2010).
  * [23] S. L. Braunstein and P. Van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * [24] A. Zeilinger, Rev. Mod. Phys. 71, S288 (1998).
  * [25] D. Bouwmeester, J. W. Pan, M. Daniell, H. Weinfurter, M. Zukowski and A. Zeilinger, Nature 394, 841 (1998).
  * [26] D. Bouwmeester, A. Ekert and A. Zeilinger, The Physics of Quantum Information (Berlin: Springer), (2000).
  * [27] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge: Cambridge University Press), (2000).
  * [28] A. Zavatta, S. Viciani and M. Bellini, Science 306, 660 (2004).
  * [29] J. Wenger, T. Brouri and P. Grangier, Phys. Rev. Lett. 92, 153601 (2004).
  * [30] G. S. Agarwal and K. Tara, Phys. Rev. A 43, 492 (1991).
  * [31] A. Zavatta, S. Viciani and M. Bellini, Phys. Rev. A 72, 023820 (2005).
  * [32] A. Zavatta, V. Parigi, M. S. Kim, H. Jeong and M. Bellini, Phys. Rev. Lett. 103, 140406 (2009).
  * [33] L. Y. Hu, X. X. Xu, Z. S. Wang and X. F. Xu, Phys. Rev. A 82, 043842 (2010).
  * [34] S. Y. Lee and H. Nha, Phys. Rev. A 82, 053812 (2010).
  * [35] A. Chatterjee, H. S. Dhar and R. Ghosh, J. Phys. B: At. Mol. Opt. Phys. 45, 205501 (2012).
  * [36] N. Alam, K. Thapliyal, A. Pathak, B. Sen, A. Verma and S. Mandal, arXiv:1708.03967, (2017).
  * [37] N. Alam, A. Verma and A. Pathak, Phys. Lett. A 382(28), 1842 (2018).
  * [38] K. Thapliyal, A. Pathak, B. Sen and J. Peřina, Phys. Rev. A 90(1), 013808 (2014).
  * [39] K. Thapliyal, N. L. Samantray, J. Banerji and A. Pathak, Phys. Lett. A 381(37), 3178 (2017).
  * [40] G. S. Agarwal, Quantum Optics, Cambridge University Press, (2013).
  * [41] N. Gisin et. al., J. Phys. B: At. Mol. Opt. Phys. 45(12), 124016 (2011).
  * [42] W. H. Louisell, Quantum Statistical Properties of Radiation, Wiley Classics Publications, (1990).
  * [43] L. Mandel, Opt. Lett. 4(7), 205 (1979).
  * [44] S. Dey and V. Hussin, Phys. Rev. A, 93, 053824 (2016).
  * [45] J. M. Vargas Martínez, H. Moya-Cessa and M. F. Guasti, Rev. Mexi. De Física E 52(1), 13 (2006).
  * [46] M. Abramowitz and I. A. Stegun, Handbook of Mathematical functions (Dover, New York, 1968).
  * [47] C. T. Lee, Phys. Rev. A, 41(3), 1721 (1990).
  * [48] K. Husimi, Proc. Phys. -Math. Soc. Japan. 3rd Ser., 22(4), 264–314 (1940).
  * [49] A. Datta et. al., Phys. Rev. A 83, 063836 (2011).
  * [50] J. Zhang et. al., Opt. Exp., 25(21), 24907 (2017).


# Self-Interacting Dark Matter and Small-Scale Gravitational Lenses in Galaxy
Clusters

Daneng Yang  ID<EMAIL_ADDRESS>Department of Physics, Tsinghua
University, Beijing 100084, China  Department of Physics and Astronomy,
University of California, Riverside, California 92521, USA  Hai-Bo Yu  ID
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of California,
Riverside, California 92521, USA

###### Abstract

Recently, Meneghetti et al. reported an excess of small-scale gravitational
lenses in galaxy clusters. We study its implications for self-interacting dark
matter (SIDM), compared with standard cold dark matter (CDM). We design
controlled N-body simulations that incorporate observational constraints. The
presence of early-type galaxies in cluster substructures can deepen
gravitational potential and reduce tidal mass loss. Both scenarios require a
relatively high baryon concentration in the substructure to accommodate the
lensing measurements, and their tangential caustics are similar. The SIDM
substructure can experience gravothermal collapse and produce a steeper
density profile than its CDM counterpart, leading to a larger radial galaxy-
galaxy strong lensing cross section, although this effect is hard to observe.
Our results indicate SIDM can provide a unified explanation to small-scale
lenses in galaxy clusters and stellar motions in dwarf galaxies.

Introduction. Strong gravitational lensing is characterized by the existence
of giant arcs, rings, and multiple images caused by the deflection of lights
by massive foreground galaxies, groups, or galaxy clusters Schneider et al.
(1992); Kneib and Natarajan (2011). It provides a powerful tool for testing
cosmological models Cao et al. (2012); Merten et al. (2015), determining the
mass distribution of clusters Newman et al. (2013); Annunziatella et al.
(2017); Bonamigo et al. (2018), probing substructures Mao and Schneider
(1998); Dalal and Kochanek (2002); Keeton et al. (2003); Xu et al. (2009);
Grillo et al. (2015); Natarajan et al. (2017); Nierenberg et al. (2017);
Gilman et al. (2020a); Minor et al. (2020a) and dark matter properties Vegetti
and Vogelsberger (2014); Minor et al. (2020b); Gilman et al. (2020b); Andrade
et al. (2020); He et al. (2020); Nadler et al. (2021). Recently, Meneghetti et
al. reported that observed substructures in galaxy clusters are more efficient
lenses than those predicted in simulations of standard cold dark matter (CDM)
Meneghetti et al. (2020), indicating that the former are more dense and
compact. Other studies also show strong lensing clusters contain more
substructures with high maximum circular velocities than predicted in CDM
simulations Grillo et al. (2015); Munari et al. (2016). The tension could be
resolved in CDM with a different implementation of baryonic feedback,
resulting in dense substructures Bahé (2021); Robertson (2021).

We study self-interacting dark matter (SIDM) Spergel and Steinhardt (2000);
Kaplinghat et al. (2016); Tulin and Yu (2018) in light of small-scale
gravitational lenses observed in Meneghetti et al. (2020). In this scenario,
dark matter collisions thermalize the inner halo over cosmological timescales
Dave et al. (2001); Vogelsberger et al. (2012); Rocha et al. (2013); Zavala et
al. (2013). Compared with CDM, SIDM better explains diverse dark matter
distributions inferred in a wide range of galactic systems from ultra-diffuse
galaxies to galaxy clusters Firmani et al. (2000); Kaplinghat et al. (2016);
Kamada et al. (2017); Creasey et al. (2017); Salucci (2019); Ren et al.
(2019); Kaplinghat et al. (2019, 2020); Yang et al. (2020); Sagunski et al.
(2021). It may also explain the origin of supermassive black holes Balberg and
Shapiro (2002); Pollack et al. (2015); Choquette et al. (2019); Feng et al.
(2020). The observations of dense substructures in galaxy clusters Meneghetti
et al. (2020) seem to challenge SIDM, as the thermalization could lead to
shallow density cores, which are prone to tidal disruption. In contrast to
this general expectation, we will show in this work that the SIDM substructure
can be as efficient as its CDM counterpart in producing small-scale lenses,
after taking into account realistic stellar distributions and tidal
environments.

We design N-body SIDM and CDM simulations to model the MACS J1206.2-0847
(MACSJ1206) cluster, one of the examples studied in Meneghetti et al. (2020),
and construct four benchmark cases, covering a representative range of halo
concentration. The presence of early-type galaxies in the substructures could
significantly accelerate the onset of SIDM gravothermal collapse Balberg et
al. (2002); Koda and Shapiro (2011); Essig et al. (2019); Huo et al. (2020).
For all benchmarks, our simulated SIDM subhalos experience collapse after
$6~{}{\rm Gyr}$ of tidal evolution and become more dense than their CDM
counterparts, assuming a self-scattering cross section per mass of
$\sigma/m=1~{}{\rm cm^{2}/g}$, which is relatively conservative Tulin and Yu
(2018). For both SIDM and CDM scenarios, the simulated substructures require a
high baryon concentration to be consistent with the observations in Meneghetti
et al. (2020).

We further model strong lensing observables and compute galaxy-galaxy strong
lensing (GGSL) cross sections for the benchmarks. The tangential caustic
predicted in the SIDM substructure is comparable to its CDM counterpart, while
the former has a much larger radial GGSL cross section, with details depending
on the source redshift and initial halo concentration. We will also show mock
lensing images and discuss observational implications of a radial GGSL cross
section.

Modeling the cluster system. We model the host cluster using a static
spherical potential characterized by a Navarro-Frenk-White (NFW) density
profile Navarro et al. (1997) and fix its corresponding scale density and
radius as $\rho_{s}=1.82\times 10^{6}~{}{\rm M_{\odot}/kpc^{3}}$ and
$r_{s}=442~{}{\rm kpc}$, respectively. It well reproduces the projected total
mass profile of MACSJ1206 Bergamini et al. (2019); Caminha et al. (2017);
Bonamigo et al. (2018). We set the apocenter of our simulated substructures to
be $400~{}{\rm kpc}$, with a tangential velocity of $1000~{}{\rm km/s}$.
During the tidal evolution, their distance to the host center oscillates in
the range $\sim 140\textup{--}{370}~{}{\rm kpc}$. The orbit controls the
significance of tidal stripping, and hence the mass loss. In our setup, the
mass of the simulated substructures at $6~{}{\rm Gyr}$ is ${\cal
O}(10^{10})~{}{\rm M_{\odot}}$. The strong lensing analysis in Meneghetti et
al. (2020) focuses on substructures within $15\%$ of the virial radius of the
host cluster, which is $300~{}{\rm kpc}$, and most of them have a mass in the
range $10^{10}\textup{--}10^{11}~{}{\rm M_{\odot}}$. Thus our simulated
cluster system well represents those studied in Meneghetti et al. (2020).

Figure 1:  _Left:_ Evolution of $V_{\rm max}\textup{--}M_{\rm sub}$ for the
simulated SIDM (solid thick) and CDM (solid thin) benchmarks with different
initial subhalo concentrations, denoted using deviations from the cosmological
median at $z=0$, i.e., $+3\sigma$ (magenta), $+2\sigma$ (red), $+1\sigma$
(green) and $+0\sigma$ (blue). The arrows denote the direction of the
evolution and the final snapshot is at $t=6~{}{\rm Gyr}$. CDM simulations
without including the star component are also shown (dotted). For comparison,
the average relation from strong lensing observations (black dashed)
Meneghetti et al. (2020), the range (gray band) and best-fit model relation
(gray dashed) from their CDM simulations are displayed.  _Right:_ SIDM
tangential (dashed) and radial (solid) GGSL cross sections vs source redshift,
normalized to their corresponding CDM counterparts for $t=6~{}{\rm Gyr}$. The
color scheme is the same as the one in the left panel. The simulated CDM
substructures without stars have a low surface density and their lensing
effect is negligible.

For the subhalos, we use an NFW profile to model their initial dark matter
distribution. We fix the initial virial halo mass to be $M_{200}=3\times
10^{12}~{}{\rm M_{\odot}}$, and choose four benchmark values for the
concentration, i.e., $c_{200}=7.49$, $9.65$, $12.4$ and $16.0$, corresponding
to $0\sigma$, $1\sigma$, $2\sigma$ and $3\sigma$ higher than the cosmological
median at $z=0$ Dutton and Macciò (2014), respectively. For each of the
benchmarks, we convert their ($M_{200}$, $c_{200}$) to ($\rho_{s}$, $r_{s}$)
to specify the initial NFW density profile in our simulations. Note for a
given set of $\rho_{s}$ and $r_{s}$, the interpretation of $c_{200}$ depends
on the redshift. Consider $z=2$, at which infall is expected to occur, the
concentration of the benchmarks, from low to high, is $-1.5\sigma$,
$-0.39\sigma$, $+0.73\sigma$, and $+1.8\sigma$ away from the median ($z=2$)
Dutton and Macciò (2014), which are representative.

We fix the initial stellar mass as $M_{\star}=6\times 10^{10}~{}{\rm
M_{\odot}}$, expected from the stellar-to-halo mass relation Moster et al.
(2013), and model its distribution with a truncated singular isothermal
profile as in Meneghetti et al. (2020), $\rho_{\star}(r)=\rho_{0}r^{4}_{\rm
cut}/[r^{2}(r^{2}_{\rm cut}+r^{2})]$, where $\rho_{0}$ is the density
normalization factor and $r_{\rm cut}$ is the cutoff radius. This is
consistent with observations of early-type galaxies Bolton et al. (2006); Treu
et al. (2006); Koopmans et al. (2006); Gavazzi et al. (2007). We take $r_{\rm
cut}=6.23~{}{\rm kpc}$ following the size-mass relation Carleton et al.
(2019), and $\rho_{0}=1.26\times 10^{7}~{}{\rm M_{\odot}/kpc^{3}}$. We use
live particles for the subhalo and stellar components and perform both SIDM
and CDM simulations. For the former, we choose $\sigma/m=1~{}{\rm cm^{2}/g}$,
approximately the lower limit that could explain observations on galactic
scales Tulin and Yu (2018). For comparison, we also perform CDM simulations
without including stars. As we will show that the stellar component is
important in producing strong lensing observables.

We use the public GADGET-2 code Springel (2005); Springel et al. (2001), and
extend it with a module modeling dark matter self-interactions Huo et al.
(2020); Yang et al. (2020), which has been validated in both gravothermal
expansion and collapse regimes with results from Ren et al. (2019); Essig et
al. (2019); Pollack et al. (2015); Koda and Shapiro (2011). We use the code
SpherIC Garrison-Kimmel et al. (2013) to generate initial conditions for the
simulated substructures. The mass of the simulated particle is $10^{6}~{}{\rm
M_{\odot}}$ for both subhalo and stellar components, and the softening length
is $0.2~{}{\rm kpc}$. The resolution is high enough to avoid numerical
artifacts concerning disruption of substructures van den Bosch (2017); van den
Bosch et al. (2018); van den Bosch and Ogiya (2018). We let the simulated
substructures evolve for $6~{}{\rm Gyr}$ in the tidal field of MACSJ1206.

Gravothermal collapse. The left panel of Fig. 1 shows the evolution of the
maximum circular velocity $V_{\rm max}$ vs the total substructure mass $M_{\rm
sub}$ for the benchmarks with $c_{200}=16.0$ (magenta) $12.4$ (red) $9.65$
(green) and $7.49$ (blue) from our SIDM (solid thick) and CDM (solid thin)
simulations, where we include both subhalo and stellar components as in
Meneghetti et al. (2020). The arrow on each curve denotes the direction of the
evolution. The subhalos lose the majority of their mass after tidal evolution,
while the stellar mass is only reduced by an ${\cal O}(1)$ factor. Our final
total stellar and subhalo masses are consistent with those of cluster
substructures from the Illustris simulations Niemiec et al. (2019);
Vogelsberger et al. (2014). For comparison, our CDM simulations without
including stars are shown (dotted).

The maximum circular velocities of the CDM substructures decrease
continuously, aside from oscillatory features due to tidal interactions. For
those with stars, the final $V_{\rm max}$ values are close to the high end of
the range predicted in cosmological hydrodynamical simulations Meneghetti et
al. (2020) (gray band) and the average value from the strong lensing
observations (black dashed). We also see that the $V_{\rm max}$ values
predicted in our CDM simulations without stars are still within the gray band.
It implies that a large population of simulated substructures with AGN
feedback in Meneghetti et al. (2020) has diffuse baryon distributions and high
dark matter fractions. For our simulated substructures with stars, the baryon
mass fraction is $\sim 40\%$ after $6~{}{\rm Gyr}$ of tidal evolution, which
is reasonable compared to that of observed cluster galaxies Annunziatella et
al. (2017).

The simulated SIDM substructures follow a similar trend for most of the
evolution time, but their $V_{\rm max}$ values spike toward the measured ones
at late stages, as gravothermal collapse occurs and their central densities
increase. At $t=6~{}{\rm Gyr}$, all four SIDM benchmarks, even the one with a
median concentration ($z=0$), are denser than their CDM counterparts with
stars. The collapse occurs earlier if $c_{200}$ is higher, leading to a higher
density at $6~{}{\rm Gyr}$, as its timescale is extremely sensitive to the
concentration Essig et al. (2019). We find the presence of stellar particles
deepens potential and accelerates gravothermal evolution, as in the isolated
case Elbert et al. (2018); Sameie et al. (2018); Feng et al. (2020). Without
stars, the collapse would not occur within $6~{}{\rm Gyr}$ unless $c_{200}$ is
$6\sigma$ higher than the median, and a subhalo with median $c_{200}$ would be
nearly destroyed Yang et al. (2020). In the cluster environment, tidal
stripping could also speed up the onset of gravothermal collapse Nishikawa et
al. (2020); Sameie et al. (2020); Kahlhoefer et al. (2019); Correa (2020);
Turner et al. (2020).

Strong lensing observables. To further see implications for strong lensing
observations, we compute GGSL cross sections for the simulated substructures.
We adopt thin-lens approximation and project the mass distribution of the host
cluster and substructure, assumed to be spherical, onto the lens plane, which
is perpendicular to the line of sight. The distance between the substructure
and the host center is fixed to be $300~{}{\rm kpc}$. We denote the angular
positions as $\boldsymbol{\theta}$ and $\boldsymbol{\beta}$ on the lens and
source planes, respectively. It is convenient to introduce effective lensing
potential as Schneider et al. (1992); Kneib and Natarajan (2011):

$\displaystyle\Psi(\boldsymbol{\theta})$ $\displaystyle\equiv$
$\displaystyle\dfrac{1}{\pi}\int
d^{2}\theta\ln|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}|\kappa(\boldsymbol{\theta}^{\prime}),$

where
$\kappa(\boldsymbol{\theta})=\Sigma(D_{L}\boldsymbol{\theta})/\Sigma_{cr}$ is
the scaled projected density and $\Sigma_{cr}={c^{2}D_{S}}/{4\pi
GD_{L}D_{LS}}$ is the critical density. $D_{L}$ and $D_{S}$ are lens and
source angular diameter distances, respectively, and $D_{LS}$ the distance
between the two. We calculate these quantities in a flat universe with matter
energy density $\Omega_{m}=0.3175$ and $h=0.671$ Ade et al. (2014). The lens
equation is
$\boldsymbol{\beta}=\boldsymbol{\theta}-{D_{LS}}\hat{\boldsymbol{\alpha}}/{D_{S}}$,
where $\hat{\boldsymbol{\alpha}}$ is the deflection angle.

Figure 2:  Total (solid) and dark matter (dashed) density profiles for the
SIDM (red) and CDM (gray) benchmarks with $c_{200}=12.4$ ($+2\sigma$) at
$t=6~{}{\rm Gyr}$. The red and gray arrows denote their Einstein radii,
respectively, assuming $z_{s}=3$. The dark matter density profile from CDM
simulations without stars is also shown (dotted); its lensing effect is
negligible. _Insert:_ Evolution of $V_{\rm max}\textup{--}r_{\rm max}$ for the
benchmarks with stars.  Figure 3:  Left: Tangential (dashed) and radial
(solid) caustics for the simulated SIDM and CDM substructures with
$c_{200}=12.4$ ($+2\sigma$), together with four mock sources (circles). The
corresponding critical lines and lensed images are shown in the Middle and
Right panels, respectively, assuming $z_{s}=3$. For the lens redshift
$z_{l}=0.439$ of MACSJ1206, one arcsec corresponds to $5.76~{}{\rm kpc}$.

For each simulated substructure, we model its total density profile with a
numerical interpolation function and construct its surface density assuming
spherical symmetry. We have checked that the substructure slightly deviates
from spherical symmetry after tidal evolution, but the deviation is minor and
it has negligible effects on the lensing observables. For a substructure plus
its host cluster, we determine the lensing potential by solving the Poisson
equation
$\nabla_{\boldsymbol{\theta}}^{2}\Psi(\boldsymbol{\theta})=2\kappa(\boldsymbol{\theta})$.
We implement the fast Fourier transformation method on a $2000\times 2000$
grid, and use the software Mathematica 12 for numerical computation.

After obtaining $\Psi$, we calculate the shear matrix as ${\cal
A}\equiv{\partial\beta_{i}}/{\partial\theta_{j}}=\left(\delta_{ij}-\Psi_{ij}\right)$
and the pseudo-vector shear
$\boldsymbol{\gamma}=\sqrt{[({\Psi_{11}-\Psi_{22}})/{2}]^{2}+\Psi_{12}^{2}}$,
where
$\Psi_{ij}\equiv{\partial^{2}\Psi}/{\partial\theta_{i}\partial\theta_{j}}$ and
$i,j$ are indices of the two spatial coordinates. The tangential and radial
critical lines are contours of $\lambda_{t}=1-\kappa-\gamma=0$ and
$\lambda_{r}=1-\kappa+\gamma=0$, respectively. We obtain their corresponding
caustic lines by mapping them onto the source plane using the lens equation,
and compute tangential and radial GGSL cross sections defined as the area
enclosed by the secondary caustic Meneghetti et al. (2013). We have further
performed convergence tests and confirmed lensing results discussed below are
robust.

The right panel of Fig. 1 shows ratios of SIDM to CDM tangential (dashed) and
radial (solid) GGSL cross sections as a function of the source redshift
$z_{s}$. The tangential cross sections are comparable for both cases and the
differences are within order unity. Interestingly, the SIDM substructure has
larger a radial cross section than its CDM counterpart, by a factor of $\sim
3\textup{--}7$ for $z_{s}\gtrsim 1$, and the difference increases with the
concentration. But this effect is hard to observe as we will discuss later.
For the CDM substructures without stars, their surface density is low and the
lensing effect is negligible.

We see the SIDM substructure can be as compact and dense as its CDM
counterpart in reproducing small-scale lenses of galaxy clusters. In contrast
to the general expectation, the lensing excess reported in Meneghetti et al.
(2020) does not rule out SIDM. Instead, it could be an indicator of
gravothermal catastrophe induced by dark matter self-interactions. We also
find a compact stellar density profile is necessary in reproducing observed
strong lenses in both SIDM and CDM; see also Meneghetti et al. (2020); Bahé
(2021); Robertson (2021) for related discussions.

Although our simulations are based on a particular choice of the cross
section, i.e., $\sigma/m=1~{}{\rm cm^{2}/g}$, the lower limit to be relevant
for dwarf galaxies, the benchmarks capture all the essential features
predicted in SIDM. For example, a larger cross section in the subhalo would
further shorten the collapse timescale, resulting in an even high central
density. Since there is a strong constraint of $\sigma/m\lesssim 0.1~{}{\rm
cm^{2}/g}$ for cluster halos with masses $\sim 10^{15}~{}{\rm M_{\odot}}$
Kaplinghat et al. (2016); Sagunski et al. (2021); Andrade et al. (2020), we
consider a velocity-dependent SIDM scenario, which can be naturally realized
in many particle physics models, see Tulin and Yu (2018).

Density profiles. We take the benchmark with $c_{200}=12.4$ ($+2\sigma$) and
perform a detailed case study. In Fig. 2, we show its dark matter (dashed) and
total (solid) density profiles with (red) and without (gray) dark matter self-
interactions. After $6~{}{\rm Gyr}$ of evolution, the collapse leads to an
overdense region within $1.3~{}{\rm kpc}$, and a less dense region $r\gtrsim
1.3~{}{\rm kpc}$, compared to the CDM subhalo. For the CDM substructure
without including stars (dotted), the density is significantly lower. The
inset displays the evolution of $V_{\rm max}$ and $r_{\rm max}$ for SIDM and
CDM substructures with stars, their $r_{\rm max}$ values decrease overall due
to tidal mass loss. The SIDM one becomes further smaller at late stages, as
the collapse occurs and the central density increases; see also Turner et al.
(2020).

The left panel of Fig. 3 shows the corresponding secondary tangential (dashed)
and radial (solid) caustics, assuming a source at $z_{s}=3$ for the SIDM (red)
and CDM (gray) benchmarks ($c_{200}=12.4$). For SIDM, the enclosed area of the
radial caustic is much larger than the tangential one. For CDM, they are
comparable, but both are smaller than the area of the SIDM radial caustic. For
a collapsed SIDM substructure, the total density profile is steeper than the
isothermal profile $r^{-2}$, and hence its associated radial caustic is larger
than the tangential one.

Mock lensed images. The left panel of Fig. 3 displays four mock sources at
four representative locations; the middle (SIDM) and right (CDM) panels show
their corresponding lensed images, together with the critical lines mapped
from the caustics using the lens equation. The innermost blue source is inside
radial and tangential caustics predicted in SIDM and CDM, and it has four
images in both cases Blandford and Narayan (1986); Schneider et al. (2006).
The source in orange sits on the second fold caustic in SIDM while it only
crosses the tangential caustic in CDM, thus it has one more image in the
former case. Similarly, the sources in green and magenta have one more image
in SIDM than in CDM. Our example explicitly illustrates that the collapsed
SIDM substructure has a higher capability of producing multiple images.
Observationally, the predicted images in the central region are difficult to
detect because they are highly demagnified and obscured by bright objects.
Thus it remains challenging to differentiate a collapsed SIDM substructure
from its CDM counterpart with lensing observations.

We compute the Einstein radii as $r_{E}=\theta_{E}D_{L}\approx 4.3~{}{\rm
kpc}$ and $4.8~{}{\rm kpc}$ for the SIDM and CDM benchmarks with stars,
respectively, denoted in Fig. 2 with arrows, where $r_{E}$ is defined as the
radius of the circle with the same area as that enclosed by the critical line
Meneghetti et al. (2013). It’s not surprising that they are comparable. In our
setup, the stellar component dominates the inner region and both cases have
similar stellar distributions after tidal evolution. In addition, during the
collapse process, the SIDM central density increases, but the total mass does
not change. As indicated in Fig. 2, $r_{\rm max}\sim r_{E}$ for CDM, while
$r_{\rm max}<r_{E}$ for SIDM. Thus the mass distribution induced by
gravothermal collapse does not produce an appreciable change in the enclosed
mass within $r_{E}$.

Discussion and Conclusions. Our simulations assume a compact stellar
distribution motivated by observations of early-type galaxies. Hydrodynamical
simulations show that for field SIDM halos with masses $\gtrsim 10^{12}~{}{\rm
M_{\odot}}$ stars could dominate the inner region Robertson et al. (2018); Di
Cintio et al. (2017); Despali et al. (2019); Robertson et al. (2020); Sameie
et al. (2021). Thus our assumption is well justified. In addition, the SIDM
halo structure is more resilient to feedback than its CDM counterpart, because
of rapid energy redistributions induced by the self-interactions; see Robles
et al. (2017); Fitts et al. (2019). We expect our overall predictions are
robust, and it would be interesting to further test them with cosmological
simulations. In summary, we have shown that a collapsed SIDM substructure
could have a steep density profile, to be consistent with small-scale
gravitational lenses observed in galaxy clusters. The features of strong
lensing observables predicted in the SIDM scenario could be further tested
using existing data Talbot et al. (2020), and upcoming observations Drlica-
Wagner et al. (2019).

## Acknowledgments

We thank Haipeng An, Seong Chan Park, and Massimo Meneghetti for useful
discussions. DY was supported by NSFC under Grant No. 11975134 and the
National Key Research and Development Program of China under Grant
No.2017YFA0402204. HBY was supported by the U.S. Department of Energy under
Grant No. de-sc0008541 and the John Templeton Foundation under Grant ID#
61884. The opinions expressed in this publication are those of the authors and
do not necessarily reflect the views of the John Templeton Foundation.

## References

  * Schneider et al. (1992) P. Schneider, J. Ehlers, and E. E. Falco, _Gravitational Lenses_ (Springer-Verlag Berlin Heidelberg New York, 1992).
  * Kneib and Natarajan (2011) J.-P. Kneib and P. Natarajan, Astron. Astrophys. Rev. 19, 47 (2011), eprint 1202.0185.
  * Cao et al. (2012) S. Cao, Y. Pan, M. Biesiada, W. Godlowski, and Z.-H. Zhu, JCAP 03, 016 (2012), eprint 1105.6226.
  * Merten et al. (2015) J. Merten et al., Astrophys. J. 806, 4 (2015), eprint 1404.1376.
  * Newman et al. (2013) A. B. Newman, T. Treu, R. S. Ellis, D. J. Sand, C. Nipoti, J. Richard, and E. Jullo, Astrophys. J. 765, 24 (2013), eprint 1209.1391.
  * Annunziatella et al. (2017) M. Annunziatella et al., Astrophys. J. 851, 81 (2017), eprint 1711.02109.
  * Bonamigo et al. (2018) M. Bonamigo, C. Grillo, S. Ettori, G. B. Caminha, P. Rosati, A. Mercurio, E. Munari, M. Annunziatella, I. Balestra, and M. Lombardi, Astrophys. J. 864, 98 (2018), eprint 1807.10286.
  * Mao and Schneider (1998) S.-d. Mao and P. Schneider, Mon. Not. Roy. Astron. Soc. 295, 587 (1998), eprint astro-ph/9707187.
  * Dalal and Kochanek (2002) N. Dalal and C. S. Kochanek, Astrophys. J. 572, 25 (2002), eprint astro-ph/0111456.
  * Keeton et al. (2003) C. R. Keeton, B. Gaudi, and A. Petters, Astrophys. J. 598, 138 (2003), eprint astro-ph/0210318.
  * Xu et al. (2009) D. Xu, S. Mao, J. Wang, V. Springel, L. Gao, S. White, C. Frenk, A. Jenkins, G. Li, and J. Navarro, Mon. Not. Roy. Astron. Soc. 398, 1235 (2009), eprint 0903.4559.
  * Grillo et al. (2015) C. Grillo et al., Astrophys. J. 800, 38 (2015), eprint 1407.7866.
  * Natarajan et al. (2017) P. Natarajan et al., Mon. Not. Roy. Astron. Soc. 468, 1962 (2017), eprint 1702.04348.
  * Nierenberg et al. (2017) A. M. Nierenberg, T. Treu, G. Brammer, A. H. G. Peter, C. D. Fassnacht, C. R. Keeton, C. S. Kochanek, K. B. Schmidt, D. Sluse, and S. A. Wright, Mon. Not. Roy. Astron. Soc. 471, 2224 (2017), eprint 1701.05188.
  * Gilman et al. (2020a) D. Gilman, X. Du, A. Benson, S. Birrer, A. Nierenberg, and T. Treu, Mon. Not. Roy. Astron. Soc. 492, L12 (2020a), eprint 1909.02573.
  * Minor et al. (2020a) Q. E. Minor, M. Kaplinghat, T. H. Chan, and E. Simon (2020a), eprint 2011.10629.
  * Vegetti and Vogelsberger (2014) S. Vegetti and M. Vogelsberger, Mon. Not. Roy. Astron. Soc. 442, 3598 (2014), eprint 1406.1170.
  * Minor et al. (2020b) Q. E. Minor, S. Gad-Nasr, M. Kaplinghat, and S. Vegetti (2020b), eprint 2011.10627.
  * Gilman et al. (2020b) D. Gilman, S. Birrer, A. Nierenberg, T. Treu, X. Du, and A. Benson, Mon. Not. Roy. Astron. Soc. 491, 6077 (2020b), eprint 1908.06983.
  * Andrade et al. (2020) K. E. Andrade, J. Fuson, S. Gad-Nasr, D. Kong, Q. Minor, M. G. Roberts, and M. Kaplinghat (2020), eprint 2012.06611.
  * He et al. (2020) Q. He, A. Robertson, J. Nightingale, S. Cole, C. S. Frenk, R. Massey, A. Amvrosiadis, R. Li, X. Cao, and A. Etherington (2020), eprint 2010.13221.
  * Nadler et al. (2021) E. O. Nadler, S. Birrer, D. Gilman, R. H. Wechsler, X. Du, A. Benson, A. M. Nierenberg, and T. Treu (2021), eprint 2101.07810.
  * Meneghetti et al. (2020) M. Meneghetti et al., Science 369, 1347 (2020), eprint 2009.04471.
  * Munari et al. (2016) E. Munari, C. Grillo, G. De Lucia, A. Biviano, M. Annunziatella, S. Borgani, M. Lombardi, A. Mercurio, and P. Rosati, Astrophys. J. Lett. 827, L5 (2016), eprint 1607.01023.
  * Bahé (2021) Y. M. Bahé, Mon. Not. Roy. Astron. Soc. (2021), eprint 2101.12112.
  * Robertson (2021) A. Robertson, Mon. Not. Roy. Astron. Soc. 504, L7 (2021), eprint 2101.12067.
  * Spergel and Steinhardt (2000) D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84, 3760 (2000), eprint astro-ph/9909386.
  * Kaplinghat et al. (2016) M. Kaplinghat, S. Tulin, and H.-B. Yu, Phys. Rev. Lett. 116, 041302 (2016), eprint 1508.03339.
  * Tulin and Yu (2018) S. Tulin and H.-B. Yu, Phys. Rept. 730, 1 (2018), eprint 1705.02358.
  * Dave et al. (2001) R. Dave, D. N. Spergel, P. J. Steinhardt, and B. D. Wandelt, Astrophys. J. 547, 574 (2001), eprint astro-ph/0006218.
  * Vogelsberger et al. (2012) M. Vogelsberger, J. Zavala, and A. Loeb, Mon. Not. Roy. Astron. Soc. 423, 3740 (2012), eprint 1201.5892.
  * Rocha et al. (2013) M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Onorbe, and L. A. Moustakas, Mon. Not. Roy. Astron. Soc. 430, 81 (2013), eprint 1208.3025.
  * Zavala et al. (2013) J. Zavala, M. Vogelsberger, and M. G. Walker, Mon. Not. Roy. Astron. Soc. 431, L20 (2013), eprint 1211.6426.
  * Firmani et al. (2000) C. Firmani, E. D’Onghia, V. Avila-Reese, G. Chincarini, and X. Hernandez, Mon. Not. Roy. Astron. Soc. 315, L29 (2000), eprint astro-ph/0002376.
  * Kamada et al. (2017) A. Kamada, M. Kaplinghat, A. B. Pace, and H.-B. Yu, Phys. Rev. Lett. 119, 111102 (2017), eprint 1611.02716.
  * Creasey et al. (2017) P. Creasey, O. Sameie, L. V. Sales, H.-B. Yu, M. Vogelsberger, and J. Zavala, Mon. Not. Roy. Astron. Soc. 468, 2283 (2017), eprint 1612.03903.
  * Salucci (2019) P. Salucci, Astron. Astrophys. Rev. 27, 2 (2019), eprint 1811.08843.
  * Ren et al. (2019) T. Ren, A. Kwa, M. Kaplinghat, and H.-B. Yu, Phys. Rev. X 9, 031020 (2019), eprint 1808.05695.
  * Kaplinghat et al. (2019) M. Kaplinghat, M. Valli, and H.-B. Yu, Mon. Not. Roy. Astron. Soc. 490, 231 (2019), eprint 1904.04939.
  * Kaplinghat et al. (2020) M. Kaplinghat, T. Ren, and H.-B. Yu, JCAP 06, 027 (2020), eprint 1911.00544.
  * Yang et al. (2020) D. Yang, H.-B. Yu, and H. An, Phys. Rev. Lett. 125, 111105 (2020), eprint 2002.02102.
  * Sagunski et al. (2021) L. Sagunski, S. Gad-Nasr, B. Colquhoun, A. Robertson, and S. Tulin, JCAP 01, 024 (2021), eprint 2006.12515.
  * Balberg and Shapiro (2002) S. Balberg and S. L. Shapiro, Phys. Rev. Lett. 88, 101301 (2002), eprint astro-ph/0111176.
  * Pollack et al. (2015) J. Pollack, D. N. Spergel, and P. J. Steinhardt, Astrophys. J. 804, 131 (2015), eprint 1501.00017.
  * Choquette et al. (2019) J. Choquette, J. M. Cline, and J. M. Cornell, JCAP 07, 036 (2019), eprint 1812.05088.
  * Feng et al. (2020) W.-X. Feng, H.-B. Yu, and Y.-M. Zhong (2020), eprint 2010.15132.
  * Balberg et al. (2002) S. Balberg, S. L. Shapiro, and S. Inagaki, Astrophys. J. 568, 475 (2002), eprint astro-ph/0110561.
  * Koda and Shapiro (2011) J. Koda and P. R. Shapiro, Mon. Not. Roy. Astron. Soc. 415, 1125 (2011), eprint 1101.3097.
  * Essig et al. (2019) R. Essig, S. D. Mcdermott, H.-B. Yu, and Y.-M. Zhong, Phys. Rev. Lett. 123, 121102 (2019), eprint 1809.01144.
  * Huo et al. (2020) R. Huo, H.-B. Yu, and Y.-M. Zhong, JCAP 06, 051 (2020), eprint 1912.06757.
  * Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), eprint astro-ph/9611107.
  * Bergamini et al. (2019) P. Bergamini, P. Rosati, A. Mercurio, C. Grillo, G. B. Caminha, M. Meneghetti, A. Agnello, A. Biviano, F. Calura, C. Giocoli, et al., Astronomy and Astrophysics 631, A130 (2019), eprint 1905.13236.
  * Caminha et al. (2017) G. B. Caminha, C. Grillo, P. Rosati, M. Meneghetti, A. Mercurio, S. Ettori, I. Balestra, A. Biviano, K. Umetsu, E. Vanzella, et al., Astronomy and Astrophysics 607, A93 (2017), eprint 1707.00690.
  * Bonamigo et al. (2018) M. Bonamigo, C. Grillo, S. Ettori, G. B. Caminha, P. Rosati, A. Mercurio, E. Munari, M. Annunziatella, I. Balestra, and M. Lombardi, Astrophysical Journal 864, 98 (2018), eprint 1807.10286.
  * Dutton and Macciò (2014) A. A. Dutton and A. V. Macciò, Mon. Not. Roy. Astron. Soc. 441, 3359 (2014), eprint 1402.7073.
  * Moster et al. (2013) B. P. Moster, T. Naab, and S. D. M. White, Mon. Not. Roy. Astron. Soc. 428, 3121 (2013), eprint 1205.5807.
  * Bolton et al. (2006) A. S. Bolton, S. Burles, L. V. Koopmans, T. Treu, and L. A. Moustakas, Astrophys. J. 638, 703 (2006), eprint astro-ph/0511453.
  * Treu et al. (2006) T. Treu, L. Koopmans, A. Bolton, S. Burles, and L. Moustakas, Astrophys. J. 640, 662 (2006), eprint astro-ph/0512044.
  * Koopmans et al. (2006) L. V. Koopmans, T. Treu, A. S. Bolton, S. Burles, and L. A. Moustakas, Astrophys. J. 649, 599 (2006), eprint astro-ph/0601628.
  * Gavazzi et al. (2007) R. Gavazzi, T. Treu, J. D. Rhodes, L. V. Koopmans, A. S. Bolton, S. Burles, R. Massey, and L. A. Moustakas, Astrophys. J. 667, 176 (2007), eprint astro-ph/0701589.
  * Carleton et al. (2019) T. Carleton, R. Errani, M. Cooper, M. Kaplinghat, J. Peñarrubia, and Y. Guo, Mon. Not. Roy. Astron. Soc. 485, 382 (2019), eprint 1805.06896.
  * Springel (2005) V. Springel, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005), eprint astro-ph/0505010.
  * Springel et al. (2001) V. Springel, N. Yoshida, and S. D. M. White, New Astron. 6, 79 (2001), eprint astro-ph/0003162.
  * Garrison-Kimmel et al. (2013) S. Garrison-Kimmel, M. Rocha, M. Boylan-Kolchin, J. Bullock, and J. Lally, Mon. Not. Roy. Astron. Soc. 433, 3539 (2013), eprint 1301.3137.
  * van den Bosch (2017) F. C. van den Bosch, Mon. Not. Roy. Astron. Soc. 468, 885 (2017), eprint 1611.02657.
  * van den Bosch et al. (2018) F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, Mon. Not. Roy. Astron. Soc. 474, 3043 (2018), eprint 1711.05276.
  * van den Bosch and Ogiya (2018) F. C. van den Bosch and G. Ogiya, Mon. Not. Roy. Astron. Soc. 475, 4066 (2018), eprint 1801.05427.
  * Niemiec et al. (2019) A. Niemiec, E. Jullo, C. Giocoli, M. Limousin, and M. Jauzac, Mon. Not. Roy. Astron. Soc. 487, 653 (2019), eprint 1811.04996.
  * Vogelsberger et al. (2014) M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu, G. F. Snyder, D. Nelson, and L. Hernquist, Mon. Not. Roy. Astron. Soc. 444, 1518 (2014), eprint 1405.2921.
  * Elbert et al. (2018) O. D. Elbert, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, A. S. Graus, and M. Rocha, Astrophys. J. 853, 109 (2018), eprint 1609.08626.
  * Sameie et al. (2018) O. Sameie, P. Creasey, H.-B. Yu, L. V. Sales, M. Vogelsberger, and J. Zavala, Mon. Not. Roy. Astron. Soc. 479, 359 (2018), eprint 1801.09682.
  * Nishikawa et al. (2020) H. Nishikawa, K. K. Boddy, and M. Kaplinghat, Phys. Rev. D 101, 063009 (2020), eprint 1901.00499.
  * Sameie et al. (2020) O. Sameie, H.-B. Yu, L. V. Sales, M. Vogelsberger, and J. Zavala, Phys. Rev. Lett. 124, 141102 (2020), eprint 1904.07872.
  * Kahlhoefer et al. (2019) F. Kahlhoefer, M. Kaplinghat, T. R. Slatyer, and C.-L. Wu, JCAP 12, 010 (2019), eprint 1904.10539.
  * Correa (2020) C. A. Correa (2020), eprint 2007.02958.
  * Turner et al. (2020) H. C. Turner, M. R. Lovell, J. Zavala, and M. Vogelsberger (2020), eprint 2010.02924.
  * Ade et al. (2014) P. A. R. Ade et al. (Planck), Astron. Astrophys. 571, A16 (2014), eprint 1303.5076.
  * Meneghetti et al. (2013) M. Meneghetti, M. Bartelmann, H. Dahle, and M. Limousin, Space Sci. Rev. 177, 31 (2013), eprint 1303.3363.
  * Blandford and Narayan (1986) R. Blandford and R. Narayan, Astrophysical Journal 310, 568 (1986).
  * Schneider et al. (2006) P. Schneider, C. Kochanek, and J. Wambsganss, _Gravitational lensing: strong, weak and micro: Saas-Fee advanced course 33_ , vol. 33 (Springer Science & Business Media, 2006).
  * Robertson et al. (2018) A. Robertson et al., Mon. Not. Roy. Astron. Soc. 476, L20 (2018), eprint 1711.09096.
  * Di Cintio et al. (2017) A. Di Cintio, M. Tremmel, F. Governato, A. Pontzen, J. Zavala, A. Bastidas Fry, A. Brooks, and M. Vogelsberger, Mon. Not. Roy. Astron. Soc. 469, 2845 (2017), eprint 1701.04410.
  * Despali et al. (2019) G. Despali, M. Sparre, S. Vegetti, M. Vogelsberger, J. Zavala, and F. Marinacci, Mon. Not. Roy. Astron. Soc. 484, 4563 (2019), eprint 1811.02569.
  * Robertson et al. (2020) A. Robertson, R. Massey, V. Eke, J. Schaye, and T. Theuns (2020), eprint 2009.07844.
  * Sameie et al. (2021) O. Sameie, M. Boylan-Kolchin, R. Sanderson, D. Vargya, P. Hopkins, A. Wetzel, J. Bullock, and A. Graus (2021), eprint 2102.12480.
  * Robles et al. (2017) V. H. Robles, , et al., Mon. Not. Roy. Astron. Soc. 472, 2945 (2017), eprint 1706.07514.
  * Fitts et al. (2019) A. Fitts et al., Mon. Not. Roy. Astron. Soc. 490, 962 (2019), eprint 1811.11791.
  * Talbot et al. (2020) M. S. Talbot, J. R. Brownstein, K. S. Dawson, J.-P. Kneib, and J. Bautista (2020), eprint 2007.09006.
  * Drlica-Wagner et al. (2019) A. Drlica-Wagner et al. (LSST Dark Matter Group) (2019), eprint 1902.01055.