text:
# A Quantum Measurement Scenario which Requires Exponential Classical Communication for Simulation

## 1 . Introduction

Quantifying the power of shared entanglement is a fundamental goal of quantum information theory. In this paper, we compare the power of shared entanglement to the power of classical correlations with communication by considering the following restricted type of quantum protocol, which we call sampling-using-entanglement: Alice and Bob are two parties which are not allowed to communicate with each other, but they share $`k`$ entangled quantum bits. Each party makes a generalized measurement on his quantum bits.

We compare this protocol to a classical simulation of it, where Alice and Bob use classical correlations while communicating classical bits to sample from a distribution that is $`ϵ`$-close of the distribution resulting from the quantum protocol.

The cost of the classical simulation is the sum of the number of shared random bits and the number of classical bits communicated between the two parties. The number of shared random bits is the classical analogue to the amount of shared entanglement, while classical communication serves to augment the limited power of classical correlations.

Our main result is that $`2^k`$ classical bits are needed to simulate the distribution that arises from the quantum protocol’s measurement on $`k`$ entangled qubits. Our result is based on the analysis of sampling protocols for the set-disjointness problem DISJ(x,y) \[ASTS+98\]. In that paper, the authors demonstrate that the set-disjointness problem can be sampled using $`\mathrm{O}\left(\mathrm{log}(n)\mathrm{log}(1/ϵ)\right)`$ qubits of communication between two parties, while classical sampling of set-disjointness requires $`\mathrm{\Theta }(\sqrt{n})`$ classical bits, whether those bits are communicated on-line during the protocol, or are a cache of prior correlations. We cast these results in a different light, showing how this result holds for a protocol in which the parties share entangled bits. In particular, we demonstrate an efficient algorithm for sampling the set-disjointness function DISJ(x,y) using $`O\left(\mathrm{log}(n)\mathrm{log}(1/ϵ)\right)`$ entangled bits and no other communication, which is again an exponential improvement over the classical sampling complexity of set-disjointness. These results show that there exist problems for which classical simulation of an entangled state would require an exponential amount of communication, demonstrating that classical correlations are less powerful than quantum correlations.

Another intruiging element of this result is that in this quantum protocol each party maintains some secret information about their respective subset, so that parties do not have full information about the other parties’ output. Such a protocol is fundamentally different than known classical algorithms, where cryptographic schemes are necessary for $`A`$ and $`B`$ to stay ignorant about each other’s outputs.

Within the community of quantum communication complexity research, there is growing body of work addressing the cost of simulation of entanglement using classical bits\[BCT99\], \[MBCC01\], \[CP01\], \[Ste99\], \[CGM00\], and similar work in the communication complexity of entanglement \[BCvD97\]. The work of Brassard et. al. \[BCT99\] showed that 8 classical bits of communication suffice to exactly simulate a Von Neumann measurement on a Bell pair if an infinite random string is shared between two parties. That paper also shows a partial function for which an exact classical simulation of an entangled system of $`n`$ qubits would require $`2^n`$ classical bits of communication. However, they left open the question of whether any classical simulation existed on a more robust quantum measurement scenario when that simulation allowed small error, or was limited in the amount of bits available in prior correlation. Massar et. al. \[MBCC01\] showed that no protocol restricted to a finite amount of prior correlation could simulate the Bell pair, but that Bell states can be simulated without an infinite random string if the expected classical communication is finite.

Our work differs from the models studied above because we consider the the communication complexity of sampling rather than of functions or partial functions. Further, we demonstrate a protocol that is robust to small errors. Lastly, as shown in Massar et. al.\[MBCC01\] an infinite amount of shared randomness is too powerful to be used as a resource. As such, we quantify the total amount of classical information that is used in the protocol, counting the amount of prior correlations as well as the amount of communicated bits. Finally, we believe that the sampling model is more appropriate for investigation of quantum communication complexity problems since measurement of a quantum mechanical system is inherently a sampling process.

In 2 we discuss classical and quantum sampling protocols and introduce our restricted quantum protocol for sampling-using-entanglement. In 4 we state the problem of set-disjointness, and introduce an sampling-with-entanglement protocol for it. Section 5 contains the main theorem proof of the upper bound on the communication complexity of sampling-using-entanglement for the set-disjointness problem, while 8 presents the proof of the main lemma for this result. In 6 we present the classical sampling results for set-disjointness, demonstrating an exponential gap between the communication complexity of sampling-using-entanglement and communication complexity of classical sampling.

## 2 . Preliminaries

Consider a scenario in which two parties, Alice and Bob wish to output samples according to some known joint probability distribution. They share $`|\psi _{AB}`$, a bipartite entangled state of $`k`$ quantum bits, but they cannot communicate classically with each other. A simple protocol, labelled $`\{U_A,U_B,|\psi _{AB}\},`$ samples a probability distribution as follows: each party performs a unitary operation, $`U_A`$ and $`U_B`$, respectively, on their portions of the entangled system $`|\psi _{AB},`$ and additional ancilla. Then, each party performs a measurement in the computational basis and outputs the results. Since a measurement of a quantum system results in a classical probability distribution, this protocol samples from the classical probability distribution $`\mathrm{\Pi }\{U_A,U_B,|\psi _{AB}\}`$ induced by performing a measurement $`U_AU_B`$ on $`|\psi _{AB}`$. For shorthand, we will refer to a unitary operation and subsequent measurement in the computational basis as a generalized measurement, and say that the protocol samples a state $`\left(U_AU_B\right)|\psi _{AB}`$ to mean that the protocol outputs samples according to the classical probability distribution $`\mathrm{\Pi }\{U_A,U_B,|\psi _{AB}\}.`$

Now imagine Alice and Bob collude; that is, they claim to share an entangled state, but in reality they have only classical communication and prior classical correlations available to them. They will again be given unitary operators and asked to output measurements according to $`\mathrm{\Pi }\{U_A,U_B,|\psi _{AB}\}`$. In order for them to succeed in convincing a third party that they share entanglement, they must “simulate” this quantum protocol by outputting results according to a probability distribution that is $`ϵ`$-close to the probability distribution $`\mathrm{\Pi }\{U_A,U_B,|\psi _{AB}\}`$.

We now define formally a two-party quantum protocol for the sampling-using-entanglement of $`|\psi _{AB}`$ and the subsequent classical simulation of such a protocol. To do so, we include an intermediate definition of q-sampling \[ASTS+98\] where the protocol that samples $`|\psi _{AB}`$ not only outputs values according to the distribution of $`|\psi _{AB}`$ but creates a pure state that is close to the $`|\psi _{AB}`$ desired.

###### Definition 1

Let $`|\psi _{AB}`$ be a bipartite entangled state, and let $`U_A`$ and $`U_B`$ be unitary operators. A two-party sampling-using-entanglement protocol $`Q`$ is a quantum protocol executed by Alice and Bob as follows: Alice and Bob perform their respective unitary operators on their portions of the $`|\psi _{AB}`$ and some additional ancilla. Then, each party performs a generalized measurement in the computational basis. The two parties do not communicate in any other way. This protocol outputs samples according to the classical probability distribution $`\mathrm{\Pi }\{\mathrm{𝒰}_\mathrm{𝒜}\mathrm{}\mathrm{𝒰}_{\mathrm{}}\mathrm{}\mathrm{}\psi _\mathrm{𝒜}\mathrm{}\}`$, defined as probability distribution induced by performing a measurement $`U_AU_B`$ on $`|\psi _{AB}`$. Further, we say that $`Q`$ q-samples a state $`|\mathrm{\Psi }`$ if $`\mathrm{\Psi }|(U_AU_B)|\psi _{AB}>1ϵ`$. The communication complexity of sampling-using-entanglement is the number of quantum bits shared between the two parties.

###### Definition 2

Let $`Q`$ be a sampling-using-entanglement protocol as above that samples the distribution $`\mathrm{\Pi }\{\mathrm{𝒰}_\mathrm{𝒜}\mathrm{}\mathrm{𝒰}_{\mathrm{}}\mathrm{}\mathrm{}\psi _\mathrm{𝒜}\mathrm{}\}`$. Let $`P`$ be a two-party classical protocol $`P`$ where each party is computationally unbounded. We say $`P`$ is a classical simulation of $`Q`$ if $`P`$ outputs samples $`x`$ and $`y`$ according to a distribution $`\mathrm{\Pi }^{}\mathrm{}\mathrm{𝒫}\mathrm{}`$, such that that $`|\mathrm{\Pi }^{}\mathrm{}\mathrm{𝒫}\mathrm{}\mathrm{}\mathrm{\Pi }\{\mathrm{𝒰}_\mathrm{𝒜}\mathrm{}\mathrm{𝒰}_{\mathrm{}}\mathrm{}\mathrm{}\psi \}\mathrm{}_{\mathrm{𝒯}\mathrm{𝒱}\mathrm{𝒟}}\mathrm{}ϵ`$, where the norm is the Total Variational Distance. The communication complexity of $`P`$ is $`k+m`$, the sum of the number of classical bits of communication $`(k)`$ and the number of classical bits of correlation $`(m)`$.

The classical simulation does not need to simulate the specific quantum state or its overall phase information; it only matches the overall probability distribution. However, the quantum sampling protocol we describe is more restrictive, and matches the appropriate quantum state, as well as the resulting distribution.

In previous literature \[ASTS+98\], the communication complexity of classical sampling is defined as the number of classical bits of communication needed to sample a function. This definition is equivalent to the number of bits of correlation needed to sample $`f`$.

## 3 . Discussion

The main result of this paper is to answer the question of whether all quantum measurement scenarios can be simulated with only a polynomial number of classical bits. We demonstrate that the answer is no, by showing that already established results in quantum communication complexity can be applied directly to this problem and provide the answer.

The main technical theorems of this paper were first proved by Ambainis et. al. \[ASTS+98\] in a paper comparing the classical and quantum communication complexities of sampling. While we provide the details of those theorems again for completeness, our main result is to show that those results apply directly to a new communication protocol of sampling-using-entanglement. A sampling-using-entanglement protocol is equivalent to a quantum measurement scenario: both protocols involve each party performing a generalized measurement on their portion of the entangled state.

Prior work in quantum information theory has failed to answer the question of whether all larger non-separable entangled systems could be simulated with even a polynomial number of bits. However, by applying the communication complexity results for sampling, we find that if such a scenario could be simulated with less than an exponential number of bits relative to the quantum scheme, then we could beat the classical sampling bound proven in Ambainis et. al. Therefore, any classical simulation of such a quantum measurement scenario must match the classical communication lower bound on sampling. We use the results of Ambainis et. al. to demonstrate a quantum measurement scenario for the set-disjointness problem on $`n`$ entangled qubits cannot be simulated with less than $`2^n`$ classical bits of communication or correlation.

## 4 . Results for sampling of set-disjointness

We now show that there exists a sampling-using-entanglement quantum protocol for the set-disjointness problem that uses only $`O(\mathrm{log}n\mathrm{log}1/ϵ)`$ qubits. This result provides an exponential speedup over the classical sampling result for set-disjointness, where $`\mathrm{\Theta }(\sqrt{n})`$ classical bits of communication and shared randomness are required. The treatment provided here is based on the treatment in Ambainis et. al, but differs greatly in its organization and presentation.

Formally, the problem of sampling set-disjointness (for size $`\sqrt{n}`$) is defined as follows: Given a set $`\mathrm{\Omega }=\{1,\mathrm{},n\}`$, two parties $`A`$ and $`B`$ sample the set-disjointness function if they output uniformly random disjoint subsets $`S`$ and $`T`$ respectively, $`S,T\mathrm{\Omega }`$, $`|S|=|T|=\mathrm{\Theta }(\sqrt{n})`$. Let $`|\chi `$ be the superposition representing all possible disjoints subsets of size $`\sqrt{n}`$: $`|\chi =_{i,j\mathrm{}}|i|j`$. A protocol can quantum sample $`|\chi `$ if it creates a state $`|\varphi `$ which is close to $`|\chi `$. The trivial sampling-using-entanglement protocol quantum samples this function: it simply measures $`|\chi `$ in the computational basis. However, the number of qubits necessary to express $`|\chi `$ is $`\mathrm{log}\left(\genfrac{}{}{0pt}{}{n}{\sqrt{n}}\right)\sqrt{n}\mathrm{log}n`$, which does not improve upon the classical sampling complexity of this problem. Instead, we will show that there exists a smaller bipartite system that contains nearly all of the necessary information about $`|\chi `$. Specifically, we will construct a pure state $`|\psi `$ of only $`\mathrm{O}(\mathrm{log}(n)\mathrm{log}(1/ϵ))`$ entangled qubits which when operated upon by local unitary operations $`U_A`$ and $`U_B`$ yields a state which is an $`ϵ`$-approximation for $`|\chi `$.

###### Theorem 1

There exists a basis $`V`$ for $`|\chi `$, the superposition ranging over all possible disjoint subsets of size $`\sqrt{n}`$ with the following property: L et $`|\psi `$ be the projection of $`|\chi `$ onto the subspace spanned by the largest $`n^{\mathrm{log}1/ϵ}`$ eigenvectors of $`V`$. Then $`\chi |\left(VV\right)|\psi >1ϵ`$.

## 5 . Proof of Main Theorem

Proof:  We begin by using a matrix representation for states and the unitary operators. Consider $`|\chi `$ as above, written as

