text:
# Surface relaxations, current enhancements, and absolute distances in high resolution STM.

\[

## Abstract

We have performed the most realistic simulation to date of the operation of a scanning tunneling microscope. Probe-sample distances from beyond tunneling to actual surface contact are covered. We simultaneously calculate forces, atomic displacements, and tunneling currents, allowing quantitative comparison with experimental values. A distance regime below which the probe becomes unstable is identified. It is shown that the real distance differs substantially from previous estimates because of large atomic displacements on the surface and at the probe-tip.

68.37.Ef, 71.15.Mb, 73.40.Gk

\]

It is hard to avoid the scanning tunneling microscope (STM) in the imaging, fabrication, or characterisation of nanoscale structures today . It is certainly the most ubiquitous tool in surface science. But while the formation of atomic wires from metal leads and their final breaking is well researched , the physics of the reverse process - an STM tip approaching the surface at very close range - is far less well documented. This is so even though the phenomena that occur - the forces, atomic relaxations, the sudden changes in image at a critical point - are the everyday experience of every STM experimenter. Even more importantly, several routes are currently being explored to fabricate and operate electronic devices on the nanoscale . In chemical methods of nanofabrication  a high precision of atomic assembly is attainable, even though the exact position of every atom is not controlled. By contrast, to initiate such a process, or to craft structures by moving each element into place, the exact position of the STM tip and the target atoms must be known. Furthermore, nanoscale research requires a high degree of knowledge about the electronic, chemical, and transport properties of potential constituents - information often obtained by STM. This information may be misleading if the tip-sample distance is not correctly estimated. Therefore, as soon as the work moves from qualitative to quantitative, distance estimates are of critical importance. The question of distances has previously been analysed to some extent for STM , and for the related scanning force microscope (SFM) . What is missing, so far, is a picture of the interplay between forces, relaxations, and currents.

In this paper we analyze the interactions between tip and sample and their effect on distance estimates and STM images from first principles. We demonstrate that our simulation technique allows us to keep track of all key quantities in this situation. We choose the Au(111) surface since it is widely used in nanotechnology and was studied in recent extensive experiments .

The computational method is based on first principles density functional theory. As is standard practice in calculations of this type only a small part of the periodic surface is used to determine the physical properties of the whole system. Moreover we employ periodic boundary conditions both parallel and perpendicular to the surface. The sample surface was mimicked by a five layer Au(111) film with nine atoms per layer. The high number of atoms per layer was necessary for a complete decoupling of the interactions between adjacent tip pyramids. The central layer of the Au(111) film was our plane of reference. The tip apex was mimicked by a tungsten tetrahedron mounted on the reverse of the Au(111) surface slab; the tip-sample separation was altered by varying the size of our unit cell normal to the surface. Tungsten is the most commonly used tip material in experiments; the low reactivity and the similar size of tungsten and gold atoms makes the support of the tip apex on gold realistic. We also checked its plausibility by a careful analysis of the forces during our simulations. In this case the force gradients are mainly due to the position of an atom, but not its chemical nature. In Fig. 1 we display the forces on the apex atom and the subsurface atoms of the tip, it can be seen that the forces on the apex atom are about one order of magnitude higher.

Even for the pseudopotential VASP (Vienna ab initio Simulation Program ) code the simulations were rather demanding and at the limit of realistic time scales. This fact ruled out a check of the results including additional gold or tungsten layers. For the same reason we did not include long-range forces. But as established by Perez et al. in their analysis of atomic force microscopy , they only cause significant relaxations in the very low distance regime (below 5 Å), and in this range short-range forces are substantially higher, as we find in our simulations. The same point was made in an earlier publication by Dürig et al. . We also neglect bias-induced surface relaxations in our treatment. The technical details of the calculation are the following: The cutoff of the ultrasoft pseudopotentials  was 230 eV, exchange and correlation effects where computed with generalized gradiant correction and using the formulation of Perdew et al. . The Brillouin zone was sampled with (4 $`\times `$ 4 $`\times `$ 1) k-points of a Monkhorst-Pack grid , forces were converged to values of less than 0.01 eV/Å.

The reading of a piezoscale in the experiments is given by the position of the plane of reference in the supercell setup of our system. From an initial tip-sample separation of 10 Å the distance was gradually reduced, initially in steps of 0.5 Å. At every step the system was fully relaxed. The onset of significant mechanical interactions between tip and sample occurs at about 5.0 Å. To determine the precise effect of these interactions in the low distance range we reduced the distance in steps of 0.1 Å from 5.0 to 4.0 Å for one tip position (STM tip on top of a surface atom), and from 4.5 to 3.5 Å for the second tip position (STM tip between the surface atoms). We computed the surface wavefunctions for the fully relaxed sample and tip systems at every distance and for two tip positions. In this part of the calculation the initial system was separated into two subsystem, the sample and the tip surface, and the surface wavefunctions computed for each subsystem separately. These wavefunctions were used as the input for our simulation of the tunnel currents . The procedure leads to separate constant current contours of the surface, depending on the position of the tip. Interpolating between the two simulations we obtain the change of corrugation due to system relaxations from first principles.

It has been shown by a number of authors  that substantial relaxations are confined to the low distance regime. Perez et al. obtained their result with an ab initio method, but they calculated the changes on a semiconductor surface and did not relate them to measured currents. Olesen et al., and Clarke et al. analysed the events on metal surfaces, but employed empirical potentials; for this reason they too could not simultaneously evaluate the current. Their results show, not surprisingly, that the onset of forces and relaxations begins when the interatomic distances are close to the typical bondlength of the atoms. However, from experiments we know that tip-sample forces on Au(111) are appreciable for distances of at least 5 Å . Di Ventra and Pantelides employed density functional theory to determine the relaxation of surface atoms for Al(110) surfaces, mimicking the tip by a single Al atom . Within the range of distances from the point contact to the tunneling regime (1.3 to 3.3 Å) they predict no corrugation due to atomic relaxations. The present calculation is in a sense complementary to their work. It is shown here that experimental conditions in STM scans do not allow such a close approach, because the system becomes unstable already at larger distances (4.6 Å). Our calculations establish an early onset of relaxations, which depend strongly on the position of the tip. For the on-top position of the STM tip we observe an onset of forces and relaxations at about 4.7 Å, about 1.0 Å higher than suggested by earlier simulations . The surface atom reaches its furthest elongation from the surface at a tip-sample separation of 4.5 Å; it is then 1.3 Å above its position on the isolated surface. This point marks the beginning of atomic transfer from the surface to the tip, the origin of the hysteresis observed in single approach/retraction cycles . Neighboring atoms on the surface are equally relaxed, their maximum displacement is 0.1 Å. The second nearest neighbors remain close to their original position (vertical displacement less than 0.04 Å). Even though the boundary conditions at our two dimensional repeat unit constrain the relaxations of individual atoms, the very small relaxation of second nearest neighbors indicates a rapid decay of the forces in lateral direction. Therefore we do not expect substantial changes in a calculation performed with a larger unit cell. If the tip is in the threefold hollow, the onset of forces and relaxations is substantially retarded and occurs at 4.3 Å. The maximum outward relaxation of the three gold atoms is somewhat lower and about 1.0 Å (reached when the tip-sample separation is 4.1 Å).

The tungsten tetrahedron of the STM tip relaxes about the same amount in both positions: the apex atom approaches the sample by about 0.5 Å at the point of closest approach. The change is due to an increase of the interlayer distance of the tungsten atoms by about 0.1 to 0.2 Å, and a buckling of the layers underneath, which extends about two layers into the lead. The forces on the tungsten tip are shown in Fig. 1. The values for tip-sample relaxations are given in Fig. 2. Owing to the strong interactions between tip and sample at close range, the core-core distance between tip and sample is decreased. The net effect of relaxations is a decrease of the distance between tip and sample. A deviation of about 2 Å between the distance inferred from simulations, and the distance inferred from measurements of the decay, has actually been observed . In the critical range of 4.6 to 4.1 Å the actual distance changes by about 2 Å. The real distances are shown in Fig. 3 (a). This points to an effect frequently observed in experiments with atomic resolution: it is difficult to obtain atomically-resolved images, because sizable corrugations in most cases require a distance below 5 Å, but it is very easy to destroy an STM tip by small increases of the current. We now understand that the actual range of imaging is less than 1.0 Å, or less than an order of magnitude in the current scale. Stability is destroyed, when the apex atom of the tip and an atom of the sample form a chemical bond. To check, whether a different metal surface would lead to different results and thus a different interval in actual measurements we computed the forces between a tungsten atom and a Cu, Ag, and Au atom in a dimer configuration. The result of our calculation is shown in the inset of Fig. 1. It can be seen that the forces are equal for all three noble metals in the range above 4 Å. Short of performing identical simulations for all three surfaces this seems the best confirmation that our result is of general validity.

Currents were calculated in a perturbation approach  and based on the fully relaxed sample and tip systems. The current was calculated for a sample bias of -0.1 V. The curves in Fig. 3 (b) show the current for the hollow position. We obtained the change due to surface and tip relaxations by calculating the full range of currents in a frozen atomic configuration. The logarithmic graphs are different for every distance, the z-labels give the tip-sample separation before relaxation. The true increase of the current is determined by taking only a single point of each graph. In practice we interpolated between the computed points to obtain the smooth black curve shown in the figure, giving an enhanced increase in the current. For an experimenter who estimates the potential barrier between tip and sample from the exponential decay of the current, the apparent barrier height increases in this range . The current calculation is based on the Baardeen approach , where sample and tip are considered decoupled. In a previous paper we studied the effects of a tip potential on the surface electronic structure , and it was established that these effects are limited to very low distances (below 4 Å). Here we find that the onset of relaxations begins about 1 Å before this point. The main effect, not included in perturbation theory, is therefore due to relaxations. But this indicates that a more refined treatment of STM currents e.g. by a scattering approach is no substantial improvement compared to the Bardeen approximation, if it does not simultaneously treat relaxations. And this, in turn, is still computationally untractable because it requires a fully ab initio treatment for every single position of the tip.

All the results presented here are derived from ab initio treatments of currents and relaxations for two lateral positions of the tip: the on-top position and the threefold hollow. The tunnel currents in our simulations were assumed to be 5.1 nA in both cases, but owing to the site-dependent relaxation of the systems the two constant-current contours differ. The piezoscale would give a distance of 4.7 Å on top of a gold atom and 4.5 Å in the fcc hollow site. These distances are at the lowest limit of stability for the tip-sample system. The difference of 0.2 Å describes the experimental corrugation of the Au(111) surface, when measured with a tungsten terminated tip . The simulated scans with the tips and surfaces frozen in the two relaxed configurations are shown in Fig. 4 (a), and (b). We stress that these images are idealizations, because they were computed with the wavefunctions of frozen atomic positions of tip and sample. In Fig. 4 (c) we show two simulated linescans across a single atom. The two separate curves describe the linescan for the two setups used (on-top, and hollow). It can be seen that the relaxation leads in effect to an enhancement of the corrugation by about 0.13 Å. Only with this relaxation do we obtain agreement between experiments and simulations .

