Dr Wei Zhang

Dr Hanna Enriquez

Dr Yongfeng F Tong

Dr Andrew J Mayne

Dr Azzedine Bendounan

Dr Yannick J Dappe

Abdelkader Kara

Prof Gérald H Dujardin

Hamid Oughaddou
email: hamid.oughaddou@u-psud.fr

Prof H Oughaddou

Prof A Kara

Phosphorus Pentamers: Floating Nanoflowers form a 2D Network

Keywords: phosphorene, 2D materials, low temperature scanning tunneling microscopy, X-ray photoelectron spectroscopy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The great success of graphene has inspired an explosive increase in the exploration of novel two-dimensional (2D) materials. [1][2][3][4] Examples include silicene, [5][6][7][8] transition-metal dichalcogenides, [9][10][11] germanene, [12][13][14] borophene [15,16] , stanene [17] , and phosphorene. [18,19] Among them, phosphorene, made of one or few atomic thick phosphorus layers, has received particular attention recently due to its exotic properties of high carrier mobilities, [19,20] inherent tunable direct band gap (from bulk ~0.3 eV to monolayer ~1.8 eV), [21][22][23][24][25][26] high on/off ratios, [19,26] and high in-plane anisotropy. [27,29] Black phosphorus as the most stable allotrope of phosphorus, [30] is composed of layers stacked together by weak van der Waals forces. [31] The exfoliation method from a 3D crystal can be used to produce a few layers of black phosphorene. [18,19,23] While each phosphorus atom is covalently bonded to three adjacent atoms, these phosphorene layers prefer to form a corrugated sheet structure unlike the planar graphene. [32] The molecular beam epitaxy (MBE) technique has proven to be advantageous in producing phosphorene on metal substrates. Blue phosphorene can be grown on Au(111), [33][34][35][36][37][38] while only phosphorus clusters are reported on Cu(111). [39,40] Thus, the effect of the substrate is an important factor determining the nature of the phosphorene growth. A moderate interaction strength has been found for a P layer adsorbed on the Ag(111) surface. [41,42] However, until now the actual production of epitaxial phosphorus sheets on Ag (111) has not yet been reported. Theoretical calculations predict that phosphorene can form many polymorphic structures that do not necessarily have a hexagonal structure. [43] Indeed, earlier calculations indicate that pentameric structures are stable. [44] In this study, for the first time, we report the epitaxial growth of single layer phosphorus on Ag(111) surface under ultrahigh-vacuum (UHV) conditions. Our results demonstrate that the phosphorus is organized in the form of pentamers presenting a band gap of 1.2 eV. This is the first observation of such a structure made by phosphorus. The detailed morphology is characterized by in situ scanning tunneling microscopy (STM) operated at 78 K. The chemical composition is determined by in situ Auger electron spectroscopy (AES) and X-ray photoelectron spectroscopy (XPS). The electronic structure of the phosphorus pentamers on Ag(111) is probed by scanning tunneling spectroscopy (STS). The adsorption and diffusion of P-pentamers on Ag(111) were calculated using DFT with self-consistent van der Waals interactions.

Results and Discussion

In situ AES was carried out during the deposition process. Typical AES spectra of Ag(111) surface before and after 0.5 monolayer (ML) of phosphorus deposition are presented in Figure 1. The characteristic spectrum of Ag MNN acquired on clean Ag(111) presents sharp peaks at 351 and 356 eV and broadened secondary peaks at 260 and 302 eV. [START_REF] Childs | Handbook of Auger electron spectroscopy: a book of reference data for identification and interpretation in Auger electron spectroscopy[END_REF] After deposition of phosphorus, the P LMM peak is clearly observed at the expected energy position of 120 eV.

Moreover, none of the features show any characteristic shifts or satellite peaks that would appear if surface alloying or compound formation had occurred. Ag(111). Spectra collected from clean (black line) and after 0.5 ML of phosphorus deposition (red line). The LEED pattern recorded at room temperature (RT) did not show any extra diffraction spots related to the phosphorus deposition (Figure S1a-b, Supporting Information). This suggests that the phosphorus layer is disordered. STM experiments performed at RT to investigate the structure of the phosphorus layer reveal a flat 2D phosphorus layer in the large-scale STM images (Figure S1c, Supporting Information). The smooth round edges to the layers have no clear crystallographic orientation with respect to the Ag(111) substrate reminiscent of a liquid-like growth phase. This liquid behavior has been observed in the growth of Perlyene on Ag(110). [START_REF] Kalashnyk | [END_REF] Furthermore, a close-up STM image of the phosphorus layer clearly shows no order even at the nanometer scale (Figure S1d, Supporting Information). The RT-STM experiments confirm the initial analysis of the LEED pattern.

Figure 1. AES spectra of

An investigation of the growth of the Phosphorus layer was made at low temperature. After P deposition, the sample was cooled down to 78 K and analyzed in the STM. A typical largearea STM image is shown in Figure 2a. In comparison with the image recorded at RT, the phosphorus appears to have undergone a phase change from a liquid-like layer to a crystalline structure. The phosphorus deposition results in a single layer of pentamers that are one atomic layer thick and forms different rotational domains covering the terraces (highlighted with different colors). Comparison of the orientation of the pentamer domains (Figure 2a) relative to the metal substrate (inset in Figure 2a) yields an angle of 19°. All distinct domains have sizes in the 50-300 nm range, exhibit a striped structure rotated by 19° relative to the three main crystallographic directions of the Ag(111) substrate, leading to six different domains.

At RT, we observe that thermal diffusion of the pentamers is facile, owing to a smooth potential energy surface with shallow binding sites. As the temperature is lowered, selfassembly becomes energetically favorable so that an ordered layer is observed in the STM images. This competition between thermal diffusion and ordering occurs at the molecular level, for example in benzene on silicon. [47] All these observations point to a weakly corrugated potential energy surface between the phosphorus layer and the Ag(111) substrate.

The atomically resolved STM image (Figure 2b) shows that the phosphorus sheet is composed of a series of rows. Careful inspection of the STM image reveals that each small structure is composed of 5 atoms in the form of a pentamer. The average distance between the nearestneighbor pentamers is around 0.76 nm from the line profile in Figure 2c (corresponding to the blue line in Figure 2b). This value is √7 × the unit cell of Ag(111) (√7 × 0.289 nm = 0.765 nm). In addition, Figure 2b shows that pentamers are close-packed into rows. The pentamers within the same row have the same orientation while a neighboring row has a different orientation giving a periodic pattern to the stripes (see Figure 3 for a detailed description).

The line profile measured across these stripes (red line in Figure 2b) gives a periodicity of 4.53 nm which is very close to 6√7 times the unit cell of Ag(111) (6√7 × 0.289 = 4.6 nm).

Hence The phosphorus pentamer structure presents a (√7 × 6√7)R19° unit cell with respect to Ag(111). The unit cell is indicated by black parallelogram.

The periodicity relationship is confirmed in the fast Fourier transform (FFT) of the STM images. Figure 2e presents the FFT pattern of the clean Ag(111) surface, where the six symmetric bright spots related to the pristine Ag(111) substrate are highlighted in white circles. Figure 2f shows the FFT pattern of Figure 2b and gives only the periodicity of the pentamers structure as only pentamers are shown in the STM image. The spots can be assigned to the phosphorus superstructure; the green and yellow circles highlight diffraction spots corresponding to √7 and 6√7, respectively. The FFT of the two STM images gives independent confirmation as the rotation of 19° is clearly visible in reciprocal space (see Figures 2e and2f) in agreement with the analysis of the STM images (see Figure 2a). The high-resolution STM image in Figure 3a reveals the finer details of the pentamer structure of phosphorus on Ag(111), where different orientations of pentamers are observed.

One of the unit cells is indicated by the black parallelogram. Along the direction a, the pentamers are similar and have the same relative alignment. However, along the direction b, each pentamer is rotated by 36° relative to its neighbor as highlighted by the red and blue pentagons in Figure 3a. The line profile (Figure 3b) measured along the lines A and B in Figure 3a gives a buckling of the pentamers in the range 0.09 -0.23 Å (very close to the calculated buckling of about 0.3 Å, see Supplemental Material), even the STM technique is sensitive to both topology and electron density. The average distance measured between the P atoms within each pentamer is around 0.27 nm. This value is larger than the 0.224 nm reported for the phosphorene sheet. [20] Furthermore, the five protrusions of each pentamer show an asymmetry in height; three bright and two dark indicating a clear buckling of the P atoms in each pentamer. STS measurements were performed to investigate the electronic structure of this new structure.

The differential conductivity (dI/dV) spectrum recorded on the phosphorus pentamers is shown in Figure 3c using the lock-in technique collected at random positions on top of pentamers. The representative dI/dV spectrum shows a semiconductor character of the phosphorus structure. In the spectrum the peaks are located at -2.0 eV and -1.0 eV below the Fermi level. The dI/dV spectrum reveals a band gap of ~1.2 eV for the new phase of phosphorus pentamers. The finite intensity observed between -0.5 and -1 V is due to the contribution of the silver substrate beneath pentamers to the total density of states probed by the STM tip. Because the metal states are closer to the Fermi level, they leak into the band gap of the phosphorus pentamer layer. This has been observed in thin insulating films on metals. [48,49] As already reported for black phosphorene the band gap is related to the thickness of the layer. [26] The band gap value of the pentamer structure is close to the one obtained experimentally for blue phosphorene, [33,35,37] but smaller than the one obtained for an exfoliated black phosphorene monolayer (2.0 eV). [25] However, like black phosphorene, the phosphorus pentamers show a high in-plane anisotropy as the five protrusions within the pentamer show a height asymmetry. The temporary formation of pentamer vacancies is observed in the STM images. Figure 4a-d depicts a series of STM images taken on the same area of the phosphorus layer. Some vacancies having the form of a single missing phosphorus pentamer remain stationary during the tip scan (highlighted by the red circle). The fact that vacancy has the size and shape of a pentamer indicates that the pentamer is the fundamental building block of this structure.

Vacancies are observed to undergo diffusion within the layer (marked by the blue dashed circle). This diffusion process points to two interesting properties of the pentamers: the pentamer entity is stable on its own, and the diffusion is possible even at 78 K. The presence of vacancies and their diffusion even at low temperature shows that the pentamer-silver potential energy surface is smooth with very weak corrugation because vacancies can appear and disappear within the layer. At RT, thermal diffusion dominates the STM images because not only the structure is disordered and liquid-like but the atomic resolution of the pentamers cannot be obtained. The long-range order is determined by the interaction between pentamers.

There is no evidence of individual pentamers on top of the P-layer, indicating that diffusion is confined to within the layer. This has been observed in the case of atomic vacancies, [50] and molecular vacancies in a molecular layer. [51,52] This is in stark contrast to metal surfaces where diffusion proceeds via adatoms on the surface. [53] Further analysis of the nature of the chemical interaction between the phosphorus and Ag(111) surface, was obtained from XPS measurements. All spectra are fitted with a Doniach-Sunjic line shape. [54] The experimental data are displayed with black dots, while the fitted curves are in black solid lines. Figure 5a-b shows the characteristic XPS Ag 3d core-level spectra acquired on the clean Ag(111) surface and after formation of phosphorus pentamers. For the clean substrate, the characteristic signals of Ag 3d 3/2 and 3d 5/2 are located at binding energies of 374.2 eV and 368.2 eV, respectively. The best fit was obtained with only one component (red curve) with a 140 meV Gaussian profile and a 280 meV Lorentzian profile, while the spin-orbit splitting is 6 eV.

