N Pauget

F Pistolesi

M Houzet

Coulomb blockade for an oscillating tunnel junction

Keywords: numbers: 73.23.-b, 85.85.+j, 74.50.+r

We consider a tunnel junction formed between a fixed electrode and an oscillating one. Accumulation of the charge on the junction capacitor induces a force on the nano-mechanical oscillator. The junction is voltage biased and connected in series with an impedance Z(ω). We discuss how the picture of Coulomb blockade is modified by the presence of the oscillator. Quantum fluctuations of the mechanical oscillator modify the I-V characteristics particularly in the strong Coulomb blockade limit. We show that the oscillator can be taken into account by a simple modification of the effective impedance of the circuit. We discuss in some details the case of a single inductance Z(ω) = iLω and of a constant resistance Z(ω) = R. With little modifications the theory applies also to incoherent transport in Josephson junctions in the tunneling limit.

I. INTRODUCTION

The improved ability to build electronic devices on the nanometer scale opens the perspective to exploit and study the coupling between electronic transport and mechanical degrees of freedom. 1,2 The elementary system typically considered in nano-electromechanics is the harmonic oscillator coupled in some way to the electronic degrees of freedom. A widely studied device is the singleelectron transistor which consists in a small metallic island connected to the leads through two tunnel barriers. The oscillating part of the device can be either a nearby gate electrode, or the island itself. The motion of the device modifies mainly two quantities: the capacitances and the (bare) tunneling rates. According to which dependence dominates, several effects have been predicted or observed. When only the capacitance is modified by the oscillation, it has been shown in the weak coupling limit that the variance of the oscillator position satisfies the equipartition theorem with the voltage bias replacing the temperature. 3 Under specific conditions, the system may also undergo dynamical instabilities. 4,5,6 In the strong coupling limit, a new kind of current blockade (Frank-Condon blockade) has been predicted. 7,8,9,10 This effect also persists in the classical limit. 11,12,13 When the oscillation modifies the distance between the metallic leads, it has been shown that a dynamical instability can occur. 14,15,16,17,18,19,20 This instability (called shuttle instability) is characterized by a synchronization of the electronic transport with the mechanical oscillations, as is clearly seen from the full counting statistics of charge transport. 21,22 Observation of shuttling is difficult, since the distance between the metallic leads must be comparable with the tunneling length in order to produce a measurable current. 23,24,25,26 Very recently indications of shuttling have been found by the authors of Ref. 27 probably in an intermediate regime between tunneling and field emission. The wide interest for the single electron transistor comes from its high sensitivity to a small variation of the gate voltage. Thus, it can be used, for instance, as a sensitive displacement detector 28 and it has been shown to reach the quantum back action limit. 29 An even simpler device is the single tunnel junction where the distance between the two metals constitutes the mechanical degree of freedom. The effect of the position dependence of the tunneling resistance has been already considered in the literature, 30,31,32 both for the current and the noise of the device. In this case, the force acting on the oscillator originates from the electron-momentum transfer. Recently, the position fluctuation has been detected by measuring the current fluctuations. 33 This experiment probed the back action of the current crossing the tunnel junction on the oscillator. Apparently, the position dependence of the tunneling matrix elements is not sufficient to explain the intensity of the back action. To our knowledge, the effect of the force coming from the variation of the capacitance C as a function of the distance between the leads has not been considered in this context so far. A particularly interesting case is when the device is voltage biased in series with an impedance Z(ω) leading, under certain conditions, to the Coulomb blockade physics. [START_REF] Ingold | Single Charge Tunneling[END_REF] Then, similarly to what happens for the single electron transistor, one can expect that a small position dependence of the capacitance can affect the current-voltage characteristics of the device.

In this paper, we consider the effect of the mechanical oscillator on the I-V curve. We will work within the assumptions of the standard Coulomb blockade theory. [START_REF] Ingold | Single Charge Tunneling[END_REF] Specifically, we assume that the tunneling resistance is sufficiently large so that the oscillator and the electromagnetic environment have time to relax to equilibrium between two tunneling events. This is the opposite limit with respect to the one considered in Refs. 31,32 for the tunneling dependent Hamiltonian. There, the tunneling is so frequent that the system may reach a stationary off-equilibrium state. Restricting to the linear coupling of the oscillator with the electromagnetic modes, we will discuss the effect of the quantum fluctuations of the mechanical oscillator on the current.

The paper is organized as follows. In Sect. II we describe the model and its range of validity. We show that within a linear approximation the mechanical oscillator can be accounted by an effective impedance. In Sect. III we consider two specific cases of Z(ω) and discuss the resulting conductance at zero temperature. In Sect. IV we consider the case of an oscillating Josephson junction by extending the results of previous sections to the superconducting case. Section V gives our conclusions.

II. MODEL

We consider the electronic transport through a circuit composed of a tunnel junction in series with an arbitrary impedance Z(ω) and biased with an ideal voltage source V . We assume that one of the two metallic leads forming the tunnel junction is free to oscillate and can thus modify its distance with respect the other one. The single degree of freedom describing the oscillation is given by x, the oscillator position. Its effective mass and spring constant are m and k, respectively. The tunneling surface is very small. Thus, it is important to account for the effective capacitor forming at the tunnel barrier. Its capacitance C(x) depends on the position of the oscillator. The tunneling resistance R T (x) may also depend on x. In the following we will neglect this dependence by showing that its effect is negligible with respect to that of C(x) in the regime of rare tunneling events that we are considering. The device is depicted schematically in Fig. 1.

A. Time scales

As anticipated in the Introduction we will work within the same hypothesis that are commonly used to describe Coulomb blockade. Specifically, we assume that the relaxation times for the electromagnetic environment, 1/γ RC = RC, and for the oscillator, 1/γ, are shorter than the average time between two tunneling events, τ T ∼ eR T /V . Here, R ∼ Re[Z] is the typical resistance associated with the impedance Z and γ is the damping coefficient of the mechanical oscillator with frequency ω m = k/m. If 1/γ RC , 1/γ ≪ τ T , the system has time to come back to the equilibrium state after an electron has crossed the tunnel barrier. This allows us to calculate

Z V m k C(x), R T (x) x I L1 L2 L3 Li C1 C2 C3 Ci ∼ = FIG.
1: Schematic representation of the system. A tunneling contact is formed between an oscillating electrode (here depicted as a STM tip) and a fixed surface. The oscillator is modelled by a spring constant k and an effective mass m.

The junction between the tip and the plate is characterized by an effective capacitance C(x) and resistance RT (x), where x denotes the position of the tip. The circuit is closed by an impedance Z and a voltage source V .

the average current on the equilibrium quantum state of the environment. It gives the following limitations on the voltage: eV ≪ E C (R T /R) and eV ≪ e 2 γR T , where

E C = e 2 /2C
is the Coulomb energy. Since we are interested in the range of voltages around E C /e we assume that R T ≫ R and R T ≫ 1/Cγ. The first condition is the standard tunneling condition, the second involves the mechanical damping. One should remember that the damping of the oscillator also depends on the impedance of the circuit. Thus, even if the oscillator is not damped by a mechanical source, its quality factor would not be infinite. In the limit of large R T the above conditions can thus always be satisfied. Let us now discuss the effect of the position dependence of the tunneling resistance. To evaluate the importance of this effect for the dynamics of the junction we consider the effective friction coefficient generated by the position dependence of R T :

