text:
\section{Introduction} \label{sec:intro}
Current sheets are spatially localized plasma structures that play an essential role in various space plasma systems: solar flares \citep{Syrovatskii81,bookParker94,Fleishman&Pevtsov18},  solar wind turbulence \citep{Servidio11:npg,Borovsky10:solarwind,Vasko22:apjl}, boundaries of planetary magnetospheres – magnetopauses \citep{deKeyser05}, magnetotails of planets \citep{Jackman14,Achilleos18GMS:jupiter,Lui18} and comets \citep{Cravens&Gombosi04,Volwerk18A&A,Volwerk18GMS}. Strong currents flowing within current sheets are subject to various instabilities resulting in magnetic field line reconnection, which converts magnetic energy to plasma heating and particle acceleration \citep[e.g.,][]{book:Gonzalez&Parker,bookBirn&Priest07}.

The properties of magnetic reconnection and the dissipation rate of magnetic energy strongly depend on the current sheet configurations. The simplest magnetic field geometry is the one-dimensional (1D) current sheet, which is either a tangential discontinuity separating two plasmas with different properties or a rotational discontinuity having a finite magnetic field component, $B_n$, normal to the current sheet. The stress balance in 1D tangential discontinuities is established by the gradients of plasma and magnetic field pressures [see Figure \ref{fig:scheme}(a) and \citet{Allanson15,Neukirch20,Neukirch20:jpp}], whereas the stress balance in 1D rotational discontinuities requires a contribution from the plasma dynamic pressure in order to balance the magnetic field line tension force $\propto B_n$ [see Figure \ref{fig:scheme}(b) and \citet{Hudson70}].

In two-dimensional (2D) current sheets, $B_n\ne 0$ naturally appears \citep{bookSchindler06}; thus, the stress balance requires 2D plasma pressure gradients [see Figure \ref{fig:scheme}(c) and \citet{YL05}] or 2D dynamic pressure gradients [see Figure \ref{fig:scheme}(d) and \citet{Birn92,Nickeler&Wiegelmann10,Cicogna&Pegoraro15}].

