text:
# Bursts in the Chaotic Trajectory Lifetimes Preceding the Controlled Periodic Motion

## Abstract

The average lifetime ($`\tau (H)`$) it takes for a randomly started trajectory to land in a small region ($`H`$) on a chaotic attractor is studied. $`\tau (H)`$ is an important issue for controlling chaos. We point out that if the region $`H`$ is visited by a short periodic orbit, the lifetime $`\tau (H)`$ strongly deviates from the inverse of the naturally invariant measure contained within that region ($`\mu _N(H)^1`$). We introduce the formula that relates $`\tau (H)/\mu _N(H)^1`$ to the expanding eigenvalue of the short periodic orbit visiting $`H`$.

Controlling chaos by stabilizing one of the many unstable periodic orbits embedded within a given chaotic attractor is attainable with small, time-dependent changes in an accessible system parameter. The idea is to observe a typical trajectory of the uncontrolled system for some transient time, until it falls sufficiently close to the desired periodic orbit, and then to activate the control mechanism. An important issue related to the utilization of this method is the average lifetime of chaotic transients that precede the controlled periodic motion .

Suppose that the uncontrolled chaotic attractor $`A`$ describes the asymptotic behavior of the dynamical system $`O:DD,D𝐑^m`$, also referred to as the original system. Let $`\stackrel{}{\xi }=O^k(\stackrel{}{\xi })`$ be a point on a particular unstable periodic orbit that we wish to stabilize. Furthermore, let the vicinity of the orbit $`HH_ϵ(\stackrel{}{\xi })`$ be an $`m`$-dimensional ball of radius $`ϵ<<1`$ centered at $`\stackrel{}{\xi }`$. The probability that a randomly started trajectory does not reach $`H`$ up to time $`t`$ is $`e^{t/\tau (H)}`$. The average lifetime $`\tau \tau (H)`$ is strongly correlated with the visitation frequency of typical trajectories to the region $`H`$, which is described in terms of the naturally invariant measure ($`\mu _N`$) contained within $`H`$ \- $`\mu _N(H)`$. Obviously, if a certain region on a given chaotic attractor is visited more frequently by typical trajectories, the average lifetime it takes for an orbit to land in that region will be smaller. In the present study we address the following question: What is the deviation of $`\tau `$ from $`\mu _N(H)^1`$ as a function of $`\stackrel{}{\xi }`$ and $`ϵ`$?

We will demonstrate the existence of bursts in the lifetimes, i.e. significant deviations of $`\tau `$ from $`\mu _N(H)^1`$, which appear when the $`H`$ region encompasses a point on a short periodic orbit. In contrast to the overall $`\tau \mu _N(H)^1`$ behavior, at these exceptional positions, the lifetime $`\tau `$ is considerably prolonged as compared to $`\mu _N(H)^1`$. As the length of the shortest cycle visiting $`H`$ increases, the parameter of this deviation, $`S(H)\tau /\mu _N(H)^1`$, decreases rapidly towards $`1`$. We will introduce a formula that relates the parameter $`S(H)`$ to the repelling properties (expanding eigenvalue) of the shortest cycle within $`H`$. Furthermore, we will demonstrate that $`S(H)`$ is independent of $`ϵ`$ (for $`ϵ<<1`$). This is consistent with the previously reported scaling $`\tau \mu _N(H)^1`$(see e.g. ).

The present study is motivated by the previous investigation of the logistic map with a hole. In this paper we present a theoretical explanation for the phenomenological result reported in Ref.  and generalize it to 1D noninvertible and 2D invertible chaotic maps. (From our considerations a conjecture follows that similar phenomena occur generally in chaotic systems.)

It will be useful to define an auxiliary modified map

$$M(\stackrel{}{\xi }^{})=\{\begin{array}{cc}O(\stackrel{}{\xi }^{}),\hfill & \stackrel{}{\xi }^{}D\backslash H\hfill \\ \text{outside of the basin of }A,\hfill & \stackrel{}{\xi }^{}H.\hfill \end{array}$$
(1)
A typical trajectory of the map $`O`$ remains on the chaotic attractor forever, while the same trajectory in the map $`M`$ eventually escapes through the region $`H`$, from now on also referred to as the hole. The average lifetime of chaotic transients created by the map $`M`$ is equal to $`\tau (H)`$, which we have defined above. Similar maps with a forbidden gap region arise in the context of communicating with chaos, and in calculation of the topological entropy.

To illustrate the concept of bursts, we consider two chaotic 1D maps: (i) the asymmetric tent map $`O(x)=k_1x,x<k_1;O(x)=k_2(1x),x>k_1`$, $`k_1,k_2>0`$, $`k_1^1+k_2^1=1`$, and (ii) the sinusoidal map $`O(x)=\mathrm{sin}\pi x`$. Fig. 1 displays $`\tau `$ as a function of the position of the hole $`\xi `$ ($`H=(\xi ϵ,\xi +ϵ)`$), for the two paradigmatic maps (see also Ref. ). The width of the hole is kept constant $`(ϵ=0.005)`$. Both graphs exhibit some common features: (i) the overall behavior of lifetimes follows the $`\mu _N(H)^1`$ pattern; (ii) strong local deviations from the $`\mu _N(H)^1`$ behavior - the bursts, observed as leaps in the lifetimes, occur when the hole interval sweeps across a short periodic orbit; (iii) the bursts are more significant if the length of the short periodic orbit is smaller.

The explanation of the burst phenomenon requires the comparison of two concepts: (i) the conditionally invariant measure  (also referred to as the $`c`$-measure) - the concept associated with the modified system, and (ii) the naturally invariant measure  of the original system. In order to define these measures, imagine that we cover the chaotic attractor with cells $`(I)`$ from a very fine grid. Then we randomly distribute a large number $`(N)`$ of points on the grid, and evolve them under the dynamics $`O`$ for a long time $`T`$. Suppose that all initial points are colored blue, and that a point irretrievably changes color from blue to red immediately after its first entrance into the region $`H`$. Thus, the point $`\stackrel{}{x}_T=O^T(\stackrel{}{x}_0)`$ at time $`T`$ is blue if $`O^t(\stackrel{}{x}_0)H`$ for $`t\{0,1,\mathrm{},T1\}`$, and red otherwise. In the limit $`T\mathrm{}`$, the fraction of points found in a given cell $`I`$, is just the natural measure contained within that cell: $`\mu _N(I)=N_T(I)/N;N=_IN_T(I)`$, where $`N_T(I)=b_T(I)+r_T(I)`$ denotes the total number of points in a given cell $`I`$ at time $`T`$ . The number of blue (red) points within a given cell $`I`$ at time $`T`$ is denoted as $`b_T(I)(r_T(I))`$. The points which change color from blue to red under the action of the map $`O`$, are those which would escape the attractor under the action of the map $`M`$. Hence, if we were to evolve exactly the same initial conditions using the map $`M`$ for the same time $`T`$, the number of surviving points (blue points) would be $`B_T=_Ib_T(I)N\mathrm{exp}(T/\tau )`$. In the limit $`N\mathrm{}`$, $`T\mathrm{}`$, the distribution of blue points converges to the $`c`$-measure of the modified system . The fraction of blue points in a given cell $`I`$ is simply the $`c`$-measure contained within that cell: $`\mu _C(I)=b_T(I)/B_T`$.

The blue point at any time $`t>0`$, was certainly not in the hole $`H`$ at time $`t1`$. Therefore, $`M^1(I)O^1(I)\backslash H`$ and $`b_T(I)=b_{T1}(M^1(I))`$, which divided by $`B_T=B_{T1}\mathrm{exp}(1/\tau )`$ yields the well known relation for the $`c`$-measure

$$\mu _C(I)=e^{\frac{1}{\tau }}\mu _C(M^1(I)).$$
(2)
By summing the equation above over all the cells $`I`$ we obtain

$$\frac{1}{\tau }=\mathrm{ln}[1\mu _C(H)]\mu _C(H).$$
(3)
We emphasize the importance of this observation. The average lifetime it takes for a typical trajectory to reach the small region $`H`$ on the attractor is an inverse of the $`c`$-measure contained within that region ($`\mu _C(H)^1`$), which may significantly differ from the inverse of the natural measure ($`\mu _N(H)^1`$).

As an illustration, in Fig. 2 we display the $`c`$-measure for the modified version of the map $`\mathrm{sin}\pi x`$, in comparison to the natural measure of the original map. For the $`c`$-measure in Fig. 2 a), the hole has been positioned at an arbitrary point, but not on the short periodic orbit. In this case, we observe that $`\mu _C(H)\mu _N(H)`$, i.e., $`\mu _N(H)^1`$ is a good approximation for the lifetime. In contrast, in Fig. 2 b) we display $`\mu _C`$ for the modified map $`\mathrm{sin}\pi x`$, with the hole positioned on the fixed point. We notice that $`\mu _C(H)`$ strongly deviates from $`\mu _N(H)`$. This case corresponds to the burst labeled $`1`$ in Fig. 1 b). The overall agreement of the two measures is evident in both Figures 2 a) and b). However, at locations above the first few images of the hole, $`\mu _C`$ takes the shape of a well, with values which are considerably lower than $`\mu _N`$. When the hole lies on the fixed point (Fig. 2 b)), the wells are just above the hole itself. This results in a pronounced deviation of $`\tau =\mu _C(H)^1`$ from $`\mu _N(H)^1`$, which manifests as a burst.

