text:
# Nonlinear Quantum Transport and Current Noise

## Abstract

We study nonequilibrium transport and noise in a generic dissipative tight-binding model. Within a real-time path integral approach, we derive formally exact series expressions in the number of tunneling events for the noise valid for arbitrary bias, frequency and temperature. At zero temperature, the low-frequency noise can be summed in analytic form. The resulting Shiba-like $`|\omega |`$-singularity is a consequence of the $`1/t^2`$ decay law for the current correlation function. At finite temperature, this singularity is smoothed out.

Low-energy excitations of 1D conduction electrons and of edge currents in fractional quantum Hall (FQH) systems  are bosonic and constitute prime realizations of a Luttinger liquid (LL) phase of interacting fermions . A sensitive probe of this collective state is the tunneling conductance of electrons through an impurity or point contact in a 1D quantum wire , and of Laughlin quasiparticles through a narrow constriction for FQH edge currents . In these systems, recent attention has been paid to current fluctuations with emphasis on dc noise. To second order in the tunneling amplitude, the dc noise $`S`$ at $`T=0`$ is classical shot noise, $`S=2e^{}I`$ . Thus upon measuring dc current and dc noise in the weak-tunneling limit, one can determine the charge $`e^{}`$ of the tunneling entity. In this way, the fractional charge of the Laughlin quasiparticles has been manifested . In the strong tunneling regime, higher order tunneling processes cause a deviation of the dc noise from shot noise . Furthermore, correlations between tunneling events lead to colored noise. Singular behavior $`S(\omega )|\omega |`$ has been found within a nonperturbative analysis of the equivalent integrable Boundary Sine-Gordon (BSG) model . In this treatise, the prefactor of the $`|\omega |`$-law depends on the particular treatment of the reflection matrix for quasiparticle scattering at the boundary coming from the Fermi sea (cf. the debate in Refs. ).

In this paper, we study the low-frequency noise upon employing real-time path integral techniques to a generic quantum transport model introduced by Schmid ,

$`H`$ $`=H_{\mathrm{TB}}+{\displaystyle \underset{\alpha }{}}\left[{\displaystyle \frac{p_\alpha ^2}{2m_\alpha }}+{\displaystyle \frac{1}{2}}m_\alpha \omega _\alpha ^2\left(x_\alpha {\displaystyle \frac{c_\alpha }{m_\alpha \omega _\alpha ^2}}q\right)^2\right].`$ (1)
The Hamiltonian (1) describes a particle which moves in the tilted tight-binding (TB) lattice $`H_{\mathrm{TB}}=\frac{1}{2}\mathrm{}\mathrm{\Delta }_n(c_n^{}c_{n+1}+c_{n+1}^{}c_n)\mathrm{}ϵq/q_0`$ with lattice spacing $`q_0`$, tunneling matrix element $`\mathrm{\Delta }`$, and biasing force $`F=\mathrm{}ϵ/q_0`$. The particle is bilinearly coupled via its position $`q=q_0_nc_n^{}c_n`$ to a bath of delocalized harmonic bosons. In the reduced density matrix description for the Brownian quantum particle, all effects of the environmental coupling are captured by the spectral density $`J(\omega )=(q_0^2/2\mathrm{})_\alpha (c_\alpha ^2/m_\alpha \omega _\alpha ^2)\delta (\omega \omega _\alpha )`$. The TB model (1) is dual to a corresponding weak-binding (WB) model, in which the particle moves in a tilted cosine potential. In the duality, the power $`s`$ in the generic form $`J(\omega )\omega ^s`$ is transformed into $`2s`$. Thus, super-Ohmic friction ($`2>s>1`$) is mapped on sub-Ohmic friction ($`0<s<1`$), and vice versa . In the Ohmic case, $`J(\omega )=2K\omega `$, the TB and WB models are selfdual to each other. The dimensionless Kondo parameter $`K`$ is mapped onto the reciprocal value $`1/K`$, and the bias $`ϵ`$ onto $`ϵ/K`$, as first observed by Schmid .

The TB and WB model describe many transport problems in solid state physics. Regarding applications of the model (1) to the beforementioned case of tunneling of a charge in a LL phase, the spectral density is strictly Ohmic for $`\omega \omega _c`$. In the scaling limit, $`\omega _c`$ is by far the largest frequency of the problem. The interaction parameter $`g`$ of the Luttinger model corresponds to $`1/K`$, and the bias in the TB model (1) is related to the voltage drop $`V`$ across the weak link or constriction by $`\mathrm{}ϵ=eV`$, where $`e`$ is the unit charge. Another important example where the model (1) applies is a voltage-biased Josephson junction . Here, the WB and TB limits correspond to the charge and phase representation, respectively.

Fluctuations are measured relative to a stationary state. Therefore, we proceed as follows. We first presuppose that the global system is in a factorized system-bath state in the infinite past. Then we switch on the bath coupling and the tilting force, and let the system evolve under reign of the full Hamiltonian (1). By the time $`t=0`$, the system has reached a stationary nonequilibrium state, indicated by the index “ne”. The fluctuations of the current $`j(t)=e\dot{q}(t)/q_0`$ are reflected in the noise spectrum, which is the Fourier transform of the symmetrized current correlation function,

$$S(\omega )=_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega t}j(t)j(0)+j(0)j(t)_{\mathrm{ne}}.$$
(2)
Observing that the current correlation function is related to the mean square displacement $`C_{\mathrm{ne}}(t)=[q(t)q(0)]^2_{\mathrm{ne}}`$ by

$`j(t)j(0)+j(0)j(t)_{\mathrm{ne}}={\displaystyle \frac{e^2}{q_0^2}}{\displaystyle \frac{^2}{t^2}}C_{\mathrm{ne}}(t),`$ (3)
we can express the noise $`S(\omega )`$ in terms of the spectral function of $`C_{\mathrm{ne}}(t)`$ as

$`S(\omega )`$ $`={\displaystyle \frac{e^2}{q_0^2}}\omega ^2\stackrel{~}{C}_{\mathrm{ne}}(\omega ),`$ $`\stackrel{~}{C}_{\mathrm{ne}}(\omega )`$ $`={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{i\omega t}C_{\mathrm{ne}}(t).`$ (4)
In particular parameter regimes, the noise can directly be expressed in terms of transport quantities, e. g. the shot noise $`S=2eVG(V)`$ and Johnson-Nyquist noise $`S=4k_\mathrm{B}TG(V=0)`$, where $`G(V)=I(V)/V`$ is the nonlinear conductance. For later convenience, we also introduce the differential conductance $`G_\mathrm{d}(V)=I(V)/V`$. The conductances are related to the mobilities of the model (1), $`\mu (F)=\dot{q}(t\mathrm{})_{\mathrm{fc}}/F`$ and $`\mu _\mathrm{d}(F)=\dot{q}(t\mathrm{})_{\mathrm{fc}}/F`$, by $`G=e^2\mu /q_0^2`$ and $`G_\mathrm{d}=e^2\mu _\mathrm{d}/q_0^2`$. The index “fc” indicates that we may choose for the system-bath state at initial time zero a factorized form. This is possible since the leading contribution to mean values at asymptotic times $`t\mathrm{}`$ does not depend on system-bath correlations in the initial state.

Within nonequilibrium linear response theory, the nonequilibrium susceptibility

$`\chi _{\mathrm{ne}}(t)`$ $`={\displaystyle \frac{i}{\mathrm{}}}\mathrm{\Theta }(t)q(t)q(0)q(0)q(t)_{\mathrm{ne}},`$ $`\stackrel{~}{\chi }_{\mathrm{ne}}(\omega )`$ $`={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{i\omega t}\chi _{\mathrm{ne}}(t)`$ (5)
comprises the differential mobility as

$`\mu _\mathrm{d}(F)`$ $`=\underset{\omega 0}{lim}\left[i\omega \stackrel{~}{\chi }_{\mathrm{ne}}(\omega )\right]=\underset{t\mathrm{}}{lim}\chi _{\mathrm{ne}}(t).`$ (6)
Formally exact series expressions in $`\mathrm{\Delta }`$ for the above quantities can be derived upon employing a real-time path integral method for the reduced density matrix (RDM). Because of the discreteness of the TB model, a simple picture of a grandcanonical gas of charges can be developed. The charges are stringed at the flip times $`\{t_j\}`$, and the effects of the environment are in the charge interactions. For a path with $`n`$ flips in the negative-time branch ($`t_1,t_2,\mathrm{},t_n<0`$) and $`m`$ flips in the positive-time branch ($`t_{n+1},t_{n+2},\mathrm{},t_{n+m}>0`$) labelled by charges $`\{\xi _j=\pm 1\},\{\eta _j=\pm 1\}`$, the interaction factor or influence function is $`F_{n,m}=G_{n,m}H_{n,m}`$, where

$`G_{n,m}`$ $`=\mathrm{exp}\left\{{\displaystyle \underset{k=2}{\overset{n+m}{}}}{\displaystyle \underset{j=1}{\overset{k1}{}}}\xi _k\xi _jQ^{}(t_kt_j)\right\},`$ $`H_{n,m}`$ $`=\mathrm{exp}\left\{i{\displaystyle \underset{j=1}{\overset{n+m1}{}}}\eta _j\theta _{j,n+m}\right\}`$ (7)
with $`\theta _{j,n+m}=_{k=j+1}^{n+m}\xi _kQ^{\prime \prime }(t_kt_j)`$. The charges $`\{\eta _j\}`$ describe propagation on the infinite RDM in diagonal direction, whereas the charges $`\{\xi _j\}`$ indicate moves perpendicular to the diagonal, i. e. represent quantum fluctuations. The pair interaction $`Q(t)`$ is the second integral of the autocorrelation function of the random force $`f(t)=_\alpha c_\alpha x_\alpha (t)`$ in thermal equilibrium, $`\ddot{Q}(t)=(q_0/\mathrm{})^2f(t)f(0)_{\mathrm{eq}}`$ . In the Ohmic scaling limit, the pair interaction reads

$`Q(t)`$ $`=Q^{}(t)+iQ^{\prime \prime }(t)=2K\mathrm{ln}\left[(\mathrm{}\beta \omega _c/\pi )\mathrm{sinh}(\pi |t|/\mathrm{}\beta )\right]+i\pi K\mathrm{sgn}(t).`$ (8)
The influence of the bias $`ϵ`$ or voltage drop $`V`$ is included in the bias phase factor

$`b_{n,m}`$ $`=\mathrm{exp}(i\varphi _{n,m}),`$ $`\varphi _{n,m}`$ $`=ϵ{\displaystyle \underset{l=1}{\overset{n+m}{}}}\xi _lt_l=ϵ{\displaystyle \underset{l=1}{\overset{n+m1}{}}}g_{l,n+m}\rho _l.`$ (9)
Here, we have introduced as an off-diagonal measure the cumulative charge after $`l`$ jumps $`g_{l,n+m}=_{k=l+1}^{n+m}\xi _k=_{k=1}^l\xi _k`$, and the time intervals $`\rho _l=t_{l+1}t_l`$.

Following the lines of Ref. , it is straightforward to arrive at formally exact series expressions for the Laplace transforms of the above quantities in which the summation over all possible sequences $`\{\eta _j\}`$ has been performed already.

The first moment in the stationary nonequilibrium state with initial condition $`q(t=0)_{\mathrm{ne}}=0`$ is given by

$`\widehat{q}(\lambda )_{\text{ne}}`$ $`=\widehat{q}(\lambda )_{\text{fc}}+\widehat{R}_{\text{ne}}(\lambda ).`$ (10)
The term $`\widehat{q}(\lambda )_{\text{fc}}`$ is the first moment with factorized system-bath initial state at time zero. This term has dynamics only in the positive-time branch. The residual term $`\widehat{R}_{\text{ne}}(\lambda )`$ describes the corrections due to the correlated initial state,