The most remarkable result of our calculations, and the one that experimenters may find most important, is the very low range of stability for atomically resolved STM scans on metal surfaces. This range is only about 0.5 Å. Below this point, at about 4.5 Å from the surface, the distance between tip and sample changes drastically by 2 Å (within 0.5 Å of measured displacement), as the two separate systems jump into contact. In this range the current may even decrease as the tip approaches, if one atom of the tip or sample system becomes separated.

The work was supported by the British Council and the National Research Council. Computing facilities at the UCL HiPerSPACE center were funded by the Higher Education Funding Council for England. AJF was also supported by an Advanced Fellowship from the Engineering and Physical Sciences Research Council. PG was supported by the National Science and Engineering Research Council of Canada and Le Fonds pour la Formation de Chercheurs et l’Aide a la Recherche de la Province de Quebec.


format: text

text:
# Full Frequency Back-Action Spectrum of a Single Electron Transistor during Qubit read-out

## Abstract

We calculate the spectral density of voltage fluctuations in a Single Electron Transistor (SET), biased to operate in a transport mode where tunneling events are correlated due to Coulomb interaction. The whole spectrum from low frequency shot noise to quantum noise at frequencies comparable to the SET charging energy $`(E_C/\mathrm{})`$ is considered. We discuss the back-action during read-out of a charge qubit and conclude that single-shot read-out is possible using the Radio-Frequency SET.

The Single Electron Transistor (SET) has recently been suggested as a read-out device for solid-state charge qubits Aassime et al. (2001); Kane et al. (2000); Lea et al. (2000); Averin (2000); Makhlin et al. (2000). Presenting new state of the art figures for sensitivity in charge measurements, Aassime et al. Aassime et al. (2001) made probable that the Radio-Frequency SET Schoelkopf et al. (1998) (RF-SET) can be used for single shot read-out of the Single Cooper-pair Box (SCB) qubit Bouchiat et al. (1998); Nakamura et al. (1999). This is possible if the measurement time $`t_{ms}`$ needed to resolve the two states of the qubit is much shorter than the time $`t_{mix}`$ required to destroy the initial qubit state because of back-action due to voltage fluctuations on the SET.

The mixing rate $`1/t_{mix}`$ is proportional to the noise at the frequency $`\omega `$ corresponding to the energy splitting between the qubit states $`\mathrm{\Delta }E`$. Moreover, in order to have well defined charge states to achieve single shot read-out, $`\mathrm{\Delta }E`$ should be comparable to the qubit charging energy $`E_{qb}`$. The ratio $`E_{qb}/E_C`$ may be varied between approximately 1/10 and 10 in order to optimize the performance of the system.

The voltage fluctuations on the SET island have previously been studied in the experimentally accessible low frequency limit ($`\mathrm{}\omega E_C`$). The noise is here fully understood both in the classical transport regime of sequential tunneling, described by a master equation with frequency independent tunneling rates Korotkov et al. (1992); Hanke et al. (1993); Korotkov (1994), as well as in the Coulomb blockade cotunneling regime Averin (2000). In Ref. Aassime et al. (2001) the back-action of the SET was estimated using simple interpolation between the low frequency shot noise limit and the high frequency limit of Nyquist noise from two independent tunnel junctions.

In this paper we evaluate the finite frequency noise on a fully quantum mechanical basis, taking into account the frequency dependence of the tunneling rates in the presence of fluctuations, separating processes that absorb or emit energy. We present a simple expression for the voltage noise of a SET biased in the transport mode, valid in the whole frequency range from low frequency classical shot noise to high frequency quantum noise. By comparing $`t_{mix}`$ with $`t_{ms}`$ for reading out a SCB qubit with the RF-SET of Ref. Aassime et al. (2001), we conclude that $`t_{ms}t_{mix}`$, i.e. single-shot read-out should indeed be possible.

Consider a small metallic SET island coupled via low transparency tunnel barriers to two external leads, and coupled capacitively to the SCB, which is controlled by a gate voltage (see Fig.1). First we calculate the voltage noise on the SET island taking the voltage $`V_g`$ on the SCB to be constant. Then the back-action of these fluctuations on the state of the SCB is estimated. This approach is appropriate in the considered limit of weak SET-qubit coupling ($`\kappa C_c/C_{qb},\kappa 1;C_L,C_RC_{qb}`$).

We follow the outline of Ref. Schoeller and Schön (1994) and model the SET by the Hamiltonian

$$H=H_L+H_R+H_I+V+H_T=H_0+H_T,$$
(1)
where

$$H_r=\underset{kn}{}ϵ_{kn}^ra_{krn}^{}a_{krn},H_I=\underset{ln}{}ϵ_{ln}c_{ln}^{}c_{ln}$$
(2)
describe noninteracting electrons in the left/right lead ($`H_r,r\{L,R\}`$) and on the island ($`H_I`$). The quantum numbers $`n`$ denote transverse channels including spin, and $`k,l`$ denote momenta. The Coulomb interaction on the island is described by

$$V(\widehat{N})=E_C(\widehat{N}n_x)^2,$$
(3)
where $`\widehat{N}`$ denotes the excess number operator, $`E_C=e^2/2C`$ the charging energy ($`C=C_L+C_R+C_c`$), $`n_x`$ the fractional number of electrons induced by the external voltages ($`n_x`$ is the fractional part of $`(C_LV_L+C_RV_R+C_cV_g)/e`$) and $`e`$ the electron charge. The tunneling term is

$$H_T=\underset{r=L,R}{}\underset{kln}{}(T_{kl}^{rn}a_{krn}^{}c_{ln}e^{i\widehat{\mathrm{\Phi }}}+T_{kl}^{rn}c_{ln}^{}a_{krn}e^{i\widehat{\mathrm{\Phi }}}),$$
(4)
where the operator $`e^{\pm i\widehat{\mathrm{\Phi }}}`$ changes the excess particle number on the island by $`\pm 1`$ and $`T_{kl}^{rn}`$ are the tunneling matrix elements. $`\widehat{\mathrm{\Phi }}`$ is the canonical conjugate to $`\widehat{N},([\widehat{\mathrm{\Phi }},\widehat{N}]=i)`$. In this case of a metallic island, containing a large number of electrons, the charge degree of freedom $`N=0,\pm 1,\mathrm{}`$ is to a very good approximation independent of the electron degrees of freedom $`l,n`$.

The Cooper-pair box is a small superconducting island coupled to a superconducting reservoir, via a tunnel junction with Josephson energy $`E_J`$ (see Fig. 1). The box is also characterized by its charging energy $`E_{qb}=e^2/2C_{qb}`$, where $`C_{qb}`$ is the total capacitance of the box. For $`E_JE_{qb}`$ two nearby charge states, corresponding to zero and one excess Cooper pair on the island, can be used as the $`|0`$ and $`|1`$ states of the qubit. The energy splitting $`\mathrm{\Delta }E`$ between these states is controlled by the gate voltage $`V_g`$.

For a given charge measurement sensitivity $`\delta q`$ and measurement time $`t_{ms}`$, the uncertainty in charge is given by $`\mathrm{\Delta }q=\delta q/\sqrt{t_{ms}}`$. In order to separate the two qubit states we need the two intervals $`0\pm \mathrm{\Delta }q`$ and $`2e\pm \mathrm{\Delta }q`$ not to overlap, which defines the measuring time $`t_{ms}=\delta q^2/(\kappa e)^2`$.

During the measurement, the voltage fluctuations on the SET island will induce transitions between the qubit states. The rate for excitation/relaxation of the qubit is proportional to the SET noise spectral density $`S_V(\omega )`$ at the negative/positive frequency corresponding to the transition $`\mathrm{\Delta }E/\mathrm{}`$. The information of the initial state is destroyed at the combined rate Makhlin et al. (2001); Devoret and Schoelkopf (2000)

$$\mathrm{\Gamma }_1=\frac{1}{t_{mix}}=\frac{e^2}{\mathrm{}^2}\kappa ^2\frac{E_J^2}{\mathrm{\Delta }E^2}\left[S_V(\frac{\mathrm{\Delta }E}{\mathrm{}})+S_V(\frac{\mathrm{\Delta }E}{\mathrm{}})\right],$$
(5)
here given in the relevant limit $`\mathrm{\Delta }EE_J`$ used below.

The spectral density of voltage fluctuations on the SET island is described by the Fourier transform of the voltage-voltage correlation function

$$S_V(\omega )=\frac{e^2}{C^2}_{\mathrm{}}^{\mathrm{}}𝑑\tau e^{i\omega \tau }Tr\{\rho _{st}(t_0)\widehat{N}(\tau )\widehat{N}(0)\}.$$
(6)
Here $`\rho _{st}(t_0)`$ is the density matrix of the system in steady-state, which is assumed to have been reached at some time $`t_0`$ before the fluctuation occurs ($`t_0<\mathrm{min}\{0,\tau \}`$). $`\rho _{st}`$ is the tensor product of the equilibrium (Fermi distributed) density matrix $`\rho _{eq}^e`$ for the electron degrees of freedom in each reservoir ($`L,R,I`$) and a reduced density matrix $`\rho _{st}^c`$, describing the charge degrees of freedom, which we assume to be diagonal Schoeller (1997) with elements $`P_N^{st}`$ denoting the probability of being in charge state $`N`$.

To evaluate $`S_V(\omega )`$ we make a perturbation expansion of the forward and backward time evolution operators in terms of the tunneling term $`H_T`$, as shown in Fig. 2a. The reservoir degrees of freedom are traced out using Wick’s theorem. In the diagrammatic language of reference Schoeller and Schön (1994) (Fig. 2) we have forward and backward propagators of the Keldysh contour (horizontal lines), with internal charge transfer vertices $`e^{\pm \widehat{\mathrm{\Phi }}}`$ (small dots) connected by reservoir lines. The $`\widehat{N}(0)`$ and $`\widehat{N}(\tau )`$ operators form external vertices (big dots). Figure 2a shows the diagrammatic expression for $`S_V(\omega )`$ in our approximation, neglecting processes involving other charge states than $`N\{0,1\}`$. Here $`\mathrm{\Pi }_{NN^{}}(\omega )`$ is the frequency dependent rate for transitions between the charge states $`N`$ and $`N^{}`$.

One may graphically write down a Dyson type of equation for $`\mathrm{\Pi }_{NN^{}}(\omega )`$ (Fig. 2b). In matrix notation this reads