Although there are many types of current sheet configurations, they can all be categorized by two plasma parameters: the Alfv\'en Mach number $M_A$ (i.e., the ratio of plasma flow speed to Alfv\'en speed) and the plasma beta $\beta$ (i.e., the ratio of plasma thermal pressure to magnetic field pressure). Current sheet configurations from different space environments are located within different domains in this $(\beta, M_A)$ space. Tangential discontinuities ($B_n=0$) do not require contributions from the plasma dynamic pressure. These discontinuities exist in systems with either large $\beta$ (if the plasma pressure balances the magnetic field pressure) or small $\beta$ [if the current sheet configuration is magnetically force free without cross-field currents, ${\bf J}\times{\bf B}=0$; see Figure \ref{fig:scheme}(a)]. Such current sheets are often observed in the solar wind \citep[see discussion in][]{Neugebauer06,Artemyev19:grl:solarwind}, where they propagate with plasma flows and have $M_A \approx 0$ in the current sheet reference frame \citep[see examples of kinetic models of such current sheet configurations in][]{deKeyser96,Harrison09:prl,Allanson16:jpp,Neukirch20}. Rotational discontinuities are also observed in the solar wind ($\beta\sim1$) and are characterized by plasma flow gradients with $M_A\sim 1$ in the current sheet reference frame [\citet[e.g.,][]{deKeyser97,Haaland12,Paschmann13:angeo,Artemyev19:grl:solarwind}; see Figure \ref{fig:scheme}(b)]. Figure \ref{fig:scheme}(e) shows the parameter regime of such current sheets in $(\beta, M_A)$ space, whereas Figure \ref{fig:space-observations}(a) shows an example of a low-$\beta$, force-free current sheet in the solar wind.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=1\textwidth]{schematic.png}
    \caption{Panels (a)-(d) show magnetic field lines and the main components of stress balance for different current sheet configurations. The coordinate system consists of the $l$ component along the main magnetic field direction, reversing sign at the current sheet neutral plane ($B_l=0$), the $m$ component along the main current density direction corresponding to variations of $B_l$, and the $n$ (normal) component along the spatial gradient of $B_l$ (i.e., $4\pi j_m/c=\partial B_l/\partial r_n$).  Panel (e) shows the typical parameter regimes of current sheets observed by THEMIS and ARTEMIS \citep{Angelopoulos08:ssr,Angelopoulos11:ARTEMIS} in the solar wind (blue), Earth's magnetosheath (shocked solar wind; black), near-Earth magnetotail (red), and lunar-distance magnetotail (green). The black dashed box defines the parameter regime where 2D force-free current sheets are expected. Our dataset includes $\sim300$ solar wind current sheets, $\sim100$ magnetosheath current sheets, $\sim 100$ current sheets in the near-Earth magnetotail, and $\sim 500$ current sheets in the lunar-distance magnetotail. All current sheets were identified from THEMIS and ARTEMIS fluxgate magnetometer measurements \citep{Auster08:THEMIS}. The selection procedure and data processing for solar wind current sheets are described in \citet{Artemyev19:jgr:solarwind}: ion temperatures are calculated using the OMNI dataset \citep{King&Papitashvili05,Artemyev18:jgr:report}, plasma flows in the current sheet reference frame are calculated as the change of the most variable flow component across the current sheet \citep[see details in][]{Artemyev19:jgr:solarwind}. The same criteria and data processing methods were applied to the magnetosheath (dayside) current sheets, but without using the OMNI data because THEMIS accurately measures magnetosheath ion temperatures in that region \citep{McFadden08:THEMIS}. The dataset of the near-Earth magnetotail current sheets (radial distances $\sim 10-30$ Earth radii) is described in \citet{Artemyev16:jgr:pressure}. The same criteria and data processing methods were applied to the lunar-distance current sheets. Details of the calculation of $\beta$ and $M_A$ are described in \citet{Artemyev19:jgr:solarwind,Artemyev17:jgr:mars}.  }
    \label{fig:scheme}
\end{figure}

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.9\textwidth]{observation.png}
    \caption{Four examples of force-free current sheets in (a) solar wind, (b) lunar-distance magnetotail, (c) Martian magnetotail, and (d) Jovian magnetotail. The top panels show the magnetic field in the local coordinate system \cite{Sonnerup68}. Grey curves show the magnetic field magnitudes; $B \approx \mathrm{constant}$ indicates the force-free current sheets. Middle panels show electron and ion $\beta$ profiles. Bottom panels show profiles of the field-aligned currents with an indication on the dominant ion species and estimations of the current sheet thickness in the ion inertial length. For the solar wind and Earth's lunar-distance magnetotail, we used ARTEMIS magnetic field \citep{Auster08:THEMIS} and plasma \citep{McFadden08:THEMIS} measurements \citep[see details of data processing procedure in][]{Artemyev19:jgr:solarwind}. For the Martian magnetotail, we used MAVEN magnetic field \citep{Connerney15:ssr} and plasma \citep{McFadden15:ssr,Halekas15:ssr} measurements \citep[see details of data processing procedure in][]{Artemyev17:jgr:mars}. For the Jovian magnetotail, we used Juno magnetic field \citep{Connerney17:ssr,Connerney17:science} and plasma \citep{McComas17:ssr,Kim20:juno:PS} measurements \citep[see details of data processing procedure in][]{Artemyev23:juno:arxiv}.  }
    \label{fig:space-observations}
\end{figure}

Current sheets in the shocked (magnetosheath) solar wind are characterized by lower $M_A$ [Figure \ref{fig:scheme}(e)] compared to that of the solar wind. The stress balance in these current sheets is significantly influenced by plasma pressure gradients, i.e., compressional current sheets with both parallel and transverse current density components \citep[e.g.,][]{Chaston&Travnicek20:magnetosheath_currents,Webster21}. Another $(\beta, M_A)$ domain, with similar $M_A$ and larger $\beta$, is found in Earth's distant magnetotail, where strong plasma flows may reach $M_A\geq 1$ and $\beta\in[10,100]$ \citep[e.g.,][]{Vasko15:jgr:cs,Artemyev17:jgr:mars}. These current sheets are mostly balanced by plasma pressure gradients. However, in the distant magnetotail and for $\beta \lesssim 10$ there have been previous observations of force-free current sheets with ${\bf J}\times{\bf B}=0$ \citep{Xu18:artemis_cs}. Figure \ref{fig:scheme}(e) shows the parameter domain of distant magnetotail current sheets, whereas Figure \ref{fig:space-observations}(b) shows an example of an almost force-free current sheet observed in the distant magnetotail. Both the magnetosheath and distant magnetotail current sheet configurations are characterized by a small but finite $B_n$, and thus the magnetic field line tension force $\propto B_n$ may be balanced by weak plasma flows or by plasma anisotropy. There is a series of models on such quasi-1D current sheets, with $B_n\ne 0$, that describe both large-$\beta$ \citep[see][]{Burkhart92TCS,Sitnov00,Sitnov06,Mingalev07,Zelenyi11PPR} and force-free \citep{Artemyev11:pop,Mingalev12,Vasko14:angeo_by} current sheets.

The distant magnetotail $(\beta, M_A)$ domain smoothly extends towards the lower $M_A$ with the decrease of the radial distance from the Earth. Near-Earth magnetotail current sheets are characterized by slightly higher $\beta>100$ and are quantified by the plasma pressure contribution to the stress balance \citep{Runov06,Artemyev11:jgr,Petrukovich15:ssr}. The important difference between solar wind/magnetosheath and magnetotail current sheets is their magnetic field line configuration; solar wind current sheets may be considered as 1D discontinuities, whereas the magnetotail current sheets are 2D [Figures \ref{fig:scheme}(a-d)]. Therefore, current sheet models with $B_n\ne0$ and low $M_A$ should include 2D plasma pressure gradients balancing the magnetic field line tension force \citep[see examples of such models in][]{SB02,Birn04,YL05,Sitnov&Merkin16}.

In the $(\beta, M_A)$ space, the large-$\beta$ domain can be described by kinetic models of current sheets with $B_n=0$ (1D models) and $B_n\ne 0$ (2D models with plasma pressure gradients). Conversely, the low-$\beta$ domain with dominant field-aligned currents (force-free current sheets) can be described only by 1D kinetic models with $B_n=0$. Large $M_A$, low-$\beta$ (1D rotational discontinuities with $B_n\ne 0$) and low $M_A$, low-$\beta$ (2D force-free current sheet) domains have been analyzed only by fluid models \citep[e.g.,][]{Cowley78,Hilmer87,Tassi08,Lukin18}. The class of such 2D force-free current sheets is not limited to observations in the solar wind and distant Earth's magnetotail, but also includes current sheets observed in the cold plasma of the Martian magnetotail \citep[e.g.,][]{DiBraccio15,Artemyev17:jgr:mars} and in the low-density Jovian magnetotail \citep[e.g.,][]{Behannon81,Artemyev14:pss}. Figures \ref{fig:space-observations}(c,d) show examples of force-free current sheets observed in the Martian and Jovian magnetotails. There are not enough known integrals of motion to describe distribution functions of charged particles in such current sheet configurations \citep[see discussion in][]{Lukin22:arxiv}. Consequently, the absence of kinetic models of force-free current sheets with $B_n \ne 0$ significantly hinders the analysis of their stability, the process responsible for the magnetic reconnection onset and particle acceleration.

In this study, we develop a 2D kinetic force-free current sheet model. In Section \ref{sec:mhd-equilibrium}, we describe a magnetohydrodynamic (MHD) model of 2D force-free current sheets. In Section \ref{sec:setup}, we initialize self-consistent, particle-in-cell (PIC) simulations using currents and magnetic fields from the MHD equilibrium. We load particle distributions using drifting Maxwellians, which initially do not necessarily satisfy the time-stationary Vlasov equation. In Section \ref{sec:results}, we obtain kinetic equilibria by evolving these particle distribution functions using the self-consistent PIC simulations and demonstrate the existence of kinetic equilibria for 2D force-free current sheets. We describe the equilibrium distribution functions and the difference between ion- and electron-dominated 2D current sheets. In Section \ref{sec:conclusion}, we summarize the results and discuss their applications.

\section{MHD equilibria of force-free current sheets} \label{sec:mhd-equilibrium}
We consider the general case of 2D force-free equilibrium in which all quantities vary with coordinates $x$ and $z$ only (i.e., $\partial / \partial y = 0$), and, as is customary in this definition, the pressure gradient force contribution to the force balance is negligible, i.e., $\nabla P \sim 0$. The divergence-free magnetic field has the form \citep[e.g.,][]{low1988spontaneous}
\begin{equation} \label{eq:b-vec}
    \mathbf{B} = \left(-\frac{\partial A}{\partial z}, B_y, \frac{\partial A}{\partial x}\right) ,
\end{equation}
where $A \mathbf{e}_y$ is the vector potential in the $y$ direction, and $B_y$ is the magnetic field in the $y$ direction. The current density is
\begin{equation} \label{eq:j-vec}
    \mathbf{J} = \frac{c}{4 \pi} \nabla \times \mathbf{B} = \frac{c}{4 \pi} \left(-\frac{\partial B_y}{\partial z}, -\nabla^2 A, \frac{\partial B_y}{\partial x}\right) ,
\end{equation}
where $c$ is the speed of light, and $\nabla^2 = \partial^2/\partial x^2 + \partial^2/\partial z^2$ is the Laplacian. The force-free condition, $\mathbf{J} \times \mathbf{B} = 0$, is equivalent to $\mathbf{J} = \alpha \mathbf{B}$, where $\alpha = \alpha(x, z)$ is a scalar function. Comparing Equations \eqref{eq:b-vec} and \eqref{eq:j-vec}, the force-free condition requires that $B_y$ is a function of $A$ only, and that $B_y$ satisfies
\begin{equation} \label{eq:vec-pot}
    \nabla^2 A + B_y(A) \frac{\mathrm{d} B_y}{\mathrm{d} A} = 0 ,
\end{equation}
and the coefficient $\alpha$ is $(c / 4\pi) \mathrm{d}B_y / \mathrm{d}A$. Once $B_y(A)$ is given, we can solve Equation \eqref{eq:vec-pot} for $A(x, z)$ with appropriate boundary conditions.

As shown in Figure \ref{fig:thermal-pressure-to-force-free}, for balancing force-free current sheets with low thermal and dynamic pressures (e.g., low $\beta$/low Mach number current sheets in the solar wind, magnetosheath, and lunar-distance magnetotail; see Figure \ref{fig:scheme}), the magnetic pressure $B_y^2$ plays an role analogous to the thermal pressure. Motivated by the functional form of pressure $p \propto \exp{\left[ 2 A / (\lambda B_0) \right]}$ describing thermal-pressure balanced current sheets \citep{LP82}, we adopt
\begin{equation}
	B_y(A) = B_0 \exp\left(\frac{A}{\lambda B_0}\right),
\end{equation}
for force-free current sheets, so that Equation \eqref{eq:vec-pot} describing magnetic field lines in the $x$-$z$ plane (i.e., contours of $A$) resembles that of the Lembege-Pellat current sheet \citep[e.g.,][]{Tassi08,Lukin18}:
\begin{equation}\label{eq:Lembege-Pellat-general}
	\nabla^2 A = -\frac{B_0}{\lambda} \exp\left(\frac{2 A}{\lambda B_0}\right) ,
\end{equation}
where $\lambda$ is the current sheet thickness, and $B_0$ is the magnetic field at large $\vert z \vert$.

\begin{figure}[tphb]
    \centering
    \begin{tikzpicture}[
		signal arrow/.style={decorate,decoration={show path construction,
				lineto code={
					\path let \p1 = ($(\tikzinputsegmentlast)-(\tikzinputsegmentfirst)$),
					\n1 = {int(mod(scalar(atan2(\y1,\x1))+360, 360))},
					\n2 = {veclen(\x1,\y1)}
					in
					(\tikzinputsegmentfirst) -- (\tikzinputsegmentlast)
					node[signal,midway,sloped,left color=red!30,right color=blue!30,draw,
					signal from=west, signal to=east,minimum width=\n2-\pgflinewidth,
					inner xsep=0pt,inner ysep=5pt,shading angle=\n1+90,
					anchor=center,#1]{};
		}}},
		signal arrow/.default=
		]
		\path (0,0) node[process1] (A) {$\begin{aligned}
				J_y B_z / c = \frac{\partial p}{\partial x} \\
				J_y B_x / c = \frac{\partial p}{\partial z}
			\end{aligned}$}
		(8,0) node[process2] (B) {$\begin{aligned}
				J_y B_z / c = J_z B_y / c = \frac{\partial}{\partial x} \frac{B_y^2}{8 \pi} \\
				J_y B_x / c = J_x B_y / c = \frac{\partial}{\partial z} \frac{B_y^2}{8 \pi}
			\end{aligned}$};
		\path[signal arrow] (A) node[anchor=south, yshift=32pt, color=red]{Near-Earth magnetotail} to node[anchor=south, yshift=5pt]{Mixed states}  (B) node[anchor=south, yshift=32pt, color=blue]{Solar wind; Lunar-distance magnetotail};
	\end{tikzpicture}
    \caption{Force balance of current sheets in different space environments. At one end of the spectrum, representing current sheets in the near-Earth magnetotail, magnetic tension force $J_y B_z /c$ and magnetic pressure force $J_y B_x / c \approx \partial (B_x^2 / 8 \pi) / \partial z$ are balanced by thermal pressure gradients $\partial p / \partial x$ and $\partial p / \partial z$, respectively. At the other end of the spectrum, representing current sheets in the solar wind and lunar-distance magnetotail, $J_y B_z /c$ and $J_y B_x / c$ are balanced by shears of $B_y^2$, i.e., $\partial (B_y^2 / 8 \pi) / \partial x$ and $\partial (B_y^2 / 8 \pi) / \partial z$, respectively. There may be a continuous spectrum of current sheets between these two ends \citep{yoon2023equilibrium}, in which $J_y B_z /c$ and $J_y B_x / c$ are balanced by a mixture of thermal pressure gradients and shears of $B_y^2$.}
    \label{fig:thermal-pressure-to-force-free}
\end{figure}

Because the scale length along current sheets is much larger than that across current sheets, we assume $A = A(\varepsilon x, z)$ to be weakly nonuniform in the $x$ direction ($\varepsilon$ being a small parameter). Equation \eqref{eq:Lembege-Pellat-general} can be approximated as
\begin{equation}\label{eq:Lembege-Pellat-magnetotail}
	\frac{\partial^2 A}{\partial z^2} = -\frac{B_0}{\lambda} \exp\left(\frac{2 A}{\lambda B_0}\right) ,
\end{equation}
which is accurate up to order $\varepsilon$. The boundary condition is
\begin{equation}\label{eq:bc-magnetotail}
	\frac{\partial A}{\partial z}\bigg\vert_{z = 0} = 0, \hspace{15pt} A\big\vert_{z = 0} = \varepsilon B_0 x ,
\end{equation}
which is equivalent to $B_x (z = 0) = 0$ and $B_z (z = 0) = \varepsilon B_0$. The solution of $A$ is
\begin{equation}
	A = - \lambda B_0 \ln\left[\cosh\left(\frac{z}{\lambda F(x)}\right) \cdot F(x)\right] ,
\end{equation}
where $F(x) = \exp\left(-\varepsilon x / \lambda\right)$. Thus the three components of the magnetic field are
\begin{equation} \label{eq:b-vec-LP}
	\begin{split}
		B_x &= B_0 \tanh\left(\frac{z}{\lambda F(x)}\right) F^{-1}(x) , \\
		B_y &= B_0 \left[\cosh\left(\frac{z}{\lambda F(x)}\right) \cdot F(x)\right]^{-1} , \\
		B_z &= \varepsilon B_0 \left[1 - \frac{z}{\lambda F(x)} \tanh\left(\frac{z}{\lambda F(x)}\right)\right] ,
	\end{split}
\end{equation}
where $B_z \neq 0$ describes rotational discontinuities. The current density is
\begin{equation} \label{eq:j-vec-LP}
	\mathbf{J} = \alpha \mathbf{B} = \frac{c \mathbf{B}}{4 \pi \lambda} \left[\cosh\left(\frac{z}{\lambda F(x)}\right) \cdot F(x)\right]^{-1} .
\end{equation}
The direction of magnetic field rotates from pointing in $+x$ at the $z>0$ half-space ($z / \lambda \gg 1$) to pointing in $+y$ around the neutral (equatorial) plane ($-1 \lesssim z / \lambda \lesssim 1$), and further to pointing in $-x$ at the $z<0$ half-space   ($z / \lambda \ll -1$). The magnitude of magnetic field, $(B_x^2 + B_y^2 + B_z^2)^{1/2} = B_0 [F^{-1}(x) + \mathcal{O}(\varepsilon)]$, has a weak dependence on $x$ and is roughly a constant at a given $x$ location. The magnitude of current density, $(J_x^2 + J_y^2 + J_z^2)^{1/2} \propto \cosh^{-1}[z / (\lambda F(x))]$, is concentrated in a layer $\vert z \vert < \lambda F(x)$.

\section{Computational setup} \label{sec:setup}
Searching for kinetic equilibria of 2D force-free current sheets, we initialize our simulated current sheets from the results of an MHD equilibrium. Importantly, this may not satisfy the time-stationary Vlasov equation because the initial particle distributions are simply drifting Maxwellians. We simulate their relaxation toward a certain kinetic equilibrium (if any) based on a massively parallel PIC code \citep{Pritchett01,Pritchett05:driven}. Our simulations have two dimensions $(x, z)$ in configuration space and three dimensions $(v_x, v_y, v_z)$ in velocity space. The results are presented in normalized units: magnetic fields are normalized to $B_0$, lengths to the ion inertial length $d_i = c / (4 \pi n_0 e^2 / m_i)^{1/2}$, time to the reciprocal of the ion gyrofrequency $\omega_{ci}^{-1} = m_i c / (e B_0)$, velocities to the Alfv\'en velocity $v_\mathrm{A} = B_0 / (4 \pi n_0 m_i)^{1/2}$, electric fields to $v_\mathrm{A} B_0 / c$, and energies to $m_i v_\mathrm{A}^2$. Here $n_0$ is the plasma density, $m_i$ is the ion mass, and $e$ is the elementary charge. The computational domain is $[-64 \leq x / d_i \leq 0] \times [-8 \leq z / d_i \leq 8]$ with a cell length $\Delta x = d_i / 32$. The time step is $\Delta t = 0.001\, \omega_{ci}^{-1}$. The ion-to-electron mass ratio is $m_i / m_e = 100$. The normalized speed of light is $c / v_\mathrm{A} = 20$, which gives the ratio of electron plasma frequency to electron gyrofrequency $\omega_{pe} / \omega_{ce} = 2$. The reference density $n_0$ is represented by $1929$ particles. The total number of particles is $1.6 \times 10^9$ including ions and electrons in each simulation.

Figure \ref{fig:init-config} shows the initial magnetic field configuration and current density, which are determined by Equations \eqref{eq:b-vec-LP} and \eqref{eq:j-vec-LP} using the input parameters $\varepsilon = 0.04$ and $\lambda = 2 d_i$. There are a few degrees of freedom in choosing the initial particle distribution functions: (1) While the current density is known a priori, the proportion of ion and electron currents is left unspecified; (2) Electron and ion temperatures and plasma $\beta$ vary across different space environments (e.g., Figure \ref{fig:space-observations}); (3) Specific forms of the velocity distribution functions are left undetermined. To this end, we choose initial electron and ion distribution functions as drifting Maxwellians
\begin{equation}
    f_{s} (\mathbf{v}, x, z) = \frac{n_0}{(2 \pi T_{s} / m_s)^{3/2}} \exp\left[-\frac{m_s \left(\mathbf{v} - \mathbf{v}_{d s}(x, z)\right)^2}{2 T_s}\right] ,
\end{equation}
where the subscript $s = i, e$ represents ions and electrons, $\mathbf{v}_{d s}(x, z)$ is the position-dependent drift velocity, and $T_s$ is the temperature. Note that $T_i$, $T_e$, and $n_0$ are constants throughout the domain so that the initial plasma pressure $n_0 (T_i + T_e)$ is a constant, as required by the force-free condition. The simulations allow either ions or electrons to be the main current carrier, while keeping the relative drift between them commensurate with the current density [Figures \ref{fig:init-config}(c), \ref{fig:init-config}(d), and \ref{fig:init-config}(e)]. Additionally, we vary particle temperatures to search for kinetic equilibrium in different plasma beta regimes, where the plasma beta is $\beta = 2(T_i + T_e) /m_i v_{\mathrm{A}}^2$. The detailed parameters for particle distribution functions in the six runs are listed in Table \ref{tab:init-dist-params}. These parameters cover the $\beta$ and $T_i/T_e$ ranges of force-free current sheets observed in the solar wind and planetary magnetotails (see Figure \ref{fig:space-observations}).

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.6\textwidth]{config.png}
    \caption{MHD equilibrium of a force-free current sheet. (a) Magnetic field in the $x$ direction $B_x$. (b) Magnetic field in the $y$ direction $B_y$. (c) Current density in the $x$ direction $J_x$. (d) Current density in the $y$ direction $J_y$. (e) Current density in the parallel direction $J_\parallel$.}
    \label{fig:init-config}
\end{figure}

\begin{table}[]
    \centering
    \begin{tabular}{c|c|c|c|c|c}
    \tableline
        Run NO. & $\mathbf{v}_{di}$      & $\mathbf{v}_{de}$       & $\beta_i = 2 T_i / m_i v_\mathrm{A}^2$      & $\beta_e = 2 T_e / m_i v_\mathrm{A}^2$      & $\beta = \beta_i + \beta_e$ \\
        \tableline
        1A      & $\mathbf{J} / (n_0 e)$ & $0$                     & $0.05$         & $0.05$         & $0.1$ \\

        1B      & $0$                    & $-\mathbf{J} / (n_0 e)$ & $0.05$         & $0.05$         & $0.1$\\
        \tableline
        2A      & $\mathbf{J} / (n_0 e)$ & $0$                     & $5/6$          & $1/6$          & $1$ \\

        2B      & $0$                    & $-\mathbf{J} / (n_0 e)$ & $5/6$          & $1/6$          & $1$ \\
        \tableline
        3A      & $\mathbf{J} / (n_0 e)$ & $0$                     & $25/3$         & $5/3$          & $10$ \\

        3B      & $0$                    & $-\mathbf{J} / (n_0 e)$ & $25/3$         & $5/3$          & $10$ \\
        \tableline
    \end{tabular}
    \caption{Parameters for particle distribution functions in the six PIC runs. The first digit in Run NO., `1', `2', `3', denotes three plasma betas, $\beta = 0.1, 1, 10$, respectively. The second digit in Run NO., `A', `B', denotes current carriers as ions and electrons, respectively. The current density $\mathbf{J}$ comes from Equation \eqref{eq:j-vec-LP}, and is visualized in Figures \ref{fig:init-config}(c) and \ref{fig:init-config}(d).}
    \label{tab:init-dist-params}
\end{table}

In our simulations, the electromagnetic fields are advanced in time by integrating Maxwell's equations using a leapfrog scheme. These fields are stored on the Yee grid. The relativistic equations of motion for ions and electrons are integrated in time using a leapfrog scheme with a standard Boris push for velocity update. The conservation of charge is ensured by applying a Poisson correction to the electric fields \citep{marder1987method,langdon1992enforcing}.

For particles crossing the $x$ boundaries, we take advantage of the symmetry between $z>0$ and $z<0$ [Figure \ref{fig:init-config}], so that particles exiting the system at a location $z$ with velocity $(v_x, v_y, v_z)$ are reinjected into the system at the conjugate location $-z$ with velocity $(-v_x, v_y, v_z)$. This is equivalent to an open boundary condition for particles because the injected particle distribution matches that at one cell interior to the boundary at all simulation times. At the $z$ boundaries, particles striking the boundary are reflected into the system with $v_z = - v_z$.

For fields at the $x$ boundaries, the guard values of the tangential magnetic fields are determined by
\begin{equation}
    \delta B^n_{y, g} = \delta B^n_{y, i1} , \hspace{15pt} \delta B^n_{z, g} = \delta B^{n-1}_{z, i1} (2 - \Delta x / c \Delta t) + \delta B^{n-2}_{z, i2} (\Delta x / c \Delta t - 1) ,
\end{equation}
where the superscript indicates the time level, and the subscripts $g, i1, i2$ indicate the guard point, first interior point and second interior point, respectively. This boundary condition for $\delta B_z$ ensures that the magnetic flux can freely cross the $x$ boundaries \citep{pritchett2005newton}. The guard values of the normal electric field $\delta E_x$ at the $x$ boundaries are determined by $\delta E_{x, g}^n = \delta E_{x, i1}^n$. Similarly, at the $z$ boundaries, the two components of the tangential magnetic fields in the guard point are determined by $\delta B_{x, g}^n = \delta B_{x, i1}^n$ and $\delta B_{y, g}^n = \delta B_{y, i1}^n$; The normal electric field $\delta E_z$ in the guard point are determined by $\delta E_{z, g}^n = \delta E_{z, i1}^n$. Unlike the driven simulations that adds magnetic flux to the simulation box, no external $E_y$ is applied at the $z$ boundaries.

\section{Results} \label{sec:results}
We track the evolution of the initially force-free current sheets in the six runs up to $t = 180\, \omega_{ci}^{-1}$, at which time the macroscopic states (e.g., electric and magnetic fields, currents, pressures) of the system are quasi-stationary. Below we examine if the configurations satisfy $\mathbf{J} \times \mathbf{B} = 0$ at the end of the simulations. We further investigate how the particle velocity distributions deviate from the initial Maxwellians and whether or not the system reaches a kinetic equilibrium.

\subsection{Initially ion-dominated current sheets}
When ions initially carry entirely the (field-aligned) currents (Runs 1A, 2A, and 3A), electrons are accelerated by transient parallel electric fields to form field-aligned currents [e.g., Figure \ref{fig:current-case1}(b)], while ions are slightly decelerated by such electric fields [e.g., comparing Figures \ref{fig:current-case1}(a) and \ref{fig:init-config}(e)]. These electron currents can be a significant fraction (e.g., $\gtrsim 1/3$) of ion currents in the off-equatorial region ($1 < \vert z \vert < 5$). Almost all currents remain field aligned [see $J_\perp \sim 0$ in Figure \ref{fig:current-case1}(c)], implying that a force-free configuration may indeed be obtained.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.6\textwidth]{prod_current_case1.png}
    \caption{Current density for the ion-dominated force-free current sheet with plasma beta $\beta = 1$ in Run 2A. The snapshot is taken at $t = 180\, \omega_{ci}^{-1}$ in the simulation. (a) Ion field-aligned currents. (b) Electron field-aligned currents. (c) Total perpendicular currents.}
    \label{fig:current-case1}
\end{figure}

During the relaxation of ion-dominated current sheets, electrostatic fields are generated and remain present in the late equilibrium states. These electrostatic fields are perpendicular to the magnetic field: At $\vert z \vert > 0$, they points away from the equatorial plane [Figure \ref{fig:efield}(a)]; Near the equator, they points in the $-x$ direction (i.e., tailward) [Figure \ref{fig:efield}(b)]. Here the $E_x$ component is much weaker than the $E_z$ component. The electrostatic fields arise due to the decoupling of the unmagnetized ions and magnetized electrons \citep{Schindler12}, which can be derived from the ordering among ion thermal gyroradius, current sheet thickness, and electron thermal gyroradius
\begin{equation}
    \rho_i : \lambda : \rho_e  = d_i \sqrt{\frac{\beta_i}{2}} : 2 d_i : d_i \sqrt{\frac{\beta_e}{2} \frac{m_e}{m_i}} = \sqrt{\frac{\beta_i}{8}} : 1 : \sqrt{\frac{\beta_e}{8} \frac{m_e}{m_i}} .
\end{equation}
Taking Run 2A for example, we have $\rho_i : \lambda : \rho_e = 0.32 : 1 : 0.01$. As the plasma beta is lowered, the electrostatic field decreases -- due to stronger magnetization of ions, which is shown in Figures \ref{fig:fbx} and \ref{fig:fbz} below.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.6\textwidth]{prod_efield_case1.png}
    \caption{Electrostatic field for the ion-dominated force-free current sheet with plasma beta $\beta = 1$ in Run 2A. The snapshot is taken at $t = 180\, \omega_{ci}^{-1}$ in the simulation. (a) Electric field in the $z$ direction $E_z$. (b) Electric field in the $x$ direction $E_x$. The magnitude of $E_x$ is much smaller than that of $E_z$.}
    \label{fig:efield}
