text:
# Single electron measurements with a micromechanical resonator

## I Introduction

As devices for processing and storing information become smaller the demands on the readout technology become ever greater. This is especially true for proposed solid state quantum computers which store information in various quantum degrees of freedom (qubits): in quantum dotsld98 , nuclear spink98 ; pvk98 , superconducting islandsnpt99  and persistent currentsfpctl2000 , to cite just a small sample.

Kane has proposed storing a qubit in the spin of a single phosphorous nucleus implanted in silicon. In his original readout scheme, this was coupled by the hyperfine interaction to the spin of the donor electron bound weakly to the nucleus. A surface gate would then draw the electron towards an adjacent ancilla donor, to which it might tunnel, producing a doubly charged $`D^{}`$ state. Under appropriate bias conditions, this transfer can only occur if the nuclear spin of the qubit is oriented opposite the ancilla.

A spin measurement is thus reduced to detecting the transfer of a single electron charge to the ancilla. This can be done by a sensitive electroscope such as a single electron transistorkmdcmsw2000 . However, the techniques used for fabricating microelectronics have recently been adapted to build mechanical structures at micron and even nanometre scalesr2000 , and mechanical electroscopes sensitive to small numbers of electrons have been constructedcr98 . We will consider how effectively such devices might perform the measurements required for quantum information processing.

Classical treatments of measurement sensitivity assume that the observable being measured has a definite value, which influences the measuring instrument in a definite way. The only question is how much data we must gather to reliably distinguish this effect from other influences on the apparatus, which produce noise. Once we know the size of the effect we wish to distinguish, and the level of noise in the system, some elementary statistics tell us the integration time required for a reliable measurement.

This assumption does not hold when we measure an observable of a quantum system. If the system is in a superposition state, the observable will not have a definite value until some sort of measurement is carried out. Any interesting quantum information device will produce such superpositions. The process by which the superposition is reduced so that the observable has a certain value imposes a minimum level of noise in the measurement, which might be increased by the same sources of technical noise that affect measurements of classical systems.

In the proposed readout scheme for the Kane computer, a donor electron is induced to tunnel between two phosphorous nuclei, depending on the state of the nuclear spins. In general, the nuclei are in a superposition of a state which would permit tunneling, and one which would prevent it. After this tunneling has occured, the electon is left in a superpostion of two position states, each localised on one nucleus. It then interacts with the electroscope, and in general with other degrees of freedom in the crystal lattice, with the result that we see it become localised on one nucleus or the other, so that the electroscope gives a definite signal that the charge is present or absent.

Note that we are not discussing an ensemble of quantum systems subject to a single measurement, but rather a single quantum system subject to a dynamical measurement process. In such a situation we need to be able to describe the instantaneous conditional state of the measured system as the measurement results accumulate. This is quite different from the usual situation that prevails in condensed matter systems, where typically a measurement is made on a large number of (almost) identical constituents undergoing quantum dynamics, and the measurement results are already an average over an ensemble. Fortunately mathematical techniques (known as quantum trajectory methods) are available to describe the conditional dynamics of a single quantum system subject to measurement with added noise, and these methods have been applied with considerable success to experiments in quantum optics and ion trapspk98 . Recently such methods have been applied to mesoscopic electronic systemswusmkdc2001 ; gmws2001 ; k99a .

## II The Mechanical Electroscope

The operation of a micro-mechanical electroscope is shown schematically in Figure 1. The active part is an electrode, mounted on a cantilever no longer than $`1\mu \text{m}`$, which is set in motion near the charge to be measured. The electrode is held at constant potential, so that its motion with respect to the unknown charge induces a flow of charge between it and its voltage source. The induced charge gives the electroscope electric potential energy as well as elastic, and changes its resonant frequency. If we envision the electroscope being used to readout a qubit in a quantum computer, there will be two charge states we wish to distinguish. We will denote the difference between the resonant frequencies of the cantilever in these two states by $`\delta \omega `$; it is determined by geometry and the mutual capacitance between the electrode and the measured charge distribution.

We will assume the mechanical motion of the cantilever is elastic and treat it as a simple harmonic oscillator. Then its motion, including the capacitive coupling to the target charge, is described by a harmonic oscillator Hamiltonian

$$H=\mathrm{}(\omega _0+\delta \omega n_1)c^{}c$$
(1)
where $`\omega _0`$ is the resonant frequency of the cantilever in the absence of surface charge, and $`c`$ the annihilation operator for its oscillation. The observable $`n_1`$ will be defined shortly.

During readout of a Kane computer, a single donor electron may occupy a bound state around either of two adjacent nuclei. We will denote these distinct spatial states by $`|\psi `$ and $`|\varphi `$. Only one state (suppose $`|\psi `$) couples to the electroscope - this is how we can distinguish them.

During readout, the surface gates will be configured to produce tunneling between the two nuclei, depending on the state of the nuclear spin qubits. This entangles the charge states with the qubit states $`|`$ and $`|`$. We will denote the combined states by $`|\text{0}=||\varphi `$, and $`|\text{1}=||\psi `$, according to the number of electrons interacting with the electroscope, which we will represent by the operator $`n_1=|\text{1}\text{1}|`$. In general the measured qubit will be in a superposition state, so the total state will take the form

$$|\mathrm{\Psi }=a|\text{0}+b|\text{1}$$
(2)
Table 1 gives numerical parameters for a cantilever electroscope fabricated in 1998. The frequency and operating temperature of this electroscope meant that themal noise completely dominated any quantum effects. Besides lowering the temperature, this could be changed by using a cantilever with a higher resonant frequency, and such devices have been fabricated. However, the sensitivity of the electroscope depends on the frequency changing significantly when change is present, and this might not be the case in higher frequency cantilevers.

We note that the interaction Hamiltonian commutes with the number operator $`\widehat{n}_1`$. Furthermore, in the absence of tunnelling, the free Hamiltonian for the charge state itself is proportional to the square of the charge (capacitive electrostatic energy) and itself commutes with the charge number operator. In the presence of the measurement the number operator is thus a constant of motion. Such a measurement is known as a quantum nondemolition measurementwam . Number eigenstates are not changed by the coupling to the apparatus, and moments of the number operator are constant in time. On the other hand any state that is initially a coherent superposition in the number basis will be reduced to a mixture diagonal in this basis, a process known as decoherence. In an ideal quantum nondemolition measurement, the probability distribution for observed results at the conclusion of the measurement should accurately reflect the intrinsic probability distributions of the quantum nondemolition variable in the quantum state at the start of the measurement.

This model, where the electroscope performs a QND measurement of the coupled charge, is idealised. If such an electroscope was used to measure any interesting device, the motion of the cantilever would disrupt the distribution of charge being measured. The nature and extent of this disruption would depend on the electrical properties of the system being measured; for the Kane computer, determining these is an unsolved problem in atomic physics. In general, back action (and interference from sources unrelated to the measurement) imposes a time limit on the measurement, after which the charge state will have been disrupted and the results will be meaningless. The results of this paper determine whether the electroscope can measure the charge with the necessary precision within that time.

To detect the change in resonant frequency, we must set the cantilever in motion with some driving mechanism. In the device described in Table 1, this was supplied by driving an alternating current through a wire on the cantilever in the presence of a magnetic field. The current induced by the field in another wire was used to monitor the response of the cantilever to the driving.

However, the details of the driving are not important. As long as the cantilever is coupled weakly to the driving system and is not damped so strongly that its state changes significantly over the period of its vibration (in other words, it has high finesse), the effect can be described by a Hamiltonian. In the interaction picture this takes the form $`\mathrm{}(\widehat{c}+\widehat{c}^{})`$, where $``$ is the strength of the driving in units of frequency. If the finess of the cantilever is low, noise from the driving system affects its motion significantly, and the dynamics due to the driving can not be approximated by a Hamiltonian.

The frequency shift could be detected in a number of ways. We could sweep the driving frequency and monitor the amplitude of the oscillations. Or else we could drive the oscillator at a constant frequency $`\omega `$, and then detect the change in phase of the oscillation due to the shift in resonance frequency when a small charge is coupled; this is the method analysed in this paper. We will assume that if the charge state is $`|\text{0}`$, the cantilever will be driven on resonance; if it is $`|\text{1}`$, the change $`\delta \omega `$ in its resonant frequency will cause its phase to differ from that of the driving force. The rate of change of the phase of the output current with frequency of the driving is greatest when the cantilever is driven near its resonant frequency.

We will measure time by the inverse damping rate $`\gamma ^1`$. Then, defining a dimensionless driving strength $`E=/\gamma `$ and a detuning $`\mathrm{\Delta }=\delta \omega /\gamma `$, the Hamiltonian for the coherent driving, in the interaction picture, is