Now we compare the two measures globally. A chaotic repeller is a set of points on the attractor that never visit the hole. A trajectory that starts close to the repeller, does not escape the attractor for a long time. Therefore, the blue points at a large time $`T`$ are located along the unstable manifold of the repeller. Their distribution along this manifold defines the $`c`$-measure. Thus, the natural measure is constructed from all the points (at time $`T`$) on all the unstable manifolds, whereas the $`c`$-measure results only from points on parts of these manifolds. The parts which extend from the repeller up to the hole. As we reduce the size of the hole $`ϵ`$, the repeller and its unstable manifold grow. Consequently, $`\mu _C`$ gradually approaches $`\mu _N`$, and for sufficiently small $`ϵ`$, the two measures are practically identical. (For $`ϵ=0`$, $`\mu _C`$ becomes $`\mu _N`$.)

However, the deviation of the lifetime $`\tau =\mu _C(H)^1`$ from $`\mu _N(H)^1`$ depends only on the values of the two measures within the hole, and therefore is a local quantity. In order to make a more accurate comparison of $`\mu _C`$ and $`\mu _N`$, we introduce the following definitions. Consider a set $`PD`$ such that $`\mu _N(P)>0`$. We define the quantity

$$a(P)=\mu _C(P)/\mu _N(P),$$
(4)
which describes the relation between $`\mu _C`$ and $`\mu _N`$ within $`P`$. We also define the influence $`i(P)`$ of the hole on the set $`P`$ as

$$i(P)=\frac{\mu _N(O^{l(P)}(P)H)}{\mu _N(O^{l(P)}(P))}.$$
(5)
$`l(P)`$ denotes the smallest integer for which the section $`O^{l(P)}(P)H`$ becomes nonempty. The natural invariant measure within $`O^{l(P)}(P)`$ is mapped to $`P`$ in $`l(P)`$ iterates. The influence is just a fraction ($`0i(P)1`$) of $`\mu _N(P)`$ that is mapped from the hole in the last $`l(P)`$ time steps.

Let $`P_ϵP_ϵ(\stackrel{}{x})D`$ be an $`m`$-dimensional ball of radius $`ϵ`$ (the same radius as the hole) centered at $`\stackrel{}{x}`$. We ask the following question: Given a chaotic attractor and choosing the hole region, what is the behavior of $`a(P_ϵ)=\mu _C(P_ϵ)/\mu _N(P_ϵ)`$ as the position of $`P_ϵ`$ on the attractor is changed?

By using Eq. (2) and the identity $`\mu _N(P_ϵ)=\mu _N(O^1(P_ϵ))`$, we can write

$$a(P_ϵ)=e^{\frac{l}{\tau }}(1i(P_ϵ))a(M^l(P_ϵ)),$$
(6)
where $`ll(P_ϵ)`$ (in what follows, $`ll(P_ϵ)`$).

Concerning the first factor in Eq. (6), note that the average lifetime typically scales like $`\tau 1/ϵ^{D_p(\stackrel{}{\xi })}`$ ($`D_p(\stackrel{}{\xi })`$ denotes the pointwise dimension at $`\stackrel{}{\xi }`$) , whereas the minimal number of iterates $`l`$ for which $`O^l(P_ϵ)H\mathrm{}`$ scales like $`l\mathrm{ln}(1/ϵ)`$. Therefore, $`\mathrm{exp}(l/\tau )1+l/\tau 1`$.

If the influence $`i(P_ϵ)`$ is small, the second factor in Eq. (6) is $`1`$. We argue that $`i(P_ϵ)`$, the influence of the hole on the region $`P_ϵ`$ decreases exponentially with $`l`$. For the 2D original map, $`O^l(P_ϵ)`$ is a narrow region which is stretched along the stable direction and squeezed along the unstable one. The intersection of $`O^l(P_ϵ)`$ with the hole $`H`$ is roughly a rectangle of length $`ϵ`$ and width $`ϵ\mathrm{exp}(\lambda _1l)`$. For the 1D map, $`O^l(P_ϵ)H`$ is an interval of width $`ϵ\mathrm{exp}(\lambda _1l)`$. In both cases, $`\lambda _1`$ denotes the positive Lyapunov exponent obtained for typical initial conditions on the attractor. Since the natural measure is concentrated along the unstable manifolds, we can relate $`\mu _N(O^l(P_ϵ)H)\mathrm{exp}(\lambda _1l)`$. Thus, due to the chaoticity of the map $`O`$ we obtain $`i(P_ϵ)\mathrm{exp}(\lambda _1l)`$.

Concerning the third factor in Eq. (6), we consider the set $`M^l(P_ϵ)`$ and the value $`a(M^l(P_ϵ))`$ in dependence of $`l`$. For the 2D maps, the set $`M^l(P_ϵ)O^l(P_ϵ)\backslash H`$ is stretched exponentially fast with increasing $`l`$ along the stable manifolds, and thus crosses many of the unstable manifolds that carry both the natural and the $`c`$-measure. For the 1D maps, the number of disjoint intervals that make the $`lth`$ preimage of $`P_ϵ`$ grows exponentially with $`l`$. Furthermore, they are scattered all over the attractor. Due to the chaoticity of the map $`O`$, in both cases $`M^l(P_ϵ)`$ becomes more democratic with larger $`l`$, in the sense that the value $`a(M^l(P_ϵ))`$ reflects the global agreement between the two measures. Thus, insofar as $`l`$ is not small, $`a(M^l(P_ϵ))1`$.

We conclude that for $`ll(P_ϵ)`$ larger than some critical value (call it $`l_c`$), all of the three factors in Eq. (6) are $`1`$ and therefore $`\mu _C(P_ϵ)\mu _N(P_ϵ)`$. This is consistent with the global agreement between the two measures. The critical value $`l_c`$ depends on the chaoticity of the original map. For example, we may take $`l_c`$ to be the smallest integer for which $`e^{\lambda _1l_c}<0.1`$ (e.g. for the $`\mathrm{sin}(\pi x)`$ map this gives $`l_c34`$). When $`H`$ maps to $`P_ϵ`$ in just a few iterates, so that $`l<l_c`$, a significant difference is observed between $`\mu _C(P_ϵ)`$ and $`\mu _N(P_ϵ)`$ (this explains the wells in Fig. 2).

Coming back to the average lifetimes, if the hole does not map back to itself in just a few iterates, i.e., if the shortest periodic orbit within $`H`$ has a period larger than $`l_c`$, then $`\mu _C(H)\mu _N(H)`$, or simply $`\tau \mu _N(H)^1`$. This explains the overall behavior of lifetimes (see Fig. 1). On the other hand, if $`H`$ encompasses a short periodic orbit (period$`l(H)<l_c`$), the two measures differ within the hole. Quantitatively, we substitute $`P_ϵH`$ in Eq. (6) and approximate $`e^{l(H)/\tau }1`$ and $`\mu _C(M^{l(H)})\mu _N(M^{l(H)})`$. This results in

$$\tau (1i(H))^1\mu _N(H)^1.$$
(7)
Although $`l(H)`$ is small, the approximation $`a(M^{l(H)})1`$ is justified if $`l(M^{l(H)})>l_c`$, or simply, if the period of the second shortest orbit within $`H`$ exceeds $`l_c`$. We have tested relation (7) and consequently the approximations that lead to it in a number of systems. We have compared $`\tau (H)`$ with $`\mu _N(H)^1`$ by changing the position of $`H`$ from ”the most exceptional” point, the fixed point, to longer cycles. In Fig. 3 we display a test of Eq. (7) for the the generalized baker’s map (see Ref. , p. 75, $`\lambda _a=0.35,\lambda _b=0.40,\alpha =0.40,\beta =0.60`$), and for the Hénon map (see Ref.  $`a=1.4,b=0.3`$). Recalling that $`i(H)`$ decreases exponentially with $`l(H)`$, and considering Eq. (7), we see that the parameter $`S(H)=(1i(H))^1`$ decreases rapidly towards $`1`$ with the increase of $`l(H)`$ (see Fig. 3). Eq. (7) is robust and can be applied for holes of different shapes, as long as $`\tau >>1`$. If (for the 2D maps) we tailor the hole as a rectangle with sides of length $`ϵ`$ parallel to the stable and unstable manifold segments, and center it on a short periodic orbit, then

