Imen Nasr
email: nasr.imen@supcom.rnu.tn

Leïla Najjar Atallah

Benoît Geller
email: benoit.geller@ensta-paristech.fr

Sofiane Cherif
email: sofiane.cherif@supcom.rnu.tn

CRB Derivation and New Code-Aided Timing Recovery Technique for QAM Modulated Signals

In this paper, we propose a maximum likelihood based Code-Aided (CA) timing recovery algorithm for square-QAM modulated signals. We also theoretically derive the analytical expression of the CA Cramer-Rao Bound for time delay estimation. Our simulations show that the proposed CA approach realizes a performance equivalent to the Data-Aided (DA) approach over a large interval of signal to noise ratio (SNR) values.

I. INTRODUCTION

High spectrum efficiency modulations are required in modern communication systems in order to transmit high data rate information over a limited bandwidth [START_REF] Zhou | Crosslayer rate allocation for multimedia applications in pervasive computing environment[END_REF]. This is why quadrature amplitude modulation (QAM) is commonly used. This modulation requires fine synchronization which is at the origin of the development of many time delay estimators. Some DA and non-data-aided (NDA) time delay estimation techniques have been employed in real systems for decades such as [START_REF] Mueller | Timing recovery in digital synchronous data receivers[END_REF] and [START_REF] Gardner | A bpsk/qpsk timing-error detector for sampled receivers[END_REF]. In DA estimation, reference signals are introduced in the transmitted frame to assist the synchronization process. Even if this technique performs well, it reduces the whole system throughput; this is avoided in the NDA mode where only useful information data, unknown to the receiver, is transmitted. This operating mode optimizes spectral efficiency, but NDA synchronization has low performance compared to the DA mode.

