P Kervella
email: pierre.kervella@obspm.fr

T Verhoelst

S T Ridgway

G Perrin

S Lacour

J Cami

X Haubois

The close circumstellar environment of Betelgeuse ⋆

Keywords: individual: Betelgeuse, Stars: imaging, Stars: supergiants, Stars: circumstellar matter, Techniques: high angular resolution, Methods: observational

Context. Betelgeuse is one the largest stars in the sky in terms of angular diameter. Structures on the stellar photosphere have been detected in the visible and near-infrared as well as a compact molecular environment called the MOLsphere. Mid-infrared observations have revealed the nature of some of the molecules in the MOLsphere, some being the precursor of dust. Aims. Betelgeuse is an excellent candidate to understand the process of mass loss in red supergiants. Using diffraction-limited adaptive optics (AO) in the near-infrared, we probe the photosphere and close environment of Betelgeuse to study the wavelength dependence of its extension, and to search for asymmetries. Methods. We obtained AO images with the VLT/NACO instrument, taking advantage of the "cube" mode of the CONICA camera to record separately a large number of short-exposure frames. This allowed us to adopt a "lucky imaging" approach for the data reduction, and obtain diffraction-limited images over the spectral range 1.04 -2.17 µm in 10 narrow-band filters.

Results. In all filters, the photosphere of Betelgeuse appears partly resolved. We identify an asymmetric envelope around the star, with in particular a relatively bright "plume" extending in the southwestern quadrant up to a radius of approximately six times the photosphere. The CN molecule provides an excellent match to the 1.09 µm bandhead in absorption in front of the stellar photosphere, but the emission spectrum of the plume is more difficult to interpret. Conclusions. Our AO images show that the envelope surrounding Betelgeuse has a complex and irregular structure. We propose that the southwestern plume is linked either to the presence of a convective hot spot on the photosphere, or to the rotation of the star.

Introduction

Betelgeuse (α Ori, HD 39801, HR 2061) is one of the brightest stars at infrared wavelengths and has one of the largest apparent photospheric diameters (∼44 mas in the K band, from Perrin et al. 2004). It is a red supergiant of spectral type M2Iab with irregular flux variations. Because of its brightness and of its angular size, it has been an easy target for interferometers. It can even be resolved with large telescopes at short wavelengths, as demonstrated by Gilliland & Dupree (1996), who resolved its extended chromosphere using the Hubble Space Telescope at ultraviolet wavelengths. Several mysteries remain to be solved for such an object among which the structure of its convection and the mechanism of its mass loss. Similarities with Mira stars have been noted, including a significant thermal infrared excess, indicative of the presence of circumstellar dust. An important difference is that Betelgeuse is not subject to regular large amplitude pulsations which could levitate material high enough above the atmosphere where dust could form and be blown away by radiation pressure. The reason for the extension of the atmosphere is still unclear and another mechanism for dust formation and mass loss needs to be invoked. A careful scrutiny of the object from the photospheric scale up to the distance in its circumstellar envelope (CSE) where dust can condense is the path to understanding this prototypical and enigmatic star. It is also known from interferometric observations in the thermal infrared domain by Bester et al. (1996) that Betelgeuse experiences episodic mass loss. Our new NACO observations aim at probing the close environment of Betelgeuse to etablish a link between the photosphere and mass loss. After a presentation of our observations in Sect. 2, we detail our analysis of the NACO images in Sect. 3, and discuss the observed features in Sect. 4.

Observations and data reduction

Instrumental setup and observation log

