\section{Introduction}

The phase space of a dynamical system represents all its possible states. Starting from an initial condition, trajectories along the phase space points represent an evolution of the system in time \cite{lichtenberg,ozorio}.

In particular, for two degrees of freedom systems, the time progression of the flow of trajectories in phase space can be studied within a bidimensional surface as the phase paths intersect it successively \cite{lichtenberg}. These surfaces are known as Poincaré sections and can be generalized for a $2N$ degree of freedom system. The symplectic maps, explained in detail in the following section, can be approximately depicted as Poincaré sections of two degrees of freedom quasi-integrable systems.

This article is for students of physics, with the goal of inviting them to delve deeper into the subject of area-preserving maps by introducing fundamental concepts of symplectic maps using two examples: the standard and standard non-twist maps. They are paradigmatic systems, serving the purpose of introducing symplectic maps.

The standard map is a family of area-preserving maps with a single parameter $k$ that expresses the disturbance of the system. It was introduced by Chirikov \cite{chirikov} establishing its universality and numerous applications, as the description of the dynamics of magnetic field lines, and independently by Taylor \cite{taylor}, modelling magnetic fields in plasmas. Another significant result for the standard map was obtained by Greene \cite{greene1979}, involving a method to determine the transition to chaos.

First introduced by del Castillo Negrete and Morrison in 1993 \cite{del-Castillo-Negrete1993} while studying the chaotic transport in shear flow and described with further detail by Wurm, Apte, Fuchss, and Morrison in 2005 \cite{wurm2005}, the standard non-twist map is defined as a map that locally violates the twist condition.

This violation of the twist condition leads to some unique properties of physics interest, concerning transport barriers and the map's topology, such as the twin islands chain appearance and the reconnection process between others. The curve where the violation happens, the shearless curve, is destroyed for specific sets of parameters. This set can be obtained with the standard non-twist map breakup diagrams, where each point represents a pair of parameters.

The standard and standard non-twist maps can be applied to model a large variety of physical theories, e.g., magnetic fields in tokamaks \cite{portela2008, morrison2000, oda1995, caldas2012} and stellarators \cite{hayashi1995} (plasma physics), planetary orbits \cite{kyner1968,moser2001} and stellar pulsations \cite{munteanu2002} (astronomy), sheared flows \cite{del-Castillo-Negrete1993, del-Castillo-Negrete2000}.

This paper is organized as follows. In section 2, we start with a review of the main concepts of quasi-integrable conservative systems in order to introduce the symplectic maps. Sections 3 and 4 are dedicated to the example maps, exploring the general properties of the symplectic maps as well as some unique from the standard and standard non-twist maps. Section 5 presents the parameter space as a tool to investigate the shearless curve breakup.

\section{Symplectic Maps}