Nevertheless, with the development of channel coding techniques [START_REF] Gallager | Low-density parity-check codes[END_REF]- [START_REF] Vanstraceele | A low complexity block turbo decoder architecture -[transactions letters[END_REF], iterative decoding has recently been extended to other receiver functions, especially time synchronization. Observing that turbo receivers usually operate at low SNR values, many researches have focused on CA synchronization algorithms [START_REF] Berrou | Turbo-Information processing: algorithms, implementations & applications[END_REF] which consist in using the available soft information provided by the turbo detector to assist the synchronizer. This idea has already been applied in a number of contributions. For instance, it was applied for phase estimation in turbo decoding receivers [START_REF] Lottici | Carrier phase recovery for turbo-coded linear modulations[END_REF]. It was also exploited in [START_REF] Noels | Turbo synchronization: an em algorithm interpretation[END_REF] with an expectation maximization algorithm in the Maximum Likelihood (ML) synchronization framework. Moreover, the timing synchronizer and the channel decoder can improve each other progressively by exchanging information using iterative techniques [START_REF] Zhang | APPA symbol timing recovery scheme for turbo codes[END_REF]. In [START_REF] Nayak | Joint timing recovery and turbo equalization for coded partial response channels[END_REF] the authors have proposed to introduce the timing recovery inside the turbo equalizer to jointly perform these tasks. In [START_REF] Barry | Iterative timing recovery[END_REF], an iterative timing recovery technique where the soft information from the decoder is fed into the Mueller and Muller (M&M) timing error detector (TED) [START_REF] Mueller | Timing recovery in digital synchronous data receivers[END_REF] is jointly performed with turbo equalization. A soft maximum likelihood timing recovery technique has also been proposed in [START_REF] Nasr | A soft maximum likelihood technique for time delay recovery[END_REF] where a priori information on the transmitted symbols is used at the demapper and is injected to the M&M TED.

The variance of the timing error is usually compared to a theoretical bound corresponding to the optimum performance. The Cramer-Rao Bound (CRB) is one of the most commonly used lower bounds on variance relative to the class of unbiased estimators and is the easiest to derive [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. Closed-form expressions of the CRB for CA synchronization parameter estimation were specifically derived for carrier phase and frequency estimation. In [START_REF] Yang | Bayesian and hybrid cramer-rao bounds for the carrier recovery under dynamic phase uncertain channels[END_REF]- [START_REF] Bay | Analytic and asymptotic analysis of bayesian cramer-rao bound for dynamical phase offset estimation[END_REF], the authors derived the Bayesian and hybrid Cramer-Rao bounds (BCRB and HCRB) for the CA, the DA, and the NDA dynamical phase estimation of QAM modulated signals and theoretically showed the possible improvement using a CA technique. For timing recovery, only closed-form expressions of the DA CRB [START_REF] Meyr | Wiley Series in Telecommunications and Signal Processing[END_REF] and the NDA CRB for linearly modulated signal [START_REF] Masmoudi | Closed-form expressions for the exact cramer-rao bounds of timing recovery estimators from bpsk, msk and square-qam transmissions[END_REF] have been derived. The modified CRB (MCRB) is easier to derive and averages the bound calculus over nuisance parameters, which in our case correspond to unknown symbols; the MCRB has been presented in [START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF], [START_REF] D'andrea | The modified cramerrao bound and its application to synchronization problems[END_REF], however, it is looser and moves away from the true CRB, especially for low SNR.

Differently from [START_REF] Meyr | Wiley Series in Telecommunications and Signal Processing[END_REF], [START_REF] Masmoudi | Closed-form expressions for the exact cramer-rao bounds of timing recovery estimators from bpsk, msk and square-qam transmissions[END_REF], the goal and contribution of this paper is to propose a CA timing recovery technique and a closed-form expression of the CRB for CA delay estimation in the case of square-QAM modulated signals.

This paper is organized as follows. In section II, the system model and the CA time delay estimation based on the ML approach are presented. In section III, the CA time synchronization algorithm is proposed. In section IV, closedform expressions of the CRB relative to CA mode are derived. Simulation results are provided in section V and validate our analysis. The last section concludes our work.

II. SYSTEM MODEL

Let us consider the linearly modulated transmitted signal s(t) written as:

s(t) = ∑ i a i h(t -iT ), (1)
where a i denotes the zero mean i.i.d. transmitted symbols drawn from a square-QAM constellation, h(t) is the impulse response of the transmission filter and T is the symbol period. The received signal is:

r(t) = s(t -τ ) + n(t), ( 2
)
where τ is an unknown delay introduced by the channel to the transmitted signal s(t) and n(t) is an additive white Gaussian noise (AWGN) of zero mean and of variance σ 2 .The maximum likelihood estimator of the time delay is given by the following equation [START_REF] Oerder | Derivation of gardner's timing-error detector from the maximum likelihood principle[END_REF]:

τ = arg max u {Λ(u, a)} , ( 3
)
where a is the vector of the transmitted symbols,

Λ(u, a) = exp ( - 1 2σ 2 ∫ T0 |r(t) -s(t -u)| 2 dt ) , ( 4
)
is the likelihood function and T 0 is the observation interval. We can equivalently use the log-likelihood function Λ L (u, a) instead of Λ(u, a). Thus from ( 1)-( 4) the time delay is estimated in the maximum likelihood sense according to:

τ = arg max u {Λ L (u, a)} , (5)
where:

Λ L (u, a) = ∑ k ( ℜ {a * k x k (u)} σ 2 - |a k | 2 2σ 2 ) , (6)
and

x k (u) = y k (u) + v k (u), (7)
y k (u) = ∑ i a i g ((k -i)T -(τ -u)) , ( 8
)
v k (u) = ∫ T0 h * (t -kT -u)n(t)dt, ( 9
)
g(t) = h ⊗ h * -(t), (10)
where ⊗ denotes the convolution operation, ℜ{z} and ℑ{z} are the real and imaginary parts of z, z * is the conjugate of z and for any function f , f -(t) = f (-t).

In equations ( 7)- [START_REF] Noels | Turbo synchronization: an em algorithm interpretation[END_REF], x k (u) is the matched filter output of the received signal, y k (u) is its useful information part, whereas v k (u) is a colored gaussian noise of zero mean and variance σ 2 and g(t) is classically assumed to be a Nyquist pulse. Since it is difficult in practice to maximize Λ L (u, a) and to find u such that ∂ΛL(u,a) ∂u = 0, we propose to implement an adaptive algorithm whose main objective aims at minimizing the derivative of the log-likelihood function towards zero. The updating equation of the adaptive algorithm is given by:

τk = τk-1 + µe k (a k , τk-1 ), ( 11
)
where µ is the step-size and e k (a k , τk-1 ) is the updating error such that:

e k (a k , τk-1 ) = ℜ { a * k ∂x k (τ ) ∂τ | τ =τ k-1 } . ( 12
)
Such approach is called the ML Detector (MLD) [START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF]. The stepsize can be optimized [START_REF] Geller | Equalizer for video rate transmission in multipath underwater communications[END_REF] [25] but this is beyond the scope of our paper.

According to equation [START_REF] Nayak | Joint timing recovery and turbo equalization for coded partial response channels[END_REF], in order to estimate τ we need to know the transmitted symbols a k . This can be made by sending pilot symbols in the case of a DA mode or by making possibly unreliable hard decisions at the receiver in the case of a NDA mode. We hereafter propose to use more reliable softinformation from the turbo-decoder block to derive a code aided estimate of the symbols a k . Also, this optimizes the spectral efficiency as it avoids to send any non informative reference signals.

III. PROPOSED CODE AIDED TIMING DETECTOR

Based on the maximum likelihood approach we propose in the following paragraph to derive a code aided adaptive algorithm for square-QAM modulated signals.

A. Decoder's soft output for Gray-coded constellation

We assume the constellation to be Gray-coded taking values in the alphabet set V = {v 1 , v 1 , ..., v M } where M is the size of the constellation. The chosen i th Gray-coded constellation element v i and the k th transmitted symbol a k are denoted as:

v i ⇔ c i 1 c i 2 ...c i log2(M ) , ( 13
)
a k ⇔ b k 1 b k 2 ...b k log2(M ) , ( 14
)
where c i j (resp. b k j ) corresponds to the j th binary information of v i (resp. a k ).

We can assume that the coded bits in a symbol are statistically independent by using a large-size interleaver [START_REF] Berrou | Near optimum error correcting coding and decoding: turbo-codes[END_REF], [START_REF] Vucetic | Turbo Codes: Principles and Applications[END_REF] so that:

P [a k = v i ] = P [b k 1 = c i 1 , b k 2 = c i 2 , ..., b k log2(M ) = c i log2(M ) ] = log2(M ) ∏ m=1 P [b k m = c i m ] (15)
Let us consider, λ k m the soft output of the decoder at any time index k such that:

λ k m = ln ( P [b k m = 1] P [b k m = 0] ) , m = 1, . . . , log 2 (M ). (16)
The soft-information λ k m is generally obtained after several decoding iterations by a soft decoder. A general expression of the a priori probability of the coded bit b k m , considering that the symbol v i is sent at time k, is given by:

P [b k m = c i m ] = exp ( (2c i m -1) λ k m 2 ) 2 cosh ( λ k m 2 ) , ( 17
)
where c i m ∈ {0, 1}. In order to compute P (a k = ṽi ) we use a Gray-mapping technique for 2 2p -QAM modulated signals.

According to [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] and [START_REF] Bay | On the hybrid cramer-rao bound and its application to dynamical phase estimation[END_REF], for independent symbol bits, the symbol probability is given by:

P (a k = ṽi) = 2p ∏ l=1 1 2 cosh ( λ k l 2 ) β p k 2p ∏ l=1 exp ( (2c i l -1) λ k l 2 ) ξ p k (ṽ i ) . ( 18
)
For notation convenience, we briefly describe the recurrence based construction of the Gray-mapping constellation starting from a 2 2p-2 -QAM Gray coded constellation:

In the upper right quarter of the complex plan, ∀ṽ i ∈ Ṽ2 2p the 2p -2 most significant bits (MSB) are those of the symbols from V 2 2p-2 . Then, ∀ṽ i ∈ -Ṽ2 ) are obtained by symmetry on the y-axis (resp. x-axis) of those of Ṽ2 2p .

Finally, the two least significant bits (LSB) of the the symbols from V 2 2p are those of a Gray-coded QPSK constellation. In other terms:

∀ṽ i ∈ Ṽ2 2p c i 2p-1 c i 2p = 11, ( 19
)
∀ṽ i ∈ (-Ṽ2 2p ) c i 2p-1 c i 2p = 00, ( 20
)
∀ṽ i ∈ Ṽ * 2 2p c i 2p-1 c i 2p = 01, ( 21
)
∀ṽ i ∈ (-Ṽ * 2 2p ) c i 2p-1 c i 2p = 10. (22)
The symbol probability (18) then becomes:

P (a k = ṽi ) = β p k γ p k (ṽi) 2p-2 ∏ l=1 exp ( (2c i l -1) λ k l 2 ) × ( 23
)
exp

( (2c i 2p -1)
λ k 2p 2 + (2c i 2p-1 -1) λ k 2p-1 2
)

.

In [START_REF] Oerder | Derivation of gardner's timing-error detector from the maximum likelihood principle[END_REF], γ p k (ṽ i ) deals with the first 2p -2 most significant bits

(MSB) of ṽi ⇔ c i 1 c i 2 . . . c i 2p-1 c i 2p .
Regarding the construction symmetry of the Gray-code ṽi , -ṽ i , ṽ * i and -ṽ * i have the same 2p -2 MSBs. Accordingly:

γ k p (ṽ i ) = γ p k (-ṽ i ) = γ p k (ṽ * i ) = γ p k (-ṽ * i ). (24)

B. ML based derivation of an adaptive CA delay estimation

In this paragraph, we derive the adaptive CA algorithm equation in the case of a 2 2p -QAM modulated signal (M = 2 2p ) by exploiting the structure of the Gray-mapping technique and the general properties of the constellation alphabet. In fact, an 2 2p -QAM modulated signal takes values into the set

V 2 2p = {v m = ±(2i-1)d p ±j(2n-1)d p ; i, n = 1, 2, . . . , 2 p },
where d p is the intersymbol distance and j 2 = -1. Since we do not know a k , the time delay estimate in the maximum likelihood sense is given by:

τ = arg max u {Λ(u)} . ( 25
)
Λ(u) is obtained by averaging the likelihood function Λ(u, a) over the possible values of the transmitted symbols a k . Thus, by utilizing the independence property between the coded symbols and averaging over the set of possible transmitted symbols we find that:

Λ(u) = ∏ k ∑ vi∈V 2 2p P [a k = v i ] exp ( ℜ{v * i x k (u)} σ 2 - |v i | 2 2σ 2
) .

(

) Let us consider Ṽ2 2p = {ṽ m = (2i-1)d p +j(2n-1)d p , i, n = 1, . . . , 2 p }, so V 2 2p = Ṽ2 2p ∪ (-Ṽ2 2p ) ∪ ( Ṽ2 2p ) * ∪ (-Ṽ * 2 2p ). 26
Then, the likelihood function can be written as:

Λ(u) = ∏ k ∑ ṽi∈ Ṽ2 2p P (a k = ṽi ) exp ( ℜ {ṽ * i x k (u)} σ 2 - |ṽ i | 2 2σ 2 ) + P (a k = -ṽ i ) exp ( - ℜ {ṽ * i x k (u)} σ 2 - |ṽ i | 2 2σ 2 ) + P (a k = ṽ * i ) exp ( ℜ {ṽ i x k (u)} σ 2 - |ṽ i | 2 2σ 2
)

+ P (a k = -ṽ * i ) exp ( - ℜ {ṽ i x k (u)} σ 2 - |ṽ i | 2 2σ 2 ) . ( 27
)
Replacing the symbol probability ( 23) into ( 27) yields the following expression of the likelihood function:

Λ(u) = ∏ k 2β p k σ 2 ∑ ṽi∈ Ṽ2 2p exp ( - |ṽ i | 2 2σ 2 ) γ p k (ṽ i ) × [ cosh ( ℜ {ṽ * i x k (u)} σ 2 + λ k 2p-1 2 + λ k 2p 2 ) (28)
cosh

( ℜ {ṽ i x k (u)} σ 2 - λ k 2p-1 2 + λ k 2p 2 ) ]
.

A similar development of the likelihood function has been made in [START_REF] Bellili | Closed-form cramer rao lower bound for snr estimation from turbo-coded bpsk-, msk-, and square-qammodulated signals[END_REF] to compute the Cramer Rao Bound for the CA SNR estimation and it has been shown that:

γ p k (ṽ m ) = θ k,2p (i)θ k,2p-1 (n), (29)
where ṽm = ±(2i -1)d p ± j(2n -1)d p and θ k,2p (i) and θ k,2p-1 (n) are recursively obtained according to the following equations:

θ k,2p (i) = θ k,2p-2 ( |2i -1 -2 p-1 | + 1 2 ) × exp ( ( 2⌊ i -1 2 p-2 ⌋ -1 ) λ k 2p-2 2
) , (30)

θ k,2p-1 (n) = θ k,2p-3 ( |2n -1 -2 p-1 | + 1 2 ) × exp ( ( 2⌊ n -1 2 p-2 ⌋ -1 ) λ k 2p-3 2
)

, ( 31
)
where ⌊x⌋ is the integer part of x, θ k,1 (1) = 1 and θ k,2 (1) = 1.

Using the identity cosh(a) + cosh(b) = 2 cosh ( a+b 2 ) cosh ( a-b 2 )
, the likelihood function (28) can then be written as:

Λ(u) = ∏ k 4β p k σ 2 H 2p k (r k (u))H 2p-1 k (i k (u)),
where

H l k (x) = 2 p-1 ∑ i=1 θ k,l (i) exp ( - (2i -1) 2 d 2 p 2σ 2 ) × cosh ( (2i -1)d p x σ 2 + (-1) l λ k l 2 ) , ( 32
)
r k (u) = ℜ {x k (u)} and i k (u) = ℑ {x k (u)}.
Then, the derivative of the log-likelihood function is:

∂Λ L (u) ∂u = ∑ k Ḣ2p k (r k (u)) H 2p k (r k (u)) ∂r k (u) ∂u + Ḣ2p-1 k (i k (u)) H 2p-1 k (i k (u)) ∂i k (u) ∂u ,
(33) where for any function f (.), we note ḟ (u) = ∂f (u) ∂u . If we consider the expressions of the real and the imaginary parts of the soft symbol:

ℜ{ã k } = Ḣ2p k (r k (u)) H 2p k (r k (u)) , ( 34
)
ℑ{ã k } = Ḣ2p-1 k (i k (u)) H 2p-1 k (i k (u)) , ( 35
)
then, we find that:

∂Λ L (u) ∂u = ∑ k ℜ { ã * k ∂x k (u) ∂u } . ( 36
)
In practice, we propose to estimate τ , using the updating equation ( 11) by substituting a k with the soft symbol ãk ; the updating equation becomes:

τk = τk-1 + µe k (ã k , τk-1 ), (37)
where the real and imaginary parts of ãk are given by ( 34) and ( 35)

for u = τk-1 and e k (ã k , τk-1 ) = ℜ { ã * k ∂x k (τ ) ∂τ | τ =τ k-1
} . The performance of these estimators can also be compared to the corresponding theoretical Cramer-Rao Bound. The analytical expression of this bound is evaluated in the next section for a code aided time synchronization algorithm.

IV. CRAMER-RAO BOUND

In this section, we derive the Cramer-Rao Bound (CRB) for a code aided delay estimation in the case of square-QAM modulated signals.

Suppose that we are able to produce an unbiased estimate τ of the delay τ from the received signal. Thus, the CRB which verifies E [ (τ -τ ) 2 ] ≥ CRB(τ ) for any unbiased estimator τ of τ is defined as [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]:

CRB(τ ) = I -1 (τ ), (38)
where I(τ ) is the Fisher information matrix (FIM) [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] expressed as:

I(τ ) = E [ ( ∂Λ L (τ ) ∂τ ) 2 ] . ( 39
)
From (33) we have:

E [ ( ∂ΛL(τ ) ∂τ ) 2 ] = ∑ k ∑ n ∑ l∈{2p,2p-1} E [ Ḣl k (z k,l (τ )) H l k (z k,l (τ )) × Ḣl n (z n,l (τ )) H l n (z n,l (τ )) ∂z k,l (τ ) ∂τ ∂z n,l (τ ) ∂τ ] (40)
+ ∑ k ∑ n E [ Ḣ2p k (z k,2p (τ )) H 2p k (z k,2p (τ )) ∂z k,2p (τ ) ∂τ × Ḣ2p-1 n (zn,2p-1(τ )) H 2p-1 n (zn,2p-1(τ )) ∂zn,2p-1(τ ) ∂τ ] ,
where

z k,l (τ ) = { r k (τ ) if l = 2p i k (τ ) if l = 2p -1 (41)
We note that

∂z k,2p (τ ) ∂τ (resp. ∂z k,2p-1 (τ ) ∂τ
) and z k,2p (τ ) (resp. z k,2p-1 (τ )) are uncorrelated (see [START_REF] Masmoudi | Closed-form expressions for the exact cramer-rao bounds of timing recovery estimators from bpsk, msk and square-qam transmissions[END_REF]), therefore the second term of (40) is equal to 0. Deriving the CRB then resorts to compute the first term of (40), which will be carried hereafter distinguishing the cases k = n and k ̸ = n. Case 1:

k = n: E [ ( Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ) 2 ( ∂z k,l (τ ) ∂τ ) 2 ] = E [ ( Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ) 2 ] E [ ( ∂z k,l (τ ) ∂τ ) 2 ] . ( 42
)
On the other hand,

E [ ( ∂z k,l (τ ) ∂τ ) 2 ] = ∑ i E[a 2 i ] ġ2 ((k -i)T ) - σ 2 2 g(0), (43)
where for any function f , we note f (τ ) = ∂ 2 f (τ ) ∂τ 2 . We have:

E [ ( Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ) 2 ] = ∫ +∞ -∞ ( Ḣl k (z k,l (τ )) H l k (z k,l (τ )
)

) 2 P (z k,l (τ ))dz k,l (τ ). (44)
g being a Nyquist filter, the received matched filtered signal can be written as:

x k (τ ) = a k + v k (τ ), (45)
so that from ( 23) and (45):

P (x k (τ )) = ∑ vi∈V exp ( - |x k (τ ) -v i | 2 2σ 2 ) P (a k = v i ) 2πσ 2 = 4β k 2πσ 2 H 2p k (r k (τ ))H 2p-1 k (i k (τ )) × exp ( - r 2 k (τ ) + i 2 k (τ ) 2σ 2 ) . ( 46
)
Given that r k (τ ) and i k (τ ) are respectively the real and imaginary parts of x k (τ ), which are two independent random variables identically distributed, therefore:

P (x k (τ )) = P (r k (τ ))P (i k (τ )), (47)
where

P (r k (τ )) = 2ω k,2p √ 2πσ 2 H 2p k (r k (τ )) exp ( - r 2 k (τ ) 2σ 2
) (48)

P (i k (τ )) = 2ω k,2p-1 √ 2πσ 2 H 2p-1 k (i k (τ )) exp ( - i 2 k (τ ) 2σ 2 ) ( 49
)
and

ω k,2p = p ∏ l=1 1 2 cosh ( λ k 2l 2 ) (50) ω k,2p-1 = p ∏ l=1 1 2 cosh ( λ k 2l-1 2 ) ( 51
)
As a result we obtain:

E [ ( Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ) 2 ] = 2ω k,l √ 2πσ 2 ∫ +∞ -∞ Ḣl k (x) 2 H l k (x) exp ( - x 2 2σ 2
) dx.

(52)

Case 2: k ̸ = n: Given that ∂z k,l (τ ) ∂τ and ∂z n,l (τ ) ∂τ are dependent on z k,l (τ ) and z n,l (τ ), in order to derive the desired expectation in the first term of (40), we first average by conditioning on z k,lk (τ ) and z n,l (τ ) , then we average the resulting expression with respect to these two random variables.

We have:

z k,l (τ ) = a k,l + α k,l (53)
where

a k,l (resp. α k,l ) is equal to ℜ{a k } (resp. ℜ{v k (τ )}) if l = 2p and to ℑ{a k } (resp. ℑ{v k (τ )}) if l = 2p -1. ∂z n,l (τ ) ∂τ = a k,l ġ((n -k)T ) + ∑ j,j̸ =k,n a n ġ((n -j)T ) + αn,l = (z k,l (τ ) ġ((n -k)T ) -ġ((n -k)T )αk, l) + ∑ j,j̸ =k,n a j,l ġ((n -j)T ) + αn,l (54)
Evaluating a conditional mean of (54) we then find that:

E [ ∂z n,l (τ ) ∂τ × ∂z k,l (τ ) ∂τ |z k,l (τ ), z n,l (τ ) ] = z k,l (τ ) ġ((n -k)T )z n,l (τ ) ġ((k -n)T ) = -z k,l (τ )x n,l (τ ) ġ((k -n)T ) 2 , ( 55
)
so that:

E [ ∂z n,l (τ ) ∂τ ∂z k,l (τ ) ∂τ Ḣl k (z k,l (τ )) H l k (z k,l (τ )) × Ḣl n (z n,l (τ )) H l n (z n,l (τ )) |z k,l (τ ), z n,l (τ ) ] (56) = - Ḣl k (z k,l (τ )) H l k (z k,l (τ )) Ḣl n (z n,l (τ )) H l n (z n,l (τ )) z k,l (τ ))z n,l (τ )) ġ((k -n)T ) 2 .
Finally, we obtain:

E [ ∂z n,l (τ ) ∂τ ∂z k,l (τ ) ∂τ Ḣl k (z k,l (τ )) H l k (z k,l (τ )) Ḣl n (z n,l (τ )) H l n (z n,l (τ )) ] = -E [ z k,l (τ )) Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ] E [ z n,l (τ )) Ḣl n (z n,l (τ )) H l n (z n,l (τ )) ] × ġ((k -n)T ) 2 = -E [ z k,l (τ )) Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ] 2 ġ((k -n)T ) 2 , ( 57
)
where

E [ z k,l (τ ) Ḣl k (z k,l (τ )) H l k (z k,l (τ )) ] = 2ω k,l √ 2πσ 2 ∫ +∞ -∞ x Ḣl k (x) exp ( - x 2 2σ 2
) dx.

(58)

We mention that the integrand functions involved in (52) and (58), respectively, decrease rapidly as x increases. Therefore, the integrals over [-∞, +∞] can be accurately approximated by a finite integral over an interval [-A, +A] and the Riemann integration method can be adequately used. The evaluation of the CRB is thus possible as all implied expressions in equation (40) are derived. Results will be displayed in the next section.

V. SIMULATION RESULTS In this section, we display the simulation results of the proposed CA time delay estimation algorithms and the proposed expression of the CRB. We consider the case of 16-QAM and 64-QAM signals with an up-sampling factor equal to 8, passed through a raised cosine filter, with a roll-off factor δ. Results are given for 200 symbol blocks averaged over 1000 Monte Carlo iterations. The LLRs are provided from the encoder which is composed of two identical RSCs concatenated in parallel with systematic rate r = 1 2 and generator polynomials (1, 0, 1, 1) and (1, 1, 0, 1). A large interleaver is placed between the two RSCs. τk is initialized to 0 and its estimated value is depicted at the end of the block when the steady state is achieved. The value of x k (τ k-1 ) is obtained via a quadratic interpolation. In Fig. 1 and 2, we display simulation results of the mean square error (MSE) on (τ /T ) as a function of the SNR for 16-QAM and 64-QAM respectively. For each figure, one can compare the estimation performance of the DA, CA and NDA mode. Compared to the NDA mode, the MSE is decreased by exploiting the soft-information from the turbo-decoder. The DA mode still achieves the best performance, however, it leads to a higher loss of spectral efficiency. Besides, over a wide range of SNRs, the CA mode performance is almost equivalent to that of the DA mode with no need for pilot symbols. For instance, for an SNR= 0 dB the MSE is equal to 6.10 -3 for the DA mode, 9.10 -3 for the proposed technique and 1.8.10 -2 for the NDA mode for 16-QAM signals and it is equal to 1.5.10 -2 for the DA mode, 1.8.10 -2 for the proposed technique and 2.7.10 -2 for the NDA mode for 64-QAM signals. The saturation of the MSE at the right side of Fig. 1 and 2 is due to the self noise of the updating error [START_REF] Nayak | Joint timing recovery and turbo equalization for coded partial response channels[END_REF]. Similar results are obtained for other roll-off factors.

