Nicolas Lutz

Basile Sauvage

Jean-Michel Dischler

Cyclostationary Gaussian noise: theory and synthesis

Keywords:

Figure 1: Cyclostationary noise is a procedural model for cyclostationary textures, which are patterns exhibiting a periodicity in their statistics. We convey existing stationary noises to a cyclostationary context, enabling the synthesis of cyclostationary textures controlled by spectra (left) and by an exemplar (right).

Introduction

Procedural noise has many applications in computer graphics, ranging from geometric modeling and animation to texture synthesis, as it is able to compactly represent continuous stochastic stationary processes. This explains the many research works that noise generated during the past decades. However, the physical nature of some natural or man-made phenomena often cause the assumption of stationarity to be violated. Modeling non-stationary stochastic processes is more difficult, and such processes have not been well studied yet.

A key application of noise is procedural texture modeling. Purely random patterns can be directly generated with noise; an important class of such patterns are textures that are realizations of stationary Gaussian processes. Recent noise synthesis approaches attempt to enlarge the range of patterns that noise is able to reach directly, that is, with no additional procedural modeling technique, like local random-phase [GSV + 14], high-performance [START_REF] Heitz | High-performance by-example noise using a histogram-preserving blending operator[END_REF] and Phasor [TEZ + 19] noise. Our goal is to follow on from these works and further increase the range of patterns reachable with a stand-alone procedural noise.

Therefore, we study a new class of non-stationary stochastic processes, called cyclostationary stochastic processes, which have not been studied yet in the field of computer graphics and procedural noise synthesis to the best of our knowledge. Such processes, characterized by the periodicity of their statistics, are well known in signal processing. They are used to model periodically-correlated random signals such as signals resulting from rotating machines, radar systems, or weather patterns. We introduce novel noises, called cy- clostationary Gaussian noises, which we define as models and algorithms implementing cyclostationary Gaussian processes. These noises can compactly and procedurally represent a class of textures we subsequently call cyclostationary noise textures, defined by their periodically random pattern. The periodicity and repetition of these textures is not obvious and it can be easier to perceive their random nature. We showcase examples of cyclostationary textures in Figure 2: they exhibit strong stochastic variations, best observable when zooming in, while having a repetitive spatial organization, best observable when zooming out.

Previous techniques introduced a class of textures, called nearregular textures [START_REF] Liu | Near-regular texture analysis and manipulation[END_REF], defined as periodic tilings with geometry and color variations. With this respect, our work provides a mathematically grounded, novel and complementary type of color variations, which can be strongly stochastic and detailed. Moreover, our cyclostationary textures can be synthesized in real-time and are unbounded, i.e. they have an infinite expanse.

First, we explore the mathematical properties of cyclostationary Gaussian processes, their statistics, and their link with cyclostationary textures in section 3. We then present various ways of generating patterns using classic stationary processes, notably with a new spatially-varying spectrum model, put in practice by adapting realtime procedural noises known as local random phase noise, Gabor noise and phasor noise in section 4. Then, we show how we achieve a high-quality by-example synthesis of cyclostationary textures by adapting the real-time procedural noises known as spot noise and high-performance noise as well as a cyclostationary version of histogram transfer in section 5. We then show how to estimate the period of cyclostationary textures in section 6; finally, we discuss details, comparisons and limitations in section 7.

Related work

Our work is related to noise and Gaussian processes, to the theory of cyclostationary processes and to near-regular textures. We briefly review these topics.

Near-regular textures. Near-regular textures have been introduced to address a specific class of repetitive textures that could not be sufficiently well generated with aperiodic tilings [Sta97, show that controlling the position of these spots helps controlling the visual aspect of the output. High-performance noise [HN18, DH18, Bur19] is a very fast method approximating the ADSN model. It blends hexagonal tiles taken randomly from a discrete exemplar and arranged on a regular lattice. The small number of blended values enables the preservation of non-Gaussian features, but diminishes the variety of the outputs. These noise models can often be adapted to the by-example case, either by using the spectrum of the exemplar [GGM11, GLLD12, GSV + 14], or, in the case of the ADSN model, by using the exemplar to build a spatial kernel [GGM11, GLM17, HN18].

All these models share issues related to the stationarity of the underlying Gaussian process that they build upon. They are thus unable to synthesize periodic statistics, and subsequently unable to model color transformations in near-regular textures. Works attempting to control the spatial organization of the underling stationary process do so by either incidentally altering the stationarity [GSDC17,PGDG16,CGG19], by varying the sparse convolution using specific point processes; or, by altering their Gaussian nature through the control of the phase [GSV + 14,TEZ + 19]; or by reducing the convergence towards a Gaussian process [START_REF] Heitz | High-performance by-example noise using a histogram-preserving blending operator[END_REF]. Conversely, we propose to alterate stationarity and keep the Gaussian nature by exploiting a generalization to Gaussian cyclostationary processes.

Cyclostationary processes. Cyclostationary processes are periodic processes that can be seen as interleaved stationary processes (and, as such, are a generalization of stationary processes). They were first presented in [Ben58], who noted that the synchronization algorithms used in communication systems exhibited nonstationary signals with a periodic mean and a periodic autocovariance function. Since then, they have been used in various contexts (such as gear systems, mechanical machinery, econometric, or astronomy), but, to our knowledge, have never yet been exploited in texture synthesis in an explicit way. Our presentation of cyclostationary processes and their properties is based on the recent work of Napolitano [Nap19], a book summarizing 60 years of research on cyclostationary processes and their properties. The work of Kipnis et al. [START_REF] Kipnis | The distortion rate function of cyclostationary gaussian processes[END_REF] was also used, as it establishes the link between our spatially-varying spectrum and the autocovariance function of cyclostationary noises.

It is important to note that even though cyclostationary random fields have been defined, the cyclostationary theory has been mostly elaborated for one-dimensional processes. As such, some of the cyclostationary theory presented in this paper is a generalization of the 1D model.

Cyclostationary model

The main idea of this paper is to consider cyclostationary (CS) Gaussian processes as generators of procedural noise, characterized by the periodicity of its statistics. This section presents their respective properties and the link between them.

The term "noise" is well defined as a stationary random function characterized by its power spectral density [LLC + 10]. In computer graphics, the terms random process, noise and textures are sometimes merged; moreover, despite similarities with classic noise, our model is characterized by additional parameters. To avoid confusions between these three terms, we use them as follows:

• "Process" refers to a general random process, defined through its mathematical properties, as in section 3.1. • "Noise" refers to a specific process, defined by a computable formula, along with algorithms to compute it. We explore several noises in Sections 4 and 5. • "Texture" refers to an actual image, one result of an algorithm, and one realization of a process.

Cyclostationary processes

Processes. A stochastic process Y is a set of random variables Y (x) defined on a probability space for which the sample space is a subset of R d , and x lie in Ω. We focus on 2D random fields, thus Ω = R 2 . Their realization can produce grayscale (d = 1) or color (d = 3) textures.

Gaussian processes. A Gaussian process Y is a stochastic process for which every linear combination of Y (x 1 ), ...,Y (x N ) such that (x 1 , ...x N ) ∈ Ω N is a random variable with a normal distribution. The first order moment and second order moment of these processes are necessary and sufficient to define their N-th order statistics for all N, that is, the probability distribution of the entire process.

Stationarity and cyclostationarity. The invariance of the sets of probabilities P under a translation t ∈ Ω is defined by

P(Y (x 1 ), .. Y (x N )) = P(Y (x 1 + t), .. Y (x N + t)). (1)
If equation ( 1) is satisfied for any t, then Y is said stationary. Y is said to be a cyclostationary (CS) or periodically correlated process with a period vector t 0 if equation ( 1) is satisfied for t = kt 0 , ∀k ∈ Z, where t 0 represents the period of the CS process, and there are no vectors t 0 with any smaller coordinate larger than 0 for which this invariance is true [Nap19]. Most of the time throughout this paper, CS processes are assumed to be periodic with not one, but two vectors t 0 and t 1 which form a basis of Ω.

As shown in Figure 3, t 0 and t 1 define a parallelogram Ω -and a lattice tiling the plane Ω. We denote as P the operator Ω → Ω - that projects any x modulo (t 0 , t 1 ):

x = P(x) + kt 0 + lt 1 , with P(x) ∈ Ω -and (k, l) ∈ Z 2 (2)
A stationary process can be seen as a CS process where t 0 and t 1 are infinitely small, making CS processes a generalization of the stationary model. Likewise, a CS process with only one period vector t 0 can be seen as one with an orthogonal t 1 infinitely small if it is stationary in this direction, or infinitely large if it has arbitrary statistics.

Ergodicity. A stationary process is said to be ergodic when its statistical properties can be estimated on a single, sufficiently large realization, rather than several realizations of the process. This property typically enables an accurate by-example texture synthesis, because the statistics of the process can be estimated from the example. The extension to CS processes is known as cycloergodicity. We assume all the processes to be cycloergodic, and we will make use of this property to estimate the CS statistics.

Cyclostationary statistics

As explained previously, the statistics of a cyclostationary process Y are periodic. All statistics presented here have an echo in stationary processes, as they only take an additional argument x ∈ Ω, due to the spatially-varying nature of CS processes.

The first order moment M(x) E[Y (x)] is the expectation of the process Y in x. The cyclostationarity implies M(x) = M(P(x)). The cycloergodicity implies that it can be estimated from an exemplar E by an empirical mean

M(x) = 1 #samples ∑ k,l E(x + k t 0 + l t 1 ) (3)
where #samples is the number of samples x + k t 0 + l t 1 in E. This equation is later used for our cyclostationary ADSN model in section 4.3 and can be visualized in the bottom row of figure 6.

The second order central moment is the autocovariance function of Y , defined by R(x, t) E[(Y (x) -M(x)) (Y (x + t) -M(x + t))] for any position x and any translation t in Ω. The cyclostationarity implies R(x, t) = R(P(x), t). The cycloergodicity implies that it can be estimated from an exemplar E by

R(x, t) = 1 #samples ∑ k,l E(y) -M(x) E(y + t) -M(x + t) , (4
) where y = x + k t 0 + l t 1 and #samples is the number of samples for which both y and y + t fall inside E. It is computed independently for each color channel.

In the stationary case, the autocovariance is often controlled in the spectral domain by its Fourier transform, i.e. the power spectral density (PSD). In the cyclostationary case, the most straightforward formula is the instantaneous power spectral density (IPSD), i.e. the Fourier transform of R with respect to t. However, another, more effective approach is possible using the polyphase components.