\end{figure}

To show the detailed force balance of these current sheet configurations, we decompose the forces in parallel and perpendicular directions with respect to local magnetic fields. The force balance in the parallel direction is rapidly established because particles move freely along field lines. The unit vectors in the two perpendicular directions are defined as $\mathbf{e}_{\perp 1} = \hat{\mathbf{z}} \times \hat{\mathbf{b}}$ and $\mathbf{e}_{\perp 2} = \hat{\mathbf{b}} \times (\hat{\mathbf{z}} \times \hat{\mathbf{b}})$, where $\hat{\mathbf{b}}$ and $\hat{\mathbf{z}}$ are the unit vectors along the magnetic field and $z$ direction, respectively. It is worth noting that $\mathbf{e}_{\perp 1}$ is roughly along the $\pm y$ direction far above/below the equator, and along the $-x$ direction at the equator, while $\mathbf{e}_{\perp 2}$ is roughly along the $+z$ direction all over the domain.

At the equator, the electrostatic fields, $\mathbf{E}_{\perp 1}$ (i.e., along the $-x$ direction), cause a drift $c \mathbf{E}_{\perp 1} \times \mathbf{b} / \vert \mathbf{B} \vert$ of both ions and electrons [Figures \ref{fig:fbx}(b-c), \ref{fig:fbx}(h-i), and \ref{fig:fbx}(n-o)], which does not produce net currents [Figures \ref{fig:fbx}(a), \ref{fig:fbx}(g), and \ref{fig:fbx}(m)]. In addition, the pressure gradient $(\nabla \cdot \mathbf{P})_{\perp 1}$ is about zero in the two cases, with $\beta = 0.1$ [Run 1A; see Figure \ref{fig:fbx}(a)] and $\beta = 1$ [Run 2A; see Figure \ref{fig:fbx}(g)]. Thus, we obtain $(\mathbf{J} \times \mathbf{B})_{\perp 1} / c = (\nabla \cdot \mathbf{P})_{\perp 1} = 0$ in the $\mathbf{e}_{\perp 1}$ direction. In the case with $\beta = 10$, however, the current sheet gets compressed and rarefied periodically in the $x$ direction. Correspondingly, $(\nabla \cdot \mathbf{P})_{\perp 1}$ oscillates around zero and acts as a restoring force [Run 3A; see Figure \ref{fig:fbx}(m)]. Such oscillations are still localized around the equilibrium $(\mathbf{J} \times \mathbf{B})_{\perp 1} / c = (\nabla \cdot \mathbf{P})_{\perp 1} = 0$ when performing a spatial average over $x$.