Essentially, a dynamical system can be described by its phase space using the Hamiltonian framework. For example, a $4$-dimensional phase space will be formed of the generalized positions $\{q_1,q_2\}$ and momentum $\{p_1,p_2\}$ showing all possible states of the system. These two degrees of freedom \emph{Hamiltonian dynamical systems} are fundamentally symplectic\footnote{The term symplectic comes from the greek word \emph{sumplektikos} and it means intertwined.} \cite{abud2013}, this is, they preserve both volume and infinitesimal areas in phase space, according to the \emph{Liouville's Theorem} \cite{meiss1992}.

Dynamical systems with two degrees of freedom may become hard to deal with, given their four-dimensional phase space. In order to simplify the analysis of the phase spaces, Poincaré introduced the concept of \emph{Poincaré sections}, as they are known nowadays. These sections are cross-sectional surfaces to the flux of trajectories in phase space \cite{poincare1893}. Thereby, for a time-independent Hamiltonian system of two degrees of freedom, the Poincaré sections allow us to simplify the problem, returning us to a two-dimensional system.

\begin{figure}[htb]
    \centering
    \includegraphics[width=6.82cm]{draw.png}
    \caption{Phase diagrams passing through a 2-dimension Poincaré Section $\sigma$ and its iterated intersection points $\mathcal{P}$.}
    \label{draw}
\end{figure}

For the nth intersection of the trajectory of a dynamical system with a Poincaré section in phase space being $\mathcal{P}_n$, the next point of intersection given the time evolution will be $\mathcal{P}_{n+1}$. We define a \emph{Poincaré map} as the discrete dynamics given by a matrix $\mathcal{M}$ in such a way that $\mathcal{P}_{n+1} = \mathcal{M}\mathcal{P}_n$ \cite{lichtenberg}, being each point $\mathcal{P}$ placed only in a single direction. Given this framework, we define a discrete \emph{orbit} as a sequence of points $\{\mathcal{P}_0,\dots,\mathcal{P}_n,\mathcal{P}_{n+1},\dots\}$ such that $\mathcal{P}_{n+1} = \mathcal{M}\mathcal{P}_n$ \cite{meiss1992}.

An indispensable concept for the understanding of the symplectic maps is the \emph{integrability}. A two-degrees-of-freedom system is said integrable if there exists exactly two independent \emph{integrals of motion} \cite{lichtenberg,reichl} (essentially the degrees of freedom of the system must match the number of integrals of motion). These so-called integrals of motion are functions $\mathcal{F}(q_n,p_n), \ n=1,2$ that are constants along a trajectory in phase space \cite{nivaldo}.

Let's now consider the following perturbed Hamiltonian \cite{lichtenberg}
\begin{equation}
\label{Perturbed Hamiltonian}
    \mathcal{H}(\theta,\mathcal{J}) = \mathcal{H}_0(\mathcal{J}) + \epsilon\mathcal{H}_1(\theta,\mathcal{J}),
\end{equation}
where ($\theta$,$\mathcal{J}$) represents the action-angle variables of the integrable Hamiltonian $\mathcal{H}_0$\footnote{The nth action variable is defined as $\mathcal{J}_n = \frac{1}{2\pi}\oint p dq, \ n = 1,\dots,N$ and the nth angular variable as $\theta_n = \pdv{W}{\mathcal{J}_{n}}, \ n = 1,\dots,N$ with the line integral along one period and $W$ being the complete integral of the Hamiltonian system \cite{nivaldo}}. The term $\mathcal{H}_1$ expresses the non-integrable part of the Hamiltonian and $\epsilon$ is the amplitude of perturbation.

From Kolmogorov-Arnold-Moser's (KAM) theorem, for frequencies $\Omega = \pdv{\mathcal{H}_0}{\mathcal{J}}$ sufficiently irrational and for $\epsilon\ll 1$, most of the invariant tori will become slightly deformed but preserving their former properties \cite{lichtenberg,reichl,ozorio}. In other words, for small perturbations, the system has just a slight change in most of its phase space. However, this kind of system has also chaotic orbits around the so-called hyperbolic points, which will be described later.

For these quasi-integrable systems, the regular orbits of irrational frequencies are presented in Poincaré sections by curves or \emph{invariant tori} \cite{lichtenberg} (topologically they match with the definition of a torus). These curves are usually referred to as orbits or simply tori. However, in the limit $\epsilon\gg1$ (overwhelming KAM's theorem validity), more tori become broken implying global chaos. Their presentation in Poincaré sections is just of densely distributed points.

In the same way, the KAM theorem states the quasi-integrability for irrational frequencies, but it fails to provide the same conclusion for the rational ones. Complementing KAM's theorem, \emph{Poincaré-Birkhoff theorem} states that, considering a torus with rational frequency $\Omega = r/s, \ r,s\in \mathbb{Z}$, from the perturbation $H_1$ creates an even number $2ms,\ m\in \mathbb{N}$, of periodic points. From these points, half of them are stable, expressing the labeled \emph{elliptical points} and the other half unstable, called \emph{hyperbolic points} \cite{birkhoff1913,aguiar}. Both of them will be described in the following section.

A \emph{periodic orbit} of period $s$ is defined as $\mathcal{M}^s(\theta_n,\mathcal{J}_n) = (\theta_{n} + r,\mathcal{J}_{n}), \ \forall\ n$, as a consequence of this statement, their frequencies are rational numbers $\Omega = r/s, \ r,s\in \mathbb{Z}$ \cite{lichtenberg}. Different from the quasi-periodic orbits, from Poincaré-Birkhoff's theorem the rational frequency implies periodic distributions of the iterated points, presenting themselves in Poincaré sections in closed ellipses.

The time evolution in \emph{quasi-integrable conservative dynamical systems} can be discrete (map) or continuous (differential equations). The first one gives rise to the \emph{symplectic maps}, defined as area-preservative conservative systems, described by Poincaré sections of time-independent bidimensional Hamiltonians \cite{lichtenberg}, given by recurrence relations of the approximated form \cite{lichtenberg}:
\begin{equation}
\label{General representation of the iteration form of symplectic maps}
\begin{cases}
    \mathcal{J}_{n+1} = \mathcal{J}_{n} + \epsilon f(\theta_n,\mathcal{J}_{n+1}) \\
    \theta_{n+1} = \theta_{n} + \Omega(\mathcal{J}_{n+1}) + \epsilon g(\theta_n,\mathcal{J}_{n+1})
\end{cases}
\end{equation}
for $n \in \mathbb{N}$ being the nth iteration, $f$ and $g$ analytic periodic functions in $\theta$ and $\Omega$ the unperturbed frequency.

\begin{figure*}[htb]
    \centering
    \begin{subfigure}[b]{.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{k10-4}
        \caption{$k = 1\times 10^{-4}$.}
        \label{k = 10(-4)}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{k0.0318}
        \caption{$k = 0.0318$.}
        \label{k = 0.0318}
    \end{subfigure}
    \vskip\baselineskip
    \begin{subfigure}[b]{.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{k0.151}
        \caption{$k = 0.151$.}
        \label{k = 0.151}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.48\textwidth}
        \centering
        \includegraphics[width=\textwidth]{k4}
        \caption{$k = 4$.}
        \label{k = 4}
    \end{subfigure}
    \caption{Example of a sequence of standard maps increasing the value of the disturbance parameter $k\geq 0$.}
    \label{Standard Map: maps}
\end{figure*}

Since these maps are symplectic, the \emph{area-preserving} property must hold. For a standard symplectic map $(\mathcal{J}_{n+1},\theta_{n+1}) = \mathcal{M}(\mathcal{J}_{n},\theta_n)$ (matrix form) this characteristic is expressed as the determinant of its Jacobian matrix being unitary \cite{meiss1992}.

\begin{equation}
\mathrm{det} \left[\pdv{(\mathcal{J}_{n+1},\theta_{n+1})}{(\mathcal{J}_{n},\theta_n)}\right] = 1
\end{equation}

This is an important statement about these maps. Since it is a unique property of them, the symplectic maps can also be called \emph{area-preserving} \emph{maps}.

In the following sections, we will work with two examples of symplectic maps, the standard and the standard non-twist maps.

\section{The Standard Map}

Also known as  \emph{Chirikov-Taylor's map} the standard map is described by the recurrence equations \cite{chirikov,taylor},
\begin{equation}
\label{Standard Map Equations}
\begin{cases}
    y_{n+1} = y_n + k\sin(2\pi x_n)\\
    x_{n+1} = x_n - y_{n+1},\ \ \ \ \ \ \ \ \ \ \ \ \  \mbox{mod } 1
\end{cases}
\end{equation}
in which $x_n$ and $y_n$ represents the nth iteration of the variables $\theta_n$ and $\mathcal{J}_n$ respectively and $k\geq 0$ is a perturbation parameter. The value of $k$ is directly connected with the constant $\epsilon$ in Equation (\ref{Perturbed Hamiltonian}) and therefore as we increase it, the non-linearity of the

\begin{figure}[htb]
 \centering
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{RN_k0.0318}
     \caption{Winding number profile for $k = 0.0318$.}
     \label{RN k = 0.0318}
 \end{subfigure}
 \hfill
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{RN_k0.151}
     \caption{Winding number profile for $k = 0.151$.}
     \label{RN k = 0.151}
 \end{subfigure}
 \caption{Examples of winding number profiles increasing the disturbance parameter $k$.}
\label{Standard Map: winding number profiles}
\end{figure}

\noindent map raises, revealing eventually chaotic structures \cite{chirikov}.

The system is a twist map since it satisfies the \emph{twist condition} \cite{lichtenberg,reichl}:
\begin{equation}
    \label{Twist-Condition}
    \pdv{x_{n+1}}{y_n} \neq 0
\end{equation}
everywhere. This circumstance is the analog of the \emph{nondegeneracy condition} from Hamiltonian systems \cite{del-Castillo-Negrete1996}:
\begin{equation}
    \label{Non-Degeneracy Condition}
    \pdv[2]{H_0}{\mathcal{J}} \neq 0
\end{equation}
and it's a crucial statement for KAM's theorem applicability to the map.

From Fig. (\ref{k = 10(-4)}) we can see the standard map for the small value of the parameter $k = 1\times 10^{-4}$. This choice makes the non-linear term $\sin(2\pi x_n)$ a diminutive contribution to the Equations (\ref{Standard Map Equations}). This implies a near-linear behavior on the map thus making the phase space with predominant regular orbits.

The almost horizontal lines in Fig. (\ref{k = 10(-4)}) are called \emph{invariant curves} or \emph{KAM curves} and each one of them represents quasi-periodic orbits in phase space. Each invariant curve is related to an initial condition $(x_0,y_0)$, in Fig. (\ref{k = 10(-4)}) only a few invariants are shown in order to keep the image clean.

If we increase the magnitude of the parameter to $k=0.0318$, as we can see in Fig. (\ref{k = 0.0318}), the non-linearity starts to take place. The invariant curves start to deform as stated by the KAM's theorem. It is also very prominent in the same figure the arising of \emph{islands} on the map that represent the periodic orbits in phase space, as predicted by the Poincaré-Birkhoff theorem. The central point of an island is a periodic elliptical point, surrounded by stable quasi-periodic motion around it \cite{reichl,greene1968}. The points on both sides of the islands in Fig.(\ref{k = 0.0318}) are the periodic  hyperbolic points, being unstable, with a small chaotic region around them \cite{abud2013,greene1968}.

Another important concept to extract from Fig. (\ref{k = 0.0318}) is the \emph{separatrix} curve. The separatrices are the invariant lines that almost connect the hyperbolic points. These curves, for quasi-integrable maps, are along the hyperbolic points, marking a boundary between two different behaviors from the phase space paths. In fact, in Fig. (\ref{k = 0.0318}), the separatrix traces the boundary between the invariant curves inside the islands and others outside.

For $k = 0.151$ in Fig. (\ref{k = 0.151}), a large number of islands appear, revealing \emph{chaotic regions} and thus the non-integrability of the system. For a high value of the parameter, as $k = 4$ for example (Fig. (\ref{k = 4})), the map gets completely chaotic, showing no sign of KAM curves and islands. The chaotic regions

\begin{figure}[htb]
 \centering
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{k0.509_redsquare}
     \caption{$k = 0.509$.}
     \label{k = 0.509}
 \end{subfigure}
 \hfill
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{k0.509_zoom}
     \caption{$k = 0.509$ map with zoom.}
     \label{k = 0.509 (zoom)}
 \end{subfigure}
 \caption{Example of the standard map fractal behavior. The red rectangle indicated in (\ref{k = 0.509}) is shown amplified in Fig. (\ref{k = 0.509 (zoom)}).}
\label{Standard Map: fractal behavior}
\end{figure}

\noindent shown in Figs. (\ref{k = 0.151}) and (\ref{k = 4}) originates from the hyperbolic points, having an unpredictable behavior in the long run due to the high sensitivity to the initial conditions. During the transition from Fig. (\ref{k = 0.151}) to (\ref{k = 4}), there is a critical value of $k$ from which the chaotic regions are no longer separated. The numeric calculation of this threshold value was made for the first time by John Greene \cite{greene1979}.

Chaotic maps can also have their topology and monotonicity obtained from \emph{winding number profiles}. These profiles are based on the definition of \emph{winding number} \cite{meiss1992,reichl} (also called \emph{rotation number})
\begin{equation}
\label{winding Number Definition}
    \omega = \lim_{n\to\infty}\frac{x_{n+1} - x_0}{n},
\end{equation}
provided the limit exists. Invariant curves and islands have a dense behavior as we can see in the following figures, thus having an existing winding number corresponding to each one of them.

For $\epsilon = 0$, invariant curves (quasi-periodic orbits) have irrational winding numbers, while islands (periodic orbits) have their rotation numbers given by a rational number $\omega = r/s, \ r,s \in \mathbb{Z}$ (where $r$ is the number of revolutions in $\theta$ and $s$ is the period) \cite{lichtenberg}.

The winding number profiles are $\omega\times y$ graphics that show how the winding numbers within an arbitrary curve act with respect to the coordinates of the vertical $y$-axis (For the standard map we will always be using the line $x = 0.5$ as the reference curve.).

Looking at the winding number profile in Fig. (\ref{RN k = 0.0318}), by comparing the major island interval in the $y$-axis with the correspondent one in the map (\ref{k = 0.0318}), it becomes clear that the larger island in the middle has a constant winding number, forming the central plateau. Thus, having the winding number profile of a map enable us to identify islands along the chosen curve without looking at the map itself and obtain its rotation number.

 As we can see in Fig. (\ref{RN k = 0.151}), a large amount of plateaus in the figure actually represents the numerous islands from the map (Fig. (\ref{k = 0.151})). However, a new element presents in the graphic, the irregular (clearly not continuous) dark oscillating lines in the winding number profile which can be identified as places where the winding number does not exist (as the limits from Equation (\ref{winding Number Definition})). If we compare the $y$-axis intervals of these regions within the reference curve $x=0$ with the ones from its correspondent map, these regions evidently are chaotic. Thus from the winding number profile, irregular regions (places where $\omega$ does not exist) as the ones from Fig. (\ref{RN k = 0.151}) represent

\begin{figure*}[htb]
 \centering
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{a0.35_b0.35_blueline}
     \caption{$a = 0.35$ and $b = 0.35$.}
     \label{a = 0.35}
 \end{subfigure}
 \hfill
 \begin{subfigure}[b]{.48\textwidth}
     \centering
     \includegraphics[width=\textwidth]{RN_a0.35.png}
     \caption{Winding number profile for a = 0.35 and b = 0.35.}
     \label{RN a = 0.35}
 \end{subfigure}
 \caption{The standard non-twist map and its nonmonotonic winding number profile as an example.}
\label{Non-Twist Map: general example}
\end{figure*}

\noindent chaotic regions on the map.

In Fig. (\ref{Standard Map: winding number profiles}) the \emph{monotonicity} of the standard map profiles is evident. A curve is said monotonic if is either strictly increasing or decreasing in an interval, thus matching the behavior seen in Fig. (\ref{Standard Map: winding number profiles}).

The standard map has an interesting property, it has a fractal behavior \cite{lichtenberg}. Generally, fractals are geometric patterns that keep their topological structures for arbitrarily small scales \cite{mandelbrot} In Fig. (\ref{k = 0.509}) (standard map for $k = 0.509$), if we choose to take a closer look inside the red rectangle, we'll get Fig. (\ref{k = 0.509 (zoom)}).

This smaller portion of Fig. (\ref{k = 0.509}) has its own central island with a chain of smaller ones around it, similar to its larger-scale version (Fig. (\ref{k = 0.509 (zoom)})). The second map figure thus represents analogous information in comparison with the first one, showing the fractal behavior of the standard map. Further successive amplifications would confirm the fractal structure \cite{lichtenberg}.

\section{The Non-Twist Map}

The so-called \emph{Standard Non-Twist map} is defined by the following recurrence equations \cite{del-Castillo-Negrete1996}.
\begin{equation}
\label{Non-Twist Map Equations}
\begin{cases}
    y_{n+1} = y_n - b\sin(2\pi x_n)\\
    x_{n+1} = x_n + a(1 - y_{n+1}^2) \ \ \ \ \  \mbox{mod } 1,
\end{cases}
\end{equation}
in which the variables $x$ and $y$, as in the Equations (\ref{Standard Map Equations}), represent the action-angle variables $\theta$ and $\mathcal{J}$ respectively. The two parameters $a,b\in \mathbb{R}$ regulate the map topology.

A non-twist map is defined as a map that violates the twist condition (Equation (\ref{Twist-Condition})) at least in one invariant curve \cite{wurm2005,del-Castillo-Negrete1996,carvalho1992, corso1998}. In this section, we will only work with the standard non-twist map, in which the twist condition is violated (even for $\epsilon = 0$), at $y=0$.

Distinguishing from the standard map, the standard non-twist map has in the definition Equations (\ref{Non-Twist Map Equations}) a quadratic function $a(1-y_{n+1}^2)$. For $b=0$ the derivative of the quadratic function with respect to $x$ is known as the \emph{shear function} and it provides an extremum point on the winding number profile. The value of $a$, thus, has a direct influence on the nonmonotonicity of the winding number profile for the non-twist map.

In Fig. (\ref{a = 0.35}) we have the standard non-twist map for the parameters $a = b = 0.35$. The curve referent to the invariant that violates the twist condition is called \emph{shearless curve} \cite{wurm2005,del-Castillo-Negrete1996,carvalho2013} and it can be seen enhanced in red in Fig. (\ref{a = 0.35}). In the same figure, this curve is centrally located on the map. As we shall see, its position can be comprehended from its winding number profile.

Figure (\ref{a = 0.35}) also shows hyperbolic points that connect the periodic orbits from the \emph{twin chains of islands} in the center region. The twin chains are the islands connected by separatrices in Fig. (\ref{a = 0.35}) on both sides of the shearless curve \cite{wurm2005}.

Figure (\ref{RN a = 0.35}) has the winding number profile of Fig. (\ref{a = 0.35}), being the reference curve the dashed blue line $x = 0$. Each point in a winding number profile shows the winding number of a specific orbit, the red one referent to the shearless curve. From the same figure, it becomes clear, in this case, that the shearless curve point is related to the maximum of the profile.

The fact that the shearless curve refers to the maximum point in Fig. (\ref{RN a = 0.35}) is a direct consequence of the twist condition (Equation (\ref{Twist-Condition})) violation, since the winding number's derivative is
\begin{equation*}
    \pdv{\omega}{y_n} = 0
\end{equation*}
 there. From calculus, we know that if the slope of a continuous function is zero, there is an extremum or inflection on the curve of the function at that point. Since an extremum can be a local maximum or a local minimum, both possibilities represent a shearless curve point in the winding number profile.

In Fig. (\ref{RN a = 0.35}) we can also notice that the winding number profile is clearly nonmonotonic, as we expected. Non-twist maps have a couple of remarkable properties, different from the standard map ones. In the following subsections, we shall explore some of them.

\subsection{Even Reconnection}

When we have a winding number profile with an extremum, the shearless curve point, there will be two orbits with the same winding number, since the curve is nonmonotonic. This fact is a unique property from the non-twist maps since it originates from the violation of the twist condition (Equation (\ref{Twist-Condition})) and it is known as \emph{degeneracy}.

If the winding number profile starts to decrease its extremum vertically for a change in parameters $a$ and $b$, becoming flat, the periodic orbits at some point will have the same rotation number (a plateau) \cite{wurm2005}. This sequence of events is called \emph{periodic orbit's collision} and it is followed by its annihilation.

The periodic orbit's collision involves another process, the reconnection of the separatrices. This reconnection is a process of global bifurcation (large stability changes in an equilibrium \cite{yorke}) that changes the topology of the separatrices, during this event the separatrices collide \cite{wurm2005}.

The way the reconnection happens depends on the parity of the orbit's period that can be obtained from the number of islands on the map. Notice that the maps are periodic in the $x$-axis, so half islands on it are the same ones from the opposite side. In this section, we take a look at the even-period scenario.

Fig. (\ref{even reconnection}) shows an example of an even reconnection sequence. Initially, in Fig. (\ref{a = 0.532}), we have two twin island chains. Decreasing the parameter $a$ give us Fig. (\ref{a = 0.505}). Now the twin island chains start to merge in a dipole topology \cite{scholarpediaWurm}, the separatrices get closer to each other, becoming one and the islands and hyperbolic points collide. In the same image, the shearless curve vanished during the reconnection process.

Decreasing again the value of $a$ even more, we see that the islands finally collided, vanishing the twin island chains (Fig. (\ref{a = 0.495})). The center has

\begin{figure*}[t]
    \centering
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.532_b0.28.png}
        \caption{$a = 0.532$ $b = 0.28$.}
        \label{a = 0.532}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.505_b0.28.png}
        \caption{$a = 0.505$ $b = 0.28$.}
        \label{a = 0.505}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.495_b0.28.png}
        \caption{$a = 0.495$ $b = 0.28$.}
        \label{a = 0.495}
    \end{subfigure}
    \vskip\baselineskip
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.532.png}
        \caption{Winding number profile for \\ $a = 0.532$ $b = 0.28$}
        \label{RN a = 0.532}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.505.png}
        \caption{Winding number profile for \\ $a = 0.505$ $b = 0.28$}
        \label{RN a = 0.505}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.495.png}
        \caption{Winding number profile for \\ $a = 0.495$ $b = 0.28$}
        \label{RN a = 0.495}
    \end{subfigure}
    \caption{Example of even period reconnection sequence of maps (upper row $y$ vs $x$) together with their respective winding number profiles (upper row $\omega$ vs $y$), referenced as WNP.}
    \label{even reconnection}
