text:
# Branching Fractions of 𝜏 Leptons to Three Charged Hadrons

CLEO Collaboration

## Abstract

From electron-positron collision data collected with the CLEO detector operating at CESR near $`\sqrt{s}`$=10.6 GeV, improved measurements of the branching fractions for $`\tau `$ decays into three explicitly identified hadrons and a neutrino are presented as $`(\tau ^{}\pi ^{}\pi ^+\pi ^{}\nu _\tau )=(9.13\pm 0.05\pm 0.46)\%,`$ $`(\tau ^{}K^{}\pi ^+\pi ^{}\nu _\tau )=(3.84\pm 0.14\pm 0.38)\times 10^3,`$ $`(\tau ^{}K^{}K^+\pi ^{}\nu _\tau )=(1.55\pm 0.06\pm 0.09)\times 10^3,`$ and $`(\tau ^{}K^{}K^+K^{}\nu _\tau )<3.7\times 10^5`$ at 90% C.L., where the uncertainties are statistical and systematic, respectively.

preprint: CLNS 03-1818preprint: CLEO 03-04

For hadronic $`\tau `$ decays, final states with kaons provide a powerful probe of the strange sector of the weak charged current. In this Letter, we present the improved measurements of branching fractions for $`\tau ^{}\pi ^{}\pi ^+\pi ^{}\nu _\tau `$, $`K^{}\pi ^+\pi ^{}\nu _\tau `$, and $`K^{}K^+\pi ^{}\nu _\tau `$ decays, and an improved upper limit for the phase space suppressed and Cabibbo-suppressed decay $`\tau ^{}K^{}K^+K^{}\nu _\tau `$. (Charge conjugate decays are implied). The decay mode $`\tau ^{}K^{}\pi ^+\pi ^{}\nu _\tau `$ has a significant contribution to the overall strange spectral function which can provide a direct determination of the strange quark mass and, potentially, precision extraction of the CKM element $`V_{us}`$ msvus . The decay $`\tau ^{}K^{}K^+\pi ^{}\nu _\tau `$ can proceed via both the vector and axial-vector currents prd47  and, therefore, has a sensitivity to the Wess-Zumino anomaly WZ .

The data used for this analysis were collected with the CLEO III detector CLEOIII  located at the symmetric $`e^+e^{}`$ Cornell Electron Storage Ring (CESR). The data sample consists of 3.26 fb-1 taken on or near the $`\mathrm{{\rm Y}}(4S)`$, corresponding to 2.97$`\times 10^6`$ $`\tau ^+\tau ^{}`$ pairs. The CLEO III detector configuration features a new four-layer silicon strip vertex detector, a new wire drift chamber, and, most importantly for this analysis, a ring imaging Cherenkov (RICH) particle identification system. Particle trajectories, momenta, and charges are measured with the tracking system, which uses hits from both the silicon detector and drift chamber. The specific ionization, $`dE/dx`$, measured in the drift chamber’s 16 inner axial and 31 outer stereo layers is sensitive to the traversing particle’s mass. The RICH detector tom  surrounds the drift chamber, and uses 1 cm-thick LiF radiators to generate Cherenkov photons from the incident charged particles. These photons then propagate through a 20 cm-thick expansion volume of gaseous nitrogen at atmospheric pressure, and are detected by multi-wire proportional chambers filled with a mixture of triethyleamine (TEA) and CH4 gases, and localized by readout from 8 mm-square cathode pads located on a cylinder of $``$1 m radius. RICH particle identification is available over the central region of the detector, $`|\mathrm{cos}\theta |<0.83`$, where $`\theta `$ is the polar angle with respect to the direction of the incident positron beam. A detailed description of the RICH performance can be found in Ref. RICH . The electromagnetic calorimeter surrounds the RICH detector and measures the energy, position, and lateral shape of showers induced by charged and neutral particles. It contains 7784 16-radiation-length-long CsI(Tl) crystals, arranged in a barrel section ($`|\mathrm{cos}\theta |<0.83`$) and two endcaps ($`0.83<|\mathrm{cos}\theta |<0.95`$). These components all operate inside a superconducting solenoid coil which creates a uniform magnetic field of 1.5 T. A muon detection system surrounds the solenoid coil with 1 m of iron absorber interspersed with Iarocci tube wire chambers operated in proportional mode.

We combine RICH and $`dE/dx`$ information to determine if a track appears more consistent with a pion or kaon identity. The deviation of the measured $`dE/dx`$ for any track from that expected under particle hypothesis $`i`$ ($`i=e`$, $`\mu `$, $`\pi `$, $`K`$, $`p`$) at that momentum, in units of the expected Gaussian width of the distribution, is defined as $`\delta _i`$. The $`dE/dx`$ resolution is about 6% for hadrons. For each detected charged particle track, the RICH detector response is condensed into a $`\chi _i^2`$-like variable for each particle hypothesis. The value of $`\chi _i^2`$ is derived from the number $`N_\gamma ^i`$ of detected Cherenkov photons and their locations relative to the Cherenkov cone expected for a particle with that momentum and mass. Pions and kaons are distinguished from one another with the combined RICH-$`dE/dx`$ variable $`\mathrm{\Delta }\chi ^2=\chi _\pi ^2\chi _K^2+\delta _\pi ^2\delta _K^2`$. We identify a track as a pion (kaon) if $`|\mathrm{cos}\theta |<0.83`$, $`N_\gamma ^\pi (N_\gamma ^K)3`$, $`\mathrm{\Delta }\chi ^2<(>)0`$, and apply these criteria to three-prong-side tracks. For the pion in $`\tau ^{}K^{}K^+\pi ^{}\nu _\tau `$, $`|\delta _\pi |<3.0`$ alone is used, yielding $``$10% higher efficiency from a looser criterion as well as the region outside the RICH acceptance, $`0.83<|\mathrm{cos}\theta |<0.93`$. The consequently higher $`K`$-fakes-$`\pi `$ rate is unimportant for this mode due to the dearth of three-kaon events.

We take advantage of the copious cascade decays $`D^+D^0\pi ^+`$ with $`D^0K^{}\pi ^+`$, in which the charged kaon and pion can be tagged kinematically without reference to RICH or $`dE/dx`$ information, to measure the probability for pions and kaons to be correctly or incorrectly identified. This method measures the area above background of the $`K^{}\pi ^+`$ mass peak with and without the particle identification criteria to determine the efficiency and fake rates for both pions and kaons. We find that the efficiency for pion and kaon identification ranges from 80$``$95% over most of the momentum range of interest for this analysis (0.5$``$3 GeV/$`c`$), and that the probability for a pion to fake a kaon or vice versa ranges from 1$``$2% for momentum values of 0.5$``$2 GeV/$`c`$ and rises to about 10$``$15% for momenta around 3 GeV/$`c`$, as shown in Fig. 1.

Events are required to have four well-measured charged particles in a one-versus-three topology (the two hemispheres being defined by the event thrust axis) and zero net charge. We use 1-prong $`\tau `$ leptonic ($`e/\mu `$) or hadronic decays ($`\rho /\pi `$) to “tag” the event. An electron tag requires the lone track to have $`\delta _e>2`$ and $`E/p`$, the ratio of energy deposited in the calorimeter to track momentum, to satisfy $`0.85<E/p<1.1`$. These criteria attain $`95\%`$ efficiency and a pion-faking-electron probability less than 1%. A muon tag must have the lone track penetrate at least three (if $`p<`$2 GeV/$`c`$) or five (if $`p>`$2 GeV/$`c`$) absorption lengths of material, which achieves 80$``$90% efficiency and a few percent pion-faking-muon probability for $`p>`$1 GeV/$`c`$.

If the tag-side track is not identified as an electron or a muon, tag-side calorimeter showers are then examined. If the invariant mass of the lone track and showers in the tagging hemisphere not matched to that track is below 0.5 GeV/$`c^2`$, the event is classified as $`\pi `$-tag. To establish a $`\rho ^{}(\pi ^{}\pi ^0)`$-tag, we require a $`\pi ^0\gamma \gamma `$ candidate on the 1-prong side. For photon candidates, we select showers which do not match the projection of any charged track trajectories into the calorimeter, have energies exceeding 30 (50) MeV in the barrel (endcaps), and have compact lateral profiles consistent with those expected from electromagnetic showers. At least one photon from any $`\pi ^0`$ candidate decay must be located in the barrel. We take the most energetic $`\gamma \gamma `$ pair which has invariant mass within $`(3.0,2.5)`$ standard deviations of $`m_{\pi ^0}`$ (typically a window $``$30 MeV/$`c^2`$ wide). After kinematically constraining the $`\gamma \gamma `$ pair to $`m_{\pi ^0}`$, the invariant mass of this fitted $`\pi ^0`$ candidate and the associated charged track, assigned the pion mass, is required to be consistent with the $`\rho `$ mass: $`0.5<m(\pi ^{}\pi ^0)<1.0`$ GeV/$`c^2`$.

