The quirk trajectory

Abstract

We for the first time obtain the analytical solution for the quirk equation of motion in an approximate way. Based on it, we study several features of quirk trajectory in a more precise way, including quirk oscillation amplitude, number of periods, as well as the thickness of quirk pair plane. Moreover, we find an exceptional case where the quirk crosses at least one of the tracking layers repeatedly. Finally, we consider the effects of ionization energy loss and fixed direction of infracolor string for a few existing searches.

1 Introduction

Solutions to the gauge hierarchy problem of the Standard Model (SM) of particle physics such as supersymmetry and composite Higgs models usually predict a colored top partner with mass around TeV scale. They have been challenged by the null results of LHC searches so far. Theories of neutral naturalness [1] aim to address the gauge hierarchy problem without introducing colored states, thus relieve the tension with the LHC searches. This class of models include folded supersymmetry [2, 3], quirky little Higgs [4], twin Higgs [5,6,7], minimal neutral naturalness model [8] and so on. In those models, some new SU(N) gauge symmetries are introduced in addition to the SM gauge group. A particle, which is charged under both the SM electroweak gauge group and the new confining SU(N) gauge group and has mass much larger than the confinement scale ($\Lambda$) of the SU(N), is dubbed as quirk. At colliders, the quirk can only be produced in pairs due to the conserved SU(N) symmetry. The infracolor force ($F_{s}$, interaction induced by the SU(N) gauge bosons) between two quirks will lead to non-conventional signals in the detector. The manifestation of the quirk signal is strongly dependent on $\Lambda$ due to $F_s \propto \Lambda ^2$ [9].

Throughout the work, we focus on quirk with mass around the EW scale, motivated by the gauge hierarchy problem. For $\Lambda \gtrsim \mathcal {O}(10)$ MeV, the strong infracolor force will lead to intensive oscillations in quirk motion. The quirk-pair system will lose kinetic energy quickly via photon and hidden glueball radiation. Subsequently, they will annihilate almost promptly into the SM particles after production. Such quirk signals can be searched through resonances in the SM final states [10,11,12,13,14,15]. When $\Lambda \in [10~\text {keV}$, $10~\text {MeV}]$, the quirk pair oscillation amplitude is microscopic ($\lesssim \upmu$m), while the energy loss due to photon and hidden glueball radiation is not efficient. The electric neutral quirk-pair system will leave a straight line inside the tracker, which will be reconstructed as a single ultra-boosted charged particle with a high ionization energy loss (different from conventional heavy stable charged particle). This signal was looked for at the Tevatron [16]. As for $\Lambda \lesssim \mathcal {O}(10)$ eV, the infracolor force is too small to push the quirk out of its helical trajectory, because finite spacial resolution is considered and typical track reconstruction allows the $\chi ^2/\text {DOF}$ in fitting as large as 5 [17]. Then the trajectory of each quirk can be reconstructed as normal track in detector. The signal can be constrained by conventional heavy stable charged particle searches at the LHC [17,18,19].

The scenario with $\Lambda \in [100~\text {eV}$, $10~\text {keV}]$ is most interesting, since the quirk pair oscillation amplitude can be macroscopic ($\sim$ cm) and the tracker can resolve hits of two quirks on each tracking layer. In this case, the quirk trajectory is no longer a helix. The hits caused by the quirks on tracking layers will be completely ignored in conventional event reconstruction at the LHC. Meanwhile, as we will show later, the quirk energy deposit in electromagnetic/hadronic calorimeter (ECal/HCal) is usually small within the time period of bunch crossing (25 ns) at the LHC. As a result, the quirk simply behaves as missing transverse energy and will be constrained by mono-jet searches at the LHC [19,20,21] if the quirk pair is produced recoiling against an energetic initial state radiated (ISR) jet. In fact, there are still many features of quirk trajectory that we could use to help identifying the quirk signal in detector. Ref. [22] innovatively pointed out that if the quirk-pair system is relatively boosted and the infracolor force is much stronger than Lorentz force ($B=4$ T inside the CMS detector), the hits of quirk pair in the detector will almost lie on a plane with deviation less than $\sim \mathcal {O}(100)~\upmu$m. They found that the coplanar hits search suffers little background while maintaining very high signal efficiency. What is more, the quirk carries electric charge and moves slowly when two quirks are widely separated (the typical velocity of quirk-pair system is $\sim 0.1$ for quirk pair transverse momentum $p_T>100$ GeV and quirk mass ${m} \sim 100$ GeV). The information of relatively large ionization energy loss at each hit in the tracker [23] can be used to further improve the coplanar hits search [24]. On the other hand, if the kinetic energy of the quirk-pair system is small, the ECal/HCal of detector will be able to stop the quirk pair. Then, after a long time oscillation inside the calorimeter, the quirk pair annihilates eventually. If the annihilation happens at a time when there are no active pp collisions, the signal would be captured by stopped long-lived particles searches at the LHC [25,26,27].

In this paper we will present an improved understanding of the quirk trajectory inside detector. In particular, we provide an approximate analytical solution for the quirk equation of motion (EoM) for the first time. Based on this solution, we can obtain a more precise expression for the quirk pair oscillation amplitude as well as the number of periods of quirk’s motion inside the detector. The thickness of quirk pair plane will be discussed in details. Its dependence on kinetic variables, the confinement scale as well as the quirk charge will be given. We also find that there is great possibility for each quirk traveling though the same tracking layer more than once, leading to multiple resolvable hits on the layer. Note that the coplanar search proposed in Ref. [22] becomes less efficient in this situation. Finally, we briefly discuss how the existing searches for the quirk pair with macroscopic oscillation amplitude will change if one includes the ionization energy loss in solving quirk EoM or assumes the direction of infracolor string is fixed (The latter assumption seems to be taken in Ref. [25]).

In the following, we will consider the infra-gauge group as SU(3). Without specification, the quirk ($\mathcal {Q}$) is fermion with SM quantum numbers (3, 1, 2/3). The colored quirk contradicts the principle of neutral naturalness, but it has large production rate at the LHC and may manifest itself first. The SM quantum numbers only affect the production channel of the quirk, which will lead to different distributions of quirk initial momenta. The analysis proposed in this work is fully applicable to quirks with other quantum numbers. Note that due to color confinement, only the quirk-quark bound state is observable in experiment and its electric charge can be either $\pm 1$ or zero. It was found by Pythia8 [28] simulation that around 30% [22] of quirk-quark bound states have charge $\pm 1$. Only the quirk bound states with non-zero electric charge are considered in this work. In our Monte Carlo simulation, the events are generated with MG5_aMC@NLO [29] framework, where the simplified quirk model is written in UFO format by FeynRules [30], the parton shower and hadronization of the ISR jet are implemented by Pythia8. However, the QCD parton shower and hadronization of colored quirk are ignored, since they will not significantly change the kinetic energy of the quirk. We adopt the CMS detector configuration, simulated the ionization energy loss of both SM particles and the quirk inside the detector, as well as including the pile-up events. The detailed introduction about the simulation can be found in our previous work [24].

The paper is organized as following. The Sect. 2 is devoted to solve the quirk EoM analytically, either with or without including the external forces. The oscillation amplitude, period number as well as the plane thickness of quirk pair trajectory in the tracker are discussed in the Sects. 3, 4 and 5, respectively. In Sect. 6, we find an exceptional case where the quirk crosses at least one of the tracking layers repeatedly. We study the effects of ionization energy loss and fixed ${\hat{s}}$ in obtaining quirk bounds from existing searches respectively in the Sects. 7 and 8, and conclude our work in Sect. 9. Moreover, the technical details are given in Appendices A, B, C and D.

2 Quirk motions inside detector

The quirk equation of motion (EoM) inside detector is given by [9]

$$\begin{aligned} \frac{\partial ({m} \gamma \vec {v})}{\partial t}&=\vec {F}_{s}+\vec {F}_{ext}, \end{aligned}$$