\end{figure*}

 \noindent become indistinguishable compared with the first image and the shearless curve has reappeared.

 From a winding number profile perspective, the  maximum point from the profile in Fig. (\ref{RN a = 0.532}) starts to lower with the decrease of the parameter $a$ as we can see from Fig. (\ref{RN a = 0.532}) to Fig. (\ref{RN a = 0.505}), becoming flat (no extremum, implying the absence of a shearless curve). The last profile Fig. (\ref{RN a = 0.495}), with the further decrease of $a$, shows a new maximum showing the shearless curve on the map once more.

 \subsection{Odd Reconnection}

In a correspondingly way, for an odd period scenario, the reconnection process occurs. However, its development is different from the even period's case as we shall see. As it is evident from Fig. (\ref{a = 0.346}), we have an odd number of islands in each chain, setting an odd reconnection scenario. We decrease the parameter $a$ originating Fig. (\ref{a = 0.3407}). Now the shearless curve has already vanished and the separatrices have connected themselves originating a contour-invariant curve around the islands.

If we further decrease the value of $a$, the separatrices from Fig. (\ref{a = 0.3407}) take a meandering shape around the islands until they are suppressed as in Fig. (\ref{a = 0.3382}). The separatrices, forming a \emph{homoclinic} behavior (joins saddle equilibrium points) in the middle, are then replaced by the \emph{heteroclinic} (connects different equilibrium points) curve exhibited in Fig. (\ref{a = 0.3382}) in which the shearless curve shows up again.