We observed Betelgeuse on the nights of 2 and 3 January 2009 using the Nasmyth Adaptive Optics System (NAOS, [START_REF] Rousset | Proc. SPIE[END_REF] of the Very Large Telescope (VLT), coupled to the CONICA infrared camera (Lenzen et al. 1998), abbreviated as NACO. NAOS is equipped with a tip-tilt mirror and a deformable mirror controlled by 185 actuators, as well as two wavefront sensors: one for visible light and one for the infrared domain. For our observations, we exclusively used the visible light wavefront sensor. The detector is a 1024 × 1024 pixel ALADDIN InSb array. As its name indicates, NACO is installed at the Nasmyth focus of the Unit Telescope 4 (Yepun), located on the eastern side of the VLT platform. The CONICA camera was read using the newly implemented cube mode that allows to record very short exposures (down to 7.2 ms) over part of the detector. As discussed in Sect. 2.2, this new capability is essential to reach truly diffraction limited imaging.

We selected the smallest available pixel scale of 13.26 ± 0.03 mas/pix (Masciadri et al. 2003). Table 1 gives the list of the observations of Betelgeuse and the PSF calibrator stars Aldebaran (HD 29139), 31 Ori (HD 36167) and δ Phe (HD 9362). Due to the brightness of Betelgeuse and Aldebaran, we employed narrow-and intermediate-band filters together with a neutral density filter (labeled "ND2 short") with a transmission of about 1.5%. The observations of our fainter, secondary PSF calibrators (31 Ori and δ Phe) were obtained without neutral density filter. The transmission curves of the 10 narrow-band filters we used for our observations are shown in Fig. 1. They are also available through the NACO instrument web page at ESO1 .

Data processing

The individual raw images were pre-processed (bias subtraction, flat-fielding, bad pixel masking) using the Yorick2 and IRAF3 software packages in a standard way, except that we did not subtract the negligible sky background.

The quality of adaptive optics images is affected by a residual jitter of the star on the detector ("tip-tilt"), and by high-order residual wavefront distorsions. The star jitter is most probably caused by a combination of instrumental vibrations and atmospheric tip-tilt residuals, while the high-order distorsions are simply due to the limitations of the instrument performance. In order to obtain truly diffraction-limited image quality, especially at short wavelengths in the J band, it is essential to eliminate the perturbations from these two effects. Taking advantage of the newly implemented cube mode of CONICA, we obtained extremely short exposures that allowed us to freeze the image jitter. After a precentering at the integer pixel level, the images were sorted based on their maximum intensity, used as a proxy of the Strehl ratio. The 10% best images of each cube were then resampled up by a factor 4 using a cubic spline interpolation, and the star image was subsequently centered using gaussian fitting. The resulting cubes were eventually averaged to obtain the master images of each star used in the image analysis process described below. The initial oversampling of the image followed by the precise centering of the images gives an improved sampling of the image, as the residual atmospheric "drizzle" allows to retrieve sub-pixel sampling. This is particularly important in the J band, where the resolution of the telescope corresponds to only two NACO pixels. The selection of the best frames allowed us to reject the frames affected by the residual high-order distorsions. The result of this processing is presented in Fig. 2 for Betelgeuse and its primary PSF calibrator Aldebaran.

Efficiency of NACO full-pupil cube vs. standard mode

In order to check the efficiency of the image recentering and selection, we processed a series of 10 data cubes obtained on the PSF calibrator Aldebaran (#001 to #010 in Table 1) with and without recentering the individual frames. The resulting images for three wavelengths are shown in Fig. 3. Already clearly visible on these images, we hereafter quantify the gain in image quality by processing the same data cubes using three methods:

-recentering, selection of the best 10% frames based on the peak flux in the image and co-addition, -recentering and co-addition of all frames (no selection), -coaddition without recentering of all frames (equivalent to a single long exposure of 72 s).

From this exercise, it appears that the recentering of the frames before co-addition results in a spectacular increase of the Strehl ratio and a very significant decrease of the full width at half maximum (FWHM) of the star image (Fig. 4). As expected, this effect is stronger at shorter wavelengths, where the adaptive optics performance is comparatively lower. In the 1.04-1.09 µm range, the improvement reaches a factor three in Strehl ratio and a factor two in FWHM. The improvement in FWHM brought by the selection of the frames with the best Strehl appears very small compared to the effect of the recentering. Selecting the 10% best frames improves the average Strehl by approximately 10% (relatively to centering without selection) for the shortest wavelengths. For our very bright targets, we chose to keep this selection criterion as the signal-to-noise ratio of the average images is essentially limited by the AO residuals, and not by the number of frames.

As shown in Fig. 4, we can also compare the image FWHM with the theoretical diffraction limit of the telescope. Using a very simplified model of an 8 m primary mirror with a central 1.12 m obscuration, and a 20 mas uniform disk (UD) for Aldebaran, we obtained the theoretical FWHM shown as a thick grey curve. We chose this UD angular size as the angular diameter of Aldebaran was measured in the K band by Richichi & Roccatagliata (2005), who found θ UD (K) = 19.96±0.03 mas and θ LD (K) = 20.58 ± 0.03 mas. As a remark, the FWHM of our images is measured over a large sub-window that includes part of the adaptive optics photometric incoherent residuals (also called "pedestal"). The resulting values are thus slightly overestimated, but have the advantage of being directly comparable between real and theoretical images.

From these results, the frame-by-frame recentering of adaptive optics observations appears as a very efficient way to im- prove the angular resolution compared to single long exposures. The spectacular increase in Strehl ratio observed at short wavelengths should also result in an improved sensitivity for faint object, but the necessity of combining a large number of frames results in a proportionally higher detector readout noise. For this reason, this data reduction strategy may not be optimal for K band observations where the gain is "only" 50% in Strehl ratio. However, in the J band, whenever a bright enough source is available to recenter the frames, the cube mode of NACO with short exposures will in most cases be more sensitive than classical long exposures.

Image analysis

Point spread function

We obtained images of two unresolved and relatively fainter stars (31 Ori and δ Phe, #022-025 and #051-054 in Table 1) without neutral density filter, to check the influence of the ND filter on the PSF shape. We do not find evidence of an effect of this filter on the PSF shape. However, considering the faintness of the CSE of Betelgeuse, we chose to strictly keep the same instrumental setup between Betelgeuse and the PSF to avoid any contamination in the deconvolved images. We therefore used Aldebaran as the PSF reference. To retrieve the true PSF of the telescope from the images of Aldebaran, we took into account the angular extension of the star. Although this star is significantly smaller (θ LD = 20.6 mas) than Betelgeuse (θ LD = 43.7 mas), it is not negligible compared to the resolution of the telescope in the J band. We therefore computed synthetic images of Aldebaran using its known LD size and model limb darkening coefficients in the J, H, and K bands from Claret's (2000) tables using Teff = 4 000K and log g = 1.5. We then computed a classical Lucy-Richardson (L-R) deconvolution of the NACO images of Aldebaran (single iteration) using these synthetic images as the "dirty beams". This procedure gives us the true PSF of the telescope for our ten narrow-band filters, i.e. the image of a point source.

Photometry of Betelgeuse

For the photometric calibration of our images, we used a classical aperture photometry approach using Aldebaran as the calibrator. As we used non-standard narrow-band filters, we have to refer our measurements to a synthetic spectral energy distribution (SED) of Aldebaran, normalized by its interferometrically measured angular diameter.

The shape of the SED of Aldebaran was read from the theoretical spectra table by [START_REF] Castelli | Proc. IAU Symposium 210, Modelling of Stellar Atmospheres[END_REF], for an effective temperature of 4 000 K, solar metallicity, V T = 2 km.s -1 and log g = 1.5 (parameters approximated from Cayrel de [START_REF] Cayrel De Strobel | [END_REF]. We then used the limb darkened disk angular diameter of 20.6 mas measured by Richichi & Roccatagliata (2005) to normalize this synthetic spectrum. We checked that the existing broadband photometric measurements are consistent with the resulting SED, but we did not used them to fit the SED. Considering the true profiles of the NACO filters we used for our observations (Fig. 1), we integrated the flux from Aldebaran in W.m -2 .µm -1 in each band. We then normalized the flux integrated on the Betelgeuse images by the flux measured on Aldebaran using the same relatively narrow aperture (50 to 70 mas in radius, slightly increasing with wavelength). This approach allows us to integrate the full flux from the resolved Betelgeuse disk, while avoiding the adaptive optics halo. We did not consider the airmass difference between the two stars and the negligible infrared background.

The resulting photometric measurements are shown in Fig. 5 for the ten narrow-band flters. In this figure, we plotted the synthetic spectrum of a normal supergiant with T eff = 3 500 K, solar metallicity, V T = 2 km.s -1 and log g = 0.0 from the [START_REF] Castelli | Proc. IAU Symposium 210, Modelling of Stellar Atmospheres[END_REF] library. We normalized the expected flux using the 43.7 mas angular diameter measured by Perrin et al. (2004).

The overall agreement with our photometry is satisfactory, but a strong absorption feature around 1.09-1.24 µm is present in our data and absent from the model spectrum. We also plotted in Fig. 5 a calibrated spectrum of Betelgeuse (Cushing et al. 2005) taken from the IRTF database4 , scaled by a factor 1.2 to match the brightness of Betelgeuse at the epoch of our observations. This spectrum also shows the absorption feature that is particularly strong at 1.09 µm.

Image deconvolution

We deconvolved the images of Betelgeuse using the PSF images (Sect. 3.1) as the dirty beams and the Lucy-Richardson (L-R) algorithm. The images corresponding to the three imaging epochs of Betelgeuse were processed separately with the relevant PSF calibrator observation. The resulting average deconvolved images in the ten NACO narrow-band filters are presented in Fig. 6. We stopped the L-R deconvolution after 30 iterations, as the resulting images in the K band (2.12 and 2.17 µm) show a full width at half maximum of the photosphere of Betelgeuse compatible with the LD size measured by interferometry (Fig. 6, left column), while no significant deconvolution artefact is present.

In this process, we took into account the difference in rotational position of the pupil of the telescope between the two stars. We rotated the images of Aldebaran obtained in each filter to precisely match the telescope pupil position angle (hereafter abridged as "PA") of the corresponding images of Betelgeuse. In other words, the secondary mirror support "spikes" were brought to the same PA in Aldebaran's and Betelgeuse's images. This allowed us to avoid artefacts in the deconvolution of the Betelgeuse images. Although the presence of the spikes has a noticeable effect on the deconvolved images, we checked that the straight deconvolution of the Betelgeuse images by the Aldebaran images (without rotation) gives very similar results. Photometry on deconvolved images is a tricky task, but McNeil & Moody (2005) showed that the Lucy-Richardson algorithm preserves the photometric accuracy relatively well compared to other classical algorithms. We therefore normalized photometrically the deconvolved images presented in Fig. 6 (407 × 407 mas field of view) by the total flux of Betelgeuse in each filter (in W.m -2 .µm -1 ) as estimated in Sect. 3.2, and by the size of the resampled image pixel (3.315 mas/pix). The resulting physical unit of the images is therefore W.m -2 .µm -1 .sr -1 (surface brightness).

Circumstellar envelope shape

In the deconvolved images presented in Fig. 6, we observe the photosphere of Betelgeuse as a disk whose size appears mostly constant with wavelength (left columns). We also detect an extension of its envelope in the southwestern quadrant, at a PA of 200 ± 20 degrees (counted positively east of North), that we will subsequently refer to as a "plume". The bright part of this plume extends up to a minimum radius of six times the photosphere of the star (≈ 90 -100 mas), and is already visible in the non-deconvolved images of Betelgeuse. At a distance of 4 R * , its surface brightness in the 1.09 µm band is ≈ 0.4% of the center of the stellar disk, and falls down to ≈ 0.1% at 6 R * (Fig. 8). Two other faint extensions appear in the northeastern (PA ≈ 75 • ) and northwestern (PA ≈ 310 • ) quadrants. Although faint, their presence at all ten wavelengths (processed independently) indicates that they are likely real features.

Photospheric and CSE spectral energy distributions

Based on the deconvolved and photometrically calibrated images, we can analyze locally the SED of the CSE of Betelgeuse. Fig. 7 shows the pseudo-slit we used for this analysis, positioned in a north-south direction to include the brightest plume. The right part of this Figure shows the average surface brightness spectrum of Betelgeuse and its CSE over the pseudo-slit. The individual spectra observed at different distances from the center of the disk of Betelgeuse along the pseudo-slit are shown in Fig. 8 (they correspond to horizontal cuts in the right part of Fig. 7). We used the RMS dispersion of the pseudo-slit measurements on the three deconvolved and photometrically calibrated images of Betelgeuse obtained in each filter to estimate the systematic uncertainty of our measurements. This gave us the error bars shown in Fig. 8. The two fainter plumes visible in the Betelgeuse images in the northeastern and northwestern quadrants have similar spectra as the southern plume, with a maximum flux in the J band between 1.08 and 1.24 µm.

Discussion

Identification of the spectral features

Both the photospheric spectrum and the spectra of the plume at different radii show a spectral feature from roughly 1.1 µm up to 1.3 µm 5 . Against the photosphere, it is seen in absorption, and without a strong background source, i.e. offset from the central star but at the location of the plume, it is observed in emission.

The near-IR spectral range of red supergiant stars such as α Ori and µ Cep is the topic of a long-standing debate: several features have alternately been attributed to CN or H 2 O (Woolf et al. 1964, Wing & Spinrad 1970, Tsuji 2000). The current consensus is that features due to photospheric CN are present, as 5 As there is a gap in the NACO spectral coverage from 1.16 to 1.22 µm (see Fig. 1), it is conceivable that we are in fact dealing with two separate features. predicted by the self-consistent hydrostatic models, and that water is detected in the poorly-understood close circumstellar environment (e.g. Jennings & Sada 1998, Verhoelst et al. 2006, Perrin et al. 2007). As we are limited by the transmission of Earth's atmosphere, our spectra do not cover the water bands which were discussed most (those with bandheads around 1.4 and 1.9 µm.)

Regarding the identification of the feature observed here, H 2 O is not a candidate since this molecule does not exhibit a clear band at this particular wavelength range. CN provides an excellent match to the bandhead, and by a careful choice of temperature (2 500 K) and column density (6×10 18 cm -2 ), it is possible to model this absorption feature with a single slab of CN, without the introduction of unobserved features at the other wavelengths of CN opacity (see the upper panel of Fig. 9 for the model fit).

Unfortunately, the same basic modelling approach does not work for the emission spectra: the prediction of the bluest feature is, in the optically thin case, far too narrow (lower panel of Fig. 9) and any combination of temperature and column density Fig. 9. Upper panel: The photospheric model spectrum (black), the model + CN slab spectrum (green), the predicted inband intensities (blue) and the observations (red). Lower panel: The emission spectrum of optically thin CN at 2 500 K in green, the observations in red at a radius of 6 R * and the predicted inband intensities in blue. Both plots use the same temperature and column density for the CN slab.

that matches this broad feature, predicts a strong emission feature in the K band starting around 1.9 µm, but is not observed. Neither clumpiness nor a difference between excitation temperature and source-function temperature or a source of continuous background emission help to solve this problem. Obscuration of the feature at 1.9 µm is in theory possible with cold water vapour (cfr. the discussion on the identification of the observed absorption features at those wavelengths), but that would require a large water column density and a very particular geometry.

An important aspect of our spectro-imaging observations is that the angular resolution decreases by a factor two in the K band compared to the J band (where the plume is most visible). In the L-R deconvolution process, the angular size of the bright photosphere is recovered properly (Fig. 6), but the detection of the plume is made more difficult by its much lower surface brightness (typically several hundred times fainter than the photosphere). We therefore cannot exclude that our marginal detection of the plume in the K band is caused by a lack of resolving power, that results in a lower sensitivity of the L-R deconvolved images to extended sources. As a consequence, as we do not observe the expected emission feature around 1.9 µm, we cannot confirm our idenfication of CN, but cannot propose a viable alternative either.

Previous high angular resolution observations

Thanks to its large angular diameter and correlatively high brightness, Betelgeuse has been an excellent target for high resolution observations using a broad palette of techniques: direct full-pupil imaging, aperture masking, speckle interferometry and long-baseline interferometry. This short summary of the existing observational reports is in no way complete, but is intended to give an overview of the detections of asymmetries and surface features on Betelgeuse. Such asymmetries could point at a link between the photosphere and the envelope of the star as observed in our NACO images.

Early observations in the visible using the rotation shearing interferometer (Roddier & Roddier 1983, 1985;Roddier et al. 1986) revealed that in November 1980 and February 1982, asymmetric emission was present around Betelgeuse, and they proposed that it was caused by scattering in a dusty layer close to the star. The presence of two companions was suggested by Karovska et al. (1986), but our NACO observations do not show any point-like source down to a magnitude difference 5 in the JHK band over the field of view of our observations (0.4×0.4 ′′ ). Observations in the visible obtained at the same epoch (February 1981) with the Mayall telescope by Christou et al. (1988) show the presence of PSF-subtraction residuals in the southwestern quadrant with a typical radius of the order of 100 mas, a size and position comparable to that of the plume detected with NACO. [START_REF] Hebden | [END_REF] observed Betelgeuse in the Hα line with the MMT in speckle interferometry at two epochs: December 1983 and November 1985. They identify mostly no asymmetry at both epochs up to a radius of ≈ 95 mas, while slightly larger extensions appear in the east-west and north-south directions at larger radii. In February 1989, Buscher et al. (1990) observed Betelgeuse using the WHT equipped with a non-redundant aperture mask and three narrow-band filters in the visible around λ = 0.7 µm. These authors detected an asymmetry on the surface of the star, that they interpret as a hot spot located at a PA of ≈ 280 • , with a relative flux contribution of 10-15%. Using the same instrument and technique, Wilson et al. (1992) re-observed Betelgeuse in January 1991, and also detected asymmetries on the disk of the star and a northeast-southwest elongation different from to the 1989 epoch. Gilliland & Dupree (1996) resolved the chromosphere of Betelgeuse at ultraviolet wavelengths using the Hubble Space Telescope equipped with the Faint Object Camera (FOC), and also obtained resolved spectra with the Goddard High Resolution Spectrograph (GHRS). These data were obtained in March 1995 and clearly show the resolved chromosphere in UV continuum, extending to more than twice the radius of the visible photosphere. The deconvolved FOC images consistently show an unresolved bright region located in the southwestern quadrant of the disk, where we observed the plume with NACO. Gilliland & Dupree (1996) estimate that the spot temperature is 200 K hotter than the surrounding chromosphere. The chromosphere observed at λ = 255 nm with HST/FOC is on average much hotter (≈ 5000 K) than the overall effective temperature of the star (≈ 3600 K, Perrin et al. 2004).

In October 1995, Burns et al. (1997) used the COAST interferometer at λ = 0.80 µm to observe Betelgeuse and found no asymmetry of the stellar disk itself, but their reconstructed image (their Fig. 4) shows an extension along the northeastsouthwest direction (PA ≈ 200 • ), aligned with the NACO plume. Further observations obtained in November 1997 with COAST and the WHT telescope at visible and near infrared wavelengths (0.70, 0.91 and 1.29 µm) by Young et al. (2000) also showed an asymmetry of the light distribution on the disk, that tends to decrease as wavelength increases. Interestingly, their 1.29 µm observations do not show any asymmetry, while our observations with NACO at 1.28 µm still show the presence of the southwestern plume. Tuthill et al. (1997) also used the WHT telescope in aperture masking mode around λ = 0.70 µm between January 1992 and December 1993. Depending on the epoch, they identified two to three hot spots with changing positions, and relative flux contributions of up to 20%, on the surface of Betelgeuse.

Betelgeuse was observed in early 1998 by Hinz et al. (1998) using nulling interferometry in the mid-infrared (λ = 10 µm) with two mirrors of the Multiple Mirror Telescope separated by a baseline of 5 meters. This interferometric technique allows to cancel out the light from the unresolved star, leaving as a residual the emission from the circumstellar environment. Their resulting nulled image (their Fig. 3) shows a very extended (≈ 5 ′′ in diameter) asymmetric dusty envelope, with extensions in the southwestern, northwestern and northeastern directions, similar to the positions of the three plumes we observed with NACO. Smith et al. (2009) recently found a mildly clumpy, spherical shell around Betelgeuse in CO emission, from spatially resolved spectroscopy at λ ≈ 4.6 µm. Some interesting structures were detected ≈ 1 ′′ west of the star, but their identification to a mass loss event or envelope asymmetry is unclear.

The single-mode two-telescope interferometric observations obtained by Perrin et al. (2004) and Chagnon et al. (2002) respectively in the K and L bands provided high accuracy photospheric angular diameters of Betelgeuse (θ LD K = 43.7 ± 0.10 mas, θ UD L = 42.17 ± 0.05 mas). However, due to the absence of closure phase measurements, it was not possible to check for the presence of asymmetry on the stellar disk. New observations of Betelgeuse with the three-telescope beam combiner IONIC installed at the IOTA interferometer were obtained in October 2005 and March 2006 by [START_REF] Haubois | Proc. SF2A conf[END_REF] in the infrared H band. They point at the presence of one very faint hot spot (relative flux 0.5%) in the northwestern quadrant of Betelgeuse (through model-fitting), or alternatively at two hot spots close to the center of the disk (through image reconstruction). Spectrally dispersed interferometric observations of Betelgeuse in the K band (obtained in January 2008) are reported by Ohnaka et al. (2009). These authors find a centrally symmetric emission on the apparent disk of the star in the continuum, but an inhomogeneous velocity field in the 2.303 µm CO line, with a North-South preferential direction. This direction appears consistent with the preferential axis of the main plume we detected with NACO.

At thermal infrared wavelengths (λ = 11.15 µm), Tatebe et al. (2007) identified an asymmetry on the disk of Betelgeuse from three-telescope interferometric observations with the ISI instrument obtained in November and December 2006. They interpreted this asymmetry as the signature of a hot spot located on the southern limb of the star at a PA of 170 -180 • . This position is compatible with the PA of the plume we detect with NACO. Perrin et al. (2007) observed Betelgeuse with the two-telescope instrument VLTI/MIDI in the N band (λ = 7.5 -13.5 µm), de-tecting a molecular layer (the MOLsphere) containing H 2 O, SiO and Al 2 O 3 at a radius of ≈ 1.4 R * from the star.

Rotation

The rotation of a supergiant star as Betelgeuse is expected to be very slow, but it can have non-negligible effects on the stellar atmosphere structure. Uitenbroek et al. (1998) presented a detailed analysis of the HST/GHRS spectra first reported by Gilliland & Dupree (1996). They showed that chromospheric emission in the UV emission lines of Mg II extends to ≈ 270 mas from the center of the star, twice further than the UV continuum. In addition, thanks to a scanning of the instrument over the stellar disk, they could determine the projected rotational velocity of (v sin i = 5 km/s), the PA of its polar axis (PA ≈ 235 • ) and the inclination of the polar axis on the line of sight of (i = 20 • ). From spectroscopy and theoretical considerations, Gray (2000) obtained a similar rotational velocity. This position of the polar cap is compatible with the hot spot visible in the FOC images and the PA of the southwestern plume observed in our NACO images (PA= 200 ± 20 • )

Origin of the southwestern plume

The presence of the bright southwestern plume detected with NACO implies that the spherical symmetry observed on the infrared photosphere of Betelgeuse is not preserved in its close environment. Among the plausible mechanisms to explain its presence, two appear as promising: convection and rotation. Following the proposal by [START_REF] Schwarzschild | [END_REF] and the simulation results by Freytag & Höfner (2008), a first hypothesis is that the very large convection cells enhance locally the massloss from the star, creating the observed plumes of molecules. However, although many studies have shown that the photosphere of Betelgeuse presents a significant degree of asymmetry, variable with time, the link with the chromosphere and the envelope remains to be established. A second hypothesis is that mass loss may be enhanced above the polar cap of the star, creating the southwestern plume of Betelgeuse. This possibility is supported by the HST imaging and spectroscopy by Uitenbroek et al. (1998), as it reveals the presence of a hot spot that could be located at the pole of the star. An interaction between this hot spot and the extended chromosphere of the star could result in a preferential production of molecules above the pole and the creation of a plume.

Conclusion

From diffraction-limited observations in ten narrow-band filters spread over the JHK spectral domain, we find evidence for the presence of a complex CSE around Betelgeuse. The most conspicuous feature of this envelope is a bright plume extending in the southwestern quadrant up to a minimum distance of six photospheric radii. Its SED in emission appears complementary to a strong absorption feature present in the J band (1.08-1.24 µm) on the photospheric spectrum. This points at the presence of a molecular layer, possibly CN, although the formal identification of the spectral signature of this molecule in the K band is not established. Two hypotheses could explain the formation of the plume: convection causing an enhanced mass loss above a large upwards moving convective cell, or rotation, through the presence of a hot spot at the location of the polar cap of the star. Further spectro-imaging observations resolving the photosphere of Betelgeuse (e.g. by interferometry) should allow to discriminate between these possibilities, if a link can be established between photospheric hot spots (characteristic of convection) and the detected plume. For each wavelength, the left image is the full-dynamic linear color scale image normalized to the maximum and minimum values. The 10 contours plotted in red correspond to flux levels of 1 000, 2 000, 4 000, 8 000, 16 000, 32 000, 64 000, 128 000 and 256 000 W.m -2 .µm -1 .sr -1 . The right image uses a square root color scale with narrower color cuts (including 95% of the pixels) to emphasize the faint extensions of the CSE. The green circle corresponds to the 43.7 mas photospheric angular diameter of the star as measured by Perrin et al. (2004) from K band interferometric observations.

Fig. 1 .

 1 Fig. 1. Transmissions of the NACO narrow-band filters (black curves) and of the atmosphere (from Lord 1992, thin grey curve).

Fig. 3 .

 3 Fig. 3. Average frames of the NACO cubes #001 (1.04 µm), #007 (1.64 µm) and #010 (2.17 µm) obtained on the PSF calibrator Aldebaran without centering or frame selection (left) and with centering and 10% best frame selection (right). The white disk in the 1.04 µm images represents the angular size of Aldebaran's photosphere (20 mas). The color scales are proportional to the square root of the flux, and are identical for the two images at each wavelength.

Fig. 4 .

 4 Fig. 4. Upper panel: Relative Strehl ratio of the recentered and co-added images compared to the classic mode of NACO, without frame selection (dashed curve) and keeping only the 10% best peak flux frames (thin solid curve). Lower panel: Full width at half maximum of the NACO images of Aldebaran for the classic mode (black curve), frame centering without selection (dashed curve) and 10% best frames selection and centering (thin solid curve), compared to the theoretical diffraction pattern for a UT observing a uniform disk of 20 mas (grey curve).

Fig. 5 .

 5 Fig. 5. Calibrated photometry of Betelgeuse. The dots represent our photometry in the 10 narrow-band filters of NACO. The thin black curve represents the model spectrum of a T eff = 3 500 K normal supergiant (Castelli & Kurucz 2003) normalized to a limb darkened disk angular diameter of 43.7 mas. The grey curve is a spectrum of Betelgeuse from the IRTF database (Cushing et al. 2005), scaled by a factor 1.2.

Fig. 7 .

 7 Fig. 7. Left panel: Positions of the 3D pseudo-slit, plotted over the image of Betelgeuse at 1.04 µm. The color scale is logarithmic to better show the extensions of the CSE. The central circle represents the photosphere of the star. The graduations from 0 to 6 are distances from the photosphere center in stellar radii for the spectra shown in Fig. 8. Right panel: Reconstructed long-slit spectrum of Betelgeuse (linear interpolation from the pseudo-slit cut in the Betelgeuse image cube). The color scale is logarithmic. The vertical lines indicate the central wavelengths of the NACO narrow-band filters (see also Fig. 1), and the 43.7 mas photospheric angular diameter is shown with green lines.

Fig. 8 .

 8 Fig. 8. Surface brightness spectra of the close environment of Betelgeuse, averaged over the pseudo-slit width (44 mas). The vertical scales are in W.m -2 .µm -1 .sr -1 , and the horizontal scale is the wavelength in µm. The radius (in R * ) from the center of the photosphere is indicated in the upper right corner of each plot.

Fig. 2 .

 2 Fig. 2. Average NACO images of Betelgeuse (left image of each column) and the PSF calibrator Aldebaran (right image of each column) in 10 narrow-band filters from 1.04 to 2.17 µm. The color scale is a function of the square root of the flux, normalized to the maximum and minimum value of each image, in order to emphasize the fainter parts of the images. The white disks in the lower right corners of the 1.04 µm images represent the photospheric angular diameters of Betelgeuse and Aldebaran.

Fig. 6 .

 6 Fig.6. Lucy-Richardson deconvolved images of Betelgeuse (30 iterations). For each wavelength, the left image is the full-dynamic linear color scale image normalized to the maximum and minimum values. The 10 contours plotted in red correspond to flux levels of 1 000, 2 000, 4 000, 8 000, 16 000, 32 000, 64 000, 128 000 and 256 000 W.m -2 .µm -1 .sr -1 . The right image uses a square root color scale with narrower color cuts (including 95% of the pixels) to emphasize the faint extensions of the CSE. The green circle corresponds to the 43.7 mas photospheric angular diameter of the star as measured byPerrin et al. (2004) from K band interferometric observations.

http://www.eso.org/sci/facilities/paranal/instruments/naco/

http://yorick.sourceforge.net/

IRAF is distributed by the NOAO, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.

http://irtfweb.ifa.hawaii.edu/∼spex/IRTF Spectral Library/

Acknowledgements. We thank Dr. C. Lidman for his help in the preparation of these observations, and the ESO Paranal staff for their perfect execution in visitor mode at the telescope. STR acknowledges partial support by NASA grant NNH09AK731. This work also received the support of PHASE, the high angular resolution partnership between ONERA, Observatoire de Paris, CNRS and University Denis Diderot Paris 7. This research took advantage of the SIMBAD and VIZIER databases at the CDS, Strasbourg (France), and NASA's Astrophysics Data System Bibliographic Services.

The "•" symbol indicates that the neutral density filter was used ("•" indicates it was not used). b DIT is the detector integration time. c σ is the DIMM seeing at Paranal in the visible and AM is the airmass.


Konstantia Nathanael

Sibo Cheng

Nina M Kovalchuk

Rossella Arcucci

Mark J H Simmons

Optimization of microfluidic synthesis of silver nanoparticles: a generic approach using machine learning

Keywords: Reaction kinetics, microfluidic synthesis, silver nanoparticles, decision tree, machine learning Highlights AgNPssilver nanoparticles, AIartificial intelligence, BO -Bayesian optimization, DCAdicarboxyacetone, PTFEpolytetrafluoroethylene, PVApolyvinyl alcohol

come

Introduction

The precise and controlled fabrication of silver nanoparticles (AgNPs) is required for a range of applications [1]. Optimisation of the chemical synthesis process is a very labour-intensive, costly and time-consuming task because it includes multiple reagents and different independent experimental conditions including the types and concentrations of reactants, temperature, reactor design and mixing conditions. For example, the size of AgNPs strongly depends on both the type and concentration of stabilizing agent [2], [START_REF] Patel | Role of stabilizing agents in the formation of stable silver nanoparticles in aqueous solution: Characterization and stability study[END_REF], [START_REF] Tejamaya | Stability of Citrate, PVP, and PEG Coated Silver Nanoparticles in Ecotoxicology Media[END_REF] and often depends on pH of the solution [START_REF] Anigol | Effect of Concentration and pH on the Size of Silver Nanoparticles Synthesized by Green Chemistry[END_REF], [START_REF] Dong | Shape Control of Silver Nanoparticles by Stepwise Citrate Reduction[END_REF], [START_REF] Qin | Size control over spherical silver nanoparticles by ascorbic acid reduction[END_REF]. Ionic stabilizers such as citrate used in this study produce a charged layer around particles which inhibit agglomeration.

Several studies showed that both decreasing [START_REF] Guardia | Controlled Synthesis of Iron Oxide Nanoparticles over a Wide Size Range[END_REF] and increasing [START_REF] Shevchenko | Study of Nucleation and Growth in the Organometallic Synthesis of Magnetic Alloy Nanocrystals: The Role of Nucleation Rate in Size Control of CoPt3 Nanocrystals[END_REF] particles sizes can be obtained by varying the amount of stabilizer. Das, Bandyopadhyay [START_REF] Das | Effect of stabilizer concentration on the size of silver nanoparticles synthesized through chemical route[END_REF] and Henglein and Giersig [START_REF] Henglein | Formation of Colloidal Silver Nanoparticles: Capping Action of Citrate[END_REF] found there is an optimum concentration of trisodium citrate in AgNP synthesis. They observed that smaller particles with a narrow size distribution can be formed for a specific range of citrate concentrations and larger particles at concentrations higher than the optimal, due to the high ionic strength of the solutions.

It was observed that the reducing ability of the completely hydrolysed citrate species and dicarboxyacetone (a by-product of citrate reaction with silver) is higher at pH 12, causing a fast reduction rate of the precursor and therefore smaller particles [START_REF] Marciniak | The Effect of pH on the Size of Silver Nanoparticles Obtained in the Reduction Reaction with Citric and Malic Acids[END_REF]. However, it was also found in Liu, Kozlovskaya [START_REF] Liu | Encapsulation of anticancer drug by hydrogen-bonded multilayers of tannic acid[END_REF] that alkaline pH in the presence of high salt concentrations can lead to thicker tannic acid-poly(N-vinylpyrrolidone) (TA-PVPON) multilayers as a result of screening of the negatively charged tannic acid. Furthermore, the hydrodynamic conditions play a significant role in the nucleation and growth steps for the synthesis of nanoparticles [START_REF] Baber | An engineering approach to synthesis of gold and silver nanoparticles by controlling hydrodynamics and mixing based on a coaxial flow reactor[END_REF], [START_REF] Pal | Tuning of Particle Size in a Helical Coil Reactor[END_REF]. Relevant works are represented in Khan, Günther [START_REF] Khan | Microfluidic Synthesis of Colloidal Silica[END_REF] and Wu, De Varine Bohan [START_REF] Wu | Synthesis of narrow sized silver nanoparticles in the absence of capping ligands in helical microreactors[END_REF]. The former demonstrated that the sizes and size distributions of colloidal silica particles can be tuned in laminar flow reactors and segmented flow reactors by varying their linear flow velocity and mean residence time. The latter showed that decreasing the helix diameter increases the mixing in helical reactor (due to generation of Dean vortices) which then leads to a controlled size distribution in continuous synthesis of AgNPs in the absence of capping agents. A similar effect was observed at increasing flow rates when the Dean number (𝐷𝑒) was above 5.

Microfluidic technology can improve the optimisation process by providing better control over the reactions and reduced reagent consumption [START_REF] Kulkarni | Microfluidic devices for synthesizing nanomaterials-a review[END_REF]. Machine learning combined with microfluidic synthesis offers a promising approach to tackle the repeated and extensive experimental tests and accelerate the development of efficient protocols [START_REF] Chen | Intelligent control of nanoparticle synthesis on microfluidic chips with machine learning[END_REF].

Recently, artificial intelligence (AI) has become a valuable aid for many scientific and engineering fields, including nanotechnology, due to its contribution in data acquisition and processing improvements [START_REF] Hrvat | Artificial Intelligence in Nanotechnology: Recent Trends, Challenges and Future Perspectives[END_REF], [START_REF] Zhang | Nanotechnology and artificial intelligence to enable sustainable and precision agriculture[END_REF], [START_REF] Adir | Integrating Artificial Intelligence and Nanotechnology for Precision Cancer Medicine[END_REF]. Machine learning (ML) algorithms which are a subset of AI have been employed specifically in the synthesis of AgNPs to predict their characteristics [START_REF] Sun | Machine Learning for Silver Nanoparticle Electron Transfer Property Prediction[END_REF], [START_REF] Khayati | A predictive model on size of silver nanoparticles prepared by green synthesis method using hybrid artificial neural network-particle swarm optimization algorithm[END_REF], [START_REF] Findlay | Machine learning provides predictive analysis into silver nanoparticle protein corona formation from physicochemical properties[END_REF] or perform other kinds of decision making under uncertainty [START_REF] Liu | Cytotoxicity of phytosynthesized silver nanoparticles: A meta-analysis by machine learning algorithms[END_REF]. Mekki-Berrada, Ren [START_REF] Mekki-Berrada | Two-step machine learning enables optimized nanoparticle synthesis[END_REF] combined a deep neural network (DNN) and Gaussian process-based Bayesian optimization (BO) to synthesize AgNPs with a desired absorbance spectrum in a droplet microfluidic device. They used the flow rates of silver seeds (𝑄 𝑠𝑒𝑒𝑑 ), of silver nitrate, SN (𝑄 𝐴𝑔𝑁𝑂 3 ), of trisodium citrate, TC (𝑄 𝑇𝐶 ) and polyvinyl alcohol, PVA (𝑄 𝑃𝑉𝐴 ) as input parameters and efficiently predicted the desired plasmon resonance for the reaction synthesis. Sattari and Khayati [START_REF] Sattari | Prediction of the size of silver nanoparticles prepared via green synthesis: A gene expression programming approach[END_REF] applied a Gene Expression Programming (GEP) to predict the size of AgNPs prepared by a green synthesis route. The proposed predictive model with coefficient of determination , 𝑅 2 = 0.9961, mean absolute error, MAE = 0.2545, and root relative square error, RSME= 0.0668 showed that the initial concentration of silver precursor and plant extract were the most influential parameters for the final particle size. Other relevant work used to predict the final particle size of AgNPs is shown in Shabanzadeh, Senu [START_REF] Shabanzadeh | Artificial intelligence in numerical modeling of silver nanoparticles prepared in montmorillonite interlayer space[END_REF].

Decision trees are an established keystone in machine learning literature [START_REF] Myles | An introduction to decision tree modeling[END_REF] as a powerful tool for prediction, interpretation and data manipulation. More specifically, a decision tree is a supervised learning algorithm which sorts a population into segments. This approach looks like an inverted tree with a root node, some internal nodes and leaf nodes. The root node is always on the top of the tree structure and represents the variable that best splits the data, whilst the leaf nodes at the bottom display the final outcome of a combination of decisions. A decision pathway is presented by a line connecting the root node with one of the leaf nodes [START_REF] Song | Decision tree methods: applications for classification and prediction[END_REF]. Compared with deep learning approaches [START_REF] Lecun | Deep learning[END_REF], the advantages of a tree-based algorithm are the good interpretability of the constructed models and the low computational cost. It can identify, in a linear way, information and relationships that can be used to design future experiments and analyse data [START_REF] Myles | An introduction to decision tree modeling[END_REF]. A limitation is that it can display overfitting or underfitting within a small data set [START_REF] Song | Decision tree methods: applications for classification and prediction[END_REF].

The literature review revealed that previous ML studies for the synthesis of AgNPs focused on input parameters which are related to specific chemical compositions and process conditions. This means that a model which is derived for a certain experimental configuration cannot be applied for another set of reagents or for a different reactor design. The aim of this study is to develop a generalised workflow to design experiments for fast prediction of the size of AgNPs synthesized in microfluidic systems. It is suggested that the chemistry of the process can be described in the terms of nucleation and growth constants, whereas the hydrodynamics of the process can be accounted for by Reynolds number, Re, and the ratio of Dean number to Reynolds number 𝐷𝑒 𝑅𝑒 ⁄ . Storage temperature is included to enable consideration of particle stability after collection.

A machine learning-guided design of experiments based on the decision tree method has been applied for the production of AgNPs in a continuous flow microreactor. The proposed strategy, which used optimal and flexible experiment designs based on uncertainty analysis, can provide more flexibility compared to existing approaches such as design of experiments [START_REF] Sethuramiah | Chapter 6 -Statistics and Experimental Design in Perspective[END_REF], since the exact values of the initial parameters are not fixed and the training features could be parameters which have been calculated based on physical parameters, such as nucleation and growth constants. The nucleation and growth constants were derived from an independent set of experiments carried out in a beaker using the Finke-Watzky (F-W) two-step mechanism [START_REF] Watzky | Transition Metal Nanocluster Formation Kinetic and Mechanistic Studies. A New Mechanism When Hydrogen Is the Reductant: Slow, Continuous Nucleation and Fast Autocatalytic Surface Growth[END_REF] which has been applied previously to describe well a wide range of kinetic processes such as protein aggregation [START_REF] Iashchishyn | Finke-Watzky Two-Step Nucleation-Autocatalysis Model of S100A9 Amyloid Formation: Protein Misfolding as "Nucleation" Event[END_REF], [START_REF] Ohgita | Mechanisms of enhanced aggregation and fibril formation of Parkinson's disease-related variants of α-synuclein[END_REF] and metal nanoparticle formation including silver [START_REF] Sandoe | Experimental probes of silver metal nanoparticle formation kinetics: Comparing indirect versus more direct methods[END_REF], gold [START_REF] Pestovsky | Formation of Aggregate-Free Gold Nanoparticles in the Cyclodextrin-Tetrachloroaurate System Follows Finke&ndash[END_REF], palladium [START_REF] Wojnicki | Kinetic studies of nucleation and growth of palladium nanoparticles[END_REF] and rhodium [START_REF] Yao | Insights into the Formation Mechanism of Rhodium Nanocubes[END_REF]. Only one set of chemicals was used in this study providing several pairs of nucleation and growth constants. These kinetic constants were used as input parameters in the decision tree to build a ML model which combines parameters that can be applied in the future for different sets of chemicals. By applying this strategy, a general model for the synthesis of AgNPs can be obtained in a more direct way by providing the region of interest for further experiments. The suggested approach can also enable a considerable reduction of ecological impact of process optimization by using microfluidics, which manipulates small samples under well-controlled conditions resulting in a considerable reduction of materials and energy consumption.

2.Materials and methods

AgNPs synthesis and characterisation

The silver nanoparticles were synthesized in a T-junction microfluidic device composed of 0.5 mm inner diameter polytetrafluoroethylene (PTFE) cylindrical tubes (Cole-Parmer). Two inlets (length of 0.6 m each) were used to supply silver nitrate, SN (0.92 mM) and the mixture of tannic acid, TA (0.123 mM) / trisodium citrate, TC (1.91 mM -3.82 mM.) respectively and an outlet tube of 2 m length was used as the reaction channel (Fig. 1). The inlet flow rates were regulated by syringe pumps (World precision instrument-AL-4000) equipped with 5 mL syringes (Fisher), while the three tubes were connected with a Tee tubing junction (0.020" (0.5 mm) Thru-Hole) from Upchurch Scientific. The outlet channel has the same cross section as the inlet channels, therefore the superficial velocity in the in the outlet channel is twice that observed in each of the two inlet channels. The concentrations of reagents were determined based on a previous study of Kašpar, Koyuncu [START_REF] Kašpar | Governing factors for preparation of silver nanoparticles using dropletbased microfluidic device[END_REF] with modifications in the concentrations of reducing agent (TA) and stabilizing agent (TC). The solution of TA and TC had a pH of 7 or 12. The reagents including silver nitrate, 99+% (𝐴𝑔𝑁𝑂 3 ), tannic acid (𝐶 7 𝐻 52 𝑂 46 ), trisodium citrate dihydrate (𝑁𝑎 3 𝐶 6 𝐻 5 𝑂 7 • 2𝐻 2 𝑂) and sodium hydroxide (𝑁𝑎𝑂𝐻) were purchased from Alfa Aesar.

A

The outlet tube was either straight or coiled onto a 3D-printed helical shape devices with diameters of 3 mm and 5 mm to improve the mixing of the reagents. The particle size was measured using dynamic light scattering (DLS) by Zetasizer Nano series (Malvern) and confirmed by transmission electron microscopy (TEM) (JEOL JEM-1400). The concentration of silver nitrate at different flow rates was measured to investigate the completeness of the reaction. A silver electrode from EDT direction connected to mV meter from Mettler Toledo-FP20 was used.

The nucleation, 𝑘 1 , and growth, 𝑘 2 , rate constants were derived from the absorbance intensity at 400 nm, characteristic for AgNPs using F-W two step mechanism [START_REF] Watzky | Transition Metal Nanocluster Formation Kinetic and Mechanistic Studies. A New Mechanism When Hydrogen Is the Reductant: Slow, Continuous Nucleation and Fast Autocatalytic Surface Growth[END_REF] which combines a homogenous nucleation reaction and an autocatalytic growth process in which nuclei and growing particles play the role of auto-catalysts. The time dependence of the absorbance of AgNPs (kinetic curve) for different experimental conditions was measured by UV-vis spectrophotometer (Jenway-6300). Reactions were carried out in a beaker and mixing was set at a point where further increase of mixing intensity does not affect the kinetic curve. Thus, AgNPs synthesis was independent of mass-transfer and the rate was determined by true chemical kinetics. To quantify the amount of silver in AgNPs at any time during synthesis, the UV-vis spectra was recorded every 5 min. All measurements were performed in triplicate.

The synthesis of AgNPs involves three different steps [START_REF] Thanh | Mechanisms of Nucleation and Growth of Nanoparticles in Solution[END_REF] including the reduction of silver ions to silver atoms, the formation of silver nuclei and their subsequent growth to generate AgNPs (see Fig. 2).

In this study, the reaction was performed through the reduction of SN in the presence of TA which is a weak reducing agent and TC which is both reducing and stabilizing agent. As described in Eq. ( 1) and (2), nuclei represented as 𝐵 are produced uniformly from the precursor represented as 𝐴. This reaction follows an autocatalytic growth of nuclei leading to the formation of 𝐵 particles (see Fig. 2).

𝑁𝑢𝑐𝑙𝑒𝑎𝑡𝑖𝑜𝑛: 𝐴

𝑘 1 → B ; 𝑛 0 = 𝑘 1 [𝐴],
(1)

𝐺𝑟𝑜𝑤𝑡ℎ: 𝐴 + 𝐵 𝑘 2 → 2𝐵; 𝐺 = 𝑘 2 [𝐴][𝐵], (2)
where, 𝑛 0 and 𝐺 are the nucleation and growth rates, 𝑘 1 and 𝑘 2 are the nucleation and growth kinetic constants and [𝐴] and [𝐵] are the molar concentrations of the precursor and silver in the nuclei/particles.

The overall reaction rate can be expressed as shown:

-

𝑑[𝐴] 𝑑𝑡 = 𝑑[𝐵] 𝑑𝑡 = 𝑘 1 [𝐴] + 𝑘 2 [𝐴][𝐵] (3)
The F-W mechanism can be expressed though the integrated form (Eq. 4) which was also used to fit all experimental kinetic data sets [START_REF] Sandoe | Experimental probes of silver metal nanoparticle formation kinetics: Comparing indirect versus more direct methods[END_REF]. It was considered that the concentration of material inside the particles changes with time 𝑡 as

[𝐵] 𝑡 = [𝐴] 0 -[𝐴] 𝑡
, where [𝐴] 0 is the concentration of the precursor solution at t = 0 and [𝐴] 𝑡 the concentration of precursor at time 𝑡.

[𝐵](𝑡) = [𝐴] 0 -

𝑘 1 𝑘 2 + [𝐴] 0 1 + 𝑘 1 𝑘 2 [𝐴] 0 * 𝑒𝑥𝑝(𝑘 1 + 𝑘 2 [𝐴] 0 ) 𝑡 . ( 4
)
For curve fitting procedures, the value of [𝐴] was estimated using Eq. ( 5), where 𝑎 is 𝐵 𝑡 𝐵 ∞ ⁄ and 𝐵 𝑡 and 𝐵 ∞ are the maximum absorbance at 𝑡 and ∞, respectively [START_REF] Amirjani | Modified Finke-Watzky mechanisms for the two-step nucleation and growth of silver nanoparticles[END_REF].

[𝐴] 𝑜 -[𝐴]

[𝐴] = 𝑎 1-𝑎 [START_REF] Anigol | Effect of Concentration and pH on the Size of Silver Nanoparticles Synthesized by Green Chemistry[END_REF] The nucleation 𝑘 1 and growth 𝑘 2 constants were estimated using the linearized form below under assumptions that nucleation happens more slowly than growth (i.e. 𝑘 1 ≪ 𝑘 2 [𝐴]) and [𝐴] < [𝐴] 0 [START_REF] Amirjani | Modified Finke-Watzky mechanisms for the two-step nucleation and growth of silver nanoparticles[END_REF].

𝑙𝑛 ( [𝐴] 0 -[𝐴] [𝐴] ) = 𝑙𝑛 ( 𝑘 1 𝑘 2 [𝐴 0 ] ) + 𝑘 2 [𝐴] 0 𝑡, (6)
By plotting, ln (𝑎 1 -𝑎) ⁄ versus 𝑡, a straight line with a slope and intercept were obtained as shown in 1.91, 2.87 and 3.82 mM. This provided six different pairs of nucleation and growth constants to be used in the ML models.

Machine learning modelling

The machine learning analysis was performed using Python programming language with the Scikitlearn package [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. Tree-based algorithms, including Decision Tree (DT), Random Forest (RF), and Extreme Gradient Boosting (XGBoost) were applied to test the prediction performance for the sizes of AgNPs. These methods have been extensively used in a wide range of engineering problems [START_REF] Natekin | Gradient boosting machines, a tutorial[END_REF], [START_REF] Cheng | Parameter Flexible Wildfire Prediction Using Machine Learning Techniques: Forward and Inverse Modelling[END_REF], [START_REF] Gong | An efficient digital twin based on machine learning SVD autoencoder and generalised latent assimilation for nuclear reactor physics[END_REF]. After implementation, all the proposed models were compared based on the accuracy of their predictions. For each experiment, the average size was recorded. Each dataset was divided into five inputs and one predicted output. The training features (inputs) and their investigated range are summarized in Table 1.

𝑅𝑒 = 𝑝𝑢𝑑 𝜇 (7)
𝐷𝑒 = 𝑅𝑒√ 𝑑 2𝑅 𝑐 ⁄ (8)
Eq. ( 7) and ( 8) define Reynolds, 𝑅𝑒, and Dean, 𝐷𝑒, numbers respectively, where 𝑑 is the inner diameter of the microchannel, 𝑅 𝑐 is the radius of helix curvature, 𝜌 is the fluid density, 𝑢 is the velocity of the fluid, and 𝜇 is the viscosity of the fluid.

Table 1: Input parameters and their investigated range.

Feature

Range investigated

Nucleation constant (min -1 ) 0.0011 -0.1

Growth constant (M -1 min -1 ) 13.66-77.97

Storage temperature (℃) 0 -20

Dean number/Reynolds number 0 -0.41

Reynolds number 0.0849-16.96

The effect of the chosen input parameters on the size of AgNPs was investigated through the Pearson correlation [START_REF] Navlani | Python Data Analysis[END_REF] shown in Eq.( 9) [START_REF] Profillidis | Chapter 5 -Statistical Methods for Transport Demand Modeling[END_REF].

𝑝 𝑋𝑌 = ∑(𝑥 𝑖 -𝑥̅ )(𝑦 𝑖 -𝑦 ̅) √∑(𝑥 𝑖 -𝑥̅ ) 2 • ∑(𝑦 𝑖 -𝑦 ̅) 2 • (9)
In Eq. ( 9) 𝑥 𝑖 is for the values of the independent variables, 𝑥̅ is the mean of the values of x-variable, 𝑦 𝑖 is for the values of the dependent variables and 𝑦 ̅ is the mean of the values of y-variable.

DT-guided design of experiments

The iterative process of prediction and experiment design carried out is shown in Fig. 3, with the key focus being to perform experiments where the predictive model has the highest uncertainties. In particular, DT is chosen to guide the experiments thanks to its good interpretability and computational efficiency [START_REF] Gilpin | Explaining Explanations: An Overview of Interpretability of Machine Learning[END_REF]. The first step in the proposed approach described previously was focused on obtaining the datasets required for the training of the algorithm, i.e., data collection based on experimental measurements. The second and third steps included data sorting or any other required data processing and the development of the most appropriate algorithm. For example, during these steps, the independent and dependent variables were identified, and a DT based algorithm was built. To ensure a good explainability of the DT model for the uncertainty evaluation, the maximum tree depth for this stage was set to be 4.

The fourth stage comprised carrying out additional experiments where the obtained DT algorithm has the highest mean squared error (MSE). The challenge of this step was to ensure that the developed DT model based on the designed experiments can improve the performance regarding different metrics. As shown in Fig. 3, once the experimental design is performed using the preliminary DT model, other ML classifiers such as RF and XGBoost can be applied to achieve a more accurate prediction based on the new experimental data.

Results and discussion

Kinetic constants

The kinetic curves describing time dependence of concentration of silver included in nanoparticles (Fig. 4A,4B) were analysed by F-W mechanism [START_REF] Watzky | Transition Metal Nanocluster Formation Kinetic and Mechanistic Studies. A New Mechanism When Hydrogen Is the Reductant: Slow, Continuous Nucleation and Fast Autocatalytic Surface Growth[END_REF] described in detail in Section 2.1. At pH 7, nucleation is rather slow, and the kinetic curves have a typical sigmoidal shape. The values of kinetic constants found from the linear fitting are given in Table 2. Using these values, it is easy to find that condition 𝑘 1 ≪ 𝑘 2 [𝐴] is valid for [A] >> 0.03 mM, which is only 3 % of initial concentration of SN, 0.92 mM. The obtained values of rate constants are validated by recalculating the full kinetic curve using Eq. ( 4) with these constants. Excellent agreement between the theoretical and experimental kinetic curves is shown in Fig. 4A (pH 7).

Multiple studies such as Anigol, Charantimath [START_REF] Anigol | Effect of Concentration and pH on the Size of Silver Nanoparticles Synthesized by Green Chemistry[END_REF], Jebakumar Immanuel Edison and Sethuraman [START_REF] Immanuel Edison | Electrocatalytic Reduction of Benzyl Chloride by Green Synthesized Silver Nanoparticles Using Pod Extract of Acacia nilotica[END_REF],

and Dong, Ji [START_REF] Dong | Shape Control of Silver Nanoparticles by Stepwise Citrate Reduction[END_REF] have shown that pH changes the kinetics in the formation of AgNPs, using Capparis Moonii fruit extract, Pod Extract of Acacia nilotica and citrate as reducing agents, respectively. At low pH, there is a slow reduction rate of the precursor while at high pH there is an improved reducing ability of the reducing agent that results in smaller size particles.

At pH 12, it was observed that the silver nuclei formation proceeded very fast, without a noticeable lag phase for all studied stabilizing agent concentrations, in agreement with literature data [START_REF] Iashchishyn | Finke-Watzky Two-Step Nucleation-Autocatalysis Model of S100A9 Amyloid Formation: Protein Misfolding as "Nucleation" Event[END_REF]. Only the right-hand side part of sigmoidal curve can thus be measured, which mostly accounts for particle growth. Although the linear fitting in Fig. 4D (pH 12) gives a reasonable fit, using values of kinetic constants provided in Table 2 shows that this linear fitting is valid only for [A] >> 0.5 mM for

[TC] = 3.82 mM and [A] >> 1 mM for [TC] = 1.91 mM. Using these constants with the full F-W model, Eq. ( 4), does not provide as good agreement with experimental data as for the lower value of pH, but it is still acceptable as can be seen from Fig. 4D (pH 12). Note, non-linear fitting directly using Eq. ( 6)

does not improve the fitting presented in Fig. 4B (pH 12). Thus, the values of nucleation and growth constants presented in Table 2 were used as the best possible fitting. For the investigated case, under alkaline pH, there is a full dissociation of citrate species [START_REF] Dong | Shape Control of Silver Nanoparticles by Stepwise Citrate Reduction[END_REF] and dicarboxyacetone (DCA) [START_REF] Marciniak | The Effect of pH on the Size of Silver Nanoparticles Obtained in the Reduction Reaction with Citric and Malic Acids[END_REF], generated as a by-product of the decarboxylation of citrate during silver ions reduction, which makes the chemistry of the system more complex. Therefore, it is possible that several parallel processes are involved with different nucleation and growth constants making fitting of the kinetic curves much more challenging.

Pearson Correlation and the effect of input parameters

The results of Pearson correlation calculations are presented in Fig. 5 which shows that all input variables impact the particle size significantly with a correlation larger than 0.12 in absolute value. The estimation of the Pearson correlation value was obtained using all of the data generated in this study.

The data used for this study can be found in the supporting information (SI). For example, it was observed that the variables including 𝑘 1 , 𝑘 2 and 𝐷𝑒 𝑅𝑒 ⁄ impact the particle size negatively, whilst the storage temperature and 𝑅𝑒 impact the size of particles positively. The Pearson correlation calculations showed that 𝑘 2 is less important compared to 𝑘 1 for the prediction of the size of AgNPs. In this study, the effect of kinetic constants was investigated for only one set of chemicals, so that k1 and k2 cannot be considered as completely uncoupled. Therefore, it cannot be concluded that faster growth results in smaller particles, because for the chosen set of chemicals all growth constants are of the same order of magnitude, whereas nucleation constants vary over 2 orders of magnitude. Thus, the particle size is mostly defined by the nucleation constant. Similar studies, such as the works of Mansouri and Ghader [START_REF] Mansouri | Experimental study on effect of different parameters on size and shape of triangular silver nanoparticles prepared by a simple and rapid method in aqueous solution[END_REF] and Liu, Zhang [START_REF] Liu | Effect of temperature on the size of biosynthesized silver nanoparticle: Deep insight into microscopic kinetics analysis[END_REF] showed that size of particles is a result of a significant effect of 𝑘 1 instead of 𝑘 2 .

Figure 5: Pearson correlation between input variables and the size of AgNPs

The samples of AgNPs were kept at 20 ℃ and 0 ℃ and Fig. 5 shows that the temperature has a positive effect on the particle size, i.e., that smaller particles were observed at low storage temperatures. The full set of data on the particles size depending on process parameters, including storage temperature is provided in SI. Measurements of the silver nitrate concentration for different flow rates of reagents show that the reduction reaction is completed before collection. For example, at 𝑅𝑒 = 5.66 with a 2 m outlet tube, the residual concentration of silver nitrate was around 0.017 ± 0.00053 mM (less than 2 % of initial concentration) whilst at highest flow rate used, 𝑅𝑒 = 16.96 and 2 m length, the residual concentration of silver nitrate was 0mM (~ 3% of initial concentration). The small residual concentrations of silver nitrate showed that the reaction was practically complete for all the investigated values of Re including highest one where residence time was smaller. Izak-Nau, Huk [START_REF] Izak-Nau | Impact of storage conditions and storage time on silver nanoparticles' physicochemical properties and implications for their biological effects[END_REF] and Peng, Krauss [START_REF] Peng | Aging induced Ag nanoparticle rearrangement under ambient atmosphere and consequences for nanoparticle-enhanced DNA biosensing[END_REF] have described in detail the ageing of AgNPs samples stored at room temperature, proposing that growth of the particle size is due to agglomeration or oxidation of silver.

When helical coils are employed, improved mixing due to the formation of Dean vortices leads to smaller particles, and this effect is one of the largest (Fig. 5). It can be assumed that a decrease in particle size due to faster mixing is the result of more intensive nuclei formation at the initial stage of reaction. Of course, faster mixing results also in faster particle growth, but comparing the relative effect of nucleation and growth constant in Fig. 5, it can be suggested that effect of faster nucleation is larger than the effect of faster growth. The effect of the Reynolds number is positive, yet the least significant.

Prediction results with DT-guided experiments

The performance of the developed models (DT, RF, XGBoost) was described through MSE, MAE and 𝑅 2 . In the following Eqs( 10)-( 12), 𝑛 is the number of data within training and testing data sets, 𝑋 𝑖 is the predicted 𝑖 𝑡ℎ value, 𝑌 𝑖 is the actual 𝑖 𝑡ℎ value and 𝑌 ̂ is the mean of true values [START_REF] Chicco | The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation[END_REF].

𝑀𝑆𝐸 = 1 𝑛 ∑ (𝑋 𝑖 -𝑌 𝑖 ) 2 𝑛 𝑖=1 (10)
𝑀𝐴𝐸 = 1 𝑛 ∑ |𝑋 𝑖 -𝑌 𝑖 | 𝑛 𝑖=1 (11)
𝑅 2 = 1 - ∑ (𝑋 𝑖 -𝑌 𝑖 ) 2 𝑛 𝑖=1 ∑ (𝑌 ̂-𝑌 𝑖 ) 2 𝑛 𝑖=1 (12)
The decision tree shown in Fig. 6 examined the synthesis of AgNPs based on the combination of five conditions. The decision tree shows different routes from up to down for preparing AgNPs. Each box presents the value of the parameter, the estimated size, the MSE and the sample sizes. The tree showed that there is a high MSE (Eq. ( 10)) in some synthesis routes e.g., when the growth constant, 𝑘 2 was smaller than 32, the storage temperature was higher than 10 ℃, and the values of 𝐷𝑒/𝑅𝑒 and 𝑅𝑒 were smaller and higher than 0.2 and 0.47 respectively. To improve the performance of the DT developed, synthesis routes were designed under the conditions with high MSE (orange routes on DT). These designed experiments were used not only to improve the performance of the developed DT model but also of other models such us RF and XGBoost. Additionally, the designed synthesis routes (set of parameters identified in the DT with high MSE) were applied to confirm the validity of this approach.

The closer the value of 𝑅 2 is to 1 and the values of MAE and MSE are to zero, the better is the fit of the developed models. The statistical characteristics of the three appropriate models for the original and designed extra data sets (10 extra designed experiment conditions where high uncertainty was found) are given in detail in Table 3.

These results are obtained on the unseen test data. In particular, the designed data sets are obtained by uncertainty analysis (on the training data) of preliminary DT predictions as illustrated in Fig. 6. As

shown, the values of 𝑅 2 are always increasing with the designed experiments for all the models while MAE and MSE are decreasing. These results confirmed that the designed experiments based on DT approach can improve the final predictive models with less cost and time.

As shown in Fig. 7 and Table 3, the computed predictive models, including DT, RF, and XGBoost, which were assisted by DT-based designed of experiments, outperform considerably the ones trained on random extra experiments. These results also confirmed that the proposed DT-based on designed experiments can significantly enhance the performance of other machine learning algorithms with a

Figure 1 :

 1 Figure 1: A) Experimental set up B) Schematic diagram of T-junction microfluidic device.

Figure 2 :

 2 Figure 2: Formation of AgNPs through three different steps, i) the reduction of silver ions to silver

Fig. 4B

  Fig. 4B in Section 3.1. The values of rate constants were found by varying pH of reducing solution at

The

  DT algorithm used to guide the experiments was built based on 20 samples, each being repeated three times to obtain the averaged '(output quantity)' as model output. For fair comparison, 10 extra experiments (30 experiments in total, including three replications) were carried out either randomly or following the DT-guided design of experiments to enhance the performance of the predictive models.

Figure 3 :

 3 Figure 3: Workflow chart and description of methods used in the building of machine learning-guided

Figure 4 :

 4 Figure 4: Application of the F-W mechanism for the analysis of AgNPs formation kinetics. Time evolution of [𝐵] at 400 nm for A) pH 7 and B) pH 12 and different concentrations of the stabilizing

(

  20 initial conditions used to build the DT model), random (10 random extra experiment conditions)

Figure 7 :

 7 Figure 7: Predicted versus real measurements for the size of AgNPs using DT. A) Original training data (20 initial conditions, 60 training samples in total including replications), B) Random extra experiments (10 experiments, 30 validation samples in total including replications) and C) Designed extra experiments (10 experiments, 30 validation samples in total including replications). The conditions that have been used for training, random extra experiments and designed extra experiments are described in detail in SI. The average size of AgNPs on X, Y axes is given in nm.