$$\widehat{H}_D=\mathrm{}E\left(\widehat{c}+\widehat{c}^{}\right)+\mathrm{}\mathrm{\Delta }n_1\widehat{c}^{}\widehat{c}$$
(3)
In reality of course the mechanical oscillations of the cantilever will be subject to frictional damping, and accompanying mechanical noise. The rate of energy dissipation is specified by the quality factor, $`Q`$, which is the ratio of the resonance frequency to the width of the resonance. For linear response, this gives $`Q=\omega _0/\gamma _M`$, where $`\gamma _M`$ is the decay rate of energy due to mechanical dissipation. Roukes et al.cr99  have measured quality factors up to $`2\times 10^4`$. With such quality factors and resonance frequencies approaching GHz, these devices are approaching low quality optical resonators. So we will treat the effect of mechanical damping with the master equation methods of quantum optics. These methods assume that the coupling of the resonator to the dissipative degrees of freedom is sufficiently weakhr85 ; hrsw86 . Specifically we assume that $`\gamma _M\omega _0,kT/\mathrm{}`$.

Under these assumptions the coupling between the oscillator and the thermal mechanical reservoir isgardiner

$$H_M=\sqrt{\gamma _M}(ca^{}(t)+c^{}a(t))$$
(4)
where $`a(t),a^{}(t)`$ are bosonic reservoir operators. The state of the reservoir will be taken to be that of a Planck thermal equilibrium density operator with temperature $`T_M`$.

We now consider in more detail the mechanism by which the small changes in resonance frequency induced by the proximity of a target charge are transduced. This may be donecr98  by fabricating a wire loop on the mechanical oscillator and placing the whole apparatus in a strong magnetic field. As the mechanical oscillator moves, an induced EMF is set up in the loop and we may measure the induced current. When the current for the driving circuit is such as to drive the mechanical oscillator at its resonance frequency, the induction current is out of phase with the driving current. However when a small target charge shifts the resonance frequency of the oscillator, the induced current shifts in phase with respect to the driving current. We can detect this phase shift by an electrical comparison of the driving current and induction current. This is essentially homodyne detection in which the driving current plays the role of a local oscillator. Unfortunately this electrical transduction of the mechanical motion introduces another source of noise for the measurement.

The induction current is coupled into an external amplifier circuit which can be treated as a bosonic reservoir, with some non zero noise temperaturelouisell , $`T_E`$. The readout circuit variable coupled to the cantilever is the current operator $`i(t)`$ in the readout circuit. We will assume that the coupling is linear in the current and coordinate degree of freedom of the cantilever. Under standard assumptions the interaction between the mechanical oscillator and the readout circuit is described by the interaction picture Hamiltonian,

$$H_R=i\sqrt{\gamma _E}\left(c^{}\mathrm{\Gamma }(t)c\mathrm{\Gamma }^{}(t)\right)$$
(5)
where $`\mathrm{\Gamma }(t)=b(t)e^{i\omega _0t}`$ with the actual current in the circuit given by $`i(t)=\sqrt{\mathrm{}\omega _0/2Lz_0}(b(t)+b^{}(t))`$ , $`L`$ being the inductance per unit length of the transmission line, and $`z_0`$ the quantisation length. We will assume that the readout circuit reservoir is bosonic and also in thermal equilibrium at some temperature $`T_E`$.

Using the interaction Hamiltonians for the reservoir coupling (Equations 4 and 5), we may obtain the Heisenberg equations of motion for the oscillator and reservoir variables. Using standard techniquesgardiner , the reservoir variables may be eliminated to give a quantum Langevin stochastic differential equation describing the dynamics of the oscillator amplitude

$$\frac{da}{dt}=i\delta \omega ai\frac{\gamma _M}{2}a\frac{\gamma _E}{2}a+\sqrt{\gamma _M}a_{\text{in}}(t)+\sqrt{\gamma _E}b_{\text{in}}(t)$$
(6)
where $`a_{\text{in}}(t)`$, $`b_{\text{in}}(t)`$ are the quantum noise sources for the mechanical and electrical reservoirs respectively. These noise terms are defined by correlation functions, which are Fourier transforms of

$`a_{\text{in}}(t)`$ $`=`$ $`b_{\text{in}}(t)=0`$ (7)
$`a_{\text{in}}^{}(\omega )a_{\text{in}}(\omega )^{}`$ $`=`$ $`\overline{n}(\omega ,T_M)\delta (\omega \omega ^{})`$ (8)
$`a_{\text{in}}(\omega )a_{\text{in}}^{}(\omega )^{}`$ $`=`$ $`(\overline{n}(\omega ,T_M)+1)\delta (\omega \omega ^{})`$ (9)
$`b_{\text{in}}^{}(\omega )b_{\text{in}}(\omega )^{}`$ $`=`$ $`\overline{n}(\omega ,T_E)\delta (\omega \omega ^{})`$ (10)
$`b_{\text{in}}(\omega )b_{\text{in}}^{}(\omega )^{}`$ $`=`$ $`(\overline{n}(\omega ,T_E)+1)\delta (\omega \omega ^{})`$ (11)
where

$$\overline{n}(\omega ,T)=\frac{1}{2}\left(\mathrm{coth}\left(\mathrm{}\omega /2k_BT\right)1\right)$$
(12)
Note the equation explicitly includes a friction term (proportional to $`\gamma _E`$) that arises form the electrical coupling to the readout circuit. The steady state average amplitude $`\alpha _n=a(t)_t\mathrm{}`$, is given by

$$\alpha _n=\frac{2i}{(\gamma _M+\gamma _E)+2i\delta \omega n_1}$$
(13)
The actual measured quantity is the current in the readout circuit, that is to say the readout variable is an electrical bath variable, $`b_{\text{out}}`$ at the output from the system interaction. The output amplitudes for both the mechanical and electrical baths are related to the input variables for these two baths and the amplitude of the mechanical oscillator bywam

$`a_{\text{out}}(t)`$ $`=`$ $`\sqrt{\gamma _M}a(t)a_{\text{in}}(t)`$ (14)
$`b_{\text{out}}(t)`$ $`=`$ $`i\sqrt{\gamma _E}a(t)b_{\text{in}}(t)`$ (15)
The average value of the electrical readout amplitude in the steady state is then found using equations Equation 6 and Equation 15.

$$b_{\text{out}}=\sqrt{\gamma _E}\alpha _n$$
(16)
where $`\alpha _n`$ is given in Equation 13. We see that the steady state amplitude of the cantilever, and hence the output electrical signal undergoes a change in phase and amplitude, see Figure 2. If we monitor the component in the imaginary direction (that is, in quadrature with the driving signal, $`E`$) we will have maximum sensitivity to this change in phase. Furthermore it is desirable to have $`E`$ as large as possible so that small changes in phase translate into large changes in the quadrature.

We can now proceed to calculating the noise power spectrum for the measured current. The calculation is analogous to that for a double-sided cavity given in referencewam . We now do not work in the rotating frame but return to the laboratory frame. The Fourier component of the output operator for the current is given by

$$b_{\text{out}}(\omega )=\frac{\left[\left(\frac{\gamma _E\gamma _M}{2}\right)i(\omega _0\omega )i\delta \omega n_1\right]b_{\text{in}}(\omega )i\sqrt{\gamma _E}(\omega )+\sqrt{\gamma _E\gamma _M}a_{\text{in}}(\omega )}{\left[\left(\frac{\gamma _E+\gamma _M}{2}\right)+i(\omega _0\omega )+i\delta \omega n_1\right]}$$
(17)
where $`(\omega )`$ is the Fourier component of the driving amplitude. If the driving is noiseless and monochromatic, $`(\omega )=\delta (\omega \omega _d)`$. However in reality there would be some noise in the driving amplitude derived from the electrical noise in the driving circuit. We will treat this as entirely classical.

Equations 13 and 16 suggest that the signal will appear in the quadrature of the current out of phase with the driving force, defined by

$$X_{2,out}(t)=i(b_{\text{out}}^{}(t)b_{\text{out}}(t))$$
(18)
with Fourier components $`X_{2,out}(\omega )`$. The measured power spectrum is then given by the correlation function,

$$S_{2,out}(\omega ,\omega ^{})=X_{2,out}(\omega ),X_{2,out}(\omega ^{})$$
(19)
Using the specified states for the electronic and mechanical noise operators, we find,

$$S_{2,out}(\omega ,\omega ^{})=\left[|(\omega )|^2(2\overline{n}(\omega ,T_E)+1)+|𝒜(\omega )|^2(2\overline{n}(\omega ,T_M)+1)\right]\delta (\omega \omega ^{})$$
(20)
where