From a winding number profile perspective, similarly to the even reconnection, from Fig. (\ref{RN a = 0.346}) to Fig. (\ref{RN a = 0.3407}) the maximum point becomes a plateau, and in the subsequent step (Fig. (\ref{RN a = 0.3382})) the reinstate of the shearless curve comes from the arising of a new maximum.

\begin{figure*}[t]
    \centering
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.346_b0.30.png}
        \caption{$a = 0.346\ b = 0.30$.}
        \label{a = 0.346}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.3407_b0.30.png}
        \caption{$a = 0.3407\ b = 0.30$.}
        \label{a = 0.3407}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{a0.3355_b0.30.png}
        \caption{$a = 0.3355\ b = 0.30$.}
        \label{a = 0.3382}
    \end{subfigure}
    \vskip\baselineskip
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.346.png}
        \caption{Winding number profile for $a = 0.346$\ $b = 0.30$}
        \label{RN a = 0.346}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.3407.png}
        \caption{Winding number profile for $a = 0.3407$ $b = 0.30$}
        \label{RN a = 0.3407}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{.32\textwidth}
        \centering
        \includegraphics[width=\textwidth]{RN_a0.3355.png}
        \caption{Winding number profile for $a = 0.3355$ $b = 0.30$}
        \label{RN a = 0.3382}
    \end{subfigure}
    \caption{Example of odd period reconnection sequence of maps (upper row $y$ vs $x$) together with their respective winding number profiles (upper row $\omega$ vs $y$), referenced as WNP.}
    \label{odd reconnection}
