Matthieu Martel

Amine Najahi

Guillaume Revy

Toward the synthesis of fixed-point code for matrix inversion based on Cholesky decomposition

Keywords: Cholesky decomposition, matrix inversion, fixed-point arithmetic, automated code synthesis, certified numerical accuracy

Matrix inversion is a computationally intensive basic block of many digital signal processing algorithms. To decrease the cost of their implementations, programmers often prefer the fixed-point arithmetic. This arithmetic requires less resources and runs faster than the floating-point arithmetic, but all the arithmetical details must be handled by the programmer. In this article, we overcome this drawback by presenting an automated approach to synthesize fixed-point code for matrix inversion based on Cholesky decomposition. First we rigorously define the square root and division operators especially in terms of rounding error, and we implement them in the CGPE library. This allows us to provide accuracy certificates for the generated code. Second we propose a workflow based on Cholesky decomposition that carefully uses these operators to produce accurate code for the basic blocks of matrix inversion. Finally we illustrate the efficiency of our approach on some benchmarks, and show how it allows us to synthesize accurate code in a few seconds and thus to reduce the development time of fixed-point matrix inversion.

I. INTRODUCTION

Matrix inversion is known to be numerically unstable. As a consequence, numerical analysts advise against using it for a large set of problems. Yet, as stated by Higham [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]Ch. 14], cases exist where the inverse conveys useful information. In the wireless communication field, matrix inversion is used in equalization algorithms [START_REF] Zhou | A novel adaptive equalization algorithm for MIMO communication system[END_REF] as well as detection estimation algorithms in space-time coding [START_REF] Chen | Layered turbo space-time coded mimo-ofdm systems for time varying channels[END_REF]. This article does not discuss these cases and leaves it to the programmer to decide whether computing the inverse is indeed justified. The goal of this article is to introduce techniques that help the programmer in the process of writing fixed-point code for matrix inversion. The originality of our work relies on the fact that certified error bounds are also produced together with the code. In a typical design, DSP programmers prototype and simulate their algorithms in high level environments like MATLAB. These environments work with the floating-point arithmetic [START_REF][END_REF] to ease and speedup the prototyping phase. However, when mapping the design to hardware, constraints on silicon area, power consumption, or throughput frequently force the implementer to convert this design to the more efficient fixed-point arithmetic [START_REF] Yates | Fixed-Point Arithmetic: An Introduction[END_REF]. This conversion is known to be a tedious and time consuming process [START_REF] Menard | Floating-to-fixed-point conversion for digital signal processors[END_REF] that may be split into two phases:

1) Range analysis: This phase allows to find the integer wordlength of each variable in the design. In a finite wordlength environment, minimizing the integer word length allows one to allocate more digits for the fractional part, thus obtaining more accuracy.

2) Precision analysis: In this phase, the number of bits to allocate to the fractional part is decided. This phase must take into account the precision needs of the application.

Over the last years, authors have suggested different strategies to tackle these conversion phases. These contributions fit into two categories:

1) Simulation based strategies [START_REF] Nikolic | Design and Implementation of Numerical Linear Algebra Algorithms on Fixed-Point DSPs[END_REF], [START_REF] Sung | Simulation-based word-length optimization method for fixed-point digital signal processing systems[END_REF]: The information that allows to estimate the required range and precision are inferred from intensive simulations carried out using a precise arithmetic, typically floating-point arithmetic. 2) Analytic strategies [START_REF] Lee | Accuracy-guaranteed bit-width optimization[END_REF], [START_REF] Fang | Fast, Accurate Static Analysis for Fixed-Point Finite-Precision Effects in DSP Designs[END_REF]: The information is obtained using formal methods such as interval arithmetic, affine arithmetic, and norm computation for digital filters. The precision analysis phase relies on optimization techniques.

In this work, we use an analytic approach based on interval arithmetic to bound the range of the variables in our design and to give strict bounds on the rounding errors. Previous works applied similar techniques but only to a small range of problems involving only additions and multiplications. This article describes algorithms to compute square roots and divisions, and the means to bound their rounding errors. By adding these operators to the CGPE tool, we are able to treat problems, previously considered intractable in fixed-point arithmetic, like Cholesky decomposition and triangular matrix inversion.

This article is organized as follows: Section II gives some insights on matrix inversion algorithms. Section III introduces some basics on fixed-point arithmetic, while Section IV details our new square root and division operators. Section V details the workflow we propose to map matrix inversion into the fixed-point arithmetic. Then some experimental results are exhibited in Section VI, before concluding in Section VII.

II. BACKGROUND ON MATRIX INVERSION

A survey of floating-point matrix inversion and linear systems solving shows that there are many algorithms in use: LU and QR decompositions, Cramer's rule, . . . [START_REF] Golub | Matrix Computations[END_REF]. A common pattern to the efficient algorithms is the decomposition of the input matrix into a product of easy to invert matrices (triangular or orthogonal matrices). LU decomposition for instance, proceeds by Gaussian elimination to decompose an input matrix A into two triangular matrices L and U such that A = LU . Inverting triangular matrices being straightforward, solving the associated linear system is equally simple and so is obtaining the inverse A -1 by the formula

A -1 = U -1 L -1 .
Almost the same chain of reasoning is applicable to QR decomposition. Clearly, the computationally intensive step of these algorithms is the decomposition part [START_REF] Golub | Matrix Computations[END_REF].

The decomposition part is also the missing link in fixedpoint arithmetic. Indeed, some works exist on fixed-point code synthesis for matrix multiplication [START_REF] Martel | Code Size and Accuracy-Aware Synthesis of Fixed-Point Programs for Matrix Multiplication[END_REF], but, to the authors knowledge, no published material suggests a rigorous methodology of code synthesis for matrix decomposition or triangular matrix inversion. The authors of the Gusto tool [START_REF] Irturk | GUSTO: An Automatic Generation and Optimization Tool for Matrix Inversion Architectures[END_REF] do provide benchmarks for the accuracy of their matrix inversion algorithms. However, the methodology for generating the fixed-point implementation is not explicit in the article and the treated matrices do not exceed size 8. Besides, the evaluation of the rounding error is an a posteriori estimation process. Therefore, no absolute error bounds are provided. In [START_REF] Nikolic | Design and Implementation of Numerical Linear Algebra Algorithms on Fixed-Point DSPs[END_REF], Frantz et al. use an approach based on the ideas introduced by [START_REF] Sung | Simulation-based word-length optimization method for fixed-point digital signal processing systems[END_REF] to study the mapping of many linear algebra routines to fixed-point DSPs. They treat matrices of sizes up to 35, but uniquely with simulation based methods and without providing certificates on the error bounds. In introducing a rigorous error model, our work is analogous to [START_REF] Lee | Accuracy-guaranteed bit-width optimization[END_REF] and [START_REF] Fang | Fast, Accurate Static Analysis for Fixed-Point Finite-Precision Effects in DSP Designs[END_REF], where affine arithmetic is used to assert properties on the error bounds. However, both approaches target discrete transformation algorithms and only [START_REF] Lee | Accuracy-guaranteed bit-width optimization[END_REF] treats briefly of small size matrix multiplication. In our approach, we suggest an interval arithmetic based formalization of all the arithmetic operators that intervene in the computation of classical linear algebra decompositions. This allows us to compute explicit bounds on the rounding errors of the codes synthesized by our tools. In the following of this section, we recall the formulas for triangular matrix inversion and Cholesky decomposition, which are useful to solve linear systems and invert matrices.

A. Triangular matrix inversion

Using the so called backward and forward substitution techniques, inverting a triangular matrix is a straightforward process in floating-point arithmetic. Indeed, for a lower triangular matrix A, its inverse N is given by the following equations:

n i,j =        1 a i,i if i = j -c i,j a i,i if i = j where c i,j = i-1 k=j a i,k • n k,j . (1)
In floating-point arithmetic, these equations are easy to implement. It is less so in fixed-point arithmetic, since the coefficient n i,j depends on other coefficients of N , namely all the n k,j with k ∈ {j, . . . , i -1}. This implies that the synthesis tool, when generating code that computes n i,j must know the ranges and formats of all the n k,j with k ∈ {j, . . . , i -1}. It is clear that such a tool must follow a determined order in synthesizing code and that it must keep track of the formats and ranges of the computed coefficients so as to reuse them. Besides, similarly to the work in [START_REF] Martel | Code Size and Accuracy-Aware Synthesis of Fixed-Point Programs for Matrix Multiplication[END_REF], this process involves multiple trade-offs between code size and accuracy.

B. Cholesky decomposition

If A is a symmetric positive definite matrix, its Cholesky decomposition is defined. Since this method exploits the structure of the matrix, it is more efficient than Gaussian elimination. It is also known for its numerical stability. The aim of Cholesky's method is to find a lower triangular matrix L such that

A = L • L T . ( 2
)
By equating the coefficients in (2) and using the properties of A, the following two formulas for the general term of L are deduced:

i,j =      √ c i,i if i = j c i,j j,j if i = j with c i,j = a i,j - j-1 k=0 i,k • j,k . (3)
Generating code for c i,j is easily achievable in fixed-point arithmetic since it involves a subtraction and a size-j dot product. Synthesizing code that computes i,j from c i,j is less obvious since the square root and division operators are not well studied in fixed-point arithmetic. To conclude this section, note that Cholesky decomposition can be used to invert any non symmetric positive definite matrix A by decomposing the following matrix: M = AA T which is guaranteed to be symmetric to obtain M = LL T . From the decomposition, A -1 can be recovered thanks to the formula:

A -1 = A T L -T L -1 .

III. BASICS ON FIXED-POINT ARITHMETIC