After phosphorus deposition, the shapes and energy positions of these two representative peaks remain unchanged. We also observed no splitting of these Ag 3d peaks, indicating the weak chemical interaction between the phosphorus and the underlying Ag substrate. The XPS P 2p core-level spectra obtained from phosphorene recorded at normal and grazing emission are respectively shown in Figure 5c, 5d. The best fit was obtained with a 200 meV Gaussian profile and a 80 meV Lorentzian profile, while the spin-orbit splitting is 0.865 eV. The P 2p core-level spectrum in Figure 5c is composed of two components S1 and S2 located respectively at 129.72 and 129.04 eV. These two components are assigned to two types of phosphorus atoms with different chemical environments. This can be explained by the buckled structure of phosphorus pentamers observed by STM. One component corresponds to P atoms in contact with the substrate and the second corresponds to P atoms in contact only with P atoms. The buckling in the pentamer can be explained by the strong physisorption (van der Waals interaction) of the pentamers on Ag(111), where the calculated binding energy was found to be 0.78 eV/atom (see Supplemental Material).

With increasing emission angle, the XPS surface sensitivity is enhanced to selectively detect the surface and sub-surface contributions. The relative intensities of the surface component peak will be increased. Thus, from the Figure 5d we can conclude that the S1 component is assigned to the phosphorus atoms buckled above the plane and the S2 component is assigned to the phosphorus atoms buckled below the plane. These results confirm the buckled structure observed in the STM images. In addition, from the XPS data we can estimate the ratio of the area of the two components S1 and S2 of the P 2p core level. In normal emission, the ratio S1/S2 is equal to 0.64 while for the grazing emission it is equal to 1.36. These values are consistent with 3 P atoms up and 2 P down in each pentamer, in agreement with the STM images showing three bright and two dark protrusions. The P 2p spectra taken at (c) normal emission and (d) grazing emission (70° off-normal), showing two asymmetric components (S1) and (S2). For all spectra, the dots correspond to data and black line overlapping the data corresponds to the best fit.

Conclusion

We have synthesized successfully a high-quality single layer of phosphorus on the Ag(111) substrate using MBE; this new 2D phase is composed of phosphorus pentamers. STM measurements demonstrated a characteristic array of buckled pentamer structure of phosphorus atoms on Ag(111) with anisotropic corrugation. STS investigations determined the semiconductor character of phosphorus pentamers with a band gap of 1.20 eV close to that of silicon, making this material interesting for electronic and opto-electronic devices. XPS measurements and AES analysis reveal a weak chemical interaction between phosphorus and Ag(111). However, the buckling in the pentamer points to a strong physisorption that is confirmed by DFT+vdW calculations. Moreover, XPS further confirmed the buckled conformation of phosphorene. This work reveals a new 2D polymorphic structure of phosphorus, opening new avenues for this novel 2D material.

Experimental Section

The experiments were performed in a commercial ultra-high vacuum (UHV) system with a base pressure of below 1.0 × 10 -10 mbar for surface preparation and characterization. The single-crystal Ag(111) substrate (99.999%) was cleaned by several cycles of Ar + sputtering (650 eV, 1 × 10 -5 mbar) followed by annealing at the temperature of 500°C for 40 minutes.

During growth process, a Knudsen source loaded with black phosphorus was used to deposit P atoms on Ag(111) surface, maintained the substrate at 150°C by radiative heating from a filament located behind the sample. STM measurements were conducted in an Omicron Scanning Tunneling Microscope, operating at 78 K or 300 K. All STM experiments were acquired in constant current mode, using electrochemically etched W tips for surface characterization. STS experiments were acquired using the lock-in technique, with an applied modulation voltage of 10 mV at a frequency of 5127.7 Hz. The photoemission measurements were performed at the TEMPO beam-line of Synchrotron SOLEIL-France using a Scienta SES-2002 electron spectrometer at RT with energy resolution better than 50 meV. [START_REF] Polack | AIP Conf. Proc[END_REF] The calibration of the XPS measurements was made with the Ag 3d core level corresponding to the clean Ag substrate. The photon energy of XPS experiments were recorded at 470 eV for Ag 3d and 280 eV for P 2p, respectively.

Computational Section

To study the adsorption of P-pentamers on Ag(111) we used the Vienna Ab initio Simulation Package (VASP) version 5.4.4. [START_REF] Kresse | [END_REF][57][58] , which contains the projector augmented wave (PAW) method [59][60] . To model the exchange-correlation interaction we employed the optB88-vdW functional 61 . This particular functional includes a nonlocal implementation of the correlation energy term with the aim of seamlessly and self-consistently incorporating the long-range van der Waals (vdW) interaction. In the past, several studies demonstrated the importance of vdW interactions (including the use of optB88) in a myriad of systems [62][63][64][65] , including organic molecules on transition metal substrates. To model the gold surface, we created a 4-layer slab with each layer containing 12 atoms. VASP uses periodic boundary conditions so we placed at least 20 Å of vacuum in between slabs. A lattice constant of 4.147 Å was taken. To achieve structural relaxation, we utilized the conjugate gradient method. The force criterion was set to 0.02 eV/Å, the plane wave energy cutoff was chosen at 400 eV, and the Brillouin zone was sampled using a 4x4x1 k-points mesh. To simulate the diffusion, we used the Nudged Elastic Band method, with Climbing Image (NEB+CI), using 5 intermediary images 66 .

Figure 2 .

 2 Figure 2. Formation of 2D phosphorus sheets on Ag(111). a) Topographic STM overview map revealing that the phosphorus covers the Ag(111) terrace and forms different domains highlighted by different colors (U = -200 mV, I = 0.80 nA). The inset shows an atomically resolved STM image of clean Ag(111) surface, the black arrows indicate [1-10] crystal directions of the Ag lattice. b) Zoom of a showing a striped-phase atomic-scale structure of the phosphorus adlayer (U = -500 mV, I = 1.00 nA). The unit cell is indicated in black. c) Line profile taken along the black line in (a), revealing the distance between nearest-neighbor pentamers of 0.76 nm. d) Line profile measured across the stripes (green line in image (b)), revealing the periodicity of about 4.53 nm. e) Fast Fourier transform (FFT) of clean Ag(111). The white circles represent the (1 × 1) Ag(111) lattice. f) FFT of (b). The green and yellow circles correspond to × √7 in one direction and × 6√7 in the other direction relative to Ag. The rotation between the unit cell of silver and the pentamer superstructure (19°) is indicated. The unit cell in the reciprocal space is also indicated in white.

Figure 3 .

 3 Figure 3. Structure of phosphorus pentamers. a) Atomically resolved filled state STM imageshowing buckled pentagonal lattice structure (U = -100 mV, I = 3.0 nA). b) Line profile along black lines in (a), revealing a corrugation in the phosphorus pentamers. c) STS dI/dV spectrum acquired on phosphorene sheet, reveals a semiconductor characteristic with a band gap of 1.2 eV. The dI/dV curve for the clean Ag(111) surface is included for comparison.

Figure 4 .

 4 Figure 4. a-d) Consecutive STM images of the same area of the single layer of phosphorus on Ag(111) surface. One vacancy (red circle) remains at the same position in all images, others (blue circle) diffuse on the surface (U = -700 mV, I = 1.10 nA).

Figure 5 .

 5 Figure 5. XPS analysis for the single layer of phosphorus on Ag(111) substrate. a) The Ag 3d core-level spectra of clean Ag(111) surface. b) After formation of phosphorus pentamers on Ag(111). In both spectra, the best fit is obtained with only one component (red curves). c,d)

Acknowledgements W.Z. would like to thank the China Scholarship Council (CSC) for the PhD financial support (scholarship). The computational work of AK is supported by the U.S. Department of Energy Basic Energy Science under Contract No DE-FG02-11ER16243. Resources of the National Energy Research Scientific Computing Center (NERSC) were used for the computations.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.


B Oger

P Vismara

B Tisseyre

Combining target sampling with within field route-optimization to optimise on field yield estimation in viticulture

Keywords: Yield estimation, Sampling, NDVI, Constraint programming, Simulation, Spatial data, Viticulture

This paper describes a new approach for yield sampling in viticulture. It combines approaches based on auxiliary information and path optimization to offer more consistent sampling strategies, integrating statistical approaches with computer methods. To achieve this, groups of potential sampling points, comparable according to their auxiliary data values are created. Then, an optimal path is constituted that passes through one point of each group of potential sampling points and minimizes the route distance. This part is performed using constraint programming, a programming paradigm offering tools to deal efficiently with combinatorial problems. The paper presents the formalization of the problem, as well as the tests performed on nine real fields were high resolution NDVI data and medium resolution yield data were available. In addition, tests on simulated data were performed to examine the sensitivity of the approach to field data characteristics such as the correlation between auxiliary data and yield, the spatial auto-correlation of the data among others. The approach does not alter much the results when compared to conventional approaches but greatly reduces sampling time. Results show that, for a given amount of time, combining model sampling and path optimization can give estimation error up to 30% lower for a given amount of time compared to previous methods.

Introduction

In order to optimize harvest organization and quality management, the wine industry needs to know the yield of each vine field. Ideally, yield has to be estimated a few days before harvest with a relative error of less than 10% [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF]. Although models have 1 3 been developed to forecast the yield at the regional level [START_REF] Cristofolini | Concentration of airborne pollen of Vitisvinifera L. and yield forecast: A case study at S.Michele all'Adige[END_REF], their results were not precise enough to manage logistic issues in relation to harvest operations at the farm or at the winery level. Therefore, precise estimation of vine field yield always requires fruit sampling and counting [START_REF] Clingeleffer | Crop development, crop estimation and crop control to secure quality and production of major wine grape varieties. A national approach: Final report to grape and wine research & development corporation[END_REF]. This estimation must be carried out quickly (few minutes per field) at a time when the workload at harvest or for the preparation of the harvest is critical. Practical constraints, like the time available to visit all the fields before harvest, limit the number of sampled sites per field. Therefore, yield estimation is based on a low number of sites sampled (4 to 5 sites per field) where yield components (number of clusters, number of berries per cluster, mean berry weight) are manually measured by a practitioner. Due to these practical constraints and the high within-field variability of grape yield usually observed [START_REF] Taylor | A comparison of the spatial variability of vineyard yield in European and Australian production systems[END_REF], the small number of observations results in high errors in yield estimation (generally around 20 to 30%).

Recent works [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF][START_REF] Uribeetxebarria | Assessing ranked set sampling and ancillary data to improve fruit load estimates in peach orchards[END_REF][START_REF] Arnó | Comparing efficiency of different sampling schemes to estimate yield and quality parameters in fruit orchards[END_REF]) have shown the interest of integrating auxiliary data to improve sampling strategies and yield estimation for perennial crops. Among possible auxiliary information, vegetation index derived from multispectral airborne images is of great interest since they permit to characterise the spatial variability of several fields; in one acquisition, with a high spatial resolution (< 1 m) and at an optimal date. In viticulture, [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] shows the potential of NDVI to drive target sampling of the main grape yield components (bunch number, berry weight, etc.) to improve yield estimation. Although spatial patterns of yield and vegetative expression, estimated by vegetation indices (i.e. NDVI) may not match systematically in all the situations [START_REF] Bramley | Spatio-temporal variability in vine vigour and yield in a Marlborough Sauvignon Blanc vineyard[END_REF]. [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] demonstrated, in a dry vineyard of southern France, the value of using NDVI information to determine relevant within field sampling sites selection based on the distribution of NDVI values.

