Julien Moreau

Vincent Loriette
email: loriette@optique.espci.fr.

Albert-Claude Boccara

Full-field birefringence imaging by thermal-light polarization-sensitive optical coherence tomography. I. Theory

A method for measuring birefringence by use of thermal-light polarization-sensitive optical coherence tomography is presented. The use of thermal light brings to polarization-sensitive optical coherence tomography a resolution in the micrometer range in three dimensions. The instrument is based on a Linnik interference microscope and makes use of achromatic quarter-wave plates. A mathematical representation of the instrument is presented here, and the detection scheme is described, together with a discussion of the validity domain of the equations used to evaluate the birefringence in the presence of white-light illumination.

Introduction

The presence of point defects inside multilayer optical coatings is a major source of losses in optical systems that require low levels of scattering. The performance of laser gyros and interferometric detectors of gravitational waves is directly affected by the amount of scattered light that is present; usually scattering levels larger than a few parts in 10 6 are not tolerated. [START_REF] Mackowski | Different approaches to improve the wavefront of low-loss mirrors used in the VIRGO gravitational wave antenna[END_REF][START_REF] Vinet | Scattered light noise in gravitational wave interferometric detectors: a statistical approach[END_REF] Lowering the number of scattering defects inside optical coatings is a technical challenge that requires the development of diagnostic tools in parallel with optimization of the manufacturing processes. The exact localization of point defects inside a multilayer structure, on the substrate-coating interface, at the interface between layers, or inside layers could produce useful information on the source of contamination. However, a single point defect may scatter only a few parts in 10 6 of the flux incident upon the optical component, so a simple tomographic image of the component's reflectance would not permit the detection of such a small signal, which would be completely masked by the specular reflectances of the various interfaces of the coating. Recently it was proposed to get rid of the specular reflectance signals by performing optical coherence tomography ͑OCT͒ measurements on tilted optical components. [START_REF] Vabre | Imagery of local defects in multilayer components by short coherence length interferometry[END_REF] Another solution could be to detect the change of polarization of the incident light caused by the presence of nonspherical scattering defects. Birefringence imaging has already proved to be a highly efficient way of revealing structures inside various media such as biological and geological specimens. Birefringence signals may be affected by chemical, mechanical, or structural causes, which are usually barely detectable with polarization-insensitive instruments if the causes do not modify the reflectance of the structure. A combination of polarization measurement and OCT has been successfully used to image birefringent structures in biological tissues. 4 -12 In these experiments the axial resolution is determined by the coherence length of the source. Usually superluminescent diodes are used, with a coherence length in the 10 -20-m range. The search for micrometer-scale resolution in classic ͑non-polarization-sensitive͒ OCT has been achieved either with femtosecond lasers [START_REF] Drexler | In vivo ultrahighresolution optical coherence tomography[END_REF] or with the addition of more original superluminescent light sources with bandwidths ten times larger than those of classic superluminescent diodes. [START_REF] Kowalevicz | Ultrahigh resolution optical coherence tomography using a superluminescent light source[END_REF] Recently, similar resolution was obtained with thermal light. [START_REF] Vabre | Thermal-light fullfield optical coherence tomography[END_REF] Thermal light may not be so effective as spatially coherent light sources in producing high-sensitivity systems, but it has the advantages of low cost and ease of implementation on a standard microscope, and, above all, permits the use of detector arrays to perform full-field microscopy. 16 -18 To achieve a phasesensitive ͑PS͒ OCT instrument with micrometer resolution in three dimensions we chose to use a thermal light source in a polarization-sensitive Linnik microscope. The detector is a CCD array, and the detection is based on sine phase-modulation interferometry. [START_REF] Dubois | Sinusoidally phasemodulated interference microscope for high-speed highresolution topographic imagery[END_REF][START_REF] Dubois | Phase-map measurements by interferometry with sinusoidal phase modulation and four integrating buckets[END_REF] In this paper we present a mathematical description of our instrument. We show that, by making two measurements with different orientations of the polarizing elements, one can characterize the magnitude and direction of birefringence or of dichroism and simultaneously obtain phase OCT images, with limited crossover among the various parameters. Crossover is induced by the use of a polychromatic light source. We show how the spectrum's shape affects the measurements of topography and birefringence and under what hypotheses the equations that are valid for monochromatic illumination can be corrected to be used with polychromatic illumination.

Setup

The Linnik microscope ͑a Michelson interferometer with microscope objectives in both arms͒ is the instrument of choice to fulfill the requirements of threedimensional imaging capabilities, lateral and axial resolution in the micrometer range, sensitivity to polarization, [START_REF] Vabre | DNA microarray inspection by interference microscopy[END_REF] and full-field imaging capabilities. The hardest requirement to fulfill is to achieve axial resolution, because it requires working either with highnumerical-aperture objectives [START_REF] Dubois | Phase measurements with wide-aperture interferometers[END_REF] or with a broadband source. The second solution is preferred because only standard sets of microscope objectives with relatively low numerical apertures and classic illumination sources can be used. However, in using a broadband source one has to overcome two difficulties: First, perfectly symmetrical arms are necessary to ensure that the resolution is not spoiled by dispersion mismatch. Then, as is the case in most polarization-sensitive interferometers, quarter-wave plates are present inside the arms. The best achromatic quarter-wave plates are Fresnel rhombs, but their use is difficult in this setup because they shift the beams laterally. The only other solution is to use thin achromatic quarter-wave plates, but those plates are achromatic only over a few hundred nanometers. Thermal light sources, such as xenon lamps, have spectral widths larger than 800 nm, whereas typical thin achromatic quarter-wave plates working in the visible spectrum have a spectral width of 300 nm, which limits the width of the usable light spectrum and thus reduces the achievable resolution. Our system is shown schematically in Fig. 1.

Measurement Method

In this section we show how we calculate the intensity of the field on the detector and extract information from the sample. We start by recalling a few basic mathematical tools that we shall use throughout our calculations. We perform the calculations first for monochromatic and then for broadband illumination. In each case the calculation is divided into two steps: First we calculate a Jones matrix that represents the whole interferometer and show that by making two measurements we can obtain the necessary information with which to evaluate both topographic and birefringence signals. Then we present a measurement method, based on a phasemodulation technique, that permits the topographic and birefringence signals to be recovered from a combination of intensity measurements. We always calculate the field intensity on a pixel that images a finite transverse section of the sample. As the measurement is local ͑i.e., pixels are independent͒, we do not explicitly write the dependence on transverse coordinates in the equations. In the following calculations all phase variables are denoted by the Greek letters for the phase terms induced by the sample and for propagation phase terms; all length variables are denoted by either ␦ for lengths that characterize the sample or ⌬ for lengths of the arms.

A. Basic Mathematical Tools

The basic tools that we use to calculate the state of polarization of the electromagnetic field in an interferometer are Jones matrices. We do not present the theory behind Jones matrix calculus or calculations of elementary Jones matrices, as they can be found in many textbooks. [START_REF] Collett | Polarized Light: Fundamental and Applications[END_REF] Jones matrices allow the amplitude of the transverse optical field emerging from a polarizing component to be calculated if the amplitude of the incident transverse field is known. The fields are represented as vectors on an orthonormal basis ͑e S , e P ͒. The directions of the transverse vectors are chosen depending on the setup geometry. In the following calculations we refer to S and P linear polarization states. In the instrument those directions are fixed by the orientation of a beam splitter cube. If the beam splitter were a perfectly non- polarizing component, any other directions could be used; but in practice a so-called nonpolarizing beam splitter always exhibits some polarization properties, and those imperfections can be used to align any other optical axis in the instrument. For example, if a linearly polarized field is incident upon the beam splitter, the reflected and transmitted fields will be slightly elliptically polarized unless the incident field is perfectly parallel or perfectly perpendicular to the plane of incidence on the beam splitter's interface. It is thus natural to define the S polarization state as an electric field vector perpendicular to the plane of incidence on the beam splitter's interface. The basic Jones matrices that we need are the following. M S and M P are the Jones matrices of vertical and horizontal polarizers:

M S ϭ ͫ 1 0 0 0 ͬ , (1)
M P ϭ ͫ 0 0 0 1 ͬ ; ( 2
)
M Q is the Jones matrix of a quarter-wave plate with its fast axis vertical:

M Q ϭ exp ͩ i 4 ͪͫ 1 0 0 Ϫi ͬ ; (3)
M is the rotation matrix of angle in the ͑e S , e P ͒ plane:

M ϭ ͫ cos sin Ϫsin cos ͬ ; (4)
M Q is the Jones matrix of a quarter-wave plate with its fast axis making an angle with the vertical axis:

M Q ϭ M Ϫ M Q M ϭ exp ͩ i 4 ͪͫ cos 2 Ϫ i sin 2 cos sin ͑1 ϩ i͒ cos sin ͑1 ϩ i͒ sin 2 Ϫ i cos 2 ͬ ; ( 5
)
M Q45 is the Jones matrix of a quarter-wave plate with its fast axis at 45°from the vertical axis:

M Q45 ϭ ͱ2 2 ͫ 1 i i 1 ͬ ; ( 6
)
M R is the Jones matrix of a nonbirefringent and nondichroic surface with a reflection factor r:

M R ϭ r ͫ 1 0 0 1 ͬ . (7)
We have chosen to use positive values for the components of M R ; this means that the reflected field is described on the same ͑e S , e P ͒ basis as the incident field, although the direction of the wave vector is reversed on reflection. This choice has no influence on the results, but it simplifies the calculation because, in this case, the Jones matrix that represents a polarizing element with axes that make an angle with the basis vectors retains its form, independently of the direction of the field incident upon it. As the goal of the measurement is the detection of changes in electromagnetic field properties of a completely polarized beam on reflection on or inside a sample, we should treat the most general case of a birefringent, dichroic, and absorptive sample, and this would lead to the use of a general Jones matrix. [START_REF] Mosin ˜o | Propagation of partially polarized light in dichroic and birefringent media[END_REF] However, we are interested chiefly in the detection of pure birefringence. We thus simplify our calculations by fully treating only this case. M B is the Jones matrix of a birefringent sample with reflection factor r on the basis of its proper axes:

M B ϭ r ͫ exp͓Ϫi͑ B ͞2͔͒ 0 0 exp͓ϩi͑ B ͞2͔͒ ͬ . ( 8
)
Angle B is the phase shift between two orthogonal states of polarization of a field reflected by a structure inside the sample. Its value depends not only on the phase shift induced by the reflection on the structure but also on the birefringence integrated along the round trip of the beam inside the sample. Thus B can be seen not as a local value of birefringence inside a sample but only as an integrated signal. Of course this remark also holds for dichroism and topography of buried structures. M B␤ is the Jones matrix of a birefringent sample on the ͑S, P͒ basis; the proper axes make an angle ␤ with the ͑S, P͒ axes:

The Jones matrix of a dichroic sample on the basis of its proper axes is

M D ϭ ͫ r x 0 0 r y ͬ , (10)
and on the ͑S, P͒ basis

M D␤ ϭ ͫ r x cos 2 ␤ ϩ r y sin 2 ␤ ͑r x Ϫ r y ͒sin ␤ cos ␤ ͑r x Ϫ r y ͒sin ␤ cos ␤ r x sin 2 ␤ ϩ r y cos 2 ␤ ͬ . (11)
In Eq. ͑11͒ the reflection factors may be complex; in this case the sample is also birefringent but the axes of birefringence are the same as the axes of dichroism. We use this matrix to represent our sample and, de-

M B␤ ϭ r ͫ cos 2 ␤ exp͓Ϫi͑ B ͞2͔͒ ϩ sin 2 ␤ exp͓ϩi͑ B ͞2͔͒ Ϫi sin 2␤ sin͑ B ͞2͒ Ϫi sin 2␤ sin͑ B ͞2͒ cos 2 ␤ exp͓ϩi͑ B ͞2͔͒ ϩ sin 2 ␤ exp͓Ϫi͑ B ͞2͔͒ ͬ . (9)
pending on the relative phase and magnitude of r x and r y , we shall present either a dichroic or birefringent sample.