$$|\chi =\underset{i,j}{}m_{ij}|i|j\text{ where }m_{ij}=\{\begin{array}{cc}\hfill 1:& ij=\mathrm{},\hfill \\ \hfill 0:& ij\mathrm{}\hfill \end{array}$$

Let $`M_\chi `$ be the linear operator corresponding to the matrix representation for $`|\chi `$: $`M_\chi =\left[m_{ij}\right]`$, where $`M_\chi `$’s rows are indexed by $`|i`$ and columns by $`|j`$.

###### Remark 1

Given a state $`|\varphi =_{i,j}m_{i,j}|i|j`$, a matrix representation $`M_\varphi =[m_{i,j}]`$, and unitary operators $`U_A`$ and $`U_B`$, the matrix representation for the state $`\left(U_AU_B\right)|\varphi `$ is  $`U_AM_\varphi U_B^{}`$. 111To see the proof of this remark, please refer to the Appendix.

Using the representation, a rotation of a state $`|\psi `$ into the $`V`$ basis can be written as $`VM_\psi V^{}`$. This notation is particularly useful because it allows us to compare the norms of pure states by comparing the trace norms of their matrix representations. We define the trace norm of a matrix $`M`$ as $`M=\sqrt{\mathrm{Tr}\left(\mathrm{M}^{}\mathrm{M}\right)}`$.

###### Remark 2

Given states $`|\psi `$ and $`|\chi `$, their corresponding matrix representations $`M_\psi `$ and $`M_\chi `$, and a basis $`V`$, then for every $`ϵ`$, if $`VM_\psi V^{}M_\chi ^2<2ϵ`$, then $`\chi |\left(VV\right)|\psi >1ϵ`$.  222To see the proof of this remark, please refer to the Appendix.

With the above two facts in place, it remains to prove that there exists an appropriately-sized basis for the matrices representing our bipartite quantum states.

###### Lemma 1

There exists a basis $`V=\{v_1,\mathrm{}v_n\}`$ for $`M_\chi `$ such that the linear operator $`M_\psi `$, the projection of $`M_\chi `$ onto the subspace spanned by $`n^{\mathrm{log}(1/ϵ)}`$ eigenvectors of $`V`$ corresponding to the largest $`n^{\mathrm{log}(1/ϵ)}`$ eigenvalues, is within $`2ϵ`$ of $`M_\chi `$:

$$VM_\psi V^{}M_\chi ^2<2ϵ$$

The proof then follows.

Intuitively, this proof follows from the fact that we need only a small number of basis vectors to closely approximate any vector in the vector space of the linear operator $`M_\chi `$. We can then rotate that system of basis vectors back into a larger space and achieve a linear operator that is close to $`M_\chi `$. In subsequent sections we prove the associated lemma.

## 6 . Classical Bound comparison

We now return to the issue of simulation of quantum measurement. Since the sampling-using-entanglement protocol is equivalent to a quantum measurement scenario that samples set-disjointness, any classical simulation of this quantum measurement scenario would need to classically sample set-disjointness. But by Ambainis et. al. \[ASTS+98\], the lower bound for classical sampling of set-disjointness is $`\mathrm{\Omega }(\sqrt{n})`$. Therefore, any classical protocol simulating the the quantum measurement scenario of set-disjointness must use $`\mathrm{\Omega }(\sqrt{n})`$ bits of communication. This demonstrates an exponential gap between sampling-using-entanglement and classical sampling, and shows that the result of Brassard et. al. \[BCT99\] that a Von Neumann measurement on a Bell pair could be simulated with a constant number of classical bits of communication does not scale to larger quantum measurement scenarios on systems of qubits.

## 7 . A cryptographic protocol

This protocol shows that given a small initially entangled state, by performing unitary operations on that state and some ancilla, we can rotate into a larger space. But because the original entangled state resided on a smaller subspace, parties $`A`$ and $`B`$ do not have full knowledge of each other’s final output. Indeed, the are only aware of $`O(\mathrm{log}(n)\mathrm{log}(1/ϵ))`$ bits of information of each others output, rather than $`\sqrt{n}`$ bits. This is a fundamental difference over classical information theory, where the existence of cryptographic one-way functions is necessary to achieve the same “hidden” bits in a sampling protocol.

We believe that this amount of hidden information could be used to create future cryptographic protocols for property testing. Those result will be shown in a later version of this paper. We also believe that this imbalance between hidden information and classical simulation results is an interesting area for future research into the power of entanglement. This idea addresses some of the concerns expressed in Collins and Popescu \[CP01\], who raised the question of what particular behavior of entanglement is quantum-mechanical. We believe that the hidden bits here may directly answer this question.

## 8 . Proof of Lemma 1

Proof:  We begin by assuming that such a basis exists, and consider the $`t`$ eigenvectors that correspond to the largest $`t`$ eigenvalues. We derive conditions on the value of $`t`$ that succeeds in achieving an approximation bound of $`2ϵ`$.

Let $`V`$ be a basis of eigenvectors for $`M_\chi `$: $`M_\chi =V\mathrm{\Lambda }V^{}`$, where $`\mathrm{\Lambda }`$ is the matrix of eigenvalues. Let $`V_t`$ denote the subspace defined by the first $`t`$ eigenvectors of $`V`$. We define $`M^{}`$ as the linear operator created by projecting $`M`$ onto the subspace of the first $`k`$ eigenvectors in $`V`$.

$`|VM^{}V^{}M_\chi ^2`$ $`=`$ $`|VM^{}V^{}V\mathrm{\Lambda }V^{}^2`$
$`=`$ $`V(M^{}\mathrm{\Lambda })V^{}^2`$
$`=`$ $`M^{}\mathrm{\Lambda }^2`$
$`=`$ $`V_tM_\chi V_t^{}\mathrm{\Lambda }^2`$
$`=`$ $`{\displaystyle \underset{i>t}{\overset{|V|}{}}}\lambda _{i}^{}{}_{}{}^{2}`$
Therefore, $`VM^{}V^{}M_\chi ^2<2ϵ`$ if $`_{i>t}^{|V|}\lambda _{i}^{}{}_{}{}^{2}<2ϵ`$. We continue by characterizing the eigenvalues of $`M_\chi `$, and use this information to derive a bound on $`t`$ to ensure $`_{i>t}^{|V|}\lambda _{i}^{}{}_{}{}^{2}<2ϵ`$.

$`M_\chi `$ can be written as a linear combination of the ones matrix and a related matrix $`B=[b_{ij}]`$ of size $`N=\left(\genfrac{}{}{0pt}{}{n}{k}\right)`$ by $`N=\left(\genfrac{}{}{0pt}{}{n}{k}\right)`$ , whose rows are indexed by $`|i`$, columns by $`|j`$, and $`b_{ij}=1\text{ if }ij=\mathrm{},1\text{ if }ij\mathrm{}`$. Since $`M_\chi `$ is a linear combination of $`B`$ and the ones matrix, $`B`$ has the same eigenvectors as $`M_\chi `$.

$`B`$ is real and symmetric, and therefore has an orthonormal set of eigenvectors $`W=\{w_1,\mathrm{},w_{N1}\}`$. Using a result of Lovasz for $`B`$ \[Lov79\], we characterize the eigenspaces of $`B`$ as follows:

###### Lemma 2

\[Lov79\] $`B`$ has $`k+1`$ eigenspaces. Eigenspace $`E_0`$ is of dimension $`1`$ and contains the all $`1^{}s`$ vector. $`E_i`$ has dimension $`\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)`$. The typical eigenvector in $`E_i`$ is indexed by $`x_1,\mathrm{},x_{2i1},x_{2i}\{1,\mathrm{},n\}`$. The corresponding eigenvector $`e`$ (unnormalized) is given by: $`e_S=0`$ if there is an index $`j:|S\{x_{2j1},x_{2j}\}|1`$, otherwise $`e_S=\mathrm{\Pi }_j(1)^{|Sx_{2j}\}|}`$. The corresponding eigenvalues are $`\lambda _0=\frac{2\left(\genfrac{}{}{0pt}{}{nk}{k}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}`$, and $`\lambda _i=\frac{2\left(\genfrac{}{}{0pt}{}{nki}{ki}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}`$.

The first $`g`$ eigenspaces contain less than $`g\left(\genfrac{}{}{0pt}{}{n}{g}\right)n^{g+1}`$ eigenvectors. Setting $`t=n^{g+1}`$, we wish to show that we can pick $`g`$ small enough to satisfy $`_{i>t}^{|V|}\lambda _i^2<2ϵ`$.

A projection of a row of $`B`$ onto subspace $`E_i`$ is the sum of the row of $`B`$ onto each vector $`w`$ spanning $`E_i`$. Let $`q_i`$ be the square of the length of the projection of a row $`b`$ onto the subspace $`E_i`$. Since the value of such a projection is $`\lambda _i`$, then $`q_i=_{wE_i}\lambda _{i}^{}{}_{}{}^{2}`$. Therefore,

$$\underset{i>t}{\overset{|V|}{}}\lambda _i^2\underset{i>g}{\overset{k+1}{}}\underset{wE_i}{}\lambda _i^2=\underset{i>g}{\overset{k}{}}q_i.$$

Since each $`E_i`$ has dimension $`\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)`$, $`E_i`$ has $`\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)`$ eigenvectors, each with value $`\lambda _i=\frac{2\left(\genfrac{}{}{0pt}{}{nki}{ki}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}`$. Therefore,

$$q_i=\left(\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)\right)\left(\frac{2\left(\genfrac{}{}{0pt}{}{nki}{ki}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}\right)^2.$$

Next, we show that the projections decay rapidly, by showing the relative change between $`q_i`$ and $`q_{i+1}`$.

$`{\displaystyle \frac{q_{i+1}}{q_i}}`$ $`=`$ $`{\displaystyle \frac{\left(\left(\genfrac{}{}{0pt}{}{n}{\left(i+1\right)}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)\right)\left(\frac{2\left(\genfrac{}{}{0pt}{}{nk\left(i+1\right)}{k\left(i+1\right)}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}\right)^2}{\left(\left(\genfrac{}{}{0pt}{}{n}{i}\right)\left(\genfrac{}{}{0pt}{}{n}{i1}\right)\right)\left(\frac{2\left(\genfrac{}{}{0pt}{}{nki}{ki}\right)}{\left(\genfrac{}{}{0pt}{}{n}{k}\right)}\right)^2}}`$
$``$ $`{\displaystyle \frac{2n}{i+1}}{\displaystyle \frac{\left(ki\right)^2}{\left(nki\right)^2}}`$
$``$ $`{\displaystyle \frac{2n}{i+1}}\left({\displaystyle \frac{2k}{n}}\right)^2`$
If $`k=\mathrm{\Theta }(\sqrt{n})`$, then $`\frac{q_{i+1}}{q_i}=O(\frac{1}{i+1})`$. Therefore, $`q_g\frac{c^g}{g}`$.

Choosing $`g`$ to be $`k\mathrm{log}n\mathrm{log}(1/ϵ)/\mathrm{log}\mathrm{log}(1/ϵ))`$,

$$\underset{i>g}{\overset{k}{}}q_i=\underset{i>g}{\overset{k}{}}c^i/i<c^g/g!<ϵ.$$

Since $`_{i>t}^{|V|}\lambda _i^2_{i>g}^{k+1}q_i`$,

$$\underset{i>t}{\overset{|V|}{}}\lambda _i^2ϵ.$$

With $`g=O(\mathrm{log}(1/ϵ)/\mathrm{log}\mathrm{log}(1/ϵ))`$, the total number of eigenvectors needed is bounded above by $`n^{g+1}=O(n^{\mathrm{log}1/ϵ})`$, and therefore only $`O(\mathrm{log}n\mathrm{log}1/ϵ)`$ entangled qubits are needed to perform the protocol. This completes the proof.

## 9 . Acknowledgements

The author wishes to thank Umesh Vazirani, Richard Josza, and Jennifer Sokol for numerous helpful discussions. This research is supported by a Siebel Scholars Fellowship and DARPA Grant F30602-00-2-060.

## Appendix A . Proof of Remark 1

We prove that given a bipartite state $`|\psi `$ and local unitary operations $`U`$ and $`V`$, we can write the result of $`\left(UV\right)|\psi `$ as $`UMV^{}`$ where $`M`$ is the matrix representing $`|\psi `$.

Proof:

A and B want to perform $`UV`$ on a bipartite state $`|\varphi =_{ij}\varphi _{ij}|i|j`$. The coefficient $`\varphi _{ij}`$ refers to the probability amplitude that a measurement of the state $`|\varphi `$ will yield $`|i`$ to party $`A`$ and $`|j`$ to party $`B`$. Performing the unitary operation $`UV`$ on $`|\varphi `$ changes the coefficients to $`|\varphi ^{}`$ according to the rule

$`\varphi _{ij}^{}`$ $`=`$ $`(i\text{th row of }Uj\text{th row of }V|\varphi `$
$`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \underset{m}{}}U_{ik}V_{jm}\varphi _{km}`$
$`=`$ $`{\displaystyle \underset{k}{}}U_{ik}{\displaystyle \underset{m}{}}V_{jm}\varphi _{km}`$
$`=`$ $`{\displaystyle \underset{k}{}}U_{ik}{\displaystyle \underset{m}{}}\varphi _{km}V_{jm}`$
Now, we show that in the matrix formulation, $`UV|\varphi =UMV^{}`$. Beginning with a matrix $`M`$, we show the values of the matrix $`UMV^{}`$.

$`(UM)_{ij}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{n}{}}}U_{ik}M_{kj}`$
$`(AV^T)_{ij}`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{n}{}}}A_{im}V_{jm}`$
$`(UAV^T)_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}U_{ik}\left({\displaystyle \underset{m}{}}A_{km}V_{jm}\right)`$
$`(UMV^T)_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}U_{ik}\left({\displaystyle \underset{m}{}}M_{km}V_{jm}\right)`$
Hence, these two formulations are equal.

## Appendix B . Proof of Remark 2

Proof:

Note that $`\mathrm{Tr}\left(\mathrm{M}_\chi ^{}\mathrm{M}_\chi \right)=\chi |\chi `$, since $`\mathrm{Tr}\left(\mathrm{M}_\chi \right)^{}\mathrm{M}_\chi =_{\mathrm{ij}}\left(\mathrm{M}_\chi \right)_{\mathrm{ij}}.`$ But since this matrix consists of all the elements of $`|\chi `$, this result is $`_{ij}\left(M_{\chi _{ij}}\right)=\chi |\chi .`$

We rewrite $`VM_\psi V^{}M_\chi ^2`$ in terms of the trace norm:

$`VM_\psi V^{}M_\chi ^2=`$
$`\mathrm{Tr}\left(\left(\mathrm{VM}_\psi \mathrm{V}^{}\right)^{}\left(\mathrm{VM}_\psi \mathrm{V}^{}\right)\right)`$
$`2\mathrm{T}\mathrm{r}\left(\mathrm{M}_\chi ^{}\mathrm{VM}_\psi \mathrm{V}^{}\right)+\mathrm{Tr}\left(\mathrm{M}_{\chi }^{}{}_{}{}^{}\mathrm{M}_\chi \right)`$
Since the trace is invariant under basis transformations, and $`V^{}V=I`$,

$`\mathrm{Tr}\left(\left(\mathrm{VM}_\psi \mathrm{V}^{}\right)\right)^{}\left(\mathrm{VM}_\psi \mathrm{V}^{}\right)`$ $`=`$ $`\mathrm{Tr}\left(\left(\mathrm{VM}_\psi ^{}\mathrm{V}^{}\mathrm{VM}_\psi \mathrm{V}^{}\right)\right)`$
$`=`$ $`\mathrm{Tr}\left(\left(\mathrm{M}_\psi ^{}\mathrm{M}_\psi \right)\right)`$
$`=`$ $`{\displaystyle \underset{ij}{}}\left(M_{\psi _{ij}}\right)=\psi |\psi `$
Continuing, we see that

$`2\mathrm{T}\mathrm{r}\left(\mathrm{VM}_\psi ^{}\mathrm{V}^{}\mathrm{M}_\chi \right)`$ $`=`$ $`2\mathrm{T}\mathrm{r}\left(\mathrm{VM}_\psi ^{}\mathrm{V}^{}\mathrm{V}\mathrm{\Lambda }\mathrm{V}^{}\right)`$
$`=`$ $`2\mathrm{T}\mathrm{r}\left(\mathrm{M}_\psi ^{}\mathrm{\Lambda }\right)`$
$`=`$ $`2\psi ,000||\chi `$
where $`000|`$ represents an ancilla, so that the $`|\chi `$ and $`|\psi ,000`$ lie in the same vector space.

Therefore,

$`VM_\psi V^{}M_\chi ^2`$ $`=`$ $`\psi |\psi +\chi |\chi 2\psi ,000|\chi `$
$`=`$ $`22\psi ,\mathrm{\hspace{0.17em}000}|\chi `$
$`<`$ $`2ϵ`$
Therefore, $`\chi |\psi ,\mathrm{\hspace{0.17em}000}>1ϵ`$. Since unitary operations on pure states cannot change their length, $`\chi |\psi ,\mathrm{\hspace{0.17em}000}=\chi |VV|\psi ,\mathrm{\hspace{0.17em}000}>1ϵ`$.


format: text

text:
# Untitled Document

ABOUT THE SIMULTANEOUS CO-EXISTENCE OF INSTANTANEOUS AND RETARDED INTERACTIONS IN CLASSICAL ELECTRODYNAMICS

LjUBA ŠKOVRLJ and TOMISLAV IVEZIĆ†,‡

Faculty of Textile Technology, University of Zagreb, Pierottieva 6, 10000

Zagreb, Croatia

Rud er Bošković Institute, P.O.B. 180, 10002 Zagreb, Croatia

E-mail: ivezic@rudjer.irb.hr

In this paper it is proved that, contrary to the results found by A.E. Chubykalo and S.J. Vlaev (Int. J. Mod. Phys. 14, 3789 (1999)), the retarded electric and magnetic fields for an uniformly accelerated charge exactly satisfy Maxwell equations (ME). Furthermore it is shown that ME are correctly written in the usual form with the partial derivatives and thus not, as proposed by Chubykalo and Vlaev, with the total derivatives.

1. Introduction

Recently Chubykalo and Vlaev1 argued that the electromagnetic fields (obtained from the Liénard-Wiechert potentials (LW)) of a charge moving with a constant acceleration along an axis do not satisfy Maxwell equations (ME) if one considers exclusively a retarded interaction. According to their interpretation the retarded interaction is mathematically described by the implicit functional dependence of the electric $`𝐄`$ and magnetic $`𝐁`$ fields on the observation time $`t.`$ Further they stated that ME will be satisfied ”if and only if” one takes into account an *explicit* functional dependence of $`𝐄`$ and $`𝐁`$ on the observation time $`t`$ *together* with the mentioned implicit dependence. This explicit dependence of fields on $`t`$ then mathematically describes an instantaneous interaction. The same assertions were also reported in series of papers, e.g., Refs. 2-4.

In the paper5 we have argued that the ”electrodynamics dualism concept”4 with the simultaneous coexistence of instantaneous and short range interaction is not correct. In this paper we show that the opinion presented in1-4 about the incorrectnes of the retarded solutions emerges, in fact, from an incorrect mathematical procedure in treating the partial derivatives of continuos functions of several variables. Their1-4 physical interpretation with the ”electrodynamics dualism concept” is then a consequence of such an incorrect mathematical treatment. Further it will be proved that, contrary to the assertions from,1-4 the ME

$`𝐄`$ $`=`$ $`\rho /ϵ_0`$
$`𝐁`$ $`=`$ $`0`$
$`\times 𝐁`$ $`=`$ $`\mu _0𝐣+\mu _0\epsilon _0𝐄/t`$ (1)
$`\times 𝐄`$ $`=`$ $`𝐁/t`$
are correctly written in the usual form, Eqs. (1), i.e., with the partial derivatives and not with the total derivatives, and that the retarded fields of a charge $`q`$ moving with a constant acceleration along the $`X`$ axis, Eqs. (22)-(27) from Ref. 1, exactly satisfy ME (1).

2. Derivation of the Fields $`𝐄`$ and $`𝐁`$ Taking into Account the Retarded Interaction Only

First it will be discussed the calculation of the LW fields $`𝐄`$ and $`𝐁`$ from a charge $`q`$ for which the retarded position $`x_{0i}(t_0)`$ as a function of the retarded time $`t_0`$ is assumed given. These fields are

$`𝐄`$ $`=`$ $`Kq\left\{(𝐑𝐕R/c)(1V^2/c^2)/s^3\right\}_{t_0}+`$ (2)
$`+`$ $`Kq\left\{\left[𝐑\times \left((𝐑R𝐕/c)\times \dot{𝐕}\right)\right]/c^2s^3\right\}_{t_0},`$
$$𝐁=\left\{\left[𝐑\times 𝐄\right]/Rc\right\}_{t_0},$$
(3)
where $`K=1/4\pi \epsilon _0`$ and $`s=R\mathrm{𝐕𝐑}/c.`$ (The notation is similar to that one in Ref. 1; $`\left\{\mathrm{}\right\}_{t_0}`$ means that all functions of $`x,y,z,t`$ in the braces $`\left\{\right\}`$ are taken at the moment of time $`t_0,`$ that is determined from the retardation condition (6), see below.) $`𝐄`$ and $`𝐁`$ are determined in the usual way

$$𝐄=\phi 𝐀/t,𝐁=\times 𝐀$$
(4)
from the LW potentials

$$\phi (𝐫,t)=\left\{Kq/s\right\}_{t_0},𝐀(𝐫,t)=\left\{Kq𝐕/s\right\}_{t_0}.$$
(5)
We note that the field and source point variables are connected by the retardation condition

$$R(x_i,x_{0i})=\left[\underset{i}{}\left(x_ix_{0i}\right)^2\right]^{1/2}=c\left(tt_0\right).$$
(6)
This equation is Eq. (8) from Ref. 1, see also Ref. 7 Sec. 20-1, or Ref. 6 Secs. 62 and 63. Since $`x_{0i}`$ is asumed given as a function of $`t_0`$ we can write

$$R(x_i,x_{0i}(t_0))=f(x,y,z,t_0)=c\left(tt_0\right),$$
(7)
which is a functional relation between $`x_i`$, $`t`$ and $`t_0,`$ see7 Eq. (20-4). The LW potentials are expressed in terms of $`R`$ (6), which means that $`\phi (𝐫,t)`$ and $`𝐀(𝐫,t)`$ are not only implicit functions of $`x,y,z`$ through $`t_0`$ but *the explicit functions as well*. Although Landau and Lifshitz6 p. 161 stated that the LW potentials (5) are only implicit functions of $`x,y,z`$ through $`t_0`$ they really calculated the fields $`𝐄`$ and $`𝐁`$ taking into account that R, (6) or (7), is also the explicit function of the field coordinates $`x,y,z`$ in the same way as in Ref. 7. This can be clearly seen from the calculation of $`t_0`$ when $`R`$ is written as

$$R=ct_0=_1R+(R/t_0)t_0,$$
(8)
which is Eq. (20-9) from7 and the same is written in the unnumbered equation preceding the numbered equation (63.7) in Ref. 6. Then the equations which constitute the transformation of the differential operators from the coordinates of the field point to those of the radiator are given as

$$/t=(R/s)/t_0,=_1+(t_0)/t_0,$$
(9)
Eqs. (20-8) and (20-11) in Ref. 7. The differential operator $`_1`$ refers to the differentation with respect to, e.g., the first argument of the function $`f`$ in Eq. (7), that is, differentation at constant retarded time $`t_0`$. Expressing this statement in another way we can say that $`_1`$ denotes the differentation with respect to the arguments $`x,y,z`$ of the function $`f,`$ i.e., $`R,`$ (6) or (7), (and thus of the potentials $`\phi `$ and $`𝐀`$ and the fields $`𝐄`$ and $`𝐁`$ as well), on which $`f`$ (consequently $`\phi ,`$ $`𝐀,`$ $`𝐄`$ and $`𝐁`$) *explicitly depend* (not only implicitly through $`t_0`$). This consideration reveals that the relations (11) from Ref. 1 have to be correctly written as

$`\phi /x_i`$ $`=`$ $`(\phi /x_{i_1})+(\phi /t_0)(t_0/x_i),`$
$`𝐀/t`$ $`=`$ $`(𝐀/t_0)(t_0/t),`$ (10)
$`A_k/x_i`$ $`=`$ $`(A_k/x_{i_1})+(A_k/t_0)(t_0/x_i),`$
where $`/x_{i_1}`$ denotes, as above $`_1,`$ the differentation with respect to the argument $`x_i`$ (i.e., $`x`$ or $`y`$ or $`z`$) on which $`\phi `$ or $`A_k`$ explicitly depend. The partial derivative, e.g., $`/x_i`$ on the l.h.s. of (10) can be called - the *complete* partial derivative - which takes into account both the explicit and the implicit dependence of the functions $`\phi `$ or $`A_k`$ on the argument $`x_i`$. The partial derivative $`/x_{i_1}`$ in the first and the third equation on the r.h.s. of (10) takes into account only the *explicit* dependence of the functions $`\phi `$ or $`A_k`$ on the argument $`x_i.`$ For the sake of brevity we shall, from now on, call such type of the partial derivative - the *explicit* partial derivative.

When such procedure is applied to the calculation of $`𝐄`$ and $`𝐁`$ by Eq. (4) from the LW potentials $`\phi `$ and $`𝐀`$ (5) then one finds the exact expressions for the usual LW fields (2) and (3) without need for addition of any new terms. We note that in Ref. 1 Eqs. (11) the partial derivatives which refer to the explicit dependence of $`\phi `$ and $`𝐀`$ on $`x_i`$, i.e., $`(\phi /x_{i_1})`$ and $`(A_k/x_{i_1}),`$ were not taken into account. As a consequence the authors of1 had to add an additional term in their Eq. (21) to find the result (3).

3. The Retarded Electromagnetic Fields of a Charge Moving with a Constant Acceleration Do Satisfy Maxwell Equations

It is argued in Ref.1 that the electric and magnetic fields for a charge $`q`$ moving with a constant acceleration along the $`X`$ axis, their Eqs. (22)-(27), do not satisfy ME (1). These fields are the retarded fields (2) and (3) but written in components, and taking into account that the velocity and the acceleration of that charge $`q`$ have only $`x`$-components, i.e., $`𝐕(V,0,0)`$ and $`𝐚(a,0,0).`$ We shall, for brevity, quote here only $`E_x`$ and $`B_y`$ components, Eqs. (22) and (26) from Ref. 1,

$`E_x(x,y,z,t)`$ $`=`$ $`Kq\left\{\left((xx_0)VR/c\right)\left(1V^2/c^2\right)/\left(RV(xx_0)/c\right)^3\right\}_{t_0}`$ (11)
$`+Kq\left\{a\left((xx_0)^2R^2\right)/c^2\left(RV(xx_0)/c\right)^3\right\}_{t_0},`$
$`B_y(x,y,z,t)`$ $`=`$ $`Kq\left\{Vz\left(1V^2/c^2\right)/c^2\left(RV(xx_0)/c\right)^3\right\}_{t_0}`$ (12)
$`Kq\left\{aRz/c^3\left(RV(xx_0)/c\right)^3\right\}_{t_0}.`$
To prove the above mentioned assertion the authors of1 supposed, as in their Sec. 2., that $`𝐄`$ and $`𝐁`$ fields, their Eqs. (22)-(27), are functions of $`x,y,z,`$ $`t`$ only through $`t_0,`$ and consequently used the differentiation rules as in their Eqs. (11). But from our discussion in Sec. 2., see the relation for $`R`$ (7), and from the explicit expressions for $`𝐄`$ and $`𝐁`$ fields, (11) and (12) or Eqs. (22)-(27) from Ref.1, one unambiguously concludes that $`𝐄`$ and $`𝐁`$ fields are also the explicit functions of the field point $`x,y,z.`$ Thence the relations (28) and (30)-(32) from Ref.1 have to be changed in accordance with our discussion given in Sec. 2. Thus, e.g., the second relation in Eq. (28) in Ref. (1) has to be changed and written as

$$E\left\{orB\right\}_k/x_i=(E\left\{orB\right\}_k/x_{i_1})+(E\left\{orB\right\}_k/t_0)(t_0/x_i),$$
(13)
where the partial derivative on the l.h.s. of (13) is (in our terminology) - the *complete* partial derivative - which takes into account both the explicit and the implicit dependence of the functions $`E\left\{orB\right\}_k`$ on the argument $`x_i`$, while the partial derivative in the first term on the r.h.s. of (13) is - the *explicit* partial derivative, which takes into account only the explicit dependence of the functions $`E\left\{orB\right\}_k`$ on the argument $`x_i.`$ Of course the relations (30)-(32) from1 have to be changed in the same way as explained above.

Then one finds that the term $`_1\times 𝐄`$ (due to the explicit dependence of $`𝐄`$ on the arguments $`x,y,z`$) *exactly cancels* the terms on the r.h.s. of Eqs. (35) and (36) from Ref. 1. In such a way we prove, as expected, that the Faraday law is exactly satisfied. It can be easily shown that not only the Faraday law but all ME (1) are exactly satisfied by the LW fields (2) and (3), and thus also by the fields from uniformly accelerated charge, (11) and (12) or Eqs. (22)-(27) from Ref. 1, when all partial derivatives are correctly treated.

At this place it is worth to discuss the results from Sec. 4 in Ref. 1, where the similar mathematically incorrect step has been made. The authors of1 argued that in the calculations of $`𝐄`$ and $`𝐁`$ from (4) and in ME (1) as well all derivatives have to be considered as *total* ones. To explain the need for such change they calculate $`t_0/t`$ and $`t_0/x_i`$, Eqs. (14) in,1 using two different expressions for $`R`$, our Eq. (6), or Eqs. (37) and (38) from Ref. 1. Then to calculate the above mentioned expressions they claim that the commonly used partial derivatives $`R/t`$ and $`R/x_i`$ have to be replaced by the total derivative $`dR/dt`$ and $`dR/dx_i`$, see Eqs. (39) in Ref. 1. Thus they also objected the use of the partial derivatives by Landau and Lifshitz6 (in the calculation of $`R/t`$ and $`R/x_i`$) considering that Landau and Lifshitz6 used the symbol $``$: ”in order to emphasize that $`R`$ depends also on other independent variables $`x,y,z`$,$`\mathrm{}`$ .” The justification for such assertions they find in the following statement: ”The point is that if a given function is expressed by two different types of functional dependences, then exclusively *total* derivatives of these expressions with respect to a *given* variable can be equated (contrary to the *partial* ones).” However this statement is not true. Namely (using our terminology) it is true that one cannot equate - the *explicit* partial derivatives, but one can equate - the *complete* partial derivatives of the mentioned expressions. Moreover when one calculates, e.g., $`R/x,`$ then it has to be written as

$$R/x=(R/x_1)+(R/t_0)(t_0/x),$$
(14)
since according to (6) and (7) $`R`$ depends explicitly not only on $`x`$ but also on $`y`$ and $`z`$ (thence $`R/x_1`$). Further, we explain what is the most important for the use of the operator $`/x`$ on the l.h.s. of (14) (also (8), (9), (10) and (13)) and not the total derivative operator $`d/dx.`$ It is the fact that $`t_0`$ in $`R`$ depends not only on the $`x`$ variable but also on other independent variables $`y`$ and $`z`$. Of course $`R/x`$ on the l.h.s. of (14) is (in our terminology) - the *complete* partial derivative of $`R`$, while $`R/x_1`$ (the first term on the r.h.s. of (14)) is - the *explicit* partial derivative of $`R`$. To see that this explanation is correct and to understand the difference between the total derivative and the partial derivative of the function (composite) of several variables one can consult some mathematical book, e.g., Ref. 8 Eqs. (8) and (9) in Sec. 153. Also one can see the nice explanation given in Ref. 7 Sec. (20-1). This means that the objections from1 to the Landau and Lifshitz6 derivations and results are unfounded. Hence in Eqs. (39) from Ref. 1 the *complete* partial derivatives $`R/t`$ and $`R/x_i`$ have to be written and not the *total* derivatives $`dR/dt`$ and $`dR/dx_i`$. Obviously the *partial* derivatives have to be retained in all other expressions and they must not be replaced by the *total* derivatives. The proposition from Ref. 1 that partial derivatives have to be replaced by the total derivatives has been also extensively used for numerous conclusions in Refs. 2-4, see particularly Eqs. (43)-(46) and the discussion in Ref. 2. We see that such changes of ME are unnecessary and, in fact, incorrect.

4. Conclusion

The consideration presented in this paper unambiguously reveals that Maxwell equations (1) are correctly written in the usual way with the partial derivatives and not with the total derivatives, as argued in Ref. 1. Furthermore, as proved in our Sec. 3., the retarded Liénard-Wiechert fields for an uniformly accelerated charge, (11) and (12) or Eqs. (22)-(27) from Ref. 1, exactly satisfy the usual form of ME (1). Hence it is not true that: ”there is a *simultaneous* and *independent* coexistence of *instantaneous* and *retarded interactions* which cannot be reduced to each other.” The retarded solution is a full-value solution, i.e., it is a complete solution of ME, which is sufficient to describe the whole electromagnetic phenomenon. In addition, it has to be mentioned (see also Ref. 5) that it is possible to write the complete solution of ME in the present-time formulation, i.e., as an action-at-a-distance formulation, and such a solution is completely equivalent to the retarded (and advanced) solution. This is explicitly shown in the general case in Ref. 9 by means of Lagrange series expansion. We shall present the potentials and fields in the present-time formulation in a *closed form* for the case of an uniformly accelerated charge, but it will be reported elsewhere.

References

1. A.E. Chubykalo and S.J. Vlaev, Int. J. Mod. Phys. 14, 3789 (1999).

2. A.E. Chubykalo and R. Smirnov-Rueda, Mod. Phys. Lett. A12, 1 (1997).

3. A.E. Chubykalo, Mod. Phys. Lett. A13, 2139 (1998).

4. A.E. Chubykalo and R. Smirnov-Rueda, Phys. Rev. E53, 5373 (1996);

E55, 3793(E) (1997); E57, 3683(1998).

5. T. Ivezić and Lj. Škovrlj, Phys. Rev. E57, 3680 (1998).

6. L.D. Landau and E.M. Lifshitz, The Classical Theory of Field (Pergamon,

Oxford, 1975).

7. W.K.H. Panofsky and M. Phillips, Classical Electricity and Magnetism

(Addison-Wesley, Reading, MA, 1962).

8. V.I. Smirnov, Kurs Visšei Matematiki, Tom I (Nauka, Moscow, 1967),

(English translation: A Course of Higher Mathematics, Vol I (Pergamon,

Oxford, 1964)).

9. R.A. Villecco, Phys. Rev. E48, 4008 (1993).


text:
# The rms peculiar velocity of galaxy clusters for different cluster masses and radii

## 1 INTRODUCTION

One of the interesting unknowns in cosmology is the large-scale peculiar velocity field in the Universe. The peculiar velocity field can be studied by using galaxies or clusters of galaxies. However, there are some advantages in studying the peculiar velocity field by using galaxy clusters. One of these advantages comes from the fact that, on scales probed by galaxy clusters, velocity fluctuations are largely in the quasi-linear regime and close to the initial state from which large scale structures developed. In addition, peculiar velocities of clusters can be determined more accurately than peculiar velocities of galaxies since the distance to each cluster can be obtained from a large number of member galaxies, thus considerably reducing the velocity uncertainties of clusters. Cluster motions could therefore provide an important tool in probing the large-scale peculiar velocity field.

Peculiar velocities of clusters of galaxies have been studied in several papers (e.g. Bahcall, Gramann & Cen 1994; Lauer & Postman 1994; Bahcall & Oh 1996; Moscardini et al. 1996; Borgani et al. 1997; Watkins 1997; Dale et al. 1999; Hudson et al. 1999; Borgani et al. 2000; Colless et al. 2001). Watkins (1997) developed a likelihood method for estimating the rms peculiar velocity of clusters from line-of sight velocity measurements. This method was applied to two observed samples of cluster peculiar velocities: the SCI sample (Giovanelli et al. 1997) and a subsample of the Mark III catalogue (Willick et al. 1997). Watkins (1997) found that the rms one-dimensional cluster peculiar velocity is $`256_{75}^{+106}`$ km s-1, which corresponds to the three-dimensional rms velocity $`459_{130}^{+184}`$ km s-1 . Dale et al. (1999) obtained Tully-Fisher peculiar velocities for 52 Abell clusters distributed over the whole the sky between $`50`$ and $`200h^1`$ Mpc. They found that the rms one-dimensional cluster peculiar velocity is $`341\pm 93`$ km s-1, which corresponds to the three-dimensional rms velocity $`591\pm 161`$ km s-1.

Radial peculiar velocities of clusters can be determined to large distances by measuring the Sunyaev-Zeldovich (1980) (SZ) effect. Rephaeli & Lahav (1991) made one of the first estimates of the possibility of measuring the peculiar velocities by using the SZ effect for a selected sample of galaxy clusters. However, most convincing measurements for individual clusters have been done only recently, using the new generation of sensitive bolometers (Holzapfel et al. 1997; Lamarre et al. 1998). The accuracy of SZ measurements for determining peculiar velocities of clusters have been studied in several papers (Haehnelt & Tegmark 1996; Aghanim, Gorski & Puget 2001; Diego et al. 2002; Holder 2002; Nagai, Kravtsov & Kosowsky 2002). With $`\mu K`$ sensitivity on arcminute scales at several frequencies it will be possible to measure peculiar velocities to an accuracy of $``$130 km s-1.

In this paper we study the rms peculiar velocity of clusters, $`v_{rms}`$, in different cosmological models assuming that the initial density fluctuation field is a Gaussian field. To investigate the nonlinear regime, we use N-body simulations. We examine cluster peculiar velocities for different cluster masses. Do cluster velocities depend on their masses? The rms peculiar velocity of peaks in the initial Gaussian field does not depend on the height of peaks (Bardeen et al. 1986). However, as an initially Gaussian density field evolves gravitationally in the nonlinear regime, it becomes non-Gaussian.

The evolution of peculiar velocities of galaxy clusters in different N-body models with an initially Gaussian density field has been examined in several papers (e.g. Bahcall, Gramann, Cen 1994; Croft & Efstathiou 1994; Suhhonenko & Gramann 1999, Colberg et al 2000; Sheth & Diaferio 2001). Groth & Efstathiou (1994) studied the cluster peculiar velocities for two different cluster richnesses, described by the mean intercluster separations $`d_{cl}=30h^1`$ Mpc and $`d_{cl}=55h^1`$ Mpc. They found that the peculiar velocity distributions are almost independent of cluster richness. Sheth & Diaferio (2001) studied rms peculiar velocity of clusters for different masses and found that in N-body simulations the rms velocity of clusters depends weakly on cluster mass, with a small trend that for more massive clusters the rms peculiar velocity decreases. On the other hand, Suhhonenko & Gramann (1999) investigated the properties of clusters using N-body simulations and found that the rms peculiar velocities of clusters increase with cluster richness.

In this paper we study the dependence of $`v_{rms}`$ on cluster masses in more detail. We also examine the rms peculiar velocities of clusters for different cluster radii. We use the N-body simulations published by the Virgo Consortium and described in detail by Jenkins et al. (1998). These simulations were carried out using a parallel, adaptive particle-particle/particle-mesh (AP3M) code (Couchman, Thomas & Pearce 1995; Pearce & Couchman 1997). In this paper we analyze the velocities in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model (see Jenkins et al. (1998) and Section 2 for the description of the cosmological parameters in these models). We also follow the evolution of particles in a similar $`\mathrm{\Lambda }`$CDM model with the same cosmological parameters but using a particle-mesh (PM) code. In this PM code we use a traditional two-point finite-difference approximation to calculate the forces on the grid. We also examine the shifted-mesh scheme (see eq. (13) below) and show that this scheme artificially boosts the velocities of clusters.

To identify clusters in the simulations we use two methods: the standard friends-of-friends (FOF) method and the method, were the clusters are defined as the maxima of the density field smoothed on the scale $`R1h^1`$ Mpc (DENSMAX). To determine the velocity of DENSMAX clusters, we use the same smoothing scale as for the density field. The velocity of FOF clusters is defined to be the mean velocity of all the particles in the cluster.

This paper is organized as follows. In Section 2 we describe the cosmological models, N-body simulations and cluster selection algorithms used. We also examine the linear theory predictions for the peculiar velocities of peaks in the Gaussian field. In Section 3 we examine the rms peculiar velocity of clusters for different cluster masses and radii, and compare the cluster velocities with the linear theory predictions. In Section 4 we briefly discuss the effect of the shifted-mesh scheme on the cluster velocities. A summary and discussion are presented in Section 5.

## 2 MODELS

### 2.1 Simulations

We analyze peculiar velocities in N-body simulations carried out by the Virgo consortium for two cosmological models, $`\mathrm{\Lambda }`$CDM ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }=0.7`$) and $`\tau `$CDM ($`\mathrm{\Omega }_0=1)`$ as described by Jenkins et al. (1998). In these cold dark matter (CDM) models the power spectrum of the initial conditions was chosen to be in the form given by Bond & Efstathiou (1984),

$$P(k)=\frac{Ak}{[1+(aq+(bq)^{3/2}+(cq)^2)^\nu ]^{2/\nu }},$$
$`(1)`$
where $`q=k/\mathrm{\Gamma }`$, $`a=6.4h^1`$ Mpc, $`b=3h^1`$ Mpc, $`c=1.7h^1`$ Mpc, $`\nu =1.13`$ and $`\mathrm{\Gamma }=0.21`$. The normalization constant, $`A`$, was chosen by fixing the value of $`\sigma _8`$ (the linearly extrapolated mass fluctuation in spheres of radius $`8h^1`$ Mpc); $`\sigma _8=0.9`$ and $`\sigma _8=0.51`$ for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model, respectively.

We investigated the linear theory predictions for peculiar velocities of peaks in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model. The linear rms velocity fluctuation on a given scale $`R`$ can be expressed as

$$\sigma _v(R)=H_0f(\mathrm{\Omega }_0)\sigma _1(R),$$
$`(2)`$
where the spectral moments $`\sigma _j`$ are defined for any integer $`j`$ by

$$\sigma _j^2=\frac{1}{2\pi ^2}P(k)W^2(kR)k^{2j+2}𝑑k,$$
$`(3)`$
$`W(kR)`$ is a window function and $`f(\mathrm{\Omega }_0)`$ is the dimensionless growth rate. The function $`f(\mathrm{\Omega }_0)=0.51`$ and $`f(\mathrm{\Omega }_0)=1.0`$ in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model, respectively. (We note that the approximation $`f(\mathrm{\Omega }_0)=\mathrm{\Omega }_0^{0.6}`$ underestimates the dimensionless growth rate by $`5`$% in the flat $`\mathrm{\Omega }_0=0.3`$ model).

Bardeen et al. (1986) showed that in the linear approximation the rms peculiar velocity at peaks of the smoothed density field differs systematically from $`\sigma _v(R)`$, and can be expressed as

$$\sigma _p(R)=\sigma _v(R)\sqrt{1\sigma _0^4/\sigma _1^2\sigma _1^2}.$$
$`(4)`$
In this approximation, the rms velocities of peaks do not depend on the height of the peaks.

Fig. 1 shows the rms peculiar velocity of peaks, $`\sigma _p(R)`$, for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model. We have used the top-hat window function. For comparison, we show also the rms peculiar velocity $`\sigma _v(R)`$ for the same models. For the radius $`R=1h^1`$ Mpc, $`\sigma _p=509`$ km/s and $`\sigma _p=562`$ km/s in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, respectively. At this radius, $`\sigma _p`$ is lower than $`\sigma _v`$ about $`2`$ per cent for the models studied. On larger scales, the difference between $`\sigma _p`$ and $`\sigma _v`$ increases.

Next we study peculiar velocities in N-body simulations. The Virgo simulations were created using a parallel adaptive particle-particle/particle-mesh (AP3M) code as described by Couchman, Thomas & Pearce (1995) and Pearce & Couchman (1997). It supplements the standard P3M algorithm (Efstathiou et al. 1985) by recursively placing higher resolution meshes, ‘refinements‘, in heavily clustered regions. The Virgo simulations were done on two large Cray T3D parallel supercomputers at the Edinburgh Parallel Computing Center and at the Computing Center of the Max Plank Society in Garching. These simulations are publicly available at http://www.mpa-garching.mpg.de/Virgo/virgoproject.html.

In the simulations used here, the evolution of particles was followed in the comoving box of size $`L=239.5h^1`$ Mpc. The number of particles was $`N_p=256^3`$. Therefore, the mean particle separation $`\lambda _p=L/N_p^{1/3}=0.9355h^1`$ Mpc. The mass of the particle $`m_p=\rho _b\lambda _p^3=6.82\times 10^{10}h^1M_{}`$ and $`m_p=2.27\times 10^{11}h^1M_{}`$ in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, respectively (here $`\rho _b`$ is the mean background density). The gravitational softening length is $`r_{soft}=25h^1`$ kpc and $`r_{soft}=36h^1`$ kpc , respectively (see Jenkins et al. (1998) for a detailed description of the force calculation scheme used in the Virgo simulations). We denote the Virgo $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models as the model $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM.

We also investigated the evolution of $`256^3`$ particles on a $`256^3`$ grid using a particle-mesh (PM) code described by Gramann (1988) and Suhhonenko & Gramann (1999). The PM code achieves the force resolution close to the mean particle separation $`\lambda _p`$. The PM method is discussed in detail by Hockney & Eastwood (1981) and Efstathiou et al.(1985). The cosmological parameters for this PM simulation were chosen similar to the $`\mathrm{\Lambda }`$CDM1 model. We chose the flat $`\mathrm{\Omega }_0=0.3`$ model and used the initial power spectrum $`P(k)`$ given in eq. (1) (with $`\sigma _8=0.9`$). The comoving box size was $`L=239.5h^1`$ Mpc. Therefore, the mean particle separation, $`\lambda _p`$, and the mass of the particle, $`m_p`$, in this model are the same as used in the $`\mathrm{\Lambda }`$CDM1 model. We denote this model as the $`\mathrm{\Lambda }`$CDM2 model.

We examined the rms velocity of particles in the simulations studied. In the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, the rms velocity of particles was $`648\mathrm{km}\mathrm{s}^1`$ and $`636\mathrm{km}\mathrm{s}^1`$, respectively. In the $`\mathrm{\Lambda }`$CDM2 model, the rms velocity was $`575\mathrm{km}\mathrm{s}^1`$. Due to the small-scale smoothing inherent to the PM method, the intrinsic velocity dispersions of clusters in the $`\mathrm{\Lambda }`$CDM2 model are smaller than in the $`\mathrm{\Lambda }`$CDM1 model. The velocity field of clusters in these models is expected to be similar.

In the $`\mathrm{\Lambda }`$CDM2 model we used the traditional two-point approximation to calculate the acceleration on the grid. We also studied the evolution of particles using a PM code with the shifted-mesh scheme. The effect of the shifted-mesh scheme on the cluster velocities is discussed in Section 4.

### 2.2 Selection of clusters

We used two different algorithms to identify clusters in simulations: the standard friends-of-friends (FOF) algorithm, and the algorithm, where clusters are defined as maxima of the density field smoothed on the scale $`R1h^1`$ Mpc (DENSMAX).

The friends-of-friends group finder algorithm was applied using the program suite developed by the cosmology group in the University of Washington. These programs are available at http://www-hpcc.astro.washington.edu.

The FOF cluster finder depends on one parameter $`b`$, which defines the linking length as $`b\lambda _p`$. The conventional choice for this parameter is $`b=0.2`$ (see e.g. Götz, Huchra & Brandenberger 1998; Jenkins et al. 2001). In this paper we also define clusters by using the value $`b=0.2`$. We also study velocities of the clusters defined by the parameters $`b=0.15`$ and $`b=0.3`$. In the limit of very large numbers of particles per object, FOF approximately selects the matter enclosed by an isodensity contour at $`1/b^3`$.

To test our FOF output data, we found the mass function of clusters in the Virgo simulations. The cluster mass function in these simulations has been studied in detail by Jenkins et al. (2001). Fig. 2 shows the mass function of clusters determined for $`b=0.164`$ in the $`\mathrm{\Lambda }`$CDM model and for $`b=0.2`$ in the $`\tau `$CDM model. For comparison we show the mass functions given by the approximations obtained by Jenkins et al. (2001) (eq. (B2) and (B1) for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM model, respectively). Jenkins et al. (2001) studied the mass function at the high mass end up to the point where the predicted Poisson abundance errors reach $`10`$%. In the simulations studied here, this limit is reached when the number density of clusters is $`n=7.3\times 10^6h^3\mathrm{Mpc}^3`$. Fig. 2 shows that if $`n`$ is larger than this value, the agreement between our results and these obtained by Jenkins et al. (2001) is very good.

We studied clusters that contained at least ten particles. The three-dimensional velocity of each cluster was defined as

$$\stackrel{}{v}_{cl}=\frac{1}{N_p}\underset{i=1}{\overset{N_p}{}}\stackrel{}{v}_i,$$
$`(5)`$
where $`N_p`$ is the number of particles in the cluster and $`\stackrel{}{v}_i`$ is the velocity of the particle $`i`$ in the cluster. To characterize the size of the cluster, we use the effective radius defined as

$$R_{eff}=\frac{1}{N_p}\underset{i=1}{\overset{N_p}{}}\left[(x_i\overline{x})^2+(y_i\overline{y})^2+(z_i\overline{z})^2\right]^{1/2},$$
$`(6)`$
where $`x_i,y_i,z_i`$ are the particle coordinates in the cluster and ($`\overline{x},\overline{y},\overline{z}`$) are the coordinates of the cluster centre. If we use the FOF method, the mean size for the high-mass clusters is larger than the mean size for the low-mass clusters.

We also selected clusters using the DENSMAX method. In this case clusters were identified in the simulations as maxima of the density field that was determined on a $`256^3`$ grid using the cloud-in-cell (CIC) scheme. To determine peculiar velocities of clusters, $`v_{cl}`$, we calculated the peculiar velocity field on a $`256^3`$ grid using the CIC-scheme, and found peculiar velocities at the grid points where the clusters had been identified. If we use the DENSMAX method, the size of the selected clusters is similar for all clusters, given by the cell size, $`\lambda _g`$. In our case, the cell size $`\lambda _g=0.9355h^1`$ Mpc.

To determine the rms peculiar velocities of clusters, we used the equation

$$v_{rms}^2=v_s^2+v_L^2=\frac{1}{N_{cl}}\underset{i=1}{\overset{N_{cl}}{}}v_{cli}^2+v_L^2,$$
$`(7)`$
where the parameter $`v_s`$ describes the dispersion of cluster velocities, $`v_{cli}`$, derived from simulation and the parameter $`v_L`$ is the linear contribution from velocity fluctuations on scales greater than the size of the simulation box $`L`$ and is given by

$$v_L^2=\frac{f^2(\mathrm{\Omega }_0)H_0^2}{2\pi ^2}_0^{\frac{2\pi }{L}}P(k)𝑑k.$$
$`(8)`$
$`N_{cl}`$ is the number of clusters studied. By using eq. (4), the linear rms peculiar velocity of peaks can be written as

$$\sigma _p^2(R)=\sigma _v^2(R)H_0^2f^2(\mathrm{\Omega }_0)\frac{\sigma _0^4(R)}{\sigma _1^2(R)}.$$
$`(9)`$
The second term in this expression is not sensitive to the amplitude of large-scale fluctuations at wavenumbers $`k<2\pi /L`$. Therefore, the linear rms velocity of peaks can be expressed as

$$\sigma _p^2(R)=\sigma _p^2(R)+v_L^2,$$
$`(10)`$
where $`\sigma _p^{}(R)`$ is determined by the power spectrum at wavenumbers $`k>2\pi /L`$ and $`v_L`$ is given by eq. (8). For $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, we found that $`v_L=220\mathrm{km}\mathrm{s}^1`$ and $`v_L=245\mathrm{km}\mathrm{s}^1`$, respectively.

If the one-dimensional velocities of clusters, $`v_{xi}`$, follow a Gaussian distribution with a mean $`\overline{v}_x=0`$ and a dispersion $`\sigma ^2`$, then the sum

$$\chi ^2=\frac{1}{\sigma ^2}\underset{i=1}{\overset{N_{cl}}{}}v_{cli}^2$$
$`(11)`$
is distributed as a $`\chi ^2`$ distribution with the number of degrees of freedom $`\nu =3N_{cl}`$. In this case, the rms error for the variable $`v_{rms}^2`$ can be determined as

$$\mathrm{\Delta }v_{rms}^2=\sqrt{\frac{2}{3N_{cl}}}v_s^2.$$
$`(12)`$
As a first step, the one-dimensional distribution of cluster velocities can be approximated as a Gaussian distribution (see e.g. Bahcall, Gramann, Cen (1994) for the study of the velocity distribution of clusters in different cosmological models). We used eq. (12) to estimate the error bars for the rms velocities of clusters.

## 3 RESULTS

First, we investigated the rms velocity of clusters in different density and mass intervals. The results are presented in Fig. 3. The rms velocities are shown for the intervals, where the number of clusters is $`N_{cl}>10`$.

The clusters defined with the DENSMAX method were divided into subgroups according to their density. We studied the rms velocity of clusters in seven subgroups, where the density $`\rho /\rho _b`$ was in the range $`15`$, $`510`$, $`1050`$,…,$`10005000`$. Table 1 shows the number of clusters and the upper right panel in Fig. 3 shows the rms velocity of clusters in different density intervals. We see that the rms velocities of clusters in the $`\mathrm{\Lambda }`$CDM1 and $`\mathrm{\Lambda }`$CDM2 models are similar. In the models studied, the rms velocity of clusters is almost independent of the density. The rms velocity somewhat increases at smaller densities. However, this increase is small ($`10`$ per cent). For the range $`\rho /\rho _b=100500`$, the rms velocity is $`505`$ km s-1 and $`570`$ km s-1 in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, respectively. These values are similar to the linear theory expectations. At the radius $`R=1h^1`$ Mpc, the rms peculiar velocity of peaks $`\sigma _p=509`$ km s-1 and $`562`$ km s-1 in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, respectively. The rms velocity for the low-density clusters is somewhat smaller than predicted by the linear theory.

The clusters determined with the FOF method were divided into subgroups according to their mass. We studied the rms velocity of clusters in eight subgroups, where the mass was in the range $`(5\times 10^{11}10^{12})h^1M_{}`$, …, $`(10^{15}5\times 10^{15})h^1M_{}`$. Table 3 shows the number of clusters and Fig. 3 demonstrates the rms peculiar velocity of clusters in different mass intervals. The lower left panel in Fig. 3 shows the results for the clusters determined by $`b=0.2`$. The right panels show the rms velocities of clusters defined by $`b=0.15`$ and $`b=0.3`$.

We see that in the $`\mathrm{\Lambda }`$CDM1 and $`\mathrm{\Lambda }`$CDM2 models the rms velocities are similar. The rms velocity of FOF clusters decreases with cluster mass. The rms velocity of massive FOF clusters is smaller than the rms velocity of high-density clusters determined with the DENSMAX method. In the $`\mathrm{\Lambda }`$CDM model for $`b=0.2`$, the rms velocity of clusters is $`525\mathrm{km}\mathrm{s}^1`$ in the mass interval $`(5\times 10^{11}10^{12})h^1M_{}`$ and $`430\mathrm{km}\mathrm{s}^1`$ in the mass interval $`(10^{15}5\times 10^{15})h^1M_{}`$. For $`b=0.15`$ and $`b=0.3`$, this effect is similar. These results are in good agreement with the results obtained by Sheth & Diaferio (2001). They studied the rms velocities of clusters in different mass intervals and found that the rms cluster velocity decreases with mass.

We also studied the rms velocity of clusters with density (or mass) higher than a given threshold density (or mass). The results are presented in Fig. 4. The DENSMAX clusters were ranked according to their density, and we selected $`N_{cl}=(L/d_{cl})^3`$ highest ranked clusters to produce cluster catalogues with a mean intercluster separation $`1050h^1`$ Mpc. Similarly, the FOF clusters were ranked according to their mass. Table 2 shows the density and mass thresholds used to produce cluster catalogues for different values of the mean cluster separation. For comparison, the number density of observed APM clusters and Abell clusters is $`n_{cl}3.4\times 10^5h^3`$ Mpc-3 ($`d_{cl}31h^1`$Mpc) and $`n_{cl}2.5\times 10^5h^3`$ Mpc-3 ($`d_{cl}34h^1`$ Mpc), respectively (Dalton et al. 1994, Einasto et al. 1997).

Fig. 4 demonstrates the rms peculiar velocity of clusters with a mean separation $`d_{cl}=1050h^1`$ Mpc. In this range, the rms velocity of DENSMAX clusters is almost independent of the density of clusters. For the clusters with a mean separation $`d_{cl}=30h^1`$ Mpc, $`v_{rms}=510\mathrm{km}\mathrm{s}^1`$ and $`v_{rms}=580\mathrm{km}\mathrm{s}^1`$ in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, respectively. These values are similar to the linear rms velocity of peaks at the radius $`R1h^1`$Mpc.

The rms velocity of FOF clusters decreases with cluster richness. For rich clusters, the rms velocity of FOF clusters is smaller than the rms velocity of clusters determined with the DENSMAX method. In the $`\mathrm{\Lambda }`$CDM model for $`b=0.2`$, the rms velocities are $`500\mathrm{km}\mathrm{s}^1`$ and $`430\mathrm{km}\mathrm{s}^1`$, if $`d_{cl}=10h^1`$ Mpc and $`d_{cl}=50h^1`$ Mpc, respectively. For $`b=0.15`$ and $`b=0.3`$, this effect is similar. For the clusters with $`d_{cl}=30h^1`$ and $`b=0.2`$, $`v_{rms}=475\mathrm{km}\mathrm{s}^1`$ and $`v_{rms}=565\mathrm{km}\mathrm{s}^1`$ in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, respectively.

In Table 3 we compare the rms velocity of FOF clusters, $`v_{rms}`$, with $`\sigma _p(R)`$ for the radius $`R=1h^1`$Mpc. We analyzed the rms velocity of clusters for different values of the cluster separation. The results are given for the clusters determined by $`b=0.2`$. The rms peculiar velocities $`v_{rms}`$ for the clusters determined by $`b=0.15`$ and $`b=0.3`$ are similar. In the $`\mathrm{\Lambda }`$CDM1 model, the rms peculiar velocity of small clusters is close to the linear theory expectations, while the rms peculiar velocity of rich clusters is smaller ($``$15% per cent for clusters with a mean intercluster separation $`d_{cl}=50h^1`$ Mpc).

Our results are in good agreement with the results obtained by Colberg et al. (2000). They also studied the rms velocity of clusters in the Virgo $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, but used a slightly different method to select clusters. High-density regions were located using a FOF method with $`b=0.05`$ and their barycentres were considered as candidate cluster centers. Any candidate centre for which mass within $`1.5h^1`$ Mpc exceeded the threshold mass $`M_t`$ was identified as a candidate cluster. The final cluster list was obtained by deleting the lower mass candidate in all pairs separated by less than $`1.5h^1`$ Mpc. The peculiar velocity of each cluster was defined to be the mean peculiar velocity of all the particles within the $`1.5h^1`$ Mpc sphere. In this method, the size of the selected clusters is same for all clusters and in this sense, this method is similar to the DENSMAX method.

Colberg et al. (2000) used the value $`M_t=3.5\times 10^{14}h^1M_{}`$. For this value, the number of clusters was $`70`$ in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models ($`d_{cl}58h^1`$Mpc). They found that the rms cluster velocities derived from simulation, $`v_s`$ are $`439\mathrm{km}\mathrm{s}^1`$ and $`535\mathrm{km}\mathrm{s}^1`$ in the $`\lambda `$CDM1 and $`\tau `$CDM models, respectively (they did not include the dispersion $`v_L^2`$). If we use the DENSMAX method, we find that the rms velocity of clusters is almost independent of the number density of clusters and for $`d_{cl}=50h^1`$ Mpc, the velocities are $`v_s=450\mathrm{km}\mathrm{s}^1`$ and $`549\mathrm{km}\mathrm{s}^1`$ in the $`\mathrm{\Lambda }`$CDM1 and $`\tau `$CDM models, respectively. These values are very close to the values found by Colberg et al. (2000), only slightly larger ($`2`$ per cent). This small difference is probably caused by the fact that in the DENSMAX method we use the smoothing length $`1h^1`$ Mpc, which is smaller than $`1.5h^1`$Mpc used by Colberg et al. (2000)\]. For comparison, the rms velocities $`v_s`$ for the FOF clusters for $`b=0.2`$ and $`d_{cl}=50h^1`$ Mpc, are $`368\mathrm{km}\mathrm{s}^1`$ and $`468\mathrm{km}\mathrm{s}^1`$ in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, respectively.