$$\widehat{\mathrm{\Pi }}(\omega )=\widehat{1}\frac{i}{\omega }+\widehat{\mathrm{\Pi }}(\omega )\frac{i}{\omega }\frac{i\widehat{\mathrm{\Sigma }}(\omega )}{\omega }\widehat{\mathrm{\Pi }}(\omega )=\frac{i}{\omega }\left[\widehat{1}\frac{i}{\omega }\widehat{\mathrm{\Sigma }}(\omega )\right]^1,$$
(7)
where $`\mathrm{\Sigma }_{NN^{}}(\omega )`$ represents the sum of irreducible diagrams, i.e. containing no free propagators, for transitions between $`N`$ and $`N^{}`$.

Since the SET is biased in transport mode, we consider only the dominating lowest order diagrams in $`\mathrm{\Sigma }_{NN^{}}(\omega )`$, corresponding to single tunneling events. In this approximation we get Schoeller and Schön (1994)

$$\mathrm{\Sigma }_{N,N\pm 1}(\omega )=\gamma _N^\pm (\omega )+\gamma _N^\pm (\omega ),$$
(8)
where $`\gamma _N^\pm (\omega )`$ denotes the rate, in second order pertubation theory, to go from charge state $`N`$ to $`N\pm 1`$ while the SET is absorbing ($`\omega `$) or emitting ($`\omega `$) a quantum of energy $`|\mathrm{}\omega |`$. The well-known expressions for these rates are

$$\gamma _N^+(\omega )=\frac{\pi }{\mathrm{}}\underset{r}{}(\mathrm{}\omega +\mathrm{\Delta }_N^r)\alpha _0^rn(\mathrm{}\omega +\mathrm{\Delta }_N^r),$$
(9)
$$\gamma _N^{}(\omega )=\frac{\pi }{\mathrm{}}\underset{r}{}(\mathrm{}\omega \mathrm{\Delta }_{N1}^r)\alpha _0^rn(\mathrm{}\omega \mathrm{\Delta }_{N1}^r)$$
(10)
where the dimensionless conductivity

$$\alpha _0^r=\underset{n}{}\left|T^{rn}\right|^2\rho _r^n\rho _I^n=\frac{R_K}{4\pi ^2R_T^r}$$
(11)
is the ratio of the quantum resistance $`R_K=h/e^2`$ and the resistance of barrier $`r\{L,R\}`$, and $`\beta =1/k_BT`$. $`n(E)=1/\{1\mathrm{exp}[\beta E]\}`$ is the Bose-Einstein distribution of electron-hole excitations and comes from the convolution of the Fermi distributions for filled initial states and empty final states. $`\mathrm{\Delta }_N^r=V(N)+\mu _rV(N+1)\mu _I`$ is the energy gained by an electron tunneling from the chemical potential of lead $`r`$ ($`\mu _r`$) to the chemical potential of the island ($`\mu _I`$), thus changing the charge state from $`N`$ to $`N+1`$ (see Figs. 3a-c). We also assume that the relaxation processes on the island are fast on the timescales we are looking at and neglect the energy dependence of both tunneling matrix elements $`T_{kl}^{rn}T^{rn}`$ and reservoir densities of states $`\rho _r^n(ϵ),\rho _I^n(ϵ)`$. To avoid renormalization effects Schoeller and Schön (1994) we choose not too small bias, i.e. $`\mathrm{max}_r\alpha _0^rln(E_c/|\mathrm{\Delta }_0^r|)1`$.

We also choose the bias not too high, so that the higher charge states are inaccessible at low frequency, i.e. $`\gamma _1^+(0)=\gamma _0^{}(0)=0`$ giving $`P_N^{st}=0`$ for $`N\{0,1\}`$. At low frequencies we may thus restrict ourselves to the two lowest charge states $`N\{0,1\}`$. The matrix inverse in Eq. (7) may then easily be performed analytically. For high enough frequencies also the charge states $`N\{1,2\}`$ are accessible. In this regime we may expand Eq. (7) to first order in $`\widehat{\mathrm{\Sigma }}/\omega `$. The results in both these regimes combine to the following expression for the noise valid, within our approximations, for all frequencies:

$$S_V(\omega )=\frac{2e^2}{C^2}\frac{P_0^{st}\left[\gamma _0^+(\omega )+\gamma _0^{}(\omega )\right]+P_1^{st}\left[\gamma _1^{}(\omega )+\gamma _1^+(\omega )\right]}{\omega ^2+\left[\gamma _0^+(\omega )+\gamma _1^{}(\omega )+\gamma _0^+(\omega )+\gamma _1^{}(\omega )\right]^2},$$
(12)
where $`P_0^{st}=\gamma _1^{}(0)/(\gamma _0^+(0)+\gamma _1^{}(0))`$ and $`P_1^{st}=1P_0^{st}`$.

Note the simple structure of this expression: a sum over the probabilities of being in the states $`N=0`$ and $`N=1`$ times the rates for possible transitions from these states, normalized by a denominator containing their finite lifetime due to the bias.

At zero frequency Eq. (12) coincides with the classical shot noise result Korotkov et al. (1992); Korotkov (1994); Devoret and Schoelkopf (2000), which is recovered by approximating all rates with their zero frequency value. In the high-frequency limit, where the biasing energy is negligible, one recovers the high-frequency Nyquist noise from two resistances connected in parallell to ground.

It is informative to analyze Eq. (12) at zero temperature where the tunneling rates are linear functions of the possible energy gain (since $`n(E)=\mathrm{\Theta }(E)`$, where $`\mathrm{\Theta }(x)`$ is the unit step function). To be specific we assume that $`\mu _L>\mu _I>\mu _R`$ and $`\mathrm{\Delta }_0>0`$ so that electrons tunnel from the left lead to the right and $`N=0`$ is the lowest charge state. We also assume that the left bias is larger than the right one, $`|\mathrm{\Delta }_0^L|>|\mathrm{\Delta }_0^R|`$ (see Fig. 3a).

For large negative frequencies $`\mathrm{}\omega <|\mathrm{\Delta }_0^L|`$ the fluctuations try to extract more energy than the bias can provide, thus there is no contribution to the noise from this frequency range (since we neglect higher order cotunneling effects). Between $`|\mathrm{\Delta }_0^L|<\mathrm{}\omega <|\mathrm{\Delta }_0^R|`$ the SET may emit energy $`|\mathrm{}\omega |`$ while tunneling through the left barrier, while for $`|\mathrm{\Delta }_0^R|<\mathrm{}\omega <0`$ tunneling through both barriers contribute. For positive frequencies the SET may absorb energy $`\mathrm{}\omega `$ while tunneling through either barrier. For frequencies $`\mathrm{}\omega >|\mathrm{\Delta }_0^{L/R}|`$ there is also a contribution from electrons tunneling backwards, against the bias, at the (L/R) junction. In the region $`|\mathrm{\Delta }_0^{L/R}|\mathrm{}\omega <\{|\mathrm{\Delta }_1^{L/R}|,|\mathrm{\Delta }_1^{L/R}|\}`$ the higher charge states $`N\{1,2\}`$ are still inaccessible and the noise approaches exactly one half of the Nyquist noise. This may be explained by the strong correlation of fluctuations at the two barriers induced by the Coulomb interaction. At high frequencies $`\mathrm{}\omega >\{|\mathrm{\Delta }_1^{L/R}|,|\mathrm{\Delta }_1^{L/R}|\}`$ also processes at barrier $`L/R`$ involving the charge states $`N\{1,2\}`$ contribute, giving the real high-frequency limit of Nyquist noise from two completely uncorrelated tunnel junctions. Figure 4 illustrates the above discussion of which processes contribute in which frequency region.

In the low frequency regime $`|\mathrm{}\omega |<|\mathrm{\Delta }_0^R|`$, where no backward processes are allowed, the noise spectral density may be written as

$$S_V(\omega )=\frac{e^2}{C^2}\frac{2I/e+2\pi \omega \left[P_0^{st}\alpha _0^L+P_1^{st}\alpha _0^R\right]}{\omega ^2+4\left[\gamma _0^+(0)+\gamma _1^{}(0)\right]^2},$$
(13)
where $`I=2e\gamma _0^+(0)\gamma _1^{}(0)/(\gamma _0^+(0)+\gamma _1^{}(0))`$ is the DC current through the SET. This is the expression for classical shot noise plus a linear $`\omega `$-term in the numerator, indicating that the probability of processes where the SET absorbs/emits energy ($`\pm \omega `$) is raised/lowered compared to the classical expression (see Figs. 3c,d and the right inset of Fig. 5). The total frequency dependent decay rate of the states, found in the denominator, is constant in this regime since the linear increase/decrease in rates for processes absorbing/emitting energy $`|\mathrm{}\omega |`$ cancel exactly, at zero temperature.

The sum of positive and negative frequency SET noise, proportional to the back-action on the qubit, is plotted in Fig. 5, using the parameters of the RF-SET in Ref. Aassime et al. (2001) with sensitivity $`\delta q=6.3\mu e/\sqrt{Hz}`$. One may define the signal-to-noise ratio ($`SNR`$) of a single shot measurement as $`SNR=\sqrt{t_{mix}/t_{ms}}`$, which is plotted as a function of qubit level splitting $`\mathrm{\Delta }E`$ in the left inset of Fig. 5. Using realistic qubit parameters for an SCB made with the same aluminum technique as the RF-SET, $`E_{qb}/k_B=1K`$, $`\mathrm{\Delta }E=2.4E_{qb}`$, $`E_J/E_{qb}=0.1`$, assuming read-out at $`\mathrm{\Delta }E/E_{qb}=2.4`$ we find that $`SNR=\sqrt{t_{mix}/t_{ms}}5`$, which indicates that single-shot measurement is possible. If one increases the qubit charging energy, e.g. by using a niobium qubit, the $`SNR`$ increases even further.

In conclusion, we have obtained a simple expression for the finite frequency voltage fluctuations of a SET island biased in the transport mode. Using this information about back-action we can conclude that single shot read-out of a SCB qubit using an RF-SET is possible.

###### Acknowledgements.

We would like to acknowledge fruitful discussions with Per Delsing, Yuriy Makhlin, Vitaly Shumeiko, Tomas Löfwander and Abdel Aassime. We are especially grateful to Michel Devoret and Daniel Esteve for the formulation of the problem. This work was supported by NFR and the European Union under the SQUBIT project.


text:
# Fluctuation formula for nonreversible dynamics in the thermostated Lorentz gas

## Abstract

We investigate numerically the validity of the Gallavotti-Cohen fluctuation formula in the two and three dimensional periodic Lorentz gas subjected to constant electric and magnetic fields and thermostated by the Gaussian isokinetic thermostat. The magnetic field breaks the time reversal symmetry, and by choosing its orientation with respect to the lattice one can have either a generalized reversing symmetry or no reversibility at all. Our results indicate that the scaling property described by the fluctuation formula may be approximately valid for large fluctuations even in the absence of reversibility.