(2.1)

$$\begin{aligned} \vec {F}_{s}&=-\Lambda ^2\sqrt{1-\vec {v}_{\perp }^{2}} {\hat{s}}-\Lambda ^2 \frac{v_{ \Vert } \vec {v}_{\perp }}{\sqrt{1-\vec {v}_{\perp }^{2}}}, \end{aligned}$$

(2.2)

$$\begin{aligned} \vec {F}_{\text {ext}}&=q\vec {v}\times \vec {B}-\left\langle \frac{dE}{dx}\right\rangle {\hat{v}}, \end{aligned}$$

(2.3)

where $\gamma =1/\sqrt{1-\vec {v}^2}$, $v_{ \Vert }=\vec {v}\cdot {\hat{s}}\nonumber$ and $\vec {v}_{\perp }=\vec {v}-v_{ \Vert }{\hat{s}}\nonumber$ with ${\hat{s}}$ being a unit vector along the string pointing outward at the endpoints. $\vec {F}_s$ corresponds to the infracolor force and is described by the Nambu–Goto action, where $\Lambda$ is the confinement scale. $\vec {F}_{\text {ext}}$ represents the external forces including Lorentz force and the effects of ionization energy loss for charged quirk propagating in magnetic field and through materials, respectively. Note that we have ignored several sub-dominating energy loss effects such as infracolor glueball and photon radiation, as well as hadronic interaction with detector.

To solve Eq. 2.1, we have to consider both the quirk pair centre of mass (CoM) frame and the lab frame. In the CoM frame, ${\hat{s}}$ is approximately parallel to the vector difference between positions of the two quirks (this is only true for $\Lambda ^2 \gg F_{\text {ext}}$, see Ref. [9]). However, the CoM frame itself is changing all the time due to effects of $\vec {F}_{\text {ext}}$, which is related to quirk velocity in the lab frame. The procedures of numerically solving the EoM by slowly increasing the time with small steps were introduced in Ref. [24]. In what follows, we will provide a detailed analytical solution for the quirk EoM which is found to be accurate in a wide range of parameter space.

With $|\vec {F}_{\text {ext}}|\ll \Lambda ^2$ (which corresponds to $\Lambda > \mathcal {O}(100)$ eV in the CMS detector), it is reasonable to first ignore $\vec {F}_{\text {ext}}$ if we do not plan to collect the information of ionization energy loss, and take the contribution of $\vec {F}_{\text {ext}}$ as a correction.

2.1 Kinematics without $\vec {F}_{\text {ext}}$

Taking $\vec {F}_{\text {ext}} =0$ in Eq. 2.1, the trajectories of two quirks will lie exactly on the plane (denoted by A ) constructed by $\vec {P}_{o1}$ and $\vec {P}_{o2}$, with $\vec {P}_{oi}$ ($E_{oi}$) corresponding to the ith quirk initial momentum (energy) in the lab frame. Denoting ${\hat{e}}_x$ and ${\hat{e}}_y$ as the unit vectors of $\left( \vec {P}_{o1}/{E_{o1}}-\vec {P}_{o2}/{E_{o2}}\right)$ and $\left( \vec {P}_{o1}-\vec {P}_{o1}\cdot {\hat{e}}_x {\hat{e}}_x\right)$ respectively, the ${\hat{s}}_i~(i=1,2)$ in Eq. 2.2 is the same as either of $\pm {\hat{e}}_x$ all the time and reverses for each time when two quirks cross each other during the oscillation. We can use $\rho _i=\frac{\vec {P}_{oi}\cdot {\hat{e}}_x}{{m}}$, $h=\frac{\vec {P}_{o1}\cdot {\hat{e}}_y}{E_{o1}}$ and $l=\frac{\sqrt{1-h^2}\Lambda ^2t}{{m}(\rho _1-\rho _2)}$ to define

$$\begin{aligned} \vec {r}_{o1}(l)=\,&\frac{{m}}{\Lambda ^2}\left( \sqrt{1+\rho _1^2}-\sqrt{1+\left( \rho _1-l(\rho _1-\rho _2) \right) ^2}\right) {\hat{e}}_x\nonumber \\&+ \frac{{m}}{\Lambda ^2}\frac{h(\rho _1-\rho _2)l}{\sqrt{1-h^2}}{\hat{e}}_y, \end{aligned}$$

(2.4)

$$\begin{aligned} \vec {r}_{o2}(l)=\,&\frac{{m}}{\Lambda ^2}\left( \sqrt{1+\left( \rho _2-l(\rho _2-\rho _1) \right) ^2}- \sqrt{1+\rho _2^2}\right) {\hat{e}}_x\nonumber \\&+\frac{{m}}{\Lambda ^2}\frac{h(\rho _1-\rho _2)l}{\sqrt{1-h^2}}{\hat{e}}_y. \end{aligned}$$

(2.5)

According to the Eq. 2.1, the trajectories of two quirks in the lab frame will be

$$\begin{aligned} \vec {r}_{1}(l)&=[l]\vec {r}_{o1}(1)+\frac{1+(-1)^{[l]}}{2}\vec {r}_{o1}({\bar{l}})+\frac{1-(-1)^{[l]}}{2}\vec {r}_{o2}({\bar{l}}), \end{aligned}$$

(2.6)

$$\begin{aligned} \vec {r}_{2}(l)&=[l]\vec {r}_{o2}(1)+\frac{1-(-1)^{[l]}}{2}\vec {r}_{o1}({\bar{l}})+\frac{1+(-1)^{[l]}}{2}\vec {r}_{o2}({\bar{l}}), \end{aligned}$$

(2.7)

where [l] is floor of l and ${\bar{l}} \equiv l-[l]$. It is straight forward to find that $\vec {r}_{1}(l)=\vec {r}_{2}(l)$ when l is an integer and the time interval between l and $l+1$ is

$$\begin{aligned} {t_p=0.0658\frac{(\rho _1-\rho _2)}{\sqrt{1-h^2}}\frac{{m}}{\text {[100~GeV]}}\frac{\text {[keV]}^2}{\Lambda ^2}~\text {[ns]}}. \end{aligned}$$

(2.8)

2.2 Kinematics with uniform magnetic field

With an uniform magnetic field $\vec {B}$ in the lab frame, the Lorentz forces on two quirks will push them out of the plane A. Then the motions of two quirks can be decomposed into two parts: one is parallel to plane A and the other is perpendicular to plane A. In the limit of $F_{\text {Lorentz}} \ll \Lambda ^2$, the motions parallel to plane A can be approximated by Eqs. 2.6 and 2.7. By adding the motion perpendicular to plane A, the trajectories of two quirks in the lab frame can be expressed as

$$\begin{aligned} \vec {r}_{i}^{~\prime }(l)&=\vec {r}_{i}(l)+z_i(l){\hat{e}}_z, \end{aligned}$$

(2.9)

where ${\hat{e}}_z={\hat{e}}_x\times {\hat{e}}_y$ with ${\hat{e}}_x$ and ${\hat{e}}_y$ defined as before.

It will be easier to calculate $z_i$ in the approximately invariant CoM frame. Each of the initial quirks momenta ($\vec {P}_{oi}^\prime$) in the CoM frame is given by $|\vec {P}_{oi}^\prime |={m}\rho$ with

$$\begin{aligned}&\rho =\sqrt{\frac{(E_{o1}+E_{o2})^2-(\vec {P}_{o1}+\vec {P}_{o2})^2}{4{m}^2}-1}. \end{aligned}$$

(2.10)

The quirks motions described by Eqs. 2.6 and 2.7 in the lab frame can be boosted into the CoM frame and expressed as

$$\begin{aligned} \vec {v}_{c1}(g)&=-\vec {v}_{c2}(g)=(-1)^{[g]} v({\bar{g}}){\hat{e}}_{xc}, \end{aligned}$$

(2.11)