Let us now consider the rms velocities of clusters for different cluster radii. We studied the rms velocity of clusters with the effective radius $`R_{eff}`$ larger than a given threshold radius. The clusters were ranked according to their effective radius and we selected $`N_{cl}=(L/d_{cl})^3`$ highest ranked clusters to produce cluster catalogues with mean separations $`1050h^1`$ Mpc. Table 4 shows the threshold radii used for different values of the mean cluster separation.

Fig. 5 illustrates the dependence of $`v_{rms}`$ on cluster radii. We see that the effect of the cluster radius on $`v_{rms}`$ is similar to the effect of the cluster mass on $`v_{rms}`$ (compare Fig. 4 and Fig. 5). In the models studied, the rms velocity of small clusters is higher than the rms velocity of large clusters. For example, in the $`\mathrm{\Lambda }`$CDM1 model for $`b=0.2`$, the rms velocities are $`520\mathrm{km}\mathrm{s}^1`$ and $`450\mathrm{km}\mathrm{s}^1`$ for clusters with $`d_{cl}=10h^1`$ Mpc and $`d_{cl}=50h^1`$ Mpc, respectively.

## 4 THE SHIFTED-MESH SCHEME

When we started the study of cluster velocities we first used a PM code with a shifted-mesh scheme. In this scheme, the acceleration on the grid was calculated by using the approximation