The Lorentz gas (LG) thermostated by a Gaussian isokinetic (GIK) thermostat is one of the most popular models in the study of the relationship between transport properties and chaotic behaviour in nonlinear dynamical systems. Since the microscopic dynamics of the LG is chaotic, and on a sufficiently long time scales it possesses a well defined macroscopic transport coefficient, it can be used to study the connection of microscopic chaos and macroscopic nonequilibrium behaviour.

The so-called fluctuation formula (FF) has first been observed numerically in a system of thermostated fluid particles undergoing shear flow shearflow . In that model, trajectory segments violating the second law of thermodynamics were found with probabilities exponentially smaller than those of trajectory segments associated with normal thermodynamical behaviour. More precisely, let $`\xi _\tau (t)`$ denote the entropy production rate $`\xi `$ averaged over a time interval of length $`\tau `$ centered around time $`t`$: $`\xi _\tau (t)=\frac{1}{\tau }\underset{\tau /2}{\overset{\tau /2}{}}\xi (t+t^{})dt^{}`$, and let us consider it as a probabilistic variable. Then its statistical properties in a steady state can be characterized by a probability density $`\mathrm{\Xi }_\tau (x)`$. The fluctuation formula states orig  that this probability density has the following property:

$$\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }\mathrm{ln}\frac{\mathrm{\Xi }_\tau (x)}{\mathrm{\Xi }_\tau (x)}=x.$$
(1)
One of the interesting features of the FF is that it seems to be valid in systems far from equilibrium, not just for vanishing external fields.

After discovering the formula numerically, analytical results were obtained about its validity in deterministic systems like transitive Anosov systems dynens  and special reversible maps baker ; dorfmannbook . In the proofs of these theorems, the time reversibility of the system plays a key role stocha . Nevertheless, proving fluctuation theorems under more general conditions seems to be exceedingly difficult. In this context, even relatively simple systems like the LG, with or without magnetic field, seem to be out of reach for the existing analytical techniques. It is also unclear how the FF should look like in systems with nonreversible dynamics rem-windtree . Consequently, reliable numerical results for such models may provide valuable hints in the search for more sophisticated theoretical approaches.

The field driven Lorentz gas consists of a charged particle subjected to an electric field moving in the lattice of elastic scatterers. For the sake of simplicity, we take a square or cubic lattice in our study, depending on the dimensionality of the system. Due to the applied electric field, one must use a thermostating mechanism to achieve a steady state in the system. Such a tool is the Gaussian isokinetic thermostat which preserves the kinetic energy of the particle; for a review see e.g. thermostats  and further references therein. We will also apply a constant external magentic field to control the reversibility of the dynamics.

Throughout our work we use dimensionless variables. We choose the units of mass and electric charge to be equal to the mass and electric charge of the particle, so we have $`m=q=1`$ in our model. The unit of distance is taken to be equal to the radius of scatterers ($`R=1`$), and the unit of time is chosen to normalize the magnitude of particle velocity to unity. Let $`𝐪=(q_1,\mathrm{},q_n)`$ denote the position and $`𝐩=(p_1,\mathrm{},p_n)`$ the momentum of the particle in the $`n`$-dimensional space ($`n=2`$ or 3). Due to the normalization, $`\left|𝐩\right|=1`$. The phase space variable of the system is $`𝚪=(𝐪,𝐩)`$; it is transformed abruptly at every elastic collision and evolved smoothly by the differential equation

$`\dot{𝐪}`$ $`=`$ $`𝐩`$
$`\dot{𝐩}`$ $`=`$ $`𝐄+𝐩\times 𝐁\alpha 𝐩`$ (2)
between them. Here $`\alpha `$ is called the thermostat variable, while $`𝐄`$ and $`𝐁`$ are constant vectors playing the roles of the external electric and magnetic fields, respectively. The GIK thermostat corresponds to the choice $`\alpha =\mathrm{𝐄𝐩}`$ in Eq. (Fluctuation formula for nonreversible dynamics in the thermostated Lorentz gas). For $`n=2`$, $`𝐁`$ is thought to be perpendicular to the plane of motion given by the directions of $`𝐄`$ and $`𝐩`$. We note that Eq. (Fluctuation formula for nonreversible dynamics in the thermostated Lorentz gas) is dissipative, but for $`𝐁=0`$ it has also time reversal symmetry.

Dissipation can be measured by the phase space contraction rate $`\sigma `$. It can be computed by taking the divergence of the right-hand side of Eq. (Fluctuation formula for nonreversible dynamics in the thermostated Lorentz gas):

$`\sigma =\text{div}\dot{𝚪}=\left(n1\right)\alpha ,`$ (3)
and it can be shown (see e.g. epr ) that in our case

$`\sigma (t)`$ $`=`$ $`\xi (t).`$ (4)
We note that the validity of this identity does depend on the choosen model and cannot be treated as a general property thermostats ; nosehoover ; detscatt .

The notion of reversibility revers , an extension of time reversal symmetry, can be formulated in terms of the phase space flow $`𝚽^t`$ defined by $`𝚪(t)=𝚽^t𝚪_0`$. We say that the flow is reversible, if there exists a map $`𝐆`$ which is an involution (i.e $`𝐆^2`$ is the identity) and bracketing the flow by $`𝐆`$ reverses the direction of time:

$`𝐆𝚽^t𝐆=𝚽^t.`$ (5)
Time reversal symmetry is a special case of reversibility with a particular choice of the involution: $`𝐆_0(𝐪,𝐩)=(𝐪,𝐩)`$.

In the LG, reversibility depends on the directions of the field vectors relative to each other and the lattice. It can be checked easily that our system is time reversible if $`𝐁=\mathrm{𝟎}`$, and it is not otherwise. In Ref. cpr , we have also shown that the system is still reversible for $`𝐁\mathrm{𝟎}`$ if the plane containing $`𝐄`$ and $`𝐁`$ is a symmetry plane of the lattice. Since the transformation $`𝐆=\mathrm{𝐌𝐆}_0`$ (where $`𝐌`$ is a mirroring of $`𝐪`$ and $`𝐩`$ with respect to the plane containing $`𝐄`$ and $`𝐁`$) satisfies Eq. (5), the smooth flow is always reversible. This means that the reversibility of the full dynamics including the collisions requires that the invariant plane of $`𝐌`$ be a symmetry plane of the lattice cpr . In the two dimensional case this is simplified to the condition that $`𝐄`$ has to be contained by the symmetry plane of the lattice.

The goal of our numerical simulations was to measure $`\mathrm{\Xi }_\tau (x)`$ with a precision which is sufficient to check the validity of the fluctuation formula. Due to Eq. (4), $`\mathrm{\Xi }_\tau (x)`$ could be measured by periodically computing $`\sigma _\tau `$ along a particle trajectory and making a histogram of these data. The disadvantage of this method is that the range of possible $`\sigma _\tau `$ values depends on the strength of the electric field. Instead we may introduce the quantity

$`\pi _\tau (t)={\displaystyle \frac{1}{\tau }}{\displaystyle \underset{\frac{\tau }{2}}{\overset{\frac{\tau }{2}}{}}}𝐧_𝐄𝐩(t+t^{^{}})dt^{^{}},`$ (6)
where $`𝐧_𝐄`$ denotes the unit vector paralell to $`𝐄`$. Since the magnitude of $`𝐩`$ is unity, $`\pi _\tau `$ always satisfies $`\pi _\tau [1,1]`$. By making a histogram of the periodically measured values of $`\pi _\tau `$, one gets an approximation of its probability density $`\mathrm{\Pi }_\tau (x)`$. Since the two probabilistic variables satisfy $`\sigma _\tau =(n1)E\pi _\tau `$, the connection of the probability densities is

$`\mathrm{\Xi }_\tau (x)={\displaystyle \frac{1}{(n1)E}}\mathrm{\Pi }_\tau \left({\displaystyle \frac{x}{(n1)E}}\right).`$ (7)
Then we can rewrite the fluctuation formula as

$`\underset{\tau \mathrm{}}{lim}{\displaystyle \frac{1}{(n1)E}}{\displaystyle \frac{1}{\tau }}\mathrm{ln}{\displaystyle \frac{\mathrm{\Pi }_\tau (x)}{\mathrm{\Pi }_\tau (x)}}=x.`$ (8)
At a first glance, $`\mathrm{\Pi }_\tau (x)`$ behaves similarly in all cases: as $`\tau `$ grows, $`\mathrm{\Pi }_\tau (x)`$ becomes more and more concentrated around its mean value. This typical shape is shown in Fig. 1. It can be noticed that the curve looks like a Gaussian, although it is clear that it must be different due to the finite range of $`x`$ rem-Gauss . In a separate paper next , we will deal with the properties of this distribution in more details.

In order to visualize the fluctuation formula, we introduce the quantity

$`D_\tau (x)={\displaystyle \frac{1}{(n1)E}}{\displaystyle \frac{1}{\tau }}\mathrm{ln}{\displaystyle \frac{\mathrm{\Pi }_\tau (x)}{\mathrm{\Pi }_\tau (x)}}`$ (9)
that must exactly be linear with a slope 1 in the $`\tau \mathrm{}`$ limit if the fluctuation formula is valid. We will investigate for different configurations of the LG how well $`D_\tau (x)`$ approaches this behavior in numerical simulations. Due to the fact that we have a finite number of data points coming from a numerical trajectory of finite length, our conclusions concerning $`D_\tau (x)`$ and thus the fluctuation formula are, of course, limited to an interval $`[\mathrm{\Delta }_\tau ,\mathrm{\Delta }_\tau ]`$ with $`\mathrm{\Delta }_\tau 1`$. In practice, if the extremal $`\pi _\tau `$ values in a series of $`N`$ measured data were $`\pi _{min}`$ and $`\pi _{max}`$, then we identified $`\mathrm{\Delta }_\tau `$ with $`\mathrm{min}(\pi _{min},\pi _{max})`$. It is easy to check that the probability of observing $`\pi _\tau `$ values outside this interval in another series is in the order of $`1/N`$. Fig. 2 shows the $`\tau `$ dependence of $`\mathrm{\Delta }_\tau `$ for different field strengths with $`N`$ fixed. For our simulations, $`N`$ was chosen to be $`10^9`$, which means that if $`D_\tau (x)`$ is found to be linear on $`[\mathrm{\Delta }_\tau ,\mathrm{\Delta }_\tau ]`$ with slope 1, then it can be interpreted as the fluctuation formula is valid for fluctuations with probabilities larger than $`10^9`$.

In the rest of the paper, we present our numerical results for the GIK thermostatted LG both in two and three dimensions, with various values of the external fields. We focus on the question whether nonreversible dynamics leads to different scaling in the fluctuations than the one found in reversible systems. As a general rule, we have not found any difference between time reversal symmetric cases (i.e., with $`𝐁=0`$) and reversible ones. Indeed, time reversal symmetry can be replaced in the known fluctuation theorems by general reversibility without affecting their validity, since the proofs do not make use of the special form of the involution $`𝐆_0`$. We note that we have tested the different dynamical cases with several choices for the field stregths and could not find significant deviations in the observed behavior as long as we stayed within the ergodic region of the parameter space.

