\section{\label{sec1}Introduction}

Andreev reflection is the fundamental scattering process in the metal-superconductor (SC) heterojunctions, in which an incident electron from the metal is reflected as a hole at the metal-SC interface and a Cooper pair forms in the SC.\cite{Andreev} When bias is less than the superconducting gap, the Andreev reflection dominates the conductance of the metal-SC heterojunction\cite{Blonder,Sun1,ZPNiu}.
Furthermore, the Andreev reflection is also responsible for various of basic physical phenomena such as the Josephson effects\cite{Golubov}, the proximity effects\cite{Buzdin} and the odd frequency pairings\cite{Bergeret,Linder}. For a conventional metal (CM) attached by SC,
the reflected hole moves back along the trajectory of the incident electron. This type of Andreev reflection is called the retro-Andreev reflection (RAR). Recently, with the rise of quantum materials, the so-called specular Andreev reflection (SAR) is first predicted in graphene\cite{Beenakker}. In SAR process, an incident electron will be specularly reflected as the hole.
Up to now, the SAR is discovered and studied extensively in the monolayer graphene\cite{Chengsg,Greenbaum,CXBai,addr1,addr2,Komatsu,Sahu,ChaoWang}, bilayer graphene\cite{Efetov1,Efetov2,Soori}, topological insulator\cite{Majidi}, Weyl semimetal\cite{Chen,Hou,Azizi}, nodal-line semimetal\cite{Cheng}, etc.

The scattering processes of RAR and SAR are distinct, but it is difficult to distinguish them
from the conductance of the metal-SC junction.
For example, the conductance spectra for the three dimensional nodal-line semimetal$-$SC junction possess the same characteristics with those for the graphene-SC junction\cite{Cheng}, although both RAR and SAR are present in the former junction while only RAR or SAR is possible in the latter junction. It has been reported that the shot noise and Fano factor\cite{Zhang} or the Aharonov-Bohm oscillation\cite{Schelter} can be used to characterize the transition between RAR and SAR in the graphene-based SC structures in the presence of ferromagnetic exchange or external field. However, a more universal method of SAR detection irrespective of specific materials is still lack.

The essential difference of RAR and SAR lies in whether the sign relation between the wave vector and the group velocity of the reflected hole is the same with that of the incident electron. We assume the wave vector and the group velocity of the incident electron are of the same (opposite) sign,
then if the two quantities of the reflected hole also have the same (opposite) sign, SAR will happen. On the other hand, if the two quantities of the reflected hole have the opposite (same) sign, RAR will happen. This is because the wave vector component along the junction interface is conserved in the scattering process.

Here, we propose a simple double-junction structure (the metal$-$metal$-$SC junction) to detect SAR without the help of external fields or other interactions, which is completely base on the intrinsic nature of SAR itself. The formation of bound states in the middle metal region of the double junctions is sensitive to the sign relations of the reflected hole. The quasiclassical quantization conditions for energy levels of the bound states are established. The conditions are strongly dependent on the length of the middle metal region for SAR but not for RAR, which lead to the periodic oscillations of conductance when SAR happens. The positions of oscillation peaks and their period can be predicted precisely by the quantization conditions. The establishment of the conditions is irrespective of specific materials and supported by the numerical results. However, to carry out numerical calculations, we construct the double junctions using the recently realized two-dimensional Weyl nodal-line semimetal (WNSM)\cite{Feng,Nie,Niu,Jin} in which the Andreev reflection is of the purely specular type.

The rest of this paper is organized as follows. In Sec. \uppercase\expandafter{\romannumeral2}, the SAR in WNSM$-$SC single junctions with different superconducting pairings and its dependences on the quasiparticle energy, the incident angle and the interfacial barrier are studied. The conductance spectra are presented, which show the same features with those for the RAR. Sec. \uppercase\expandafter{\romannumeral3} builds the quasiclassical quantization conditions for different pairings in the WNSM$-$WNSM$-$SC junctions. The oscillations of SAR and conductance are presented and analysed in details, which have the essential difference from those for the RAR.
Sec. \uppercase\expandafter{\romannumeral4} gives some discussions on our proposed method
and Sec. \uppercase\expandafter{\romannumeral5} concludes this paper.

\section{\label{sec2}SAR and conductance in WNSM$-$SC single junction}
\begin{figure}[!htb]
\centerline{\includegraphics[width=1\columnwidth]{fig1.eps}}
\caption{(a) The energy band of WNSM for electrons. The red arrows denote spin. The product of the wave vector and the group velocity (the orange arrow) is positive for the spin-up electron. (b) The energy band of WNSM for holes. The product of the wave vector and the group velocity (the orange arrow) is also positive for the spin-down hole. (c) The scattering processes in the WNSM$-$SC junction. The Andreev reflection is of the purely specular type. (d) The energy band of CM for electrons. The product of the wave vector and the group velocity (the orange arrow) is positive for the spin-up electron. (e) The energy band of CM for holes. The product of the wave vector and the group velocity (the orange arrow) is negative for the spin-down hole. (f) The scattering processes in the CM$-$SC junction. The Andreev reflection is of the purely retro-type.}
\label{fig1}
\end{figure}
We first study the two-dimensional WNSM$-$SC single junction in the $xy$ plane in order to demonstrate the behaviour of the SAR and its equivalent contribution to conductance as the RAR. The junction consists of the semi-infinite WNSM and SC, as well as an interface located at $x=0$ as schematically shown in Fig. \ref{fig1}. The interfacial barrier is expressed as $U(x)=U_{0}\delta(x)$ by the Dirac delta function $\delta(x)$ and the barrier magnitude $U_{0}$.

\subsection{Hamiltonian and wave functions for WNSM}
The Hamiltonian in the spin space of electrons for WNSM can be written as\cite{Jin}
\begin{eqnarray}
\hat{H}_{WNSM}(k)=\left(\begin{array}{cc}
\frac{\hbar^{2} k^2}{2m}-\mu_{W}&0\\
0&-\frac{\hbar^{2} k^2}{2m}+\mu_{W}\\
\end{array}\right).\label{WH}
\end{eqnarray}
Here $m$ is the effective mass, $\mu_{W}$ the material-dependent parameter characterizing the size of the nodal line and ${\bf{k}}=(k_{x},k_{y})$ the wave vector of electrons. Since the Hamiltonian breaks the time-reversal symmetry, the spin degeneracy is lifted. The dispersion for electrons are plotted in Fig. 1(a).

The Bogoliubov-de Gennes(BdG) Hamiltonian in the particle-hole$\otimes$spin space can be written as
\begin{eqnarray}
\check{H}_{WNSM}=\left(\begin{array}{cc}
\hat{H}_{WNSM}(k)&0\\
0&-\hat{H}^*_{WNSM}(-k)\\
\end{array}\right).
\end{eqnarray}
From the BdG Hamiltonian, it is easy to find that the spin-down (up) hole has the same sign relation of the wave vector and the group velocity with the spin-up (down) electron. For example, both wave vector and group velocity of the left-going (right-going) spin-down hole are negative (positive), and their product is positive regardless of the moving direction, as shown in Fig. \ref{fig1}(b). For the spin-up electron, the product of the wave vector and group velocity is also positive [see Fig.\ref{fig1}(a)], which is the same as that of the spin-down hole.
Due to the sign of the product for the spin-up electron and spin-down hole being the same, this determines that the Andreev reflection in WNSM is the SAR as shown in Fig. \ref{fig1}(c).

The sign relation for a CM is totally different. Consider a CM with the BdG Hamiltonian
\begin{eqnarray}
\check{H}_{CM}(k)=\left(\begin{array}{cc}
\hat{H}_{CM}(k)&0\\
0&-\hat{H}_{CM}^{*}(-k)\\
\end{array}\right),\label{CMH}
\end{eqnarray}
in which $\hat{H}_{CM}(k)=(\frac{\hbar^{2} k^2}{2m}-\mu_{N})1_{2\times2}$ with $1_{2\times2}$ the $2\times2$ identity matrix.
The product of the wave vector and group velocity for the spin-up electron is positive [see Fig.\ref{fig1}(d)], but the product for the spin-down hole is negative [Fig.\ref{fig1}(e)]. This sign difference leads to that the Andreev reflection is the RAR
as shown in Fig. \ref{fig1}(f). In fact, the same (opposite) sign relation for the incident electron
and the reflected hole is the intrinsic character of the SAR (RAR), which is independent of the specific materials\cite{Beenakker,Chengsg}.
In this paper, we will propose a method to detect the SAR and RAR based on the same and opposite sign relations.

Through solving the BdG equation $\check{H}_{WNSM}(-i\frac{\partial}{\partial x},k_{y})\psi(x)=E\psi(x)$ with the substitution of $-i\frac{\partial}{\partial x}$ for $k_{x}$, the wave functions for electrons and holes can be obtained. Due to the translation invariance of the junctions along the $y$ axis, the component $k_{y}$ is conserved in the scattering processes.
When a spin-up electron is injected from WNSM to the interface, both spin-down hole and spin-up electron are specularly reflected. The wave function in WNSM with $x<0$ is solved as
\begin{eqnarray}
\psi_{1}=\left(\begin{array}{c}
1\\
0\end{array}\right)e^{ik^{e\uparrow}_{x}x}
+a_{\downarrow}\left(\begin{array}{c}
0\\
1\end{array}\right)e^{-ik^{h\downarrow}_{x}x}
+b_{\uparrow}\left(\begin{array}{c}
1\\
0\end{array}\right)e^{-ik^{e\uparrow}_{x}x},\label{wfwnsm}
\end{eqnarray}
with the SAR coefficient $a_{\downarrow}$ and the normal reflection coefficient $b_{\uparrow}$. The wave vectors are given by $k_{x}^{e\uparrow}(k_{x}^{h\downarrow})=\sqrt{\frac{2m}{\hbar^2}(\mu_{W}+ E)-k_{y}^2}$ with $k_{y}=\sqrt{\frac{2m}{\hbar^2}(\mu_{W}+E)}\sin\theta$. The angle $\theta$ is the incident angle of the electron as shown in Fig. 1(c). The wave function $\psi_{2}$ for the spin-down electron incidence can be obtained in a similar way. We use $a_{\uparrow}$ and $b_{\downarrow}$ to denote SAR and the normal reflection in this scattering process, respectively.

\subsection{Hamiltonian and wave functions for SC}
For SC, we consider the $s$-wave, $d$-wave and chiral $p$-wave pairings. The BdG Hamiltonian for the SCs is given by
\begin{eqnarray}
H_{SC}=\left(\begin{array}{cccc}
\hat{\epsilon}({\bf{k}})&\hat{\Delta}({\bf{k}})\\
-\hat{\Delta}^{*}(-{\bf{k}})&-\hat{\epsilon}({\bf{-k}}),
\end{array}\right),
\end{eqnarray}
where the single-particle energy $\hat{\epsilon}({\bf{k}})=(\frac{\hbar^2 k^2}{2m}-\mu_{S})\hat{1}_{2\times2}$, the energy matrix $\hat{\Delta}({\bf{k}})=\Delta_{0}f({\bf{k}})i\hat{\sigma}_{y}$ for the $s$-wave or $d$-wave SC and $\hat{\Delta}({\bf{k}})=\Delta_{0}f({\bf{k}})\hat{\sigma}_{x}$ for the $p$-wave SC. Here, $\hat{\sigma}_{x}$ and $\hat{\sigma}_{y}$ are Pauli matrices in the spin space. For the $s$, $d_{x^2-y^2}$ and $d_{xy}$-wave pairings\cite{KashiwayaRPP,Hirai}, $f({\bf{k}})=1$, $\hat{k}_{x}^2-\hat{k}_{y}^2$ and $2\hat{k}_{x}\hat{k}_{y}$, respectively, while for the $p$-wave pairing\cite{Hirai,Mackenzie}, $f({\bf{k}})=\hat{k}_{x}+i\hat{k}_{y}$. The effective mass for SCs has been taken as the same with that for WNSM and $\mu_{S}$ is the chemical potential.

