text:
# Dissipation and noise in adiabatic quantum pumps

## I Introduction

A resent experiment by Switkes et al.  has stimulated increasing interest in adiabatic quantum charge pumping. Idealy in such an experiment one aims at generating a dc-current by slowly modulating the shape of a mesoscopic conductor with the help of oscillating gate voltages. A single potential oscillating at frequency $`\omega `$ does not generate a dc-current, but two potentials oscillating with the same frequency but out of phase can generate a dc-current. The effect is of interest under conditions in which electron motion is phase-coherent and is thus termed quantum pumping. The frequency of the potential modulation is small compared to the characteristic times for traversal and reflection of electrons and the pump is thus termed adiabatic. Thus carriers traversing the sample see an almost static potential. The last circumstance allows to give an elegant formulation of quantum pumping  which is based on the scattering matrix approach to low frequency ac transport in phase coherent mesoscopic systems .

Recently Avron et al.  investigated adiabatic quantum pumping with the aim to formulate criteria for an ”optimal pump”. The term ”optimal” means that such a pump is noiseless and transports integer charge in each cycle. To this extent they have investigated not only the dc-current but also the dissipation and the noise generated by a pump. Avron et al. express their results in terms of an energy shift matrix $`i\mathrm{}\widehat{s}/t\widehat{s}^{}`$ where $`\widehat{s}`$ is the time-dependend scattering matrix. This is an elegant formulation which gives a correct description of time-dependent adiabatic currents and dissipation. However, for quantities which invoke correlations at different times the approach is valid only for pump frequencies $`\mathrm{}\omega <<k_BT`$.

It is the purpose of this work to investigate the distribution function, heat flow and noise properties of an adiabatic pump for an arbitrary relation of pump frequency and temperature. To achieve this we start with the scattering matrix approach for ac-transport. This approach leads to expressions for the quantities of interest in terms of the side bands of particles exiting the pump. The side bands correspond to particles which have gained or lost a modulation quantum $`\mathrm{}\omega `$. In particular, the approach presented here allows a description of the quantum statistical correlation properties (noise) of an adiabatic quantum pump.

The adiabatic quantum pump  of interest here should be distinguished from a variety of other pumping mechanisms. For certains pumps  the charge transferred in each cycle is quantized. Quantized charge pumping is most easily achieved in devices based on the Coulomb blockade effect  where the charge on a quantum dot is quantized. This is of considerable metrological interest . Other effects which lead to pumping are the photovoltaic effect  and the acoustoelectric effect .

The paper is organized as follows. In Sec.II the essential assumptions we make are described. In Sec.III we calculate the nonequilibrium distribution function for the outgoing particles produced by the pump. In Sec.IV we formulate the condition which is necessary to pump dc-current. In Sec.V we calculate the heat flows produced by an oscillating mesoscopic scatterer. In Sec.VI we consider the shot noise produced by the pump and analyze the noise in terms of uncorrelated movement of nonequilibrium quasi-particles (quasi-electrons and holes) generated by the pump and correlations between them . In Sec.VII we present explicit results for the particular case of a two-leads scatterer with the time-reversal symmetry.

## II The main assumptions

To describe the response of a mesoscopic phase coherent sample to slowly oscillating (with a frequency $`\omega `$) external real parameters $`X_j(t)`$ (gate potential, magnetic flux, etc.)

$$X_j(t)=X_j+X_{\omega ,j}e^{i(\omega t\phi _j)}+X_{\omega ,j}e^{i(\omega t\phi _j)},$$
(1)
we will use the scattering matrix approach . The sample is connected via leads (which we will number via Greek letters $`\alpha `$, $`\beta `$, $`\gamma `$, etc.) to $`N_r`$ reservoirs. The scattering matrix $`\widehat{s}`$ being a function of parameters $`X_j(t)`$ depends on time. Two main assumption will be used. First, we suppose that the external parameter changes so slowly that we can apply an ”instant scattering” description using the scattering matrix $`\widehat{s}(t)`$ frozen at some time $`t`$. Physically this means that the scattering matrix changes only a little while an electron is scattered by the mesoscopic sample (i.e., the frequency $`\omega `$ is much smaller than the inverse Wigner time delay ). In this sense we use the term ”adiabatic” pump.

Second, we assume that the amplitude $`X_{\omega ,j}`$ is small enough to keep only the terms linear in $`X_{\omega ,j}`$ in an expansion of the scattering matrix

$$\widehat{s}(t)\widehat{s}+\widehat{s}_\omega e^{i\omega t}+\widehat{s}_{+\omega }e^{i\omega t}.$$
(2)
In the limit of small frequencies the amplitudes $`\widehat{s}_{\pm \omega }`$ can be expressed in terms of parametric derivatives of the on-shell scattering matrix $`\widehat{s}`$,

$$\widehat{s}_{\pm \omega }=\underset{j}{}X_{\omega ,j}e^{\pm i\phi _j}\widehat{s}/X_j.$$
(3)
The expansion Eq.(2) is equivalent to the nearest sidebands approximation  which implies that a scattered electron can absorb or emit only one energy quantum $`\mathrm{}\omega `$ before it leaves the scattering region.

The kinetic properties (charge current, heat current, etc.) which are of interest here depend on the values of the scattering matrix within the energy interval of the order of $`\mathrm{max}(k_BT,\mathrm{}\omega )`$ near the Fermi energy. In the low frequency ($`\omega 0`$) and low temperature ($`T0`$) limit we assume the scattering matrix to be energy independent.

## III Outgoing distribution function

In a pump setup the mesoscopic scatterer is coupled to reservoirs $`\alpha =1,2,\mathrm{}N_r`$ with the same temperatures $`T_\alpha `$ = $`T`$ and electrochemical potentials $`\mu _\alpha `$ = $`\mu `$. Thus electrons with the energy $`E`$ entering the scatterer are described by the Fermi distribution function

$$f_\alpha ^{(in)}(E)=f_0(E)=\frac{1}{1+e^{\frac{E\mu }{k_BT}}}.$$

Due to the interaction with an oscillating scatterer an electron can absorb or emit an energy quantum $`\mathrm{}\omega `$ that changes the distribution function. Our aim is to find the distribution function for outgoing particles (i.e., for electrons leaving the mesoscopic sample and entering the reservoir) far from the scatterer.

Let us consider a single transverse channel in one of the leads. We introduce two kinds of carriers . First, incoming particles which are going from the reservoir to the scatterer. And, second, outgoing particles which are leaving the scattering region. We can express the operators $`\widehat{b}_\alpha `$ which annihilate outgoing carriers in the lead $`\alpha `$ in terms of operators $`\widehat{a}_\beta `$ annihilating incoming electrons  in lead $`\beta `$. Applying the hypothesis of an instant scattering we can write

$$\widehat{b}_\alpha (t)=\underset{\beta }{}s_{\alpha \beta }(t)\widehat{a}_\beta (t).$$
(4)
Here $`s_{\alpha \beta }`$ is an element of the scattering matrix $`\widehat{s}`$; the time dependent operator is $`\widehat{a}_\alpha (t)`$ = $`𝑑E\widehat{a}_\alpha (E)e^{iEt/\mathrm{}}`$, and the energy dependent operators obey the following anticommutation relations

$$[\widehat{a}_\alpha ^{}(E),\widehat{a}_\beta (E^{})]=\delta _{\alpha \beta }\delta (EE^{}).$$

Note that above expressions correspond to single (transverse) channel leads and spinless electrons. For the case of many-channel leads each lead index ($`\alpha `$, $`\beta `$, etc.) includes a transverse channel index and any repeating lead index implies implicitly a summation over all the transverse channels in the lead. Similarly an electron spin can be taken into account.

Using Eq.(2) and Eq.(4) we obtain

$`\widehat{b}_\alpha (E)={\displaystyle \underset{\beta }{}}s_{\alpha \beta }\widehat{a}_\beta (E)`$ (5)
$`+s_{\omega ,\alpha \beta }\widehat{a}_\beta (E+\mathrm{}\omega )+s_{+\omega ,\alpha \beta }\widehat{a}_\beta (E\mathrm{}\omega ).`$ (6)
The distribution function for electrons leaving the scatterer through the lead $`\alpha `$ is $`f_\alpha ^{(out)}(E)`$ = $`<\widehat{b}_\alpha ^{}(E)\widehat{b}_\alpha (E)>`$, where $`<\mathrm{}>`$ means quantum-mechanical averaging. Substituting Eq.(6) we find

$`f_\alpha ^{(out)}(E)={\displaystyle \underset{\beta }{}}|s_{\alpha \beta }|^2f_0(E)`$ (7)
$`+|s_{\omega ,\alpha \beta }|^2f_0(E+\mathrm{}\omega )+|s_{+\omega ,\alpha \beta }|^2f_0(E\mathrm{}\omega ).`$ (8)
Note that the distribution function for outgoing carriers is a nonequilibrium distribution function generated by the nonstationary scatterer. The above expression gives a simple physical interpretation for the Fourier amplitudes of the scattering matrix. $`|s_{\omega ,\alpha \beta }|^2`$ ($`|s_{+\omega ,\alpha \beta }|^2`$) is the probability for an electron entering the scatterer through the lead $`\beta `$ and leaving the scatterer through the lead $`\alpha `$ to emit (to absorb) an energy quantum  $`\mathrm{}\omega `$. Note that $`|s_{\alpha \beta }|^2`$ is the probability for the same scattering without the change of an energy. Below we will use Eq.(8) to analyze the kinetics of a pump.

## IV dc-current

To be definite we take currents from the scatterer to the reservoirs to be positive. Using the distribution functions $`f_0(E)`$ for incoming electrons and $`f_\alpha ^{(out)}(E)`$ for outgoing electrons we find for the dc-current $`I_\alpha `$ in the lead $`\alpha `$ far from the scatterer

$$I_\alpha =\frac{e}{h}_0^{\mathrm{}}𝑑E[f_\alpha ^{(out)}(E)f_0(E)].$$
(9)
Substituting Eq.(8) we get

$$I_\alpha =\frac{e\omega }{2\pi }[T_{+\omega ,\alpha }T_{\omega ,\alpha }].$$
(10)
Here we have introduced the total probabilities for electrons scattered into the lead $`\alpha `$ (irrespective of the lead through which they entered the scattering region) to absorb $`T_{+\omega ,\alpha }`$ or to emit $`T_{\omega ,\alpha }`$ an energy quantum $`\mathrm{}\omega `$

$$T_{\pm \omega ,\alpha }=\underset{\beta }{}|s_{\pm \omega ,\alpha \beta }|^2.$$
(11)
We see that only a scatterer with the property

$$T_{+\omega ,\alpha }T_{\omega ,\alpha },$$
(12)
can pump current into the lead $`\alpha `$.

It is useful to express these probabilities in terms of a bare scattering matrix $`\widehat{s}`$. To this end we introduce a generalized parametric emissivity matrix $`\widehat{\nu }[X]`$