$$g_x(i+\frac{1}{2},j,k)=\phi (i,j,k)\phi (i+1,j,k),$$
$`(13)`$
where $`\phi (i,j,k)`$ is the gravitational potential on the grid. A similar scheme for one-dimensional systems was proposed by Melott (1986). The traditional two-point finite-difference approximation for the acceleration is given by

$$g_x(i,j,k)=\frac{\phi (i1,j,k)\phi (i+1,j,k)}{2}.$$
$`(14)`$
Gramann (1987) studied the evolution of an one-dimensional plane-wave perturbation by using a PM code with differential operators (13) and (14) and found that the approximation (13) leads to smaller deviations from the exact solution than the approximation (14).

Fig. 6 shows the rms peculiar velocities of clusters in a PM simulation, where we used the approximation (13) for a three-dimensional system. The cosmological parameters in this simulation were chosen similar to the $`\mathrm{\Lambda }`$CDM1 and $`\mathrm{\Lambda }`$CDM2 models. We investigated the evolution of $`256^3`$ particles on a $`256^3`$ grid in the comoving box of size $`L=239.5h^1`$Mpc. For comparison we show the rms peculiar velocities in the $`\mathrm{\Lambda }`$CDM2 model with a standard finite-difference approximation (14). The clusters are defined by the FOF method with $`b=0.2`$. The rms cluster velocities in the $`\mathrm{\Lambda }`$CDM2 model are similar to the $`\mathrm{\Lambda }`$CDM1 model, which achieves a force resolution smaller than the mean particle separation $`\lambda _p`$ (see Fig. 3 and Fig. 4).