$$\tau (1\mathrm{\Lambda }_u^1)^1\mu _N(H)^1.$$
(8)
$`\mathrm{\Lambda }_u`$ denotes the magnitude of the expanding eigenvalue of that orbit. Eq. (8) also applies to 1D maps. Note that the approximation $`i(H)\mathrm{\Lambda }_u^1`$ assumes that the natural measure is smooth along the unstable direction within $`H`$. We observe that $`\tau /\mu _N(H)^1`$ is independent of $`ϵ`$. This is in accordance with the statement that the lifetime $`\tau `$ scales with $`ϵ`$ just like $`\mu _N(H)^1`$ .

Let us consider an application of Eq. (8). Suppose that we wish to control a chaotic system around an unstable fixed point. In order to obtain the position of the fixed point, its unstable eigenvalue, and other information required for the control, an observation of the free running system is needed. From this observation, we can also evaluate the visitation frequency to the $`ϵ`$-vicinity of the fixed point, i.e. $`\mu _N(HH_ϵ(\stackrel{}{\xi }))`$. $`ϵ`$ is determined by the maximally allowed deviation of the control parameter from its nominal value. The question of interest is how many iterates are needed on the average ($`\tau `$), before a chaotic trajectory enters the region $`H`$, when the control becomes attainable. The prediction given by $`\mu _N(H)^1`$ is an underestimate, since we are on the fixed point. For example, if the underlying dynamics of the system is the asymmetric tent map (with the same parameters as in Fig. 1), and if $`H=(0.62111801\mathrm{}0.002,0.62111801\mathrm{}+0.002)`$, the estimate for the lifetime $`\mu _N(H)^1`$ gives $`250`$ iterates. On the other hand, the numerically calculated lifetime is $`627`$ iterates, which is more than twice as long. The lifetime obtained from formula (8) is $`641`$ iterates, which is very close to the numerically calculated lifetime. Thus, Eq. (8) can be utilized to easily and accurately obtain $`\tau `$ from an observation of the free running system.

In summary, we have studied the average lifetime ($`\tau `$) it takes for a randomly started orbit to land in a small region ($`H`$) on a chaotic attractor. That problem was introduced in Ref.  as an important issue for controlling chaos. Our main result is that if a low-period unstable periodic orbit visits the region $`H`$, then the lifetime $`\tau `$ significantly deviates from the inverse of the natural measure contained within $`H`$ ($`\mu _N(H)^1`$). The parameter of this deviation, $`\tau /\mu _N(H)^1`$, is a function of the expanding eigenvalue of that low-period orbit.


format: text

text:
# 1 Introduction

## 1 Introduction

Search for SUSY particles is one of the primary tasks of the upgraded Fermilab Tevatron and the upcoming CERN Large Hadron Collider (LHC). Because of the unknown masses of the sparticles and other free parameters, various possibilities and strategies have to be considered in the search. Among the plethora of sparticles, superpartners of the top quark, i.e., the stops, especially the lighter of the two mass eigenstates, denoted as $`\stackrel{~}{t}_1`$, is of particular interest. This is because that this stop has color interactions and is likely to be the lightest sfermion. Therefore, it could be produced in the gluon-rich environment of high energy hadron colliders. The lightness of the stop is usually argued for the following reasons. Firstly, the large top quark Yukawa coupling can lead to a large negative one-loop contribution to the stop masses. Hence, the stops could be significantly lighter than other sfermions at the electroweak scale due to the renormalization group evolution, even if all sfermions have an universal mass at the unification scale. Secondly, since the mixing between the sfermions corresponding to the left- and right-handed states of a given fermion is proportional to the mass of the fermion, the large top quark mass can lead to a large mixing of the two stops. This in turn causes a sizable mass splitting between the two mass eigenstates to make the lighter one, i.e., $`\stackrel{~}{t}_1`$, even lighter so as to be accessible to the current and future hadron colliders. Thirdly, the existence of a light stop is preferred by electroweak baryogenesis . Finally, on the theoretical ground, the scenario that the first two-generation sfermions are as heavy as 10 TeV while the third generation sfermions are significantly lighter conforms with the naturalness principle .

In the framework of minimal supersymmetric model (MSSM) with $`R`$-parity conservation, several possibilities of light stop searches from top quark decay have been considered in specific scenarios in the literature. We recapitulate them briefly below.

If the stop is the next-to-the-lightest super particle (NLSP), its only two-body decay mode is $`\stackrel{~}{t}_1c\stackrel{~}{\chi }_1^0`$ via loops , where the lightest neutralino, $`\stackrel{~}{\chi }_1^0`$, is assumed to be the lightest sparticle (LSP). In the case that the stop is sufficiently light, one can consider the exotic top decay $`t\stackrel{~}{t}_1\stackrel{~}{\chi }_1^0`$ followed by $`\stackrel{~}{t}_1c\stackrel{~}{\chi }_1^0`$. Studies showed that this decay chain in the $`t\overline{t}`$ pair events, if realized, can be observable in a large part of the SUSY parameter space at the future runs of the Tevatron collider .

Another possible decay mode of the stop in the case of $`R`$-parity conservation, if the stop is light but heavier than the NLSP, is $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+`$ through the tree-level coupling, where $`\stackrel{~}{\chi }_1^+`$ denotes the lightest chargino. This decay mode will be the dominant decay channel of $`\stackrel{~}{t}_1`$ whenever it is allowed kinematically. The phenomenology of the $`t\overline{t}`$ production followed by the decay chain $`t\stackrel{~}{t}_1\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+`$ has been studied soon after the observation of the top quark . But now the significant higher values of the lower bounds of the masses of $`\stackrel{~}{t}_1`$ and $`\stackrel{~}{\chi }_1^0`$, which are given by about 122 GeV  and 45 GeV  respectively, albeit under certain assumptions, make these top decay chains discussed above less likely.

From the above discussion we see that the discovery of the light stop through the top quark decay is not a promising possibility. However, the light stop offers a more direct route for its discovery at the Tevatron and LHC, i.e., the direct production of the stop pair, assuming the production cross section is sufficiently large and there exits a suitable decay channel for its identification. We note that since the stop pair production is a QCD process, the only uncertainty in the production cross section is the mass of the stop and how it decays.

If the stop is the NLSP and thus its only two-body decay mode is $`\stackrel{~}{t}_1c\stackrel{~}{\chi }_1^0`$, then the signal of the stop pair production at hadron colliders can only be two jets plus missing energy . The large QCD background renders the signal impossible to uncover. In this article, we examine, in the MSSM with $`R`$-parity conservation, the case of a stop with a mass close to that of the top quark and is heavier than the lightest chargino so that it decays dominantly through $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+`$, followed by $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^0f\overline{f}^{}`$ via a real or virtual $`W`$-boson intermediate state, where $`\stackrel{~}{\chi }_1^+`$ is the lightest chargino and $`\stackrel{~}{\chi }_1^0`$ the LSP. Then a stop pair event will look like a top quark pair event and can be easily masked by the latter . Through a detailed Monte Carlo simulation, the possibility of uncovering the possible stop pair events from the top quark counting experiment at the Tevatron and LHC colliders is investigated in detail. As also being demonstrated below, for a stop as light as the top quark, the stop sample will generally be buried in the top sample at Run 1 because of the small statistics. But at the future runs of the upgraded Tevatron and at the LHC, such a stop sample can be revealed through a series of suitable selection cuts.

## 2 Stop pair production and signatures

Similar to the top quark pair production, in hadron collisions the stop pair can be produced in the $`q\overline{q}`$ annihilation and gluon-gluon fusion due to the $`g\stackrel{~}{t}_1\overline{\stackrel{~}{t}}_1`$ coupling. The lowest-order matrix elements will be used in our Monte Carlo simulation. The absolute values squared of the two processes are given by

$`||^2(q\overline{q}\stackrel{~}{t}_1\overline{\stackrel{~}{t}}_1)`$ $`=`$ $`16g_s^4{\displaystyle \frac{\widehat{t}_1\widehat{u}_1m_{\stackrel{~}{t}_1}^2\widehat{s}}{\widehat{s}^2}},`$ (1)
$`||^2(gg\stackrel{~}{t}_1\overline{\stackrel{~}{t}}_1)`$ $`=`$ $`2g_s^4\left[24\left(12{\displaystyle \frac{\widehat{t}_1\widehat{u}_1}{\widehat{s}^2}}\right){\displaystyle \frac{8}{3}}\right]\left[12{\displaystyle \frac{m_{\stackrel{~}{t}_1}^2\widehat{s}}{\widehat{t}_1\widehat{u}_1}}\left(1{\displaystyle \frac{m_{\stackrel{~}{t}_1}^2\widehat{s}}{\widehat{t}_1\widehat{u}_1}}\right)\right],`$ (2)
where $`\widehat{s}`$ is the center-of-mass energy squared of the parton process, $`\widehat{t}_1=\widehat{t}m_{\stackrel{~}{t}_1}^2`$ and $`\widehat{u}_1=\widehat{u}m_{\stackrel{~}{t}_1}^2`$ with $`\widehat{t}`$ and $`\widehat{u}`$ being Mandelstam variables. For parton distribution functions, we use CTEQ5L with $`\mu =\sqrt{\widehat{s}}`$ . For a stop with a mass close to the top quark, the QCD corrections enhance the total cross section of stop pair by a factor of $`1.2`$ at the Tevatron and $`1.4`$ at the LHC energies . These enhancements are taken into account in our calculation.

Although the QCD coupling in stop production processes is as strong as the top quark production, the production rate of stop pair at a given energy is much smaller than the top pair for similar masses. The suppression of the stop pair production is caused largely by the fact that they are spin-0 particles: (1) There is no sum over the spin projects of the final states that can enhance the production rate by several fold. (2) The $`P`$-wave coupling in the $`q\overline{q}`$ annihilation process give rise to a $`\beta ^3`$-dependence  that caused the cross section to be suppressed strongly near the threshold. While the suppression factors work at all collider energies, the stop production at the Tevatron is suppressed more severely because the dominant production of stop at the Tevatron is through the $`q\overline{q}`$ annihilation.

As stated in the Introduction, we focus on the possibility of the $`\stackrel{~}{t}_1`$ decay chain $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+`$ and $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^0f\overline{f}^{}`$. We assume a SUSY spectrum in which all the sparticles involved in the decay chain are on-shell. We will specify the relevant mass values of the sparticles below. As pointed out early, the two-body decay channel $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+`$ will be dominant if the final state particles are on shell. Thus in our calculation we approximate the branching ratio of this mode as 100%. For the subsequent 3-body decay of the chargino, we use its full matrix element in our Monte Carlo simulation and take the total width of chargino to be the sum of these 3-body decay channels of all allowed $`f\overline{f}^{}`$. We also assume that the charged Higgs boson, sleptons and squarks are much heavier than the $`W`$-boson so that these three-body decays proceed dominantly through the $`W`$-boson intermediate state .