$$\widehat{\nu }[X]=\frac{1}{2\pi i}\frac{d\widehat{s}}{dX}\widehat{s}^{}.$$
(13)
with matrix elements

$$\nu _{\alpha \beta }[X]=\frac{1}{2\pi i}\underset{\gamma }{}\frac{ds_{\alpha \gamma }}{dX}s_{\beta \gamma }^{}.$$
(14)
The diagonal element $`\nu _{\alpha \alpha }[X]`$ of the parametric emissivity matrix  is the charge that leaves the sample through contact $`\alpha `$ in response to a variation of the parameter $`X`$. The non-diagonal element $`\nu _{\alpha \beta }[X]`$ ($`\alpha \beta `$) of the parametric emissivity matrix determines the correlations between current amplitudes generated in the contacts $`\alpha `$ and $`\beta `$ due to a variation of the parameter $`X`$ (see Eq.(29)).

Using Eq.(3) and Eq.(17) we express the probabilities $`T_{\pm \omega ,\alpha }`$ given by Eq.(11) in terms of the matrix elements of the parametric emissivity matrix (to lowest order in $`X_{\omega ,j}`$)

$$T_{\pm \omega ,\alpha }=4\pi ^2\underset{\beta }{}\left|\underset{j}{}X_{\omega ,j}e^{\pm i\phi _j}\nu _{\alpha \beta }[X_j]\right|^2.$$
(15)
The quantities $`T_{\pm \omega ,\alpha }`$ admit a simple interpretation in the quasi-particle picture. Due to scattering the electron system gains an energy from the nonstationary (oscillating) scatterer. Absorption of an energy quantum $`\mathrm{}\omega `$ leads to creation of a nonequilibrium (quasi-)electron-hole pair. Note that at any temperature $`T0`$ equilibrium electron-hole pairs exist. This nonequilibrium pair is neutral but transfers an energy $`\mathrm{}\omega `$. From Eq.(8) it follows that $`T_{+()\omega ,\alpha }`$ is proportional to the number of nonequilibrium quasi-electrons (holes) leaving the scattering region through the lead $`\alpha `$. The electron and hole (belonging to the same pair) can be scattered either into one lead (see Fig.1.a) or into different leads (see Fig.1.b). If they are scattered into the same lead they do not contribute to the current. But if they are scattered into different leads they do contribute. In any case they contribute to the heat transfer from the oscillating scatterer into the reservoirs.

The dc-current $`I_\alpha `$ in the lead $`\alpha `$ can be represented as a sum of two contributions $`I_\alpha `$ = $`I_\alpha ^{(e)}`$ \+ $`I_\alpha ^{(h)}`$, where $`I_\alpha ^{(e)}`$ = $`e\omega T_{+\omega ,\alpha }/(2\pi )`$ and $`I_\alpha ^{(h)}`$ = $`e\omega T_{\omega ,\alpha }/(2\pi )`$ are currents carried by nonequilibrium quasi-electrons and holes, respectively (here e (-e) is an electron (a hole) charge).

Now we will show that current is conserved, i.e., $`_\alpha I_\alpha =0`$. To this end we use the fact that the scattering matrix is unitary

$$\widehat{s}(t)\widehat{s}^{}(t)=1.$$
(16)
For the expansion Eq.(2) this leads to the relations

$$\underset{\gamma }{}[s_{\alpha \gamma }s_{\beta \gamma }^{}+s_{\omega ,\alpha \gamma }s_{\omega ,\beta \gamma }^{}+s_{+\omega ,\alpha \gamma }s_{+\omega ,\beta \gamma }^{}]=\delta _{\alpha \beta },$$
(17)
$$\underset{\gamma }{}s_{\alpha \gamma }s_{\omega ,\beta \gamma }^{}=\underset{\gamma }{}s_{\beta \gamma }^{}s_{+\omega ,\alpha \gamma },$$
(18)
$$\underset{\gamma }{}s_{\beta \gamma }^{}s_{\omega ,\alpha \gamma }=\underset{\gamma }{}s_{\alpha \gamma }s_{+\omega ,\beta \gamma }^{}.$$
(19)
Multiplying Eq.(18) and Eq.(19) by parts, summing the result over $`\alpha `$ and taking into account Eq.(11) we obtain (neglecting the higher powers of $`s_{\pm \omega ,\alpha \gamma }`$)

$$\underset{\alpha }{}T_{\omega ,\alpha }=\underset{\alpha }{}T_{+\omega ,\alpha }.$$
(20)
Using Eqs.(10) and (20) we see that the scatterer does not produce any current $`_\alpha I_\alpha `$ = 0 but it can only push a current from some reservoir to another reservoir.

An alternative (but equivalent) way to find the dc-current is to average the time-dependent current

$$I_\alpha =\underset{\mathrm{\Delta }t\mathrm{}}{lim}\frac{1}{\mathrm{\Delta }t}\underset{0}{\overset{\mathrm{\Delta }t}{}}𝑑t<\widehat{I}_\alpha (t)>.$$

The current operator is

$$\widehat{I}_\alpha (t)=\frac{e}{h}[\widehat{b}_\alpha ^{}(t)\widehat{b}_\alpha (t)\widehat{a}_\alpha ^{}(t)\widehat{a}_\alpha (t)].$$
(21)
Substituting Eqs.(2) and (4) into Eq.(21) and performing quantum mechanical and time averaging we obtain Eq.(10).

Note that in a pump setup where the external reservoirs are at the same macroscopic conditions (electrochemical potential, temperature, etc.) and a periodic in time perturbation is applied directly to the mesoscopic conductor there is no linear regime for dc transport (only ac currents are linear in perturbation). The dc-currents (charge, heat, etc.) are of a quantum mechanical nature and arise because of a nonlinear (quadratic) dependence on the quantum-mechanical (scattering) amplitudes (see Eq.(11)).

## V Heat flow

Particles traversing the sample absorb energy from a time dependent scatterer and carry it into the reservoirs. We assume that the reservoirs are large enough to absorb this energy and to remain still in thermal equilibrium. In the leads the energy is transferred by electrons only (we neglect any inelastic processes in the leads). Thus to calculate  an energy flow $`I_{E,\alpha }`$ entering the reservoir $`\alpha `$ we can use an electron distribution function and write

$$I_{E,\alpha }=\frac{1}{h}_0^{\mathrm{}}𝑑E(E\mu )[f_\alpha ^{(out)}(E)f_0(E)].$$
(22)
Substituting Eq.(8) we obtain

$$I_{E,\alpha }=\frac{\mathrm{}\omega ^2}{4\pi }[T_{+\omega ,\alpha }+T_{\omega ,\alpha }].$$
(23)
Comparing Eq.(10) and Eq.(23) we see that the time dependent scatterer always generates heat flows (because $`T_{\pm \omega ,\alpha }`$ are positively defined) and can be considered as a mesoscopic (phase-coherent) heat source which can be useful, for instance, for studying various thermoelectric phenomena in mesoscopic structures. In contrast the existence of a dc-current Eq.(10) requires a special condition (see Eq.(12)). Another difference is that the heat flow is directed (at any lead) from the scatterer to the reservoir (if all the reservoirs are at the same temperature) but the charge flow, if it exists, can be directed either from the reservoir to the scatterer (at some lead) or vice versa (at another lead) because of charge conservation.

The quasi-particle description gives a simple physical interpretation of Eq.(23). We can say that the heat is transported by two kinds of quasi-particles, the quasi-electrons and holes. Each quasi-particle has an energy $`\mathrm{}\omega /2`$ (on average). This is because the absorption of each energy quantum $`\mathrm{}\omega `$ creates two quasi-particles, a quasi-electron and a hole. Thus the heat (energy) transferred by quasi-electrons and holes is $`I_{E,\alpha }^{(e)}`$ = $`(\mathrm{}\omega /2)(I_\alpha ^{(e)}/e)`$ and $`I_{E,\alpha }^{(h)}`$ = $`(\mathrm{}\omega /2)(I_\alpha ^{(h)}/(e))`$, respectively (the quasi-electron $`I_\alpha ^{(e)}`$ and hole $`I_\alpha ^{(h)}`$ currents are defined in the previous section after Eq.(15)). The sum of these contributions gives Eq.(23)

## VI Current fluctuations

The problem of current noise in a quantum pump is closely connected with the problem of quantization of the charge pumped in one cycle . On the other hand the noise in mesoscopic phase coherent conductors is interesting in itself  because it is very sensitive to quantum-mechanical interference effects and can give additional information about the scattering matrix.

To describe the current-current fluctuations we will use the correlation function

$$S_{\alpha \beta }(t,t^{})=\frac{1}{2}<\mathrm{\Delta }\widehat{I}_\alpha (t)\mathrm{\Delta }\widehat{I}_\beta (t^{})+\mathrm{\Delta }\widehat{I}_\beta (t^{})\mathrm{\Delta }\widehat{I}_\alpha (t)>,$$
(24)
where $`\mathrm{\Delta }\widehat{I}`$ = $`\widehat{I}`$ \- $`<\widehat{I}>`$ and $`\widehat{I}_\alpha (t)`$ is the quantum-mechanical current operator in the lead $`\alpha `$ given by Eq.(21). Note that in the case of a time-dependent scatterer the correlation function depends on two times $`t`$ and $`t^{}`$.

Here we are interested in the noise averaged over a long time  ($`\mathrm{\Delta }t2\pi /\omega `$) and we investigate

$$S_{\alpha \beta }(t)=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{}}𝑑t^{}S_{\alpha \beta }(t,t^{}).$$

In addition we restrict our consideration to the zero-frequency component of the noise spectra $`S_{\alpha \beta }`$ = $`𝑑tS_{\alpha \beta }(t)`$. Substituting the current operator Eq.(21) and taking into account Eq.(2) and Eq.(4) we can write the zero-frequency noise power

$`S_{\alpha \beta }={\displaystyle \frac{2e^2}{h}}{\displaystyle _0^{\mathrm{}}}𝑑E<\widehat{S}_{\alpha \beta }(E,E)`$ (25)
$`+\widehat{S}_{\alpha \beta }(E,E\mathrm{}\omega )+\widehat{S}_{\alpha \beta }(E,E+\mathrm{}\omega )>.`$ (26)
Here

$$\widehat{S}_{\alpha \beta }(E,E^{})=\frac{1}{2}[\mathrm{\Delta }\widehat{I}_\alpha (E)\mathrm{\Delta }\widehat{I}_\beta (E^{})+\mathrm{\Delta }\widehat{I}_\beta (E^{})\mathrm{\Delta }\widehat{I}_\alpha (E)];$$

$`\mathrm{\Delta }\widehat{I}_\alpha (E)`$ = $`\widehat{I}_\alpha (E)`$ \- $`<\widehat{I}_\alpha (E)>`$ and