31 γ T ≃ R Q R T x o λ T 2 ω m (1)
where R Q = 2π /e 2 is the quantum of resistance, x 2 o = /2mω m is the (square of the) zero point motion amplitude of the oscillator and λ T is the tunneling length defined by λ -1 T = (dR T /dx)/R T . The explicit expression of the friction coefficient induced by the electromagnetic dissipation is given below for a specific case [cf. Eq. ( 28)]. Independently on its explicit form, what matters is that it does not depend on the tunneling resistance since the damping is due to the environment and the capacitive coupling. Thus the ratio γ T /γ can always be made small by choosing a large enough tunneling resistance. Then we may safely neglect the position dependence of the tunneling matrix elements. We will come back to this point at the end of Section III.

B. Derivation of the general expression for the current

The Hamiltonian describing the system is the following:

H = H j + H osc + H EM .
(

) 2
where

H j = k ξ k c † k c k + q ξ q c † q c q + H T with H T = A + A † and A = k,q t kq e iφ(t) c † k c q .
Here, k and q label the femionic eigenstates in the fixed and oscillating electrode, respectively, c k is a destruction operator in state k, with energy ξ k , and t kq the tunneling matrix element. The phase φ(t) is obtained by integrating the voltage difference at the junction:

φ(t) = t -∞ (eV J (t ′ )/ )dt ′ .
For a constant voltage biased junction φ(t) = eV t/ . However, due to the impedance in series with the junction, the voltage fluctuates and, in general, V J is different from V . The Coulomb blockade physics stems from the behavior of the phase φ and, specifically, from its correlation function. The current through the junction can be obtained by perturbation theory in H T . The operator for the current from the mobile electrode to the fixed one reads I = -ie(A -A † )/ . Linear response theory gives:

I = e 2 kq |t kq | 2 +∞ -∞ dte -i(ξ k -ξq)t/ × e -iφ(t) e iφ(0) f k (1 -f q ) -e iφ(0) e -iφ(t) f q (1 -f k ) (3)
where f k,q = 1/(e βξ k,q + 1) are the Fermi distributions in the leads, and β is the inverse temperature. In order to include the quantum fluctuation of the electromagnetic fields, the phases in this expression must be regarded as operators whose dynamics is determined by the rest of the Hamiltonian (2). We will restrict to the quadratic part of the Hamiltonian. Within this approximation, the fields φ satisfy Wick theorem and expression (3) can be simplified by making use of the exact relation:

e i φ(t) e -i φ(0) = e J(t)-J(0) (4)
where φ = φ(t) -eV t/ and J(t) = φ(t) φ(0) . The current can then be expressed in terms of the Fourier transform of e J(t)-J(0) :

P (E) = +∞ -∞ dt 2π e J(t)-J(0)+iEt/ . (5)
The function P (E) gives the probability that the environment absorbs a quantity of energy equal to E during the tunneling process. Finally the current in terms of P (E) reads: 34

I = 1 eR T +∞ -∞ EdEP (eV -E) (1 -e -βeV ) 1 -e -βE , (6)
where R T is the tunneling resistance in the absence of the environment [i.e. when P (E) = δ(E)].

In order to find the behavior of the device we need to determine the function P (E) in the presence of the oscillator. Since we will restrict to the quadratic part of the Hamiltonian, we can obtain the quantum phase correlation function J(t) from the classical response function of the phase to a current source at the junction. Specifically we add to the total Hamiltonian (2) the term -I ext (t) φint (t) /e, where φint (t) is the phase in the interacting picture. The linear response to the external field

I ext gives φ(t) = - +∞ -∞ dt ′ χ(t -t ′ )I ext (t ′ ) e (7)
with

χ(t) = i [ φint (0), φint (t)] θ(t).
Using the fluctuation-dissipation theorem one can show that

J(ω) = -2 Im[χ(ω)] 1 -e βω . ( 8
)
Since the Hamiltonian is quadratic, χ(ω) can be calculated from the classical response function (see for instance Ref. 35). The problem is now reduced to the determination of the classical retarded response function. This can be readily obtained by solving the classical equations of motion by Laplace transform.

We determine now the equations of motion for the phase, the position of the oscillator, and the degrees of freedom of the environment. To describe the dissipation we introduce a bath of harmonic oscillators coupled to the mechanical oscillator and a bath of electromagnetic oscillators coupled to the current through the tunnel junction. The coupling of the mechanical and electrical part comes from the charge energy Q 2 /2C(x) of the capacitor. The equations of motion read:

mẍ = -kx - p k p (x -x p ) + 1 2 C ′ (x) φ2 e 2 (9)
C(x) φ = -C ′ (x) ẋ φ + n 1 L n (φ e -φ n -φ) + I ext (t) e (10) m p ẍp = k p (x -x p ) (11)
C n φn = 1 L n (φ e -φ n -φ) . (12)
Here x p and φ n are the degrees of freedom of the environment, φ e (t) = eV t/ is due to the constant external voltage source and we choose the zero of x as the equilibrium position of the oscillator when V = 0. We also define C ′ (x) = dC/dx.

For I ext = 0, this system of equations has a stationary solution given by x = x eq = C ′ (x eq )V 2 /2k, φ = eV / , x p = x eq and φ n = 0. This solution is unique for a given gauge. We thus expand the equations around this solution and consider the quadratic fluctuations of the fields only: x(t) = x eq + x(t), x p (t) = x eq + xp (t), and φ(t) = eV t/ + φ(t). Let us introduce the Laplace transform for the fluctuating fields. For example, for x it reads: x(s) = +∞ 0 dte -st x(t), with Re[s] > 0 and x(t) = +i∞+a -i∞+a (ds/2πi)e st x(s) with a > 0. The system of differential equations then becomes an algebraic system of equations. We are interested in the response function. Thus, we can neglect transient terms. Keeping only linear terms in the fluctuating fields, and solving first for the degrees of freedom of the environment one finds:

s 2 + γ(s)s + ω 2 m x = C ′ V em s φ ( 13
) [Cs + Y (s)] s φ = -C ′ sx eV + I ext (t) e , (14)
where we have introduced the memory function of the oscillator and the impedance of the circuit:

γ(s) = 1 m p k p s s 2 + ω 2 p , Z(s/i) -1 = n 1 L n s s 2 + ω 2 n .
The two coupled equations ( 13) and ( 14) describe the response of the phase φ to the external current I ext flowing through the tunnel junction. Solving for φ and letting s → iω + 0 + we obtain the expression for χ(ω). Using Eq. ( 8), the expression for J(ω) reads:

J(ω) = 4π ω Re[Z t (ω)] R Q 1 1 -e βω (15)
with

Z t (ω) = iωC + Z(ω) -1 + C ′ 2 V 2 iω/m ω 2 m + γ(iω)iω -ω 2 -1 . (16)
The usual result of the Coulomb blockade can be obtained by setting C ′ = 0 in this expression. Then, only the electromagnetic environment would contribute to the response function. The presence of the oscillator simply modifies the response of the system; for γ constant, it would be equivalent to a RLC circuit in parallel with the tunnel junction.