Polyphase components. It is possible to decompose a cyclostationary signal into stationary components [Nap19]. One possibility is the translation time series decomposition, or its discrete form, called the polyphase decomposition. For any x ∈ Ω -, the collection of positions {x + kt 0 + lt 1 , (k, l) ∈ Z 2 } congruent to x define a discrete sub-process Y x called the x-th polyphase component of the process Y (also called x-th decimated component), as shown in Figure 3 left. This decomposition is helpful because Y x are stationary processes: their histogram H x capture the first order statistics (Figure 3 right), which we may want to control (see section 5.3); their power spectral densities S x are intuitive local spectra (see sections 3.3 and 4.2), that control the second order moment R [START_REF] Kipnis | The distortion rate function of cyclostationary gaussian processes[END_REF].

Since our processes are Gaussian, wide-sense and strict-sense cyclostationarity are equivalent [Lap17], which implies that Y is entirely defined by M and R, or by all S x .

Textures

We model a texture as a function I : Ω → R d , with d the dimension of a random vector returned by noise (d = 1 for grayscale, d = 3 for color).

Cyclostationary textures. We define cyclostationary textures (CS textures) as realizations of arbitrary cyclostationary processes, implying they have periodic statistics, including but not limited to their mean and autocovariance function; a plausible set of period vectors and lattice by which the CS textures can be generated by a CS process also appear clearly when examining CS textures, as Figure 3 shows. The parameters of the underlying CS process can be estimated as long as the CS texture is sufficiently large, and the CS process is assumed cycloergodic.

Stationary random phase textures. Random phase textures are the result of a random phase noise. They are also commonly called "Gaussian textures", along with the realizations of other types of noises, because they often have statistics similar to those of the realizations of a Gaussian process, although this depends on the spectrum of the noise. They are defined by Galerne et al. [GGM11] as textures having a random phase in the Fourier domain:

