Rebecca Salles
email: rebeccapsalles@acm.org

Esther Pacitti
email: esther.pacitti@lirmm.fr

Eduardo Bezerra
email: ebezerra@cefet-rj.br

Celso Marques

Carla Pacheco
email: carla.pacheco@eic.cefet-rj.br

Carla Oliveira
email: carla.oliveira@ibge.gov.br

Fábio André

Machado Porto
email: fporto@lncc.br

Eduardo Ogasawara
email: eogasawara@ieee.org

Cefet / Rj

Carla Olveira

Fabio Porto

TSPredIT: Integrated Tuning of Data Preprocessing and

Keywords: Time Series Prediction Models time series, prediction, data preprocessing, machine learning, hyperparameter optimization

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 prediction of time series has gained more attention in the last decades. Many time series prediction methods have been developed and can be found in the literature [START_REF] Cheng | Time series forecasting for nonlinear and non-stationary processes: A review and comparative study[END_REF]. An adequate prediction method is mandatory for building the right model [START_REF] Salles | TSPred: A framework for nonstationary time series prediction[END_REF], especially for data-driven models. They are generally organized between statistical and machine learning [START_REF] Lindemann | A survey on long short-term memory networks for time series prediction[END_REF]. These types of methods usually have to set hyperparameters. In this sense, hyperparameter optimization is a fundamental step since it can influence the predictive performance of the resulting models [START_REF] Lim | Time-series forecasting with deep learning: A survey[END_REF][START_REF] Khalid | A survey on hyperparameters optimization algorithms of forecasting models in smart grid[END_REF].

Additionally, these models might be improved by adequate data preprocessing activities. Most of these methods tend to be optimistic regarding their assumptions over the time series and are not ready to handle nonstationarity [START_REF] Cheng | Time series forecasting for nonlinear and non-stationary processes: A review and comparative study[END_REF][START_REF] Salles | Nonstationary time series transformation methods: An experimental review[END_REF]. They also suffer from the presence of concept drift [START_REF] Lu | Learning under Concept Drift: A Review[END_REF] or lack of data [START_REF] Talavera | Data Augmentation techniques in time series domain: A survey and taxonomy[END_REF]. These two cases are related. When concept drift occurs, generally, there are few samples to support model building. Usually, nonstationarity demands transformation methods to address this issue [START_REF] Salles | Nonstationary time series transformation methods: An experimental review[END_REF]. Besides, while working with small samples, data augmentation techniques are also needed [START_REF] Mumuni | Data augmentation: A comprehensive survey of modern approaches[END_REF][START_REF] Wen | Time Series Data Augmentation for Deep Learning: A Survey[END_REF][START_REF] Talavera | Data Augmentation techniques in time series domain: A survey and taxonomy[END_REF].

preprint For the algorithm to make predictions with greater accuracy, optimizing the hyperparameters is necessary [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. Hyperparameters are values that make up the initial configuration of the learning algorithm [START_REF] Garcia | Data Preprocessing in Data Mining[END_REF]. Hyperparameters are also present in data preprocessing methods. Several factors influence the predictive performance of time series models, mainly choosing and tuning the right methods for data preprocessing and hyperparameters.

Regarding statistical learning, some methods seamlessly integrate hyperparameters optimization of data processing techniques. It includes the autoregressive integrated moving average (ARIMA) algorithm that optimize parameters (p, d, q) for ARIMA [START_REF] Hyndman | Automatic time series forecasting: The forecast package for R[END_REF]. The d parameter represents the Integrated part of ARIMA and performs the differentiation of observations internally as a preprocessing step for the series to be stationary. The p parameter corresponds to the AR part of ARIMA, the number of autoregressive terms. The MA model works with the size of the moving average window and is represented by the q parameter. This method tunes altogether differentiation, autoregressive, and moving average models [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF]. Conversely, there are many frameworks for machine learning, such as Scikit learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF][START_REF] Hao | Machine Learning Made Easy: A Review of Scikit-learn Package in Python Programming Language[END_REF], which provides (i) a broad range of prediction methods, (ii) an extensive set of preprocessing methods, (iii) hyperparameter optimization features for machine learning. However, directly optimizing data processing and machine learning is left for users to program according to their needs.

In this context, this paper presents TSPredIT, an evolved version of TSPred [START_REF] Salles | TSPred: A framework for nonstationary time series prediction[END_REF] that seamlessly integrates the tuning of data preprocessing and time series prediction models for univariate time series. It only concerns regression models and is specialized integrating data transformation methods and data augmentation to aid in building machine learning methods (MLM) prediction models. TSPredIT is made available as an R-package. It is the first tool to seamlessly integrate a broad range of data transformation and preprocessing methods and state-of-the-art statistical and machine learning prediction methods for addressing nonstationary time series. The package automates the time series prediction process and parameterization while enabling user-defined prediction methods and data transformations, including code built in other languages like Python. Due to that, the features provided by TSPredIT are shown to be competitive regarding time series prediction accuracy.

Besides this introduction, the paper is organized into five more sections. Section 2 presents the background, while Section 3 presents the related work. Section 4 presents the TSPredIT, which evolves from the previous version of TSPred [START_REF] Salles | TSPred: A framework for nonstationary time series prediction[END_REF]. Section 5 provides a clear example of using TSPredIT's main features. Conversely, Section 6 characterizes the effect of choosing data preprocessing and different MLM during prediction. Finally, Section 7 concludes the work.

Background

Time series prediction is commonly associated with the scenario of regression. For simplicity, the paper may refer to prediction and regression interchangeably. Relevant models adopted for time series prediction generally fall into the categories of statistical or machine learning models [START_REF] Salles | A framework for benchmarking machine learning methods using linear models for univariate time series prediction[END_REF]. The accuracy of the predictions depends on the quality of the historical data, the appropriateness of the model, and the assumptions made about the underlying processes driving the time series [START_REF] Hyndman | Forecasting: principles and practice[END_REF][START_REF] Gujarati | Essentials of Econometrics[END_REF]. Figure 1 presents a general time series prediction process. It encompasses five main activities. It provides a general framework for predicting a time series based on a particular setup of preprocessing methods and prediction models. They are briefly described here, and some parts are detailed in the following sections. In such a case, the training slice is again split into a novel training and validation set. Alternative models exploring hyperparameters values are built using the novel training and evaluated using the validation set. Once hyperparameters are fixed, a single model is built using the entire training dataset. From this moment, the model is available for use. Activity 3, also in blue, refers to the model prediction. It is worth mentioning that the predicted values are not in the time series domain. In this sense, they can not be directly evaluated. Data postprocessing is needed in this case. Activity 4, also in purple, corresponds to the postprocessing of predictions, reversing transformations applied to the time series data in Activity 1. In a macro view, the data is normalized (scaled) and given as input to an algorithm. After the forecast, a denormalization process maps back the predicted values into the original scale of the time series. An example of postprocessing can be seen in previous work [START_REF] Ogasawara | Adaptive Normalization: A novel data normalization approach for non-stationary time series[END_REF], where the data is denormalized to measure the error in the same scale for comparison purposes.

Finally, Activity 5, in pink, is the evaluation of prediction errors yielded by the model, as well as model fitness metrics. If the results are inadequate, this process can be revised and repeated to refine models. This process iteratively improves the quality of predictions (for time series prediction) or model adjustment (for time series modeling). The prediction can be evaluated in several ways, mostly measuring the errors between prediction and actual observation, such as Mean Square Error (MSE) and symmetric MAPE (sMAPE) in a test set. Alternatively, they can be measured by the level of model adjustment, such as Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) [START_REF] Salles | TSPred: A framework for nonstationary time series prediction[END_REF].

This process provides a systematic way of predicting a time series based on particular preprocessing and prediction methods. It also focuses on prediction and model evaluation, that is, evaluating the accuracy of prediction and the fitness of a model. Such evaluation may indicate a demand for refining and perfecting the preprocessing-modeling setup and its parameters to obtain a more accurate model. This process may be repeated if the evaluated time series prediction model does not reach the desired accuracy. This process enables benchmarking different preprocessingmodeling setups.

Related Work

Several authors focused on the task of exploring different models. Ramey [START_REF] Ramey | sortinghat: sortinghat[END_REF] and Lessmann et al. [START_REF] Lessmann | Benchmarking classification models for software defect prediction: A proposed framework and novel findings[END_REF] developed frameworks for evaluating classification models and algorithms. Moreover, Bischl et al. [START_REF] Bischl | Mlr: Machine learning in R[END_REF] and Eugster and Leisch [START_REF] Eugster | Bench Plot and Mixed Effects Models: First Steps toward a Comprehensive Benchmark Analysis Toolbox[END_REF] developed the R-packages mlr and benchmark, respectively, which provide tools for executing automated experiments when benchmarking a set of models for data mining tasks such as classification and regression. These packages are designed to support tabular data and focus on benchmarking based on plot visualization.