$`(\omega )`$ $`=`$ $`{\displaystyle \frac{\frac{\gamma _E\gamma _M}{2}i\left((\omega \omega _0)+\delta \omega n_1\right)}{\frac{\gamma _E+\gamma _M}{2}+i\left((\omega \omega _0)+\delta \omega n_1\right)}}`$
$`𝒜(\omega )`$ $`=`$ $`{\displaystyle \frac{\sqrt{\gamma _E\gamma _M}}{\frac{\gamma _E+\gamma _M}{2}+i\left((\omega \omega _0)+\delta \omega n_1\right)}}`$
To estimate the signal to noise ratio (SNR) we evaluate the spectrum at the driving frequency (that is to say, at the central Fourier component of the coherent driving);

$$S(\omega _0)=\frac{\left[\left(\frac{\gamma _E\gamma _M}{2}\right)^2+(\delta \omega n_1)^2\right](2\overline{n}(\omega ,T_E)+1)+\gamma _E\gamma _M(2\overline{n}(\omega ,T_M)+1)}{\left(\frac{\gamma _M+\gamma _E}{2}\right)^2+(\delta \omega n_1)^2}$$
(21)
Equations 13, 16 and 18 show that the magnitude of the Fourier component of the mean signal at the driving frequency is given by

$$|X_{2,out}(\omega _D)|=\frac{8\sqrt{\gamma _E}\delta \omega n_1}{(\gamma _M+\gamma _E)^2+4\delta \omega ^2n_1}$$
(22)
The signal is a sharp peak at $`\omega =\omega _d=\omega _0`$, in which there is a noise power $`S(\omega _0)`$ per root Hertz. So the SNR per root Hertz is $`|X_{2,out}(\omega _D)|^2/S(\omega _0)`$, or

$$\text{SNR}=\frac{16\gamma _E^2\delta \omega ^2n_1}{\left[(\gamma _M+\gamma _E)^2+4\delta \omega ^2n_1\right]\left[(\gamma _E\gamma _M)^2+4\delta \omega ^2n_1\right](2\overline{n}(\omega ,T_E)+1)+\gamma _E\gamma _M(2\overline{n}(\omega ,T_M)+1)}$$
(23)
If the SNR required for the measurement is $`\text{SNR}_r`$, then we must average over noise for a time $`t`$ such that $`\text{SNR}_r=\text{SNR}/\sqrt{t}`$. If we set $`n_1=1`$, so we are measuring the charge on one electron, the sensitivity is then $`e\sqrt{t}=e\text{SNR}_r/\text{SNR}`$.

## III Unconditional description of the measurement.

When we measure a quantum system, we bring an extremely large set of independent observables of our instrument and its environment into correlation with the measured system observable. The environment of the electroscope has two distinct components. Firstly there is the environment associated with the mechanical oscillator, which is responsible for mechanical damping and noise. Secondly there is the environment associated with the electrical readout, which is responsible for Johnson-Nyquist noise in the electrical circuit, and ultimately provides the measured result. However, we are interested in what the measurement tells us about the system, not in the exact quantum state of the instrument and its environment. Useful instruments must operate independently of the detailed state of their environments.

There are two ways to describe the partial state of the charge and oscillator. Firstly we can ignore the results of the measurement and average over states of the environment completely. In this case the evolution of the charge and oscillator is described by a master equation. Effectively we are averaging over the ensemble of partial states distinguished by different measurement records

Secondly, we can ask for the conditional states of the charge and oscillator, given a particular measurement record. Each member of the ensemble of partial states is associated with a distinct measurement record of the instrument. For it to be an effective measurement, observers must be able to distinguish the states of the instrument. In other words the charge must end up correlated with some simple macroscopic quantity, like the current in a wire or the position of a pointer on a scale. It is then possible to ask for the particular partial state of the measured system that is correlated with a known pointer value. In other words we need to be able to specify the conditional state of the system given a readout of the instrument variable that distinguishes different charge states. This is the conditional, or selective, description of the measured system. Of course if we average over the readout variables, we must obtain the unconditional description of the system.

We begin with the unconditional description of the measurement. The dominant sources of excess noise that limit the quality of the measurement are the thermal mechanical noise and thermal electrical noise on the readout circuit. Under certain Markoff and rotating wave assumptionscarmichael1 ; gardiner , the explicit states of the mechanical and electrical reservoirs may be traced out. This leaves the following master equation for the density operator of the composite system of charge and cantilever,

$`\dot{\rho }(t)`$ $`=`$ $`i[\mu (en_1)^2\widehat{c}^{}\widehat{c},\rho ]i[\left(\widehat{c}+\widehat{c}^{}\right),\rho ]`$ (25)
$`+{\displaystyle \underset{i=M,E}{}}\gamma _i(\overline{n}_i+1)𝒟[c]\rho +\gamma _i\overline{n}_i𝒟[c^{}]\rho `$
where the superoperator $`𝒟`$ is defined by

$$𝒟[c]\rho =c\rho c^{}\frac{1}{2}(c^{}c\rho +\rho c^{}c)$$
(26)
This can be written in a more standard form

$`\dot{\rho }`$ $`=`$ $`i[\mu (en_1)^2\widehat{c}^{}\widehat{c},\rho ]i[\left(\widehat{c}+\widehat{c}^{}\right),\rho ]`$ (28)
$`+\gamma (\overline{n}+1)𝒟[c]\rho +\gamma \overline{n}𝒟[c^{}]\rho `$
where $`\gamma \gamma _M+\gamma _E`$, and $`\overline{n}(\gamma _M\overline{n}(\omega ,T_M)+\gamma _E\overline{n}(\omega ,T_E))/\gamma `$.

We will begin solving this master equation by separating the dynamics of the cantilever and the charge. As before, we assume there is only one charge in the system, and consider the charge states $`|\text{0}`$ and $`|\text{1}`$. We can decompose $`\rho `$ into a $`2\times 2`$ matrix of cantilever operators

$$\rho =\widehat{A}|\text{0}\text{0}|+\widehat{B}|\text{1}\text{1}|+\widehat{Z}|\text{0}\text{1}|+\widehat{Z}^{}|\text{1}\text{0}|$$
(29)
Since $`\rho `$ is Hermititan, we need only three cantilever operators, $`\widehat{A}`$, $`\widehat{B}`$ and $`\widehat{Z}`$. We can now decompose Equation 28 into three independent equations involving only cantilever operators:

$$\frac{d\widehat{A}}{dt}=i[E(\widehat{c}+\widehat{c}^{}),\widehat{A}]+(\overline{n}+1)𝒟[c]\widehat{A}+\overline{n}𝒟[c^{}]\widehat{A}$$
(30)
$$\frac{d\widehat{B}}{dt}=i[E(\widehat{c}+\widehat{c}^{})+\mathrm{\Delta }\widehat{c}^{}\widehat{c},\widehat{B}]+(\overline{n}+1)𝒟[c]\widehat{B}+\overline{n}𝒟[c^{}]\widehat{B}$$
(31)
$$\frac{d\widehat{Z}}{dt}=i[E(\widehat{c}+\widehat{c}^{}),\widehat{Z}]+i\mathrm{\Delta }\widehat{Z}\widehat{c}^{}\widehat{c}+(\overline{n}+1)𝒟[c]\widehat{Z}+\overline{n}𝒟[c^{}]\widehat{Z}$$
(32)
As before, we are now measuring time relative to the damping time $`1/\gamma `$.

If we measured the state of the charge by means other than the cantilever, the state of the cantilever immediately after the measurement would be $`\widehat{B}`$ if the charge were present, or $`\widehat{A}`$ if it were absent. Hence $`\widehat{A}`$ and $`\widehat{B}`$ must be density operators, and Equations 30 and 31 have the form of master equations for a damped harmonic oscillator. Such equations, and their solutions, are familiar to quantum opticians. The stable solution is a displaced thermal state, which can be written

$$\rho =\left(1e^{\lambda (t)}\right)D(\alpha (t))e^{\lambda (t)\widehat{c}^{}\widehat{c}}D^{}(\alpha (t))$$
(33)
, where $`D(\alpha )`$ is a displacement operator $`\mathrm{exp}(\alpha \widehat{c}^{}\alpha ^{}\widehat{c})`$, and in the steady state $`\lambda =\mathrm{}\omega _0/k_bT`$. In the limit of low temperature, $`kT\mathrm{}\omega _0`$, this becomes a coherent state $`|\alpha \alpha |`$. In the steady state, the cantilever has as many thermal phonons as a resevoir mode with the same frequency, i.e. $`e^\lambda =\overline{n}/(\overline{n}+1)`$. Its coherent amplitude $`\alpha _0`$ reaches a balance with the driving and damping after a time around $`2/\gamma `$:

$$\alpha (t)=\alpha _0e^{\kappa t/2}\frac{2iE}{\kappa }(1e^{\kappa t/2})$$
(34)
$$\kappa =\{\begin{array}{cc}1& n=0\\ 1+2i\mathrm{\Delta }& n=1\end{array}$$
(35)
During measurement, the cantilever states $`\widehat{A}`$ and $`\widehat{B}`$ are displaced thermal states with distinct coherent amplitudes.

As the measurement proceeds, we expect the charge state to evolve from a coherent superpostion of $`|\text{0}`$ and $`|\text{1}`$ to an incoherent mixture; in terms of our decomposition, we expect the off-diagonal term $`Z`$ to decay with time. An operator of the form

$$Z=z(t)D(\alpha )\mathrm{exp}(\lambda \widehat{c}^{}\widehat{c})D^{}(\beta )$$
(36)
where $`z(t)`$ is a (possibly complex) amplitude, solves Equation 32 if $`\alpha `$, $`\beta `$, $`\lambda `$ and $`z`$ obey the following differential equations:

$$\frac{dl}{dt}=(\overline{n}+1)l^2(2\overline{n}+1+i\mathrm{\Delta })l+\overline{n}$$
(37)
$$\frac{da}{dt}=\left(i\mathrm{\Delta }+(\overline{n}+1)l\overline{n}\frac{1}{2}\right)aiE(1l)$$
(38)
$$\frac{db}{dt}=\left((\overline{n}+1)l\overline{n}\frac{1}{2}\right)b+iE(1l)$$
(39)
$$\frac{dk}{dt}=iE(ab)+(\overline{n}+1)(l+ab1)+1$$
(40)
Here $`l=\mathrm{exp}(\lambda )`$, $`a=\alpha l\beta `$, $`b=\beta ^{}l\alpha ^{}`$, and $`k=\mathrm{log}z+l\alpha ^{}\beta \frac{1}{2}(|\alpha |^2+|\beta |^2)`$

In general, these equations can be solved numerically. However, there are some special cases where we can get interesting information analytically. First we consider the zero temperature limit, where the off diagonal term $`\widehat{Z}`$ is a projector $`z|\alpha \beta |`$. The amplitudes $`\alpha `$ and $`\beta `$ are the amplitudes of the diagonal terms given by Equation 35, and $`z`$ is a complex amplitude. Once $`\alpha `$ and $`\beta `$ have reached their steady state, the trace of the off-diagonal term, decays exponentially with a rate $`|\alpha \beta |^2/2`$.

If we assume the detuning $`\mathrm{\Delta }`$ is small, and hence $`|\beta ||\alpha |=2E`$, The difference between the steady state amplitudes of Equation 35 is

$$|\alpha \beta |^2=\frac{16E^2\mathrm{\Delta }^2}{1+4\mathrm{\Delta }^2}4|\alpha |^2\mathrm{\Delta }^2$$
(41)
Cleland and Roukes give enough information about their devices for us to calculate this explicitly cr98 . Using the data in Table 1, we can calculate $`\alpha `$ from the definition of the annihilation operator for a torsional pendulum

$$\alpha =\widehat{c}=\sqrt{\frac{\kappa }{2\mathrm{}\omega _0}}\theta _{\text{max}}=5.3\times 10^6$$
(42)
The normalised detuning can be calculated from the frequency shift per electron and the measured quality factor:

$$\mathrm{\Delta }=\delta \omega /\gamma =\frac{2\pi \delta \nu Q}{\omega _0}=2.4\times 10^4$$
(43)
The decoherence rate is then $`3.2\times 10^6\gamma `$, or $`8.1\times 10^9\text{s}^1`$.

As $`\overline{n}`$ increases from zero, the amplitudes $`\alpha `$ and $`\beta `$ for the off diagonal operator $`\widehat{Z}`$ are reduced, as shown in Figure 2. The initial decay of $`z(t)`$ is shown in Figure 4, and the steady decay rate, i.e. the limit of $`|z^{}(t)/z(t)|`$ when $`t1/\gamma `$, in Figure 3. At low temperatures (below $`130\text{mK}`$), the increased thermal noise from the bath causes $`\widehat{Z}`$ to decay more rapidly as the temperature of the bath is increased. Contrary to expectations, the steady decoherence rate of the charge superposition decreases as the bath temperature increases above $`130\text{mK}`$. The extra thermal noise increases the overlap between the oscillator states corresponding to the presence and absence of charge.

## IV Conditional description

We now turn to the correlations between the charge and the resevoir system. These are important because we must be able to distinguish the results corresponding to different charges to make a measurement of the charge at all. They can be studied most simply using quantum trajectory theory, which associates charge states with possible observed states of the apparatus carmichael-lect .

We will assume we monitor the current in the electrical resevoir; this is equivalent to an optical homodyne measurement sa96 . The inferred state of the charge as such a measurement proceeds is governed by a Wiener process, which is generated by a stochastic increment $`dW`$. The average of $`dW`$ over the ensemble of possible measurement results is zero. Since the deviation of the Wiener process represented by $`dW`$ increases proportional to $`\sqrt{t}`$, the average of $`(dW)^2`$ is $`dt`$. The simplest way to manipulate such differentials is to modify the chain rule, to give what is know as Ito calculus.

Given a particular measurement result, labelled by a Wiener increment $`dW`$, the evolution of the charge and cantilever is

$$d|\psi =\left(\frac{1}{i\mathrm{}}\widehat{H}dt\frac{\gamma }{2}\left(c^{}c2\frac{x}{2}c+\frac{x}{2}^2\right)dt+\sqrt{\gamma }\left(c\frac{x}{2}\right)dW\right)|\psi $$
(44)
When we insert the charge and cantilever Hamiltonian, and normalise time by the damping rate as before, this becomes

$$d|\psi =\left(i(E(c+c^{})+\mathrm{\Delta }nc^{}c)dt\frac{1}{2}\left(c^{}c2\frac{x}{2}c+\frac{x}{2}^2\right)dt+\left(c\frac{x}{2}\right)dW\right)|\psi $$
(45)
When a particular function $`dW`$ is selected from the Wiener ensemble, this can be solved to show the evolution of a pure state $`|\psi `$. These states form an ensemble with density operator $`\rho `$. Of course $`\rho `$ can be decomposed into many ensembles, so the evolution generated by Equation 45 is not unique. The details are given in Carmichaelcarmichael-lect .

Mixed states of the cantilever and charge must be written in the form of Equation 29. However pure states can always be written as

$$|\psi =|A|\text{0}+|B|\text{1}$$
(46)
as before we will assume the state of the cantilever is initially coherent, so that

$$|\psi =p|\alpha 0+q|\beta 1$$
(47)
The differential of a scaled coherent state $`q(t)|\beta (t)`$ is

$$d(q|\beta )=\left(dq\frac{1}{2}qd|\beta |^2\right)|\beta +q\dot{\beta }dtc^{}|\beta $$
(48)
comparison with Equation 45 gives Equations 34 for the evolution of $`\alpha `$ and $`\beta `$ as before. Some Ito calculus manipulations show that

$$d|q|^2=|pq|^2(x_\alpha x_\beta )dW$$
(49)
where $`x_\alpha `$ is the expectation value of the amplitude quadrature $`x`$ in a coherent state $`|\alpha `$, which is just $`2\text{Re}\alpha `$. The normalisation of $`|\psi `$ requires that $`d|p|^2=d|q|^2`$.

We need to compare the gain in knowledge shown by this trajectory picture to the decay of coherence modelled by the master equation. The results of the measurement are the probabilities $`|p|^2`$ and $`|q|^2`$; the pure state which the observer will infer from these has a density operator

$$\rho =|p|^2|00|+|q|^2|11|+|pq|(|10|+|01|)$$
(50)
The off-diagonal terms in this have magnitude $`|pq|`$; we can average over $`dW`$ to see the behaviour of the density operator for the ensemble of measurement results.

Some more routine Ito calculus gives the evolution of this:

$$d|pq|=|pq|\left(\frac{1}{8}\left(x_\beta x_\alpha \right)^2dt+\frac{1}{2}\left(|p|^2|q|^2\right)\left(x_\beta x_\alpha \right)dW\right)$$
(51)
Since the average of $`dW`$ over different measurement results is zero, on average

$$d|pq|=\frac{1}{8}(x_\beta x_\alpha )^2|pq|dt$$
(52)
If the difference between the charges associated with states $`|0`$ and $`|1`$ is $`e`$, then in the state $`p|\alpha 0+q|\beta 1`$, the uncertainty in the charge is given by

$$\left((ne)^2ne^2\right)^{{\scriptscriptstyle \frac{1}{2}}}=e|pq|$$
(53)
From Equation 52, this decreases exponentially as the measurement progresses, at a rate

$$\frac{1}{8}(x_\beta x_\alpha )^2=\frac{8\mathrm{\Delta }^2E^2}{(1+4\mathrm{\Delta }^2)^2}$$
(54)
This differs from the square root decay of classical uncertainty as measurements are averaged over time, but exponential decay is what we would expect for decay of coherencewam . For the device described in cr98 , this is almost equal to the decoherence rate. In real devices, thermal noise will cause the trajectory states to be mixed, however the evolution of such mixed states is much harder to calculate.

## V Discussion

To estimate the time required for our measurement, we have calculated how long it takes for an initially pure superposition of charge states to be reduced to a mixture, and how long (in some sense) it takes us to find out which charge eigenstate we have been left with. While these questions are interesting in their own right (they composed the deepest mystery of physics for the best part of a century), it could be argued that they don’t reflect the way measurements would be used in a real computer.

The most that we could do with measurements on pure states is state preparation. In a coherent quantum computer this would be rather pointless though, since if we know the initial state we could just rotate it into the eigenstate we want. We carry out measurements to find out something we don’t know: in other words we apply them to mixed states, with a view to finding out which of the possibilities is real.

Information theory provides tools to quantify this, such as conditional entropy and mutual information. Unfortunately calculating any of these requires knowledge of the ensemble of trajectories generated by each component of the mixture, and the overlaps between them. In general it is hard to find the probability distribution of trajectories; we usually just calculate averages. It might be worth doing this numerically, however.

There is a more straightforward limitation to our analysis: in present day devices the thermal effects that we have neglected in the trajectory treatment utterly dominate the vacuum noise we have considered. Hence the measurement time will be limited by the need to average classical fluctuations. It is possible that future devices operating at higher frequencies will reduce the level of thermal noise so that quantum effects will be important. This presents the remarkable prospect of a solid cantilever with position and momentum known to the limit allowed by the uncertainty principle.


format: text

text:
# Shearing a Glassy Material: Numerical Tests of Nonequilibrium Mode-Coupling Approaches and Experimental Proposals

\[

## Abstract

The predictions of a nonequilibrium schematic mode-coupling theory developed to describe the nonlinear rheology of soft glassy materials have been numerically tested in a sheared binary Lennard-Jones mixture. The theory gives an excellent description of the stress/temperature ‘jamming’ phase diagram of the system both at the microscopic and the macroscopic levels. In the present paper, we focus more particularly on the issue of an effective temperature $`T_{\mathrm{eff}}`$ for the slow modes of the fluid, as defined from a generalized fluctuation-dissipation theorem. As predicted theoretically, many different observables are found to lead to the same value of $`T_{\mathrm{eff}}`$, suggesting several experimental procedures to measure $`T_{\mathrm{eff}}`$. New, simple experimental protocols to access $`T_{\mathrm{eff}}`$ from a generalized equipartition theorem are also proposed, and one such experiment is numerically performed. These results give strong support to the thermodynamic interpretation of $`T_{\mathrm{eff}}`$ and make it experimentally accessible in a very direct way.

\]

Glassy materials are usually defined by the fact that their relaxation time is larger than the experimental time scale. In simple molecular systems, the associated glass transition temperature corresponds to very high viscosities, making it difficult to investigate experimentally their rheological properties. In complex fluids (e.g. colloids, emulsions) it is however possible to reach a glassy situation, in the sense of large relaxation times, with systems having viscosities or shear moduli that allow for rheological investigations . Such materials have been described as ‘soft glassy materials’ .

In its glassy state, a material is by definition out of equilibrium. Physical properties are then a function of the time $`t_w`$ spent in the glassy phase, a behavior called aging . Interestingly, recent experiments on various complex fluids have demonstrated striking similarities with other, more standard, glassy systems . Upon imposing a steady, homogeneous shear flow, a different kind of nonequilibrium situation is obtained. The flow, characterized by the shear rate $`\gamma `$, creates a nonequilibrium steady state, in which time translation invariance is recovered . This situation can therefore be used to probe the glassy state, with the convenient feature of having the shear rate $`\gamma `$ rather than the waiting time $`t_w`$ as a control parameter. Moreover, this way of probing the nonequilibrium properties of glassy systems is probably more relevant experimentally than the aging approach, at least in the case of soft glassy materials.

Recently, a general scenario was proposed for glassy systems subject to an external forcing , based on the study of mean-field models. The rationale of this approach is that the equilibrium dynamics of these models is equivalent to the ‘schematic’ mode-coupling approach of slowing down in supercooled liquids . The study of their nonequilibrium dynamics can thus be seen as a nonequilibrium schematic mode-coupling approach . To our knowledge, the mode-coupling theory of supercooled fluid has not yet been extended to fluid under shear beyond the linear response regime. However, in analogy to what was done for aging or supercooled systems, it is sensible to bypass this aspect and to carry out a direct comparison between mean-field predictions and experimental or numerical results.

The aim of this work is thus to check on a realistic model of a fluid the predictions that emerged from the theoretical approach of Ref. . Several earlier studies have been devoted to simulating glassy sheared fluids. Ref.  focused on the dynamics at the molecular level. Our study is devoted to more global aspects, with the aim of providing experimentally testable predictions. Ref.  also proposed to use the shear rate as a control parameter for jamming systems. Some dynamic aspects of a model similar to ours, with the difference of being athermal (zero temperature) have been investigated, as discussed in a companion paper .

We have investigated the stress/temperature jamming phase diagram (Fig. 1) of the ‘standard model’ for supercooled liquids, namely a 3D 80:20 binary Lennard-Jones mixture, in a simple shear flow defined by $`𝒗=\gamma y𝒆_𝒙`$. The system consists of $`N=2916`$ particles in a cubic simulation box. It has been characterized in much details at the reduced density $`\rho =1.2`$, where we carry out our simulations, both in supercooled  and aging  regimes. Our simulations follow the protocol detailed in Ref. : the system is first made stationary on a timescale of a few $`\gamma ^1`$. Then we perform our measurements in a range $`T[0.15,0.6]`$ and $`\gamma [10^4,10^1]`$. Standard Lennard-Jones units  are used.

Following the theoretical analysis of Ref. , our approach consists in three main steps. We first investigate macroscopic rheological aspects (flow curves). We then focus on microscopic properties by an analysis of the density fluctuations: structure factor, density-density correlation functions. Both aspects will be discussed in much more detail in a longer publication . Last, we analyse the main new feature for the field of rheology which emerged in Ref. , namely the existence, behavior and properties of an effective temperature, defined through a generalized fluctuation-dissipation theorem . This point, together with its experimental consequences are the main object of this paper.

The macroscopic behavior of the system under shear is summarized in a stress/temperature phase diagram in Fig. 1. Our results are qualitatively similar to the one reported in Ref.  on a different glass-former. A Newtonian behavior is observed in the high-$`T`$, low-$`\sigma `$ part of the phase diagram. This corresponds roughly to situations where $`\gamma ^1`$ is larger than the relaxation time of the fluid, which is thus only weakly affected by the flow. Outside this region, we find that the viscosity of the fluid decreases when the shear rate increases, a shear-thinning behavior reported in various ‘soft glassy materials’ . We used the phenomenological relation

$$\sigma \sigma _0+a\gamma ^n.$$
(1)
to extract the behavior of the yield stress $`\sigma _0(T)`$, reported in Fig. 1. This ‘jamming transition’ line was recently experimentally investigated . Although our results cover 3 decades in shear rate, we cannot report a definitive functional form for the shear-thinning of the viscosity  which can also be satisfactorily described by $`\sigma _00`$, which amounts to describe the system as a power-law fluid, $`\eta \gamma ^{n1}`$. In the supercooled regime we find $`n1/3`$, while at lower temperatures, the exponent $`n`$ is temperature dependent with $`n0`$ when $`T0`$, indicating that a yield stress could exist in this limit only . The latter power-law behavior is precisely the one predicted theoretically in Ref. .

Equilibrium mode-coupling theory gives rise to several quantitative predictions regarding the scaling properties of the intermediate scattering function,

$$C_𝒌(t)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}\left(i𝒌\left[𝒓_𝒋(t+t_0)𝒓_𝒋(t_0)\right]\right),$$
(2)
when the glassy phase is approached by lowering $`T`$ with $`\sigma =0`$. It was shown in Ref.  that similar scaling properties are expected when the glassy phase is approached by lowering the shear stress $`\sigma `$ at constant temperature (vertical line in Fig. 1). We have tested in a detailed way these predictions in our simulation . Here, we only report the validity of the theoretically predicted ‘time-shear superposition property’  in Fig. 2, which shows that the slow decay of $`C_𝒌(t)`$ has the scaling form $`C_𝒌(t)F(t/t_{\mathrm{rel}})`$, where the relaxation time is defined, as usual, by $`C_𝒌(t_{\mathrm{rel}})e^1`$. The scaling function $`F(x)`$ is well described by a stretched exponential, $`F(x)\mathrm{exp}(x^\beta )`$, with an exponent $`\beta `$ which increases from $`\beta 0.77`$ for $`TT_c`$ to the value $`\beta =1`$ as $`T0`$.