I(x) = ∑ ξ A(ξ) cos (2πξx + Φ(ξ)) (5)
for all frequencies ξ in the frequency domain Ω, with random phases Φ uniformly distributed in [0, 2π[. These types of textures are interesting, because they can be compactly modeled in a procedural way using procedural noises.

Cyclostationary random phase textures. Similarly, we define cyclostationary random phase (CSRP) textures as the realizations of CS random phase noise processes, to model cyclostationary textures through a random phase:

I(x) = ∑ ξ A(x, ξ) cos (2πξx + Φ(ξ)). ( 6
)
The novelty is that the amplitude A : Ω × Ω → R is spatiallyvarying, i.e. it is a periodic function of x (periods t 0 and t 1 ). The phase Φ is random, while A entirely defines the statistics of the CSRP texture. In this scenario, A(x, ξ) 2 is equal to the power spectral density S x (ξ) of the x-th polyphase component.

Although not all CS textures are CSRP textures, this model approximates them well in the same way that random phase textures approximate "stationary" textures (called micro-textures in Galerne et al. [GGM11]), and establishes the possibility of synthesizing textures by defining a spatially-varying amplitude and randomizing their phase. Note that random phase textures and Gaussian textures are very similar but not equivalent [GGM11]. For simplicity, in the CS case, we call them CS noise textures too.

Cyclostationary Gaussian noises

In this section, we introduce cyclostationary Gaussian noises. They consist in different formulations of cyclostationary Gaussian processes, giving rise to different synthesis algorithms of cyclostationary noise textures.

Modulation of a stationary signal

Cyclostationary processes can be expressed as a simple aggregation of stationary components periodically defined on Ω [Nap19]. A simple way to obtain a cyclostationary signal is therefore to modulate a set of jointly stationary signals {Y 1 , ...,Y K } as

I(x) = ∑ s ws(x)Ys(x), ∀x ∈ Ω, (7)
where ws : Ω → R are (t 0 , t 1 )-periodic weighting functions.

This approach was used by Guingo et al.

[GSDC17] to synthesize heterogeneous textures from an input example, from which they extract a few stationary signals along with weighting functions. It can also be found to a lesser extent in Haindl et al. [START_REF] Haindl | Near-regular texture synthesis[END_REF], where near-regular textures are considered as a set of cells inside which a single noise is extracted and synthesized. The drawback of this approach is the definition of the weights ws and the stationary signals Ys: a manual design is tedious, while an automated extraction is difficult.

Cyclostationary random phase noise

We have established in equation 6 that cyclostationary random phase textures can be expressed as a sum of cosines with random phases, with statistics being controlled through the spatiallyvarying amplitude A. This paves the way to adapt existing stationary procedural noises based on the PSD: random phase noise [GGM11, GSV + 14], Gabor noise [LLDD09] and phasor noise [TEZ + 19]. The main difficulty is the control of the "spatiallyvarying PSD" A 2 (x, •) = S x .

We validate this idea in Figure 4 by a direct synthesis in the Fourier domain, where A 2 is linearly interpolated from a small set of spectra defining a periodic variation. In Figure 1 left, we adapted the Gabor noise, the local random-phase noise, and the phasor noise. In this example, A 2 is a periodic function, which makes a bi-lobe vary in orientation according to a single period t 0 (red vector). Note that phasor noise is actually not Gaussian; however, it is based on a Gaussian process, so we decided to present it here.

Cyclostationary asymptotic discrete spot noise

The stationary asymptotic discrete spot noise (ADSN) [GGM11] is a stationary and ergodic Gaussian process that computes an output I of arbitrary size as the sparse convolution of a discrete spot J:

I = 1 |J| (J -µ 1) * W, (8)
where W is the realization of a white noise made of uncorrelated components with distribution N (0, 1) and µ is the estimated mean of the discrete spot J, 1 is a constant function equal to 1 with the same support as J, and |J| is the number of texels of the discrete spot.

We introduce a new noise model, the cyclostationary asymptotic discrete spot noise (CS-ADSN). It is a cyclostationary and cycloergodic Gaussian process that computes an output texture I from a spot J as

I = 1 |J| (J -M) * W L , (9)
where M is the first order moment of the cyclostationary process described in section 3.2, and W L is a discrete process defined on the lattice whose nodes are congruent to 0 with respect to the periods t 0 and t 1 . W L is a sum of diracs (one per node) with uncorrelated random magnitudes distributed according to N (0, |Ω -|). I is a process with first order moment equal to 0, requiring the addition of M after the process. The second order moment of I corresponds to the second order moment of the spot, which we prove in appendix A.

This model is close to a specific kind of locally controlled spot noise presented in Pavie et al. [PGDG16] and Cavalier et al. [CGG19]

, where the results can be made regular by controlling the position where the spots are shot: we conjecture that our CS-ADSN is the underlying cyclostationary model for the perfectly regular version of the controlled spot noise they present without consideration for the cyclostationary mean and the variance of the output, which enables us to go one step further and to achieve a byexample synthesis for cyclostationary textures in the next section.

Figure 5: Exemplars (top row) synthesized with our cyclostationary spot noise (middle) and our cyclostationary high-performance noise (bottom). Our first order statistics transfer (section 5.3) was applied to all results. All textures are unbounded and continuous, and can be sampled in real time at any location of infinite space. Here, we show a sampling matching the resolution of the input.

By-example synthesis

The parameters of the noise models can be difficult to define manually. Instead, noise by-example consists in estimating the parameters from an input example image E. In our context, the idea is to find an underlying process Y which could have been likely to generate the example.

The cyclostationary random phase noise of equation 6 is unfortunately not a convenient model for by-example synthesis, because it requires an accurate estimation of the PSD S x of polyphase components; however, we dispose of only a few sparse samples in E for each polyphase component. Because of this, we choose to elaborate on the CS-ADSN model presented in section 4.3. We applied it to two procedural noises using an exemplar as a discrete spot: spot noise (section 5.1), and high-performance noise (section 5.2), with results of these algorithms shown in Figure 5.

In this section, we assume that t 0 and t 1 are already estimated from an exemplar; we show how to do this in section 6.

Cyclostationary spot noise

The by-example stationary spot noise algorithm [GGM11] is a direct implementation of the ADSN model of equation (8) that consists in shooting E (the spot J = E) on random positions. The random positions are determined by a sparse Poisson process [vW91] instead of a white noise W , which increases the space coverage. As shown in Figure 6 top, it is unable to reproduce the pattern regularity.

We define the cyclostationary spot noise from equation 9, by shooting E only on the regular lattice described by t 0 and t 1 , as shown in Figure 6 bottom. Since W L is sparse, it can be faster to compute the cyclostationary spot noise directly instead of simulating it by shooting E several times. For this process, the first order moment M needs to be estimated from the exemplar, which can 3) are the critical steps for preserving the statistics of the exemplar. The addition of µ and M required after equations (8) and (9) is omitted. be done with Equation 3. The asymptotic result of this by-example noise takes the form of a Gaussian cyclostationary process with the statistics of the CS-ADSN of E.

Cyclostationary high-performance noise

The stationary high-performance noise algorithm [HN18, DH18, Bur19] is a real-time by-example synthesis designed as a compromise between a tiling algorithm and a noise algorithm. It consists in blending hexagon-shaped spots centered on the vertices of a triangular grid. For each vertex x of the grid, an hexagonal tile is taken from the exemplar E at a random position given by a hash function h(x). The overlapping tiles are linearly blended within each triangle. As shown in Figure 7 left, it fails to reproduce CS textures. The reason is that x and h(x) are not congruent.

We define a cyclostationary variant of the high-performance noise by controlling the positions of the tiles. The hash function h(x) returns a random position congruent to x, i.e. P(h(x)) = P(x) (see yellow vectors in Figure 7). Additionally, it is compatible with the enhancement proposed by Burley [Bur19] that modifies the blending weights.

First order statistics transfer

The previous algorithms are designed for textures with cyclostationary Gaussian statistics. However, many natural stochastic textures do not have Gaussian statistics. Therefore, many works tried to push the algorithms beyond the pure Gaussian scope of application. Among them, optimal transport has been used successfully to control the first order statistics [GLR17, HN18, DH18], namely to process textures with non Gaussian histograms.

In the following, we briefly recall stationary histogram transfer, show its limitations, propose a CS histogram transfer, discuss robustness issues and filtering.

Stationary histogram transfer

Let E H be an exemplar with non Gaussian histogram H. The principle is as follows:

• Pre-compute a transfer function T that maps H onto a Gaussian histogram G. • Pre-compute T -1 in a lookup table.

• Pre-compute E G = T (E H ), a "Gaussianized" example.

• At runtime, compute the result T -1 (Synthesis(E G )) by applying the synthesis algorithm to E G and compositing by T -1 .

The results has H as non-Gaussian histogram, just as the example E H . When applied to a CS texture, this procedure produces artifacts, as shown in Figure 8. These artifacts are due to the global nature of the histogram, which does not capture spatial variations.

Cyclostationary histogram transfer

As explained in section 3, a CS process can be decomposed into a set of polyphase components (PC) which are stationary processes. In this section, the intuition is to preserve the histogram of each PC independently.

Let E be a CS texture we want to synthesize, which we assume to be the result of an underlying CS process. Let E x be the polyphase components estimates: for any x ∈ Ω -, it is an estimation, from the example E, of the x-th PC of the underlying process. In concrete terms, E x is a subsampling of E with |E|/|Ω -| samples. Let H x be the histogram of E x (Figure 3), which is an estimation of As shown in Figure 8, this cyclostationary transfer helps against false colors. It notably removes region bleeding (when the output color is incorrectly interpreted as that of a different region of the exemplar), and reduces ghosting when attempting to synthesize textures with slight geometric irregularities. The inverse transfer functions T -1 P(x) are stored for all texels x ∈ Ω -. In practice, the memory footprint is at most equivalent to the footprint of the exemplar itself.

Robust histogram estimation

The CS histogram transfer (previous section) relies on the estimations H x . However, for small examples and large periods, |E x | becomes small, causing issues in the histogram estimation and artifacts in the result (see Figure 9). To make our estimation more robust, we assume a local stationarity: in a small neighborhood of size δ around x we assume the PC to be similar. Practically, we add into E x the samples of the neighboring PC estimates: δ = 1 enlarge the sampling by a factor 9, δ = 2 by a factor 25, etc. As shown in Figure 9, extending the neighborhood yields a compromise between the stationary (wrong colors, visible blur) and cyclostationary (better colors, visible noise) transfer. In the limit (as the neighborhood covers the example), all estimates are equal and it boils down to a stationary transfer.

Real-time rendering and filtering

The synthesis algorithms presented in this paper are tailored for real-time synthesis and rendering, while keeping a low memory load. It requires to filter the resulting texture in real-time. Both the spot noise [vW91] and the high-performance noise [START_REF] Heitz | High-performance by-example noise using a histogram-preserving blending operator[END_REF] can be filtered using a pre-computed MIP-map of the example E. This technique is still valid with our cyclostationary model. We directly applied the histogram transfer on the result, as Heitz & Neyret [START_REF] Heitz | High-performance by-example noise using a histogram-preserving blending operator[END_REF]. Deliot et al. [START_REF] Deliot | Procedural stochastic textures by tiling and blending[END_REF] improve the results by computing one pre-filtered lookup table per MIP level: it can be adapted with one lookup table per level per polyphase component, however its practical implementation remains to be investigated.

Estimation of period vectors

In this section, we investigate the estimation of the period vectors (t 0 , t 1 ) from an exemplar E. As stated in the related works section, many methods have been proposed to estimate manually, semiautomatically or completely automatically repetitive structures in patterns. These approaches are useful, but sometimes lack accuracy when combined with our by-example synthesis algorithms, which require a sub-pixel precision. When using them, polyphase components may be incorrectly estimated and cause ghosting artifacts similar to those observed in Figure 8. Furthermore, automatic methods tend to be designed for harder problems and rely on complex software such as neural networks, which are heavy for the task at hand. We propose to use simple techniques such as a manual estimation, which is then refined by exploiting the characteristics of cyclostationary textures.

Our refinement algorithm is based on the observation that (t 0 , t 1 ) represent the translations to the next element in the repeated pattern, which corresponds to another sample in the same polyphase components. We expect the PC to be strongly coherent, i.e. that any two points y and z in a given PC E x are similar. Thus, we define

d PC (t 0 , t 1 , E x ) = 1 |E x | 2 ∑ y∈E x ∑ z∈E x (E(y) -E(z)) 2 , (10)
which is a measure of the internal coherence of E x . Then we measure the global relevance of the vectors (t 0 , t 1 ) by

d(t 0 , t 1 ) = x∈Ω - d PC (t 0 , t 1 , E x ) dx. ( 11
)
To refine the initial estimation (t 0 , t 1 ), we locally minimize d using a stochastic gradient descent. A visualization of d is shown in Figure 10, where an optimal point is found, even on the third exemplar exhibiting geometric irregularity. Note that d is not convex; however, the procedure robustly converges to the local minimum, as long as the the initial estimation is precise enough.

Results and discussion

Implementation When a discrete exemplar needs to be sampled at any continuous location (e.g. for the CS-spot noise), it is bilinearlyinterpolated between the pixels. The CS-spot noise is implemented with the enhancement of Galerne et al.

[GGM11] that consists in smoothing the border of the exemplar to remove linear artifacts. As in the original stationary version [START_REF] Heitz | High-performance by-example noise using a histogram-preserving blending operator[END_REF], the CS-HPN is easier to implement with a periodic examplar, otherwise hexagons that cross the border must be avoided. The performances of the stationary and cyclostationary algorithms are similar, so we invite the reader to refer to the original papers to get detailed information about timings. Indeed, our algorithms inherit the pros and cons of the original algorithms. Our codes are publicly available in the ASTex library [ast17] and on shadertoy.

Synthesis of CS random-phase textures

We have extended to the CS context several noises from the sate of the art: Figure 4 is a direct implementation of random phase noise ; Gabor noise, local random-phase noise, and phasor noise are shown in Figure 1 (left). Some additional high resolution results are available in our supplemental material.

We designed the spatially-varying PSD using functions or linear interpolation between discrete PSD. Other interpolation operators could be used, such as b-spline interpolation or Wasserstein interpolation [BPC16]. An estimation of this spectrum from an exemplar would be a powerful tool. This is however a difficult task, as samples of polyphase components estimates are sparse. Moreover, Guingo et al. [GSDC17] show that the extraction of spatial variations of the spectrum is difficult.

By-example synthesis We adapted spot noise and highperformance noise. Many results on various patterns are shown in Figure 5. The textures are unbounded and continuous: they can be sampled in real-time at any location of an infinite plane. A histogram transfer extends the application to examples with non Gaussian histograms, and helps reducing artifacts (Figure 8).

We use Figure 11 to further investigate the comparison between these two algorithms: row 1 shows the exemplar, row 2 shows our CS-spot noise with no histogram transfer (to show what unaltered Gaussian cyclostationary results look like), and row 3 shows our CS high-performance noise with histogram transfer. The major strength of the CS-spot noise is to create variety. For instance, columns A and B exhibit colored tiles, for which the colors are close to the result of a normal variable: CS-spot noise (row 2) implicitly estimates the Gaussian CS process that could have generated this tiling and creates new colors that match that process. This result is only possible when the distribution of the colors of each polyphase component corresponds to the result of a Gaussian process, otherwise it creates false colors (columns C and D).

Real-life CS textures corresponding to the realization of a pure Gaussian CS process are pretty rare. The global appearance and periodicity of other CS textures is well kept by the CS-spot noise, but the fine details and complex structures of these textures can be lost by the process. As can be seen in columns E and F, our CS high-performance noise (row 3) better reproduces such features. This is due to the underlying tiling and blending algorithm, which, in counterpart, may create slight visual repetitions (red circles). On some other examples such as colored tiles (columns A, B, C and D), unwanted patterns are created by the linear blending.

To summarize, our CS-spot noise should be used when variety is especially important; our CS-HPN should instead by used when fine details are more important; and neither process is capable of preserving the appearance of CS textures for which the distribution of the colors on each polyphase component is not a Gaussian histogram.

Our cyclostationary histogram transfer may be used to correct the colors, as it is able to correct the first order statistics of the exemplar. However, it does not act on higher-order statistics, and neighboring transfers are independent, causing potentially unexpected color differences between neighbors, such as in column D of Figure 11.

Among other existing by-example noise algorithms, it would be worth investigating texton noise [GLM17]. However, a robust estimation of the cyclostationary texton is not trivial.

Comparison with near-regular textures Cyclostationary textures represent a subset of near-regular textures for which there are no, or very little, geometric alterations, but important color alterations. Liu et al [START_REF] Liu | Near-regular texture analysis and manipulation[END_REF] synthesize color transformations using principal component analysis (PCA). In our experiments, we observed that PCA often induces important losses of contrast and a higher loss of details compared to our CS-spot noise, as shown in Figure 12. Because they separates individual tiles, the output tiles need to be stitched together, making real-time rendering harder to obtain, unless a lower quality blending of tiles is used. Haindl et al. [START_REF] Haindl | Near-regular texture synthesis[END_REF] propose another offline synthesis algorithm for a subclass of CS textures, which are characterized by a stationary noise within a periodic structure. We show in Figure 13 that CS-spot noise better preserves finely structured and contrasted details. CS-HPN even better preserves complex details.

Comparison with optimization and CNN techniques Until now, optimization and CNN techniques are computationally expensive and rarely suited for the rendering of unbounded textures, whereas speed and unbounded size are two core properties of our methods.

We show in Figure 14 a side by side comparison with "Deep correlations" [SCO17], which attempts to recover automatically repetitive structures in exemplars. Our methods better preserve fine details. More comparisons are given in the supplementary materials. Our model has the advantage of offering some guarantees: exemplars that fall in the class of cyclostationary noise textures are faithfully synthesized, and geometric regularity is perfectly reproduced. . This synthesis technique is not specifically addressing cyclostationary textures, and sometimes fails to capture the periodicity of statistics. Additionally, it is offline and works for small bounded textures: the opposite of our approach.

We believe CNN techniques could be helpful for detecting optimal period vectors according to more complex criteria, such as visual criteria, and thus replace the method described in section 6, or they could possibly provide a more robust way of estimating cyclostationary histograms.

Geometric variations

As of now, geometric variations are not handled with our approach. Liu et al [START_REF] Liu | Near-regular texture analysis and manipulation[END_REF] apply both geometric and color transformations, which we could a priori also do. Spatial warping of grids and barycentric coordinates are used to apply geometric alterations. However, we experienced that this works only well if colors have low contrasts and low frequencies. When colors are characterized by high frequencies and strong contrasts, spatial warping generates visual artifacts such as overstretching of details, unnatural distorted patterns and sampling artifacts. Devising an elaborated cyclostationary-like process, which considers both alterations simultaneously, seems like a difficult task that could be an interesting topic for future research.

Conclusion

In this paper, we introduced cyclostationary Gaussian noise, a new type of noise able to synthesize cyclostationary noise textures, characterized by the periodicity of their statistics. As a generalization of stationary noise, it broadens the range of patterns that noise is able to address. Thereby, a key application is the more efficient management of stochastic color transformations in near-regular textures. We defined a frame for our models by presenting cyclostationary processes and noises, their statistics and estimators, and cyclostationary textures. We provided a model to synthesize a cyclostationary texture by modulating a stationary process and by designing a spatially-varying power spectral density; we also provided a cyclostationary version of the ADSN model. This allowed us to generalize five procedural noises to the cyclostationary case: local random phase noise, Gabor noise, phasor noise, spot noise and highperformance noise. We showed how to execute a first order statistics correction on outputs generated with by-example synthesis in a cyclostationary context to extend the scope of Gaussian cyclo-stationary noise to a wider class of cyclostationary textures. We showed how to estimate the periods of an exemplar, required for by-example synthesis. Finally, we compared our synthesis to state of the art methods, and discussed limits and future works.

Figure 2 :

 2 Figure 2: Examples of cyclostationary (CS) textures, which exhibit a periodic organization when zooming out and, at the same time, strong stochastic color variations when zooming in. These types of near-regular textures are difficult to model, because colors are neither repetitive nor purely random: color randomness is periodic.

CSHD03,

  Wei04, LHW + 04] or patch-based synthesis [EF01, LLX + 01, KSE + 03, WY04, LH05]. Liu et al. [LLH04, LTL05] define near-regular textures as a periodic tiling that undergoes spatial, intensity and color transformations. Various works followed: evaluations and comparative studies [WHC + 06], analysis methods to detect regularities [DED05,CB11,PBCL09,CB13,LPVV17, LNS + 15], and improvements of synthesis quality [NMMK05, HH09, RHE11]. But managing simultaneously regularity and randomness of such textures still remains a challenging issue. When color variations are dominant (e.g. the resulting pattern is close to Gaussian noise), synthesis quality quickly breaks down, PCA [LLH04] not being an appropriate tool and patch-based methods not providing sufficient variety [RHE11]. We remark that none of these methods proposed a procedural noise model for handling color transformations more efficiently and in real-time. In the last years, texture synthesis has been carried through convolutional neural networks (CNN) [GEB15]. Some of the recent approaches attempt to automatically recover regularities and self-similarities in near-regular exemplars. Liu et al. [LGX16] use spectrum constraints and constraints learned from a CNN, while deep correlations [SCO17] use a structural energy, based on correlations and deep learning. These methods are not adapted to "on-the-fly" synthesis, and are restricted to medium size textures. Faster and unbounded texturing models have been proposed [ULVL16, GRGH20], but remain computationally expensive and are not suited for near-regular texture synthesis. Gaussian processes. Noise synthesis, carried through Gaussian processes, finds its roots in sparse convolution noise [Lew84]. Since then, different variations of this model have been proposed: Gabor noise [LLDD09, LLD11, TNVT19] sums a Gabor kernel with random weights to generate noise with controlled spectrum. Random phase noise [GGM11] uses Fourier series and a random phase to generate noise. Gilet et al. [GSV + 14] shows how fixing some phases can help preserve the spatial organization of an exemplar; however, freezing phases makes the output tend towards a perfectly periodic signal, which comes at the expense of variety. Guingo et al. [GSDC17] explores the possibility of using several spectra to provide a higher control over the output spatial organization, but without synthesizing near-regular textures. Phasor noise [TEZ + 19] generate patterns with strong variations of intensity, by controlling the instantaneous phase and a wave profile. Spot noise [vW91] is a notable implementation of the ADSN model presented later by Galerne et al. [GGM11], which consists in blending a discrete kernel at random positions. These works were improved with texton noise [GLM17], which computes a summarized version of the kernel, optimized for the synthesis. The work of Pavie et al. [PGDG16] and Cavalier et al. [CGG19]