The current voltage characteristics can now be obtained by combining Eqs. ( 5), ( 6), (15), and ( 16). In the following section we consider two representative cases for the circuit impedance and discuss the expected conductances.

III. RESULTS IN SPECIFIC CASES

We consider now two simple cases for the impedance of the circuit: Z = iωL and Z = R. In the first case (Sect. III A), the electromagnetic environment is an LCoscillator. Thus, the system reduces to two coupled oscillators: one mechanical and the other electromagnetic. This example is interesting for his simplicity. It clearly shows how the mechanical and electric part becomes coupled and how this coupling appears in the current voltage characteristics. By contrast, its experimental realization is non-trivial, since it is in general not easy to get rid of parasitic capacitances and resistances. An inductance and a capacitance are present in the environment for the oscillating tunneling junction in the experiment of Ref. 33. The environment forms in this way a radio frequency resonator that is used to improve the band-width of the detection scheme. However, one should note that experiment of Ref. 33 is not in the regime of parameters considered in the present paper.

In the second case (Sect. III B), the environment provides the dissipation and the conditions for the standard Coulomb blockade with the suppression of the current at low voltages. This case is more relevant from the experimental point of view, since high impedance environment can be more easily devised. For electronic transport without mechanical motion an accurate experimental study of the effect of a purely resistive environment has been carried out recently by tuning the environment resistance for a given tunnel junction. [START_REF] Altimiras | [END_REF]

A. Inductive case Z(ω) = iLω

We consider the case of undamped harmonic oscillator. Strictly speaking, the hypothesis of thermal equilibrium for the oscillator and for the electromagnetic modes is not valid. We will assume that an infinitesimal dissipation is present and that the tunneling resistance is sufficiently large to let the time for the system to relax between two tunneling events. The electromagnetic part of the circuit is characterized by the oscillator frequency ω LC = 1/ √ LC. The pure electromagnetic case is discussed in the review paper Ref. 34. Due to the coupling with the mechanical mode two resonant frequencies appear, as can be seen by finding the poles of Z t (ω) [cfr. Eq. ( 16)]:

ω 2 ± = ω 2 m + ω 2 LC + ω 2 I ± (ω 2 LC + ω 2 m + ω 2 I ) 2 -4ω 2 LC ω 2 m 2 .
(17) Here we introduced the "coupling frequency"

ω 2 I = C ′2 V 2 /Cm. The real part of Z t reduces to a sum of four delta functions Re[Z t (ω)] R Q = σ=± ρ σ ω σ 2 [δ(ω -ω σ ) + δ(ω + ω σ )] (18) with ρ ± = (E C / ω ± )|ω 2 ± -ω 2 m |/|ω 2 ± -ω 2 ∓ |.
It is convenient to express the coupling of mechanical and electronic degrees of freedom in terms of a parameter that depends only on the device and does not depend on the voltage bias. The reason is that it is easy to vary experimentally the voltage for a given device. A physically relevant parameter is the variation of the elastic energy E e when the oscillator is displaced by a distance ∆x = E C (C ′ /C)/k in response to the Coulomb force generated by a single electron on the capacitor. With this definition we have

E e = E 2 C (C ′ /C) 2 /2k and ω 2 I = ω 2 m E e (eV ) 2 /E 3 C . ( 19
) For a junction made of a STM tip on top of a suspended carbon nanotube, an electromechanical coupling with the radial breathing mode of the tube was reported 37,38 . In this context, ω m and E c were of the same order of magnitude (around 10meV ). Thus, we estimate the coupling parameter E e /E c ∼ (x 0 /L c ) 2 , where x 0 is the zeropoint motion of the mode (a fraction of Angström) and L c = C/C ′ (a fraction of nm). In the following, we assume a (rather optimistic) ratio of E e /E c 0.1

For small value of E e /E C and ω m = ω LC the new frequencies ω ′ m and ω ′ LC are only weakly modified by the interaction with respect to their bare values:

ω ′2 m /ω 2 m = 1 + ω 2 I /(ω 2 m -ω 2 LC ) (20) ω ′2 LC /ω 2 LC = 1 -ω 2 I /(ω 2 m -ω 2 LC ) (21)
(note that ω ′ m = ω + if ω m > ω LC and ω ′ m = ω -if ω m < ω LC
) and the weights of the poles simplify:

ρ m ≈ E C ω m ω 2 I ω 2 m (ω 2 LC -ω 2 m ) 2 ≪ ρ LC ≈ E C ω LC . ( 22
)
The inequality holds if the bare frequencies ω m and ω LC are of the same order of magnitude. On the other hand, when ω m = ω LC ≫ ω I , we find ω ± = ω m ± ω I /2 and ρ + = ρ -= ρ LC /2. ¿From Eq. ( 15), one can obtain J(t) -J(0) = σ ρ σ j σ (t) with: 34

j σ (t) = coth β ω σ 2 (cos(ω σ t) -1) -i sin(ω σ t) . ( 23
)
We consider now the zero temperature limit. In this case, the function P (E) obtained from Eq. ( 5) has a simple expression in terms of an infinite series:

P (E) = ∞ n,n ′ =0 p n (ρ m )p n ′ (ρ LC )δ(E -n ω ′ m -n ′ ω ′ LC ) (24) where p n (ρ) = e -ρ ρ n /n!.
The differential conductance G = dI/dV can be obtained from the current (6). Neglecting the weak Vdependence of ρ σ and ω σ in the derivative we obtain:

R T G(V ) = nω ′ m + n ′ ω ′ LC <eV p n (ρ m )p n ′ (ρ LC ) . ( 25
)
It leads to a double series of steps spaced by ω m and ω LC . If ρ m = 0, one recovers the usual sequence of steps with a suppression of the conductance at low voltage and R T G = 1 for eV ≫ E C . The interesting situation is when ρ LC is of the order of 1. In this case steps are clearly visible for eV ∼ E C (cf Fig. 2). If ρ m = 0 and ω m < ω LC the mechanical oscillator manifests itself as a series of additional steps of width ω ′ m particularly visible at the beginning of each step of the LC circuit. In practice the picture is slightly more complex due to the dependence on V of the four parameters ρ σ and ω σ . The renormalization of the frequency and the change of ρ introduces a shift of the steps and a weak smooth Vdependence of the conductance plateaux, as shown in Fig. 2. The presence of the mechanical oscillator thus gives a similar result to what has been predicted for a single electron transistor 7 with Frank-Condon steps and a reduction of the current at low voltage. This result is found here in a much simpler structure: a single tunnel junction coupled in series with an external inductance.

B. Ohmic case Z(ω) = R

Let us now consider the pure Ohmic case for the external impedance Z(ω) = R. The first step is again to find the poles of Z t (ω). For γ(ω) constant, three poles are present: one is pure imaginary, iγ ′ RC , and the two other ones are complex, ω ± = ±ω ′ m + iγ ′ /2. In the relevant limit of small coupling ω I ≪ ω m , γ RC and low mechani-cal damping γ ≪ ω m , γ RC we get:

γ ′ RC = γ RC - ω 2 I γ RC γ 2 RC + ω 2 m + . . . , (26)
ω ′ m = ω m + 1 2 ω m ω 2 I γ 2 RC + ω 2 m - γ 2 8ω m + . . . , (27)
γ ′ = γ + γ RC ω 2 I γ 2 RC + ω 2 m + . . . . (28)
The mechanical oscillator acquires a damping even if γ = 0. This intrinsic damping is proportional to the coupling. ¿From here, we drop the bare mechanical damping since it does not induce new effects compared to those already introduced by the presence of R. An important parameter for this device is the dimensionless conductance of the tunnel junction g = R Q /R. For g ≫ 1 the Coulomb blockade is no more observable, since the charge can vary continuously. In the absence of the oscillator the appearance of the Coulomb blockade is signaled by the change of the low-energy behavior of P (E), that goes like E 2/g-1 . 34 Thus for g < 2 the current is suppressed at low voltage V < E C /e. If we consider the interesting case of ω m ∼ E C , then the condition g ≪ 1 implies also that γ RC ≪ ω m . In this limit the real part of Z t is simply the sum of three Lorentzian functions:

Re[Z t (ω)] R Q = E C [ǫL(ω, 2γ ′ RC ) + γ ′ γ RC [L(ω -ω ′ m , γ ′ ) + L(ω + ω ′ m , γ ′ )] (29)
with ǫ = γ ′ RC /γ RC , and L(ω, γ) = (γ/2)/(ω 2 + γ 2 /4)/π. Since the coupling is weak, the damping of the mechanical mode is also very small and the Lorentzian function can be approximated by a Dirac delta-function [L(ω, γ → 0) = δ(ω)]. At zero temperature this gives for the phase correlation function:

J(t) -J(0) = E C π γ RC j o (γ ′ RC t) + ρ m [e -iω ′ m t -1] (30)
where

ρ m = 2EC ωm γ ′ γRC and j o (τ ) = 2 ∞ 0 e ixτ -1
x(x 2 +1) dx is the correlation function of the circuit in the absence of the mechanical oscillator. P (E) thus reads:

P (E) = ∞ n=0 p n (ρ m ) ǫ P o [(E -n ω ′ m )/ǫ] (31)
where

P o (E) = 1 2π +∞ -∞ dte iEt/ +(EC /π γRC)jo(tγRC ) (32)
is the function P (E) in the absence of the oscillator. For the case at hand of small g the function P o (E) is peaked at E ≈ E C with a width ∼ (E C /π) 2g| log g| that vanishes for g → 0. Thus, the resulting P (E) is a sequence of peaks shifted by ω m . Figure 3 shows the dependence of P (E) on E and V as obtained by numerical integration. For eV /E C ≪ E C /E e (Eqs. ( 21) and ( 29)) the mechanical coupling is negligible and one recovers the usual picture of Coulomb blockade, with a single peak at E ≈ E C . For larger values of eV /E C new peaks appear and their positions slowly drift as a function of V towards higher energies. The first peak instead drifts towards smaller energies, due to the ǫ factor in (31). Since γ ′ RC reduces for larger values of V , the width of the peaks also reduces.

In Fig. 4 we show the same function, but for a higher mechanical frequency ω m /E C = 0.5. It shows the appearance of the peaks and their shrinking.

An important case is g → 0, for which P o (E) = δ(E -E C ). This is clearly the most favorable case to observe the quantum fluctuations of the oscillator directly in the I-V characteristics. Indeed, even in the weak electromechanical coupling the sharpness of the peaks gives rise to steps in G(V ), like in the case of a purely electromagnetic resonator. The conductance has thus the form (again neglecting the small voltage dependence of ρ m and ω m in performing the derivative):

R T G(V ) = nω ′ m <eV -ǫEC p n (ρ m ) ( 33
)
the result is plotted in the insert of Fig. 5 for E e /E C = 0.1 and g = 10 -4 . For the general case the conductance can be obtained numerically by integrating P (E). It is shown for three different values of the ratio E e /E C in Fig. 5. The presence of the mechanical coupling introduces oscillations of the conductance on the scale eV ∼ ω m that are precursors of the step behavior in the extreme g = 0 case. The small increase of the conductance at low voltage is the consequence of the renormalisation of E C due to the electromechanical coupling (the ǫ factor).

A last comment is in order for the resistive case. Even 4) (dot-dashed line) and for vanishing coupling (full line). Inset: G(V ) for g = 10 -4 from the analytical expression (33) for ωm = 0.1EC , Ee = 0 (full line) and Ee = 0.1EC (dashed line). The conductance is much more sensitive to the electromechanical coupling in presence of a high impedance environment.

if this system is described by a very similar electric circuit to that used for a single electron transistor without gate, some differences should not be overlooked. In the single electron transistor the two resistance are tunnel junctions, while here one is a tunnel junction and one a ohmic resistor. A tunnel junction is substantially different from a ohmic resistor, since in the first case the charge can be transferred only by a sudden tunneling event, while in the second the charge slowly leaks through the resistor. The charge quantization is thus not enforced by the Ohmic resistor, and this has consequences also on the time dependence of the force acting on the mobile part. For these reasons our results complement those obtained by others authors for suspended single electron transistors. 3,4,5,7,8,9,11,12,13

C. Range of validity

Let us return to the comparison of position dependence of the tunneling and capacitive terms. Using the explicit expression (28) for γ ′ , we find that the condition to neglect the effect of the tunneling force (case γ RC ≪ ω m ) is:

1 ≫ γ T γ ′ = x o λ T 2 ω m R Q E 3 C γ RC R T E e e 2 V 2 . ( 34
)
This gives the explicit condition on the voltage in order to neglect the position dependent part of the Hamiltonian. Formally, it is clear that it can be always fulfilled for R T sufficiently large.

A second possible source of deviation from the behavior we found comes from the cubic term of of the interaction energy. We neglected the term ( φ) 2 with respect to eV φ/ . This is justified when ( φ) 2 1/2 ≪ eV / . We thus obtain the condition of validity of the quadratic theory at zero temperature:

eV E C ≫ g| ln g| . (35)
IV. JOSEPHSON JUNCTION Observing a Josephson supercurrent through a STM tunnel junction has been shown to be feasible. 39 Particularly appealing is the possibility of studying an oscillating Josephson junction. The effect of the position dependence of the electronic tunneling amplitude on the dynamics of the Josephson junction has been considered recently. 40 However the effect of the position dependence of the capacitance (Coulomb interaction) has not been explored so far. The theory of Coulomb blockade for normal metallic electrodes presented above can be generalized with little change to the superconducting case. The main difference is that the phase difference is now a truly observable quantity and it is conjugated with the number of Cooper pairs that cross the junction.

The Josephson junction is characterized by a critical current I c that defines the energy scale E J = I c /2e. For E J /E C ≪ 1 a tunneling approach can be used to obtain the dissipative current of the device. We remind the main steps of the procedure. [START_REF] Ingold | Single Charge Tunneling[END_REF] The current operator is I S = I c sin(2φ) while the Josephson part of the Hamiltonian reads H J = E J cos(2φ) (we keep the same definition of the phase, the factor two intervenes in the usual Josephson relation).

Calculating the current in the linear response limit we get As long as ωm remains smaller than the width of the peak, the coupling manifests itself in a broadening of the Coulomb peak. For ωm/EC = 0.5, contributions of the mechanical oscillator clearly appear.