Hyndman and Khandakar [START_REF] Hyndman | Automatic time series forecasting: The forecast package for R[END_REF] and Hyndman et al. [START_REF] Hyndman | A state space framework for automatic forecasting using exponential smoothing methods[END_REF] present frameworks for automatic forecasting using mainly statistical models such as ARIMA and exponential smoothing state space model (ETS). Hyndman and Khandakar [START_REF] Hyndman | Automatic time series forecasting: The forecast package for R[END_REF] produced the well-known R-package named forecast, which can be used for automatic time series prediction. The R-package of Moreno, Rivas, and Godoy [START_REF] Moreno | predtoolsTS: Time Series Prediction Tools[END_REF] also facilitates time series prediction using simple differencing (diff) and Box-Cox transform (BCT). Furthermore, we observed three works worth mentioning. Diebold and Mariano [START_REF] Diebold | Comparing predictive accuracy[END_REF] propose various tests to compare the predictive accuracy of two different prediction models. Diebold and Lopez [START_REF] Diebold | 8 Forecast evaluation and combination[END_REF] propose an ensemble approach using different prediction models. Kumar et al. [START_REF] Kumar | Model Selection Management Systems: The Next Frontier of Advanced Analytics[END_REF] propose a class of analytics systems to manage model selection using key ideas from data management research.