In Figs. 3 and 4, we plot $`D_\tau \left(x\right)`$ for reversible dynamics in two and three dimensions (2D and 3D), respectively. The fluctuation formula appears to be valid in both cases; the convergence to the linear limit, however, seems to be different in them. For 2D, the $`D_\tau \left(x\right)`$ curve has deviations, decreasing in size with $`\tau `$ increasing, from the linear shape, while for 3D, $`D_\tau \left(x\right)`$ exhibits strongly linear behavior with slopes approaching 1 as $`\tau `$ increases. It is worth noting that the latter convergence can also be observed in the 2D random LG unpub .

Our results for the nonreversible versions are shown in Figs. 5 and 6.

The most striking difference compared to the reversible cases is the fact that there seems to be a cubic term present in $`D_\tau \left(x\right)`$ that does not disappear for larger $`\tau `$ values. This term leads to a breakoff from the diagonal line for $`|x|x_c0.3`$, which means that there can be deviations from the fluctuation formula for large fluctuations. The slope of the linear part, however, is still 1 in the large $`\tau `$ limit, so the FF can be a good approximation for small to moderate size fluctuations. The fact that the region of validity of the FF does not shrink considerably for larger $`\tau `$ values suggests that the coefficient of the cubic term in $`D_\tau (x)`$ may have only weak dependence on $`\tau `$. This also means that as the distribution $`\mathrm{\Pi }_\tau (x)`$ is concentrating around its mean value for increasing $`\tau `$ values, the total statistical weight of the large fluctuations that are not covered by the linear regime is decreasing. In other words, the FF becomes more and more valid in a probabilistic sense as $`\tau \mathrm{}`$, since the larger fluctuations become less and less likely in that limit.

We may conclude that the FF appears to be valid in the GIK thermostated LG with reversible dynamics, both in two and three dimensions. For nonreversible dynamics, we have found indications that the FF may still describe the scaling properties of fluctuations in a moderate size regime, although for large fluctuations there are clear deviations from it due to higher order terms in $`D_\tau \left(x\right)`$. The fact that the slope of the linear part in the scaling behavior is the same in reversible and nonreversible cases suggests a kind of robustness for the FF in the thermostated LG. It would be interesting to see if this remains valid in other nonequlibrium systems with nonreversible dynamics.

###### Acknowledgements.

The authors are grateful to Tamás Tél for fruitful discussions and a careful reading of the manuscript. This work was supported by the Bolyai János Research Grant of the Hungarian Academy of Sciences and by the Hungarian Scientific Research Foundation (Grant No. OTKA T032981).


text:
# THE MUON ANOMALOUS MAGNETIC MOMENT - SIGNIFICANCE OF THE NEW MEASUREMENT

## 0.1 INTRODUCTION

The beautiful experiment measuring the muon anomalous magnetic moment (g-2) now in progress at Brookhaven  has already attained a precision that challenges our ability to calculate the expected result of the experiment. The uncertainty of the Brookhaven g-2 measurement is now about 1.5 in units of $`10^{10}a_\mu `$ ( where $`a_\mu =\frac{g2}{2}`$), with 1 year of data yet to be analyzed according to reference 1. Weak and electromagnetic contributions to g-2 can be calculated according to a well-understood theory with about 1/30’th of the uncertainty of the measurement, as is noted in the reference. The difference between that calculation and the measurement is 71.8 with the uncertainty of the experimental result. The difference is ascribed to hadronic effects (see  for a recent review).

The significance of the new result is, in my view, that it measures the hadronic contribution to the muon g-2 with a precision of the order of 1 %. The precision may increase by about an order of magnitude when the experiment is completed. There is no theory of hadrons that predicts hadronic interactions so precisely. The appropriate theory should presumably be grounded in quantum chromodynamics (QCD), which does not yet lend itself to such precise predictions. There are phenomenological arguments which relate the hadronic contribution to g-2 to other experiments and provide the basis for current estimates which range  from 69.2 to 72.5, each with a quoted uncertainty of the order of a percent.

The quoted uncertainties of current estimates reflect only the uncertainties of the experimental data used for phenomenological estimates and the computational uncertainties of the estimates. They do not, in my view, adequately reflect the uncertainties in the theories underlying the phenomenological estimates, as I shall explain in the remainder of this note.

## 0.2 DESCRIPTION OF THE HADRONIC CONTRIBUTION

In the notation of Sommerfield , the muon propagator $`G_{mu}`$ relates to a mass function $`M`$ as

$$G_{mu}^1=\gamma \mathrm{\Pi }+M$$
(1)
where $`\mathrm{\Pi }_\mu =(peA)_\mu `$, $`p`$ is the muon momentum and $`A`$ is the potential for a weak constant external field. The anomalous g-factor occurs as a term proportional to $`\sigma _{\mu \nu }F_{\mu \nu }`$ in the mass function $`M`$. That function is calculated from the equation

$$M=m_0+ie^2\gamma G_{mu}\mathrm{\Gamma }𝒢$$
(2)
involving the photon Green’s function $`𝒢`$ and the “dressed photon-muon vertex ” $`\mathrm{\Gamma }`$.

The photon Green’s function is assumed  to obey a Källen-Lehman representation . According to that representation each element $`𝒢_{\mu ,\nu }`$ may be written as a dispersion integral

$$𝒢_{\mu ,\nu }(k^2)=\frac{\delta _{\mu ,\nu }}{k^2}[\delta _{\mu ,\nu }\frac{k_\mu k_\nu }{k^2}]_0^{\mathrm{}}\frac{s(x)}{x+k^2iϵ}𝑑x.$$
(3)
The weight function $`s(x)`$ is proportional to the sum over possible intermediate states produced by the electromagnetic current, according to the equation

$$k^4s(k^2)\delta _{\mu \nu }=\frac{1}{3}\underset{\alpha }{}<0j_\mu (0)\alpha ><\alpha j_\nu (0)0>\delta ^4(\alpha k)$$
(4)
where $`j`$ is the electromagnetic current operator. The hadronic contribution to $`𝒢`$ is then given by that part of $`s(x)`$ for which the states $`\alpha >`$ consist of hadrons. The matrix elements $`<0j_\mu \alpha >`$ are the amplitudes for production of the state $`\alpha >`$ from the vacuum by the electromagnetic current. The corresponding contributions to $`s(x)`$ are therefore proportional, to lowest order in the fine-structure constant, to the cross-sections for producing states $`\alpha `$ in $`e^+e^{}`$ collisions. It would accordingly seem that such measured cross-sections could be used to calculate the hadronic contributions to $`𝒢`$.

The components $`\mathrm{\Gamma }_\mu `$ of the vertex function $`\mathrm{\Gamma }`$ contain terms corresponding to the production of a virtual hadron that couples to the muon. The lowest order such terms (in powers of the charge $`e`$) - the so-called “light-light scattering” terms - are of order $`e^6`$. Power counting suggests that such terms contribute to the anomalous muon moment at a level of about $`\frac{1}{10}`$th of the uncertainty in the Brookhaven measurement. There is, however, no theory comparable to electroweak theory for calculating these terms, as Hayakawa and Kinoshita have emphasized  in connection with their own estimates. Those estimates are not severely tested at the present level of experimental accuracy.

## 0.3 UNCERTAINTIES IN HADRONIC-CONTRIBUTION ESTIMATES

The Källen-Lehman representation  for the photon propagator is of the nature of a dispersion relation. The dispersion relation is a statement that $`𝒢_{\mu ,\nu }(k^2)`$ is a real analytic function of the variable $`k^2`$, cut along the negative (that is, timelike) $`k^2`$ axis . It is also, however, a statement that $`𝒢`$ is a polynomially-bounded function for large, complex $`k^2`$. I contend that the latter statement is unprovable . There does not, in any event, appear to be any existing proof that $`𝒢`$ obeys such a boundedness property. Any attempt to assign a numerical uncertainty to a calculation based upon a Källen-Lehman representation for the photon propagator must therefore be regarded as speculative.

There are recent estimates of the hadronic contributions to the muon g-2 that are within a few per cent of values that one would deduce from the new experimental result . Is this near agreement evidence for the polynomial boundedness of $`𝒢`$? I present in the next section models, by way of counter-examples, where such near agreement can occur even with functions that grow exponentially. Polynomial boundedness cannot, therefore, be deduced from the present data. Uncertainty estimates based upon dispersion-relation calculations of $`𝒢`$ must be understood as lower limits because the theories underlying the calculations are incompletely defined.

## 0.4 Dispersion Relations for Functions Having Exponential Growth

I present two models to show that an attempt to calculate a function by a dispersion relation, when the assumption of polynomial boundedness is incorrect, need not lead to a result that is wildly different from the correct result. The fact that dispersive calculations of the hadronic component of the muon g-2 are close to the values deduced from experiment is, therefore, not evidence that such calculations can give correct results.

Both of the models happen to involve functions related to Bessel functions. The essence of the demonstration, however, may be made with the two elementary functions natural logarithm and exponential. Consider, for this purpose the following function $`Q(z)`$ of a complex variable $`z`$:

$$Q(z)=\mathrm{ln}(z)+e^{bz}$$
(5)
where $`b`$ is some real positive constant. $`Q(z)`$ is an analytic function of the complex variable $`z`$ except for logarithmic branch points at $`0`$ and $`\mathrm{}`$ and an essential singularity (of the exponential) at $`\mathrm{}`$. $`Q`$ is made single-valued by a cut along the positive real axis, giving $`Q`$ an imaginary part equal to $`i\pi `$ as the real axis is approached from above.

The point of the models is that they involve functions that, like $`Q`$, have an essential singularity and are therefore not polynomially bounded. One may nevertheless try write a dispersion relation for $`Q`$ by integrating around the cut, because the imaginary part, being a constant in this case, is well-behaved. The result of the dispersion integral will be to recover only the $`\mathrm{ln}`$ part of $`Q`$. The remaining part of $`Q`$, the exponential, will be negligible (compared with the logarithm) almost everywhere along the positive real axis, depending upon the size of the parameter $`b`$. The erroneous assumption that $`Q`$ satisfies a dispersion relation therefore makes only a small, model-dependent error almost everywhere along the positive real axis.

### 0.4.1 A Hankel-Function Model

In the first model I assume that the correct function that I am trying to calculate is the Hankel function $`H_0^{(1)}(z)`$. The imaginary part of the Hankel function is the Neumann function $`Y_0(z)`$. Both functions emulate a coulomb-like singularity at $`z=0`$, albeit only logarithmically, so a numerical integration of a dispersion relation must be cut off at small $`z`$. The function $`H_0^{(1)}(z)`$ behaves for large, complex $`z`$ like $`\frac{e^{iz}}{\sqrt{z}}`$, and is not, therefore, polynomially bounded.