\end{figure*}

\begin{figure*}[t]
\centering
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[width=\textwidth]{a0.354_b0.4}
    \caption{$a = 0.354$ $b = 0.4$.}
    \label{b = 0.4}
\end{subfigure}
\hfill
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[width=\textwidth]{a0.354_b0.56}
    \caption{$a = 0.354$ $b = 0.56$.}
    \label{b = 0.56}
\end{subfigure}
\vskip\baselineskip
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[width=\textwidth]{a0.354_b0.7}
    \caption{$a = 0.354$ $b = 0.7$.}
    \label{b = 0.7}
\end{subfigure}
\hfill
\begin{subfigure}[b]{.49\textwidth}
    \centering
    \includegraphics[width=\textwidth]{a0.354_b0.92}
    \caption{$a = 0.354$ $b = 0.92$}
    \label{b = 0.92}
\end{subfigure}
\caption{Sequence of maps serving as an example of the shearless curve destruction process.}
\label{shearless curve destruction}
\end{figure*}

\subsection{Shearless Curve Destruction}

Looking at the standard non-twist map definition, Equations (\ref{Non-Twist Map Equations}), the perturbation term $b$ evidently has an important role in the disturbance of the map. If we input a large value for $b$, the number of deformed and destroyed tori will increase expressing regular and chaotic regions \cite{del-Castillo-Negrete1996,del-Castillo-Negrete1997}. Eventually, with the increase of the perturbation in the map, the shearless curve can be destroyed.

In this section, we shall dissect this Shearless curve breakup process step by step. Starting from similar parameters in comparison with Fig. (\ref{a = 0.346}), we have Fig. (\ref{b = 0.4}). Now, instead of changing $a$ as we did to show the odd reconnection scenario, in order to grow disturbance in the map, we increase $b$. Within an odd reconnection process, Fig. (\ref{b = 0.56}) shows the standard non-twist map for a larger parameter $b$. On the image, we can see that the chaotic regions of the map have grown and the KAM curves and periodic orbits within the separatrices start to reconnect.

Increasing the perturbation parameter to $b = 0.7$ brings to Fig. (\ref{b = 0.7}). From this point, there are no longer invariant curves, islands, and separatrices, instead, we have a darker region of \emph{stickiness}, around the still-standing shearless curve. Stickiness is the temporary confinement of chaotic orbits in a specific region of the phase space before they diffuse to a larger region after, usually, a large number of iterations \cite{contopoulos2010,szezech2009}.

Even with the vanishment of the majority of the KAM tori in the center region of the map (\ref{b = 0.7}), the shearless curve is kept intact. This illustrates the shearless curve's resilience to perturbation, being it the most resistant torus \cite{abud2013,del-Castillo-Negrete1996,scholarpediaWurm}. In the rupture limit to global chaos, the shearless torus acts as a \emph{transport barrier}. This name comes from one of its traits: the barrier prevents transport between the chaotic regions \cite{szezech2009,fonseca2014}. Even after the shearless curve destruction, the stickiness remains, keeping the transport barrier a little longer \cite{szezech2009}.

In the last image of the process (Fig. (\ref{b = 0.92})), the chaotic sea takes over the last torus (the shearless curve) and the darker stickiness region becomes homogeneous chaos. This expresses the shearless curve destruction, also known as \emph{Transition to Global Chaos in Non-Twist Maps} \cite{del-Castillo-Negrete1996}. Thus by increasing the parameter $b$, the perturbation on the map increases. Since this parameter is directed related to the disturbance parameter $\epsilon$ in Equation (\ref{Perturbed Hamiltonian}), for $b\gg 1$, KAM's theorem implies a larger number of chaotic tori, the same behavior exposed in Fig. (\ref{odd reconnection}).

\section{Parameter Space}

The shearless breakup shown in Fig. (\ref{shearless curve destruction}) occurs for a specific set of parameters $(a,b)$. In order to understand which ones lead to this scenario, there is a method proposed by Shinohara and Aizawa \cite{shinohara1997}, using the \emph{indicator points} to create a \emph{parameter space} or \emph{shearless breakup diagram}. Indicator points can be generally described as a set of four points
\begin{equation}
    \label{Indicator Points}
    z_1^{(\pm)} = \left(\pm \frac{1}{4},\pm \frac{b}{2}\right), \ \ z_2^{(\pm)} = \left(\frac{a}{2}\pm\frac{1}{4},0\right),
\end{equation}
that necessarily are contained in the shearless curve \cite{shinohara1998}, if it exists.

The parameter space of the shearless curve

\begin{figure}[htb]
    \centering
    \includegraphics[width=0.7\textwidth]{paramSpace.png}
    \caption{Shearless breakup diagram for a mesh of $2000\times 2000$ points, intervals $0<a<1$ and $0<b<1$, threshold range $\abs{y} < 10$ and $10^6$ iterations.}
    \label{Parameter Space}
\end{figure}

\noindent consists of a mesh of points whose axis are the parameters $a$ (x-axis) and $b$ (y-axis), usually limited to the same interval. One way of building a parameter space is checking, for a large number of iterations if the $y$ value of the indicator points for a determined point $(a,b)$ is restrained to a specified threshold range.

If we set the points inside the threshold range to be colored red and the ones outside white, the resulting graphic is an approximation of the parameter space in the limited interval. For a mesh of $2000\times 2000$ points, the intervals $0<a<1$, $0<b<1$ and the threshold range $\abs{y} < 10$, Fig. (\ref{Parameter Space}) shows the resulting parameter space after $10^6$ iterations of the indicator points.

The white points (out of the threshold range) represent the parameters $(a,b)$ for which the shearless curve has been broken, while in the red ones it is maintained. Fig. (\ref{Parameter Space}) shows that the points for which the shearless curve does exist and the ones it does not are separated by a sharped boundary. Taking a closer look at Fig. (\ref{Parameter Space}) we can visualize fractal and continuous boundaries \cite{mathias2019,mugnaine2020}.

From Fig. (\ref{Parameter Space}), when $b$ increases, it is noticeable the shearless breakup tendency grows. However, this can be compensated by decreasing $a$. Thus, the parameter space acts as a useful tool in order to investigate which set of parameters $(a,b)$ preserves the shearless curve and to visualize how a change in the value of one of them or both influences the possibility of this scenario. The continuous boundaries of the parameter space are related to the collision of the twin islands \cite{wurm2005}.

\section{Conclusion}

In this review, we introduced the standard map and the non-twist standard map. These maps can be derived as approximations of Poincaré maps of Hamiltonian of quasi integrable systems of two degrees of freedom. Invariant curves, islands, and chaotic trajectories are introduced. The route to global chaos in phase space is shown by varying control parameters. Special features of the standard non-twist systems contained in the non-twist map are presented. Among these features, are the onset and break up of shearless invariants curves, as well as their dependence on the control parameters.