Along the $\mathbf{e}_{\perp 2}$ direction (i.e., roughly the $z$ direction), both the pressure gradient $(\nabla \cdot \mathbf{P})_{\perp 2}$ and Lorentz force $(\mathbf{J} \times \mathbf{B})_{\perp 2} / c$ vanish [Figures \ref{fig:fbz}(a), \ref{fig:fbz}(g), \ref{fig:fbz}(m)]. The dominant component of the electrostatic fields, $\mathbf{E}_{\perp 2}$ (along the $\pm z$ directions), leads to a drift $c \mathbf{E}_{\perp 2} \times \mathbf{b} / \vert \mathbf{B} \vert$ (roughly in the $+y$ direction) of both ions and electrons [Figures \ref{fig:fbz}(b-c), \ref{fig:fbz}(h-i), and \ref{fig:fbz}(n-o)], which does not give net currents. All aforementioned perpendicular drift velocities are small, when compared to ion and electron parallel flow velocities that carry currents [e.g., Figure \ref{fig:current-case1}].

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.9\textwidth]{prod_fb_xstrip_perp1_avg.png}
    \caption{Force balance in the $\mathbf{e}_{\perp 1} = \hat{\mathbf{z}} \times \hat{\mathbf{b}}$ direction at the equator for different plasma betas. Note $\mathbf{e}_{\perp 2}$ is along the $-x$ direction at the equator. All the snapshots are taken at $t = 180\, \omega_{ci}^{-1}$ in the simulations. The volume forces are averaged in the range $-1 \leq z / d_i \leq 1$. The rows from top to bottom represent Runs 1A, 1B, 2A, 2B, 3A, and 3B, respectively. The columns from left to right represent the force balance from the point view of single fluid, ion fluid, and electron fluid, respectively.}
    \label{fig:fbx}
\end{figure}

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.9\textwidth]{prod_fb_zstrip_perp2_avg.png}
    \caption{Force balance in the $\mathbf{e}_{\perp 2} = \hat{\mathbf{b}} \times (\hat{\mathbf{z}} \times \hat{\mathbf{b}})$ direction for different plasma betas. Note $\mathbf{e}_{\perp 2}$ is approximately along the $+z$ direction. All the snapshots are taken at $t = 180\, \omega_{ci}^{-1}$ in the simulations. The volume forces are averaged in the range $-16 \leq x / d_i \leq -12$. The format is the same as Figure \ref{fig:fbx}.}
    \label{fig:fbz}
\end{figure}

\subsection{Initially electron-dominated current sheets}
When electrons initially carry entirely the (field-aligned) currents (Runs 1B, 2B, and 3B), there is no change in the macroscopic states of the system at the end of simulations [e.g., Figure \ref{fig:current-case2}]: The currents remain field aligned, and are carried by electrons only. Compared to the ion-dominated current sheets, the relaxation of electron-dominated current sheets does not alter the macroscopic states: neither the proportion of electron and ion currents is changed, nor do the electrostatic fields develop. In all three cases with different plasma betas, the system remains in a force-free configuration, with no plasma pressure gradients [Figures \ref{fig:fbx}(d-f), \ref{fig:fbx}(j-l), \ref{fig:fbx}(p-r), \ref{fig:fbz}(d-f), \ref{fig:fbz}(j-l), \ref{fig:fbz}(p-r)]. Despite the system being identical to the initial ones from the fluid point of view, the underlying velocity distribution functions evolve substantially from the initial Maxwellians to reach the Vlasov equilibrium, shown below.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.6\textwidth]{prod_current_case2.png}
    \caption{Current density for the electron-dominated force-free current sheet with plasma beta $\beta = 1$ in Run 2B. The snapshot is taken at $t = 180\, \omega_{ci}^{-1}$ in the simulation. (a) Ion field-aligned currents. (b) Electron field-aligned currents. (c) Total perpendicular currents.}
    \label{fig:current-case2}
\end{figure}

\subsection{Reaching kinetic equilibrium or not?}
Since the force-free equilibria are reestablished toward the end of simulations for both ion- and electron-dominated current sheets from the fluid point of view, a logical next step would be to further examine if such current sheets reach the kinetic equilibrium. Because all the currents are field aligned, we integrate the full distribution function $f_s(x, z, v_\parallel, v_\perp, \phi, t)$ of species $s$ to obtain the reduced distribution $g_s(x, v_\parallel, t)$ at the equator as a function of $x$, $v_\parallel$, and $t$:
\begin{equation}
    g_s(x, v_\parallel, t) = \int_{-d_i}^{d_i} \mathrm{d} z \int_{0}^{\infty} \mathrm{d} v_\perp \int_{0}^{2 \pi} \mathrm{d} \phi\, f_s(x, z, v_\parallel, v_\perp, \phi, t) ,
\end{equation}
where $v_\parallel$ and $v_\perp$ are the parallel and perpendicular velocities, respectively, and $\phi$ is the gyrophase.

Figure \ref{fig:vphase_space_case1} shows the difference between the reduced distributions at the end of simulations and the initial Maxwellians, $\Delta g_s = g_s(x, v_\parallel, t= 180\, \omega_{ci}^{-1}) - g_s(x, v_\parallel, t= 0)$, from Run 2A as an example. The ion phase space density is transported toward velocities both above and below the mean ion flow velocity [Figure \ref{fig:vphase_space_case1}(a)]. Consequently, the ion velocity distribution becomes wider and less peaked, compared to the initial Maxwellian [Figure \ref{fig:vphase_space_case1}(b)]. The mean ion flow velocity is slowed more on the earthward side than the tailward side [see the inset of Figure \ref{fig:vphase_space_case1}(a)]. Correspondingly, electrons are accelerated on average in the $v_\parallel < 0$ direction to develop an electron current to compensate the reduction of the ion current [see the inset of Figure \ref{fig:vphase_space_case1}(c)]. Such changes of mean flow velocities are not apparent in the phase space density plot because the thermal velocities are much larger than the mean flow velocities. The peak of the electron phase space density is transported toward larger velocities in both $v_\parallel > 0$ and $v_\parallel < 0$ [Figures \ref{fig:vphase_space_case1}(c) and \ref{fig:vphase_space_case1}(d)].

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.9\textwidth]{prod_vphase_space_case1.png}
    \caption{Phase space distributions for the ion-dominated force-free current sheet with plasma beta $\beta = 1$ in Run 2A. (a) The difference between the reduced ion distribution and the initial Maxwellian $\Delta g_i$. The inset plot shows the mean ion flow velocity as a function of $x$ at $t=0$ (gray line) and $t = 180\, \omega_{ci}^{-1}$ (black line). (b) A cut of the final ion distribution $g_i (x, v_\parallel, t = 180\, \omega_{ci}^{-1})$ (black line), $g_i (x, v_\parallel, t = 0)$ (gray line), and their difference $\Delta g_i$ multiplied by $5$, for better visibility (red line) between $-16 \leq x / d_i \leq -12$. (c) The different between the reduced electron distribution and the initial Maxwellian $\Delta g_e$. The inset plot shows the mean electron flow velocity as a function of $x$ at $t=0$ (gray line) and $t = 180\, \omega_{ci}^{-1}$ (black line). (d) A cut of the final electron distribution $g_e (x, v_\parallel, t = 180\, \omega_{ci}^{-1})$ (black line), $g_e (x, v_\parallel, t = 0)$ (gray line), and their difference $\Delta g_e$ multiplied by $5$, for better visibility (red line) between $-16 \leq x / d_i \leq -12$.}
    \label{fig:vphase_space_case1}