$$\widehat{I}_\alpha (E)=\widehat{b}_\alpha ^{}(E)\widehat{b}_\alpha (E)\widehat{a}_\alpha ^{}(E)\widehat{a}_\alpha (E).$$

For the energy independent sattering matrix in the lowest order in $`\widehat{s}_{\pm \omega }`$ we obtain $`S_{\alpha \beta }`$ = $`S_{\alpha \beta }^{(th)}`$ \+ $`S_{\alpha \beta }^{(pump)}`$. Here the thermal (or Nyquist-Johnson) noise is

$$S_{\alpha \beta }^{(th)}=2e^2k_BT/h[2\delta _{\alpha \beta }|s_{\alpha \beta }|^2|s_{\beta \alpha }|^2].$$

The noise power produced by the pump is

$$S_{\alpha \beta }^{(pump)}=\frac{2e^2}{h}F(\mathrm{}\omega ,k_BT)\left(\delta _{\alpha \beta }[T_{\omega ,\alpha }+T_{+\omega ,\alpha }]T_{\alpha \beta }^{(cor)}\right),$$
(27)
where

$$T_{\alpha \beta }^{(cor)}=\left|\underset{\gamma }{}s_{\beta \gamma }s_{\omega ,\alpha \gamma }^{}\right|^2+\left|\underset{\gamma }{}s_{\beta \gamma }s_{+\omega ,\alpha \gamma }^{}\right|^2,$$
(28)
and $`F(\mathrm{}\omega ,k_BT)`$ = $`\mathrm{}\omega \mathrm{coth}[\mathrm{}\omega /(2k_BT)]2k_BT`$.

In addition to the probabilities $`T_{\pm \omega ,\alpha }`$ which determine the dc current Eq.(10) and heat flow Eq.(23) there appears a third key quantity $`T_{\alpha \beta }^{(cor)}`$, which describes the effect of correlations between (quasi-)particles. Similarly to $`T_{\pm \omega ,\alpha }`$ (see Eqs.(15)) this probability can be expressed in terms of a generalized emissivity matrix $`\widehat{\nu }`$ (see Eq.(13))

$$T_{\alpha \beta }^{(cor)}=4\pi ^2\underset{\eta =+1,1}{}\left|\underset{j}{}X_{\omega ,j}e^{i\eta \phi _j}\nu _{\alpha \beta }[X_j]\right|^2.$$
(29)
Note that there is no summation over $`\alpha `$ or $`\beta `$: Consequently this probability which determines the current cross-correlation is directly proportional to the off-diagonal element of the emissivity matrix.

Now we will analyze the noise power Eq.(27). We can see that the current cross-correlations $`S_{\alpha \beta }^{(pump)}`$ ($`\alpha \beta `$) produced by the pump are negative: that is quite general for nonequilibrium noise in the system of fermions . The noise generated by the pump obeys the following sum rules  $`_\alpha S_{\alpha \beta }^{(pump)}`$ = $`_\beta S_{\alpha \beta }^{(pump)}`$ = 0. This is a straightforward consequence of an instant scattering description applied here. Indeed, using Eqs.(4),(16) and (21) we can see that the conservation law holds not only for a quantum averaged current $`<\widehat{I}_\alpha (t)>`$ but for a current operator as well : $`_\alpha \widehat{I}_\alpha (t)=0`$. Thus the current correlations $`S_{\alpha \beta }<\widehat{I}_\alpha \widehat{I}_\beta >`$ must obey the same sum rule.

The function $`F(\mathrm{}\omega ,k_BT)`$ describes the effect of thermal fluctuation on shot noise and determines the dependence of the noise on the pump frequency $`\omega `$. At sufficiently high temperature $`\mathrm{}\omega k_BT`$ the noise Eq.(27) is quadratic in $`\omega `$. This is in agreement with Ref.. But at low temperature $`k_BT\mathrm{}\omega `$ the noise is linear in $`\omega `$ and this is in agreement with the counting statistics calculations of Levitov .

Next consider the three terms in the brackets of the r.h.s. of Eq.(27). Consider the low temperature limit $`k_BT\mathrm{}\omega `$ and devide the expression for noise into two parts $`S_{\alpha \beta }^{(pump)}`$ = $`\delta _{\alpha \beta }S_\alpha ^{(pump),(P)}`$ \+ $`S_{\alpha \beta }^{(pump),(cor)}`$. The first part

$$S_\alpha ^{(pump),(P)}=\frac{e^2\omega }{\pi }[T_{\omega ,\alpha }+T_{+\omega ,\alpha }],$$

is due to an uncorrelated movement of nonequilibrium quasi-electrons and holes. To verify this we apply the Schottky formula  for shot noise $`S_{\alpha ,q}^{(Sch)}`$ = $`2qI_\alpha ^{(q)}`$ (here $`q`$ is a particle charge and the index ”q” means that the current is carried by the particles with the charge q). Substituting the current carried by the quasi-electrons $`I_\alpha ^{(e)}`$ = $`e\omega T_{+\omega ,\alpha }/(2\pi )`$ and by the holes $`I_\alpha ^{(h)}`$ = $`e\omega T_{+\omega ,\alpha }/(2\pi )`$ (here e(-e) is an electron (hole) charge) into Schottky’s formula we obtain $`S_\alpha ^{(pump),(P)}`$ = $`S_{\alpha ,e}^{(Sch)}`$ \+ $`S_{\alpha ,h}^{(Sch)}`$ (in the literature the Schottky result is referred as the Poisson value of shot noise that we indicate by the upper index ”P”).

The second part

$$S_{\alpha \beta }^{(pump),(cor)}=\frac{e^2\omega }{\pi }T_{\alpha \beta }^{(cor)},$$

is due to correlations between quasi-electrons and holes. These correlations are a consequence of the common origin of the electron and hole forming a pair and their subsequent scattering into different leads Fig.1.b. Thus we can say that the cross-correlations are exclusively due to dissolving (neutral) electron-hole pairs. Because of charge conservation this gives a simple explanation of a negative sign of cross-correlations in our case Eq.(27). Note that Schottky’s result gives no correlations between currents at different leads. Due to $`S_{\alpha \alpha }^{(pump),(cor)}`$ the current correlation at the same lead $`S_{\alpha \alpha }^{(pump)}`$ is below the Poisson value $`S_\alpha ^{(P)}`$ and the Fano factor characterizing the deviation of the actual shot noise from the Poisson noise (see, e.g., ) $`F`$ = $`S_{\alpha \alpha }^{(pump)}/S_\alpha ^{(P)}`$ is, in general, less than unity. We would like to emphasize that when we calculate the Poisson value of shot noise we do not use the total current $`I_\alpha `$ in the lead $`\alpha `$ but we calculate the sum of the Poisson noises produced by both the quasi-electrons (the current is $`I_\alpha ^{(e)}`$) and the holes (the current is $`I_\alpha ^{(h)}`$).

## VII Applications

In this section we consider a simple but a quite generic case of a two-terminal mesoscopic conductor with a time-reversal symmetry (without magnetic fluxes). In addition we assume that the external, time-dependent parameters $`X_j(t)`$ do not change the total charge on a sample (if this is not the case we need to take into account the self-consistent internal potential ). In this case we can choose the scattering matrix as follows (we consider spinless electrons)

$$\widehat{s}=\left(\begin{array}{cc}re^{i\theta }& it\\ it& re^{i\theta }\end{array}\right).$$
(30)
Here $`r^2`$ and $`t^2`$ are the reflection and the transmission probability, respectively ($`r^2+t^2=1`$). We assume the quantities $`r,t,\theta `$ to be the functions of external parameters $`X_j`$ (see Eq.(1)).

### A Heat flow and noise in one-parameter ”pumps”

If only one external parameter $`X`$ is varied then we get from Eq.(3) $`|s_{\omega ,\alpha \beta }|^2`$ = $`|s_{+\omega ,\alpha \beta }|^2`$ for any $`\alpha `$ and $`\beta `$. Thus the one-parameter adiabatic ”pump” gives the same probability for absorption and emission at both leads ($`\alpha =1,2`$)

$`T_{\omega ,\alpha }^{(1)}=T_{+\omega ,\alpha }^{(1)}=T_{\omega ,\alpha }^{(1)},`$ (31)
(32)
$`T_{\omega ,\alpha }^{(1)}[X]=X_\omega ^2{\displaystyle \underset{\beta }{}}\left|{\displaystyle \frac{ds_{\alpha \beta }}{dX}}\right|^2,`$ (33)
and does not produce a dc-current: $`I_\alpha =0`$ (see Eq.(10)). In contrast it does produce the heat flows Eq.(23) and the noise Eq.(27). The contribution to the noise Eq.(28) due to correlations between quasi-particles is

$$T_{\alpha \beta }^{(cor)}=8\pi ^2X_\omega ^2\left|\nu _{\alpha \beta }[X]\right|^2.$$
(34)
Here $`\nu _{\alpha ,\beta }[X]`$ is a matrix element of a generalized parametric emissivity matrix $`\widehat{\nu }[X]`$ Eq.(13). With the scattering matrix Eq.(30) we find a heat flow

$`I_{E,1}^{(1)}=I_{E,2}^{(1)}`$ (35)
$`={\displaystyle \frac{\mathrm{}\omega ^2X_\omega ^2}{2\pi }}\left(r^2\left({\displaystyle \frac{d\theta }{dX}}\right)^2+\left({\displaystyle \frac{dr}{dX}}\right)^2+\left({\displaystyle \frac{dt}{dX}}\right)^2\right),`$ (36)
and a noise $`S_{11}^{(pump)}`$ = $`S_{22}^{(pump)}`$ = $`S_{12}^{(pump)}`$ = $`S_{21}^{(pump)}`$ = $`S^{(1)}`$

$`S^{(1)}={\displaystyle \frac{4e^2X_\omega ^2}{h}}F(\mathrm{}\omega ,k_BT)`$ (37)
$`\times \left(r^2t^2\left({\displaystyle \frac{d\theta }{dX}}\right)^2+\left({\displaystyle \frac{dr}{dX}}\right)^2+\left({\displaystyle \frac{dt}{dX}}\right)^2\right).`$ (38)
We see that the noise produced by the one-parameter ”pump” Eq.(38) gives us direct information on the dependence of the scattering matrix on the varying external parameter $`X`$. This dependence $`\widehat{s}(X)`$ is important for calculating the current produced by the pump if two external parameters are varied Eq.(48). In a real experimental situation the dependence $`\widehat{s}(X)`$ is unknown and can not be calculated in a simple way. Thus the possibility to obtain this dependence from the experimental data seems useful.

Note that for some particular conditions the noise and the heat flow at the same lead are related by a simple relation. For instance, if the amplitude $`r`$ of a reflection coefficient is independent of $`X`$ the ratio $`S^{(1)}`$/$`I_E^{(1)}`$ at low temperature ($`k_BT\mathrm{}\omega `$) is $`8\pi G/\omega `$, where $`G=e^2t^2/h`$ is the conductance of our mesoscopic sample.