The model is unrealistic in the sense that the imaginary part of the “true function” is oscillatory, unlike the strictly positive cross-section data that is used as input in the photon-propagator calculations. That particular unrealistic feature is corrected in the second model.

The game, then, is to pretend that we know from experiment the imaginary part of the “correct” function. We then put the “experimental” data into an unsubtracted dispersion relation under the erroneous assumption that the “correct” function vanishes sufficiently rapidly at large $`|z|`$, cutting off the integration at some small value of $`z`$. We want to determine how much our answer deviates from the correct answer which is, in this case, the Bessel function $`J_0(z)`$.

The answer can be determined analytically, because

$$\frac{1}{\pi }_0^{\mathrm{}}\frac{Y_0(x)dx}{xz}=J_0(z)\frac{2}{\pi ^2}S_{1,0}(z)$$
(6)
where the principal value of the integral is intended and $`S_{1,0}(z)`$ is a Lommel function whose explicit form has been given by Watson . It is a solution of the inhomogeneous Bessel Equation

$$S^{\prime \prime }+\frac{1}{z}S^{}+S=\frac{1}{z^2}$$
(7)
subject to the condition that it behaves asymptotically like $`\frac{1}{z^2}`$. It grows for small $`z`$ like $`ln^2(z)`$, but rapidly becomes small as $`z`$ increases. Fig. 1 compares the right hand side of Eq. 6, shown by the solid line, with the “correct” answer $`J_0(z)`$, shown by the dotted line. It is apparent that the error from using misusing the dispersion integral can be made quite small if small values of $`z`$ are excluded, as they are in “real life” when coulomb and infra-red effects play a role.

### 0.4.2 A Strictly Positive Imaginary Part

A more realistic model would require that the simulated experimental data be given by a function that is strictly positive, like a cross-section. A convenient choice is the modified Hankel function $`K_0(z)`$ which approaches zero asymptotically for large $`z`$ like $`\frac{e^z}{\sqrt{z}}`$ and grows logarithmically for small $`z`$ like $`ln(z)`$. For $`z`$ approaching the positive real axis from above, there is a function $`f(z)`$ given by

$$f(z)=\frac{1}{\pi }_0^{\mathrm{}}K_0(x)𝑑x(\frac{1}{xz}\frac{1}{x+z})=\pi (𝐋_\mathrm{𝟎}(z)I_0(z))+iK_0(z)$$
(8)
Here $`𝐋_\mathrm{𝟎}`$ is the modified Struve function and $`I_0`$ is the modified Bessel function, each of order zero.

The complex function $`f(z)`$ does not satisfy a dispersion relation. This is because the integral of the function around a closed loop at large $`z`$ cannot be made arbitrarily small, since its imaginary part, $`K_0(z)`$, grows without bound when $`\frac{\pi }{2}<arg(z)<\frac{3\pi }{2}`$. The imaginary part of $`f(z)`$ again represents the simulated experimental data.

The game, as before, is to put the simulated experimental data into an unsubtracted dispersion relation and compare with the correct answer, given in Eq 8. The result is shown in Fig. 2 which displays the real part of the dispersion result as a solid line and the real part of $`f(s)`$ as a dotted line (the imaginary parts are, of course, identical). It is evident that for values of $`z`$ greater than about unity, the dispersion result is not greatly different from the “correct” value.

## 0.5 Conclusion

Hadronic contributions to g-2 originate, according to our present understanding, from electromagnetic excitation of virtual quark pairs from the vacuum. The task of calculating the effects of such excitations therefore falls within the realm of Quantum Chromodynamics (“QCD”). Such calculations have not yet been done. It has been tempting to suppose, as an alternative, that dispersion relations using empirical data could be used to predict the results of QCD calculations, at least for the vacuum polarization corrections to the photon propagator.

The difficulty with the supposition is that a dispersion relation involves an assumption about the behavior of a function for large values of its argument, namely, that the function is at least polynomially bounded. In the case of the photon propagator there is no way to justify that assumption.

It might be tempting to try to evade this lack of justification by reversing the burden of proof, asking: ”What is the physical origin of the analogous terms in the physically relevant case? In particular, if Exp\[-z\] —$`>`$ Exp\[-E/E0\], what is the meaning and value of the energy scale E0?”

Because the presence of an essential singularity at infinite energy, as I have modeled with an exponential term, is purely conjectural, any mass, such as the pion mass, the QCD scale or a spontaneously generated mass  is a possible candidate for the energy scale E0. One doesn’t see such terms in the usual discussions of analytic properties of amplitudes because those discussions are customarily guided by perturbation theory  , where such terms cannot occur. But we know of models where perturbation theory breaks down. . And essential singularities are a common feature of solutions of differential equations such as the relativistic Schroedinger equation, cf. .

Proof of the existence of an essential singularity would, of course, end the discussion; dispersion relations would not provide an algorithm for calculating the hadronic contribution to vacuum polarization. That is not the intent of this note, however. The intent is merely to point out that there is, as yet, no mathematical justification for the use of dispersion relations for that purpose, and the accuracy of the resulting predictions must therefore remain uncertain.

What I have shown here is that violation of the assumption of polynomial boundedness can, in some circumstances, lead to results that are close to the correct results. There is, however, no apparent a priori way to estimate the size of the error resulting from the violation. There is consequently no way to estimate the error involved in using a dispersion relation to calculate the hadronic contribution in question.

The Brown, et al. measurement of the muon g-2 is therefore significant because it measures the hadronic contribution to g-2 with unprecedented precision. That measurement stands as a challenge to our understanding of the quark structure of the vacuum.

I am indebted to Geoffrey Bodwin, Tim Tait, Carlos Wagner, Alan White, and Cosmas Zachos for some helpful discussions, to Professor T. T. Wu for his assistance in identifying reference , and to the (anonymous) referee, whose question led me to become aware of references , ,  and several related papers. This work was supported by the U.S. Department of Energy, Division of High Energy Physics, Contract W-31-109-ENG-38.


text:
# VISCOUS COSMOLOGIES WITH EXTRA DIMENSIONS

## Abstract

We present an analysis of a n-dimensional vacuum Einstein field equations in which 4-dimensional space-time is described by a Friedmann Robertson-Walker (FRW) metric and that of the extra dimensions by a Kasner type Euclidean metric. The field equations are interpreted as four dimensional Einstein equations with effective matter properties. The effective matter is then treated as viscous fluid. We consider the theories of imperfect fluid given by Eckart, truncated Israel-Stewart and full Israel-Stewart theories and obtain cosmological solutions for a flat model of the universe. The imperfect fluid assumption admits here power law inflation for the 3-physical space in some cases.

PACS numbers : 04.20 Jb, 98.80 H

Model building in higher dimensions was initiated by Kaluza and Klein1 who tried to unify gravity with electromagnetic interaction by introducing an extra dimension. Kaluza-Klein theory is basically an extension of Einstein general relativity in 5D which is of much interest in Particle Physics and Cosmology. In the last few decades the study of a higher dimensional theory has been revived and considerably generalized after realizing that many interesting theories of particle interactions need more than four dimensions for their formulation. Attempts have been made to build cosmological models2 in higher dimensions which may undergo a spontaneous compactification leading to a product space $`M^4\times M^d`$, with $`M^d`$ describing the compact $`\mathrm{𝑖𝑛𝑛𝑒𝑟}\mathrm{𝑠𝑝𝑎𝑐𝑒}`$. In the usual approach one uses Einstein’s field equation

$$G_{ab}=8\pi GT_{ab}$$
(1)
which match the Einstein tensor and the energy momentum tensor. Recently another approach3 has been developed taking $`G_{ab}=0`$, the extra terms which appear in these equations due to extra dimensions are interpreted as induced or effective properties of matter in ordinary 4D space-time. Ibánez and Verdaguer \[henceforth IV\]4 obtained a set of solutions of Einstein’s equation in an n-dimensional vacuum. They found homogeneous solutions with expanding three dimensional isotropic space. The solutions are then identified with the observable four dimensional subspace with perfect fluid. Later a class of solutions in flat universe model have been obtained by Gleiser and Diaz.5 Krori et al.6 using a higher dimensional anisotropic cosmology (Bianchi-I) obtained 4D perfect fluid solutions. It was shown that the perfect fluid solutions obtained by them are compatible with contraction of all the extra dimensions. The solutions are identified with a 4 dimensional anisotropic cosmological model. However, it is known that the matter distribution in the early universe may not be as simple as predicted by perfect fluid due to a number of dissipative processes. Hence it is essential to study cosmological solutions in the framework of imperfect fluid description. The simple theory is the Eckart theory,7 however, it suffers from serious drawbacks concerning causality and stability.8,9 These difficulties can, however, be removed by considering higher order theories i.e., extended irreversible thermodynamics (EIT).10,11 In this paper we intend to study a higher dimensional vacuum Einstein equations which leads to an observed 4 dimensional universe with imperfect fluid. We obtain exact solutions in the framework of the Eckart, truncated Israel-Stewart (TIS) and full Israel-Stewart (FIS) theories.

We now consider a higher dimensional metric in the form

$$ds^2=ds_{FRW}^2+\underset{i=1}{\overset{d}{}}b_i^2(t)dx_i^2$$
(2)
where $`d`$ is the number of extra dimensions $`(d=n4)`$ and $`ds_{FRW}^2`$ represents the line element of the FRW metric in four dimensions which is given by

$$ds_{FRW}^2=dt^2+a^2(t)\left[\frac{dr^2}{1kr^2}+r^2\left(d\theta ^2+sin^2\theta d\varphi ^2\right)\right]$$
(3)
where $`a(t)`$ is the scale factor of the 4 dimensional spacetime and $`b_i`$’s are the scale factors for the extra dimensions and $`k=0,+1,1`$ represents flat, closed and open universe respectively.

The vacuum Einstein’s equation in n-dimensions using the metric (2) can be written as

$$3\left(\frac{\dot{a}^2+k}{a^2}\right)=\frac{\ddot{D}}{D}\frac{1}{8}\left(\underset{i=1}{\overset{d}{}}\frac{\dot{b}_i}{b_i}\right)^2+\frac{1}{8}\underset{i=1}{\overset{d}{}}\left(\frac{\dot{b}_i}{b_i}\right)^2,$$
(4)
$$2\frac{\ddot{a}}{a}+\frac{\dot{a}^2+k}{a^2}=\frac{\dot{a}\dot{D}}{aD}+\frac{1}{8}\left(\underset{i=1}{\overset{d}{}}\frac{\dot{b}_i}{b_i}\right)^2\frac{1}{8}\underset{i=1}{\overset{d}{}}\left(\frac{\dot{b}_i}{b_i}\right)^2,$$
(5)
$$\frac{\ddot{b_i}}{b_i}+3\frac{\dot{a}\dot{b}_i}{ab_i}+\frac{\dot{D}\dot{b}_i}{Db_i}\frac{\dot{b_i}^2}{b_i^2}=0$$
(6)
where we denote $`D^2=\mathrm{\Pi }_{i=1}^db_i(t)`$.