$`\widehat{q}(\lambda )_{\text{fc}}`$ $`={\displaystyle \frac{i}{\lambda ^2}}{\displaystyle \underset{m\mathrm{even}}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{0,m}(\lambda ){\displaystyle \underset{\{\xi _l\}}{}}a_{0,m}^{(1)}G_{0,m}\mathrm{sin}\varphi _{0,m},`$ (11)
$`\widehat{R}_{\text{ne}}(\lambda )`$ $`={\displaystyle \frac{i}{\lambda }}{\displaystyle \underset{n,m=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{n,m}(\lambda ){\displaystyle _0^{\mathrm{}}}𝑑\tau _n𝑑s_ne^{\lambda s_n}{\displaystyle \underset{\{\xi _l\}}{}}a_{n,m}^{(1)}G_{n,m}\mathrm{sin}\varphi _{n,m}.`$ (12)
The sum is over all charge sequences $`\{\xi _l\}`$ which form a neutral cluster, $`_{l=1}^{n+m}\xi _l=0`$. The integration symbol includes the $`\mathrm{\Delta }`$-factors associated with each jump and the integrations over the interval lengths except for the interval $`\rho _n`$,

$$_0^{\mathrm{}}\stackrel{~}{𝒟}_{n,m}(\lambda )\mathrm{}=\mathrm{\Delta }^{n+m}\left(\underset{l=n+1}{\overset{n+m1}{}}_0^{\mathrm{}}𝑑\rho _le^{\lambda \rho _l}\right)\underset{k=1}{\overset{n1}{}}_0^{\mathrm{}}𝑑\rho _k\mathrm{}.$$
(13)
The interval $`\rho _n`$ around $`t=0`$ is divided into the negative-time part $`\tau _n=t_n`$ and the positive-time part $`s_n=t_{n+1}`$ and is treated separately. Finally, the phase factor

$`a_{n,m}^{(1)}`$ $`=i{\displaystyle \frac{q_0}{2}}(1)^{(n+m2)/2}{\displaystyle \underset{j=1}{\overset{n+m1}{}}}\mathrm{sin}(\theta _{j,n+m})`$ (14)
subsumes the interaction between the $`\{\eta _j\}`$ and $`\{\xi _j\}`$ charges. In the Ohmic scaling limit, the phase $`\theta _{j,n+m}`$ is independent of the flip times and is $`\theta _{j,n+m}=\pi Kg_{j,n+m}`$.

Similarly, the series expression for the nonequilibrium susceptibility reads

$`\widehat{\chi }_{\text{ne}}(\lambda )`$ $`={\displaystyle \frac{q_0}{i\mathrm{}\lambda }}{\displaystyle \underset{n,m=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{n,m}(\lambda ){\displaystyle _0^{\mathrm{}}}𝑑\tau _n𝑑s_ne^{\lambda s_n}{\displaystyle \underset{\{\xi _l\}}{}}g_{n,n+m}a_{n,m}^{(1)}G_{n,m}\mathrm{cos}\varphi _{n,m}.`$ (15)
The differential mobility can be obtained either from this form using the relation (6) or from its definition with eq. (11). In any event, we arrive at the series expression

$`\mu _\mathrm{d}(ϵ)`$ $`={\displaystyle \frac{q_0^2}{2\mathrm{}}}{\displaystyle \underset{m\mathrm{even}}{\overset{\mathrm{}}{}}}(1)^{m/21}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{0,m}(\lambda 0){\displaystyle \underset{\{\xi _l\}}{}}\left({\displaystyle \underset{l=1}{\overset{m1}{}}}\mathrm{sin}\theta _{l,m}\right)G_{0,m}B_{0,m},`$ (16)
$`B_{0,m}`$ $`={\displaystyle \frac{}{ϵ}}\mathrm{sin}\varphi _{0,m}=\left({\displaystyle \underset{l=1}{\overset{m1}{}}}g_{l,m}\rho _l\right)\mathrm{cos}\varphi _{0,m}.`$
As well, the mean square displacement can be divided into the second moment with factorized initial state at time zero, $`\widehat{q}^2(\lambda )_{\text{fc}}`$, and the correction term $`\widehat{K}_{\text{ne}}(\lambda )`$,

$`\widehat{C}_{\text{ne}}(\lambda )`$ $`=\widehat{q}^2(\lambda )_{\text{fc}}+\widehat{K}_{\text{ne}}(\lambda ),`$ (17)
$`\widehat{q}^2(\lambda )_{\text{fc}}`$ $`={\displaystyle \frac{1}{\lambda ^2}}{\displaystyle \underset{m\mathrm{even}}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{0,m}(\lambda ){\displaystyle \underset{\{\xi _l\}}{}}a_{0,m}^{(2)}G_{0,m}\mathrm{cos}\varphi _{0,m},`$ (18)
$`\widehat{K}_{\text{ne}}(\lambda )`$ $`={\displaystyle \frac{1}{\lambda }}{\displaystyle \underset{n,m=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{n,m}(\lambda ){\displaystyle _0^{\mathrm{}}}𝑑\tau _n𝑑s_ne^{\lambda s_n}{\displaystyle \underset{\{\xi _l\}}{}}a_{n,m}^{(2)}G_{n,m}\mathrm{cos}\varphi _{n,m}.`$ (19)
The phase factor resulting from the $`\{\eta _j\}`$ summation differs from $`a_{n,m}^{(1)}`$. We have

$`a_{n,m}^{(2)}`$ $`=iq_0a_{n,m}^{(1)}{\displaystyle \underset{k=n+1}{\overset{n+m1}{}}}\mathrm{cot}(\theta _{k,n+m}).`$ (20)
The formally exact series expression for the noise follows from (4) and (17) - (19) as

$`S(\omega )`$ $`=2(e^2/q_0^2)\text{Re}\left[\lambda ^2\widehat{C}_{\text{ne}}(\lambda )\right]_{\lambda =i\omega +0^+}.`$ (21)
Consider now the $`\lambda `$-expansions of eqs. (10) - (12) and eqs. (17) - (19) about $`\lambda =0`$, which determine upon Laplace transformation the behavior of the respective quantities at long times. We first observe that the phase factor $`a_{n,m}^{(1)}`$ vanishes if the system reaches a diagonal state $`g_{l,n+m}=0`$. Therefore, all paths that contribute to the first moment hit the diagonal of the RDM only in the initial and final state. Thus each interval $`\rho _l`$ separates two charged clusters, and for $`\lambda =0`$, the interaction factor $`G_{n,m}`$ at $`T>0`$ or the bias factor $`B_{n,m}`$ at $`T=0`$ ensures the convergence of the respective integration. For the first moment with factorized initial conditions, eq. (11), we find

$`\lambda ^2\widehat{q}(\lambda )_{\mathrm{fc}}`$ $`=F\mu +\lambda q_{\mathrm{}}+𝒪(\lambda ^2),`$ $`q_{\mathrm{}}`$ $`=\underset{\lambda 0}{lim}(/\lambda )\left[\lambda ^2\widehat{q}(\lambda )_{\mathrm{fc}}\right].`$ (22)
The series expression for $`\mu `$ is as in (16), but with $`B_{0,m}=\mathrm{sin}(\varphi _{0,m})/ϵ`$.

The series (12) can be rearranged upon using the relation (we put $`\rho _n=\tau _n+s_n`$)

$`{\displaystyle _0^{\mathrm{}}}𝑑\tau _n𝑑s_ne^{\lambda s_n}A(\tau _n+s_n)`$ $`={\displaystyle \frac{1}{\lambda }}{\displaystyle _0^{\mathrm{}}}𝑑\rho _n[\mathrm{\hspace{0.17em}1}e^{\lambda \rho _n}]A(\rho _n),`$ (23)
where $`A(\rho _n)`$ subsumes all correlations across the intervall $`\rho _n`$. It is then straightforward to see that the $`\lambda `$-dependence of the expression $`\lambda ^2\widehat{R}_{\mathrm{ne}}(\lambda )`$ is fully cancelled by $`\lambda ^2\widehat{q}(\lambda )_{\mathrm{fc}}`$. In the end, we find using eq. (10) for all $`\lambda `$

$`\lambda ^2\widehat{q}(\lambda )_{\mathrm{ne}}`$ $`=F\mu .`$ (24)
This verifies that the mean position with initial value $`q(t=0)_{\mathrm{ne}}=0`$ evolves in the stationary nonequilibrium state for all times as $`q(t)_{\mathrm{ne}}=F\mu t`$.

After these preliminaries, consider now the mean square displacement, eq. (17). In the phase factor $`a_{n,m}^{(2)}`$, the factor $`\mathrm{cos}(\pi Kg_{l,n+m})`$ can be allocated to any interval in the positive-time branch. However, upon employing the relation (23), the series for $`\widehat{C}_{\mathrm{ne}}(\lambda )=\widehat{q}^2(\lambda )_{\text{fc}}+\widehat{K}_{\text{ne}}(\lambda )`$ can be rearranged so that most contributions of the sum in eq. (20) cancel out. The only remaining contribution is the one in which the cosine factor covers the first interval in the positive-time branch, labelled by $`\rho _n`$. Thus, this is the only interval which is allowed to be occupied by a diagonal state, $`g_{n,n+m}=0`$. Since $`\rho _n`$ separates two neutral $`\{\xi _j\}`$-clusters, the attractive bias and charge interaction factors are missing, and hence this interval gets very long as $`\lambda 0`$. It is convenient to factorize the interaction term $`G_{n,m}`$ ($`n`$, $`m`$ even) into the intra-cluster interaction factors $`G_{n,0}`$ and $`G_{0,m}`$ and the inter-cluster interaction factor $`G_{\mathrm{int}}`$, $`G_{n,m}=G_{n,0}G_{0,m}G_{\mathrm{int}}`$. Expanding the inter-cluster interaction stretching over the long interval $`\rho _n`$, we get

$`G_{\mathrm{int}}`$ $`=1\ddot{Q}^{}(\rho _n){\displaystyle \underset{l=1}{\overset{n1}{}}}g_{l,n+m}\rho _l{\displaystyle \underset{k=n+1}{\overset{n+m1}{}}}g_{k,n+m}\rho _k.`$ (25)
The first term is the contribution in the absence of the cluster interaction, and the second term represents the dipole-dipole interaction. Accordingly, we get two contributions to $`\widehat{C}_{\mathrm{ne}}(\lambda )`$, called $`\widehat{C}_{\mathrm{ne}}^{(a)}(\lambda )`$ and $`\widehat{C}_{\mathrm{ne}}^{(b)}(\lambda )`$. In the first case, the cluster to the left of $`\rho _n`$ is just $`F\mu =\lambda ^2\widehat{q}(\lambda )_{\text{ne}}`$, whereas the cluster to the right of $`\rho _n`$ represents $`\lambda ^2\widehat{q}(\lambda )_{\text{fc}}`$ \[Note that each cluster is actually an infinite series in $`\mathrm{\Delta }^2`$\]. Since the two clusters are noninteracting, the integration over $`\rho _n`$ gives a factor $`1/\lambda `$. Thus we find

$`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(a)}(\lambda )`$ $`=2\left[F\mu \right]\left[\lambda ^2\widehat{q}(\lambda )_{\mathrm{fc}}\right]/\lambda =2F^2\mu ^2/\lambda +2F\mu q_{\mathrm{}}+𝒪(\lambda ).`$ (26)
In the second form, we have used the result (22).

Consider next the dipole contribution $`\widehat{C}_{\mathrm{ne}}^{(b)}(\lambda )`$. Upon inserting the second term of eq. (25) into the series for $`\widehat{C}_{\mathrm{ne}}(\lambda )`$, we observe that the two sums can be generated by suitable differentiations with respect to the bias. Thus we find