\section*{Acknowledgments}

We acknowledge the financial support from the scientific agencies: São Paulo Research Foundation (FAPESP) under Grant Nos. 2022/04251-7, 2022/05667-2 and 2018/03211–6, and Conselho Nacional de Desenvolvimento Cientíıfico e Tecnológico (CNPq) under Grant No. 304616/2021-4.

\printbibliography

\end{document}


\section{Introduction}
Observations of large-scale cosmic phenomena, galaxy clusters,  the matter power spectrum and cosmic microwave background radiation are evidence for the existence of non-baryonic matter that constitutes about 85\% of the total matter content of the Universe~\cite{Kadribasic:2017, Ade:2014}. Non-baryonic matter, known as dark matter, is one of the most fundamental topics in cosmology, astronomy, and high-energy physics~\cite{Cushman:2013,aghanim:2020,bertone:2005}. Most dark matter laboratories are designed and implemented based on direct and indirect detection approaches.
The indirect detection approach corresponds to the annihilation or decay of dark matter particles~\cite{Jungman:1996}. In contrast, direct detection is via the elastic collision of a dark matter particle with the nuclei of atoms~\cite{Gaitskell:2004} and generates typical recoil energy of
${\mathcal O({\mathrm{ keV}})}$~\cite{Goodman:1985,OHare:2015}. This energy operates to create an electron-hole pair inside the detector. In recent  decades, many studies for the direct search of dark matter were employed, such as XENON10~\cite{Xe10}, XENON100~\cite{Xe100:2012}, SIMPLE~\cite{Girard:2012}, CoGeNT~\cite{Aalseth:2013} and DAMA/LIBRA~\cite{Bernabei:2010, Bernabei:2018}.
In most direct dark matter experiments, a dark matter particle has been considered with a mass of ${\mathcal O(10-100)}$ GeV/c$^{2}$
~\cite{Cushman:2013}. Recently, with the development of convincing theoretical models, dark matter masses less than 10 GeV/c$^{2}$ have also been introduced~\cite{Kadribasic:2017, Cushman:2013, Romani:2018, Crisler:2018}.

In many experiments, the target material is either low-pressure gas or scintillating liquids. These experiments typically feature a rather high threshold energy for detecting dark matter based on nuclear  recoil~\cite{Kadribasic:2017, Mayet:2016}. Also, in a gaseous detector, the large volume of gas required is a deterrent. Due to the low energy of nuclear recoils for dark matter with low mass, high-mass detectors with low detection thresholds are more desirable \cite{Agnese:2014, Mirabolfathi:2017,Iyer:2021}. Semiconductor targets, such as silicon or germanium, can be utilized as cryogenic solid-state ionization detectors, where a sensitivity to single electron events has been demonstrated~\cite{Kadribasic:2017, Mirabolfathi:2017, Agnese:2018,Heikinheimo:2022}. When a particle colliding with a nucleus of the detector material transfers enough energy to the nucleus, it causes a defect in the configuration of the crystal lattice~\cite{Kadribasic:2017}. The effective threshold energy in generating the defect depends on the recoil angle in silicon and germanium crystals. Holmstrom et al.~\cite{Holms:2008, Holms:2010} obtained the threshold energy required to create a defect in germanium and silicon crystals for different directions using molecular dynamics simulations.  It has been argued~\cite{Heikinheimo:2019,Kadribasic:2017} that this threshold displacement energy is equivalent to the ionization threshold, which implies that the ionization signal will also be sensitive to the recoil direction.
Like most spiral galaxies, the current models assume that the Milky Way is immersed in a halo of dark matter. Inside the halo, the solar system moves toward the constellation Cygnus~\cite{OHare:2015}, and the direction of motion of dark matter particles in the lab frame is opposite to this. The Earth's rotation and motion around the Sun cause the laboratory frame velocity of the dark matter particles to change over time. Due to the anisotropy of the ionization threshold, this results in daily and annual modulation effects in the dark matter interaction rate~\cite{Heikinheimo:2019,sassi:2021,Coskuner:2019,Trickle:2019,Coskuner:2022}.

In this paper, we study the daily and annual rate of dark matter interaction with silicon nuclei, taking into account the energy threshold and the direction of the recoil induced by the collision with a dark matter particle. We use threshold energy data for silicon obtained via molecular dynamics simulations. We investigate the modulation of the ionization rate for dark matter interactions due to the directional dependence of the threshold energy and the motion of the laboratory frame with respect to the galactic rest frame.

The layout of this paper is as follows: in Section~\ref{sec2} we calculate the directional event rate of a dark matter in the presence of a direction-dependent energy threshold, and in Section~\ref{sec3} we present the results.

\section{Dark matter rate }\label{sec2}
To calculate the interaction rate of dark matter with the nuclei of the detecting material, we need the distribution function of the dark matter velocity in the galactic halo. We consider the Maxwell-Boltzmann distribution~\cite{lewin:1996, sassi:2021} for the velocity of dark matter in the halo as follows:
\begin{equation}
f(\textbf{v})=\frac{1}{N_{\rm esc}(2\pi\sigma_v^2)^\frac{3}{2}}\exp\left(-\frac{\textbf{v}^2}{2\sigma_v^2}\right)\Theta(v_{\rm esc}-|\textbf{v}|),
\end{equation}
where $v_{\rm esc}=544 $ kms$^{-1}$~\cite{smith:2007} is the escape velocity. The dark matter velocity dispersion $\sigma_v=\frac{v_0}{\sqrt{2}}$, where $v_0=220$ kms$^{-1}$ is the local circular speed~\cite{kerr:1986}. The normalization constant ${N_{\rm esc}}$ is given by
\begin{equation}
N_{\rm esc}={\rm erf}\left(\frac{v_{\rm esc}}{\sqrt{2}v_0}\right)-\sqrt{\frac{2}{\pi}}\frac{v_{\rm esc}}{\sigma_v}\exp\left(-\frac{v_{\rm esc}^2}{2\sigma_v^2}\right).
\end{equation}
Applying a Galilean transformation on the velocity distribution in the Galactic frame, we obtain the distribution in the laboratory frame as
\begin{equation}
 f_{\rm lab}(\textbf{v})=f_{\rm gal}(\textbf{v}+\textbf{v}_{\rm lab}).
\end{equation}
 The velocity of the laboratory frame $\textbf{v}_{\rm lab}$ with respect to the galactic rest frame is given by