Fig. 1 .Fig. 2 .

 12 Fig. 1. MSE vs SNR for 16QAM modulated signal

Fig. 3 . 3 Fig. 3

 333 Fig. 3. Comparison between the empirical CRB and the analytical expression for different modulations, δ = 0.3

  2p the 2p -2 MSB are obtained by central point symmetry of those of Ṽ2 2p . The 2p-2 MSB of the elements in Ṽ * 2 2p (resp. -Ṽ * 2 2p

VI. CONCLUSION

In this paper, we proposed a new code aided estimation algorithm for time delay recovery of square QAM modulated signals. The proposed CA algorithm performs better than the NDA algorithm and reaches the performance of the DA one over a large interval of SNR values with no need for reference signals. We also established the analytical expression of the CRB for code aided time delay estimation for square-QAM signals. This expression matched the empirical results obtained by simulation. Some future works aim at transposing our approach to the Bayesian canvas [28].

This work has been supported by the Greencocom project of the ANR program


M F Orihuela

M Ortuño

A M M Somoza

J Colchero

E Palacios-Lidón

T Grenet

J Delahaye

Charging of highly resistive granular metal films

Keywords: numbers: 64, 60, ah, 71, 23, An 72, 20, Ee

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

During the last decades, slow conductance relaxation has been studied in many disordered insulators and, in particular, in granular metals by means of field effect measurements. After a quench at low temperatures, a change in the gate voltage is accompanied by a sudden increase in the conductivity, which subsequently slowly decreases with a roughly logarithmic dependence on time [START_REF] Vaknin | Nonequilibrium field effect and memory in the electron glass[END_REF][START_REF] Grenet | Anomalous electric-field effect and glassy behaviour in granular aluminium thin films: electron glass?[END_REF][START_REF] Ovadyahu | Relaxation Dynamics in Quantum Electron Glasses[END_REF][START_REF] Havdala | Ultra-slow relaxation in discontinuous-film based electron glasses[END_REF] . Memory effects and aging are often seen in the same type of experiments [START_REF] Grenet | Manifestation of ageing in the low temperature conductance of disordered insulators[END_REF] . All these glassy effects have been interpreted in different ways, but there is a growing tendency to explain them in terms of electron glasses, i.e., systems with states localized by the disorder and long-range Coulomb interactions between carriers [START_REF] Pollak | A model for the electron glass[END_REF][START_REF] Pollak | The Electron Glass[END_REF][START_REF] Amir | On relaxations and aging of various glasses[END_REF] .

These glassy effects were first limited to low temperature T , but they have recently been observed at room temperature in discontinuous Au [START_REF] Eisenbach | Glassy Dynamics in Disordered Electronic Systems Reveal Striking Thermal Memory Effects[END_REF] , amorphous NbSi 10 and granular Al films 11 . In addition, high resolution techniques such as scanning force microscopy (SFM) open the possibility of studying these systems at room temperature allowing a detailed real space analysis of the problem. Some of us have applied the SFM technique to study slow relaxation in the surface potential of conducting polymers [START_REF] Ortuño | Conducting polymers as electron glasses: surface charge domains and slow relaxation[END_REF] . We observed logarithmic relaxation over four decades of time, as well as full aging in terms of the time of application of the gate voltage. The technique proved appropriate to monitor slow relaxation at the nanoscale.

The use of local probe techniques, such as scanning Kelvin probe microscopy (SKPM), presents two advantages as compared with the conductance measurements performed so far: i) it allows a study of the phenomena at the nanometer scale, which can shed some light on the mesoscopic effects at work in the conductance relaxations, previously observed in indium oxide [START_REF] Orlyanchik | Electron glass in samples approaching the mesoscopic regime[END_REF] and granu-lar Al films [START_REF] Delahaye | Coexistence of anomalous field effect and mesoscopic conductance fluctuations in granular aluminium[END_REF] , but not fully understood; ii) samples with larger resistances can be measured, which is interesting because we know that in the electron glass experiments, the larger the resistance, the larger the conductance relaxations (in % of conductance).

In this problem the first question to answer is the mechanism of charge injection in a strongly disordered insulator. It is an interesting general question and has rarely been addressed. The most thorough study on this subject was the work by Adkins' group [START_REF] Benjamin | Resisitive inhomogeneity in discontinuous metal films[END_REF][START_REF] Adkins | Potential disorder in granular metal systems: The field effect in discontinuous metal films[END_REF] . Analysing non local electrical measurements, they found that inhomogeneities play a crucial role in the charging process and that the scale of inhomogeneities grows as the proportion of metal increases and grains coalesce. They also found strong evidence for the percolative nature of both transport and the charging process. Finally, they monitored the total charge accumulated in the sample as a function of time and observed that the corresponding time dependence was clearly broader than the predictions of diffusion models. In discontinuous metal films and in many other strongly insulating systems, the basic assumption of the conduction models is the exponential-like distribution of the microscopic hopping rates [START_REF] Pollak | The Electron Glass[END_REF] . Its implications in macroscopic transport properties have been systematically confirmed. It gives rise, for example, to the widely observed variable-range hopping [START_REF] Tsigankov | Variable Range Hopping in Two-Dimensional Systems of Interacting Electrons[END_REF][START_REF] Somoza | Collective variable-range hopping in the Coulomb gap: Computer simulations[END_REF] , but no direct microscopic evidence have been produced so far.

In this paper, we study with the SFM technique the time evolution and the spatial dependence of the surface potential of granular Al thin films when a gate voltage is suddenly applied and later on switched off, after a relatively long period of time. We measure samples that are qualitatively similar to those previously studied by some of us in the context of slow relaxation in the conductivity, but that have been grown with the appropriate resistance for this experiment, that is, with a resistance high enough to ensure that the charging processes take place in the scale from milliseconds to hours. As our experimental results are difficult to understand with a model of homogeneous conductivity, we have constructed a more general network resistance-capacitance model able to handle strong inhomogeneities [START_REF] Ortuño | Electronic transport in two-dimensional high dielectric constant nanosystems[END_REF] .

In the next section we describe the SFM technique employed and how the samples have been grown. In section III we present the experimental results, first for the charging experiments and later for the discharging experiments. In section IV we introduce our network model and describe the results of the numerical simulations. Finally, we extract some conclusions.

II. EXPERIMENTAL METHOD A. Sample preparation

The horizontal and vertical schemes of the sample are shown in Fig. 1 (a). The granular Al films used in this study were prepared as described elsewhere by the ebeam evaporation of pure Al under a partial pressure of oxygen 2 . The parameters here are an Al evaporation rate of 1.8 A/s and a partial oxygen pressure around 3 • 10 -5 mbar, for a film thickness of 10 nm. Beyond an oxygen pressure of about 3.2 • 10 -5 mbar, the resistance of the films cannot be measured with our experimental set-up at room temperature (their resistance per square is larger than ≈ 100 GΩ). Previous X rays and TEM investigations have revealed that such films are made of crystalline Al grains of a few nanometers size, dispersed in an amorphous alumina matrix. For the SFM Kelvin probe investigations, the granular Al films were deposited on top of heavily doped Si wafers (the gate) coated with a 100 nm layer of thermally grown SiO 2 (the gate insulator). Al contacts 30 nm thick and evaporated on the granular Al film without opening the e-beam evaporation chamber define channels 30-100 micrometers wide (see the topographic picture). No leaking current to the gate was detected (smaller than 1 pA).

B. SFM microscopy