Figure 3 :

 3 Figure3: A cyclostationary signal characterized by two periods t 0 and t 1 which define a regular lattice (yellow). Left: the signal can be decomposed into polyphase components (collection of blue dots, or magenta dots), which are discrete stationary processes. The xth polyphase component (PC) is made of all values on {x + kt 0 + lt 1 , (k, l) ∈ Z 2 } congruent to x. Right: histograms of the PCs are estimators of the statistics of the underlying process (provided it is cycloergodic).

Figure 4 :

 4 Figure 4: Cyclostationary random phase noise generated by Fourier series with a spatially-varying amplitude A (computed from an array of PSD) and a random phase. The array was constructed by bilinear interpolation of 9 spectra and defined periodically on the parallelogram described by t 0 and t 1 .

Figure 6 :

 6 Figure 6: Comparison between the stationary (top) and our cyclostationary (bottom) spot noise, using an exemplar texture as the spot. The estimation of the lattice and the mean M (Equation 3) are the critical steps for preserving the statistics of the exemplar. The addition of µ and M required after equations (8) and (9) is omitted.

Figure 7 :

 7 Figure 7: High-performance noise: stationary (left) and cyclostationary (right) version. For the latter, the positions of the hexagonal tiles in the example (h(x)) and in the result (x) are congruent.

Figure 8 :

 8 Figure 8: Histogram transfer: comparison between stationary (left) and cyclostationary (right) transfer with the HPN synthesis algorithm [HN18]. Our CS transfer helps reducing color artifacts.

Figure 9 :

 9 Figure 9: Robust histogram estimation, by collecting the polyphase component estimates in a small neighborhood of size δ. Extending δ trades noisy for blurry artifacts.

Figure 10 :

 10 Figure 10: Distance function d applied around the local minimum of d (corresponding to ground truth period vectors of t 0 and t 1 ) on various textures. d was only evaluated with orthogonal period vectors to produce this visualization.

Figure 11 :

 11 Figure 11: Exemplars (row 1) synthesized with a CS-spot noise and no histogram transfer (row 2), and with a CS-HPN with cyclostationary histogram transfer (row 3). The Gaussian nature of our CS-spot noise makes it work especially well with exemplars having cyclostationary Gaussian statistics (A,B). Our CS-HPN instead better preserves complex, non-Gaussian structures (E,F), but can produce repetitions (circled in red) and new aberrant tiles (A,B,C).

Figure 12 :

 12 Figure 12: Comparison with the tile synthesizer of Liu et al. [LLH04]. Tiles produced with their PCA technique (middle) show a loss of contrast and detail, and further require a restitching of tiles to produce output textures. Our real-time CS spot noise (right) avoids these drawbacks.

Figure 13 :

 13 Figure 13: Comparison with the offline synthesis algorithm of Haindl et al. [HH09] which reproduces a sub-class of CS textures. Our CS spot and CS high-performance noises are real-time, and better preserve finely structured and contrasted random patterns.

Figure 14 :

 14 Figure 14: Comparison with the CNN synthesizer of Sendik et al. [SCO17]. This synthesis technique is not specifically addressing cyclostationary textures, and sometimes fails to capture the periodicity of statistics. Additionally, it is offline and works for small bounded textures: the opposite of our approach.

c 2021 The Author(s) Computer Graphics Forum c 2021 The Eurographics Association and John Wiley & Sons Ltd.

Acknowledgments

This work has been partially funded by the project HDWorlds from the Agence Nationale de la Recherche (ANR-16-CE33-0001). We thank Xavier Chermain for his help with shadertoy. All textures are provided by https://textures.com, except for those of Figure 13 which are taken from Haindl et al. [HH09].

Appendix A: Autocovariance of the cyclostationary ADSN We prove the claim of section 4.3 that the second order moments of the spot J and its CS-ADSN I are identical. We prove it similarly to Galerne et al. [GGM11] for the stationary case. To enable an accurate estimation of the autocovariance on a spot, we consider that J is periodic. Let us denote by L the lattice described by the vectors t 0 and t 1 as in Figure 3. We call |L| the number of nodes of the lattice on the spot J.

The autocovariance of the spot J (which is the autocovariance of section 3.2 with the mean subtracted to each term) can thus be estimated by

∈ L or v / ∈ L, and |Ω -| otherwise. I is therefore a centered cyclostationary Gaussian process with second order moment R J (x, t).


Chieu Quang

Bui

Bassem Salem

Hervé Roussel

Xavier Mescot

Youssouf Guerfi

Carmen Jiménez

Vincent Consonni

Gustavo Ardila

Effects of Thermal Annealing on the Structural and Electrical Properties of ZnO Thin Films for Boosting their Piezoelectric Response

Keywords: Zinc Oxide, Metal-Organic Chemical Vapor Deposition, Annealing, Coalescence, Nanostructures, Piezoelectricity

The growth process and post-deposition treatment have strong effects on the structure of ZnO thin films, which can also significantly affect their electrical and piezoelectric properties. In this report, the effect of high temperature annealing under oxygen atmosphere on ZnO stacked thin films grown on silicon by pulsed-liquid injection metal-organic chemical vapor deposition is investigated in detail. By varying the annealing temperature from 600 to 1000 °C, the structural analysis reveals the occurrence of c-axis oriented grain coarsening and growth processes associated with the improvement of crystallinity and removal of carbon impurities. The electrical resistivity also strongly varies with thermal annealing. Those variations are accompanied by an enhancement of the piezoelectric amplitude measured by piezo-response force microscopy (PFM). Interestingly, the PFM measurements reveal the predominant development of O-and Zn-polar domains at low and high temperature annealing, respectively, originating from the anisotropic grain boundary velocity. The present findings show that the high temperature annealing under oxygen atmosphere in ZnO thin films can enhance not 2 only its piezo-amplitude but also its Zn-polarity uniformity, which unambiguously leads to the improvement of its piezoelectric performance.

Introduction