$`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(b)}(\lambda )`$ $`={\displaystyle \frac{q_0^2}{2\lambda ^2}}{\displaystyle \underset{n\mathrm{even}}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m\mathrm{even}}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}𝑑\rho _ne^{\lambda \rho _n}\ddot{Q}^{}(\rho _n){\displaystyle _0^{\mathrm{}}}\stackrel{~}{𝒟}_{n,m}(\lambda )`$ (27)
$`\times {\displaystyle \underset{\{\xi _l\}_n}{}}{\displaystyle \underset{\{\xi _k\}_m}{}}G_{n,0}G_{0,m}\left({\displaystyle \frac{}{ϵ}}\mathrm{sin}\varphi _{n,0}\right)({\displaystyle \frac{}{ϵ}}\mathrm{sin}\varphi _{0,m}){\displaystyle \underset{l=1,ln}{\overset{n+m1}{}}}\mathrm{sin}\theta _{l,n+m}.`$
The sums $`_{\{\xi _l\}_n}`$ and $`_{\{\xi _k\}_m}`$ cover all neutral charge clusters, i.e., they are constrained by $`_{l=1}^n\xi _l=0`$ and $`_{l=n+1}^{n+m}\xi _l=0`$.

At $`T=0`$, the dipole-dipole-interaction is found from eq. (8) as $`\ddot{Q}^{}(t)=2K/t^2`$. Substituting this form, the integration over $`\rho _n`$ in the limit $`\lambda 0`$ yields

$`{\displaystyle _0^{\mathrm{}}}𝑑\rho _ne^{\lambda \rho _n}\ddot{Q}^{}(\rho _n)2K\lambda \mathrm{ln}(\lambda /\lambda ^{}).`$ (28)
Now, in leading order of $`\lambda `$ both the cluster on the left and the cluster on the right of $`\rho _n`$ can be identified with the differential mobility $`\mu _\mathrm{d}`$, eq. (16). Thus we find

$`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(b)}(\lambda )=(4\mathrm{}^2K/q_0^2)\mu _\mathrm{d}^2\lambda \mathrm{ln}(\lambda /\lambda ^{})[1+𝒪(\lambda )].`$ (29)
So far, we have considered the contributions to $`\widehat{C}_{\mathrm{ne}}(\lambda )`$ in which the period $`\rho _n`$ the system dwells in a diagonal state is very long. However, there are also contributions in which the diagonal state is either occupied for a short period only or not occupied at all, $`g_{n,n+m}0`$. In both cases, the charges form a single neutral cluster with regular behavior in the limit $`\lambda 0`$. These terms are combined in the expression

$`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(c)}(\lambda )=2D^{}+𝒪(\lambda ).`$ (30)
Putting all contributions together, the mean square displacement takes the form

$`\lambda ^2\widehat{C}_{\mathrm{ne}}(\lambda )={\displaystyle \frac{2}{\lambda }}F^2\mu ^2+(2F\mu q_{\mathrm{}}+2D^{}){\displaystyle \frac{4\mathrm{}^2K}{q_0^2}}\mu _\mathrm{d}^2\lambda \mathrm{ln}(\lambda /\lambda ^{})+𝒪(\lambda ).`$ (31)
Turning to the time regime, the expression (31) gives the behavior at long times,

$`C_{\mathrm{ne}}(t)`$ $`=F^2\mu ^2t^2+2Dt+(4\mathrm{}^2K/q_0^2)\mu _\mathrm{d}^2\mathrm{ln}(t/t^{}).`$ (32)
Here, all terms linear in $`t`$ are combined in the diffusion constant

$`D`$ $`={\displaystyle \frac{1}{2}}\underset{t\mathrm{}}{lim}{\displaystyle \frac{}{t}}\left[C_{\mathrm{ne}}(t)q(t)_{\mathrm{ne}}^2\right]=D^{}+F\mu q_{\mathrm{}}.`$ (33)
The current correlation function at long times is found from eqs. (3) and (32) as

$`j(t)j(0)+j(0)j(t)_{\mathrm{ne}}=2(eF\mu /q_0)^2{\displaystyle \frac{2\mathrm{}e^2\mu _0}{\pi q_0^2}}{\displaystyle \frac{[\stackrel{~}{\mu }_\mathrm{d}]^2}{t^2}}.`$ (34)
Here, we have normalized the differential mobility on the mobility of a free Brownian particle, $`\stackrel{~}{\mu }_\mathrm{d}=\mu _\mathrm{d}/\mu _0`$ with $`\mu _0=q_0^2/2\pi \mathrm{}K`$.

Using (21), it is easy to commute the $`\lambda `$-expansion (31) into the low-frequency noise spectrum. The $`1/\lambda `$-term in the expansion (31) does not contribute to the noise. The term of order $`\lambda ^0`$ directly relates the dc noise $`S(\omega =0)`$ to the diffusion constant,

$`S(\omega =0)=4e^2D/q_0^2.`$ (35)
In the integrable BSG model, a basis of interacting quasiparticles which scatter off the impurity or point contact one by one can be introduced. Upon generalizing a Landauer-type approach familiar from transport of free electrons to this interacting case, it has been argued that the dc noise at $`T=0`$ is connected with the conductance by differentiation with respect to the voltage  or the tunneling amplitude

$`S(\omega =0)=eV\mathrm{\Delta }G/\mathrm{\Delta }=[1/(K1)]eV^2G/V.`$ (36)
In order $`\mathrm{\Delta }^2`$, this formula reproduces shot noise, $`S(\omega =0)=2eVG=2eI`$. The higher order terms describe corrections due to quantum coherence for strong tunneling. The relation (36) can be proven on the basis of (16) - (19) and is a consequence of (i) a generally valid detailed balance relation (DBR) that follows from analytic properties of the pair interaction $`Q(t)`$ and (ii) of a special DBR which is restricted to the scaling limit, i.e. relies on the particular form (8) for $`Q(t)`$ . The general DBR states that at $`T=0`$ there is no real occupation of TB states higher than the initial state. The special DBR states that only those paths for which all $`\{\eta _j\}`$ charges have equal sign, i.e., which are free of backward hops parallel to the diagonal, contribute to the noise.

The low-frequency correction $`\mathrm{\Delta }S(\omega )`$ is captured by the term of order $`\lambda \mathrm{ln}(\lambda /\lambda ^{})`$ in (31). Using (21), we obtain for the noise the nonperturbative universal form

$`\mathrm{\Delta }S(\omega )`$ $`=2\mathrm{}\mu _0(e^2/q_0^2)\stackrel{~}{\mu }_\mathrm{d}^2|\omega |=2\mathrm{}G_0(g)\stackrel{~}{G}_\mathrm{d}^2|\omega |.`$ (37)
In the second form, we have used transport quantities of the related Luttinger model, i. e. the dimensionless differential conductance $`\stackrel{~}{G}_\mathrm{d}=G_\mathrm{d}/G_0(g)`$, normalized with the static limit of the microwave conductance of a quantum wire or the Hall conductance, $`G_0(g)=ge^2/h`$. This result agrees with the findings of Ref. .

The algebraic $`1/t^2`$ decay law of the current correlation function (34) and the resulting $`|\omega |`$-singularity in the noise spectrum (37) are signatures of the unscreened dipole-dipole interaction. These forms are universal since they do not depend explicitly on the parameter $`K`$ or $`g`$ and the applied bias. The parameters are only in the mobilities or conductances which determine the prefactor of the $`1/t^2`$\- and $`|\omega |`$-law. The forms (34) and (37) are analogous to the $`1/t^2`$ decay law of the position correlation function  in the Ohmic two-state system and the equivalent Shiba relation in frequency space, respectively . The difference, however, is in the prefactor. For the position autocorrelation function, each of the two sets of dipoles corresponds to the static susceptibility $`\chi _0`$, whereas in the present case each dipole set represents the differential mobility, which is related to the nonequilibrium susceptibility by eq. (6). In fact, in transport problems with a stationary nonequilibrium state, the static susceptibility diverges, and the differential mobility $`\mu _\mathrm{d}(F)`$ is the proper quantity here.

At finite temperature, the dipole-dipole interaction is given by

$`\ddot{Q}^{}(t)=2K(\pi /\mathrm{}\beta )^2/\mathrm{sinh}^2(\pi t/\mathrm{}\beta ).`$ (38)
Performing the integration over the interval $`\rho _n`$ between the dipoles, eq. (28), and switching to the noise, the low-frequency contribution from $`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(b)}(\lambda )`$ takes the form

$`\mathrm{\Delta }S(\omega ,T)`$ $`={\displaystyle \frac{e^2}{\pi K}}\left[\stackrel{~}{\mu }_\mathrm{d}\right]^2\omega \mathrm{coth}\left(\mathrm{}\beta \omega /2\right){\displaystyle \frac{2e^2}{\pi K\mathrm{}\beta }}\left[\stackrel{~}{\mu }_\mathrm{d}\right]^2+{\displaystyle \frac{e^2\mathrm{}\beta }{6\pi K}}\left[\stackrel{~}{\mu }_\mathrm{d}\right]^2\omega ^2+𝒪(\omega ^4).`$
Thus, the $`|\omega |`$-singularity is smoothed out at finite temperature. We get a contribution to the dc noise which is proportional to $`T`$ and a correction term of order $`\omega ^2`$. However, there are additional contributions $`\omega ^0`$ and $`\omega ^2`$ from $`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(a)}(\lambda )`$ and $`\lambda ^2\widehat{C}_{\mathrm{ne}}^{(c)}(\lambda )`$. However, the corresponding prefactors can not be given in analytic form for all $`K`$.


format: text

text:
# The Gluon Propagator without lattice Gribov copies

## 1 Introduction

Over the last twenty years, widely different conjectures have been proposed for the infrared behaviour of the gluon propagator. Although it is a gauge dependent quantity, it can be discussed in a given gauge. Even within the same gauge, the proposals for the infrared dependence differ drastically . We mainly summarize here the results that are given in the literature within the Landau gauge, since that gauge is widely used in studies of Dyson Schwinger equations (DSE) as well as in lattice QCD. Early predictions were obtained by solving approximately the DSE. Mandelstam  obtained a solution of a set of truncated DSE equations with an infrared behaviour of the form $`(q^2)^2`$ for the gluon propagator. Such an infrared enhancement was shown, if obtained in any gauge, to lead to an area law for the Wilson loop  and thus to be sufficient for confinement. Infrared enhancement was assumed in various phenomenological studies  and corroborated by later studies of DSE with refined approximations . A different perspective was taken by Gribov , who showed that avoiding gauge copies one would obtain a gluon propagator which vanishes in the infrared in the Landau and Coulomb gauges, of the form

$$D(q^2)\frac{q^2}{q^4+m^4}.$$
(1.1)
An infrared suppressed behaviour was advocated by Stingl , and recently by others , as a possible solution to DSE. Following a procedure similar to that by Gribov, Zwanziger  gave arguments to show that, on the lattice, for any finite spacing in the limit of infinite volume, $`D(q^2=0)=0`$.

We will also consider in this work the parametrization deduced by Cornwall  using a resummation of Feynman graphs which leads to gauge-invariant amplitudes. The gluon propagator is obtained as a solution to this special set of DSE where the claim is that the only gauge dependence appears in the free part. Cornwall’s solution, in addition to fulfilling the Ward identities, allows a dynamical mass generation. Thus this formulation has the additional attractive feature that the gluon mass vanishes in the ultraviolet as required perturbatively. Since the self-energy obtained by Cornwall is claimed to be gauge-independent, we will use his model to fit the propagator both in the Landau and in the Laplacian gauge.