This section presents some basics on fixed-point arithmetic. It explains the means to implement and bound the error of fixedpoint addition/subtraction, multiplication, and shift, the most useful operators in DSP algorithms. The reader is encouraged to refer to [START_REF] Yates | Fixed-Point Arithmetic: An Introduction[END_REF] or [START_REF] Mouilleron | Automated Synthesis of Target-Dependent Programs for Polynomial Evaluation in Fixed-Point Arithmetic[END_REF] for a more detailed introduction to fixedpoint arithmetic.

A. Fixed-point number representation

Fixed-point arithmetic allows to represent a real value by means of an integer associated to an implicit scaling factor. Let X be a k-bit signed integer in radix 2, encoded using two's complement notation. Combined with a factor f ∈ Z, it represents the real value x defined as follows:

x = X • 2 -f .
In the sequel of this article, Q i.f denotes the format of a given fixed-point variable v represented using a k-bit integer associated with a scaling factor f , with k = i + f . Here i and f denote the number of bits in the integer and fraction part of v, respectively, while k represents its wordlength. Hence v is such that:

v ∈ {V • 2 -f } with V ∈ Z ∩ [-2 k-1 , 2 k-1 -1]
by step of 2 -f .

B. Arithmetic and error model

Given an arithmetic operation ∈ {+, -, ×, <<, >>}, to bound the range of the resulting value and the rounding error entailed by its computation, we keep track for each fixed-point variable v, of the three following intervals: 1) Val(v) enclosing the values of v computed at run-time with finite precision, 2) Math(v) enclosing the values taken by v if computations were done with infinite precision, 3) Err(v) enclosing the rounding errors occurring while computing v, such that:

Err(v) = Math(v) -Val(v).
Notice that all the computations involving Val(v) and Err(v) are done using interval arithmetic [START_REF] Moore | Interval Analysis[END_REF]. And thanks to the formula above, keeping track of Val(v) and Err(v) suffices to deduce Math(v).

To show how Val(v) and Err(v) are computed, let us define three fixed-point variables v, v 1 , and v 2 along with their formats Q i.f , Q i 1 .f 1 , and Q i 2 .f 2 , respectively. When ∈ {+, -, ×} and v = v 1 v 2 , we have:

Val(v) = Val(v 1 ) Val(v 2 ) -Err .
For physical shifts, that is, for ∈ {<<, >>}, we simply have:

Val(v) = Val(v 1 ) -Err ,
since these operators do not modify the value of the input v 1 but both the bit string of V 1 and its fixed-point format. Interpreting v 1 in a different fixed-point format without modifying the bit string of V 1 is called virtual shift1 in [START_REF] Mouilleron | Automated Synthesis of Target-Dependent Programs for Polynomial Evaluation in Fixed-Point Arithmetic[END_REF]: this operation is error-free, and it is the only shift operation that modify the value of v 1 :

Val(v) = Val(v 1 ) × 2 f1-f .
Addition and subtraction. In absence of overflow, addition and subtraction are error-free. Hence for ∈ {+, -} we have:

Err(v) = Err(v 1 ) Err(v 2 ) and, i = max(i 1 , i 2 ) + 1 and f = max(f 1 , f 2 ).
Multiplication. If is a multiplication, we have:

Err(v) = Err × + Err(v 1 ) • Err(v 2 ) + Err(v 1 ) • Val(v 2 ) + Val(v 1 ) • Err(v 2 ),
where Err × is the error entailed by the multiplication itself.

Remark that exact fixed-point multiplication results in a number having a fraction of f 1 + f 2 bits. And usually in fixedpoint arithmetic the error Err × is due to the truncation of the exact result of the multiplication to fit in a smaller format with f < f 1 + f 2 fraction bits. Hence we have:

Err × = [0, 2 -f -2 -(f1+f2) ].
Most 32-bit DSP processors provide a 32 × 32 multiplier that returns the 32 most significant bits of the exact result, which is the multiplier considered in this work. In that case,

i = i 1 + i 2 and f = 32 -i.
Left and right shift. If ∈ {<<, >>}, we have:

Err(v) = Err(v 1 ) + Err .
Left shifts of s bits entail no error but only a possible overflow:

Err << = [0, 0] and (i, f ) = (i 1 -s, f 1 + s).
However right shifts of s bits may be followed by a truncation to fit the result in a smaller format with f < f 1 fraction bits. Thus, we have:

(i, f ) = (i 1 + s, f 1 -s) and Err >> = [0, 2 -f1+s -2 -f1 ].

IV. FIXED-POINT SQUARE ROOT AND DIVISION

Many research articles have been published that deal with the four fixed-point operations presented in Section III [START_REF] Mouilleron | Automated Synthesis of Target-Dependent Programs for Polynomial Evaluation in Fixed-Point Arithmetic[END_REF], [START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF], [START_REF] Lopez | Sum-of-products evaluation schemes with fixed-point arithmetic, and their application to IIR filter implementation[END_REF], [START_REF] Lee | Optimized Custom Precision Function Evaluation for Embedded Processors[END_REF]. Square root and division did not receive as much attention, and prior to this work, we are not aware of any published material that formalizes rigorously square root and division in fixed-point arithmetic, by suggesting an implementation together with an error model. Yet square root is useful when computing euclidean norms. In a linear algebra context, Cholesky decomposition and triangular inverse involve square roots and divisions, as shown in Section II. Division may also be used to evaluate rational polynomial approximants [START_REF] Cheung | Automating custom-precision function evaluation for embedded processors[END_REF]. The following of the section presents our new fixed-point square root and division operators.

A. Fixed-point square root

This first part presents the operator of fixed-point square root we have implemented. Assuming

v 1 ≥ 0, for v = √ v 1 , we have: Val(v) = Val(v 1 ) -Err √
since the computed value is truncated, while Err(v) is:

Err(v) = Math(v 1 ) -Val(v 1 ) + Err √ ,
where Err √ is the error entailed by the square root operation itself. The error term is given by the following formula:

Err(v) = Val(v 1 ) + Err(v 1 ) -Val(v 1 ) + Err √ (4) = Val(v 1 ) • 1 + Err(v 1 ) Val(v 1 ) -1 + Err √ .
Here the factorization is used to remedy the interval dependency phenomenon inherent to interval arithmetic.

Notice that this formula does not yield tight error bounds as soon as Val(v 1 ) smallest elements are of the same order of magnitude than Err(v 1 ). To overcome this issue, we may use the subadditivity property of the square root function, which holds as long as x and x + y are both positive:

√ x -|y| ≤ √ x + y ≤ √ x + |y|.
Hence we deduce the following bounds on Err(v):

Err √ -|Err(v 1 )| ≤ Err(v) ≤ Err √ + |Err(v 1 )|. (5)
In practice, we compute the intersection of the enclosures (4) and [START_REF] Yates | Fixed-Point Arithmetic: An Introduction[END_REF].

As for Err √ , it depends on the square root algorithm. To explicit such an algorithm, let us remember that we have

v 1 = V 1 •2 -f1 .
In the following, we propose two different algorithms together with a piece of C code to perform fixed-point square root, and we compare them by exhibiting their error bounds. Algorithm 1. A first attempt to implement square root is to start from the following rewriting:

√ v 1 = √ V 1 • 2 -f1/2 if f 1 is even, V 1 /2 • 2 -(f1-1)/2 if f 1 is odd.
An implementation of this approach would compute √ v 1 using one of the following:

V 1 • 2 -f1/2 or V 1 /2 • 2 -(f1-1)/2
where √ • is the integer square root operation. This operation may be implemented in hardware or in software using multiple techniques such as digit-recurrence, and Newton-Raphson or Goldschmidt iteration [START_REF] Ercegovac | Improving Goldschmidt division, square root, and square root reciprocal[END_REF], [START_REF] Ercegovac | Digital Arithmetic[END_REF]. We deduce that the value √ v 1 has ≈ f 1 /2 fraction bits, and that

Err √ = [0, 2 -f 1 2 ] or Err √ = [0, 2 -f 1 -1 2 ]
depending on the parity of f 1 . Compared to the error bound of the previous operations, an error bound of ≈ 2 -f1/2 for the square root is not acceptable in practice.

Algorithm 2. To overcome the accuracy issue of Algorithm 1 above, let use rewrite v 1 as

v 1 = 2 η • V 1 • 2 -(f1+η)
with the integer η being a parameter of the algorithm chosen at synthesis-time such as f 1 + η is even. Using this scaling factor, it follows that

√ v 1 = 2 η • V 1 • 2 -(f 1 +η) 2 . ( 6
)
The cases η = 0 and η = -1 correspond to the even and odd cases of Algorithm 1, respectively. An algorithm that exploits (6) shifts the integer representation V 1 of v 1 by η bits to the left and computes its integer square root. The result of this algorithm is a fixed-point variable with (f 1 + η)/2 bits of fraction part. Hence using this approach, we conclude that

i = i 1 /2 , f = f 1 + η 2 , and √ v 1 = 2 η • V 1 •2 -(f 1 +η) 2
, where √ 2 η • V 1 is computed using an integer square root operation. It is clear now that with such an algorithm, we obtain the following bound on Err √ :

Err √ = 0, 2 -(f 1 +η) 2 .
Notice that it would not make sense to choose η < 0, since it would result in an increase of Err √ . Hence in the following of the section, we assume η ≥ 0.

Listing 1 gives a C implementation for a 32 → 32bit square root operation using Algorithm 2 based on digitrecurrence iteration when 0 ≤ η < 32. This is by far the most frequent case in practice. In this code, the 32 least significant