We now focus on the issue of an effective temperature for the slow modes of the sheared fluid. This quantity, naturally included in nonequilibrium mode-coupling theories , is defined through a nonequilibrium generalization of the fluctuation-dissipation theorem (FDT) . Consider two physical observables, $`O(t)`$ and $`O^{}(t)`$, their connected cross-correlation function $`C_{OO^{}}(t)O(t+t_0)O^{}(t_0)O(t_0)O^{}(t_0)`$, and the response function $`R_{OO^{}}(t)\frac{\delta O(t+t_0)}{\delta h_O^{}(t_0)}`$, where $`h_O^{}`$ is the field thermodynamically conjugated to $`O^{}(t)`$. At equilibrium, both quantities satisfy the FDT, $`TR_{OO^{}}(t)=\frac{\mathrm{d}C_{OO^{}}(t)}{\mathrm{d}t}`$. The susceptibility $`\chi _{OO^{}}(t)_0^tdt^{}R_{OO^{}}(t^{})`$ can be measured by applying a small (to ensure linear response), constant field $`h_O^{}`$ between times $`0`$ and $`t`$ and FDT implies a simple linear relation $`T\chi _{OO^{}}(t)=\left(C_{OO^{}}(0)C_{OO^{}}(t)\right)`$. In the sheared fluid, an effective temperature is defined by