The field eqs. (4)-(6) can be rewritten as the Einstein equations with imperfect fluid in 4 dimensions. Using $`b_i=D^{2p_i}`$ with $`_{i=1}^dp_i=1`$ the field equations reduce to

$$3\left(\frac{\dot{a}^2+k}{a^2}\right)=\rho $$
(7)
$$2\frac{\ddot{a}}{a}+\frac{\dot{a}^2+k}{a^2}=P_{eff}$$
(8)
where $`\rho =\frac{\ddot{D}}{D}+\alpha \left(\frac{\dot{D}}{D}\right)^2`$ , $`P_{eff}=p+\mathrm{\Pi }=\frac{\dot{a}\dot{D}}{aD}+\alpha \left(\frac{\dot{D}}{D}\right)^2`$ and $`\alpha =\frac{1}{2}\left[_{i=1}^dp_i^21\right]`$. In the above $`p`$ represents the thermodynamic pressure and $`\mathrm{\Pi }`$ represents the bulk viscous stress. Using the relation $`p=(\gamma 1)\rho (1\gamma 2)`$ one obtains the bulk viscous stress which is given by

$$\mathrm{\Pi }=\left(\gamma \frac{4}{3}\right)\frac{\ddot{D}}{D}\alpha (\gamma 2)\left(\frac{\dot{D}}{D}\right)^2.$$
(9)
It is evident above that for $`\alpha =0`$ and $`\gamma =\frac{4}{3}`$ one recovers the cosmological solutions for a radiation dominated universe i.e. without viscosity. The solutions are obtained by IV in the context of perfect fluid in 4 dimensions. We study here cosmological solutions for $`\alpha 0`$, which leads to interesting results. The evolution of the bulk viscous stress in (3 + 1) dimensions is given by10

$$\mathrm{\Pi }+\tau \dot{\mathrm{\Pi }}=3\zeta H\frac{ϵ}{2}\tau \mathrm{\Pi }\left(3H+\frac{\dot{\tau }}{\tau }\frac{\dot{\zeta }}{\zeta }\frac{\dot{T}}{T}\right),$$
(10)
where $`\zeta `$ is the coefficient of bulk viscosity, $`\tau `$ is the relaxation time and $`T`$ is the temperature and H ( = $`\frac{\dot{a}}{a}`$ )

is the Hubble parameter. Here the parameter $`ϵ`$ can take the value either 0 or 1; $`ϵ=0`$ for TIS and $`ϵ=1`$ for FIS theory. One obtains Eckart theory for $`\tau =0`$. The system of eqs. (7)-(10) are not closed as number of unknowns are more than the number of equations. It is, therefore, necessary to assume the following adhoc but commonly chosen relations for the bulk viscosity co-efficient and relaxation time :

$$\zeta =\beta \rho ^qand\tau =\beta \rho ^{q1}$$
(11)
where $`q0`$ and $`\beta `$ is a constant parameter. The behavior of temperature $`T`$ of a universe with imperfect fluid is obtained in FIS theory using eq.(10).

We now look for cosmological solutions. Let us consider a power law model of the expansion of the universe given by

$$a=a_ot^m$$
(12)
where $`a_o`$ and m are constants. We consider a flat universe ( $`k=0`$ ) and obtain solutions for Eckart, TIS and FIS theories for two cases $`\gamma =\frac{4}{3}`$ and $`\gamma =2`$ respectively.

Case I : Eckart theory :

Eckart theory corresponds to $`\tau =0`$. The eq.(10) determines the viscous stress, using $`q=\frac{1}{2}`$ i.e., $`\zeta =\beta \rho ^{1/2}`$ we get the following solutions :

(i) For $`\gamma =\frac{4}{3}`$ one gets power law expansion with $`\beta =\frac{2(2m1)}{\sqrt{3}m}`$. The coefficient of bulk viscosity $`\beta `$ is a positive quantity for $`m>\frac{1}{2}`$. The model admits a large power law inflation ($`m\mathrm{}`$) when $`\beta \frac{4}{\sqrt{3}}`$.

(ii) For $`\gamma =2`$, we get $`\beta =\frac{2(3m1)}{3\sqrt{3}m}`$. In this regime a large power law inflation is obtained when $`\beta \frac{2}{\sqrt{3}}`$.

Case II : TIS theory ( $`ϵ=0`$ ) :

(i) For $`\gamma =\frac{4}{3}`$, eq.(9) determines the bulk viscous stress which is given by

$$\mathrm{\Pi }=\frac{2\alpha (13m)^2}{3t^2}.$$
(13)
The bulk viscous stress depends both on $`\alpha `$ and on the exponent $`m`$. For $`m=\frac{1}{3}`$ , the co-efficient of bulk viscosity vanishes. Thus one requires $`\alpha <0`$ and $`m\frac{1}{3}`$ for an acceptable solution. The first constraint in the model is satisfied for $`_{i=1}^dp_i^2<1`$. For simplicity we consider $`q=\frac{1}{2}`$ and obtain $`\beta =\frac{2\sqrt{3}(2m1)m}{9m^2+8m4}`$ from eq.(10). In this case a realistic scenario is obtained for $`m>\frac{1}{2}`$. One gets a large power law inflation when $`\beta \frac{4\sqrt{3}}{9}`$.

(ii) For $`\gamma =2`$, the bulk viscous stress is given by

$$\mathrm{\Pi }=\frac{2m(3m1)}{t^2},$$
(14)
which remains negative for $`m>\frac{1}{3}`$. The bulk viscous stress is independent of $`\alpha `$. For $`q=\frac{1}{2}`$, we get $`\beta =\frac{2\sqrt{3}m(3m1)}{9m^2+12m4}`$, which allows a realistic scenario for $`m>\frac{1}{3}`$. In this case a large power law expansion is derived when $`\beta \frac{2\sqrt{3}}{3}`$.

Case III : FIS theory ( $`ϵ=1`$ ) :

In this case temperature of the universe can be evaluated from the eq.(10).

(i) For radiation ( $`\gamma =\frac{4}{3}`$ ) and viscous fluid model, the bulk viscous stress is evaluated which is given by

$$\mathrm{\Pi }=\frac{2\alpha }{3}\frac{(13m)^2}{t^2}.$$
(15)
For simplicity we take $`q=\frac{1}{2}`$ and obtain the temperature of the universe which is given by

$$T=T_{01}t^{a_1}$$
(16)
where $`a_1=3m4+\frac{2\sqrt{3}m}{\beta }+\frac{27m^3}{\alpha (13m)^2}`$ and $`T_{01}`$ is an integration constant. One gets a decreasing mode of temperature in this regime for $`\beta >\frac{2\sqrt{3}m(2m1)}{3m^2+11m4}`$,

which in turn determines the relaxation time and coefficient of bulk viscosity. Thus a lower bound on viscosity is obtained depending on the exponent of the scale factor. A realistic scenario is obtained for $`m\frac{1}{2}`$. However, for a large power law inflation one requires a large value of the bulk viscosity (as $`\beta >\frac{4\sqrt{3}}{3}`$).

(ii) For a stiff matter ( $`\gamma =2`$ ) and viscous fluid model, the bulk viscous stress is

$$\mathrm{\Pi }=\frac{2m(3m1)}{t^2},$$
(17)
which is independent of $`\alpha `$. Here $`\mathrm{\Pi }<0`$ is obtained for $`m>\frac{1}{3}`$. Using $`q=\frac{1}{2}`$ one gets temperature of the universe which varies as

$$T=T_{02}t^{a_2}$$
(18)
where $`a_2=\frac{415m}{3m1}+\frac{2\sqrt{3}m}{\beta }`$ and $`T_{02}`$ is an integration constant. In this case a decreasing mode of temperature is found for $`\beta >\frac{2\sqrt{3}m(3m1)}{15m4}`$. For a realistic scenario one requires $`m\frac{1}{3}`$ in this case. One gets a large power law expansion for $`\beta >\frac{2\sqrt{3}}{5}`$.

Thus we obtain cosmological solutions of a higher dimensional universe considering matter sector of the theory which arises from the extra dimensions. The evolution of the extra dimensions scale factor is described by Kasner-type behavior. The solutions obtained here are identified with the (3 + 1) dimensional flat (k = 0) universe with imperfect fluid. The solutions are obtained in Eckart, TIS and FIS theories respectively. The causal fluid distribution considered here admit cosmological models with power law inflation as well as models with slow power law expansion of the early universe. The overall evolution of the internal space is contracting in nature ( as $`_{i=1}^dp_i^2<1`$ and $`_{i=1}^dp_i=1`$ ) whereas the 3-physical space expands. We note the following :

$``$ In the TIS or FIS theory one may get a universe with comparatively large power law expansion when $`\gamma =\frac{4}{3}`$ than that when $`\gamma =2`$ for a given bulk viscosity. However for $`\gamma =\frac{4}{3}`$ with $`m=\frac{1}{2}`$ one gets back the solution corresponding to perfect fluid distribution. The usual radiation dominated universe $`\left(a(t)\sqrt{t}\right)`$ is recovered in the case.4

$``$ A large viscosity drives a power law inflation. In this model the viscosity parameter $`\beta `$ is least in TIS theory and greatest in FIS theory for a given expansion of the universe for $`\gamma =\frac{4}{3}`$ as well as for $`\gamma =2`$. Thus in the TIS theory, one gets a large power law inflation for a comparatively lower bulk viscosity. However one gets a range of values of $`\beta `$ in FIS theories which requires a higher bulk viscosity to begin with for $`\gamma =\frac{4}{3}`$ compared to the other models. The viscosity of the universe decreases as the universe evolves.

$``$ The temperature of the universe can be determined explicitly in the FIS theory which is found to have a decreasing mode determined by the viscous fluid distribution. The bulk viscosity of the universe decreases with time with the choice $`q=\frac{1}{2}`$ considered here in eq. (11). One interesting aspect of the solution is that for a given expansion of the universe $`\beta `$ is fixed in Eckart and TIS theories with values $`\beta (Eckart)>\beta (TIS)`$. However one finds a lower bound on the value of viscosity in the FIS theory. In this case a hot universe cools rapidly when $`\beta `$ i.e., viscosity of the universe is large.

Acknowledgements :

I would like to thank the Inter-University Centre for Astronomy and Astrophysics ( IUCAA ) Pune for awarding a Visiting Associateship and providing a facility where this work was carried out. I am grateful to Professor N. Dadhich for his suggestion and cordial discussion. Finally, I like to thank Professor S. Mukherjee for providing $`\mathrm{𝐼𝑈𝐶𝐴𝐴}\mathrm{𝑅𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒}\mathrm{𝐶𝑒𝑛𝑡𝑟𝑒}`$ (IRC) facilities at North Bengal University for completing the work.