On the other hand if the phase $`\theta `$ of a reflection coefficient is independent of the varying parameter $`X`$ the contribution of quasi-particle correlations to the current correlations at the same lead vanishes, i.e., $`T_{\alpha \alpha }^{(cor)}=0`$. In this case the noise $`S_{\alpha \alpha }^{(pump)}`$ reaches the Poisson value (the Fano factor is $`F=1`$) and the ratio of the noise to the heat flow is a universal function of the temperature and the pump frequency

$$S^{(1)}/I_E^{(1)}=F^{(1)}(\mathrm{}\omega ,k_BT)=4e^2F(\mathrm{}\omega ,k_BT)/(\mathrm{}\omega )^2,$$

which is independent of individual features of a scatterer. Comparing Eq.(23) and Eq.(27) we see that this conclusion is a quite general feature of a weak amplitude pump: if $`T_{\alpha \alpha }^{(cor)}=0`$ the ratio $`S_{\alpha \alpha }^{(pump)}/I_{E,\alpha }`$ = $`F^{(1)}(\mathrm{}\omega ,k_BT)`$. Note that even if the phase $`\theta `$ is independent of $`X`$, $`T_{\alpha \beta }^{(cor)}0`$ ($`\alpha \beta `$) and the quasi-particle correlations are important for the current cross-correlations. To illustrate this fact in the next subsection we consider a particular case of a multiterminal (three-terminal) conductor. We investigate a case when the phase of the transmission (reflection) amplitude is unimportant but the current cross-correlations are present.

### B Noise and heat flow of an oscillating wave splitter

Let us consider a wave splitter in which one lead $`\alpha =1`$ couples via a tunnel barrier with transparency $`ϵ`$ symmetrically to two leads $`\alpha =2,3`$. We assume that this three lead structure is described by the single parameter scattering matrix

$$\widehat{s}=\left(\begin{array}{ccc}(a+b)& \sqrt{ϵ}& \sqrt{ϵ}\\ \sqrt{ϵ}& a& b\\ \sqrt{ϵ}& b& a\end{array}\right).$$
(39)
where $`a=(\sqrt{12ϵ}1)/2`$ and $`b=(\sqrt{12ϵ}+1)/2`$. For $`ϵ=0`$ carriers incident from lead $`1`$ are completely reflected, for $`ϵ=1/2`$ carriers incident from lead $`1`$ are transmitted (without reflection) with equal probability into leads $`2`$ and $`3`$.

We choose the transparency $`ϵ`$ as an external parameter and assume that it is subject to small amplitude oscillations $`ϵ(t)=ϵ+2ϵ_\omega \mathrm{cos}(\omega t)`$ ($`ϵ_\omega ϵ`$). With the scattering matrix Eq.(39) the parametric emissivity matrix $`\widehat{\nu }[ϵ]`$ Eq.(13) is

$$\widehat{\nu }[ϵ]=\frac{1}{4\pi \sqrt{(12ϵ)ϵ}}\left(\begin{array}{ccc}0& i& i\\ i& 0& 0\\ i& 0& 0\end{array}\right).$$
(40)
Using Eqs.(15) and (29) we can calculate $`T_{\pm \omega ,\alpha }`$ and $`T_{\alpha \beta }^{(cor)}`$. Substituting these probabilities into Eqs.(10), (23) and (27) we obtain the quantities of interest here. We see that the dc charge current $`I_\alpha `$ is zero at all leads (as it is expected because only one parameter is varied). However the heat flow $`I_{E,\alpha }`$ and zero-frequency noise power show interesting features.

Despite the fact that lead $`\alpha =1`$ is coupled only weakly, the heat flow $`I_{E,\alpha }`$ (for any $`ϵ>0`$) and the shot noise power $`S_{\alpha \alpha }^{(pump)}`$ are the same in all three leads

$$I_{E,1}=I_{E,2}=I_{E,3}=\frac{\mathrm{}\omega ^2ϵ_\omega ^2}{4\pi ϵ(12ϵ)}.$$
(41)
The heat flow is related to the noise in the simple way discussed at the end of the previous subsection $`S_{\alpha \alpha }^{(pump)}/I_{E,\alpha }`$ = $`4e^2F(\mathrm{}\omega ,k_BT)/(\mathrm{}\omega )^2`$.

The asymmetry between leads appears in cross-correlations $`S_{\alpha \beta }^{(pump)}`$ ($`\alpha \beta `$). The cross-correlation of current fluctuations at the two symmetrically coupled leads vanishes $`S_{23}^{(pump)}=0`$ but the cross-correlations invoking the fluctuating current of the weakly coupled lead $`1`$ are non-vanishing, $`S_{12}^{(pump)}`$ = $`S_{13}^{(pump)}`$ = $`\frac{1}{2}S_{\alpha \alpha }^{(pump)}`$.

### C Heat flow and noise in two parameter pumps

Now we return to the scattering matrix Eq.(30) and assume that it depends on two parameters $`X_1(t)`$ and $`X_2(t)`$. In this case we can represent the probabilities for absorption and emission in terms of a symmetric and antisymmetric contribution

$$T_{\pm \omega ,\alpha }^{(2)}=T_{\omega ,\alpha }^{(s)}\pm T_{\omega ,\alpha }^{(a)}.$$
(42)
Here the symmetric (with respect to absorption and emission of a modulation quantum $`\mathrm{}\omega `$) $`T_{\omega ,\alpha }^{(s)}`$ and antisymmetric $`T_{\omega ,\alpha }^{(a)}`$ parts which determine the heat flow Eq.(23) and the dc-current Eq.(10), respectively, are

$`T_{\omega ,\alpha }^{(s)}=T_{\omega ,\alpha }^{(1)}[X_1]+T_{\omega ,\alpha }^{(1)}[X_2]`$ (43)
$`+2X_{\omega ,1}X_{\omega ,2}\mathrm{cos}(\phi _2\phi _1)Re[\mathrm{\Pi }_\alpha (X_1,X_2)],`$ (44)
$$T_{\omega ,\alpha }^{(a)}=2X_{\omega ,1}X_{\omega ,2}\mathrm{sin}(\phi _2\phi _1)Im[\mathrm{\Pi }_\alpha (X_1,X_2)],$$
(45)
with the quantity $`\mathrm{\Pi }_\alpha (X_1,X_2)`$ being

$`\mathrm{\Pi }_\alpha (X_1,X_2)={\displaystyle \underset{\beta }{}}\left[{\displaystyle \frac{s_{\alpha \beta }^{}}{X_1}}{\displaystyle \frac{s_{\alpha \beta }}{X_2}}\right]`$ (46)
$`=4\pi ^2{\displaystyle \underset{\beta }{}}\nu _{\alpha \beta }^{}[X_1]\nu _{\alpha \beta }[X_2].`$ (47)
In Eq.(44) the quantity $`T_{\omega ,\alpha }^{(1)}[X]`$ is given by Eq.(32).

Substituting Eq.(42) into Eq.(10) we immediately reproduce the result obtained by Brouwer  for the pumped current (at small amplitudes $`X_{\omega ,j}`$)

$$I_\alpha =\frac{2e\omega }{\pi }X_{\omega ,1}X_{\omega ,2}\mathrm{sin}(\phi _2\phi _1)Im[\mathrm{\Pi }_\alpha (X_1,X_2)].$$
(48)
We see that varying only one (no matter which) parameter we can not generate a dc-current. But if at least two parameters $`X_1`$ and $`X_2`$ are varied periodically but out of phase $`\phi _1\phi _2`$ then the mesoscopic sample can continuously pump charge between reservoirs with the same chemical potentials. This is a consequence of quantum-mechanical interference effects. These effects manifest themselves not only in the dc-current but also in the heat flow and in the noise.

To characterize the contribution of interference effects to the heat flow we consider the difference $`\mathrm{\Delta }I_{E,\alpha }`$ between the heat flow $`I_E^{(2)}`$ when the two parameters are varied simultaneously and the sum of heat flows when only one parameter oscillates $`I_{E,\alpha }^{(1)}[X_1]`$ \+ $`I_{E,\alpha }^{(1)}[X_2]`$. This difference is

$$\mathrm{\Delta }I_{E,\alpha }=\frac{\mathrm{}\omega ^2}{\pi }X_{\omega ,1}X_{\omega ,2}\mathrm{cos}(\phi _2\phi _1)Re[\mathrm{\Pi }_\alpha (X_1,X_2)].$$
(49)
We see that the additional heat production $`\mathrm{\Delta }I_{E,\alpha }`$ and the dc-current give a full description of the quantity $`\mathrm{\Pi }_\alpha (X_1,X_2)`$, i.e., they determine the real and the imaginary parts.

For the scattering matrix Eq.(30) we get

$`Re[\mathrm{\Pi }_1(X_1,X_2)]=Re[\mathrm{\Pi }_2(X_1,X_2)]`$ (50)
$`=r^2{\displaystyle \frac{\theta }{X_1}}{\displaystyle \frac{\theta }{X_2}}+{\displaystyle \frac{r}{X_1}}{\displaystyle \frac{r}{X_2}}+{\displaystyle \frac{t}{X_1}}{\displaystyle \frac{t}{X_2}},`$ (51)
$`Im[\mathrm{\Pi }_1(X_1,X_2)]=Im[\mathrm{\Pi }_2(X_1,X_2)]`$ (52)
$`={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{r^2}{X_1}}{\displaystyle \frac{\theta }{X_2}}{\displaystyle \frac{\theta }{X_1}}{\displaystyle \frac{r^2}{X_2}}\right).`$ (53)
As we did for the heat production we calculate an additional noise $`\mathrm{\Delta }S_{\alpha \beta }^{(pump)}`$ generated by two simultaneously oscillating parameters $`X_1`$ and $`X_2`$ over the sum of noises produced by each of them separately,

$`\mathrm{\Delta }S_{\alpha \beta }^{(pump)}={\displaystyle \frac{8e^2}{h}}X_{\omega ,1}X_{\omega ,2}F(\mathrm{}\omega ,k_BT)`$ (54)
$`\times \mathrm{cos}(\phi _2\phi _1)Re[N_{\alpha \beta }(X_1,X_2)],`$ (55)
where

$$N_{\alpha \alpha }(X_1,X_2)=4\pi ^2\underset{\gamma \alpha }{}\nu _{\alpha \gamma }^{}[X_1]\nu _{\alpha \gamma }[X_2],$$
(56)
$$N_{\alpha \beta }(X_1,X_2)=4\pi ^2\nu _{\alpha \beta }^{}[X_1]\nu _{\alpha \beta }[X_2],\alpha \beta .$$
(57)
Note the cosine dependence of the additional heat and noise on the phase difference $`\mathrm{\Delta }\phi `$ = $`\phi _2\phi _1`$. In contrast the pumped current Eq.(48) is determined by the sine of the phase difference $`\mathrm{\Delta }\phi `$. As a consequence if the pumped current is large (as a function of $`\mathrm{\Delta }\phi `$) the additional noise and heat flow are small.