\end{figure}

The macroscopic states of electron-dominated current sheets do not differ from their initial configurations, as seen in Figure \ref{fig:current-case2}. The electron reduced distribution function $g_e(x, v_\parallel)$, however, shows systematic deviations from the initial Maxwellians, as displayed in Figure \ref{fig:vphase_space_case2}. Such deviations are similar to those in ion-dominated current sheets. The peak of the electron velocity distributions is reduced and redistributed toward larger velocities of $v_\parallel > 0$ and $v_\parallel < 0$ (but without changes in mean electron flow velocities) [Figures \ref{fig:vphase_space_case2}(c) and \ref{fig:vphase_space_case2}(d)]. The deviation of the ion velocity distribution from the initial Maxwellian is small [Figure \ref{fig:vphase_space_case2}(b)], but is required to reach the Vlasov equilibrium.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.9\textwidth]{prod_vphase_space_case2.png}
    \caption{Phase space distributions for the electron-dominated force-free current sheet with plasma beta $\beta = 1$ in Run 2B. The format is the same as Figure \ref{fig:vphase_space_case1}.}
    \label{fig:vphase_space_case2}
\end{figure}

To quantify the convergence of the system toward kinetic equilibrium, we define the metric
\begin{equation}
    G_s(t) = \int_{-\infty}^{\infty} \mathrm{d} v_\parallel \left\langle \left(g_s(x, v_\parallel, t) - g_s(x, v_\parallel, 0)\right)^2 \right\rangle_x ,
\end{equation}
where $\langle \cdot \rangle_x$ denotes the average over spatial coordinate $x$. This metric is chosen in such a way that the non-Maxwellian features of $g_s(x, v_\parallel, t)$ are properly accounted for. During the evolution of the system, it is expected that $G_s(t)$ experiences significant changes at the beginning of the simulation and eventually reaches a steady state, if a kinetic equilibrium is established. Figure \ref{fig:reaching_kequibrium} shows the absolute values of the time derivative $\vert \partial G_s /\partial t \vert$ for both ion- and electron-dominated current sheets with different plasma betas. All runs end up with $\vert \partial G_s /\partial t \vert \lesssim 0.01$, except for Run 1A (i.e., the ion-dominated current sheet with $\beta = 0.1$). For the same plasma beta, electron-dominated current sheets reach the kinetic equilibrium faster than ion-dominated current sheets. In the latter case, the system simply takes more time to redistribute the currents between unmagnetized ions and magnetized electrons. For both types (ion- and electron-dominated) of current sheets, those with higher plasma $\beta$ evolve toward kinetic equilibrium faster, because the current sheets with higher plasma $\beta$ can support larger transient electric fields to redistribute particles more rapidly in phase space toward the equilibrium.

At present, we cannot simulate further in time the relaxation of the low-$\beta$, ion-dominated current sheet in Run 1A toward the final kinetic equilibrium, because the current sheet becomes unstable at a later time (most likely due to boundary conditions). However, the progression of the relaxation for different beta values and the behavior of Run 1A up to $t \cdot \omega_{ci} \sim 10$ suggests that it also conforms to the explanatory model presented herein. The distribution functions at that time provide a sufficiently good representation of the Vlasov equilibrium to use in future modeling.

\begin{figure}[tphb]
    \centering
    \includegraphics[width=0.6\textwidth]{prod_kin_eqbrm_combine.png}
    \caption{Evolution of force-free current sheets toward kinetic equilibrium with different plasma betas. The $x$-axis has two scales, a linear scale in $0 \leq t \omega_{ci}^{-1} \leq 10$ and a log scale in $10 \leq t \omega_{ci}^{-1} \leq 180$. The $y$ axis denotes $\vert \partial G / \partial t \vert$. The rows from top to bottom represent Runs 1A, 1B, 2A, 2B, 3A, and 3B, respectively. The grey line $\vert \partial G / \partial t \vert = 0.01$ in each panel is drawn as a reference.}
    \label{fig:reaching_kequibrium}
\end{figure}

\section{Conclusion and discussion} \label{sec:conclusion}
In summary, we demonstrate that kinetic equilibria of 2D force-free current sheets exist for different proportions of ion and electron currents in various plasma betas. When initial currents are carried purely by ions, field-aligned electron currents are developed by transient parallel electric fields, while perpendicular electrostatic fields are generated due to unmagnetized ions and magnetized electrons. When initial currents are carried exclusively by electrons, the macroscopic state of the system remains unchanged from its initial state, described by the MHD model. In both scenarios, the electron and ion distribution functions at the late equilibrium states show systematic deviations from the initial drifting Maxwellians. These deviations occur in order to satisfy the time-stationary Vlasov equation.

The existence of a kinetic equilibrium for 2D force-free current sheets indicates that there is at least one hidden symmetry in the system, implying the existence of an additional integral of motion (i.e., an additional invariant). Our system has five dimensions $N_D = 5$ in the phase space (2D in the coordinate space and 3D in the velocity space). The two known invariants are the total energy $H$, and the $y$ component of the canonical momentum $P_y = m v_y + e A /c$. The degrees of freedom of the system are $N_f = N_D - N_I$, where $N_I$ is the number of invariants. The particle phase space densities at equilibrium are constructed as a function of such invariants of motion. To fully describe such kinetic equilibrium, which we now know it does exist, we must have $N_I > N_f$, or equivalently, $N_I > N_D / 2$. Thus, $N_I$ is at least $3$ (in our case) and there is at least one hidden symmetry. Starting from particle trajectory data in the equilibrium electromagnetic fields, future research on this topic could make use of machine learning models to help find the number of invariants or even discover the analytic formula of these invariants \citep[e.g.,][]{Liu&Tegmark21:ML,Liu&Tegmark22:Poinkare2}.

\section{Acknowledgments}.
This work was supported by NASA grants 80NSSC20K1788, 80NSSC22K0752, and NAS5-02099. We acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation \citep{cheyenne}.

redpajama_set_name: RedPajamaArXiv


text:
\section{Summary:} We present Rust-Bio, the first general purpose bio\-informatics library for the \change{innovative} \emph{Rust} programming language.
Rust-Bio leverages the unique combination of speed, memory safety and high-level syntax offered by Rust to provide a fast and safe set of bioinformatics algorithms and data structures with a focus on sequence analysis.