M bsR and M bsT are the reflection and transmission matrices, respectively, of a perfect nonpolarizing beam splitter cube:

M bsR ϭ r bsͫ 1 0 0 1 ͬ , (12)
M bsT ϭ t bsͫ 1 0 0 1 ͬ . (13)
The Jones vector of a P-polarized field of intensity I is

E P ϭ ͱI ͩ 0 1 ͪ . (14)
B. Monochromatic Source

Jones Representation of the Arms

To calculate the field amplitude on the detector we build a single Jones matrix that represents the whole interferometer, using the elementary Jones matrices presented in Subsection 3.A. This matrix is the sum of two matrices, M ref and M sam , that represent, respectively, the reference arm and the sample arm. Although we use a broadband source to obtain the required resolution, it is useful to start the calculations with a monochromatic source. In this way we can deal with much simpler formulas and get better insight into the problem. The broadband source is treated next. The Jones matrix of a reference arm of length ⌬ ref is

M ref ϭ M bsT M Q M R M Q M bsR , ( 15
)
M ref ϭ i ͫ cos 2 sin 2 sin 2 Ϫcos 2 ͬ r bs t bs r ref exp͑i ref ͒. (16)
The global phase term ref represents the phase shift of the field after a round trip in the reference arm.

Although is should be sensitive to the topography of the reference mirror as well as to the modulation of its position, as is explicitly detailed below, the reference mirror is assumed to be perfectly flat. This assumption is justified by the fact that in practice the reference mirror is a superpolished flat with a rms roughness less than 0. The global phase term sam is the phase shift of the field after a round trip in the sample arm.

Field on the Detector for the Monochromatic Source

The field incident upon the beam splitter is P polarized. To extract the amplitude of the birefringence ͑or the dichroism͒ and the direction of the axes we must make two measurements with two orientations of the analyzer. The first measurement is performed with a P-oriented analyzer. The amplitude of the field on the detector is given by

E det ͑P͒ ϭ M P ͑M ref ϩ M sam ͒ E P . ( 20
)
Using Eqs. ͑14͒, ͑18͒, and ͑19͒ in Eq. ͑20͒ and using a perfectly symmetrical beam splitter with r bs ϭ t bs ϭ ͌ 2͞2, we obtain

E det ͑P͒ ϭ ͱI 2 ͫ Ϫr ref cos 2 exp͑i ref ͒ Ϫ r x Ϫ r y 2 exp͑Ϫ2i␤͒exp͑i sam ͒ ͬ . ( 21
)
The intensity is thus

I det ͑P͒ ϭ I 4 ͑r ref cos 2͒ 2 ϩ ͉r x Ϫ r y ͉ 2 4 ϩ r ref cos 2͓Re͑r x Ϫ r y ͒cos͑ Ϫ 2␤͒ Ϫ Im͑r x Ϫ r y ͒sin͑ Ϫ 2␤͔͒ . ( 22
)
Phase ϭ sam Ϫ ref is the phase shift between the two beams. It is different for each pixel. We can write it by making the topography appear explicitly as

͑ x, y͒ ϭ sam Ϫ ref ϩ z ͑ x, y͒. ( 23
)
M ech ϭ M bsR M Q45 M S␤ M Q45 M bsT , (17)
M sam ϭ 1 2 ͫ ͑r x Ϫ r y ͒exp͑2i␤͒ i͑r x ϩ r y ͒ i͑r x ϩ r y ͒ Ϫ͑r x Ϫ r y ͒exp͑Ϫ2i␤͒ ͬ r bs t bs exp͑i sam ͒.

(18)

M sam ϭ i ͫ Ϫsin͑ B ͞2͒exp͑2i␤͒ cos͑ B ͞2͒ cos͑ B ͞2͒ sin͑ B ͞2͒exp͑Ϫ2i␤͒ ͬ r bs t bs r exp͑i sam ͒. (19)
The ͑x, y͒ coordinates are lateral coordinates on the sample. z ͑x, y͒ is the classic topography phase signal, z ͑x, y͒ ϭ 4␦ z ͑x, y͒͞c, where ␦ z ͑x, y͒ is the sample topography. sam Ϫ ref ϭ 4͑⌬ sam Ϫ ⌬ ref ͒͞c is the phase shift averaged over the ͑x, y͒ image field.

It is important to insist on the fact that, if an interference term is generated by birefringence or dichroism, then its phase depends not only on the direction of the principal axis ␤ but also on the sample's topography ␦ z . Unless a pure topography map of the sample is available, it is not possible to extract the value of the direction of the birefringence-dichroism axes independently of the value of ␦ z with only one measurement.

The second measurement is performed with an S-oriented analyzer. In this case we obtain

E det ͑S͒ ϭ ͱI 2 ͫ r ref sin 2 exp͑i ref ͒ ϩ i r x ϩ r y 2
exp͑i sam ͒ ͬ [START_REF] Mosin ˜o | Propagation of partially polarized light in dichroic and birefringent media[END_REF] and, for the intensity,

I det ͑S͒ ϭ I 4 ͭ ͑r ref sin 2͒ 2 ϩ ͉r x ϩ r y ͉ 2 4 Ϫ r ref sin 2͓Re͑r x ϩ r y ͒sin͑͒ ϩ Im͑r x ϩ r y ͒cos͔͑͒ ͮ . (25)
For a dichroic sample we write r x ϩ r y ϭ 2r sam and r x Ϫ r y ϭ ␦r; Eqs. ͑22͒ and ͑25͒ can be rewritten as

I det ͑P͒ ϭ I 4 ͫ ͑r ref cos 2͒ 2 ϩ ␦r 2 4 ϩ r ref ␦r cos 2 cos͑ Ϫ 2␤͒ ͬ , (26)
I det ͑S͒ ϭ I 4 ͓͑r ref sin 2͒ 2 ϩ r sam 2 Ϫ 2r ref r sam sin 2 sin͔͑͒ . ( 27
)
When the sample is birefringent, r x ϩ r y ϭ 2r cos͑ B ͞2͒ and r x Ϫ r y ϭ Ϫ2ri sin͑ B ͞2͒, Eqs. ͑22͒ and ͑25͒ can be rewritten as

I det ͑P͒ ϭ I 4 ͫ ͑r ref cos 2͒ 2 ϩ ͩ r sam sin B 2 ͪ 2 ϩ 2r ref r sam cos 2 sin B 2 sin͑ Ϫ 2␤͒ ͬ , ( 28
)
I det ͑S͒ ϭ I 4 ͫ ͑r ref sin 2͒ 2 ϩ ͩ r sam cos B 2 ͪ 2 Ϫ 2r ref r sam sin 2 cos B 2 sin͑͒ ͬ . ( 29
)
Using Eq. ͑23͒, we write the results of the two measurements as

I det ͑P͒ ϭ I 0 ͑P͒ ϩ A ͑P͒ cos͑ sam Ϫ ref ͒ ϩ B ͑P͒ sin͑ sam Ϫ ref ͒, (30)
I det ͑S͒ ϭ I 0 ͑S͒ ϩ A ͑S͒ cos͑ sam Ϫ ref ͒ ϩ B ͑S͒ sin͑ sam Ϫ ref ͒, (31)
with

I 0 ͑P͒ ϭ 1 4 I ͫ ͑r ref cos 2͒ 2 ϩ ͩ r sam sin B 2 ͪ 2 ͬ , ( 32
)
A ͑P͒ ϭ ϩ I 2 r ref r sam cos 2 sin B 2 sin͑ z Ϫ 2␤͒, ( 33
)
B ͑P͒ ϭ ϩ I 2 r ref r sam cos 2 sin B 2 cos͑ z Ϫ 2␤͒, ( 34
)
I 0 ͑S͒ ϭ 1 4 I ͫ ͑r ref sin 2͒ 2 ϩ ͩ r sam cos B 2 ͪ 2 ͬ , ( 35
)
A ͑S͒ ϭ Ϫ I 2 r ref r sam sin 2 cos B 2 sin͑ z ͒, (36)
B ͑S͒ ϭ Ϫ I 2 r ref r sam sin 2 cos B 2 cos͑ z ͒. (37)
To measure the amplitude of the birefringence signal we must first extract the four values A ͑P͒ , B ͑P͒ , A ͑S͒ , and B ͑S͒ . Then the birefringence amplitude can be calculated from

tan 2 ͩ B 2 ͪ ϭ tan Ϫ2 ͑2͒ A ͑P͒2 ϩ B ͑P͒2 A ͑S͒2 ϩ B ͑S͒2 ; ( 38
)
the topography from

tan͑ z ͒ ϭ A ͑S͒ B ͑S͒ ; (39)
and the birefringence direction from

tan͑2␤͒ ϭ A ͑S͒ B ͑P͒ Ϫ A ͑P͒ B ͑S͒ B ͑S͒ B ͑P͒ ϩ A ͑S͒ A ͑P͒ . ( 40
)
The value of B is of course independent of and, in most papers that treat polarization-sensitive OCT experiments, ϭ ͞8. This configuration is useful if the birefringence amplitude is unknown and can a priori take any value. We prefer to leave the value of unspecified because in some cases, for example, when one is testing nonbiological samples, the birefringence may be small everywhere. In these cases, choosing a lower value for enhances the weight of the sin͑ B ͞2͒ terms and thus the signal-to-noise ratio.

Modulation

FWHM ϭ 0.6 1 Ϫ cos͑arcsin N.A.͒ , ( 41
)
where N.A. is the numerical aperture of the objective and is the wavelength. In practice the modulation amplitude is always smaller than this length because one great advantage of a broadband source is the possibility of using objectives with a depth of field large compared to the average wavelength, because the axial resolution is limited by the coherence length. By using objectives with numerical apertures smaller than 0.6 we can assume that the contrast does not depend on the position of the reference mirror.

We write sam Ϫ ref as

sam Ϫ ref ϭ 0 sin͑t ϩ ͒, ( 42
)
where 0 is the modulation depth and phase term is the phase shift between the reference mirror's modulated position and the CCD acquisition sequence. These two parameters are in practice tunable by electronic means. The constraint that 0 Ͻ 2 (43) ensures that the modulation amplitude is small compared to the depth of field of the objective:

sin͑ sam Ϫ ref ͒ ϭ 2 ͚ kϭ0 ϱ J 2kϩ1 ͑ 0 ͒sin͓͑2k ϩ 1͒͑t ϩ ͔͒, ( 44
) cos͑ sam Ϫ ref ͒ ϭ J 0 ͑ 0 ͒ ϩ 2 ͚ kϭ1 ϱ J 2k ͑ 0 ͒cos͓2k͑t ϩ ͔͒. ( 45
)
The CCD array used to image the sample integrates the signal, so we cannot directly extract the signal at the modulation frequency. Instead, we use the integration property of the detector to acquire four con-secutive images, labeled S 1 -S 4 , integrated during one quarter of modulation period T. Figure 2 is a simulated signal recorded on one pixel of the CCD array:

S q ͑P͒ ϭ ͐ ͑qϪ1͒T͞4 qT͞4 I det ͑P͒ ͑t͒dt, S q ͑S͒ ϭ ͐ ͑qϪ1͒T͞4 qT͞4 I det ͑S͒ ͑t͒dt. (46)
Substituting Eqs. ͑42͒, ͑44͒, and ͑30͒ into Eq. ͑46͒, we have

S q ͑P͒ ϭ T 4 ͓I 0 ͑P͒ ϩ J 0 ͑ 0 ͒ A ͑P͒ ͔ ϩ TA ͑P͒ ͚ kϭ1 ϱ J 2k ͑ 0 ͒ 2k ( sin ͫ 2k ͩ q 2 ϩ ͪͬ Ϫ sin ͭ 2k ͫ ͑q Ϫ 1͒ 2 ϩ ͬͮ ) Ϫ TB ͑P͒ ͚ kϭ0 ϱ J 2kϩ1 ͑ 0 ͒ 2k ϩ 1 ( cos ͫ ͑2k ϩ 1͒ ͩ q 2 ϩ ͪͬ Ϫ cos ͭ ͑2k ϩ 1͒ ͫ ͑q Ϫ 1͒ 2 ϩ ͬͮ ) , (47)
and, using Eq. ͑31͒, we obtain a similar expression for S q ͑S͒ . To extract A ͑P͒ and B ͑P͒ we use two different linear combinations of the four values of S q ͑P͒ : Substituting Eq. ͑47͒ into Eqs. ͑48͒ and ͑49͒ leads to