uint32_t isqrt32hu(uint32_t V1, uint16_t eta) { uint64_t V1_extended = ((uint64_t)V1) << eta; // eta >= 0 uint64_t V = 0; int64_t one = 0x4000000000000000ll; // 2ˆ(62) while (one != 0){ if (V1_extended >= V + one){ V1_extended = V1_extended -(V + one); V = V + (one << 1); } V >>= 1; one >>= 2; } return (uint32_t)V; } Listing 1. C code of 32 → 32-bit fixed-point square root operation.
bits of the variable V are returned. Hence the parameter η must be carefully chosen, to ensure that f 1 + η is even, and that the wordlength of the result is at most 32, otherwise an overflow may occur.

B. Fixed-point division

This second part presents our fixed-point division operator. Let v = v 1 /v 2 . When the quotient is defined, i.e. when v 2 = 0, that is, 0 / ∈ Val(v 2 ), we have:

Val(v) = Val(v 1 ) Val(v 2 ) -Err /
while Err(v) is defined as:

Err(v) = Math(v 1 ) Math(v 2 ) - Val(v 1 ) Val(v 2 ) + Err / ,
where Err / is the error entailed by the division itself. It follows that the error term is defined as:

Err(v) = Val(v 2 ) • Err(v 1 ) -Val(v 1 ) • Err(v 2 ) Val(v 2 ) • (Val(v 2 ) + Err(v 2 )) + Err / .
From a theoretical point of view, the quotient v 1 /v 2 when defined must be a fixed-point variable in the format Q i.f with

i = i 1 + f 2 and f = f 1 + i 2 since:
1) The largest possible dividend is -2 i1-1 while the smallest divisor is 2 -f2 . Thus the quotient could be as large as -2 i1+f2-1 . 2) As for the opposite case, the smallest dividend is 2 -f1 while the largest divisor is -2 i2-1 . To be precise, the fractional part must be of size

f 1 + i 2 .
Remark that a special care must be taken when 0 ∈ Val(v 2 ) to avoid division by zero. In this case, we first compute both error and value bounds twice, using the two intervals Val(v 2 ) and Val(v 2 ), such that:

Val(v 2 ) ∪ Val(v 2 ) = Val(v 2 ) \ {0}.
Then we compute the union of the resulting intervals.

As for the square root, the fixed-point format of v depends on the algorithm implemented. In the following, we suggest two different algorithms to perform fixed-point division together with a piece of C code, and we analyze the difference between these algorithms by exhibiting their error bounds. To do this, let us remember that we have

v 1 = V 1 • 2 -f1 and v 2 = V 2 • 2 -f2 .
Algorithm 1. A first attempt to implement division is to start from the following rewriting:

v 1 v 2 = V 1 • 2 -f1 V 2 • 2 -f2 = V 1 V 2 • 2 -(f1-f2) .
This algorithm computes the quotient of the two integers V 1 and V 2 , and then considers f 1 -f 2 as the new implicit scale factor. The C standard requires integer division to be computed by discarding the fractional part of the exact division result, even for nonpositive results, that is, by rounding the exact result toward zero [22, § 6.5.5]. Since our implementations are intended to be compliant with the C language, we chose to implement a division operator satisfying this requirement. We denote by trunc(•) this operation. Hence:

v 1 v 2 = trunc V 1 V 2 • 2 -(f1-f2) .
Below we give an implementation of this method using the C standard integer division. If this option is not available or is too costly, this operation may also be implemented in hardware or in software using digit-recurrence, and Newton-Raphson or Goldschmidt iteration. In this case, the error Err / is as follows:

Err / = [-2 -(f1-f2) , 2 -(f1-f2) ].
It follows that the result has no more fraction bits, as soon as v 1 and v 2 have the same fixed-point format. This leads to a huge loss of accuracy since

Err / = [-1, 1].
Algorithm 2. To overcome the accuracy issue of Algorithm 1 above, we use a technique similar to the one presented for square root, by introducing a scaling factor η. Hence

v 1 v 2 = V 1 • 2 -f1 V 2 • 2 -f2 = V 1 • 2 η V 2 • 2 -(f1-f2+η) .
Usually, an implementation of this algorithm computes

v 1 v 2 = trunc V 1 • 2 η V 2 • 2 -(f1-f2+η)
which results in a variable having the fractional part

f = f 1 -f 2 + η.
In our context, we consider that all the variables have the same wordlength, particularly i + f = i 1 + f 1 . It follows that the integer part of the result is:

i = i 1 + f 2 -η. (7)
Then we deduce that the error Err / is as follows:

Err / = [-2 -(f1-f2+η) , 2 -(f1-f2+η) ].
Remark that even when v 1 and v 2 have the same format, Err / = [-2 -η , 2 -η ] remains tight as long as η is large enough. Also as for square root, it would not make sense to choose η < 0, since it would result in an increase of Err / . Again in the following of the section, we assume η ≥ 0. Listing 2 gives a C implementation for a 32 × 32 → 32-bit division operation using Algorithm 2 when 0 ≤ η < 64. This code returns the int32_t div32hs(int32_t V1, int32_t V2, uint16_t eta) { int64_t t1 = ((int64_t)V1) << min(32, eta); // eta >= 0 int64_t t2 = ((int64_t)V2); int64_t V = (t1 / t2) << max(0, eta-32); return (int32_t) V; } Listing 2. C code of 32 × 32 → 32-bit fixed-point division operation.

32 least significant bits of the variable V. It must be clear to the reader that here again the parameter η must be chosen carefully, since it greatly influences the result by impacting its integer part i. Indeed picking a large η leads to a smaller value i than the theoretical one and it minimizes the error bound and ensures more accuracy on the result. However, by doing so, we suppose that the result is not large enough. More precisely, this means that the largest values in magnitude eventually taken by the result are ignored and discarded. This approach is equivalent to discarding the smallest values eventually taken by the variable v 2 in Val(v 2 ), that is, the values around 0 in Val(v 2 ). This may also be taken in consideration when computing Val(v) and Err(v). In our experiments, we have implemented various heuristics to decide the parameter η.

V. SYNTHESIS TOOL FOR MATRIX INVERSION

In this section, we detail the synthesis tool we have implemented for the particular case of matrix inversion.

A. Cholesky decomposition based approach

As mentioned earlier, we choose to implement matrix inversion through Cholesky decomposition. Given an input matrix A, it is done in three steps: 1) matrix decomposition: L such as A = LL T , 2) triangular matrix inversion: L -1 , 3) and matrix multiplication:

A -1 = L -T L -1 .
In fixed-point arithmetic, we can frequently reduce the ranges of the inputs to a range included in [-1, 1]. Indeed, for Cholesky decomposition, if a matrix A does not satisfy this condition, instead of decomposing A, we can decompose

B = 2 -k • A with k ∈ Z, (8)
to obtain B = RR T , where B is a matrix with coefficients in

[-1, 1]. It follows that A = 2 k • RR T = LL T where L = 2 k/2 • R.
Notice that k must be chosen even. For triangular matrix inversion, we can still compute A -1 as follows:

A -1 = 2 -k • B -1
, where B is as in [START_REF] Sung | Simulation-based word-length optimization method for fixed-point digital signal processing systems[END_REF].

In fixed-point arithmetic, these scalings are just a matter of manipulating the fixed-point format of the coefficients.

We have automated this in the FPLA tool2 which relies on the CGPE library. 3 These tools are detailed below. [START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF], CGPE is dedicated to the automated synthesis of fast and accurate fixed-point codes to evaluate mathematical expressions, like polynomials, dot products, and summations. It takes as input intervals that bound the range of each variable that appears in the expression, and produces a code which is able to evaluate the input expression for any instance of its arguments as long as they are in the original interval of values. Initially devoted to evaluating expressions involving only addition/subtraction, multiplication, and shift, CGPE has been here enhanced by adding the square root and division operators presented in Section IV.

CGPE (Code Generation for Polynomial Evaluation). Developed by Revy and Mouilleron

FPLA (Fixed-Point Linear Algebra). We developed this tool with the aim of generating fixed-point code for the most frequently used linear algebra routines. It handles the aspects peculiar to each class of input problems and relies on the CGPE library for the low level and code synthesis details. For instance, when prompted to generate code for matrix multiplication, it is capable of finding good trade-offs between code size and numerical accuracy as in [START_REF] Martel | Code Size and Accuracy-Aware Synthesis of Fixed-Point Programs for Matrix Multiplication[END_REF], by using the same code for output coefficients obtained from input rows and columns with elements having close enough fixed-point formats. For triangular matrix inversion and Cholesky decomposition, since some of the output coefficients are dependant on previously generated code, this strategy is harder to apply. Hence in these cases, FPLA takes as input the intervals that bound the ranges of the coefficients of the input matrix and works by making successive calls to CGPE with the appropriate parameters. It iterates on the output matrix and asks for a code to compute a given output coefficient, defined by a specific expression as in [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] or [START_REF] Chen | Layered turbo space-time coded mimo-ofdm systems for time varying channels[END_REF]. Once this code generated by CGPE, FPLA stores the range and error of its output. Note that FPLA also takes care of correctly ordering the calls to CGPE in such a way that each coefficient's code is generated only after all the information on which it depends has been collected. In this sense, FPLA's mode of operation greatly depends on the linear algebra problem to be solved.

B. Order of code synthesis in FPLA

In triangular inversion, the diagonal elements do not depend on any generated code. Therefore they may be computed in any order. The non diagonal coefficients depend only on the coefficients that precede them on the same column. FPLA must therefore follow a row major, column major or even a diagonal major approach. The latter consists in generating the elements on the diagonal, followed by those on the first sub-diagonal and so on. The last code generated in this fashion would be that of the bottom left coefficient n-1,0 . For Cholesky decomposition, a diagonal element i,i depends on the generated coefficients that precede it on row i. A non diagonal element i,j depends on the first j elements of row i as well as the first j + 1 elements of row j. FPLA may satisfy these dependencies by following either a row major or column major synthesis strategy but not a diagonal major strategy.