Although interesting, the methodology proposed by [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] presents a significant drawback. Indeed, it does not take into consideration the relative position of the sites to be sampled and the fact that vine fields are structured in rows. This peculiarity implies that rows cannot be crossed, leading to sampling plans optimized in terms of prediction but potentially unrealistic in terms sampling routes and resulting travelled distance (and time) for the operator. This paper proposes a new approach to optimally design withinfield sampling routes which takes into account the spatial organisation of the crop (rows) and spatial location of sampling sites. The originality of this approach, called constrained sampling, is to combine statistical and computer methods. It can be decomposed into two steps. In the first step, potential sampling sites are sorted into different groups according to their auxiliary data value in a similar way to traditional targeted sampling. The second step finds an optimal route that passes through one sampling site from each group. A constraint programming solver is used to build an optimal route in terms of travelled distance. This kind of solver has already been used in precision viticulture to solve the differential harvest problem [START_REF] Briot | Integration of Operational Constraints to Optimize Differential Harvest in Viticulture[END_REF].

The objective of this work is therefore to propose the resolution of a sampling problem by combining two methodologies from two very different scientific fields: a stochastic approach aimed at considering a large number of representative candidates (the sampling sites) with a constraint programming approach whose objective is to search for the optimal solution among a number of identified candidates. This combination is an interesting scientific question since it is necessary to simplify the exploration of possible candidates in order to be able to calculate an optimal route, while preserving the abundance of information provided by the variability provided high spatial resolution data. The scientific hypothesis of this work is therefore to test and validate the contribution of a combination of methods for spatial sampling in agriculture. In particular, the aim is to verify to what extent the introduction of an optimisation approach decreases the quality of the estimates by taking conventional sampling approaches as a reference. It will also examine whether gains in sampling time compensate for decreases in the quality of the estimates. Tests are performed either on experimental data from France or on simulated plots to better characterize its performance under different spatial structures.

Materials and methods

Sampling sites and selection principles

Overview

The purpose of constrained sampling (CS) is to select N sampling sites constituting a sampling route in the vine field. Accounting for classical sampling practices in viticulture, N will vary between 5 and 10. It is assumed that there is a finite number of sites on the plot where sampling can be carried out, these sites are called potential sampling sites. For instance, considering the plant as a potential sampling site, a one-hectare vineyard plot planted at 4000 vines/ha consists in 4000 potential sampling sites. For each potential sampling site (PSS), the method assumes that: (i) the coordinate of the PSS, (ii) the row that the PSS belongs to and (iii) the corresponding auxiliary data value are known.

The method requires computation of a distance matrix. This matrix gives the distance between each couple of PSSs. Distance must correspond to the shortest walking distance between the two PSSs. It must take into account the structure of the vineyard. Each PSS, located on one vine rows, can be accessed through two inter-rows (Fig. 1). As a result, each site appears twice in the distance matrix. If the sites are in the same inter-row, it corresponds to the classical Euclidian distance. If they belong to different inter-rows, the distance is computed considering that the practitioner has to get out the row, to reach the desired rows passing by all ends of intermediate rows and finally to reach the targeted sampling site (Fig. 1). As rows have two extremities, two different distances can be computed and the shortest one is kept.

Representative sampling based on auxiliary information

The same approach as [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] was used to consider auxiliary data. It is summarised hereafter. The sampling approach aims at calibrating a linear regression which relates the yield (sampled) to auxiliary data. [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] showed that yield components, especially berry weight, are linearly related to auxiliary data (NDVI) in non-irrigated vineyard of south of France. The sampling approach aims at selecting sampling sites (SSs) representatively to build this linear model. This model is then used to estimate yield using all available high-resolution auxiliary data.

The approach proposed in this paper relies on the following principle. PSSs are split into N homogeneous groups (N corresponding to the number of samples) according to their auxiliary data values. Once groups are formed, one PSS from each group is selected in order to optimize the length of the route connecting all selected PSSs. This ensures the PSSs are representative of the field as selected points are spread across the auxiliary data distribution. Quantiles were used to make groups. Related work concluded that this approach was adapted to yield estimation [START_REF] Oger | Combining target sampling with route-optimization to optimise yield estimation in viticulture[END_REF]. Quantiles method guarantees that the groups have the same number of PSS, K being the number of PSSs and N the number of samples required, then each group has K N elements (Fig. 2).

Route optimisation

The second step of the approach consists in selecting N sampling sites (one sampling site per group). These N sampling sites (SSs) must be all different and have to form the shortest possible sampling circuit. There are many possible choices to select these SSs and many ways to order them to form a circuit. It is therefore a highly combinatorial optimization problem. Constraint Programming is one of the programming paradigms able to deal with such problems. It aims at solving a problem expressed as a set of variables and a set of constraints on these variables. Such a problem is called a Constraint Satisfaction Problem (CSP). A Constraint Solver is used to find a solution to the problem that satisfies all the constraints. The efficiency of these solvers relies on the implementation of many methods such as filtering, which allows a quick detection of combinations of values that do not lead to an optimal solution. The interest of constraint solvers lies in their ability to address many types of constraints. Without going into detail, let S = {1, … , K} be the set of PSSs and G i i∈(1,…,N) the set of quantile groups covering S, formed in the previous step from the auxiliary data. Since all these groups are disjoint subsets, G i i∈(1,…,N) is a partition of S . P i is defined as the selected site for group G i . Hence, set S Selected (the set sampling sites) will be equal to

P i i∈(1,…,N) .
The first constraint imposes that all P i must be different; because all quantile groups are disjoint subsets, it is immediately satisfied in the case of the presented approach. P 0 represents the point of departure and arrival, it is a fixed parameter representing the initial position of the practitioner. The length of the optimum route passing through all the P i∈(0,…,N) must be minimum. This is a particular case of vehicle routing problem (VRP) where the goal is not to find a Hamiltonian tour (visiting once every site) but a tour covering only a subset of sites. Recent work about the WeightedSubCircuit constraint [START_REF] Vismara | A circuit constraint for multiple tours problems[END_REF] has proposed a filtering algorithm that is well adapted to address this type of situation. All these constraints and variables constitute the constraint satisfaction problem. An instance of this problem is built from each dataset and solved with the solver in order to get an ordered set of sites that form a sampling circuit. The program returns the list of sampling sites, the order in which they are visited and the associated walking distance.

Yield estimation

The aim is to estimate Y , the average yield of the field. For each selected sampling site s ( s ∈ S Selected ), GW(s) is the grape weight per vine value. A linear model linking the aux- iliary data ( AD ) to GW is built from these sites (Eq. 1):

For a given site s , ĜW(s) represents an estimate of GW(s) . The parameters a and b are obtained from a linear regression on the N sites selected by the sampling method.

With S = {1, … , K} being the full set of PSSs available, ŶCS , the yield estimate with Constrained Sampling (CS) approach, can be computed from the model using all these PSSs (Eq. 2):

Estimation error

Regardless of the method used, the estimation error is a deviation from the actual yield value ( Y ) and its estimation ( Ŷ ), expressed as a percentage (Eq. 3).

(1) ĜW

(s) = a × AD(s) + b (2) ŶCS = mean s∈S ( ĜW(s))

Reference methods

The method is compared to two references:

A conventional random sampling (RS) approach where the N sampling sites are randomly selected among all the PSSs. ŶRS , the yield estimated with RS, is therefore directly calculated from the mean of observed GW values. RS represents what is generally done in practice in terms of yield estimation.

The second approach is model sampling (MS) whose method principles have been described by [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF]. SSs are chosen according to NDVI values. One site is randomly selected for each of the N NDVI quantiles. MS uses a model based on the NDVI/yield relationship, as described in Eq. 1. Unlike constrained sampling method, the selection of SS does not consider their position on the plot.

To compare the length of the routes between the different methods, the optimal route between the selected sites was computed for both reference methods. This was done with the R TSP package [START_REF] Hahsler | TSP-infrastructure for the traveling salesperson problem[END_REF].

Theoretical fields

Methodology

Simulated data were used to study the properties and limitations of the approach. Theoretical fields are intended to compare CS to reference methods in a wide range of known situations. Each simulation aims at providing two variables for each theoretical field: auxiliary observation (i.e. NDVI) and variable of interest (i.e. the yield), both spatially auto-correlated.

The simulation assumed that a main underlying phenomenon (i.e. environmental factors like soil, climate, elevation, etc.) drives the within field variability of the plant response. The simulation process therefore starts by generating a theoretical auto-correlated variable (noted G ), representing the spatial variability of the underlying factor. G is simulated as a spatialized Gaussian field with no nugget effect (Fig. 3a). Two new variables, respectively V1 (Fig. 3b) and V2 (Fig. 3e), were derived from G by adding a non-auto-correlated noise following a nor- mal centred distribution of respective variances V 1 and V 2 (Eqs. 5 and 6). V1 and V2 are therefore intrinsically correlated with each other, the level of correlation depends on V1 and V2 . In the followings, V1 will be used for the variable of interest while V2 will be used for the auxiliary variable.

(3) Four parameters were considered to vary across the theoretical dataset: (i) the distance of auto-correlation of G, V1 and V2, defined by the range of the semi-variogram of each variable; (ii) the ratio of nugget effect/sill of the V1 variable; (iii) the degree of correlation between V1 and V2 and (iv) the number of outliers on V2. These parameters fit the diversity of theoretical data by controlling the link between the two variables and their semi-variograms. Hereafter is described the process used to set up these parameters during the simulation process:

Error(%) = | | | Y - Ŷ| | | Y (4) ŶRS = mean s∈S Selected (GW(s)) (5) V1 i = G i + V1 i with V1 i ∼ N 0, 2 V1 (6) V2 i = G i + V2 i with V2 i ∼ N 0, 2 V2 With Cor (V1, V2) [0,1]
2 V1
, the variance of the noise added to V1 was chosen to obtain the expected nugget effect/sill ratio on V1. G having no nugget effect, 2

V1 is therefore directly equal to the nugget effect of V1 . It can be directly deduced from the ratio and 2 s (Eq. 7):

2 V2 , the variance of the noise added to V2 was used to calibrate the degree of correlation between V1 and V2 . This was deduced from the covariance and Pearson correlation formu- las (Eqs. 8 and 9). The variable of interest will be the yield and the auxiliary data, the NDVI. They were respectively derived from V1 and V2 by a linear change of scale ( ) to be centred around

(7) Ratio = nugget(Y) sill(Y) = 2 V1 2 G ⇒ 2 V1 = Ratio × 2 G (8) Cov(V1, V2) = Cov G + V1 , G + V2 Cov(V1, V2) = Cov(G, G) + Cov G, V1 + Cov G, V2 + Cov V1 , V2 Cov(V1, V2) = Cov(G, G) = Var(G) = 2 G Cor(V1, V2) = Cov(V1, V2) √ Var(V1) × √ Var(V2) (9) Cor(V1, V2) = 2 G √ 2 sG + 2 V1 × √ 2 G + 2 V2 2 V2 = ⎛ ⎜ ⎜ ⎜ ⎝ 2 G � 2 G + 2 V1 × Cor(V1, V2) ⎞ ⎟ ⎟ ⎟ ⎠ 2 -2 G
the desired average yield ( Y theo ) and the desired average NDVI ( NDVI theo ) with the appropriate dispersion (Eqs. 10 and 11) This transformation does not affect the correlation or any of the three other parameters. (Fig. 3c, f and Eqs. 10 and 11):

An optional step consists of adding outlier zones on the NDVI simulated maps. Outlier zones intended to represent abnormal phenomenon like weed patches (abnormally strong NDVI) or local diseased vines (abnormally low NDVI) who may locally alter the correlation between yield and NDVI.

The number of outlier zones will vary from 0 to 6 (Table 1). The Location of each outlier zones was randomly chosen. Their size varies from 10 to 30 pixels and all these pixels have the same NDVI value taken from [0.1,0.25] ∪ [0.75,0.9] . All these parameters are drawn randomly for each desired outlier zones. Pixels around the outlier zone were smoothed in order to simulate a short gradient with "normal" surrounding NDVI values (Fig. 3g). Introduction of outlier zones may lead to a slight decrease for the correlation parameter.