For the scattering matrix Eq.(30) we have $`Re[N_{11}]`$ = $`Re[N_{22}]`$ = -$`Re[N_{12}]`$ = -$`Re[N_{21}]`$ = $`Re[N^{(2)}]`$

$`Re[N^{(2)}(X_1,X_2)]=r^2t^2{\displaystyle \frac{\theta }{X_1}}{\displaystyle \frac{\theta }{X_2}}`$ (58)
$`+{\displaystyle \frac{r}{X_1}}{\displaystyle \frac{r}{X_2}}+{\displaystyle \frac{t}{X_1}}{\displaystyle \frac{t}{X_2}}.`$ (59)
From Eq.(51) and Eq.(59) we can see that the additional heat production and the additional noise vanish if the amplitude of the reflection coefficient $`r`$ and its phase $`\theta `$ depend on a single parameter only (not necessary the same). However there remains a heat production Eq.(36) and noise Eq.(38) owing to independently oscillating parameters. As it is evident from Eq.(32) this unavoidable heat production (see Eq.(23)) and noise (see Eq.(27)) are present always if only the ”pump” is working. On the other hand if the phase $`\theta `$ of the reflection coefficient depends on only one varying parameter the ratio of additional noise Eq.(55) to additional heat production Eq.(49) does not depend on the scattering matrix and is equal to $`4e^2F(\mathrm{}\omega ,k_BT)/(\mathrm{}\omega )^2`$.

Under some conditions the additional noise $`\mathrm{\Delta }S_{\alpha \alpha }^{(pump)}`$ can be related to the dc-current $`I_\alpha `$ at the same lead. If $`r=r(X_1)`$ (or $`r=r(X_2)`$) then their ratio

$$\frac{\mathrm{\Delta }S_{\alpha \alpha }^{(pump)}}{I_\alpha }=()4et^2\mathrm{cot}(\phi _2\phi _1)\frac{\theta }{r^2}\frac{F(\mathrm{}\omega ,k_BT)}{\mathrm{}\omega },$$

is independent of the varying parameters. On the other hand if $`\theta =\theta (X_2)`$ (or $`\theta =\theta (X_1)`$) then we get

$$\frac{\mathrm{\Delta }S_{\alpha \alpha }^{(pump)}}{I_\alpha }=()\frac{e}{t^2r^2}\mathrm{cot}(\phi _2\phi _1)\frac{r^2}{\theta }\frac{F(\mathrm{}\omega ,k_BT)}{\mathrm{}\omega }.$$

In the low temperature limit this ratio becomes independent of frequency.

## VIII Conclusion

In this work we have developed the approach to the kinetics of an adiabatic quantum pump for an arbitrary relation of pump frequency and temperature. Our consideration is based on the scattering matrix approach for ac-transport. This approach takes into account the existence of the side bands of particles exiting the pump and thus allows a description of the quantum statistical correlation properties (e.g., noise) of an adiabatic quantum pump. The side bands correspond to particles which have gained or lost a modulation quantum $`\mathrm{}\omega `$. We find that our results for the pump current, the heat flow and the noise can all be expressed in terms of a parametric emissivity matrix. In particular we find that the current cross-correlations of a multiterminal pump are directly related a to a non-diagonal element of the parametric emissivity matrix.

Using the quasi-particle picture we have given a simple physical interpretations of processes leading to charge and energy transfer in the system. Due to the oscillations of the scatterer the electron system gains energy (the side bands arise). Absorption of an energy quantum $`\mathrm{}\omega `$ leads to the creation of a nonequilibrium (quasi-)electron-hole pair. These quasi-particles carry energy from the scatterer to the reservoirs. On average the electron-hole pair is neutral thus the pump is not a source of a charge current but under some conditions  can only push charge from some reservoirs to others. These conditions can be realized if the quasi-electron and hole (belonging to the same pair) leave the scatterer through different leads (say, $`\alpha `$ and $`\beta `$). In this case these quasi-particles contribute to the charge transfer between the reservoirs $`\alpha `$ and $`\beta `$. These quasiparticles are correlated since they are created in the same event. This is also the source of the correlations between the currents at the leads $`\alpha `$ and $`\beta `$ (cross-correlations). Thus we conclude that the existence of the dc-current in a weak amplitude pump is always accompanied by current correlations (shot noise). This type of a pump can not be optimal (in particular, noiseless) in the sense of Ref..

To assess the possibility of an optimal adiabatic quantum pump further investigations are necessary. In particular large amplitude variations of the external parameter have to be considered. We hope that by taking into account many photon processes (see, e.g., ) the approach developed in the present paper can be generalized to the case of a strong amplitude adiabatic quantum pump.

ACKNOWLEDGMENTS

M. B. thanks G. M. Graf for stimulating discussions. Both of us acknowledge useful discussions with A. Alekseev. M.M. gratefully acknowledges the support of the University of Geneva where this work was done. This work was supported by the Swiss National Science Foundation.


format: text

text:
# Synchronization in large populations of limit cycle oscillators with long-range interactions

## 1 Introduction

Coupled limit-cycle oscillators have been extensively used to model the collective behavior of systems in physics , chemistry and biological sciences . The existence of a self organized periodic motion is a direct consequence of the energy exchange between the system and its surroundings keeping the system far from equilibrium. The mathematical description of those systems relies on differential equations, or its discrete versions, since no Hamiltonian exists, and the appearance of a macroscopically synchronized state has been treated as non-equilibrium phase transition . The coupling among oscillatory units is usually considered to be uniform and of infinite range or with first neighbors only. The reason for these choices is simplicity, uniform coupling makes possible analytical calculations, and nearest neighbor coupling is a computationally efficient form of short range coupling. Here we would like to consider a long range distance dependent coupling of the form $`r^\alpha `$, where $`r`$ is the distance between oscillators and $`\alpha >0`$ defines the range: $`\alpha =0`$ corresponds to uniform coupling, and $`\alpha \mathrm{}`$ to nearest neighbor coupling. For this we consider two different approaches: the Kuramoto model for the phases and a cell dynamical system (CDS) model for phases and amplitudes along the limit cycle. We observe that, as in equilibrium systems, the critical value for synchronization will depend on the system size and on $`\alpha `$. For both models we show that a scaling proposed for equilibrium nonextensive spin systems works well for the synchronization process.

## 2 Long Range Coupling in Equilibrium Spin Systems

It well known that, depending on the dimensionality $`d`$ of the system, interacting potentials of the type $`J/r^\alpha `$ may lead to nonextensive energies. This problem is usually avoided in the case of uniform or mean field coupling, for which $`\alpha =0`$, by the redefinition of the coupling term with the division by the size of the system, $`N`$. This correction may be extended to other values of $`\alpha `$ in the nonextensive regime with the definition of a $`N^{}(\alpha ,d)`$ that may be calculated from the requirement that the total energy of the system must be finite. The equilibrium properties of chains of ferromagnetically coupled spins have been considered in the nonextensive regime for the Ising model  and XY model , the finiteness condition for the energy leads to a common correction factor

$$N^{}(\alpha ,N)=_1^{N^{1/d}}𝑑rr^{d1}r^\alpha =\frac{N^{1\alpha /d}1}{1\alpha /d}.$$
(1)
The extensive behavior is recovered with the definition of the scaled interaction strength $`J^{}=J/(2^\alpha N^{})`$. For $`\alpha >d`$, (1) provides a correction for finite size effects, and for $`\alpha <d`$ it removes the dependency of the critical coupling strength on $`\alpha `$ and $`N`$. In the limit $`N\mathrm{}`$, $`N^{}=N`$ for mean-field coupling ($`\alpha =0`$) and, $`N^{}=\mathrm{ln}N`$ for $`\alpha =d`$. The behavior of thermodynamical quantities like magnetization and energy as a function of temperature $`T`$ for different values of $`N`$ and $`\alpha `$ are found to collapse with the mean field curve when plotted as a function of $`T/(N^{}2^\alpha )`$, indicating that this regime contains all the information regarding the nonextensive behavior.

## 3 Kuramoto Model

One of the simplest models for systems of coupled limit-cycle oscillators is the Kuramoto model for the time evolution of the phases. As shown in , it possible to reduce any system of $`N`$ weakly coupled limit cycle oscillators near the Hopf bifurcation, to the set of equations

$$\dot{\theta }_i=\omega _i+\underset{j}{\overset{N}{}}\mathrm{\Gamma }_{ij}(\theta _j\theta _i),(i=1,2,\mathrm{},N)$$
(2)
where $`\theta _i`$ and $`\omega _i`$ are the phase and natural frequency of the $`ith`$ oscillator respectively. The coupling term $`\mathrm{\Gamma }_{ij}(\theta _i\theta _j)`$ is a periodic function that may be obtained from the original equations of motion, regardless of the dimensionality of the oscillating field. The case of uniform, or mean field, coupling has been considered with detail by Kuramoto  for the simplest form of the coupling term,

$$\mathrm{\Gamma }_{ij}=N^1K\mathrm{sin}(\theta _j\theta _i).(K>0)$$
(3)
Kuramoto considered the possibility of a collective oscillatory motion when the natural frequencies are picked from a given distribution $`g(\omega )`$, as would happen in a real system. As the coupling modifies the frequencies, acting as an external oscillatory force, it is necessary to define the asymptotic frequency of the $`ith`$ oscillator after the transient regime as

$$\stackrel{~}{\omega }_i=\underset{T\mathrm{}}{lim}\frac{1}{T}\left[\theta _i(t+T)\theta _i(t)\right].$$
(4)
The existence of a collective behavior with phase and frequency entrainment depends on the coupling strength and on the choice of $`g(\omega )`$ . In what follows we consider only the case where $`g(\omega )`$ is Gaussian. Kuramoto characterized the degree of synchronization in the asymptotic regime by means of two order parameter like quantities: $`R`$ for the frequencies, and $`\sigma `$ for the phases. Frequency entrainment may be quantified simply by analyzing the formation of frequency clusters and defining

$$R=\underset{N\mathrm{}}{lim}\frac{N_s}{N},$$
(5)
where $`N_s`$ is the size of the largest cluster of oscillators with a common asymptotic frequency $`\stackrel{~}{\omega }`$. The complex quantity

$$\sigma =\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}(i\theta _j)$$
(6)
is appropriate for the phase ordering as it measures the concentration of phases on a certain value. $`R`$ and $`\sigma `$ behave as equilibrium order parameters in the sense that, for $`K`$ below some critical value $`K_c`$, oscillations are incoherent leading to $`R`$ and $`\sigma `$=0. For strong enough coupling, or $`K>K_c`$, collective oscillation takes place and $`R,\sigma >0`$.