Through solving the BdG equation, we can obtain the wave functions in the SC region with $x>0$. For the spin-up electron incidence, the function is written as
\begin{eqnarray}
\Psi_{1}=c_{\uparrow}\left(\begin{array}{c}
u_{+}e^{i\phi_{+}}\\
v_{+}\end{array}\right)e^{ik_{x}x}
+d_{\downarrow}\left(\begin{array}{c}
v_{-}e^{i\phi_{-}}\\
u_{-}\end{array}\right)e^{-ik_{x}x}.
\end{eqnarray}
Here, $c_{\uparrow}$ and $d_{\downarrow}$ represent the transmissions as electron-like quasiparticle and hole-like quasiparticle, respectively. The wave vector is $k_{x}=\sqrt{\frac{2m}{\hbar^2}\mu_{S}-k_{y}^2}$
under the Andreev approximation\cite{Andreev}. For the $s$-wave and $p$-wave SCs, $u_{+/-}=\sqrt{\frac{E+\Omega}{2E}}$ and $v_{+/-}=\sqrt{\frac{E-\Omega}{2E}}$ with $\Omega=\sqrt{E^2-\Delta_{0}^2}$. For the $d$-wave SC, $u_{\pm}=\sqrt{\frac{E+\Omega_{\pm}}{2E}}$ and $v_{\pm}=\sqrt{\frac{E-\Omega_{\pm}}{2E}}$ with $\Omega_{\pm}=\sqrt{E^2-\Delta_{\pm}^2}$. The direction dependent energy gap $\Delta_{\pm}=\Delta_{0}\cos{[2(\theta_{s}\mp\beta)]}$ in which $\theta_{s}$ is the transmission angle for electron-like quasiparticles and $\beta$ is the angle between the interface normal and the crystallographic $a$-axis of the $d$-wave SC. The angle $\theta_{s}$ can be expressed by $\theta$ under the conservation of $k_{y}$. For the $d_{x^2-y^2}$-wave SC, $\beta=0$ and for the $d_{xy}$-wave SC, $\beta=\frac{\pi}{4}$. The internal phase factors $e^{i\phi_{\pm}}$ are $1$, $\frac{\cos(2\theta_{s}\mp2\beta)}{|\cos(2\theta_{s}\mp2\beta)|}$ and $\frac{\pm\cos{\theta_{s}}+ i\sin{\theta_s}}{|\cos{\theta_s}+i\sin{\theta_s}|}$ for the $s$-wave pairing, the $d$-wave pairing and the $p$-wave pairing, repsectively. For the spin-down electron incidence, the wave function $\Psi_2$ for SC can be derived in a similar way.

\subsection{Boundary conditions and conductance}
The boundary conditions which ensure the probability conservation for the WNSM-SC single junction are given by
\begin{eqnarray}
\begin{split}
&\psi_{1(2)}|_{x=0^-}=\Psi_{1(2)}|_{x=0^+},\\
&\Psi'_{1}|_{x=0^+}-\hat{\tau}_{z}\psi'_{1}|_{x=0^-}=\frac{2mU}{\hbar^2}\psi_{1}|_{x=0},\\
&\Psi'_{2}|_{x=0^+}+\hat{\tau}_{z}\psi'_{2}|_{x=0^-}=\frac{2mU}{\hbar^2}\psi_{2}|_{x=0},\label{WSBC}
\end{split}
\end{eqnarray}
with $\hat{\tau}_{z}$ being the Pauli matrix in the particle-hole space. Under the boundary conditions, the reflection coefficients $a_{\uparrow}$, $a_{\downarrow}$, $b_{\uparrow}$ and $b_{\downarrow}$ can be solved (see Appendix).

We define the probabilities for SAR and the normal reflection as
\begin{eqnarray}
A_{\uparrow}=\text{Re}\left[\frac{k_{x}^{h\uparrow}}{k_{x}^{e\downarrow}}\right]|a_{\uparrow}|^2,~~~B_{\uparrow}=|b_{\uparrow}|^2,\\
A_{\downarrow}=\text{Re}\left[\frac{k_{x}^{h\downarrow}}{k_{x}^{e\uparrow}}\right]|a_{\downarrow}|^2,~~~B_{\downarrow}=|b_{\downarrow}|^2.
\label{prob}
\end{eqnarray}
According to the Blonder-Tinkham-Klapwijk theory\cite{Blonder}, the conductance can be expressed as
\begin{eqnarray}
\sigma_{\uparrow}=1+A_{\downarrow}-B_{\uparrow},\\
\sigma_{\downarrow}=1+A_{\uparrow}-B_{\downarrow},\label{fencond}
\end{eqnarray}
which are caused by the spin-up electron incidence and the spin-down electron incidence, respectively.
The total conductance normalized by the normal conductance is given by
\begin{eqnarray}
\sigma(eV)=\frac{\int_{-\pi/2}^{\pi/2}(\sigma_{\uparrow}+\sigma_{\downarrow})\cos{\theta}d\theta}
{\int_{-\pi/2}^{\pi/2}(\sigma_{n\uparrow}+\sigma_{n\downarrow})\cos{\theta}d\theta},\label{totalcond}
\end{eqnarray}
where $\sigma_{n\uparrow}$ and $\sigma_{n\downarrow}$ are the conductance when SCs are in the normal state and $V$ is the bias of the junction. Note, the zero temperature conductance is considered in this paper. In this situation, we have the relations $E=eV$ and $\sigma(E)=\sigma(eV)$ in accordance with the Blonder-Tinkham-Klapwijk theory\cite{Blonder}.

\subsection{\label{sec3}Numerical results}
For simplicity, we set $\mu_{W}=\mu_S=\mu$ in our calculations. The mismatch between $\mu_{W}$ and $\mu_{S}$ will suppress the conductance but will not fundamentally change our physical results. Since the conductance are independent on the spin of the incident electron, we only show the numerical results of probability $A_{\downarrow}$ for the spin-up electron incidence. We define the effective interfacial barrier as $z=\frac{2mU_{0}}{\hbar^2k_{F}}$ with $k_{F}=\sqrt{\frac{2m\mu}{\hbar^2}}$.

\begin{figure}[!htb]
\centerline{\includegraphics[width=0.9\columnwidth]{SAR.eps}}
\caption{The incident angle $\theta$ and the incident energy $E$ dependences of the SAR probability at $z=1$ for (a) the $s$-wave pairing, (b) the $d_{x^2-y^2}$-wave pairing, (c) the $d_{xy}$-wave pairing and (d) the $p$-wave pairing.}
\label{SAR}
\end{figure}

Now, we discuss the incident angle $\theta$ and the quasiparticle energy $E$ dependences of SAR. For the transparent junction with $z=0$, the subgap SAR probability $A_{\downarrow}$ for the $s$-wave and the $p$-wave SCs is always $1$, which is irrespective of $\theta$ and $E$. This will lead to the normalized conductance of the value $2$ in the whole gap. However, for the $d$-wave SCs, the $\theta$-independent SAR with $A_{\downarrow}=1$ only happens at $E=0$. As a result, the zero-bias conductance (ZBC) is still $2$ but the conductance for $eV>0$ is suppressed.

For the nonzero interfacial barrier with $z\ne0$, we present the SAR probabilities for $z=1$ in Fig. {\ref{SAR}}. For the $s$-wave pairing in Fig. \ref{SAR}(a), it is found that the probability $A_{\downarrow}$ is dramatically weakened for the small energy near $E=0$ compared with the SAR of $z=0$. When the incident energy $E$ is increased from 0 to $\Delta_{0}$, the SAR probability tends to $1$. This means the conductance increases with the bias in the gap. The similar thing happens for the $d_{x^2-y^2}$-wave pairing in Fig. {\ref{SAR}}(b). The conductance spectra at $z=0$ and $z=1$ for the $s$-wave and the $d_{x^2-y^2}$-wave pairings are shown in Figs. {\ref{cond}}(a) and (b). The behaviours of the conductance are consistent with the above analyses on SAR.
The situations for the $d_{xy}$-wave and the $p$-wave pairings become very different. The SARs are reduced near $E=\Delta_{0}$ when one increases $z$ from $0$ to $1$, but they keep the value $A_{\downarrow}\sim 1$ in a large range of the incident angle near $E=0$ as shown in Figs. \ref{SAR}(c) and (d). The conductance will decrease as the bias is increased in the gap. The conductance spectra for $z=0$ and $z=1$ are shown in Figs. \ref{cond}(c) and (d).

If the interfacial barrier continues to increase, for example to $z=3$, the features of SAR presented in Fig. \ref{SAR} will become more obvious. The resulting conductance spectra are plotted in Figs. \ref{cond}(a)-(d) for different pairings. For the $s$-wave and the $d_{x^2-y^2}$-wave pairings, the shapes of conductance evolve towards the bulk density of states of SCs\cite{Blonder,Tanaka}. For the $d_{xy}$-wave and the $p$-wave pairings, the conductance possesses the ZBC peaks, which characterize the existence of the bound states at the surface of SCs\cite{Kashiwaya}. It can be concluded that the SAR can also well reflect the anisotropic properties of SCs as RAR.

\begin{figure}[!htb]
\centerline{\includegraphics[width=0.9\columnwidth]{cond.eps}}
\caption{The normalized conductance spectra at different interfacial barriers for (a) the $s$-wave pairing, (b) the $d_{x^2-y^2}$-wave pairing, (c) the $d_{xy}$-wave pairing and (d) the $p$-wave pairing.}
\label{cond}
\end{figure}

The SAR in the WNSM$-$SC single junction possesses the same features with the RAR in the CM$-$SC single junction. Beneficially, we can obtain the same amplitude of SAR in WNSM as RAR in CM. Unfortunately, as an observable quantity, the conductance for the WNSM$-$SC junction is of the same form with that for the CM$-$SC junction\cite{Blonder,Tanaka,Yamashiro}. This indicates the contributions of the SAR and RAR to the conductance are equivalent. It is difficult to distinguish the two types of reflections by measuring the electric transport of the single-junction structure.

\section{\label{sec4}Detection of SAR in WNSM$-$WNSM$-$SC junctions}

In order to distinguish SAR and RAR, we consider the WNSM$-$WNSM$-$SC double-junction structure as shown in Fig. \ref{wws}. The two interfaces are located at $x=0$ and $x=L$, which potential can be described by $U(x)=U_{1}\delta(x)+U_{2}\delta(x-L)$. When an electron is injected from the left WNSM, it will transmit into the middle WNSM through the interface at $x=0$. The transmitted electron impings on the interface at $x=L$ and is specularly reflected as the left-going hole. The left-going hole is normally reflected as the right-going one at the interface with $x=0$ and then specularly reflected as the left-going electron at the interface with $x=L$. The left-going electron will be normally reflected as the right-going one at the interface $x=0$. After twice SARs and twice normal reflections, a closed path of the quasiparticle motion is formed as depicted in Fig. \ref{wws}(a).

\begin{figure}[!htb]
\centerline{\includegraphics[width=0.7\columnwidth]{wws.eps}}
\caption{(a) Schematic illustration of the WNSM-WNSM-SC double junctions and the scattering processes in the middle WNSM for the electron incidence from the left WNSM. The group velocities are denoted by blue lines which form a closed path. Along the path, the wave vectors are denoted by red lines, correspondingly. (b) The same closed path (blue lines) and the corresponding wave vectors (red lines) along the path in the CM-CM-SC double junctions. In the figure, the solid lines for electron and the dashed ones for hole.
Here the directions of the wave vectors for holes in (a) and (b) are
opposite.}
\label{wws}
\end{figure}
\subsection{$z_{2}=0$}
We first consider the junctions with $z_{2}=0$ to clarify the basic physics.
For the $s$-wave and $d_{x^2-y^2}$-wave SCs, if the phase accumulated along the closed path satisfies the following quasiclassical quantization condition,
\begin{eqnarray}
-\sum_{\alpha=\pm}\arccos{\frac{E}{|\Delta_{\alpha}|}}+2(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})L=2n\pi,\label{sdqqc}
\end{eqnarray}
with $n$ an integer number, the energy levels of bound states will be formed. For the $s$-wave pairing, $\Delta_{\pm}=\Delta_{0}$ while for the $d_{x^2-y^2}$-wave pairing, $\Delta_{\pm}=\Delta_{0}\cos(2\theta\mp2\beta)$ with $\beta=0$. The wave vectors $k_{x}^{e\uparrow}$ and $k_{x}^{h\downarrow}$ are the ones in Eq. (\ref{wfwnsm}). The first term in Eq. (\ref{sdqqc}) is acquired from the Andreev reflection (see Appendix for details). The second term originates from the motion of quasiparticles along the closed path. For the zero energy bound state with $E=0$ and the normal incidence with $k_{y}=0$, the quantization condition degenerates into
\begin{eqnarray}
k_{F}L=\frac{(2n+1)\pi}{4},\label{p1}
\end{eqnarray}
with $k_{F}=\sqrt{\frac{2m\mu}{\hbar^2}}$. This indicates ZBC of the junctions will oscillate with the length $L$ of the middle WNSM. The peaks of the ZBC appear at $k_{F}L=\frac{\pi}{4}, \frac{3\pi}{4}, \frac{5\pi}{4}\cdots.$ The oscillation period defined as the spacing between the neighbouring peaks is $\Delta k_{F}L=\frac{\pi}{2}$.