ZnO is not only an abundant, biocompatible material, but also a semiconductor compound possessing many interesting properties that can be used for various applications, [START_REF] Özgür | A Comprehensive Review of ZnO Materials and Devices[END_REF][START_REF] Kołodziejczak ; Radzimska | Zinc Oxide-From Synthesis to Application: A Review[END_REF] such as UV light emitting diodes and lasers, [START_REF] Willander | Zinc Oxide Nanorod Based Photonic Devices : Recent Progress in Growth , Light Emitting Diodes and Lasers[END_REF] UV photodetectors, [START_REF] Tian | Nanoscale Ultraviolet Photodetectors Based on One-Dimensional Metal Oxide Nanostructures[END_REF] electron transporting material and transparent conductive oxides in solar cells, [START_REF] Consonni | Nanowires for Solar Cells : A Comprehensive Review[END_REF] or piezoelectric components. [START_REF] Wang | Towards Self-Powered Nanosystems : From Nanogenerators to Nanopiezotronics **[END_REF] The use of ZnO in piezoelectric applications has attracted a lot of attention, in particular because ZnO is an eco-friendly, sustainable material and does not contain lead in contrast to the widely used lead zirconate titanate material. [START_REF] Wang | Towards Self-Powered Nanosystems : From Nanogenerators to Nanopiezotronics **[END_REF][START_REF] Wei | An Overview of Lead-Free Piezoelectric Materials and Devices[END_REF] The intrinsic piezoelectric property of ZnO results from its non-centrosymmetric wurtzite structure, which is its most stable phase under standard temperature and pressure conditions. [START_REF] Özgür | A Comprehensive Review of ZnO Materials and Devices[END_REF] ZnO can be integrated in the form of thins films or nanowires into piezoelectric devices, such as surface acoustic wave (SAW) applications, [START_REF] Yildirim | Comprehensive Characterization of ZnO Thin Films for Surface Acoustic Wave Applications[END_REF] nano-generators and pressure/strain sensors. [START_REF] Briscoe | Piezoelectric Nanogenerators -a Review of Nanostructured Piezoelectric Energy Harvesters[END_REF][START_REF] Tao | Performance of ZnO Based Piezo-Generators under Controlled Compression[END_REF][START_REF] Parmar | A Demonstration of the Mechanical Sensing Capability of Individually Contacted Vertical Piezoelectric Nanowires Arranged in Matrices[END_REF] However, the main drawback of ZnO for piezoelectric applications is related to its residual n-type doping. [START_REF] Wang | Towards Self-Powered Nanosystems : From Nanogenerators to Nanopiezotronics **[END_REF][START_REF] Liu | Interface Engineering on P-CuI/n-ZnO Heterojunction for Enhancing Piezoelectric and Piezo-Phototronic Performance[END_REF] The incorporation of many residual impurities (i.e. Al, Ga, In, H…) is typically achieved during the growth of ZnO, regardless of the deposition techniques used. [START_REF] Özgür | A Comprehensive Review of ZnO Materials and Devices[END_REF][START_REF] Janotti | Fundamentals of Zinc Oxide as a Semiconductor[END_REF] These residual impurities can act as shallow donors, generating a large amount of free electrons inside its structure. [START_REF] Özgür | A Comprehensive Review of ZnO Materials and Devices[END_REF][START_REF] Janotti | Fundamentals of Zinc Oxide as a Semiconductor[END_REF] Those free electrons can freely move and screen the piezoelectric potential, in turn reducing drastically its amplitude. [START_REF] Wang | Piezoelectric Nanogenerators Based on Zinc Oxide Nanowire Arrays[END_REF][START_REF] Tao | Unveiling the Influence of Surface Fermi Level Pinning on the Piezoelectric Response of Semiconducting Nanowires[END_REF] It is thus essential to reduce the amount of free electrons in order to increase the piezoelectric efficiency. Among the vast number of deposition techniques available to form ZnO thin films and nanowires, metal-organic chemical vapor deposition (MOCVD) offers the great advantages of forming ZnO deposits with a high crystallinity and of monitoring their morphology by simply adjusting the growth temperature, [START_REF] Park | Temperature-Induced Morphological Changes of ZnO Grown by Metalorganic Chemical Vapor Deposition[END_REF][START_REF] Malandrino | Relationship between Nanostructure and Optical Properties of ZnO Thin Films[END_REF][START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF] the flow rate or the ratio of zinc precursor and oxidizing agent. [START_REF] Fanni | Increasing Polycrystalline Zinc Oxide Grain Size by Control of Film Preferential Orientation[END_REF][START_REF] Aebersold | Quantifying Competitive Grain Overgrowth in Polycrystalline ZnO Thin Fi Lms[END_REF] This opens several approaches to tune the properties of ZnO deposits to be suitable for different applications. The morphology transitions from stacked thin films to nanowires through columnar thin films have been revealed in Ref. [START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF] by increasing growth temperature from 400 to 700 °C. ZnO thin films grown by MOCVD have further been shown to exhibit a high potential for piezoelectric applications. [START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF][START_REF] Emanetoglu | Epitaxial ZnO Piezoelectric Thin Films for Saw Filters[END_REF][START_REF] Khomchenko | Structure and the Emission and Piezoelectric Properties of MOCVD Grown ZnS , ZnS -ZnO , and ZnO Films[END_REF] However, further optimizations of the structural and electrical properties of ZnO thin films to enhance the piezoelectric amplitude are still required.

The effects of thermal annealing on the structural and electrical properties can lead to the performance enhancement of the ferroelectric materials such as PZT, [START_REF] Smith | PZT-Based Piezoelectric MEMS Technology[END_REF] (Bi 0.5 Na 0.5 )TiO 3 -(Bi 0.5 K 0.5 )TiO 3 -BaTiO 3 , [START_REF] Zhan | Effect of Annealing Temperature on Piezoelectric and Mechanical Properties of (Bi0.5Na0.5)TiO3-(Bi0.5K0.5)TiO3-BaTiO3 Thin Film[END_REF] and LiNbO 3 , [START_REF] Bhagavannarayana | Enhancement of Crystalline , Piezoelectric and Optical Quality of LiNbO3 Single Crystals by Post-Growth Annealing and Poling[END_REF] or the piezoelectric materials such as AlN, [START_REF] Vergara | Effect of Rapid Thermal Annealing on the Crystal Quality and the Piezoelectric Response of Polycrystalline AlN Films[END_REF] and ZnO. [START_REF] Lim | Crystal Growth and Piezoelectric Characterization of Mechanically Stable ZnO Nanostructure Arrays[END_REF] By using post-deposition thermal annealing under oxygen atmosphere, not only the crystallinity and caxis orientation of ZnO nano-pyramidal structures were improved, but also some defects were removed, resulting in better piezoelectric response as reported in Ref. [START_REF] Lim | Crystal Growth and Piezoelectric Characterization of Mechanically Stable ZnO Nanostructure Arrays[END_REF] . Many research groups have studied the effects of thermal annealing on the structural, optical, and electrical properties of ZnO thin films grown by different techniques such as sputtering, [START_REF] Fang | Influence of Post-Annealing Treatment on the Structure Properties of ZnO Films[END_REF][START_REF] Daniel | Effect of Annealing Temperature on the Structural and Optical Properties of ZnO Thin Films Prepared by RF Magnetron Sputtering[END_REF] molecular beam epitaxy, [START_REF] Ogata | Effects of Thermal Annealing of ZnO Layers Grown by MBE[END_REF] wet chemistry [START_REF] Richardson | Thermally Induced Pore Formation in Epitaxial ZnO Films Grown from Low Temperature Aqueous Solution[END_REF] or MOCVD. [START_REF] Gritsenko | Effect of Thermal Annealing on Properties of Polycrystalline ZnO Thin Films[END_REF] They have shown that the thermal annealing can improve the crystallinity of ZnO thin films and reduce the density of extended defects in their center. Under oxygen atmosphere, the thermal annealing causes oxygen chemisorption onto grain boundaries, which could lead to the decrease in free charge carrier density and mobility. [START_REF] Gritsenko | Effect of Thermal Annealing on Properties of Polycrystalline ZnO Thin Films[END_REF] These consequences would be favorable to improve the piezoelectric performance. However, the effects of thermal annealing on the piezoelectric coefficients, and especially the crystal polarity of ZnO thin films as well as their correlation with the structural and electrical properties have not been investigated yet. The uniformity of the crystal polarity in ZnO thin films is a critical condition for piezoelectric applications as the coexistence of domains with two opposite polarities can neutralize the piezoelectric output of each other. Unlike ferroelectric materials in which the polarity can be reversed and aligned by poling, ZnO is a piezoelectric material and cannot be poled. Therefore, it is crucial to form ZnO thin films not only with a high c-axis orientation, but also with a uniform crystal polarity in the final structure to enhance its piezoelectric performance.

In this work, we systematically study the effects of thermal annealing under oxygen atmosphere on the structural, electrical, and piezoelectric properties of ZnO thin films grown by pulsed-liquid injection MOCVD (PLI-MOCVD). By varying the annealing temperature from 600 to 1000 °C, the recrystallization and grain coarsening processes accompanied by the change of electrical resistivity, piezoelectric amplitude, and polarity distribution are described in detail. The structural properties of ZnO thin films are thoroughly investigated by using field-emission scanning electron microscopy (FESEM), X-ray diffraction (XRD), and Raman spectroscopy. Their electrical resistivity and piezoelectric properties are measured by transmission line measurement (TLM) and piezo-response force microscopy (PFM), respectively.

Experiment

ZnO thin film growth and thermal annealing ZnO thin films were grown by PLI-MOCVD in an Annealsys MC-200 MOCVD system.

Heavily doped p-type Si (100) wafer were used as substrates and as bottom electrode during PFM measurements. The zinc chemical precursor solution used was composed of diethylzinc (Zn(C 2 H 5 ) 2 , DEZn, Sigma-Aldrich) diluted in cyclohexane at 0.4 mol/l. It was injected to the reactor chamber at a flow rate of 0.5 g/min with 500 sccm Ar as gas carrier. The injection rate for zinc precursor was set at 3 Hz. The O 2 flow rate was set to 500 sccm, corresponding to an O/Zn ratio of about 100. The temperature of the Si substrate and the chamber pressure were maintained at 500 °C and 3 mbar during the growth, respectively. The ZnO thin film grown in 10 minutes has its thickness at of 1100 nm.

The as-grown ZnO thin films were then annealed under oxygen atmosphere at a temperature ranging from 600 to 1000 °C. The thermal annealing was carried out in a horizontal quartz tubular furnace under pure oxygen atmosphere. The pressure in furnace was close to 1 atm. The temperature was ramped up from room temperature to the annealing temperature with a rate of 25 °C/min. Then, the temperature was maintained for 1 hour. Finally, the temperature was cooled down to room temperature naturally.

Structural property measurements

The morphology of ZnO thin films before and after being annealed was investigated by Quanta 250 FEI FEG-SEM. The XRD patterns were recorded by Bruker D8 Advance diffractometer using CuKα1 radiation according to the Bragg-Brentano configuration. Their structural properties were also examined by Raman scattering measurements at room temperature using a Jobin Yvon/Horiba Labram spectrometer with an excitation Ar + laser at 488 nm. The cross-polarization (VH) configuration was used in order to reduce the signals from Si substrate.

Electrical property measurements

The resistivity of ZnO thin films was analyzed by using TLM. The metal electrodes, made by 50 nm-thick nickel (Ni) and 120 nm-thick gold (Au) layers, were prepared on the ZnO thin film surface using photolithography, e-beam evaporator and lift-off processes (Figure S1). The I-V curves between adjacent electrodes were measured using Karl Suss PM8 system. The total resistance (R Total ) obtained from the slope of I-V curves can be expressed as follows:

𝑅 = 2. 𝑅 + 2. 𝑅 + . 𝑑 (1)
where R m is the metallic electrode resistance, R c is the contact resistance at the electrode/thin film interface, and R s is the sheet resistance of the thin film; W is the width of the electrode, and d is the distance between electrodes (Figure S1). By recording the variation of R Total in function of d and fitting a linear curve, the sheet resistance (R s ) as well as the resistivity of the thin film can be deduced.

Piezoelectric and crystal polarity measurements

The piezoelectric properties of ZnO thin films were analyzed by PFM measurement using Bruker Dimension Icon atomic force microscope equipped with PtSi-NCH tips from Nanosensors, whose spring constant lies in the range of 43 -50 N/m. The metallic tip was kept in contact with the surface of the thin film during the scanning. Simultaneously, a 5 V AC was applied on the thin film to induce its deformation owing to its piezoelectric property. The amplitude and phase of the deformation were detected by the AFM tip, corresponding to the piezoelectric coefficient and polarity of the thin film, respectively. The frequency of AC voltage was set to 7 kHz to avoid the contact resonance frequency and high background signal contribution as reported in Ref. [START_REF] Jungk | Quantitative Analysis of Ferroelectric Domain Imaging with Piezoresponse Force Microscopy[END_REF]