So the stop pair $`\stackrel{~}{t}_1\overline{\stackrel{~}{t}_1}`$ production followed by the decay chain $`\stackrel{~}{t}_1b\stackrel{~}{\chi }_1^+b\stackrel{~}{\chi }^0f\overline{f}^{}`$ gives rise to top-like signatures except for two extra neutralinos which escape detection. There are three possible observing channels for the stop-pair event: dilepton+2-jet, single lepton+4-jet and all(six)-jet. All three channels are associated with a significant amount of missing energies. The all-jet channel has the largest rate but will subject to a very large QCD background, and thus is not suitable for isolating the stop signal. The dilepton channel has the lowest rate and, furthermore, it is difficult to find a mechanism to enhance the stop/top rate to find the ”smoking gun” of the stop pair production. So we will not use it, either. In the remaining single lepton+4-jet channel, the best signal is $`\mathrm{}+4j/b+\overline{)E_T}`$ for the purpose of distinguishing the stop event from a top pair event. Here $`4j/b`$ represents a 4-jet event with at least one of the jets passing the $`b`$-tagging criterion. As is shown below, we can find very effective selection cuts to enhance the stop/top ratio for this signal.

## 3 Relevant SUSY parameters

There are several SUSY parameters involved in our calculation. First of all, the stop mass is the most important parameter. We will fix it in the range of 170 GeV in most of our numerical examples. But we will vary it to find out how heavy it can be for the stop signal to be observable. For the neutralino and chargino masses, as well as their couplings, there are four independent parameters: $`M`$, $`M^{}`$, $`\mu `$ and $`\mathrm{tan}\beta `$. $`M`$ is the $`SU(2)`$ gaugino mass and $`M^{}`$ the hypercharge $`U(1)`$ gaugino mass, $`\mu `$ the coefficient of the Higgs mixing term in the superpotential, $`\mu H_1H_2`$, and $`\mathrm{tan}\beta =v_2/v_1`$ the ratio of the vacuum expectation values of the two Higgs doublets. We work in the framework of the MSSM and assume the grand unification of the gaugino masses, which gives the relation $`M^{}=\frac{5}{3}M\mathrm{tan}^2\theta _W0.5M`$. This reduces the independent parameters needed to three: $`M`$, $`\mu `$ and $`\mathrm{tan}\beta `$. For the three independent parameters, the chargino-neutralino sector can be divided into two regions, the gaugino-like region ($`M<|\mu |`$) and the higgsino-like region ($`M>|\mu |`$).

The gaugino-like region is favorable for the discovery of the stop signal. In this region the lightest neutralino $`\stackrel{~}{\chi }_1^0`$ is mainly composed of the hypercharge $`U(1)`$ gaugino (bino), and the lightest chargino $`\stackrel{~}{\chi }_1^+`$ is mainly composed of the charged $`SU(2)`$ gaugino (wino). So the $`\stackrel{~}{\chi }_1^0`$ mass is about half of that of the $`\stackrel{~}{\chi }_1^+`$. The large mass splitting between $`\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\chi }_1^+`$ is needed to produce the required energetic jets or lepton in the decay $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^0f\overline{f}^{}`$, so that they can pass the necessary kinematic cuts.

On the contrary, the higgsino-like region is unfavorable for the stop signal. In this case, both the lightest neutralino $`\stackrel{~}{\chi }_1^0`$ and the lightest chargino $`\stackrel{~}{\chi }_1^+`$ are mainly composed of the higgsino fields. As a result, $`\stackrel{~}{\chi }_1^0`$ is almost degenerate with (but lighter than) $`\stackrel{~}{\chi }_1^+`$. Then the lepton or jets produced in the decay $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^0f\overline{f}^{}`$ will be too soft to pass our selection cuts. So in this case the stop signal will be significantly reduced and likely hidden under the top events even stop pairs are produced.

In our calculation we choose the following representative set of values for the parameters in the gaugino-like region

$`M=100\mathrm{GeV},\mu =200\mathrm{GeV},\mathrm{tan}\beta =1.`$ (3)
The chargino and neutralino masses in units of GeV are then given by

$`m_{\stackrel{~}{\chi }_1^+}=120,m_{\stackrel{~}{\chi }_2^+}=220,`$
$`m_{\stackrel{~}{\chi }_1^0}=55,m_{\stackrel{~}{\chi }_2^0}=122,m_{\stackrel{~}{\chi }_3^0}=200,m_{\stackrel{~}{\chi }_4^0}=227.`$ (4)
As expected, $`m_{\stackrel{~}{\chi }_1^0}`$ is about half of $`m_{\stackrel{~}{\chi }_1^+}`$.

It should be remarked that SUSY parameters are generally not well-constrained experimentally at the present time. The only robust constraints are the LEP and Tevatron lower bounds on some of the sparticle masses . In addition, the intermediate value of $`\mathrm{tan}\beta `$ is favored by low energy experiments . Therefore, the above SUSY parameter values used in our calculation are not the only choice. They are a set of representative values which are allowed by the current experimental bounds and often applied for simulation.

## 4 Selection of $`\mathrm{}+4j/b+\overline{)E_T}`$ events

In making the analyses, we simulate the energy resolution of the detector by assuming a Gaussian smearing of the energy of the final state particles,

$`\mathrm{\Delta }E/E`$ $`=`$ $`30\%/\sqrt{E}1\%,\mathrm{for}\mathrm{leptons},`$ (5)
$`=`$ $`80\%/\sqrt{E}5\%,\mathrm{for}\mathrm{hadrons},`$ (6)
where $`E`$ is in GeV, and $``$ indicates that the energy-dependent and energy-independent terms are added in quadrature.

The basic selection cuts are chosen as follows. For the Tevatron, the cuts are

$`p_T^{\mathrm{}}`$ $``$ $`20\mathrm{GeV},`$
$`p_T^{\mathrm{miss}}`$ $``$ $`20\mathrm{GeV},`$
$`p_T^{jet}`$ $``$ $`15\mathrm{GeV},`$
$`\eta _{jet},\eta _{\mathrm{}}`$ $``$ $`2.0,`$
$`\mathrm{\Delta }R_{jj},\mathrm{\Delta }R_j\mathrm{}`$ $``$ $`0.5.`$ (7)
For the LHC, the cuts are chosen to be

$`p_T^{\mathrm{}}`$ $``$ $`20\mathrm{GeV},`$
$`p_T^{\mathrm{miss}}`$ $``$ $`30\mathrm{GeV},`$
$`p_T^{jet}`$ $``$ $`20\mathrm{GeV},`$
$`\eta _{jet},\eta _{\mathrm{}}`$ $``$ $`3.0,`$
$`\mathrm{\Delta }R_{jj},\mathrm{\Delta }R_j\mathrm{}`$ $``$ $`0.4.`$ (8)
Here $`p_T`$ denotes the transverse momentum, $`\eta `$ is the pseudo-rapidity, and $`\mathrm{\Delta }R`$ is the separation in the azimuthal angle-pseudo rapidity plane $`(\mathrm{\Delta }R=\sqrt{(\mathrm{\Delta }\varphi )^2+(\mathrm{\Delta }\eta )^2})`$ between a jet and a lepton or between two jets.