For the $d_{xy}$-wave and the $p$-wave SCs, the quasiclassical quantization condition becomes
\begin{eqnarray}
\pi-\left[\sum_{\alpha=\pm}\arccos{\frac{E}{|\Delta_{\alpha}|}}+2\delta\theta\right]
+2(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})L=2n\pi,\label{pdqqc}
\end{eqnarray}
where $\Delta_{\pm}=\Delta_{0}$ and $\delta=1$ for the $p$-wave SC while $\Delta_{\pm}=\Delta_{0}\cos[2(\theta\mp\beta)]$ with $\beta=\frac{\pi}{4}$ and $\delta=0$ for the $d_{xy}$-wave SC. The difference between Eq. (\ref{sdqqc}) and Eq. (\ref{pdqqc}) is caused by the anisotropy of the $p$-wave and the $d$-wave SCs (see Appendix for details). The phase $(-\arccos{\frac{E}{|\Delta_{+}|}}-\delta\theta)$ is acquired by the SAR with the conversion from electron to hole while the phase $(\pi-\arccos{\frac{E}{|\Delta_{-}|}}-\delta\theta)$ is acquired by the SAR with the conversion from hole to electron as shown in Fig. {\ref{wws}}(a).
For the $s$-wave and the $d_{x^2-y^2}$-wave SCs, the two phases are the same. According to Eq. (\ref{pdqqc}), the zero energy levels for the normal incidence ($\theta=0$) can be formed when
\begin{eqnarray}
k_{F}L=\frac{2n\pi}{4},\label{p2}
\end{eqnarray}
which indicates the peaks of the ZBC appear at $k_{F}L=0,\frac{\pi}{2},\pi\cdots$. The oscillation period is $\Delta k_{F}L=\frac{\pi}{2}$.

The expressions in Eqs. (\ref{sdqqc}) and (\ref{pdqqc}) are our critical results which can be used to distinguish SAR from RAR. The key point is the presence of the term $2(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})L$ derived from the integration of wave vectors along the closed path. This term embodies the intrinsic character of the specularly Andreev reflected hole. For the spin-up  electron in the middle WNSM, its wave vectors and its group velocities have the same sign, i.e., they are in the same direction, see Figs. \ref{fig1}(a) and \ref{wws}(a).
The specularly reflected spin-down hole also possesses the relation. Its wave vectors and group velocities are also in the same direction as shown in
Figs. \ref{fig1}(b) and \ref{wws}(a).
As a result, along the closed path, the phase acquired by the motion of the hole is $2k_{x}^{h\downarrow}L$ when the phase for the electron is $2k_{x}^{e\uparrow}L$.

However, it is not the case for RAR in CM-CM-SC double junctions since the retro-Andreev reflected hole has the opposite sign relation with the incident electron as shown in Figs. \ref{fig1}(d), \ref{fig1}(e) and \ref{wws}(b).
For the spin-up electron in the middle CM, its wave vectors and its group velocities have the same sign, but the wave vectors and group velocities of the retro-reflected hole have the opposite sign\cite{Zagoskin,Chengphyc,ChengJLTP}.
Along the closed path, the phase acquired by the motion of the hole is $-2k_{x}^{h\downarrow}L$ when the phase for the electron is $2k_{x}^{e\uparrow}L$. For $E=0$, we have $k_{x}^{e\uparrow}=k_{x}^{h\downarrow}$ exactly. The phase from the quasiparticle motion is zero. The quantization condition is independent on the length of CM. There is no oscillation effect for ZBC. Even if $E\ne0$, the conclusion still holds as long as $E$ is much smaller than the chemical potential.

\begin{figure}[!htb]
\centerline{\includegraphics[width=0.9\columnwidth]{osi.eps}}
\caption{The length $k_{F}L$ and the incident angle $\theta$ dependences of SAR at $E=0$ for (a) the $s$-wave pairing and (b) the $d_{xy}$-wave pairing. (c) The oscillation of ZBC as the length $k_{F}L$ for the $s$-wave pairing and the $d_{xy}$-wave pairing in the WNSM-WNSM-SC junctions. The effective interfacial barriers have been taken as $z_{1}=5$ and $z_{2}=0$. The width-independent straight dashed lines denote the ZBCs for the $s$-wave pairing (the lower one) and the $d_{xy}$-wave pairing (the upper one) in the CM-CM-SC junctions under the same parameters.}
\label{osi}
\end{figure}

To support our physical analyses, we calculate numerically SAR and conductance in WNSM-WNSM-SC double junctions. We define the effective interfacial barriers as $z_{1}=\frac{2mU_{1}}{\hbar^2k_{F}}$ and $z_{2}=\frac{2mU_{2}}{\hbar^2k_{F}}$. In the calculations, $z_{1}=5$ and $z_{2}=0$ are taken. Fig. {\ref{osi}}(a) and (b) show the length and the incident angle dependences of SAR at $E=0$ for the $s$-wave pairing and the $d_{xy}$-wave pairing, respectively. For $k_{F}L \rightarrow 0$, SAR is almost completely suppressed for the $s$-wave pairing while it can happen for the $d_{xy}$-pairing. The former corresponds to the ZBC near to zero and the latter corresponds to the ZBC peak. They are consistent with the results of the WNSM-SC single junction as given in Figs. \ref{cond}(a) and (c).

The oscillation of the SAR probabilities with the length of the middle WNSM is very obvious for a given incident angle $\theta$, which can be well explained from the quantization conditions Eqs. (\ref{sdqqc}) and (\ref{pdqqc}). The bright fringes demonstrate the formation of the zero energy bound states. Each bright fringe corresponds to an integer number $n$. From the bottom to the top, $n$ is $0,1,2\cdots$. For a certain bright fringe with the number $n$, it bends in the direction of the length increase as the incident angle is raised. This is because the wave vectors $k_{x}^{e\uparrow}$ and $k_{x}^{h\downarrow}$ are reduced when one raises the incident angle $\theta$. According to Eqs. (\ref{sdqqc}) and (\ref{pdqqc}), we need to enlarge the length $L$ to ensure the establishment of the equations with the fixed number $n$.

Fig. \ref{osi}(c) presents the conductance integrated about $\theta$ for the $s$-wave pairing and the $d_{xy}$-wave pairing. Although the integration about $\theta$ will bring about broadening of the ZBC peaks, the positions of peaks still fit Eqs. (\ref{p1}) and (\ref{p2}) very well. The $\frac{\pi}{4}$ shift of the peaks for the $d_{xy}$-pairing relative to those for the $s$-wave pairing is derived from the extra $\pi$ phase in Eq. (\ref{pdqqc}) due to the anisotropy of the $d$-wave SC. This shift also shows itself in SAR as given in Fig. \ref{osi}(a) and (b). The SAR and ZBC for the $d_{x^2-y^2}$-wave and the $p$-wave pairings are not presented here since the former is similar to the $s$-wave case and the latter is similar to the $d_{xy}$-wave case. In addition, the measurement on the period of the discrete ZBC peaks can provide the immediate information about the nodal-line size in WNSM.

For comparision, we also present the ZBCs of the CM-CM-SC junctions for the $s$-wave pairing (the lower straight dashed line) and the $d_{xy}$-wave pairing (the upper straight dashed line) in Fig. \ref{osi}(c). They are both independent on the width of the middle CM as we have discussed. Actually, the width-independent RAR for the zero energy in the CM-CM-SC junctions are heavily suppressed for the $s$-wave pairing by the high barrier with $Z_{1}=5$, while the probability of RAR is still $1$ for the $d_{xy}$-wave pairing. Accordingly, the ZBC for the $s$-wave pairing is almost zero and that for the $d_{xy}$-pairing has a peak with a large value. In Fig. \ref{osi}(c), the value for the lower dashed line has been magnified eighty times when the value for the upper dashed line has been reduced five times. These behaviors of RAR and ZBC for the CM-CM-SC double junctions are the same as those for the CM-SC single junction but totally different from the features for the WNSM-WNSM-SC double junctions.

\subsection{$z_{2}\ne0$}
\begin{figure}[!htb]
\centerline{\includegraphics[width=0.9\columnwidth]{osi2.eps}}
\caption{The length $k_{F}L$ and the incident angle $\theta$ dependences of SAR at $E=0$ for (a) the $s$-wave pairing and (b) the $d_{xy}$-wave pairing. (c) The oscillation of ZBC as the length $k_{F}L$ for the $s$-wave pairing and the $d_{xy}$-wave pairing. The effective interfacial barriers have been taken as $z_{1}=5$ and $z_{2}=2$.}
\label{osi2}
\end{figure}
Finally, we consider the more realistic junctions with $z_{2}\ne0$. Fig. \ref{osi2} shows SAR probabilities and ZBC for $z_{1}=5$ and $z_{2}=2$ at $E=0$. For the $d_{xy}$-wave pairing, both SAR and ZBC oscillations remain virtually unchanged, which include the positions of peaks and the period. However, there are significant changes for the $s$-wave pairing. First, the ZBC peaks are greatly weakened, which have been magnified ten times in Fig. \ref{osi2}(c). Second, the distribution of peaks is no longer uniform. The spacing between the $2n$th and the $(2n+1)$th peaks with $n\geq1$ becomes small. This can also be seen from the SAR probability in Fig. \ref{osi2}(a). Actually, the spacing will get smaller and smaller as $z_{2}$ increases. Two neighbor peaks will form a main peak with split when $z_{2}$ is large enough, for example $z_{2}=5$. The two changes result from the suppression of SAR and the enhancement of normal reflection of the $s$-wave pairing when the effective barrier $z_{2}$ is raised at $E=0$.

The normal reflection in the middle WNSM also leads to conductance oscillation. Twice normal reflections, one at $x=0$ and the other at $x=L$, will form a closed path. The quasiclassical quantization condition is simply $k_{F}L=n\pi$ for the normal incidence with $\theta=0$. Hence, the positions of main peaks are located at $k_{F}L=0,\pi,2\pi,\cdots$. Each main peak possesses two split peaks from SAR. The similar thing also happens for the $d_{x^2-y^2}$ wave and the $p$-wave pairings. It is just that the peaks for the $p$-wave pairing are not weakened so much since there is still a range of incident angle for strong SAR at $E=0$ as given in Fig. \ref{SAR}(d). The situation is different for the $d_{xy}$-wave pairing which ZBC peaks from SAR can survive even if $z_{2}=5$. The total SAR at $E=0$ is unaffected by the interfacial barrier as shown in Fig. \ref{SAR}(c). In brief, the SAR induced oscillation is completely distinguishable for all pairings at least when $z_{2}\leq2$. For more large $z_{2}$, the split of main peaks is still a strong signal for SAR.

\section{\label{secdiscussion} Discussions}