\section{Availability:}
Rust-Bio is available open-source under the MIT license at \href{https://rust-bio.github.io}{https://rust-bio.github.io}.

\section{Contact:} \href{koester@jimmy.harvard.edu}{koester@jimmy.harvard.edu}
\end{abstract}

\section{Introduction}

With ever increasing amounts of experimental data being generated, their computational analysis becomes \change{increasingly} challenging.
For novel or custom problems where carefully engineered high-performance standalone tools (like read mappers) are not yet available, general purpose bioinformatics libraries can help to minimize the coding effort.
Bioinformatics libraries are published for many popular programming languages, e.g., SeqAn for C++, Biopython, Bioperl and BioRuby \citep{Doring2008,Cock2009,Stajich2002,Goto2010}.
Choosing the programming language for a specific task usually entails a tradeoff between execution and development speed.
Low-level system programming languages like C or C++ provide optimal performance at the cost of increased complexity.
Higher level languages like Python or Perl \change{provide a more concise} syntax while leading to computational overhead introduced by online memory management (e.g. reference counting or garbage collection), type inference and not being compiled but interpreted during execution.
Often, the combination of a high-level language with some carefully engineered implementations of a bioinformatics library is a good choice to quickly solve a problem with reasonable performance.
\change{However, the amounts of data the bioinformatics community is facing in the coming years and the need to handle nature's resources carefully implies that using a high-performance, compiled language is still beneficial for certain problems.}

Recently, \emph{Rust}\footnote{\href{http://www.rust-lang.org}{http://www.rust-lang.org}} has gained attention as a new programming language combining speed with memory safety and high-level syntactical features.
Being compiled with LLVM \citep{Lattner2004}, Rust has many advantages of low-level, system programming languages, such as speed and a small memory footprint.
Supporting automatic type inference, it's code is often less verbose than C or C++ code.
With Rust, type inference happens at compile time, such that runtime overhead (appearing with scripting languages like Python) can be avoided.
The key feature of Rust is a concept of ownership and borrowing of variables, that enables the compiler to automatically decide about lifetime of objects during compile time, making an online memory management superfluous without requiring manual freeing of resources.
At the same time, this concept prevents common sources of errors with low-level languages like accessing invalid memory regions.
\change{Finally, the ownership concept enforces thread-safety, such that race conditions cannot occur.}
These features make Rust a promising solution to above tradeoff problem.

In this work, we present Rust-Bio, the first general purpose bioinformatics library for the Rust programming language.
Rust-Bio provides a high-level, fast and safe API for many state-of-the-art data structures and algorithms used in bioinformatics.

\section{Library}

Rust-Bio is built with the following principles in mind.
Where possible, iterators are returned. This allows to process streams of data with minimal memory footprint. On top, using the extensive set of iterator tools available in Rust, iterators can be e.g. filtered, modified, chained or combined in an easy way.
If a language data type appears suitable, we avoid to enclose data into a custom object. This mimimizes memory usage and increases flexibility when handling the data: e.g. biological sequences are represented as vectors or slices of bytes in ASCII encoding. This allows to use sequences with all algorithms and functions in e.g. the Rust standard library that work with byte vectors or slices.
Each implemented algorithm is automatically tested via continuous integration\footnote{\href{https://travis-ci.org}{https://travis-ci.org}}.
For each algorithm and data structure, we provide complexities in the documentation. Where more than one alternative is available, the documentation tries to highlight distinguishing use cases.
So far, Rust-Bio is focused on algorithms and data structures for biological sequences.
A central component of Rust-Bio are \emph{alphabets}, which, e.g., allow to check in linear time whether a given sequence is a word over the alphabet, transform symbols to their lexicographical ranks and perform bit-encoding to save memory or iterate over q-grams.
Rust-Bio can read and write common file formats like FASTA, FASTQ and BED.  For SAM/BAM, CRAM, and VCF/BCF support it is complemented by Rust-HTSlib.

Especially when considering sequencing data, many problems can be solved with a set of well established data structures like suffix arrays \citep{Manber1990}, the Burrows-Wheeler Transform \citep{Burrows1994}, rank/select data structures \citep{Jacobson1988} and $q$-gram indices.
In line with that, Rust-Bio implements induced sorting for suffix array construction \citep{Nong2009}, the FM-Index \citep{Ferragina2000} for  pattern matching on top of the Burrows-Wheeler Transform, a practical variant of a rank/select data structure \citep{Gonzalez2005} and a $q$-gram index for arbitrary alphabets and $q \leq 32$.
Further, Rust-Bio implements the FMD-Index \citep{Li2012}, that allows to find supermaximal exact matches in DNA sequences and their reverse complements in linear time.

Implementations for many classical pattern matching algorithms are provided, including the algorithm of Knuth, Morris and Pratt, Backward Nondeterministic DAWG Matching, Backward Oracle Matching, the algorithm of Horspool, and the Shift-And algorithm \citep{Knuth1977,GonzaloNavarro,Allauzen1999,Horspool1980,Wu1992}.
\change{In the supplement, we compare the speed of these algorithms against the C++ based Seqan, which is among the fastest bioinformatics libraries \citep{Doring2008}. The benchmarks exemplify that the speed of Rust-Bio is comparable to that of C++ based implementations.
For approximate pattern matching}, Ukkonen's dynamic programming based algorithm \citep{Ukkonen1985} and Myer's bit-parallel algorithm \citep{Myers1999} are provided.
Finally, Rust-Bio implements local, global and semi-global pairwise sequence alignment as variants of the Smith-Waterman and Needleman-Wunsch algorithms \citep{Needleman1970, Smith1981}.
An example for using the Rust-Bio API can be seen in Listing \ref{lst_example}.

\begin{lstlisting}[caption={Creating an FM-Index for a given text with an occurence table sampling rate of 3. Here, the alphabet is used to provide guarantees for being able to limit memory usage during FM-Index construction. Afterwards, we iterate over a FASTQ file, use the alphabet to validate read sequences and search for exact matches in the FM-Index.},float=t,label=lst_example,belowskip=-1.5\baselineskip]
let alphabet = alphabets::dna::iupac_alphabet();
let pos = suffix_array(text);
let bwt = bwt(text, &pos);
let fmindex = FMIndex::new(&bwt, 3, &alphabet);

let reader = fastq::Reader::from_file("reads.fastq");
for record in reader.records() {
    let seq = record.seq();
    if alphabet.is_word(seq) {
        let interval = fmindex.backward_search(seq.iter());
	let positions = interval.occ(&pos);
    }
}
\end{lstlisting}

\section{Conclusion}
Rust-Bio is a general purpose bioinformatics library.
Building on the innovative Rust programming language, Rust-Bio combines memory safety with speed, complemented by rigorous continuous integration tests.
So far, a wide set of algorithms and data structures for biological sequences is provided, ranging from index data structures to pattern matching and alignment, complemented by readers and writers for common file formats.

\section*{Acknowledgements}
We thank Christopher Schr\"oder and Peer Aramillo Irizar for code contributions as well as Sven Rahmann, Dominik Kopczynski, Tobias Marschall and Marcel Martin for their inspiring lecture notes "Algorithms on Sequences".

\bibliographystyle{natbib}

\section*{Rust-Bio---a fast and safe bioinformatics library: Supplement}
\subsection*{Benchmarks of pattern matching algorithms}
Since Rust-Bio is based on a compiled language, similar performance to C/C++ based libraries can be expected. Indeed, we find the pattern matching algorithms of Rust-Bio to perform in the range of the C++ library SeqAn\footnote{\url{http://www.seqan.de}}:

\begin{center}
\begin{tabular}{lrr}
\hline
Algorithm & Rust-Bio & Seqan \\
\hline
BNDM  &	77ms &	80ms \\
Horspool &	122ms &	125ms \\
BOM &	103ms &	107ms \\
Shift-And &	241ms &	545ms \\
\hline
\end{tabular}
\end{center}
We measured $10,000$ iterations of searching pattern
$$GCGCGTACACACCGCCCG$$
in the sequence of the human MT chromosome (assembly hg38). Initialization time of each algorithm for the given pattern was included in each iteration. Benchmarks were conducted with \emph{Cargo bench} for Rust-Bio and \emph{Python timeit} for SeqAn on an Intel Core i5-3427U CPU. Benchmarking SeqAn from Python timeit entails an overhead of around 1.46ms for calling a C++ binary.
This overhead was subtracted from above run times.

Note that this benchmark only compares the two libraries to exemplify that Rust-Bio has comparable speed to C++ libraries: all used algorithms have advantages depending on text and pattern structures and lengths. Details about when to use which pattern matching algorithm can be found in the documentation of Rust-Bio's pattern matching module.
\end{document}


text:
\section{INTRODUCTION}

Galactic black hole binaries are classified by their X-ray features,
such as the strength and temperature of the soft multi color
blackbody component, the hard X-ray emission, the X-ray luminosity
and   timing properties \citep{mitsuda84,vander95,remillard06}.
Black hole X-ray transient outbursts are known to pass through a
number of X-ray spectral states over the course of their outbursts;
these spectral states are also associated with the specific jet
behaviors \citep{vander95,vander00,vander05,mcclintock06}. Several
states (very high state, high/soft state, intermediate state,
low/hard state, quiescent state etc.) have been identified to
characterize black hole accretion
\citep{tanaka95,vander94,poutanen99,mcclintock06}. The soft state,
commonly described as $\sim 1$ keV thermal emission, was usually
observed when the source was   at higher luminosity, thereby
prompting the name``high/soft state''. The hard state, with a
typical photon index $\Gamma \sim 1.7$, is generally observed when
the source was in low luminosity, hence named "low/hard state". The
quiescent state, appearing at low luminosity levels, characterizes
the long periods of quiescence of transient systems. In the very
high state, black hole becomes exceedingly bright ($L_x$ $\geqslant$
0.2 $L_{Edd}$), and its X-ray spectrum displays a substantial
nonthermal radiation, which may constitute 40--90\% of the total
flux. In such a case, the photon index is typically $\Gamma \geq
2.4$, which is steeper than  that ($\Gamma \sim 1.7$) observed  in
the hard state \citep{mcclintock06}. More than 130 soft X-ray
transients (SXTs) have been found \citep{brocksopp04} up to now, and
their number is growing every year. In particular, there are a
number of supposed black hole SXTs which do not show a soft
component in their X-ray spectra. Instead they remain  in the
low/hard state (the X-ray emission was dominated by the power-law
spectra) throughout the outburst \citep{brocksopp04}. In addition,
there are very similar X-ray spectral properties in this subset, but
the light curve morphologies of these sources are very different,
including the X-ray, optical and radio light curves. As declaimed,
{Swift~J1753.5$-$0127 {}}, a newly detected black hole X-ray transient, likely
belongs to the subset of BHCs \citep{miller06}.

A new, bright, variable  $\gamma$-ray and X-ray transient source,
{Swift~J1753.5$-$0127 {}}, was detected firstly by the {\it Swift} BAT on  May 30,
2005 \citep{palmer05}, then its position,  $R.A.=17^h53^m28.3^s$,
$Dec.=-1^{o}27^{'}9.3^{''}$ (J2000.0) with an estimated 90\%
confidence uncertainty of $6^{''}$, was given by X-ray
Telescope(XRT) \citep{burrows05}. {Swift~J1753.5$-$0127 {}} was also clearly
detected in soft X-ray with the {\it Swift} X-ray Telescope (XRT)
 \citep{morris05}. Later on, its optical counterpart   was observed
by the MDM 2.4 m telescope on July 2, 2005 \citep{halpern05}, and
its radio counterpart was detected at 1.7 GHz with the Multi-Element
Radio-Linked Interferometer Network (MERLIN) on  July 3, 2005.
Further radio observations on July 4 and 5 indicate that this source
is variable, strengthening the likely association with the X-ray
transient \citep{fender05}. These radio observations suggest that
the source might in its low/hard state during the outburst.

To fully understand the BHC behavior  of X-ray transient, it is
instructive to outline the source properties for the whole  range of
X-ray flux from the maximum of the outburst to quiescence, so as for
this purpose {\sl RXTE} has proved to be an extremely unique tool
for studying the transients. Using the  {\sl RXTE} and {\sl
XMM-Newton} data, \cite{miller06} have analyzed the spectra of
{Swift~J1753.5$-$0127 {}} when the source approached the quiescence state. Their
results tend to support a cool (kT $\sim 0.2$ keV) accretion disk
around the central compact object, which may locate  at/near  the
innermost stable circular orbit (ISCO). Furthermore, the detection
of   a 0.6 Hz quasi-periodic oscillation (QPO) by {\sl RXTE} and the
measurement of a hard power-law spectrum by {\sl Swift/XRT}
\citep{morris05} suggest a hard state for this source
\citep{morgan05} during the outburst.

In this paper, we have analyzed the  TOO observations carried out by
Proportional Counter Array (PCA) and High-Energy X-ray Timing
Experiment (HEXTE) onboard {\sl RXTE}, and investigated in detail
the evolutions of the outburst, with the focus on the spectral and
the timing properties. The paper is organized to show the
observations and data reduction in Sec.2, the presentation of the
main results in Sec.3, and finally the discussion and conclusion are
performed in Sec.4 and Sec.5.

\section{OBSERVATION AND DATA REDUCTION}
Onboard the {\sl RXTE} satellite there are three detectors. The PCA
and HEXTE are the two co-aligned spectrometers with  narrow fields
of view. They provide a broad energy coverage from 2 to $\sim$ 250
keV. Another detector, the so-called ASM, is used to track the
long-term behavior of the source in the energy band 2-13 keV. The
PCA consists of 5 non-imaging Xe multiwire proportional counter
units (PCUs). Each unit covers the energy  2-60 keV and has the time
resolution $\sim 1\,\mu{\rm s}$ ($2^{-20}$ s). PCA has
$1^o\times1^o$ field of view and provides in sum the collecting area
of $\sim 6500\,{\rm cm}^2$. The HEXTE comprises two clusters, each
with 4 scintillation detectors sensitive to photons in the range
15-250 keV, collimated to view a common $1^o$ field, and these eight
detectors in the clusters provide a total collecting area of
1600$cm^2$ \citep{gruber96}.

The {\sl RXTE} TOO data of {Swift~J1753.5$-$0127 {}} are available for observations
from July 6, 2005 to March 11, 2006 (Fig.1). The software package
HEASOFT version 6.0 was adopted for reduction of the PCA data. We
get rid of the data for which the angle of the source above the
Earth limb is less than 10 degrees or for which the pointing offset
is greater than 0.02 degrees. We take the Standard-2  mode data ( 16
s time resolution and 129 energy channels between 2-60 keV) for
energy spectral analysis,  the Event mode E\_125us\_64M\_0\_1s
($\sim122$ $\mu$s time resolution and 64 energy channels between
2-60 keV) and the Standard-1 mode (0.125 s time resolution and 1
energy channel between 2-60 keV) for timing analysis. Only
proportional counter units (PCUs) 0 and 2 were always active during
our observation, but the propane layer on PCU-0 was lost when a
micro-meteorite created a small hole in the front window of PCU-0.
Therefore, we take the Standard-2 events from the PCU-2 and extract
the energy spectra from the upper xenon layer.

The extracted PCA spectrum covers the 2-60 keV
bandpass in 129 channels, each sampled in every 16 seconds. The
energy spectra at 18-180 keV  are derived from the HEXTE-A data.
Because the detector 2 of the HEXTE cluster B   lost the spectral
capability, only the  data from cluster A are
analyzed. We use the software package XSPEC, version 12.2.1, to
analyze the broadband (3-180 keV) spectra. The background spectra
are made by the FTOOLS {\em pcabackest}, with the latest `bright
source' background model enclosed. An instrument response files are
obtained with the tool {\em pcarsp}. An additional  0.5\% systematic
errors are added to the PCU-2 spectra because of the calibration
uncertainties.

\section{ANALYSIS and RESULTS}

\subsection{ASM and PCA Light Curves, Color-Color Diagram}
 ASM monitored the light curve of the entire outburst in the
1.3--12.1 keV energy band. As shown in  Fig.\ref{asmlc} is the ASM
light curve with each bin averaged over one day. The outburst
started on June 29, 2005, reached the maximum  of $\sim 210$ mCrab
(1 Crab $\sim$ 75 counts/s
\footnote{http://xte.mit.edu/ASM\_lc.html}) on July 06, 2006, and
finally decayed  exponentially to quiescence. The fast-rise and
exponential-decay (hereafter FRED) light curve is considered as the
`classic' signature of a soft X-ray transient \citep{chen97}. To
investigate the  evolution of the energy spectrum, two colors are
investigated (see, e.g. \cite{hansiger89}): the soft color is the
count ratio of 3.0--5.0 keV to  1.3--3.0 keV, and the hard color is
the count ratio of 5.0--12.1 keV to 3.0--5.0 keV. In Fig.\ref{asmlc}
the spectral evolution is  obvious only for soft color during the
formation of the outburst.

As shown in panel A of Fig.\ref{lc_color_qpo_time} is the available
PCA observations of {Swift~J1753.5$-$0127 {}} during the decay of outburst. The
residuals resulted from  fitting this light curve with two
exponential functions  indicate either the fluctuation of the flux
or several  additional small bursts. The light curve is subdivided
into three parts according to the evolution of the  colors and the
OPQ, as will be shown in the following.

The PCA colors are defined as the soft ( 6.1-9.4 keV/2.1-5.7 keV)
and the hard ( 13.5-16.9 keV/9.8-13.1 keV). Both colors increase
linearly when the burst  decay from the flux maximum to where the
linear relationship becomes rather weak, and the corresponding  time
period  is defined  as the  time zone 'part 1' (panel C/D of
Fig.\ref{lc_color_qpo_time}). Later on, both colors keep likely
constant during the rest of the decay. This trend is rather
straightforward  in the color-color diagram ( Fig.\ref{coco}).
Similar linear relationship between the  soft  and hard color was
reported in GRS 1915 + 105  by \citep{belloni00} in its low/hard
state. The second time zone (part 2) is ended when the QPO signal
vanishes. One sees from the Fig.\ref{lc_color_qpo_time} that in the
third time zone (part 3) the QPO signal appears occasionally around
the time when  the source is likely to have additional flare in flux
and evolution in hard color.

\subsection{Timing Analysis}

The power density spectra are produced with FTOOLS {\em powspec} for
the Event mode (E\_125us\_64M\_0\_1s) data in the energy band 2-60
keV. As shown in Fig.\ref{qpo_PDS} is one sample of the  PDS
obtained from the PCA observation ID 9143-01-01-00. The PDS consists
of several components  typical for the low-hard state in accreting
BHs \citep{cui99}: the   white noise at the low frequencies and the
red noise at the high frequencies. These noises are generally fitted
by a broken power law function  and the QPO in PDS is fit with a
Lorentzian function. The errors are estimated by  setting $\Delta
\chi^2$ = 1. For {Swift~J1753.5$-$0127 {}}, 39 of 62 PCA observations are found to
have low-frequency quasi-periodic oscillations (LFQPOs) at the
frequencies $\sim$0.1--0.9 Hz (see Table \ref{qpotable}). These
results are obtained without taking into account  the harmonics in
the spectral fitting. The LFQPOs are a common feature for the X-ray
BH binary while staying in their low/hard state, and elsewhere
detected already in V404 Cyg, GRO J1719-24, GRO J0422+32, GS 1354-64
and XTE J1118+480 \citep{brocksopp04}.

During the decay of the outburst, the QPO has a rough trend of
evolving to the lower frequencies (Fig.\ref{lc_color_qpo_time}).
Such a trend is obvious in the plot of the flux against  QPO
frequency (Fig.\ref{qpoflux}), where   a  linear relationship is
well established. Similar feature was reported in the other BHC like
GRS 1915 + 105 \citep{markwardt99,reig00}. It is worth noting that
 several QPO frequencies in Fig.\ref{qpoflux} strongly deviate
from the overall linear relationship. These QPOs are detected  when
the outburst either reached to the flux maximum or decayed to rather
low flux levels.

\subsection{Spectral Analysis}

The energy spectra of the black hole X-ray binaries are generally
described by the joint components of a blackbody and a power-law
shape\citep{remillard06}. In such a model the soft X-rays are
thought to be produced in the accretion disk while the hard X-rays
from the inverse Compton scattering of the soft photons off the
electrons in the hot corona residing in the inner disk.

To find the proper components in the energy spectrum we take the
{\sl RXTE} observation 91432-01-01-00 for spectral fitting. This
observation was carried out on July 7, 2005 (indicated in
Fig.\ref{lc_color_qpo_time} by the arrow symbol), and the data have
almost the best statistics. After standard screening, the remained
exposures for PCA and HEXTE are  6.4 ksec and 2.1 ksec,
respectively. The {\sl RXTE} PCU-2 (3.0--23.0 keV) and HEXTE-A
(18.0--180.0 keV) spectra are combined. A scaling factor is
introduced to account for the different normalization between PCA
and HEXTE. This parameter is set to unit for PCA and free for HEXTE.
The hydrogen column density is fixed at  $\rm N_H$ =$1.6\times
10^{21} cm^{-2}$. To account for the major component in the energy
spectrum, we have performed several trials. The fitting   results in
a reduced $\chi^2$ $\sim$ 10.9 with a simple power law model, a
reduced $\chi^2$ $\sim$ 6.14 with a compTT model and a reduced
$\chi^2$ $\sim$ 5.10 with a broken power-law model (see
Fig.\ref{phamodel}). Therefore we take for the first step the broken
power-law model. For the improvement of fittings, the excesses
showing up in the residuals of Fig.\ref{phamodel}  require the extra
models. By adding  a disk blackbody model ("diskbb"in XSPEC) the fit
is improved with the reduced $\chi^2$ $\sim$ 2.21. The excess
residing at the energies between 5 and 8 keV could in principle be
the contribution from the disk reflection.  By taking the smeared
edge model (\cite{ebisawa94}, ``smedge''in XSPEC) into account, the
reduced $\chi^2$ goes down to 0.97. However, the parameters of disk
blackbody and smeared edge cannot be constrained. A substitution of
the smeared edge model by a Gaussian with the center fixed at 6.4
keV results  in  the reduced $\chi^2$ of 1.16, but the errors for
disk blackbody parameters are not constrained as well by the data.
Accordingly,  the disk blackbody model is replaced by a blackbody
model (``bbodyrad'' in XSPEC). This gives a reduced $\chi^2$ of
1.19, and the parameters are well constrained as the follows: the
blackbody temperature kT = $0.66\pm 0.04$ keV, normalization $\rm
Norm_{bb}$ = $330^{+140}_{-100}$, equivalent blackbody radius
$21.3^{+4.2}_{-3.4}$ (d/8.5  kpc) km, break energy BreakE = $36.4^{+
1.8}_{- 1.7}$ keV, photon index $\Gamma_1$ = $1.64\pm 0.01$ and
$\Gamma_2$ = $2.04\pm 0.03$, equivalent width of Fe K$\alpha$ line
EW = 0.152 keV,  an unabsorbed flux of $4.7\times 10^{-9} \ \rm erg
\ cm^{2}\ s^{-1}$(2-10 keV) and a luminosity of $4.1\times 10^{37}\
\rm erg \ s^{-1}$ $(d/\rm8.5\ kpc)^2$ (or $L_X/L_{Edd} = 3.1\times
10^{-2} \ (d/8.5\ \rm kpc)^2 \ (M/10M_\odot )$). The comparison of
the trials with different models
 refers to Table \ref{tab:model}. By applying the
bbodyrad+bknpower+gauss model (hereafter BBG model), the whole data
are analyzed and the results are presented in Table \ref{tab:bknpo}.

\section{DISCUSSION}
Comparing the light curve morphology of the 2005-2006 outburst of
{Swift~J1753.5$-$0127 {}} with that of the 2000 first outburst of XTE J1118+480
\citep{wood01}, we find that  both exhibit  the  very similar
FRED-shaped form, which is often considered as a characteristic for
the `soft X-ray transients'(SXTs) \citep{lasota01}, then the account
for the  SXT outbursts is usually  referred to  the thermal-viscous
disk instability model (DIM) \citep{lasota01,brocksopp04} or to a
diffusion model \citep{wood01}. The FRED-shaped light curve of
{Swift~J1753.5$-$0127 {}} can be well described using above two models.

The high and low spectral states firstly identified in Cyg X-1
\citep{tananbaum} have   been observed in a number of X-ray
binaries, where the ``high'' and ``low'' terminology is originally
chosen based on the 2-10 keV X-ray flux. It is later found that the
low state corresponds to a hard spectrum and the high state to a
soft spectrum. The spectra and fluxes (see Table \ref{tab:bknpo})
indicate that {Swift~J1753.5$-$0127 {}} is in its low/hard state. Strikingly, the
observed state of {Swift~J1753.5$-$0127 {}} is similar to those of a few BHCs (for
instance, XTE J1118+480 and GRO J0422+32), remaining in the low/hard
state during outburst.  For explaining the behavior of the SXTs
which can remain in the low/hard state during outburst, a model was
proposed by \cite{meyer04}. In their model, a transition to the
high/soft state depends on the peak mass flow rate in outburst, and
the critical mass flow rate range is 0.02 to 0.05 $L_X/L_{\rm Edd}$.
The thermal component comes from the blackbody emission of a
truncated thin disk and the non-thermal one from the Comptonization
of the soft X-ray. The highest flux of Swift~J1753.5$-$0127 {} is $L_X/L_{Edd} =
3.1\times 10^{-2} \ (d/8.5\ \rm kpc)^2 \ (M/10M_\odot )$, which lies
within the declaimed critical mass flow rate. We also find that the
thermal component  presented in the spectra from the outburst to
quiescence  might come from the disk. But the parameters of the
thermal component cannot be constrained by the PCA  data.
\cite{miller06} analyzed the {\sl XMM-Newton} and {\sl RXTE}
observations when the source is near its quiescence. A thermal
component which comes from a thin disk was also found, however, a
smaller disk radius is not consistent with Meyer-Hofmeister's model.
\cite{miller06} suggest that the {Swift~J1753.5$-$0127 {}} may have a low mass black
hole than the typical black hole candidates. The low mass black hole
has a less evaporation efficiency to allow the disk to grow near the
ISCO, which  accounts  for the thermal component\citep{meyer04}.

The X-ray PDS of {Swift~J1753.5$-$0127 {}} displays LFQPOs in the range 0.1 to 0.9
Hz (see panel E of Fig.\ref{lc_color_qpo_time}), which is
extraordinarily important for understanding  the behaviors of
accreting black holes \citep{vander05,belloni02}, and usually they
are associated with the nonthermal or low/hard states and
transitions, implying the QPO emitting regimes to be far away from
ISCO, perhaps $\sim$100 Schwarzschild radius for the $\sim$10 solar
mass central object if the Keplerian frequency (maybe not necessary)
is taken into account. Two models are often mentioned to describe
the LFQPO \citep{brocksopp04}. \cite{titarchuk00} suggest that the
LFQPO of X-ray transients may be caused by the vertical oscillation
of the inner edge of accretion disk, and they predict a proportional
correlation between QPO frequency and the color parameter.
Furthermore, they also suggest  that the spectrum gets harder and
the QPO frequency increases when the disk mass accretion rate drops.
However, our results would expect the inner edge of disk to move
outwards with decreasing the QPO frequency, while moving back to its
``quiescent state'', which has been hinted in panel E of
Fig.\ref{lc_color_qpo_time}. With the exponential-decay light curve,
the frequencies of LFQPOs are decreasing, as shown in
Fig.\ref{lc_color_qpo_time}. So, generally, there exists an
approximate linear relation between the QPO frequency and  flux, and
the similar phenomena have been noted in GRS 1915 + 105
\citep{markwardt99,reig00}. In Fig.\ref{qpo_color}, a plot of the
soft/hard color versus the QPO frequency, we can find that both the
soft and hard colors drop with the QPO frequency, and a similar
anti-correlation between QPO frequency and hardness ratios has been
found in  BHC GRS 1915+105 \citep{reig00}. Therefore, LFQPOs of
{Swift~J1753.5$-$0127 {}} might be unsatisfactorily explained by the model of
\cite{titarchuk00}.

We would like to note that in Fig.\ref{qpoflux} there are several
data points staying away from the overall trend of a linear
relationship between the flux and QPO. This phenomenon in
  could not be well explained  by the models of
\cite{markwardt99}, \cite{reig00}, and \cite{titarchuk00}. As shown
in Fig.\ref{lc_color_qpo_time}, an additional small flare occurred
in the third time region is correlated to the three QPO signals
standing out of the lowest flux  in Fig.\ref{qpoflux}. During this
small flare, the spectrum becomes soft (panels D and E in
Fig.\ref{lc_color_qpo_time}). This might suggest  that the QPO
origins from  the perturbation at somewhere of accretion disk, since
otherwise the high QPO frequency would require larger luminosity
than observed according to the main trend of a linear relationship.
The QPO frequency (Fig.\ref{qpoflux}) detected at the peak flux of
the outburst  drops  from 0.906 Hz to 0.654 Hz within one day, while
the  flux decreases by only 9 percent and is largely deviated from
the overall linear trend. This is in general  not consistent with
the arguments of Meyer-Hofmeister (2004) for explaining the
relationship between the QPO frequency and  flux. In the paper of
Meyer-Hofmeister (2004), the transition in spectral state  is the
result of the balance between the evaporation of the inner edge of
  disk and the inward accretion flow. The closer  to the black
hole for the inner disk, the larger the flux and the higher value
the  QPO frequency. In the  region very close to the black hole, the
coronal evaporation becomes dominant, and instead of being
completely truncated the inner edge of the disk might constitute a
soft extension toward the corona. As a result, the source might
remain hard state due to the evaporation of the accretion flow, and
presents the larger frequency of QPO in the power spectrum. The
deviation of this point from the overall trend might provide the
evidence from the observational point of view that the inner edge of
the disk might have a soft extension and evaporate rapidly to
prevent from establishing effectively the high/soft state.

An alternative explanation might be the radial inflow
effect\citep{vander01}. By connecting the deviated QPO points in
Fig.\ref{qpoflux} one gets a line which is roughly parallel to the
track of the main linear trend. Similar parallel tracks are observed
in some kHz-QPO sources in plots of QPO frequency vs. X-ray
luminosity. According to \cite{vander01}, in the so-called radial
inflow scenario, the QPO frequency is determined by some balance
between the accretion rate through disk and radiative stresses, and
it can be expressed as: $\nu(t)\propto\left({\dot M_d(t)\over
L_x(t)} \right)^\beta,(\beta>0)$, where $\dot M_d(t)$ is the
instantaneous accretion rate and $L_x(t)$ is the luminosity.
$L_x(t)$ has both an immediate response to $\dot M_d(t)$ and the
filtered one $\av{\dot M_d}(t)$ averaged over a special time period:
$L_x(t)\propto\dot M_d(t)+\alpha\av{\dot M_d}(t),(\alpha>0)$. The
adjustments on the values of $\dot M_d(t)$ and $\av{\dot M_d}(t)$
could lead to the parallel tracks  in the plot of QPO vs. flux, as
might be the case as well for {Swift~J1753.5$-$0127 {}}.

In panel a \& b of Fig.\ref{qpo1_vs_parameter}, we plot the
relations between the QPO frequency and the flux of blackbody and
broken power-law respectively, and find that both of them have a
positive linear relation, implying  the LFQPO  correlation  with the
spectral parameters for both the thermal and nonthermal components.
Similarly,  \cite{markwardt99} and \cite{trudolyubov99} have shown
that the frequency of QPO between 0.5 and 10 Hz of GRS 1915 + 105 is
positively correlated with both the thermal and nonthermal flux.
However, in imagination, the oscillation could  originate in the
accretion disk if the power-law spectrum comes from the inverse
Compton scattering of disk photons,  and so the thermal and
nonthermal components would have some close correlation. Thus this
QPO appears to link both the accretion disk and density of Compton
scattering electrons.  Using the normalization of blackbody model,
we calculate the equivalent radius of blackbody area, and plot its
relation with the QPO frequency in Fig.\ref{qporadius}, showing  a
linear relation. When the flux drops, the equivalent radius gets
smaller, i.e. the radiation area of blackbody in the disk area
decreases, and the QPO frequency decreases because the radius of
inner disk increases.

In the work by \cite{titarchuk04}, they declaim that the observed
low QPO frequency---spectral-index correlation is a natural
consequence of an adjustment of the Keplerian disk flow to the
innermost sub-Keplerian boundary conditions near the central object.
This ultimately leads to the formation of the sub-Keplerian
transition layer between the adjustment radius and the innermost
boundary (the horizon for BH). In panel c and d of
Fig.\ref{qpo1_vs_parameter}, we present the correlation of photon
index $\Gamma$ versus QPO frequency. For QPO frequency values above
0.4 Hz, the QPO frequency and power-law index appear to be
correlated each other positively. However, for the frequency $\leq$
0.4 Hz, the photon power-law index is determined less accurate, and
remains almost constant. The similar properties of the index--QPO
correlation in GRS 1915+105 and XTE J1550-546 have been found
\citep{vignarca03,titarchuk04,shaposhnikov06}. Then in GRS 1915+105,
\cite{titarchuk04} stress that the hard state is related to an
extended Compton cloud (cavity) characterized by a photon index of
around 1.7 and low QPO frequencies below 1 Hz.
 This is the regime where thermal Comptonization dominates the
upscattering of soft disk photons and the spectral shape (index) is
almost independent of mass accretion rate. Thus, panel c and d of
Fig.\ref{qpo1_vs_parameter} suggest that the   source was in the
low/hard state. In addition, \cite{shaposhnikov06} argues that the
index saturation level is determined by the temperature of the
converging flow where the soft (disk) photons are upscattered by
electrons to the energy of falling electrons. In principle, one can
evaluate the mass of the central black hole using the index-QPO
relation because QPO frequencies are inversely proportional to the
BH mass \citep{titarchuk04,shaposhnikov06}.

\section{ CONCLUSION}

Almost all the lights under the current analysis point to the fact
that {Swift~J1753.5$-$0127 {}} kept staying in the low/hard state during the
outburst, along with the evolutions of  the spectral and the timing
features similar to those of  other BHCs like XTE J1118+480 and GRO
J0422+32. By far all these sources belong to this small subset of
BHC were observed to have almost only the low/hard state while in
outburst. The corresponding scenario could be that, due to the low
peak luminosity, the coronal evaporation truncates the thin disk
even in outburst and the spectra have to stay hard \citep{meyer04}.
During the outburst, {Swift~J1753.5$-$0127 {}} was detected to have its QPO
frequency varying  linearly with the flux for most of the time. An
interesting deviation from this trend as observed at the peak flux
might provide the first evidence of 'seeing' a rapid evaporation for
the inner edge of the disk, which prevents from effectively
establishing the high/soft state during the early stage of the
flare.  The relationship between the QPO frequency and spectral
color is anti-correlation. Although on a theoretical ground the
mechanism of the LFQPO of {Swift~J1753.5$-$0127 {}} is still open, this seems not
likely to be satisfied  by the global disk oscillation model
\citep{titarchuk00}.

\acknowledgments This work has been partially supported by the
National Nature Science Foundation of China , Hebei Province Nature
Science Fundamental Project No.A2006000128 and made use of data
obtained through the HEASARC online service, provided by the
NASA/GSFC.