Rémy Ul
email: remy.ul@insa-cvl.fr

Rémi Rouffaud

Louis-Pascal Tran-Huu-Hue

Franck Levassort

Mai Pham-Thi

Claire Bantignies

Complete Electroelastic Set of co doped Barium Titanate for Transducer Applications

Keywords: barium titanate, piezoelectricity, functional characterization, piezocomposite

BaTiO 3 ceramic doped with cobalt and calcium (BTCaCo), processed by a conventional technique with a reduced sintering temperature, is characterized. Particular attention is paid to the delivery of complete and accurate set of elastic, dielectric, and piezoelectric constants. These databases are relatively rare in the literature but are essential for numerical simulations of new devices, integrating these materials. The processed piezoelectric material exhibits an electromechanical thickness coupling factor (k t ) of 45% and a relative dielectric constant of 1180, at constant strain. A 1-3 piezo-composite is successfully fabricated to improve the performance (k t =47%). Finally, BTCaCo disc is used to fabricate a transducer with a center frequency at 4.3 MHz, showing competitive performance compared to a standard PZT-based transducer.

I. INTRODUCTION

Barium titanate (BaTiO 3 ) is a well-known ferroelectric material that was intensively studied in the early 50's [1]. Indeed, BaTiO 3 has been rapidly substituted by lead zirconate titanate (PZT), which exhibits higher piezoelectric properties. However, PZT material contains lead and their increasing success is associated to health and environmental problems. Therefore, efforts have been devoted to the development of competitive lead-free counterparts and consequently, BaTiO 3 compositions have comeback on center stage as promising candidates. Various papers have been reported on BaTiO 3 based materials, processed by employing conventional solidstate techniques, with improved piezoelectric properties [START_REF] Li | Enhanced ferroelectric properties in (Ba1-xCax)(Ti0.94Sn0.06)O3 lead-free ceramics[END_REF]. Here, barium titanate co-doped with cobalt/calcium (BTCaCo) is prepared by conventional solid-state reaction method. Three samples with different geometries including pellets with a thickness of 0.6mm and a diameter of 13.5mm are fabricated. In order to obtain the complete set of material constants, resonance-antiresonance method is first used. In this case, measured and calculated values are mixed. This can lead to violations in the interrelations [START_REF] Ikeda | Fundamentals of Piezoelectricity[END_REF] between particular groups of electroelastic moduli and non-consistent database is provided. In this study, a method based on a genetic algorithm is used to optimize the consistency of the full set of material constants through the minimization of a defined criterion. This method has already been applied to KN single crystal characterization [START_REF] Rouffaud | Complete electroelastic set for the (YXt)-45° cut of a KNbO3 single crystal[END_REF]. With a BTCaCo disc, a single element tranducer is fabricated and the measured performances are compared with the theoretical behavior to evaluate the efficiency of this leadfree composition. Then, 1-3 piezocomposite is fabricated to evaluate the possible improved properties in this configuration for transducer applications. The following section describes the fabrication of ceramic and composite samples. Section III is devoted to the database of BTCaCo material. Finally, the last part (section IV), gives the electromechanical performance of the 1-3 piezocomposite and acoustic properties of the transducer.

II. FABRICATION AND ELECTRICAL MEASUREMENTS

A. Ceramic samples fabrication

Co-doped barium titanate is fabricated using oxide and carbonate powder as raw materials: BaCO 3 , TiO 2 , Co 3 O 4 , CaCO 3 and Li 2 O. The calcium dopant is introduced in substitution to the barium and the cobalt in substitution to the titanium.

According to the chemical formula (Ba 0.99 Ca 0.01 )(Ti 0.9925 Co 0.0075 )O 3-δ the powders are weighed, mixed, and then ball-milled by attrition in water for 1 h. The mixture is then dried calcined at 1100 °C for 2 h. Thereafter, calcined powders are ball-milled by attrition for 0.5 h, dried, and then mixed with 1 mol% of Li 2 O powder and with polyvinyl butyral and dibutyl phthalate used as a binder and a plasticizer. The resulting powder is then pressed by uniaxial pressing into a cylinder of 16.3 mm in diameter and 10-10.5 mm in thickness. After burning the binder and plasticizer at 600 °C for 0.5 h, the pressed cylinder is sintered at 1100 °C for 4 h in air. After sintering, the ceramic is cut into either several 0.6 mm thickness pellets or bars of 9 mm² in surface and 9 mm in height or platelets of 3 mm width 12 mm length and 0.55 mm thick. The geometries of the samples are designed to ensure different modes of vibration of the ceramic, such as thickness or planar (Fig. 1) and to give insight to the different constants associated with that modes.

The samples are coated with silver paint and fired at 500 °C for 1 h. The coated samples are then poled in a silicone oil bath by Field Cooling, consisting with the applying of a DC electric field of 600 V/mm during the cooling process of the samples from 140-150 °C to 40-50 °C.

B. 1-3 Piezocomposite fabrication

Dice and fill method [START_REF] Savakus | PZT-epoxy piezoelectric transducers: a simplified fabrication procedure[END_REF] is used for the fabrication of 1-3 piezocomposite from a BTCaCo disc. Dicing conditions are tuned for this material which has slightly higher elastic parameters than conventional lead-based ceramics. A ceramic volume fraction of 50% is chosen with a kerf around 30 µm. Samples are lapped to reduce the thickness to reach 400 µm. Finally, gold electrodes are sputtered on both faces (Fig. 2).

C. Electrical measurements