Listing 3 gives the global code produced by FPLA for a size-3 triangular matrix inversion. The coefficients of the upper triangular part are explicitly set to zero. Then compute_i_i computes the coefficient n i,i as n i,i = 1/d while compute_i_j computes the coefficient n i,j as

n i,j = (a 0 • b 0 + • • • + a i-j • b i-j )/d where a 0 =A[i][j], b 0 =N[j][j], a i-j =A[i][i-1], b i-j =N[i-1][j], and d =A[i][i].
The six compute_x_y functions of this example are each generated by CGPE in a specific C file. Listing 4 shows the code of compute_2_0 in Listing 3.

// A: input matrix --N: inverse matrix of A N[0][0] = compute_0_0( A[0][0] ); N[0][1] = 0; N[0][2] = 0; N[1][0] = -compute_1_0( A[1][0], N[0][0], A[1][1] ); N[1][1] = compute_1_1( A[1][1] ); N[1][2] = 0; N[2][0] = -compute_2_0( A[2][0], A[2][1], N[0][0], N[1][0], A[2][2] ); N[2][1] = -compute_2_1( A[2][1], N[1][1], A[2][2] ); N[2][2] = compute_2_2( A[2][2] );
Listing 3. FPLA ouput code for a 3 × 3 triangular matrix inversion.

C. How to use correctly fixed-point division?

Division in Section IV is the trickiest among the arithmetic operators treated. It requires one to provide a parameter η that impacts the format of the output. In practice, instead of choosing the parameter η, we compute it according to an expected acceptable output integer part using [START_REF] Nikolic | Design and Implementation of Numerical Linear Algebra Algorithms on Fixed-Point DSPs[END_REF]. This integer part can be set using multiple ways: 1) Set to a constant, 2) Using a function of the formats of the operands. For instance, if we want division results to have an integer part two bits larger than the integer part of its left operand denoted i 1 , we will use the function f (i 1 ) = i 1 + 2 to compute i. In this example, using [START_REF] Nikolic | Design and Implementation of Numerical Linear Algebra Algorithms on Fixed-Point DSPs[END_REF], we deduce that the parameter η = f 2 -2.

At first sight, the first solution seems to be either too restrictive or too unsafe. Indeed with long chains of computations, the format of the intermediate operands tend to grow and choosing a small enough output integer part is a good idea to bring the result coefficient to a manageable range. Conversely choosing a small output integer part increases the chance that overflows occur at run-time. In Section VI, we illustrate the interest of using the first solution on some cases, and we show experimental evidence of the problems caused by the different methods of deciding η.

VI. EXPERIMENTAL RESULTS

In this section, we conduct some experiments to show the interest of using the operators specified in Section IV as well as the Cholesky decomposition and triangular inversion basic blocks of Section V. First we investigate the impact of the output format of division. Second we study the speed of the generation and the sharpness of the error bounds of the result. Finally we show the impact of the matrix condition number on the accuracy of the generated code.

A. Impact of the output format of division on accuracy

As mentioned in Section V-C, we need to explicitly fix the output format of each division. More particularly, we define each output integer part using a function of the formats of the operands, the output fraction part being determined so as each result fits on 32 bits. This first experiment shows the impact of this function on the output accuracy. To do so, we defined four functions:

f1(i1, i2) = t, f2(i1, i2) = min(i1, i2) + t, f3(i1, i2) = max(i1, i2) + t, and f4(i1, i2) = (i1 + i2)/2 + t,
where t ∈ Z is a user defined parameter, and i 1 and i 2 are the integer parts of the numerator and denominator, respectively. The function f 1 consists in fixing all the division results to the same fixed-point format. The experiment consists in computing the Cholesky decomposition and the triangular inversion of matrices of size 5 and 10, respectively. Using FPLA, we synthesize codes for each problem and each function in {f 1 , f 2 , f 3 , f 4 }, for t ranging from -2 to 8. Then for each synthesized code, 10000 example instances are generated and solved both in fixed and floating-point arithmetics. Each example input is a matrix having 32-bits coefficients ranging from -1 to 1. Then the error considered is obtained by comparing the results to floating-point computations and by considering the maximum errors among the 10000 samples. The results are shown in Figure 1. No curve in the figures means that all of the examples overflew, and none succeeded. Obviously, we can observe that function used to determine the output format of division has a great impact on the accuracy of the generated code. For example, if we consider the case t = 0 on 5 × 5 Cholesky decomposition on Figure1(a), using f 1 leads to an error of ≈ 2 -28 , while using f 3 gives an error ≈ 2 -15 , that is, twice larger than f 1 . More particularly, we can observe that a good function choice is one that minimizes the output integer part. Indeed, as long as t ≥ -1, using the function f 1 always leads to better maximum error than using function f 3 . In addition surprisingly, as long as t ≥ -1, the function that gives the best results is f 1 (i 1 , i 2 ) = t, namely the function that fixes explicitly all the division results of a resulting code to the same fixed-point format independently of the input formats. Indeed the problem of using a function that depends on the input formats comes from the fact that it quickly leads to a growth of the integer part of each computed coefficient, since it relies on the previously computed coefficient themselves. Hence the interest of f 1 is that it avoids this fast growth, and leads to result coefficients having a fixed and relatively small integer part, thus to tighter errors than the other functions. This remark is particularly true on 5 × 5 Cholesky decomposition on Figure1(a) when t ≥ 0 where f 1 is the only function that leads to successful results. This phenomenon becomes obvious as the matrix size increases.

However we cannot restrain ourselves to function f 1 when implementing matrix inversion basic blocks. Indeed, cases occur where f 1 leads to unsuccessful results. This occurs when the result of some division has a fixed-point format with an integer part greater than t. This is the case on 10×10 triangular inversion on Figure1(b), where the diagonal coefficients of the inverse are 1/a i,j . Since a i,j may be arbitrarily small, then 1/a i,j may be too large to fit into a format with an integer part of t bit. In these cases, other functions should be preferred.

B. Sharpness of the error bounds and generation time

The originality of our approach is the automatic generation of certified error bounds along with the synthesized code. This enables the generated code to be used in critical and precision sensitive applications. However, it is equally important that these bounds be sharp, at least for a large class of matrices. To investigate their sharpness, we compare in this second experiment the error bounds for the case of triangular matrix inversion with the experimental errors obtained from inverting 10000 sample matrices. This experiment is carried out using the f 4 function introduced in the previous experiment with t = 1. For each matrix size from 4 to 40, C code and error bounds are obtained by running FPLA. The generation time for each size is shown by the third curve with triangle-shaped dots and the right ordinate of Figure 2. It does not exceed 14 seconds on an Intel Core i7-870 2.93 GHz and for 40 × 40 matrices, an improvement of several orders of magnitude over a hand written fixed-point code.

The remaining two curves show the evolution of the error bounds and the experimental errors. The bounds vary from 2 -26 to 2 2 while the experimental errors vary from 2 -28 to 2 -6 . The difference between the error bounds and experimental errors is less than 2 bits for size-4 matrices and is inferior to 5 bits for size-15 matrices, and it grows as the size of the input matrices grows. Although the bounds obtained for matrices of size larger than 35 are too large to be useful in practice, the experimental errors are still tight enough and do not exceed 2 -6 . These issues may be tackled by considering other means to handle division that are more suited to large matrices. Indeed, our experiments tend to show that the output format of division impacts heavily the accuracy of the result and that there is no way to determine this format that is adapted to all matrix sizes. We also argue that a bound of 2 -12 on the inversion of size-20 matrices is satisfying for a broad range of applications, and this is a large improvement over hand-written fixed-point codes or codes whose numerical quality is asserted uniquely by simulations and a posteriori processes.

C. Impact of the matrix condition number on accuracy

The sample matrices considered in the previous experiments were randomly drawn in the input intervals. In this third experiment, we consider some rather known matrices namely, Hilbert, Cauchy, Kms, Lehmer, and Prolate matrices. Among these, Hilbert and Cauchy matrices and to a lower extent Prolate are ill conditioned. Nonetheless, with a fixed-point code generated for matrices in the input format Q 1.31 , we were able to check that the fixed-point results, whenever computable, are quite accurate as shown in Figure 3. For sizes larger than 8 and 9, respectively, overflows occur when computing the decompositions of Cauchy and Hilbert matrices. But this fact does not invalidate our approach. Indeed, these matrices are very ill conditioned and are difficult to decompose accurately even in floating-point arithmetic.

On the other hand, Lehmer and Kms matrices have a linearly growing condition number and are therefore very well suited to our approach. Indeed as shown by the two bottom curves, the code generated by FPLA decomposes these matrices with a precision of up to 25 bits.

VII. CONCLUSION

In this article, we presented an automated approach to help in writing codes in fixed-point arithmetic for the particular case of matrix inversion based on Cholesky decomposition. First we formalized rigorously the fixed-point square root and division, especially in of error bounds, and implemented both of them in the CGPE tool. Second we introduced FPLA, dedicated to the synthesis of linear algebra subroutines. Our work focused particularly on the correct handling of fixedpoint division, which is not as straightforward as the other operators and for which the output format must be carefully chosen to ensure a certain accuracy on the result. Enhancing our tool chain with these operators allowed us to tackle Cholesky decomposition and triangular matrix inversion. We finally showed that accurate fixed-point codes accompanied by bounds on the rounding errors can be automatically generated in a few seconds to invert and to decompose matrices of sizes up to 40.

In addition our further research direction is twofold: As a first direction, we intend to investigate the synthesis of other linear algebra basic blocks, like LU or QR decompositions, that may also be used to implement matrix inversion, and to implement them in FPLA. As a second direction, we aim at investigating, similarly to the work done in [START_REF] Martel | Code Size and Accuracy-Aware Synthesis of Fixed-Point Programs for Matrix Multiplication[END_REF], the different trade-offs involved in the code synthesis process and especially the one between code size and accuracy.