For the signal, we require to tag at least one $`b`$ in $`\mathrm{}+4j/b+\overline{)E_T}`$. The tagging efficiency is 53% at Run 1 and expected to reach 85% at Run 2 and Run 3 . For the LHC we assume the tagging efficiency to be the same as the Tevatron Run 2. Under the above basic selection cuts and $`b`$-tagging, the ratio of the top events $`\mathrm{}+4j/b+\overline{)E_T}`$ to the QCD backgrounds is about 12:1 , which we will use to evaluate the QCD backgrounds.

We noticed that for the top events and $`W`$+jets background events the missing energy comes only from the neutrino of the W decay, while for the stop events the missing energy contains two extra neutralinos. From the transverse momentum of the lepton, $`\stackrel{}{P}_T^{\mathrm{}}`$, and the missing transverse momentum, $`\stackrel{}{P}_T^{\mathrm{miss}}`$, we construct the transverse mass

$$m_T(\mathrm{},p_T^{\mathrm{miss}})=\sqrt{(|\stackrel{}{P}_T^{\mathrm{}}|+|\stackrel{}{P}_T^{\mathrm{miss}}|)^2(\stackrel{}{P}_T^{\mathrm{}}+\stackrel{}{P}_T^{\mathrm{miss}})^2}.$$
(9)
As is well-known, if $`P_T^{\mathrm{}}`$ and $`p_T^{\mathrm{miss}}`$ are from the decay products of a parent particle, the transverse mass is bound by the mass of the parent particle. For the top quark and $`W`$+jets background events, where the only missing energy is from the neutrino of the W decay, $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ is always less than $`M_W`$ and peaks just below $`M_W`$, although kinematic smearings can push the bound and the peak above $`M_W`$. For the stop events, there is no such peak due to the extra missing energies of the neutralinos. The transverse mass distributions of the stop and top quark events are shown in Fig.1. The transverse mass distribution of the top quark events indeed conforms with the expectation, i.e., it peaks just below 80 GeV and significant distribution appears above 80 GeV due to the smearing. In order to substantially enhance the ratio stop/top, Fig. 1 suggests that we apply the following cut,

$$m_T(\mathrm{},p_T^{\mathrm{miss}})50100\mathrm{G}\mathrm{e}\mathrm{V}.$$
(10)
To further enhance the stop/top ratio, we construct four different invariant masses, denoted as $`M(3j)`$ by using three jets out of the four jets in the event. We define the one which is closest to 175 GeV as the reconstructed top quark and denote its value as $`M_{\mathrm{top}}(3j)`$. The $`M_{\mathrm{top}}(3j)`$ distribution at the upgraded Tevatron energy is shown in Fig.2. As expected, for the top quark events, there is a peak at the top quark mass, $`M_t=175`$ GeV. To enhance the stop/top ratio further, we suppress the top quark events by selecting $`M_{\mathrm{top}}(3j)`$ to be 20 GeV away from the top quark mass (175 GeV), i.e.,

$`|M_{\mathrm{top}}(3j)M_t|`$ $``$ $`20\mathrm{GeV}.`$ (11)
## 5 Numerical results

For the parameter values specified in Sec. 3, the $`\mathrm{}+4j/b+\overline{)E_T}`$ cross section from the stop and top quark events at the Tevatron and LHC, under various cuts, are presented in Table 1. The basic selection cuts Eqs. (7) and (8), which are necessary to reduce the QCD background, affect the stop cross section more than that of the top quark and therefore lower the stop/top ratios. This is because the lepton and jets from the stop events are relatively softer and thus harder to pass the selection cuts, resulting in the stop/top ratios below 5%. This makes the discovery of a stop impossible if the systematic uncertainty is considered. However, the cut of the reconstructed top quark mass suppresses the top quark and stop cross sections respectively by factors of about 12 and 1.6, drastically enhancing the stop/top ratio. The transverse mass cut can further increase the stop/top ratio by suppressing the top cross section by a factor of about 3.3 and that of the stop by 1.9.

In extracting the new physics signal from the top quark events, various uncertainties have to be taken into account besides the experimental statistical and systematic errors. The present uncertainty of the standard model $`t\overline{t}`$ cross section is at the 5% level . Additional uncertainties that come from the error in $`m_t`$ will be much reduced with the expected more precise determination of $`m_t`$ (within 2.8 and 0.8 GeV as quoted for Run 2 and 3 ). For illustration we use the total systematic error of 5%. Combined with the statistical error for each run, we obtain the total errors for the top quark events as listed in Table 2. In estimating the statistical errors, the QCD backgrounds which are mainly W+jets have also been taken into account. Under the basic selection cuts they are reduced to about 1/12 of the top events in the channel of $`\mathrm{}+4j/b+\overline{)E_T}`$ . The cut on the reconstructed top mass $`M_{\mathrm{top}}(3j)`$ cannot suppress the QCD backgrounds significantly and, therefore, to be conservative, we neglect such suppressions. However, the transverse mass cut is expected to suppress the QCD W+jets backgrounds significantly because the missing energy is only from the neutrino of the $`W`$ decay just like in a top event. So we assumed the transverse mass cut suppress the QCD backgrounds by a factor of about 3 as in the top quark case.

As indicated in Tables 1 and 2, although the stop/top ratio is enhanced significantly by the suitable cuts, Run 1 of the Tevatron is unable to observe the stop events because of the small statistics. So even for the favorable case (gaugino-like region) under consideration, the stop pair events will still be hidden in the top pair samples.

For Run 2A with a luminosity of 2 fb<sup>-1</sup>, the statistical error is still large. As showed in Table 2, after combined with the 5% systematic error, the total error is 6%, 14% and 24% under three different selection cuts. Comparing with the stop/top ratio in Table 1, which is 4%, 30% and 49% under the corresponding cuts, we see that the $`M_{\mathrm{top}}(3j)`$ and $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cuts drive the sensitivity to the $`2\sigma `$ level. This is still below the discovery limit which is usually required to be a $`5\sigma `$ deviation or more.

For Run 2B (15 fb<sup>-1</sup>) and Run 3 (30 fb<sup>-1</sup>), the statistical errors are significantly reduced, as shown in Table 2. For example, comparing the total errors (5%, 6%, 8% under the three selection cuts) at Run 3 with the stop contributions (4%, 30%, 49% under the three corresponding cuts), we conclude that the stop event under the $`M_{\mathrm{top}}(3j)`$ and $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cuts is observable ($`5\sigma `$).

For LHC, because of the large production rates, even for the low luminosity run (say $`10`$ fb<sup>-1</sup>), the statistical error is reduced to be negligible. So the total error under each selection cut is dominated by the systematical error which is assumed to be 5%. Comparing with the stop contributions (5%, 40%, 62% under the three selection cuts ), one sees that the stop sample after the $`M_{\mathrm{top}}(3j)`$ and $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cuts will undoubtedly be observable.

In the above results of stop events we fixed stop mass to be 170 GeV. In Figs.3 and 4 we present the stop/top ratio versus stop mass under the basic plus $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ plus $`M_{\mathrm{top}}(3j)`$ cuts. The horizontal dotted lines are the limits required by the discovery ($`5\sigma `$), evidence ($`3\sigma `$) and (if not observed) exclusion ($`2\sigma `$) of the production of stop pairs. We see that the LHC (10 fb<sup>-1</sup>) is able to discover a $`135215`$ GeV stop, while Run 2B (15 fb<sup>-1</sup>) of the Tevatron can discover a $`135175`$ GeV stop. If not discovered, a stop lighter than $`245`$ ($`200`$) GeV will be excluded by LHC (Run 2B of the Tevatron) at $`95\%`$ C.L.. Of course, such results are valid only for the gaugino-like scenario with the specific parameter values we considered as given in Sec. 3

The peaks in Figs. 3 and 4 are the artifact of the cuts applied and can be understood as follows. As the stop mass decreases, the stop pair production rate increases. However, the $`b`$-jet from $`\stackrel{~}{t}_1\stackrel{~}{\chi }_1^+b`$ becomes softer and thus harder to pass the selection cuts so as to decrease the ratio stop/top. At low values of stop mass, this latter effect is stronger and thus the net effect is to decrease the ratio stop/top for decreasing stop mass. At the high end of the stop mass, the phase space suppression of the production rate leads to decreasing stop/top for increasing stop mass. The balance of these opposite effects give rise to the peaks. The peak will shift to higher value for higher parton center of mass energy as shown in Figs. 3 and 4.

## 6 Summary and discussion