We see that the shifted-mesh scheme artificially boosts cluster velocities. The effect is particularly strong for rich clusters. The rms velocities are $`560`$ km s-1 and $`735`$ kms-1, if $`d_{cl}=10h^1`$ Mpc and $`d_{cl}=50h^1`$ Mpc, respectively. This effect is probably caused by the fact that using the scheme (13), we calculate the components ($`g_x,g_y,g_z`$) at different locations on the grid and this can lead to artificial enhancement of particle acceleration in one dimension. This effect does not arise in an one-dimensional system, where $`g_y=g_z=0`$.

In the $`\mathrm{\Lambda }`$CDM2 model we used the standard two-point finite-difference approximation (14). In this model, the rms velocities of clusters are similar to the $`\mathrm{\Lambda }`$CDM1 model. Therefore, we can use the PM code to study peculiar velocities on scales that are close to the mean particle separation and force resolution. But it is important to use a correct difference operator in simulations.

## 5 SUMMARY AND DISCUSSION

In this paper we have examined the rms peculiar velocities of galaxy clusters, $`v_{rms}`$, for different cluster masses and radii. We analyzed clusters in the Virgo simulations for two cosmological models, $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM (Jenkins et al. 1998). These simulations were carried out using the AP3M code. We used the simulations where the mean particle separation $`\lambda _p1h^1`$ Mpc. We also analyzed clusters in an N-body simulation where the evolution was followed using a PM code with the same mass resolution that was used in the Virgo simulations. The cosmological parameters for this simulation were chosen similar to the Virgo $`\mathrm{\Lambda }`$CDM model. We found that the rms velocities of clusters in the PM simulation are similar to the rms velocities in the AP3M simulation. We can use the PM code to study peculiar velocities on scales that are close to the mean particle separation $`\lambda _p`$.