Measurements of the complex electrical impedance around the fundamental resonance of the favor mode (i.e. IEEE standard [START_REF]ANSI/IEEE Standard on Piezoelectricity[END_REF] are performed with an Agilent 4395A or 4294A impedance spectrum analyzer (Agilent Technologies Inc., Palo Alto, CA) and an impedance test kit. For the three BTCaCo samples, resonance and antiresonant frequencies are measured from maximum values of electrical admittance and impedance, respectively. According to relations given in [START_REF] Ikeda | Fundamentals of Piezoelectricity[END_REF], thickness (k t ) and planar (k p ) coupling factors are deduced with sample 1, k 33 with sample 2 and k 31 with sample 3 (Fig. 1). For the characterization of the 1-3 piezocomposite, the KLM equivalent electrical circuit [START_REF] Leedom | New equivalent circuits for elementary piezoelectric transducers[END_REF] is used to compute the theoretical behavior of the electrical impedance and a fitting process to deduce the thickness-mode parameters of the sample (Fig. 2). In this case, five parameters are evaluated: the loss factors (mechanical, δm , and electrical, δe ), the effective thickness coupling factor kt, the longitudinal wave velocity cL and the dielectric constant at constant strain.

III. COMPLETE SET OF CO DOPED BT

For modeling ultrasonic transducers with numerical methods such as Finite Element Method, all components from elastic, dielectric and piezoelectric tensors are required. Depending on the calculation performed, only one subset of tensors of the four possible ones is chosen (for example, s E , d and ε T ). Regardless of this choice, the calculation has to give a unique and similar result. If the result depends on the choice of the subset, the complete set is inconsistent [START_REF] Yu | Inconsistencies of the complete sets of electromechanical constants of relaxor-ferroelectric single crystals[END_REF][START_REF] Rouffaud | Complete electroelastic set for the (YXt)-45° cut of a KNbO3 single crystal[END_REF]. The complete set provided in this paper avoids this important defect.

A. Consistency and quantification criterion

In order to provide a complete and consistent set, all of the 16 intercoefficient relations [START_REF] Ikeda | Fundamentals of Piezoelectricity[END_REF] that allow switching from a subset to another, must be respected. Because of the measured parameters, experimental inconsistencies can be introduced in the complete set due to the experimental method chosen. So, intercoefficient relations cannot be strictly respected but, in the method reported here, particular attention will be paid to reduce these errors. Among the 40 constants that make a 6mm structure database, 18 constants are measured or calculated in a first step (part II.C) and are indicated in Table I with stars. So, 22 values have to be determined to complete the material set. The genetic algorithm (GA) [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF] is used to examine the wide search space created by these 22 parameters. In order to limit the search space without damaging the quality of the solution, upper and lower boundaries are defined for the investigated values between 0 and 2 times the reference values. In this work, reference values are taken from the Bechmann's results [START_REF] Bechmann | Elastic, Piezoelectric, and Dielectric Constants of Polarized Barium Titanate Ceramics and Some Applications of the Piezoelectric Equations[END_REF] about the complete electroelastic set of barium titanate ceramic because of the similar compositions. The GA attempts to reduce the fitness value p as it is defined in [START_REF] Rouffaud | Complete electroelastic set for the (YXt)-45° cut of a KNbO3 single crystal[END_REF]. Relative distances are calculated between both right-and left-hand sides of each intercoefficient relation. The value p represents

B. Determination of the electroelastic moduli

The complete set of obtained material parameters is given in Table II (last page). As mentioned before, 18 parameters are measured and are marked with a star. The maximum p value is found at 7.4% for the following intercoefficient relation (only for k=3 and q=3):

g kq = h kp s pq D (1)
More precisely, the value of the left member (g 33 ) of the relation ( 1) is equal to 0.0123 V.mN -1 while the right member has the value 0.0130 V.mN -1 .

This p-value obtained is in agreement with those that can be found in the literature calculated from the available database. As example, one of the best complete electroelastic sets for a lead-free single crystal [11] exhibits a p-value of 8.5%.

IV. 1-3 PIEZOCOMPOSITE AND TRANSDUCERS

Table II summarizes the electromechanical properties of 1-3 piezocomposite for the thickness mode. The resonant frequency in air is at 6.1 MHz for a thickness around 400 μm. According to the design of the piezocomposite (pitch and kerf values), the thickness mode is not disturbed by the first lateral mode, as shown in Fig. 3. The corresponding coupling factor slightly increases compared to the ceramic (from 45% to 47%) but do not reach the k 33 value of the ceramic given in Table I. With BTCaCo disc (thickness of 588μm and area of 140 mm 2 ), a single-element transducer is fabricated. The backing is made with epoxy resin loaded with tungsten powder (10% volume fraction) to deliver an acoustical impedance of 5.2 MRa. No matching layer is added. Two electrical contacts are made on each electrode with thin copper wires and conductive epoxy resin. The transducer is inserted into an insulating housing (Fig. 4). A coaxial cable with the characteristic impedance of Zc = 50 Ω and a length of around 40 cm is used for the electrical connection to the pulser/receiver (PR5900, panametrics). The transducer is placed in water in front of a metallic target to measure the electroacoustic response. It is recorded with an oscilloscope (Lecroy Waverunner 64XI). Fig. 5 shows this experimental response which has a center frequency of 4.8 MHz. According to the low backing acoustical impedance used with no matching layer, fractional bandwidth is logically limited (20% at -6dB and 59% at -20 dB). With the KLM scheme [START_REF]ANSI/IEEE Standard on Piezoelectricity[END_REF], this response is modeled with good accuracy (Fig. 5) and same sensitivity. Finally, the BTCaCo disc is replaced by a standard PZT (Pz27, MEGGITT A/S [START_REF] Levassort | A matrix method for modeling electroelastic moduli of 0-3 piezo-composites[END_REF]) with the same area. Only the thickness is adapted to keep the same center frequency. This last theoretical electroacoustic response is superimposed with the two previous ones (Fig. 5). Differences of sensitivity and bandwidth (20.7% at -6dB) are low.

V. CONCLUSIONS

Co doped barium titanate (BTCaCo) was fabricated with an objective of a preindustrial process in particular with a low sintering temperature while keeping good electromechanical performance. Full consistent database (40 elastic, dielectric and piezoelectric parameters) of BTCaCo was obtained with, for example a thickness coupling factor (k t ) of 45%.

Piezocomposites with a ceramic volume fraction of 50% were successfully fabricated, thanks to a dense microstructure and low grain size. The BTCaCo-based transducer, with a low acoustical impedance backing (5.2 MRa), possess a fractional bandwidth (-6dB) at 20% and a resonant frequency at 4.8 MHz. These results were in good agreement with simulated transducer behaviors using KLM scheme. They were compared with a transducer based on a standard PZT pellet to deliver similar sensitivities and bandwidth showing the efficiency of co-doped barium titanate.

978-1-5386-3383-0/17/$31.00 ©2017 IEEE

Fig. 1 .

 1 Fig. 1. Schematic representations of the different geometries needed to access to a complete electroelastic moduli.

Fig. 2 .

 2 Fig. 2. Optical photographs of (a) fabricated 1-3 piezocomposite disc (top view) and (b) zoom of the structure (squares : ceramic rods surrounded by epoxy resin).

Fig. 3 .

 3 Fig. 3. Complex electrical impedance (Z) of the BTCaCo based 1-3 piezocomposite (gray solid lines : theoretical curves with KLM scheme, dashed black lines : experimental curves).

  5386-3383-0/17/$31.00 ©2017 IEEE the maximum values of these relative distances and gives an indication of the consistency of the delivered full database.

Fig. 4 .

 4 Fig. 4. Photography of the fabricated transducer including BTCaCo disc.

Fig. 5 .

 5 Fig. 5. Experimental and theoretical electroacoustic responses of BTCaCo based transducer (dashed black line: experimental; thick solid black line: theoretical) and Pz27 theoretical thin solid black line).

TABLE I

 I

		.	FUNCTIONAL PROPERTIES OF BTCACO-BASED 1-3
			PIEZOCOMPOSITE
	t (μm)	A (mm 2 )	ε ε33 S /ε0	kt (%)	δ δm (%)	δ δe (%)	Z (MRa)
	404	95	580	47	1.5	1	17.1

a. t: thickness; A: area of the sample; ε33 S /ε0: relative dielectric constant at constant strain; kt: effective thickness coupling factor; δm: mechanical losses; δe: dielectric losses; Z: acoustical impedance.

TABLE II .

 II COMPLETE SET OF BTCACO CERAMIC ACCORDING TO FUNDAMENTAL PIEZOELECTRIC RELATIONS FOR VARIOUS INDEPENDENT SETS. STARS

			INDICATE MEASURED VALUES
		( ) , S E -type fundamental relations : E ij c (10 9 Pa), ij e (Cm -2 ),	S ij ε ε	0
	11 E c *	12 E c	13 E c	33 E c *	44 E c	66 E c
	150	66.94	54.26	138.9	40.84	41.53
	31 e	33 e *	15 e		ε	11 S	33 S ε *
	-8.90	19.31	13.14		1215	1179
		( ) , S D -type fundamental relations : D ij c (10 9 Pa), ij h (10 9 Vm -1 ), S ij β (10 9 mF -1 )
	11 D c	12 D c	13 D c	33 D c *	44 D c	66 D c
	157.59	74.58	38.89	174.70	56.42	41.50
	31 h	33 h	15 h		β	11 S	β	33 S
	-0.85	1.73	1.21		0.093	0.089
		( ) , T E -type fundamental relations : E ij s (10 -12 m 2 N -1 ), ij d (10 -12 CN -1 ),	T ij ε ε	0
	11 E s *	12 E s *	13 E s *	33 E s *	44 E s	66 E s
	8.89	-3.15	-2.26	8.97	24.49	24.08
	31 d *	33 d *	15 d		ε	11 T	33 T ε *
	-92.3	205	321.73		1690	1824
		( ) , T D -type fundamental relations : D ij s (10 -12 m 2 N -1 ), ij g (VmN -1 ), T ij β (10 9 mF -1 )
	11 D s *	12 D s *	13 D s *	33 D s *	44 D s	66 D s
	8.35	-3.69	-1.09	6.47	17.72	24.08
	31 g *	33 g *	15 g		β	11 T	β	33 T
	-0.0058	0.0123	0.021		0.067	0.059
		Coupling factors k (%)					Density (kg/m 3 )
	t k *	p k *	31 k *	33 k *				ρ *
	45.2	42.9	24.4	53.9				5810

This work was supported by the French Research Agency (HEcATE, ANR 14-CE07-0028-01)


C Bantignies
email: c.bantignies@vermon.com

E Filoux

P Mauchamp

R Dufait

M Pham

R Rouffaud

J M Grégoire

F Levassort

Lead-free high-frequency linear-array transducer (30 MHz) for in vivo skin imaging

Keywords: lead-free, 1-3 piezocomposite, linear array, high frequency, in-vivo skin imaging I

This work presents the fabrication of a 30 MHz, linear-array transducer based on a KN, 1-3 piezocomposite. Performances of the transducer were characterized and compared to a PZT-based linear array with similar structure. The composites were designed to minimize lateral modes of vibration which can severely degrade imaging performances. Fabrication steps were optimized to achieve a 40 MHz resonant frequency in air with a composite thickness of 69 microns. The measured thickness coupling factor was around 50 %. A 128element, linear array was then fabricated with 100 m pitch and 1,5 mm elevation aperture. The structure of the transducer (backing, matching layers, and electric components) was optimized to deliver good fractional bandwidth and sensitivity. The final probe was integrated in a prototype, real-time, 128channel scanner to acquire high-resolution images of the human skin in vivo. Results showed that, compared to PZT ceramics, KN single crystals provide low density and high acoustic velocity, both highly desirable for the manufacturing of HF transducers. The central frequency of the linear-array transducer was 30 MHz despite the KN composite being 20% thicker than equivalent PZT-based composites and the relative bandwidth was about 50%. High-resolution images of the human skin were acquired. A large ultrasound penetration due to good signal sensitivity was obtained and detailed features could be visualized.

INTRODUCTION