To suppress backgrounds with one or more neutral particles, we do not allow any shower in the event unassociated with a charged track or a $`\pi ^0`$ participating in $`\rho `$-tagging with energy exceeding 100 MeV. After the event selection, the backgrounds with extra neutrals contribute 4.1%, 2.1% and 2.6% of the total observed events for the channels with 0, 1 and 2 kaons, respectively. In order to reject two-photon backgrounds (characterized by the missing momentum along the beam direction and low visible energy), the missing momentum must point into the CLEO detector ($`|\mathrm{cos}\theta _{miss}|<0.95`$) and the visible energy must exceed 20% of the $`e^+e^{}`$ center-of-mass energy. For all modes, we further demand that the three-prong mass, calculated using mass assignments determined from particle identification, lie below the $`\tau `$ mass. For the decay $`\tau ^{}K^{}\pi ^+\pi ^{}\nu _\tau `$, a veto of $`K_S^0`$ candidates with $`K_S^0\pi ^+\pi ^{}`$ within 10 MeV/$`c^2`$ of the $`K_S^0`$ nominal mass and a production vertex more than 1 cm from the primary interaction point is applied. Any $`K_s^0`$ content in both modes with two or more pions is treated as background, not signal.

Efficiencies and the remaining backgrounds are evaluated with Monte Carlo (MC) events from the KORALB korb  ($`\tau `$-pairs) and JETSET jetset  ($`q\overline{q}`$ hadrons) generators passed through the GEANT-based GEANT  CLEO detector simulation. Any MC event which has a tag-side particle which is truly $`e^\pm `$, $`\mu ^\pm `$, $`\pi ^\pm `$, $`K^\pm `$, or $`\rho ^\pm `$ is considered signal for purposes of efficiency; other tags are considered background even if the three-prong side is a signal decay mode. For $`\tau ^{}K^{}K^+K^{}\nu _\tau `$ decay, we generate signal MC samples according to phase space. The particle identification simulation models the general features of the efficiency and fake rates reasonably well (at the few percent level). We correct the Monte Carlo efficiencies and fake rates by hand, using momentum-dependent scale factors derived from measured and MC particle identification rates.

The corrected efficiencies for events in which one $`\tau `$ decays into a tag mode and the other into the relevant signal mode (integrated over all tag modes, which constitute about 72.4% of all $`\tau `$ decays), observed numbers of events in the data, and the expected numbers of $`\tau `$ and $`q\overline{q}`$ background events are given in Table 1. There are substantial signals evident in all but the three-kaon mode. The largest backgrounds come from $`\tau `$ decays to three charged hadrons with a single misidentification of a pion or kaon. Many such events are eliminated by the restriction of three-prong mass to be below the $`\tau `$ mass. Each decay in Table 1 contributes the dominant background for the channel below it due to the particle misidentification. For the $`\tau ^{}K^{}K^+K^{}\nu _\tau `$ decay, the $`\tau `$ background is not taken from the Monte Carlo sample, but rather the measured $`\pi `$-fakes-$`K`$ rate is applied to the $`K^{}K^+\pi ^{}`$ events in the data (many of which subsequently fail the three-prong mass restriction).

In Fig. 2, we show the three-hadron masses for the four modes considered here along with the results from the overlaid MC simulation. The two-particle substructure is presented in Fig. 3. Good agreement between the data and MC is observed mytune ; pi2pi0 . The $`\tau ^{}K^{}K^+\pi ^{}\nu _\tau `$ mode can be used to measure the contribution from the Wess-Zumino anomaly, offering an alternative channel to $`\tau ^{}\eta \pi ^{}\pi ^0\nu _\tau `$. A quantitative estimate requires determination of the structure functions prd47 , for which a substantially larger data sample would be needed.

Most sources of systematic uncertainty affect the four modes similarly. MC statistics for signal events add 1% uncertainty. The uncertainty in the number of $`\tau `$-pair events produced, and therefore in the normalization of the branching fractions, derives from the uncertainty of luminosity lumins  (2%) and $`\tau `$-pair cross section (2%) korb . Studies comparing track-finding efficiencies in the MC simulation to data result in a 0.5% per track uncertainty, added linearly to become 2% for the entire event. Several effects are probed individually by tightening and/or relaxing an event selection criterion, recalculating the branching fractions, and assigning an uncertainty in how well the MC models reality based on the changes induced: allowing for an extraneously reconstructed track in signal events adds 2% uncertainty; one systematic error arising from the assumed absence of two-photon physics background is 0.5% which can be tested by imposing a harsh net transverse momentum of the three-prong side requirement of 3 GeV/$`c`$; uncertainty due to the background levels from decay modes with extra neutrals allowed in is estimated by relaxing the veto energy on extra calorimeters showers to as high as 250 MeV and found to be 1%; and variations of the remaining non-particle-identification selection criteria contribute an error of 1%. Uncertainty in $`\tau `$ branching fractions results in 1% error due to changes in $`\tau `$ backgrounds. The efficiencies of particle identification were studied exhaustively for variations with time (0.5%) or charged-track multiplicity (0.5%), bias of the efficiency-measuring technique (0.5%), the data/MC efficiency correction procedure (0.4%), and statistics of the $`D^{}`$ data samples (0.3%), which all sum in quadrature to 1% per pion or kaon, which then are added linearly for the event as a whole to 3% overall. When the branching fractions are broken out by individual tags, consistent results are obtained, verifying both the modeling of tag efficiencies and the lack of unexpected non-$`\tau `$ backgrounds which sometimes show up in the non-leptonic tags. Together these errors sum in quadrature to 5%, which applies to the branching fractions of all four modes.

Other systematics vary from mode to mode, which will be referred to in the order they appear in Table 1. To estimate uncertainties from $`q\overline{q}`$ backgrounds, we compare event totals above the $`\tau `$ mass between the data and $`q\overline{q}`$ MC, and assign an additional $``$50% (relative) uncertainty from the observed differences. This results in uncertainties of 0.2%, 2%, 1%, and 3% for the four modes, respectively. An overall scale uncertainty in the fake rate per pion or kaon of 0.25% (absolute) is estimated by comparing measured results to those predicted by the RICH simulation, cross-checking with MC for possible bias in the fake-rate measurement technique, examining the rate of the forbidden $`\tau ^{}K^+\pi ^{}\pi ^{}\nu _\tau `$ (which is dominated by $`\tau ^{}\pi ^+\pi ^{}\pi ^{}\nu _\tau `$ with a pion faking a kaon), and estimating the error from the momentum-reweighting correction procedure. This fake rate uncertainty propagates to the branching fractions differently because of different levels of feed-down: 0.1%, 9%, 2%, and 12%, respectively. The $`K^{}K^+\pi ^{}`$ mode substructure is not well understood, and variations in the possibilities result in a 2% uncertainty in the efficiency. The three kaon mode substructure is completely unknown and assumed to occur via phase space with no assigned error. The remaining modes appear to be adequately described by the models used. Adding these contributions together in quadrature with the common sources results in total relative systematic uncertainties of 5%, 10%, 6%, and 14%, respectively.

For all decay modes, the yields are obtained from the numbers of the observed events with the backgrounds subtracted, the measured branching fractions (with $`K_S^0`$ excluded from the modes with two or more pions) are