In contrast to all the approaches described above, lattice QCD provides a framework for the calculation of the gluon propagator starting directly from the QCD Lagrangian and can thus yield a conclusive result. Attempts to calculate the gluon propagator started more than ten years ago  on rather small lattices. These early results could be interpreted in terms of a massive scalar propagator, but confirmed the expectation that a Lehmann-Källen representation is not applicable: positivity of the transfer matrix is lost after the non-local gauge fixing. Results on larger lattices were accounted for by assuming a positive anomalous dimension . Recently, a detailed study of the gluon propagator on very large lattices  has been performed, which makes an impressive effort towards bringing under control errors due to the finite lattice spacing and to the finite lattice volume. However, up to now, all lattice studies have used a similar implementation of the Landau gauge on the lattice. Gauge-fixing is accomplished by using a local iterative procedure which identifies local stationarity, but in general fails to determine the global extremum. Which local extremum (“lattice Gribov copy”) is selected depends on the starting condition. These lattice Gribov copies cannot be eliminated. In this situation, their effect has repeatedly been claimed to be small . As discussed in Section 3, we are not convinced by such claims. Results obtained so far will have a safer foundation if the effects of lattice Gribov copies are better understood.

In this work we address the problem of Gribov copies. We use a different gauge condition, which produces a smooth gauge field like the Landau gauge, but which specifies the gauge uniquely: no ambiguity arises due the lattice gauge fixing procedure. This is accomplished by using the Laplacian gauge . The motivation and implementation of this gauge are given in section 3.

We calculate the gluon propagator in quenched QCD on lattices of sizes $`8^4,16^4`$ and $`16^3\times 32`$ at $`\beta =5.8`$ and $`6.0`$, in an attempt to study its zero-temperature behaviour. Our procedure can be extended straightforwardly to finite temperature where the infrared behaviour of the propagator yields the chromo-electric and chromo-magnetic screening masses. The results that we obtain, within the Laplacian gauge, show the same ultraviolet behaviour as in Landau gauge. However, there are significant modifications in the infrared. In particular we find that the zero-momentum propagator is finite, obeys scaling, and becomes volume independent for large enough volumes. It should not however be used as a definition of the gluon mass, since the zero-momentum limit of the propagator is gauge dependent. It is simply a measure of the susceptibility of the gauge-fixed field $`A_\mu `$ in the Laplacian gauge. A quantity which instead can be shown to be gauge independent to all orders in perturbation theory is the pole mass of the transverse part $`D(q^2)`$ of the propagator . To determine this pole if it exists at all, an extrapolation to negative $`q^2`$ is necessary. We compare the inverse propagator $`D^1(q^2)`$ in the Laplacian and the Landau gauges. Using a variety of extrapolation ansätze, in particular a fit to Cornwall’s model  which describes the momentum dependence of our results rather well, we find that data in the Laplacian gauge give support for the existence of a pole at a mass of $`640(140)`$ MeV. Data at smaller momenta are needed to consolidate this result.

Section II introduces our notation; Section III motivates and describes our choice of the Laplacian gauge; Section IV presents our results. They are summarized in Section V.

## 2 Definition of the gluon propagator

The gluon propagator in the continuum is given by

$$𝒟_{\mu \nu }^{ab}(q)=id^4x0|T[A_\mu ^a(x)A_\nu ^b(0)]|0e^{iq.x}$$
(2.2)
This tensor can be decomposed into a transverse and a longitudinal part:

$$𝒟_{\mu \nu }^{ab}(q)=\left(\delta \mu \nu \frac{q_\mu q\nu }{q^2}\right)\delta ^{ab}D(q^2)+\frac{q_\mu q_\nu }{q^2}\delta ^{ab}\frac{F(q^2)}{q^2}$$
(2.3)
For a covariant gauge $`F(q^2)`$ reduces to a constant and corresponds to the gauge fixing parameter $`\xi `$ which in the Landau gauge is zero. Since we want to make a comparison with the recent results  obtained in the Landau gauge, we study the transverse scalar function $`D(q^2)`$ which can be extracted from $`𝒟_{\mu \nu }^{ab}(q)`$:

$$D(q^2)=\frac{1}{3}\left\{\underset{\mu }{}\frac{1}{8}\underset{a}{}𝒟_{\mu \mu }^{aa}(q)\right\}\frac{1}{3}\frac{F(q^2)}{q^2}.$$
(2.4)
$`F(q^2)`$ is determined by projecting the longitudinal part of $`𝒟_{\mu \nu }^{aa}(q)`$ using the symmetric tensor $`q^\mu q^\nu `$.

On the lattice the dimensionless gluon field can be defined by

$$A_\mu (x+\widehat{\mu }/2)=\frac{1}{2ig_0}\left\{\left[U_\mu (x)U_\mu ^{}(x)\right]\frac{1}{3}\mathrm{Tr}\left[U_\mu (x)U_\mu ^{}(x)\right]\right\}+𝒪(a^3)$$
(2.5)
where $`a`$ is the lattice spacing. One may consider different definitions for the gluon field $`A_\mu `$, accurate to higher order in $`a`$. It has been found  that these different definitions give rise to modifications that can be absorbed in the multiplicative field renormalization constant.

The gluon propagator in momentum space is constructed by taking the discrete Fourier transform of $`A_\mu `$ for each colour component

$$A_\mu ^a(q)=\underset{x}{}e^{iq.(x+\widehat{\mu }/2)}A_\mu ^a(x+\widehat{\mu }/2)$$
(2.6)
where the discrete momentum $`q=(q_\mu ,\mu =1,..,4)`$ takes values

$$q_\mu =\frac{2\pi }{a\mathrm{L}_\mu }n_\mu ,n_\mu =(\frac{1}{2}\mathrm{L}_\mu 1),\mathrm{},(\frac{1}{2}\mathrm{L}_\mu )$$
(2.7)
and the momentum-space gluon propagator $`D_{\mu \nu }^{ab}(q)`$ is defined by

$$\mathrm{V}\delta (qq^{})D_{\mu \nu }^{ab}(q)=A_\mu ^a(q)A_\nu ^b(q^{}).$$
(2.8)
with $`V`$ the lattice volume. In the ultraviolet the gluon propagator is expected to behave like $`1/q^2`$. Since on the lattice the free massless propagator behaves as

$$D(q)=\frac{1}{_\mu \left(\frac{2}{a}\mathrm{sin}\frac{q_\mu a}{2}\right)^2},$$
(2.9)
to reduce errors due to the finite lattice spacing we take as our momentum variable the usual

$$\widehat{q}_\mu =\frac{2}{a}\mathrm{sin}\frac{q_\mu a}{2}$$
(2.10)
To relate the bare lattice propagator to the renormalized continuum propagator $`D_\mathrm{R}(q;\mu )`$ one needs the renormalization constant $`Z_3(\mu ,a)`$:

$$a^2D(qa)=Z_3(\mu ,a)D_\mathrm{R}(q;\mu ).$$
(2.11)
Imposing a renormalization condition such as

$$D_\mathrm{R}(q)|_{q^2=\mu ^2}=\frac{1}{\mu ^2}$$
(2.12)
at a renormalization scale $`\mu `$ allows a determination of $`Z_3(\mu ,a)`$. Connection to other continuum renormalization schemes can then be made.

## 3 Gauge fixing procedure

### 3.1 Motivation

The gluon propagator is normally considered in Landau gauge, $`_\mu A_\mu (x)=0x`$. On the lattice, this condition becomes:

$$F(\mathrm{\Omega })\underset{x,\mu }{}\mathrm{Re}\mathrm{Tr}(\mathrm{\Omega }(x)^{}U_\mu (x)\mathrm{\Omega }(x+\widehat{\mu }))\mathrm{maximum}$$
(3.13)
The gauge-fixing functional $`F`$ has many local maxima. To specify the gauge uniquely, the gauge condition above refers to the global maximum. This defines the Fundamental Modular Region (FMR) Landau gauge. In practice however, the gauge transformation $`\mathrm{\Omega }`$ is found by an iterative local maximization of $`F`$, which terminates when any local maximum has been reached. A different gauge condition is thus implemented, which one might call the random Landau gauge, and which depends on the details of the maximization procedure.

It is commonly believed that the effect of choosing a local maximum of (3.13) rather than the global maximum is small, so that the “random” Landau gauge is a good approximation to the FMR Landau gauge. The following argument is often presented to support this view. A given gauge configuration is gauge-fixed $`n`$ times, each time after performing a random gauge transformation; this procedure generates many gauge copies, each corresponding to the local maximum nearest to the random starting point along the gauge orbit. It is observed  that the difference between gluon propagators measured on copies corresponding to the largest and the smallest values of (3.13) is found to be statistically insignificant. A possible problem with this argument however is that the number $`n`$ of gauge copies considered in such comparisons (typically 30 or less) is extremely small compared to the total number of local extrema of (3.13): for simple entropic reasons, all copies considered miss the global maximum by similar amounts, and no reliable information can be extracted about the gluon propagator in the global maximum configuration. It is therefore possible, and we believe quite likely, that the “random” Landau gauge and the FMR Landau gauge are significantly different.

Further evidence for this situation has recently been provided in another gauge, the Direct Maximal Center (DMC) gauge . Although the functional $`F_{DMC}(\mathrm{\Omega })`$ to be maximized differs from (3.13), a similar approach of local iterative maximization is taken, leading to the “random” DMC gauge, with similar problems. In this case however, it is also possible to converge to a large value $`\stackrel{~}{F}_L`$ of $`F_{DMC}`$ by starting from a Landau gauge copy (“Landau” DMC gauge). This value $`\stackrel{~}{F}_L`$ can then be compared with the values obtained from $`n`$ random starting points. One may fit the maximum value among $`n`$ copies, $`\stackrel{~}{F}(n)`$, by a reasonable ansatz like a series in $`1/n`$, and extrapolate to $`n\mathrm{}`$. It turns out that the extrapolated value falls well below $`\stackrel{~}{F}_L`$, which is itself below the global maximum . Furthermore, the properties of the gauge-fixed field are qualitatively different between the “random” and the “Landau” DMC gauges: the former confines after center projection, while the latter does not .

Since in Landau gauge as in DMC gauge, the number of local maxima is expected to grow exponentially with the lattice volume, we expect a similar situation in Landau gauge, leading to large differences between the “random” (local maximum) and the FMR (global maximum) gauges. One might argue that this is not a problem, and that the local maximization of (3.13) implements in the thermodynamic limit a well-defined, but stochastic gauge condition. The relationship between that gauge condition and its perturbative version $`_\mu A_\mu (x)=0`$ is unclear however. Therefore, one should consider the possible effects of selecting a local rather than the global maximum of (3.13) with a great deal of caution. This is the motivation for our study of the gluon propagator in a well-defined, unambiguous gauge.

### 3.2 $`SU(3)`$ Laplacian gauge fixing

In , Vink and Wiese proposed a simple method to fix the gauge unambiguously in $`SU(N)`$. It uses $`N`$ auxiliary Higgs fields, which are chosen as the $`N`$ lowest-lying eigenvectors $`v^{(i)}`$ of the covariant Laplacian. Under a local gauge transformation $`\mathrm{\Omega }(x)`$, these eigenvectors transform covariantly: $`v^{(i)}(x)\mathrm{\Omega }(x)v^{(i)}(x)`$. Therefore, the gauge can be fixed by requiring, at each space-time point $`x`$, $`\{\mathrm{\Omega }(x)v^{(i)}(x),i=1,..,N\}`$ to have some predefined orientation in color space. Specifically, each eigenvector $`v^{(i)}(x)`$ has $`N`$ complex color components, so that the $`N`$ eigenvectors form a complex $`N`$ by $`N`$ matrix $`M`$. Ref. projects this matrix onto $`SU(N)`$ by polar decomposition: $`M=WP,WU(N),P=(M^{}M)^{1/2}`$. The required gauge transformation is then $`\mathrm{\Omega }(x)=e^{i\alpha }W^{}`$, where $`\alpha =\frac{1}{N}\mathrm{arg}(detW)`$. $`\mathrm{\Omega }(x)`$ rotates $`M`$ “parallel” to the identity $`\mathrm{𝟏}_N`$ at each space-time point. The gauge is unambiguously defined, except for these gauge configurations where some of the $`N`$ lowest eigenvalues are degenerate. Such configurations are genuine Gribov copies; they never occur in practice. This approach has been tested for $`SU(2)`$ and $`U(1)`$  and it was shown to reduce to the Landau gauge in the continuum limit aside from exceptional configurations (e.g. an instanton background). Here, we use a slightly modified procedure which requires only $`(N1)`$ eigenvectors (2 for $`SU(3)`$), as follows .