Besides, hyperparameters optimization is also a deeply studied subject [START_REF] Izaú | Towards Robust Cluster-Based Hyperparameter Optimization[END_REF][START_REF] Khalid | A survey on hyperparameters optimization algorithms of forecasting models in smart grid[END_REF][START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF]. The studies may focus on exploring the hyperparameter search space using a certain heuristic. Conversely, some approaches target the right establishing of the hyperparameters to explore using either grid search or a more advanced search strategy. The Grid Search approach is commonly adopted to explore a broad range of hyperparameter settings. It consists of repeatedly training the learning algorithm with different possible hyperparameter settings combinations. At the end of the process, the hyperparameter setting that resulted in the lowest prediction errors (measured in a separate validation set) is chosen [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. Such optimized hyperparameter settings can then be used to fit the learning model [START_REF] Ran | Parameter identifiability in statistical machine learning: A review[END_REF]. All these approaches try to lower global prediction error in machine learning but are resilient to the problem of data overfitting. Such an issue occurs when the fitted model is too dependent on the training dataset. One consequence is the fitted model's inability to generalize to unseen data observations [START_REF] Sarwar Murshed | Machine Learning at the Network Edge: A Survey[END_REF].

All in all, several works present frameworks and tools for MLM performance assessment. Nonetheless, to our knowledge, no work proposes and implements a framework for the seamless integration of hyperparameter optimization of data preprocessing and time series prediction methods.

The main modules of the TSPredIT framework are depicted in Figure 2 as a UML class diagram. TSPredIT has five main functionality modules: Preprocessing (in purple), Modeling (in blue), Sampling (in yellow), Evaluating (in pink), and Tuning (in green). Together, they are used to support the time series prediction process. The colors of the classes are associated with their participation in the time series prediction process as depicted in Figure 1.

All classes are inherited from TSBase. It provides a basic fit method and some attributes for introspection (to support provenance). It also includes a fit analysis of the data to adjust basic parameter values. A TSData class also provides a uniform perspective for time series data and its transformation to sliding windows. These two types are a specialization of TSData.

The first module is responsible for preprocessing (transform) and postprocessing (inverse transform) a time series. The model groups two main features. The first is related to data transformation. Especially it includes the implementation of the main nonstationary time series transformation methods [START_REF] Salles | Nonstationary time series transformation methods: An experimental review[END_REF], being either mapping-based, namely the logarithmic transform (LT), BCT, percentage change transform (PCT), moving average smoother (MAS), and diff, or splitting-based, such as empirical mode decomposition (EMD) and wavelet transform (WT). All these methods are implemented as a specialization of the Transformation class. The basic implements fit, transform, and inverse transform to provide the desired behavior.

Furthermore, it also groups a set of methods related to data augmentation. All methods include warping, flipping, and jittering [START_REF] Mumuni | Data augmentation: A comprehensive survey of modern approaches[END_REF][START_REF] Talavera | Data Augmentation techniques in time series domain: A survey and taxonomy[END_REF]. These methods are used during tuning, which is explained later. The Modeling module, in blue, is responsible for modeling (fit) and predicting (predict) a time series based on a particular time series prediction method. These tasks are specialized for either statistical or machine learning models. For the latter, the framework is prepared to perform any necessary machine learning life-cycle tasks during the training and prediction steps, including coercing data into sliding windows, normalizing and transforming the input data. This module includes the implementation of the statistical models: ARIMA, Holt-Winter's exponential smoothing (HW), theta forecasting (TF), and ETS. MLM models include multilayer perceptron network (MLP), random forest regression (RFR), preprint support vector machine (SVM), MLP, and extreme learning machines network (ELM). Furthermore, the module provides deep learning models available in the PyTorch library, namely convolutional neural network (Conv1D) and long shortterm memory neural network (LSTM). Models are associated with zero or one data preprocessing transformation and data augmentation technique. However, it might be set up with multiple candidate options chosen during tuning.

The Tuning module is designed to provide hyperparameter optimization. It is invoked during the fitting of a model. Hyperparameter optimization occurs whenever the modeling of time series or preprocessing transformations has a degree of freedom to adjust. The default Tuning applies time series cross-validation using the training set [START_REF] Hyndman | Forecasting: principles and practice[END_REF]. Data augmentation might be applied in each partition, and the model is trained after applying the data transformation. The fittest model and data preprocessing method were discovered using a grid-search, i.e., the one that leads to better prediction during cross-validation, is chosen for training using the entire training set. Since this is not the only way of conducting Hyperparameter optimization, the default Tuning class can be specialized to provide other ways of enhancing this feature.

Finally, the Evaluating module, in pink, is responsible for assessing the model fitness and quality of predictions. These tasks are specialized for computing either prediction accuracy (error) measures or model fitting criteria. The available prediction accuracy measures include MSE, sMAPE, and maximal error. It also includes model fitness criteria such as AIC, BIC, and log-likelihood [START_REF] Davydenko | Measuring Forecasting Accuracy: The Case Of Judgmental Adjustments To Sku-Level Demand Forecasts[END_REF].

TSPredIT can integrate the described modules in a workflow, connecting the five modules described. The package provides the means to perform the benchmarking of several prediction models. It is important to remark that although providing several pre-implemented options, TSPredIT design enables the user to define and apply customized time series prediction methods.

Moreover, the package provides several automatized features for any time series prediction application. Among them, some of the main features are (i) seamless recursive combination of two or more transformation methods; (ii) seamless integration of transformation methods to the prediction process [START_REF] Salles | Nonstationary time series transformation methods: An experimental review[END_REF], which demands the combination of predictions for each component resulting from data decomposition (first package to include this approach); (iii) transformation and model parameter selection; (iv) multistep-ahead or one-step-ahead predictions; (v) rolling origin evaluation [START_REF] Hyndman | Forecasting: principles and practice[END_REF] for both statistical and machine learning models, and (vi) machine-learning life-cycle tasks performed during training and prediction steps. Data normalization and sliding window transformation are seamlessly conducted during machine learning model training.

The framework is implemented in R using the S3 class system [START_REF] Wickham | Advanced R[END_REF]. TSPredIT is currently available on GitHub1 . It is an ongoing evolution of TSPred [START_REF] Salles | TSPred: A framework for nonstationary time series prediction[END_REF], built on top of the DAL Toolbox2 .

Usage of TSPredIT

This section gives examples of TSPredIT usage. The first example corresponds to a time series prediction using a wrapper for the MLP model using sliding windows min-max normalization. The hyperparameter tuning applies a grid search using time series cross-validation.

The Listing 1 the TSPredIT R-package processes a time series. The lines (1-3) of code target and load the installation of TSPredIT. The components for the time series process can be defined separately to enable reuse. Besides, the dataset used is made available in the R-package. It is loaded using the data function (line 6).

The time series is converted into sliding windows (line 8). All sliding windows are shifted with overlap with step 1 by default. In the example, the size of the sliding windows is 8. Besides, the last 4 windows are reserved for testing (line 10), and the complement is used for testing in order to control more precisely which observations is being considered, as a fine tuning. Finally, the training data is separated into input and output (line 12).

The hyperparameter setup is established in lines 15-17. It indicates the data preprocessing option of min-max sliding windows. The input size for model building varies between 3 and 7. The base model is related to MLP. No data augmentation method is used in this example: ts augment(). Also, some specific parameters for MLP are indicated in lines 18-19. It provides ranges for the number of neurons in the hidden layer, the rate of decay during training, and the maximum number of iterations (fixed). In lines 20-21, the actual tuning is executed using the training set. Internally, it splits the data using time series cross-validation. The build model hyperparameters that work better during crossvalidation are used to build the final model using the entire training set. In this example, 500 configurations were explored, each one ten times due to the default ten-fold cross-validation. # Loading f e r t i l i z e r s d a t a s e t

6

> data ( f e r t i l i z e r s ) 7 # C o n v e r t i n g t o s l i d i n g windows

8

> t s <t s data ( f e r t i l i z e r s $ b r a z i l n , sw = 8 ) 9 # P a r t i t i o n i n g i n t o t r a i n i n g and t e s t i n g 10 > samp <t s sample ( ts , t e s t s i z e = 4 ) 11

# S e p a r e t i n g i n p u t and o u t p u t f o r t r a i n i n g

12 > i o t

5

> ranges <l i s t ( n h i d = 1 : 2 0 , 6 a c t f u n =c ( ' s i g ' , ' radbas ' , ' t r i b a s ' , ' r e l u ' , ' p u r e l i n ' ) )

7 > model <-f i t ( tune , x = i o t r a i n $ i n p u t , 8 y = i o t r a i n $ o u t p u t , ranges )
Lines 23-26 present the level of adjustment for the time series in the training set. This aspect is important since the error level in training is commonly higher during testing. It provides an expected entry error. Lines 28-31 present the prediction for testing. It applies a rolling origin with one step-ahead prediction, leading to four predictions from previously known observations. The sMAPE is presented in line 33.

To clarify how extensible is TSPredIT in providing alternative MLM, Listing 2 changes four lines of code to switch the MLP to ELM, with different ranges of hyperparameters to explore. This feature is possible due to the wrapper classes provided by TSPredIT that integrate state-of-the-art methods. Additionally, novel methods can be wrapped. Writing the fit and predict methods is needed to incorporate a novel method at TSPredIT. TSPredIT was experimentally evaluated to expose the main features of the framework. For that, it was derived a time series dataset from public data available at the International Fertilizer Association (IFA)3 as a proof of concept of the present framework. This dataset of fertilizers was explored in deep in previous work [START_REF] Pacheco | Exploring Data Preprocessing and Machine Learning Methods for Forecasting Worldwide Fertilizers Consumption[END_REF]. It contains data on the annual consumption of three fertilizers (K 2 O, N, P 2 O 5 ) among the top ten main consumer countries. Each time series contains 60 observations from 1961 to 2020. Observations from 1961-2016 are used for training, and observations from 2017-2020 are used for testing. For the evaluation, we selected Brazil, the third major fertilizer consumer. All coding for the experimental evaluation is available 4 .

The goal of this paper is to explore multiple facets of TSPredIT. The first experiment evaluated the effect of data transformation (swminmax, diff, an, gminmax) during prediction using MLP as MLM. A second evaluation explored the adoption of data augmentation techniques while fixing both MLP as MLM and swminmax as a data preprocessing technique. Figure 4 compares the performance of not applying data augmentation (none), using jittering (jitter) and warping stretching (stretch). As it can be observed, none was better both in K 2 O and P 2 O 5 . However, for N, both jitter and stretch were worth value. The prediction performance increased by more than 1%. The data augmentation technique was seamlessly applied during time series cross-validation for training during hyperparameter optimization.

Finally, the third evaluation studied the adoption of different MLM for predicting P 2 O 5 using an. Figure 5 presents the evaluation of using Conv1D, ELM, MLP, RFR, SVM, LSTM for this scenario. Table 1 presents the hyperparameters explored.

All these methods, except for Conv1D and LSTM, explored a similar amount of hyperparameter combinations (about 250 options each). The methods mostly could not improve the performance of the MLP. The exception was LSTM, which improved prediction by more than 1%.

As a proof of concept, these results can explore the capability of TSPredIT in combing a broad range of data preprocessing techniques and state-of-the-art MLM. Besides, the framework can provide hyperparameter optimization exploring both features, aiding the selection choice for these methods.

Figure 1 :

 1 Figure 1: Time series prediction process [29] Activity 1, depicted in Figure 1 in purple, refers to acquiring the time series and performing data preprocessing. It is generally associated with data cleaning, normalization, and transformation but might include other techniques, such as data augmentation. The transformations commonly change the time series domain values, and their parameters must be stored to support later detransformation to the original domain. For time series prediction, splitting the time series into a training and test set is also important during data preprocessing. All data preprocessing parameters should be computed during training and reapplied from the tune-values of training during the test. The model is built using the training slice and evaluated using the unseen test set, always ahead. However, when the goal is to adjust a model for the time series, the model does not need to be partitioned into a training and test set.

  The second module is related to Sampling. It is responsible for converting a time series to sliding windows. It is also responsible for separating the data (TSData) into training and testing. The test size uses only recent observations to avoid introducing new data during training. Like augmentation, sampling is used during tuning. Other relevant preprocessing methods for time series prediction are specialized from Preprocessing class, in purple. It includes support for handling missing values and data normalization such as Sliding Windows min-max normalization (swminmax), Min-max normalization (gminmax), diff, and Adaptive normalization (an) methods [24].

Figure 2 :

 2 Figure 2: TSPred-IT functionality modules and pre-implemented algorithms

preprint Listing 1 : 1 > l i b r a r y ( d a l t o o l b o x ) 2 #

 112 Example of time series prediction process in TSPredIT h t t p s : / / c e f e t -r j -d a l . g i t h u b . i o / t s p r e d i t /3> l i b r a r y ( t s p r e d i t ) 4 5

Figure 3

 3 compares these methods during testing for the three fertilizers K 2 O, N, and P 2 O 5 in Brazil. It presents the sMAPE error during testing. The swminmax outperformed other methods for K 2 O and N. However, in P 2 O 5 , an was better followed close by swminmax.

Figure 3 :

 3 Figure 3: Comparison of data transformations applied (swminmax, diff, an, gminmax)

  preprint

Figure 4 :

 4 Figure 4: Comparison of data augmentation applied (none, jitter, stretch)

Figure 5 :

 5 Figure 5: Comparison of evaluation of using different MLM for predicting P 2 O 5 using an

https://github.com/cefet-rj-dal/tspredit

https://cran.r-project.org/web/packages/daltoolbox/index.html

http://www.fertilizer.org

https://eic.cefet-rj.br/ ˜dal/tspredit/

Acknowledgements

The authors thank CNPq, CAPES (finance code 001), and FAPERJ for partially sponsoring this research.

The combination of time series transformation methods in prediction with decomposed time series, transformation and model parameter selection, multi-step or one-step-ahead prediction, rolling origin evaluation, and the management of sliding windows is a key differentiation for TSPredIT. Several benchmark datasets from time series prediction competitions come bundled with TSPredIT. This new version enables users to practice data transformation and prediction methods, gaining confidence in the developed prediction models. Besides, the framework was designed to enable users to implement their customized methods. For example, both LSTM and Conv1D were added to TSPredIT as customized methods implemented in Python. Future updates will expand the range of implemented preprocessing methods, MLM, and evaluation metrics, especially empowering the hyperparameter selection targeting choosing combinations of data preprocessing and MLM that led to more stable models. This work leaves room for future implementation including multivariate time series, cost of computation time, and short and long term prediction in classical problems.


Chloé Andrade

Saly Hawila

Ahmad Abdallah

Jean-Luc Rukemampunzi

Adel Mesbah

Nathalie Guillou

Florent Perret

Stefan Wuttke

Thomas Niehaus

Régis Debord

Olivier Boisron

Stéphane Pailhès

Aude Demessence

A p-type semi-conducting copper( i )-1,3-benzenedithiolate 2D coordination polymer with high Seebeck coefficient

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A new p-type semi-conducting and 2D coordination polymer (CP) based on thiolate ligand and Cu(+I) exhibits a record Seebeck coefficient among CPs, equivalent to that obtained for inorganic compounds. Such high value is originating from the 1D {Cu-S} n tubular inorganic network, that is connected through the bridging 1,3-benzenedithiolate ligand.

Considering the actual global warming and our over-reliance on fossil fuels, there is a tremendous need in new energy technologies with low environmental and societal impacts. Thermoelectricity, which is the ability of a material to convert heat into electricity and vice versa, is an ideal technology that would allow transforming heat waste from industries, engines, power plants, household appliances and even human bodies into power; providing sustainable energetic cycles. To reach this goal, a thermoelectric material should have a good performance and its efficiency is evaluated by a dimensionless figure of merit, ZT = S 2 σT / κ, where S, σ, T and κ are the Seebeck coefficient, the electrical conductivity, the absolute temperature and the thermal conductivity, respectively. Bearing in mind the Carnot cycle, the ZT value of a material should be greater than 3 for a device efficiency. 1 From the discovery of thermoelectricity phenomenon, extensive studies have been carried out on inorganic materials, mostly the chalcogenides such as Bi 2 Te 3 , Sb 2 Te 3 , PbTe and SnSe. Those thermoelectric (TE) materials have high Seebeck coefficient, up to around 500 µV.K -1 , but they have high thermal conductivity, which results in a modest ZT, still below 1. 2 The record ZT value has been reported in 2019 for an alloy, Fe 2 V 0.8 W0 .2 Al, with a ZT value of 5 to 6. 3 In addition, these thermoelectric inorganic materials have the drawbacks of poor processability, high cost and toxicity, preventing today further sustainable applications. In the meantime, -delocalized organic polymers, such as poly (3,4-ethylenedioxythiophene) polystyrene sulfonate (PEDOT:PSS) or polyaniline (PANI), have appeared as good alternatives as TE materials with good shaping opportunity, light weight, non-toxicity, low cost and more importantly low thermal conductivity ( 0.2 W.m -1 .K -1 ). 4 Nevertheless, the organic materials have low Seebeck coefficient ( 30 µV.K -1 ) resulting in still modest ZT. Thus, there is huge need to further improve the ZT value of materials to achieve the Carnot efficiency, and significant efforts have to be devoted into the development of TE materials exhibiting high Seebeck coefficient, high electrical conductivity, and low thermal conductivity. An alternative is to use hybrid materials to take advantages of the inorganic and organic networks. Thus, TE composite materials which are the mixture of two compounds, are a good strategy, but the synergy between the two entities is difficult to control and requires multistep syntheses. 5 Metal Organic Frameworks (MOFs) and Coordination Polymers (CPs), constructed from metal ions bonded to organic moieties, offer a diversity of structural topologies ranging from one-dimensional (1D) to three-dimensional (3D) arrangements, with possible porosity coupled with other different functionalities, which allow fine tuning of their electronic structures, resulting in intrinsic physical properties. 6 Those MOFs and CPs have recently appeared as potential TE materials (Table S1). 7 Indeed, the 2D copperbenzenehexathiolate, Cu x BHT, exhibits a record electrical conductivity up to 2000 S.cm -1 , due to the delocalized charge through the π-d conjugated system. 8 The TE properties have been evaluated for different formula of Cu x BHT with x = 3, 4 and 5.5 and the best ZT of 0.08 at RT has been obtained with a Seebeck coefficient of 80 µV.K -1 and a thermal conductivity of 0.59 W.m -1 .K -1 . 9 The lowest thermal conductivity of 0.2 W.m - 1 .K -1 has been measured in the two Ni 3 (2,3,6,7,10,11hexaiminotriphenylene) 2 and [Ni 3 (perthiolated coronene)] n CPs, nevertheless their low Seebeck coefficients of -12 and 47 µV.K -1 imply low ZT of 1.19.10 -3 and 0.003, respectively. 10 Zn-HAB (HAB : hexaaminobenzene), a 3D porous network, shows the best intrinsic Seebeck coefficient of 200 µV.K -1 at RT and under vacuum, but this MOF presents a low conductivity and suffer of instability under O 2 , due to the weak interactions between Zn(II) and the amino groups. 11 Consequently, there is still a tremendous need to explore more MOF/CP systems as new TE materials to understand the impacts of the metal ions, their oxidation state, the coordinating functions, the topology of the ligands and the dimensionality of the network on the electric and thermal properties. Based on the record properties of Cu x BHT and the high stability of copper-thiolate compounds, due to the strong interactions between the thiolate ligands (soft base) and the Cu(+I) (soft acid), copperthiolate CPs are very promising systems as TE materials. 12 Apart from Cu x BHT, few more conducting coppermonothiolate CPs have been reported, such as [Cu(6-mna)] n (5.4 × 10 -5 S.cm -1 at RT), 13 {[Cu 2 (6-Hmna)(6-mna)]•NH 4 } n (11 S.cm -1 at RT on single crystal) 14 (6-Hmna = 6-mercaptonicotinic acid, 6-mna = 6-mercaptonicotinate), [Cu(p-SPhOH)] n (120 S.cm -1 ), 15 and Cu[Cu(pdt) 2 ], (pdt = 2,3-pyrazinedithiolate) (6.10 -4 S.cm -1 at RT), pointing out their ability as semiconducting CPs. 16 Here we report, the copper(I)-1,3benzenedithiolate, [Cu 2 (1,3-BDT)] n , a new 2D CP that is made of aligned 1D Cu-S tubes connected together through the bridging ligands. This CP exhibits intrinsic and bulk electrical conductivity of 1.5 mS.cm -1 under air at RT and a Seebeck coefficient of 420 µV.K -1 , a record value among CPs and MOFs, and equivalent to inorganic materials.

[Cu 2 (1,3-BDT)] n is obtained by mixture of two equivalents of copper(I) chloride with one equivalent of 1,3-benzenedithiol in an acidic aqueous solution, heated in a sealed vial at 120°C. As shown by the scanning electron microscopy image (Fig. S1), thin needles of few micrometers length are formed. Elemental and thermo-gravimetric analyses lead to the chemical formula of [Cu 2 (1,3-BDT)] n , which confirms the S / Cu(I) = 1 ratio and the charge balance (see ESI and Fig. S2). The infrared spectrum also confirms the coordination of the thiolate functions with the absence of thiol function, usually observed at 2500 cm -1 (Fig. S3). The powder X-ray diffraction (PXRD) data, collected at the CRISTAL beamline of the SOLEIL synchrotron facility (Gif-Sur-Yvette, France), show good crystallinity of the compound and allow its ab initio structural determination. The final Rietveld refinement (Fig. S4) leads to satisfactory model indicator and profile factors (Table S2). [Cu 2 (1,3-BDT)] n crystallizes in the monoclinic P2 1 /c space group. Its structure is made of inorganic Cu-S chains, where copper and sulphur atoms are three coordinated, forming distorted Cu 3 S 3 hexagons packed together along b axis (Fig. 1a andb, S5). These Cu-S tubular chains are then connected together through the bridging 1,3-BDT ligands to form a 2D network (Fig. 1c), totally different to that observed for the 3D [Ag 2 (1,3-BDT)] n analogue. 17 Such Cu(I)-S tubular chains made of distorted Cu 3 S 3 hexagons have been already observed in few [Cu(SR)] n CPs with monotopic thiolate ligands, with R = Me, 18 Ph, p-PhMe, p-PhOMe 19 and o-PhCO 2 H. 20 In this 2D CP, the Cu-S distances range from 2.20(2) Å to 2.26(3) Å, the Cu-S-Cu angles from 75.0(7) ° to 125(1) ° and the S-Cu-S angles from 109.1(7) ° to 127.2(8) ° (Table S3). However, [Cu 2 (1,3-BDT)] n is the first 2D copper-thiolate CP bearing with these 1D Cu-S tubular chains, thus making it possible to combine in one compound, a 1D inorganic chain and a 2D metal--conjugated path. The Cu-Cu distances of 2.947(4) Å and 3.797(3) Å point out the absence of cuprophilic interactions, that can sometimes be observed in d 10 coinage metal compounds. 21 The X-ray photoelectron spectroscopy (XPS) of [Cu 2 (1,3-BDT)] n also confirms the oxidation state of Cu(I). Indeed, the 2pcopper core levels exhibit dominant peaks at 932.8 eV and 952.6 eV, corresponding to Cu(I) 2p3/2 and Cu(I) 2p1/2, respectively (Fig. S6). The direct band gap obtained from the Tauc plot of the UVvisible absorption spectroscopy curve is 2.26 eV (Fig. S7), which is of the same order of magnitude as that of the other reported 1D copper(I)-thiolate CPs (2.23-2.68 eV). [19][20] The thermoelectric properties of [Cu 2 (1,3-BDT)] n were evaluated from dense pellets with a diameter of 4 mm and a thickness of 1 mm (Fig S8), prepared at room temperature Please do not adjust margins Please do not adjust margins under pressure, with typically, around 40 mg of crushed sample pressed at 0.6 MPa for 15 minutes. PXRD data of these pellets show that the sample remains crystalline and stable under this mechanical pressure, and that no decomposition is observed when it is stored under air for few months (Fig. S8). The electrical transport properties were investigated first by four-probe method at the surface of the pellet from RT to 390 K. At RT, the current-voltage measurements on the pellets in the four-probe configuration, show a linear behaviour following the Ohm's law (Fig. S9, which reveals an electrical conductivity of 8.3.10 -4 S.cm -1 . This modest value is close to the one of the 1D [Cu(SPh)] n analogue, which has an electrical conductivity of 3.7.10 -4 S.cm -1 , pointing out the implication of the Cu-S chains in the electron transport rather than the aromatic benzene ring. 22 The variable-temperature measurements with a four-probe configuration using metallic tips applied on the surface of a pressed pellet show that the electrical conductivity increases with the rise in temperature, which is consistent with bulk transport dominated by a thermally activated process (Fig. 2a andS10). Fitting the data with the Arrhenius equation yields an activation energy (E a ) of 0.13 eV (Fig. S11). 23 In order to ensure that the electrical conductivity is not a surface effect on the pressed pellet, bulk conductivity also was measured. For that, a corded saw was used to cut the sintered pellet into a rectangular cuboid of 1.4  2  2.5  4 mm (Fig. S12) and two electrodes have been welded on the metallised sides of this cuboid to apply a uniform current. The voltage is then measured between differently spaced electrodes, alloying the measurements of the voltage drop as function of distances, in order to extract the resistance and the contact resistance (see insert of Fig. 2a andS13). A low contact resistance of 0.017 kΩ was determined (two orders of magnitude lower than the one of the sample). The bulk electrical conductivity of the CP measured in this configuration is 15.10 -4 S.cm -1 , close to the value obtained with the previous surface four-probe experiment, proving that the electrons can move through the bulk pressed pellet prepared without any care of orientation. The Seebeck coefficient (S) was measured on bulk samples of [Cu 2 (1,3-BDT)] n , inserted in a primary vacuum chamber, at different temperatures ranging from RT to 400 K, below the decomposition temperature of 520 K (Fig. S2). The S value at RT is positive and high: +420 ± 10 µV.K -1 (Fig. S14). The positive sign indicates that the dominant charge carriers are holes with a Fermi level close to the valence band. This S is the highest reported among conducting MOFs, CPs, composite materials and organic polymers (Table S1). 7 Interestingly, the S increases upon heating up to a certain temperature, around 350 K, from which it decreases reversibly up to 380 K (see Fig. 2b). This unusual existence of temperature maximum will be correlated below to the particularities of the electronic band structure. 24 The power factor (PF), expressed as PF = , is 0.026 µW.m - 𝑆 2 𝜎 1 .K -2 at RT and reaches 0.034 µW.m -1 .K -2 at 349 K. Due to the modest electrical conductivity of [Cu 2 (1,3-BDT)] n , this PF is below the best values reported for MOFs and CPs with a record of 88.2 µW.m -1 .K -2 for the nanosheets of Cu x BHT (Table S1). 9b Nevertheless, as mention earlier, the electrical conductivity of [Cu 2 (1,3-BDT)] n is an average over all the directions of the needle-like crystallites constituting the powder, that has been simply pressed to form the pellets. Further studies on aligned crystallites or on a single crystal would allow extracting the intrinsic conductivity of the 1D-CuS chains that might be much higher. In order to understand the origin of the high Seebeck values and its anomalous temperature dependence in [Cu 2 (1,3-BDT)] n , the electronic band structure has been calculated and the Seebeck coefficient using ab initio DFT calculations determined. In the optimized [Cu 2 (1,3-BDT)] n unit cell, the sulfur atoms are negatively charged, while the copper is in the Cu(I) oxidation state. The electronic band structure along the main symmetry directions is reported in the SI (S15a). A direct band gap of 1.5 eV is extracted from the band structure. This value is underestimated at the current theoretical level, but is sufficient to analyse the nature of the electronic states which mainly contribute to the electrical conductivity and the Seebeck coefficient. Interestingly, the band structure exhibits very flat bands along the directions in momentum space perpendicular to the Cu-S chains, confirming the 1D like character of the electrical conduction. These flattened bands result in peaks when integrated in momentum as shown in the total, eDOS, and projected electron density of states, p-eDOS (see black arrows on Fig. 3a andS15b). Moreover, at the Fermi Please do not adjust margins Please do not adjust margins level, near the top of the valence band, the p-eDOS shows that the dominant character of the propagating electronic states is provided by sulfur and copper orbitals, indicating that there is no conduction between the Cu-S chains, which also confirms that the conduction takes place in the 1D-like Cu-S channels. 24 On the other hand, at the bottom of the conduction band, a more significant contribution of the C-backbone appears, which indicates an inter-chains conduction.

It is well known in the eDOS that sharp features close to the Fermi level result in large values of the Seebeck coefficient. 25 From the transport integrals, the Seebeck is enhanced by sharp and thermally activated electronic singularities in the eDOS. As a direct consequence of the 1D-like character of the valence bands close to the Fermi level, the eDOS of [Cu 2 (1,3-BDT)] n exhibits such sharp peaks (see black arrow on Fig. 3a). This explains the large and abrupt changes of the Seebeck coefficient according to energy close to the Fermi level as seen in Fig. 3b, which can also be observed by varying the temperature as shown in the insert of Fig. 3b (and Fig. S16). This qualitatively explains the high positive Seebeck values observed experimentally and also the maximum observed in the temperature dependence of the Seebeck (insert in Fig. 3b). Thus, the reasons for both the high Seebeck values and a maximum value around 350 K are due the strong quasi-1D character of the electronic transport. This can be seen from the ab-initio simulated electronic band structures (Fig. S16), the flatness of the bands along the Cu-S chains is responsible for peaks (Van Hove singularities) in the eDOS and the temperature dependence of S(T).

Conclusions

A neutral 2D copper-thiolate CP, [Cu 2 (1,3-BDT)] n , has been synthesized under hydrothermal conditions. This CP shows good air stability, due to the strong soft base -soft acid interactions between the thiolate functions and the copper(I) ions, and it remains crystalline upon pressure for pellet formation. This Cu(I)-thiolate network exhibits at RT a mild conductivity of 1.5 mS.cm -1 and a record Seebeck coefficient of 420 µV.K -1 among CPs, equivalent to inorganic solids. The resulting PF is 0.03 µW.m -1 .K -2 . Theoretical studies showed that the conductivity is originating from the 1D inorganic Cu-S chains and explains the high Seebeck coefficient by the presence of a significant peak in the total electronic DOS at the Fermi level, near the top of the valence band. This study highlights that CPs, and more especially copper-thiolates, hold great promises as emerging TE materials. In addition, compared to the highly conducting 2D copperbenzenehexathiolate, which required multistep synthesis of the ligand, 1,3-benzenedithiol is commercially available, reducing considerably the time and cost of the final material production, and allowing formation of large pressed pellets under air. Finally, systematic syntheses of different copperthiolate CPs and in-depth electronic and theoretical studies have to be done in order to increase the electrical conductivity while keeping a good Seebeck coefficient. Changing ligand functionalities, 12b mixing sulfur, using selenium or telluriumbased linkers, 26 carrying out metal doping or synthesising heterometallic composition CPs are some paths to improve the electrical conductivity. 27

Figure 1 .

 1 Figure 1. Structure of [Cu 2 (1,3-BDT)] n : a. the Cu-S chain along the b axis; b. the two deformed Cu 3 S 3 hexagons and c. the central projection along the b axis. Orange, yellow and gray spheres correspond to copper, sulfur and carbon atoms. Hydrogen atoms were omitted for clarity.

Figure 2 .

 2 Figure 2. a. Temperature dependence of the electrical conductivity measured on a [Cu 2 (1,3-BDT)] n pressed pellet by surface four-probe configuration and photo of the device to apply the bulk measurement of the electrical conductivity; b. Seebeck coefficient with temperature (errors come from the linear fits of the raw data) and c. the power factor with temperature.

  journal is © The Royal Society of Chemistry 20xx

Figure 3 .

 3 Figure 3. a. Calculated density of states (DOS) of [Cu 2 (1,3-BDT)] n : projected electronic DOS (p-eDOS) with grey line for C, black for H, orange for S and blue for Cu electrons and the total electronic DOS (eDOS) in red. The insert is a zoom close to the Fermi level and with switched axes. The black arrows show the shard peaks associated to high Seebeck coefficient. b. Calculated Seebeck coefficient at 300 K. The insert shows its temperature dependence at constant energy represented by the dotted red line on the Figure b.

Acknowledgements

The authors acknowledge SOLEIL synchrotron (Gif-sur-Yvette, France) for provision radiation facilities (BAG for X-ray diffraction 20201440) and Erik Elkaïm for assistance in using CRISTAL beamline. This work was supported by the French National Agency (MOTIC ANR-21-CE08-0045). SH acknowledges the CNRS for her PhD grant and AA the Fondation de la Maison de la Chimie for his postdoctoral position. The European Commission is acknowledged by AD for her Marie Skłodowska-Curie Individual Fellowship (101031503 -AniMOC -H2020-MSCA-IF-2020). SP acknowledges V. Giordano and S. Le Floch for their support in using the Seebeck apparatus and the pellet press respectively. The pellets, the thermoelectric and XPS measurements were performed using the facilities at the platforms "PLECE, Transport and PLYRA@ILMTech", respectively.

Author Contributions

CA, SH, AA and JLR performed the synthetic experiments and routine characterizations; AM and NG solved the structure; FP was involved in the purification; SW hosted AD; TN performed the DFT, RD and CA carried out the conductivity; OB executed the XPS; SP and AD conceived the ideas, designed the experiments and wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.


Martin Hamann

Sébastien Andrieux

Markus Schütte

Daniel Telkemeyer

Meik Ranft

Wiebke Drenckhan

Directing the pore size of rigid polyurethane foam via controlled air entrainment

Keywords: PUR foam, air bubbles, nucleation, pore size, morphology, double-syringe mixing

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the wake of increasing environmental awareness, reducing carbon dioxide emission has become a major declared goal of politics and society. Among others, heating and refrigeration significantly contribute to the global annual carbon dioxide output. [START_REF] Nejat | A global review of energy consumption, CO2 emissions and policy in the residential sector (with an overview of the top ten CO2 emitting countries)[END_REF][START_REF] González-Torres | A review on buildings energy information: Trends, end-uses, fuels and drivers[END_REF] In this context, efficient thermal insulation of buildings, refrigerators, storage tanks etc. has become an important topic for both academic research and industry. As a consequence, interest in the optimization of polyurethane rigid (PUR) foams for thermal insulation applications has been revived in the last years, [START_REF] Eling | Current and Future Trends in Polyurethanes: An Industrial Perspective[END_REF][START_REF] Reignier | The morphology of rigid polyurethane foam matrix and its evolution with time during foaming -New insight by cryogenic scanning electron microscopy[END_REF][START_REF] Brondi | The effect of organofluorine additives on the morphology, thermal conductivity and mechanical properties of rigid polyurethane and polyisocyanurate foams[END_REF][START_REF] Brondi | Competing bubble formation mechanisms in rigid polyurethane foaming[END_REF][START_REF] Merillas | Nanoparticles addition in pu foams: The dramatic effect of trapped-air on nucleation[END_REF][START_REF] Liu | Effect of microcrystalline cellulose on the preparation and performance of rigid polyurethane foam[END_REF][START_REF] Omrani | High performance biobased pour-in-place rigid polyurethane foams from facile prepared castor oil-based polyol: Good compatibility with pentane series blowing agent[END_REF][START_REF] Zarzyka | Rigid polyurethane foams modified with borate and oxamide groups -Preparation and properties[END_REF][START_REF] Brondi | Role of Air Bubble Inclusion on Polyurethane Reaction Kinetics[END_REF] since they excel with very low thermal conductivities. [START_REF] Ball | Conductivity Rigid Urethane Foams[END_REF][START_REF] Engels | Polyurethanes: Versatile materials and sustainable problem solvers for today's challenges[END_REF] Key to successful optimization of the physical properties of PUR foams for thermal insulation is their porous morphology [START_REF] Ashby | Mechanical Properties of Cellular Solids[END_REF][START_REF] Schuetz | A basic study of heat transfer through foam insulation[END_REF][START_REF] Mills | Polymer Foams Handbook[END_REF] . For thermally insulating rigid foams in particular, a desirable foam morphology is characterized by small, closed pores with a narrow pore size distribution. [START_REF] Schuetz | A basic study of heat transfer through foam insulation[END_REF][START_REF] Mills | Polymer Foams Handbook[END_REF][START_REF] Gibson | Cellular solids[END_REF][START_REF] Glicksman | Radiation heat transfer in foam insulation[END_REF] To develop strategies allowing for efficient control over the porous foam morphology, it is crucial to understand the main mechanisms involved in PUR foam formation. In general, formulations for the preparation of PUR foams consist of two components. The so-called "Acomponent" consists mainly of polyols but typically also contains additives such as surfactants, catalysts, and blowing agents. [START_REF] Ashida | Polyurethane and Related Foams[END_REF][START_REF] Szycher | Szycher's Handbook of Polyurethanes. 2nd Ed[END_REF] The "B-component" usually only consists of (poly-)isocyanates. Shortly after both components are mixed, a crosslinked PU network via the polyaddition [START_REF] Tieke | Makromolekulare Chemie: Eine Einführung[END_REF] of polyols and (poly-)isocyanates starts to form. In parallel, foam blowing takes place as inert low boiling point blowing agents evaporate due to the strongly exothermic nature of the polyaddition reaction. Besides physical blowing agent evaporation, additional blowing is achieved chemically by adding water which reacts with isocyanates to liberate CO2. [START_REF] Ashida | Polyurethane and Related Foams[END_REF] Once the urethane network formation has sufficiently progressed, the initially liquid matrix becomes solid, hence "freezing" the cellular structure of the blown liquid foam.

Before the onset of PU network formation and foam blowing, micrometric air bubbles are entrapped into the system due to the (pre-)processing procedure, e.g. during the homogenization of the A-component and particularly during the vigorous mixing of the A-and B-component, which results in a bubbly, reactive mixture that eventually forms the PUR foam. The importance of these pre-dispersed air bubbles in fixing the final foam morphology was predicted decades ago. [START_REF] Kanner | Urethane Foam Formation-Role of the Silicone Surfactant[END_REF][START_REF] Baumhäkel | Influence of Stirring Velocity and Air Loading on the Formation of Flexible Urethane Foams[END_REF] It was suggested by Kanner et al. [START_REF] Kanner | Urethane Foam Formation-Role of the Silicone Surfactant[END_REF] that these dispersed air bubbles may serve as "precursors" for the pores of the final PUR foam, since the released blowing agent vapors may diffuse into them without any energy barrier to overcome (type IV non-classical nucleation [START_REF] Jones | Bubble nucleation from gas cavitiesa review[END_REF] ) and that there are typically sufficient air bubbles to account for all the pores found in the final PUR foam. [START_REF] Kanner | Urethane Foam Formation-Role of the Silicone Surfactant[END_REF] Peculiarly, this topic has received little experimental coverage for a long time and related studies were published only in the past few years. [START_REF] Reignier | The morphology of rigid polyurethane foam matrix and its evolution with time during foaming -New insight by cryogenic scanning electron microscopy[END_REF][START_REF] Brondi | Competing bubble formation mechanisms in rigid polyurethane foaming[END_REF][START_REF] Merillas | Nanoparticles addition in pu foams: The dramatic effect of trapped-air on nucleation[END_REF] For instance, Merillas et al. [START_REF] Merillas | Nanoparticles addition in pu foams: The dramatic effect of trapped-air on nucleation[END_REF] demonstrated very recently that the addition of solid nanoclay particles as nucleating agents hardly affects PUR foam morphology if they are dispersed within the system under vacuum.

The authors argued that the vacuum prevents the particles from incorporating tiny air bubbles into the liquid reactive mixture while being dispersed and that only the air bubbles affect PUR foam morphology, not the nanoclay particles themselves. This is in line with findings by Brondi et al. [START_REF] Brondi | Competing bubble formation mechanisms in rigid polyurethane foaming[END_REF] who showed that there are two types of bubble nucleation regimes in nascent PUR foam systems depending on whether dispersed air bubbles are present in the reactive mixture or not.

Without dispersed bubbles, classical homogeneous bubble nucleation caused by gas supersaturation of the liquid phase is observed, while in the presence of dispersed bubbles, the existing bubbles take up all the released vapor of the blowing agent and the CO2 without the nucleation of new bubbles.

If pore formation in PUR foams is in fact mostly driven by the inflation of entrapped air bubbles, the air bubble density within the initially liquid reactive mixture should dictate the pore size of the final PUR foam. To prove this hypothesis, we developed a new double-syringe processing technique which not only allows for controlling the number of micrometric air bubbles that are entrapped into the reactive mixture (Figure 1a and SI Figure S1), but also prevents further air entrainment during the blending step of the A-and B-component. Additionally, the loss of volatile blowing agents during processing is inhibited, since the setup is sealed. The doublesyringe mixing technique as such has already been reported in the literature for the generation of liquid-and solid foams [START_REF] Zowada | Polyfoam: Foam-Templated Microcellular Polymers[END_REF][START_REF] Gaillard | Controlled foam generation using cyclic diphasic flows through a constriction[END_REF] as well as emulsions [START_REF] Nastasa | Moderately stable emulsions produced by a double syringe method[END_REF][START_REF] Soullem | A Novel Low-Pressure Device for Production of Nanoemulsions[END_REF][START_REF] Yu | Development of an Elongational-Flow Microprocess for the Production of Size-Controlled Nanoemulsions : Batch Operation[END_REF] . However, to the best of our knowledge, this study is the first report on the exploitation of this technique in connection with PUR foam systems. Employing this processing and mixing technique together with a simplified model PUR foam formulation, we show that dispersed micrometric air bubbles indeed affect PUR foam morphology crucially and, most importantly, that there is a direct relation between the initial air bubble density in the reactive mixture forming the PUR foam and the mean pore size of the final PUR foam.

Experimental approach, results and discussion

Air entrainment

The first step of our process is the controlled generation of bubbles via air entrainment into the A-component only. For all experiments conducted in this study, the A-component consists of tripropylene glycol (TPG; 93.95 wt%), an amine-based reaction catalyst (1.55 wt% of the Acomponent), a polysiloxane-based surfactant (2.00 wt% of the A-component), water (2.50 wt% of the A-component) as chemical blowing agent, and cyclopentane as physical blowing agent (13.5 wt% with respect to 100% A-component added on top). Detailed information on the PUR foam formulation and preparation is provided in the supporting information (SI, Section 1). As sketched in Figure 1a, 17 ml of the A-component are poured into a 60 ml plastic syringe and attached to another 60 ml plastic syringe filled with 17 ml of air via a narrow silicone tubing (internal diameter dint = 3.2 mm ; external diameter dext = 4.8 mm) and a three-way valve. Then, the connected syringes are fixed within a home-built double-syringe mixing device [START_REF] Gaillard | Controlled foam generation using cyclic diphasic flows through a constriction[END_REF] which drives a periodical motion of a two-armed crankshaft to push the pistons of the connected syringes in an alternating fashion. Thus, the content of both syringes is pushed repeatedly through the narrow tubing which leads to the formation of air bubbles via air entrainment into the polyol-based A-component. A more detailed description of the setup is provided in the supporting information (SI, Section 2). This polyol-based A-component containing air bubbles is what we refer to as "premix" in the following. The double-syringe technique allows for an explicit control of the premix generation for several reasons:

1.

The ratio between liquid and gas in the syringes can be set precisely.

2.

The velocity of the crankshaft movement, i.e. the energy input, can be varied.

3.

The number of mixing cycles, i.e. the mixing time, can be varied.

4.

The system is sealed. Thus, uncontrolled air entrainment from the environment into the system is prevented as well as loss of volatile blowing agents due to evaporation. To alter the number of air bubbles formed during the air entrainment step (see step I in Figure 1b), we varied the number of mixing cycles Nc between air and the A-component while keeping the average speed of the crankshaft arms constant at 〈v〉 = 9 cm s -1 . To be more precise, we chose six different conditions of the air entrainment step with Nc being [0 ; 1 ; 2 ; 5 ; 10 ; 50].

For each Nc, a sample of the premix is collected directly after the air entrainment step and analyzed using transmitted light microscopy. Automatic image analysis via ImageJ is used to measure the number of entrained air bubbles per unit volume of the premix, i.e. the bubble density ωpre, and the mean bubble radius 〈Rb〉. For detailed information on the microscopy analysis, the reader is referred to the supporting information (SI, Section 3.1). An overview over the observed bubble densities ωpre and the mean bubble radii 〈Rb〉 obtained in the premix directly after the air entrainment step for different numbers of air entrainment cycles Nc is provided in Figure 2. Note that for each cycle number Nc, the air entrainment experiments were repeated and analyzed at least three times. The microscopy images (Figure 2a) as well as the diagram in Figure 2b reveal that increasing the number of air entrainment cycles Nc leads to a steady increase of the air bubble density ωpre in the premix. Premixes with air bubble densities ωpre ranging between 0 and roughly 4000 bubbles per mm 3 can hence be obtained in a highly reproducible manner. In contrast, the mean bubble radius 〈Rb〉 is hardly affected by the number of air entrainment cycles Nc since 〈Rb〉 remains constant between 13 -15 µm. These findings confirm that the chosen double-syringe process is very well-suited for studying the relationship between the air bubble density ωpre in the premix (and therefore also in the reactive mixture) prior to PUR foam formation and the final PUR foam morphology, as it provides good control over the bubble formation during air entrainment.

PUR foam preparation

Directly after the air entrainment step, the whole premix is collected in only one of the two syringes (step II in Figure 1b), while the opposite empty syringe is replaced as quickly as possible (typically 40-50 seconds) by a syringe filled with the B-component (oligomeric methylene diphenyl diisocyanate, 198 wt% with respect to A-component; NCO-Index = 120).

Blending of the reactive components, i.e. the now bubbly A-component (premix) and the Bcomponent, is then conducted using three cycles of pushing back and forth both components through the connecting constriction with the double-syringe device. Note that three blending cycles proved to be appropriate to obtain sufficiently homogeneously dispersed "reactive mixtures" of the A-and B-component (see supporting information, Section 3.2), while being short enough to prevent premature foam blowing, i.e. foam blowing inside the syringes during the blending step. As the system is sealed, additional air entrainment and evaporation of volatile components during the blending step of the reactive components is avoided. This is another major advantage of the double-syringe technique, since for most other blending approaches uncontrolled evaporation takes place and air is trapped into the matrix arbitrarily while the reactive compounds are being blended. [START_REF] Baumhäkel | Influence of Stirring Velocity and Air Loading on the Formation of Flexible Urethane Foams[END_REF] The double-syringe process therefore allows to trace back potential changes in PUR foam morphology exclusively to variations of the initial air bubble density ω in the liquid reactive mixture, given that all other process parameters are kept constant (e.g. PUR foam formulation, blending step conditions). After blending, the liquid reactive mixture is ejected via the outlet channel of the three-way valve into a paper cup (V = 735 ml) where foam blowing and solidification take place. Finally, 3-5 days after complete foam solidification, photographs of the whole foam and of cubes cut out of the center of the foam are taken. Moreover, the porous morphology is analyzed via scanning electron microscopy (SEM) (see SI, Section 4). Photographs and SEM images of the PUR foams that were obtained using reactive mixtures of different air bubble density ω are provided in Figure 3. Note that "ω" refers to the air bubble density in the reactive mixture, while "ωpre" refers to the bubble density of the premix. To obtain ω from ωpre, we consider the volumetric dilution of the bubble density in the premix upon addition of the B-component during the blending step, i.e. Figure 3 demonstrates that with increasing number of air entrainment cycles Nc, i.e. with increasing air bubble density ωpre in the premix and, as a consequence, in the liquid reactive mixture, the mean pore radius 〈Rp〉 of the corresponding PUR foam decreases. This is supported by Figure 4a, which plots the mean pore radius 〈Rp〉 as a function of air bubble density ω in the liquid reactive mixture.

𝜔 = 𝑉

A 𝑉 A +𝑉 B 𝜔 pre = 0.38 𝜔 pre . [1]

Correlation between entrained air bubble density and PUR foam pore size

For small ω, the mean pore radius remains close to 1600 µm (colored zone in Figure 4a), whereas it decreases rapidly beyond an air bubble density of ω  250 -300 mm -3 (blank zone in Figure 4a). This behavior is in good agreement with recent findings of Brondi et al. [START_REF] Brondi | Competing bubble formation mechanisms in rigid polyurethane foaming[END_REF] which can be summarized as follows: without dispersed air bubbles in the premix, nucleation and growth in reactive PUR foam systems can be described by type I classical homogeneous nucleation, where gas bubbles nucleate only due to high levels of gas supersaturation in the liquid phase. [START_REF] Jones | Bubble nucleation from gas cavitiesa review[END_REF] Thus, the size of the nucleated bubbles and their growth are determined by physico-chemical properties of the system such as reaction heat and gas supersaturation. In contrast, in the presence of pre-dispersed air bubbles, the released blowing agent vapor has no energy barrier to be overcome in order to nucleate a gas bubble and can therefore readily diffuse into the pre-existing air bubbles (type IV non-classical nucleation). [START_REF] Jones | Bubble nucleation from gas cavitiesa review[END_REF] In analogy, we argue that mostly classical homogeneous nucleation due to gas supersaturation (type I) is observed for very low ω (colored zone in Figure 4a). Increasing the air bubble density ω, type I classical nucleation and type IV non-classical nucleation are likely to happen in parallel. However, surpassing this "transition regime", type IV non-classical nucleation becomes the dominating mechanism (blank zone in Figure 4a). In this regime, where the pre-dispersed air bubbles dominate pore formation, a direct link between the air bubble density ω in the liquid reactive mixture and the final PUR foam size is conceivable. Assuming that all the bubbles entrained into the liquid reactive mixture are preserved during PUR foam blowing and take up the liberated blowing agent volume quantitatively to form the pores of the PUR foam, a minimum possible pore size should therefore be accessible. Thus, to estimate a lower bound for the final pore size, we propose a model based on the following approximations:

1.

All the released blowing agent volume (cyclopentane vapor + CO2) is taken up quantitatively by the pre-dispersed air bubbles.

2.

The total volume of the pre-dispersed air bubbles is negligible with respect to the total liberated blowing agent volume, i.e. the mean volume 〈Vb〉 of entrained bubbles is negligible with respect to final, mean pore volume 〈Vp〉 of the foam (〈Vb〉 ≪ 〈Vp〉).

3.

The number of air bubbles is preserved as soon as PU foam blowing and PU network formation begin (Nb = const).

4.

The density of the initially liquid reactive mixture is approximately the same as the density of the polyurethane polymer that forms the solid continuousphase of the foam.

Hence, we neglect the expansion or shrinkage of the fcontinuous phase of the foam during polymerization (Vlrm ≈ Vpol). We estimate that this introduces an error of maximally 10%, which is of the order of other experimental errors.

We can then establish the relationship

1 〈𝑅 p 3 〉 = 𝑐 • 𝜔, [2]
with c being

𝑐 = 4 3 𝜋 • 1 -𝜑 𝜑 . [3]
where 〈Rp 3 〉 is a measure of the mean pore volume and φ is the porosity of the PUR foam. Note that a detailed derivation of Equations [2] and [3] is provided in the supporting information (SI Section 5). This simple model predicts an inverse relation between the final pore volume 〈Rp 3 〉 and the initial air bubble density ω in the liquid reactive premix. The blue points correspond to the raw data, whereas the pink diamonds take into account a correction for bubble ageing during the blending step. The red line represents the proposed model given in Eq. [2].

To check the validity of the model introduced in Eq. [2], we replotted Figure 4a using 1/〈Rp 3 〉.

The resulting diagram is shown in Figure 4b. The red line represents Eq. [2], while the blue data points correspond to the raw experimental data. While both show the same trend predicted by Eq. [2], i.e. 1/〈Rp 3 〉 ~ ω, they also show a strong discrepancy in the pre-factor. However, this first comparison overlooks a crucial issue: rapid bubble aging through bubble dissolution and

Ostwald ripening is affecting the bubble density ω as soon as the air entrainment step is stopped until foam blowing is triggered. Coalescence between air bubbles in the premix is considered negligible at this stage due to the small probability of bubble/bubble encounters thanks to their low density. Using microscopy of the bubbly premix, we established quantitative estimations of the bubble ageing, which are presented in more detail in the Supplementary Information (SI, Figure S4 and S5). We found that roughly 39% of the air bubbles present in the premix directly after the air entrainment step (see Figure 2) vanish within first 100 seconds due bubble dissolution and gas exchange, 100 s being approximately the delay between the air entrainment step and the beginning of PU foam blowing. We therefore apply this correction for bubble ageing to all the experimental raw data, leading to the pink diamonds shown in Figure 4b.

As can be seen, the slope of line fitted through the ageing-corrected data in the bubbledominated regime is closer to the slope of the red line representing the model proposed in Eq.

[2]. However, the difference in slope remains a factor of 5: while the slope of Eq. [2] is predicted as 11 x 10 -11 µm -3 /mm -3 , we find 2.2 x 10 -11 µm -3 /mm -3 for the experiment. An explanation for this result could be that we underestimated the extent of bubble ageing prior to the onset of PU foam blowing. This may be due to the fact that the bubble ageing was studied in the premix, i.e. in the A-component only, while the full bubbly reactive mixture also contains 62 vol% of B-component into which air can dissolve. Unfortunately, studying this bubbly reactive mixture via transmitted light microscopy is complicated due to fast reaction kinetics. Moreover, our bubble ageing studies were conducted under static conditions in a quasi-two-dimensional setup (SI Section 3.1) which is very different from the situation present in the double-syringe device.

Measuring the air bubble density in-situ during double-syringe processing remains a challenge that we have yet to address.

Other than that, our proposed model is based on simple, idealized approximations that may not fully match with the experimental reality. For instance, we calculated that the volume of the reactive system should expand by a factor of ~ 40 upon PUR foam blowing given that the blowing agent vapor and CO2 are quantitatively integrated into the PUR foam. However, experimentally we find expansion factors between 35 -37 in our obtained foams (SI Table S2).

Last but not least, our model does not take into account potential coarsening of the bubbles while the reactive system expands due to foam blowing. In the same context, coalescence of inflating air bubbles being pushed against each other during foam blowing cannot be ruled out.

The monitoring and quantification of both coarsening and coalescence in-situ in reacting systems is a delicate challenge that needs further investigation.

Another interesting observations is that the experimentally found line shows an offset on the bubble density-axis unlike the model proposed in Eq. [2]. We hypothesize that this is due to the physical dissolution of a fixed amount of air bubbles in the isocyanate during the blending step of the premix and the B-component since air is highly soluble in isocyanate. Again, quantifying the actual extent of air bubble dissolution in the B-component during blending remains a challenge due to the hazardous nature of the isocyanate and the fast reaction kinetics. Moreover, a certain proportion of the liberated blowing agent vapor volume might be unavailable for air bubble inflation since a certain solubility of cyclopentane vapor and CO2 in both the A-and the B-component cannot be denied [START_REF] Caprio | Polyether polyol/CO2 solutions: Solubility, mutual diffusivity, specific volume and interfacial tension by coupled gravimetry-Axisymmetric Drop Shape Analysis[END_REF][START_REF] Caprio | Mass transport and physical properties of polymeric methylene diphenyl diisocyanate/CO2 solutions[END_REF] .

Due to these different deviations between the experimental reality and the idealized assumptions made in our model, Eq. [2] should therefore be taken as a lower limit for the estimation of the final pore size. However, we emphasize that our data shows that the air bubble density ω in the liquid reactive mixture can in fact be used to alter the mean pore radius 〈Rp〉 of PUR foams and that, as in the model, there is an inverse relation between the mean pore volume 1/〈Rp 3 〉 and the air bubble density ω in the air bubble-dominated nucleation regime.

Conclusion

In summary, we describe for the first time the use of a double-syringe mixing device for the preparation of PUR foams. The main advantage of this processing approach is that the entrainment of micrometric air bubbles and the mixing of the reactive components can be systematically decoupled. Moreover, the double-syringe processing technique allows for a controlled variation of the air bubble density in the liquid reactive mixture prior to the onset of foam blowing and PU network formation. We profit from these advantages to analyze the relationship between the air bubble density in the initially liquid reactive mixture and the pore size of the final PUR foam. Our results confirm that the bubble nucleation and growth mechanisms fundamentally depend on the amount of pre-existing air bubbles in the reacting PUR foam system. In absence of dispersed air bubbles (and for very low air bubble densities) in the liquid reactive mixture, bubble nucleation and growth are dominated by supersaturation of the liquid phase with released blowing agent vapor (type I classical homogeneous nucleation). In this case very coarse PUR foams are obtained and the pore size is independent of the initial air bubble density. For high air bubble densities in the liquid reactive mixture, bubble growth is dominated by the inflation of pre-existing air bubbles which serve as reservoirs for released blowing agent vapor which can diffuse into the reservoirs without surpassing any energy barrier (type IV non-classical nucleation). We demonstrated that in the air bubbledominated regime, there is indeed a distinct correlation between the air bubble density in the liquid reactive mixture and the mean pore radius of the final PUR foam. For instance, we could reduce the mean pore radius of the foam from 1600 µm to 400 µm only by increasing the amount of pre-dispersed air bubbles within the system. Moreover, we proposed a simple model which predicts the experimentally observed inverse proportionality between the mean pore volume and the air bubble density in the liquid reactive mixture. However, to reach full quantitative agreement between the model and the experiment, all contributions that lead to air bubble and foam ageing need to be addressed quantitatively in future studies. To conclude, we demonstrate how the pore size of PUR foams can be controlled simply via air bubble inclusion into the system prior to foam blowing and PU network formation. These insights highlight the crucial impact of pre-dispersed air bubbles on PUR foam morphology and suggest a new path for foam optimization that focusses primarily on maximizing air bubble entrainmentrather than on the formulationto reduce the final pore size for optimized thermal insulation. We are confident that the presented approach may not only provide a valuable model experiment to scan formulations in R&D laboratories, but it may also provide concrete suggestions as to how to improve current industrial processes to reach smaller pore sizes at identical formulation for better thermal insulation.

Figure 1 :

 1 Figure 1: (a) Schematic drawing of the double-syringe device used for bubble generation via air entrainment and blending of reactive PUR foam formulations. (b) Three-step process used to generate the premix and the liquid reactive mixture with a controlled number of micrometric air bubbles before foam blowing.

Figure 2 :

 2 Figure 2: (a) Transmitted light microscopy photographs of bubbly, polyol-based Acomponents (premix) directly after the air entrainment step for different numbers of air entrainment cycles Nc. (b) Measured air bubble density in the premix ωpre after air entrainment (green diamonds) and the corresponding mean bubble radii 〈Rb〉 (blue circles) for different numbers of air entrainment cycles Nc in the double-syringe device.

Figure 3 :

 3 Figure 3: Overview of PUR foams obtained using the double-syringe process with increasing number of air entrainment cycles. (a) Images of the entire cup foams. (b) Images of cubic foam samples cut out of the center of the cup foams. (c) Scanning electron microscopy (SEM) images of slices cut out of the cubic foam samples. Note that for the foams corresponding to Nc = 0 and 1, the SEM images are missing, since the pores were too large. For the case of Nc = 0 and 1, the photographs shown in (b1) and (b2) were used to determine the mean pore size.

Figure 4 :

 4 Figure 4: (a) Mean pore radius 〈Rp〉 of the obtained solid PUR foams versus initial air bubble density ω in the liquid reactive mixture. The green labels represent the number of cycles that was used for the air entrainment step respectively. The colored area represents the regime of homogeneous bubble nucleation, whereas the clear area stands for the non-classical nuclearion regime dominated by the pre-dispersed air bubbles. (b) Mean inverse pore volume 1/〈Rp 3 〉 of solid PUR foams versus the initial air bubble density ω in the corresponding reactive mixture.The blue points correspond to the raw data, whereas the pink diamonds take into account a correction for bubble ageing during the blending step. The red line represents the proposed model given in Eq.[2].

Acknowledgement

The authors gratefully acknowledge Marc Bottinot, Vincent Klein, Jérémy Sanchez, and Rémy Brauge at the Laboratoire de Physique des Solides of the University Paris-Saclay for building the double-syringe device used in this study. The authors would also like to thank Marie-Pierre Krafft for discussions concerning the interaction of fluorocarbons with interfaces, and Manon Jouanlanne, Aurélie Hourlier-Fargette, and Leandro for fruitful exchanges concerning PU formulations and experimental details. Last but not least, the authors gratefully acknowledge the PLAMICS facility of the Institute Charles Sadron for providing access to the SEM. supported by IdEx Unistra (ANR-10-IDEX-0002) and SFRI (STRATUS project, ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program. Moreover, the authors gratefully acknowledge the financial support granted by BASF SE for this study.

ORCID ID

Martin Hamann https://orcid.org/0000-0001-8573-6217

Funding

This work has been financed by an ERC Consolidator Grant (agreement 819511 -METAFOAM) and profited also from an IdEx Unistra "Attractivity grant" (Chaire W. Drenckhan). It was conducted in the framework of the Interdisciplinary Institute HiFunMat, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was

Conflict of interest

The authors declare that there is no conflict of interest.

Supporting information

Detailed supporting information for this article is available online.