Let us give short discussions on the detection of SAR in other materials.
For WNSM here, we study the behaviors of ZBC, which are related to the formation of zero energy Andreev levels. In this situation, the conductance oscillates for SAR while it does not for RAR.
However, the zero energy is not the necessary requirement to distinguish the two Andreev reflections.
Even if the finite bias conductances with $E\ne0$ are considered,
which will oscillate both for RAR and SAR, their oscillation behaviors are still distinct.
The key is that the quantization conditions for the levels of Andreev bound states are different for RAR and SAR due to their intrinsic characters.
The conditions will form different Andreev levels for RAR and SAR in the junction parameter space,
which can show themselves in conductances and be detected.
For materials such as the monolayer graphene\cite{Beenakker}, bilayer graphene\cite{Efetov2}, topological insulator\cite{Majidi} and nodal-line semimetal\cite{Cheng}, the establishment of quantization conditions could be a little more complicated, but our method is still effective,
which can still predict some distinguishable features for RAR and SAR.

Finally, we discuss the difference between the CM$-$CM$-$SC (WNSM$-$WNSM$-$SC) junctions
and the SC$-$CM$-$SC (SC$-$WNSM$-$SC) junctions
to show that the proposed method to distinguish the RAR and SAR
can only work in the WNSM$-$WNSM$-$SC junctions.
Actually, for clean SC$-$CM(WNSM)$-$SC junctions with transparent interfaces,
the quantization condition for the Andreev levels can be simply written as
(here we take the $s$-wave SC as an example)\cite{Zagoskin}
\begin{eqnarray}
-2\arccos{\frac{E_{n}^{\pm}}{\Delta_{0}}}\pm\varphi
+(k_{x}^{e\uparrow}-k_{x}^{h\downarrow})L=2n\pi ,\label{R1}
\end{eqnarray}
for the RAR with the superconducting phase difference $\varphi$,
and the quantization condition is
\begin{eqnarray}
-2\arccos{\frac{E_{n}^{\pm}}{\Delta_{0}}}\pm\varphi
+(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})L=2n\pi ,\label{R2}
\end{eqnarray}
for the SAR.
Although the quantization conditions in Eqs.(\ref{R1}) adn (\ref{R2}) are similar with
that in Eq.(\ref{sdqqc}),
this dose not mean that the Josephson current from the RAR in the SC$-$CM(WNSM)$-$SC junctions
is not dependent on the length $L$ because the Josephson current includes contributions from all Andreev bound states with $E_{n}^{\pm}<0$
even at zero temperature.
From Eq.(\ref{R1}), the formation of the zero energy level $E=0$ is irrespective of the length $L$,
but the levels with $E\ne 0$ is dependent on the length $L$, leading to the Josephson current
depending on the length $L$ also.

\section{\label{sec5}Conclusions}
We propose a more universal method to detect SAR by using
the semimetal-semimetal-SC double-junction structure.
The method is completely base on the intrinsic character of
the specularly reflected hole in SAR which has the same sign relation
about the group velocity and the wave vector with the incident electron.
However, the sign relation for the retro-Andreev reflected hole in RAR
is opposite. This essential difference between SAR and RAR
can be seized subtly by the accumulated phase of the quasipartilce
motion along a close path in the middle semimetal.
Through establishing the quasiclassical quantization conditions
for the energy levels of bound states, one finds the ZBC from SAR
periodically oscillates with the length of the middle region
while that from RAR does not.
Further, the positions and the period of ZBC peaks can be predicted
precisely with the quantization conditions.
This provides strong distinguishable signatures of SAR.
Actually, even though the phases acquired by electron and hole
do not cancel for RAR at the non-zero incident energy,
the oscillation features of ZBC, including positions and period of peaks,
are still recognizably different from those for SAR.

\section*{\label{sec5}ACKNOWLEDGMENTS}
This work was financially supported by National Key R and D Program of China (2017YFA0303301),
NSF-China under Grants Nos. 11921005 and 11447175,
the Strategic Priority Research Program of Chinese Academy of Sciences (XDB28000000)
and the Natural Science Foundation of Shandong Province under Grants No. ZR2017QA009.

\section{Appendix}
\setcounter{equation}{0}
\setcounter{subsection}{0}
\renewcommand{\theequation}{A.\arabic{equation}}
\renewcommand{\thesubsection}{A.\arabic{subsection}}
\subsection{Derivation of the quasiclassical quantization conditions}
Using the boundary conditions in Eq.(\ref{WSBC}), the SAR and normal reflecton coefficients for the spin-up electron incidence in WNSM-SC junction are solved as
\begin{eqnarray}
\begin{split}
a_{\downarrow}=\frac{4k_{x}k_{x}^{e\uparrow}u_{-}v_{+}}{k_{x}^2\xi_{-}+k_{x}(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})\xi_{+}
+(z-i k_{x}^{e\uparrow})(z+i k_{x}^{h\downarrow})\xi_{-}},\\
b_{\uparrow}=\frac{k_{x}(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})\xi_{+}
-k_{x}^2\xi_{-}+( k_{x}^{e\uparrow}-iz)(k_{x}^{h\downarrow}-iz)\xi_{-}}{k_{x}^2\xi_{-}+k_{x}(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})\xi_{+}
+(z-i k_{x}^{e\uparrow})(z+i k_{x}^{h\downarrow})\xi_{-}},
\end{split}
\end{eqnarray}
with $\xi_{\pm}=u_{+}u_{-}e^{i\phi_{+}}\pm v_{+}v_{-}e^{i\phi_{-}}$.

We consider the transparent junction with $z=0$.
Under the Andreev approximation, we take $k_{x}=k_{x}^{e\uparrow}=k_{x}^{h\downarrow}$. The SAR coefficient degenerates into
\begin{eqnarray}
a_{\downarrow}=\frac{v_{+}}{u_{+}}e^{-i\phi_{+}}.
\end{eqnarray}
For the $s$-wave and the $d_{x^2-y^2}$-wave pairings, $a_{\downarrow}=e^{-i\arccos{\frac{E}{|\Delta_{+}|}}}$ with $E<|\Delta_{+}|$. The expression implies a phase of $(-\arccos{\frac{E}{|\Delta_{+}|}})$ is acquired by the specularly reflected hole in the SAR process.
For the $d_{xy}$-wave and the $p$-wave pairings, the coefficient $a_{\downarrow}=e^{-i\arccos{\frac{E}{|\Delta_{+}|}}}e^{-i\delta\theta}$. This implies the specularly reflected hole will acquire a phase of $(-\arccos{\frac{E}{|\Delta_{+}|}}-\delta\theta)$. Here, $\Delta_{+}=\Delta_{0}$ for the $s$-wave and $p$-wave pairings while $\Delta_{+}=\Delta_{0}\cos(2\theta-2\beta)$ for the $d_{x^2-y^2}$-wave pairing $(\beta=0)$ and the $d_{xy}$-wave pairing $(\beta=\frac{\pi}{4})$. In addition, $\delta=1$ for the $p$-wave SC and $\delta=0$ for the $d_{xy}$-wave SC.

On the other hand, when a spin-down hole is injected from WNSM, the coefficients $\tilde{a}_{\uparrow}$ for the specularly reflected electron and $\tilde{b}_{\downarrow}$ for the normally reflected hole can also be solved as
\begin{eqnarray}
\begin{split}
\tilde{a}_{\uparrow}&=\frac{4k_{x}k_{x}^{h\downarrow}u_{+}v_{-}e^{i\phi_{+}e^{i\phi_{-}}}}
{k_{x}^2\xi_{-}+k_{x}(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})\xi_{+}
+(z-i k_{x}^{e\uparrow})(z+i k_{x}^{h\downarrow})\xi_{-}},\\
\tilde{b}_{\downarrow}&=-\frac{k_{x}(k_{x}^{e\uparrow}-k_{x}^{h\downarrow})\xi_{+}
+k_{x}^2\xi_{-}-(k_{x}^{e\uparrow}+iz)(k_{x}^{h\downarrow}+iz)\xi_{-}}
{k_{x}^2\xi_{-}+k_{x}(k_{x}^{e\uparrow}+k_{x}^{h\downarrow})\xi_{+}
+(z-i k_{x}^{e\uparrow})(z+i k_{x}^{h\downarrow})\xi_{-}}.
\end{split}
\end{eqnarray}
Under the Andreev approximation, $\tilde{a}_{\uparrow}$ becomes
\begin{eqnarray}
\tilde{a}_{\uparrow}=\frac{v_{-}}{u_{-}}e^{i\phi_{-}}.
\end{eqnarray}
For the $s$-wave and $d_{x^2-y^2}$-wave pairings, $\tilde{a}_{\uparrow}=e^{-i\arccos{\frac{E}{|\Delta_{-}|}}}$ with $E<|\Delta_{-}|$, which means a phase of $(-\arccos{\frac{E}{|\Delta_{-}|}})$ is obtained by the specularly reflected electron. For the $d_{xy}$-wave and $p$-wave pairings, $\tilde{a}_{\uparrow}=-e^{-i\arccos{\frac{E}{|\Delta_{-}|}}}e^{-i\delta\theta}$ which means a phase of $(\pi-\arccos{\frac{E}{\Delta_{0}}}-\delta\theta)$ is obtained by the specularly reflected electron. Here, $\Delta_{-}=\Delta_{0}$ for the $s$-wave and $p$-wave pairings while $\Delta_{-}=\Delta_{0}\cos(2\theta+2\beta)$ for the $d_{x^2-y^2}$-wave pairing $(\beta=0)$ and the $d_{xy}$-wave pairing $(\beta=\frac{\pi}{4})$. In addition, we still have $\delta=1$ for the $p$-wave SC and $\delta=0$ for the $d_{xy}$-wave SC.

Along the closed path in the middle WNSM in Fig. \ref{wws}(a), there are twice SAR processes. One is for the right-going electron which is specularly reflected as hole at $x=L$. The other is for the right-going hole which is specularly reflected as electron . Combing the acquired phases, the total phase obtained by the two SAR processes can be derived as given in Eqs. (\ref{sdqqc}) and (\ref{pdqqc}). Note, we have assumed $\theta>0$ for the $d_{xy}$-wave SC in the above discussions. The situation of $\theta<0$ can be considered in a similar way, which will not change the final results in Eqs. (\ref{sdqqc}) and (\ref{pdqqc}).

\subsection{SAR probability and conductance for the WNSM-WNSM-SC double junctions}
The wave function in the left WNSM ($x<0$) for the spin-up electron incidence is given by
\begin{eqnarray}
\psi_{LW}=\left(\begin{array}{c}
1\\
0\end{array}\right)e^{ik^{e\uparrow}_{x}x}
+a_{\downarrow}\left(\begin{array}{c}
0\\
1\end{array}\right)e^{-ik^{h\downarrow}_{x}x}
+b_{\uparrow}\left(\begin{array}{c}
1\\
0\end{array}\right)e^{-ik^{e\uparrow}_{x}x},
\end{eqnarray}
with $a_{\downarrow}$ and $b_{\uparrow}$ still the SAR and the normal reflection coefficients. The wave function in the middle WNSM ($0<x<L$) is
\begin{eqnarray}
\begin{split}
\psi_{MW}=f_{1}\left(\begin{array}{c}
1\\
0\end{array}\right)e^{ik^{e\uparrow}_{x}x}
+f_{2}\left(\begin{array}{c}
1\\
0\end{array}\right)e^{-ik^{e\uparrow}_{x}x}\\
+
f_{3}\left(\begin{array}{c}
0\\
1\end{array}\right)e^{ik^{h\downarrow}_{x}x}
+f_{4}\left(\begin{array}{c}
0\\
1\end{array}\right)e^{-ik^{h\downarrow}_{x}x}
\end{split}
\end{eqnarray}
where $f_{1}$ and $f_{2}$ are coefficients for the right-going electron and the left-going electron and $f_{3}$ and $f_{4}$ are coefficients for the right-going hole and the left-going hole. The wave function in SC ($x>L$) is written as
\begin{eqnarray}
\Psi_{S}=c_{\uparrow}\left(\begin{array}{c}
u_{+}e^{i\phi_{+}}\\
v_{+}\end{array}\right)e^{ik_{x}x}
+d_{\downarrow}\left(\begin{array}{c}
v_{-}e^{i\phi_{-}}\\
u_{-}\end{array}\right)e^{-ik_{x}x},
\end{eqnarray}
with $c_{\uparrow}$ and $d_{\downarrow}$ still the coefficients for quasiparticle transmissions. For the spin-down electron incidence, the wave functions can also be given in a similar way.