Table 2 :

 2 Nucleation and growth constants for varied pH and concentrations of trisodium citrate.

	A		C
	pH 7
	B		D
	pH 12
	pH	[TC]	Nucleation constant,𝒌 𝟏	Growth constant, 𝒌 𝟐
		(mM)	(𝒎𝒊𝒏 -𝟏 )	(𝑴 -𝟏 𝒎𝒊𝒏 -𝟏 )
			(average)	(average)
	7	1.91	0.0011	55.51
	7	2.87	0.0013	37.86
	7	3.82	0.0015	44.81

Table 3 :

 3 Decision tree (DT), Extreme Gradient boosting (XGBoost) and Random Forest (RF) resultsThe decision tree based on designed experiments is an effective and low-cost strategy to improve a system in a much more directed manner. The obtained model showed the areas of interest for further experiments (regions with high MSE) and numerical results showed that well designed synthesis routes can improve significantly the performance not only of DT itself but also of other regression methods such as RF and XGBoost compared to randomly extra synthesis routes. In particular, considerable reduction of MSE compared to the original model, 26% for XGBoost and 38% for RF, was found.

		for original, random and designed experiments.
	A. Decision Tree
		MSE	MAE	𝑹 𝟐
	Original	9.43	2.46	0.45
	Random	8.65	2.24	0.49
	Designed	6.41	2.01	0.62
	B. Gradient Boosting
		MSE	MAE	𝑹 𝟐
	Original	8.88	2.27	0.47
	Random	7.26	2.01	0.57
	Designed	6.56	2.16	0.61
	C. Random Forest