$$R_{OO^{}}(t)=\frac{1}{T_{\mathrm{eff}}^{OO^{}}(C_{OO^{}})}\frac{\mathrm{d}C_{OO^{}}(t)}{\mathrm{d}t},$$
(3)
where the $`T_{\mathrm{eff}}^{OO^{}}(x)`$ are a priori arbitrary functions of their argument, which may in general depend on the observables $`O`$ and $`O^{}`$ under study. For any pair of observables, $`T_{\mathrm{eff}}^{OO^{}}(x)`$ can be measured in the sheared system, by following the same linear response procedure as above, so that $`\chi _{OO^{}}(t)=_{C_{OO^{}}(t)}^{C_{OO^{}}(0)}\frac{\mathrm{d}x}{T_{\mathrm{eff}}^{OO^{}}(x)}.`$ The existence of an effective temperature is thus demonstrated if a straight line is obtained in a susceptibility-correlation plot parameterized by the time. Obviously, the introduction of an effective temperature is of ‘thermodynamic’ interest only if this quantity is actually independent of the observables under consideration, $`T_{\mathrm{eff}}^{OO^{}}T_{\mathrm{eff}}`$. This crucial feature is true at the mean-field level , and we shall prove that it is nicely satisfied in our model.

We have already shown the existence of such an effective temperature for the slow modes of the system for a single correlation function, Eq. (2), at a given wavevector . We also found that the essential phenomenological idea that a system sheared more vigorously has a higher $`T_{\mathrm{eff}}`$ is indeed captured by this definition .

Here, we go much further in our investigations and prove that several different observables lead to the same value of $`T_{\mathrm{eff}}`$. Since the numerical measurement of $`T_{\mathrm{eff}}`$ is very demanding, our strategy has been to compute for $`T=0.3`$ and $`\gamma =10^3`$ the value of $`T_{\mathrm{eff}}`$ from the observables of Ref.  with a great accuracy. We found $`T_{\mathrm{eff}}0.65`$. We then computed $`T_{\mathrm{eff}}`$ using different observables, and checked that the value so obtained was compatible with $`T_{\mathrm{eff}}=0.65`$. Our results are summarized in Fig. 3 which shows 10 among the 14 different susceptibility-correlation measurements we have performed. All our data are well compatible with the single value of $`T_{\mathrm{eff}}=0.65`$ for the slow modes of the fluid.

We first investigated the density fluctuations as an observable, taking $`O(t)=\frac{1}{N}_{j=1}^N\epsilon _j\mathrm{exp}\left(i𝒌𝒓_𝒋(t)\right)`$, and $`O^{}(t)=2_{j=1}^N\epsilon _j\mathrm{cos}\left(𝒌𝒓_𝒋(t)\right)`$. For some wavevectors, this was done separately for both types ($`A`$ and $`B`$) of Lennard-Jones particles (noted ‘$`A+B`$’ in Fig. 3). Taking $`\epsilon _j=1`$ selects the coherent part of the intermediate scattering function (noted ‘coh’ in Fig. 3), while $`\epsilon _j=\pm 1`$ selects the incoherent one. Such response-correlation plots could be obtained experimentally in soft condensed matter systems. While the correlation functions are reasonably easily obtainable through light scattering experiments, the same is not true of response functions. To obtain the latter, one has to manipulate the particles through some externally applied potential, modulated at the same wavevector as used in the light scattering experiment. One suggestion would be to use some non-index matched tracer particles in an index matched colloidal suspension. The tracer particles would then be sensitive to the intensity of the local electric field, as in optical tweezers. An interference pattern would actually realize the modulated external potential considered here. Reading of the response could then be obtained from the scattering at the wavelength corresponding to this pattern.

We also used as a correlation the mean square displacement of a tagged particle. The associated response function is the displacement induced by applying a small, constant external force to this tagged particle . Both quantities are linked by a FDT, the Einstein relation. We computed both quantities separately for particles of type $`A`$ and $`B`$ (noted ‘Self $`A`$’ and ‘Self $`B`$’ respectively). Its interest, especially in view of experimental realizations, lies in the fact that the full time dependence of the correlation and response functions is not needed to extract $`T_{\mathrm{eff}}`$. Indeed, at large times, both quantities become proportional to the time which defines the diffusion constant and the mobility. One may therefore define $`T_{\mathrm{eff}}`$ simply as the ratio of diffusion to mobility. Again, experiments could be considered if tracer particles sensitive to an external force field (e.g. magnetic particles) could be introduced into the system, and their mobility measured together with their diffusion constant.

A completely different observable, relevant to flow situations, is the stress $`\sigma `$. We have studied the case in which the observables $`O`$ and $`O^{}`$ are equal to the diagonal stress in the direction transverse to the flow, $`\sigma _{zz}(t)`$. To add a field conjugated to $`\sigma _{zz}`$, a compression $`\delta L_z`$ of the simulation box is realized by rescaling all particle coordinates at time $`t=0`$, in analogy to stress relaxation experiments. The corresponding curve is labelled ‘Stress’ in Fig. 3. Experimentally, the off-diagonal component of the stress would be used as the observable. Preliminary results in this direction, using an extremely sensitive rheometer, have been obtained in aging colloidal systems .

Definition (3) of $`T_{\mathrm{eff}}`$ implies that, if the fluid is used as a thermal bath to equilibrate a subsystem (a ‘thermometer’) of typical time scale $`t_st_{\mathrm{rel}}`$, the thermometer does not measure the microscopic temperature, but rather the effective temperature associated with its characteristic time scale . We propose here to use tracers of mass $`m_{\mathrm{tr}}`$ as a thermometer, since tuning $`m_{\mathrm{tr}}`$ allows to control their vibration time scale $`t_s`$. We have considered 10 massive tracers with $`m_{\mathrm{tr}}[1,10^7]`$, but being otherwise identical to $`A`$ particles. Since $`t_s\sqrt{m_{\mathrm{tr}}}`$, heaviest particles have a frequency typically $`10^3`$ times smaller than the light ones, implying $`t_st_{\mathrm{rel}}`$. Reading of the temperature is done by measuring the average mean square velocity of the tracers in the direction $`z`$. Results are shown in Fig. 4, which shows that light particles measure the bath temperature while heaviest particles measure $`T_{\mathrm{eff}}`$. This implies that a generalized equipartition theorem holds,

$$\frac{1}{2}m_{\mathrm{tr}}v_z^2=\frac{1}{2}T_{\mathrm{eff}}.$$
(4)
This result could be tested against experiments involving for instance colloidal particles, the tracers, in a complex fluid, e.g. polymeric, or by investigating rotational degrees of freedom of non-spherical tracer particles.