int32_t compute_2_0 (

 compute_2_0 int32_t a0 / * Q1.31 in [-1,1] * / , int32_t a1 / * Q2.30 in [-2,2] * / , int32_t b0 / * Q1.31 in [-1,1] * / , int32_t b1 / * Q2.30 in [-2,2] * / , int32_t d / * Q1.31 in [0.88,0.99] * / ){ int32_t r0 = mul(a0, b0); // Q2.30 in [-1,1] int32_t r1 = r0 >> 2; // Q4.28 in [-1,1] int32_t r2 = mul(a1, b1); // Q4.28 in [-4,4] int32_t r3 = r1 + r2; // Q4.28 in [-5,5] int32_t r4 = div32hs(r3, d, 31); // Q4.28 in [-8,8] return r4; } Listing 4. C code of the compute_2_0 function.

  Triangular 10 × 10.

Figure 1 .

 1 Figure 1. Maximum error of Cholesky decomposition and triangular inverse with various functions used to determine output formats of division.

Figure 2 .

 2 Figure 2. Comparison of the error bounds and experimental errors together with generation time, for the inversion of triangular matrices of size 4 to 40.

Figure 3 .

 3 Figure 3. Maximum errors measured when computing the Cholesky decomposition of various kinds of matrices for sizes varying from 4 to 14.

This operation is not detailed in our article, but we encourage the reader to refer to[START_REF] Mouilleron | Automated Synthesis of Target-Dependent Programs for Polynomial Evaluation in Fixed-Point Arithmetic[END_REF] for more details.

An archive containing the code of the FPLA tool is available upon request.

CGPE is freely available at http://cgpe.gforge.inria.fr/.


Ahmad Sinno

Pascal Ruaux

Luc Chassagne

Suat Topçu

Yasser Alayli

Gilles Lerondel

Pascal Royer

Aurélien Bruyant

Sylvain Blaise

Large area sample holder unit for enhanced near field microscopy applications

Keywords: Dimensional nanometrology, long range of travel, laser interferometer, closed loop control, near-field microscopy sample-holder, high resolution displacement system

Near-field microscopy applications sometimes require scanning the sample over the millimetre range. This paper describes a home-made displacement system that can replace the limited range XY scanner of a classical near-field microscope. The system lays on a coarse/fine displacements configuration to attain the millimetre. Nanometre scale repeatability and resolution are met thanks to a control loop based on phase shifting techniques. The system was tested under an Atomic Force Microscope. The setup of the experiments and the results obtained are herein presented.

INTRODUCTION

Mapping the properties of a macro system to the nanometre world is currently the area of study for a lot of industries. Nevertheless, many research institutes are seeking to push the limits of the nanometre to a larger scope [START_REF] St | Nanomanipulation in a scanning electron microscope[END_REF]. Imaging small objects at a large scale is one of the obstacles nanometrology faces today. State of the art Scanning Probe Microscopes (SPMs) have an imaging area of about a hundred micrometres. This is mainly due to the limited range integrated scanner. Piezoelectric actuators, in addition to presenting nonlinearity defects, do not allow millimetre scans. Furthermore, commercial displacement systems that attain the millimetre do not achieve the desired nanometre repeatability and accuracy. Yet, many applications require scanning the sample over the millimetre range: Near field characterization of optical integrated structures, electronic lithography for designing highly resolved patterns, and microelectronics for example. Scanning optical integrated structures with a combined Atomic Force Microscope (AFM) and a Scanning Near Field Optical Microscope (SNOM) may reveal internal losses of the light that is propagating along the waveguides of the structure. Thus investigating a waveguide along its light propagation axis turns out to be necessary [START_REF] Bachelot | Probing Photonic and Optoelectronic Structures by Apertureless Scanning Near-Field Optical Microscopy[END_REF][START_REF] Blaise | Apertureless scanning near-field optical microscopy for ion exchange channel waveguide characterization[END_REF]. Such a waveguide sometimes spreads over a millimetre and even centimetre long. Electronic lithography using a Scanning Electron Microscope (SEM) often shows degraded resolutions due to the positioning system that does not manage to accurately stitch the sample [START_REF] Rai-Choudhury | SPIE Handbook of Microlithography, Micromachining and Microfabrication[END_REF]. This causes the operator to widen the pattern after the sample is moved. High resolution versus a reduced image size is what traditional Atomic Force Microscopes (AFMs) offer. The number of pixels or data points that can be acquired per a line of scan is usually limited to 1024 whatever the size of the image is [START_REF]Nano-R atomic force microscope specifications[END_REF]. If a full range image is needed, the operator will find the resolution reduced. For a 100 µm square scan for example, the resolution is found to be around 100 nm square per pixel. An alternative solution that allows to overcome this limitation and even to cross the 100 µm range barrier is stitching [START_REF] Blaise | Apertureless scanning near-field optical microscopy for ion exchange channel waveguide characterization[END_REF]. Stitching is a time consuming technique during which the surface of the sample is partitioned into small scan areas imaged each apart and then stitched together. These elementary scans called windows must overlap to allow accurately match the overlapping sections afterwards. Therefore, a sample surface that is not perfect is needed in order to facilitate the matching process by means of image recognition techniques. This makes the stitching technique not always applicable: firstly when the surface of the sample does not show distinguishable topography patterns, and secondly when the measurement parameters such as temperature are time dependent. Moreover, traceability in addition to accuracy is sometimes hard to ensure. Accreditation institutes are currently developing displacement and measurement systems. This is the case of Metrological Atomic Force Microscopes (MAFM) used as reference machines for characterizing and standardizing samples [START_REF] Lee | An ultraprecision stage for alignment of wafers in advanced microlithography[END_REF][START_REF] Schneir | Design of an atomic force microscope with interferometric position control[END_REF][START_REF] Bienias | A metrological scanning force microscope used for coating thickness and other topographical measurements[END_REF][START_REF] Picotto | A sample scanning system with nanometric accuracy for quantitative SPM measurements[END_REF][START_REF] Gonda | Real-time, interferometrically measuring atomic force microscope for direct calibration of standard[END_REF][START_REF] Haycocks | A metrological atomic force microscope for calibration of transfer standards[END_REF][START_REF] Koops | Nanometrology standards in the Netherlands: the traceable scanning probe microscope[END_REF]. *ahmad.sinno@lisv.uvsq.fr; phone 01 39253023; fax 01 39253025; www.lisv.uvsq.fr These often suffer bad repeatability and accuracy due to the defects of their displacement system. Moreover, Coordinate Measuring Machines (CMM) are intended to make measurements over a large scale. However, the highest resolution ever reached for a 1 mm x 100 µm scan is 100 nm [START_REF] Wilkening | Nanoscale Calibration Standards and Methods[END_REF].

We present a 2D displacement system with a 2-stage configuration to attain the millimetre range. The first stage performs coarse displacements while the second rectifies the defects of the first stage. The system is used as a sample holder as it was integrated into an Atomic Force Microscope to allow scans beyond the microscope's scanning area limit.

EXPERIMENTAL SETUP

Sketch of the system

The system developed is a Sample Holder Unit (SHU) that can be integrated to most of current commercial near-field microscopes. To achieve this, the SHU performs 2D lateral displacements relative to a fixed microscope's tip while carrying a sample (Fig. 1). Scans are carried out on a coarse and fine displacements basis. Two stages are used for each displacement axis. As a coarse stage, the SHU integrates two millimetre range motorized tables, Aerotech ANT50 [14], that allow a resolution of 10 nm and speeds up to 2 mm/s over 50 mm with a straightness of 1 µm. The fine stage comprises a 2D piezoelectric actuator (PZT), PiezoSystem Jena, used for accurate positioning thanks to its subnanometric resolution and high bandwidth. Two laser interferometers, Zygo ZMI2001 [15], permit real time position sensing with a resolution of 0.31 nm. The laser source is a frequency stabilized HeNe laser [15] that ensures measurements traceable to the Metre. A Zerodur mirror cube acting as a moving mirror for the laser interferometers is placed on top of the displacement stages. This cube carries the sample to scan on its top face. A home made high frequency electronic board ensures the nanometre performances [START_REF] Topçu | Heterodyne interferometric technique for displacement control at sub-nanometric scale[END_REF]. The board can in fact order PZT displacement steps resolved at 0.26 nm [START_REF] Chassagne | Highly accurate positioning control method for piezoelectric actuators based on phase-shifting optoelectronics[END_REF]. It also allows locking the system in its position thanks to a servo loop controller associated to a Phase Locked Loop (PLL) scheme. The electronic board uses phase shifting techniques which allow nanometric resolution and repeatability. The integrity of the system is controlled through a Labview interface.

AFM integration

To allow topography imaging, the SHU was integrated to an Atomic Force Microscope (AFM). For this, an AFM (APHB-0100 Park Scientific Instruments) was placed above the SHU in such a way the tip of the microscope can sense the sample surface (Fig. 2). The XY PZT actuators of the microscope were inhibited. The SHU ensures the lateral displacements on a line per line basis. The SHU does not integrate z-axis control. It is the AFM that handles the force feedback loop for z-position sensing and imaging.

The microscope presents a major limitation: the number of data points acquired at each line of scan is limited to 1024. This yields to an image of 100 µm square resolved at 100 nm per pixel. The SHU operates over a millimetre range. A resulting resolution of 1 µm per pixel would then not be suitable. For this, we shunted the acquisition module of the AFM and a Labview Real Time based data acquisition board (NI PCI 6221) took hand over the acquisition. The aim is to have a resolution that is not correlated to the size of the image. The 16-bit registers of the acquisition board and its high sampling frequency (up to 100 KHz) agree with the desired resolution and ensure handling the high flow of data when the SHU runs at high speeds. The board acquires a signal that is proportional to the elongation of the Z piezoelectric actuator, that is, the topography signal. The whole is placed on an active anti-vibration table thus decoupling external noise.

Principle of the control loop

The displacements of the sample holder are based on double stage actuators: the linear motors and the piezoelectric actuators. These stages cooperate together to explore the surface of the sample. In fact, it is the PZT stage that rectifies the defects encountered by the motors along the path. The motors have mainly straightness defects. However, the PZT itself presents some limitations due to intrinsic defects such as hysteresis. For this, an upper control layer was introduced to help the PZT actuator track the desired path. The interferometers observe the position along both X and Y axes.

The control loop is based on a Phase Locked Loop (PLL) scheme. The interferometers deliver at their output a phase modulated position signal. We designed a high frequency optoelectronic board that controls the phase of the position signal with respect to a reference signal in order to control the position [START_REF] Chassagne | Highly accurate positioning control method for piezoelectric actuators based on phase-shifting optoelectronics[END_REF]. If the PZT actuators are to be moved, the position signal is phase-shifted with respect to the reference signal. If the PZT actuators are to be held at a specific position, the phase of the position signal is locked on that of the reference signal. This phase locking based principle ensures repeatable and highly resolved displacements. The optoelectronic board uses ultra-high stability and high frequency integrated circuits. The position resolution can be programmed down to 0.26 nm.

Contrarily to classical model-based control schemes [START_REF] Dirscherl | online correction of scanning probe microscopes[END_REF], this PLL control loop tries to track the desired path by locking to a reference signal. No PZT model is needed. This has the advantage that the control loop bypasses any change that arises on the PZT behaviour. Consequently, the PZT actuator does not have to be calibrated previous to operation.

Another advantage is that the PLL is hardware implemented in the optoelectronic board. This may be advantageous when the SHU runs at a high speed; the bandwidth lays on the effectiveness of the servo loop controller that is implemented. A software control loop would require some processing time which could be limiting.

Mechanical adjustments

The SHU ensures long range displacements by stacking two stages: the coarse stage and the fine stage. Above these, the Zerodur mirror cube stands as a moving mirror for the laser interferometers. These three stages must have their axes aligned (Fig. 3); measurement errors might be encountered otherwise. These axes are the axes of the linear motors (Xm,Ym), the axes of the piezoelectric actuators (Xp,Yp) and the axes of the optical beams of the laser interferometers (Xo,Yo). For this to be achieved, three tilt adjustment elements are inserted (Fig. 4) in order to ensure three degrees of freedom (3-DOF) according to the yaw, pitch and roll angles.

Furthermore as a 2D system, the SHU requires its axes to be orthogonal [START_REF] Bobroff | Critical alignment in plane mirror interferometry[END_REF]. For this, the axes of the motors were aligned with an uncertainty of 5 arcsec. The mirror cube has its faces orthogonal at 1 arcsec. All these alignments are necessary before operating the SHU. The SHU is in fact operating over a range of one millimetre which magnifies any residual mechanical misalignment. Any remaining misalignment is processed online using the closed loop control as long as the misalignment remains below the range of the PZT actuator.

  Y p X p  Z D      Y m X m Y m Y p Y o X m X p X o L in
Another constraint to deal with is the parallelism between the upper surface of the sample and the displacement plane of the stages. Such a tilt misleads to an artefact on the resulting image giving the impression that the upper face of the sample is tilted with respect to its bottom face. Moreover, this tilt will make the Z-PZT actuator of the AFM contract/retract to track the surface of the sample. When scanning over a millimetre range, the Z actuator might not have enough range to recover the tilt, which might cause the AFM tip to unlock. Fig. 4. Photography of the displacement system. It shows (bottom to up): the X and Y long range linear motors, the pitch and roll adjustment elements, the yaw adjustment element, the XY PZT actuator, and the mirror cube that moves with respect to the reference mirrors of the interferometers. Note also the head of the microscope that holds the tip and that is kept close to the top face of the mirror cube where the sample is placed.

Environment supervision

Besides mechanical adjustments that are critical in a multi-stage millimetre range scanner, care should be taken for the surrounding environment. Thermal drift will act to expand the metrological chain of the system. This chain includes the body of the microscope in addition to the multi-stacked elements of the SHU. A thermal expansion of any of the mechanical parts will be rectified by the servo control loops, leading to an artefact. For this, the mechanical parts of the system are made of low expansion coefficient materials. In addition to this, temperature of the room is controlled, and the system is turned on several hours previous to operation for heat up. Temperature control is also required to prevent any drift on the laser measurements [START_REF] Demarest | High resolution, high speed, low data age uncertainty, heterodyne displacement measuring interferometer electronics[END_REF].

Furthermore, the interferometers used feature immunity to refractive index gradients due to their Differential Plane Mirror Interferometer [21] configuration in which both arms of the interferometer are aligned. In addition to this, a small distance (1 cm) separates the optical beams of each of the interferometers. The reference mirrors and the moving mirror cube are kept close to each other to reduce any effect of the surrounding environment.

As for vibrations, the SHU and the microscope lay on an active anti-vibration table that attenuates external noise.

MEASUREMENTS AND DISCUSSIONS

Parameters of the experiment

We chose to image a highly integrated silicon-on-insulator (SOI) substrate. The sample features waveguides on its surface. The waveguides are parallel but irregularly distanced. Their height is around 160 nm and their width 2.8 µm.

The sample was placed on the top face of the mirror cube, and then the tip of the AFM approached. The AFM was operating in tapping mode. A speed of scan of 40 µm.s -1 was chosen. The SHU allows higher speeds but the microscope's operation regulations prohibit excessive speeds. The temperature of the room was controlled to prevent the limited range Z scanner (around 7 µm) from saturating. These 7 µm must be enough to handle thermal drifts over time, pitch and roll misalignments over the millimetre scanning range, in addition to abrupt tip overshoots in case the tip encounters high topography patterns. Instead of normal operation of the microscope where the Z actuator is centred, we chose to bias the actuator so it can operate over its full range (from -3.5 to +3.5 µm) rather half its range.

Results for long topography imaging

The image presented in Fig. 5 is 1 mm long and 725 nm wide. The long scan axis is labelled Y while X is used for the short scan axis. So the SHU was travelling perpendicularly to the waveguides along the Y axis, and then changing lines. The resolutions used are 20 nm along Y and 25 nm along X, so 29 lines were achieved. The number of data points acquired is then 29 times 50000. At a speed of 40 µm.s-1, each line took 25 seconds to be scanned; the whole image took around 25 minutes due to the repositioning process when the SHU was changing lines.

Figure 5 is a raw image displayed using Matlab. It presents progressive shading off colours meaning that it is not perfectly horizontal. The image spreads from negative values of the topography signal (due to biasing the actuator) till positive values. The image is in fact tilted due to different probable factors. Firstly, although regulated, the temperature may have altered which may have caused a thermal drift on the resulting image due to the long measurement duration of 25 minutes. This tilt appears along the short scan axis X. Secondly, the pitch and roll angles were adjusted at best in order to prevent the Z scanner from saturating along a tilted surface. However, we had no means to measure the absolute tilt with respect to a zero reference. This tilt appears along the long scan axis Y. Thirdly, this tilt may also be due to the sample itself which may have been tilted during fabrication (sample's upper and lower faces not parallel).

Fig. 5. 3D Atomic Force Microscope topography of a substrate comprising on its surface a set of parallel but unequally distanced waveguides. The waveguides spread along the X axis. The image is 1 mm long and 750 nm wide. Shading off colours reflect the fact that the image is tilted. This tilt is greater along X than it is along Y. This is because the SHU spent more time to travel along X than it did along Y (temperature effect).

For an in-depth view, a cross section is pointed out on Fig. 6. This figure shows a line travelled by the SHU along the long scan axis Y, that is, temperature had a minor effect on the tilt of this figure. It is the misalignment errors (pitch misalignment) that induced this tilt of around 1.5 µm. Fig. 6. Cross section of the previous image showing one of its lines along Y. A tilt appears meaning that the surface of the sample was not perfectly parallel to the plane of motion of the SHU. This made the piezoelectric actuator of the microscope contract of about 1500 nm when the SHU was moving forward (positive Y direction), and then retract of the same amount during the backward motion. A larger tilt or a longer move along Y would cause the tip of the AFM to unlock if the tilt exceeds the range of the Z actuator. Fig. 7 shows the rectified cross section. This figure reveals more information: it presents an un-straight profile which means that the surface of the sample is wavy. On the hand, some peaks appear on the lower right part of each waveguide. These are overshoots that may be due to abrupt tip-sample interactions which may themselves be due to scanning at a relatively high lateral speed of 40 µm.s -1 . Overshoots appear at the right bottom of each waveguide, and may be due to abrupt tip-sample interactions. On the other hand, the cross section has a non straight profile meaning that the surface of the sample is wavy.

Fig. 8 shows a zoom of a single waveguide around the Y position 217 µm. The figure shows that the resolution is far higher from the dimensions of the waveguide. The dimensions of the waveguide can be measured on this cross section. The 157 nm height and 2.85 µm width are compatible with those expected.

Fig. 1 .

 1 Fig. 1. Bloc diagram of the system. It shows the displacement system serving as a 2D scanner under an Atomic Force Microscope (AFM). The displacement system comprises 2 linear motors to execute coarse moves, and a 2D piezoelectric actuator for precise positioning. Two laser interferometers monitor the X and Y positions which are fed back to a computer for control. The computer processes the information and sends the necessary control commands to a homemade optoelectronic board able to accurately control the position of the XY PZT actuators, and to the controller of the linear motors for long range motion control. Simultaneously, a Motion Synchronized Pulse Firing Controller (MSPFC) senses the displacements. It allows acquiring the AFM signal (topography signal) synchronously with the lateral displacements using a Data Acquisition Board (DAB).

  The acquisition of the topography signal is synchronized with the displacements of the SHU. A Motion Synchronized Pulse Firing Controller (MSPFC) ensures this job. Each time the SHU steps, the MSPFC triggers an acquisition point. The MSPFC sampling signals are derived from signals of the optical encoders of the coarse stage and internal signals of the PLL control board.

Fig 2 .

 2 Fig 2. Top view photography showing the displacement and imaging systems: we can see the mirror cube which ispositioned on top of the displacement system, the AFM with its head that was turned, and the X and Y reference mirrors. The sample to scan is placed on the top of the mirror cube allowing the tip of the microscope sense its surface. The whole is placed on an active anti-vibration table thus decoupling external noise.

  ie z o e le c tr ic a c tu a to r Z e r o d u r c u b e w ith m ir o r fa c e s R o ta tin g s ta g e w ith p ie z o e le c tr ic c o n tr o l

Fig. 3 .

 3 Fig. 3. (Left) this figure shows the Computer Aided Design of the SHU and highlights its degrees of freedom. The linear motors' axes (Xm,Ym), the 2D piezoelectric actuator's axes (Xp,Yp), and the optical axes (Xo,Yo) are to be aligned (forcing α=β=0). This can be achieved thanks to tilt and rotation adjustment elements. These allow adjusting the pitch, roll and yaw angles denoted θ 1 , θ 2 and θ 3 . (Right) ZD stands for the Zerodur mirror cube.

Fig. 7 .

 7 Fig. 7. Rectified cross section along the Y axis; the tilt was subtracted for better understanding the topography encountered by the SHU as it travelled along the long scan axis. Note the different scales used for the axes of this figure. This makes the waveguides seem distorted and appear that high. The figure reveals two properties: overshoots and waviness.Overshoots appear at the right bottom of each waveguide, and may be due to abrupt tip-sample interactions. On the other hand, the cross section has a non straight profile meaning that the surface of the sample is wavy.

ACKNOWLEDGMENTS

The authors are grateful towards Mr. E. Oseret from the Laboratoire de Parallélisme, Réseaux, Systèmes, et Modélisation (PRiSM) for his useful discussions about handling high density images, and would like to thank Mr. P. Juncar from the LNME (Laboratoire National de métrologie et d'Essais) for his collaboration in dimensional metrology. This project is co-funded by the French Government (ACI Nanosciences) and by the Ile-de-France region.

These consecutive zooms show one of the advantages of the SHU: no need to scan once again in order to zoom; the entire information can be retrieved from the original image. With a classical AFM, it would have been necessary to scan once again for a zoom. Moreover, the original image and its zooms have the same resolution because the resolution is dissociated from the size of the image. Large images can be achieved with resolutions not accessible with a classical microscope.

Discussions

We reported the results of a millimetre long scan achieved under an Atomic Force Microscope. The width of the image is 725 nm, limited by the range of the AFM actuator (Z axis actuator). In fact, the image presents a tilt that occupies the full range of the actuator. For wider scans, this tilt has to be better controlled. This can be done by better controlling the temperature in addition to finely aligning the surface of the sample with respect to the moving stages. Moreover, temperature control would offer more accurate measurements. Furthermore, the SHU does not integrate Z displacements. This causes the measurements to be limited by the range of the Z actuator of the AFM. A Z control shall help the SHU adapt to large altitude variations.

The SHU was integrated in this experimentation to an AFM but it may be integrated to other Scanning Probe Microscopes like the Scanning Near Field Optical Microscope (SNOM) for near field characterization of optical integrated structures.


M Wakim

L Chassagne

S Topsu

Y Alayli

J P Wallerand

P Juncar

G Geneves

Novel concept of interferometric control of the moving coil in watt balance

Keywords: watt balance, kilogram, precision motion control, nanotechnology, heterodyne laser interferometer

Novel concept of interferometric control of the moving coil in watt balance

INTRODUCTION

The kilogram is the last remaining unit of the International System of Units which is still defined by a material artefact. This artefact is made from an alloy of 90 % platinum, 10% iridium. In spite of the good choice of the alloy, Surface contamination and bulk outgassing phenomena can vary the amount of matter in the artefact. At the end of the twentieth century international reference and all national references were compared between them. Variations of 3.10 -8 appeared between the international reference and the mean value of all national prototypes [START_REF] Girard | The Third Periodic Verification of National Prototypes of the Kilogram (1988-1992[END_REF]. Since it is very difficult to know which mass has varied or if all masses have varied, actual reference is maintained till finding a new reference. Several laboratories of metrology develop methods to link the unit of mass to fundamental or atomic constants hoping to find a new definition of the kilogram. One of the most promising methods seems to be the watt balance. The"Laboratoire National de métrologie et d'Essais" (LNE) is developing its own watt balance [START_REF] Geneves | The BNM Watt Balance Project[END_REF]. The principle was proposed by B. P. Kibble (NPL) in 1975 [START_REF] Kibble | Atomic Masses and Fundamental Constants[END_REF]. The aim of the project is to compare a mechanical power with an electromagnetic power in order to link mass to the Planck's constant [START_REF] Fujii | [END_REF][5][6]. The measurement system is mainly composed of a force comparator, a coil immersed in a 1T radial induction field suspended below one of its arms. The coil, weighing 2kg, is made from a 260mm diameter epoxy cylindrical support on which 600 spires are enrolled. The watt balance uses two-mode procedure [START_REF] Eichenberger | Tracing the Planck's constant to the kilogram by electromechanical methods[END_REF]: static mode and dynamic mode. The Laplace's force counterbalances the mass of the etalon in the static mode. It's generated when the coil is supplied by a known current. In the dynamic mode the coil moves vertically in the permanent magnetic field with a velocity of 2mm/s and induces Electromotive force (ε) of 1V. The vertical displacement is insured thanks to a set of six flexure strips. By combining the two-mode results and measuring the electrical quantities using the Josephson Effect [START_REF] Hamilton | Evaluating the uncertainty of Josephson voltage standards[END_REF] and the quantum hall effect [START_REF] Jeckelmann | The quantum Hall effect as an electrical resistance standard[END_REF], mass can be expressed as follows:

h V g f K m J   2 (1)
Where K a constant, f J is the Josephson frequency, g the local acceleration of gravity, V the velocity of the moving-coil, and h the Planck's constant. To obtain "m" with an uncertainty of 10 -8 , all physical parameters of the above formula should be measured and controlled with an uncertainty better than 10 -8 and notably the velocity (V). We propose in this paper a novel method to measure and control velocity dedicated to the French watt balance. The method is based on the use of a heterodyne Michelson's interferometer, a two level translation stage, and a home-made high frequency phase shifting electronic circuit. Our method allows a centimetric displacement with an intrinsic relative accuracy on the velocity of 1.9x10 -9 over 100s. Linear translation stage displaces the coil in the magnetic field with a velocity of 2mm/s. Even with the best high precision translation stages, it is impossible to obtain a stability of the velocity below to 10 -6 in relative value over centimetric displacement. In the other hand Constraints of the experience force the use of a motor with stability on the velocity of 10 -4 (Motor is used in vertical and in a vacuum enclosure). Moreover vertical movement of the coil must be perfectly vertical, horizontal component of the velocity induces an undesired electromotive force and disturb the electrical measurements. For this purpose an intermediate vertical translation stage was developed in the interest to guide the coil according to the vertical direction. The guide stage is decoupled of the linear translation stage and insures straightness within a few micrometers over 80mm of travel. Moving the coil on an 80mm of travel with a relative uncertainty of 10 -9 requests a second stage. Our method is based on the use of a double pass heterodyne Michelson's interferometer, a two-level translation stage, and a homemade high frequency phase-shifting electronic circuit (fig. 1). The two level translation stage is made by superposing a linear motor and a piezoelectric actuator. The piezoelectric actuator, placed between the translation stage and the coil, corrects the defect of the first stage. The coil is replaced by a simple mirror that consist the movable mirror of the interferometer. Commercial laser head generates two orthogonal polarized beams with 20 MHz of difference between frequencies. The first beam constitutes the measurement signal (reflect on MM) and the second one constitutes the reference signal (reflect on SM). Information on the position and the velocity is coded in the output signal of the interferometer. Wavelength of the Laser head is carefully chosen to have a traceable measurement with regard of the meter and the second definition. For this purpose we choose one of the wavelengths recommended by the "comité international des poids et des mesures". The home made high frequency electronic circuit generates three signals s1, s'1, and s2 at the same frequency [10][START_REF] Topçu | High accuracy velocity control method for the French moving coil watt balance[END_REF]. The first one (S1) feeds the Bragg cell and creates a synchronous heterodyning with the electronic signals. Signals (S'1, S2) are phase compared to the signal (S3) at the exit of the interferometer. Error signal resulting of the mixing of S2 and S3 creates the input of the proportional integrator derivative module that controls the piezoelectric actuator length. In static, imagine that vibration on the system moves the mirror away or closer to the reference mirror. On the mixer level appears a non null signal error due to the Doppler Effect. Doppler Effect results in a positive or negative phase shift by an amount of 2π for a λ/4 of displacement, neglecting refractive index of the air. The displacement can be deduced according with the formula:

PRINCIPLE AND EXPERIMETAL SETUP

   2 4     x (2)
Where λ is the wave length of the measurement beam and ΔΦ the phase shift. When the position of the mirror is controlled, the noise level is about 1nm of amplitude. Let us handle the problem at reverse. Suppose that S1 or S2 undergo a phase shift. Here again error  becomes non null and the mirror moves to compensate the phase shift imposed by the electronic circuit. Mirror can move in the 2 ways of the vertical direction depending on where the phase shifts are imposed with a step whose value is defined by the angle of the phase shift. Electronic circuit is able to impose phase shifts up to 2π/32 (5nm) in normal mode and 2π/612 in differential mode (260pm) [START_REF] Topcu | Heterodyne interferometric technique for displacement control at the nanometric scale[END_REF]. Differential mode consists to impose phase shifts in the same time on the two signals (S1 and S2) Knowing that the resolution of our interferometer doesn't exceed 0.31 nm so we work in normal mode to control the velocity with steps of 5nm which is sufficient to get 10 -9 of relative stability on the velocity. When the movable mirror moves at 2mm/s the phase shift becomes too high and shifts the value of the frequency at the exit of the interferometer. Signal demodulated with a reference signal s'1 to reconstitute the Doppler frequency (12,64khz for 2mm/s of velocity). By imposing the same rhythm of phase shift on S1 or S2 (depends on the direction) and by controlling with the servo loop the piezoelectric will change the real time position of the mirror in relation with the motor referential and compensate defects of the motor velocity.

RESULTS AND DISCUSSIONS

The optical signal at the exit of the interferometer is demodulated to obtain the Doppler frequency that contains the information of the velocity and its stability. The resulting signal is sent into a FFT apparatus (SR785, Stanford) and a high stability frequency counter HP53132A with a relative accuracy of 10 -10 . The first apparatus monitors the control quality by examining the peak width. The second one measures the frequency and collects the real time value in a computer. The software uses a non classical statistic method to plot the Allan standard deviation curve [START_REF] Chassagne | Highly accurate positioning control method for piezoelectric actuators based on phase-shifting optoelectronics[END_REF]. This method is a graphical data analysis technique used like a metrological tool to make an estimation of the stability of primary frequency standards [START_REF] Jeckelmann | The quantum Hall effect as an electrical resistance standard[END_REF]. The mirror moves on an 80mm travel at 2 mm/s of velocity corresponding to 40 seconds by travel. We consider that the motor and the translation system take 10 seconds to achieve the transient state (acceleration, deceleration and vibration) at the beginning and the end of the travel leading 20 seconds of acquisition time in steady state. In this last state a measure of the Doppler frequency shift is took by averaging measurements on each second. 20 points are recorded by travel.

Step is repeated 100 times to get 2000 points. Figure 2 represents a standard deviation of the Doppler frequency shift corresponding to a velocity of the moving mirror of 2mm/s in the vertical direction. The results show that a stability of 2x10 -8 is occurred over 1 second. The same figure makes an estimation on stability of the velocity and shows a stability of 1.9x10 -9 over 100 seconds. Results are obtained in an air conditioned room where temperature is controlled to 0.1 degrees Celsius. The watt balance experiment is planned to be placed under vacuum, which will ameliorate results.

ALTERNATIVE METHOD

In spite of good results acquired with the system described above some strength can be harmful to the good functioning of the watt balance. Since the horizontal component of the velocity annoys the measurement of the Josephson Effect it's interesting to measure the pitch and yaw of the coil. Controlling velocity and measuring simultaneously pitch and yaw demands a system able to work on three channels. Thus the laser beam will be divided into three beams. The small power of the laser head used in this method (500microwatt) will be a limitation. The laser source contains two orthogonal linear polarized beams with different frequencies. Practically polarisations present some elipticity and 0.5 degrees of uncertainty on the perpendicularity. Knowing that the watt balance will be installed in a vacuum enclosure, polarizations risk to be lost by the changing of the refractive index of air and by the crossing of the optical window.

Using fiber optic to carry the signal in the vacuum enclosure is not the best solution, Polarizations risk to be mixed in the fiber. Moreover laser head can not be installed in the vacuum enclosure because of enormous thermal dissipation.

Fig. 3: OL, local oscillator at 640MHz; ÷8, ÷128, very low phase noise synthesiser divider; AOM, acousto-optic modulator; AOM Ampli, electronic amplifier for the signal feeds the acousto-optic modulator; PZT piezoelectric actuator; PID, proportional integrator derivative module; FO, fiber optic cable.

Optical isolator

In order to resolve those problems a new concept for a new heterodyne laser source is setup in order to adapt the method to the watt balance (fig. 3). It's very important to find a system that insures: power, linearity, good resolution, small thermal dissipation, traceability and to be fibered. Nowadays no commercial system is able to insure those criteria. We have developed an alternative method able to answer to all our demands. The principle is based on the use of a Nd:yag doubled laser with an accordable wavelength that generates a monomode longitudinal laser beam on 532nm with 20mw of power. Faraday isolator is used to prevent unwanted feedback to the laser and made the optical oscillator more stable. 1mw of the total power is used to lock on laser frequency to a saturation absorption line of the iodine hyperfine structure transition. Locking the laser frequency on the iodine absorption peak permit stability on the frequency of 10 -12 . To reduce the absorption of the laser beam in the iodine cell and to have the best frequency stability, the temperature of the cold finger is fixed on -15 degree Celsius in order to reduce pressure to 0. In the vacuum enclosure the two beams are divided by a beam-splitter (fig. 4). The first two beams are beating together thanks to a Kösters Interference double prism and are phase compared to electronic signal (S1) at 5 MHz. The phase of the beating signal is locked on the electronic signal (R1) via an analogical servo-loop. The servo-loop acts on the piezoelectric actuator to move mirror and change the length of the optical path and than to fix the beating frequency at 5 MHz. system allowing to overcome all frequency changing due to thermal and all other effects in fibers. The other part of the two beams beats at the exit of a homemade double pass heterodyne interferometer. The interferometer is composed from a polarizing beam splitter (PBS) a quarter waveplate (λ/4) and a cone reflector (CR). The beams follow the path like showed in the figure 4. The first beam constitutes the reference beam and is reflected on the holed mirror (SM). The second one constitutes the measurement beam; it passes through the holed mirror and is reflected on the movable mirror (MM). The two beams are beating together again in a Kösters Interference double prism and are phase compared to the electronic signal (R2) at 5 MHz and another analog servo-loop lock-in the phase consequently the position of the movable mirror. As seen previously phase shifts can be imposed on the signal R1 or R2 to move the movable mirror in two directions but this time the resolution is better thanks to a smaller wavelength and difference between the two frequencies.

Where p  the minimum is the displacement achieved in normal mode and λ the wavelength of the laser beam.

CONCLUSION

We have developed a new method to measure and control the velocity of the moving coil in the watt balance with a relative uncertainty of 1.9x10 -9 over 100 seconds. The developed system achieved satisfactory results, however it's not enough adapted to the watt balance. In this perspective a new method is conceived to overcome drawbacks and eliminates weakness points. The system was conceived to be fibered insuring good resolution, high power to measure simultaneously velocity, pitch and yaw. Moreover the system insures linearity, small thermal dissipation, and traceability. The new system will be placed under vacuum; we expect better results and more stability on the velocity. The high power generated by our system will allow in the future the control of the velocity of three lateral points of the coil and thus enables us to surpass the pitch and yaw effects.

Fig: 1 .

 1 Fig: 1. experimental setup. DDS, digital direct signal generator; F(p), proportional integral derivative system; HF, high frequency electronic circuit, MM, movable mirror, SM stationary mirror. FFT, fast Fourier transform apparatus; HPTS, high precision translation stage; DPHMI, double pass heterodyne Michelson's interferometer; PZT piezoelectric actuator;

Fig. 2 :

 2 Fig.2: Allan standard deviation of the Doppler frequency shift for a velocity of 2mm/s in the vertical direction. Each point is averaged over 1 second.

8

 8 Pa in the cell thanks to Peltier Effect Thermoelectric Cooler. The other part of the laser beam passes trough a two acousto-optic modulators. The two first order beams are injected into two optical fibers. Frequencies of the two first order beams at the exit of the acoustooptical modulators are shifted of the value of the frequency that supplies it. The two acousto-optical modulators are fed respectively with two electrical signals S1 and S2 with 80MHz and 85 MHz of frequency. The two first order beams generated from the two acousto-optical modulators with a difference of frequency of 5 MHz are injected in two optical fibers to bring signals into the vacuum enclosure. The first signal is directly injected. The second one passes through a polarizing beam splitter (PBS) and a quarter waveplate before being reflected by the mirror (M). After the passing in the quarter waveplate (λ/4) the polarization is converted from horizontally linear to right circular polarization. Reflected beam with left circular polarization pass again in the same quarter waveplate (λ/4) to obtain a vertically linear polarization and reflects on the polarizing beam splitter (PBS) and than it will be injected in a fiber optic. The mirror is fixed on a piezoelectric actuator (PZT). The system is conceived to lock on the phase at the input of the vacuum enclosure. A home made high frequency electronic card generates two signals at 80 MHz and 85 MHz by dividing the signal of the local oscillator OL (640MHz) by 8 and 128 thanks to a very low phase noise synthesiser divider. The signal at 80 MHz (S1) supplies the first acousto-optical modulator. The other one (R1) is mixed with the signal (S1) at 80 MHz. A band-pass filter isolates all frequencies and passes the signal (S2) at 85 MHz that supplies the second acousto-optical modulator. Another very low phase noise synthesiser divider divides the signal of the local oscillator to obtain the reference signal R2.

Fig. 4 :

 4 Fig.4: BS: beam splitter, KP: Kösters Interference double prism; PBS, polarizing beam splitter; CR, cone reflector; λ/4, quarter waveplate; SM, stationary mirror; MM movable mirror.