$`(\tau ^{}\pi ^{}\pi ^+\pi ^{}\nu _\tau )`$ $`=`$ $`(9.13\pm 0.05\pm 0.46)\%,`$
$`(\tau ^{}K^{}\pi ^+\pi ^{}\nu _\tau )`$ $`=`$ $`(3.84\pm 0.14\pm 0.38)\times 10^3,`$
$`(\tau ^{}K^{}K^+\pi ^{}\nu _\tau )`$ $`=`$ $`(1.55\pm 0.06\pm 0.09)\times 10^3,`$
$`(\tau ^{}K^{}K^+K^{}\nu _\tau )`$ $`<`$ $`3.7\times 10^5\mathrm{at}90\%\mathrm{CL}.`$
The errors are statistical and systematic, respectively. For the decay $`\tau ^{}K^{}K^+K^{}\nu _\tau `$, the upper limit is obtained according to the procedure described in the PDG PDG  with both statistical and systematic errors taken into accout. This represents the first direct measurement of $`(\tau ^{}\pi ^{}\pi ^+\pi ^{}\nu _\tau )`$ with three pions explicitly identified, and our result is in good agreement with the PDG fitted value derived from the PDG constrained fit to the measurements of $`\tau ^{}(3h)^{}\nu _\tau `$ decays with and without identified charged kaons PDG . The measurement of $`(\tau ^{}K^{}K^+\pi ^{}\nu _\tau )`$ is consistent with the previous value PDG  and improves upon the precision previously attained by a factor of two. Our result $`(\tau ^{}K^{}\pi ^+\pi ^{}\nu _\tau )`$ is more precise than and consistent with both the CLEO cleoK2pi  and OPAL opalK2pi  results, while it is larger than the ALEPH alephK2pi  value by nearly three standard deviations. The upper limit for the Cabibbo-suppressed decay $`\tau ^{}K^{}K^+K^{}\nu _\tau `$ has been improved over the PDG value PDG  by a factor of five. Taken together, these measurements significantly enhance our understanding of the kaon content of three-prong $`\tau `$ decays.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. This work was supported by the National Science Foundation, the U.S. Department of Energy, the Research Corporation, and the Texas Advanced Research Program.


format: text

text:
# Integrable Markov processes and quantum spin chains

## Abstract

A set of Markov processes corresponding to systems of hard-core particles interacting along the line are shown to be solvable via a dynamic matrix product ansatz (DMPA). We show that quantum spin Hamiltonians can be treated by the DMPA as well, and demonstrate how the DMPA, originally formulated for systems with open ends, works for periodic systems.

MCS numbers: 60J25, 37K10, 17B69

Recent advances in nonequilibrium statistics have allowed for the clarification of many important issues such as exact description of nonequilibrium phase transitions, microscopic structure of shocks, dynamical scaling etc. . Dynamical systems of stochastically hopping hard core particles (the so-called exclusion processes) have been at the heart of this development, in particular due to remarkable observation that one can find exactly the probabilities of all system configurations in a stationary state. It was subsequently noticed that for the totally symmetric exclusion process one can as well solve the evolution dynamics, namely to find conditional probabilities $`P(\nu ,t|\nu ^0,t^0)`$ of the system being at state $`\nu `$ at time $`t`$ provided that it was in state $`\nu ^0`$ at time $`t^0`$ via coordinate Bethe Ansatz . The dynamical matrix product ansatz (DMPA) proposed in  has recovered the results of  in a more straighforward manner. But for technical reasons to be explained later, application of the DMPA was incomplete in periodic systems.

In this communication, we show how the DMPA is applied for periodic systems. Alongside, we demonstrate that the DMPA can be used for testing the stochastic processes and/or one-dimensional spin chains for integrability. The integrability restrictions come from the consistency conditions on the algebra generated by the DMPA, which has to satisfy the Yang-Baxter equation. The important difference from the conventional approach is that we start directly from the Hamiltonian and do not need a parameter-dependent transfer matrix, commuting with it.

The continuous time dynamics of Markov process is defined via the master equation for the probability $`P(\nu ,t)`$ of the system being at state $`\nu `$ at time $`t`$,

$$P(\nu ,t)/t=\underset{\nu ^{}}{}P(\nu ^{},t)w_{\nu ^{}\nu }P(\nu ,t)\underset{\nu ^{}}{}w_{\nu \nu ^{}}$$
(1)
where $`w_{\nu \nu ^{}}`$ is the rate with which the system can change its configuration from $`\nu `$ to $`\nu ^{}`$.

For definiteness we consider a system of identical particles on a lattice of $`L`$ sites. Each site $`k`$ can be either occupied by a particle (local occupation number $`n_k=1`$) or be empty ( $`n_k^0=1`$). The set of occupation numbers for all sites constitutes the configuration of the system. The system can change its configuration via hopping of any particle to an empty nearest neighbour site. At any infinitesimally small moment of time $`dt`$ at most one such process can happen, namely the particle can hop to the empty site to the right/to the left with the rates $`D_R`$ and $`D_L`$ respectively. Now, consider a vector space $`𝒲=^LC^2`$. Notice that one can establish one-to-one correspondence between configurational space of our particle system and canonical basis of $`𝒲`$; $`|\eta 𝒲`$, with the scalar product $`\eta ^{}|\eta =\delta _{\eta \eta ^{}}`$. The local particle number operators have form $`n_k=I\mathrm{}Ie_{22}\mathrm{}I`$ , and $`e_{11}`$ respectively for holes. $`e_{ij}`$ denotes a matrix with $`(e_{ij})_{kl}=\delta _{ik}\delta _{jl}`$, and it stays in the $`k`$-th place in the tensor product expression above. $`I`$ is the unit $`2\times 2`$ matrix.

With the help of above notations the set of equations (1) can be rewritten as single vector equation ( we ask the reader to refer e.g. to  for the details of the scheme).

$$\frac{d|P(t)}{dt}=H|P(t)$$
(2)
where $`|P(t)=_{\eta 𝒲}P_\eta (t)|\eta `$, $`H_{\eta \eta ^{}}=w_{\eta ^{}\eta },H_{\eta \eta }=_\eta ^{}w_{\eta \eta ^{}}`$.

Due to the fact that the hopping is restricted to the nearest-neighbouring sites, the Hamiltonian above is the sum of local two-site Hamiltonians, $`H=_{n=1}^Lh_{n,n+1}`$ given by

$$4h_{n,n+1}=D_R\left((1\sigma _n^z)(1+\sigma _{n+1}^z)\sigma _n^+\sigma _{n+1}^{}\right)+D_L\left((1+\sigma _n^z)(1\sigma _{n+1}^z)\sigma _n^{}\sigma _{n+1}^+\right)$$
(3)
where $`\sigma ^z,\sigma ^\pm =\sigma ^x\pm i\sigma ^y`$ are Pauli matrices.

Note that (2) has a form of a Schrödinger equation in imaginary time, and the Hamiltonian (3) the form of a nonhermitian spin Hamiltonian with nearest neighbour exchange interaction.

The global Hamiltonian $`H`$ possesses an U(1) symmetry, $`[H,_n\sigma _n^z]=0`$. The general hermitian spin Hamiltonian satisfying the latter relation is $`XXZ`$ Hamiltonian,

$$H=J\underset{n}{}\sigma _n^x\sigma _{n+1}^x+\sigma _n^y\sigma _{n+1}^y+\mathrm{\Delta }\left(\sigma _n^z\sigma _{n+1}^zI\right).$$
(4)
Thus Schrödinger quantum problem and Markov process for interacting particles are closely connected. In the following we show how to find the spectrum of $`H`$ and how to single out integrable models.

To demonstrate how the DMPA works, we start from the Hamiltonian of the form

$$H=\underset{n=1}{\overset{L}{}}h_{n,n+1};h=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& f& u& 0\\ 0& v& g& 0\\ 0& 0& 0& 0\end{array}\right)$$
(5)
with arbitrary entries $`f,g,u,v`$, which contains both the stochastic (3) and quantum spin system (4). We shall consider a system on a ring with $`L`$ sites. One proposes to search for time-dependent probability distribution vector $`|P(t)`$ in the following form:

$$|P(t)\text{Tr }\left[\left(\begin{array}{c}E\\ D\end{array}\right)^{L1}\left(\begin{array}{c}EQ\\ DQ\end{array}\right)\right]$$
(6)
where $`E,D,Q`$ are some time-dependent matrices. The probability of a given configuration $`\nu `$ is given by the scalar product of the above vector with the basis vector $`|\nu `$ of the configuration, $`P(\nu ,t)=\nu |P(t)`$. For example, probability $`P_\eta (t)`$ of a state $`\eta =\{\mathrm{0\; 0}A0A\}`$ for a system of 5 sites is sought in the form $`P_\eta (t)\text{Tr }\left(EEDEDQ\right)`$, up to normalization.

One requires matrices $`E,D,Q`$ such that (6) satisfies the master equation (2). To simplify matters, first take $`Q`$ in (6) to be a unit matrix $`Q=I`$ (this would allow to retrieve the translationally invariant part of the $`H`$ spectrum as we shall see below). Given that, and that the Hamiltonian in (2) splits into local ones, $`H=h_{1,2}+h_{2,3}+\mathrm{}h_{L1,L}+h_{L,1}`$, it is sufficient for the purpose that the matrices $`E,D`$ satisfy quadratic relation