⌺ B ͑P͒ ϭ ϪS 1 ͑P͒ ϩ S 2 ͑P͒ ϩ S 3 ͑P͒ Ϫ S 4 ͑P͒ , ( 48
)
⌺ A ͑P͒ ϭ ϪS 1 ͑P͒ ϩ S 2 ͑P͒ Ϫ S 3 ͑P͒ ϩ S 4 ͑P͒ . ( 49
)
⌺ B ͑P͒ ϭ 4T B ͑P͒ ⌫ B ͑ 0 , ͒, (50)
⌺ A ͑P͒ ϭ 4T A ͑P͒ ⌫ A ͑ 0 , ͒, (51)
with

⌫ B ͑ 0 , ͒ ϭ ͚ kϭ0 ϱ ͑Ϫ1͒ k J 2kϩ1 ͑ 0 ͒ 2k ϩ 1 sin͓͑2k ϩ 1͔͒, ( 52
)
⌫ A ͑ 0 , ͒ ϭ ͚ kϭ1 ϱ ͓1 Ϫ ͑Ϫ1͒ k ͔ J 2k ͑ 0 ͒ 2k sin͑2k͒. ( 53
)
We use the same combinations to extract A ͑S͒ and B ͑S͒ . Equation ͑38͒ can thus be rewritten as

Evaluation of Modulation Depth 0 and Modulation Phase

To calculate the value of ⌫ A ͑ 0 , ͒ and ⌫ B ͑ 0 , ͒ we need to measure 0 and . We could rely on a precise calibration of the reference mirror's displacement system to obtain 0 , but we have no easy way to measure . To overcome this difficulty we need to use other combinations of images to fix those parameters to known values. Once those values are known, we can modify them to optimize the signal-to-noise ratio. As they are independent of the orientation of the analyzer, we can choose to measure them with a P-oriented analyzer. We use the following notation for the image background:

S 0 ͑P͒ ϭ T 4 ͓I 0 ͑P͒ ϩ J 0 ͑ 0 ͒ A ͑P͒ ͔. ( 55
)
The four images can be written as

S 1 ͑P͒ ϭ S 0 ͑P͒ ϩ T ͓ϪA ͑P͒ ⌫ A ͑ 0 , ͒ Ϫ B ͑P͒ f 2 ͑ 0 , ͔, ( 56
)
S 2 ͑P͒ ϭ S 0 ͑P͒ ϩ T ͓A ͑P͒ ⌫ A ͑ 0 , ͒ Ϫ B ͑P͒ f 3 ͑ 0 , ͔͒, (57)
S 3 ͑P͒ ϭ S 0 ͑P͒ ϩ T ͓ϪA ͑P͒ ⌫ A ͑ 0 , ͒ ϩ B ͑P͒ f 2 ͑ 0 , ͔͒, (58)
S 4 ͑P͒ ϭ S 0 ͑P͒ ϩ T ͓A ͑P͒ ⌫ A ͑ 0 , ͒ ϩ B ͑P͒ f 3 ͑ 0 , ͔͒, (59)
with

f 2 ͑ 0 , ͒ ϭ ͚ kϭ0 ϱ J 2kϩ1 ͑ 0 ͒ 2k ϩ 1 ͕cos͓͑2k ϩ 1͔͒ ϩ ͑Ϫ1͒ k sin͓͑2k ϩ 1͔͖͒, (60)
f 3 ͑ 0 , ͒ ϭ ͚ kϭ0 ϱ J 2kϩ1 ͑ 0 ͒ 2k ϩ 1 ͕cos͓͑2k ϩ 1͔͒ Ϫ ͑Ϫ1͒ k sin͓͑2k ϩ 1͔͖͒. (61)
It is easy to show that

⌫ B ͑ 0 , ͒ ϭ 1 ⁄2 ͓ f 2 ͑ 0 , ͒ Ϫ f 3 ͑ 0 , ͔͒. (62)
A practical way to find a specific value of is to tune it until the two images S 2 and S 4 become identical. At this point the value of f 3 ͑ 0 , ͒ vanishes. If this holds for any value of 0 , which we check by changing the modulation amplitude, then necessarily ϭ ͞4.

We can then compute

⌺ A ͑P͒2 ϩ ⌺ B ͑P͒2 : ⌺ A ͑P͒2 ϩ ⌺ B ͑P͒2 ϭ ͩ 4T ͪ 2 ͕͓A ͑P͒ ⌫ A ͑ 0 , ͔͒ 2 ϩ ͓B ͑P͒ ⌫ B ͑ 0 , ͔͒ 2 ͖. (63)
Then to fix the value of 0 we use the fact that, when 0 ϭ 2.0759 rad ͓when ⌫ A ͑ 0 , ͒ ϭ ⌫ B ͑ 0 , ͔͒, we obtain

⌺ A ͑P͒2 ϩ ⌺ B ͑P͒2 ϭ ͫ 4T ⌫ A ͑ 0 , ͒ I 2 r ref r sam cos 2 sin B 2 ͬ 2 , (64)
which is obviously independent of the topography signal. In practice one generates a topography signal by tilting the reference mirror such as to get many fringes in the images and then adjusts 0 to eliminate this topography signal by making the fringes vanish. So, although we do not directly measure 0 and , we can calculate B by from

tan 2 ͩ B 2 ͪ ϭ tan Ϫ2 ͑2͒ ⌺ A ͑P͒2 ϩ ⌺ B ͑P͒2 ⌺ A ͑S͒2 ϩ ⌺ B ͑S͒2 . (65)

C. Broadband Source

We can now study the case of a broadband source, with an apparent spectrum f ͑͒. In what follows, all numerical examples are calculated with a sample spectrum that is flat from 650 to 950 nm. The birefringence amplitude depends on wavelength, and we write B ͑͒ ϭ 2␦ B ͞c. The direction of birefringence axes ␤ is a geometrical characteristic of the

tan 2 ͩ B 2 ͪ ϭ tan Ϫ2 ͑2͒ ͓⌫ B ͑ 0 , ͒⌺ A ͑P͒ ͔ 2 ϩ ͓⌫ A ͑ 0 , ͒⌺ B ͑P͒ ͔ 2 ͓⌫ B ͑ 0 , ͒⌺ A ͑S͒ ͔ 2 ϩ ͓⌫ A ͑ 0 , ͒⌺ B ͑S͒ ͔ 2 . ( 54
)
sample, which depends mainly on the local orientation of the structures that compose it. We make the hypothesis that within a given structure the directions of the birefringence axes do not change significantly with wavelength. Our goal is to obtain an equation, similar to Eq. ͑54͒, that will allow us to evaluate ␦ B by using combinations of the eight images. As for the monochromatic source, we first have to find a way to measure the modulation parameters. When the beam splitter, the reference mirror, the quarter-wave plates, and the polarizers are perfectly achromatic, we can write Eqs. ͑28͒ and ͑29͒ with an explicit dependence on as

I det ͑P͒ ͑͒ ϭ I 4 ͭ ͑r ref cos 2͒ 2 ϩ ͫ r sam sin B ͑͒ 2 ͬ 2 ϩ 2r ref r sam cos 2 sin B ͑͒ 2 sin͓͑͒ Ϫ 2␤͔ ͮ , (66)
I det ͑S͒ ͑͒ ϭ I 4 ͭ ͑r ref sin 2͒ 2 ϩ ͫ r sam cos B ͑͒ 2 ͬ 2 Ϫ 2r ref r sam sin 2 cos B ͑͒ 2 sin͓͔͑͒ ͮ . (67)
Modulation phase does not depend on wavelength but the modulation depth does, as 0 ϭ 4␦ 0 ͞c, where ␦ 0 is the oscillation amplitude of the reference mirror. The four images are now

S 1 ͑P͒ ϭ S 0 P ϩ T ͫ Ϫ ͐ 0 ϱ f ͑͒ A ͑P͒ ͑, ␦ B ͒⌫ A ͑, ␦ 0 , ͒d Ϫ ͐ 0 ϱ f ͑͒ B ͑P͒ ͑, ␦ B ͒ f 2 ͑, ␦ 0 , ͒d ͬ , (68) S 2 ͑P͒ ϭ S 0 P ϩ T ͫ ͐ 0 ϱ f ͑͒ A ͑P͒ ͑, ␦ B ͒⌫ A ͑, ␦ 0 , ͒d Ϫ ͐ 0 ϱ f ͑͒ B ͑P͒ ͑, ␦ B ͒ f 3 ͑, ␦ 0 , ͒d ͬ , (69) S 3 ͑P͒ ϭ S 0 P ϩ T ͫ Ϫ ͐ 0 ϱ f ͑͒ A ͑P͒ ͑, ␦ B ͒⌫ A ͑, ␦ 0 , ͒d ϩ ͐ 0 ϱ f ͑͒ B ͑P͒ ͑, ␦ B ͒ f 2 ͑, ␦ 0 , ͒d ͬ , (70) S 4 ͑P͒ ϭ S 0 P ϩ T ͫ ͐ 0 ϱ f ͑͒ A ͑P͒ ͑, ␦ B ͒⌫ A ͑, ␦ 0 , ͒d ϩ ͐ 0 ϱ f ͑͒ B ͑P͒ ͑, ␦ B ͒ f 3 ͑, ␦ 0 , ͒d ͬ . (71)
We cannot at this point use the equations that are valid for the monochromatic source. The result would not be a measurement of the birefringence weighted by the source spectrum, even if there were no topography signal ͑␦ z ϭ 0͒, because the functions ⌫ A and ⌫ B vary significantly with wavelength, as can be seen from Fig. 3, and thus cannot be extracted from the integrals of Eqs. ͑68͒-͑71͒.

D. Evaluation of Oscillation Amplitude ␦ 0 and Modulation Phase for the Broadband Source

To fix the values of ␦ 0 and we perform the same operations as for the monochromatic source. For the problem is simply solved because f 3 ͑, ␦ 0 , ͞4͒ ϭ 0@, ␦ 0 . So, by making S 2 ͑P͒ ϭ S 4

͑P͒

, we obtain ϭ ͞4. It is more difficult to handle ␦ 0 , and in fact there is no way to cancel S 2 ͑P͒ Ϫ S 1 ͑P͒ for any wavelength or for ⌺ A ͑S͒2 ϩ ⌺ B ͑S͒2 to be independent of the topography for any wavelength. We establish the setting by using a nonbirefringent sample and working with an S-oriented analyzer: Equation ͑72͒ cannot be made independent of ␦ z whatever its value, but we can try to make our notation insensitive to ␦ z near a particular value, for example, ␦ z ϭ 0. We want to find ␦ 0opt that satisfies the equation (75) or, substituting Eqs. ͑73͒ and ͑74͒ into Eq. ͑72͒,