Over the last decades, lead zirconate titanate (Pb(Zr, Ti)O 3 , PZT) based ceramics and single crystals (PMN-PT) have been extensively used to manufacture high-frequency (HF, > 20MHz) ultrasound transducers for medical imaging as they have high piezoelectric coefficients and high dielectric permittivity. Due to health and environmental concerns, the European legislation recently adopted the RoHS (Restriction of the use of certain Hazardous Substances in electrical and electronic equipment) directive that requires lead-based components to be substituted by safer materials in most electrical devices. Even though medical devices are exempted from this legislation for the moment, it is necessary to look for new, safer piezoelectric components. For the past ten years, a lot of effort has been put into the development of lead-free piezoelectric materials that provide high performances and are suitable for industrial manufacturing [START_REF] Rödel | Perspective on the development of lead-free piezoceramics[END_REF]. Shrout et al. [START_REF] Shrout | Lead-free piezoelectric ceramics: alternatives for PZT?[END_REF] classified lead-free materials in three groups: BaTiO3, (Na,Bi)Ti-BaTiO3-(K,Bi)TiO3 (BNBK), and (K,Na)NbO3 (KNN) and several lead-free materials (ceramic and single crystal) were already tested for high-frequency transducer applications [START_REF] Chen | High-frequency ultrasonic transducer fabricated with lead-free piezoelectric single crystal[END_REF][START_REF] Lee | High-frequency ultrasonic transducer based on lead-free BSZT piezoceramics[END_REF][START_REF] Yan | Lead-free BNT Composite film for high-frequency broadband ultrasonic transducer applications[END_REF] and compared with lead-based transducers [START_REF] Jadidian | 25 MHz ultrasonic transducers with lead-free piezoceramic, 1-3 PZT fiver-epoxy composite, and PVDF polymer active elements[END_REF]. In 2012, Zhou et al. [START_REF] Zhou | Lead-free piezoelectric single crystal based 1-3 piezocomposites for ultrasonic transducer applications[END_REF] fabricated a 1-3 lead-free piezocomposite (NBT-BT/epoxy) for single transducer and linear array in the frequency range 3-5 MHz. Potassium niobate (KN) single crystal (KNN group) was identified as a good alternative to PZT because of its high thickness coupling coefficient [START_REF] Davis | Large and stable thickness coupling coefficients of [100]c-oriented KNbO3 and Li-modified (K,Na)NbO3 single crystal[END_REF] Nakamura et al. [START_REF] Nakamura | Domain structures in KNbO3 crystals and their piezoelectric properties[END_REF][START_REF] Nakamura | Orientation Dependence of Electromechanical Coupling Factors in KNbO3[END_REF] showed that a piezoelectric coupling coefficient of up to 69% can be achieved for specifics cuts and orientations of the KN crystal. Other interesting properties of KN are its relatively low mechanical impedance compared to PZT, for better impedance matching to tissues, and a high phase transition temperature (225°C) suitable for industrial processing and imaging applications. One drawback of single crystals is that they are usually more brittle than lead-based ceramics, which makes conventional processing such as dicing and lapping more difficult. In previous work [START_REF] Bantignies | Focused 20MHz single-crystal piezocomposite ultrasound array[END_REF] we demonstrated the feasibility of a 20 MHz, linear-array probe based on a single-crystal (PMN-PT), 1-3 piezocomposite. Here, a 30 MHz, 128-element probe was manufactured using a KN, 1-3 piezocomposite. The composite and the final probe were characterized, and in-vivo images of the human skin were acquired using a real-time, high-resolution imaging system developed at Tours University [START_REF] Bantignies | 40 MHz piezo-composite linear array for medical imaging and integration in a high resolution system[END_REF] .

II. ARRAY DESIGN AND CHARACTERIZATION

A. KN piezocomposite

Large polydomain plates, cut at 45° from an as-grown, [001]coriented KN crystal, were used to fabricate a 1-3 piezocomposite with fine structure. The composite was designed to minimize lateral modes of vibration which can severely degrade imaging performances. Previously, spatial homogeneity (of initial plate sample) of the electromechanical properties and as a function of the frequency was measured [START_REF] Rouffaud | Electro-elastic Moduli and Frequency Dependence of KN Single Crystal[END_REF] . Dicing conditions had to be finely tuned because the KN material has a higher Young's modulus than conventional lead-based single crystals, and it breaks more easily. However, the superior hardness of KN is an advantage when it comes to high frequencies due to the higher sound velocity in the crystal, allowing for high resonant frequencies to be achieved with relatively thick KN plates. The 1-3 piezocomposite had a kerf of about 10µm, a volume fraction of 64%, and a 69 µm thickness (Fig. 1).

Fig. 1 : Photograph of the KN single crystal 1-3 piezocomposite

The piezocomposite plate was machined to fabricate a linear array made of 128 elements. The complex electrical impedance of each element was measured using a HP4395 spectrum analyzer with its impedance test kit and specific spring-clip fixture. An equivalent electrical circuit scheme was used (KLM model [START_REF] Krimholtz | New equivalent circuit for elementary piezoelectric transducers[END_REF] ) to model the electrical behavior of the piezocomposite elements around their thickness-mode resonant frequency, and deduce their electrical, acoustic and electromechanical properties. Fig. 2 : shows the experimental impedance curve and the numerical fit used to determine the effective thickness coupling factor, the dielectric constant (at constant strain), the longitudinal wave velocity, and the electrical and mechanical losses. The resonant frequency in air was around 36MHz and the thickness coupling factor was 0.47. A low relative dielectric constant of 14 was obtained, which is likely due to a degradation of the crystal integrity during the dicing and lapping processes. B. 30 MHz, 128-element, lead-free probe A 128-element probe was manufactured using the KN-based 1-3 piezocomposite and the same processing steps than described in [START_REF] Michau | Piezocomposite "Piezocomposite 30MHz linear array for medical imaging: design challenges and performances evaluation of a 128 elements array[END_REF] The tip of the probe was sealed by an acoustic lens providing an elevation focusing of 8 mm.

Specifications of the transducer:

Center frequency: 30MHz Type of array: linear Number of elements: 128 Pitch: 100µm Elevation: 1.5 mm This acoustic head was integrated into housing and interconnected to a coaxial cable. This cable had 80 Ohms impedance and its length was 2 m (Fig. 3). A connector was mounted on the cable for interconnecting with the imaging system [START_REF] Bantignies | 40 MHz piezo-composite linear array for medical imaging and integration in a high resolution system[END_REF] The 128 elements were functional.

Fig. 3 : Photograph of the complete 128 elements array

Electro-acoustic characterization

The electro-acoustical properties of the complete array were characterized. A Panametrics 5073R pulser/receiver was used to excite each element using 50 settings with 20dB total gain on receive. Pulse-echo measurements were performed in water using a stainless steel, planar reflector positioned 4 mm away from the probe. The measured pulse-echo and corresponding spectrum of one element is shown in Fig. 4. Fréquence (MHz)

Im(Z) (Ω)
Frequency (MHz) Frequency (MHz)

µ µ µ µ F c =27 MHz
LCF@-6dB=20.3 MHz BW@-6dB=48.5% HCF@-6dB=32.3MHz Pulse length@-20dB = 143 ns Fc: center frequency; BW@-6dB: fractional bandwidth at -6dB; LCF:lower frequency; HCF: higher frequency at -6dB

Considering those measurements, it was estimated that the signal sensitivity was more than sufficient at 30 MHz to obtain a good penetration on in-vivo images.

2. Impedance measurement Electrical impedance of the final transducer without the cable was measured using an Agilent E5100A impedance analyzer. Fig. 5 shows the real and imaginary parts of the impedance for element 64. Measured values at 30MHz (R+jX with R=144 Ohms and X=-460 Ohms) were less favorable than those of a PZT-based probe previously reported in [START_REF] Michau | Piezocomposite "Piezocomposite 30MHz linear array for medical imaging: design challenges and performances evaluation of a 128 elements array[END_REF]. This result is a direct consequence of the low permittivity of the KN, 1-3 piezocomposite.

Fig. 5 : Real and imaginary parts of electrical impedance

for element 64

3. Acoustical field measurement A measurement of the acoustic field was performed in water using a 240 µm-diameter steel wire mounted on 2-axes, motorized stage. A representative element was excited using a 5900 Panametrics pulser-receiver with 40dB gain in receive. Results are presented in Fig. 6 and confirm that good sensitivity was achieved over a large depth-of-field. A custom, multi-channel, real-time, imaging system (ECODERM [START_REF] Bantignies | 40 MHz piezo-composite linear array for medical imaging and integration in a high resolution system[END_REF]) was used to perform in-vivo, highfrequency imaging of the human skin. Linear scanning was performed using 15 elements in transmit and 32 in receive. A fixed focus at 8 mm was used in transmit. A frame-rate of 40 Hz was achieved with this configuration. A 20 MHz PZTbased probe was also used for comparison. The 20 MHz probe had identical element count, pitch and kerf, and the same beamforming configuration was used. Fig. 7 shows images of a human forearm skin performed successively with the two probes. First, one can see that image resolution of the PZT probe was lower than that of the leadfree probe, which is mainly due to its lower central frequency.

The higher resolution of the 30 MHz probe allowed for clearly distinguishing the inner wall of the artery. A large depth-of field of over 6 mm could be obtained with the 30 MHz probe.

Fig. 2 :

 2 Fig. 2 : complex electrical impedance of 1-3 piezocomposite measured in air (blue line : experimental; redline : theoretical)

Fig. 4 :

 4 Fig. 4 : Typical measured time and frequency responses of one element of the final transducer.The characteristics deduced from the electro-acoustic response of one representative element (element 64) are summarized below:

Fig. 6 :

 6 Fig. 6 : Radiation pattern with an excited representative element