Results and discussion

Structural characterization Figure 1 shows the top-view (left) and cross sectional-view (right) FESEM images of ZnO thin films before and after being annealed at different temperatures. From the cross sectional-view image, the as-grown ZnO thin film is found to be composed of stacked grains, as reported in our previous work [START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF] for that growth temperature of 500 °C. From the top-view images, it can be seen that the grains in all annealed ZnO thin films exhibit a more pronounced hexagonal shape evidencing the wurzite structure, which is the most stable phase of ZnO in these conditions. The average grain size constantly increases from 44 nm in as-grown thin film to 319 nm in 1000 °C-annealed thin film, as presented in

Figure S2

. This shows that grains coalesce and coarsen during an annealing-induced grain growth process to form large-sized grains when the annealing temperature is increased from 600 to 1000 °C.

Along with these phenomena, the root mean square roughness measured by AFM also increases from 13 nm in as-grown thin film to 76 nm in 1000 °C-annealed thin film. Grain coarsening and growth processes are the result of the migration of zinc and oxygen atoms to reduce the density of grain boundary with a higher energy. [START_REF] Thompson | Structure Evolution during Processing of Polycrystalline Films[END_REF] Simultaneously, the pore coalescence also occurred after thermal annealing, as indicated by the presence of some pores observed in top-view images and some dots with a lower contrast in cross-sectional view images, especially in 900 °C-and 1000 °C-annealed thin films. The similar pore formation in ZnO thin films after thermal annealing has been reported in Refs. [START_REF] Richardson | Thermally Induced Pore Formation in Epitaxial ZnO Films Grown from Low Temperature Aqueous Solution[END_REF][START_REF] Yi | Effect of Annealing on Performance of ZnO Thin Film Transistors[END_REF][START_REF] Nijikovsky | Microstructure of ZnO Films Synthesized on MgAl2O4 from Low-Temperature Aqueous Solution : Growth and Post-Annealing[END_REF] . These pores could affect the electrical properties of ZnO thin films. 3, 47.5, 56.6, 62.9, and 81.4°, respectively. Among them, the intensity of the 002 diffraction peak is much higher than the others, indicating that the ZnO thin films are strongly oriented with the polar c-axis normal to the Si substrate surface. As the annealing temperature is raised, the intensity of all ZnO diffraction peaks is significantly increased as revealed in Figure 2b, showing the improvement of the crystallinity of the thin films. Concomitantly, the full-width-at-half-maximum (FWHM) of the diffraction peaks is decreased and its position slightly varies. To get a more reliable and accurate information on the evolution of the grain orientation as a function of the annealing temperature, the texture coefficients C hkl of ZnO thin films for a given (hkl) plane were determined in percent unit by the following equation: 37

𝐶 (%) = , × Ʃ , (2)
where I hkl is the hkl peak intensity, I 0,hkl is the reference hkl peak intensity from the 00-036-1451 file of the ICDD, and N is the number of peaks considered (N = 6 here).

The FWHM and position of diffraction peaks are also further exploited to assess the homogeneous and inhomogeneous strains, as well as the crystallite size of ZnO thin films and their dependence on the annealing temperature. The homogeneous strain represents the mean residual strain in ZnO thin films at the end of the growth and annealing processes, which lead to the shift of the diffraction peak from its theoretical value in the XRD pattern. The homogeneous strain ε hs in ZnO thin film was calculated as follows:

𝜀 = (3)
where c o is the theoretical c-lattice parameter of wurtzite ZnO taken from the 00-036-1451 ICDD file and equal to 5.2066 Å, whereas c is the value deduced from the 002 diffraction peak position by using Bragg's law equation:

𝜆 = 2𝑑 sin 𝜃 ( 4
)
where λ is the wavelength of the Kα1(Cu) source equal to 0.15406 nm, θ is the Bragg angle of the 002 diffraction peak, d is the spacing between (002) planes and equal to half the c-lattice parameter in the wurtzite structure. The texture coefficient of the 002 diffraction peak is around 98.9 % in as-grown thin film and increases up to 99.9 % in 900 °C-annealed thin film, as shown in Figure 3a. Upon thermal annealing, grains coarsen and the ones with the polar c-axis normal to the Si substrate preferentially grow to form large-sized grains at the expense of the grains with other orientations that are progressively consumed.

This phenomenon is called abnormal grain growth in polycrystalline thin films [START_REF] Thompson | Structure Evolution during Processing of Polycrystalline Films[END_REF] and strongly depends on the processing parameters. [START_REF] Consonni | Thickness Effects on the Texture Development of Fluorine-Doped SnO2 Thin Films : The Role of Surface and Strain Energy[END_REF][START_REF] Consonni | Preferential Orientation of Fluorine-Doped SnO2 Thin Films : The Effects of Growth Temperature[END_REF] The grain boundary velocity is dependent upon the grain boundary mobility and the difference in the total free energy per unit volume representing the driving force for grain growth, both of them being strongly affected by the grain orientation. As the development of caxis oriented grains exhibiting the highest surface energy is promoted, [START_REF] Wander | Stability of Polar Oxide Surfaces[END_REF] an anisotropy in the grain boundary mobility is expected here as the dominant mechanism. A slight decrease in the 002 texture coefficient to 99.8% in the 1000 °C-annealed thin film could be due to the deterioration of the structural quality through the formation of large pores, as reported in Ref. [START_REF] Fang | Influence of Post-Annealing Treatment on the Structure Properties of ZnO Films[END_REF] . In Figure 3b, the compressive homogeneous strain of -0.15 % in as-grown thin film is measured. For that growth temperature of 500 °C, the homogeneous strain originates from intrinsic stress that is generated during the Volmer-Weber growth mode: [START_REF] Nix | Crystallite Coalescence: A Mechanism for Intrinsic Tensile Stresses in Thin Films[END_REF][START_REF] Floro | The Dynamic Competition between Stress Generation and Relaxation Mechanisms during Coalescence of Volmer -Weber Thin Films[END_REF] island nucleation and coalescence at the onset of growth generates high biaxial tensile stress that is relieved and switched to biaxial compressive stress by interface and inclined shears or impurity segregation at grain boundaries. The compressive homogeneous strain decreases down to -0.02 % in 700 °C-annealed thin film owing to grain coarsening and growth. The homogeneous strain, however, suddenly increases to -0.13% in 800 °Cannealed thin films. In that annealing temperature range, the homogeneous strain comes from extrinsic stress, basically thermal stress induced by the thermal expansion coefficient mismatch between ZnO thin film and Si substrate. This can cause a strong residual homogeneous strain after the sample was cooled down from high annealing temperature to room temperature. [START_REF] Fang | Influence of Post-Annealing Treatment on the Structure Properties of ZnO Films[END_REF][START_REF] Daniel | Effect of Annealing Temperature on the Structural and Optical Properties of ZnO Thin Films Prepared by RF Magnetron Sputtering[END_REF] Indeed, the thermal expansion coefficient of hexagonal ZnO in a-axis is around 3.2 -6.05 x 10 -6 °C-1 , while the thermal expansion coefficient of Si is 2.5 x 10 -6 °C-1 at room temperature. [START_REF] Jou | Substrate Dependent Internal Stress in Sputtered Zinc Oxide Thin Films[END_REF] As the thermal expansion coefficient of ZnO is higher as compared to the Si substrate, ZnO shrinks faster than Si during the cooling down process, resulting in a residual tensile strain in the a-axis of ZnO (ε xx ) at the interface with the Si substrate, which can be expressed as follows:

𝜀 = (𝛼 -𝛼 )∆𝑇 (5)
where α ZnO and α Si are the thermal expansion coefficients of ZnO and Si, respectively, ΔT is temperature difference. From this expression, the residual tensile strain in the a-axis of ZnO (ε xx ) is calculated to be in the range of 0.05 -0.28 % when the temperature is cooled down from 800 to 25 °C (room temperature). The tensile strain in the a-axis can be converted into the compressive strain in the c-axis (ε zz ) as follows: [START_REF] Wright | Elastic Properties of Zinc-Blende and Wurtzite AlN, GaN, and InN[END_REF][START_REF] Wagner | Properties of Strained Wurtzite GaN and AlN : Ab Initio Studies[END_REF] 𝜀 = -2 𝜀 (6)

where C 13 and C 33 are elastic constants of ZnO, which are equal to 104.2 and 213.8 GPa, respectively. [START_REF] Hong | Influence of Buffer Layer Thickness on the Structure and Optical Properties of ZnO Thin Films[END_REF] Thus, the compressive strain along the c-axis (ε zz ) is estimated to range from -0.27 to -0.05 %. This compressive strain value is comparable to the homogenous strain deduced from XRD results, indicating that the thermal expansion coefficient mismatch can cause notable residual strain after cooling down process, especially at higher annealing temperature. However, the residual strain decreases again when annealing temperature increased to 900 and 1000 °C. The appearance of pores in these thin films as seen in FESEM images in Figure 1 could be responsible for this strain relaxation.

Furthermore, the inhomogeneous strain caused by the microstructural defects in the ZnO structure along with the crystallite size were evaluated by using the Williamson-Hall method. The inhomogeneous strain and crystallite size both contribute to the broadening of the diffraction peaks.

By considering that the peaks are mainly described by a Gaussian function, their contribution to the FWHM denoted β of diffraction peaks can be expressed as follows: 47

𝛽 = 𝛽 + 𝛽 (7)
where β is the FWHM of 100, 002, 101, 102, 110, 103, 004 and 104 diffraction peaks, 𝛽 = 𝐶𝜀 tan 𝜃 represents the peak broadening contribution from the inhomogeneous strain and 𝛽 = represents the peak broadening contribution from the crystallite size, for which θ is the Bragg angle, C is a constant equal to 4, K is a constant depending on the shape of the particles and typically is equal to 0.9, 47 ε IS is the inhomogeneous strain, and L is the average crystallite size. Equation ( 7) can be rearranged as follows:

(𝛽 cos 𝜃) = (𝐶𝜀 sin 𝜃) + (8)
By plotting (βcosθ) 2 vs (sinθ) [START_REF] Kołodziejczak ; Radzimska | Zinc Oxide-From Synthesis to Application: A Review[END_REF] and fitting linear curves, we can deduce the inhomogeneous strain (ε IS ) from the slope and the crystallite size (L) from the intercept, respectively. Since the Williamson-Hall method includes several simplifying assumptions, it is a reliable method to qualitatively assess the evolution of the inhomogeneous strain and crystallite size with the annealing temperature.