⌺ A ͑S͒2 ϩ ⌺ B ͑S͒2 ϭ ͩ 4T ͪ 2 ( ͭ ͐ 0 ϱ f ͑͒ A ͑S͒ ϫ ͑, ␦ B ͒⌫ A ͓͑, ␦ 0 , ͔͒d ͮ 2 ϩ ͭ ͐ 0 ϱ f ͑͒ B ͑S͒ ϫ ͑, ␦ B ͒⌫ B ͓͑, ␦ 0 , ͔͒d ͮ 2 ) , (72)
d d␦ z ͫ ⌺ A ͑S͒2 ϩ ⌺ B ͑S͒2 ͬͯ␦ z Ϸ 0 ϭ 0,
ͫ ͐ 0 ϱ f ͑͒⌫ A ͑, ␦ 0opt ͒d ͬ 2 Ϫ ͫ ͐ 0 ϱ f ͑͒⌫ B ͑, ␦ 0opt ͒d ͐ 0 ϱ f ͑͒⌫ B ͑, ␦ 0opt ͒ 2 d ͬ ϭ 0. (76)
We can verify that, when f ͑͒ is the Dirac distribution ␦͑ Ϫ 0 ͒, we end with ⌫ A ͑S͒ ϭ ⌫ B ͑S͒ . For an arbitrary spectrum the equation must be solved numerically. The left-hand side of Eq. ͑76͒ is plotted in Fig. 4. It was calculated from the flat spectrum from 650 to 950 nm. The first zero of this function gives the value of ␦ 0 , which makes Eq. ͑72͒ independent of ␦ z near ␦ z ϭ 0. One can perform this calculation experimentally as in the monochromatic case by tilting the reference mirror to obtain Fizeau fringes and then finding the value of ␦ 0 that make the fringes at the center of the interferogram disappear. The width of this region, ⌬ z , is nearly equal to the coherence length of the source. Figure 5 shows a calculation of such an interferogram, in which ␦ z varies from Ϫ2 to ϩ2 m and ␦ 0 varies from 0.85 to 1.15 times its optimal value. We can clearly see the central fringes disappear as ␦ 0 reaches its optimal value. From the value of ␦ 0opt we can calculate an equivalent wavelength 0 , defined by

⌫ Aͩ 4 ␦ 0opt 0 , 4 ͪ ϭ ⌫ Bͩ 4 ␦ 0opt 0 , 4 ͪ . ( 77
)
The equivalent wavelength can be thought of as the wavelength that should be used in the monochromatic case with the same modulation amplitude ␦ 0opt .

Once the values of ␦ 0 and are fixed we can try to generalize Eq. ͑54͒ to encompass a broadband source. We define ⌰ A and ⌰ B by

⌰ A ͑␦ 0opt ͒ ϭ ͐ 0 ϱ f ͑͒⌫ Aͫ ͩ , ␦ 0opt , 4 ͪͬ d, ⌰ B ͑␦ 0opt ͒ ϭ ͐ 0 ϱ f ͑͒⌫ Bͫ ͩ , ␦ 0opt , 4 ͪͬ d (78)
and estimate birefringence ␦ B by using

tan 2 ͩ ␦ B 0 ͪ ϭ tan Ϫ2 ͑2͒ ϫ ͓⌰ B ͑␦ 0opt ͒⌺ A ͑P͒ ͔ 2 ϩ ͓⌰ A ͑␦ 0opt ͒⌺ B ͑P͒ ͔ 2 ͓⌰ B ͑␦ 0opt ͒⌺ A ͑S͒ ͔ 2 ϩ ͓⌰ A ͑␦ 0opt ͒⌺ B ͑S͒ ͔ 2 , ( 79
) topography ␦ z with tan ͩ 2 ␦ z 0 ͪ ϭ ⌰ B ͑␦ 0opt ͒⌺ A ͑S͒ ⌰ A ͑␦ 0opt ͒⌺ B ͑S͒ , (80)
and the birefringence direction as

tan͑2␤͒ ϭ ⌰ A ͑␦ 0opt ͒⌰ B ͑␦ 0opt ͓͒⌺ A ͑S͒ ⌺ B ͑P͒ Ϫ ⌺ A ͑P͒ ⌺ B ͑S͒ ͔ ⌰ A 2 ͑␦ 0opt ͒⌺ B ͑S͒ ⌺ B ͑P͒ ϩ ⌰ B 2 ͑␦ 0opt ͒⌺ A ͑S͒ ⌺ A ͑P͒ . (81)
For a monochromatic source, only the measurement of the direction of birefringence axes ␤ was affected by the topography. For the broadband source, neither ␤ nor ␦ B can be made rigorously independent of ␦ z . Figure 6 shows the calculated value of the topography plotted relative to the true topography. Figure 7 shows the minimum and maximum values of the estimated birefringence as functions of the true birefringence for two values of the birefringence direction, as ␦ z varies from ␦ z ϭ 0 to ␦ z ϭ 0.75 ϫ 0 . The estimation error depends on the value of ␦ z . For values of ␦ z Ͼ ⌬ z ͞2, Eqs. ͑79͒-͑81͒ cannot be used. This simply means that, outside the coherence length of the source, the measurement cannot be made, which is a standard limitation of phase detection methods. [START_REF] Hariharan | White-light phase stepping interferometry: measurement of the fractional interference order[END_REF] A way to overcome this limitation is first to perform a test on the interferometer envelope [START_REF] Harasaki | Improved verticalscanning interferometry[END_REF] to discriminate between points at which the evaluation birefringence can or cannot give accurate results. One could of course use efficient combinations of more than four images [START_REF] Larkin | Efficient nonlinear algorithm for envelope detection in white light interferometry[END_REF][START_REF] Sandoz | An algorithm for profilometry by white-light phase-shifting interferometry[END_REF] to obtain a precise measurement of the coherence envelope, but this test should not be used directly to measure topography, as it would give poor results. Let function ⍀ be defined by

⍀͑␦ B , ␤, ␦ z ͒ ϭ ⌺ A ͑S͒2 ϩ ⌺ A ͑P͒2 ϩ ⌺ B ͑S͒2 ϩ ⌺ B ͑P͒2 . ( 82
)
Figure 8 is a plot of the maximum and minimum values of ⍀ as functions of ␦ z when B takes the values 0 and rad. We see that ⍀ is only slightly dependent on the birefringence. ⍀ can thus be used as an estimate of the coherence envelope. By calculating the value of this function at each point on the image field we roughly estimate the topography, independently of the birefringence. In practice, for recording full-field images, one fixes the sample upon an axial translation stage that permits step-by-step movement with step size d z smaller than the FWHM of ⍀. For a given position of the sample, the birefringence is evaluated only inside the region where the rough topography has a value smaller than d z ͞2. This condition determines the threshold level that is used with the function ⍀.

Discussion

We have presented a method for combining phasesensitive OCT and thermal-light OCT to obtain polarization-sensitive tomographic images. We have shown that, within determined limits, birefringence or dichroism and a tomographic phase can be measured nearly independently of each other. A combination of coherence envelope detection and tomographic phase measurement permits discrimination between points where the evaluation of birefringence is strongly dependent on topography and points where simple formulas can be used to obtain accurate results. We have presented formulas that explicitly include the spectrum of the light source and allow us to evaluate birefringence by use of a thermal-light polarization-sensitive optical coherence tomography instrument. Micrometer axial resolution is then possible with such a white-light source. These formulas ͓Eqs. ͑79͒-͑81͔͒ are simple combinations of eight images; they hide the experimental difficulties introduced by the use of imperfectly achromatic components and by the need to combine two different arrangements of polarizing elements. The polarization scheme used rejects light whose polarization state has not been modified by the sample. This property, combined with micrometer resolution, makes this instrument a valuable tool for the study of small-scale birefringent defects inside optical components, in particular, inside optical multilayer coatings. Although our primary motivation for developing this measurement method was to use it for the study of optical components, biological applications could of course also benefit from it.

Fig. 1 .

 1 Fig. 1. Schematic of the instrument: NPBS, nonpolarizing beam-splitter cube; AQWPs, achromatic quarter-wave plates; PZT, piezoactuated translation stage.

Fig. 2 .

 2 Fig. 2. Typical interferometric signal I det recorded by a CCD pixel as a function of time. This signal is integrated over four consecutive quarters of the modulation period to give four images, S 1 -S 4 .

withAFig. 3 .

 3 Fig. 3. Values of ⌫ A and ⌫ B as a function of wavelength for ␦ 0 ϭ ␦ 0opt .

Fig. 4 .

 4 Fig. 4. Value of Eq. ͑76͒ as a function of the reference mirror's modulation amplitude: circle, point for calculating Fig. 5͑b͒; triangles, points used in Figs. 5͑a͒ and 5͑c͒.

Fig. 5 .

 5 Fig. 5. Simulation of the procedure proposed for fixing the value of ␦ 0 : interferograms with micrometer times micrometer dimensions: ͑a͒ ␦ 0 ϭ 0.85 ϫ ␦ 0opt , ͑b͒ ␦ 0 ϭ ␦ 0opt , and ͑c͒ ␦ 0 ϭ 1.15 ϫ ␦ 0opt .

Fig. 6 .

 6 Fig. 6. Estimated versus true topography ␦ z for a source spectrum that is flat from 650 to 950 nm.

Fig. 7 .

 7 Fig. 7. Magnitudes of estimated versus true birefringence ␦ B for various values of topography ␦ z ranging from 0 ͑dashed curves͒ to 0.75 ϫ 0 ͑solid curves͒ for the flat source spectrum: ͑a͒ ␤ ϭ ͞4, ͑b͒ ␤ ϭ 0.

Fig. 8 .

 8 Fig. 8. Values of ⍀ as a function of ␦ z expressed in units of 0 for B ϭ 0 ͑solid curve͒ and B ϭ ͑dashed curve͒.

  1 nm. Therefore ref depends only on the reference mirror's position, ref ϭ 4⌬ ref͞c . The Jones matrix of the sample arm is For a birefringent sample r x ϩ r y ϭ 2r cos͑ B ͞2͒ and r x Ϫ r y ϭ Ϫ2ri sin͑ B ͞2͒. The Jones matrix for the sample arm of length ⌬ sam reads as

  ͑30͒ and ͑31͒ by use of a different set of values for sam Ϫ ref and by calculation of different linear combinations of those results. Four values are necessary to produce a set of two parameters, A ͑P͒ and B ͑P͒ or A ͑S͒ and B ͑S͒ . One might think that only two values would be necessary, but, as we make clear below, two extra parameters appear, a phase shift and a modulation depth, that need to be evaluated too. We can change sam Ϫ ref , for example, by moving the reference mirror to any of four positions.

	We
	chose to modulate the reference arm's length sinusoi-
	dally with angular frequency . This signal extrac-
	tion method can be used only if the reference mirror's
	position span is shorter than the depth of field of the
	microscope objectives; otherwise different positions of
	the reference mirror will lead not only to a different
	state of interference but also to a different contrast of
	the interference signal. An expression of the
	FWHM of the contrast function can be found in Ref.
	22:

In this last step of the calculation we describe a practical way to compute the values of A ͑P͒ , B ͑P͒ , A ͑S͒ , and B ͑S͒ . These values can be obtained from Eqs.

July

͞ Vol. 42, No. 19 ͞ APPLIED OPTICS

APPLIED OPTICS ͞ Vol. 42, No. 19 ͞ 1 July 2003

The authors thank Arnaud Dubois and Laurent Vabre for helpful discussions.


Julien Moreau

Vincent Loriette
email: loriette@optique.espci.fr.

Albert-Claude Boccara

Full-field birefringence imaging by thermal-light polarization-sensitive optical coherence tomography. II. Instrument and results

We describe an instrument for measuring the magnitude of birefringence of tomographic images and the principal directions of axes that uses thermal-light polarization-sensitive optical coherence tomography. The instrument permits full-field measurements with an axial resolution of 1.5 m and a transverse resolution limited by diffraction. We obtained a sensitivity of 84 dB, limited by shot noise, when we integrated the signal for 1 s. We verified the validity of the measurement by measuring the birefringence of a variable phase shifter. We present typical results obtained with optical samples.

Introduction