Fig. 7 :Fig. 8 :

 78 Fig. 7 : Ultrasound image of a human forearm skin performed with (a) 30 MHz lead-free probe and (b) 20 MHz PZT probe

ACKNOWLEDGMENT

Authors thank Leong Ratsimandresy (for acoustical measurements), Frédéric Ossant and Tony Mateo (for skin images). This work was funded by the French Research Agency (ANR HYPERCAMPUS Mat&Pro) and European Union (FEDER Funds).


Franc ¸ois Orieux

Raphael Chinchilla
email: raphael@chinchillas.net

Semi-Unsupervised Bayesian Convex Image Restoration with Location Mixture of Gaussian

Convex image restoration is a major field in inverse problems. The problem is often addressed by hand-tuning hyperparameters. We propose an incremental contribution about a Bayesian approach where a convex field is constructed via Location Mixture of Gaussian and the estimator computed with a fast MCMC algorithm. Main contributions are a new field with several operator avoiding crosslike artifacts and a fallback sampling algorithm to prevent numerical errors. Results, in comparison to standard supervised results, have equivalent quality in a quasi-unsupervised approach and go with uncertainty quantification.

I. INTRODUCTION

Image restoration, or deconvolution, is a major problem with an abundant literature and applications found in, for instance, optics, radio astronomy, microscopy [START_REF] Chakrova | Deconvolution methods for structured illumination microscopy[END_REF]- [START_REF] Giovannelli | Positive deconvolution for superimposed extended source and point sources[END_REF]. Deconvolution is also related to the resolution of ill-posed inverse problems where the likelihood of data presents defects like missing data or instability [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF].

Optimization methods, based on the minimization of mixed criteria, are common [START_REF] Komodakis | Playing with Duality: An Overview of Recent Primal-Dual Approaches for Solving Large-Scale Optimization Problems[END_REF], [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] and popular in high-dimension, thanks to efficient algorithms, especially in the convex case [10]. However, the hyper-parameters that determine the balance between the data and prior adequacy are usually handtuned (supervised) and uncertainty about the optimum estimate are not well defined or studied.

Bayesian methods, on the contrary, present a natural way to estimate the hyper-parameters (unsupervised) by considering them as nuisance parameter and marginalization of the joint a posteriori law [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF], [START_REF] Scarpa | Hierarchical Multiple Markov Chain Model for Unsupervised Texture Segmentation[END_REF]. Moreover, a posteriori law analysis brings confidence interval analysis and quantification of uncertainty about the estimation.

Within the Bayesian methods, variationnal approaches had interest by providing apparently fast algorithms [START_REF] Zheng | Efficient Variational Bayesian Approximation Method Based on Subspace optimization[END_REF]- [START_REF] Babacan | Variational Bayesian Super Resolution[END_REF]. However, variationnal approaches build an approximation of the a posteriori law and the uncertainty is known to be underestimated, induced by the separability and the noncorrelation model.

The approach of this work is the use of Monte Carlo Markov Chains (MCMC) algorithms for exact a posteriori law exploration and computation of the Expectation a posteriori estimator (EAP). However, results on unsupervised method are sparser with some existing work [START_REF] Pereyra | A Survey of Stochastic Simulation and Optimization Methods in Signal Processing[END_REF], for instance based on marginal likelihood [START_REF] Descombes | Estimation of Markov Random Field prior parameters using Markov Chain Monte Carlo Maximum Likelihood[END_REF], quadratic prior [START_REF] Gilavert | Efficient Gaussian Sampling for Solving Large-Scale Inverse Problems Using MCMC[END_REF], [START_REF] Orieux | Sampling high-dimensional Gaussian distributions for general linear inverse problems[END_REF], per variables Metropolis-Hastings [START_REF] Chaari | A Hierarchical Bayesian Model for Frame Representation[END_REF] or more recently Moreau approximation [START_REF] Pereyra | Maximum-a-posteriori estimation with Bayesian confidence regions[END_REF].

In [START_REF] Orieux | Bayesian estimation of regularization and psf parameters for wiener-hunt deconvolution[END_REF], authors present a fast MCMC deconvolution method limited to quadratic prior. Work of J.-F. Giovannelli [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF] presents an unsupervised convex deconvolution approach. However, the use of multiple regularization operator, like horizontal and vertical gradient, is missing. In addition, some numerical instabilities are present due to special functions. Finally, divergence of the MCMC chains may be observed in practice when the full set of hyper-parameters is estimated.

We propose an incremental work based on the method of [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF]. We show that multiple operators can be used, providing better results and no crosslike effect of the Laplacian around the edges [START_REF] Chinchilla | Efficiently Sampling from High-Dimensional Mixture of Gaussians Distributions and Applications to Inverse Problems[END_REF]. To prevent numerical instability of the simulation we present a fallback solution based on a simple and efficient Metropolis-Hastings algorithm. Finally, by just considering a known or estimated noise level, the method becomes robust and estimate the a priori shape and precision. The proposed algorithm does not pretend to achieve the best possible results but rather to complement ongoing work for unsupervised convex image restoration.

The paper is organized as follows : Sec. II presents the notation and the methodological context. Sec. III is devoted to the convex potential, the unsupervised approach and the MCMC algorithm. Finally, Sec. IV shows some results and comparison to standard algorithm of the same class.

II. NOTATIONS AND PROBLEMS

We consider N pixels images x ∈ R N . The data model is y = Hx + n where y ∈ R N collects the data, H ∈ R N ×N is the convolution operator and n ∈ R N an unknown random noise. The convolution are circulant and computed via Discrete Fourier Transform F with a diagonalizable operator H = F † Λ H F , where Λ H = diag( h0 , . . . , hN-1 ).

The a priori noise law is supposed white Gaussian with a known noise level

σ n = γ -1/2 n > 0 leading to the likelihood p(y | x; γ n ) = (2π) -N 2 γ N 2 n exp - γ n 2 y -Hx 2 . (1)
We also consider a priori image law using Markov Random Field

p(x | θ) = K(θ) -1 exp - 1 2 E θ (x) (2)
where K(θ) = exp(-E θ (x)/2) dx is the partition function.

In the general case, K is dependant wrt. regularization parameters θ which is a major difficulty that blocks θ estimation [START_REF] Descombes | Estimation of Markov Random Field prior parameters using Markov Chain Monte Carlo Maximum Likelihood[END_REF], [START_REF] Orieux | Bayesian estimation of regularization and psf parameters for wiener-hunt deconvolution[END_REF], [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF].

For the well-posedness of the inverse problem we consider the energy E θ (x) = γ x c∈C φ (d t c x; θ) with the set C of cliques c and neighborhoods d c [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF] and γ x > 0. For edge preserving consideration and uniqueness of the supervised solution, the potential φ is chosen as convex.

A common estimator in image processing community is the maximum a posteriori (MAP) obtained by minimization of the anti logarithm of the a posteriori law

x MAP = arg min x -ln(p(x | y, γ n , θ)) (3)
= arg min

x y -Hx 2 + γ x γ n E θ (x). (4)
In that case of convex optimization, an abundant literature is available to compute it efficiently [10], [START_REF] Idier | Convex Half-Quadratic Criteria and Interacting Auxiliary Variables for Image Restoration[END_REF].

III. BAYESIAN FRAMEWORK AND ALGORITHM

A. LogErf potential and LMG

The LogErf potential φ [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF], [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF] is the convex function

φ(u) = -2 ln(χ(u) + χ(-u)) (5)
with

χ(u) = exp γ b 2 u erfc γ x 2 γ b 2γ x + u (6)
with γ x > 0 and γ b > 0. This potential, liken to Huber potential, brings several advantages when used for Gibbs field: 1) φ is convex and, therefore, E θ (x) is convex wrt. x.

2) The distribution p(u) can be expressed as a Location Mixture of Gaussian (LMG)

p(u) ∝ exp - 1 2 φ(u) ∝ R p(u | b) p(b) db (7) = P -1 R exp - γ x 2 (u -b) 2 - γ b 2 |b| db (8)
with a Laplace distribution p(b) on the mean of p(u | b).

3) The partition function is tractable and does not depend on b. For the LogErf potential, we have

P = √ 2πγ -1 2 x γ b 4 -1 . (9)
A Gibbs field on the image x can be constructed as

p(x | γ x , γ b ) = K -1 exp - γ x 2 Dx -b 2 c∈C exp γ b 2 |b c | db (10)
that is, an LMG for each clique d t c x associated with an auxiliary variable b c . Initial work [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF] states the necessity that the number of clique equals the number of pixels. The operator was circulant convolution with second order difference (Laplacian), leading to a priori independent clique in image and Fourier space

p(x) ∝ N c=1 R exp -γ x |x c -b c | 2 exp γ b 2 |b c | db c (11)
where x = Dx, and xc the c-element of x. However, this model introduces crosslike artifacts around the edges [START_REF] Chinchilla | Efficiently Sampling from High-Dimensional Mixture of Gaussians Distributions and Applications to Inverse Problems[END_REF]. Moreover, this limitation is not necessary and a field build with more clique and auxiliary variables than pixels is possible, at the price of no longer orthogonal cliques.

For image, we consider two group of cliques (horizontal and vertical gradient for instance) D h and