First, apply a gauge transformation $`\mathrm{\Omega }^{(1)}(x)`$ which rotates $`v^{(1)}(x)`$ to $`\left(\begin{array}{c}|v^{(1)}(x)|\\ 0\\ 0\end{array}\right)`$. Five real components of the rotated $`v^{(1)}(x)`$ must vanish, which specifies five constraints. Therefore $`\mathrm{\Omega }^{(1)}(x)`$ is not fully specified, but has $`85=3`$ degrees of freedom. Any satisfactory $`\mathrm{\Omega }^{(1)}`$ can be used.

To completely fix the gauge, we use the second eigenvector $`v^{(2)}`$, already rotated by $`\mathrm{\Omega }^{(1)}`$ to $`\left(\begin{array}{c}v_1^{(2)}\\ v_2^{(2)}\\ v_3^{(2)}\end{array}\right)`$. Three additional constraints are obtained by requiring $`v^{(2)}`$ to be rotated to $`\left(\begin{array}{c}v_1^{(2)}\\ \sqrt{|v_2^{(2)}|^2+|v_3^{(2)}|^2}\\ 0\end{array}\right)`$. This fixes the gauge completely and uniquely.

Note that the second rotation is in an $`SU(2)`$ subgroup, since it leaves $`v_1^{(2)}`$ untouched. This indicates how to generalize this construction to $`SU(N)`$: the first rotation fixes $`(2N1)`$ constraints, which leaves $`(N^21)(2N1)=((N1)^21)`$ degrees of freedom, forming a subgroup $`SU(N1)`$. The next step reduces the gauge freedom to $`SU(N2)`$, etc… down to $`SU(2)`$. It is easily seen that, in this recursive procedure, the matrix $`M`$ is reduced to upper triangular form (with real positive diagonal elements) by the rotation $`\mathrm{\Omega }(x)`$. This is why the $`N^{th}`$ eigenvector needs not be computed: it is only transformed by a phase, which is separately determined by the requirement that $`\mathrm{\Omega }(x)SU(N)`$. Our procedure can thus be viewed as a $`QR`$ decomposition of $`M`$. The gauge, which is globally well-defined (provided the $`N`$ eigenvalues are distinct), may be ill-defined on a sub-manifold of points $`x`$ where our recursive process breaks down. It can be seen that such local gauge defects occur at isolated points, where for $`SU(3)`$, $`\sqrt{|v_2^{(2)}|^2+|v_3^{(2)}|^2}=0`$. The correlation of these points with instantons is studied in .

The Laplacian gauge so defined has the great virtue of being unambiguous. Hence it is the appropriate tool to address our concern about the effect of local extrema of the usual Landau gauge. It also has strong similarities with Landau gauge: it is smooth, Lorentz-symmetric, and gauge-fixes a pure gauge lattice configuration (gauge-transformed from $`U_\mu (x)=\mathrm{𝟏}x,\mu `$) back to $`U=\mathrm{𝟏}`$. Nevertheless, it is a different gauge: its perturbative definition is under consideration ; it differs from Landau gauge most strongly where the magnitude of the eigenvectors $`|v^{(1,2)}(x)|`$ becomes small.

## 4 Results

Since most of the previous studies were performed in the Landau gauge, it is important to compare our Laplacian-gauge propagator with the Landau-gauge one. For this purpose, we have taken, for our analysis, lattice configurations available on the Gauge Connection database , which had already been gauge fixed to Landau gauge with the usual local Over-Relaxation method . These are 200 configurations of a $`16^3\times 32`$ lattice, at $`\beta =5.8`$ and $`6.0`$ each.

The transverse gluon propagator is shown in Fig. 1 for the two gauges at $`\beta =6.0`$. As expected the ultraviolet behaviour is identical in the two gauges, whereas in the infrared, which is the region of interest, significant differences are visible. Since we use a different gauge, this should not come as a surprise. We show the usual quantity $`\widehat{q}^2D(q^2)`$. The Laplacian propagator is clearly not as large as the Landau propagator at low momenta.

The difference between Landau and Laplacian gauge can also be seen in the deviation of $`F(q^2)`$ from zero. Whereas in Landau gauge we find that

$$\widehat{q}^\mu \widehat{q}^\nu 𝒟_{\mu \nu }^{aa}1$$
(4.14)
as expected, in the Laplacian gauge $`F(q^2)`$ is not small, and has a maximum at low momenta. The behaviour of $`F(q^2)`$ is shown in Fig. 2 for $`8^4`$ and $`16^3\times 32`$ lattices at $`\beta =6.0`$. Since $`F(q^2=0)`$ can not be obtained by our projection, we only have one point, at the smallest momentum $`2\pi /32`$ on the larger lattice, to ascertain that $`F(q^2)`$ really has a maximum and does not keep diverging as $`q^20`$. But since the data are systematically higher for the smaller lattice than for the larger one, it seems unlikely that increasing the lattice size further would bring the infrared data up and remove the maximum.

It is interesting to examine the volume dependence of the zero-momentum propagator $`𝒟(0)\frac{1}{4}_{\mu \mu }_\alpha 𝒟_{\mu \mu }^{\alpha \alpha }(q^2=0)`$. We note that in order to determine the transverse part of the propagator at zero-momentum, $`D(0)`$, one must subtract from $`𝒟(q^2)`$ $`F(q^2)/q^2|_{q^2=0}`$, which we can only obtain as $`\mathrm{lim}_{q^20}F(q^2)/q^2`$.

¿From Fig. 2 it can be seen that to extract the limit of $`F(q^2)/q^2`$ as $`q^20`$ reliably, one needs more data in the infrared. Therefore, with the volumes at our disposal we can only examine the zero-momentum limit $`𝒟(0)`$ of the full propagator. In the Landau gauge, Zwanziger has argued that $`D(0)=𝒟(0)`$ should vanish in the infinite lattice volume limit . A recent lattice study in $`SU(2)`$ at finite temperature  seems indeed to indicate such a behaviour. With the data on our present volumes the needed subtraction in the Laplacian gauge cannot be reliably performed, and thus we cannot extract $`D(0)`$. What we find is that $`𝒟(0)`$, the zero-momentum propagator, is finite and volume independent for large enough volumes. The volume dependence and scaling of the renormalized zero-momentum propagator in physical units is displayed in Fig. 3 where we collected results from $`\beta =5.8,6.0,6.2`$ and $`6.5`$. To obtain the renormalized propagator we impose the renormalization condition given in eq. (2.12) where we choose the renormalization point to be $`\mu =a^1`$ for $`\beta =6.0`$ i.e. $`\mu =1.943`$ GeV. This determines $`Z_3(\mu ,a_{\beta =6.0})2.312`$. We then use eq. (2.11) to find the ratio of the factors $`Z3`$ for different $`\beta `$ values at the same physical momentum $`q`$, e.g. for $`\beta =5.8`$ at $`q=\mu `$ we have:

$`{\displaystyle \frac{Z_3(\mu ,a_{\beta =5.8})}{Z_3(\mu ,a_{\beta =6.0})}}`$ $`=`$ $`{\displaystyle \frac{a_{\beta =5.8}^2D(qa_{\beta =5.8})}{a_{\beta =6.0}^2D(qa_{\beta =6.0})}}`$ (4.15)
$`=`$ $`0.97(4)`$ (4.16)
In this way we obtain the renormalization factors at all $`\beta `$ values. For $`\beta =6.2`$ and $`\beta =6.5`$ we find 1.04(4) and 1.07(6) respectively as compared to the value at $`\beta =6.0`$. We also obtained consistent results by fitting our data in the ultraviolet regime using the asymptotic one-loop result for $`D_RZ/q^2(1/2\mathrm{ln}(q^2/\mathrm{\Lambda }^2))^{d_D}`$, with $`d_D=13/22`$ as in the Landau gauge since in this regime the results in the Laplacian and Landau gauges are the same. As can be seen from Fig. 3 the renormalized $`𝒟(0)`$ displays reasonable scaling, and appears quite volume-independent for volumes larger than $`1/2`$ fm<sup>4</sup>. We find a value of $`𝒟(0)=5.020(16)`$ GeV<sup>-2</sup>, or $`𝒟(0)^{1/2}=445(3)`$ MeV, corresponding to a length scale of $`0.5`$ fm. Since the zero-momentum propagator measures the susceptibility of the $`A_\mu ^a`$ field, the length associated with it determines the domain over which the gluon field remains correlated in the Laplacian gauge. If the lattice dimensions become of the order of this characteristic length, then one expects finite size effects to become appreciable. This is indeed what is observed, as shown in Fig.3, with an approximate volume dependence of $`exp(V/V_0)`$ with $`V`$ the lattice volume and $`V_0𝒟(0)^2`$.

On the lattice, the Lorentz symmetry is only approximately restored. Lattice artifacts cause some dependence of $`D(q)`$ on the orientation of the vector $`q`$ rather than just on $`q^2`$. To minimize these discretization effects, we filter our data by making a cylindrical cut in momentum along a reference direction $`\widehat{n}=\frac{1}{2}(1,1,1,1)`$, in the same manner as in Ref. . Namely, we only consider momenta obeying the criterion $`|\mathrm{\Delta }\widehat{q}|<2\pi /L_s`$, where $`L_s`$ is the number of sites in the spatial direction, and $`\mathrm{\Delta }\widehat{q}`$ is the momentum transverse to $`\widehat{n}`$ ($`\mathrm{\Delta }\widehat{q}=\widehat{q}\widehat{q}.\widehat{n}\widehat{n}`$). Using these filtered data which allow a direct comparison with , we examine the various proposals discussed in the Introduction for the infrared behaviour of the propagator. We find that Gribov type parametrizations  as well as infrared enhancement of the type $`(q^2)^2`$  are excluded . The ansatz of Marenzoni et al. ,

$$D(q^2)=\frac{Z}{(q^2)^{1+\alpha }+M^2},$$
(4.17)
with a non-perturbative anomalous dimension $`\alpha `$, gives a better description of the lattice data than the aforementioned parametrizations, but, as seen in Fig. 4, underestimates the peak of the propagator. On the other hand, Cornwall  allows for a dynamically generated gluon mass which vanishes at large momentum in accord with perturbation theory. Using a special set of DSE referred to as a gauge invariant “pinch technique”, he obtains the following solution for the gluon propagator

$$D(q^2)=Z\left[\left(q^2+M^2(q^2)\right)\mathrm{ln}\frac{q^2+4M^2(q^2)}{\mathrm{\Lambda }^2}\right]^1$$

with

$$M(q^2)=M\left\{\frac{\mathrm{ln}\left[(q^2+4M^2)/\mathrm{\Lambda }^2\right]}{\mathrm{ln}\left[4M^2/\mathrm{\Lambda }^2\right]}\right\}^{6/11}$$
(4.18)
Cornwall’s proposal provides a reasonable fit to the data over the whole momentum range (with $`\chi ^2/\mathrm{n}.\mathrm{d}.\mathrm{f}=2.5`$). The quality of this fit can be seen in Fig. 4. For comparison we also fitted our data to the form suggested by Leinweber et al.  where two terms were used, one to describe the ultraviolet behaviour of the form $`D_{\mathrm{UV}}\frac{1}{q^2+M^2}L(q^2,M)`$, and one the infrared of the form $`D_{\mathrm{IR}}1/(q^2+M^2)^{1+\alpha }`$. The exact form, referred to as model A, as taken from ref. , is