To identify clusters in the simulations we used two methods: the standard friends-of-friends (FOF) method and the method where the clusters are defined as maxima of the density field smoothed on the scale $`R1h^1`$ Mpc (DENSMAX). The velocity of DENSMAX clusters was determined using the same smoothing scale as for the density field. If we use the DENSMAX method, the size of the selected clusters is similar for all clusters. The velocity of FOF clusters was defined to be the mean velocity of all the particles in the cluster. If we use the FOF method, the size of the high-mass clusters is larger than the size of the low-mass clusters. We studied the velocities of FOF clusters defined by the parameters $`b=0.15`$, $`b=0.2`$ and $`b=0.3`$ (the parameter $`b`$ defines the linking length as $`b\lambda _p`$).

We found that the rms velocity of clusters defined with the DENSMAX method is almost independent of the density of clusters. The rms velocity of FOF clusters decreases with the cluster mass and radius. The effect of the cluster radius on $`v_{rms}`$ is similar to the effect of the cluster mass on $`v_{rms}`$ (see Fig. 4 and Fig. 5). The rms velocity of small clusters is higher than the rms velocity of large massive clusters. For different values of $`b`$, this effect is similar. In the $`\mathrm{\Lambda }`$CDM model, the rms peculiar velocity of massive clusters with an intercluster separation $`d_{cl}=50h^1`$ Mpc is $`15`$% smaller than the rms velocity of clusters with a separation $`d_{cl}=10h^1`$ Mpc.