The final step consisted in extracting the values of both information (Yield and NDVI) corresponding to the rows of the vine. The rows take the values of the nearest pixel. (Fig. 3d andh).

Implementation

Theoretical fields were designed with an area of one hectare (100 m × 100 m) with a 2.5 m distance between rows which corresponds more or less to typical fields area and plantation density of 4000 vines/ha found in south of France. The resolution was 1 pixel/m 2 . For theoretical yield data, parameters of the simulation were defined so that average yield corresponds to common average yield of the region ( -Y theo = 1000g∕vine) and Coefficient of Vari- ation (CV) to previous works: CV = 30%; 2 G = (1000 × 0.3) 2 (Krstic et al. 1998, Dunn and Martin 2000). For theoretical NDVI, the average value was set at -NDVI theo = 1∕2 , and the factor in Eq. 11 at 1 18 G to ensure that all the NDVI values lay within the range of [0,1].

The choice of the different possible values for the four parameters (Table 1) was based on observations from literature in precision viticulture [START_REF] Bramley | Spatio-temporal variability in vine vigour and yield in a Marlborough Sauvignon Blanc vineyard[END_REF][START_REF] Bramley | Understanding variability in winegrape production systems. 1. Within vineyard variation in yield over several vintages[END_REF][START_REF] Hall | Within-season temporal variation in correlations between vineyard canopy and winegrape composition and yield[END_REF][START_REF] Li | Spatial variation of winegrape yield and berry composition and their relationships to spatiotemporal distribution of soil water content[END_REF][START_REF] Taylor | A comparison of the spatial variability of vineyard yield in European and Australian production systems[END_REF][START_REF] Acevedo-Opazo | The potential of high spatial resolution information to define within-vineyard zones related to vine water status[END_REF]. For each parameter, three possible levels were defined (Table 1), encompassing a range of variability allowing to account for the existing diversity in vineyard fields. Each parameter will ( 10) 1. This procedure will ensure to test the effect of each parameter on the sampling results.

Y i = Y theo + V1 i (11) NDVI i = NDVI theo + V2 i ×

Real data

Real fields used to test the method come from INRA Pech-Rouge (Narbonne, France). The experiment and the data base was detailed by [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF]. It is briefly summarised hereafter. NDVI values from nine different vine fields were considered. All of them are non-irrigated and exposed to Mediterranean climate with precipitation occurring during spring with hot and dry summer. The characteristics of each plot are shown in Table 2. NDVI values were derived from a 1 m. resolution multi-spectral image taken the 31th of August 2008 by Avion Jaune (Montpellier, Hérault, France). The spectral regions captured in the images were: blue (445-520 nm) green (510-600 nm), red (632-695 nm) and nearinfrared (757-853 nm). From these 1 m square image pixels, aggregation method described in [START_REF] Acevedo-Opazo | The potential of high spatial resolution information to define within-vineyard zones related to vine water status[END_REF]) was used to obtain 9 m square image pixels reducing the effect of canopy discontinuity and bare soil on measured values. NDVI was finally computed from processed images according to [START_REF] Rouse | Monitoring vegetation systems in the great plains with ERTS[END_REF]. Mechanical or chemical weeding was performed over the inter-row spacing; therefore, row cover crop did not affect much NDVI values. PSSs were selected regularly over the fields with measurement made on each node of a 15 m 2 width sampling grid (Fig. 4). At each node, yield was measured on five consecutive vines in the row and average yield was affected to the location corresponding to the central vine. The final data base was a set of 313 sites over the nine different fields. It is noteworthy that, unlike simulated fields, the number of PSSs is reduced for real fields. Indeed, since it was not practically possible to measure the yield on each vine stock, the consequence is that the number of PSSs depends on the number of available measurement sites, therefore PSSs per site varied from 19 to 45. Each PSS was then characterized by a Grape Weight per vine value (GW) and a NDVI value.

For each field, the average of all available measured GW values was used to estimate the average yield of the field ( Y).

Implementation

The core of the sampling approach was written in Java, the program used the Choco solver [START_REF] Prud'homme | Choco documentation[END_REF]. The calculations to obtain the distance matrix are made with Python. Theoretical fields, quantile classification according to their value for auxiliary data, estimation errors, and route distance were computed with R. Packages "gstats" and "sp" were used to generate Gaussian fields for the G variable.

As explained in the description of the constraints, the approach presented here takes into account the starting site of the practitioner to include it in the sampling circuit. Varying the starting site thus changes the result that will be obtained. In order to increase the number of situations tested for real data, this starting site is positioned on different ends of row across the vineyards. The approach was then applied to 86 different situations instead of 9. Results for the different starting sites were averaged for each field. For theoretical data, thirty simulations per set of parameters are tested (270 in total). CS was applied once to each simulated plot. The starting site is located on one of the corners of the plot. The two outer rows on each side and the first three vines of each row cannot be selected as, in practice, they are subject to border effects. For both theoretical and experimental data, RS and MS were applied 1000 times. These repetitions were possible as they rely on random or partially random selection of SSs. The following figures are based on the average of the results obtained on each field, the result of each field being the average across repetitions with different starting site.

Results and discussion

Sampling and vine diversity

Number of sampling sites

Figure 5 presents the results obtained for the three sampling approaches, constrained sampling (CS), model sampling (MS) and random sampling (RS) for simulated fields with parameters set at their default values (Table 1). The different sampling approaches were tested on each simulated plot with a number of measurement sites ranging from N = 5 to N = 10. Results represent the averages obtained with thirty simulated plots.

All the sampling methods follow the same trend with a decreasing error as the number of sampled sites increases. This result is logical, and consistent with the literature. As [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] have already shown, taking into account auxiliary data, MS approach slightly improves the quality of yield estimation compared to a RS (Fig. 5a). With default values, it seems that the MS approach proposed by [START_REF] Carrillo | Use of multi-spectral airborne imagery to improve yield sampling in viticulture[END_REF] gives the best results in terms of estimation error. CS and RS both present similar errors whatever the number of SSs. Observed errors with CS are higher than for MS (i.e. without constraints). This result may be logical considering that the addition of the constraints may limit possibilities when choosing among the PSS.

Figure 5b clearly illustrates the gains brought by CS in terms of travel distance across the vineyard. Logically, the travelled distance within the plot increases linearly with N, the number of SSs. Travel distances with CS are at least 85% better than with other sampling methods. Distances are also less sensitive to the increase in the number of SSs with CS. This is explained by the gain brought by considering this criterion when selecting SSs. In view of these first results, CS offers a compromise between MS and distance criterion optimization, with a higher estimation error in favour of a significant reduction in distance.

With the hypothesis of a walking speed of 0.9 m/s and 60 s needed per SS, Table 3 shows that MS and CS perform better than RS. For a given amount of time, CS allows a higher number of observations to be made compared to other sampling approaches and therefore the lowest estimation error.

Impact of the semi-variogram range

Figure 6 shows the estimation error and the distance results obtained with the 3 sampling approaches tested on theoretical fields with varying range (10 m, 20 m, 40 m) of semi-variograms. Varying the range does not affect significantly the estimation error in function of the number of SSs compared to previous results. MS still presents the best estimation error compared to CS and RS (Fig. 6a). For a range of 40 m, CS seems more erratic with the highest error for a low number of SSs (N > 8) and an estimation error which tends to the error observed with MS for a higher number of SSs (N > 8) As for estimation error, results of distances associated with RS and MS do not seem to be affected by the range parameter either (Fig. 6b). The increase in range is nevertheless associated with a slight increase in distance with CS.

Figure 7 gives an illustration of the sampling route obtained with the method for plots with different ranges (10, 20 and 40 m).

The method always promotes the measurement sites in the immediate vicinity of the starting site (0.0 coordinates). This is expected as the method aims to minimize travel time. A robustness study (result not shown) showed that the results were similar regardless of the starting site, the method always promote measurement sites close to the starting point whatever the plot characteristics. Figure 7 shows that the average distance between sampling sites tends to increase with the plot range. On these plots, the maximum distance between sampling points corresponds approximately to the range of the plot; it is of 40 m, 30 m and 10 m respectively for plot range of 40 m, 20 m and 10 m. This result is expected since as the range increases, the distance to find a higher diversity of values also increases leading to longer sampling routes for larger ranges. This result is consistent with that of Fig. 6b. One surprising aspect, at first glance, is the close proximity of the sampling points. This characteristic is related to the nugget effect, which introduces erratic variance and high variability over short distances. This erratic variance makes it possible to find a high variability of values in the immediate vicinity of the starting point and other measurement sites. The sampling method takes advantage of this variability to minimize travel time. It should be remembered that in these simulations, the nugget effect was defined according to literature figures in precision viticulture. It is rather high since it represents more than 30% of the plot variability. A decrease in this nugget effect leads to longer sampling route. For example, if the nugget effect is set to 0 (no erratic variability), then the sampling distance is longer since the distance needed to find representative values necessarily increases. Conversely in the case of no spatial autocorrelation, the method chooses contiguous independent measurement sites on a same row and the sampling distance is in this case very short (result no shown). Figure 7 therefore represents the result of optimal sampling that takes into account the combined effect of the nugget effect and the range considering realistic figures of spatial variability in viticulture.

Impact of the ratio nugget/sill

Figure 8 shows the impact of the ratio for the three sampling strategies. Whatever the sampling method and the number of sampling points, an increasing ratio affect the estimation error (Fig. 8a). This result is logical since the increase in the ratio corresponds to an increase in the proportion of erratic (non-autocorrelated) variance in the total variance of the NDVI. Thus, sampling methods tend to select SSs whose NDVI value is not necessarily correlated to the expected yield value. Note that for a high ratio (ratio = 0.5), the estimation error with CS is more affected compared to other approaches. CS approach appears to be more sensitive to high ratio values. Indeed, the short-range variability introduced by higher erratic variance increase the chance to have sites in different quantile groups close to each other. The CS approach that minimizes the route from the starting site might select close SSs. Close SSs often provide redundant information, which leads to an increase in estimation error. On the other hand, the ratio does not seem to have any significant influence on the length of the sampling circuits (Fig. 8b).

Impact of the correlation level with auxiliary data

As expected, RS is not affected by a low level of correlation since SSs selection is not based on the NDVI data (Fig. 9a). This is not the case for MS and CS which both show lowest estimation errors when correlation between yield and NDVI is high while they show high estimation errors decreases when the correlation decreases. Since, MS and CS use the relationship between the two variables, it was expected that the quality of prediction decreases when this relationship is weakened. When the correlation is close to 0, MS tends to have the same results as RS while CS presents the worst estimation error. For the theoretical dataset, a very low correlation corresponds to a strongly noisy NDVI (Eq. 9). The resulting short range erratic variability in NDVI could explain this result for CS. It follows the same phenomenon as already observed in the previous section for an increasing ratio. These results highlight the sensitivity of CS to noise in NDVI for the selection of SSs aiming at optimising the distance. As for previous result, sampling distance does not seem much affected by the correlation, apart for CS for which it seems to slightly decrease with the correlation (Fig. 9b). This result supports the idea of a decreasing correlation allowing the CS approach to find SSs closer to each other.

Impact of outlier zones

The addition of local outlier (Fig. 10a) seems to slightly affect all methods, including the RS which should not be affected by changes in auxiliary data. The effect of outlier zones is still more important on the CS. It is difficult to draw conclusions as the effect is not proportional to the number of outlier zones. As it could be expected, the length of the sampling circuits is not affected by these local outlier zones (Fig. 10b). other, while this was possible for theoretical fields. Overall, this method offers a good compromise between the quality of the estimate and the travel constraint on the plot.

Table 4 compares the performance of the different approaches on sampling time and estimation error for real data. Walking speed is set at 0.9 m/s and 60 s are required for each SS. As for simulated data, CS perform better (up to 30%) than other approaches for a given amount of time.