\begin{equation}
\textbf{v}_{\rm lab}=\textbf{v}_{\rm GalRot}+\textbf{v}_{\rm Solar}+\textbf{v}_{\rm Earth}+\textbf{v}_{\rm EarthRot},
\end{equation}
where $\textbf{v}_{\rm GalRot}$, $\textbf{v}_{\rm Solar}$, $\textbf{v}_{\rm Earth}$, and $\textbf{v}_{\rm EarthRot}$ are the galactic rotation, the peculiar velocity of the Solar system with respect to the galactic halo, Earth’s revolution, and Earth’s rotation velocities respectively. Even though the rotation velocity of the Earth is small compared to the other components, it plays an essential role in the daily modulation of the interaction rate, due to the effect on the direction of the dark matter velocity in the laboratory frame. The relations and algebraic calculations for converting the velocity and coordinates from the galactic frame to the laboratory frame are given in the reference~\cite{bozorgnia:2011}.
Since the direction of the laboratory velocity is towards the constellations of Cygnus~\cite{OHare:2015}, the direction of the velocity of the dark matter can be considered as the inverse of the direction of Cygnus.

Fig.~\ref{cyg&Sun} shows the position of the Sun (yellow line) and the inverse of the Cygnus's position (blue line) in the sky for a laboratory located at latitude 45.2 and longitude 6.7 on June 24, 2020. The red dots indicate the position of the Sun and Cygnus at 6, 12, and 18 o'clock.
Since the direction of the dark matter wind is from the Cygnus and the dominant low energy neutrino background is from the Sun, it is useful to notice that these directions never overlap in the figure.
Fig.~\ref{cygnus} shows the hourly changes in the Cygnus's direction relative to the laboratory's location for latitude 45.2 and longitude 6.7 in one year. Each red dot corresponds to one hour from January 1, 2020, 0:00, until December 31, 2020, 23:00. According to the Earth's rotation around the Sun, the direction of the Cygnus changes during the  year, affecting the daily modulation of the event rate.
\begin{figure}
 \centerline{\includegraphics[width=9cm]{Fig_1.pdf}}