This last result opens the way for new, simple determinations of $`T_{\mathrm{eff}}`$ in out of equilibrium glassy materials. Indeed, no ‘complex’ dynamic functions such as correlations or susceptibilities are needed here. In this context, it would be interesting to reproduce Perrin’s experiment on barometric equilibrium of colloidal suspensions . We expect indeed that the barometric equilibrium of heavy particles inside a horizontally sheared fluid should be ruled by $`T_{\mathrm{eff}}`$ instead of the room temperature.

All these results give strong support to the theoretical scenario elaborated from mean-field theories to describe the rheology of soft glassy materials . Note in particular that $`T_{\mathrm{eff}}`$, which we have shown to be a physically relevant quantity, is a natural outcome of the theory. While numerical simulations provide a test of the theory on short time scales, systematic experiments that could test quantitatively existing theories of nonequilibrium glassy dynamics are still needed. In that sense, the situation is more or less similar to the mid-eighties when schematic equilibrium mode-coupling were already derived, but with little experimental confirmation of its main features. This is why we tried, as much as possible, to suggest experimental counterparts to our numerical measurements. We hope that our findings and suggestions will motivate further experimental work in the field.

We thank L. Bocquet, J. Kurchan and W. Kob for discussions. This work was supported by the PSMN at ENS Lyon and the CDCSP at Université de Lyon.


text:
# Status of the QCDOC projectPresented by NHC and TW at Lattice 2001, Berlin.

## 1 INTRODUCTION

The debate over which computing platforms are suitable for the simulation of lattice gauge theories has continued at this year’s lattice conference. It has been shown that recent off-the-shelf PC processors by Intel or AMD are capable of delivering impressive floating-point performance on the order of 1 GFlops if the vector processing units are employed and cache management is optimized . Unfortunately, this single-node performance cannot be fully utilized in a PC cluster if the local lattice volume becomes small. In this case, the surface-to-volume ratio is large, and the communications latency inherent in standard network solutions such as Ethernet or Myrinet slows down the individual processors. Therefore, to achieve reasonable efficiencies with a PC cluster, the local lattice volume must not be too small. This implies that for a given total lattice volume, the number of nodes working on a single problem is limited.

However, to finish a dynamical fermion calculation on a moderately large lattice in a reasonable amount of time, it is essential to distribute the total volume onto as many nodes as possible, which implies very small local lattice volumes, perhaps as small as $`2^4`$. This requires communications between neighboring nodes with extremely low latencies that cannot be achieved using off-the-shelf networking components. Massively parallel machines with custom-designed communications hardware (with latencies 10–100 times smaller than in the case of Myrinet) appear to be the only viable alternative. This is especially true in the field of lattice gauge theory with its very regular communications requirements. In this contribution, we describe the design of the QCDOC architecture which is capable of delivering computing power in the tens of TFlops range at a price/performance ratio of 1 US-$ per sustained MFlops.

## 2 OVERVIEW OF QCDOC

A first description of the QCDOC project has been given at last year’s lattice conference . In the nomenclature of parallel computers, this is a multiple-instruction, multiple-data (MIMD) machine with distributed memory. The name QCDOC stands for “QCD On a Chip”. An essential element of this project is the collaboration with IBM Research and the resulting ability to utilize IBM’s system-on-a-chip technology. Nowadays it is possible to integrate most of the components that make up a single processing node on a single chip, creating an application-specific integrated circuit, or ASIC. The QCDOC chip is such an ASIC, consisting of

* 500 MHz, 32-bit PowerPC 440 processor core
* 64-bit, 1 GFlops floating-point unit
* 4 MBytes on-chip memory (embedded DRAM)
* controllers for embedded and external memory
* nearest-neighbor serial communications unit with estimated latencies of 120ns (300 ns) for send (receive), overlapped between the 12 independent directions and an aggregate bandwidth of 12 Gbit/s
* other components such as Ethernet controller, interrupt controller, etc.

The chip will occupy a die size of about $`(12\mathrm{m}\mathrm{m})^2`$, and the power consumption is expected to be in the range of 1–2 W.

Two such ASICs will be mounted on a daughterboard, together with two industry-standard double data rate (DDR) SDRAM modules (one per ASIC) whose capacity will be determined by the price of memory when the machine is assembled. (At the time of this writing, 256 MByte registered DIMMs with ECC were available for $40.) 32 daughterboards will be mounted on a motherboard, and 8 motherboards in a crate with a single backplane. The final machine then consists of a certain number of such crates.

There are two separate networks: the physics network and an Ethernet-based auxiliary network. The physics network consists of high-speed serial links between nearest neighbors with a bandwidth of 2$`\times `$500 Mbits/s per link. The nodes are arranged in a 6-dimensional torus which allows for an efficient partitioning of the machine in software, as described in detail in Ref. . The serial-communications unit in the QCDOC ASIC provides direct memory access, single-bit error detection with automatic resend, and a low-latency store-and-forward mode for global operations. The auxiliary network is used for booting, diagnostics, and I/O over Ethernet, with an Ethernet controller integrated on the ASIC. The Ethernet traffic to and from the ASIC will run at 100 Mbit/s. Hubs or switches on the motherboard will provide a bandwidth of 1.6 Gbit/s off a motherboard to commercial switches and the host workstation.

## 3 THE QCDOC ASIC

A schematic picture of the QCDOC ASIC is shown in Fig. 1. In this section we describe the various components in more detail. Many of the ingredients of the ASIC are provided by IBM, see Secs. 3.13.2. The major custom-designed components are the serial communications unit (SCU) and the prefetching eDRAM controller (PEC), see Secs. 3.33.4.

### 3.1 IBM Library Components

* There are three main busses:

+ The processor local bus (PLB) is a 128-bit wide, fully synchronous bus that runs at 1/3 of the CPU frequency. It has separate read and write busses, supports pipelining (up to 4 deep), and has a large number of sophisticated features such as split transactions, burst transfers, etc.
    + The on-chip peripheral bus (OPB) is a 32-bit wide, synchronous bus running at 1/2 of the PLB frequency. Its basic purpose is to off-load slower devices from the PLB bus.
    + The device control register bus (DCR) is a slow and simple bus used to read and write control registers. The devices on this bus are connected in a ring-like fashion. The 440 is the only master on this bus.
* The 440 PowerPC processor core is a dual-issue, superscalar, 32-bit implementation of the IBM Book-E architecture (the E stands for embedded). The core has a 7 stage pipeline. There are instruction and data caches of size 32 kBytes each. These are 64-way associative, partitionable, and lockable. Hardware memory protection is provided through a translation-lookaside buffer (TLB). The 440 has three PLB master interfaces, one each for instruction read, data read, and data write. The 440 will run at 500 MHz and is connected to a 64-bit, 1 GFlops IEEE floating-point unit (for details of the FPU see Ref. ). An important feature of the 440 is that it contains a JTAG interface. JTAG (Joint Test Action Group) is an industry-standard protocol that allows an external device to take complete control of the processor. This functionality will be used for booting and debugging, see below.
* The PLB arbiter provides programmable arbitration for the up to eight allowed masters that can control PLB transfers. In our design we have seven masters: the three PLB master interfaces of the 440 (two of which, namely data read and data write, will be channeled through the PEC, see below), the MCMAL, two master interfaces for the send and receive operations of the SCU and one for the DMA unit.
* The 4 MBytes of embedded DRAM (eDRAM) are a standard IBM library component that will be accessed with low latency and high bandwidth through a custom-designed controller, the PEC, see Sec. 3.4 below.
* The Universal Interrupt Controller (UIC) processes the various interrupts generated on- and off-chip and provides critical and non-critical interrupt signals to the 440.
* The DDR controller is a slave on the PLB, capable of transferring data to/from external DDR SDRAM at a peak bandwidth of 2.6 GBytes/s. It supports an address space of 2 GBytes and provides error detection, error correction, and refresh of the off-chip SDRAM.
* The PLB-OPB bridge is used to transfer data between the two busses. It is the only master on the OPB and a slave on the PLB.
* The Ethernet media access controller (EMAC) provides a 100 Mbit/s Ethernet interface. The media-independent interface (MII) signals at the ASIC boundary are connected to a physical layer chip on the daughterboard. The EMAC is a slave on the OPB and has sideband signals to the MCMAL on the PLB.
* The DMA-capable Memory Access Layer (MCMAL) loads/unloads the EMAC through the sideband signals. It is a master on the PLB.
* The inter-integrated circuit (I2C) controller is a slave on the OPB used to read the configuration of the DDR DIMM. This information is then used to configure the DDR controller.
* The general purpose I/O (GPIO) unit is another slave on the OPB whose 32-bit wide data bus will be taken out to the ASIC boundary. It can be used, e.g., to drive LEDs.
* The high-speed serial links (HSSL) used by the SCU each contain eight independent ports (four each for send/receive) through which bits are clocked into/out of the ASIC at 500 Mbit/s per port. The bits are converted to bytes (or vice versa) in the HSSL. The HSSL input clocks are phase-aligned by another IBM macro, the phase-locked loop (PLL).