The Williamson-Hall plots of ZnO thin films before and after being annealed are shown in Figure 3c. The inhomogeneous strain and average crystallite size are presented in Figures 3d and3e, respectively. In contrast to the evolution of the homogeneous strain, the inhomogeneous strain constantly decreases from 0.49 to 0.08 % (Figure 3d) when the annealing temperature increases to 1000 °C. This suggests that a number of extended defects like dislocations and stacking faults are removed at high annealing temperature. The average crystallite size also increases from 37 nm in asgrown thin film up to 138 nm in 1000 °C-annealed thin film (Figure 3e). Grain coarsening and growth processes thus occur both in the plane of the film and all along its thickness, as shown by top-view and cross-sectional view FESEM images in Figure 1, progressively replacing the stacked thin films to form an almost columnar thin film. 18 Figure 4 shows the Raman scattering spectra of ZnO thin films before and after annealing at different temperatures. The two dominant Raman lines located at 99 and 438 cm -1 are assigned to the E 2 low and E 2 high modes of the ZnO wurtzite structure, respectively. [START_REF] Cuscó | Temperature Dependence of Raman Scattering in ZnO[END_REF] Similarly to the diffraction peaks in the XRD patterns, their intensity significantly increases when the annealing temperature is increased, which again indicates that the crystallinity is largely improved with annealing temperature.

The E 2 high line also shifts to the bulk value at 438 cm -1 after annealing as shown in Figure S3, demonstrating again the strain relaxation. Furthermore, the relative intensity of the E 2 high line associated with the oxygen lattice with respect to the E 2 low line associated with the zinc lattice continuously increases as the annealing temperature is raised to 1000 °C. This reveals that the oxygen vacancies are progressively filled and vanish as the thermal annealing proceeds. Other Raman lines coming from the ZnO wurtzite structure such as E 2 high -E 2 low , A 1 (TO), E 1 (TO) and A 1 (LO) located at 333, 378, 410 and 574 cm -1 , respectively, are hardly seen in the Raman scattering spectrum of asgrown thin film, but gradually appear when annealing temperature is increased. [START_REF] Cuscó | Temperature Dependence of Raman Scattering in ZnO[END_REF] Notably, two large broad bands centered at around 1360 cm -1 and 1580-1600 cm -1 are observed in the Raman scattering spectra of the as-grown thin film, which can be assigned to C-C bonds. [START_REF] Ferrari | Interpretation of Raman Spectra of Disordered and Amorphous Carbon[END_REF] This indicates that there exist a contamination of carbon in ZnO thin films grown by PLI-MOCVD. The signals from C-C bonds are gradually decreased when annealing temperature is increased and completely vanish at 800 °C, indicating that those carbon contaminations are removed after annealing at sufficient high temperature.

Electrical characterization 3 Ω.cm in as-grown thin film to 50.7 Ω.cm in 800 °C-annealed thin film, but then decreases to 31.3 Ω.cm when annealing temperature was further increased to 900 °C. The electrical resistivity increases again to 81.0 Ω.cm at 1000 °C. The resistivity measurement was repeatedly performed within several months and showed similar results, indicating that the electrical resistivity of ZnO thin films is stable with time.

The variation of the electrical resistivity after annealing can be related to the variation of the charge carrier density and mobility, which are inversely proportional to the electrical resistivity. By increasing annealing temperature, some reports have shown that the resistivity of ZnO increases, [START_REF] Gritsenko | Effect of Thermal Annealing on Properties of Polycrystalline ZnO Thin Films[END_REF][START_REF] Asghar | Post-Annealing Modification in Structural Properties of ZnO Thin Films on p-Type Si Substrate Deposited by Evaporation[END_REF] while other reports have shown that it decreases. [START_REF] Ogata | Effects of Thermal Annealing of ZnO Layers Grown by MBE[END_REF][START_REF] Yi | Effect of Annealing on Performance of ZnO Thin Film Transistors[END_REF] The increase in the electrical resistivity was explained by i) the atoms relocated during thermal annealing in such a way that a reduced amount of defects that act as shallow donors in ZnO occurs, thus leading to a decrease in the charge carrier density, [START_REF] Asghar | Post-Annealing Modification in Structural Properties of ZnO Thin Films on p-Type Si Substrate Deposited by Evaporation[END_REF] or by ii) the oxygen adsorbed at the grain boundary after thermal annealing under oxygen atmosphere, creating depletion zones that decrease the charge carrier density and hamper the movement of charge carriers. [START_REF] Gritsenko | Effect of Thermal Annealing on Properties of Polycrystalline ZnO Thin Films[END_REF] In other reports, the decrease in the electrical resistivity was explained by the improvement of crystallinity, the decrease of defect density or the decrease of grain boundary density, which helps to reduce the scattering of charge carriers and thus enhance the charge carrier mobility. [START_REF] Ogata | Effects of Thermal Annealing of ZnO Layers Grown by MBE[END_REF][START_REF] Yi | Effect of Annealing on Performance of ZnO Thin Film Transistors[END_REF] In our case, the electrical resistivity increases when the annealing temperature is increased to 800 °C, which could be due to more oxygen adsorption creating the depletion zone at the grain boundary at higher annealing temperature. Besides, the decrease in the density of defects that can act as shallow donors such as carbon impurities [START_REF] Mohammadbeigi | Carbon Related Donor Bound Exciton Transitions in ZnO Nanowires[END_REF] could also be one of the reasons, which is suggested by the decrease of C-C bond signals in Raman spectra (Figure 4). These combined effects could cause a large decrease of charge carrier density, leading to an increase in the electrical resistivity when the annealing temperature is increased to 800 °C. Although the charge carrier mobility could also be increased thanks to the grain coarsening and growth processes and to crystallinity improvement with the increase of annealing temperature, it could not overcome the decrease in the charge carrier density.

From 800 to 900 °C, the electrical resistivity, however, decreases. Along with this observation, the Raman spectra of 800 °C-and 900 °C-annealed thin films are similar since the C-C bonds signals are completely removed in 800 °C-annealed thin film (Figure 4). Because the carbon impurity can act as a shallow donor, [START_REF] Mohammadbeigi | Carbon Related Donor Bound Exciton Transitions in ZnO Nanowires[END_REF] no reduction of C-C bond signals could state that there is no significant decrease of charge carrier density. On the other hand, the charge carrier mobility could continue increasing thanks to the grain coarsening and growth processes along with the crystallinity improvement, as shown in Figure 3e. Thus, the electrical resistivity decreases when annealing temperature increased from 800 to 900 °C. At 1000°C, the increase in the electrical resistivity is due to a significant impact of pores that interrupted the charge carrier motion. [START_REF] Yi | Effect of Annealing on Performance of ZnO Thin Film Transistors[END_REF] Piezoelectric characterization that the annealed thin films have a higher piezoelectric amplitude than the as-grown thin film. Also, the thin films annealed at higher temperatures have a larger number of area with higher piezoelectric amplitude than the thin films annealed at lower temperatures. Correlatively, their piezoelectric amplitude histograms shift and stretch to a higher value (Figure 7a-f). As the piezoelectric amplitude is proportional to the piezoelectric coefficient, this indicates an increase in the piezoelectric coefficient with the annealing temperature.

In the piezoelectric phase images as presented in Figure 6m-r, the area with the phase at around 70 ° (bright area) is attributed to the Zn-polar domains, while the area with the phase at around -90 ° (dark area) is assigned to the O-polar domains. [START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF] The phases corresponding to the two opposite piezoelectric polarities of ZnO were confirmed by performing PFM on ZnO monocrystals with an identified polarity. [START_REF] Bui | Morphology Transition of ZnO from Thin Film to Nanowires on Silicon and Its Correlated Enhanced Zinc Polarity Uniformity and Piezoelectric Responses[END_REF] Both dark and bright areas can be seen clearly in the phase response images of the thin films before and after annealing at 600 °C (Figure 6m-n). However, it becomes more difficult to distinguish between dark and bright areas in the phase response image of the 700 °C-annealed thin film (Figure 6o). At higher annealing temperatures (Figure 6p-r), the borders between dark and bright areas can be clearly seen again. Interestingly, the bright area becomes more dominant when the annealing temperature increases. Correspondingly, the negative peak around -90 ° and positive peaks around 70 ° in their phase histograms (Figure 7g-l) are assigned to Zn-and O-polar domains, respectively. Both negative and positive peaks can be seen in histograms of as-grown and 600 °Cannealed thin films (Figure 7g-h), which relate to Zn-and O-polar domains shown in their phase images (Figure 6m-n). Compared to the as-grown thin film, the thin film annealed at 600 °C has its phase histogram distribution reduced on the positive side and increased on the negative side. This reveals that the O-polar domains in the thin film are more developed after being annealed at 600 °C. In 700 °C-annealed thin film, a broad peak centered at nearly 0 ° is shown in its histogram (Figure 7i), corresponding to the dark-bright mixed areas that are observed in its phase image (Figure 6o). This could be explained by the coarsening process between grains with Zn-and O-polarity upon thermal annealing. As they are of opposite polarities, the coarsening and growth processes result in the formation of domains with compensated polarity, in which there is no piezoelectric signal. The phase response of non-piezoelectric area, which has a large area in 700 °C-annealed thin film, can be randomly distributed from negative to positive values, leading to a broad peak centered at 0 ° in its histogram, as seen in Figure 7i. In the thin films annealed at higher temperature, their phase responses are distributed more on the positive side (Figure 7j-l), exhibiting that the Zn-polar domains are predominant over the O-polar domains. This indicates that the grains with Zn-polarity prevailed over the grains with O-polarity after the grain coarsening and growth processes at high annealing temperature. Given that the surface energy of Zn-polar c-plane is significantly larger than the surface energy of O-polar c-plane, [START_REF] Tang | Activity of ZnO Polar Surfaces: An Insight from Surface Energies[END_REF] this supports again that the grain coarsening and growth processes at high temperatures is driven by an anisotropy in the grain boundary mobility rather than in the minimization of the free energy per unit volume governed by the surface energy contribution. [START_REF] Thompson | Structure Evolution during Processing of Polycrystalline Films[END_REF] Furthermore, the PFM data was treated to separate the piezoelectric responses of the Zn-and Opolar domains, and to remove noise from non-polarity areas. The data treatment is described as follows: data points with the phase response higher than 50 ° were assigned to Zn-polar domains and kept without modification; in contrast, data points with the phase response lower than -50 ° were assigned to O-polar domains and had their amplitude value multiplied by -1. Other data points were assigned to non-polarity areas and were removed. Figure 8 shows piezoelectric amplitude histograms after data treatment of ZnO thin films before and after annealing at different temperatures. These results show that the area of Zn-polar domains decreases in thin films annealed at 600 and 700 °C (Figure 8a-c), but increases again in thin films annealed at 800, 900, and 1000 °C (Figure 8d-f). The piezoelectric amplitude distribution of the Zn-polar domains also stretches to higher value in those high temperature annealed thin films (Figure 8d-f). On the other hand, the O-polar domains have its piezoelectric amplitude and area increased in the 600 °C-annealed thin film, but then its area decreases as being compensated by the Zn-polar 20 domains in 700 °C-annealed thin film and being consumed at higher annealing temperature (Figure 8d). The average piezoelectric amplitude response of the Zn-polar domains retains a value close to 9.5 pm in as-grown and annealed thin films up to 700 °C, and then increases significantly to 13.2 pm at 800 °C (Figure 9a). At the same time, the average piezoelectric amplitude response of the O-polar domains increases from a smaller value of 7.5 pm in the as-grown thin film to 11.3 pm in the thin film annealed at 600 °C. It then decreases slightly to 10.8 pm at 700 °C, and eventually increases slightly to 12 pm at 800 °C (Figure 9a). The average piezoelectric amplitude of both polar domains does not change at the annealing temperatures above 800 °C. It should be noted that the bottom electrode used in this report is highly doped Si substrate. We suspect that, after annealing at 1000 °C, the substrate was partially oxidized, and ZnO thin films started being peeled off the substrate. This could cause poor electrical contact between the ZnO thin film and bottom electrode, leading to an underestimation of the piezoelectric amplitude of this sample. The area percentage of the Zn-polar domains decreases from 71.7 % in the as-grown thin film to 30.5 % in the thin film annealed at 600 °C, but then increases constantly up to 93.2 % at 1000 °C (Figure 9b). These results again indicate that i) the piezoelectric amplitude in both polar domains increases with annealing temperature, and that ii) the Zn-polar domains becomes dominant after annealing at high temperature. The increase in the piezoelectric amplitude with annealing temperature could be due to the larger vertical alignment of the c-axis oriented domains in the ZnO thin films annealed at higher temperature, as suggested in the XRD patterns (Figure 3a). [START_REF] Lim | Crystal Growth and Piezoelectric Characterization of Mechanically Stable ZnO Nanostructure Arrays[END_REF] Another explanation could be related to the decrease in the charge carrier density after annealing at higher temperature as suggested in the TLM method (Figure 5), leading to a reduction of the screening effect. [START_REF] Lim | Crystal Growth and Piezoelectric Characterization of Mechanically Stable ZnO Nanostructure Arrays[END_REF] These effects appear to be more pronounced in the case of the Znpolar domains than in the case of the O-polar domains. The predominance of the O-and Zn-polar domains after annealing at low and high temperatures can be explained as follows. In the low annealing temperature range, a faster coarsening and growth rate of the O-polar grains as compared to the Zn-polar grains occurs, because the development of the O-polar domains with the lower surface energy [START_REF] Tang | Activity of ZnO Polar Surfaces: An Insight from Surface Energies[END_REF] is favorable for minimizing the free energy of the system and hence boosting the driving force for grain growth. In contrast, in the high annealing temperature range, a faster coarsening and growth rate of the Zn-polar grains as compared to the O-polar grains takes place owing to the predominant anisotropy in the grain boundary mobility. A comparison between the migration rates of grain boundary along Zn-and O-polarity was demonstrated in Ref. [START_REF] Lee | Effects of Polarity on Grain-Boundary Migration in ZnO[END_REF] . Their results showed that the migration rate is significantly faster in the Zn-polar domains than in the O-polar domains, which validates the present considerations. As a consequence, the development of polar domains upon thermal annealing in the ZnO stacked thin films is governed by the annealing temperature-dependent competitive process between the difference in the free energy per unit volume and the grain boundary mobility, both of them being strongly anisotropic. Overall, the prevalence of Zn-polar domains with a larger piezoelectric amplitude after high temperature annealing can lead to the polarity uniformity in ZnO stacked thin films grown by MOCVD, which significantly enhances its piezoelectric performance.