Using the following boundary conditions,
\begin{eqnarray}
\begin{split}
&\psi_{LW}|_{x=0}=\psi_{MW}|_{x=0},\\
&\psi_{MW}|_{x=L}=\Psi_{S}|_{x=L},\\
&\psi'_{MW}|_{x=0}-\psi'_{LW}|_{x=0}=\frac{2mV_1}{\hbar^2}\hat{\tau}_{z}\psi_{LW}|_{x=0},\\
&\Psi'_{S}|_{x=L}-\hat{\tau}_{z}\psi'_{MW}|_{x=L}=\frac{2mV_2}{\hbar^2}\psi_{MW}|_{x=L},
\end{split}
\end{eqnarray}
$a_{\downarrow}$ and $b_{\uparrow}$ will be solved. The similar boundary conditions for the spin-down electron incidence can lead to the solutions of $a_{\uparrow}$ and $b_{\downarrow}$. The probabilities for SAR and the normal reflection are defined according to Eq. (\ref{prob}). The conductance of the double junctions can also be written as Eqs. (\ref{fencond}) and (\ref{totalcond}).

\section*{REFERENCES}


\section{Introduction}

In the past decades the $\Lambda$CDM model has successfully passed many tests \citep{Komatsu2011,Planck2014,
DelPopolo2013}, deserving to be dubbed "the standard model of cosmology". One of the assumption on which modern cosmology is based, is that General Relativity is the correct theory of gravitation. In this case, observations indicate the existence of a larger content of mass-energy than predicted \cite{DelPopolo2007,DelPopolo2014aa,Bull2016}.
The mass-energy of the universe is dominated by non-baryonic and non-relativistic particles, indicated as "cold dark matter" \cite{DelPopolo2014}, and a second component, dubbed "dark energy" (DE), a fluid with exotic properties, like that of having negative pressure, and giving rise to the accelerated expansion of the universe.
In its simplest form, and in the $\Lambda$CDM model, DE is represented by the cosmological constant
$\Lambda$.

Nevertheless the success of the $\Lambda$CDM model \cite{Spergel2003,Komatsu2011,DelPopolo2007}, precision data are revealing drawbacks, and tensions both at large scales \citep{Eriksen2004,Schwarz2004,Cruz2005,Copi2006,Macaulay2013,Planck2014_XVI,Raveri2016}, and at small ones
\cite{Moore1999,deBlok2010,Ostriker2003,BoylanKolchin2011,DelPopolo2014a,DelPopolo2014d,DelPopolo2017a}.

To start with, the particles that should constitute the DM has never been observed \cite{Klasen2015}, despite a large number of indirect evidences from small to large scales  \cite{Einasto2001,Bertone2005,Bouchet2004,Kilbinger2015}, and a large campaign of direct and indirect searches \cite{Bertone2005,Klasen2015, DelPopolo2014}.
{In particular, in a recent direct-detection experiment, XENON1T \citep{Aprile2018}, no significant excess over the background was found. For indirect searches, a large parameter space for annihilating dark matter was ruled out based on the radio data of the Andromeda galaxy \citep{Chan2019} and the gamma-ray data of the A2877 and Fornax clusters \citep{Chan2017}.}

The so called "small scale problems" of the $\Lambda$CDM \cite{DelPopolo2017a} are plaguing the model, and  several recipes, cosmological \citep{Zentner2003}, different nature of the dark matter particles \citep{Colin2000,Goodman2000,Hu2000,Kaplinghat2000,Peebles2000,SommerLarsen2001}, MG theories, e.g., $f(R)$ \citep{Buchdahl1970,Starobinsky1980}, $f(T)$ \citep{Bengochea2009,Linder2010,Dent2011,Zheng2011}, MOND \citep{Milgrom1983}
and astrophysical \cite{Brooks2013,Onorbe2015,DelPopolo2017a}, have been proposed to solve those problems.

Apart the issues related to the DM component of the $\Lambda$CDM model, the cosmological constant $\Lambda$
suffers from the ``cosmological constant fine tuning problem", and the ``cosmic  coincidence problem" \cite{Astashenok2012,Velten2014,Weinberg1989}.