Further reflections

The method presented in this paper aims to select SSs accounting for variability highlighted by auxiliary data and the practitioner constraints simultaneously. It is intended to be general enough to be applicable to various combinations of auxiliary data and variable of interest. This type of strategy could be applied to any type of crop where the travel route of operators is constrained by the organization of the crop (trellised structure). It may be necessary to keep in mind the assumptions on which the approach is based on: the fact that there is a correlation between the variable of interest and the auxiliary variable and the relevant auxiliary variable is available with a high spatial resolution.

The choice of the auxiliary variable depends on the variable of interest to be estimated and the pedo-climatic context of the crop. Returning to the application case presented in this paper, it is based on preliminary studies that indicate that the NDVI measured at veraison was a relevant auxiliary variable to guide yield sampling in the specific context of non-irrigated Mediterranean vineyard. However, this same auxiliary variable may not be relevant in other soil and climatic contexts. As a result, the choice of the variables used requires prior knowledge and expertise for the correct implementation of the approach.

Overall, it seems that CS offers a strategy with relevant compromise between estimation error and sampling distance. It significantly improves the performance of the latter criterion for a small increase of the estimation error. For a given amount of time, CS presents better results than reference methods. For a given number of SSs, the gap for the average error between CS and MS, which also uses a model for estimating the variable of interest, could be explained by two points. First, the optimization of the distance tends to regroup the SSs around the starting point (Fig. 6), and it has been shown that this could be a disadvantage when the auxiliary data is noisy, indeed in this case, CS favours the choice of different NDVI values very close to each other and not necessarily correlated with the variable of interest. Conversely, the other methods can be more representative as they are more likely to choose SSs anywhere over the plot. The second reason is that choosing points close to each other reduces the distance of the circuit but might increase the autocorrelation between measurements. Two close sites will provide more or less redundant information depending on the distance between them. It should be noted that on the plot, practitioners generally rely on a random sampling limited to a fraction of the plot (often a pair of rows) as they cannot cover the whole field. The results obtained would then correspond to a degraded version of the random sampling due to the same auto-correlation problem.

A potential solution to minimize the gap between MS and CS would be to set a minimum distance between two SSs. This minimum distance could be derived from spatial auto-correlation of available auxiliary data or historical yield maps. This distance would be determined with the range of the experimental semi-variogram obtained with the auxiliary data. Depending on the value, this could make it possible to avoid or limit the autocorrelation issue between the SSs and to better results of CS for slightly longer sampling times. Another practical option when the auxiliary data is particularly noisy (ratio > 0.33), would be to smooth the data using, for example, a moving window average. This was done on real fields of this study and gave particularly interesting results. The methodology proposed by [START_REF] Tisseyre | How to define the optimal grid size to map high resolution spatial data[END_REF] could be used to decide on the size of the smoothing window to be used according to the erratic variance that is to be eliminated. Each of these areas for consideration could improve the sampling method presented here with reduced estimation errors and more widely spaced measurement sites.

Other areas for improvement are possible, such as, for example, the ability to integrate multiple auxiliary data or to use more complex models.

Conclusion

The methodology presented in this paper describes a new approach, Constrained Sampling (CS), for yield sampling in viticulture. The originality of the approach comes from the association of method from Carrillo et al. based on auxiliary data and optimisation algorithms to propose relevant sampling routes in term of estimation error and travelled distance. While the model sampling principle guides sampling points choice considering auxiliary information, optimisation through constraint programming ensures the relevancy of the chosen route in term of walking distance for the practitioner. CS appears however sensitive to unfavourable situations (low correlation, poor spatial structure, erratic variance of the auxiliary data) while other methods relying on random aspect may fare better. In favourable situations (good correlation between auxiliary data and yield and strong spatial structure), CS gives very good results. The estimation error is close to what is proposed by Carrillo et al. However, CS makes it possible to obtain much shorter sampling in distances and times. This saved time can then be used to increase the number of measurements and the reliability of the estimation.

Fig. 1

 1 Fig. 1 Distances across vineyards. The left illustrates how the vineyard structure affect the moving from point a to point b. The others illustrate the four different ways of going from point c to point b

Fig. 2

 2 Fig. 2 NDVI values for groups with quantile approach and N = 5. Each group contains 20% (100/N) of PSS according to their NDVI values

Fig. 3

 3 Fig. 3 Workflow of the theoretical fields simulation process. a Variable G is generated as a fully spatialized Gaussian field. b Variable V1 derived from G by adding a random noise. c Yield variable derived from V1(linear transformation). d Final theoretical yield data after row extraction. e Variable V2 derived from G by adding a random noise. f Variable NDVI derived from V2 (linear transformation). g NDVI variable complemented with outlier zones. h Final NDVI data after row extraction

  Hence:G , V1 and V2 being independent random variables, their covariance are equal to zero. Finally:With the Pearson correlation formula: it results in: And finally:

Fig. 4

 4 Fig. 4 INRA Pech Rouge plot with row edges in blue and potential sampling sites in red (Color figure online)

Fig. 5

 5 Fig. 5 Results for theoretical data with default parameters in function of the number of sampling sites; a estimation error, b sampling route distance

Fig. 6

 6 Fig. 6 Effect of the semi-variogram range on sampling strategies. a Estimation error & b sampling route distance

Fig. 7

 7 Fig. 7 Illustration of sampling routes for N = 9 and three different ranges: a range = 40 m; b range = 20 m; c range = 10 m

Fig. 8

 8 Fig. 8 Effect of the proportion of the ratio nugget/sill on sampling strategies. a Estimation error & b sampling route distance

Fig. 9

 9 Fig. 9 Effect of the correlation between NDVI and yield on sampling strategies. a Estimation error & b sampling route distance

Table 1

 1 Values for theoretical field parameters setting the others to their default values, indicated in bold in Table

	Parameter	Low	Medium	High
	Range (m)	10	20	40
	Ratio Nugget effect/Sill	1/10	1/3	1/2
	Pearson correlation coefficient	0.1	0.4	0.7
	Number of Outlier zones	0	3	6
	Default values are in bold font

Table 2

 2 Description of the experimental fieldsNugget effect could not be estimated because of the resolution of yield data

	Field	Area (ha)	Variety	Number of potential	Range of semi vari-	Pearson correlation
				sampling sites	ogram yield (m)	coefficient (NDVI/
						yield)
	P22	1.72	Syrah	45	21.14	0.13
	P63	1.33	Syrah	42	7.37	0.28
	P65	0.69	Syrah	33	27.71	0.86
	P76	1.14	Carignan	37	22.20	0.39
	P77	1.24	Syrah	19	9.25	0.48
	P80	0.54	Syrah	40	20.34	0.63
	P82	1.15	Syrah	53	21.69	0.47
	P88	0.85	Syrah	21	14.08	-0.04
	P104	0.81	Carignan	23	12.96	0.18

Table 3

 3 Sampling time and estimation error for RS, MS and CS on simulated data

	N Random sampling	Model sampling	Constrained sam-
	(RS)		(MS)		pling (CS)
	Error (%) Time (s) Error (%) Time (s) Error (%) Time (s)
	5 14,6	767	10,8	773	15,2	347
	6 13,5	877	9,8	881	13,9	411
	7 12,8	993	9,1	994	12,9	474
	8 12,1	1111	8,5	1111	11,8	537
	9 11,5	1192	7,9	1223	10,6	603
	10 11,1	1387	7,6	1371	9,5	660

Acknowledgements This work was supported by the French National Research Agency under the Investments for the Future Program, referred as ANR-16-CONV-0004 (#Digitag).

1 3

Evaluation of sampling strategies on real data

Figure 11 shows the averaged result of the three sampling methods on real fields. The same logical decrease in error as the number of sampled sites increases is observed (Fig. 11a.) The irregularity of the curves associated with CS can be explained by a smaller number of experiments compared to results obtained with theoretical fields. With N = 6 put aside, CS estimation errors are just halfway between those of RS and MS. Considering the characteristic of real fields in terms of ratio and correlation between yield and NDVI (Table 1), these results are consistent with the results obtained on the simulated data. CS has been shown to be sensitive to the correlation between yield and NDVI and noise in NDVI values. The average correlation between yield and NDVI being quite large for real data (0.38 on average, Table 1) and NDVI values being smoothed, real fields present average characteristics close to the default values for theoretical fields. The reduction of erratic local variability in NDVI may have favoured results observed with real fields with CS. Figure 11b illustrates the gains in terms of travel distance across the vineyard when using CS compared to other sampling approaches. Distance is reduced by approx. 50% with CS compared to MS and RS. This result is again consistent with those observed with theoretical data. However, the gain here is much smaller than with theoretical data because of the lower number of PSSs available for real fields Considering the minimum distance between two PSSs is 15 m, the distance required to travel through the selected SS is necessarily higher for real fields because of the impossibility for the algorithm to find SSs closer than 15 m to each


Augustin Mouze
email: augustin.mouze@univ-lille.fr

Augustin Mouze Universal Taylor

UNIVERSAL TAYLOR SERIES WITH RESPECT TO A PRESCRIBED SUBSEQUENCE

Keywords:

For a holomorphic function f in the open unit disc D and ζ ∈ D, Sn(f, ζ) denotes the n-th partial sum of the Taylor development of f at ζ. Given an increasing sequence of positive integers µ = (µn), we consider the classes U(D, ζ) (resp. U (µ) (D, ζ)) of such functions f such that the partial sums {Sn(f, ζ) : n = 1, 2, . . . } (resp. {Sµ n (f, ζ) : n = 1, 2, . . . }) approximate all polynomials uniformly on the compact sets K ⊂ {z ∈ C : |z| ≥ 1} with connected complement. We show that these two classes of universal Taylor series coincide if and only if lim sup n µ n+1 µn < +∞. In the same spirit, we prove that, for ζ = 0, we have the equality U (µ) (D, ζ) = U (µ) (D, 0) if and only if lim sup n µ n+1 µn < +∞. Finally we deal with the case of real universal Taylor series.

Introduction

As usual N, Q denote the sets of positive integers and rational numbers respectively. Let D := {z ∈ C : |z| < 1} be the open unit disc of the complex plane. Throughout the paper, H(D) denotes the vector space of all holomorphic functions on D endowed with the topology of uniform convergence on all compact subsets of D. Also for a compact set K of C we denote by A(K) the set of all functions which are holomorphic in the interior K o of K and continuous on K. As usual, for a holomorphic function f in the unit disc and ζ ∈ D, S n (f ) or S n (f, ζ) stands for the n-th partial sum of the Taylor development of f with center at 0 or at ζ respectively. In 1996, Nestoridis proved the following result [START_REF] Nestoridis | Universal Taylor series[END_REF].

Theorem 1.1. [START_REF] Nestoridis | Universal Taylor series[END_REF] There exist Taylor series f = n≥0 a n z n such that, for every compact set K ⊂ {z ∈ C; |z| ≥ 1} with connected complement and for every function h ∈ A(K) there exists a subsequence (λ n ) ⊂ N such that S λn (f ) converges to h, as n → +∞, uniformly on K.will

In the sequel, such Taylor series will be called universal Taylor series and we will denote by U (D, 0) the set of such universal Taylor series. In the same spirit, for ζ ∈ D, we can replace S λn (f ) by S λn (f, ζ) in the previous theorem to obtain the class U (D, ζ) of universal Taylor series with center ζ. The sets U (D, 0) or U (D, ζ) enjoy very strong properties. For example, these sets are G δ dense subsets of H(D) and contain, apart from 0, a dense vector subspace of H(D) [START_REF] Bayart | Abstract theory of universal series and applications[END_REF][START_REF] Nestoridis | Universal Taylor series[END_REF]. This last property means that the sets U (D, 0) and U (D, ζ) are algebraically generic. A very nice theorem asserts that for all ζ ∈ D, U (D, 0) = U (D, ζ). We refer the reader to [START_REF] Gehlen | On the existence of O-universal functions[END_REF] and [START_REF] Melas | Universality of Taylor Series as a Generic Property of Holomorphic Functions[END_REF]. A crucial tool for the proof is a result initiated by Gehlen, Luh and Müller [START_REF] Gehlen | On the existence of O-universal functions[END_REF] which asserts that every universal Taylor series actually possesses Ostrowski-gaps, in the sense of the following definition. Definition 1.2. Let ζ ∈ C. Let +∞ j=0 a j (zζ) j be a complex power series with radius of convergence r ∈ (0, +∞). We say that it has Ostrowski-gaps (p m , q m ) if (p m ) and (q m ) are sequences of natural numbers with (1) p 1 < q 1 ≤ p 2 < q 2 ≤ . . . and lim m→+∞ qm pm = +∞, (2) for I = ∪ ∞ m=1 {p m + 1, . . . , q m }, we have lim j∈I |a j | 1/j = 0. The fact that every universal Taylor series possesses Ostrowski-gaps is at the core of many beautiful results (see for instance [START_REF] Bayart | Boundary behavior and Cesàro means of universal Taylor series[END_REF][START_REF] Charpentier | On countably universal series in the complex plane[END_REF][START_REF] Gehlen | On the existence of O-universal functions[END_REF][START_REF] Katsoprinakis | Coincidence of some classes of universal functions[END_REF][START_REF] Melas | Universality of Taylor Series as a Generic Property of Holomorphic Functions[END_REF][START_REF] Mouze | Polynomial inequalities and universal Taylor series[END_REF]). Thus the combination of the existence of Ostrowski-gaps with a result of Luh (see Theorem 1 of [START_REF] Luh | Universal approximation properties of overconvergent power series on open sets[END_REF] and Lemma 2.3 below) allows to obtain the aforementioned equality U (D, 0) = U (D, ζ). Notice that we give a new proof of [START_REF] Luh | Universal approximation properties of overconvergent power series on open sets[END_REF]Theorem 1] in Section 2.1. Next, in order to prove the algebraic genericity of the class of universal Taylor series, the following subclass of universal series was introduced. Definition 1.3. Let ζ ∈ D. Let µ = (µ n ) be an increasing sequence of positive integers with µ n → +∞ as n tends to infinity. A holomorphic function f ∈ H(D) belongs to the class U (µ) (D, ζ) if for every compact set K ⊂ {z ∈ C; |z| ≥ 1} with connected complement and for every function h ∈ A(K) there exists a subsequence (λ n ) ⊂ N such that S µ λn (f, ζ) converges to h, as n → +∞, uniformly on K.