Optical coherence tomography ͑OCT͒ is an attractive technique for obtaining information about buried structures with high spatial resolution. OCT is oriented mainly toward the study of diffusive samples, but it can be used as well for testing transparent samples with low optical contrast. In this case polarization-sensitive OCT, which brings the advantages of dark-field microscopy, proves to be an effective technique for revealing structures that otherwise would hardly be detectable. [START_REF] De Boer | Two-dimensional birefringence imaging in biological tissue by polarization-sensitive optical coherence tomography[END_REF][START_REF] De Boer | Imaging thermally damaged tissue by polarization sensitive optical coherence tomography[END_REF][START_REF] Everett | Birefringence characterization of biological tissue by use of optical coherence tomography[END_REF][START_REF] Schoenenberger | Mapping of birefringence and thermal damage in tissue by use of polarization-sensitive optical coherence tomography[END_REF][START_REF] De Boer | Determination of the depth-resolved Stokes parameters of light backscattered from turbid media by use of polarization-sensitive optical coherence tomography[END_REF][START_REF] Yao | Two-dimensional depth-resolved Mueller matrix characterization of biological tissue by optical coherence tomography[END_REF][START_REF] Saxer | High-speed fiber-based polarization-sensitive optical coherence tomography of in vivo human skin[END_REF][START_REF] Hitzenberger | Measurement and imaging of birefringence and optic axis orientation by phase resolved polarization sensitive optical coherence tomography[END_REF][START_REF] Roth | Simplified method for polarization-sensitive optical coherence tomography[END_REF] The combination of a polarization-sensitive instrument with a broadband source allows one to obtain polarizationdependent information with micrometer-scale resolution in three dimensions. With a CCD array, full-field microscopy is possible. 10 -12 The coherence length of the broadband source determines the axial resolution of the instrument, so microscope objective lenses with large depths of field or low numerical apertures can be used without spoiling the axial resolution. Unless one uses immersion objectives and samples with low index variations, the use of low-numerical-aperture objectives is obligatory when one tests samples that exhibit large index changes because aberrations inside the sample put severe limitations on the depth that can be probed. However, the use of a broadband source raises a major technical problem: Although achromatic polarizers can be manufactured easily, making a fairly achromatic phase-shifting component is a difficult task. The most common achromatic phase-shifting component is the Fresnel rhomb. Standard commercially available Fresnel rhombs exhibit phase errors of a few degrees over a span of a few hundred nanometers. Unfortunately the standard design of Fresnel rhombs makes them unsuitable for OCT microscopy because they shift the emerging beam laterally and are highly sensitive to misalignment. The use of multilayer phaseshifting plates permits the problem of beam shift to be overcome, but phase-shifting plates are far less efficient than rhombs in the blue part of the spectrum, and their performance varies greatly from one manufacturer to the other. The present authors have chosen to use achromatic quarter-wave plates, chiefly because of their ease of use, although, for some experimental configurations, this choice is made to the detriment of accuracy and sensitivity. In Section 2 of this paper we describe the experimental arrangement, and in Section 3 we discuss its performance. We focus on elimination of the spurious signal induced by phase errors in the quarter-wave plates. In Section 4 we present some results for multilayer optical samples.

Experimental Arrangement

The experimental arrangement is sketched in Fig. 1. The instrument is based on a Linnik interference microscope, i.e., a Michelson interferometer with a microscope objective lens in each arm. The reference surface should have topographic features much smaller than the sample, and to optimize the contrast its reflectance should be as close as possible to the reflectance of the sample. Uncoated glass, which has a low reflection factor of ϳ4% and can be superpolished, is a good choice when one is testing biological or transparent optical samples. Light is emitted by a xenon lamp with 300-nm bandwidth centered about 800 nm. It is polarized with a Glan-Taylor polarizer that has an extinction ratio of 2 ϫ 10 Ϫ5 . The beam is divided by a nonpolarizing beam-splitter cube. Each arm contains an achromatic quarter-wave plate ͑Newport Corporation Model ACWP͒ and a 20ϫ, microscope objective with a N.A. of 0.4 ͑Olympus 20ϫ plan achromat͒. The emerging light is collected by a CCD array ͑256 ϫ 256 pixel, 8-bit, maximum 200-Hz readout frequency; Dalsa Semiconductors Model CAD-1͒ after it has passed through a thin plate polarizer with an extinction ratio of 4 ϫ 10 Ϫ5 . The plate polarizer is fixed upon a motorized rotation stage to permit angular positioning with a precision of 0.17 mrad. Both reference mirror and sample are fixed upon motorized translation stages. The reference is attached to a piezoelectric actuator oscillating at 45 Hz, and the sample is translated with a 74-nm step motorized translation stage. The amplitude of piezo-oscillation can be adjusted to produce a particular modulation depth. The piezoactuator driver and the CCD are phase locked by two waveform generators ͑Hewlett-Packard Model 33120A͒. The total acquisition time depends on the number of averages for each position of the sample and on the total volume scanned; a good order of magnitude is 4 s͞section for a 100-image average.

Setup and Measurement Procedure

The orientations of the various polarization components are adjusted by use of extinction configurations. At the end of this step the polarizer is oriented to produce a P-polarized field incident upon the beam splitter. The quarter-wave plate in the sample arm is oriented such that its fast axes make an angle of 45°with the P axis, and the quarter-wave plate in the reference arm is oriented at 11°. These values allows the signal-to-noise ratio to be maximized, as explained below. During this step of the setup the CCD is replaced by a silicon photodetector for better sensitivity. With digital waveform generators, both modulation depth and phase are adjusted independently by use of the method discussed in Ref. 13. By measuring four images with a tilted nonbirefringent test sample and calculating S 2 ͑P͒ -S 4 ͑P͒ we can adjust the modulation phase to ͞4. By making the fringes disappear at the center of the interferogram, as shown in Fig. 2, we can fix the modulation depth close to its theoretical value of 2.0759 rad with an estimated precision of 98%. Once the modulation parameters are fixed, the nonbirefringent test sample is replaced by the sample that we wish to study. We then record two sets of four images ͕S i ͑S͒ , S i ͑P͒ ͖, i ϭ 1-4, acquired with two different orientations, S and P, of the analyzer. For a given analyzer position the corresponding four images are acquired consecutively; each image is integrated during one quarter of the modulation period. This eightimage acquisition is performed for each position of the sample. We form four linear combinations: and calculate birefringence magnitude ␦ B , orientation ␤, and signal topography ␦ z with the formulas

⌺ B ͑S͒ ϭ ϪS 1 ͑S͒ ϩ S 2 ͑S͒ ϩ S 3 ͑S͒ Ϫ S 4 ͑S͒ , (1)
⌺ A ͑S͒ ϭ ϪS 1 ͑S͒ ϩ S 2 ͑S͒ Ϫ S 3 ͑S͒ ϩ S 4 ͑S͒ , (2)
⌺ B ͑P͒ ϭ ϪS 1 ͑P͒ ϩ S 2 ͑P͒ ϩ S 3 ͑P͒ Ϫ S 4 ͑P͒ , (3)
⌺ A ͑P͒ ϭ ϪS 1 ͑P͒ ϩ S 2 ͑P͒ Ϫ S 3 ͑P͒ ϩ S 4 ͑P͒ ; (4)
tan 2 ͩ ␦ B 0 ͪ ϭ tan Ϫ2 ͑2͒ ͓⌰ B ⌺ A ͑P͒ ͔ 2 ϩ ͓⌰ A ⌺ B ͑P͒ ͔ 2 ͓⌰ B ⌺ A ͑S͒ ͔ 2 ϩ ͓⌰ A ⌺ B ͑S͒ ͔ 2 , ( 5
) tan ͩ 2 ␦ z 0 ͪ ϭ ⌰ B ⌺ A ͑S͒ ⌰ A ⌺ B ͑S͒ (6) tan͑2␤͒ ϭ ⌰ A ⌰ B ͓⌺ A ͑S͒ ⌺ B ͑P͒ Ϫ ⌺ A ͑P͒ ⌺ B ͑S͒ ͔ ⌰ A 2 ⌺ B ͑S͒ ⌺ B ͑P͒ ϩ ⌰ B 2 ⌺ A ͑S͒ ⌺ A ͑P͒ , (7)
where, knowing the apparent spectrum of the source, we calculate ⌰ A and ⌰ B , and 0 is an equivalent wavelength of the polychromatic source.

Results

We calibrated our instrument by inserting a Babinet compensator into the sample arm between the quarter-wave plate and the microscope objective. We generated a well-defined birefringence signal with known axial directions and varied its amplitude from Ϫ͞2 to ϩ͞2. Figure 3 shows the measured birefringence as a function of the true retardation imposed by the compensator. The best linear fit is shown ͑solid line͒; its coefficient of linear regression is 0.99985, slope is 0.993 Ϯ 0.005, and intercept is Ϫ3 Ϯ 5 mrad, demonstrating the good agreement of the measurement with the true values. However, it must be noted that to obtain this result we found it necessary to take into account the imperfections of the quarter-wave plates. This correction is described in Section 5 below. We validated our setup with a multilayer optical coating sample whose structure is sketched in Fig. 4. For an optical system with a scattering specification level of a few parts in 10 6 ͑laser gyros or interferometric detectors of gravitational waves [START_REF] Mackowski | Different approaches to improve the wavefront of low-loss mirrors used in the VIRGO gravitational wave antenna[END_REF][START_REF] Vinet | Scattered light noise in gravitational wave interferometric detectors: a statistical approach[END_REF] ͒, point defects inside the multilayer coating are a major source of loss. Knowing the locations of these scattering defects, and furthermore using a nondestructive technique, will help in finding the source of contamination. This detection cannot be made by use of a tomographic image of the reflectance of the structure because the specular light from the interfaces would mask the small signal coming from the defect. However, if the scattering defect is nonspherical, it can change the polarization of the incident light. Therefore, in a tomographic image of the birefringence the point defect will appear as a bright spot in a dark field, as the interfaces are not birefringent. Our sample is an infrared interference coating designed for a CO 2 laser ͑ ϭ 10.6 m͒, with five layers of alternate high-and low-index dielectric material deposited upon a ZnSe substrate. The layers are alternately half-wave and quarter-wave plates. Figure 5 shows the interferogram envelope: ⍀͑z͒, [START_REF] Moreau | Full-field birefringence imaging by thermal-light polarization-sensitive optical coherence tomography. 1. Theory[END_REF] which in this case is similar to the classic OCT to-Fig. 3. Relation of estimated to true retardation of a Babinet compensator: squares, uncorrected data; crosses, data corrected for quarter-wave plate error following to Eq. ͑15͒. mography signal. It must be noted that the scale of the Z axis is given by the translation of the sample upon its motorized stage. To obtain true optical thickness values one should correct the Z scaling by an obliquity factor that depends on the objective's N.A. and on the sample's index of refraction. The axial resolution of our instrument enables us to resolve the six interfaces. We made a threedimensional birefringence map of this mirror ͑60 images; step, 0.4 m͒ to enable us to detect the presence of point defects inside the structure. Figures 6 and7, respectively, show XY cuts of the magnitude and direction of the birefringence axes about such defects. Figures 8 and9, respectively, show XZ and YZ cuts. In these cuts the noise was measured to be 2.9 ϫ 10 Ϫ4 rad. Figure 7 shows a measurement of the direction of the birefringence axes in the region where the defect was detected. To ensure quantitative measurement we evaluated the birefringence only when ⍀ took a significant value. A defect in birefringence can be seen and precisely located on the second interface. Its size is equal to or less than the resolution of our instrument in three dimensions. The spurious birefringence signal coming from the interfaces ͑ϳ0.1 rad͒ is due to a retardation error of the achromatic quarter-wave plate ͑see Subsection 5.D below͒. The time needed to acquire a complete three-dimensional map of a sample depends on the number of XY cuts as well as on the number of images averaged at each position of the sample. In the measurement presented, the sample was moved to 60 different Z positions, and an average of 400 sequences was performed at each position, so the total number of images acquired was 60 ϫ 8 ϫ 400. The time needed for the sample to move from one position to the next was 200 ms, and the total acquisition time was 17 min.

Performance and Limitations

A. Resolution and Position Accuracy

The axial resolution in air is evaluated with a single interferogram recorded by use of a mirror sample. From the measurement presented in Fig. 10 one can find the best Gaussian envelope of the interferogram and estimate the full-width at half-maximum. We found axial resolution ⌬ z ϭ 1.5 m in air. The axial resolution is limited by the apparent spectrum width of the light source presented in Fig. 11 ͑i.e., which takes into account the spectral response of the CCD array and of the optical components͒ and by uncompensated dispersion inside the interferometer and inside the sample. The uncompensated dispersion can in principle be corrected if one knows the dispersion curves and the thicknesses of the various optical elements and sample. The accuracy of the axial position depends on the algorithm that uses one to find the position of the maximum of the interferogram envelope and is limited by the axial step size of the sample displacement. In our setup the displacement step size is 74 nm, which can be taken as a good order of magnitude of the accuracy of the axial position.