$$\left(\frac{1}{2}\frac{d}{dt}+𝐡\right)\left(\genfrac{}{}{0pt}{}{E}{D}\right)\left(\genfrac{}{}{0pt}{}{E}{D}\right)=\left(\genfrac{}{}{0pt}{}{X^0}{X}\right)\left(\genfrac{}{}{0pt}{}{E}{D}\right)\left(\genfrac{}{}{0pt}{}{E}{D}\right)\left(\genfrac{}{}{0pt}{}{X^0}{X}\right)$$
(7)
conveniently written in short form

$$\left(\frac{1}{2}\frac{d}{dt}+h\right)|𝐀|𝐀=|𝐗|𝐀|𝐀|𝐗.$$
(8)
where $`𝐀|=(ED)`$, $`𝐗|=(X^0X)`$.

For $`QI`$ which is generic case one gets two additional relations in addition to (8), connected with the action of the terms $`h_{L1,L},h_{L,1}`$ of the Hamiltonian:

$`\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}+h\right)|𝐀|𝐀Q`$ $`=`$ $`|𝐗|𝐀Q|𝐀|𝐗Q`$ (9)
$`\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}+h\right)|𝐀Q|𝐀`$ $`=`$ $`|𝐗Q|𝐀|𝐀Q|𝐗`$ (10)
Equation (8) gives the relations for the matrices $`E,D`$ while Eqs. (9,10) impose restrictions on $`Q`$.

Equation (7) has four components:

$`{\displaystyle \frac{1}{2}}\left(\dot{E}E+E\dot{E}\right)`$ $`=`$ $`X^0EEX^0`$ (11)
$`{\displaystyle \frac{1}{2}}\left(\dot{E}D+E\dot{D}\right)+fED+uDE`$ $`=`$ $`X^0DEX`$ (12)
$`{\displaystyle \frac{1}{2}}\left(\dot{D}E+D\dot{E}\right)+gDE+vED`$ $`=`$ $`XEDX^0`$ (13)
$`{\displaystyle \frac{1}{2}}\left(\dot{D}D+D\dot{D}\right)`$ $`=`$ $`XDDX`$ (14)
Using the reduction

$$\dot{E}=0;X^0=0,$$
(15)
Eq. (11) is satisfied trivially. Supposing additionally that $`E^1`$ exists, we obtain from (12), (13):

$$2X=(gf)D+vEDE^1uE^1DE;$$
(16)
$$\dot{D}+(g+f)D+vEDE^1+uE^1DE=0$$
(17)
Substituting $`X,\dot{D}`$ into (14), one gets:

$$(g+f)DD+uDE^1DE+vEDE^1D=0$$
(18)
Now, consider the two equations (9,10) with $`Q`$ present. First, suppose that $`Q`$ is time-independent, and $`Q^1`$ exists. Then, Eq.(9) is satisfied with the choices (16-18) since multiplication of (9) by $`Q^1`$ from the right reduces (9) to (8). In order to satisfy Eq.(10) with the choice (15-18),the matrix $`Q`$ must obey

$$[Q,DE^1]=[Q,E^1D]=0.$$
(19)
Note that in general we assume that $`Q`$ does not commute neither with $`D`$ nor with $`E`$ separately, $`[Q,E]0;[Q,D]0`$.

Proceeding along the lines , we can get rid of time-dependence of $`D`$ via a formal Fourier Transform

$$𝒟_p(t)=\underset{k}{}\text{e}^{ipk}E^{k1}D(t)E^k$$
(20)
Indeed, it follows from (20) that

$$E^1𝒟_p(t)E=\text{e}^{ip}𝒟_p(t);$$
(21)
Substituting (20), (21) into (17), one gets

$$\dot{𝒟}_p+ϵ_p𝒟_p=0;𝒟_p(t)=e^{ϵ_pt}𝒟_p(0)$$

with the “dispersion relation”

$$ϵ_p=g+f+v\text{e}^{ip}+u\text{e}^{ip}$$
(22)
The inverse Fourier Transform is

$$D(t)=E\frac{1}{2\pi }_\pi ^\pi \text{d}p𝒟_p(t)=E\frac{1}{2\pi }_\pi ^\pi \text{d}p\text{e}^{ϵ_pt}𝒟_p(0)$$
(23)
In what follows we shall omit the time argument in initial matrices $`𝒟_p(0)`$. Inserting (23) into (14), and using (21), one obtains:

$$\frac{E^2}{2\pi }_\pi ^\pi _\pi ^\pi \text{e}^{(ϵ_{p_1}+ϵ_{p_2})t}𝒟_{p_1}𝒟_{p_2}a(p_1,p_2)\text{d}p_1\text{d}p_2=0$$
(24)
where

$$a(p_1,p_2)=u\text{e}^{ip_1+ip_2}+(f+g)\text{e}^{ip_1}+v$$
(25)
Now, integral in (24) may be split in two parts as

$$_\pi ^\pi _\pi ^\pi \mathrm{}\text{d}p_1\text{d}p_2=_\pi ^\pi _\pi ^{p_1}\mathrm{}+_\pi ^\pi _{p_1}^\pi \mathrm{}=I_{p2<p1}+I_{p2>p1}$$
(26)
By changing the order of integration and interchanging $`p_1p_2`$ in the last term in (26) we obtain from (24):

$$_\pi ^\pi _\pi ^{p_1}\left(a(p_1,p_2)𝒟_{p_1}𝒟_{p_2}+a(p_2,p_1)𝒟_{p_2}𝒟_{p_1}\right)\text{e}^{(ϵ_{p_1}+ϵ_{p_2})t}\text{d}p_1\text{d}p_2=0$$

Since this is satisfied at all times $`t`$, we must require the expression inside the integrand to vanish

$$𝒟_{p_1}𝒟_{p_2}=\frac{a(p_2,p_1)}{a(p_1,p_2)}𝒟_{p_2}𝒟_{p_1}$$
(27)
Finally, we require the matrix $`Q`$ from (9,10) to satisfy $`[Q,E^{k1}DE^k]=0`$, which contains (19) as a special case and yields

$$[Q,𝒟_p(t)]=0$$
(28)
To obtain Bethe Ansatz equations for the spectrum, consider the quantity

$$\text{Tr }(𝒟_{p_1}𝒟_{p_2}\mathrm{}𝒟_{p_N}E^LQ),$$
(29)
which is the building block of the expression for expectation values. Indeed, for a system with just two particles the expectation value $`n_x(t)n_y(t)`$, according to (6), (21),(23), can be written as

$`n_x(t)n_y(t)=\text{Tr }(E^{x1}DE^{yx1}DE^{Ly}Q)/Z_L=`$ (30)
$`{\displaystyle \text{d}p_1𝑑\text{p}_2\text{e}^{(ϵ_{p_1}+ϵ_{p_2})t}\text{e}^{ip_1xip_2y}\text{Tr }(𝒟_{p_1}𝒟_{p_2}E^LQ)/Z_L},`$ (31)
where $`Z_L`$ is an appropriate normalization factor, $`Z_L=_{x<y=1,2,\mathrm{}L}n_x(t)n_y(t)`$. Generally, one is interested in a system with $`N`$ particles which leads to (29).

Using (27), and commuting $`𝒟_{p_k}`$ around the circle, one obtains using (21), (28) and the cyclic invariance property of the trace:

$$\text{Tr }(𝒟_{p_1}𝒟_{p_2}\mathrm{}𝒟_{p_N}E^LQ)=(1)^{N+1}\underset{n=1}{\overset{N}{}}\frac{a(p_n,p_k)}{a(p_k,p_n)}e^{ip_kL}\text{Tr }(𝒟_{p_1}𝒟_{p_2}\mathrm{}𝒟_{p_N}E^LQ)$$

and consequently

$$(1)^{N+1}\underset{n=1}{\overset{N}{}}\frac{a(p_n,p_k)}{a(p_k,p_n)}\text{e}^{ip_kL}=1$$
(32)
One recognizes in (32), (22) the well-known expressions for the spectrum of the $`XXZ`$ spin model, for $`f=g=2J\mathrm{\Delta }`$, $`u=v=2J`$.

It may seem that the proposed line of argument is not changed if the auxiliary matrix $`Q`$ that we have introduced is chosen to be a unity matrix $`Q=I`$. This trivial choice however would bring about an additional relation. Namely let us rewrite (29) as

$$\text{Tr }(𝒟_{p_1}𝒟_{p_2}\mathrm{}𝒟_{p_N}E^L)=\text{Tr }(𝒟_{p_1}𝒟_{p_2}\mathrm{}𝒟_{p_N}EE^{L1})$$
(33)
and then commute single matrix E around the circle. We shall obtain