Obviously we have

U (µ) (D, 0) ⊂ U (D, 0) (or U (µ) (D, ζ) ⊂ U (D, ζ)).
In [START_REF] Mouze | Polynomial inequalities and universal Taylor series[END_REF], as a consequence of a more general result relating to the weighted densities of subsequences along which the partial sums of universal Taylor series realize the universal approximation, the authors exhibit non-trivial subsequences µ of N such that U (µ) (D, 0) = U (D, 0). For example, the sequences (µ n ) = (n 2 ), (µ n ) = (2 n ) or the sequence of prime numbers satisfy this property. In [START_REF] Mouze | Polynomial inequalities and universal Taylor series[END_REF]Section 5], it is asked to characterize the subsequences µ for which the equality U (µ) (D, 0) = U (D, 0) holds. In this paper, we are going to answer this question by establishing the following result.

(i) U (D, ζ) = U (µ) (D, ζ) (ii) lim sup n→+∞ µ n+1 µ n < +∞.
To prove the implication (2) ⇒ (1) we use in an essential way the fact that all universal Taylor series possess Ostrowski-gaps. To obtain the converse implication, we employ a constructive method based on a Bernstein-Walsh type theorem given by Costakis and Tsirivas (see Theorem 2.4 below), when they studied the phenomenon of disjoint universal Taylor series [START_REF] Costakis | Doubly universal Taylor series[END_REF]. As consequence of Theorem 1.4 we show the independence of the class U (µ) (D, ζ) with the center of expansion ζ provided that lim sup n→+∞ µ n+1 µ n < +∞. This phenomenon was already noticed in specific cases in [START_REF] Vlachou | Disjoint universality for families of Taylor-type operators[END_REF][START_REF] Vlachou | Subclasses of universal Taylor series and center independence[END_REF]. Furthermore, using a constructive method similar to that of the proof of Theorem 1.4 and the ideas of the new proof of [8, Theorem 1], we obtain the following characterization. The following assertions are equivalent:

(i) U (µ) (D, ζ) = U (µ) (D, 0) (ii) lim sup n→+∞ µ n+1 µ n < +∞.
We refer the reader to Corollary 2.6 and Theorem 2.7. Finally in the last section we deal with the case of real universal Taylor series.

2. Universal Taylor series versus universal Taylor series with respect to a prescribed subsequence 2.1. Preliminary results. In this subsection we state some results that we will use for the proof of the main theorem. On one hand we are interested in the fact that all universal Taylor series possess Ostrowski-gaps. Actually a slightly more precise result holds (se [9, Theorem 9.1]). Combined Lemma 2.1 with the definition of universal Taylor series we immediately deduce the following useful lemma. Lemma 2.2. Let ζ ∈ D. Let f be in H(D) and suppose that the Taylor series of f at 0 has Ostrowski-gaps (p m , q m ). Then for every sequence (r m ) with p m < r m ≤ q m the difference between partial sums S rm (f, ζ)(z) -S pm (f, ζ)(z) converges uniformly to zero (as m → +∞) on compact sets of C.

In connection with the Ostrowski-gaps, we also have the following result ([8, Theorem 1] or [START_REF] Melas | Universality of Taylor Series as a Generic Property of Holomorphic Functions[END_REF]Lemma 9.2]). The published proof uses the Hadamard three-circle theorem. Here we give an elementary proof of the result.

Lemma 2.3. Let f be in H(D) and ζ 0 , ζ ∈ D. Suppose that the Taylor series of f at ζ 0 has Ostrowski-gaps (p m , q m ). Then the difference

S p k (f, ζ)(z) -S p k (f, ζ 0 )(z) converges to zero (as k → +∞) uniformly on compact sets of D × C (ζ ∈ D, z ∈ C).
Proof. Without loss of generality, we can suppose that ζ 0 = 0 (and ζ = ζ 0 !). Let K ⊂ C and L ⊂ D be fixed compact sets. We define

r = sup ζ∈L |ζ| and M > max(2, sup{|z -ζ| : z ∈ K, ζ ∈ L}). Let us choose 0 < ε < 1 such that r(1 + ε) < 1 and 2εM < 1. We write, for all |z| < 1, f (z) = +∞ k=0 a k z k . Thus, for all j ≥ 0, we have f (j) (ζ) = j! +∞ k=j a k k j ζ k-j .
Using the equality

z k = k l=0 a k k l (z -ζ) l ζ k-l , we get, for all ζ ∈ L, (1)
S pm (f, ζ)(z) -S pm (f )(z) = pm j=0   +∞ k=j a k k j ζ k-j   (z -ζ) j - pm j=0   pm k=j a k k j ζ k-j   (z -ζ) j = pm j=0   qm k=1+pm a k k j ζ k-j   (z -ζ) j + pm j=0   +∞ k=1+qm a k k j ζ k-j   (z -ζ) j := A 1,m (ζ, z) + A 2,m (ζ, z)
Now we are going to estimate the series A 1,m (ζ, z)and A 2,m (ζ, z). On one hand, by the triangle inequality and the inequality k j ≤ 2 k , we get

sup ζ∈L sup z∈K |A 1,m (ζ, z)| ≤ pm j=0   qm k=1+pm |a k |2 k r k-j   M j .
Using the Ostrowski-gaps, one can find m 1 such that, for all m ≥ m 1 and 1 + p m ≤ k ≤ q m , |a k | < ε k . We deduce, using the property M > 2r,

(2) sup

ζ∈L sup z∈K |A 1,m (ζ, z)| ≤ pm j=0   qm k=1+pm (2εr) k   M r j ≤ 1 1 -2εr (2εM ) 1+pm .
On the other hand, we need the following classical inequality

(3) k j ≤ k j j! ≤ ek j j
and the fact that, for t < k, the function t → ek t t is increasing. It follows

(4) sup ζ∈L sup z∈K |A 2,m (ζ, z)| ≤ pm j=0   +∞ k=1+qm |a k |e pm k p m pm r k-j   M j .
Since p m /q m → 0, as m tends to infinity, we have, for all k ≥ 1 + q m and 0 ≤ j ≤ p m ,

e pm k-j e pm k-j log( k pm ) ≤ e pm 1+qm-pm e pm k(1-j/k) log( k pm ) ≤ e pm 1+qm-pm e - pm k(1- pm 1+qm ) log( pm k ) → 1 as m → +∞.
Therefore, since lim sup |a k | 1/k ≤ 1 and p m /q m → 0, as m tends to infinity, there exists a positive integer m 2 such that for all m ≥ m 2 , k ≥ 1 + q m and 0

≤ j ≤ p m , ke p m pm/(k-j) |a k | 1/(k-j) r ≤ (1 + ε)r.
Since the choice of ε ensures (1 + ε)r < 1 and M > (1 + ε)r, we deduce that for all m ≥ m 2

(5)

sup ζ∈L sup z∈K |A 2,m (ζ, z)| ≤ pm j=0   +∞ k=1+qm (1 + ε) k-j r k-j   M j ≤ ((1 + ε)r) 1+qm 1 -(1 + ε)r pm j=0 M (1 + ε)r j ≤ 1 + p m 1 -(1 + ε)r M (1 + ε)r pm ((1 + ε)r) 1+qm .
Finally taking into account p m /q m → 0 again and combining (1) with ( 2) and ( 5), we derive On the other hand, we will need a specific version of Bernstein-Walsh theorem. It is a polynomial approximation theorem which allows in some sense to control both the degree and the valuation of the polynomials. This elegant statement was given in [START_REF] Costakis | Doubly universal Taylor series[END_REF]. For given sequence (x n ), (y n ) of positive real numbers, the notation

sup ζ∈L sup z∈K |S pm (f, ζ)(z) -S pm (f )(z)| → 0, as m → +∞.
y n = O(x n ) means that the sequence (y n /x n ) is bounded. Theorem 2.4. [4, Theorem 2.1] Let (σ n ), (τ n ) be strictly increasing sequences of positive integers, let K ⊂ C \ D be a compact set with connected complement and let r ∈ (0, 1). If 1 ≤ τ n /σ n → +∞ as n → +∞ and if U is open in C with K ⊂ U ,
then there is θ ∈ (0, 1) so that for every h ∈ H(U ) there exists a sequence of polynomials (P n ) of the form Theorem 2.5. Let µ = (µ n ) be a strictly increasing sequence of positive integers. The following assertions are equivalent:

P n (z) =
(i) U (D, 0) = U (µ) (D, 0) (ii) lim sup n→+∞ µ n+1 µ n < +∞. Proof. (ii) ⇒ (i): assume that lim sup n→+∞ µ n+1 µ n < +∞. Since the inclusion U (µ) (D, 0) ⊂ U (D, 0)
is obvious, it suffices to prove U (D, 0) ⊂ U (µ) (D, 0). Let also f be in U (D, 0). Thus according to Lemma 2.1 for all compact subset K ⊂ C \ D with connected complement and for all h ∈ A(K), there exists two sequences of positive integers (p m ), (q m ) such that (1) the Taylor series of f at 0 has Ostrowski-gaps (p m , q m ), (2) and sup