D v with conditional mean b h ∈ R N and b v ∈ R N p(x | b) ∝ exp - γ x 2 Dx -b 2 (12) = (2π) -N 2 |Q| 1 2 exp - (x -µ) t Q(x -µ) 2 (13)
with [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF], the distribution on x is not a product of independent distributions on each clique, or two a priori law on each clique group [START_REF] Chinchilla | Efficiently Sampling from High-Dimensional Mixture of Gaussians Distributions and Applications to Inverse Problems[END_REF]. The prior mean µ is therefore a mix of the auxiliary variables. Eq. ( 13) is sufficient statistics form but Eq. ( 12) is more natural where each clique have an unknown mean and, consequently, each clique has a LMG. The full marginal a priori law for the image writes

D t = [D t h , D t v ], b t = [b t h , b t v ], Q = γ x (D t h D h + D t v D v ) and µ = Q -1 D t b, if Q is not singular. Contrary to
p(x) = p(x | b) 2N c=1 p(b c ) db, (14)
remains an LMG and, by choosing a Laplace distribution for b c , the partition function is expressed wrt. hyper-parameters as

K(γ x , γ b ) ∝ γ -N 2 x γ -2N b . (15)
Finally, the existence condition is not that the number of clique equals the number of pixels but that Q is not singular. Since the gradient operator leads to singular Q and improper prior with a null eigen-value for the mean level, we refer to papers [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF] and [START_REF] Orieux | Bayesian estimation of regularization and psf parameters for wiener-hunt deconvolution[END_REF] and the change to

K(γ x , .) ∝ γ -(N -1)/2 x .

B. Posterior Law and Estimators

A full joint a posteriori law of the unknown can now be expressed, with the horizontal and vertical gradient,

p(x, b, γ x , γ b | y; γ n ) ∝ p(y | x; γ n ) p(x | b, γ x ) p(b | γ b ) p(γ x ) p(γ b ). (16)
The a priori laws is correlated Gaussian for x, Laplacian for b and conjugate for γ x and γ b , that is the non-informative Jeffrey's distribution [START_REF] West | Outlier Models and Prior Distributions in Bayesian Linear Regression[END_REF]. The a posteriori law writes

p(x, b, γ x , γ b | y; γ n ) ∝ γ - (N -1) 2 x γ -2N b exp - γ n 2 y -Hx 2 exp - γ x 2 Dx -b 2 exp - γ b 2 b 1 . ( 17
)
Several estimators are studied:

1) The Supervised Expectation a posteriori (SEAP)

x SEAP = x p(x | γ x , γ b , y) dx (18) = x p(x, b | γ x , γ b , y) dx db. ( 19
)
This estimator is the main reference for the proposed unsupervised EAP estimator. It allows finding the best hyper parameter values γ x and γ b given the true signal x and a measure. Then the image

x SEAP = x p(x | γ x , γ b , y) dx (20)
is considered as a reference, being the best possible reconstructed image given that model.

2) The supervised Maximum a posteriori (MAP)

x MAP = arg min x -ln(p(x | γ x , γ b , y)) (21)
= arg min

x -ln p(x, b | γ x , γ b , y) db . ( 22
)
Champagnat et. al. shown in [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF] that this estimator can be computed by an EM algorithm equivalent to the efficient HQ Geman & Yang optimization scheme. The MAP also allow to qualify and quantify the differences, if any, with the SEAP, and can also be used to determine best γ x and γ b .

3) The semi-unsupervised Expectation a posteriori (EAP)

x EAP = x p(x | y) dx (23)
= x p(x, b, γ x , γ b | y) dx db dγ x dγ b . ( 24
)
The EAP estimator is the objective and will be compared to supervised reconstructions. The algorithm, presented Sec. III-C, also provides the EAP γx and γb that can be compared to γ x and γ b . 4) The unsupervised Expectation a posteriori with quadratic prior, that is with prior model

p(x | γ x ) ∝ exp -γ x Dx 2 /2 . (25)
This estimator x QUAD is known to produce near optimal hyper-parameter estimation when a quadratic penalization is used [START_REF] Orieux | Bayesian estimation of regularization and psf parameters for wiener-hunt deconvolution[END_REF].

C. MCMC Algorithm for EAP and SEAP

The computation of the EAP estimator is based on a Gibbs sampler that successively simulate, after an initialization k = 0, b (0) = 0 and γ (0)

x = γ (0) b = 1, the conditional a posteriori laws as 1) x (k+1) ∼ p x | b (k) , γ (k) x , y , 2) b (k+1) ∼ p b | x (k+1) , γ (k)
x , γ

(k) b , 3) γ (k+1) x ∼ p γ x | x (k+1) , b (k+1) , 4) γ (k+1) b ∼ p γ b | b (k+1) , 5) k ← k + 1.
For the SEAP estimator, step 3 and 4 are not undertaken as γ x , γ b are fixed, eventually to γ x and γ b .

The conditional a posteriori law of x is Gaussian. Since all operators are circulant convolution, this law can be simulated very efficiently with diagonal matrix in Fourier space

p x | b, γ x , y ∝ exp - γ n 2 ẙ -Λ H x 2 - γ x 2 Λ h x -bh 2 - γ x 2 Λ v x -bv 2 . ( 26
)
A sample is x (k+1) = F † x(k+1) where

x(k+1) = γ n ΣΛ † H F y + Σ 1 2 F (27) with Σ -1 = γ n |Λ H | 2 +γ x (|Λ h | 2 +|Λ v | 2 )
, and ∼ N (0, I).

The conditional a posteriori law of auxiliary variables b is more delicate but is a posteriori conditionally independent

p(b | x, γ x , γ b ) ∝ 2N c=1 exp - γ x 2 |x c -b c | 2 -γ b |b c | (28)
where x = Dx. The choice of the LogErf allows fast simulation by inversion of the cumulative density function (icdf) as described in [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF]. However, practical use of the icdf shows that numerical error and floating point arithmetic overflow may arise, even if small fraction of b is concerned. In that case we propose a fallback on a random-walk Metropolis-Hastings (MH) step. This fallback algorithm rarely occurs, is for independent scalar variables, and is used only for the elements of b where the icdf has failed. In practice, the MH appears not to be difficult to tune, is efficient, and produce satisfactory results.

Lastly, the conditional a posteriori laws for hyperparameters are Gamma distribution

p(γ x | x, b) ∝ γ - (N -1) 2 x exp -γ x Dx -b 2 2 (29) and p(γ b | b) ∝ γ -2N b exp -γ b 2N c=1 |b c | 2 . ( 30
)
These scalar distribution are simulated with standard toolbox available in Matlab or Numpy for instance.

The SEAP and EAP are approximated by the empirical mean of K samples, after discarding samples of the burn-in period,

x EAP ≈ 1 K K k=1 x (k) , b EAP ≈ 1 K K k=1 b (k) . (31)
The pixel variance, diagonal of the posterior covariance, is also approximated by

σ 2 EAP = diag(Σ x,EAP ) ≈ 1 K -1 K k=1 x (k) -x EAP 2 . (32) IV. RESULTS
The proposed method is tested on the "cameraman" image, a N = 256×256 pixels image with strong discontinuities. The observation operator has a squared impulsionnal response (IR) of size 5 × 5 pixels. The convolution is done by filtering in frequencies space as y = F † Λ H F x +n, where x stands for the true image. The tested noise inverse variance γ n are 100, 1 and 0.1. The horizontal and vertical gradients are computed via the frequency space

x h = F † Λ h F x and x v = F † Λ v F x.
The MCMC algorithm produces 1000 samples [x, b, γ x , γ b ], with a burn-in period of 500 samples. These numbers are chosen to have visually a sufficiently good exploration after the burn-in period. The algorithm, implemented with Python and Scipy has been run on a 1.9 GHz processor. The total time to produce 1000 samples is around 80 seconds, with 0.08 seconds per iteration. The majority of time (60%) is spent inside the icdf simulation and special functions evaluation.

Fig. 1 shows the true image x and the data y with γ n = 1. The best Supervised Expectation a posteriori x SEAP is illustrated in Fig. 1d, where the hyper-parameters γ x , γ b are determined by minimization of the 1 error x SEAP (γ x , γ b ) -x 1 . The minimization is done by exhaustive search to reach the near optimal SEAP solution. The Fig. 1c shows our proposition where hyper-parameters are automatically estimated from the data without hand-tuning and without knowledge of the true image, except the noise level.

For comparison, the unsupervised quadratic solution, x QUAD Fig. 1e, has some noise residual, well visible inside flat region, and Gibbs effect is present near the edges. These defects are no more visible on the proposed solution x EAP , see Fig. 1c, as for the other convex solutions. To assess the good results of the semi-unsupervised EAP estimator, Fig. 1f shows the best supervised convex x MAP estimator, known to provide good results [START_REF] Idier | Convex Half-Quadratic Criteria and Interacting Auxiliary Variables for Image Restoration[END_REF]. In that case, the hyper-parameters are also found via the minimization of the 1 error x MAP (γ x , γ b ) -x 1 . Almost no differences is visible by eyes between supervised x SEAP , x MAP and our unsupervised proposition x EAP .