I S = I c E J 4 +∞ -∞ dt ′ e 2iφ(t) , e -2iφ(t ′ ) e 2ieV t . (36
One can thus express the current again in terms of the Fourier transform of the phase correlator:

I(V ) = πE 2 J e [P s (2eV ) -P s (-2eV )] (37)
where the definition of P s (E) differs from the definition of P (E) simply by a factor of 4 in front of J(t):

P s (E) = +∞ -∞
dte 4J(t)-4J(0)+iEt/ .

The function J(t) depends on the mechanical response function as shown in Section II. The results of Sects. II and III can thus be applied to the superconducting case to obtain P s (E) and hence I(V ). In the Figure 6 we show the expected current-voltage characteristics for an oscillating Josephson junction in the pure resistive case of Sect. III. The superconducting device is much more sensitive to the presence of the mechanical oscillator, since the voltage dependence is directly related to P s (E), and not to its integral. One can then observe directly the Coulomb blockade peaks in the I-V characteristics.

V. CONCLUSION

In this work, we considered the current-voltage characteristics of an oscillating tunnel junction. At low tunnel coupling and moderately low voltage bias, we found that the coupling between the electronic and mechanic degrees of freedom mainly originates from the position dependence of the junction capacitance. That is, the electrons crossing the tunnel junction couple to an electromagnetic environment formed by the rest of the circuit which is itself coupled to the mechanical oscillator. In our study, we included both the electromagnetic and mechanical environments in the Coulomb blockade theory. We found that the mechanical properties of the junction result in an additional mechanism for the current suppression that is similar to the Frank-Condon blockade predicted in an oscillating single-electron transistor. As in the conventional Coulomb blockade, similar, though more pronounced, effects would also appear in a voltage-biased oscillating Josephson junction.

At larger tunneling rates, two new ingredients should be added to our study. First, the electromagnetic and mechanical degrees of freedom would not have time to relax to equilibrium between two tunneling events. Second, the position dependence of the tunneling matrix element may not give anymore negligible effects compared to those arising from the position dependence of the capacitance. Considering both, it may be interesting to investigate whether a shuttling instability takes place in a single oscillating tunnel junction.

FIG. 2 :

 2 FIG.2: Differential conductance of the device as a function of the bias voltage for the case Z = iωL. The bare oscillating frequencies are ωm = 0.5EC and ωLC = EC. The electromechanical coupling is Ee = 0 (full line) and Ee = 0.1EC (dashed line) (in this case, signatures of mechanical coupling are not visible for ωm below 0.5EC ).

FIG. 3 :

 3 FIG.3: P (E) is displayed for different values of V . We take g = 0.015, ωm = 0.1EC and Ee = 0.1EC .

FIG. 4 :FIG. 5 :

 45 FIG. 4: P (E) is obtained numerically for different values of V for g = 0.015, ωm = 0.5EC and Ee = 0.1EC .

  FIG.6: I(V ) is plotted for g = 0.065. We distinguish the case Ee = 0 (full line), Ee = 0.1EC with ωm = 0.1EC (dashed line) and Ee = 0.1EC with ωm = 0.5EC (dot-dashed line). As long as ωm remains smaller than the width of the peak, the coupling manifests itself in a broadening of the Coulomb peak. For ωm/EC = 0.5, contributions of the mechanical oscillator clearly appear.

Acknowledgements

We acknowledge useful discussions with D. Bagrets. This work has been supported by the French Agence Nationale de la Recherche under contract ANR-06-JCJC-036 NEMESIS. F. P. thanks A. Buzdin and his group for hospitality at the Centre de Physique Moleculaire Optique et Hertzienne of Bordeaux (France) where part of this work has been completed.


Loïc Denis
email: loic.denis@univ-st-etienne.fr

Thierry Fournel

Corinne Fournier

Dominique Jeulin

Reconstruction of the rose of directions from a digital micro-hologram of fibers

Digital holography makes it possible to acquire quickly the interference patterns of objects spread in a volume. The digital processing of the fringes is still too slow to achieve on line analysis of the holograms. We describe a new approach to obtain information on the direction of illuminated objects. The key idea is to avoid reconstruction of the volume followed by classical three-dimensional image processing. The hologram is processed using a global analysis based on autocorrelation. A fundamental property of diffraction patterns leads to an estimate of the mean geometric-covariogram (MGC) of the objects projections.

The rose of directions is connected with the MGC through an inverse problem.

In the general case, only the 2D rose of the object projections can be reconstructed. The further assumption of unique-size objects gives access with the knowledge of this size to the 3D direction information. An iterative scheme is suggested to reconstruct the 3D rose in this special case. Results are provided on holograms of paper fibers.

Introduction

Non-intrusive on-line control is of great importance in most production processes. In the case of 3D moving objects, designing this control is challenging. To track the directions of small objects spread in a volume one requires both a fast 3D imaging facility and an efficient image processing scheme.

Holography is a technique that allows a description of a three-dimensional distribution of objects to be stored on a two-dimensional detector. It ensures non intrusive control of still or moving objects. The recording of the hologram on a digital camera suppresses the wet chemical processing step and makes it possible to acquire and analyze volume objects within a short time, the limiting factor being the numerical processing speed. It gives rise to a whole range of potential algorithms for the processing of the holograms. Digital holography has been widely studied these last years (see recent contributions published in (Poon et al., 2006)) and applied to many different fields such as microscopy, particle image velocimetry and deformation analysis. It plays an important role for the study of small objects and is used in fluid mechanics for 3D flow visualization.

Since fast acquisition of 3D objects is possible with digital holography, we focus in this article on the hologram processing. We describe a method to study the direction of fibers randomly distributed in a volume from a single hologram. We follow an original stereological approach that differs from classical analysis of the volume reconstructed from the hologram. This work follows a sizing method based on the hologram-autocorrelation analysis [START_REF] Denis | Direct extraction of mean particle size from a digital hologram[END_REF].

To quantify anisotropy of the fiber distribution, a statistical representation is often used: the rose of directions R [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. It is the distribution of directions weightened by the mean fiber-length in each direction. If the fiber process is isotropic, then R is the uniform distribution. By contrast, the fact that R is uniform does not imply that the underlying fiber process is isotropic [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Stoyan | Anisotropy analysis for particle systems[END_REF][START_REF] Kiderlen | Algorithms to estimate the rose of directions of a spatial fibre system[END_REF]. A non-uniform rose of directions clearly indicates anisotropy. In particular, preferred directions of a fiber system can be detected [START_REF] Kiderlen | Algorithms to estimate the rose of directions of a spatial fibre system[END_REF]. The estimation of R requires either direct measurement of spatial directions or an indirect approach. In the latter case, R is derived from another quantity carrying the same directional information while being directly estimated from the data. The mean geometric covariogram (MGC) [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], defined as the autocorrelation of an object, averaged over the objects population, is such a quantity. Another commonly used quantity to derive R is the rose of intercepts γ, defined as the diametral variation in a given direction [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF].

After a short introduction to digital holography (section 2), we derive the rose of directions of the fibers projections from their mean geometric covariogram computed from the hologram. The mean geometric covariogram of the objects projections is estimated directly from the hologram (section 3). In the general case, only the 2D rose of directions can be derived from the estimated MGC (section 4). With a further assumption one can reconstruct the 3D rose of the objects. An iterative reconstruction algorithm is suggested in this case and the method is illustrated on simulated and experimental holograms (section 5).

In-line digital holography

Let's give a preliminary description of in-line digital holography to introduce the basics and the notations we will follow in the next sections.

The in-line holography recording setup, also named "Gabor holography" after its inventor, involves no imaging lens. The intensity of the incident coherent beam, diffracted when crossing the objects and free propagated beyond is recorded on the sensor (a digital camera) (see Figure 1). This setup, used to record holograms of small opaque objects such as microparticles or fibers, is robust with respect to vibrations thanks to the unique beam used.

Increasing the concentration of objects yields a degradation of the hologram quality due to speckle noise. This technique therefore is restricted to collections of objects whose projected surface on the hologram plane remains under a limit (1% of the hologram surface) [START_REF] Royer | An application of high-speed microholography: the metrology of fogs[END_REF].

The recorded hologram can be explained as the interference pattern between two beams: the reference beam and the objects beam. Unlike with two-beams setups, in Gabor holog-raphy there is a single beam which can be interpreted past the objects as the superposition of two waves: a reference wave identical to the incident wave and a object wave emerging from the objects [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. In all-optic holography, the reconstruction is led by reilluminating the hologram with the reference beam. The hologram then diffracts the beam and gives rise to the real and virtual images. Diffraction of a plane wave by a plane object of given transmittance t can be modeled under Fresnel approximation as a spatial convolution [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. The complex amplitude A z of the diffracted wave in a plane located at a distance z related to the objects is given by:

A z (x, y) = t * h z (x, y) (1)
with the following kernel:

h z (x, y) = 1/(iλz) exp iπ(x 2 + y 2 )/(λz)
where i stands for the imaginary unit, λ is the wavelength and * denotes the two-dimensional convolution. The complex amplitude A z is called the Fresnel transform of transmittance function t with parameter λz (scale parameter). Since diffraction governs both the hologram formation and its reconstruction, Fresnel transform is a fundamental tool for both hologram modeling and digital reconstruction.

Fresnel transform is a totally redundant scaling transform [START_REF] Onural | Diffraction from a wavelet point of view[END_REF]. By changing the scale parameter, images can be reconstructed at various depths z. The redundancy of the transform results in the presence of out-of-focus objects in the reconstructed images. These out-of-focus objects appear as low contrast concentric fringes. The numerical transform is performed using Fast Fourier transforms and a digital focusing is obtained by modifying the parameter z of the convolution kernel.

Estimating the mean geometric covariogram of the objects projections

We derive in this section an estimate of the mean geometric covariogram (MGC) of the objects projections from the expression of the hologram intensity.

In the context of in-line holography of small opaque objects, the hologram intensity is related with the objects through a diffraction phenomenon. If we consider the incoming wave to remain quasi-plane for each object, which is justified for holograms of low-concentration collections of objects, the resulting diffracted wave in the recording plane is [START_REF] Denis | Direct extraction of mean particle size from a digital hologram[END_REF]:

A ≈ 1 - N j=1 ϑ j * h z j * δ x j ,y j
where z j is the distance between object j and the recording plane, N is the number of objects, ϑ j denotes the aperture of object j, δ x j ,y j stands for Dirac's distribution centered in (x j , y j ).

In this expression, the first term can be interpreted as the reference wave (with normalized energy and null phase in the hologram plane) ; the other term is a sum of the contribution of each object. These contributions are complex valued and are summed up.

The recorded signal in the hologram plane corresponds to the intensity I H of the wave, related to the complex amplitude via the following relation:

I H = |A zr | 2 .
The hologram intensity can therefore be expanded as follows:

I H = 1 -2 N j=1 ϑ j * ℜ(h z j ) * δ x j ,y j + N j=1 N k=1 ϑ j * h z j * δ x j ,y j × ϑ k * h -z k * δ x k ,y k , (2)
where ℜ denotes the real part of a complex number. The first term is constant and corresponds to the reference wave, the second term is the sum of all diffraction patterns, and the last term expresses both inter-object interferences (j = k) and the non-linear component of a diffraction pattern (j = k) . Under the hypothesis of both small objects and lowconcentration holograms, the third term can be neglected. This gives the following linear hologram-formation approximation:

I H ≈ 1 -2 N j=1 ϑ j * ℜ(h z j ) * δ x j ,y j . (3)
When considering objects with a random distribution in space, implanted on points of a 3D homogeneous Poisson point process with intensity c (i.e. c is the average number of points per unit volume), the hologram can then be described in the context of random set models as a dilution random function (Serra, 1968;[START_REF] Jeulin | Multivariate random image models[END_REF]. The primary functions ϑ j * ℜ(h z j ) depend on the distance z j to the hologram and on the object aperture ϑ j (defined as the indicator function of the projection of the j-th object). They are implanted on points of a 2D

Poisson point process resulting from the projection with intensity c • e (e being the thickness of the volume) and whose realization is the set of (x j , y j ) j objects locations.

The

ℜ(h z ) * ℜ(h -z ) = 1 2 [δ 0,0 +ℜ(h 2z )].
As Fresnel functions are symmetrical and ℜ(h -z ) = ℜ(h z ), the last equation can be rewritten in terms of cross-correlation:

ℜ(h z ) ⋆ ℜ(h z ) = 1 2 [δ 0,0 + ℜ(h 2z )
] (with ⋆ denoting the cross-correlation operation, connected to the convolution through the relation f ⋆ g = f * ǧ with f and g any distribution and ǧ the symmetric of g with respect to the origin). The aperture of object j is nearly recovered when the jth diffraction pattern ϑ j * ℜ(h z j ) * δ x j ,y j is cross-correlated with ℜ(h z j ). This is the basic principle classically used in digital in-line holography. The distance z j has to be known in order to get an in-focus image of the object.

Instead of using ℜ(h z j ) to perform the cross-correlation, the diffraction pattern can be autocorrelated. The (x j , y j ) location is then lost and the remaining terms are: 1 2 ϑ j ⋆ϑ j + 1 2 ϑ j ⋆ϑ j * ℜ(h 2z j ). The autocorrelation of the centered hologram intensity ĨH (i.e. its autocovariance) can be expanded into three terms [START_REF] Denis | Direct extraction of mean particle size from a digital hologram[END_REF]:

ĨH ⋆ ĨH ∝ 2 N j=1 ϑ j ⋆ ϑ j + 2 N j=1 ϑ j ⋆ ϑ j * ℜ(h 2z j ) + 4 N j=1 N k=1 k =j ϑ j ⋆ ϑ k * ℜ(h z j ) * ℜ(h z k ) * δ x j -x k ,y j -y k . (4)
The first term corresponds to the mean geometrical covariogram of the apertures (up to a proportionality factor accounting for the number of objects and experimental parameters such as the camera gain and the laser power); the second to a virtual image of the geometric covariogram; the third to slightly out-of-focus and largely out-of-focus cross-correlation terms. It has been shown in the quantitative study carried out in reference [START_REF] Denis | Direct extraction of mean particle size from a digital hologram[END_REF] that the MGC dominates the other terms on most of its range. Thus, in the case of holograms of low-concentration volumes and under the assumption of a homogeneous Poissonian object process, the hologram autocovariance leads to a good estimate of the MGC.

Reconstruction of the rose of directions

We have shown in the previous section that the mean geometric-covariogram (MGC) can be estimated easily from a hologram. The MGC carries information about the direction, size and shape of the objects projections. However the contributions of shape(s) and size(s) mix with that of directions. The resulting MGC can therefore be difficult to interpret quantitatively.

We describe in this section how the MGC can be processed to get a reconstruction of a rose of directions.

In the general case, only a projection of the rose of directions of the objects can be reached from the estimated MGC. We describe in section 4.A how the 2D rose of directions can be derived form the MGC. We then describe in section 4.B how the MGC of the objects projections can be connected with the theoretical geometric-covariogram of a parametrized model. The MGC can easily be computed for given shapes, sizes and directions of the objects.

The reconstruction problem of the 3D rose of directions is an inverse problem. Due to the projection operation that connects the objects to the MGC, additional information is required to reconstruct the 3D rose of directions. We suggest a reconstruction algorithm in the case of short fibers under the assumption of constant and known object size. This description is of more general scope and the derivation can easily be carried out for other types of shapes.

4.A. Derivation and inversion of the rose of intercepts

We review in this paragraph the first possible approach for MGC processing: the derivation of the rose of intercepts.

The rose of intercepts γ(α) defines the diametral variation in the direction orthogonal to α. It can be easily measured on an image of the objects. In our case, objects are not directly accessible and the rose of intercepts must be derived from the MGC expressed in polar coordinates k(r, α) using the following relation [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]:

γ(α) = - ∂ k(r, α) ∂r r=0
(5)

The rose of directions can then be reconstructed from the rose of intercepts using one of the inversion schemes described in the literature (see reference [START_REF] Kiderlen | Algorithms to estimate the rose of directions of a spatial fibre system[END_REF] for an overview of available methods).

This approach is however limited by the difficulty to get a robust estimation of the zero derivative of the MGC as the presence of uncorrelated noise in the hologram tends to interfere with it. Furthermore, the MGC is underused which is regrettable as the underlying inverse problem is ill-posed. To circumvent these drawbacks, a more general approach making use of the whole MGC is described in the following paragraphs. The objects under study to illustrate the reasoning are micro-fibers.

4.B. Formulation of the general inverse problem

We derive in this section the connection between the MGC, estimated directly from the hologram, and the distribution of directions of the underlying fiber population.

Let a short fiber be defined by its direction (θ, φ) (given in spherical coordinates), its length L and its thickness t (see Figure 2). Its projection along the z axis can be approximated by a rectangle of height h = t and width w = L sin φ + t cos φ.

The geometric covariogram K(r, α) of a w × h rectangle directed along θ is the common area between the rectangle and its translate by a vector of length r and direction α (see Figure 3):

         K t,L,θ,φ (r, α) = (w -r| cos(α -θ)|) (h -r| sin(α -θ)|) if w -r| cos(α -θ)| ≥ 0 and h -r| sin(α -θ)| ≥ 0; K t,L,θ,φ (r, α) = 0 if not. (6)
A graphical representation of the geometric covariogram is given in Figure 4.

As the contribution of each fiber is added to form the mean geometric covariogram (see section 3), the overall MGC k is connected to the density f (t, L, θ, φ) of sizes and directions through the integral relation:

k(r, α) data = K t,L,θ,φ (r, α) model f (t, L, θ, φ) unknown dtdLdθdφ + n(r, α) noise . (7)
This equation links together the data and the unknown and underlines the inverse problem nature of the reconstruction. The geometric covariogram K t,L,θ,φ (r, α) represents the model of the shape/size/direction dependency in the MGC and is known analytically. Since the MGC is obtained from experimental data and with an approximation, we introduce a noise term that represents the difference between the modelization and the actual data.

Equation 7 expresses both a projection operation and an integration over all possible objects. The reconstruction of the density f requires to invert both operations. The projection can not be inverted without prior information. Therefore, the scheme to recover the 3D density information described in the following subsections rely on necessary assumptions.

4.C. Discrete problem

We describe here how equation 7 can be sampled and we suggest an iterative algorithm to reconstruct the unknown distribution f under additional assumptions.

The number of parameters must first be reduced in order to allow the reconstruction. Two problems of interest can be studied:

(P1) reconstruct the 3D direction distribution under the assumption that t and L are known and constant;

(P2) reconstruct the 2D direction and projected length distribution of fibers under the assumption that t is known and constant (only the projection of the directions is studied).

These two problems are of equivalent difficulty since the projected length and the out-ofplane direction play similar roles in the expression of the geometric covariogram K. We therefore describe the reconstruction of f (θ, φ) in the case of fibers of constant length. The results obtained can directly be adapted to problem P2.

The MGC is estimated on a disc that includes the bounded support of the MGC. The sampling used is Cartesian as the estimation is based on Fourier transforms. The noise n(r, α) associated with the data is not stationary as it is derived from an autocorrelation.

To take into account the difference in noise level between the center of the autocorrelation and farther regions, we introduce a weighting matrix W that reduces the influence of data points with low confidence due to the noise level. W is defined in the general case as the inverse data-covariance matrix W = cov( k) -1 .

The unknown distribution f is computed at chosen points (θ, φ). To compare easily the reconstructed density in two different directions one has to use a uniform sampling in θ. A uniform sampling in spherical coordinates would lead to an over-sampling effect for directions close to φ = 0. This region however is badly defined as the directions are nearly collinear to the projection direction. A form of regularization can be introduced by using a looser sampling than spherical coordinates. Since the direction discrimination is connected with the difference between the projections, the sampling can be defined uniformly in the projected space (i.e. in a disc). For this purpose, we chose a sampling based on a spiral distribution inspired from the paper [START_REF] Rakhmanov | Minimal discrete energy on the sphere[END_REF].

The discretized form of equation 7 can be written as a matrix product:

            k(r 1 , α 1 ) . . . k(r n , α 1 ) . . . k(r n , α m )             =             K θ 1 ,φ 1 (r 1 , α 1 ) • • • K θp,φq (r 1 , α 1 ) . . . . . . K θ 1 ,φ 1 (r n , α 1 ) • • • K θp,φq (r n , α 1 ) . . . . . . K θ 1 ,φ 1 (r n , α m ) • • • K θp,φq (r n , α m )                         f (θ 1 , φ 1 ) . . . f (θ p , φ 1 ) . . . f (θ p , φ q )             +             n(r 1 , α 1 ) . . . n(r n , α 1 ) . . . n(r n , α m )             k = K f + n , (8)
where (r k , α l ) represents a sampling point in the autocorrelation space, (θ k , φ l ) a sampling point in the directions space (S 2 ) and bold faced symbols k, K, f and n are compact notations for their vector and matrix counterparts.

4.D. Reconstruction algorithm

There exists a wide range of reconstruction algorithms to solve inverse problems. The choice of a specific approach requires to analyze the characteristics of the problem and especially to identify general properties common to all possible solutions: a priori knowledge. In the case of the rose of directions of fibers, a reconstruction of the rose is of special interest only if there exists one or a few significant directions. The distribution of directions f then is close to zero for all directions but that of the fibers. This "mainly zero" characteristic together with the positivity constraint (f (θ, φ) ≥ 0 for all (θ, φ)) lead us to use a constrained maximum likelihood approach: the solution f (ML) is that which maximizes the likelihood under the positivity constraint:

f (ML) = arg max f Pr( k|f ) = arg max f φ ML (f ), subject to ∀ (θ, φ), f (θ, φ) ≥ 0, (9)
with φ ML (f ) the log-likelihood:

φ ML (f ) ∝ log[Pr( k|f )]+constant.
Well known and easy to implement iterative algorithms have been proposed to solve such problems. They are referred to as multiplicative algorithms [START_REF] Lantéri | Penalized maximum likelihood image restoration with positivity constraints: multiplicative algorithms[END_REF] and differ depending on the noise nature. In the case of Poissonian noise, two identical algorithms have been introduced: Richardson-Lucy Algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF] and the expectation-maximization procedure [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF]. In the case of Gaussian noise, the Image Space Reconstruction Algorithm (ISRA) builds a solution through the following recursion formula (Daube-Witherspoon & Muehllehner, 1986):

f (k+1) j = f (k) j K T W k j (K T WKf (k) ) j , (10)
where the subscript • j denotes element j of the vector, and the superscript • T is used for the transpose operator. As long as the initial value f (0) is positive, the reconstructed densities f (k) remain positive (or null) as the iterative definition of equation 10 involves only positive reals.

Let us notice that the mean geometric covariogram is known up to a proportionality factor, f is therefore retrieved up to that factor and the actual density is obtained by normalization by j f j .

W clearly appears as a weighting factor that reduces the influence of data values with high noise variance. Even if the Gaussian model does not in our case appear as a natural model, this weighting capability (i.e. non-stationarity) is of great interest and overcomes the validity considerations on the Gaussian hypothesis. The approximations described in section 3 indeed lead to a bias that depends on (r, α) parameters. To take into account this bias, one can weight the data to give less influence to data points subject to higher bias.

Application to holograms of paper fibers

The model of short fibers introduced in the previous section can be applied to different kind of objects. The main constraint is that the fibers are required to be straight which is often true for short enough fibers (less than a millimeter). We worked on paper fibers as the on-line study of fiber directions in the very first step of paper making can have significant impact for paper industry.

The principles presented in section 4 are illustrated here on simulated and experimental data. The scheme of the proposed method to reconstruct a distribution of 3D directions in the case of fibers of constant size is depicted on Figure 5. A hologram is processed in two steps: first, the MGC is estimated as described in section 3; then the MGC is inverted to reconstruct the distribution of directions, based on section 4. Results are given on simulated data in the following subsection. In the case of fibers with variable lengths, the derivation of the rose of intercepts is illustrated on an experimental hologram in section 5.B.

5.A. Numerical simulation

The diffraction approximation (equation 1) makes it possible to simulate holograms. Our method can then be confirmed and the precision of the reconstruction assessed by comparing the retrieved distribution with the underlying exact distribution of the random process used to generate numerically the hologram.

Fiber population A set of 100 fibers was generated according to a random uniform distribution in a 7mm × 7mm × 7mm volume (fig. 6). The fibers are oriented along one of two mean directions: (θ 1 , φ 1 ) = (57 • , 30 • ) and (θ 2 , φ 2 ) = (18 • , 83 • ). Two thirds of the fibers have a direction close to (θ 1 , φ 1 ), with a standard deviation of 4 • ; the last third is directed following to (θ 2 , φ 2 ) with a lower standard deviation (1.5 • ). The fibers are 20µm thick and 500µm long. Reconstruction of the distribution of 3D directions The MGC can then be inverted using equation 10. Figure 9 displays the reconstructed distribution after a thousand iterations. A comparison of the retrieved distribution with the exact distribution (Figure 6(b)) shows good agreement. The whole reconstruction process takes about 30 seconds using Matlab R language on a Pentium 4 3GHz PC. The convergence speed can be improved by using a C language implementation in combination with techniques such as the non-linear conjugate gradients or the limited-memory variable metric methods described in [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF].

5.B. Experimental hologram of fibers

We illustrate the general case of fibers with varying sizes and directions in this paragraph.

Objects The objects under study are short paper fibers spread on a glass plate with a brush. The choice of a planar distribution of the fibers was caused by the ease to depict the actual distribution. A reconstructed plane can be computed to show an image of the fibers, for pedagogical reasons.

Hologram The hologram is recorded at 50mm from the glass plate with a 632nm HeNe laser. The camera used has a 1024 × 1280 pixel array with 6.7µm pixel width. The hologram is displayed on Figure 10.

Estimated rose of intercepts

The MGC is estimated from the hologram using the autocovariance technique described in section 3. The derivative at the origin of the MGC is then estimated using a mean square approach [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] to get the rose of intercepts as pointed out in paragraph 4.A. The resulting rose of intercepts is shown on Figure 11. Its shape is characteristic of an anisotropic image. This can be confirmed with the reconstructed image computed in Figure 12.

The 2D rose of directions of the fibers projections can be computed by inverting the rose of intercepts (see section 4).

Conclusion

We have described and illustrated in this paper an original stereological approach for inline digital holograms. The use of an imaging technique capable of fast 3D acquisition in association with direct processing of the data leads to a system suitable for on-line control.

In-line digital holography gives access to moving objects in a deep field (edges a few centimeters long). The hologram can be approximated as a sum of diffraction patterns corresponding to each object. The hologram auto-covariance then appears to lead to a good estimate of the MGC of the projections of the objects. The MGC carries interesting information on the shape of the projection of the objects. In the case of anisotropic objects such as fibers, the MGC can be used to estimate the main direction of a set of objects, as classically done with conventional images. Our MGC derivation from a hologram extends this method to "diffractive projections" of the objects.

In order to get quantitative results, the MGC requires further processing. The reconstruction of the distribution of directions of the objects is an inverse problem that requires additional knowledge to be inverted. We have suggested a scheme to reconstruct the 3D rose of directions of fibers under the assumption of known and unique size. If no assumption can be done on the fibers sizes or directions, only the 2D rose of the projected directions can be derived through the computation of the rose of intersection from the MGC.

We illustrated both cases on simulated or real data. We focused our attention on the extraction of the direction information rather than on an exhaustive study of the performance of our inversion technique. Good results are obtained with our inversion method when there are only a few main directions. Further work should focus on both the robustness of the inversion algorithm and on further improving its convergence speed.

  reconstruction of the volume recorded in a hologram is based on the diffraction of a plane wave by the hologram. The digital processing of a hologram is traditionally carried out via convolution with Fresnel functions. An in-focus image of an object is obtained when the parameter of the Fresnel function h z matches the object-hologram distance. This is made possible thanks to the inversion property of Fresnel function: h z * h -z = δ 0,0 . The physical interpretation of this property is the well-known principle of ray reversibility. The real part of Fresnel functions satisfies a weakened version:

Figure 6

 6 Figure 6(a) displays a 3D view of the fibers. The distribution of directions is shown on Fig. 6(b): directions are represented as vertices on a discretization of the unit sphere.

Fig. 1 .Fig. 3 .Fig. 6 .

 136 Fig. 1. Optical setup used in digital in-line holography

Fig. 10 .

 10 Fig. 10. Experimental hologram of paper fibers spread on a glass plate.

Fig. 12 .

 12 Fig. 12. Reconstructed plane of the fibers showing their varying sizes and directions.

Acknowledgements

This research was supported by a grant from Rhône-Alpes region.