z∈K |S pm (f )(z) -h(z)| → 0,
|S pm (f )(z) -S µ jm (f )(z)| → 0, as m → +∞.
From the triangle inequality we get

sup z∈K |S µ jm (f )(z) -h(z)| → 0, as m → +∞.
This implies f ∈ U (µ) (D, 0).

(i) ⇒ (ii): to do this, we assume that lim sup n→+∞ µ n+1 µ n = +∞ and it suffices to exhibit an universal series f ∈ U (D, 0) such that f / ∈ U (µ) (D, 0). By hypothesis there exists an increasing subsequence of positive integers (n j ) such that µ n j +1 µ n j → +∞ as j → +∞. We set, for all j ≥ 1, (6) u j = µ n j + 1, w j = µ n j +1 and v j = ⌊ √ u j w j ⌋.

Clearly there exists N 0 ∈ N such that, for all j ≥ N 0 , u j < v j < w j and

v j u j = ⌊ √ u j w j ⌋ u j ≤ µ n j +1 µ n j → +∞ as j → +∞ and w j 1 + v j = w j 1 + ⌊ √ u j w j ⌋ ≤ µ n j +1 (1 + µ n j )µ n j +1 → +∞ as j → +∞.
Let (f q ) be an enumeration of all the polynomials with coefficients in Q + iQ. Let (K m ) be a sequence of compact sets with connected complement and K m ∩ D = ∅ for every m ∈ N such that for every compact set K ⊂ C \ D with connected complement there exists n ∈ N such that K ⊂ K n (see [START_REF] Nestoridis | Universal Taylor series[END_REF]Lemma 2.1]). We consider an enumeration (K ms , f qs ), s = 1, 2, . . . , of all couples (K m , f q ), m, q = 1, 2, . . . . Let also (r l ) be an increasing sequence of real numbers with 0 < r l < 1 and r l → 1, as l → +∞. We fix z 0 ∈ C \ D. For all l ≥ 1, we set Kl = K l ∪ {z 0 } and note that Kl is a compact set with connected complement and Kl ⊂ C \ D. First we deal with Km 1 , f q 1 and r 1 . By applying Theorem 2.4, we find j 1 ∈ N and polynomials

P 1 (z) = v j 1 k=u j 1 c j 1 ,k z k , P1 (z) = w j 1 k=v j 1 +1 c j 1 ,k z k such that sup z∈ Km 1 |P 1 (z) -f q 1 (z)| < 1 2 2 , sup |z|≤r 1 |P 1 (z)| ≤ 1 2 2 and sup z∈ Km 1 | P1 (z) + f q 1 (z)| < 1 2 2 , sup |z|≤r 1 | P1 (z)| ≤ 1 2 2 .
Observe that we have, by the triangle inequality, sup

z∈ Km 1 |P 1 (z) + P1 (z)| ≤ sup z∈ Km 1 |P 1 (z) -f q 1 (z)| + sup z∈ Km 1 | P1 (z) + f q 1 (z)| < 2 2 2 .
Further we find j 2 ∈ N with j 2 > j 1 and polynomials

P 2 (z) = v j 2 k=u j 2 c j 2 ,k z k , P2 (z) = w j 2 k=v j 2 +1 c j 2 ,k z k such that sup z∈ Km 2 |P 1 (z) + P1 (z) + P 2 (z) -f q 2 (z)| < 1 3 2 , sup |z|≤r 2 |P 2 (z)| ≤ 1 3 2 and sup z∈ Km 2 | P2 (z) + f q 2 (z)| < 1 2 2 , sup |z|≤r 2 | P1 (z)| ≤ 1 3 2 .
Observe that we have, by the triangle inequality,

sup z∈ Km 2 2 i=1 P i (z) + Pi (z) ≤ sup z∈ Km 2 |P 1 (z) + P1 (z) + P 2 (z) -f q 2 (z)| + sup z∈ Km 2 | P2 (z) + f q 2 (z)| < 2 3 2
. We argue by induction. Suppose that for a natural number s ≥ 2 we have already defined integers

j 1 < j 2 < • • • < j s
and polynomials P i , Pi , i = 1, . . . , s such that [START_REF] Katsoprinakis | Coincidence of some classes of universal functions[END_REF] P i (z) =

v j i k=u j i c j i ,k z k , Pi (z) = w j i k=v j i +1 c j i ,k z k with (8) sup z∈ Kms s-1 i=1 P i (z) + P s (z) -f qs (z) < 1 (s + 1) 2 , sup |z|≤rs |P s (z)| ≤ 1 (s + 1) 2 (9) sup z∈ Kms | Ps (z) + f qs (z)| < 1 (s + 1) 2 , sup |z|≤rs | Ps (z)| ≤ 1 (s + 1) 2 .

and (10) sup

z∈ Kms s i=1

(P i (z) + Pi (z)) < 2 (s + 1) 2 .
Using Theorem 2.4, we find j s+1 ∈ N with j s+1 > j s and polynomials

P s+1 (z) = v j s+1 k=u j s+1 c j s+1 ,k z k , Ps+1 (z) = w j s+1 k=v j s+1 +1 c j s+1 ,k z k such that sup z∈ Km s+1 s i=1 P i (z) + P s+1 (z) -f q s+1 (z) < 1 (s + 2) 2 , sup |z|≤r s+1 |P s+1 (z)| ≤ 1 (s + 2) 2 and sup z∈ Km s+1 | Ps+1 (z) + f q s+1 (z)| < 1 (s + 2) 2 , sup |z|≤r s+1 | Ps+1 (z)| ≤ 1 (s + 2) 2 .
Thus from the triangle inequality we get sup z∈ Km s+1 s+1 i=1

(P i (z) + Pi (z)) ≤ sup z∈ Km s+1 s i=1 P i (z) + Pi (z) + P s+1 (z) -f q s+1 (z) + sup z∈ Km s+1 | Ps+1 (z) + f q s+1 (z)| < 2 (s + 1) 2 .
The induction is valid. Finally we set

f (z) = i≥1 P i (z) + Pi (z) .
Thanks to the second inequalities of ( 8) and ( 9) this series converges on all compact subsets of D. So f ∈ H(D). The first inequality of (8) ensures that f ∈ U (D, 0). Indeed, let K ⊂ C \ D be a compact set with connected complement and h ∈ A(K). Set ε > 0 and s 0 ∈ N with 1/(s 0 + 1) 2 < ε/2. By hypothesis, one can find a positive integer p > s 0 such that K ⊂ K mp and sup K |h -

f qp | < ε/2. Thus from (8) we get sup z∈K |S v jp (f )(z) -h(z)| ≤ sup z∈Km p |S v jp (f )(z) -f qp (z)| + sup z∈Km p |h(z) -f qp (z)| < 1 (p + 1) 2 + ε 2 < ε.
Moreover observe that the property ( 6) implies that

{µ n ; n ≥ n j 1 + 1} ∩ (∪ s≥1 {u js , u js + 1, . . . , w js -1}) = ∅.
Therefore the equation [START_REF] Mouze | On doubly universal functions[END_REF] guarantees that, for all n ∈ N,

|S µn (f )(z 0 )| ≤ 2 s≥1 1 (s + 1) 2 = π 2 3 -2.
From this last inequality we easily deduce that f / ∈ U (µ) (D, 0). , the distance between ζ and the unit circle. By hypothesis there exists an increasing sequence of integers (n j ) such that µ n j +1 µ n j → +∞, as j → +∞.

We set, for all j ≥ 1, ( 11)

u j = µ n j + 1, w j = µ n j +1 and v j = ⌊ √ u j w j ⌋.
As in the proof of Theorem 2.5, we have for all j large enough u j < v j < w j and (12) v j u j → +∞ and

w j 1 + v j → +∞ as j → +∞.
Let us also consider an enumeration (K ms , f qs ), s = 1, 2, . . . , of all couples (K m , f q ), m, q = 1, 2, . . . , where (f q ) is an enumeration of all the polynomials with coefficients in Q + iQ and (K m ) is a sequence of compact sets with connected complement and K m ∩ D = ∅ such that for every compact set K ⊂ C \ D with connected complement there exists n ∈ N such that K ⊂ K n . Let also (r l ) be an increasing sequence of real numbers with 0 < r l < 1 and r l → 1, as l → +∞. For all l ≥ 1, we set Kl = K l ∪ {z 0 }. Observe that Kl is a compact set with connected complement and Kl ⊂ C \ D.

First we deal with Km 1 , f q 1 and r 1 . By applying Theorem 2.4, we find j 1 ∈ N and polynomial

P 1 (z) = w j 1 k=1+v j 1 a k z k such that, for all j ≥ j 1 , (13) 2|z
j 0 -(z 0 -ζ) j | = 2(1 -(1 -|ζ|) j ) > 1 and |f q 1 (z 0 )| ≤ w j 1 , sup Km 1 |P 1 (z) -f q 1 (z)| < 1 2
and sup

|z|≤r 1 |P 1 (z)| < 1 2 .
By the triangle inequality we have

(14) |P 1 (z 0 )| ≤ |P 1 (z 0 ) -f q 1 (z 0 )| + |f q 1 (z 0 )| ≤ 1 2 + w j 1 ≤ 1 + w j 1 .
We define

a 1+w j 1 = - P 1 (z 0 ) z 1+w j 1 0 -(z 0 -ζ) 1+w j 1 .
Thanks to ( 13) and ( 14), we have

|a 1+w j 1 | ≤ 2(1 + w j 1 ). Then we set R 1 (z) = a 1+w j 1 z 1+w j 1 .
By induction we construct an increasing sequence of integers (j l ), two sequences of polynomials (P l ) and (R l ) such that, for all l ≥ 2, (15)

µ n j l 1 + µ n j l-1 → +∞, as l → +∞ |f q l (z 0 )| ≤ w j l , (16)
P l (z) = w j l k=1+v j l a k z k , R l (z) = a 1+w j l z 1+w j l with (17) sup Km l P l (z) + l-1 i=1 (P i (z) + R i (z)) -f q l (z) < 1 2 l and sup |z|≤r l |P l (z)| < 1 2 l , and (18)
a 1+w j l = - P l (z 0 ) + l-1 i=1 (P i (z 0 ) + R i (z 0 )) z 1+w j l 0 -(z 0 -ζ) 1+w j l .
By the triangle inequality we have

|P l (z 0 )| ≤ 1 2 l + w j l + ≤ 1 + w j l and, taking into account (13), we deduce (19) |a 1+w j l | ≤ 2(1 + w j l ).
We also define, for all l ≥ 2, (20)

a k = 0 for k = 2 + w j l-1 , . . . , v j l .
By construction we get, for all l ≥ 1, using (18),

w j l j=0   1+w j l k=j a k k j ζ k-j   (z 0 -ζ) j = w j l j=0   w j l k=j a k k j ζ k-j   (z 0 -ζ) j + a 1+w j l w j l j=0 1 + w j l j ζ 1+w j l -j (z 0 -ζ) j = l-1 i=1 (P i (z 0 ) + R i (z 0 )) + P l (z 0 ) + a 1+w j l z 1+w j l 0 -(z 0 -ζ) 1+w j l = 0. Now we consider f (z) = +∞ i=1 (P i (z) + R i (z)) := k≥0 a k z k . (21)
Thanks to the estimates (17) and ( 19) the series f belongs to H(D). The equations ( 16) and (17) ensure that, for all l ≥ 1, sup

z∈Km l |S w j l (f )(z) -f q l (z)| < 1 2 l .
Since, for all l ≥ 1, w j l = µ 1+n j l , we deduce that f belongs to U (µ) (D, 0). To finish the proof, we are going to prove that f / ∈ U (µ) (D, ζ). Notice that we have