Table S1 :

 S1 Particles sizes of AgNPs at different conditions. The sizes were measured as a function of nucleation and growth constants, storage temperature, Reynolds number and Dean number/Reynolds number.

				Designed experiments based on decision tree
		0.0011	55.51	20	0.41	0.0849	8.932	10.141	11.35
		0.0011	55.51	0	0	0.849	11.63	12.9	15.41
		0.0128	26.01	0	0.41	0.849	15.7	14.28	17.1
		0.0013	37.86	0	0	0.0849	12.19	11.34	10.91
		0.0015	44.81	0	0	0.849	11.42	11.92	10.58
		0.0128	26.01	20	0	0.0849	20.78	24.4	23.01
		0.0013	37.86	20	0.41	0.0849	7.161	6.45	8.69
		0.0011	55.51	20	0.41	0.849	7.87	7.074	6.27
		0.0013	37.86	0	0	0.849	10.1	8.78	9.38
		0.0011	55.51	20	0	3.395	18.52	16.47	18.2
				Random experiments
		0.1	77.97	0	0.41	0.849	4.55	3.54	4.547
		0.0013	37.86	0	0	3.395	6.414	6	6.4
		0.0015	44.81	20	0	0.849	12.95	13.65	11.67
		0.0013	37.86	20	0.41	0.849	9.034	8.58	9.88
		0.1	77.97	0	0.41	1.698	3.98	3.68	4.38
		0.012	13.66	0	0.41	1.698	6.996	7.38	6.996
		0.1	77.97	20	0.41	1.698	5.22	5.27	4.94
		0.0015	44.81	20	0	1.698	13.98	11.7	9.878
		0.0011	55.51	0	0	1.698	16.2	15.65	16.11
		0.0011	55.51	20	0	16.96	17.48	19.86	19.01
	Nucleation constant (min -1 ) 0.0011 0.0011 0.0011 0.012 0.012	Growth constant (M -1 min -1 ) 55.51 55.51 55.51 13.66 13.66	Storage temperature (℃) 20 20 20 20 Original experiments Dean number / Reynolds Reynolds number 0 5.662 0.32 0.0849 0.32 3.395 number 0.41 0.0849 0 0 0.849	AgNPs size (1) 15.45 14.05 12.14 9.91 7.45	AgNPs size (2) 20.416 13.87 13.4 10.29 8.7	18.1 AgNPs 13.25 size (3) 12.614 11.31 6.1
		0.0128 0.012	26.01 13.66	20 20	0.41 0	3.395 3.395	19.05 9.24	22.08 11.6	23.65 9.7
		0.0015 0.012	44.81 13.66	0 0	0 0	3.395 1.698	11.1 7.73	13.03 6.45	12 6.28
		0.0015 0.0015	44.81 44.81	20 0	0.41 0	3.395 0.424	8.53 10.96	12.2 10.414	11.16 11.39
		0.0128 0.0015	26.01 44.81	0 0	0 0	3.395 1.698	17.62 12.48	16.71 10.78	16.99 10.95
		0.0015 0.0015	44.81 44.81	0 0	0.41 0	3.395 2.546	5.68 9.77	7.35 9.13	6.62 9.19
		0.0128 0.0013	26.01 37.86	20 0	0 0	3.395 0.424	23.38 16.84	25.52 14.7	24.8 15.67
		0.0128 0.0013	26.01 37.86	0 0	0.41 0	0.0849 2.546	19.83 13.73	19.81 13.57	23.49 13.7
		0.0015 0.0013	44.81 37.86	20 0	0 0	0.0849 1.698	12.64 15.6	13.54 13.14	13.54 13.36
		0.1	77.97	0	0.41	3.395	4.39	3.58	3.79
		0.0013	37.86	20	0	3.395	9.68	8.77	9.42
	0.0128 0.0011 5.References	26.01 55.51	0 0	0 0	0.849 3.395	19.77 15.52	17.39 15.24	18.24 11.34
		0.1	77.97	20	0	3.395	5.08	4.26	4.75
		0.0015	44.81	20	0.41	0.849	6.73	7.54	7.95
	1.	0.012 Nathanael, K., et al., Computational modelling and microfluidics as emerging approaches to 13.66 20 0.41 3.395 9.53 7.28 8.68
	2.	0.012 synthesis of silver nanoparticles -A review. Chemical Engineering Journal, 2022. 436: p. 13.66 0 0 3.395 5.3 6.097 0.0011 55.51 20 0.41 3.395 8.1 10.2 135178. 0.0015 44.81 0 0 0.0849 11.41 10.37	5.29 8.83 10.07
		0.0128	26.01	20	0.41	0.0849	23.93	23.93	23.2
		0.0128	26.01	20	0	0.849	21.52	22.95	25.47

Amir, D., et al., Effect of Stabilizers in the Synthesis of Silver Nanoparticles and Methylene Blue Oxidation. IOP Conference Series: Materials Science and Engineering, 2021. 1192(1): p. 012031.

Acknowledgements

The authors gratefully acknowledge support from the Engineering & Physical Sciences Research Council, UK, through the PREMIERE Programme Grant (EP/T000414/1). KN is funded by a PhD scholarship from the School of Chemical Engineering, University of Birmingham.

considerable reduction (-26% for XGBoost and -38% for RF) of MSE compared to the original model. Note, this paper consists of a 'proof of concept' of the proposed DT-guided design of experiments and more experimental data should be collected to obtain more reliable predictions.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supporting information