$$\begin{aligned} \vec {r}_{c1}(g)&=-\vec {r}_{c2}(g)=(-1)^{[g]} r({\bar{g}}){\hat{e}}_{xc}, \end{aligned}$$

(2.12)

$$\begin{aligned} v(g)&= \frac{\rho (1-2g)}{\sqrt{1+\rho ^2(1-2g)^2}}, \end{aligned}$$

(2.13)

$$\begin{aligned} r(g)&=\frac{{m}}{\Lambda ^2}\left( \sqrt{1+\rho ^2}-\sqrt{1+\rho ^2(1-2g)^2} \right) , \end{aligned}$$

(2.14)

where ${\hat{e}}_{xc}$ is the unit vector of $\vec {P}_{o1}^\prime$, ${\hat{e}}_{yc}={\hat{e}}_{z}\times {\hat{e}}_{xc}$, $g=\frac{\Lambda ^2 t_c}{2{m}\rho }$, ${\bar{g}}=g-[g]$, and $t_c$ is time in the CoM frame. Note that we have used subscript c for variables in the CoM frame.

Since the velocities of both quirks are approximately along $\pm {\hat{e}}_{xc}$ in the CoM frame, the ${\hat{e}}_{z}$-component of $\vec {F}_{\text {ext}}$ is induced by the ${\hat{e}}_{yc}$-component of the magnetic field and the ${\hat{e}}_{z}$-component of the electric field, which are

$$\begin{aligned} B_c&=\frac{\left( {1-\beta _o^2}\right) \vec {P}_{o1}\cdot {\hat{e}}_{x\prime }\vec {B}\cdot {\hat{e}}_{y\prime }-{\left( \vec {P}_{o1}\cdot {\hat{e}}_{y\prime }-\beta _o E_{o1}\right) \vec {B}\cdot {\hat{e}}_{x\prime }} }{{m}\rho \left( {1-\beta _o^2}\right) }, \end{aligned}$$

(2.15)

$$\begin{aligned} E_z&=-\frac{\beta _o \vec {B}\cdot {\hat{e}}_{x\prime } }{\sqrt{1-\beta _o^2}}, \end{aligned}$$

(2.16)

where $\vec {\beta }_o=\frac{\vec {P}_{o1}+\vec {P}_{o2}}{E_{o1}+E_{o2}}$, $\beta _o=|\vec {\beta _o}|$, ${\hat{e}}_{y\prime }$ and ${\hat{e}}_{x\prime }$ are the unit vectors of $\vec {\beta }_o$ and $( \vec {P}_{o1}-\vec {P}_{o1}\cdot {\hat{e}}_{y\prime }{\hat{e}}_{y\prime })$, respectively.

Defining

$$\begin{aligned} z^{+}(g)&\equiv \frac{1}{2}\left( z_1(g)+z_2(g) \right) , \end{aligned}$$

(2.17)

$$\begin{aligned} z^{-}(g)&\equiv \frac{1}{2}\left( z_1(g)-z_2(g) \right) , \end{aligned}$$

(2.18)

because electric charges of two quirks have opposite sign and their velocities satisfy Eq. 2.11, we can infer that $z^{+}(g)$ arises from the Lorentz force caused by $B_c$, and $z^{-}(g)$ arises from the electric field force caused by $E_z$. The detailed derivations for $z^{+}(g)$ and $z^{-}(g)$ are given in Appendix A. Considering $t_{ci}=\frac{t_i-\vec {\beta }_o \cdot \vec {r}_i(l_i)}{\sqrt{1-\beta _o^2}}$, we have

$$\begin{aligned} g_i(l)&=\frac{1}{\sqrt{1-\beta _o^2}}\left( \frac{(\rho _1-\rho _2)}{2\rho \sqrt{1-h^2}}l-\frac{\Lambda ^2\vec {r}_i(l)\cdot \vec {\beta }_o}{2{m}\rho }\right) . \end{aligned}$$

(2.19)

The trajectories of two quirks in the lab frame can be expressed as

$$\begin{aligned} \vec {r}_{1}^{~\prime }(l)&=\vec {r}_{1}(l)+\left( z^{+}(g_1(l))+z^{-}(g_1(l)) \right) {\hat{e}}_z, \end{aligned}$$

(2.20)

$$\begin{aligned} \vec {r}_{2}^{~\prime }(l)&=\vec {r}_{2}(l)+\left( z^{+}(g_2(l))-z^{-}(g_2(l)) \right) {\hat{e}}_z. \end{aligned}$$

(2.21)

2.3 Variable substitution and benchmark point

We have derived the quirk trajectories in the lab frame in terms of $\vec {B}$, m, $\Lambda$, $\vec {P}_{o1}$ and $\vec {P}_{o2}$. In the following discussion, we will take the configuration of the CMS detector whenever discussing the experimental measurements. This corresponds to $\vec {B}=4$ T along the ${\hat{z}}$-axis. The $\vec {P}_{o1}$ and $\vec {P}_{o2}$ parameters will be replaced by $k_i=|\vec {P}_{oi}/m|~(i=1,2)$, $\alpha$, $\theta$ and $\phi$. The $\alpha$ is the angle between $\vec {P}_{o1}$ and $\vec {P}_{o2}$. The $\theta$ stands for the angle between $\vec {P}_{o}=\vec {P}_{o1}+\vec {P}_{o2}$ and ${\hat{z}}$. The $\phi$ is related to the angle between $\vec {P}_{o1} \times \vec {P}_{o2}$ and $\vec {P}_{o} \times {\hat{z}}$. To be specific, we have