The lateral resolution in air is limited by diffraction and depends on the objective lenses that are used. We use an optical scheme in which each pixel of the CCD array receives the image of a cross section of the sample approximately equal to the diffraction spot size, so the accuracy of the lateral position is equal to the resolution. When better lateral accuracy is required, expanding imaging optics can be added such that the diffraction spot is imaged on more than a single pixel.

B. Maximum Detectable Birefringence

The short coherence length limits the maximum magnitude of birefringence that can be detected: If the retardation induced by the sample is larger than the coherence length of the source, then, after interaction with the sample, the two eigenstates of polarization will lose their mutual coherence and the field returning from the reference arm will become a superposition of two orthogonally polarized incoherent fields. The maximum value of the birefringence that can be measured is max͑␦ B ͒ ϭ ⌬ z , or, expressed in terms of phase retardation and using the equivalent wavelength of the source,

max͑͒ ϭ 2 ⌬ z 0 . ( 8
)
In our setup this value is ϳ3.

C. Sensitivity

The performance of the instrument is limited by two factors, residual noise and imperfections and misorientation of the polarizing components, especially the quarter-wave plates and the beam-splitter cube. The former puts an upper limit on sensitivity and the latter, together with uncertainties in the modulation parameters, degrade measurement precision. The fundamental sensitivity limits of the instrument are discussed in depth in the paper by Dubois. [START_REF] Dubois | Phase-map measurements by interferometry with sinusoidal phase modulation and four integrating buckets[END_REF] We measured the signal-to-noise ratios ͑SNRs͒ of individual images and found that they are shot-noise limited up to an integration duration of 20 s. During 1-s integration we measured a SNR of 70 dB; the theoretical value of the shot-noise-limited SNR for a 13 ϫ 10 4 charge per pixel CCD is 73 dB. Usually in OCT-related papers the SNR is given as a minimum sample reflectance, measured with a glass reference surface. [START_REF] Dubois | Phase-map measurements by interferometry with sinusoidal phase modulation and four integrating buckets[END_REF] In that case a correction equal to Ϫ20 ϫ log 10 ͑ ͌ R glass ͒ ϭ 14 dB is added to the SNR value, which in the study reported here increases to 84 dB. When the duration of integration becomes large, other sources of noise, possibly mechanical, acoustic, and thermal, become dominant. The shape of the measured sensitivity-versus-time plot presented in Fig. 12 suggests that this technical noise may exhibit an f Ϫ1 frequency dependency, and its presence places an upper limit of 81 dB ͑95 dB in terms of minimum reflectance with an uncoated glass reference͒ on the SNRs of individual images. However, we are interested not in the SNRs of individual images but in the SNR of the birefringence measurement. The way to compute the SNR of the birefringence measurement, SNR͑͒, starting from the SNR of individual images, SNR͑S͒, is straightforward but lengthy, and we simply outline it here. We start our calculation by assuming that each image exhibits the same average value and noise level as the other images and that the noise of the various images is uncorrelated. These assumptions allow to calculate the SNRs of the various combinations of ⌺ A ͑S͒ , ⌺ A ͑P͒ , ⌺ B ͑S͒ , and ⌺ B ͑P͒ , using Eqs. ͑1͒-͑4͒, and then the SNR of the birefringence measurement, using Eq. ͑5͒.

There is a simple relation between the SNRs of individual images SNR͑S͒ and the SNR of the birefringence measurement SNR͑͒ that depends on the magnitude of the birefringence and on the orientation of the quarter-wave plate in the reference arm:

SNR͑͒ ϭ SNR͑S͒ ϫ g͑, ͒, (9)
where g is a known function of , , and the modulation amplitude. The plot of g in Fig. 13, calculated for ⌰ A ϭ ⌰ B , shows that the value of that gives the highest SNR for the broadest range of is ϳ11°͓the optimal value is 22.5°only when one is interested in measuring tan͑͞2͔͒. Figure 14 is a plot of g͑, 11°͒, which shows that SNR͑͒ is roughly ten times less than SNR͑S͒. We thus estimate the maximum SNR for a birefringence measurement to be 61 dB ͑50 dB for 1-s integration͒.

D. Precision

The most important limitation on instrument precision comes from the various imperfections of the optical components, in particular, their behavior with respect to polarization. The requirement of achromaticity is scarcely achievable in polarization components. For polarizers it appears as a poor extinction factor; for phase-shifting components, as spurious wavelength dependent retardation. To eliminate any displacement of the image on the CCD when we rotate the analyzer, we have to use a thin polarizer with an extinction two times larger than one could obtain with a Glan-type polarizer. Then, as was found experimentally, coupling of polarization states can occur inside the beam-splitter cube if the cube is slightly misaligned with respect to the axes of the two arms. It is rather difficult to overcome this problem because it does not show in pure intensity images and is hardly distinguishable a priori for misorientations of the polarizer or the analyzer. The optical components that cause the largest imperfections are the two achromatic quarter-wave plates, as we emphasized above. A simple calculation can be performed for a monochromatic source in the presence of spurious retardation of ␦ Q waves in the quarter-wave plates. In this case the Jones matrix of the quarterwave plates reads as

M Q ϭ ΄ exp ͫ i ͩ 4 ϩ ␦ Qͪ ͬ 0 0 exp ͫ Ϫi ͩ 4 ϩ ␦ Qͪ ͬ ΅ . (10)
The apparent birefringence of the sample, app , calculated with the standard equation that is valid for the monochromatic source, is related here to the true value of the birefringence by

tan 2 ͩ app 2 ͪ ϭ sin 2 ͑͞2͒ ϩ h͑, ␤, ␦ Q ͒ cos 2 ͑͞2͒ Ϫ h͑, ␤, ␦ Q ͒ , (11)
with

h͑, ␤, ␦ Q ͒ ϭ ͓cos 2 ͑͞2͒ Ϫ sin 2 ͑2␤͒sin 2 ͑͞2͔͒sin 2 ͑2␦ Q ͒ Ϫ sin͑͒sin͑2␤͒cos͑4␦ Q ͒. (12)
The true value of ␦ Q is an unknown function of wavelength that depends on the structure of the quarterwave plate used. Its apparent value can nevertheless be measured by use of a nonbirefringent sample; then the apparent birefringence is

tan 2 ͩ app 2 ͪ ϭ sin 2 ͑2␦ Q ͒ 1 Ϫ sin 2 ͑2␦ Q ͒ Ϸ 4␦ Q 2 1 Ϫ 4␦ Q 2 . ( 13
)
We found experimentally that ␦ Q ϭ 25 mrad; the value of h can then be estimated from the apparent values of app and ␤ app . It is worth noting that h does not depend on topography ␦ z . To calculate h we need the sign of sin͑ app ͒ ϫ sin͑2␤ app ͒. It can be found because, when 0 Յ Յ ͞4, one combination of the four signals ⌺ A S , ⌺ B S , ⌺ A P , and ⌺ B P has the same sign:

sign͓sin͑͒sin͑2␤͔͒ ϭ sign͓⌺ A ͑P͒ ⌺ B ͑S͒ Ϫ ⌺ A ͑S͒ ⌺ B ͑P͒ ͔, 0 Յ Յ ͞4. (14)
The true value of can then be estimated by

tan 2 ͩ 2 ͪ Ϸ cos 2 ͑ app ͞2͒ Ϫ h͑ app , ␤ app , ␦ Q ͒ sin 2 ͑ app ͞2͒ ϩ h͑ app , ␤ app , ␦ Q ͒ . ( 15
)
The calibration plot of Fig. 3 was obtained with relation ͑15͒.

Conclusion

We have constructed a polarization-sensitive Linniktype interference microscope that works with a thermal light source. Our system uses full-field illumination with sinusoidal phase modulation and four integrating buckets. This system can produce en face ͑x-y͒ tomographic images of the birefringence with a resolution of 1 m in the transverse direction and 1.5 m in the axial direction without scanning. Quantitative measurements can be obtained for birefringence from 0.1 to 3 rad; these limits are fixed by the imperfections of the achromatic quarter-wave plates and the coherence length of the source. The maximum sensitivity of our setup, in terms of minimum detectable sample reflectance, was found to be 95 dB ͑84 dB for 1-s integration͒ for an individual image and ϳ61 dB ͑ϳ50 dB for 1-s integration͒ for a birefringence measurement. We used our system as a dark-field microscope to detect and localize micrometer-sized defects inside a multilayer optical coating infrared mirror. Our primary interest was to gain a better understanding of the locations of scattering defects in multilayer optical coatings and therefore to determine the sources of contamination. Nevertheless, we are also currently investigating ways in which to use our microscope to study biological samples.

Fig. 1 .

 1 Fig. 1. Schematic of the instrument described in this paper: NPBS, nonpolarizing beam-splitter cube; AQWPs, achromatic quarter-wave plates; PZT, piezoactuated translation stage.

Fig. 2 .

 2 Fig. 2. Procedure used to fix the value of ␦ 0 : interferograms with ͑a͒ ␦ 0 ϭ 0.85 ϫ ␦ 0opt , ͑b͒ ␦ 0 ϭ ␦ 0opt , and ͑c͒ ␦ 0 ϭ 1.15 ϫ ␦ 0opt .

Fig. 4 .

 4 Fig. 4. Schematic of the multilayer optical coating mirror used as a sample.

Fig. 6 .Fig. 7 .

 67 Fig.6. En face birefringence image of the multilayer second interface on which a defect is visible: XY cut, 0.7 m ϫ 0.7 m resolution. The noise in this image is estimated to be 2.9 ϫ 10 Ϫ4 rad.

8 .

 8 High-resolution ͑0.7 m ϫ 1.5 m, lateral ϫ axial͒ birefringence image of the multilayer; XZ cut; same noise level as in Fig.6. Dashed lines, positions of the six interfaces.

Fig. 9 .

 9 Fig. 9. YZ cut.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Measured axial response of the system and best Gaussian envelope. The FWHM is 1.5 m.

Fig. 12 .

 12 Fig.12. Image SNR as a function of the number of images averaged. Circles, measurements; solid line, shot-noise limit; dashed curve, estimate of the SNR if an f Ϫ1 noise source is present; horizontal dashed-dotted line, SNR limit imposed by the f Ϫ1 noise.

Fig. 13 .

 13 Fig. 13. Value of g as a function of birefringence magnitude and orientation of the quarter-wave in reference arm .

Fig. 14 .

 14 Fig. 14. Value of g as a function of birefringence magnitude for ϭ 11.26°.

July

͞ Vol. 42, No. 19 ͞ APPLIED OPTICS

The authors thank Arnaud Dubois and Laurent Vabre for helpful discussions.


J Eynard

S Grieu
email: grieu@univ-perp.fr

M Polit
email: polit@univ-perp.fr

Optimal control of a multi-energy district boiler: a case study ⋆

Keywords: Modelling, Identification, Prediction, Neural networks, Predictive control, Boiler, Energy storage, Renewable energy systems

The present paper deals with the optimization of a multi-energy district boiler (La Rochelle, France) which supplies domestic hot water and heats residential and public buildings, using renewable and fossil resources. First, a combination of white, grey and black boxes was used to model the plant, thanks to a modular approach. Next, a stratified thermal storage tank was modelled and completed the just-mentioned plant model. Using these models and forecasted sequences about outdoor temperature and thermal power consumption, a model predictive controller allows optimizing the use of both the tank and the wood boiler. As a result, fossil energy consumption and CO 2 emissions are minimized. Energy is stored during low-demand periods and used when demand is high, instead of consuming gas and fuel oil.

INTRODUCTION