We have investigated the potential of the detection of the top-like events of light stop pair at hadron colliders for the case that the mass of the lighter stop is close to that of the top quark, and the SUSY parameters lie in a range which allows the extraction of the stop signal. Because of the much lower production rate relative to the top quark pair production, the extraction of the stop event from the top sample requires special consideration to enhance the ratio stop/top. Through Monte Carlo simulation with suitable cuts, we found that a stop signal in the channel $`\mathrm{}+4j/b+\overline{)E_T}`$ may be detectable in the top counting experiment in the upgraded Tevatron and LHC. However, because of the small statistics, Run 1 of the Tevatron is unable to detect such a stop presence, leaving the stop event hidden in the top pair sample even if it is produced.

We note that our calculation represents the results of a limited set of numerical examples rather than the scanning of the whole SUSY parameter space allowed. Our results are dependent on the mass values of the sparticles involved: the stop $`\stackrel{~}{t}_1`$, chargino $`\stackrel{~}{\chi }_1^+`$ and the neutralino $`\stackrel{~}{\chi }_1^0`$. To validate our analysis, the stop mass must be larger than those of the lightest chargino and neutralino, and their mass spectrum has to be gaugino like. Finally, as already stated in Sec. 3, we also note that if the lightest neutralino and chargino are higgsino-like and thus their masses are close, the leptons and jets in the final states would be too soft to pass our proposed isolation cuts. Then, the stop signal would not be observable. They will remain hidden in the top pair sample even if the stops are produced at the upgraded Tevatron and LHC.

It should be pointed out that throughout our analysis we worked in the MSSM with $`R`$-parity conservation. If R-parity is violated, there are also some interesting phenomenologies in the top-stop sector at the Tevatron and LHC energies, some of which have been explored elsewhere .

## Acknowledgment

JMY thanks Y. Sumino and C.-H. Chang for helpful discussions. BLY acknowledges the hospitality extended to him by Professor Zhongyuan Zhu and colleagues at the Institute of Theoretical Physics, Academia Sinica, where part of the work was performed. This work is supported in part by a grant of Chinese Academy of Science for Outstanding Young Scholars.

Table 1: Cross sections for $`\mathrm{}+4j/b+\overline{)E_T}`$, from stop and top quark pairs at the Tevatron and LHC. The basic cuts are given in Eqs.(7) and (8). The $`M_{\mathrm{top}}(3j)`$ and $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cuts are given by $`|M_{\mathrm{top}}(3j)M_t|>20`$ GeV and $`m_T(\mathrm{},p_T^{\mathrm{miss}})50100`$ GeV. The stop events were calculated by assuming $`M_{\stackrel{~}{t}_1}=170`$ GeV, $`M=100`$ GeV, $`\mu =200`$ GeV and $`\mathrm{tan}\beta =1`$. Tagging at least one $`b`$-jet is assumed for 53% efficiency for the Tevatron (1.8 TeV), $`85\%`$ efficiency for the upgraded Tevatron (2 TeV) and LHC. The ’No cut’ column gives the result under the condition of b-tagging in the absence of any kinematic cuts. The charge conjugate channels are included.

|  |  | no cut | basic cuts | basic cuts | basic cuts |
| --- | --- | --- | --- | --- | --- |
|  |  |  |  | $`+`$ | $`+`$ |
|  |  |  |  | $`M_{\mathrm{top}}(3j)`$ cut | $`M_{\mathrm{top}}(3j)`$ cut |
|  |  |  |  |  | $`+`$ |
|  |  |  |  |  | $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cut |
|  | stop (fb) | 59 | $`12`$ | $`7.4`$ | $`3.9`$ |
| Tevatron | top (fb) | 750 | $`317`$ | $`26.0`$ | $`8.0`$ |
| (1.8 TeV) | stop/top (%) | 7.9 | $`3.8`$ | $`28.5`$ | $`48.8`$ |
|  | stop (fb) | 136 | $`27`$ | $`17`$ | $`8.8`$ |
| Tevatron | top (fb) | 1652 | $`690`$ | $`57`$ | $`18`$ |
| (2 TeV) | stop/top (%) | 8.2 | $`3.9`$ | $`30.0`$ | $`48.9`$ |
|  | stop (pb) | 26 | $`2.92`$ | $`1.94`$ | $`0.90`$ |
| LHC | top (pb) | 170 | $`60.8`$ | $`4.86`$ | $`1.46`$ |
| (14 TeV) | stop/top (%) | 15.3 | $`4.8`$ | $`40.0`$ | $`62.0`$ |

Table 2: Numbers of expected top quark events in the channel $`\mathrm{}+4j/b+\overline{)E_T}`$, together with the associated errors. The estimated total error is obtained by combing the statistical errors and a 5% systematic uncertainty.

|  |  | basic cuts | basic cuts | basic cuts |
| --- | --- | --- | --- | --- |
|  |  |  | $`+`$ | $`+`$ |
|  |  |  | $`M_{\mathrm{top}}(3j)`$ cut | $`M_{\mathrm{top}}(3j)`$ cut |
|  |  |  |  | $`+`$ |
|  |  |  |  | $`m_T(\mathrm{},p_T^{\mathrm{miss}})`$ cut |
|  | top events | $`32`$ | $`3`$ | $`1`$ |
| Run 1 ($`0.1`$ fb<sup>-1</sup>) | Stat. error (%) | $`18.5`$ | $`81.6`$ | $`141.4`$ |
|  | Total error (%) | $`19.2`$ | $`81.8`$ | $`141.5`$ |
|  | top events | $`1380`$ | $`114`$ | $`36`$ |
| Run 2A ($`2`$ fb<sup>-1</sup>) | Stat. error (%) | $`2.8`$ | $`13.3`$ | $`23.6`$ |
|  | Total error (%) | $`5.7`$ | $`14.2`$ | $`24.1`$ |
|  | top events | $`10350`$ | $`855`$ | $`270`$ |
| Run 2B ($`15`$ fb<sup>-1</sup>) | Stat. error (%) | $`1.0`$ | $`4.8`$ | $`8.6`$ |
|  | Total error (%) | $`5.1`$ | $`7.0`$ | $`10.0`$ |
|  | top events | $`20700`$ | $`1710`$ | $`540`$ |
| Run 3 ($`30`$ fb<sup>-1</sup>) | Stat. error (%) | $`0.7`$ | $`3.4`$ | $`6.1`$ |
|  | Total error (%) | $`5.1`$ | $`6.1`$ | $`7.9`$ |
|  | top events | $`6.1\times 10^5`$ | $`4.9\times 10^4`$ | $`1.5\times 10^4`$ |
| LHC ($`10`$ fb<sup>-1</sup>) | Stat. error (%) | $`0.1`$ | $`0.6`$ | $`1.2`$ |
|  | Total error (%) | $`5.0`$ | $`5.0`$ | $`5.1`$ |


text:
# Solitons in quadratic nonlinear photonic crystals

## Abstract

We study solitons in one-dimensional quadratic nonlinear photonic crystals with modulation of both the linear and nonlinear susceptibilities. We derive averaged equations that include induced cubic nonlinearities and numerically find previously unknown soliton families. The inclusion of the induced cubic terms enables us to show that solitons still exist even when the effective quadratic nonlinearity vanishes and conventional theory predicts that there can be no soliton. We demonstrate that both bright and dark forms of these solitons are stable under propagation.

The physics and applications of photonic band-gap (PBG) materials, or photonic crystals, have been active topics of research for more than a decade. The theory of linear photonic crystals is now well understood, and many of their fundamental properties and technical applications have been characterized . The next important step in the application of photonic crystals is to create tunable PBGs. Tunability is possible in linear photonic crystals through, e.g., the temperature dependence of the refractive index  or the electro-optic effect . Ultrafast dynamical tunability of the PBG can be accomplished using nonlinearity, as was first demonstrated with a constant Kerr nonlinearity .

Here we consider quadratic nonlinear photonic crystals (QNPCs) that have a linear grating (periodic dielectric constant) and/or a nonlinear grating (periodic second-order or $`\chi ^{(2)}`$ susceptibility). QNPCs are of interest for all-optical components due to the fast and strong nonlinearity they can provide through the parametric cascading effect . The efficiency of the cascading process depends critically on the phase mismatch between the fundamental and second-harmonic (SH) waves, but two powerful methods exist that use exactly a periodic photonic crystal structure to control the mismatch . In one method, a QNPC with a linear Bragg grating is used to bend the dispersion curve near the PBG . However, the short period, which is of the order of the optical wavelength, can be inconvenient. The second scheme, quasi-phase-matching (QPM), controls the phase mismatch  using a nonlinear grating with a period equal to the comparatively long beat length (typically of the order of microns). QPM is also possible with linear gratings, but this is much less effective .