$$\begin{aligned} \theta&={\left\{ \begin{array}{ll} \cos ^{-1}(\delta ),&{}\quad \delta \ge 0\\ \pi -\cos ^{-1}(\delta ),&{}\quad \delta <0 \end{array}\right. },\end{aligned}$$

(2.22)

$$\begin{aligned} \phi&={\left\{ \begin{array}{ll} \pi /2-\cos ^{-1}(\psi ),&{}\quad \delta \ge 0\\ \cos ^{-1}(\psi )-\pi /2,&{}\quad \delta <0 \end{array}\right. }, \end{aligned}$$

(2.23)

where $\delta =\frac{\vec {P}_o\cdot {\hat{z}}}{|\vec {P}_o|}$ and $\psi =\frac{\left( \vec {P}_{o1}\times \vec {P}_{o2}\right) \cdot \left( \vec {P}_{o}\times {\hat{z}} \right) }{|\vec {P}_{o1}\times \vec {P}_{o2}||\vec {P}_{o}\times {\hat{z}}|}$.

Fig. 1

Distributions of quirk initial kinematic variables for different quirk masses and quirk-pair system transverse momenta. In the upper-left panel, the solid and dashed lines correspond to $k_1$ and $k_2$, respectively

Figure 1 shows the distributions of $k_{1,2}$, $\alpha$, $\theta$ and $\phi$ with different quirk masses and quirk-pair system transverse momenta in our benchmark model. Note that the quirk-pair system is required to be relatively boosted along the transverse direction, in order to produce detectable signals inside the tracker. This is implemented by requiring a hard initial state radiated (ISR) jet that is recoiling against the quirk pair. Thus, the events of quirk production in the detector can be triggered by either the relatively large missing transverse energy or a hard jet. Meanwhile, the behaviors of quirk inside the tracker will be also affected by the hardness of the ISR jet. We can find that larger $p_T/{m}$ leads to greater $k_{1,2}$, $\theta$ and $|\phi |$, while rendering smaller $\alpha$.

Since many observables of interests contain complex coefficients, it will be more intuitive to show the results on a benchmark point, which we choose to be

$$\begin{aligned} k_1&=k_2=2.5,\quad \alpha =\pi /4, \quad \theta =0.3, \quad \phi =1.0, \nonumber \\ {\Lambda _0^2}&= 4\times 10^5~\text {eV}^2,\quad m=100~\text {GeV}. \end{aligned}$$

(2.24)

From Fig. 1, we can see that this point has high probability in most cases.

3 Oscillation amplitude in the lab frame

It is known that the quirks are traveling oscillatingly, with the characteristic amplitude of oscillation in the CoM frame [9] $\ell _{c} \sim 2~\text {cm}(\gamma -1) (\frac{{m}}{\text {100~GeV}}) (\frac{\Lambda }{\text {keV}})^{-2}$, which can be expressed via the new parameters as $\ell _{c} \sim 2~\text {cm}(\sqrt{1+\rho ^2} -1) \frac{m}{\text {[100~GeV]}}\frac{\text {[keV]}^2}{\Lambda ^2}$ with $\rho =\sqrt{\frac{\sqrt{1+k_1^2} \sqrt{1+k_2^2} - k_1 k_2 \cos (\alpha ) -1}{2}}$ from Eq. 2.10. The definitions of $k_1$, $k_2$ and $\alpha$ are given in Sect. 2.3. In fact, there is another more useful parameter relevant to the width of the quirk oscillation in the lab frame

$$\begin{aligned} L=2 \frac{\mathcal {R}}{\rho }\ell _{c}, \end{aligned}$$

(3.1)

with $\mathcal {R}=\frac{k_1 k_2 \sin (\alpha )}{\sqrt{k_1^2 +k_2^2 + 2 k_1 k_2 \cos (\alpha )}}$. L corresponds to twice the length of the projection of $\ell _{c}$ onto the plane perpendicular to $\vec {\beta }_o$. It can be calculated immediately that $L=3.8$ cm for our benchmark point.

4 Number of periods in the tracker

Since two quirks are produced at the same interaction point and traveling oscillatingly, we can define each of the second time when two quirks meet as one period of quirk motion, corresponding to the l increased by one in our previous discussion. The overall shape of the tracker system of the CMS detector is cylindrical and described by $R_{max}=118.5$ cm (3.95 ns in natural units) and $|Z_{max}|=293.5$ cm (9.78 ns) [31]. So we can calculate the number of periods ($n_T$) that the quirk pair is going through inside the tracker system in 25 ns (bunch crossing time at the LHC)

$$\begin{aligned} n_T&={\left\{ \begin{array}{ll} \frac{\text {Min}[25,~9.78/\beta _z]~\text {ns}}{t_p},\quad \beta _z/\beta _T\ge 2.4768\\ \frac{\text {Min}[25,~3.95/\beta _T]~\text {ns}}{t_p},\quad \beta _z/\beta _T<2.4768 \end{array}\right. }, \end{aligned}$$

(4.1)

where 2.4768=$|Z_{max}|/R_{max}$, $\beta _z=|\vec {\beta }_o\cdot {\hat{z}}|$, $\beta _T=\sqrt{\beta _o^2-\beta _z^2}$ and $t_p$ is given in Eq. 2.8. So we get $n_T \propto \Lambda ^2/m$.

We can find that our benchmark point is going through 19.5 periods inside the tracker. The $n_T$ is insensitive to the $\phi$ parameter. In Fig. 2, distributions of $n_T$ for different quirk masses and quirk-pair transverse momenta and the dependence of $n_T$ on the $k_1-k_2$ and $\theta -\alpha$ are shown. In the $k_1-k_2$ plane, increasing $k_1$ and $k_2$ lead to larger $t_p$ and shorter time of the quirk pair staying in tracker, and thus smaller $n_T$. In the $\theta -\alpha$ plane, $n_T$ decreases with the increasing small $\alpha$ because of the increased $t_p$. When $\alpha$ is not small, increasing $\alpha$ leads to longer time of the quirk pair staying in tracker and thus larger $n_T$. The quirk pair leaves the tracker by crossing the outermost endcap (barrel) when $\theta <(>) 0.384$ ($\tan [0.384]$=$R_{max}/|Z_{max}|$ ). Similarly, different quirk masses and quirk-pair transverse momenta lead to different distributions of $n_T$ due to the differences in the parameter space of $k_1-k_2$ and $\theta -\alpha$.

Fig. 2

Left: distributions of number of periods for different quirk masses and quirk-pair transverse momenta (from numerical simulation). Middle and right: the projected (parameters are set at benchmark points if they are not varying) relation between $n_T$ and $k_1-k_2$, and $\theta -\alpha$ (from analytical calculation). Here, we have fixed $\Lambda =\Lambda _0$ for all cases

5 Thickness of quirk pair plane

Quirk-antiquirk pair propagating through the tracker system of a detector will leave N hits located at $\vec {h}_i~(i=1,2,\ldots ,N)$. It was pointed out firstly in the Ref. [22] that in a wide range of parameter space with $\Lambda \sim \mathcal {O}(1)~\text {keV}$, these hits will largely lie on a plane. The averaged distance of hits to a virtual plane which contains the interaction point (at origin) is calculated by

$$\begin{aligned} d(\vec {n})=\sqrt{\frac{1}{N-1}\sum _{i=1}^{N}\left( \vec {n}\cdot \vec {h}_i\right) ^2}, \end{aligned}$$

(5.1)

where $\vec {n}$ is the normal vector of the plane. The plane giving the smallest $d(\vec {n})=d_{min}$ is called the quirk pair plane and the $d_{min}$ is called thickness of the quirk pair plane.

With an uniform magnetic field $\vec {B}$ in the detector, trajectories of the quirk pair can be described by Eqs. 2.20 and 2.21. The number of periods that the quirk pair is going through inside the tracker system in 25 ns is $n_T\approx l_0$, where $l_0$ is the integer closest to $n_T$.

According to the discussions in Sect. 2.2 and Appendix A, the quirk pair plane in the CoM frame can be approximately obtained by rotating the ${\hat{e}}_{xc}-{\hat{e}}_{yc}$ plane by an angle $(-j\eta )$ around ${\hat{e}}_{yc}$ and then moving it along ${\hat{e}}_{z}$ by a distance $c_z$. Then distances of two quirk trajectories to the plane obtained above are

$$\begin{aligned} d_{c1}(g)=\,&z_1(g)-c_z-(-1)^{[g]}j\eta r({\bar{g}}), \end{aligned}$$

(5.2)

$$\begin{aligned} d_{c2}(g)=\,&z_2(g)-c_z+(-1)^{[g]}j\eta r({\bar{g}}), \end{aligned}$$

(5.3)

respectively, so that thickness of quirk pair plane can be approximately expressed as

$$\begin{aligned} d_c^2&=\frac{1}{2l_0}\int _{0}^{l_0}\left( d_{c1}^2(g)+d_{c2}^2(g) \right) dg\nonumber \\&\approx \frac{1}{l_0}\sum _{i=1}^{l_0}\left( j-[\frac{i}{2}] \right) ^2\int _{0}^{1} \eta ^2r^2(g)dg \nonumber \\&\quad +\int _{0}^{1} \left( z^{+}_o(g)-c_z \right) ^2 dg+\int _{0}^{1} z^{-}_o(g)^2 dg. \end{aligned}$$

(5.4)

Finally, thickness of quirk pair plane in the lab frame can be estimated as (The detailed discussions are provided in Appendix B.)

$$\begin{aligned} d&= \sqrt{ \frac{\eta ^2r^2(0.5)}{3l_0}\sum _{i=1}^{l_0}\left( j(l_0)-[\frac{i}{2}]\right) ^2+\frac{4}{3}\left( \frac{2{z_o^{+ }}^2(1)}{15}+\frac{{z_{o}^{-}}^2(1)}{5}\right) } \end{aligned}$$

(5.5)

$$\begin{aligned}&\approx \sqrt{\left( \frac{1}{3l_0}\sum _{i=1}^{l_0}\left( j(l_0)-[\frac{i}{2}]\right) ^2+\frac{4}{15}\right) d_E^2 +\frac{8}{45} d_B^2}, \end{aligned}$$

(5.6)

where we have used the relations of Eq. A.17 in the second line, and

$$\begin{aligned} j(l_0)&=\frac{1}{l_0}\sum _{i=1}^{l_0}[\frac{i}{2}], \end{aligned}$$

(5.7)

$$\begin{aligned} d_E&={z_{o\prime }^{-}}(1)=\frac{2mqE_z}{\Lambda ^4}\rho \sinh ^{-1}[\rho ], \end{aligned}$$

(5.8)

$$\begin{aligned} d_B&={z_o^{+ }}(1)=\frac{2mqB_c}{\Lambda ^4}\left( \rho -\tan ^{-1}[\rho ] \right) . \end{aligned}$$

(5.9)

We provide the validation of Eq. 5.6 in Appendix C, where we can conclude that our analytic formula for the quirk pair plane thickness matches the numerical result within an order magnitude when quirk pair goes through $\sim \mathcal {O}(1-100)$ periods inside tracker. The main difference between the quirk pair plane thickness calculated from Eq. 5.6 ($D_A$) and that obtained from numerical simulation $(D_N)$ in Appendix C attributes to the following reasons. First of all, the former uses the shapes of overall quirk trajectories, or positions of infinite points from every parts of quirk trajectories. The latter only uses $\mathcal {O}$(10) positions of hits caused by quirk crossing detector layers in tracker. Besides, we employ the relations of Eq. A.17 to obtain Eq. 5.6 as an approximation of Eq. 5.5.

Fig. 3

Thickness of the quirk pair plane in the projected parameter space: $k_1-k_2$ (left), $\theta -\phi$ (middle), $\alpha$ (right). Parameters are set at benchmark points if they are not considered as varying. The confinement scale $\Lambda$ is fixed to be $\Lambda _0$ for all cases

In Fig. 3, we show the projected thickness of the quirk pair plane on the parameter space of $k_1-k_2$, $\theta -\phi$ and $\alpha$, respectively, by using the Eq. 5.6. The irrelevant parameters are set at benchmark point in the projection. For $l_0\sim \mathcal {O}(10)$, $d_E$ plays a more important role in d than $d_B$, since $d\approx \sqrt{\mathcal {O}(1) d_E^2+\frac{8}{45}d_B^2}$. In the $k_1-k_2$ plane, increasing $k_1$ and $k_2$ will lead to increased $E_z$ and $\rho$ thus larger $d_E$. Similarly, in the $\theta -\phi$ plane, larger $\theta$ and $|\phi |$ give greater $E_z$ thus larger $d_E$. The dependence on the $\alpha$ is more complicated. Increasing $\alpha$ will lead to larger $\rho$ parameter while smaller $E_z$. In the small $\alpha$ region, the $\rho$ parameter is dominating. So the quirk pair plane thickness is increased with $\alpha$. On the other hand, the $E_z$ becomes dominant in the large $\alpha$ region, which gives decreased plane thickness for increasing $\alpha$. Note that the non-smooth behavior of the contours are originated from the fact that $l_0$ is not a smooth function of $t_0$.

It will be more useful to predict a model or parameter space which has large thickness of quirk pair plane by using Eq. 5.6 and features in Fig. 3 instead of conducting time consuming numerical simulation. In the left panel of Fig. 4, we plot the distributions of quirk pair plane thickness with varying quirk mass (m), transverse momentum of quirk pair ($p_T$) and confinement scale ($\Lambda$). The dependence on $\Lambda$ is obvious: the thickness decreases with increasing $\Lambda$. However, for the considered quirk production process (dominated by $gg\rightarrow \mathcal {Q} \mathcal {Q}$), greater $p_T/{m}$ gives larger $k_{1,2}$, $\theta$ and $|\phi |$, but smaller $\alpha$ and $l_0$. Thus we can find the thickness dependence on m and $p_T$ is mild (we have checked with several parameter choices which are not shown in the plot.). On the other hand, if the quirk pair is produced from a heavy resonant decay $p p \rightarrow j Z'(\rightarrow \mathcal {Q} \mathcal {Q})$, the quirk pair plane thickness will be much larger for heavier $Z'$.

5.1 Charge dependence of the quirk pair plane thickness

In previous discussions, we have chosen the electric charge of quirks to be $\pm 1$ and the quirk-pair system is electric neutral. However, our conclusion can be naturally applied to the quirks with different charges as long as the quirk pair is keeping neutral. Because it is the Lorentz force rendering the quirk traveling outside the plane, the plane thickness is linearly proportional to each of the quirk charges. This feature is clearly shown by the solid lines in the right panel of Fig. 4, where the thickness is calculated precisely from numerical simulation on the samples with different charges while keeping the initial momentum and $\Lambda$ the same.

The quirk pair plane thickness will be dramatically increased for non-neutral charge quirk-pair system, since the trajectory of the quirk-pair system will be bended by the Lorentz force. In this case our discussion for the quirk pair plane thickness can only be used to roughly estimate the thickness increasement in one period of quirk’s motion. As a result, the total quirk pair plane thickness will be increased more intensely by the number of periods inside the tracker, comparing to the thickness of electric neutral quirk-pair system. For the process with ${m}=100$ GeV, $p_T>100$ GeV, $\Lambda =\Lambda _0$, the number of periods of quirk motion inside tracker is around 20, so the quirk pair with charges $\pm 1/3$ and $\pm 2/3$ should have plane thickness $\sim 20$ times larger than quirk pair with charges $\pm 1/3(2/3)$ and $\mp 1/3(2/3)$, as demonstrated in the right panel of Fig. 4.

Fig. 4

Left: quirk pair plane thickness (calculated from numerical simulation instead of using Eq. 5.6) with varying quirk mass (m), transverse momentum of quirk pair ($p_T$) and confinement scale ($\Lambda$). Right: quirk pair plane thickness for different quirk charges, where we choose ${m}=100$ GeV, $p_T>100$ GeV, $\Lambda =\Lambda _0$

Table 1 The fraction of quirk pair events that leave at most n hits on each tracking layer of the CMS detector is denoted by F$_{n^{\max }_{\text {hit}}=n}$. The different processes are characterized by fixing quirk mass m, confinement scale $\Lambda$ and least transverse momentum of quirk-pair system ($p_T$)

6 Crossing the same tracking layer more than once

The coplanar search proposed in Ref. [22] is designed to be based on the assumption that the quirk pair induces two hits on each barrel tracking layer. However, the quirk pair trajectories are highly dependent on the quirk mass, confinement scale as well as initial momentum. In many cases, the quirk pair may induce greater number of hits on tracking layers. Those intensive coplanar hits on single layer may serve as a useful handle to further suppress the backgrounds in the quirk search. We consider models with five different sets of parameters ($m, ~\Lambda ,~p_T$). For each case, 10 K events are generated (assuming QCD production of quirk pair) to characterize the initial momentum distribution. Given an event in the CMS detector, different tracking layers can collect different numbers of hits. Among them, the largest one is recorded as $n^{\max }_{\text {hit}}$. Table 1 shows the fraction of events with a certain $n^{\max }_{\text {hit}}$ and in any event there is no minimum requirement on the number of hits in each layer. Note that we use the notation F$_{n^{\max }_{\text {hit}}<2}=\text {F}_{n^{\max }_{\text {hit}}=0}+\text {F}_{n^{\max }_{\text {hit}}=1}$. In event with $n^{\max }_{\text {hit}}<2$, the quirk pair can only induce at most one hit on each tracking layer. Similarly, F$_{n^{\max }_{\text {hit}}>4}=\text {F}_{n^{\max }_{\text {hit}}=5}+\text {F}_{n^{\max }_{\text {hit}}=6}+\cdots + \text {F}_{n^{\max }_{\text {hit}}=\infty }$. We can see that there is a large fraction of events (typically around $10\%$) that will induce more than two hits in at least one tracking layer. This intensive hit fraction (F$_{n^{\max }_{\text {hit}}>2}$) is considerable for $\Lambda \gtrsim \mathcal {O}(1)$ keV. As for a given $\Lambda$, the fraction F$_{n^{\max }_{\text {hit}}>2}$ increases with increasing ${m}/p_T$ because of the increased quirk oscillation amplitude.

To be specific, we use two parameters to characterize the shape of quirk pair trajectories. One is L given in Eq. 3.1 and the other is defined as

$$\begin{aligned} D= 2\frac{m}{\Lambda ^2} \frac{\rho \beta _o}{\sqrt{1-\beta _o^2}} . \end{aligned}$$

(6.1)

L corresponds to the width of the belt which the tracks are traveling inside, and D is the distance between two consecutive crossing points of two trajectories. If we only consider the hits on barrel tracking layer, we can project the quirk trajectories onto the transverse plane. The L and D will be projected into $L^\prime =L \cos \phi$ and $D^\prime =D\sin \theta$. In particular, we demonstrate in the Appendix D that the quirk pair can induce much more than 2 hits on a single tracking layer with radius R, if R lies between $D'$ and $L'$, or R is much larger than $L'$.

Fig. 5

Left: energy deposit of quirk pair in ECal and HCal within time of 25 ns. Right: the cumulative curve for the distributions of reconstructed MET. Note that the tails are dominated by the statistical fluctuation. The dashed lines have taken into account the quirk energy deposit

Fig. 6

Left: the distributions of standard deviations of the s angles for different quirk production processes. Right: the difference of quirk pair plane thickness between taking fixed ${\hat{s}}$ and true ${\hat{s}}$ respectively in solving quirk EoM

7 Effects of ionization energy loss in mono-jet search

As pointed out in Ref. [19] the non-helical trajectory of quirk will not be reconstructed in conventional searches at the LHC. So it will simply represent as missing transverse energy (MET) in event analyses. The mono-jet search at the LHC can be used to constrain the signal of quirk production with recoiling against a hard ISR jet.

In fact, as studied in Ref. [24], the quirk pair is not fully invisible. Since the quirk usually carries electric charge and travels with speed much smaller than the speed of light due to its heavy mass, it can deposit a certain amount of its energy inside the electromagnetic calorimeter (ECal) and hadronic calorimeter (HCal). In the following, we will consider the effects of quirk ionization energy loss inside ECal and HCal on the selection efficiency of mono-jet search.

After numerical simulation as introduced in Ref. [24], we can obtain the energy deposition in calorimeters of the CMS detector for different quirk production processes. In principle, slower quirk tends to deposit more energy inside calorimeters, since the ionization energy loss is proportional to $\sim v^{-2}$(for velocity $v \gtrsim 0.1$). However, in our simulation, we only consider the energy deposition within 25 ns, which is the time interval of bunch crossing at the LHC. The slowly moving quirk-pair can not pass through the calorimeters in time. It leads to the low energy deposit for the quirk production process with ${m}=500$ GeV, $p_T>100$ GeV, and $\Lambda =\Lambda _0$, as shown in the left panel of Fig. 5. While the process with ${m}=100$ GeV, $p_T>100$ GeV, and $\Lambda =\Lambda _0$ deposits largest energy because of its longest travel distance in the calorimeters. Note that larger $\Lambda$ leads to more accelerated quirk thus relatively smaller ionization energy loss.

In experimental analyses, the MET of an event is reconstructed by using all energy deposits in the calorimeters. It means the quirk energy deposit will also be taken into account. As a result, the MET is overestimated if one ignores the quirk energy deposit, as done in Ref. [19]. In the right panel of Fig. 5, we show how much the energy deposit of quirk will change the cut efficiency on MET, where the cumulative curves for the reconstructed MET of three different processes are shown. The dashed lines have taken into account the quirk energy deposits. In the lower subplot, the ratios between the MET without and with quirk energy deposits are given. We can see that cut efficiencies of MET is typically overestimated by a factor of 1.05, if the energy deposits of quirk are not included. This turns out to be a small effect in practical analyses.

8 Variation of the ${\hat{s}}$ direction

In solving quirk EoM, one usually assumes the straight-string approximation [9], i.e., the infra-color string is straight at a given time in the CoM frame. In order to ensure the simultaneity in the CoM frame, the space-time position in the lab frame for two quirks ($t_{1,2},\vec {r}_{1,2}$) should satisfy

$$\begin{aligned} t_1-t_2=\vec {\beta }\cdot (\vec {r}_1-\vec {r}_2). \end{aligned}$$

(8.1)

It requires that the time increasing step $\epsilon _{1,2}$ in numerical simulation satisfies

$$\begin{aligned}&\epsilon _{1}[1- \vec {v}_{1} \cdot \vec {\beta }-\frac{\vec {r}_{1}-\vec {r}_{2}}{E_{1}+E_{2}}\cdot (\vec {F}_{1}-\vec {v}_{1} \cdot \vec {F}_{1} \vec {\beta })]\nonumber \\&\quad = \epsilon _{2}[1-\vec {v}_{2} \cdot \vec {\beta }-\frac{\vec {r}_{2}-\vec {r}_{1}}{E_{1}+E_{2}}\cdot (\vec {F}_{2}-\vec {v}_{2} \cdot \vec {F}_{2} \vec {\beta })], \end{aligned}$$

(8.2)

where $\vec {F}_{i}=\vec {F}_{si}+\vec {F}_{exti}$ includes the infracolor force and external forces. Then, at any time $t'_{1,2}$, the ${\hat{s}}_{1}$ and ${\hat{s}}_{2}$ used in Eq. 2.2 for two quirks in the lab frame are the unit vectors of

$$\begin{aligned} \vec {r}_{s1}&=(\vec {r}'_1-\vec {r}'_2)-(t'_1-t'_2)\vec {v}_1, \end{aligned}$$

(8.3)

$$\begin{aligned} \vec {r}_{s2}&=(\vec {r}'_2-\vec {r}'_1)-(t'_2-t'_1)\vec {v}_2, \end{aligned}$$

(8.4)

respectively.

It is clear that the directions of ${\hat{s}}_{1,2}$ are varying at each time step of numerical solution. For each event, we can calculate the standard deviation ($\sigma ({\hat{s}})$) for the angle between ${\hat{s}}$ at all time steps and $\pm {\hat{e}}_x$ (the initial ${\hat{s}}$) to characterize the varying range of ${\hat{s}}$. The distributions of $\sigma ({\hat{s}})$ for all events of several processes are shown in the left panel of Fig. 6. The deviation of ${\hat{s}}$ increases as $\Lambda$ is increased, which can be around $\mathcal {O}(0.1)$ [rad] for $\Lambda \sim$ keV. Moreover, since the deviation is mainly induced by the Lorentz force, larger $p_T/{m}$ renders more significant deviation.

Finally, we give a brief discussion on how a few existing quirk searches will change if one simply assumes the ${\hat{s}}$ is fixed (which seems taken in Ref. [25]). In the right panel of Fig. 6, we plot the changes in the quirk pair plane thickness when the fixed ${\hat{s}}_{1,2}$ are used instead of Eqs. 8.3 and 8.4. The quirk pair plane thickness is changed dramatically (Note that a true quirk pair thickness is around $\mathcal {O}(100)$ $\mu m$ for our parameter choice). The influence is sensitive to the $p_T/m$ and mildly depends on the $\Lambda$. So choosing the correct ${\hat{s}}$ direction is critical in coplanar quirk search. On the other hand, we find that fixing ${\hat{s}}$ can only lead to at most 2–3% changes in the cut efficiencies of mono-jet search. Moreover, each quirk trajectory can still be reconstructed approximately as helix with $\chi ^2 <5$ when the confinement scale $\Lambda \lesssim \mathcal {O}(10)$ eV [19]. In this case, we find that fixing ${\hat{s}}$ can only change the $\chi ^2$ by less than 1% for ${m}=100$ GeV, $p_T>100$ GeV, and $\Lambda =10$ eV.

9 Conclusion

We solved the quirk equations of motion analytically in the limit of $|\vec {F}_{\text {ext}}|\ll \Lambda ^2$, such that the external force can be treated as a correction to the infracolor force. According to the analytical solutions, the quirk pair oscillation amplitude can be expressed in a precise way in both CoM frame and laboratory frame. Meanwhile, the number of periods for quirk traveling inside tracker can be calculated immediately by using kinematic variables and model parameters, without conducting time consuming numerical simulation.

The coplanar search proposed in Ref. [22] is one of the most efficient method for searching the quirk signal at colliders, when the confinement scale $\Lambda \in [\mathcal {O}(100)~\text {eV}, \mathcal {O}(10)~\text {keV}]$. We provided an approximate expression for the thickness of quirk pair plane in terms of kinematic variables and model parameters. Comparing with the precise numerical simulation results, we found the analytically expression can be valid up-to one order of magnitude if the number of quirk periods is between 1-100. This expression is especially useful to predict a model or parameter space which has large thickness of quirk pair plane, so that the coplanar search become less efficient. Also, we studied the electric charge dependence of the quirk pair plane thickness, and found that the plane thickness is linearly proportional to the quirk charges if the quirk-pair system is electric neutral, while the plane thickness will be increased by a factor of quirk period number if the quirk-pair system carries electric charge. The coplanar search becomes less efficient if the quirk crosses at least one of the tracking layers more than once. The probability of this multi-crossing increases with increasing $\Lambda$ and ${m}/p_T$, which is typically $\sim \mathcal {O}(0.1)$ for the parameters of interest in this work.

The effect of ionization energy loss inside the detector (including tracker, electromagnetic calorimeter, hadronic calorimeter and so on) is usually ignored in quirk signal analysis. We showed that this effect will lead to an overestimated MET cut efficiency by $\sim$ 5%. Moreover, the variation of the infracolor string direction ${\hat{s}}$ is typically small (much smaller than 0.1 [rad], depending on the kinematic variables and the $\Lambda$) for relatively small $\Lambda$. So it may be assumed that the direction is fixed in some analyses for simplification. We found that the correct direction is critical in coplanar quirk search. However, the mono-jet search and the heavy stable charged particle search are quite insensitive to the true ${\hat{s}}$.

Appendices

Appendix A: Derivation for $z^{+}(g)$ and $z^{-}(g)$

1.1 a. $z^{+}(g)$

The acceleration related to $z^{+}(g)$ is

$$\begin{aligned} \frac{d^2 z^{+}(g)}{d t_c^2}{\hat{e}}_z=\frac{q \vec {v}_{c1}(g)\times (B_c{\hat{e}}_{yc})}{m}{\sqrt{1-\vec {v}^2_{c1}(g)}}, \end{aligned}$$

(A.1)

leading to

$$\begin{aligned} z^{+}(g)&=\frac{1+(-1)^{[g]}}{2}z_o^{+}({\bar{g}})+\frac{1-(-1)^{[g]}}{2}z_o^{+}(1-{\bar{g}}), \end{aligned}$$

(A.2)

$$\begin{aligned} z_o^{+}(g)&=\frac{m q B_c}{\Lambda ^4}\left( 2\rho g- \tan ^{-1}[\rho ]+ \tan ^{-1}[\rho (1-2g)]\right. \nonumber \\&\quad \left. -\rho (1-2g)\ln \left[ \sqrt{\frac{1+\rho ^2}{1+\rho ^2(1-2 g)^2}}\right] \right) , \end{aligned}$$

(A.3)

with all relevant variables defined as in the main text.

1.2 b. $z^{-}(g)$

According to Eqs. A.2 and A.3

$$\begin{aligned} z^{+}(2n)=0,\quad \frac{d z^{+}(2n)}{dt_c}=0\quad n=0,1,2\ldots , \end{aligned}$$

(A.4)

which imply the total momentum of the quirk pair is invariant at $g=2n$. Moreover, from Eq. 2.12, we know the average torque on the system of two quirks is zero when g is increased from 2n to $2n+2$ such that the total angular momentum of the quirk pair is also invariant at $g=2n$. Then there must be

$$\begin{aligned} z_{1}(2n)=-z_{2}(2n)=z^{-}(2n)=0\quad n=0,1,2\ldots . \end{aligned}$$

(A.5)

This means that two quirks meet each other and have opposite velocity at $g=2n$ in the CoM frame. We hence conclude that kinematics of quirk system at $g+2$ can be obtained by rotating the system at g with an angle $-\eta$ around ${\hat{e}}_{yc}$, leading to

$$\begin{aligned} z^{-}(2n+\Delta )&=z^{-}(\Delta )\cos (n\eta )+(-1)^{[\Delta ]} r({\bar{\Delta }})\sin (n\eta )\nonumber \\&\approx z^{-}(\Delta )+(-1)^{[\Delta ]} r({\bar{\Delta }}) n\eta , \end{aligned}$$

(A.6)

where $\Delta \in [0,2]$ and $\sin \eta =\frac{\sqrt{1+\rho ^2}}{\rho }\frac{d z^{-}(g)}{dt_c}\Big |_{g=2}\approx \eta \ll 1$. From Eq. 2.11,

$$\begin{aligned} \vec {v}_{c1}(2n+\Delta )=\vec {v}_{c1}(2n-\Delta ), \end{aligned}$$

(A.7)

we have

$$\begin{aligned} z^{-}(2n+\Delta )-r(\Delta ) n\eta = z^{-}(2n-\Delta )+r(\Delta ) n\eta , \end{aligned}$$

(A.8)

and thus

$$\begin{aligned} z^{-}(\Delta )- z^{-}(2-\Delta )=r(\Delta )\eta , \end{aligned}$$

(A.9)

where $\Delta \in [0,1]$. We infer from Eqs. A.6 and A.9 that $z^{-}(g)$ can be formally written as

$$\begin{aligned} z^{-}(g)&=\frac{1+(-1)^{[g]}}{2}\left( z_o^{-}({\bar{g}})+\sin \left( \left[ \frac{g}{2}\right] \eta \right) r({\bar{g}})\right) \nonumber \\&\quad +\frac{1-(-1)^{[g]}}{2}\left( z_o^{-}(1-{\bar{g}})-\sin \left( \left[ \frac{g+1}{2}\right] \eta \right) r({\bar{g}})\right) . \end{aligned}$$

(A.10)

For $g\in [0,1]$, we have

$$\begin{aligned} \frac{d^2 z^{-}(g)}{d t_c^2}=\frac{d^2 z_o^{-}(g)}{d t_c^2}&\approx \frac{q E_z}{m}\sqrt{1-\vec {v}^2_{c1}(g)}, \end{aligned}$$

(A.11)

leading to

$$\begin{aligned} z^{-}_o(g)&\approx z^{-}_{o\prime }(g), \end{aligned}$$

(A.12)

$$\begin{aligned} z_{o\prime }^{-}(g)&=\frac{mq E_z}{\Lambda ^4}\left( \sqrt{1+\rho ^2}-\sqrt{1+\rho ^2(1-2g)^2}\right. \nonumber \\&\quad \left. -\rho (1-2g)\left( \sinh ^{-1}[\rho ]- \sinh ^{-1}[\rho (1-2g)]\right) \right) . \end{aligned}$$

(A.13)

Results from the numerical calculations show that the difference between $z^{-}_o(g)$ and $z^{-}_{o\prime }(g)$ can not be ignored when $g\in [0.5,1]$, even though they are still of the same order of magnitude. Using

$$\begin{aligned} z^{-}(0.5)&= w_{0.5}, \end{aligned}$$

(A.14)

$$\begin{aligned} z^{-}(1)&= w_{1} \end{aligned}$$

(A.15)

from the numerical calculations, we can approximately construct $z_o^{-}(g)$ as

$$\begin{aligned} z_o^{-}(g)&= {\left\{ \begin{array}{ll} \frac{w_{0.5}}{z_{o\prime }^-(0.5)}z_{o\prime }^-(g),&{}\quad g\in [0,0.5]\\ \left( \frac{w_{0.5}}{z_{o\prime }^-(0.5)} +(\frac{w_{1}}{z_{o\prime }^-(1)} -\frac{w_{0.5}}{z_{o\prime }^-(0.5)})(2g-1)\right) z_{o\prime }^-(g),&{}\quad g\in [0.5,1] \end{array}\right. }. \end{aligned}$$

(A.16)

1.3 c. Validation of the $z^{+}(g)$ and $z^{-}(g)$

The above unknown parameters $w_{0.5}$, $w_1$ and $\eta$ can be obtained from the numerical calculations. Moreover, the magnitude order relation

$$\begin{aligned} r(0.5)|\eta | \sim |w_1| \sim |z^{-}_{o\prime }(1)|, \end{aligned}$$

(A.17)

is also obtained by the numerical results. Following the definition,

$$\begin{aligned} z_1(g)&\equiv z^{+}(g)+z^{-}(g), \end{aligned}$$

(A.18)

$$\begin{aligned} z_2(g)&\equiv z^{+}(g)-z^{-}(g), \end{aligned}$$

(A.19)

in the CoM frame, in Fig. 7, we make a comparison plot for $z_1(g)$ and $z_2(g)$ calculated from Eqs. A.2 and A.10, with those obtained by numerical simulations. It can be seen that our analytical expressions for $z_1(g)$ and $z_2(g)$ match the numerical results precisely.

Fig. 7

$z_1$ and $z_2$ with respect to time variation obtained from analytical expressions as well as numerical simulation. For illustration, we choose $\vec {B}=(0,0,4)$ T in the lab frame, $\Lambda =500$ eV, quirks have charges of $\pm e$ and $m=100$ GeV, $\vec {P}_{o1}=(-279.6,250.7,4.8)$ GeV, and $\vec {P}_{o2}=(-36.5,102.5,303.6)$ GeV

Appendix B: Thickness of quirk pair plane

To simplify Eq. 5.4, we take the following approximations

$$\begin{aligned}&r(g)\approx {\left\{ \begin{array}{ll} 2r(0.5)g&{}\quad g\in [0,0.5]\\ 2r(0.5)(1-g)&{}\quad g\in [0.5,1]\\ \end{array}\right. },\end{aligned}$$

(B.1)

$$\begin{aligned}&z_o^{+}(g)\approx {\left\{ \begin{array}{ll} 2z_o^{+}(1)g^2&{}\quad g\in [0,0.5]\\ z_o^{+}(1)\left( 1-2(1-g)^2\right) &{}\quad g\in [0.5,1]\\ \end{array}\right. },\end{aligned}$$

(B.2)

$$\begin{aligned}&z_o^{-}(g)\approx z_o^{-}(1)g^2\quad g\in [0,1]. \end{aligned}$$

(B.3)

Then the conditions of $\frac{\partial d_c^2}{\partial j}= \frac{\partial d_c^2}{\partial c_z}=0$ give

$$\begin{aligned} j(l_0)&=\frac{1}{l_0}\sum _{i=1}^{l_0}[\frac{i}{2}], \end{aligned}$$

(B.4)

$$\begin{aligned} c_z&=\frac{z_o^{+}(1)}{2}. \end{aligned}$$

(B.5)

Using Eqs. B.4 and B.5, we obtain

$$\begin{aligned} d_c^2&=C_1+C_2, \end{aligned}$$

(B.6)

$$\begin{aligned} C_1&=\frac{\eta ^2r^2(0.5)}{3l_0}\sum _{i=1}^{l_0}\left( j(l_0)-[\frac{i}{2}]\right) ^2, \end{aligned}$$

(B.7)

$$\begin{aligned} C_2&=\frac{2z_o^{+ 2}(1)}{15}+\frac{z_{o}^{- 2}(1)}{5}. \end{aligned}$$

(B.8)

Fig. 8

The ratio between the quirk pair plane thickness obtained from numerical simulation $(D_N)$ and that calculated from Eq. 5.6 ($D_A$), with respect to the quirk trajectory period ($n_T$) inside the tracker

Table 2 Number of hits in different cylindrical barrels of CMS tracker induced by the quirk pair within 25 ns. We also provide the corresponding $L'$ and $D'$ values in the left part of the table

However, we know that the ideal quirk plane in the lab frame should also contain the interaction point at which the quirk pair was produced, which means we need to correct $C_2$ in $d_c^2$ (CoM frame thickness) to get $d^2$ (lab frame thickness). For illustration, considering a rectangle which has length L and width $2\sqrt{C_2}$ ($2\sqrt{C_2}\ll L$), the diagonal line of the rectangle is the line that crosses one of the vertexes while having smallest distance square to the points on the edges. The corresponding average of the distance square is $\frac{4}{3}C_2$ when $2\sqrt{C_2}\ll L$. So we multiply $C_2$ by $\frac{4}{3}$ in Eq. B.6 to get

$$\begin{aligned} d^2&=\frac{\eta ^2r^2(0.5)}{3l_0}\sum _{i=1}^{l_0}\left( j(l_0)-[\frac{i}{2}]\right) ^2\nonumber \\&\quad +\frac{4}{3}\left( \frac{2{z_o^{+ }}^2(1)}{15}+\frac{{z_{o}^{-}}^2(1)}{5}\right) . \end{aligned}$$

(B.9)

Appendix C: Validation of the quirk pair plane thickness

On one hand, we can use Eq. 5.6 to estimate the quirk pair plane thickness approximately with initial quirk kinematics. On the other hand, a more precise but time consuming way to obtain the thickness will be simulating the tracker configuration according to a specific detector and solving the quirk EoM numerically. Then, the quirk pair plane thickness square ($(D_N)^2$) corresponds to the smallest eigenvalue of the two-tensor [22]

$$\begin{aligned} T(\vec {h}_a)_{ij} = \frac{1}{N-1} \sum _{a=1}^{N} \vec {h}^a_i \vec {h}^a_j, \end{aligned}$$

(C.1)

where $\vec {h}_a$ is the position of ath hit in the tracker caused by the quirk pair. In Fig. 8, we plot the ratio between quirk pair plane thickness obtained from numerical simulation and that calculated from Eq. 5.6. The analytical result matches the numerical result within an order of magnitude when the quirk trajectory period ($n_T$) inside the tracker is around $\mathcal {O}(1-100)$. The case with $n_T<1$ can not be described by Eq. B.9, since the derivations in Sect. 5 use the oscillatory feature of the quirk motion which means that the number of periods needs to be larger than one. The quirk plane designed in Sect. 5 does not give the smallest plane thickness when $n_T<1$ and thus the plane thickness given by Eq. B.9 is larger than the numerical result. Strictly speaking, $E_z$ and $B_c$ respectively in Eqs. 5.8 and 5.9 change with the number of periods due to the rotation of the quirk-pair system described by Eq. A.10, but are thought to be invariant because of the small $\eta$. The effect of the system rotation can not be ignored when $n_T>\mathcal {O}(100)$, so the plane thickness given by Eq. B.9 is less than the numerical result.

Appendix D: More than two hits on single tracking layer

Here, we provide a benchmark study to count the number of hits on each tracking layer to illustrate our discussions in Sect. 6. The benchmark point is chosen as ${m}=100$ GeV, $k_1=12.22, k_2=12.86$, $\alpha =2.81, \theta =1.06$ and $\phi =0.72$. The confinement scale $\Lambda$ is varying from 400 eV to 3 keV. In Table 2, we show the number of hits produced by the quirk pair on each tracking layer (characterized by the radius R where we have adopted the CMS detector configuration). On the left part of the table, we also give the corresponding values for $L'$ and $D'$ that are defined in Sect. 6. This benchmark study clearly shows that the quirk pair can induce more than 2 hits on the tracking layer with radius R lying between $D'$ and $L'$.