SFM measurements were carried out using a Nanotec SFM system with a PLL/dynamic measurement board. Once the SFM tip was placed in the granular Al channel between the two lateral contacts topography and surface potential images were acquired simultaneously by means of scanning Kelvin probe microscopy (SKPM) at room temperature and ambient conditions using Pt-coated tips (OMCL-AC240TM-R3, nominal k = 3 N/m). To assure quantitative measurements data were acquired in frequency modulation non-contact dynamic mode (FM-DSFM) with an oscillation amplitude of 2 nm while surface potential images were obtained using FM-SKPM mode with V ac = 500 mV at 7 kHz. Details of the data acquisition set-up can be found in Ref. 20. Freely available WSxM software has been used for image acquisition and processing 21 . To improve the time resolution and to minimize the topography crosstalk artifacts [START_REF] Ortuño | Conducting polymers as electron glasses: surface charge domains and slow relaxation[END_REF] , the charging/discharging experiments were performed on a single scanning spatial line, that is, the y-scanning was blocked and the V SKPM was measured along the same scanning line. This means that in the corresponding surface potential images the horizontal axis is the position along the line, while the vertical axis is time.

Fig. 1 (b) shows a topographic image of the granular aluminum microchannel and (c) the SKPM image acquired simultaneously. Since the maximum lateral scanning length of our set-up is 50x50 microns, the panoramic granular Al channel image has been built up from three images taken across the channel, by slightly displacing the tip from image one to another. The images have been acquired with symmetric potentials on source and drain (V source = V drain = 0 V) after the gate voltage had been kept at 0 V for three days, so that the sample can be assumed at equilibrium. In the topography image (Fig. 1 (b), the electric contacts can be recognized as the higher regions (40 ± 10 nm) at the left and right borders. The lower region corresponds to the very low conductivity granular Al film. In the SKPM image the contrast is reversed and those regions have different surface potential, due to their different work function. Before discussing the charging behaviour of the films, we note that the electric contacts as seen in the topographic and SKPM images do not coincide. In fact, the low conductivity channel as seen by the SKPM image (bright central region, Fig. 1 (c)) is significantly narrower than the channel as measured from the topographic image (lower topographic region, Fig. 1 (b)). As discussed in more detail elsewhere [START_REF] Palacios-Lidon | Quantitative electrostatic Force Microscopy on heterogeneous nanoscale samples[END_REF] , the resolution in the FM-SKPM mode is of the order of 20-50nm, and thus much higher than the discrepancy between the topographic and SKPM images (of the order of 10 µm). Instead, we assume that this discrepancy is due to a shadow effect of the mask used for the evaporation of the gate and drain contacts: if the mask is not sufficiently close to the Algrains, some Al may enter the region between the contacts, increasing the conductivity of the Al-grains near the contacts. Since -within the height resolution of the topographic images-no step is observed in the lower region, we believe that most probably this Al diffuses into the granular Al film. In what follows, we will assume that the true Al-grain sample is the region as seen in the SKPM image.

III. EXPERIMENTAL RESULTS

A. Charging experiments

To study charging processes in the highly resistive granular Al samples, we first let the system relax to equilibrium for a very long time (typically 3 days) by fixing the two leads and the gate voltage at V g = 0. Then, at some instant that is taken as the time origin, t = 0, the gate voltage is switched to V g = 0 and the surface potential V SKPM evolution is monitored. Figure 2 (a) shows V SKPM as a function of the lateral dimension (horizon-tal axis) and the time on a logarithmic scale (vertical axis) in a sample region close to one of the leads. Due to scanning size limitations, only one half of the channel is shown. However, it has been checked that a symmetric behavior take places at the other contact. When the gate voltage V g is switched from 0 to 3 V, the surface potential of the granular Al suddenly changes by the same amount, while the metal lead potential remains unaffected. Just after the voltage switching no current has yet been able to flow into the Al grain film. As a first approximation, at this instant the sample can therefore be considered as a dielectric. The SFM tip therefore "sees" the gate voltage through the granular Al and the SiO 2 insulating layer (see Fig. 1). As time evolves, in order to reach a new equilibrium state, charge from the source and drain contacts enters into the granular aluminum to screen the V g . This is seen best in Fig. 2 (b) where V SKPM lateral profiles at several times have been represented. These results are fairly independent of the value of the applied gate voltage. Note that while the time intervals between successive curves grow roughly exponentially, the change in slope between successive curves shown is similar. Thus, we observe an extremely slow charging behavior, which cannot be explained by a standard diffusive model with low conductivity.

To quantify the behavior of V SKPM the central region of each curve has been fitted to a straight line

V SKPM = a(x -x 0 ) . (1)
All the curves can be overlapped by plotting them as a function of the scaled variable x = (x+x 0 )/a as shown in Figure 2 (c). We first note the relatively high quality of the overlap. In second place, the symmetry between the two tails of the curves is also quite remarkable; we are dealing with an odd function with respect to the value at the center. In the inset of Figure 2 (c), we plot the time dependence of the fitted slopes on a semilogarithmic scale. The roughly logarithmic dependence of the slopes with time is noticeable, a fact difficult to understand with uniform models of the conductivity and a clear indication of hopping processes with exponential distributions of resistances.

B. Memory effects

In our experiments, we have also monitored how the system relaxes back to equilibrium after the gate voltage is switched off. The results are shown in Fig. 3a where V SKPM is plotted on a color scale versus the lateral distance (horizontal axis) and time on a logarithmic scale (vertical axis). They are also plotted in Fig 3b as a function of the lateral distance for several values of the time, with the time intervals between successive profiles increasing exponentially; the origin of time is taken as the instant in which the gate voltage is switched off.

As in the excitation phase, the curves for different times are roughly equally spaced for the exponentially spaced time intervals chosen. To quantify this fact, we measured the slope of the straight segment before the minimum value of V SKPM for the different curves and the results are plotted versus time on a logarithmic scale in panel c of Fig. 3. Within the experimental error, the data follows a straight line with some bending at long times. The continuous curve corresponds to the standard relaxation behavior observed in systems showing logarithmic time evolutions [START_REF] Grenet | Anomalous electric-field effect and glassy behaviour in granular aluminium thin films: electron glass?[END_REF][START_REF] Amir | On relaxations and aging of various glasses[END_REF] x(t) = a ln 1 + t w t .

(

) 2
where the 'waiting' time t w is the time interval during which the gate has been kept at a certain voltage in the excitation phase (here about 5 h).

IV. THEORETICAL ANALYSIS

A diffusion model (or, equivalently, a model with homogeneous conductivity) cannot explain even qualitatively the main features of these experimental results. The value of the diffusion constant D establishes the time scale, but cannot get the wide distribution of relaxation times implied by a quasi logarithmic behavior. One can solve the diffusion problem in terms of modes characterized by a wavenumber k = (πn)/L, where n are integers and L is the distance between the leads. The characteristic time of a mode is 15

τ k = 1 k 2 D , (3)
resulting in a scaling form of the type f (x/ √ t). Each mode reduces proportionally at all positions with an exponential time dependence. The shapes of the potential curves predicted with this type of model are very different from those observed experimentally in the Al-grain samples discussed above. In particular, it is difficult to explain the almost constant value of the potential away from the leads when its value at closer positions has already changed noticeably.

A. Model

The logarithmic dependence of the dynamical variables (inset in Fig. 2c, for example) is indicative of processes with rates depending exponentially on smoothly distributed random variables. Based on this hint and on the highly resistive nature of our samples we study a model suitable for hopping processes and adequate for numerical simulations. In order to keep the model as simple as possible we construct it in two dimensions. We consider a network model that consists on N sites at random on a square sample of size L × L = N and capacitors in the links joining nearest neighbor sites. Each site corresponds to a metallic grain (or set of metallic grains as explained bellow) and each capacitor to the capacitance between nearest neighbor grains [START_REF] Ortuño | Electronic transport in two-dimensional high dielectric constant nanosystems[END_REF] . A scheme of the network model is shown in Fig. 4. In order to perform a simulation as close as possible to the experimental setups and at the same time minimize finite size effects, we have included two leads at a potential V in opposite sides of the sample and periodic boundary conditions in the other two sides. Sites near the leads are connected to them by capacitors. All sites are connected to the gate by capacitances C 0 0 , which take into account the field lines that go out of the system without ending in a nearby grain. The capacitances C i,j between sites take random values extracted from the distribution

C i,j = C r 0 e η , η ∈ [-W/2, W/2] . (4)
We have taken C r = 1 (which sets our unit of energy) and W = 2. The choice of this particular distribution is not crucial, since it just introduces some small randomness in the capacitances and it has the property that the average value of the interaction between nodes is independent of W [START_REF] Ortuño | Electronic transport in two-dimensional high dielectric constant nanosystems[END_REF] . The charge Q i on each site can only take integer values and is split between the plates of the capacitors connected to this site. Each capacitor has opposite charges on its two plates

Q i = j q i,j q i,j = -q j,i = j C i,j (V j -V i ) (5)
where V i is the potential at site i and the sum on j runs over all sites, the two contacts and the gate. The charge on the leads can take any (fractional) value. We can rewrite the previous equation in matrix form

Q i = j A i,j V j (6)
where A i,j is the capacitance matrix, defined as

A i,j = C i,j -C 0 + k C i,k δ i,j . (7)
Eq. ( 6) can also take care of the leads with a proper extension of the definition of the charge vector and of the capacitance matrix. In the leads, the potential V lead is fixed and the charge is V lead C lead , where C lead is the sum of all the capacitances connecting sites of the sample with the lead. It is convenient to define the following charge vector of dimension N + 2 (see Ref. [START_REF] Ortuño | Electronic transport in two-dimensional high dielectric constant nanosystems[END_REF] for details of the model)

Q = (Q 1 , • • • , Q N , C left V left , C right V right ) ( 8
)
and to extend the definition of the capacitance matrix to include the conection with the leads. We usually have an arrangement of charges and want to know the potential at each site. To solve this problem, we have to invert the capacitance matrix F = A -1 . This new matrix F plays the role of an effective interaction between charges. The total energy of the system is directely obatained from F

H = 1 2 QFQ T . (9)
Given a set of charges, we calculate the energy through Eq. ( 9) and, following the hopping model of conductivity in localized systems, perform transitions of charge 1 between sites with probability proportional to 7

Γ i,j = exp - 2r i,j ξ - E i,j kT . (10)
where r i,j is the distance between sites i and j, ξ the localization length, E i,j the energy difference due to the transition, k Boltzmann constant and T the temperature. We consider T = 1 a temperature higher than the average charging energy of a site, but smaller than the applied voltage. Under these conditions, we expect that the factor of distance in Eq. ( 10) will be the crucial one. The energy factor would at most renormalize the percolation effects produced by the spatial factor, explained below.

The sites in our model can represent real metallic grains or effective ones formed by (possibly) many grains merging together by quantum tunneling in highly connected regions. This will just renormalize the energy scale towards smaller values and the time and spatial scales to larger values.

B. Numerical results

We have performed Monte Carlo simulations with our network model changing systematically its parameters to try to reproduce as much as possible the experimental results. The key parameter of the model turns out to be the localization length. When it is larger than the typical site separation, the transition rates between sites are relatively homogeneous and the conduction mechanism can be modeled by diffusion. In the opposite case, hops of charges between close pairs are much more likely than hops between distant sites, the sample is very heterogeneous and conduction can be modeled by percolation theory [START_REF] Pollak | The Electron Glass[END_REF] . As we will see, most experimental features can be explained with the model parameters adjusted to the regime of percolation theory.

If the transition rate between two sites is given by Eq. ( 10) and the localization length verifies ξ 1 electron jumps between close sites are exponentially more likely that between distant sites. Following the standard approach to percolation in hopping systems [START_REF] Pollak | The Electron Glass[END_REF] , we consider that two sites are connected when several hops between them have been produced on average. Then at a given time t sites with a separation smaller than

r i,j ≈ ξ ln t (11)
are connected, while those more separated remain unconnected. As explained above, we have assumed that the spatial factor in the hopping rate, Eq. ( 10), is the relevant one for our conditions. If sites are at random (or grain separations are roughly uniformly distributed), the exponential dependence of the hopping rate dominates the bond connectivity, and the proportion of bonds connected at a given time will be of the form

p(t) ≈ ξ ln t t 0 (12)
where t 0 is a constant that ensures the correct initial condition. Charge can penetrate in the sample along the clusters formed by connected pairs. According to percolation theory, the typical size of these clusters R(t) grows with time as 7

R(t) ≈ (p c -p(t)) -4/3 . ( 13
)
p c is the critical percolation probability for the specific model considered, at which there is an extended connected cluster, and 4/3 is the universal correlation length exponent for percolation in two dimensions. As in Kelvin probe microscopy one adjusts the parameters to avoid as much as possible field lines between the point and the surface, the site potential of our model is directly related to the experimental surface potential. The results of a numerical simulation of the charging process are shown in Fig. 5. The average site potential is plotted as a function of the lateral distance to the electrode for several times and measured from the instant in which the gate voltage is switch on. Note that the time interval between successive curves grow exponentially. The localization length is ξ = 0.2, so conduction is well in the percolative regime. The overall shape of the curves are similar to the experimental results, but more importantly the time dependence of the slopes of the potential is roughly logarithmic, quite similar to the experimental behavior. In the inset of Fig. 5 we represent the slope of the surface potential versus time on a logarithmic scale. According to our model, the exponentially large distribution of characteristic times is due to the exponential dependence of the hopping rate on distance (and to a roughly uniform distance distribution). At low temperatures, the energy factor in the hopping rate can contribute further to the broadening of the times involved. As our experiments are at room temperature, we do not expect a significant contribution from the energy factor in the hopping rates. The gate voltage potential (10 for the results in Fig. 5) is chosen so that the charge in each grain is clearly higher than one. In this case the results are independent of the gate voltage, as in the experiments.

After the last time for which we measure the charging of the sample, we switch off the gate potential in our simulation and monitor again the time and spatial evolution of the site potential. The results are plotted in Fig. 6 for the same parameters as in the charging process, Fig. 5. The main features of the discharging curves are again reproduced by our simulations, although there are differences in the details. The overall logarithmic behavior is well reproduced as can be seen in the inset of Fig. 6, where the slopes of the first part of the curves are plotted versus time on a logarithmic scale. The continuous line corresponds to a fit to Eq. ( 2).

We can have a more restrictive test of our percolation model. If one pays attention to the experimental curves, Fig. 2, one can observe that they seem to accelerate (in a logarithmic scale) their approach to equilibrium at large time intervals. This is clearly evident in panel a where the initially almost straight lines become bent at the bottom of the figure. We interpret this as a sign that, at the largest times involved in our experiments, we are close to reaching percolation, so that the size of the connected clusters increases faster, according to Eq. 14, as the system approaches critical percolation and the whole sample would be connected. To quantify this observation we plot the lateral distance of the point at which the voltage is equal to the mean between the lead and the gate voltages as a function of time. In Fig. 7, we show the results for both the experimental data (black points, left axis in µm) and the numerical simulations (red points, right axis in grains separation). The horizontal time axis is in seconds for the experimental points and in Monte Carlo sweeps for the simulation results, and it has been shifted for the simulation results to achieve a maximum overlap between both set of points. The line is a fit of the experimental data to the analytical expression of the typical cluster size R in percolation theory

R(t) = a ln t c t -4/3 . ( 14
)
To derive this equation, we have taken into account the time dependence of the connecting probability, Eq. ( 12), t c is the critical time needed to reach percolation, and 4/3 is the correlation exponent in two-dimensional percolation. The quality of the fit and the analogy between experimental and numerical simulation results is strong indication of the percolative character of the charging process.

V. CONCLUSIONS

We have repeated the previous measurements for other samples with different resistances. If the samples are conducting enough so that its conductivity can be measured, the processes studied here take place at a fast scale and our set up is not fast enough to detect intermediate processes. If the samples are much more resistive than the ones analyzed here, we do not see any dynamics. As resistance varies exponentially with grain separation, we can only measure a relatively narrow range of metallic concentrations.

In principle, this technique can be adapted to measure the conductivity of highly resistive samples for which standard macroscopic methods are not sensitive enough to determine it. We now intend to apply the SKPM technique to study at the nanometer scale the glassy effects reported in the conductance of several disordered insulators. [START_REF] Grenet | Anomalous electric-field effect and glassy behaviour in granular aluminium thin films: electron glass?[END_REF][START_REF] Ovadyahu | Relaxation Dynamics in Quantum Electron Glasses[END_REF][START_REF] Havdala | Ultra-slow relaxation in discontinuous-film based electron glasses[END_REF][START_REF] Grenet | Manifestation of ageing in the low temperature conductance of disordered insulators[END_REF] . Indeed in some of them, in the appropriate resistance range, these occur up to room temperature. [START_REF] Eisenbach | Glassy Dynamics in Disordered Electronic Systems Reveal Striking Thermal Memory Effects[END_REF][START_REF] Delahaye | Observation of thermally activated glassiness and memory dip in a-NbSi insulating thin films[END_REF][11] The technique can be very useful to analyze memory effects, and combined with conductance measurements, may shine light on the mechanisms involved.

FIG. 1 .

 1 FIG. 1. (a)Top and side views of the scheme of the sample. (b) Topography (∆z = 35 nm) and (c) SKPM images of the granular aluminum microchannel (z scale 400 mV).

FIG. 2 .

 2 FIG. 2. (a) Surface potential image in a region close to a contact (distance in horizontal axis and log t in the vertical axis). The horizontal dashed line indicates t = 0 when the gate voltage is switched from 0 V to 3 V. The vertical dashed line marks the position of the effective contact lead. (b) Surface potential profiles from panel (a) at different times. The profile at Vg = 0 has been included as a reference. (c) Same data as in (b) represented versus the scaled lateral distance. Inset: Slope of VSKPM versus time on a logarithmic scale

FIG. 3 .

 3 FIG. 3. (a) Color map of the data in panel a with distance in horizontal axis and log t in the vertical axis. (b) Surface potential as a function of distance at the time intervals exponentially spaced with respect to the instant when the gate voltage has been switched back to zero. (c) Slopes of the decreasing regions of the curves in panel (b) plotted versus time. The continuous curve is a fit to Eq. (2).

FIG. 4 .

 4 FIG.4. Top and side views of the scheme of the capacitors network used to model the system. Charges can jump between nodes, which is equivalent to having resistances between them (not shown).

FIG. 5 .

 5 FIG.5. Site potential of our network model as a function of distance for several times exponentially spaced. The localization length is ξ = 0.2. Inset: Slope of the site potential at a given distance versus time on a logarithmic scale.

FIG. 6 .

 6 FIG. 6. Numerical simulation of the site potential as a function of distance for several times exponentially spaced after the gate voltage has been switch back to zero. Inset: Slope of the first part of the curves of the site potential versus time on a logarithmic scale.

FIG. 7 .

 7 FIG. 7. Lateral distance of the point where the surface potential is equal to the mean value between the potential in the lead and in the middle of the sample as a function of time for the experimental results shown in Fig. 2 (black points, bottom and left axis) and for the simulations (red points, top and right axis).

ACKNOWLEDGEMENTS

The work was supported by the Spanish MINECO Grants No. FIS2015-67844-R, ENE2013-48816-C5-1-R, ENE2016-79282-C5-4-R, and CSD2010-00024, and by Fundación Séneca grant 19907/GERM/15.


Jeremy Bleyer

Jean-François Molinari

Microbranching instability using a variational phase-field model

come

Microbranching instability using a variational phase-field model

Many aspects of the microbranching instability [START_REF] Sharon | Microbranching instability and the dynamic fracture of brittle materials[END_REF] are not fully understood such as the difference of crack patterns observed in PMMA on one hand and glass or brittle gels on the other hand [START_REF] Bouchbinder | Nonuniversality in microbranching instabilities in rapid fracture[END_REF]. In this work, we simulate the microbranching instability of fast propagating cracks using a variational phase-field approach and show that three-dimensional effects are fundamental for understanding the origin of this instability.

In particular, following on a previous work regarding the branching mechanism in 2D phase-field simulations [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF], we show that microbranching does not occur in 2D simulations or for very thin 3D samples and that different microbranching regimes are obtained when increasing the sample thickness, ranging from translationally invariant across the thickness and quasi-periodic patterns to localized microbranching events driven by complex 3D dynamics (cf. Figure ). We further show that the transition between these different regimes is controlled by the value of the regularization length scale inherent to the phasefield modelling of brittle fracture.

Université Paris-Est, Laboratoire Navier, (Ecole des Ponts-IFSTTAR-CNRS UMR 8205), jeremy.bleyer@enpc.fr

Ecole Polytechnique Fédérale de Lausanne, Département de Génie Civil, {jeremy.bleyer,jean-francois.molinari}@epfl.ch


A Dantan

C Genes

D Vitali

M Pinard

Self-cooling of a movable mirror to the ground state using radiation pressure

Keywords: numbers: 42.50.Lc, 03.67.Mn, 05.40.Jc

We show that one can cool a micro-mechanical oscillator to its quantum ground state using radiation pressure in an appropriately detuned cavity (self-cooling). From a simple theory based on Heisenberg-Langevin equations we find that optimal self-cooling occurs in the good cavity regime, when the cavity bandwidth is smaller than the mechanical frequency, but still larger than the effective mechanical damping. In this case the intracavity field and the vibrational mechanical mode coherently exchange their fluctuations. We also present dynamical calculations which show how to access the mirror final temperature from the fluctuations of the field reflected by the cavity.

Cooling of mechanical resonators both at the microand at the macro-level, is an important technical challenge in various fields of physics, such as ultra-high precision measurements [1], and detection of gravitational waves [2]. It is also a prerequisite for any possible use of optomechanical systems for quantum information processing [START_REF] Mancini | [END_REF]4]. Active noise control techniques have been proposed to reduce their thermal noise and bring the oscillator motion to its ground state [5]. Recently, various experiments have demonstrated significant cooling of the vibrational mode coupled to an optical cavity [6,7,8,9,10,11,12,13,14,15]. The experiments of Refs. [8,10,12,13] in particular have adopted the socalled back-action [16], or self-cooling scheme, in which the radiation pressure of an appropriately detuned cavity interacting with the mechanical oscillator, is used. In all these experiments, however, the resulting equilibrium state of the oscillator is classical, because the new mean excitation number is still much larger than one. Therefore it is important to establish the fundamental limits of self-cooling and if it is able to cool a mechanical degree of freedom down to its quantum ground state. We will describe the system dynamics in terms of quantum Langevin equations (QLE) for the Heisenberg operators of the system and we will show that ground state selfcooling is possible in the good cavity regime, i.e. when the mechanical resonance frequency is larger than the cavity bandwidth, and provided the cavity detuning and bandwidth are appropriately chosen so as to maximize the scattering of noisy photons into the cavity mode.