One of the spectacular manifestations of nonlinearity is the soliton, a self-localized entity that can propagate unchanged over long distances. Homogeneous $`\chi ^{(2)}`$ materials support solitons in all dimensions , and gap solitons exist in QNPCs with a linear Bragg grating . In this Letter we focus on the open fundamental problem of whether solitons exist in 1D QNPCs with both a linear and a nonlinear QPM grating. Such a simultaneous linear grating is difficult to avoid when, for example, creating nonlinear QPM gratings in GaAs/AlAs semiconductors through quantum-well disordering .

Solitons exist in 1D QNPCs with a nonlinear QPM grating , but a simultaneous linear grating can reduce the effective $`\chi ^{(2)}`$ nonlinearity . Thus the global existence of solitons in such QNPCs is nontrivial. Building on previous findings that nonlinear QPM gratings induce cubic nonlinearities , we find soliton solutions that are stable under propagation. The induced cubic terms provide an elegant means of showing that, against intuition, the QNPC supports stable bright and dark solitons even when there is no effective $`\chi ^{(2)}`$ nonlinearity. This is analogous to the existence of solitons in dispersion-managed fibres with no average dispersion .

We consider the interaction of a cw beam (carrier frequency $`\omega `$) with its SH, propagating in a lossless 1D QNPC under conditions for type I second-harmonic generation (SHG), as sketched in Fig. 1. We assume that the modulation of the refractive index is weak ($`\mathrm{\Delta }n_j(z)/\overline{n}_j`$$``$1, where $`n_j(z)`$=$`\overline{n}_j+\mathrm{\Delta }n_j(z)`$ and $`j`$ refers to the frequency $`j\omega `$), and we consider only gratings for forward QPM. The grating period is then much longer than the optical period, in which case Bragg reflections can be neglected. The evolution of the slowly varying beam envelopes is then described by

$`i{\displaystyle \frac{w}{z}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2w}{x^2}}+\alpha _1(z)w+\chi (z)w^{}v\mathrm{e}^{i\beta z}`$ $`=`$ $`0,`$ (1)
$`i{\displaystyle \frac{v}{z}}+{\displaystyle \frac{1}{4}}{\displaystyle \frac{^2v}{x^2}}+2\alpha _2(z)v+\chi (z)w^2\mathrm{e}^{i\beta z}`$ $`=`$ $`0,`$ (2)
where $`w`$=$`w(x,z)`$ and $`v`$=$`v(x,z)`$ are the envelope functions of the fundamental and SH, respectively. The transverse and propagation coordinates $`x`$ and $`z`$ are in units of the input beam width $`x_0`$ and the diffraction length $`L_d`$=$`k_1x_0^2`$, respectively. The parameter $`\beta `$=$`\mathrm{\Delta }kL_d`$ is proportional to the mismatch $`\mathrm{\Delta }k`$=$`k_2`$-$`2k_1`$, $`k_j`$=$`j\omega \overline{n}_j/c`$ being the average wavenumber. Thus $`\beta `$ is positive for normal dispersion and negative for anomalous dispersion. The normalized refractive index grating is given by $`\alpha _j(z)`$=$`L_d\omega \mathrm{\Delta }n_j(z)/c`$ and the normalized nonlinear grating by $`\chi (z)`$=$`L_d\omega d_{\mathrm{eff}}(z)/(\overline{n}_1c)`$, where $`d_{\mathrm{eff}}`$=$`\chi ^{(2)}/2`$ is given in MKS units. The model (2) describes both temporal and spatial solitons .

The aim is now to average Eqs. (2) and derive accurate equations for the average field. To do so we focus on first order QPM using the conventional square gratings with 50% duty cycle, shown in Fig. 2. We expand the grating functions in Fourier series:

$$\alpha _j(z)=a_j\underset{n}{}g_n\mathrm{e}^{in\kappa z},\chi (z)=d_0+d\underset{n}{}g_n\mathrm{e}^{in\kappa z},$$
(3)
where $`g_n`$=$`2s/(i\pi n)`$ for $`n`$ odd and $`g_n`$=0 for $`n`$ even, with $`s=\mathrm{sign}(\kappa )`$. The gratings drive the system, which means that we may expand the envelope functions in Fourier series also:

$$w=\underset{n}{}w_n(z,x)\mathrm{e}^{in\kappa z},v=\underset{n}{}v_n(z,x)\mathrm{e}^{i(n\kappa \stackrel{~}{\beta })z},$$
(4)
assuming that the coefficients $`w_n(z,x)`$ and $`v_n(z,x)`$ vary slowly in $`z`$ compared to $`\mathrm{exp}(i\kappa z)`$. The residual mismatch $`\stackrel{~}{\beta }`$=$`\beta \kappa `$ is ideally zero.

Three physical length scales are in play: the diffraction length $`L_d`$, the coherence length $`L_c`$, and the grating domain length $`L_0`$. In normalized units $`L_d`$=1, $`L_c`$=$`\pi /|\beta |`$, and $`L_0`$=$`\pi /|\kappa |`$. We assume a typical QPM grating with a domain length that is much shorter than the diffraction length, $`L_0`$$``$1. Furthermore, the grating is of good quality, with the domain length being close to the coherence length, $`L_0`$$``$$`L_c`$, so the residual mismatch is small, $`|\stackrel{~}{\beta }|`$$``$$`|\kappa |`$. In this case $`|\kappa |`$$``$1 and we can use perturbation theory with the small parameter $`ϵ`$=$`1/|\kappa |1`$.

Following the approach of Ref. , we insert the Fourier expansions (3) and (4) into the dynamical equations and assume the harmonics $`w_{n0}`$ and $`v_{n0}`$ to be of order $`ϵ`$. To lowest order ($`ϵ^1`$), this gives the harmonics

$`w_{n0}`$ $`=`$ $`(a_1g_nw_0+dg_{n1}w_0^{}v_0)/(n\kappa ),`$ (5)
$`v_{n0}`$ $`=`$ $`(2a_2g_nv_0+dg_{n+1}w_0^2)/(n\kappa ).`$ (6)
Using these solutions, we obtain to first order $`ϵ`$ the averaged equations for the DC components $`w_0`$ and $`v_0`$:

$`i{\displaystyle \frac{w_0}{z}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2w_0}{x^2}}+\rho w_0^{}v_0+\gamma (|v_0|^2|w_0|^2)w_0=0,`$ (7)
$`i{\displaystyle \frac{v_0}{z}}+{\displaystyle \frac{1}{4}}{\displaystyle \frac{^2v_0}{x^2}}+\stackrel{~}{\beta }v_0+\rho ^{}w_0^2+2\gamma |w_0|^2v_0=0.`$ (8)
These equations also describe $`m`$th order QPM (where $`\stackrel{~}{\beta }`$=$`\beta m\kappa `$ is ideally zero) and any other type of periodic grating, the parameters $`\rho `$ and $`\gamma `$ being simply given as sums over the Fourier coefficients of the grating . Incorporating time or the spatial $`y`$ coordinate is also straightforward. For the square grating (3), $`\rho `$ and $`\gamma `$ can be explicitly calculated to

$$\rho =i\frac{2d}{s\pi }+i\frac{4d_0(a_1a_2)}{s\pi \kappa },\gamma =\frac{d_0^2+d^2(18/\pi ^2)}{\kappa }.$$
(9)
From Eqs. (8) follows the important result that cubic nonlinearities are induced in QNPCs by nonlinear QPM gratings. This cubic nonlinearity has the form of self- and cross-phase modulation (SPM and XPM), and is a result of non-phase-matched coupling between the wave at the main spatial frequency $`\kappa `$ and its higher harmonics. It is thus of a fundamentally different nature than the material Kerr nonlinearity, which is reflected in the fact that the SPM term is absent for the SH.

The averaged model (8) is identical to the known model for nonlinear QPM gratings with no DC component ($`a_j`$=$`d_0`$=0), for which simulations confirmed that bright solitons had properties that were not predicted by the conventional model with only quadratic terms, but were accurately described with the inclusion of the cubic terms . It was further shown that the induced cubic nonlinearity affects the phase modulation of cw waves, enabling efficient switching , and that its strength can be increased by modulating the grating .

For the more general QNPCs considered here, the induced cubic nonlinearity depends on both the DC part and the modulation part of the nonlinear grating, but is independent of the linear grating. The cubic terms may lead to either a focusing or a defocusing effect, depending on the relative intensity of the fields and the sign of the phase mismatch $`\beta `$, since $`\mathrm{sign}(\kappa )=\mathrm{sign}(\beta )`$. The strength of the effective $`\chi ^{(2)}`$ nonlinearity depends on the difference in the linear grating strengths at the fundamental and SH frequencies and on the DC component of the nonlinear grating. We thus recover the well-known effect that the interplay between the linear and nonlinear gratings can increase or decrease the effective $`\chi ^{(2)}`$ nonlinearity, depending on the physical situation .

The averaged model (8) has stationary, localised soliton solutions of the form $`w_0(x,z)`$=$`\mathrm{e}^{i\lambda z}\stackrel{~}{w}(x)/|\rho |`$ and $`v_0(x,z)`$=$`\mathrm{e}^{2i\lambda z}\stackrel{~}{v}(x)/\rho `$, which obey the equations