For the case of uniform coupling, phase and frequency entrainment appear simultaneously at the critical value

$$K_c=\frac{2}{\pi g(\omega _0)},$$
(7)
where $`\omega _0`$ is the average value of the distribution of natural frequencies .

The long range version of the coupling term is written as

$$\mathrm{\Gamma }_{ij}=\frac{K}{r_{ij}^\alpha }\mathrm{sin}(\theta _j\theta _i),K>0(r_{ij}=1,2,3\mathrm{}),$$
(8)
and was considered in the case of oscillator chains in . In Hamiltonian systems the limit of non-extensivity can be found by requiring the upper value of the energy to be finite as $`N\mathrm{}`$. Limit cycle oscillators are open systems and energy is a meaningless quantity so, instead, we establish the finiteness condition for $`(1/N)_i\dot{\theta }_i`$ in the asymptotic regime. This condition may be established by requiring that the upper value of the coupling term divided by $`N`$ must be finite, as in the mean field coupling case. With this, we end up with a condition identical to (1), but, as will be clear from the graphs below, the factor $`2^\alpha `$ does not seem to be present in this case. For any value of $`d`$ and $`\alpha `$, it is then possible to incorporate the dependence on $`\alpha `$ and $`N`$ of the coupling term into the value of $`K`$ with the definition of a scaled coupling strength

$$K^{}(\alpha ,N)=KN^{}.$$
(9)
In order to numerically study the non-extensive region, we used two dimensional square lattices of coupled oscillators with natural frequencies sorted from the Gaussian distribution $`g(\omega )=(1/\sqrt{2\pi })\mathrm{exp}(\omega ^2/2)`$ and periodic boundary conditions. Equation (2) with the coupling (8) was then integrated using the Euler method with d$`t=0.1`$. The first 3000 iterations where discarded as a transient, and after 6000 iterations past the transient, we calculated the asymptotic frequencies as defined in (4). We used the Hoshen-Kopelman method  to analyze the frequency pattern, two sites were considered to belong to same cluster if the frequency difference between them was less than $`2\pi /T`$, where $`T`$ is the time elapsed after the transient. $`R`$ was then associated to the largest cluster, and averaged over 10 initial conditions for the natural frequencies.

The non-extensivity is evident if we consider the behavior of $`R`$ as a function of $`K`$ for $`\alpha <d`$. For fixed $`\alpha `$, the critical value of $`K_c`$ decreases with increasing $`N`$ as can be seen in Fig. (1.a), as more terms contribute to the coupling. The shift in $`K_c`$ may be corrected with the use of (1). In fact, if we plot $`R`$ as a function of $`K^{}`$, there is no dependence on the system size, the curves collapse and a unique value of $`K_c^{}`$ appears. Figure (1.b) shows the behavior of $`R`$ as a function of $`K^{}`$ for fixed $`\alpha `$.

The dependence of $`K_c`$ on $`\alpha `$ is evident if we plot $`R`$ for fixed $`N`$ and different values of $`\alpha `$. Figure (2.a) shows the curves for $`\alpha =0.25`$, 0.50, 0.75, and 0 (mean-field coupling without the correction factor $`1/N`$). It is important to notice that we did not use the factor $`2^\alpha `$ as predicted for the spin systems. Consistently with dependence on system size, as $`\alpha `$ decreases, more terms are included in the coupling, leading to smaller values of $`K_c`$. Again, (1) fixes the dependence on $`\alpha `$ providing a perfect collapse even with the mean-field curve as can be seen in Figure (2.b). Notice that the value of $`K_c^{}`$ in Figures (1) and (2) is in good agreement with (7), our frequency distribution leads to $`K_c^{}=1.59`$.

## 4 Cell Dynamical System Model

A more complete description of the system includes information about the amplitudes by means of the definition of an oscillatory field

$$z(n,t)=\rho (n,t)\mathrm{exp}[i\varphi (n,t)].$$

This is an important aspect since oscillator death is possible as we effectively increase the coupling strength with the variation of system size and coupling range .

The spirit of cell dynamical modeling is the discretization of time and space from the start, as opposed to formulating a continuous model and then further adopting some discretization scheme. In this sense, $`(n,t)`$ defines the position of the $`n`$-th cell at time $`t`$, $`\rho `$ and $`\varphi `$ are then the amplitude and the phase of the field at given cell. The time evolution of the system is then given by the repeated application of the map

$$z(i,t+1)=F[z(i,t)]+\underset{ji}{}D_{ij}\{F[z(j,t)]F[z(i,t)]\}$$
(10)
where the local dynamics is regulated by $`F[z(t);A]=A\mathrm{exp}(i\omega )\frac{z(t)}{\sqrt{1+|z(t)|^2(A^21)}}`$. This form of $`F`$ corresponds to an asymptotically stable limit cycle of unit radius for all $`A>1`$ . Here the parameter $`A`$ is a measure of the attraction to the limit cycle. The long range coupling is implemented through the definition of $`D_{ij}=D_0/r_{ij}^\alpha `$.

For observation of synchronization, we considered chains of $`N=100`$ oscillators with periodic boundary conditions, $`A=2.0`$, Gaussian distribution of natural frequencies with unit mean value and variance 0.1. We iterated (10) 4000 times, the first 1000 being discarded as a transient, and the order parameter $`\sigma `$ was then calculated. This procedure was repeated for different chain sizes and values of the decay exponent $`\alpha `$ for several values of the coupling strength $`D_0`$. Figure 3(a) shows the synchronization curves averaged over 20 samples. As expected, the critical value of $`D_0`$ is dependent on $`N`$ and $`\alpha `$. Although the calculation of the upper bound of the coupling term is rather different from the one in the Kuramoto model or spin systems, the same form of scaling is observed here. Figure 3(b) shows the synchronization curves of $`\sigma `$ as a function of $`D_0^{}=D_02^\alpha N^{}`$.

It is important to verify the effect of nonextensivity on the stability of the oscillatory state. It is well known that, depending on the amplitudes and frequency spread, an increase in coupling strength may lead to stabilization of the rest state. This has been observed for systems with mean field coupling  and for pairs of oscillators . In order to observe this effect we increased the frequency spread to 0.7, then, for a fixed value of $`N`$ and $`\alpha `$, we varied $`D_0`$ and $`A`$. After the 4000 iterations we verified whether the amplitude was non zero. The same procedure was repeated for $`\alpha =0.1`$, 0.5 and 0.9. Figure 4(a) shows the $`AD_0`$ phase diagram indicating regions of synchronization, incoherent motion and death for the different values of $`\alpha `$. The scaled phase diagram can be seen in Fig. 4(b).

## 5 Conclusion

The above results bring the process of synchronization even closer to an order disorder transition, as we directly use the scaling proposed for nonextensive equilibrium systems. Also, we see that the important quantity for synchronization and stability of the oscillatory phase is the scaled coupling constant. With this in mind, one can easily tune the interaction strength to obtain the desired regime, in the case of distant dependent coupling. Some points still remain unclear. First, there is no a priori reason for the scaling of the CDS model to be the same as that of Hamiltonian spin systems. Second we cannot justify the absence of the factor $`2^\alpha `$ in the scaling of the Kuramoto model, which has a condition for extensivity identical to Hamiltonian spin systems. Actually we found by tentative and error that for the Kuramoto oscillators a better form of scaling should have a factor $`2^{\alpha /d}`$. Since the results for spin Hamiltonians and for the CDS model are for one dimensional systems, this correction would not affect their scaling, although we are not able to analytically justify this form.

## acknowledgments

The authors thank Francisco Tamarit and Celia Anteneodo for interesting discussions and the Brazilian agencies CNPq, FAPERJ and FUJB-UFRJ for financial support.


text:
# 𝐽/𝜓 production in relativistic heavy ion collisions from a multi-phase transport model

## I Introduction

A central interest in modern nuclear physics is to produce and study the Quark-Gluon Plasma (QGP) predicted by the Quantum Chromodynamics (QCD) . Studying the properties of the QGP and the deconfinement phase transition not only is important for understanding the QCD but also has astronomical implications . Experiments at the Super Proton Synchrotron (SPS) involving collisions between heavy nuclei have shown that a large amount of energy is deposited around mid-rapidity, and the energy density may already be sufficiently large for the formation of the QGP during the initial stage of collisions . More recent experiments at RHIC at Brookhaven National Laboratory, where the collision energy is much higher than at SPS, have also shown possible effects due to the formation of a partonic matter .

To find the signal for the Quark-Gluon Plasma in relativistic heavy ion collisions, Matsui and Satz proposed to study $`J/\psi `$ production in these collisions . Based on the fact that the effective potential between charm and anti-charm quarks changes in the QGP due to color screening effect, no bound states can be formed between them beyond a certain critical temperature, which is somewhat higher than the deconfinement temperature. As a result, $`J/\psi `$ production is expected to be suppressed if the QGP is formed in heavy ion collisions. There are other possible mechanisms for $`J/\psi `$ suppression  in heavy ion collisions, e.g., $`J/\psi `$ can be destroyed by collisions with incoming nucleons or with gluons in the initial partonic matter  or with comoving hadrons in the hadronic matter . The observed abnormal suppression of $`J/\psi `$ in central Pb+Pb collisions at SPS may require the formation of the Quark-Gluon Plasma in these collisions .

For heavy ion collisions at RHIC energies, unlike in previous fixed target experiments at SPS, multiple pairs of charm-anti-charm quarks can be produced in one event, and these charm and anti-charm quarks can interact and produce the $`J/\psi `$ . This mechanism is the inverse reaction of $`J/\psi `$ absorption by gluons. Using a kinetic approach, it was found that despite $`J/\psi `$ dissociation due to color screening, this reaction is important in the quark-gluon plasma with temperatures between those for $`J/\psi `$ dissociation and the deconfinement, and would lead to an enhanced production of $`J/\psi `$ in heavy ion collisions at RHIC. A possible signal for QGP at RHIC has thus been changed from $`J/\psi `$ suppression to $`J/\psi `$ enhancement. The statistical fragmentation model , which has been extensively used for studying particle production in heavy ion collisions at SPS energies, predicts a comparable $`J/\psi `$ production. The main assumption of this model is that $`J/\psi `$ formed during hadronization of the QGP is in chemical equilibrium with charm mesons.

In the kinetic formation model, the heavy ion collision dynamics is treated schematically. For a more quantitative study of this new $`J/\psi `$ production mechanism, we study $`J/\psi `$ production from both the charm-anti-charm quark interactions and the charm-anti-charm meson interactions using the AMPT model . In particular, we consider central (b=0) Au+Au collisions at the highest RHIC energy $`\sqrt{s}=200`$ AGeV. With this dynamic transport model, we find that the net number of produced $`J/\psi `$ from the parton and hadron phases is smaller than that expected from initial nucleon-nucleon collisions, contrary to the $`J/\psi `$ enhancement predicted by the kinetic formation model. The production of $`J/\psi `$ is further suppressed if one includes the color screening effect in the parton phase.