We consider a driven optical cavity coupled by radiation pressure to a micromechanical oscillator. The typical experimental configuration is a Fabry-Perot cavity with one mirror much lighter than the other (see e.g. [7,8,9,10,11,15]), but our treatment applies to other optomechanical systems, such as the silica toroidal microcavity of Refs. [12,17]. Radiation pressure typically excites several mechanical degrees of freedom of the system with different resonant frequencies. However, we consider the detection in a narrow frequency bandwidth, containing a single mechanical resonance peak. In this case the oscillator dynamics can be well approximated by that of a single harmonic oscillator with dimensionless position q(t), and momentum p(t), ([q, p] = 2i), frequency Ω m , mass M , and damping Γ. The coupled dynamics of the mirror with the intracavity field mode a (t) in the frame rotating at the frequency of the driving laser ω L , is described by the QLE

q = Ω m p (1) ṗ
= -Ω m q -Γp + G √ 2a † a + ξ (2) ȧ = -(κ + i∆ c ) a + iGaq/ √ 2 + √ 2κa in , (3)
where G = (ω c /L) /M Ω m is the opto-mechanical coupling constant, ω c and L being the cavity resonance frequency and length, respectively. a in (t) is the input field, satisfying a in (t)a in † (t ′ )a in 2 = δ(t -t ′ ), with a in the incident mean field. κ is the damping rate of the cavity mode, ∆ c = ω c -ω L the cavity detuning and ξ the noise operator accounting for the mirror Brownian motion at thermal equilibrium at temperature T [START_REF]ξ(t)ξ(t ′ ) = (Γ/πΩm) dΩe -iΩ(t-t ′ ) Ω[coth( Ω 2k B T ) + 1[END_REF].

Steady state analysis -In steady state the mean intracavity field a is given by a = √ 2κa in /(κ + i∆), where ∆ is the mean detuning of the cavity, given by ∆ = ∆ c -∆ nl = ∆ c -G 2 |a| 2 /Ω m . These two coupled equations give a third-order relation between a and ∆ which leads to the well-known bistable behavior of a cavity with a movable mirror [START_REF] Dorsel | [END_REF]. The stability condition of the system can be written as κ 2 + ∆ 2 + 2∆∆ nl > 0.

If we then consider the linearized fluctuations of the various operators around the steady state, we get linearized QLE which can be solved by moving to the frequency domain. The Fourier transform of the position fluctuations can then be simply expressed as the sum of a thermal noise term and a radiation pressure term

q[Ω] = χ[Ω] (F R [Ω] + F T [Ω]) , (4)
the response to the noise terms being given by an effective susceptibility

χ[Ω] -1 = χ[Ω] -1 -M Ω 2 m 2∆ nl ∆ κ 2 D[Ω] (5)
with

χ[Ω] = M Ω 2 m 1 -Ω 2 /Ω 2 m -iΓΩ/Ω 2 m and D[Ω] = (1 -iΩ/κ) 2 + (∆/κ) 2 . The thermal noise force is simply F T [Ω] = M Ω m ξ[Ω] and the radiation pressure force [20] F R [Ω] = √ 2M Ω 3 2 m ∆ 1 2 nl κ 3 2 D[Ω] (κ -iΩ)x in [Ω] + ∆y in [Ω]
with x in = a in + a in † and y in = i(a in † -a in ). An exact expression for the mirror variances is then obtained by integrating the contributions of these two forces to the noise spectrum (ω = Ω/Ω m )

∆q 2 = dω 2π S q (ω), ∆p 2 = dω 2π ω 2 S q (ω) (6)
where we have used the fact thatp = q/Ω m and where the noise spectrum is given by

S q (ω) =   2ω coth Ωmω 2kB T Q + 4ϕ nl 1 + ϕ 2 + b 2 ω 2 (1 -b 2 + ϕ 2 ) 2 + 4b 2 ω 2   |(1 -ibω) 2 + ϕ 2 | 2 |[(1 -ibω) 2 + ϕ 2 ][1 -ω 2 -iω/Q] -2ϕϕ nl | 2 , ( 7
)
with b = Ω m /κ, ϕ = ∆/κ, ϕ nl = ∆ nl /κ and Q = Ω m /Γ.
These two exact expressions for the variances coincide with Eq. 5 of Ref. [21]. Cooling to the ground state implies reaching ∆q 2 = 1 simultaneously with ∆p 2 = 1. Fig. 1 shows both variances as a function of b, for a cavity detuning equal to one bandwidth (ϕ = b), which shows that, starting from a mean thermal excitation number

n i T = (exp{ Ω m /k B T } -1)
-1 equal to 10 2 , it is possible to achieve a lower than unity final excitation number state in the good cavity limit: b > 1, i.e. Ω m > κ. Physical interpretation in the adiabatic limit -To get more insight into the cooling mechanism we will now derive an analytical expression for the frequency integrals of (6). The latter may be considerably simpli-fied by noting that, for high-Q cavities and under certain conditions that we will detail, the mirror response is peaked around the natural mechanical resonance frequency. One may then express the effective susceptibility as χ

[Ω] -1 ≃ M Ω2 [1-Ω 2 / Ω2 -i/ Q]
with an effective quality factor Q = Ω/ Γ and effective resonance frequency and relaxation rate

Ω = Ω m 1 -2ϕϕ nl Re[D[Ω m ] -1 1/2 (8) Γ = Γ 1 + 2ϕϕ nl Q Im[D[Ω m ] -1 (9)
The effect of the frequency-shift in the resonance is typically negligible, so that one may assume Ω ≃ Ω m . The relaxation rate of the mirror is on the contrary strongly affected by the presence of the cavity field and, depending on the sign of the cavity detuning, may either be enhanced or reduced, resulting in either cooling or heating. As pointed out in [22,23,24], the most interesting regime for self-cooling is when the effective damping rate Γ is strongly enhanced, but still less than the cavity bandwidth and when the effective mechanical Q is still larger than one, Γ ≪ Γ ≪ κ. In this case, the effective susceptibility is still peaked around ω = 1 and we can well approximate the smoothly varying function of ω inside the square brackets by taking its value at ω = 1. One then gets an analytical expression for the variance

∆q 2 ≃ Γ Γ 2n i T + 1 + 2ϕ nl Q 1 + b 2 + ϕ 2 (1 -b 2 + ϕ 2 ) 2 + 4b 2 (10)
Eq. ( 10) is the basic equation for our analysis of the quantum limits of self-cooling, and it is possible to see that it coincides with the results of [22,24], which were obtained with a different method. Ref. [21] instead gives slightly more general expressions for the two variances, which, however, reduce to that of Eq. (10) when Ω m is large enough, i.e, Ω m > Γ.

Depending on the respective frequencies involved, one can adopt different strategies to cool the mirror down to the ground state. In particular, the cavity detuning has to be well-chosen in order to optimize the self-cooling. Fig. 2 shows the variation of the variance with the cavity detuning, as calculated from the exact solution of Eq. ( 7) and from the harmonic approximation (10). The analytical result is in good agreement with the numerical result and we will use it as a basis for the cooling optimization discussion. In order to interpret Eq. ( 10) we rewrite it as

∆q 2 = 1 + 2n i T 1 + f ϕ nl Q + f ϕ nl Q 1 + f ϕ nl Q 1 + b 2 + ϕ 2 2ϕb , (11)
or, equivalently,

∆q 2 = (1 -η) ∆q 2 T + η ∆q 2 R . ( 12
)
∆q 2 T = 1 + 2n i T and ∆q 2 R = (1 + b 2 + ϕ 2 )/2ϕb represent the thermal noise-induced and the cavity field-induced contributions, respectively, to the position fluctuations of the mirror.

η = f ϕ nl Q 1 + f ϕ nl Q (13)
is the weight between those two quantities, with

f = 4ϕb (1 -b 2 + ϕ 2 ) 2 + 4b 2 . ( 14
)
The physical interpretation of Eq. ( 12) is that the cavity field and the mirror vibrational mode coherently exchange their fluctuations during the interaction, all the more so that the quantity η is close to unity. The thermal excess noise is then replaced by the cavity field fluctuations, which are those of the vacuum, hence resulting in decreasing the effective temperature of the mirror. The quantity η naturally appears as a quantum state transfer efficiency between the field and the vibrational mode, whereas the product ϕ nl Q can be interpreted as a coherent coupling strength, quite similarly to the cooperativity parameter, which appears in quantum state transfer schemes between optical fields and atomic ensembles [25,26,27]. Since the mirror is initially in a noisy thermal state, the cavity field acts as a thermal noise "eater" during the interaction, and allows to reach a much lower effective temperature for the mirror in steady state. This situation is also quite similar to cavity assisted Dopplercooling of atoms [28], for which the temperature of the atoms is decreased by enhancing the scattering of photons into the cavity mode.

If one starts with a very large thermal excess noise, Eq. ( 11) clearly shows that the thermal noise contribution will be suppressed by a factor f ϕ nl Q, which is typically large. Minimizing the mirror final temperature is then equivalent to maximizing the transfer efficiency η, and, for a given value of ϕ nl Q, maximizing f with respect to ϕ. This gives an optimal cavity detuning equal to

ϕ * = [(b 2 -1 + 2 1 + b 2 + b 4 )/3] 1 2 . ( 15
)
The most favorable conditions to reach the quantum ground state of the mechanical oscillator then occur for a cavity bandwidth smaller than the mechanical resonance. Indeed, assuming b ≫ 1, one has ϕ * ≃ b, f * ≃ 1 and for a sufficiently high Q-factor (ϕ nl Q ≫ 2n i T + 1), the resulting mirror fluctuations are then mostly given by the radiation pressure:

(∆q 2 ) * ≃ ∆q 2 R ≃ 1 + 1 2b 2 (16)
and the normalized variance tends to unity for b ≫ 1.

It is therefore possible to reach the mechanical ground state using self-cooling in the good cavity limit. However, there exists an optimal bandwidth that minimizes the mirror final temperature. Indeed, as one decreases the cavity bandwidth (for a fixed value of the mechanical resonance frequency), the effective response time of the mirror becomes comparable to that of the intracavity field and the previous adiabatic approximation is no longer valid. One can also show that this situation is detrimental to the noise transfer between the field and the mirror [25,26]. As mentioned previously and as can be seen in Fig. 1, the adiabatic approximation breaks down when ϕ nl Ω m 2κ.

Dynamical situation -We now address the dynamical behavior of the coupled field-mirror system. Under good self-cooling conditions one expect the mirror initial thermal noise to be mainly transferred to the cavity field during the cooling phase, and hence to the field leaking out of the cavity. Numerically solving the linearized Langevin equations in time allows to compute the twotime correlation functions of q(t) and p(t), as well as those of the outgoing field quadratures, x out and y out , where a out = √ 2κa -a in . To simplify, we assume a step-like injection of a coherent state cooling field into the cavity at t = 0 and start with an occupation number n i T = 10 2 for the mirror [START_REF]This analysis can straightforwardly be extended to arbitrary non-stationary interaction[END_REF]. The resulting mirror variances at time t are represented in Fig. 3, as well as the temporal evolution of the equal-time variances of the outgoing field C x,y (t, t), obtained from exact numerical calculations. It appears clearly that, during the cooling phase, thermal noise is taken from the mirror and transferred to the cavity field. The amount of thermal noise taken from the mirror can be deduced by measuring the outgoing field fluctuations. Indeed, for Ω m = ∆ and in the adiabatic limit Γ ≪ Γ ≪ κ, it can be shown that the outgoing field correlation functions are of the form [23,[START_REF]which the vacuum contribution has been subtracted, ie Cu(t, t ′ ) = u out (t)u out (t ′ ) -δ(t -t ′ ) (u = x, y)[END_REF] C x,y (t, t ′ ) = 2η Γe -Γ(t+t ′ ) + n i T η(1 -η) Γe -|t-t ′ | (17

)
Following the method developed in Ref. [START_REF] Dantan | [END_REF], one can perform a homodyne detection of the outgoing field fluctuations with a temporally matched local oscillator, E LO (t) ∼ e -Γt [32]. The measured quantity is then ∆x 2 m = 1 + η(n i T -n f T ) + η(1 -η)n i T , where n i,f T are the initial and final thermal excitation numbers. When the cooling is optimal (η ∼ 1), one indeed measures the change in temperature of the mirror: ∆x 2 m ≃ 1+n i T -n f T .

We have presented a general theory for the self-cooling of a mechanical oscillator to the ground state, which allows for deriving analytical results for the mirror final temperature and provides a simple interpretation of the cavity mediated self-cooling process. When the cavity is suitably detuned the initial thermal noise of the mirror vibrational mode is essentially transferred to the cavity field mode. This noise exchange reflects in the field leaking out of the cavity, providing useful information on the mirror temperature.

FIG. 1 :

 1 FIG. 1: Normalized variances versus b, for a detuning ϕ = b Q = 10 4 , n i T = 10 2 , ϕ nl = 0.1. The final excitation number is n f T ∼ 0.15 for b ∼ 5.

2 FIG. 2 :

 22 FIG. 2: Normalized variance versus cavity detuning [dashed:exact result of Eq. (7), plain: approximate result of Eq. (10)]. Parameters: Q = 10 4 , n i T = 10 2 , ϕ nl = 0.1, b = 10.

FIG. 3 :

 3 FIG. 3: (a) Temporal evolution of the mirror position and momentum normalized variances when the cooling field is on from t = 0 (Γ -1 units). (b) Temporal evolution of the outgoing field variances Cx,y(t, t) (in κ units). Parameters: κ = 10 3 Γ, Q = 10 4 , n i T = 10 2 , ϕ nl = 0.1, ϕ = b = 10.