$$\text{e}^{ip_1+ip_2+\mathrm{}ip_N}=1$$
(34)
then. This means putting $`Q=I`$ (or more generally allowing for $`[Q,E]=0`$) we recover only the part of spectrum with the total quasimomentum equal to zero, $`p_1+p_2+\mathrm{}p_N=0`$.

We have thus shown how from the Dynamis Matrix Product Ansatz one obtains the spectrum for spin $`1/2`$ $`XXZ`$ chain. Note that until now we have not got any restriction on parameters in (5). This is a consequence of the fact that the Hamiltonian (5) is integrable for any choice of parameters. In the following we show how the the integrability restrictions enter in a more general setting.

Consider nearest neighbour quantum spin 1 Hamiltonians. This case corresponding to three states at a site, e.g. $`n_i=0,1,2`$ comprises many integrable systems (see e.g. a list given in  where however several symmetry restrictions are imposed). With those restrictions being relaxed, there are many more, e.g. .

To be specific, we shall consider the hermitian Hamiltonian of the form

$$h=\left(\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& a& 0& g_A& 0& 0& 0& 0& 0\\ 0& 0& b& 0& 0& 0& g_B& 0& 0\\ 0& g_A& 0& a& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& d& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& c& 0& g& 0\\ 0& 0& g_B& 0& 0& 0& b& 0& 0\\ 0& 0& 0& 0& 0& g& 0& c& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& e\end{array}\right)$$
(35)
with $`8`$ independent parameters. If $`g=g_B=g_A=1`$, $`d=e=0`$ and $`a=b=c=1`$, it reduces to isotropic Hamiltonian of Sutherland  $`h=PI`$, where $`P`$ is permutation operator in $`SU(3)SU(3)`$ defined by $`P(AB)=(BA)P`$ for any matrices $`A,B`$ from $`SU(3)`$. Again, the complete Hamiltonian has the form $`H=_{n=1}^Lh_{n,n+1}`$ and we are interested in a solution of Schrödinger equation $`i\frac{d}{dt}|P=H|P`$.

Analogously to (6) one proposes an Ansatz for $`|P(t)`$ as

$$|P(t)\text{Tr }\left[\left(\begin{array}{c}E\\ D^A\\ D^B\end{array}\right)^{L1}\left(\begin{array}{c}EQ\\ D^AQ\\ D^BQ\end{array}\right)\right]$$
(36)
with matrices $`D^A,D^B`$ referring to local variable $`n_i=0,1,2`$ respectively.

Ansatz (36) leads to the equations (8,9,10) where now $`𝐀|=(ED^AD^B)`$, $`𝐗|=(X^0X^AX^B)`$. Eq.(8) has 9 components, written analogously to (11-14). The first one is satisfied with the choice (15). Those four that are linear in $`E`$, yield analogically to (16,17):

$$\dot{D}^A+2aD^A+g_A\left(E^1D^AE+ED^AE^1\right)$$
(37)
$$2X^A=g_A\left(ED^AE^1E^1D^AE\right);$$

and the Eq. for $`\dot{D}^B,X^B`$ is obtained by exchanging $`ab,AB`$. The additional Eqs.(9,10) lead to the conditions $`[Q,E^{k1}D^Z(t)E^k]=0`$, $`k=0,\pm 1,\mathrm{},Z=A,B`$ on $`Q`$-matrix.

Eliminating time derivatives in the four remaining equations that are quadratic in $`D^Z`$ we obtain:

$`(2ad)D^AD^A+g_A\left(D^AE^1D^AE+ED^AE^1D^A\right)=0`$ (38)
$`(a+bc)D^AD^BgD^BD^A+g_AED^AE^1D^B+g_BD^AE^1D^BE=0`$ (39)
$`(a+bc)D^BD^AgD^AD^B+g_AED^BE^1D^A+g_AD^BE^1D^AE=0`$ (40)
$`(2be)D^BD^B+g_B\left(D^BE^1D^BE+ED^BE^1D^B\right)=0`$ (41)
One proceeds analogically then to what we have done in the spin $`1/2`$ case. Introducing generalized Fourier Transform $`𝒟_p^A(t)=_k\alpha ^k\text{e}^{ipk}E^{k1}D^A(t)E^k`$, $`𝒟_p^B(t)=_k\beta ^k\text{e}^{ipk}E^{k1}D^B(t)E^k`$ leads to eliminating the time-dependence in (37): $`\dot{𝒟}_p^Z+ϵ_p𝒟_p^Z=0;𝒟_p^Z(t)=e^{ϵ_p^Zt}𝒟_p^Z(0);Z=A,B,`$ where

$$ϵ_p^A=2a+g_A(\alpha \text{e}^{ip}+\text{e}^{ip}/\alpha ),ϵ_p^B=2b+g_B(\beta \text{e}^{ip}+\text{e}^{ip}/\beta ).$$
(42)
The inverse Fourier Transform $`D^{A,B}(t)=E\frac{1}{2\pi }_\pi ^\pi \text{d}p𝒟_p^{A,B}(t)`$, being inserted into (38-41), yields

$`0`$ $`=`$ $`{\displaystyle \text{d}p_1\text{d}p_2a_{12}𝒟_{p_1}^A𝒟_{p_2}^A\text{e}^{(ϵ_{p_1}^A+ϵ_{p_2}^A)t}}`$ (43)
$`0`$ $`=`$ $`{\displaystyle \text{d}p_1\text{d}p_2[c_{12}𝒟_{p_1}^A𝒟_{p_2}^Bg\beta \text{e}^{ip_2}𝒟_{p_2}^B𝒟_{p_1}^A]\text{e}^{(ϵ_{p_1}^A+ϵ_{p_2}^B)t}}`$ (44)
$`0`$ $`=`$ $`{\displaystyle \text{d}p_1\text{d}p_2[d_{12}𝒟_{p_1}^B𝒟_{p_2}^Ag\alpha \text{e}^{ip_2}𝒟_{p_2}^A𝒟_{p_1}^B]\text{e}^{(ϵ_{p_1}^B+ϵ_{p_2}^A)t}}`$ (45)
$`0`$ $`=`$ $`{\displaystyle \text{d}p_1\text{d}p_2b_{12}𝒟_{p_1}^B𝒟_{p_2}^B\text{e}^{(ϵ_{p_1}^B+ϵ_{p_2}^B)t}}`$ (46)
with the functions

$`a_{12}a(p_1,p_2)`$ $`=`$ $`\alpha \text{e}^{ip_1}(2ad)+g_A(1+\alpha ^2\text{e}^{ip_1+ip_2})`$ (47)
$`b_{12}b(p_1,p_2)`$ $`=`$ $`\beta \text{e}^{ip_1}(2be)+g_B(1+\beta ^2\text{e}^{ip_1+ip_2})`$ (48)
$`c_{12}c(p_1,p_2)`$ $`=`$ $`(a+bc)\alpha \text{e}^{ip_1}+g_B\alpha \beta \text{e}^{ip_1+ip_2}+g_A`$ (49)
$`d_{12}d(p_1,p_2)`$ $`=`$ $`(a+bc)\beta \text{e}^{ip_1}+g_A\alpha \beta \text{e}^{ip_1+ip_2}+g_B`$ (50)
The trick with the $`p_1`$ and $`p_2`$ exchange in (43-46) leads consequently (see (27) and the paragraph above it) to

$`𝒟_{p_1}^A𝒟_{p_2}^A={\displaystyle \frac{a_{21}}{a_{12}}}𝒟_{p_2}^A𝒟_{p_1}^A;`$ $`𝒟_{p_1}^B𝒟_{p_2}^B={\displaystyle \frac{b_{21}}{b_{12}}}𝒟_{p_2}^B𝒟_{p_1}^B`$ (51)
$`c_{12}𝒟_{p_1}^A𝒟_{p_2}^Bg\beta \text{e}^{ip_2}𝒟_{p_2}^B𝒟_{p_1}^A`$ $`+`$ $`c_{21}𝒟_{p_2}^A𝒟_{p_2}^Bg\beta \text{e}^{ip_1}𝒟_{p_1}^B𝒟_{p_2}^A=0`$ (52)
$`d_{12}𝒟_{p_1}^B𝒟_{p_2}^Ag\alpha \text{e}^{ip_2}𝒟_{p_2}^A𝒟_{p_1}^B`$ $`+`$ $`d_{21}𝒟_{p_2}^B𝒟_{p_1}^Ag\alpha \text{e}^{ip_1}𝒟_{p_1}^A𝒟_{p_2}^B=0`$ (53)
and a restriction $`g_B=\alpha /\beta g_A;\alpha =\pm \beta `$. The latter one is a consequence of the constraint $`ϵ_{p_1}^A+ϵ_{p_2}^B=ϵ_{p_1}^B+ϵ_{p_2}^A`$, needed to equate the time-dependent factors after $`p_1p_2`$ exchange in (44,45).