Thanks to the MCMC algorithm, the EAP estimator also provides the uncertainties for every quantities, notably the image x illustrated Fig. 2a, where uncertainty is concentrated near region of strong gradient. Chains of hyper-parameters are illustrated in Fig. 3. The γ b chains converge in short time and present a small dispersion in regards to the γ x chain. An explanation is the small variation of the auxiliary variables during the iteration. On the contrary, the γ x chain presents a longer convergence time with more intra-correlation. The chain still converges after approximately 200 iterations. A possible explanation is the greater sensibility of prior adequacy wrt. the auxiliary variables changes but further investigations are needed. Finally, tests with more and less noise is also presented Fig. 4. V. CONCLUSION This paper presents an incremental contribution to the unsupervised convex image restoration within a Bayesian approach based on [START_REF] Giovannelli | Unsupervised Bayesian Convex Deconvolution Based on a Field With an Explicit Partition Function[END_REF]. However, the proposed approach allows to use several regularization operator for better image restoration avoiding crosslike effect. In addition, we propose a fallback algorithm based on a random walk Metropolis-Hastings step to avoid the numerical instability of the original algorithm. The noise level is supposed known and it's estimation jointly with the other is the main perspective.

  Fig 2b shows the Expectation a posteriori of b v . The values are naturally around zero and variables are clearly able to detect edge as awaited.

Fig. 1 .

 1 Fig. 1. Results with noise precision γn = 1. The true image x and data y are represented Fig. 1a and 1b. The Fig. 1d is the best SEAP solution and Fig. 1c the proposed unsupervised restoration. For comparison Fig. 1e shows the quadratic unsupervised solution and 1f the best supervised MAP.