\caption[]{The position of the Sun (yellow) and the direction of the dark matter flux (blue) during a day in latitude 45.2° and longitude 6.7° on June 24, 2020. The red dots indicate the positions at 6, 12, and 18 o'clock. The dashed line in the middle is the horizon. }
 \label{cyg&Sun}
\end{figure}

 \begin{figure}
 \centerline{\includegraphics[width=9cm]{Fig_2.pdf}}
 \caption[]{The hourly changes in the position of the constellation Cygnus during a year for latitude 45.2$^0$ and longitude 6.7$^0$. Each red dot corresponds to one hour.}
 \label{cygnus}
\end{figure}
 The dark matter scattering event rate as a function of recoil energy, direction in the laboratory frame, and time is given by \cite{Kadribasic:2017,Heikinheimo:2019,sassi:2021},
 \begin{equation}
\frac{d^3R}{dE_rd\Omega_r dt}=\frac{\rho \sigma_{m_D-n}}{4\pi m_D \mu_{m_D-n}^2\Delta t}A^2F^2(E_r) \hat{f}_{\rm SHM}(v_{\rm min},\textbf{q},t),
\end{equation}
Where $\rho=0.3$ GeV/cm$^{3}$ is the local dark matter density, $\textbf{q}$ the recoil direction in the detector, $\sigma_{m_D-n}$ is the spin-independent dark matter-nucleon cross-section, $A$ the mass number of the target in the detector, $m_D$ the dark matter particle mass, $\mu_{m_D-n}={m_D m_n}/({m_D+m_n})$ the dark matter-nucleon reduced mass and $F^2(E_r)$ the nuclear Helm form factor~\cite{duda:2007}. Finally, $\hat{f}_{\rm SHM}(v_{\rm min},\textbf{q},t)$ is the radon transform of the dark matter velocity distribution~\cite{Heikinheimo:2019}. The analytical formula of the radon transform given by~\cite{Heikinheimo:2019},
\begin{eqnarray}
\hat{f}_{\rm SHM}(v_{\rm min},\textbf{q},t)=&\frac{1}{N_{\rm esc}(2\pi\sigma_v^2)^\frac{1}{2}} \Bigl[\exp\left(-\frac{|v_{\rm min}+\textbf{q}\cdot \textbf{v}_{\rm lab}|^2}{2\sigma_v^2}\right)\nonumber \\
&-\exp(-\frac{|v_{\rm esc}|^2}{2\sigma_v^2})\Bigr],
\label{f_SHM}
\end{eqnarray}
where $v_{\rm min}=\sqrt{2m_n E_r}/{2\mu_{m_D-N}}$ is the minimum velocity required to produce a nuclear recoil with energy equal to $E_r$, where $\mu_{m_D-N}$ is the dark matter-nucleus reduced mass and $\sigma_v$ is the velocity dispersion.
Equation (6) clearly shows that when $v_{\rm min}<v_{\rm lab}$ , in the case of $\textbf{q}\cdot \textbf{v}_{\rm lab}=-v_{\rm min}$ or $\textbf{q}=- \textbf{v}_{\rm lab}$, the rate has a maximum. In other words, the relationship between the recoil direction and the laboratory motion in the galactic rest frame confirms the existence of a dark matter signal.

To determine the recoil direction in a detector, we consider the axes $\hat{\textbf{x}}$, $\hat{\textbf{y}}$, and $\hat{\textbf{z}}$ as the north, east and vertical directions, respectively,
\begin{equation}
\hat{q}=s\sin(\theta)\cos(\phi)\hat{\textbf{x}}+\sin(\theta)\sin(\phi)\hat{\textbf{y}}+\cos(\theta)\hat{\textbf{z}}.
\end{equation}
In an ionization detector, the total accepted signal rate is obtained by integrating the differential rate over the recoil direction and recoil energy
\begin{equation}
R(t)=\int_{4\pi}\int_{E_{\rm th}(\theta,\phi)}^{E_r^{\rm max}}\frac{d^2R}{dE_rd\Omega_r}{dE_rd\Omega_r},
\label{eq_R}
\end{equation}
where $E_{\rm th}$ is the minimum observable recoil energy, which is here assumed to be correlated with the threshold displacement energy, and therefore depends on the recoil direction with respect to the crystal lattice~\cite{Kadribasic:2017,Heikinheimo:2019}. The threshold energies are obtained from molecular dynamics simulations~\cite{Kadribasic:2017}. The silicon threshold ranges from 17.5 to 77.5 eV.
Fig.~\ref{threshold_energy} represents the angular distribution of the threshold displacement energy. The brighter points correspond to directions of higher threshold energy, indicating that more energy is needed to create a defect, and according to our assumption, to create a detectable ionization signal. The energy threshold surface shown in Fig.~\ref{threshold_energy} therefore determines the lower limit of the energy integral in equation~\ref{eq_R}.
\begin{figure}[ht
\centering
\includegraphics[width=3 in]{Fig_3.pdf}
\qquad
\caption[]{Angular distribution of the defect creation energy threshold of  silicon as a function of the recoil direction.}
\label{threshold_energy}
\end{figure}
\\
\section{Results}\label{sec3}
Here, we study the time evolution of the event rate of dark matter interacting with silicon and germanium nuclei. The dark matter event rate is obtained for each time bin by averaging 24,155 and 84,936 directions for silicon and germanium, respectively. We calculate the rate by integrating over the recoil energy for each sampled direction ($\theta_i$,$\phi_i$) with the corresponding threshold energy value ($E_\text{min}$($\theta_i$,$\phi_i$)). Throughout the analysis, we take the detector mass to be one ton. Fig.~\ref{Ge&Si_daily} shows the time evolution of the event rate of a dark matter particle of mass 300 MeV/c$^2$
interacting with silicon (top) and germanium (bottom)
via spin-independent cross-section 10$^{-39}$ cm$^2$
for January 1, 2020. As expected the
event rates of both detectors oscillate due to the Earth's rotation. As seen in Fig.~\ref{Ge&Si_daily}, for the germanium detector the overall rate is much lower than for silicon. The difference in the total rate is mostly due to the larger nuclear mass in germanium, implying a larger $v_{\min}$ for a given recoil energy, and therefore cutting out a larger fraction of the DM velocity distribution.
\begin{figure}[ht
\centering
\includegraphics[width=6cm]{Fig_4a.pdf}
\includegraphics[width=6cm]{Fig_4b.pdf}
\qquad
\caption[]{Time evolution of the event rate for a 300 MeV/c$^2$ dark matter particle with a 10$^{-39}$ cm$^2$ cross-section in interaction with silicon (top) and germanium (bottom).}
\label{Ge&Si_daily}
\end{figure}

The normalized diurnal modulation of the event rate for the dark matter mass range from 225 to 400 MeV/c$^2$ interacting with silicon nuclei is depicted in Fig.~\ref{Si_mass}. It was already noted in~\cite{Kadribasic:2017,Heikinheimo:2019} that the modulation amplitude increases for decreasing dark matter mass for a germanium target. Here we observe the same effect for silicon. For low-mass dark matter, the diurnal modulation effect can be utilized to discriminate dark matter from the background.

\begin{figure}[ht
\centering
\includegraphics[width=8cm]{Fig_5.pdf}
\qquad
\caption[]{Normalized diurnal modulation of the event rate for different masses range of dark matter from 225 to 400 MeV/c$^2$ interacting with silicon nuclei.}
\label{Si_mass}
\end{figure}
Silicon offers a better sensitivity to low mass dark matter due to its smaller atomic number. This is evident in Fig.~\ref{R_mass}, showing the daily average of the event rate for the dark matter mass range from 200 MeV/c$^{2}$ to 5 GeV/c$^{2}$ for silicon (blue line) and germanium (red line) on January 1, 2020.
The event rate for silicon is about 200 times larger than for germanium at low mass (300 MeV/c$^{2}$). The high event rate for silicon detectors provides more sensitivity for direct dark matter searches for low masses. For high-mass dark matter (1 GeV/c$^{2}$), the event rate for the silicon detector is about 1.5 times greater than the germanium rate. Our numerical studies show that the event rate for silicon detector approaches germanium at dark matter masses larger than 1 GeV/c$^{2}$.
\begin{figure}[ht
\centering
\includegraphics[width=8cm]{Fig_6.pdf}
\qquad
\caption[]{The daily average together with error bars of the event rate for the dark matter mass range from 200 MeV/c${^2}$ to 5 GeV/c$^2$ for silicon (blue line) and germanium (red line) on January 1, 2020. }
\label{R_mass}
\end{figure}

Fig.~\ref{Si_anual} represents the annual modulation of a dark matter particle with  mass of 300 MeV/c$^2$ and a cross-section of 10$^{-39}$ cm$^{2}$ scattering on silicon target over the year 2020.
Due to the rotation of Earth around the Sun, the orientation of the detector with respect to the inverse direction of Cygnus varies throughout the year.
Although the fluctuation reaches its maximum in June, Fig.~\ref{Si_anual} shows that the diurnal modulation remains significant for the low mass of 300 MeV/c$^2$ throughout the year. Since the direction of the neutrinos is from the Sun, the annual interaction rate of dark matter can be considered a helpful tool for distinguishing between dark matter signal and the neutrino background, as discussed in~\cite{sassi:2021}.
\begin{figure}[ht
\centering
\includegraphics[width=8cm]{Fig_7.pdf}
\qquad
\caption{ The annual event rate on silicon for a 300 MeV/c$^2$ dark matter particle with a $10^{-39} cm^2$ cross-section.}
\label{Si_anual}
\end{figure}
To further analyze the periodicity of the event rate, we applied the fast Fourier transformation (FFT) to the event rate time series (Fig.~\ref{Si_anual}). Table~\ref{tab1} represents the periods and amplitudes for silicon event rate. As shown in the table, beyond the leading annual modulation ($\sim$8760 hours), we observed the main periods of about 8 hours (three maxima per day) and 12 hours (two maxima per day).
The patterns with the period of 8 and 12 hours are due to silicon crystal orientation in laboratory coordinates (about 45-degree patterns of the angular distribution of the defect creation energy threshold for the silicon detector, Fig.~\ref{threshold_energy}). As the Earth rotates around its axis, the inverse direction of Cygnus scans various patterns of the minima (dark regions in Fig.~\ref{threshold_energy}) and maxima (bright regions in Fig.~\ref{threshold_energy}) of the threshold energy throughout the day. Depending on the orientation of the crystal frame, the path scanned by the inverse direction of Cygnus may cross two (12-hour period) and three (8-hour period) minima/maxima of the threshold surface during a daily cycle.
\begin{table}
\center
\caption{The period and amplitudes of event rate determined by the fast Fourier transformation (FFT) to the event rate time series.}
\begin{tabular}{ccc}
\\\hline
                     &  Period (hour)            & Amplitudes (Normalized)
 \\\hline
$\mathrm{T_1}$       &   $8756\pm 53.6$          & $(225.4\pm 2.4)\times10^{-3}$      \\
$\mathrm{T_2}$       &   $4170\pm 19.1$          & $(3.8\pm 0.07)\times10^{-3}$ \\
$\mathrm{T_3}$       &   $ 24\pm 0.0005$         & $(2.8\pm 0.04)\times10^{-3}$ \\
$\mathrm{T_4}$       &   $12\pm 0.0001 $         & $(13.1\pm 0.16)\times10^{-3}$ \\
$\mathrm{T_5}$       &   $ 8\pm  0.0008$         & $(2.9\pm 0.05)\times10^{-3}$                   \\\hline
\end{tabular}
\label{tab1}
\end{table}
\section{Conclusions}
We have studied the dark matter interaction rate in silicon assuming that the ionization threshold directional dependence follows similar functional form to the threshold displacement energy simulated using molecular dynamics methods. ~\cite{Holms:2008}.
We obtained a daily and annual modulation of the dark matter event rate for the silicon target. We found that the directional dependence of the threshold energy and the motion of the laboratory with respect to the galactic rest frame results in a modulation of the dark matter event rate in silicon, thus it provides a cosmological signature to identify dark matter from background interactions.
Therefore silicon appears as very appropriate target for identifying dark matter from backgrounds with a high event rate expected for low-mass dark matter. Using  FFT analysis, we obtained the 8 and 12 hours periodicity in simulated time series of event rates for the silicon detector that is due to the 45-degree pattern in the displacement threshold energy surface of silicon. We also observed the periods 24 hour, six months, and one year as the harmonics of the primary modes (8 and 12 hours).

\bibliographystyle{apsrev4-1}