In Section II, we give a short description of the multi-phase transport model (AMPT). Results are presented in Section III for $`J/\psi `$ production from the initial parton phase and the final hadron phase. In this section, we also discuss the dependence of $`J/\psi `$ production in the parton phase on the charm quark mass and in the hadron phase on the charm meson mass. Finally, we summarize the paper in Section IV.

## II Multi-Phase Transport model and $`J/\psi `$ production mechanisms

The multi-phase transport model  that we use for the present study is a hybrid model based on three Monte-Carlo models for the three stages of relativistic heavy ion collisions, i.e., the Heavy Ion Jet Interaction Generator (HIJING)  for the initial conditions, Zhang’s Parton Cascade (ZPC)  for the parton evolution, and A Relativistic Transport (ART) model  for the hadron evolution. In the HIJING model, nucleons are distributed according to the Woods-Saxon distribution. The parton distribution in a nucleon is adjusted according to the transverse position of the nucleon and a shadowing function based on the Mueller-Qiu parton recombination mechanism . These partons undergo hard or semi-hard collisions and generate an initial parton system after the passage of two colliding nuclei. HIJING provides the momentum space information of the partons produced in the hard or semi-hard collisions. The formation time of a parton is generated according to the Gyulassy-Wang model of parton coherent production . The position of a physical parton is then obtained from the position of its parent nucleon by adding a displacement given by its velocity times the formation time. The resulting parton phase space distribution is used as the initial conditions for starting the ZPC parton evolution. In the present version of the ZPC model, only elastic gluon-gluon scattering are included with a cross section of 3 mb, which is consistent with the estimated average screening mass in the partonic matter formed in the collisions. When the partonic matter stops interacting, partons are connected back with their parent nucleons, which also suffer soft interactions through string excitations. The Lund fragmentation model  is then called from the HIJING model to convert these wounded nucleons into hadrons. These hadrons become physical after a formation time given by the average freeze-out proper time of partons in a string plus an additional 0.7 fm/c to account for the time that is required for producing quark-anti-quark pairs. The ensuing hadron evolution is modeled by the ART model that takes into account interactions among nucleons, anti-nucleons, mesons, and resonances. For a better description of experimental hadron spectra, the original ART model has been improved by including the $`K^{}`$ meson interactions, the anti-nucleon production, and the $`\rho `$ meson interactions.

We have previously used this multi-phase transport model to study the suppression of primarily produced $`J/\psi `$ in heavy ion collisions at RHIC . These $`J/\psi `$’s are produced in nucleon-nucleon collisions during the initial passage of two colliding nuclei. Because of color screening, these $`J/\psi `$’s are not physical ones but correlated pairs of charm-anti-charm quarks close in phase space. They become physical $`J/\psi `$’s if they remain close when the partonic matter cools below the $`J/\psi `$ dissociation temperature. Otherwise, they will combine with nearby light quarks to form charm mesons at hadronization. Gluons can also destroy the phase space correlation between the charm quark and anti-quark pairs and hence destroy the $`J/\psi `$. It was found that finite space size and lifetime have significant effects on $`J/\psi `$ suppression, and different mechanisms lead to different suppression factors. Also, the dependence of $`J/\psi `$ suppression on the size of colliding nuclei was shown to be useful for differentiating the different suppression mechanisms. Similar results have also been obtained in other studies .

In the present study, we include the reactions of $`J/\psi `$ production from charm-anti-charm quark interactions in the parton phase and from charm meson interactions in the hadron phase. As the primary $`J/\psi `$ number is much smaller than the number of other particles participating in the equilibration, primary $`J/\psi `$ evolution won’t affect secondary $`J/\psi `$ production. We do not include primary $`J/\psi `$’s in the following calculations.

The charm and anti-charm quarks are generated from the PYTHIA model , with a cross section of 350 $`\mu `$b for producing a pair of charm-anti-charm pair in a nucleon-nucleon collision. This cross section does not include nuclear shadowing effect of gluons, so the present study gives an upper bound for charm meson and $`J/\psi `$ production in heavy ion collisions. In the partonic stage of heavy ion collisions, $`J/\psi `$’s are then produced from charm-anti-charm quark interactions through the reaction $`c+\overline{c}J/\psi +g`$. The $`J/\psi `$ can also be destroyed by gluons through the reaction $`g+J/\psi c+\overline{c}`$. We further include the dissociation of $`J/\psi `$ due to color screening in the partonic matter. To treat this process, we follow the method used in Ref.  by introducing a critical radius in the transverse plane of a heavy ion collision, within which the parton energy density is larger than the critical energy density for $`J/\psi `$ dissociation by color screening and the $`J/\psi `$ thus cannot be formed. The critical radius decreases as the partonic matter expands. The time evolution of this critical radius in Au+Au collisions at $`\sqrt{s}=200`$ AGeV has been parameterized in Ref. . In addition to producing the $`J/\psi `$, charm and anti-charm quarks also undergo elastic scattering with gluons. This elastic cross section is taken to be 3 mb in the present study.

After partons stop interacting, charm and anti-charm quarks are converted to $`D`$ and $`\overline{D}`$ mesons. This conversion is carried out using a delta function fragmentation scheme. In this approach, the momentum of a charm meson is taken to be the same as the momentum of its parent charm quark. A formation time of 1 fm/c is also introduced in combining the charm quark with a nearby light quark. Because of the large mass of charm quark, the delta function fragmentation scheme gives a reasonable description of charm quark fragmentation in hadron reactions. For simplicity, we do not differentiate between $`D`$ mesons and $`D^{}`$ mesons. Instead, we vary the $`D`$ meson mass to study the influence of different production thresholds for $`J/\psi `$ production from charmed hadrons.

In the hadronic phase, $`D`$ and $`\overline{D}`$ mesons can interact to produce a $`J/\psi `$ and a light meson through the reaction $`D\overline{D}J/\psi M`$, where $`M`$ denotes a light meson. The $`J/\psi `$ can also be destroyed by light mesons through the inverse reaction $`J/\psi MD\overline{D}`$. We only include reactions that involve $`\pi `$ mesons or $`\rho `$ mesons as they are the most abundant ones in the nearly baryon free matter formed in heavy ion collisions at RHIC. The cross section for light meson destruction of $`J/\psi `$ is taken to be 3 mb as in previous transport study at SPS energies . This cross section is consistent with those from recent theoretical studies based on effective hadronic Lagrangians , the quark-exchange model , and the QCD sum rules . The cross section for the inverse reaction of $`J/\psi `$ production from $`D`$ mesons is calculated according to the detailed balance relation. As we do not differentiate between $`D`$ mesons and $`D^{}`$ mesons, the spin and isospin degrees of freedom for both $`D`$ and $`D^{}`$ mesons are used in evaluating the inverse reaction cross section. This is equivalent to including all possible combination of $`D`$ and $`D^{}`$ states. We will vary the $`D`$ meson mass to study the effect due to changing $`J/\psi `$ production threshold.

In earlier studies based on kinetic models , $`J/\psi `$ production from the hadron phase of relativistic heavy ion collisions was found to be important only at LHC  but not at RHIC as the cross section for $`D\overline{D}J/\psi M`$ used in these studies was smaller than that used in present study.

Although multiple pairs of charm quarks are expected to be produced in each Au+Au collision, the number of charm quarks is small compared with other particles. For example, the rapidity density of charm quark pairs is about 1.5 at $`\sqrt{s}=200`$ AGeV. The expected number of $`J/\psi `$ per event is even much smaller. To obtain sufficient statistics for $`J/\psi `$ production, charm particles and the $`J/\psi `$ are treated in the transport model by the perturbative method , i.e., neglecting their effects on heavy ion collision dynamics. This method has been used extensively for studying rare particle production in heavy ion collisions at low energies. Explicitly, we include many charm quark events from primary nucleon-nucleon interactions in the evolution of partons from a single HIJING event. However, only charm quarks from the same event are allowed to interact among themselves, and the usual method of treating binary collisions in transport models is used in treating these collisions. For collisions between charm and uncharm particles, we keep their effects on the charm particles but neglect those on the uncharm particles. In this way, the dominant computation time for the interactions between uncharm particles remains the same as for a single charm event. Since there are much more uncharm particles than charm particles, e.g., about 300 gluons verse about 3 charm and anti-charm quarks per unit rapidity, neglecting the effects of charm particles on the heavy ion collision dynamics is expected to be a very good approximation.

## III Results

### A $`J/\psi `$ production in parton phase

We first study $`J/\psi `$ production in the parton phase. Unless otherwise indicated, the rates and numbers in the following are all rapidity densities averaged with $`|y|<1`$. Fig. 1 gives the time evolution of the production and destruction rates per unit rapidity for $`J/\psi `$’s with $`|y_{J/\psi }|<1`$. Also shown is the rate per unit rapidity for collisions between charm and anti-charm quarks which would have produced the $`J/\psi `$ in the absence of color screening effect. These unsuccessful collisions mainly occur before 0.5 fm/c when most partonic matter have effective temperatures above the $`J/\psi `$ dissociation temperature. The screening effect largely diminishes after 1 fm/c when the parton density decreases due to expansion. At this time, the $`J/\psi `$ production rate starts to increase quickly, and this is followed by an increase of the $`J/\psi `$ destruction rate. After reaching their peak values, both the production and destruction rates decrease and become comparable at about 4 fm/c.

To see more clearly the relative importance of $`J/\psi `$ production and destruction in the parton phase, we show in Fig. 2 the time evolution of the number of produced $`J/\psi `$ and the number of destructed $`J/\psi `$ as well as the net $`J/\psi `$ number. The number of $`J/\psi `$ saturates at about 3 fm/c. It decrease slowly afterwards as the production and destruction rates become comparable as shown in Fig. 1.

In the absence of color screening effect, collisions between charm and anti-charm quarks in the partonic matter that is above the $`J/\psi `$ dissociation temperature can also produce the $`J/\psi `$. This leads to a very high $`J/\psi `$ production rate before 0.5 fm/c as shown in Fig. 3. The $`J/\psi `$ destruction rate during this time also reaches its peak value of about 1/3 of the maximum production rate. In contrast to the case including color screening effect, the production and destruction rates become comparable already at about 0.5 fm/c. As shown in Fig. 4, these result in a much faster increase in the numbers of produced and destructed $`J/\psi `$, and an earlier peaking in the net $`J/\psi `$ number than in the case with color screening effect. Since the destruction rate after 0.5 fm/c is still slightly larger than the production rate, the net $`J/\psi `$ number slowly decreases from its maximum value of about 0.003 to a final value of about 0.0014 per event. This number is about a factor of two larger than the final net $`J/\psi `$ number for the case with color screening effect, which is about 0.0007 per event.