text:
# Parasitic Light in NGST instruments: the accuracy of photometric redshifts and the effect of filter leaks in the visible and near IR camera.

## 1 Introduction

The NGST Design Reference Mission (DRM) is a representative science program elaborated in sufficient detail to aid in the development of functional requirements for the NGST mission. One of the main results of the DRM has been the identification of a large fraction of subprograms (8 out of 12) requiring observations at optical wavelengths. Such an extension poses significant technical challenges. In particular it requires an efficient blocking of photons at undesired wavelengths over a range covering more than three octaves. In the following we examine the consequences of leaks in the transmission of filters for a typical programme of determination of photometric redshifts of high-redshift galaxies.

## 2 Baseline of the Instrument

Table 1 shows the main characteristics of the Visible-NIR camera as conceived by the Visible and NIR camera Subcommittee of the “Ad Hoc Science Working Group” (ASWG). Two options are envisaged: an optimal and a bare bones version. More details can be found at the URL
http://www.stecf.org/ngst/instruments.html.

## 3 Band-passes of the filters

We have assumed as reference bandpasses for the filters of the optical-NIR camera the standard $`V606I814z\mathrm{Gunn}F110F160KLM`$ set commonly used on HST with WFPC2 ($`V606I814`$) or NICMOS ($`F110F160`$) and from the ground (Fig. 1). Details about the filter properties are given in Table 2

The actual choice of the filters to be used in space with NGST may be different but the arguments developed in the following do not depend critically on this detail.

## 4 Typical SEDs of the sources and the locus of galaxies in multicolor space

To simulate a range of galaxy SEDs we have adopted the templates of Coleman et al. (1980) for a typical elliptical, Sbc, Scd and irregular galaxy plus a GISSEL (Bruzual and Charlot, 2001) template for a very blue galaxy (with a constant star formation rate and an age of 0.1 Gyr). The SEDs are shown in Fig. 2 and in the following will be referred to as extended CWW templates.

Fig. 3 shows the evolution of the apparent magnitude of typical galaxies of $`10^{11}M_{}`$ in the $`F606`$ and $`L`$ bands and gives a rough idea of the magnitudes of interest. It can be appreciated how an elliptical fades away relatively quickly in the $`F606`$ band, while other star-forming galaxies remain - as observed - significantly brighter.

The upper panel of Fig. 4 shows the evolution of the colors of the galaxies described in the previous section in a $`(F606F814)`$ vs. $`(F110K)`$ plane. The color tracks are limited to $`z<3.5`$ for an elliptical SED and to $`z<5`$ for Sbc galaxies. With this combination of bandpasses it is straightforward to identify star-forming galaxies with $`z>5`$ and ellipticals at $`z>1`$. The lower panel shows the same in a $`(F110F160)`$ vs. $`(F160K)`$ plane and demonstrates how, without optical bands, it becomes difficult to isolate star-forming galaxies at $`z<9`$.

## 5 Photometric redshifts and the effects of “(red) leaks”

In order to quantify the different performance of different sets of filters and the effects of “red-leaks” in a filter configuration extending from the visible to the infrared we have carried out simulations of the assignment of photometric redshifts to a sample of faint galaxies. To this end, we have tried to generate a simulated photometric catalog covering a very large range of galaxy spectral properties and to estimate photometric redshifts from an economic and independent a set of templates as possible. Galaxies were simulated in a redshift interval $`0<z<20`$ according to GISSEL library. This spectral synthesis model is controlled by a number of free parameters. The star formation rate for a galaxy at a given age is governed by the assumed $`efolding`$ star formation time-scale $`\tau `$. Several values of $`\tau `$ and galaxy ages are necessary to reproduce the different observed spectral types. For the shape of the initial mass function (IMF) we have assumed a Salpeter law, which turns out not to be a critical choice. In addition to the GISSEL parameterization, we have added the internal reddening of each galaxy by applying the observed attenuation law for local starburst galaxies by Calzetti et al. (1997), with the $`E_{BV}`$ normalization varying between 0 and 0.2 mag. We have also included the Lyman absorption produced by the intergalactic medium as a function of redshift. More details on this type of simulations can be found in Arnouts et al. (1999). Fig. 5 shows the galaxy counts as a function of the apparent magnitude, the age and redshift distribution, the $`zmag`$ plane characterizing the simulation.

Magnitudes and colors have been derived from the simulated SEDs by convolution with the filter passbands described in Sect.3. Measurement errors have been added, according to the NGST mission simulator (http://www.ngst.stsci.edu/nms/main/),
assuming an exposure time of $`10^4`$s.

From this set of “simulated observed colors” photometric redshifts have been derived using a standard procedure of $`\chi ^2`$ optimization, comparing the observed fluxes with errors to the extended CWW set of templates described in Sect. 4.

The $`\chi ^2`$ has been computed as:

$$\chi ^2=\underset{i}{}\left[\frac{F_{observed,i}AF_{template,i}}{\sigma _i}\right]^2$$
(1)
where $`F_{observed,i}`$ is the flux observed in a given filter $`i`$, $`\sigma _i`$ is its uncertainty, $`F_{template,i}`$ is the flux of the template in the same filter and the sum is over the various filters. The template fluxes have been normalized to the observed ones with the factor $`A`$ minimizing the $`\chi ^2`$ value :

$$A=\underset{i}{}\left[\frac{F_{observed,i}F_{template,i}}{\sigma _i^2}\right]/\underset{i}{}\left[\frac{F_{template,i}^2}{\sigma _i^2}\right]$$
(2)
Fig. 6 shows the comparison between photometric redshift estimates obtained with a $`V606I814`$ $`z\mathrm{Gunn}F110F160KLM`$ filter set and a filter set limited to the infrared bands. As suggested in the previous section in order to correctly identify star forming galaxies in the redshift range $`5<z<9`$ and avoid confusion with other SEDs it appears mandatory to have, besides the IR filters, photometric information in the optical bands. In the following we will adopt the set $`V606I814z\mathrm{Gunn}F110F160KLM`$ as standard.

In order to explore the effects of filter leaks the following approach was adopted:

1. two types of leaks have been foreseen: a constant leak at a level of $`10^4`$ of the peak transmission over the whole spectral range and a leak of Gaussian shape placed at 1.5 times the effective wavelength of the filter with an amplitude of $`10^3`$ of the peak transmission and a width of 5% of the central wavelength. The effect of Gaussian leaks on the filter transmission is shown in Fig. 9
2. The leaks have been included in the system spectral response and new magnitudes and colors have been produced with the GISSEL models.
3. The new colors have been used to produce photometric redshifts with the extended CWW templates. No knowledge of the leaks was used in this step, i.e. the system responses assumed in the z-phot estimation are those of Fig. 1, without leaks.

The typical errors deriving from photometric uncertainties and unaccounted leaks are illustrated in Fig. 7-8 and 10-11. The dot-dashed, long-dashed, short-dashed and continuous lines correspond to the 95, 90, 68 and 50 percentile of the absolute deviation $`|z_{sim}z_{phot}|`$, respectively. It can be seen that errors are in all cases relatively large at $`z5`$. This is an obvious effect due to the lack of bandpasses bluer than $`V_{606}`$, and the consequent difficulty to identify the Lyman break at relatively low redshift and distinguish it, for example, from the 4000 Å break at lower redshift.

The effect of the leaks, despite the rather pessimistic assumptions, is in general a moderate increase of the error level. For example for $`z=11`$ and $`M_{AB}<30`$ the 68 percentile corresponds to $`|z_{sim}z_{phot}|=0.23`$ with the standard photometric errors, while it increases at $`|z_{sim}z_{phot}|=0.25`$ with the unaccounted constant leak at a level of $`10^4`$ of the peak transmission over the whole spectral range. No qualitative change of the accuracy is observed when the leaks are introduced. No relevant difference is also found between the continuous and the Gaussian leak cases. Again for $`z=11`$ and $`M_{AB}<30`$ the 68 percentile corresponds to $`|z_{sim}z_{phot}|=0.27`$ when the leak of Gaussian shape is added. In the range $`5.5<z<12.5`$ a rather good accuracy is achieved with a typical error $`|z_{sim}z_{phot}|0.2`$ in all cases down to $`M_{AB}30`$. The situation deteriorates slightly around $`z15`$, but the result is due mainly to the insufficient sampling of the SEDs between the $`K_s`$ and the $`M`$ filters, in particular to the “hole” at 2.9 $`\mu `$m.

The importance of the leaks, relative to the uncertainty due to photometric errors, decreases with decreasing flux. At $`M_{AB}<30`$ and $`z=11`$ $`|z_{sim}z_{phot}|=0.18,0.23,0.23`$ for the no-leak, constant-leak and Gaussian leak case, respectively.

## 6 Conclusions

A detailed analysis of NGST multiband photometry applied to the reference case of the study of high-redshift galaxies has been carried out with simulations based on galaxy SEDs derived from the currently available empirical and model templates and on plausible standard filter-sets.

In order to correctly identify star forming galaxies in the redshift range $`5<z<9`$ and early-type galaxies above $`z>2`$ and avoid confusion with other SEDs it is mandatory to have photometric information in optical bands, besides a standard IR set like $`F110F160KLM`$. In particular by adding the $`V606,I814`$ and $`z\mathrm{Gunn}`$ filters a good discrimination is obtained above $`z>5`$ for star forming galaxies and $`z>1`$ for early-types. The case for an extension of the NGST wavelength domain to the optical range is therefore strongly supported by this analysis.

The effects of leaks in the filter blocking have also been investigated. In spite of rather pessimistic assumptions (a constant leak at a level of $`10^4`$ of the peak transmission over the whole spectral range or a leak of Gaussian shape placed at 1.5 times the effective wavelength of the filter with an amplitude of $`10^3`$ of the peak transmission and a width of 5% of the central wavelength) the effects are not dramatic: the accuracy of the determination of photometric redshifts with a standard $`V606I814z\mathrm{Gunn}F110F160KLM`$ filter set is not significantly deteriorated for a sample limited to $`M_{AB}=28`$. In the range $`5.5<z<12.5`$ a rather good accuracy a typical error $`|z_{sim}z_{phot}|0.2`$ is achieved.

As standard specifications for leaks are $`<10^4`$ and as any larger leaks could be measured once the filters are produced, it does not look as though leaks constitute a serious problem, at least for the particular application described here. The interpretation of fields containing highly reddened sources may, however, be more seriously compromised. Extended and continuous spectral coverage appears to be the driving factor to maximize the scientific output of the NGST camera.

## 7 Acknowledgments

We are grateful to A.Moorwood for enlightening discussions about the performances of IR filters and detectors and a critical reading of the manuscript.