$`D(q^2)`$ $`=`$ $`Z[{\displaystyle \frac{AM^{2\alpha }}{(q^2+M^2)^{1+\alpha }}}+{\displaystyle \frac{1}{q^2+M^2}}L(q^2,M)]`$ (4.19)
$`L(q^2,M)`$ $`=`$ $`\left\{{\displaystyle \frac{1}{2}}\mathrm{ln}\left[(\mathrm{q}^2+\mathrm{M}^2)(\mathrm{q}^2+\mathrm{M}^2)\right]\right\}^{13/22}`$ (4.20)
This parametrization, which includes one more parameter than Cornwall’s and is purely phenomenological, does fit the data best over the whole momentum range (with $`\chi ^2/\mathrm{n}.\mathrm{d}.\mathrm{f}=1.2`$).

We address the question of scaling by comparing our results at $`\beta =5.8`$ and $`\beta =6.0`$ on the largest lattice. In the scaling regime the renormalized propagator $`D_\mathrm{R}(q;\mu )`$ is independent of the lattice spacing. Therefore, as in ref. , we can use Eq. (2.11) to obtain the following expression for the ratio of unrenormalized lattice propagators at some physical momentum scale $`q`$

$$\frac{D_1(qa_1)}{D_2(qa_2)}=\frac{Z_3(\mu ,a_1)D_\mathrm{R}(q;\mu )/a_1^2}{Z_3(\mu ,a_2)D_\mathrm{R}(q;\mu )/a_2^2}=\frac{Z_1}{Z_2}\frac{a_2^2}{a_1^2}$$
(4.21)
and where the labels $`1,2`$ refer to the data at $`\beta =6.0`$ and $`\beta =5.8`$ respectively. The scaling properties of the lattice gluon propagator can now be investigated using Eq. (4.21) by adjusting the ratios $`Z_1/Z_2`$ and $`a_1/a_2`$. In Fig. 5 we show the two sets of data lying on the best scaling curve. The shifts required along the horizontal and vertical axes determine the ratios of the wavefunction renormalization constants and of the lattice spacings. We find

$$a_{\beta =6.0}/a_{\beta =5.8}=0.71\pm 0.02\mathrm{and}Z_{\beta =6.0}/Z_{\beta =5.8}=1.07\pm 0.075.$$
(4.22)
with strongly correlated errors. The ratio of lattice spacings is in agreement with the value of $`0.72(4)`$ obtained from a detailed analysis of the static potential . The ratio of the Z-factors is within what is expected from perturbation theory, and in agreement with the value of $`1.04(3)`$ of Ref. . In other words, scaling is very well satisfied for the Laplacian gauge, and performing the fits at $`\beta =6.0`$ gives the behaviour of the gluon propagator in the physical regime.

We focus now on the infrared behaviour of the transverse propagator. Figs. 6 and 7 show the inverse propagator as a function of $`\widehat{q}^2`$ in the two gauges. Two advantages of the Laplacian gauge become visible.

First, the orientation of the momentum $`q`$ has less effect than in Landau gauge: the data points at a given value of $`\widehat{q}^2`$ show less scatter, and the cylindrical cut is not as essential as in Landau gauge in the infrared region. At a given lattice spacing, the Laplacian gauge approximates better the Lorentz symmetry of the continuum. This reduction of lattice artifacts is understandable since the gauge is fixed by considering the lowest-lying eigenvectors of the Laplacian, which are the least sensitive to UV-cutoff effects. In contrast, Landau gauge comes from the iteration of a completely local, UV-dominated process. Better rotational symmetry allows for better accuracy, or for the same accuracy on coarser lattices.

Second, the inverse propagator is closer to a linear function of $`\widehat{q}^2`$ in Laplacian gauge. If it were the propagator for a free boson, it would be described by a straight line since $`1/D(q^2)=Z^1(q^2+m^2)`$. Having curvature means that one has a momentum-dependent effective mass $`\mathrm{\Pi }(q^2)`$. In particular, the infrared mass $`\mathrm{\Pi }(0)`$ and the pole mass $`\mathrm{\Pi }(q^2)`$ such that $`(q^2+\mathrm{\Pi }(q^2))=0`$ are different.

The latter is of special interest, because of its gauge independence at least to all orders in perturbation theory. Finding a pole, i.e. a zero of the inverse propagator, requires the extrapolation of our data to negative $`\widehat{q}^2`$. The less curvature in the inverse propagator in the infrared, the more reliable the extrapolation will be.

Three types of extrapolation are displayed in the figures: quadratic and cubic polynomials in $`\widehat{q}^2`$, and our fit to Cornwall’s model. The location of the pole, and even its existence, are affected by the choice of extrapolation in the Landau gauge. The coefficients $`a_1,a_2,a_3`$ of the cubic polynomial extrapolation keep increasing, indicating poor stability. Essentially, no statement about a pole can be made in that gauge. Differentiating between a cubic fit (which gives a pole) and a Cornwall-type fit (which doesn’t) will require extremely accurate data on large lattices. Extracting the pole is also difficult in the Laplacian gauge but at least one finds a pole with all the ansätze that we tried. Given the convexity of the data, a lower bound is provided by a linear fit near $`q^2=0`$, which defines the (gauge-dependent) infrared mass. Quadratic and cubic terms in the polynomial extrapolation represent small corrections of decreasing size. One can thus make some estimate of the gluon pole mass. A similar study on a lattice of double size, as was considered in Ref. , would produce more than four times as many points in the same $`\widehat{q}^2`$ interval, and should allow for an accurate determination of the pole mass.

We also measure the correlator of the gluon field averaged over a time-slice. Namely, we measure

$$C(t)=\frac{1}{L_s^3}\frac{1}{8}\underset{a=1}{\overset{8}{}}\frac{1}{3}\underset{\mu =1}{\overset{3}{}}(\underset{x}{\overset{L_s^3}{}}A_\mu ^a(\stackrel{}{x},0))(\underset{x}{\overset{L_s^3}{}}A_\mu ^a(\stackrel{}{x},t))$$
(4.23)
which is displayed in Fig.8. At large time separations $`t`$, this correlator should decay exponentially like $`exp(m_{pole}t)`$, giving us another approach to extracting the pole mass. We use this observable to perform a crosscheck on this mass, and as a further study of the systematic errors in its determination. This correlator is measured on the same configurations as $`D(q^2)`$, so it contains no additional information. But the same information is given a different weight, so that a fit to $`C(t)`$ will give different results than a fit to $`D^1(q^2)`$, especially after the cylindrical momentum cut. Therefore, we fit Cornwall’s model directly to $`C(t)`$ instead of $`D^1(q^2)`$. Remarkably, the difference is rather small, which attests again to the soundness of the model. The dashed lines in Fig.8 show the original fit of Cornwall’s ansatz to $`D^1(q^2)`$, which already provides a fair description of the data. The solid lines represent a direct fit of the same 3-parameter ansatz to $`C(t)`$, excluding the first few time-slices which otherwise completely dominate the fit. The fit started from $`t=4`$ and $`t=2`$ at $`\beta =6.0`$ and 5.8 respectively, which amounts to discarding similar intervals in physical units. Given the 3 fitted parameters, one can then solve $`D^1(q^2)=0`$ numerically, with $`D(q^2)`$ as per eq.(4.18). The corresponding pole mass varies little from one fit to the other, and remains roughly constant in physical units at $`\beta =5.8`$ and $`6.0`$. Also, a model-independent extraction of the pole mass, by measuring the effective mass $`m_{eff}(t)=Ln(C(t+1)/C(t))`$, gives a consistent value. Taking these results into account, together with the quadratic and cubic extrapolations displayed in Fig.7, we estimate the pole mass to lie in the interval $`[500,785]`$ MeV, where we used $`a^1(\beta =6.0)=1.943`$ GeV to convert to physical units  with $`\sqrt{\sigma }=440`$ MeV. The lower bound is given by the infrared mass, which corresponds to a linear extrapolation of $`D^1(q^2)`$; the upper bound is provided by the largest value obtained when fitting to our data Cornwall’s model. A reasonable central value is 640 MeV, which corresponds to Cornwall’s extrapolation in Fig.7.

We have performed a similar exercise for the Landau gauge. The fit of Cornwall’s model to $`D(q^2)`$ or $`C(t)`$ is quite satisfactory, but the equation $`D^1(q^2)=0`$ gives a complex pole, far from the real axis. Note that Ref.  also finds oscillatory behaviour for the time-slice correlator in $`3d`$ $`SU(2)`$ theory fixed to Landau gauge, reflecting a complex pole. This disagreement with the Laplacian gauge is puzzling, since one expects the pole to be gauge invariant. Possible causes include the inadequacy of the Landau gauge fixing procedure on the lattice, or finite-size effects. Larger volume studies, currently under way, should elucidate this issue.

## 5 Conclusions

We have evaluated the gluon propagator using the Laplacian gauge which avoids lattice Gribov copies. We extracted the transverse part of the gluon propagator and verified its scaling in this gauge. Examining the scaling and volume dependence of the zero-momentum propagator $`𝒟(0)`$, we reached the conclusion that it is a constant beyond a lattice size of $`0.8`$ fm. This size is consistent with the characteristic length scale determined from $`𝒟(0)`$ itself as the range beyond which the gluon field decorrelates in this gauge.

Among the various proposals for the transverse propagator which are physically founded, Cornwall’s model  provides a reasonable fit to the lattice results over the whole momentum range. We find it satisfying that the lattice data seem to favour a model with a dynamically generated mass.

By looking at the inverse propagator $`D^1(q^2)`$ at small momenta, we see that the Laplacian gauge is superior to the Landau gauge in its restoration of Lorentz symmetry on the lattice. Furthermore, it turns out that the inverse propagator is almost linear in $`\widehat{q}^2`$ in the Laplacian gauge. This allows for a more reliable extrapolation to $`\widehat{q}^2<0`$, as compared to the Landau gauge. We test a variety of extrapolation ansätze. They consistently yield a pole mass at $`640\pm 140`$ MeV.

Acknowledgements: The $`16^3\times 32`$ lattice configurations were obtained from the Gauge Connection archive . We thank A. Cucchieri, S. Dürr, M. Pepe and O. Philipsen for helpful discussions.


text:
# Radio Halos of Galaxy Clusters from Hadronic Secondary Electron Injection in Realistic Magnetic Field Configurations

## 1 Introduction

Some clusters of galaxies contain large, Mpc-sized regions of diffuse radio emission, the so called cluster radio halos and cluster radio relics. Although known for decades, these sources receive increasing attention due to the improved detection sensitivities and their important theoretical implications for the non-thermal content of the ICM. For reviews see Valtaoja Valtaoja (1984), Jaffe Jaffe (1992), Feretti & Giovannini Feretti & Giovannini (1996), Feretti Feretti (1999), Enßlin Enßlin (1999b), and Giovannini et al. Giovannini et al. (1999).

Whereas the cluster relic sources are nowadays widely believed to trace the presence of large shock waves in the ICM Enßlin et al. (1998); Roettiger et al. (1999); Enßlin & Krishna (2000) the origin of radio halos is less clear. The lifetime of the radio emitting CRe is much shorter than any reasonable transport time over cluster scales. This means that the CRe need to be continuously re-energized, or freshly injected on a cluster-wide scale. The latter would happen if the CRe result from the decay of charged pions produced in hadronic interaction of energetic CRp distributed throughout the cluster volume Dennison (1980); Vestrand (1982); Valtaoja (1984); Schlickeiser et al. (1987); Blasi & Colafrancesco (1999). The CRp have lifetimes in the IGM of the order of a Hubble time Berezinsky et al. (1997); Enßlin et al. (1997), giving them enough time to diffuse away from their source and to maintain a cluster wide distribution.