Eqs. (52,53) can be written in the form

$`𝒟_{p_1}^A𝒟_{p_2}^B=𝒟_{p_2}^B𝒟_{p_1}^A{\displaystyle \frac{d_{12}g\beta \text{e}^{ip_2}d_{21}g\beta \text{e}^{ip_1}}{c_{12}d_{12}g^2\alpha \beta \text{e}^{2ip_1}}}+𝒟_{p_2}^A𝒟_{p_1}^B{\displaystyle \frac{g^2\alpha \beta \text{e}^{ip_1+ip_2}d_{12}c_{21}}{c_{12}d_{12}g^2\alpha \beta \text{e}^{2ip_1}}}`$ (54)
$`𝒟_{p_1}^B𝒟_{p_2}^A=𝒟_{p_2}^A𝒟_{p_1}^B{\displaystyle \frac{c_{12}g\alpha \text{e}^{ip_2}c_{21}g\alpha \text{e}^{ip_1}}{c_{12}d_{12}g^2\alpha \beta \text{e}^{2ip_1}}}+𝒟_{p_2}^B𝒟_{p_1}^A{\displaystyle \frac{g^2\alpha \beta \text{e}^{ip_1+ip_2}c_{12}d_{21}}{c_{12}d_{12}g^2\alpha \beta \text{e}^{2ip_1}}}`$ (55)
The Eqs.(51,54,55) between the Fourier components $`𝒟_p^A,𝒟_p^B`$ define an algebra conveniently written via a scattering matrix $`S`$ as

$$𝒟_{p_1}^Z𝒟_{p_2}^Z^{}=S_{YY^{}}^{ZZ^{}}(p_1,p_2)𝒟_{p_2}^Y^{}𝒟_{p_1}^Y$$
(56)
Applying it twice we see that it satisfies the inversion relation ($`S^T`$ means transposition) $`S(p_1,p_2)S^T(p_2,p_1)=1`$.

Now, note that application of (56) exchange the order of indexes $`p_1`$ and $`p_2`$. For the cubic terms $`𝒟_{p_1}^X𝒟_{p_2}^Y𝒟_{p_3}^Z`$ the exchange of order of indexes $`123321`$ can be done either by subsequent series of $`12,13,23`$ or $`23,13,12`$ exchanges. Associativity of the algebra requires the result to be the same in both cases, which amounts to

$$S_{i^{\prime \prime }j^{\prime \prime }}^{ij}(p_1,p_2)S_{i^{}k^{\prime \prime }}^{i^{\prime \prime }k}(p_1,p_3)S_{j^{}k^{}}^{j^{\prime \prime }k^{\prime \prime }}(p_2,p_3)=S_{j^{\prime \prime }k^{\prime \prime }}^{jk}(p_2,p_3)S_{i^{\prime \prime }k^{}}^{ik^{\prime \prime }}(p_1,p_3)S_{i^{}j^{}}^{i^{\prime \prime }j^{\prime \prime }}(p_1,p_2).$$
(57)
This is the Yang-Baxter equation (YBE). Summing over repeating indexes is implied. The YBE guarantees the consistency of the algebra on the cubic and further levels.

The YBE impose severe restrictions on the $`S`$ matrix elements, which depend on the coefficients of the Hamiltonian (35). Nontrivial solutions of the YBE single out integrable models. We have investigated (57) with the help of Mathematica, and found the following solution:

$`\alpha =i_1\beta ,g_B=i_1g_A,`$ $`g_A=i_2g,d=2a+(i_31)g,`$ (58)
$`e=2b+(i_41)g,`$ $`c=a+b+i_5g,`$ (59)
where $`i_1,i_2,i_3,i_4,i_5=\pm 1`$.

Sutherland’s  and, more generally, solution $`\mathrm{\#}`$ 4 of  is contained in the above. To the best of our knowledge, the general manifold describes a new solution to the YBE.

To summarise, the Dynamic Matrix Product State method turns out to be an effective tool to check the integrability of both stochastic processes and quantum Hamiltonians. Applying it saves the conventional tedious step of calculating parameter-dependent transfer matrices satisfyng the the Yang-Baxter equation, commuting with the original Hamiltonian. E.g. for investigating spin $`1`$ chains conventionally, we would need first to solve YBE with in principle $`27^2`$ independent components. In the DMPA approach we needed only the original Hamiltonian to start with, and then additionally to solve at most $`8^2`$ equations (57) for the algebra associativity. (In fact the number of equations is much smaller due to symmetries). At the same time, we believe that for all the integrable models obtained the commuting transfer-matrices do exist, although we do not need them directly. Also, Bethe Ansatz equations for the spectrum can be constructed directly from the algebra, as was demonstrated for $`XXZ`$ model (see (32)). For the other models the Bethe Ansatz equations are readily constructed  by following e.g. the algebraic approach .

The DMPA approach does not make reference to the actual representations of the matrices involved. However, it would be interesting to find them. In special cases (namely, when they are not time dependent, which means only stationary (ground) state is investigated) the representations are known .


text:
# Proton number fluctuation as a signal of the QCD critical end-point

(January 2003)

## Abstract

We argue that the event-by-event fluctuation of the proton number is a meaningful and promising observable for the purpose of detecting the QCD critical end-point in heavy-ion collision experiments. The long range fluctuation of the order parameter induces a characteristic correlation between protons which can be measured. The proton fluctuation also manifests itself as anomalous enhancement of charge fluctuations near the end-point, which might be already seen in existing data.

The event-by-event fluctuations in heavy-ion collisions carry information about the degrees of freedom of the created system and their correlations ga . In particular, thermodynamic properties of QCD can be inferred from event-by-event fluctuation measurements C\_V ; srs ; srslong ; jk ; asakawa ; jeon .

Of particular interest are fluctuations originating from the QCD critical end-point barducci ; berges ; srs ; srslong ; berdnikov ; fodor ; hatta . Since the fluctuation of the order parameter induces characteristic correlations among particles, in particular pions, it is expected that the end-point affects the event-by-event fluctuations of certain observables in a nontrivial way srs ; srslong ; misha .

Here we discuss a new observable which may serve as a signal of the end-point; the event-by-event fluctuation of the net proton number, i.e., the number of the protons minus the number of antiprotons observed 111Similar conclusions apply to the fluctuations of the number of protons or antiprotons separately. The net proton number appears more naturally in the theoretical analysis, and we focus on it for clarity..

Our starting point is the fact that the baryon number susceptibility $`\chi _B`$ larry ; kunihiro ; go ; gottlieb ; bernard ; gavai  diverges at the critical end-point berges ; srs ; gavin ; hatta . $`\chi _B`$ is related to the average magnitude of the fluctuation $`\delta B`$ of the baryon number:

$`\chi _B={\displaystyle \frac{1}{VT}}(\delta B)^2`$ (1)
where $`V`$ and $`T`$ are the volume and the temperature. The divergence of $`\chi _B`$ is a consequence of the fact that the critical point is the end-point of a line of first order phase transitions, which are characterized, in particular, by a jump in the baryon number density $`B/V`$ 222Note that the rise of $`\chi _B`$ near $`T_c`$ at $`\mu _B=0`$ gottlieb ; bernard  is not due to the critical fluctuations, but to the liberation of QCD degrees of freedom kunihiro ; go ; gavai . In fact, $`\chi _B`$ is finite at $`\mu _B=0`$ and $`T_c`$..

If, in a heavy-ion collision experiment, we could measure all the baryons, the enhancement of the event-by-event fluctuation of the baryon number in a given subvolume would be a signature of the end-point. However, about one half of the emitted baryons are undetected neutrons which certainly contribute to the fluctuation of the baryon number. To what extent does the proton number fluctuation alone reflect the divergence of $`\chi _B`$?

This letter is devoted to clarifying where, in the observed quantities, the divergence occurs and advocating the proton number fluctuation as a sensible and promising observable for the search of the critical point in the heavy-ion experiments.

In this work, we confine ourselves to equilibrium thermodynamic fluctuations. Various important issues such as the nonequilibrium evolution of the fluctuations will be (and some already have been) studied separately.