$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2\stackrel{~}{w}}{x^2}}\lambda \stackrel{~}{w}+\stackrel{~}{w}\stackrel{~}{v}+\stackrel{~}{\gamma }(\stackrel{~}{v}^2\stackrel{~}{w}^2)\stackrel{~}{w}=0,`$ (10)
$`{\displaystyle \frac{1}{4}}{\displaystyle \frac{^2\stackrel{~}{v}}{x^2}}+(\stackrel{~}{\beta }2\lambda )\stackrel{~}{v}+\stackrel{~}{w}^2+2\stackrel{~}{\gamma }\stackrel{~}{w}^2\stackrel{~}{v}=0,`$ (11)
where $`\stackrel{~}{\gamma }`$=$`\gamma /|\rho |^2`$ depends only on $`\kappa `$, $`a_1a_2`$, and $`d_0/d`$. The slowly varying approximation gives valid solutions when the soliton period is longer than the grating period, i.e. when the soliton parameter $`\lambda `$ is small, $`|\lambda ||\kappa |`$.

It is important to stress that Eqs. (11) cover a much more general situation than in , which only considered a nonlinear grating. A given value of the parameter $`\stackrel{~}{\gamma }`$ represents a range of physical situations with different combinations of linear and nonlinear gratings. In Fig. 3 we illustrate representative combinations for exact phase matching ($`\stackrel{~}{\beta }`$=0) that all give the same value of $`\stackrel{~}{\gamma }`$.

The first, simple case ($`f`$=1) is typical for domain inversion in ferroelectric materials, such as LiNbO3. It has only a nonlinear grating that flips from positive to negative with no DC component, $`d_0`$=0. The second case has a DC level of $`\chi ^{(2)}`$ nonlinearity, $`d_0/d`$=3, but no linear grating, corresponding to the nonlinear part of the LiNbO3/H:LiNbO3 structure reported in . The third case is the GaAs/GaAlAs structure reported in , which has a nonlinear grating with $`d_0/d`$=5/3 and a linear grating with $`(a_1a_2)/\beta `$=$`0.07`$. For the second and third cases to give the same $`\stackrel{~}{\gamma }`$ as the first, the grating wavenumber $`\kappa `$ must be multiplied by $`f`$=48.5 and $`f`$=26.6, respectively. Physically this can be done by changing the input beam width $`x_0`$, thus maintaining $`\stackrel{~}{\beta }`$=0.

We have numerically found the bright soliton solutions of Eqs. (11) using a standard relaxation technique. Bright soliton properties were investigated in , but not for normal dispersion, $`\stackrel{~}{\gamma }>0`$. In Fig. 4 we present the properties for normal ($`\stackrel{~}{\gamma }`$=0.02) and anomalous ($`\stackrel{~}{\gamma }`$=$`0.02`$) dispersion, together with the zeroth-order solution ($`\stackrel{~}{\gamma }`$=0). The ratio $`R=\stackrel{~}{v}^2(0)/\stackrel{~}{w}^2(0)`$ of peak intensities, shown in Fig. 4(a), confirms that the zeroth-order approximation becomes increasingly inaccurate for large $`\lambda `$. Also, for a given $`\stackrel{~}{\gamma }`$, $`R`$ approaches the same limiting value as $`\lambda `$ increases, regardless of the value of $`\stackrel{~}{\beta }`$. In this limit the SH is stronger than the fundamental for $`\stackrel{~}{\gamma }>0`$ ($`R>1`$) and much weaker for $`\stackrel{~}{\gamma }<0`$ ($`R0`$). The total power $`P`$=$`_{\mathrm{}}^{\mathrm{}}(\stackrel{~}{v}^2+\stackrel{~}{w}^2)𝑑x`$ is shown in Fig. 4(b). For $`\beta >0`$ this reveals the interesting property that the power threshold for existence, decreases for $`\stackrel{~}{\gamma }>0`$ and increases for $`\stackrel{~}{\gamma }<0`$, as compared to the zeroth-order value.

The bright soliton solutions of the average model (8) were tested for the three QNPCs of Fig. 3 by mapping them back to the variables $`w`$ and $`v`$, and using them as initial conditions in simulations of the field Eqs. (2). The evolution consists of small, regular oscillations superimposed on the slow average beam. Properties of the propagating solitons were calculated by averaging over an integer number of grating periods and were then compared with the predictions of the average model. Figure 5 displays the ratio of peak intensities versus $`\lambda `$ for exact phase matching, $`\stackrel{~}{\beta }`$=0, and reveals that, for both anomalous and normal dispersion, the solutions of the average model are accurate for small $`\lambda `$ and large $`|\kappa |`$, as expected. Even when $`|\kappa |`$=4 the first-order solutions provide a much better fit than the zeroth-order solutions. Our analysis thus shows that bright solitons exist and propagate stably in QNPCs with many types and combinations of linear and nonlinear QPM gratings.

However, this is provided one is careful and does not eliminate the effective $`\chi ^{(2)}`$ nonlinearity by using a grating with $`\kappa `$=$`2(a_2a_1)d_0/d`$. This could happen in realistic QNPCs without violating the assumption $`|\kappa |1`$. For the LiNbO3 and GaP/AIP structures given in  at phase matching ($`\kappa `$=$`\beta `$), the presence of the linear grating changes the effective $`\chi ^{(2)}`$ nonlinearity by a factor of $``$=$`1+2d_0(a_1a_2)/(d\kappa )`$=1.4 and 0.3, respectively. The linear grating thus adds constructively in the LiNb03 structure and destructively in GaP/AIP. In fact, modifying the nonlinear grating in the GaP/AIP structure slightly to $`\chi _a^{(2)}`$=40pm/V (max.) and $`\chi _b^{(2)}`$=19pm/V (min.) eliminates the effective $`\chi ^{(2)}`$ nonlinearity entirely.

Conventional average models  would predict that no soliton could exist with no nonlinearity, $`\rho `$=0. However, in the model (8) the induced cubic nonlinearity predicts that solitons should still exist as solutions of nonlinear Schrödinger equations. In the case when the SH is strong ($`v_0/w_0=\sqrt{5}\mathrm{exp}[i(\stackrel{~}{\beta }\lambda /2)z]`$), a family of bright solitons, $`w_0`$=$`\sqrt{\lambda /(2\gamma )}\mathrm{sech}(\sqrt{2\lambda }x)\mathrm{e}^{i\lambda z}`$, exists for normal dispersion ($`\gamma ,\lambda `$$`>`$0), and a family of dark solitons, $`w_0=\sqrt{\lambda /(4\gamma )}\mathrm{tanh}(\sqrt{|\lambda |}x)\mathrm{e}^{i\lambda z}`$, exists for anomalous dispersion ($`\gamma ,\lambda `$$`<`$0). With no SH, bright solitons, $`w_0=\sqrt{2\lambda /|\gamma |}\mathrm{sech}(\sqrt{2\lambda }x)\mathrm{e}^{i\lambda z}`$, exist for anomalous dispersion ($`\gamma `$$`<`$0,$`\lambda `$$`>`$0), whereas dark solitons, $`w_0=\sqrt{|\lambda |/\gamma }\mathrm{tanh}(\sqrt{|\lambda |}x)\mathrm{e}^{i\lambda z}`$, exist for normal dispersion ($`\gamma `$$`>`$0,$`\lambda `$$`<`$0). We test these solutions by mapping them back to the variables $`w`$ and $`v`$ and launching them as initial conditions in simulations of the field Eqs. (2).

Figure 6 shows the evolution of a bright and dark soliton with no SH over a distance of 25 grating periods. The SH displays small, regular oscillations around the mean value zero, corresponding to the oscillations of the fundamental seen in Fig. 6. The simulations thus confirm that bright and dark solitons can indeed propagate in a stable manner in QNPCs with no effective $`\chi ^{(2)}`$ nonlinearity.

In summary, we have shown that bright solitons exist and propagate in a stable manner in 1D quadratic nonlinear photonic crystals (QNPCs) with many types and combinations of linear and nonlinear QPM gratings. By deriving first-order averaged equations, we have shown that such QNPCs have an induced cubic nonlinearity, and we have numerically found previously unknown families of bright solitons. Even with no effective quadratic nonlinearity, the QNPCs support both bright and dark solitons due to the induced cubic nonlinearity. We have found analytical expressions for these solitons and shown that they propagate in a stable manner.

Dark solitons are always unstable in homogeneous $`\chi ^{(2)}`$ media in settings for type I SHG, due to modulational instability of the back-ground plane-waves . Our results show, for the first time, a dark soliton that appears to be stable under propagation, with the stabilizing mechanism necessarily originating from the photonic crystal structure of the QNPC. An analysis of dark solitons in the general case is being carried out.

The project is supported by the Danish Technical Research Council through Talent Grant No. 9800400.