Fig. 2 . 2 ( 2 ( 1 Fig. 3 .

 22213 Fig.2. Results with noise precision γn = 1. Fig.2ashows the estimated standard deviation of each pixel. Fig.2bshows the estimated vertical auxiliary variables. Both image are slightly saturated for display purpose.

Fig. 4 .

 4 Fig. 4. Results with different noise precision γn = 0.1 and 100. Quality and absence of Gibbs effect on the camera arm is notable.


Ioane Muni Toke
email: ioane.muni-toke@ecp.fr

"Market making" in an order book model and its impact on the spread

It has been suggested that marked point processes might be good candidates for the modelling of financial high-frequency data. A special class of point processes, Hawkes processes, has been the subject of various investigations in the financial community. In this paper, we propose to enhance a basic zero-intelligence order book simulator with arrival times of limit and market orders following mutually (asymmetrically) exciting Hawkes processes. Modelling is based on empirical observations on time intervals between orders that we verify on several markets (equity, bond futures, index futures). We show that this simple feature enables a much more realistic treatment of the bid-ask spread of the simulated order book.

Introduction

Arrival times of orders: event time is not enough As of today, the study of arrival times of orders in an order book has not been a primary focus in order book modelling. Many toy models leave this dimension aside when trying to understand the complex dynamics of an order book. In most order driven market models such as [START_REF] Cont | Herd behavior and aggregate fluctuations in financial markets[END_REF][START_REF] Lux | Volatility clustering in financial markets[END_REF][START_REF] Alfi | Minimal agent based model for financial markets i -origin and self-organization of stylized facts[END_REF], and in some order book models as well (e.g. [START_REF] Preis | Multi-agent-based order book model of financial markets[END_REF]), a time step in the model is an arbitrary unit of time during which many events may happen. We may call that clock aggregated time. In most order book models such as [START_REF] Challet | Analyzing and modeling 1+ 1d markets[END_REF][START_REF] Mike | An empirical behavioral model of liquidity and volatility[END_REF][START_REF] Cont | A Stochastic Model for Order Book Dynamics[END_REF], one order is simulated per time step with given probabilities, i.e. these models use the clock known as event time. In the simple case where these probabilities are constant and independent of the state of the model, such a time treatment is equivalent to the assumption that order flows are homogeneous Poisson processes. A probable reason for the use of non-physical time in order book modellingleaving aside the fact that models can be sufficiently complicated without adding another dimension -is that many puzzling empirical observations (see e.g. [START_REF] Chakraborti | Econophysics: Empirical facts and agent-based models[END_REF] for a review of some of the well-known "stylized facts") can (also) be done in event time (e.g. autocorrelation of the signs of limit and market orders) or in aggregated time (e.g. volatility clustering). However, it is clear that physical (calendar) time has to be taken into account for the modelling of a realistic order book model. For example, market activity varies widely, and intraday seasonality is often observed as a well known U-shaped pattern. Even for a short time scale model -a few minutes, a few hours -durations of orders (i.e. time intervals between orders) are very broadly distributed. Hence, the Poisson assumption and its exponential distribution of arrival times has to be discarded, and models must take into account the way these irregular flows of orders affect the empirical properties studied on order books.

Let us give one illustration. On figure 1, we plot examples of the empirical density function of the observed spread in event time (i.e. spread is measured each time an event happens in the order book), and in physical (calendar) time (i.e. measures are weighted by the time interval during which the order book is idle). It appears that density of the most probable values of the time- weighted distribution is higher than in the event time case. Symmetrically, the density of the least probable event is even smaller when physical time is taken into account. This tells us a few things about the dynamics of the order book, which could be summarized as follows: the wider the spread, the faster its tightening. We can get another insight of this empirical property by measuring on our data the average waiting time before the next event, conditionally on the spread size. When computed on the lower one-thirdquantile (small spread), the average waiting time is 320 milliseconds. When computed on the upper one-third-quantile (large spread), this average waiting time is 200 milliseconds. These observations complement some of the ones that can be found in the early paper [START_REF] Biais | An empirical analysis of the limit order book and the order flow in the paris bourse[END_REF].

Counting processes with dynamic intensity

There is a trend in the econometrics literature advocating for the use of (marked) point processes for the modelling of financial time series. One may find a possible source of this interest in [START_REF] Engle | Forecasting the frequency of changes in quoted foreign exchange prices with the autoregressive conditional duration model[END_REF][START_REF] Robert | The econometrics of Ultra-High-Frequency data[END_REF], which introduce autoregressive conditional duration (ACD) and autoregressive conditional intensity (ACI) models. [START_REF] Hall | Modelling the buy and sell intensity in a limit order book market[END_REF] fit that type of models on the arrival times of limit, market and cancellation orders in an Australian stock market order book.

A particular class of point processes, known as the Hawkes processes, is of special interest for us, because of its simplicity of parametrization. A univariate linear self-exciting Hawkes process (N t ) t>0 , as introduced by [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF][START_REF] Hawkes | A cluster process representation of a Self-Exciting process[END_REF], is a point process with intensity:

λ(t) = λ 0 + t 0 ν(t -s)dN s , (1)
where the kernel ν is usually parametrized as ν(t) = αe -βt or in a more general way ν(t) = p k=1 α k t k e -βt . Statistics of this process are fairly wellknown and results for a maximum likelihood estimation can be found in [START_REF] Ogata | On Lewis' simulation method for point processes[END_REF]. In a multivariate setting, mutual excitation is introduced. A bivariate model can thus be written:

λ 1 (t) = λ 1 0 + t 0 ν 11 (t -s)dN 1 s + t 0 ν 12 (t -s)dN 2 s λ 2 (t) = λ 2 0 + t 0 ν 21 (t -s)dN 1 s + t 0 ν 22 (t -s)dN 2 s ( 2
)
The use of these processes in financial modelling is growing. We refer the reader to [START_REF] Bauwens | Modelling financial high frequency data using point processes[END_REF] for a review and [START_REF] Hautsch | Modelling irregularly spaced financial data[END_REF] for a textbook treatment. In [START_REF] Bowsher | Modelling security market events in continuous time: Intensity based, multivariate point process models[END_REF], a bivariate (generalized) Hawkes process is fitted to the time series of trades and mid-quotes events, using trading data of the General Motors stock traded on the New York stock Exchange. In [START_REF] Large | Measuring the resiliency of an electronic limit order book[END_REF] a ten-variate Hawkes process is fitted to the Barclay's order book on the London Stock Exchange, sorting orders according to their type and aggressiveness. It is found that the largest measured effect is the exciting effect of market orders on markets orders. [START_REF] Hewlett | Clustering of order arrivals, price impact and trade path optimisation[END_REF] fits a bivariate Hawkes model to the time series of buy and sell trades on the EUR/PLN (Euro/Polish Zlotych) FX market. Using the simplest parametrization of Hawkes processes and some (very) constraining assumptions, some analytical results of trade impact may be derived. [START_REF] Bacry | Modeling microstructure noise using point processes[END_REF] fits a bivariate Hawkes process to the trade time series of two different but highly correlated markets, the "Bund" and the "Bobl" (Eurex futures on mid-and long-term interest rates).

Organization of this paper

In this paper, we propose to enhance a basic order book simulator with arrival times of limit and market orders following mutually (asymmetrically) exciting Hawkes processes. Modelling is based on empirical observations, verified on several markets (equities, futures on index, futures on bonds), and detailed in section 1. More specifically, we observe evidence of some sort of "market making" in the studied order books: after a market order, a limit order is likely to be submitted more quickly than it would have been without the market order. In other words, there is a clear reaction that seems to happen: once liquidity has been removed from the order book, a limit order is triggered to replace it. We also show that the reciprocal effect is not observed on the studied markets. These features lead to the use of unsymmetrical Hawkes processes for the design of an agent-based order book simulator described in section 2. We show in section 3 that this simple feature enables a much more realistic treatment of the bid-ask spread of the simulated order book.

1 Empirical evidence of "market making" We use Reuters RDTH tick-by-tick data from September 10th, 2009 to September 30th, 2009 (i.e. 15 days of trading). For each trading day, we use only 4 hours of data, precisely from 9:30 am to 1:30 pm. As we are studying European markets, this time frame is convenient because it avoids the opening of American markets and the consequent increase of activity.

Our data is composed of snapshots of the first five limits of the order books (ten for the BNPP.PA stock). These snapshots are timestamped to the millisecond and taken at each change of any of the limits or at each transaction. The data analysis is performed as follows for a given snapshot:

1. if the transaction fields are not empty, then we record a market order, with given price and volume;

2. if the quantity offered at a given price has increased, then we record a limit order at that price, with a volume equal to the difference of the quantities observed; 3. if the quantity offered at a given price has decreased without any transaction being recorded, then we record a cancellation order at that price, with a volume equal to the difference of the quantities observed; 4. finally, if two orders of the same type are recorded at the same time stamp, we record only one order with a volume equal to the sum of the two measured volumes.

Therefore, market orders are well observed since transactions are explicitly recorded, but it is important to note that our measure of the limit orders and cancellation orders is not direct. In table 1, we give for each studied order book the number of market and limit orders detected on our 15 4-hour samples. On the studied period, market activity ranges from 2.7 trades per minute on the least liquid stock (LAGA.PA) to 14.2 trades per minute on the most traded asset (Footsie futures).

Empirical evidence of "market making"

Our idea for an enhanced model of order streams is based on the following observation: once a market order has been placed, the next limit order is likely to take place faster than usual. To illustrate this, we compute for all our studied assets:

• the empirical probability density function (pdf) of the time intervals of the counting process of all orders (limit orders and market orders mixed), i.e. the time step between any order book event (other than cancellation)

• and the empirical density function of the time intervals between a market order and the immediately following limit order.

If an independent Poisson assumption held, then these empirical distributions should be identical. However, we observe a very high peak for short time intervals in the second case. The first moment of these empirical distributions is significant: one the studied assets, we find that the average time between a market order and the following limit order is 1.3 (BNPP.PA) to 2.6 (LAGA.PA) times shorter than the average time between two random consecutive events. On the graphs shown in figure 2, we plot the full empirical distributions for four of the five studied assets 2 . We observe their broad distribution and the sharp peak for the shorter times: on the Footsie future market for example, 40% of the measured time steps between consecutive events are less that 50 milliseconds ; this figure jumps to nearly 70% when considering only market orders and their following limit orders. This observation is an evidence for Fig. 2. Empirical density function of the distribution of the time intervals between two consecutive orders (any type, market or limit) and empirical density function of the distribution of the time intervals between a market order and the immediately following limit order. X-axis is scaled in seconds. In insets, same data using a log-log scale. Studied assets: BNPP.PA (top left), LAGA.PA (top right), FEIZ9 (bottom left), FFIZ9 (bottom right).

some sort of market-making behaviour of some participants on those markets.

2 Observations are identical on all the studied assets.

It appears that the submission of market orders is monitored and triggers automatic limit orders that add volumes in the order book (and not far from the best quotes, since we only monitor the five best limits).

In order to confirm this finding, we perform non-parametric statistical test on the measured data. For all four studied markets, omnibus Kolmogorov-Smirnov and Cramer-von Mises tests performed on random samples establish that the considered distributions are statistically different. If assuming a common shape, a Wilcoxon-Mann-Withney U test clearly states that the distribution of time intervals between a market orders and the following limit order is clearly shifted to the left compared to the distributions of time intervals between any orders, i.e. the average "limit following market" reaction time is shorter than the average time interval between random consecutive orders.

Note that there is no link with the sign of the market order and the sign of the following limit order. For example for the BNP Paribas (resp. Peugeot and Lagardere) stock, they have the same sign in 48.8% (resp. 51.9% and 50.7%) of the observations. And more interestingly, the "limit following market" property holds regardless of the side on which the following limit order is submitted. On figure 3, we plot the empirical distributions of time intervals between a market order and the following limit order, conditionally on the side of the limit order: the same side as the market order or the opposite one. It appears for all studied assets that both distributions are roughly identical. In other words, we cannot distinguish on the data if liquidity is added where Fig. 3. Empirical density function of the distribution of the time intervals between a market order and the immediately following limit order, whether orders have been submitted on the same side and on opposite sides. X-axis is scaled in seconds. In insets, same data using a log-log scale. Studied assets: BNPP.PA (left), LAGA.PA (right).

the market order has been submitted or on the opposite side. Therefore, we do not infer any empirical property of placement: when a market order is submitted, the intensity of the limit order process increases on both sides of the order book. This effect we have thus identified is a phenomenon of liquidity replenishment of an order book after the execution of a trade. The fact that it is a bilateral effect makes its consequences similar to "market making", event though there is obviously no market maker involved on the studied markets.

1.3 A reciprocal "market following limit" effect ?

We now check if a similar or opposite distortion is to be found on market orders when they follow limit orders. To investigate this, we compute for all our studied assets the "reciprocal" measures:

• the empirical probability density function (pdf) of the time intervals of the counting process of all orders (limit orders and market orders mixed), i.e. the time step between any order book event (other than cancellation)

• and the empirical density function of the time step between a market order and the previous limit order.

As previously, if an independent Poisson assumption held, then these empirical distribution should be identical. Results for four assets are shown on figure 4. Contrary to previous case, no effect is very easily interpreted. For the three stocks (BNPP.PA, LAGA.PA and PEUP.PA (not shown)), it seems that the empirical distribution is less peaked for small time intervals, but difference is much less important than in the previous case. As for the FEI and FFI markets, the two distributions are even much closer. Non-parametric tests confirms these observations. Performed on data from the three equity markets, Kolmogorov tests indicate different distributions and Wilcoxon tests enforce the observation that time intervals between a limit order and a following market order are stochastically larger than time intervals between unidentified orders. As for the future markets on Footsie (FFI) and 3-month Euribor (FEI), Kolmogorov tests does not indicate differences in the two observed distributions, and the result is confirmed by a Wilcoxon test that concludes at the equality of the means.

Order book models with mutually exciting order flows

Following these previous observations, we enhance a basic agent-based order book simulator with dependence between the flows of limit and market orders.

The basic Poisson model

We use as base model a standard zero-intelligence agent-based market simulator built as follows. One agent is a liquidity provider. This agent submits limit orders in the order books, orders which he can cancel at any time. This agent is simply characterized by a few random distributions: Fig. 4. Empirical density function of the distribution of the time intervals between two consecutive orders (any type, market or limit) and empirical density function of the distribution of the time intervals between a limit order and an immediately following market order. In insets, same data using a log-log scale. Studied assets: BNPP.PA (top left), LAGA.PA (top right), FEIZ9 (bottom left), FFIZ9 (bottom right).

1. submission times of new limit orders are distributed according to a homogeneous Poisson process N L with intensity λ L ; 2. submission times of cancellation of orders are distributed according to homogeneous Poisson process N C with intensity λ C ; 3. placement of new limit orders is centred around the same side best quote and follows a Student's distribution with degrees of freedom ν P 1 , shift parameter m P 1 and scale parameter s P 1 ; 4. new limit orders' volume is randomly distributed according to an exponential law with mean m V 1 ; 5. in case of a cancellation, the agent deletes his own orders with probability δ.

The second agent in the basic model is a noise trader. This agent only submits market order (it is sometimes referred to as the liquidity taker). Characterization of this agent is even simpler: 6. submission times of new market orders are distributed according to a homogeneous Poisson process N M with intensity µ; 7. market orders' volume is randomly distributed according to an exponential law with mean m V 2 . For all the experiments, agents submit orders on the bid or the ask side with probability 0.5. This basic model will be henceforth referred to as "HP" (Homogeneous Poisson).

Assumptions 1, 2 and 6 (Poisson) will be replaced in our enhanced model. Assumption 3 (Student) is in line with empirical observations in [START_REF] Mike | An empirical behavioral model of liquidity and volatility[END_REF]. Assumptions 4 and 7 are in line with empirical observations in [START_REF] Chakraborti | Econophysics: Empirical facts and agent-based models[END_REF] as far as the main body of the distribution is concerned, but fail to represent the broad distribution observed in empirical studies. All the parameters except δ, which we kept exogenous, can be more or less roughly estimated on our data. In fact δ is the parameter of the less realistic feature of this simple model, and is thus difficult to calibrate. It can be used as a free parameter to fit the realized volatility.

Adding dependence between order flows

We have found in section 1.2 that market data shows that the flow of limit orders strongly depends on the flow of market order. We thus propose that in our experiment, the flow of limit and market orders are modelled by Hawkes processes N L and N M , with stochastic intensities respectively λ and µ defined as:

       µ M (t) = µ M 0 + t 0 α MM e -βMM (t-s) dN M s λ L (t) = λ L 0 + t 0 α LM e -βLM (t-s) dN M s + t 0 α LL e -βLL(t-s) dN L s (3)
Three mechanisms can be used here. The first two are self-exciting ones, "MM" and "LL". They are a way to translate into the model the observed clustering of arrival of market and limit orders and the broad distributions of their durations. In the empirical study [START_REF] Large | Measuring the resiliency of an electronic limit order book[END_REF], it is found that the measured excitation MM is important. In our simulated model, we will show (see 3.2) that this allows a simulator to provide realistic distributions of the durations of trades.

The third mechanism, "LM", is the direct translation of the empirical property we have presented in section 1.2. When a market order is submitted, the intensity of the limit order process N L increases, enforcing the probability that a "market making" behaviour will be the next event. We do no implement the reciprocal mutual excitation "ML", since we do not observe that kind of influence on our data, as explained in section 1.3.

Rest of the model is unchanged. Turning these features successively on and off gives us several models to test -namely HP (Homogeneous Poisson pro-cesses), LM, MM, MM+LM, MM+LL, MM+LL+LM -to try to understand the influence of each effect.

3 Numerical results on the order book

Fitting and simulation

We fit these processes by computing the maximum likelihood estimators of the parameters of the different models on our data. As expected, estimated values varies with the market activity on the day of the sample. However, it appears that estimation of the parameters of stochastic intensity for the MM and LM effect are quite robust. We find an average relaxation parameter βMM = 6, i.e. roughly 170 milliseconds as a characteristic time for the MM effect, and βLM = 1.8, i.e. roughly 550 milliseconds characteristic time for the LM effect. Estimation of models including the LL effect are much more troublesome on our data. In the simulations that follows, we assume that the self-exciting parameters are similar (α MM = α LL , β MM = β LL ) and ensure that the number of market orders and limit orders in the different simulations is roughly equivalent (i.e. approximately 145000 limit orders and 19000 market orders for 24 hours of continuous trading). Table 2 summarizes the numerical values used for simulation. Fitted parameters are in agreement with an assumption of asymptotic stationarity.

Model

µ0 αMM βMM λ0 αLM βLM αLL βLL HP 0.22 --1.69 ----LM 0.22 --0.79 5.8 1.8 --MM 0.09 1.7 6.0 1.69 ----MM LL 0.09 1.7 6.0 0.60 --1.7 6.0 MM LM 0.12 1.7 6.0 0.82 5.8 1.8 --MM LL LM 0.12 1.7 5.8 0.02 5.8 1.8 1.7 6.0 Common parameters: m P 1 = 2.7, ν P 1 = 2.0, s P 1 = 0.9

V 1 = 275, m V 2 = 380 λ C = 1.35, δ = 0.015
Table 2. Estimated values of parameters used for simulations.

We compute long runs of simulations with our enhanced model, simulating each time 24 hours of continuous trading. Note that using the chosen parameters, we never face the case of an empty order book. We observe several statistics on the results, which we discuss in the following sections.

Impact on arrival times

We can easily check that introducing self-and mutually exciting processes into the order book simulator helps producing more realistic arrival times. Figure 5 shows the distributions of the durations of market orders (left) and limit orders (right). As expected, we check that the Poisson assumption has to be discarded, while the Hawkes processes help getting more weight for very short time intervals. Fig. 5. Empirical density function of the distribution of the durations of market orders (left) and limit orders (right) for three simulations, namely HP, MM, LL, compared to empirical measures. In inset, same data using a semi-log scale.

We also verify that models with only self-exciting processes MM and LL are not able to reproduce the "liquidity replenishment" feature described in section 1.2. Distribution of time intervals between a market order and the next limit order are plotted on figure 6. As expected, no peak for short times is observed if the LM effect is not in the model. But when the LM effect is included, the simulated distribution of time intervals between a market order and the following limit order is very close to the empirical one.

Impact on the bid-ask spread

Besides a better simulation of the arrival times of orders, we argue that the LM effect also helps simulating a more realistic behaviour of the bid-ask spread of the order book. On figure 7, we compare the distributions of the spread for three models -HP, MM, MM+LM -in regard to the empirical measures. We first observe that the model with homogeneous Poisson processes produces a fairly good shape for the spread distribution, but slightly shifted to the right. Small spread values are largely underestimated. When adding the MM effect in order to get a better grasp at market orders' arrival times, it appears that we flatten the spread distribution. One interpretation could be that when the Fig. 6. Empirical density function of the distribution of the time intervals between a market order and the following limit order for three simulations, namely HP, MM+LL, MM+LL+LM, compared to empirical measures. In inset, same data using a semi-log scale. process N M is excited, markets orders tend to arrive in cluster and to hit the first limits of the order book, widening the spread and thus giving more weight to large spread values. But since the number of orders is roughly constant in our simulations, there has to be periods of lesser market activity where limit orders reduce the spread. Hence a flatter distribution.

Here, the MM+LM model produces a spread distribution much closer to the empirical shape. It appears from figure 7 that the LM effect reduces the spread: the "market making" behaviour, i.e. limit orders triggered by market orders, helps giving less weight to larger spread values (see the tail of the distribution) and to sharpen the peak of the distribution for small spread values. Thus, it seems that simulations confirm the empirical properties of a "market making" behaviour on electronic order books.

We show on figure 8 that the same effect is observed in an even clearer way with the MM+LL and MM+LL+LM models. Actually, the spread distri- bution produced by the MM+LL model is the flattest one. This is in line with our previous argument. When using two independent self exciting Hawkes processes for arrival of orders, periods of high market orders' intensity gives more weight to large spread values, while periods of high limit orders' intensity gives more weight to small spread values. Adding the cross-term LM to the processes implements a coupling effect that helps reproducing the empirical shape of the spread distribution. The MM+LL+LM simulated spread is the closest to the empirical one.

A remark on price returns in the model

It is somewhat remarkable to observe that these variations of the spread distributions are obtained with little or no change in the distributions of the variations of the mid-price. As shown on figure 9, the distributions of the variations of the mid-price sampled every 30 seconds are nearly identical for all the simulated models (and much tighter than the empirical one). This is due to the fact that the simulated order books are much more furnished than the empirical one, hence the smaller standard deviation of the mid price variations. One solution to get thinner order books and hence more realistic values of the variations of the mid-price would be to increase our exogenous parameter δ. But in that case, mechanisms for the replenishment of an empty order book should be carefully studied, which is still to be done.

Conclusion

We have shown the the use of Hawkes processes may help producing a realistic shape of the spread distribution in an agent-based order book simulator. We emphasize on the role of the excitation of the limit order process by the market order process. This coupling of the processes, similar to a "market making" behaviour, is empirically observed on several markets, and simulations confirms it is a key component for realistic order book models.

Future work should investigate if other processes or other kernels (ν LM in our notation) might better fit the observed orders flows. In particular, we observe very short characteristic times, which should lead us to question the use of the exponential decay. Furthermore, as pointed out in the paper, many other mechanisms are to be investigated: excitation of markets orders, link with volumes, replenishment of an empty order book, etc.

Fig. 1 .

 1 Fig.1. Empirical density function of the distribution of the bid-ask spread in event time and in physical (calendar) time. In inset, same data using a semi-log scale. This graph has been computed with 15 four-hour samples of tick data on the BNPP.PA stock (see section 1.1 for details).

1. 1

 1 Data and observation setupWe use order book data for several types of financial assets: • BNP Paribas (RIC 1 : BNPP.PA): 7th component of the CAC40 during the studied period • Peugeot (RIC: PEUP.PA): 38th component of the CAC40 during the studied period • Lagardère SCA (RIC: LAGA.PA): 33th component of the CAC40 during the studied period • Dec.2009 futures on the 3-month Euribor (RIC: FEIZ9) • Dec.2009 futures on the Footsie index (RIC: FFIZ9)

Fig. 7 .

 7 Fig.7. Empirical density function of the distribution of the bid-ask spread for three simulations, namely HP, MM, MM+LM, compared to empirical measures. In inset, same data using a semi-log scale. X-axis is scaled in euro (1 tick is 0.01 euro).

Fig. 8 .

 8 Fig.8. Empirical density function of the distribution of the bid-ask spread three simulations, namely HP, MM, MM+LM, compared to empirical measures. In inset, same data using a semi-log scale. X-axis is scaled in euro (1 tick is 0.01 euro).

Fig. 9 .

 9 Fig.9. Empirical density function of the distribution of the 30-second variations of the mid-price for five simulations, namely HP, MM, MM+LM, MM+LL, MM+LL+LM, using a semi-log scale. X-axis is scaled in euro (1 tick is 0.01 euro).

Table 1 .

 1 Number

	Code Number of limit orders Number of market orders
	BNPP.PA	321,412	48,171
	PEUP.PA	228,422	23,888
	LAGA.PA	196,539	9,834
	FEIZ9	110,300	10,401
	FFIZ9	799,858	51,020

of limit and markets orders recorded on 15 samples of four hours (Sep 10th to Sep 30th, 2009 ; 9:30am to 1:30pm) for 5 different assets (stocks, index futures, bond futures)