It is the purpose of this work to demonstrate that a simple model for hadronic electron injection in a realistic magnetic field configuration leads to radio halos which reproduce several observations: the profile of the radialy decreasing radio emission, the low radio polarization, the radio to X-ray surface brightness correlation, and the cluster radio halo luminosity-temperature relation. We do not attempt to reproduce all observational information on cluster radio halos. For example the observed spatial spectral index variations of individual clusters Giovannini et al. (1993); Deiss et al. (1997) are not reproduced by the presented model. It is our goal to investigate the observational consequences of a conceptually simple model. Sophisticated modifications, which fine-tune the model to the data, are principally possible. But first an understanding of idealized models is most helpful. The assumptions of our model are:

1. CRe are secondary particles of hadronic interactions of CRp with the background protons of the ICM. The origin of the CRp is not specified here, but a likely explanation is particle acceleration in shock waves of cluster merger events and in accretion shocks of the clusters Colafrancesco & Blasi (1998). ICM CRp might also migrate from radio plasma of radio galaxies into the ICM Valtaoja (1984); Enßlin et al. (1997); Enßlin (1999a); Blasi & Colafrancesco (1999) or result from supernova driven galactic winds Volk et al. (1996).
2. The energy density in CRp is assumed to be proportional to the energy density of the thermal ICM gas. This is reasonable if the thermal gas and the CRp were energized in the same shock waves of cluster mergers, since a constant fraction of energy dissipated in the shock should go into the CRp population. Previous work ignored the spatial profiles Dennison (1980); Vestrand (1982) or assumed one single central point source from which the CRp diffused away, resulting in a centrally very peaked CRp profile Blasi & Colafrancesco (1999).
3. The CRp spectral index is independent of position or cluster and therefore the radio spectral index is constant in this model. Although this is certainly an oversimplification since e.g. the radio halo of the Coma cluster does show spectral index variations, it allows to predict radio halo statistics.
4. No significant re-acceleration or particle diffusion of the CRe is occurring.
5. The radio spectrum is assumed to be quasi-stationary, since the electron cooling time is short and a stationary CRe population is rapidly established. The spectrum of this population results from the shape of the proton spectrum, which is assumed to be a single power-law, and the synchrotron and inverse Compton cooling. Hence electron and radio spectra are also power laws.
6. The magnetic field configuration results from amplification of pre-cluster seed fields as simulated by Dolag et al. Dolag et al. (1999). We therefore investigate very inhomogeneously magnetised clusters. Although the importance of realistic field configurations for the appearence of radio halos was stressed by Tribble Tribble (1991), even recent work on halos of hadronic origin assumed spatially constant field strength.

## 2 MHD Cluster Formation Simulation

We used the cosmological MHD code described in Dolag et al. Dolag et al. (1999) to simulate the formation of magnetised galaxy clusters from an initial density perturbation field. The evolution of the magnetic field is followed from an initial seed field. This field is amplified by the compression during the cluster collapse. Merger events and sheer flows, which are very common in the cosmological environment of large-scale structure, lead to Kelvin-Helmholtz instabilities. They further increase the field strength by a large factor. In order to reproduce the present $`\mu `$G fields an initial nG field strength is required.

### 2.1 GrapeSPH

The code combines the merely gravitational interaction of a dark-matter component with the hydrodynamics of a gaseous component. The gravitational interaction of the particles is evaluated on GRAPE boards, while the gas dynamics is computed in the SPH approximation. It was also supplemented with the magneto-hydrodynamic equations to trace the evolution of the magnetic fields which are frozen into the motion of the gas because of its assumed ideal electric conductivity. The back-reaction of the magnetic field on the gas is included. Extensive tests of the code were successfully performed and described in the previous paper. $`B`$ is always negligible compared to the magnetic field divided by a typical length scale. The code also assumes the ICM to be an ideal gas with an adiabatic index of $`\gamma =5/3`$ and neglects cooling. The surroundings of the clusters are dynamically important because of tidal influences and the details of the merger history. In order to account for this the cluster simulation volumes are surrounded by a layer of boundary SPH particles in order to represent accurately the sources of the tidal fields in the cluster neighborhood. For details of the code, the models and the obtained magnetic field structure see Dolag et al. (1999; in preparation).

### 2.2 Initial conditions

For an accurate calculation of the synchrotron emission the resolution of the magnetic field is crucial. We therefore perform simulations in a CDM Einstein-de-Sitter cosmology ($`\mathrm{\Omega }_\mathrm{m}^0=1.0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }^0=0`$, $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$), with twice the mass resolution of the simulation described in Dolag (1999). We simulate a set of ten different realisations which result in clusters of different final masses and different dynamical states at redshift $`z=0`$.

As shown in Dolag 1999, the initial field structure is irrelevant because the final field structure is mostly determined by the dynamics of the cluster collapse. The initial field strengths in this simulations where chosen so that the final field strength reproduces the observed Faraday rotation data collected by Kim et al. Kim et al. (1991). The central values of the ten clusters vary between $`0.5\mu `$G and $`2\mu `$G, whereas Sim.9 reaches values of $`1\mu `$G. The simulated magnetic field strength drops by two orders of magnitudes from the central part out to the virial radius (see Fig. 1).

Given the large uncertainties a field strength twice as high may also be consistent with the Faraday data. Therefore we also give results for the extrapolated case of field strength twice as high. We note that neither the simulation, nor the Faraday rotation measurements can exclude the existence of stronger small-scale ($`1`$ kpc) magnetic field components. Here, we assume only the presence of large-scale ($`20`$ kpc) and relatively weak fields. A comparison of the synthetic Faraday rotation measurements with the observations is shown in Fig. 2.

## 3 Hadronic Secondary Electron Injection

A CRp population may be well described by a power-law in kinetic energy $`E_\mathrm{p}`$

$$f_\mathrm{p}(x,E_\mathrm{p})dE_\mathrm{p}dV=\frac{N_\mathrm{p}(x)}{\mathrm{GeV}}\left(\frac{E_\mathrm{p}}{\mathrm{GeV}}\right)^{\alpha _\mathrm{p}}dE_\mathrm{p}dV.$$
(1)
We assume a lower cutoff at $`E_{\mathrm{p},\mathrm{min}}1\mathrm{GeV}`$. The contribution of He (assuming a 24% mass fraction) to the cosmic rays and thermal gas is implicitly taken into account in our calculations by replacing proton densities by the correct nucleon densities without further notification. We choose the normalization $`N_\mathrm{p}(x)`$ so that the (kinetic) CRp energy density $`\epsilon _{\mathrm{CRp}}(x)`$ is a constant fraction of the thermal energy density $`\epsilon _{\mathrm{th}}(x)`$ of the ICM.

$$\epsilon _{\mathrm{CRp}}=X_{\mathrm{CRp}}\epsilon _{\mathrm{th}}=\frac{N_\mathrm{p}(x)}{\alpha _\mathrm{p}2}\left(\frac{E_{\mathrm{p},\mathrm{min}}}{\mathrm{GeV}}\right)^{2\alpha _\mathrm{p}}\mathrm{GeV}$$
(2)
Since the spectral index is sufficiently steep in our model ($`\alpha _\mathrm{p}2.4`$) the upper cutoff was set to infinity.

The CRp interact hadronically with the background gas and produce pions. The charged pions decay into secondary electrons (and neutrinos) and the neutral pions into gamma-rays:

$`\pi ^\pm `$ $``$ $`\mu ^\pm +\nu _\mu /\overline{\nu }_\mu e^\pm +\nu _e/\overline{\nu }_e+\nu _\mu +\overline{\nu }_\mu `$
$`\pi ^0`$ $``$ $`2\gamma `$
The resulting electron injection spectrum $`q_\mathrm{e}(x,E_\mathrm{e})dE_\mathrm{e}dV`$ peaks at the energy of $`m_\pi c^2/435\mathrm{MeV}`$ due to rough energy equipartition between electrons and neutrinos in the charged pion decay Mannheim & Schlickeiser (1994). The radio emission is produced by electrons with an energy of several GeV (depending on observing frequency and magnetic field strength), which have to result from protons of even higher energies. For this energy range the electron production spectrum can be calculated following Mannheim & Schlickeiser Mannheim & Schlickeiser (1994) to be

$$q_\mathrm{e}(x,E_\mathrm{e})2^6\sigma _{\mathrm{pp}}cn_\mathrm{p}(x)\frac{N_\mathrm{p}(x)}{\mathrm{GeV}}\left(\frac{24E_\mathrm{e}}{\mathrm{GeV}}\right)^{\frac{4}{3}(\alpha _\mathrm{p}\frac{1}{2})},$$
(3)
where $`\sigma _{pp}=32\mathrm{mbarn}`$ is the inelastic p-p cross section, and $`n_\mathrm{p}`$ is the target proton density . We note, that the neutrino production spectrum is $`q_\nu (x,E_\nu )=3q_\mathrm{e}(x,E_\nu )`$, and the $`\pi ^0`$-decay induced gamma-ray spectrum is $`q_\gamma (x,E_\gamma )=3q_\mathrm{e}(x,E_\gamma /4)/(16\alpha _\mathrm{p}8)`$ for $`E_\gamma m_{\pi ^0}c^2/268\mathrm{MeV}`$.

The steady-state CRe spectrum is shaped by the injection and cooling process and is governed by

$$\frac{}{E_\mathrm{e}}\left(\dot{E_\mathrm{e}}(x,E_\mathrm{e})f_\mathrm{e}(x,E_\mathrm{e})\right)=q_\mathrm{e}(x,E_\mathrm{e}).$$
(4)
For $`\dot{E_\mathrm{e}}(x,p)<0`$ this equation is solved by

$$f_\mathrm{e}(x,E_\mathrm{e})=\frac{1}{|\dot{E_\mathrm{e}}(x,E_\mathrm{e})|}_{E_\mathrm{e}}^{\mathrm{}}𝑑E_\mathrm{e}^{}q_\mathrm{e}(x,E_\mathrm{e}^{}).$$
(5)
The cooling of the radio emitting CRe is dominated by synchrotron and inverse Compton losses giving

$$\dot{E_\mathrm{e}}(x,E_\mathrm{e})=\frac{4\sigma _Tc}{3m_\mathrm{e}^2c^4}\left(\frac{B^2(x)}{8\pi }+\frac{B_{\mathrm{CMB}}^2}{8\pi }\right)E_\mathrm{e}^2,$$
(6)
where $`B(x)`$ is the local magnetic field strength and $`B_{\mathrm{CMB}}^2/(8\pi )`$ is the energy density of the cosmic microwave background expressed by an equivalent field strength. The CRe have therefore a power-law spectrum

$$f_\mathrm{e}(x,E_\mathrm{e})=\frac{N_\mathrm{e}(x)}{\mathrm{GeV}}\left(\frac{E_\mathrm{e}}{\mathrm{GeV}}\right)^{\alpha _\mathrm{e}},$$
(7)
with $`\alpha _\mathrm{e}=\frac{4}{3}\alpha _\mathrm{p}+\frac{1}{3}`$ and

$$N_\mathrm{e}(x)=\frac{2^7\mathrm{\hspace{0.17em}3}^2\pi \mathrm{\hspace{0.17em}\hspace{0.17em}24}^{\frac{4}{3}(\alpha _\mathrm{p}\frac{1}{2})}}{4\alpha _\mathrm{p}5}\frac{\sigma _{\mathrm{pp}}m_\mathrm{e}^2c^4}{\sigma _T\mathrm{GeV}}\frac{n_\mathrm{p}(x)N_\mathrm{p}(x)}{B^2(x)+B_{\mathrm{CMB}}^2}.$$
(8)
The radio emissivity $`j_\nu `$ at frequency $`\nu `$ and per steradian is given by

$$j_\nu =c_2(\alpha _\mathrm{e})N_\mathrm{e}(x)B_{}(x)^{\frac{\alpha _\mathrm{e}+1}{2}}\left(\frac{\nu }{c_1}\right)^{\frac{\alpha _\mathrm{e}1}{2}}$$
(9)
Lang (1999) with $`c_1=3e\mathrm{GeV}^2/(2\pi m_\mathrm{e}^3c^5)`$,