The quoted issues motivated the investigations of other explanations, and models to clarify the universe accelerated expansion. These alternative models generate the DE effects through additional matter fields (e.g., quintessence \citep{Copeland2006}), or MG models\citep{Horndeski1974,Milgrom1983,Zwiebach1985,Moffat2006,Nojiri2005,Bekenstein2010,DeFelice2010,Linder2010,
Milgrom2014,Lovelock1971,Horava2009,Rodriguez2017,Horndeski1974,Deffayet2010,Deffayet2010}. In some cases, the quoted theories tried to explain the accelerated expansion as the manifestation of extra dimensions, or higher-order corrections effects, as in the Dvali-Gabadadze-Porrati (DGP) model \cite{Dvali2000} and in $f(R)$ gravity.
Disentangle between the plethora of models is not an easy task. The solution of the problem, or at least a better understanding of the same, may come from future surveys like:
Euclid\footnote{\url{http://www.euclid-ec.org}}, JDEM\footnote{\url{http://jdem.lbl.gov/}},
SKA\footnote{\url{https://www.skatelescope.org}}, LSST\footnote{\url{https://www.lsst.org}}, or
from new studies of the CMB \citep{Battye2018a,Battye2018c}.

DE and MG effects can be probed with structure formation. By means of hydro-dynamical simulations, one can obtain observables (e.g., splashback radius \cite{Adhikari2018}, the halo profile, the TAR \cite{Bhattacharya2017,Lopes2018}, or the mass-temperature relation (MTR) \cite{Hammami2017}) which can be directly compared with observations. In a recent paper \citep{DelPopolo2019}, we showed that the MTR relation cannot be used to put constraints on MG theories, since the $\Lambda$CDM model gives similar predictions to that of the MG theories.

Recently, the turnaround radius (TAR)\footnote{The TAR is defined as the distance from the center of a structure to the surface at which radial velocity is zero. The majority of works studies the maximum TAR, defined as the radius at which the radial acceleration is null, $\ddot{r} =0$.} has been proposed as a promising way to test cosmological models \citep{Lopes2018}, DE, and disentangle between $\Lambda$CDM model, DE, and MG models \citep{Pavlidou2014,Pavlidou2014a,Faraoni2015,Bhattacharya2017,Lopes2018,Lopes2019}. The attention on TAR increased when was shown that MG can affect the maximum TAR (see \citep{Lopes2019}).
According to some authors (e.g., \citep{Pavlidou2014}), the TAR is a well-defined, and unambiguous boundary of a structure, in the spherical collapse model (SCM), simulations, analytic calculations, and is clean from baryons physics. Concerning the last point, contrarily to \citep{Pavlidou2014}, we already showed in \citep{DelPopolo2013a,DelPopolo2013b,Pace2014,Mehrabi2017,Pace2019} that the non-linear equation driving the evolution of the overdensity contrast depends on shear, and vorticity, which conversely depends from the mass of the forming structure, and the way it forms. In the following, we will also show that the TAR also depends from dynamical friction. All the papers dealing with the determination of TAR, usually find it from some metric, and equations not taking into account the real physics of structure formation. Then, we will show in this paper, that TAR depends from baryons physics.

\citep{Pavlidou2014} calculated the TAR for $\Lambda$CDM, and \citep{Pavlidou2014a} did the same for smooth DE.
{\citep{Lee2016} proposed to use the zero velocity surface of large structures to look for
the violation of the maximum upper bound of $R_{\rm t}$.
}

In MG theories \citep{Capozziello2019} found a general relation for the maximum TAR in $f(R)$ theories, and \citep{Faraoni2015} found a method to get the same quantities in generic gravitational theories.

In the present paper, we use an extended spherical collapse model (ESCM) introduced, and adopted in   \citep{DelPopolo2013,DelPopolo2013a,Pace2014,Mehrabi2017,Pace2019}. The ESCM takes into account the effect of a non null shear and vorticity on the collapse, to show how the TAR is changed.
As we have already shown, shear and vorticity change the typical parameters of the SCM \citep{DelPopolo2013}, the mass function \citep{DelPopolo2013,DelPopolo2013a,Pace2014,Mehrabi2017}, the two-point correlation function \citep{DelPopolo1999}, and the weak lensing peaks \citep{Pace2019}. We will find a more general formula than that obtained by \citep{Pavlidou2014} for smooth DE, showing that shear, and vorticity, give rise to a smaller TAR especially at the galactic mass scales. Then, we will find the TAR - mass relation for the $\Lambda$CDM model, and several dynamical DE models, and compare our results to the prediction of the $R_{\rm t}$-$M_{\rm t}$ relation of \citep{Lopes2018} for $f(R)$ models. The result of the comparison shows that the $R_{\rm t}$-$M_{\rm t}$ relation of the $\Lambda$CDM model has a very similar behavior to that of the $f(R)$ models.
This makes impossible to disentangle between the MG results and those of GR.
A similar situation happens for the DE models.

The paper is organized as follows.
Section~\ref{sect:model} describes the model used to derive the $R_{\rm t}$-$M_{\rm t}$ relation in $\Lambda$CDM, and smooth DE cosmologies.
Section~\ref{sect:Results} is devoted to the presentation and the discussion of our results.
Section~\ref{sect:conclusions} is devoted to conclusions.

\section{The Model}\label{sect:model}

The SCM introduced by \cite{Gunn1972}, extended and improved in several following papers \citep{Fillmore1984,Bertschinger1985,Hoffman1985,Ryden1987,Subramanian2000,
Ascasibar2004,Williams2004}, is a very popular method to study analytically the non-linear evolution of perturbations of DM and DE.
Shortly, the model describes the evolution of a spherical symmetric overdensity, how it decouples from the Hubble flow, reach a maximum radius, dubbed TAR, and finally collapse and virialize\footnote{
In the virialization process, the collapse kinetic energy is converted into random motions.}.
\cite{Gunn1972} is a very simple SCM model just considering the radial collapse of the structure. It does not take account of tidal angular momentum \citep{Peebles1969,White1984}, random angular momentum \citep{Ryden1987,Ryden1988,Williams2004}, dynamical friction (\citep{AntonuccioDelogu1994,Delpopolo2009}), etc. The way to take account of angular momentum was studied in several papers\citep{Ryden1987,Gurevich1988a,Gurevich1988b,White1992,Sikivie1997,Nusser2001,Hiotelis2002,
LeDelliou2003,Ascasibar2004,Williams2004,Zukin2010}\footnote{Particles angular momentum is distributed in random directions in order the mean angular momentum at any point in space is zero \cite{White1992,Nusser2001}}, that of dynamical friction was studied in \citep{AntonuccioDelogu1994,Delpopolo2009}, while \citep{Hoffman1986,Hoffman1989,Zaroubi1993} discussed the role of shear in the gravitational collapse.

The SCM
with negligible DE perturbations was studied by
\citep[see, e.g.][]{Bernardeau1994,Bardeen1986,Ohta2003,Ohta2004,Basilakos2009,Pace2010,Basilakos2010}
while in \citep[see][]{Mota2004,Nunes2006,Abramo2007,Abramo2008,Abramo2009a,Abramo2009b,Creminelli2010,Basse2011,Batista2013}
the effects of the DE fluid perturbation were taken into account.
\citep{Pettorino2008,Wintergerst2010a,Tarrant2012} extended the SCM model to coupled DE models, and
\citep{Pettorino2008,Pace2014} to extended DE (scalar-tensor) models.

\cite{DelPopolo2013a,DelPopolo2013b} studied the effects of shear and rotation in smooth DE models.
The effects of shear and rotation were investigate in \cite{DelPopolo2013a,DelPopolo2013b} for smooth DE models, \cite{Pace2014b} in clustering DE cosmologies, and
\cite{DelPopolo2013c} in Chaplygin cosmologies.

\subsection{The ESCM}
Here, we show how the evolution equations of $\delta$ in the non-linear regime can be obtained.
Those equations were obtained by \cite{Bernardeau1994,Padmanabhan1996,Ohta2003,Ohta2004,Abramo2007,Pace2010} in the spherical and ellipsoidal collapse models scheme. We assume the equation of state $P=w\rho c^2$ for the fluid, and the Neo-Newtonian expressions for the relativistic Poisson equation, the Euler, and continuity equations \citep{Lima1997}
\begin{eqnarray}
  \frac{\partial\rho}{\partial t}+\nabla_{\vec{r}}\cdot(\rho\vec{v})+
  \frac{P}{c^2}\nabla_{\vec{r}}\cdot\vec{v} & = & 0 \label{eqn:cnpert}\;,\\
  \frac{\partial\vec{v}}{\partial
    t}+(\vec{v}\cdot\nabla_{\vec{r}})\vec{v}+
  \nabla_{\vec{r}}\Phi +{\bf \frac{c^2}{c^2 \rho+P} \nabla P} & = & 0\;, \label{eqn:enpert}\\
  \nabla^2\Phi-4\pi G\left(\rho+\frac{3P}{c^2}\right) & = & 0\;,\label{eqn:pnpert}
\end{eqnarray}
where we indicated with $\vec{r}$ the physical coordinate, with $\Phi$ the Newtonian gravitational potential, and the velocity in three-space is given by $\vec v$.

The perturbations equations, using comoving coordinates $\vec{x}=\vec{r}/a$ are given by
\begin{eqnarray}
 \dot{\delta}+(1+w)(1+\delta)\nabla_{\vec{x}}\cdot\vec{u} & = & 0\;, \label{eq:pertCont}\\
 \frac{\partial \vec{u}}{\partial t}+2H\vec{u}+(\vec{u}\cdot\nabla_{\vec{x}})\vec{u}+\frac{1}{a^2}\nabla_{\vec{x}}\phi &
 = & 0\;,\label{eq:pertEuler} \\
 \nabla_{\vec{x}}^2\phi-4\pi G(1+3w)a^2\bar{\rho}\delta & = & 0\;, \label{eq:pertPois}
\end{eqnarray}
being $\vec{u}(\vec{x},t)$ the comoving peculiar velocity, and $H(a)$ is the Hubble function.
{\bf Note that in Eq. (\ref{eq:pertEuler}) the pressure terms are zero, because
in the top-hat, with $\rho=\overline{\rho} +\delta \rho$, and $P=\overline{P} +\delta P$ the gradient is zero.}

The non-linear evolution equation is obtained combining the previous equations, which in the case $w=0$ (dust)
is given by:
\begin{equation}\label{eqn:nleq11}
\begin{split}
\ddot{\delta}+2H\dot{\delta}-\frac{4}{3}\frac{\dot{\delta}^2}{1+\delta}-
4\pi G\bar{\rho}\delta(1+\delta)-\\
(1+\delta)(\sigma^2-\omega^2) & = 0\;
\end{split}
\end{equation}
Eq. (\ref{eqn:nleq11}) is Eq.~41 of \cite{Ohta2003}, and a generalization of Eq.~7 of \cite{Abramo2007} to the case of a non-spherical configuration of a rotating fluid.

In Eq. (\ref{eqn:nleq11}), $\sigma^2=\sigma_{ij}\sigma^{ij}$, and $\omega^2=\omega_{ij}\omega^{ij}$ are the shear, and rotation term, respectively. The shear term is related to a symmetric traceless tensor, dubbed shear tensor, while rotation term is related to an antisymmetric tensor, given by
\begin{eqnarray}
 \sigma_{ij} & = & \frac{1}{2}\left(\frac{\partial u^j}{\partial x^i}+\frac{\partial u^i}{\partial x^j}\right)
 -\frac{1}{3}\theta\delta_{ij}\;, \\
 \omega_{ij} & = & \frac{1}{2}\left(\frac{\partial u^j}{\partial x^i}-\frac{\partial u^i}{\partial x^j}\right)\;.
\label{eq:tens}
\end{eqnarray}
being $\theta=\nabla_{\vec{x}}\cdot\vec{u}$ the expansion.

In terms of the scale factor, $a$,
Eq. (\ref{eqn:nleq11}) takes the form
\begin{equation}\label{eqn:wnldeq}
 \begin{split}
  \delta^{\prime\prime}+\left(\frac{3}{a}+\frac{E^\prime}{E} \right)
\delta^\prime-\frac{4}{3}\frac{\delta^{\prime 2}}{1+\delta}-
  \frac{3}{2}\frac{\Omega_{\mathrm{m},0}}{a^5 E^2(a)}\delta(1+\delta)-&\\
  \frac{1}{a^2H^2(a)}(1+\delta)(\sigma^2-\omega^2)&=0\;,
 \end{split}
\end{equation}
where $\Omega_{m,0}$ is DM density parameter at $t=0$ ($a=1$), and $E(a)$ is given by
\begin{equation}\label{eq:e}
E(a)=\sqrt{\frac{\Omega_{\mathrm{m},0}}{a^3}+\frac{\Omega_{\mathrm{K},0}}{a^2}+ \Omega_{\mathrm{Q},0}g(a)}\;,
\end{equation}
being
\begin{equation}
g(a)=\exp{\left(-3\int_1^a\frac{1+w(a')}{a'}~da^{\prime}\right)}\;.
\end{equation}

{
Recalling that $\delta=\frac{2GM_{\rm m}}{\Omega_{m,0} H^2_0}(a/R)^3-1$, where $R$ is the effective perturbation radius and inserting into
Eq.~(\ref{eqn:nleq11}), we get

\begin{eqnarray} \label{eqn:wnldeq0}
 \ddot{R} &=& -\frac{GM_{\rm m}}{R^2} - \frac{GM_{\rm de}}{R^2}(1+3w_{\rm de})-\frac{\sigma^2-\omega^2}{3}R= \nonumber\\
& &
-\frac{GM_{\rm m}}{R^2} - \frac{4 \pi G \bar{\rho_{\rm de}} R}{3}(1+3w_{\rm de})-\frac{\sigma^2-\omega^2}{3}R\,,
\end{eqnarray}
{being $M_{\rm m}= \frac{4 \pi R^3}{3} (\bar{\rho}+\delta \rho)$ and $M_{\rm de}$ is the mass of the dark-energy component enclosed in the volume, $\bar{\rho}_{\rm de}$, and $w_{\rm de}$ being respectively its background density and equation-of-state \citep{Fosalba1998a,Engineer2000,Ohta2003,Pace2019}.
$M_{\rm m}$, as shown, contains background and perturbation.}

Eq. (\ref{eqn:wnldeq0}), for $w=-1$, becomes
\begin{equation} \label{eqn:spher0}
 \ddot{R}=-\frac{GM_{\rm m}}{R^2} -\frac{\sigma^2-\omega^2}{3}R+ \frac{\Lambda}{3} R\
\end{equation}
and is similar to the usual expression for the SCM with angular momentum \citep[e.g.][]{Peebles1993,Nusser2001,Zukin2010}, and cosmological constant:
\begin{equation}
\frac{d^2 R}{d t^2}= -\frac{GM}{R^2} +\frac{L^2}{M^2 R^3}+\frac{\Lambda}{3}R= -\frac{GM}{R^2} + \frac{4}{25} \Omega^2 R+ \frac{\Lambda}{3}R,
\label{eqn:spher}
\end{equation}
The term $\frac{4}{25} \Omega^2 R$ comes from the expression of angular momentum, $L=I \Omega$, and the
momentum of inertia of a sphere, $I=2/5 M R^2$.

{
In the simple case of a uniform rotation with angular velocity $\Omega=\Omega_{\rm z} {\bf e}_{\rm z}$, we have that $\Omega=\omega/2$ (see also \cite{Chernin1993}, for a more complex and complete treatment of the interrelation of vorticity and angular momentum in galaxies). Then, the previous equations show a strict connection between vorticity, $\omega$, and angular velocity $\Omega$.

It is then convenient to define the dimensionless $\alpha$-number as the ratio between the rotational and the gravitational term in Eq. (\ref{eqn:spher}):
\begin{equation}
\alpha(M)=\frac{L^2}{M^3RG}
\end{equation}
The above ratio, $\alpha$, is mass dependent. It has maximum values for galaxies\footnote{In the case of our galay it is $\simeq 0.1$.}, and decreases going toward clusters of galaxies.

In order to solve Eq. (\ref{eqn:wnldeq}), we should know how the term $\sigma^2-\omega^2$ depends on the density contrast. This can be done using the above outlined argument for rotation, namely the connection between angular momentum and shear, and recalling that Eq. (\ref{eqn:spher}) from which $\alpha$ was obtained is equivalent to Eq. (\ref{eqn:wnldeq0}) which is also equivalent to Eq. (\ref{eqn:nleq11}), and Eq. (\ref{eqn:wnldeq}).

One may then calculate the same ratio between the gravitational and the extra term appearing in Eq. (\ref{eqn:nleq11}) or Eq. (\ref{eqn:wnldeq}) thereby obtaining
\begin{equation} \label{eq:nonlinf}
\frac{\sigma^2-\omega^2}{H^2_0}=-\frac{3}{2}\frac{\alpha \Omega_{\rm m,0}}{a^3}\delta.
\end{equation}

As shown in \citep{DelPopolo2013b} this is a reasonable assumption, and it was also used in \citep{DelPopolo2013a,DelPopolo2013b,Pace2014,Mehrabi2017}.
}

In \citep{DelPopolo2013a,DelPopolo2013b} it was obtained calculating the threshold of collapse $\delta_c$, and matching it with $\delta_c$ obtained by \citep{Sheth1999}.
}

Then the nonlinear equation to solve is obtained substituting Eq. (\ref{eq:nonlinf}) into Eq. (\ref{eqn:wnldeq})
\begin{equation}\label{eqn:wnldeqq}
  \delta^{\prime\prime}+\left(\frac{3}{a}+\frac{E^\prime}{E} \right)
\delta^\prime-\frac{4}{3}\frac{\delta^{\prime 2}}{1+\delta}-
  \frac{3}{2}\frac{\Omega_{\mathrm{m},0} (1-\alpha)}{a^5 E^2(a)}\delta(1+\delta)=0.
\end{equation}

The threshold of collapse, and the turnaround, can be obtained solving Eq. (\ref{eqn:wnldeqq}) following the method described in \citep{Pace2010}, or solving Eq. (\ref{eqn:wnldeq0}).

An important point is that Eq. (\ref{eqn:wnldeq0}), can be written in a more general form taking into account dynamical friction
 \citep{Kashlinsky1986,Kashlinsky1987,Lahav1991,Bartlett1993,
 AntonuccioDelogu1994,Peebles1993,DelPopolo1998,DelPopolo1998a,DelPopolo2006b,Delpopolo2009,DelPopolo2019}
\begin{equation}\label{eq:coll}
 \ddot{R} = -\frac{GM}{R^2} + \frac{L^2(R)}{M^{2}R^3} + \frac{\Lambda}{3}R -
 \eta\frac{{\rm d}R}{{\rm d}t}\,,
\end{equation}
being $\eta$ the dynamical friction coefficient.
Eq. (\ref{eq:coll}) can be obtained via Liouville's theorem \citep{DelPopolo1999}, and the dynamical friction force per unit mass, $\eta\frac{{\rm d}R}{{\rm d}t}$, is given in \citep{Delpopolo2009}(Appendix D, Eq. D5), and \citep{DelPopolo2006b}, Eq. 5).

In terms of shear, vorticity, and DE, the equation, similarly to Eq. (\ref{eqn:wnldeq0}), can be written as
\begin{equation} \label{eqn:wnldeqqq}
\ddot{R} = -\frac{GM_{\rm m}}{R^2} - \frac{GM_{\rm de}}{R^2}(1+3w_{\rm de})-\frac{\sigma^2-\omega^2}{3}R -\eta \frac{dR}{dt}.
\end{equation}

A similar equation (excluding the dynamical friction term) was obtained by several authors
\citep[e.g.,][]{Fosalba1998a,Engineer2000,DelPopolo2013b}) and generalized to smooth DE models in
\cite{Pace2019}.

Eq. (\ref{eqn:wnldeqqq}), and Eq. (\ref{eqn:wnldeq0}) differs for the presence of the dynamical friction term. As shown in several papers (e.g., \citep{AntonuccioDelogu1994,Delpopolo2009,DelPopolo2006b,DelPopolo2019} dynamical friction has a similar effect to that of rotation, and cosmological constant: it delays the collapse of a structure (perturbation).
As shown in Fig. 1 of \citep{DelPopolo2017}, and in Fig. 11 of \citep{Delpopolo2009}, the magnitude of the effect of rotation, dynamical friction, and cosmological constant are of the same order of magnitude with differences of a few percent.
Willingly, by means of the relation $\delta=\frac{2GM_{\rm m}}{\Omega_{m,0} H^2_0}(a/R)^3-1$, Eq. (\ref{eqn:wnldeqqq}) can be written in terms of $\delta$, similarly to Eq. (\ref{eqn:wnldeq}).

Summarizing, the typical quantities of the SCM depends from shear, vorticity, and dynamical friction.
Contrarily to \citep{Pavlidou2014}, and several other authors (e.g., \citep{Lopes2018,Bhattacharya2017}), the SCM results depend from the baryon physics, since shear, rotation, and dynamical friction depends from the mass, and the way structures formed.

In the following, we specialize on the effect shear, and rotation has on the SCM, and how the TAR is modified.
To this aim, in this paper, dynamical friction will be neglected.

\section{Results}\label{sect:Results}

\begin{figure*}[!ht]
 \centering
 \includegraphics[width=22cm,height=6cm,angle=0]{dc2_6.eps}
 \caption[justified]{The threshold of collapse $\delta_c$, as function of redshift, and mass.
In the left panel, we plot $\delta_c$ vs redshift. The red line represents the predictions of the ESCM for  $\delta_c(z)$ for the $\Lambda$CDM model. From top to bottom, the value of $\alpha$ varies from 0.05 (solid line) (galactic mass scale)
to 0.03 (long dashed line),
to 0.01 (dotted line) (cluster mass scale),
and to 0 (dashed line). Left panel: $\delta_c$ vs mass. The solid red line represents the result of the ESCM model, the red dashed line that of {the elliptictal collapse model of} \citep{Sheth2001}, and the brown line the $\Lambda$CDM expectation when no shear and rotation are taken into account.}
 \label{fig:comparison}
\end{figure*}

Several authors have studied the effect of shear \citep{Hoffman1986,Hoffman1989, Zaroubi1993}, and rotation \citep{DelPopolo1998,DelPopolo1999,DelPopolo2000,DelPopolo2001,DelPopolo2002,DelPopolo2013b,Pace2019} on the SCM and elliptical collapse  model.
The effect of shear and rotation is that of slowing down the collapse \citep{Peebles1990,Audit1997,DelPopolo2001,DelPopolo2002}, changing turnaround epoch and collapse time, as well as the typical parameters of the SCM.
As a consequence also the mass function \citep{DelPopolo1999,DelPopolo2000,DelPopolo2013b,Pace2014,Mehrabi2017,DelPopolo2017,Pace2019}, the two-point correlation function, scaling relations like the mass-temperature \citep{DelPopolo2002a,DelPopolo2019}, and luminosity-temperature relation \citep{DelPopolo2005}, are modified.

\subsection{Threshold of collapse with shear, and rotation}

\begin{figure*}[!ht]
 \centering
\includegraphics[width=25cm,angle=0]{RtAS_22.eps}
 \caption[justified]{Turnaround radius $R_{\rm t}$ vs mass, $M_{\rm t}$. Left panel: TAR predicted by the standard SCM (solid red line), and the ESCM (dashed red line)
for the $\Lambda$CDM model. Right panel: the result of the ESCM for the $\Lambda$CDM model (red dashed line), and that of the AS model (black dotted line).
}
 \label{fig:comparison}
\end{figure*}

In order to show the effects of shear and rotation on the SCM predictions, we show, in Fig. 1, how the collapse threshold,  $\delta_c$, is modified. In the left panel, we show the dependence of $\delta_c$ from the redshift. The red line represents the predictions of the ESCM for  $\delta_c(z)$ for the $\Lambda$CDM model. From top to bottom, the value of $\alpha$ varies from 0.05 (solid line) corresponding to a mass $\simeq 10^{11} M_{\odot}$, to 0.03 (long dashed line), corresponding to a mass $\simeq 10^{13} M_{\odot}$, to 0.01 (dotted line), corresponding to a mass $\simeq 10^{15} M_{\odot}$, and to 0 (dashed line).
We also show the predictions of the extended SCM in the case of one DE model, the Albrecht-Skordis \citep{Albrecht2000} (AS) DE model (black line). We chose to plot this model, because the other quintessence models considered in previous papers \citep{Pace2010,DelPopolo2013a} INV1 ($w_0=-0.4$), INV2 ($w_0=-0.79$), 2EXP ($w_0=-1$), CNR ($w_0=-1$), CPL ($w_0=-1$), SUGRA ($w_0=-0.82$) (see \citep{Pace2010,DelPopolo2013b})\footnote{$w_0$ is the value of $w$ nowadays} have a similar behavior to that of the AS model, and because the values of $\delta_c(z)$ predicted by them is located in the region between the $\Lambda$CDM and the AS model (see Fig. 4 of \citep{Pace2010}).
The inclusion of the term proportional to $\sigma^2-\omega^2$ changes only the value of the linear overdensity parameter. The ratios with the $\Lambda$CDM model does not change. This means, that
{if the $\delta_c$ relative to the $\Lambda$CDM is higher than $\delta_c$ relative to a given DE model, the same relationship will remain after taking into account the $\sigma^2-\omega^2$ term. This can be seen comparing the top left panel in Fig. 4 of \citep{Pace2010} with the top left panel in Fig. 2 of \citep{DelPopolo2013a}, and also looking at Fig. 1 of the present paper.}
The DE models similarly to the $\Lambda$CDM model are "retarded" (the collapse is slowed down) by the presence of the shear, and rotation term,
and their values for $\delta_c(z)$ are all smaller than the $\delta_c(z)$ of the $\Lambda$CDM model obtained with the ESCM. This is expected, since at early times the DE amount is larger, and consequently the value of $\delta_c(z)$ must be lower to have objects collapsing.

\begin{figure*}[!ht]
 \centering
 \includegraphics[width=12cm,angle=0]{pavl5.eps}
 \caption[justified]{Mass-radius relation of stable structures for different $w$. The solid lines from top to bottom represent $w=-0.5$ (solid green line), -1 (black solid line), -1.5 (blue solid line) ,-2 (pink solid line), -2.5 (red solid line). The dashed lines are the same of the previous lines, but they are obtained using the ESCM model. The dots with error bars, are data as in \citep{Pavlidou2014a}.}
 \label{fig:comparison}
\end{figure*}

The plot also shows the different behavior of $\delta_c(z)$ in the case of no shear, and rotation ($\alpha=0$). In this case, $\delta_c$ has a weak dependence from redshift in the $z$ range $[0,2]$, and then becomes constant, with a value equal to that of the Einstein de Sitter model. Shear, and rotation ($\alpha \neq 0$), change the behavior of $\delta_c(z)$, as follows: 1. it assumes larger values with respect the case shear and rotation are not present; 2. the larger is $\alpha$, the larger is the difference between the values of $\delta_c(z)$; 3. $\delta_c(z)$ has a decaying behavior. The value of $\delta_c(z)$ for $\alpha=0.05$, is $\simeq 30\%$ larger than in the $\alpha=0$ case. In the right panel of Fig. 1, we plot $\delta_c(M)$ versus the mass. When shear and rotation are not present the value of $\delta_c$ is constant (brown line), while when they are taken into account the threshold $\delta_c$ becomes mass dependent, and is characterized by larger values at smaller masses, converging to the standard $\Lambda$CDM model constant value moving to the largest clusters. The dashed line shows the result of \citep{Sheth1999,Sheth2001}.
The behavior of the threshold implies that less massive perturbations (e.g., galaxies) in order to form structures must cross a higher threshold than more massive ones. In terms of the peak formalism, and the peak height $\nu=\delta_c/\sigma(M)$\footnote{$\sigma(M)$ is the mass variance}, the angular momentum acquired by a peak is proportional to turnaround time, $t_{\rm ta}$, and anticorrelated with the peak height $j \propto t_{\rm ta} \propto \nu^{-3/2}$ \citep{Hoffman1986,Delpopolo2009,Polisensky2015}.
Since low peaks acquire larger angular momentum than high peaks, they need a higher density contrast to collapse and form structures \citep{DelPopolo1998,DelPopolo2001,DelPopolo2002,Delpopolo2009,Ryden1988,Peebles1990,Audit1997}.
The cosmological constant has similar effects to that of angular momentum in collapse, slowing down it, but its effects vanishes at high redshift.
\subsection{TAR, in the $\Lambda$CDM, and DE models}

\begin{figure*}[!ht]
 \centering
\includegraphics[width=18cm,angle=0]{RTm3.eps}
 \caption[justified]{Turnaround radius, $R_{\rm t}$ vs mass, $M_{\rm t}$. Left panel: the result obtained by \citep{Lopes2018} for $R_{\rm ta}$ vs mass for the $\Lambda$CDM model, the model with $\epsilon=1/3$, and $f(R)$ with $f_{\rm R0}=10^{-6}$, $10^{-5}$, and $10^{-4}$. Right panel: the same as in left panel compared with the prediction of the ESCM. The red band represents our prediction for the $\Lambda$CDM model, obtained with the ESCM, with the 68\% confidence level region.}
 \label{fig:comparison}
\end{figure*}

\begin{figure*}[!ht]
 \centering
\includegraphics[width=18cm,angle=0]{RTmde.eps}
 \caption[justified]{The turnaround radius, $R_{\rm t}$ vs mass, $M_{\rm t}$ obtained with the ESCM for the DE models. In the left panel, from top to bottom, the cyan, blue, brown, magenta, black, red, and green lines represent the INV1, INV2, SUGRA, $w_{09}$\footnote{Namely the model having $w=-0.9$}, AS, $\Lambda$CDM without shear, and rotation, and $\Lambda$CDM with shear, and rotation, respectively. In the right panel, the solid lines are the $\Lambda$CDM, and $f(R)$ models of \citep{Lopes2018} (similarly to the left panel of our previous figure), while the blue dashed line, and the green dashed lines, are the INV2 model, and the $\Lambda$CDM with shear, and rotation plotted in the left panel.}
 \label{fig:comparison}
\end{figure*}

Our main goal is to understand how shear and rotation changes the TAR. We start finding an expression for the TAR. {Just in order to compare with \citep{Pavlidou2014} and \citep{Pavlidou2014a} results, we calculate the maximum TAR, MTAR,
namely the radius of the surface where $\ddot R =0$}. Using Eq. (\ref{eqn:wnldeq0}), we obtain
\begin{equation} \label{eq:turnn}
R_{\rm ta}= \left [\frac{-3M}{4 \pi \rho_{\rm de}(1+3w)+(\sigma^2-\omega^2)} \right]^{1/3}
\end{equation}
In the case of the $\Lambda$CDM model, ($w=-1$), no shear, and rotation, Eq. (\ref{eq:turnn}) reduces to
\begin{equation} \label{eq:turn1}
R_{\rm ta}= \left [\frac{3GM}{\Lambda} \right]^{1/3}
\end{equation}
which is the same as that obtained by \citep{Pavlidou2014}. In the case of DE models, with no shear, and rotation, Eq. (\ref{eq:turnn}) reduces to
\begin{equation} \label{eq:pavl}
R_{\rm ta}= \left [\frac{-3M}{4 \pi \rho_{\rm de}(1+3w)} \right]^{1/3}
\end{equation}
which is the same as that obtained by \citep{Pavlidou2014a}.
In other terms, the turn around obtained by \citep{Pavlidou2014a} is a peculiar case of the more general form given by Eq. (\ref{eq:turnn}).

{In the paper, we will get and plot the TAR, not MTAR, since we compare with \citep{Lopes2018}, which calculated the TAR.}

When taking into account shear, and rotation, the TAR is smaller than expected from \citep{Pavlidou2014a} result, as shown in Fig. 2. In the left panel of Fig. 2, we plot the TAR predicted by the standard SCM (solid red line), and the ESCM (dashed red line)
for the $\Lambda$CDM model. As expected, similarly to the threshold of collapse, $\delta_c$, when shear and rotation are taken into account, the collapse is slowed down and the turn around is smaller.
The differences between the standard SCM predictions, and the ESCM, are larger at smaller masses, because smaller mass objects have larger rotation.
In the right panel of Fig. 2, we compare the result of the ESCM for the $\Lambda$CDM model (red dashed line) with that of the AS model (black dotted line).
The AS model has a slightly larger TAR with respect to the $\Lambda$CDM model.
\citep{Pavlidou2014a} tried to constrain the equation of state by means of Eq. (\ref{eq:pavl}), comparing the predicted TAR at different $w$ with mass and radii of small $z$ cosmic structures.
In this way is possible to constrain the equation of state parameter $w$ at $z \simeq 0$ ($w_0$)\footnote{Many of the DE models can be described by an equation of state parameter, $w(a)$, depending on its value nowadays, $w_0$, that at matter-radiation equality epoch, and some other parameters at the same epoch (see Eq. 23 of \citep{Pace2010}). In order to constrain the evolution of the equation of state parameter, $w$, it is necessary to use structures at high redshift. }. They used data for: Milky Way (MW), M81/M82 group, Local Group, Virgo cluster, Fornax-Eridanus group. \citep{Peirani2008} took into account the cosmological constant and obtained different values for M31/MW, M81, Virgo, and also obtained the R-M values for CenA/M83, IC342, NGC253. We used the data obtained by \citep{Peirani2008}, and that of \citep{Pavlidou2014a} for Local Group, and Fornax-Eridanus.

In Fig. 2, the solid lines are obtained from the equation of TAR (Eq. \ref{eq:pavl}) not taking into account shear, and rotation. The dashed lines are the corrections obtained when shear, and rotation are taken into account. The lines correspond to $w=-2.5$, -2, -1.5, -1, -0.5, from bottom to top. The region of the parameter space above each line gives the range of $w$ for which no stable structures should exist. \citep{Pavlidou2014a} discussed some constraints to $w$, based on the highest mass objects. Here we want just to stress that at masses smaller than $10^{14} M_{\odot}$ the TAR is modified by the presence of shear, and rotation. As a consequence, structures at smaller masses can give different constraints to $w$.

For example, the Local Group gives the constrain $w \geq -2.5$, while if shear, and rotation are not taken into account, $w$ is more negative.
{The other constraints are reported in Table 1}.
{Therefore, including shear and rotation can help to constrain the equation of state of DE by means of the data of galaxies. Furthermore,
can only marginally satisfy the data of NGC253 and M81, the standard $\Lambda$CDM model can also be strongly constrained by galactic data.}

\subsection{Comparison with TAR in $f(R)$ theories}

\citep{Lopes2018} investigated the evolution of the TAR, and its dependence from structure mass in General Relativity (GR), and in the $f(R)$ theories. Their Eq. (3.3) takes into account the effects of MG through the parameter $\epsilon(a,k)$, where $k$ is the (angular) wavenumber. When $\epsilon(a,k)$ is zero GR is recovered, and their Eq. (3.3) is equal to our Eq. (\ref{eqn:wnldeq}) without shear, and rotation. In other terms, \citep{Lopes2018} includes the effects of the $f(R)$ theories but not that of shear, and rotation, which changes the behavior of the main quantities of the SCM. As a consequence,
their Figs. 1, 3 predict a monotonic increase for $\delta_c(z)$ as
in the $\Lambda$CDM, and DE models without shear, and rotation, as shown in \citep{Pace2010}, and in the left panel of our Fig. 1 (case $\alpha=0$). Shear, and rotation modify this behavior.
Similarly, their Fig. 4 shows that $\delta_c(M)$, differs from the $\Lambda$CDM predictions by $\simeq 1\%$.
As shown by \citep{Sheth2001}, and by \citep{DelPopolo2013a,DelPopolo2013b,Pace2014,Mehrabi2017,DelPopolo2017},
in order to have a mass function reproducing simulations, the threshold
must be: 1. mass dependent, and 2. a decaying function of mass.
For precision's sake, in the excursion set formalism, the halo statistics, and the mass function, depends from the statistical properties of  the average overdensity, $\overline{\delta}_{\rm R_f}$, within a window of a given radius, $R_f$. $\overline{\delta}_{\rm R_f}$ vs $R_f$ is a random walk \citep{Bond1991,DelPopolo2017}, and haloes form when the random walk crosses the threshold, $\delta_c$, dubbed barrier \citep{Sheth2001,DelPopolo2007a,DelPopolo2017}. A mass function that reproduces simulations and observations must have a mass dependent, and monotonically decreasing $\delta_c$, as in \citep{Sheth2001,Sheth2002},
and \citep{DelPopolo1998,DelPopolo1999,DelPopolo2017}. This shows that the results of \citep{Lopes2018} cannot reproduce the mass function obtained in simulations, and observations, because the effects of
shear, and rotation
are not taken into account.
~\\

Another important point is the use of the TAR to disentangle between GR and MG, discussed by several authors (e.g., \citep{Lopes2018}).
In the following, we will compare \citep{Lopes2018} results with the prediction of the ESCM for the TAR in the $\Lambda$CDM, and DE models. To this aim, we compared in Fig. 4, their result concerning the $R_{\rm t}$ vs mass with the $\Lambda$CDM model predictions in the ESCM. In the figure, we plot the result obtained by \citep{Lopes2018} for $R_{\rm t}$ vs mass for the $\Lambda$CDM model, the model with $\epsilon=1/3$, and $f(R)$ with $f_{\rm R0}=10^{-6}$, $10^{-5}$, and $10^{-4}$. The red band represents our prediction for the $\Lambda$CDM model, obtained with the ESCM, with the 68\% confidence level region.
{In order to calculate the 68\% confidence level, we used a Monte Carlo simulation.
In short, this approach involves the repeated simulation of samples within the probability density functions of the input data. In the first step, one looks for the sources of uncertainty, and identifies the probability density functions (PDFs) that gives a good fit to those the data sources.
For all data inputs identified as sources of uncertainty, Monte Carlo simulations are run, and the resulting simulations are applied to models (equations), in order to get the distribution of the quantity of interest. Finally one can get the confidence intervals needed.
In our case, in order to obtain the TAR, one has to solve Eq. (\ref{eqn:wnldeq0}), or Eq. (\ref{eqn:wnldeqq}), using the condition $\dot R=0$ at turnaround or equivalently a relation between
$\delta$, $\delta^{'}_t$, and at turn around (see eq. 3.8 of \citep{Lopes2018}). If for example, we consider Eq. (\ref{eqn:wnldeq0}), we may write it as
\begin{equation}
\dot R^2=0=2 \frac{GM}{R}-\frac{4 \pi \rho_{\rm de}}{3} +\frac{2 H^2_0 \Omega_{m,0}}{2 a^3 G}
\int_{R_i}^{R_t} \frac{L^2}
{M^3} \delta dR
\end{equation}
where $R_i$ is the initial radius of the perturbation, related to $z_i$ given in the appendix.
One can get numericaly the relation between $R_t$, and $M_t$. In order to obtain the confidence interval by means of the Monte Carlo method, we observe that $R_t$ depends from mass, overdensity, and angular momentum.
Angular momentum has a lognormal distribution (\citep{Vitviska2002,Romanowsky2012}), and the overdensity has a Gaussian distribution. By means of Monte Carlo simulations we can get the 68\% confidence region. Another way, is to get $R_t$ in terms of the overdensity, since the angular momentum is proportional to $M^{5/3}$ (\citep{Ryden1987,Ryden1988,Catelan1996a,Romanowsky2012}, and $M=\frac{4 \pi}{3} R^3 \overline{\rho} (1+\delta)$. Recalling that $\delta$ has a Gaussian distribution, we may get the confidence region by means of the monte carlo method.
}
The plot shows that the $\Lambda$CDM model prediction is similar to that of the MG theories. Similarly to \citep{DelPopolo2019}, in which we showed that the mass-temperature relation {is not a good probe to test gravity theories beyond GR}, and to disentangle between GR, and $f(R)$ theories, here we have a similar result: the TAR is not a good probe to disentangle between GR, and $f(R)$ theories. In Fig. 5, we compare the predictions of some of the quintessence DE models previously cited
to the same $f(R)$ models of \citep{Lopes2018} plotted in Fig. 4. In the left panel of Fig. 5, we plot the predictions of the DE models for the TAR, obtained with the ESCM. From top to bottom, the cyan, blue, brown, magenta, black, red, and green lines represent the INV1, INV2, SUGRA, $w_{09}$\footnote{Namely the model having $w=-0.9$}, AS, $\Lambda$CDM without shear, and rotation, and $\Lambda$CDM with shear, and rotation, respectively.

In the right panel, we compare the \citep{Lopes2018} prediction with that of the DE models, plotted in the left panel. In order to have a more readable figure, we plotted just the INV2 (blue dashed line), and the $\Lambda$CDM with shear, and rotation (green dashed line) plotted in the left panel. These two lines contain all the other DE models (except INV1). The solid lines are the $\Lambda$CDM, and $f(R)$ models of \citep{Lopes2018}.
Similarly to what happens in Fig. 4, the DE models have similar predictions for the TAR to that of \citep{Lopes2018}, showing that the quoted radius cannot be used to disentangle DE, and $f(R)$ theories predictions.

\begin{table}
\caption{The ruled out ranges of $w$ based on the ESCM model.}
 \label{table1}
 \begin{tabular}{@{}lccc}
  \hline
  Stable structure & $M$ ($M_{\odot}$) & $R$ (Mpc) & range of $w$ ruled out \\
  \hline
  M31/MW & $(2.4 \pm 1.0) \times 10^{12}$ & $0.25 \pm 0.02$ & $ w >  -2.5$ \\
  M81 & $(8.9 \pm 2.2) \times 10^{11}$ & $0.90^{+0.06}_{-0.05}$ & $w \ge -1$ \\
  IC342 & $2.0^{+1.2}_{-1.0} \times 10^{11}$ & $0.51^{+0.03}_{-0.05}$ & $w \ge -1$ \\
  NGC253 & $2.34^{+0.99}_{-1.34} \times 10^{11}$ & $0.70 \pm 0.10$ & $w \ge -0.5$ \\
  CenA/M83 & $2.00^{+0.42}_{-0.33} \times 10^{12}$ & $1.37^{+0.08}_{-0.11}$ & $w \ge -1$ \\
  Local Group & $3.14^{+0.57}_{-0.54} \times 10^{12}$ & $0.96 \pm 0.03$ & $w \ge -2.5$ \\
  Fornax-Eridanus & $1.92 \times 10^{14}$ & $4.60^{+1.0}_{-0.72}$ & $w \ge -1.5$ \\
  Virgo & $1.67^{+0.31}_{-0.35} \times 10^{15}$ & $6.5^{+1.0}_{-1.5}$ & $w \ge -2.5$ \\
  \hline
 \end{tabular}
\end{table}

\section{Conclusions}\label{sect:conclusions}

In the present work, we discussed how shear, and rotation change the TAR, and
some of the parameters of the SCM. We used an extended SCM taking into account the effects of shear and vorticity to determine the $R_{\rm t}-M_{\rm t}$, in $\Lambda$CDM, and in DE scenarios. From the condition that the radial acceleration is zero, we got a formula for the maximum TAR more general than those found in literature. Shear, and rotation reduce the value of the TAR, particularly at the galactic scales where rotation is larger. As a consequence, using the $R_{\rm t}-M_{\rm t}$ relationship, and data from stable structures to constrain the $w$ parameter, one gets smaller values (absolute value) of the quoted parameter for structures having masses smaller than $10^{13} M_{\odot}$.
We also compared the $R_{\rm t}-M_{\rm t}$ relationship obtained for $\Lambda$CDM, and DE scenarios
with the prediction of the $f(R)$ theories, calculated by \citep{Lopes2018}. The result of the comparison shows that the $R_{\rm t}-M_{\rm t}$ relationship in the $f(R)$ models are practically identical to that of
the $\Lambda$CDM model, and DE scenarios. This implies that the $R_{\rm t}-M_{\rm t}$ relationship is not a good probe to disentangle between GR, and DE models predictions.