### 3.2 Ethernet-JTAG Interface

As mentioned above, the 440 core has a JTAG interface over which one can take complete control of the processor. In particular, this interface can be used to load boot code into the instruction cache and start execution. This completely eliminates the need for boot ROM. The question is how the JTAG instructions should be loaded into the 440. (There are special tools that use the JTAG interface, but it would be impractical to connect one tool per ASIC for booting.) A solution to this question has already been developed at IBM Research, implemented using a field-programmable gate array (FPGA), that converts special Ethernet packets to JTAG commands and vice versa. This logic is now part of the QCDOC ASIC and will be used not only for booting but also to access the CPU for diagnostics/debugging at run time. The unique MAC address of each ASIC is provided to the Ethernet-JTAG component by location pins on the ASIC, and the IP address is then derived from the MAC address.

### 3.3 Serial Communications Unit (SCU)

The task of the SCU is to reliably manage the exchange of data between neighboring nodes with minimum latency and maximum bandwidth. The design takes into account the particular communication requirements of lattice QCD simulations. A schematic picture of the SCU is shown in Fig. 2.

The custom protocol governing the data transfers defines packets that contain a 64-bit data word and an 8-bit header containing control and parity bits. When the receive unit receives a packet, it first interprets the header, buffers the bytes from the HSSL, and assembles the 64-bit word. It then transfers the word to the receive register or passes it on to the send unit. The receive buffer can store three 64-bit words so that the send unit (in a neighboring node) can send three words before an acknowledgment is received. The functionality of the send unit is essentially the inverse of that of the receive unit. Send and receive operations can proceed simultaneously. Each send or receive unit is controlled by a DMA engine which then transfers the data between memory and a send/receive register. Each DMA engine is controlled by block-strided-move instructions stored in SRAM in the SCU itself.

A low-latency passthrough mode is provided for global operations. Because of the latencies associated with the HSSL, the most efficient method to perform global sums is “shift-and-add”, using a store-and-forward capability built into the SCU.

The expected latency of a send (receive) operation is 120 ns (300 ns). This is one to two orders of magnitude lower than the latencies associated with Myrinet. Since a write instruction from the 440 can initiate many independent transfers on any subset of the 24 send or receive channels, the latencies associated with multiple transfers can be overlapped to some degree. The total off-chip bandwidth using all 24 HSSL ports is 12 Gbit/s. In a 4-dimensional physics application only 16 of the 24 HSSL ports will be used, resulting in a total bandwidth of 8 Gbit/s. This provides a good match for the communications requirements of typical applications. As a worst-case example, we consider a $`2^4`$ local lattice and the application of a staggered fermion conjugate gradient. In this case we obtain

| FPU alone | 75% efficiency |
| --- | --- |
| FPU + n.n. communications | 48% efficiency |
| FPU + n.n. + global sums | 38% efficiency |

The 75% efficiency in the first line is a reasonable assumption based on the performance figures quoted in Table 1 below.

### 3.4 Prefetching eDRAM Controller (PEC)

The PEC (shown schematically in Fig. 3) is designed to provide the 440 with high-bandwidth access to the eDRAM. It interfaces to the 440 data read and data write busses via a fast version of the PLB that runs at the CPU frequency and that we call processor direct bus (PDB). The PEC also contains a PLB slave interface to allow for read and write operations from/to any master on the PLB as well as a DMA engine to transfer data between eDRAM and the external DDR memory.

The maximum PDB bandwidth is 8 GBytes/s for read and 8 GBytes/s for write. The latency of the PDB itself is 1-2 CPU cycles. This very low latency eliminates the need for pipelining. The access to eDRAM (which is memory-mapped) proceeds at 8 GBytes/s. ECC is built into the PEC, with 1-bit error correct and 2-bit error detect functionality. The PEC also refreshes the eDRAM.

The read data prefetch from eDRAM occurs in two 1024-bit lines. Three read ports (PDB, PLB slave, DMA) arbitrate for the common eDRAM. The coherency between PDB, PLB slave, and DMA is maintained within the PEC. Each read port has four 1024-bit registers that are paired in two sets to allow for ping-ponging between different memory locations. There are also two 1024-bit write buffer registers each for the PDB/PLB slave/DMA write interfaces.

## 4 SOFTWARE

In this section the various aspects of the QCDOC software environment are described. Efforts are being made to make this environment user-friendly and to allow a wide range of application codes to run on the QCDOC platform without the need for extensive porting.

Boot kernel. The boot kernel will be loaded from the host work station via Ethernet using the reliable UDP protocol. The Ethernet-JTAG interface in the ASIC converts the packets to JTAG commands which in turn load the boot code into the instruction cache and start execution (see also Sec. 3.2). The boot kernel then assigns an IP address to the ASIC (based on the MAC address provided by the location pins) and loads the Ethernet driver.

Run-time kernel. The run-time kernel communicates with the host using the Ethernet driver loaded by the boot kernel. It performs complete node diagnostics and provides read, write, and execute capability. It is based on the UNIX RPC protocol.

RISCWatch debugger. This is an IBM-provided PowerPC development/debugging environment for a set of selected nodes. It provides single-stepping capability with an advanced, user-friendly graphical user interface showing a cycle-by-cycle display of the internal state of the PowerPC processor as OS or user code executes.

Extensive diagnostic programs. These will be provided to permit extensive, real-time hardware monitoring during program execution as well as targeted diagnostic tests for hardware maintenance.

Host Operating System. The host machine will be an SMP workstation with a suitable number of Gbit Ethernet ports, using a threaded transport protocol to efficiently manage tens of thousands of nodes. The host OS will manage both the boot and the run-time kernel operation. It will provide both a qcsh and a socket-based interface to the QCDOC operation, thus allowing for both csh-scripted and Perl-based machine control.

Disk software. We are planning a parallel disk system supporting a general distribution of disks. A simplified UNIX-like hierarchical file system will be used. A utility will be provided to store parallel files as a single, flattened, position-indexed file for archiving and use on other systems.

Compiler, linker, loader. The compiler and linker are GNU-based tools supporting C and C++ as well as Book-E compliant, PowerPC assembly language.

Lattice QCD application code. A code structure is being developed that will both provide the cross-platform capabilities needed by the UKQCD Collaboration and follow the QCD-API standard being developed by the U.S. SciDAC-supported lattice QCD software initiative.

## 5 STATUS AND PERFORMANCE
FIGURES

The functional design of the QCDOC ASIC is nearly finished. Single-node physics code is running on the VHDL simulator, and from these simulations we obtained the performance figures shown in Table 1. These all refer to double-precision calculations.

We have set up a test environment in which two nodes can successfully communicate in the VHDL simulator over their Ethernet interfaces and/or the HSSLs. The operating system is being developed on a system with six PowerPC 405GP boards, one PowerPC 750 board, five Ethernet-JTAG boards, and a RISCWatch probe. This hardware configuration is shown in Figure 4.

The current design has been successfully synthesized, and an analysis net list has been provided to IBM for floorplanning purposes. Anticipating a few iterations to satisfy timing constraints the complete design should be transferred to IBM by the end of 2001.

## 6 CONCLUSIONS

The QCDOC architecture combines state-of-the-art system-on-a-chip technology with custom-designed logic specially optimized for lattice QCD calculations to provide computing power in the tens of TFlops range at a price/performance ratio of 1 US-$ per sustained MFlops. Optimized QCD code is expected to run with an efficiency of roughly 50%. The communications latencies of the QCDOC design are one to two orders of magnitude smaller than those of commercial network solutions used in PC clusters, thus allowing for a large problem to be distributed onto many nodes. This fact, along with significantly lower power consumption, higher reliability, more compact packaging and a better price/performance ratio, motivates the design efforts that are required to construct such a supercomputer.

The first QCDOC machines with a combined peak speed of about 7 TFlops are projected to be finished in 2002, and by the end of 2003 the combined peak speed of QCDOC machines world-wide is targeted at 40 TFlops (10 TFlops of which will be installed in Edinburgh, 10 TFlops are planned for the RBRC and 20 TFlops for the U.S. lattice community).

## ACKNOWLEDGMENTS

This research was supported in part by the U.S. Department of Energy, the Institute of Physical and Chemical Research (RIKEN) of Japan, and the U.K. Particle Physics and Astronomy Research Council.