For simplicity and clarity we shall work in QCD with exact isospin invariance. The relevant corrections due to isospin breaking are small as we discuss below. Let us first show that in this case the isospin number susceptibility, $`\chi _I`$, is finite at the end-point. The proof is based on the fact that the singular behavior of thermodynamic quantities near the critical point is due to the divergence of a certain correlation length. It is the correlation length in the $`\sigma `$-channel, the channel with quantum numbers of the chiral condensate $`\overline{\psi }\psi `$ berges ; srs . A density-density correlator, such as $`\chi _I=(1/T)d^3xV_0(x)V_0(0)`$ can diverge only if the density can mix with the $`\sigma `$ field. For the isospin density this mixing is strictly forbidden by the SU(2)V (isospin) symmetry. The isospin density, $`V_0(x)`$, transforms as a triplet, $`\mathrm{𝟑}`$. On the other hand, $`\sigma `$ is a singlet. The mixing is forbidden and there is no singular contribution in $`\chi _I`$. 333For a similar reason the axial isospin susceptibility remains finite at $`T=T_c`$, $`\mu _B=0`$ muia .

Small explicit breaking of the SU(2)V symmetry by the quark mass difference $`m_um_d`$ or the isospin chemical potential $`\mu _I`$ will induce singularity in $`\chi _I`$, since $`m_um_d`$ and $`\mu _I`$ are SU(2)V triplets 444In the sense that they couple to SU(2)V triplet fields. The argument here is similar to kogut ., and can produce $`V_0\sigma `$ mixing. In the context of heavy ion collisions corresponding singular contributions are negligible.

We summarize by writing the singular parts of the baryon and isospin number susceptibilities:

$`\chi _B\xi ^{\gamma /\nu },\chi _I=0\text{(singular parts only)},`$ (2)
where $`\xi `$ is the divergent correlation length of the sigma field: $`\xi =1/m_\sigma `$. The zero in eq. (2) neglects small isospin breaking terms as well as finite terms. The universal values of the exponents are given by $`\gamma 1.2`$, $`\nu 0.63`$, $`\alpha 0.12\gamma `$. Note that $`\gamma /\nu =2\eta `$, where $`\eta 0.04`$.

Let us discuss the implications of (2). In particular, let us consider charge susceptibility. Study of charge fluctuations in heavy-ion collisions has attracted much attention recently. It was proposed that these fluctuations might reflect thermodynamic conditions earlier in the collision history, due to charge conservation asakawa ; jeon . The charge fluctuation magnitude per entropy is a measure of the charge per particle or degree of freedom. In the quark-gluon plasma (QGP) the charge per degree of freedom is smaller. If the charge fluctuation is equilibrating too slowly, the observed value will be smaller than the equilibrium magnitude which can be calculated using the resonance gas jk . However, the estimates of the charge diffusion shuryak  suggest that for the rapidity windows achievable in present experiments charge diffusion is very efficient in equilibrating charge fluctuations, thus practically washing out the “history” effects. Experimentally, the QGP suppression of the charge fluctuation is not seen reid ; spsdata , which is consistent with the diffusion estimates shuryak . The effect of the critical fluctuations we are discussing here is crucially different from the QGP suppression. While the latter is the “history” effect, the critical fluctuations are the equilibrium fluctuations pertaining to the freezeout point, and the diffusion is necessary to establish them.

The measure of charge fluctuations, the charge number susceptibility, $`\chi _Q`$, can be expressed in terms of $`\chi _B`$ and $`\chi _I`$ using the relation $`Q=B/2+I_3`$ and the fact that isospin symmetry requires $`\delta B\delta I_3=0`$:

$`\chi _Q={\displaystyle \frac{1}{VT}}(\delta Q)^2={\displaystyle \frac{1}{4}}\chi _B+\chi _I.`$ (3)
Eq. (2) then implies that the charge susceptibility diverges at the critical point, due to the divergence of $`\chi _B`$.

We now wish to relate the susceptibilities $`\chi _B`$, $`\chi _I`$ and $`\chi _Q`$ to observable particle number fluctuations. For simplicity, we shall limit our discussion by considering only protons, neutrons and pions. Accounting for other particles will not alter our conclusions. In the hadron language, the susceptibilities may be written as

$`\chi _B`$ $`=`$ $`{\displaystyle \frac{1}{VT}}\left(\delta N_{p\overline{p}}+\delta N_{n\overline{n}}\right)^2`$
$`\chi _I`$ $`=`$ $`{\displaystyle \frac{1}{VT}}\left({\displaystyle \frac{1}{2}}\delta N_{p\overline{p}}{\displaystyle \frac{1}{2}}\delta N_{n\overline{n}}+\delta N_{\pi ^+\pi ^{}}\right)^2`$ (4)
and

$`\chi _Q`$ $`=`$ $`{\displaystyle \frac{1}{VT}}\left(\delta N_{p\overline{p}}+\delta N_{\pi ^+\pi ^{}}\right)^2`$ (5)
where we introduced notation $`N_{p\overline{p}}N_pN_{\overline{p}}`$ for the net proton number fluctuation, with $`\delta `$ denoting event-by-event deviation from the equilibrium value. Similar notations are used for neutrons and pions.

Now we concentrate on singular parts of the susceptibilities and ask a question: what does (2) imply for the individual particle number fluctuations? It is easy to check that the following set of relations between singular parts of the particle correlators reproduces the correct singular behavior given in (2):

$`\delta N_{p\overline{p}}\delta N_{p\overline{p}}=\delta N_{n\overline{n}}\delta N_{n\overline{n}}=\delta N_{p\overline{p}}\delta N_{n\overline{n}}`$
$`\delta N_{\pi ^+\pi ^{}}\delta N_{\pi ^+\pi ^{}}=0`$
$`\delta N_{p\overline{p}}\delta N_{\pi ^+\pi ^{}}=\delta N_{n\overline{n}}\delta N_{\pi ^+\pi ^{}}=0`$
$`\text{(singular parts only)}.`$ (6)
Some of these equations follow trivially from isospin invariance, but some, for instance, the last equation on the first line and that on the second line, require a stronger condition. Such relations occur naturally if we attribute the divergences to the exchange of a sigma meson, which is an isospin singlet. Using eqs. (6) we obtain

$`\chi _B`$ $`=`$ $`{\displaystyle \frac{4}{VT}}\delta N_{p\overline{p}}\delta N_{p\overline{p}},`$
$`\chi _I`$ $`=`$ $`0,`$
$`\chi _Q`$ $`=`$ $`{\displaystyle \frac{1}{VT}}\delta N_{p\overline{p}}\delta N_{p\overline{p}}.`$ (7)
(singular parts only)
Remarkably, the singular part of the charge fluctuation comes from the protons. In other words, had we considered only contributions from charged pions in $`\chi _Q`$, the singular parts of $`\pi ^+\pi ^+`$, $`\pi ^{}\pi ^{}`$, $`\pi ^+\pi ^{}`$ correlators (all are singular at the critical point srslong ) would have canceled each other. We see also that the proton number fluctuation completely reflects the singularity of the baryon number susceptibility, which justifies its use as a sensible probe of the QCD critical end-point.

To provide a simple estimate of how large the net proton number fluctuation can become near the critical point, we begin by calculating the correlator

$`\delta n_𝒑\delta n_𝒌,`$ (8)
where $`n_𝒑`$ is the net proton number in the momentum bin labeled by the value $`𝒑`$. In addition to the usual statistical fluctuation, the correlator (8) receives a contribution from the effective interaction with the sigma field $`\sigma `$, $`_{\sigma pp}=g\sigma \overline{P}P`$, where $`g`$ is the dimensionless sigma-nucleon coupling and $`P`$ is the Dirac field of a proton. All fluctuation observables of the protons can be constructed from (8) srs ; misha .

Near the critical point, the singular term in (8) is represented by a diagram of forward proton-proton scattering. A straightforward calculation following misha  gives,

$`V\delta n_𝒑\delta n_𝒌`$ $`=`$ $`{\displaystyle \frac{g^2}{m_\sigma ^2T}}{\displaystyle \frac{4m^2}{E_𝒑E_𝒌}}\left[n_𝒑^+(1n_𝒑^+)n_𝒑^{}(1n_𝒑^{})\right]`$ (9)
$`\times `$ $`\left[n_𝒌^+(1n_𝒌^+)n_𝒌^{}(1n_𝒌^{})\right]`$
(singular part only)
where $`m=940`$ MeV is the proton mass, $`E_𝒑=\sqrt{𝒑^2+m^2}`$ and $`n_𝒑^\pm =\left[\mathrm{exp}\left\{(E_𝒑\mu _B)/T\right\}+1\right]^1`$, while $`m_\sigma =1/\xi `$ is the sigma meson (screening) mass.