The rms velocity of massive FOF clusters is smaller than the rms velocity of high-density clusters determined with the DENSMAX method. What is the reason for this difference? Do we select different objects by using different methods or do we define different velocities for the same clusters? We analyzed the fraction of DENSMAX clusters, $`F`$, which match FOF clusters (in terms of their positions). We studied the FOF clusters determined with $`b=0.2`$. In the $`\mathrm{\Lambda }`$CDM model, if we compared DENSMAX clusters for $`d_{cl}=50h^1`$Mpc with the FOF clusters for $`d_{cl}=50h^1`$Mpc, we found that the fraction of matched clusters is $`F=0.72`$. For $`d_{cl}=30h^1`$Mpc in both methods, the fraction $`F=0.74`$. In the $`\tau `$CDM model, the fraction $`F`$ was similar. We select the same objects by using different methods, but we rank them in a somewhat different way. In the DENSMAX method, the clusters are ranked according to their density and in the FOF method according to their mass. If we compared DENSMAX clusters for $`d_{cl}=50h^1`$Mpc with the FOF clusters for $`d_{cl}=40h^1`$Mpc, we found that $`F=0.98`$.

But we do assign different velocities for the same clusters by using different methods. Fig. 7 shows the peculiar velocities determined with different methods for the same clusters. We show the velocities of rich clusters that match for $`d_{cl}=50h^1`$ Mpc in both methods. We see that the velocities of clusters defined by the FOF method are systematically smaller than the velocities of clusters defined by the DENSMAX method. This is probably caused by the fact that in the FOF method, the size of rich clusters is larger than the smoothing scale $`R1h^1`$ Mpc that was used to define the cluster velocities in the DENSMAX method.

We compared the rms velocities of clusters with the linear theory predictions for the rms peculiar velocities of peaks, $`\sigma _p`$, for the smoothing radius $`R=1h^1`$ Mpc. At this radius, $`\sigma _p=509`$ kms-1 and $`\sigma _p=562`$ km s-1 in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, respectively. We analyzed the rms velocity of FOF clusters for different values of the cluster separation. In the $`\mathrm{\Lambda }`$CDM model, the rms peculiar velocity of small clusters is close to the linear theory expectations, while the rms peculiar velocity of rich clusters is smaller ($``$15 per cent for clusters with a mean separation $`d_{cl}=50h^1`$Mpc). The rms velocity of DENSMAX clusters is almost independent of the cluster density and is similar to the linear theory expectations. For the DENSMAX clusters with a mean separation $`d_{cl}=30h^1`$ Mpc, $`v_{rms}=510\mathrm{km}\mathrm{s}^1`$ and $`v_{rms}=580\mathrm{km}\mathrm{s}^1`$ in the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models, respectively. On scales probed by galaxy clusters, velocity fluctuations are in the quasi-linear regime.

## ACKNOWLEDGEMENTS

We thank J. Einasto, M. Einasto, P. Heinamäki, G. Hütsi, A. Melott and E. Saar for useful discussions. This work has been supported by the ESF grant 3601. The N-body simulations used in this paper are available at http://www.mpa-garching.mpg.de/Virgo/virgoproject.html. These simulations were carried out at the Computer Center of the Max-Planck Society in Garching and at the EPCC in Edinburgh, as part of the Virgo Consortium project. The FOF programs used in this paper are available at http://www-hpcc.astro.washington.edu. These programs were developed in the University of Washington.


text:
# Analysis of accelerations in turbulence based on generalized statistics

## Abstract

An analytical expression of probability density function (PDF) of accelerations in turbulence is derived with the help of the statistics based on generalized entropy (the Tsallis entropy or the Rényi entropy). It is revealed that the derived PDF explains the one obtained by Bodenschatz et al. in the measurement of fluid particle accelerations in fully developed turbulence at $`R_\lambda =970`$.

Key words: multifractal analysis, fully developed turbulence, PDF of fluid particle accelerations, Rényi entropy, Tsallis entropy.

PACS: 47.27.-i, 47.53.+n, 47.52.+j, 05.90.+m.

The multifractal analysis of turbulence by the statistics based on the generalized entropy of Rényi’s or of Tsallis’ has been developed by the present authors . The Rényi entropy  has the extensive character as the usual thermodynamical entropy does, whereas the Tsallis entropy  is non-extensive. The multifractal analysis belongs to the line of study based on a kind of ensemble theoretical approaches that, starting from the log-normal model , continues with the $`\beta `$-model , the p-model , the 3D binomial Cantor set model  and so on. After a rather preliminary investigation of the p-model , we developed further to derive the analytical expression for the scaling exponents of velocity structure function , and to determine the probability density function (PDF) of velocity fluctuations  and of velocity derivative  by a self-consistent statistical mechanical approach.

In this paper, we will derive the formula for the PDF of the accelerations of a fluid particle in fully developed turbulence by means of the multifractal analysis. With the theoretical PDF, we will analyze the PDF of accelerations at $`R_\lambda =970`$ (the Taylor microscale Reynolds number) obtained in the Lagrangian measurement of particle accelerations that was realized by Bodenschatz and co-workers  by raising dramatically the spatial and temporal measurement resolutions with the help of the silicon strip detectors.

We assume that the turbulent flow, satisfying the Navier-Stokes equation

$$\stackrel{}{u}/t+(\stackrel{}{u}\stackrel{}{})\stackrel{}{u}=\stackrel{}{}\left(p/\rho \right)+\nu ^2\stackrel{}{u}$$
(1)
of an incompressible fluid, consists of a cascade of eddies with different sizes $`\mathrm{}_n=\delta _n\mathrm{}_0`$ where $`\delta _n=2^n`$ $`(n=0,1,2,\mathrm{})`$. The quantities $`\rho `$, $`p`$ and $`\nu `$ represent, respectively, the mass density, the pressure and the kinematic viscosity. The acceleration $`\stackrel{}{\mathrm{a}}`$ of a fluid particle is given by the substantive time derivative of the velocity: $`\stackrel{}{\mathrm{a}}=\stackrel{}{u}/t+(\stackrel{}{u}\stackrel{}{})\stackrel{}{u}`$. At each step of the cascade, say at the $`n`$th step, eddies break up into two pieces producing an energy cascade with the energy-transfer rate $`ϵ_n`$ that represents the rate of transfer of energy per unit mass from eddies of size $`\mathrm{}_n`$ to those of size $`\mathrm{}_{n+1}`$ (the energy cascade model). The Reynolds number $`\mathrm{Re}`$ of the system is given by $`\mathrm{Re}=\delta u_0\mathrm{}_0/\nu =(\mathrm{}_0/\eta )^{4/3}`$ with the Kolmogorov scale  $`\eta =(\nu ^3/ϵ)^{1/4}`$ where $`ϵ`$ ($`=ϵ_0`$) is the energy input rate to the largest eddies with size $`\mathrm{}_0`$.111 The velocity fluctuation $`\delta u_n`$ is defined by $`\delta u_n=|u(+\mathrm{}_n)u()|`$ where $`u`$ represents a component of the velocity field $`\stackrel{}{u}`$.  Introducing the pressure (divided by the mass density) difference $`\delta p_n`$ at two points separated by the distance $`\mathrm{}_n`$, i.e., $`\delta p_n=|p/\rho (+\mathrm{}_n)p/\rho ()|`$, and the acceleration $`\mathrm{a}_n=\delta p_n/\mathrm{}_n`$ belonging to the $`n`$th step in the energy cascade, one can estimate accelerations by $`|\stackrel{}{\mathrm{a}}|=lim_n\mathrm{}\mathrm{a}_n`$. For high Reynolds number $`\mathrm{Re}1`$, or for the situation where effects of the kinematic viscosity $`\nu `$ can be neglected compared with those of the turbulent viscosity, the Navier-Stokes equation (1) is invariant under the scale transformation : $`\stackrel{}{r}\lambda \stackrel{}{r}`$, $`\stackrel{}{u}\lambda ^{\alpha /3}\stackrel{}{u}`$, $`t\lambda ^{1\alpha /3}t`$ and $`\left(p/\rho \right)\lambda ^{2\alpha /3}\left(p/\rho \right)`$. The exponent $`\alpha `$ is an arbitrary real quantity which specifies the degree of singularity in the acceleration for $`\alpha <1.5`$, i.e., $`lim_n\mathrm{}\mathrm{a}_n=lim_{\mathrm{}_n0}\delta p_n/\mathrm{}_nlim_{\mathrm{}_n0}\mathrm{}_n^{(2\alpha /3)1}\mathrm{}`$ which can be seen with the relation $`\delta p_n/\delta p_0=(\mathrm{}_n/\mathrm{}_0)^{2\alpha /3}`$.

The multifractal analysis rests on the assumption that the distribution of the exponent $`\alpha `$ is multifractal, and that the probability $`P^{(n)}(\alpha )d\alpha `$ to find, at a point in physical space, an eddy of size $`\mathrm{}_n`$ having a value of the degree of singularity in the range $`\alpha \alpha +d\alpha `$ is given by

$`P^{(n)}(\alpha )\left[1(\alpha \alpha _0)^2/(\Delta \alpha )^2\right]^{n/(1q)}`$ (2)
with $`(\Delta \alpha )^2=2X/[(1q)\mathrm{ln}2]`$. The range of $`\alpha `$ is $`\alpha _{\mathrm{min}}\alpha \alpha _{\mathrm{max}}`$ with $`\alpha _{\mathrm{min}}=\alpha _0\Delta \alpha `$ and $`\alpha _{\mathrm{max}}=\alpha _0+\Delta \alpha `$. Here, we assume that the distribution function at the $`n`$th multifractal depth has the structure $`P^{(n)}(\alpha )[P^{(1)}(\alpha )]^n`$. This is consistent with the relation  $`P^{(n)}(\alpha )\delta _n^{1f(\alpha )}`$ that is a manifestation of scale inveriance and reveals how densely each singularity, labeled by $`\alpha `$, fills physical space. Within the present model, the multifractal spectrum $`f(\alpha )`$ is given by

$$f(\alpha )=1+(1q)^1\mathrm{log}_2\left[1\left(\alpha \alpha _0\right)^2/\left(\mathrm{\Delta }\alpha \right)^2\right].$$
(3)
To make the paper self-contained, we put here its brief derivation. The distribution function (2) is derived by taking an extremum of the generalized entropy, the Rényi entropy  $`S_q^\mathrm{R}[P^{(1)}(\alpha )]=\left(1q\right)^1\mathrm{ln}𝑑\alpha P^{(1)}(\alpha )^q`$ or the Tsallis entropy  $`S_q^\mathrm{T}[P^{(1)}(\alpha )]=\left(1q\right)^1\left(𝑑\alpha P^{(1)}(\alpha )^q1\right)`$, under the two constraints, i.e., the normalization of distribution function: $`𝑑\alpha P^{(1)}(\alpha )=\text{const.}`$ and the $`q`$-variance being kept constant as a known quantity: $`\sigma _q^2=(𝑑\alpha P^{(1)}(\alpha )^q(\alpha \alpha _0)^2)/𝑑\alpha P^{(1)}(\alpha )^q`$. In spite of the different characteristics of these entropies, the distribution function giving their extremum has the common structure (2). The dependence of the parameters $`\alpha _0`$, $`X`$ and $`q`$ on the intermittency exponent $`\mu `$ is determined, self-consistently, with the help of the three independent equations, i.e., the energy conservation: $`ϵ_n=ϵ`$, the definition of the intermittency exponent $`\mu `$: $`ϵ_n^2=ϵ^2\delta _n^\mu `$, and the scaling relation:222 The scaling relation is a generalization of the one derived first in  to the case where the multifractal spectrum has negative values .  $`1/(1q)=1/\alpha _{}1/\alpha _+`$ with $`\alpha _\pm `$ satisfying $`f(\alpha _\pm )=0`$. The average $`\mathrm{}`$ is taken with $`P^{(n)}(\alpha )`$. For the region where the value of $`\mu `$ is usually observed, i.e., $`0.13\mu 0.40`$, the three self-consistent equations are solved to give the approximate equations : $`\alpha _0=0.9989+0.5814\mu `$, $`X=2.848\times 10^3+1.198\mu `$ and $`q=1.507+20.58\mu 97.11\mu ^2+260.4\mu ^3365.4\mu ^4+208.3\mu ^5`$.

Since there are two mechanisms in turbulent flow to rule its evolution, i.e., the one controlled by the kinematic viscosity that takes care thermal fluctuations, and the other by the turbulent viscosity that is responsible for intermittent fluctuations related to the singularities in acceleration, it may be reasonable to assume that the probability $`\mathrm{\Lambda }^{(n)}(y_n)dy_n`$ to find the scaled pressure fluctuations $`|y_n|=\delta p_n/\delta p_0`$ in the range $`y_ny_n+dy_n`$ has two independent origins:

$$\mathrm{\Lambda }^{(n)}(y_n)dy_n=\mathrm{\Lambda }_\mathrm{S}^{(n)}(|y_n|)dy_n+\mathrm{\Delta }\mathrm{\Lambda }^{(n)}(y_n)dy_n.$$
(4)
The singular part $`\mathrm{\Lambda }_\mathrm{S}^{(n)}(|y_n|)`$ of the PDF stemmed from multifractal distribution of the singularities, and the correction part $`\mathrm{\Delta }\mathrm{\Lambda }^{(n)}(y_n)`$ from thermal dissipation and/or measurement error.333 Needless to say that each term in (4) is a multiple of two PDF’s, i.e., the PDF for one of the two independent origins to realize and the conditional PDF for a value $`\omega _n`$ in the range $`\omega _n\omega _n+d\omega _n`$ to come out. This is of course in a generalized sense in which the second correction term may weaken the first singular contribution.  The former is derived through $`\mathrm{\Lambda }_\mathrm{S}^{(n)}(|y_n|)dy_nP^{(n)}(\alpha )d\alpha `$ with the transformation of the variables: $`|y_n|=\delta _n^{2\alpha /3}`$. The $`m`$th moments of the pressure fluctuations are given by

$$|y_n|^m_{\mathrm{}}^{\mathrm{}}𝑑y_n|y_n|^m\mathrm{\Lambda }^{(n)}(y_n)=2\stackrel{~}{\gamma }_m^{(n)}+(12\stackrel{~}{\gamma }_0^{(n)})a_{2m}\delta _n^{\zeta _{2m}}$$
(5)
with $`a_{3\overline{q}}=\{2/[\sqrt{C_{\overline{q}}}(1+\sqrt{C_{\overline{q}}})]\}^{1/2}`$, $`C_{\overline{q}}=1+2\overline{q}^2(1q)X\mathrm{ln}2`$ and

$$2\stackrel{~}{\gamma }_m^{(n)}=_{\mathrm{}}^{\mathrm{}}𝑑y_n|y_n|^m\mathrm{\Delta }\mathrm{\Lambda }^{(n)}(y_n).$$
(6)
We used the normalization: $`1=1`$. The quantity

$$\zeta _m=\frac{\alpha _0m}{3}\frac{2Xm^2}{9\left(1+\sqrt{C_{m/3}}\right)}\frac{1}{1q}\left[1\mathrm{log}_2\left(1+\sqrt{C_{m/3}}\right)\right]$$
(7)
is the so-called scaling exponent of velocity structure function, whose expression was derived first by the present authors . The formula explains quite well experimental data . Note that the formula is independent of the length $`\mathrm{}_n`$, and, therefore, independent of $`n`$.

Let us derive the PDF $`\widehat{\mathrm{\Lambda }}^{(n)}(\omega _n)`$ of the accelerations by dividing it into two parts with respect to $`\omega _n`$, i.e.,

$$\widehat{\mathrm{\Lambda }}^{(n)}(\omega _n)=\{\begin{array}{cc}\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)\hfill & \text{for }|\omega _n|\omega _n^{}\hfill \\ \widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)\hfill & \text{for }\omega _n^{}|\omega _n|\omega _n^{\mathrm{max}}\hfill \end{array}$$
(8)
where the scaled variable $`\omega _n`$ is defined by $`|\omega _n|=\mathrm{a}_n/\mathrm{a}_n^2^{1/2}=\overline{\omega }_n\delta _n^{2\alpha /3\zeta _4/2}`$ with $`\overline{\omega }_n=[2\stackrel{~}{\gamma }_2^{(n)}\delta _n^{\zeta _4}+(12\stackrel{~}{\gamma }_0^{(n)})a_4]^{1/2}`$, and $`\omega _n^{\mathrm{max}}=\overline{\omega }_n\delta _n^{2\alpha _{\mathrm{min}}/3\zeta _4/2}`$. Assuming that, for smaller accelerations $`|\omega _n|\omega _n^{}`$, the contribution to the PDF comes, mainly, from thermal fluctuations related to the kinematic viscosity or from measurement error, we take for the PDF $`\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)\omega _n=\mathrm{\Lambda }_\mathrm{S}^{(n)}(y_n)dy_n+\mathrm{\Delta }\mathrm{\Lambda }^{(n)}(y_n)dy_n`$ a Gaussian function, i.e.,

$$\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)=\stackrel{~}{\mathrm{\Lambda }}_\mathrm{S}^{(n)}\mathrm{exp}\left\{\frac{1}{2}\left[1+\frac{3}{2}f^{}(\alpha ^{})\right]\left[\left(\omega _n/\omega _n^{}\right)^21\right]\right\}$$
(9)
with $`\stackrel{~}{\mathrm{\Lambda }}_\mathrm{S}^{(n)}=3(12\stackrel{~}{\gamma }_0^{(n)})/(4\overline{\omega }_n\sqrt{2\pi X|\mathrm{ln}\delta _n|})`$. On the other hand, we assume that the main contribution to $`\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)`$ may come from the multifractal distribution of singularities related to the turbulent viscosity, i.e., $`\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)d\omega _n=\mathrm{\Lambda }_\mathrm{S}^{(n)}(|y_n|)dy_n`$:

$$\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)=\stackrel{~}{\mathrm{\Lambda }}_\mathrm{S}^{(n)}\frac{\overline{\omega }_n}{|\omega _n|}\left[1\frac{1q}{n}\frac{\left(3\mathrm{ln}|\omega _n/\omega _{n,0}|\right)^2}{8X|\mathrm{ln}\delta _n|}\right]^{n/(1q)}$$
(10)
with $`|\omega _{n,0}|=\overline{\omega }_n\delta _n^{2\alpha _0/3\zeta _4/2}`$. The point $`\omega _n^{}`$ was defined by $`\omega _n^{}=\overline{\omega }_n\delta _n^{2\alpha ^{}/3\zeta _4/2}`$ with $`\alpha ^{}`$ being the smaller solution of $`\zeta _4/22\alpha /3+1f(\alpha )=0`$, at which $`\widehat{\mathrm{\Lambda }}^{(n)}(\omega _n^{})`$ has the least $`n`$-dependence for large $`n`$. Here, $`\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)`$ and $`\widehat{\mathrm{\Lambda }}_<^{(n)}(\omega _n)`$ were connected at $`\omega _n^{}`$ under the condition that they should have the same value and the same slope there. The specific form of Gaussian function (9) comes out through the connection.

With the expression (9), we can obtain $`\mathrm{\Delta }\mathrm{\Lambda }^{(n)}(y_n)`$, and have the formula to evaluate $`\stackrel{~}{\gamma }_m^{(n)}`$ in the form

$$2\stackrel{~}{\gamma }_m^{(n)}=\left(K_m^{(n)}L_m^{(n)}\right)/\left(1+K_0^{(n)}L_0^{(n)}\right)$$
(11)
where

$`K_m^{(n)}={\displaystyle \frac{3\delta _n^{2(m+1)\alpha ^{}/3\zeta _4/2}}{2\sqrt{2\pi X|\mathrm{ln}\delta _n|}}}{\displaystyle _0^1}𝑑zz^m\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}\left[1+{\displaystyle \frac{3}{2}}f^{}(\alpha ^{})\right]\left(z^21\right)\right\},`$ (12)
$`L_m^{(n)}={\displaystyle \frac{3\delta _n^{2m\alpha ^{}/3}}{2\sqrt{2\pi X|\mathrm{ln}\delta _n|}}}{\displaystyle _{z_{\mathrm{min}}}^1}𝑑zz^{m1}\left[1{\displaystyle \frac{1q}{n}}{\displaystyle \frac{\left(3\mathrm{ln}|z/z_0^{}|\right)^2}{8X|\mathrm{ln}\delta _n|}}\right]^{n/(1q)}`$ (13)
with $`z_{\mathrm{min}}=\omega _{\mathrm{min}}/\omega _n^{}=\delta _n^{2(\alpha _{\mathrm{max}}\alpha ^{})/3}`$ and $`z_0^{}=\omega _{n,0}/\omega _n^{}=\delta _n^{2(\alpha _0\alpha ^{})/3}`$. Now, the PDF of fluid particle accelerations (8) is determined by the intermittency exponent $`\mu `$ and the number $`n`$ of steps in the energy cascade which gives the eddy size $`\mathrm{}_n`$.

The comparison between the present PDF of accelerations and that measured in the experiment  at $`R_\lambda =970`$ is plotted in Fig. 1 on log and linear scale. The intermittency exponent $`\mu =0.259`$ and the number $`n=16.4`$ of steps in the cascade are extracted by the method of least squares with respect to the logarithm of PDF’s as the best fit of our theoretical formulae (8) to the observed values of the PDF  by discarding those points whose values are less than $`2\times 10^6`$ since they scatter largely in the log scale. Substituting the extracted value of $`\mu `$ into the self-consistent equations, we have the values of parameters: $`q=0.431`$, $`\alpha _0=1.149`$ and $`X=0.307`$. With these values, other quantities are determined, e.g., $`\Delta \alpha =1.248`$, $`\alpha _+\alpha _0=\alpha _0\alpha _{}=0.7211`$ and $`\omega _n^{}=0.557`$ ($`\alpha ^{}=1.005`$). We see an excellent agreement between the measured PDF of accelerations and the analytical formula of PDF derived by the present self-consistent multifractal analysis. The value of $`\omega _n^{}`$ tells us that the contribution from thermal fluctuation and/or measurement error is restricted to smaller values of $`\omega _n`$, i.e., less than about one half of its standard deviation. Then, the values of $`\alpha `$ responsible for the intermittency due to the scale invariance turn out to be smaller than $`\alpha ^{}1`$. This is the range within the condition $`\alpha <1.5`$ in which the singularity appears in fluid particle accelerations.

The values of the scaling exponents $`\zeta _m`$, given by (7), for $`m=1\mathrm{}20`$ are listed in Table 1 for future convenience in comparison with experiments or other theories. Note that $`\mu `$ is related to $`\zeta _6`$ by the relation $`\mu =2\zeta _6`$ within the present analysis. The flatness $`F_\mathrm{a}^{(n)}\mathrm{a}_n^4/\mathrm{a}_n^2^2=\omega _n^4`$ of the PDF of accelerations has the value $`F_\mathrm{a}^{(n)}=57.8`$ which is compatible with the value of the flatness $`60`$ reported in  as it should be. It is quite attractive to see that the distance $`r=\mathrm{}_n`$ corresponding to the extracted value $`n=16.4`$ reduces to $`r=0.821\mu `$m ($`r/\eta =0.0456`$), which is close to the value of the spatial resolution $`0.5\mu `$m ($`1/40`$ of the Kolmogorov distance) of the measurement in .

We derived, in this letter, the formula of the PDF of accelerations (8) within the approach of the multifractal analysis constructed by the present authors, and analyzed successfully the beautiful experiment conducted by Bodenschatz and co-workers  in the Lagrangian frame. We expect that other data for the same system, in addition to the PDF of accelerations, will be provided in the near future, such as the intermittency exponent $`\mu `$, the scaling exponents $`\zeta _m`$, the PDF’s of the velocity fluctuations or of the velocity derivatives. Then, we can cross-check the validity of the present analysis based on the formula (8).444 We have checked, with the help of the DNS data reported by Gotoh et al. in , the accuracy of the value of the intermittency exponent extracted out of the measured PDF’s. With the formula (7), we determined in  the value $`\mu =0.240`$ for the longitudinal velocity fluctuations by fitting, with the method of least squares, the ten data of the scaling exponents $`\zeta _m`$ ($`m=1,2,\mathrm{},10`$) at $`\mathrm{R}_\lambda =381`$. On the other hand, we have extracted the value $`\mu =0.237`$ by the method of least squares as the best fit of the theoretical formula for the PDF of velocity derivatives, derived in  by the multifractal analysis, to the corresponding PDF data.  Note that the formula for the PDF of accelerations is different from the one for the PDF’s of velocity fluctuations or of velocity derivatives. The empirical PDF, $`\widehat{\mathrm{\Lambda }}_{\mathrm{emp}}(\omega )=C\mathrm{exp}\left\{\omega ^2/\left[\left(1+|\omega \beta /\sigma |^\gamma \right)\sigma ^2\right]\right\}`$ with $`\beta =0.539`$, $`\gamma =1.588`$, $`\sigma =0.508`$ and $`C=0.786`$ proposed in  for the data at $`R_\lambda =970`$, gives a line very close to the one provided by the present PDF $`\widehat{\mathrm{\Lambda }}^{(n)}(\omega )`$ of (8) for the region $`|\omega |<30`$ where the data of PDF exist. They deviates, however, for $`|\omega |>30`$, i.e., $`\widehat{\mathrm{\Lambda }}_{\mathrm{emp}}(\omega )<\widehat{\mathrm{\Lambda }}^{(n)}(\omega )`$. The extraction of the PDF of accelerations out of the DNS data obtained by Gotoh et al.  is one of the attractive investigations in order to check the validity of the present formula (8) derived by the unified approach providing various PDF’s , since the accuracy of their DNS data for PDF’s is very high up to the order of 10-10 .

The authors would like to thank Prof. C. Tsallis for his fruitful comments with encouragement. The authors are grateful to Prof. E. Bodenschatz for his kindness to show his data prior to publication.