The number of $`J/\psi `$ expected from primary nucleon-nucleon collisions in Au+Au central collisions can be estimated from the $`J/\psi `$ production cross section in $`pp`$ collisions. The cross section per unit rapidity at $`\sqrt{s}=200`$ AGeV is 0.63 $`\mu `$b . Using the nuclear overlap function $`T_{Au+Au}(b=0)=30`$ mb<sup>-1</sup>, we expect 0.019 $`J/\psi `$ per unit rapidity to be produced in Au+Au central collisions. This is much larger than the predictions of 0.0007 and 0.0014 from the dynamical multi-phase transport model. Our study thus predicts a suppression of $`J/\psi `$ production in the parton phase of heavy ion collisions at RHIC, and this is opposite to the enhancement predicted by the kinetic formation model.

### B $`J/\psi `$ production in hadron phase

In hadronic matter, the threshold for $`J/\psi `$ production is reduced in comparison with that in the partonic matter as the $`D`$ meson mass is larger than the charm quark mass. This threshold effect is expected to compensate for the lower charm density in the hadron phase than in the parton phase. To demonstrate this threshold effect, we vary the $`D`$ meson mass in studying $`J/\psi `$ production in the hadronic matter. Fig. 5 shows the time evolution of $`J/\psi `$ production and destruction rates per unit rapidity in the hadron phase for different $`D`$ meson masses. These results are obtained with the color screening effect included in the initial parton phase. The upper panel is for $`m_D=1.70`$ GeV, which is smaller than the $`D`$ meson mass in vacuum. We use this smaller mass because theoretical studies have shown that the $`D`$ meson mass may be reduced in medium . The middle panel is for $`m_D=1.87`$ GeV, which is the $`D`$ meson mass in vacuum. Here, we treat $`D^{}`$ meson mass as $`D`$ meson mass, and the results are thus a lower bound for $`J/\psi `$ production. The lower panel is for $`m_D=2.01`$ GeV, which is the $`D^{}`$ mass in vacuum. Treating $`D`$ meson mass as $`D^{}`$ meson mass in this case then gives an upper bound for $`J/\psi `$ production. The destruction rates for all three $`D`$ meson masses are seen to have similar magnitude. For the production rate, it is larger for larger $`D`$ meson masses.

In Fig. 6, we show the time evolution of the numbers of produced and destructed $`J/\psi `$ together with that of the net $`J/\psi `$ number. For $`m_D=1.87`$ GeV and $`m_D=2.01`$ GeV, the number of produced $`J/\psi `$ is always greater than the number of destructed $`J/\psi `$, leading to a net production of $`J/\psi `$ from the hadron phase of relativistic heavy ion collisions. With increasing $`D`$ meson mass, the destructed $`J/\psi `$ number increases slightly while the produced $`J/\psi `$ increases drastically. For $`m_D=1.70`$ GeV, the produced and destructed $`J/\psi `$ numbers are about equal, and the final net $`J/\psi `$ number is thus about the same as that from the parton phase. On the other hand, for $`m_D=2.01`$ GeV, the final $`J/\psi `$ number equals approximately the number of produced $`J/\psi `$ from the hadron phase. The final $`J/\psi `$ number reflects essentially the one produced from the hadron phase.

Quantitatively, the final net number of $`J/\psi `$ per unit rapidity averaged over $`|y_{J/\psi }|<1`$ is about 0.0007 per event for $`m_D=1.70`$ GeV and increases to 0.0019 for $`m_D=1.87`$ GeV and to 0.0053 for $`m_D=2.01`$ GeV. The time evolution of the net $`J/\psi `$ number in the three cases are similar as it is largely determined by the heavy ion collision dynamics. The $`J/\psi `$ number saturates at about 15-20 fm/c, which is about the lifetime of the hadronic matter.

The results for neglecting the color screening effect in the initial parton phase are shown in Fig. 7 for time evolution of the $`J/\psi `$ production and destruction rates and in Fig. 8 for the time evolution of the numbers of produced and destructed $`J/\psi `$ as well as the net $`J/\psi `$ number. Although there are more $`J/\psi `$’s destructed in the hadronic matter, particularly for $`m_D=1.70`$ GeV where the number of destructed $`J/\psi `$ is even larger than the number of produced $`J/\psi `$ in the hadron phase, the final $`J/\psi `$ number is larger than in the case with color screening in the parton phase. Their numbers are 0.0010, 0.0024, and 0.0057 for $`m_D=1.70`$, 1.87, and 2.01 GeV, respectively.

### C discussions

The relative importance of $`J/\psi `$ production from the parton and hadron phases is more clearly seen in Fig. 9, where we show the number of produced $`J/\psi `$ per unit rapidity with $`|y_{J/\psi }|<1`$ as a function of the $`D`$ meson mass. The results with color screening in the parton phase are shown by open circles for production from the parton phase and filled circles for the final $`J/\psi `$ number including production from the hadron phase. Without color screening in the parton phase, the corresponding results are given by diamonds. As expected, including color screening effect leads to a decrease in the final $`J/\psi `$ number. As the $`D`$ meson mass increases, the final $`J/\psi `$ yield increases strongly due to production from the hadron phase, and the effect due to color screening becomes less important. The results using $`m_D=2.01`$ GeV give the upper bound for $`J/\psi `$ production while those using $`m_D=1.87`$ GeV give a lower bound. We note again that the net effect of a reduced $`D`$ meson mass of $`m_D=1.70`$ GeV is the destruction of $`J/\psi `$’s produced from the parton phase in the absence of color screening.

Part of increased production of $`J/\psi `$ in midrapidity when the effective $`D`$ meson mass increases is due to increase in the charm meson rapidity density during hadronization. This is shown in Fig. 10, where we give the midrapidity densities of charm quarks before (lower dashed line) and after (lower solid line) parton evolution as well as of charm mesons before (upper dashed line) and after (upper solid line) hadron evolution. During parton evolution, the rapidity density of charm and anti-charm quarks decreases slightly due to collisions with gluons. Hadronization increases the midrapidity density of resulting charm mesons, and its effect increases with increasing charm meson effective mass. The midrapidity density of $`D`$ and $`\overline{D}`$ mesons is seen to increase by 8% when $`m_D`$ increases from 1.70 GeV to 2.01 GeV. However, this would lead to at most 20% increase in the $`J/\psi `$ midrapidity density and cannot account for the more than factor of 5 increase in $`J/\psi `$ production shown in Fig. 9. The dominant effect for such a larger increase is thus due to the reduced threshold for $`J/\psi `$ production when $`D`$ meson mass increases.

The charm quark mass strongly influences the number of $`J/\psi `$ produced from the parton phase. Fig. 11 shows the $`J/\psi `$ rapidity density as a function of the charm quark mass. To get the upper limit of $`J/\psi `$ production, the $`D`$ meson mass is set to be 2.01 GeV and color screening is not included. The $`J/\psi `$ number increases by a factor of 72% when the charm quark mass increases from the value 1.35 GeV used in PYTHIA to 1.5 GeV used in the constituent quark model. If the charm quark mass is taken to be 1.8 GeV, we observe about a factor of 4 increase in $`J/\psi `$ production from the parton phase compared with the $`m_c=1.35`$ GeV case.

In the kinetic model approach of Ref. , the midrapidity density of charm quarks in Ref.  is about 2.5 because the largest rapidity range for about 10 pairs of charm and anti-charm quarks is 4 in that study. This is larger than that shown in Fig. 10, which is about 3.48/2=1.74. Since $`J/\psi `$ production is roughly proportional to the square of the number of charm quarks, this can lead to about a factor of 2 difference between our results and those from Ref. . On the other hand, the cross section we use for $`J/\psi `$ production from the charm and anti-charm quarks is about a factor two larger than the peak value used in the kinetic model. The effects due to different midrapidity charm quark densities and $`J/\psi `$ production cross sections in these two studies are thus largely cancelled out. The different conclusions of $`J/\psi `$ suppression in our study and $`J/\psi `$ enhancement in the kinetic approach are results of different heavy ion collision dynamics used in the two studies. In our approach, the parton phase is generated from the mini-jet gluons from initial hard and semi-hard collisions between nucleons and is evolved using the transport model. In the kinetic approach, an equilibrated quark-gluon plasma is assumed to be formed in the initial stage and later expands according to the 1-D Bjorken hydrodynamical model . Further studies are required to determine the validity of these predictions.

## IV Summary

In this paper, we have studied the effects due to interactions between charm and anti-charm quarks in the initial parton phase and between charm mesons in the final hadron phase on $`J/\psi `$ production in central Au+Au collisions at RHIC. Using a multi-phase transport model, we find that because of the high density of charm quarks at early times, the $`J/\psi `$ can be produced from the interactions of charm quarks. This is consistent with the conclusion of the kinetic formation model. However, the more realistic space-time evolution given by the transport model leads to an overall suppression of $`J/\psi `$ particles in the parton phase, instead of an enhancement predicted by the kinetic model.

We have also shown that $`J/\psi `$ can be produced in the hadron phase at RHIC. This is due to the large $`D`$ meson mass compared with charm quark mass. Although we have not taken into account parton-hadron interactions in the multi-phase model, both parton phase and hadron phase can exist in the same space-time region. This allows the production of $`J/\psi `$’s in the mixed phase. At RHIC energies, the multi-phase model shows that hadron medium effect may lead to a net destruction of the $`J/\psi `$ produced in the parton phase. In addition, color screening is shown to have significant effect on the final $`J/\psi `$ yield.

In the present study, we have not included the effect due to gluon shadowing , which would reduce the charm quark and $`J/\psi `$ numbers. As the multi-phase transport model predicts an overall suppression of the directly produced $`J/\psi `$, including gluon shadowing effect on the initial charm quark distribution will not change qualitatively the conclusion from the multi-phase model. Also, $`J/\psi `$ from decays of $`\psi ^{}`$ or $`\chi _c`$ was not considered in our study. The contribution from primary $`\psi ^{}`$ or $`\chi _c`$ to final $`J/\psi `$ are negligible due to their large destruction cross sections compared with that for $`J/\psi `$. Although the production cross sections for $`\psi ^{}`$ and $`\chi _c`$ from charm quark interactions and charm meson interactions are expected to be larger than those for $`J/\psi `$, the thresholds for these reactions are also higher. In particular, the statistical fragmentation model  indicates that the feeddown should be small for central collisions. We thus expect that our conclusion of $`J/\psi `$ suppression in heavy ion collisions at RHIC will not be qualitatively modified.

###### Acknowledgements.

We thank P. Braun-Munzinger, M. Gyulassy, J. Kapusta, S. Pratt, R. Rapp, A.T. Sustich, C. Teal, R. Thews, W. Weise, and I. Zahed for helpful discussions. We also thank the Parallel Distributed System Facility at the National Energy Research Scientific Computer Center for providing computer resources. This work is supported by the U.S. National Science Foundation under Grant No. 0088934 and Grant No. 0098805, the Arkansas Science and Technology Authority under Grant No. 01-B-20, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.