By the end of the first commitment period of the Kyoto Protocol in 2012, the European Union has decided for 2020 to promote renewable energy (EU ( 2009b)) (20 %), to reduce GreenHouse Gases (GHG) emissions (EU (2009a)) (-20 %) and to increase energy efficiency (+20 %). As part of the OptiEnR research project [START_REF] Eynard | Gestion optimale de l'nergie dans un procédé multi-source pour le chauffage de bâtiments[END_REF]), the present paper deals with optimizing a multi-energy district boiler, situated at La Rochelle (west coast of France) and managed by Cofely GDF-Suez. The plant supplies domestic hot water and heats residential and public buildings. It is described in section 2 and uses renewable (wood) and fossil energies (gas and fuel oil). The main objective of the project was optimizing the district boiler functioning while minimizing the consumption of fossil energy. This has been completed according to 4 successive tasks (Fig. 1). The first task (section 3) was to model the district boiler. Because of both the complexity of the plant as a whole and the strong interactions between the sub-systems, a modular approach was proposed. To carry out that task, measurements campaigns results as well as several modelling approaches, such as black-box, parametric or knowledge modelling approaches, were considered. Because energy storage is a way to optimize the performance of a district boiler, we studied the impact of adding a stratified hot water tank to the plant of La Rochelle (section 4). Next, a Model-based Predictive Controller (MPC) was designed to optimize both the flow of the water passing through the tank and the wood boiler set-point temperature (section 6) minimizing a cost function based on fuel oil consumption. The controller uses the developed overall model and forecasted outdoor temperature and thermal power consumption sequences (section 5). Section 6 focuses on the way the MPC was tuned and on simulation results. One ⋆ Industrial partner: Cofely GDF-Suez Fig. 1. The OptiEnR research project can highlight a significant decrease of the fossil energy consumed, leading to a reduction in costs and CO 2 emissions.

DISTRICT BOILER OF LA ROCHELLE

The district boiler of La Rochelle is a complex plant (Fig. 2). The collecting hydraulic circuit connects the wood (4.5 MWh) and the gas-fuel oil (7 MWh) boilers, using two feed pumps (F p W B and F p GF B ). Due to the low price of wood, the wood boiler functions all the time while the gas-fuel oil boiler is used during very cold periods only, when the wood boiler fails to respond to the hot water demand. First, the gas-fuel oil boiler is fed using gas. When the 6000 m 3 of gas allowed every day to Cofely GDF-Suez are consumed and if the boiler must continue to run, fuel oil is used until the end of the day (in this case, gas is much more expensive than fuel oil). The primary hydraulic circuit (3000 m 3 ) supplies hot water to heat residential and public buildings, for a total of 2700 accommodations. Domestic hot water is also produced, for a total of 3500 accomodations. The circuit is composed of a hot water distribution network, a feed pump (F p DN ) used to control the network differential pressure, and a cogeneration plant which produces electricity. A part of the cold water coming back from the distribution network is warmed up by the cogeneration plant. Using measurements of both the Fig. 2. Synopsys of the district boiler with the proposed modifications (in red, green, blue and purple)

outdoor temperature and the temperature of the water entering the distribution network, a monitoring system allows tracking the functioning of the two boilers. Finally, the breaking pressure bottle pulls apart the two justmentioned hydraulic circuits, because of the difference between their respective flows.

DISTRICT BOILER MODELLING

Overall approach

Taking a look at the state of the art about the modelling of district boilers, one can highlight some interesting works. First, [START_REF] Curti | An environomic approach for the modeling and optimization of a district heating network based on centralized and decentralized heat pumps, cogeneration and/or gas furnace. Part I: Methodology[END_REF] proposed an environomic approach for modelling and optimizing a district heating network. [START_REF] Dias | Energetic and exergetic analysis in a firewood boiler[END_REF] modelled a wood boiler using a detailed energy and exergy analysis, allowing quantifying the types, causes, and locations of losses. Physical principles and artificial neural networks can be used together as it has been shown by [START_REF] Lu | Dynamic nonlinear modelling of power plant by physical principles and neural networks[END_REF] for the development of a nonlinear power plant model. Finally, [START_REF] Ghaffari | Soft computing approach for modeling power plant with a once-through boiler[END_REF] presented a soft computing approach to model electrical power generating plants and to characterize the essential dynamic behavior of the plant sub-systems. Their approach consisted of fuzzy logic, artificial neural networks and genetic algorithms. As previously mentioned and because of both the complexity of the district boiler of La Rochelle as a whole and the strong interactions between the sub-systems, a modular approach was proposed to model the plant (Eynard et al. (2011a)). According to the availability of information (as a result of measurement campaigns or taking into consideration expert knowledge about the sub-systems functioning and the boilers control systems), a combination of white, grey and black boxes was used to keep the modeling process on track. With white boxes which can be used when one can easily describe the interactions between physical parameters, experimental data serve only to validate the models. Both grey boxes [START_REF] Sjöberg | Nonlinear black-box modeling in system identification: a unified overview[END_REF]) (in this case, one needs to find appropriate model inputs and outputs, thanks to physical considerations and/or analyzing the process dynamic properties, while experimental data Fig. 3. Modular approach to model the district boiler are used to optimize both the model topology and parameters [START_REF] Ljung | System identification: theory for the user[END_REF]) as well as validation data) and black boxes (in this case, no physical considerations are taken into account, a standard topology is used and, as when using grey boxes, one needs to find appropriate model inputs and outputs and to optimize, with experimental data, the model parameters) are used in case of missing information about physical behaviours.

Modelling process

First, a data processing phase allowed eliminating abnormal data and upsampling data with a too-high sampling rate, considering a desired sampling rate of 5 minutes.

Let us also note that some relevant but not-measured parameters were estimated thanks to appropriate information. Next, for each parameter to be modelled, both the model topology (a differential, algebraic or logical equation) and inputs were defined. A first iterative process of optimisation was performed to identify the coefficents of the models. We used the trust-region reflective Newton method [START_REF] Coleman | On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds[END_REF]) for nonlinear leastsquares [START_REF] Kelley | Iterative methods for optimization[END_REF]) to solve the minimization problem.

Analyzing the results allows validating the chosen topologies (if a result is not satisfactory, one needs to think about another model topology and of course to find the right equation coefficients again). With the boiler parameters correctly described, using optimized differential, algebraic or logical equations, sub-models were defined. It can be noticed that the black-box methodology used for the modelling of the wood boiler was based on the Hammerstein-Wiener theory [START_REF] Zhu | Estimation of an N-L-N Hammerstein-Wiener model[END_REF]). A second iterative process of optimisation was carried out to adjust the sub-models coefficients. Finally, when the sub-models accuracy was as good as possible, these models were combined to obtain the overall district boiler model. Fig. 3 depicts the modular approach, leading to very good results.

ENERGY STORAGE MODELLING

Hypothesis

Adding to the district boiler of La Rochelle an energy storage unit allows using, when demand is high, the excess of energy produced by the wood boiler during low-demand periods, instead of consuming gas and fuel oil. As a result, the coverage rate of the fossil energy used can be significantly reduced thanks to a better exploitation of the available renewable resources. Many options about the materials one can use to store energy were considered but due to technical reasons, their cost, the constraints related to their integration into the district boiler, their toxicity, their flammability or their relative low power, phase change materials were not used. So, we decided for a hot water tank. Some hypothesis about its shape, size and position were made. We proposed a vertical cylindrical tank whose diameter is equal to its height with the aim of minimizing its surface and, as a consequence, thermal losses. For the same reason, we supposed the tank, whose thermal insulation is insured by 10 cm of polypropylene, to be buried in the ground, where ambiant temperature (T amb ) is about 10 ℃. The inlet fluid temperature (T in ) affects the first layers of water in the tank only, because of a grid used to protect its thermal stratification.

Modelling results

After taking a look at the review of [START_REF] Han | Thermal stratification within the water tank[END_REF] dedicated to the modelling of one-dimension hot water tanks, i.e. with vertical thermal stratification, we developed a discrete space model. The model is adapted from the model proposed by [START_REF] Alizadeh | An experimental and numerical study of thermal stratification in a horizontal cylindrical solar storage tank[END_REF] and defined by equations ( 1) and ( 2) for energy storage and destocking.

T i,k is the temperature of the i th layer of water at time index k, g is the number of layers, e is the number of layers directly affected by the inlet fluid temperature, V is the volume of one layer of water while ∆V is the volume of fluid passing through the tank during one sampling time.

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (i > g -e) ∶ T i,k+1 = V -∆V e ⋅T i,k + ∆V e ⋅T in k +βi⋅T amb k αi (i ≤ g -e) ∶ T i,k+1 = (V -∆V )⋅T i,k +∆V ⋅Ti+1(k)+βi⋅T amb k αi (1) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (i ≤ e) ∶ T i,k+1 = V -∆V e ⋅T i,k + ∆V e ⋅T in k +βi⋅T amb k αi (i > e) ∶ T i,k+1 = (V -∆V )⋅T i,k +∆V ⋅T i-1,k +βi⋅T amb k αi
(2)

Integration of the storage unit into the district boiler

The characteristics of the thermal storage unit being defined, one may think about the way the tank would be integrated into the district boiler. That is why we proposed an hydraulic modification of the plant, as shown in red on Fig. 2. We decided to place the tank between the breaking pressure bottle and the hot water distribution network, before the feed pump of the distribution network and after the cogeneration plant. A new pump (F p T S ) is used to control the flow of the water passing through the tank. To reflect the impact of the thermal storage unit on the district boiler, some subsidiary models were developed, as shown by Fig. 4. Indeed, some flows between the breaking pressure bottle and the thermal storage tank changed while some temperatures were affected by the tank functioning.

EXOGENOUS VARIABLES FORECASTING

Materials and methods

Both the outdoor temperature (T out ) and the thermal power consumption of the hot water distribution network One could speak of "MRA-ANN" methodology, described in detail in [START_REF] Eynard | Waveletbased multi-resolution analysis and artificial neural networks for forecasting temperature and thermal power consumption[END_REF]. The discrete wavelet transform allowed decomposing M sequences of past data (of l points) in subsequences, according to different frequency domains. We used a filter bank, composed of lowpass and high-pass filters, to obtain approximation (i.e. low frequency components) and detail (i.e. high frequency components) coefficients, according to a decomposition level (N ) and a wavelet order (R) (Fig. 5). From these coefficients, Multi-Layer Perceptrons (MLP) were used to estimate future subsequences of 4h30. As depicted by Fig. 6, N + 1 neural networks were used to estimate the N details and the approximation of level N of a target sequence P (the sequence to be forecasted). Future values of outdoor temperature or thermal power consumption were then obtained by simply summing up the estimated coefficients (Fig. 6). Substituting the prediction task of an original time series of high variability by the prediction, using MLP neural networks, of its wavelet coefficients on different levels of lower variability, is the main idea of the methodology. The wavelet-based multi-resolution analysis allowed isolating the global trend and the 24-hour pseudoperiod characterizing the considered time series from the variability caused by climatic phenomena. To place the model in time, the sequences of past data were completed, for each of their components, by the minute of the day and the day of the year. According to the respective values of M , l and of the sub-sampling time T SST (set initially to 30 minutes; forecasted sequences were then upsampled to be in agreement with the sampling time of all the parameters measured at the district boiler), the results can be more or less accurate. Indeed, the further the forecasting horizon (set to 4h30, as requested by the boiler operators; this leads to l=9), the more inaccurate it can be. Inaccuracy

T out and W DN forecasting results

Table 1 depicts the optimal configuration used to forecast T out and W DN (the parameters are inversely proportionnal), according to the wavelet order (R), the wavelet decomposition level (N ), the common topology of the networks used (F ), the number of considered past sequences (M ) and, finally, the sub-sampling time used (T SST ). Table 2 (lines 1 and 2) specifies, for the considered period (from February 15, 2009to March 02, 2009), the Mean Relative Error (MRE), the Mean Absolute Error (MAE) and the curve fitting (FIT) observed when forecasting T out and W DN , using the MRA-ANN methodology.

OPTIMAL CONTROL OF THE DISTRICT BOILER