Let us compare the singularity in (9) to the singularity in (2). The exponent $`\gamma /\nu =2\eta `$ in (2) is very close to 2 (the anomalous dimension $`\eta `$ of the $`\sigma `$ field is small) and is equal to 2 in the mean field approximation ($`\eta =0`$), which is the same as the power of $`1/m_\sigma `$ in (9).

In a realistic heavy-ion collision environment finiteness of the space-time volume severely prevents $`m_\sigma `$ from vanishing exactly srslong ; berdnikov . The smallest achievable value is estimated to be around $`(3\text{ fm})^1`$.

One possible concern is that the rescattering in the final hadronic stage washes out critical point fluctuations. In this respect, one should bear in mind that the rescattering in question includes the exchange of the $`\sigma `$ quanta, which, near the critical point, is the source of the fluctuations we consider. The critical fluctuations are washed out if the final (kinetic) freezeout occurs sufficiently far from the critical point. In order to see the effect, one should dial control parameters (e.g., reduce the size of the ions to raise the freezeout temperature) to bring the freezeout closer to the critical point.

In (9), $`g`$ is taken at zero momentum transfer, i.e., off the sigma mass shell. In vacuum, $`gm/f_\pi `$ ($`f_\pi =93`$ MeV is the pion decay constant) is quite large $``$10. For the quantitative estimate below, we assume that $`g`$ does not change appreciably from its vacuum value near the chiral phase transition. (See, however, hatsuda .)

First, let us assume that Au-Au collisions at RHIC at $`\sqrt{s}=130`$ GeV froze out in the vicinity of the end-point and estimate the effect of (9) in terms of the unknown mass $`m_\sigma `$, which is the measure of the proximity of the end-point. At the chemical freezeout, $`T=174`$ MeV and $`\mu _B=46`$ MeV braun . Integrating over $`𝒑`$ and $`𝒌`$, we obtain the net proton number fluctuation $`(\delta N_{p\overline{p}})^2`$ and divide it by the sum of proton and antiproton numbers $`N_{p+\overline{p}}`$:

$`{\displaystyle \frac{(\delta N_{p\overline{p}})^2}{N_{p+\overline{p}}}}|_{\mathrm{RHIC}}1.0+0.062\left({\displaystyle \frac{g}{10}}\right)^2\left({\displaystyle \frac{200\text{MeV}}{m_\sigma }}\right)^2.`$ (10)
The unity on the r.h.s. is the trivial statistical contribution. Taking $`g=10`$ and $`m_\sigma =60`$ MeV $`(3\text{ fm})^1`$ srslong ; berdnikov , we find $`(\delta N_{p\overline{p}})^2/N_{p+\overline{p}}1.7`$.

On the other hand, if the end-point were located at a value of $`\mu _B`$ of order of a few hundred MeV, as inferred from simple model estimates berges  and suggested by the recent lattice simulation fodor , it is possible that the SPS freezout is in the proximity of the critical point. Using the freezeout parameters ($`T,\mu _B`$)=(168 MeV, 266 MeV) at SPS sps  we obtain

$`{\displaystyle \frac{(\delta N_{p\overline{p}})^2}{N_{p+\overline{p}}}}|_{\mathrm{SPS}}1.0+1.5\left({\displaystyle \frac{g}{10}}\right)^2\left({\displaystyle \frac{200\text{MeV}}{m_\sigma }}\right)^2.`$ (11)
Note that the coefficient of the second term has a much bigger value than in (10). This is because the singular term given by (9) grows as the square of the net proton number and also because there is a partial cancellation between protons and antiprotons at RHIC. We stress that the main feature in (10) and (11) is the singular dependence on $`m_\sigma `$, which makes the effect large when the freezeout occurs near the critical point. There are other effects, which contribute to the r.h.s. of (10) and (11), but which are not singular near the critical point (e.g., initial volume fluctuations caused by impact parameter fluctuations).

Experimentally, separating protons and measuring proton fluctuations is a feasible task in the RHIC as well as the SPS detectors. We hope that such data analysis will be available soon.

In order to test our ideas on the existing data we can, using (10) and (11), estimate the contribution of the proton fluctuation to the total charge fluctuation characterized by $`D4(\delta Q)^2/N_{\mathrm{tot}}`$ jeon . As (7) shows, pions do not contribute to the singular part of the charge fluctuation, but they dilute such a contribution of the protons. Using $`N_{\pi ^++\pi ^{}}10N_{p+\overline{p}}`$ and $`(\delta N_{\pi ^+\pi ^{}})^2/N_{\pi ^++\pi ^{}}10.3=0.7,`$ where the negative contribution $`0.3`$ is due to the resonance decays jk , we obtain

$`{\displaystyle \frac{D}{4}}{\displaystyle \frac{(\delta Q)^2}{N_{\mathrm{tot}}}}`$ $``$ $`{\displaystyle \frac{(\delta N_{p\overline{p}})^2+(\delta N_{\pi ^+\pi ^{}})^2}{N_{p+\overline{p}}+N_{\pi ^++\pi ^{}}}}`$ (12)
$``$ $`0.8,`$
where we have neglected the cross terms between $`\delta N_p`$ and $`\delta N_{\pi ^\pm }`$. We see that the fluctuation anomaly in the proton sector can result in a larger charge fluctuation than the resonance gas value $`0.7`$ by about 10$`\%`$. At SPS this effect is even stronger. At SPS, using (11), $`g=10`$ and $`m_\sigma =200`$ MeV we get $`(\delta N_{p\overline{p}})^2/N_{p+\overline{p}}2.5`$, and with $`N_{\pi ^++\pi ^{}}5N_{p+\overline{p}}`$ we obtain: $`D/41`$.

Before these estimates can be compared to experiment, one must take into account the effect of limited acceptance of a given detector. It is easy to see that this effect reduces deviations from $`D=4`$. Its estimates range from few percent corrections jeon  to almost complete elimination of deviations from $`D=4`$ zaranek , depending on the assumptions on the rapidity correlator of fluctuations and the width of the acceptance window. We do not discuss these issues here, and refer the reader to the literature.

Experimentally, the data from RHIC suggest that the magnitude of the fluctuation is slightly larger than a thermodynamical fluctuation in a resonance gas reid . This effect is even more pronounced at SPS spsdata . There are, of course, a number of possible explanations, for example, (i) acceptance, as we have just discussed zaranek ; (ii) remnant initial state correlations; (iii) decay of multiply charge clusters; (iv) other nonequilibrium fluctuations (e.g., by a mechanism similar to gavin ); In this letter we wish to point out that, quite independently of these other effects, an equilibrium critical fluctuation due to the proximity of the end-point could explain the enhancement of charge fluctuations observed at RHIC and SPS. The independent measurement of proton fluctuations that we suggested would be necessary to confirm and sufficient to rule out this effect.

Is it possible that the light sigma effect is seen in both RHIC and SPS experiments? If this happens, the region of $`\mu _B`$ where $`m_\sigma `$ is small ($`<200`$ MeV) is rather wide – of order hundred MeV. Although unlikely, this might not be completely unnatural if one takes into account the fact that $`m_\sigma `$ is suppressed on the crossover line stretching from the end-point to $`\mu _B=0`$ axis, even though it vanishes only at the end-point (see, e.g., Fig. 5 in hatta ).

As it should be clear from our discussion, measuring the charge fluctuations is not the most efficient way to search for the end-point, although the effect may be seen in such observables too. A direct measurement of the proton number fluctuation as a function of the $`\sqrt{s}`$ of the collision is both feasible and is less afflicted by other effects. Correlation of such a measurement with other proposed signatures of the critical point (such as $`p_t`$ fluctuations srs ) would affirm the discovery of the QCD critical point.

In conclusion, protons carry both the baryon and the electric charges. They are sensitive to the fluctuation of the order parameter. Due to the peculiar nature of the end-point – isospin blindness of the sigma field – the singularity of the baryon number susceptibility is completely reflected in the proton number fluctuation. Thus the net proton number fluctuation is a very useful observable. By studying the $`\mu _B`$ dependence of this fluctuation one may discover and determine the location of the critical point on the phase diagram of QCD.

We are grateful to L. McLerran and R. Pisarski for many discussions as well as to M. Asakawa, M. Kaneta and T. Kunihiro for comments. Y. H. thanks RIKEN BNL Center for hospitality during the completion of this work. M.S. thanks RIKEN BNL Center and U.S. Department of Energy \[DE-AC02-98CH10886\] for providing facilities essential for the completion of this work. The work of M.S. is supported in part by DOE Outstanding Junior Investigator Grant and by the Alfred P. Sloan Foundation.