S m (f, ζ)(z 0 ) = m j=0   +∞ k=j a k k j ζ k-j   (z 0 -ζ) j .
First we deal with the case m = w j l . We have, thanks to the properties (20) and ( 21),

(22) S w j l (f, ζ)(z 0 ) = w j l j=0   +∞ k=1+v j l+1 a k k j ζ k-j   (z 0 -ζ) j .
Using the equation ( 15), we argue as in the estimates (3), ( 4) and ( 5) of the proof of Lemma 2.3 to obtain (23) S w j l (f, ζ)(z 0 ) → 0, as l → +∞.

Now, for all l ≥ 1 and for all positive integer m ∈ [1 + w j l , u j l+1 ], we can write, thanks to (20) and (21) again,

S m (f, ζ)(z 0 ) = w j l j=0   +∞ k=1+v j l+1 a k k j ζ k-j   (z 0 -ζ) j + m j=1+w j l   +∞ k=j a k k j ζ k-j   (z 0 -ζ) j = w j l j=0   +∞ k=1+v j l+1 a k k j ζ k-j   (z 0 -ζ) j + a 1+w j l (z 0 -ζ) 1+w j l + m j=1+w j l   +∞ k=1+v j l+1 a k k j ζ k-j   (z 0 -ζ) j := T 1,l (z 0 ) + a 1+µ n j l +1 (z 0 -ζ) 1+µ j l +1 + T 2,l (z 0 ).
By (19) we get |a 1+w j l (z 0ζ) 1+w j l | → 0, as l tends to infinity. Again inspired by the estimates (3), ( 4) and ( 5) of the proof of Lemma 2.3 we get, for i = 1, 2, T i,l (z 0 ) → 0, as l tends to infinity (let us recall that m ∈ [1 + µ j l +1 , u j l+1 ] and (1 + v j l+1 )/u j l+1 → +∞). In summary we have, for all positive integer m ∈ [w j l , u j l+1 ],

(24) S m (f, ζ)(z 0 ) → 0, as l → +∞.

By construction, we have [START_REF] Melas | Universality of Taylor Series as a Generic Property of Holomorphic Functions[END_REF][START_REF] Müller | Universal overconvergence and Ostroski-gaps[END_REF]). To see this, it suffices to note that Lemma 2.1, Lemma 2.2, Lemma 2.3 and Theorem 2.4 remain valid in this context. Thus the proofs work along the same lines.

(25) {µ n ; n ∈ N}   l≥1 ([1 + u j l+1 , w j l+1 -1] ∩ N)   = ∅.

The real case

As far as we know the first example of universal series was introduced by Fekete [14] which showed that there exists a real formal power series n≥1 a n x n satisfying the following universal property: for every continuous function g on [-1, 1] with g(0) = 0 there exists an increasing sequence (λ n ) of positive integers such that sup

x∈[-1,1] λn k=1 a k x k -g(x) → 0, as n → +∞.
Further combining this result with Borel's theorem we obtain C ∞ functions vanishing at 0 whose partial sums of its Taylor series with center 0 approximate every continuous function vanishing at 0 locally uniformly in R [START_REF] Erdmann | Universal families and hypercyclic operators[END_REF]. We denote by C ∞ 0 (R) the space of infinitely differentiable function on R vanishing at 0. Definition 3.1. Let µ = (µ n ) be an increasing sequence of positive integers with µ n → +∞ as n tends to infinity. A function f ∈ C ∞ 0 (R) belongs to the class U (µ) of universal functions with respect to µ if for every compact set K ⊂ R and every continuous functions h : R → R with h(0) = 0, there exists an increasing sequence (λ n ) of positive integers such that sup x∈K |S λn (f )(x)h(x)| → 0, as n → +∞.

For µ = N, we will denote U (µ) by U .

In this context, the analogue of Theorem 1.4 states as follows.

Theorem 3.2. Let µ = (µ n ) be a strictly increasing sequence of positive integers. The following assertions are equivalent:

(1) U = U (µ) (2) lim sup n→+∞ µ n+1 µ n < +∞.
To prove this result, we need the following results that are the C ∞ versions of Lemma 2.1 and Theorem 2.4 respectively. |P n (x)h(x)| → 0, as n → +∞.

Sketch of the proof of Theorem 3.2.

(2) ⇒ (1): thanks to Lemma 3.3, it suffices to mimic the proof (ii) ⇒ (i) of Theorem 2.5.

(1) ⇒ (2): we argue as in the proof of (i) ⇒ (ii) of Theorem 2.5. We define the same sequences (u j ), (v j ), (w j ) and we denote by (f j ) an enumeration of all the polynomials with coefficients in Q. Using Lemma 3.4, we construct step by step an increasing sequence of positive integers (j s ) and polynomials and polynomials P i , Pi , i = 1, . . . , s such that Thus f / ∈ U (µ) .

Theorem 1 . 4 .

 14 Let ζ ∈ D. Let µ = (µ n ) be a strictly increasing sequence of positive integers. The following assertions are equivalent:

Theorem 1 . 5 .

 15 Let ζ ∈ D, ζ = 0. Let µ = (µ n ) be a strictly increasing sequence of positive integers.

Lemma 2 . 1 .

 21 Let ζ ∈ D. Let f ∈ U (D, ζ). Let K ⊂ C \ D be a compact set with connected complement and let h ∈ A(K). Then there exist two sequences of positive integers (p m ), (q m ) such that (1) the Taylor series of f at ζ has Ostrowski-gaps (p m , q m ), (2) and sup z∈K |S pm (f, ζ)(z)h(z)| → 0, as m → +∞.

  Notice that the statements of [8, Theorem 1] or [9, Lemma 9.2] are given in the more general case where the open unit disc D is replaced by a simply connected domain Ω with Ω ⊂ C. Obviously the same proof does the job with easy modifications.

2 . 2 .

 22 τn k=σn c n,k z k with sup z∈K |h(z) -P n (z)| = O(θ τn ) and sup |z|≤r |P n (z)| = O(θ τn ). Proof of Theorem 1.4. In order to simplify the notations, we write the proof for the class U (D, 0). The proof works along the same lines in the case of the class U (D, ζ), ζ ∈ D.

2. 3 .Corollary 2 . 6 .

 326 Universal Taylor series and center independence. Let us recall that the classes U (D, 0) and U (D, ζ) coincide for all ζ ∈ D[START_REF] Gehlen | On the existence of O-universal functions[END_REF][START_REF] Melas | Universality of Taylor Series as a Generic Property of Holomorphic Functions[END_REF]. The proof is based on Lemma 2.1, Lemma 2.2 and Lemma 2.3 which asserts that if a Taylor series f has Ostrowski-gaps (p m , q m ) then the differenceS pm (f, ζ)(z) -S pm (f, 0)(z) converges to zero (as m → +∞) uniformly on compact sets of D × C (ζ ∈ D, z ∈ C). But if you can choose q m ∈ µ, there is no evidence that you can choose p m ∈ µ. Thus it is not clear that we have U (µ) (D, 0) = U (µ) (D,ζ). Nevertheless Theorem 1.4 immediately leads to the following corollary. Let ζ ∈ D. Let µ = (µ n ) be a strictly increasing sequence of positive integers with lim sup n→+∞ µ n+1 µ n < +∞. Then we have U (µ) (D, ζ) = U (µ) (D, 0). Proof. We know that, for all ζ ∈ D, U (D, ζ) = U (D, 0) [5, 9]. Assume that lim sup n→+∞ µ n+1 µ n < +∞. Hence Theorem 1.4 ensures that U (µ) (D, ζ) = U (D, ζ) and U (µ) (D, 0) = U (D, 0). We get U (µ) (D, ζ) = U (µ) (D, 0). Corollary 2.6 covers all the known examples of sequences µ such that U (µ) (D, ζ) = U (µ) (D, 0) [15, 16]. Moreover Corollary 2.6 is optimal in the following sense. Theorem 2.7. Let ζ ∈ D, ζ = 0. Let µ = (µ n ) be a strictly increasing sequence of positive integers with lim sup n→+∞ µ n+1 µ n = +∞. Then we have U (µ) (D, ζ) = U (µ) (D, 0). Proof. Let ζ ∈ D, ζ = 0. We are going to build an universal series f ∈ U (µ) (D, 0) such that f / ∈ U (µ) (D, ζ). Let us consider z 0 with |z 0 | = 1 such that |z 0 -ζ| = d(ζ, ∂D)

From ( 24 Example 2 . 8 .

 2428 ) and (25), we get S µn (f, ζ)(z 0 ) → 0, as n → +∞. It follows that f / ∈ U (µ) (D, ζ). This finishes the proof. Let ζ ∈ D, ζ = 0. We can apply Theorem 2.7 with the sequence (µ n ) = (n!). Therefore we obtain U ((n!)) (D, 0) = U ((n!)) (D, ζ). Remarks 2.9. (1) The combination of Corollary 2.6 with Theorem 2.7 gives Theorem 1.5. (2) Notice also that Theorem 1.4, Corollary 2.6 and Theorem 2.7 remain valid for the classes of universal Taylor series U (Ω, ζ) and U (µ) (Ω, ζ) where you replace the unit disc D by a simply connected domain Ω, with Ω ⊂ C and Ω = C, and ζ ∈ Ω, provided that the universal Taylor series possess Ostrowski-gaps (see

Lemma 3 . 3 .

 33 [START_REF] Mouze | Polynomial inequalities and universal Taylor series[END_REF] Proposition 4.4] Let f ∈ U . Let h : R → R be a continuous function, with h(0) = 0. There exist two sequences of natural numbers (λ n ), (µ n ) such that (1) the Taylor series of f around zero has Ostrowski-gaps (λ n , µ n ), (2) S λn (f ) → h, uniformly on each compact subset of R as n → +∞. Lemma 3.4. [10, Lemma 3.2] Let (l n ) and (m n ) be two strictly increasing sequences of positive integers such that l n ≤ m n and mn ln → +∞ as n → +∞. Let A > 0. For every continuous function h : R → R, with h(0) = 0, there exists a sequence (P n ) of real polynomials of the form P n (x) = mn k=ln c n,k x k , such that sup x∈[-A,A]

cP

  j i ,k x k , Pi (x) = w j i k=v j i +1 c j i ,k x k with (27) sup x∈[-s,s] s-1 i=1 P i (x) + P s (x)f s (x) < 1 (s + 1) 2 (28) sup x∈[-s,s] | Ps (x) + f s (x)| < 1 (s + 1) 2 . i (x) + Pi (x) + P s (x)f s (x) + sup x∈[-s,s] | Ps (x) + f s (x)| < 2 (s + 1) 2 .Finally let us consider the formal power series f (x) = i≥1 P i (x) + Pi (x) . By Borel's theorem one can find a function f ∈ C ∞ 0 (R) such that its Taylor development at zero is f . The inequality (27) ensures that f ∈ U . Moreover observe that the property[START_REF] Erdmann | Universal families and hypercyclic operators[END_REF] implies that∀n ≥ n j 1 + 1, ∀s ≥ 1 µ n / ∈ [u js , w js -1]and the equation (29) guarantees that, for all n ∈ N,sup x∈[-1,1]|S µn (f )(x)| ≤ 2

  as m → +∞. ∈ N, such that for all m ≥ N one can find µ jm with p m ≤ µ jm < q m . Hence we apply Lemma 2.2 to obtain sup z∈K

	Since we have both
	lim sup n→+∞	µ n+1 µ n	< +∞ and	q m p m	→ +∞, as m → +∞,
	there exists N

Acknowledgments

The author was partly supported by the grant ANR-17-CE40-0021 of the French National Research Agency ANR (project Front).

The author was partly supported by the grant ANR-17-CE40-0021 of the French National Research Agency ANR (project Front).