Control strategies

First, some criteria to be minimized were developed and used to optimize the performance of the district boiler. Equations (3-8) define J 1 as the cost related with fuel oil consumption, J 2 as the positive difference between the temperature of the water entering the distribution network (T E DN ) and the temperature setpoint (T SDN ), J 3 as the energy consumed (in LHV), J 4 as the fossil energy coverage rate, J 5 as the level of CO 2 emissions and J 6 as the overall performance of the district boiler with a controlled thermal storage unit compared with its performance without energy storage. Table 3 presents the Unitary Cost (U C), the Unitary Low Heating Value (U LHV ) and the Unitary CO 2 emissions in Life-Cycle Assessment (U LCA) for each of the fuels used (ADEME ( 2007)). Let us note that we used for wood the number of tappet strokes per five minutes (the mean weight of the wood introduced into the boiler by the tappet is about 37 kg).

J 1 = U C wood ⋅ N tappet + U C gas ⋅ V gas + U C f uel-oil ⋅ V f uel-oil (3)
J 2 = 1 2 ⋅ N ⋅ N k=1 T SDN (k) -T E DN (k) -T E DN (k) -T SDN (k) (4)
J 3 = U LHV wood ⋅ N tappet + U LHV gas ⋅ V gas + U LHV f uel-oil ⋅ V f uel-oil (5)
J 4 = U LHV gas ⋅ V gas + U LHV f uel-oil ⋅ V f uel-oil J 3 (6) J 5 = U LCA wood ⋅ N tappet + U LCA gas ⋅ V gas + U LCA f uel-oil ⋅ V f uel-oil (7)
J 6 = 1 5 ⋅ 5 n=1 100 ⋅ J n (with thermal storage) J n (without thermal storage) (8)
Minimizing one of the just-mentioned criteria defines the control strategy used for optimizing the district boiler performance. According to economic, technical, environmental or related-to-energy-efficiency considerations, one can choose one of the first five criteria we proposed. With the criterion J 6 , one can optimize the district boiler performance in an unspecific way, according to all the justmentioned considerations.

Model-based optimal predictive controller

The model predictive controller developed to optimize the district boiler performance uses the models of both the plant and the thermal storage unit as well as forecasted sequences of T out (the outdoor temperature) and W DN (the thermal power consumption of the hot water distribution network). It computes optimal command sequences to be applied to the considered process over a prediction horizon, according to an objective function to be minimzed.

Taking a look at the literature, one can note that through MPC is not inherently more or less robust than classical feedback, it can be adjusted more easily for robustness [START_REF] García | Model predictive control: Theory and practice-a survey[END_REF]). Moreover, it has good properties of stability [START_REF] Mayne | Robust timeoptimal control of constrained linear systems[END_REF]) and can be used to control nonlinear systems [START_REF] Kambhampati | A stable one-step-ahead predictive control of non-linear systems[END_REF]), hybrid systems [START_REF] Olaru | Predictive control for hybrid systems. implications of polyhedral pre-computations[END_REF]) or even fast systems [START_REF] Pannocchia | Fast, large-scale model predictive control by partial enumeration[END_REF]). Such a controller is commonly used to control industrial processes in real-time [START_REF] Qin | A survey of industrial model predictive control technology[END_REF]). The wood boiler set-point temperature (T SW B ) and the flow of the water passing through the hot water tank (Q T S ) are the two controlled parameters (Fig. 2). Let us note that information related to measured data is given to the MPC by the previously-implemented control schemes. Equation ( 9) shows the command vectors to be applied to the process and to be optimized thanks to the minimization of the chosen criterion J 1 . Equation ( 10) defines the path constraints, with p and c the prediction and control horizons, respectively. The thermal storage unit can not store hot water when the gas-fuel oil boiler is engaged (Eng GF B = 1) while thermal stratification must always be ensured. As a consequence, the temperature of the water leaving the breaking pressure bottle (T L BP B ) has to be higher than the temperature of the upper layer of the tank (T sup-T S ). The resolution of the nonlinear optimization problem was carried out using a pattern search algorithm [START_REF] Sherif | Optimization by pattern search[END_REF]). Such a global optimization method is well adapted for nonlinear problems with unknown gradient fonctions. It is commonly used for the resolution of complex industrial problems, for example a constrained optimization problem in a way that insures a secure-economic system operation (Al-Othman and El-Naggar (2008)).

Q T S (k k), . . . , Q T S (k + c -1 k) T SW B (k k), . . . , T SW B (k + c -1 k) (9)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩
Overall district boiler with thermal storage model whose inputs are T 4h30 out and W hyb DN

Q T S ≤ Q DN -10 (Eng GF B = 1) ∨ (T L BP B ≤ T sup-T S ) ⇒ Q T S ≤ 0 90 ℃ ≤ T SW B ≤ 97 ℃ ∆Q T S (k + h) = 0 ∀h ∈ c; p -1 ∆T SW B (k + h) = 0 ∀h ∈ c; p -1 (p, c) = (54, 36) (10)

Optimal predictive control results

The proposed controller was tested in simulation during 45 days, from mid-January to early March. We focused on the impact of the hot water tank volume (V tot ), ranging between 500 m 3 and 4000 m 3 . Table 4 highlights three remarkable configurations for the proposed control scheme (M P C 1 , M P C 2 and M P C 3 defined according to three different tank volumes and allowing minimizing J 1 ), taking as a reference the modelled functionning of the district boiler without thermal storage unit and MPC controller. Taking a look at Table 4, one can note that configuration M P C 1 reduces the fuel oil consumption by 37 % (-1244 l), the set-point tracking error (J 2 ) by 55 % (-0.313 ℃), the CO 2 emissions (J 5 ) by 40 % (-56.2 tCO 2 ), and the performance criterion J 6 by 36 %. With configuration M P C 2 , the wood consumption, the economic cost (J 1 ) and the primary energy used are reduced by 14 %, 24 % (-14.55 ke) and 18 % (-682 MW⋅h), respectively. Finally, configuration M P C 3 minimizes the gas consumption by 55 % (-20 794 m 3 ) as well as the fossil energy coverage rate (J 4 ) by 42.5 % (-4.93 points).

CONCLUSION

The present paper deals with optimizing the performance of a multi-energy district boiler which supplies hot water via a distribution network. First, black, grey and white boxes were used to model the considered process, according to a modular approach. Such an approach was the consequence of both the complexity of the plant as a whole and the strong interactions between the subsystems. Next, a stratified hot water tank, connected to the primary circuit, was modelled and completed the district boiler model. Because one needs, using the overall model proposed, forecasted sequences of both the outdoor temperature and the thermal power consumption of the distribution network (both parameters are the overall model inputs) to optimize the use of the tank on-line, a prediction tool was also developed. It is based on a multi-resolution analysis and artificial neural networks. Finally, a model predictive controller was developed to optimize, over a prediction horizon, the use of both the thermal storage unit and the wood boiler. It uses the overall model and the just-mentioned forecasted sequences. One can highlight that the proposed control scheme allows reducing the fossil energy consumption significantly and, as a result, functioning costs and CO 2 emissions. Future work will focus on improving both the district boiler model and the storage process, using phase change materials as well as another way to implement the hot water tank. One can also think about alternative control schemes, such as fuzzy or neural control schemes, before implementing in-situ the developed tools.

Fig. 4 .

 4 Fig.4. Integration of the model of the thermal storage unit (W DN ) are the district boiler model inputs. One needs, using the model, forecasted sequences about these two exogenous parameters to optimize on-line the use of the thermal storage tank. During the last years, schemes for temperature forecasting were mainly developed in the context of load forecasting and power utilities management. The complex and nonlinear nature of temperature variations as well as the abundance of historical data suggest that computational intelligence data-based modeling techniques would be good candidates to forecast temperature[START_REF] Neto | Comparison between detailed model simulation and artificial neural network for forecasting building energy consumption[END_REF]). The proposed methodology, adapted from the methodology developed by[START_REF] Tran | Wavelets decomposition and Box-Jenkins methodology or neural networks for shortterm electric consumptions forecast[END_REF], deals with the concept of time series, even if not only past values are considered for estimating future values, and uses a wavelet-based Multi-Resolution Analysis (MRA) and Artificial Neural Networks (ANN). One could speak of "MRA-ANN" methodology, described in detail in[START_REF] Eynard | Waveletbased multi-resolution analysis and artificial neural networks for forecasting temperature and thermal power consumption[END_REF]. The discrete wavelet transform allowed decomposing M sequences of past data (of l points) in subsequences, according to different frequency domains. We used a filter bank, composed of lowpass and high-pass filters, to obtain approximation (i.e. low frequency components) and detail (i.e. high frequency components) coefficients, according to a decomposition level (N ) and a wavelet order (R) (Fig.5). From these coefficients, Multi-Layer Perceptrons (MLP) were used to estimate future subsequences of 4h30. As depicted by Fig.6, N + 1 neural networks were used to estimate the N details and the approximation of level N of a target sequence P (the sequence to be forecasted). Future values of outdoor temperature or thermal power consumption were then obtained by simply summing up the estimated coefficients (Fig.6). Substituting the prediction task of an original time series of high variability by the prediction, using MLP neural networks, of its wavelet coefficients on different levels of lower variability, is the main idea of the methodology. The wavelet-based multi-resolution analysis allowed isolating the global trend and the 24-hour pseudoperiod characterizing the considered time series from the variability caused by climatic phenomena. To place the model in time, the sequences of past data were completed, for each of their components, by the minute of the day and the day of the year. According to the respective values of M , l and of the sub-sampling time T SST (set initially to 30 minutes; forecasted sequences were then upsampled to be in agreement with the sampling time of all the parameters measured at the district boiler), the results can be more or less accurate. Indeed, the further the forecasting horizon (set to 4h30, as requested by the boiler operators; this leads to l=9), the more inaccurate it can be. Inaccuracy

Fig. 5 .

 5 Fig. 5. Wavelet-based multi-resolution analysis

Table 2

 2 DN . Although the proposed methodology provided better results than linear curve fitting, it underestimated sometimes peaks of consumption. Both approaches are valid but have different characteristics. That is why they were used in tandem to improve the forecasts accuracy.

(line 3) also presents the results obtained when using T out to forecast W DN , via a simple linear curve fitting. Finally, a hybrid approach, combining the MRA-ANN methodology and linear curve fitting (Table

2

(line 4)), was used to fore-cast W

Table 1 .

 1 Optimal configuration

	Parameter	Symbol Optimal value
	Wavelet order	R	4
	Wavelet decomposition level	N	5
	Common number of hidden neurons	F	5
	Number of past sequences	M	4
	Sub-sampling time	T SST	30 mn

Table 2 .

 2 T out and W DN forecasting results

	Variable	FIT	MRE	MAE
	T 4H30 out	60.6 % 4.14 %	1.15 ℃
	W 4H30 DN W T 4H30 out DN W hyb DN	44.7 % 6.07 % 663.6 kW 38.4 % 6.94 % 758.3 kW 46.9 % 5.85 % 639.1 kW

Table 3 .

 3 Characteristics of the fuels used

		Wood [-] Gas m 3	Fuel oil [l]
	U C [e]	1.8648	0.378	0.40
	U LHV [kW⋅h]	133.2	10.5	9.76
	U LCA [kgCO 2 ]	1.7316	2.283 84	2.928

Table 4 .

 4 Remarkable configurations

	Criteria		Configurations
	Name	Unit	Reference	M P C 1	M P C 2	M P C 3
	Vtot	m 3	0	1000	1500	4000
	Ntappet	-	24 635	21 260	21 203	21 482
	Vgas	m 3	38 001	17 553	17 533	17 207
	V f uel	l	3287	2043	2244	2436
	J 1	ke	61.62	47.10	47.07	47.54
	J 2	℃	0.56	0.247	0.252	0.272
	J 3	MW⋅h	3712	3036	3030	3066
	J 4	%	11.6	6.73	6.80	6.67
	J 5	tCO 2	139.1	82.9	83.3	83.6
	J 6	%	100	64.0	64.3	65.2