Conclusion

In summary, we have demonstrated the effect of high temperature annealing under oxygen atmosphere on the structural, electrical, and piezoelectric properties of ZnO stacked thin films grown by PLI-

22

MOCVD. By varying the annealing temperature from 600 to 1000 °C, the c-axis oriented grain coarsening and growth processes along with the highly improvement of crystallinity has been observed. A large amount of remaining carbon impurity from the growth process has also been removed after annealing. At the same time, the electrical resistivity has significantly been affected by the annealing process. Those reasons could explain the increase of piezoelectric amplitude with the annealing temperature. Interestingly, the PFM measurement also revealed that grains with O-and Znpolarity have become predominant after the coarsening and growth processes at low and high annealing temperature, respectively, which has led to polarity uniformity. An annealing temperaturedependent competitive process between the difference in the anisotropic free energy per unit volume and the anisotropic grain boundary mobility has been discussed to account for the dominant O-and Zn-polarity, respectively. The post-deposition thermal annealing of ZnO thin films not only enhances its c-axis orientation, but also its polarity uniformity, which are two critical conditions for enhancing the performances of piezoelectric applications. These results have shown a strong correlation of the structural, electrical and piezoelectric properties of ZnO thin films and point out that high temperature annealing under oxygen atmosphere can significantly improve its piezoelectric performance.

Figure 1 .

 1 Figure 1. Top-(left) and cross-sectional-view (right) FESEM images of ZnO thin film grown by PLI-MOCVD at 500 °C (a) and ZnO thin films after being annealed at (b) 600, (c) 700, (d) 800, (e) 900, and (f) 1000 °C.

Figure 2 .

 2 Figure 2. XRD patterns of ZnO thin films plotted in the range of 2θ (a) from 20 to 90° and (b) from 29to 37° after its intensity normalization using the baseline. The intensity is plotted in logarithm scale.

Figure 3 .

 3 Figure 3. (a) Texture coefficient of the 002 diffraction peak vs annealing temperature. (b) Homogeneous strain vs annealing temperature. (c) Williamson-Hall plot of ZnO thin films. (d) Inhomogeneous strain vs annealing temperature deduced from the Williamson-Hall plot. (e) Average crystallite size vs annealing temperature inferred from the Williamson-Hall plot.

Figure 4 .

 4 Figure 4. Raman scattering spectra of ZnO thin films. The intensity is plotted in logarithm scale.

Figure 5 .

 5 Figure 5. (a) Measured resistance vs distance in transmission line measurement and (b) deduced resistivity of ZnO thin films.

Figure 6 .

 6 Figure 6. Topography of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d) 800, (e) 900 °C, and (f) 1000 °C. Piezoelectric amplitude of ZnO thin films (g) before and after being annealed at (h) 600, (i) 700, (j) 800, (k) 900 °C, and (l) 1000 °C. Piezoelectric phase of ZnO thin films (m) before and after being annealed at (n) 600, (o) 700, (p) 800, (q) 900 °C, and (r) 1000 °C.

Figure 7 .

 7 Figure 7. Raw piezoelectric amplitude histograms of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d) 800, (e) 900 and (f) 1000 °C. Raw piezoelectric phase histograms of ZnO thin films (g) before and after being annealed at (h) 600, (i) 700, (j) 800, (k) 900 and (l) 1000 °C.

Figure 6 Figure 7 .

 67 Figure 6 shows the topography, piezoelectric amplitude and phase images from PFM measurements recorded over a 6 x 6 µm 2 area of the ZnO thin films before and after annealing at different temperatures. The corresponding histograms of piezoelectric amplitude and phase are shown in Figure 7. It can be seen that the piezoelectric response in ZnO thin films is inhomogeneous in both amplitude and phase images. The piezoelectric amplitude images as presented in Figure 6g-l show

Figure 8 .

 8 Figure 8. Piezoelectric amplitude histograms of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d) 800, (e) 900 and (f) 1000 °C after the data treatments.

Figure 9 .

 9 Figure 9. Average piezoelectric amplitude on (a) Zn-polar and O-polar area, and (b) percentage of Znpolar area of ZnO thin films.

Figure 1 .

 1 Figure 1. Top-(left) and cross-sectional-view (right) FESEM images of ZnO thin film grown by PLI-MOCVD at 500 °C (a) and ZnO thin films after being annealed at (b) 600, (c) 700, (d) 800, (e) 900, and (f) 1000 °C.

Figure 2 .

 2 Figure 2. XRD patterns of ZnO thin films plotted in the range of 2θ (a) from 20 to 90° and (b) from 29to 37° after its intensity normalization using the baseline. The intensity is plotted in logarithm scale.

Figure 3 .

 3 Figure 3. (a) Texture coefficient of the 002 diffraction peak vs annealing temperature. (b) Homogeneous strain vs annealing temperature. (c) Williamson-Hall plot of ZnO thin films. (d) Inhomogeneous strain vs annealing temperature deduced from the Williamson-Hall plot. (e) Average crystallite size vs annealing temperature inferred from the Williamson-Hall plot.

Figure 4 .

 4 Figure 4. Raman scattering spectra of ZnO thin films. The intensity is plotted in logarithm scale.

Figure 5 .

 5 Figure 5. (a) Measured resistance vs distance in transmission line measurement and (b) deduced resistivity of ZnO thin films.

Figure 6 .

 6 Figure 6. Topography of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d)800, (e) 900 °C, and (f) 1000 °C. Piezoelectric amplitude of ZnO thin films (g) before and after being annealed at (h) 600, (i) 700, (j) 800, (k) 900 °C, and (l) 1000 °C. Piezoelectric phase of ZnO thin films (m) before and after being annealed at (n) 600, (o) 700, (p) 800, (q) 900 °C, and (r) 1000 °C.

Figure 7 .

 7 Figure 7. Raw piezoelectric amplitude histograms of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d) 800, (e) 900 and (f) 1000 °C. Raw piezoelectric phase histograms ofZnO thin films (g) before and after being annealed at (h) 600, (i) 700, (j) 800, (k) 900 and (l) 1000 °C.

Figure 8 .

 8 Figure 8. Piezoelectric amplitude histograms of ZnO thin films (a) before and after being annealed at (b) 600, (c) 700, (d) 800, (e) 900 and (f) 1000 °C after the data treatments.

Figure 9 .Figure S1 .

 9S1 Figure 9. Average piezoelectric amplitude on (a) Zn-polar and O-polar area, and (b) percentage of Znpolar area of ZnO thin films.

Figure S2 .

 S2 Figure S2. Average grain size of ZnO before and after being annealed at 600 -1000 °C measured from top-view SEM images.

Figure S3 .

 S3 Figure S3. The E 2 high peak position in Raman spectra of ZnO before and after being annealed at 600 -1000 °C.

Acknowledgements

This work was supported by the French National Research Agency in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02) through the project CDP NEED. Q.C.B. held a doctoral fellowship from the project CDP NEED. The authors further acknowledge the support from the CNRS Renatech Network through the "Plateforme Technologique Amont" in a cleanroom environment. This research has also benefited from some of the characterization equipments of the Grenoble INP-CMTC platform. The authors are further grateful to Odette Chaix-Pluchery, LMGP, Grenoble, France, for her assistance in the Raman spectroscopy. This work benefited from the facilities and expertise of the OPE)N(RA characterization platform of FMNT (FR 2542, fmnt.fr) supported by CNRS, Grenoble INP and UGA.