$$c_2(\alpha _\mathrm{e})=\frac{\sqrt{3}}{16\pi }\frac{e^3}{m_\mathrm{e}c^2}\frac{\alpha _\mathrm{e}+\frac{7}{3}}{\alpha _\mathrm{e}+1}\mathrm{\Gamma }\left(\frac{3\alpha _\mathrm{e}1}{12}\right)\mathrm{\Gamma }\left(\frac{3\alpha _\mathrm{e}+7}{12}\right),$$
(10)
and $`B_{}(x)`$ being the magnetic field component within the image plane. Integration of $`j_\nu (x)`$ along the $`z`$-direction gives the surface brightness of the radio halo.

Thus we find that the radio emissivity as a function of position in the hadronic secondary electron model scales with the ICM properties as

$$j_\nu X_{\mathrm{CRp}}n_\mathrm{p}\epsilon _{\mathrm{th}}\frac{B_{}^{\alpha _\nu +1}}{B^2+B_{\mathrm{CMB}}^2}\nu ^{\alpha _\nu },$$
(11)
where $`\alpha _\nu =(\alpha _\mathrm{e}1)/2=(2\alpha _\mathrm{p}1)/3`$. It follows that the emissivity of strong field regions ($`B^2B_{\mathrm{CMB}}^2`$) is nearly independent of the magnetic field strength, whereas in weak field regions ($`B^2B_{\mathrm{CMB}}^2`$) it is very sensitive to it. The factor $`n_\mathrm{p}\epsilon _{\mathrm{th}}`$ is nearly identical to the X-ray emissivity $`l_Xn_\mathrm{p}^2(kT_\mathrm{e})^{1/2}`$ so that for a radially decreasing magnetic field strength one always gets a radio luminosity which declines more rapidly than the X-ray emission:

$$\frac{j_\nu }{l_X}X_{\mathrm{CRp}}(kT_\mathrm{e})^{1/2}\frac{B_{}^{\alpha _\nu +1}}{B^2+B_{\mathrm{CMB}}^2}\nu ^{\alpha _\nu },$$
(12)
The normalization constant of Eq. 11 depends on the assumed scaling between CRp and thermal energy density. We fix this relation by comparing Coma to a similar simulated cluster (Sim. 9) and assume it to hold for all clusters in order to predict their radio luminosity.

## 4 Synthetic Radio Maps

To obtain the synthetic radio maps we calculate all quantities needed in a $`5\mathrm{Mpc}^3`$ box containing $`200^3`$ grid cells, centered on the cluster position and distribute the CRe as described before. From this data cube we calculate the maps by projecting the emission along the three spatial directions. As the radio emission depends on the angle between the magnetic field and the observing direction, this leads to different total radio luminosities in different observing directions. As the magnetic field in the simulations is nearly isotropic, the differences are small. Additionally projection effects of substructure in the underlying density distribution can lead to different radial profiles as seen in Fig. 4.

As a reference point we chose one simulated cluster (Sim. 9), which has an emission weighted temperature (9.4 keV) and a velocity dispersion (1200 km/s derived from the virial mass and radius) comparable to the values observed in the Coma cluster of galaxies (9 keV (Donelly et al., 1999), 1140 km/s Zabludoff et al. (1993)). This cluster is in a post-merger stadium, like Coma Burns et al. (1994); Bravo-Alfaro et al. (2000).

We assume $`\alpha _p=2.375`$, which gives $`\alpha _\nu =1.25`$, typical for cluster radio halos, and adjust $`X_{\mathrm{CRp}}`$ to $`0.14(E_{\mathrm{p},\mathrm{min}}/\mathrm{GeV})^{0.375}`$ so that this cluster reproduces the observed total 1.4 GHz radio emission of Coma. In the case of a twice the field strength $`X_{\mathrm{CRp}}`$ needs to be $`0.04(E_{\mathrm{p},\mathrm{min}}/\mathrm{GeV})^{0.375}`$. Note that in this case the magnetic fields and the CRp are in equipartition for $`E_{\mathrm{p},\mathrm{min}}=1\mathrm{GeV}`$.

The synthetic radio maps are shown in Fig. 3 in comparison with the radio map of Coma. In the calculation of the radial radio emission profile of Coma the so called radio bridge of Coma (see Fig. 3, top-left panel at $`x=3.5`$ Mpc, $`y=1.5`$ Mpc) is included. This can be justified because a similar bridge is also visible in the X-ray emission, which suggests that the radio bridge results from ICM substructure. This substructure is believed to be the core of a smaller merging cluster which already passed the core region of Coma Burns et al. (1994); Bravo-Alfaro et al. (2000). The radio relic in Coma (see Fig. 3, top-left panel at $`x=4.5`$ Mpc, $`y=1`$ Mpc) does not affect the comparison of the radio profiles due to its peripheral location.

Although we cannot expect this simulated cluster to match the morphology of the radio halo of Coma exactly, the radial emission profiles in the $`y`$\- and $`z`$-projection are quite comparable (Fig. 4). In the $`x`$-projection the two sub-lumps are projected on top of each other. This leads to a steeper radial radio profile (Fig. 4). The second simulated cluster with flat radio emission profile is also in a post-merger stadium with substantial substructure.

The large extension of low level radio emission of Coma is not reproduced by the model. As the simulation has its highest resolution in the center, the simulation might underestimate the small scale magnetic fields, and therefore the radio emission, in the outer regions. This could explain the less extended synthetic radio halo. Another possibility is that e.g. diffusion established a more extended CRp profile than it is assumed in this work.

The observed radio halos do not show any radio polarization. We estimate for our simulated field configuration the expected degree of radio polarization following the treatment of internal Faraday rotation of Burn Burn (1966). The synthetic polarization maps has peak values of the order of $`10^6`$. This undetectable signal would be further reduced by beam depolarization due to the small-scale polarization pattern.

It is difficult to compare morphological different clusters of galaxies. But Eq. 12 indicates a close relation between radio and X-ray emissivities within the hadronic halo model. This should be approximatively conserved in the projected fluxes, which allows for a less morphology-dependent comparison of different real or simulated clusters. For Sim. 9 we compare X-ray and radio surface brightness point by point. In order to process the data in the same way as real observational data we average over 250 kpc boxes and use the RMS as an error estimate (plus systematic errors of $`\delta L_X=10^7\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ and $`\delta L_\nu =10^{21}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{Hz}^1`$). The results are given in Fig. 5. For instrumentally detectable regions ($`L_\nu >\delta L_\nu `$, $`L_X>\delta L_X`$) the radio and X-ray emissivity can be related via a power law $`L_\nu =aL_X^b`$, where $`a=2.2310^{13}`$ and $`b=1.26`$ in cgs units. Govoni et al. Govoni et al. (2000) give for four clusters a point by point radio to X-ray count rate correlation. Their values ($`b=0.64\pm 0.07`$ for Coma, $`b=0.82\pm 0.04`$ for A2319, $`b=0.94\pm 0.04`$ for A2255, $`b=0.99\pm 0.05`$ for A2744) are significantly smaller, which is also supported by the visual impression of the more concentrated simulated halo in comparison to Coma (Fig. 3). This indicates that the radial CRp profile should be more extended than the thermal energy of the ICM, as we have assumed in this work, if radio halos are of hadronic origin.

## 5 Gamma-Rays and Neutrinos

The luminosity of the simulated cluster Sim. 9 in neutrinos and gamma-rays is $`Q_{\nu _{\mathrm{e},\mu }}(E_\nu )=5.510^{45}E_\nu ^{2.5}/\mathrm{GeV}^{1.5}`$ and $`Q_\gamma (E_\gamma )=5.910^{45}E_\gamma ^{2.5}/\mathrm{GeV}^{1.5}`$. This leads to a gamma ray flux of $`5.210^8\mathrm{cm}^2\mathrm{s}^1`$ above 100 MeV ($``$ since the spectrum is significant below the above power law near $`E_\gamma =100`$ MeV). This is comparable to the upper limit $`410^8\mathrm{cm}^2\mathrm{s}^1`$ for the flux from Coma measured by EGRET (Sreekumar et al., 1996). In the case of twice the field strength one finds $`Q_{\nu _{\mathrm{e},\mu }}(E_\nu )=1.510^{45}E_\nu ^{2.5}/\mathrm{GeV}^{1.5}`$, and $`Q_\gamma (E_\gamma )=1.610^{45}E_\gamma ^{2.5}/\mathrm{GeV}^{1.5}`$, which leads to a flux above 100 MeV of $`1.410^8\mathrm{cm}^2\mathrm{s}^1`$ well below the EGRET limit. For comparison with gamma ray and neutrino flux predictions for clusters of galaxies see Völk et al. Volk et al. (1996), Enßlin et al. Enßlin et al. (1997), Colafrancesco & Blasi Colafrancesco & Blasi (1998), and Blasi Blasi (1999).

## 6 The Temperature-Radio-Luminosity Relation

Clusters with radio halos exhibit a strong correlation between the radio luminosity (at 1.4 GHz: $`L_{1.4\mathrm{GHz}}`$) and the ICM temperature ($`T_X`$) Liang (1999); Colafrancesco (1999). This is expected because of the known correlations of cluster size and energy content with temperature Mohr & Evrard (1997); Schindler (1999); Mohr et al. (1999). But the shape of the $`T_X`$$`L_{1.4\mathrm{GHz}}`$ relation is an important touchstone of every radio halo theory. It can be seen in Fig. 6 that the hadronic halo model reproduces the observed relation if a cluster independent $`\epsilon _{\mathrm{CRp}}/\epsilon _{\mathrm{th}}`$ ratio is assumed.

## 7 Conclusion

A set of simulated clusters of galaxies in a CDM EdS cosmology, including fully MHD simulated magnetic field structures, was used to test a conceptual simple hadronic secondary electron model for radio halo formation of clusters. The seed magnetic fields were chosen so that the final magnetic fields reproduce the observed Faraday rotation measurements. It was assumed further, that in all clusters a spatially constant ratio between the cosmic ray proton and the thermal energy density exists. We normalized this ratio by comparison of a simulated radio halo of a Coma like cluster in our sample to the observed radio halo of the Coma cluster. With the further simplifying assumption of an unique power law proton energy spectrum, the resulting radio halo, gamma ray and neutrino emission could be calculated. The main results are:

1. A cosmic ray proton to thermal gas energy ratio of 4 … 14 % $`(E_{\mathrm{p},\mathrm{min}}/\mathrm{GeV})^{0.375}`$ seems to be required to reproduce the total flux at 1.4 GHz of the Coma cluster with the simulated Coma like cluster. $`E_{\mathrm{p},\mathrm{min}}`$ is the lower kinetic energy cutoff of the protons with assumed spectral index $`\alpha _\mathrm{p}2.375`$.
2. The radial profile of the simulated Coma like cluster, is comparable to that of Coma.
3. But a detailed point by point radio to X-ray comparison, reveals a relation ($`L_{1.4\mathrm{GHz}}L_\mathrm{X}^b`$, $`b=1.26`$) which seems to be steeper than the observed relation ($`b1`$). This would imply that if radio halos are of hadronic origin, the parent cosmic ray proton population needs to have a spatially wider energy density distribution than the thermal gas.
4. The expected radio polarization is of the order of $`10^6`$. This is mainly a property of the magnetic field configuration, and should also hold for other radio halo models.
5. The expected gamma ray fluxes are below present observational limits, but possibly in the range of future instruments.
6. The set of 10 simulated clusters reproduces the observed radio luminosity-temperature relation of clusters of galaxies surprisingly well.

###### Acknowledgements.

We thank Bruno Deiss, Wilhelm Reich, Harald Lesch, and Richard Wielebinski for the access to the radio map of Coma. We also want to thank Haida Liang for the access to her collected data of the radio/x-ray luminosity relation, and Henk Spruit and Sabine Schindler for carefully reading the manuscript. Finally, we like to acknowledge many helpful comments by Gabriele Giovannini, the referee.