Validation of Milner's visco-elastic theory of sintering for the generation of porous polymers with finely tuned morphology.

Sacrificial sphere templating has become a method of choice to generate macro-porous materials with well-defined, interconnected pores. For this purpose, the interstices of a sphere packing are filled with a solidifying matrix, from which the spheres are subsequently removed to obtain interconnected voids. In order to control the size of the interconnections, viscous sintering of the initial sphere template has proven a reliable approach. To predict how the interconnections evolve with different sintering parameters, such as time or temperature, Frenkel's model has been used with reasonable success over the last 70 years. However, numerous investigations have shown that the often complex flow behaviour of the spheres needs to be taken into account. To this end, S. Milner [arXiv:1907.05862] developed recently a theoretical model which improves on some key assumptions made in Frenkel's model, leading to a slightly different scaling. He also extended this new model to take into account the visco-elastic response of the spheres. Using an in-depth investigation of templates of paraffin spheres, we provide here the first systematic comparison with Milner's theory. Firstly, we show that his new scaling describes the experimental data slightly better than Frenkel's scaling. We then show that the visco-elastic version of his model provides a significantly improved description of the data over a wide parameter range. We finally use the obtained sphere templates to produce macro-porous polyurethanes with finely controlled pore and interconnection sizes. The general applicability of Milner's theory makes it transferable to a wide range of formulations, provided the flow properties of the sphere material can be quantified. It therefore provides a powerful tool to guide the creation of sphere packings and porous materials with finely controlled morphologies.


1) Introduction
Macro-porous materials with pore sizes larger than 10 µm are employed in many fields [1][2][3][4] , ranging from catalysis, absorption, acoustics, heat transfer to tissue engineering applications. This widespread use is motivated by the many different types of macro-porous materials which can be manufactured in terms of their porous structure (open-vs closed-cell, low-vs high density, etc.) 5,6 and their base materials (metals, ceramics, synthetic or natural polymers, etc.) [7][8][9] . Key parameters of the porous structure are the dimensions of the pores and of their interconnections. Their impact on the material properties has been put in evidence in different fields. For example, in acoustics, Trinh et al. 10 and Jahani et al. 11 showed how the number and size of the interconnections relates to the sound absorption coefficient. In catalysis, a high degree of interconnectivity helps to reduce the pressure drop in reactors 12 . Langlois et al. 13 demonstrated how the number and the size of the interconnections are related to the flow resistivity. Last but not least, in tissue engineering it has been shown that pore and interconnection sizes impact significantly the behaviour of cells 14,15 , fibrotic encapsulation and vascularisation within the body [16][17][18] . Of particular interest to this article are macro-porous polymeric materials. The explicit control of their pore and interconnection sizes remains an important challenge. Liquid foams or emulsions may be used as templates for this purpose [19][20][21][22] . The pore dimensions can then be controlled via the choice of the foaming/emulsification methods 23 and/or via the formulation 24 . In particular, microfluidic techniques allow to generate porous materials with highly monodisperse and even periodic pores 25,26 . However, since the pore opening mechanisms are badly understood 27,28 , it remains a major challenge to control the stability of the initially liquid template and to tune explicitly the presence and size of the interconnections. Moreover, the necessity to stabilise the initially liquid template of two immiscible fluids puts important constraints on possible formulations. One way to overcome these challenges is the so-called sphere templating approach [29][30][31] (Figure 1), which is is based on the use of solid spherical particles which are tightly packed in a mould. The interstices between the spheres are infiltrated by an initially liquid matrix (monomers, polymer solutions or melts, etc.), which is then solidified through polymerisation and/or cross-linking. Once the matrix is solidified, the spheres are selectively dissolved in order to leave an interconnected network of pores. In this process, the final pore diameter is given by the diameter of the sacrificial spheres, while the interconnections are obtained by sintering of the initial sphere template. The latter is obtained by heating the sphere template to temperatures which allow for a material transport from the spheres to the interconnections 32 , without melting the entire structure. This material transport is driven by the surface energy of the spheres, which is reduced by a progressive filling of the interconnections creating a growing "neck" between the spheres. The material transport can arise via different phenomena 32 . Of interest here are those which result from viscous flow. Such sphere templating approaches have been successfully employed in the past for the generation of macro-porous polymeric materials. For example, Somo et al. 33 and Chen et al. 34 used Poly (methyl methacrylate) PMMA spheres to produce scaffolds with micrometric pore size (20 -110 μm) for tissue engineering purposes. The spheres were sintered for variable times (0 -30 h) and temperatures (110 -175°C) to tune the interconnection sizes. Paraffin is also regularly used as a porogen, since spherical paraffin beads can be easily produced via quenching of an emulsion 35,36 . For example, Grenier et al. 37 produced a polyurethane scaffold using paraffin as sacrificial template, but without sintering. Ma et al. 30 sintered paraffin spheres at 37°C for 20 min to generate a porous scaffold. Takagi et al. 38 sintered polyethylene spheres to form the template of a porous bio-ceramic. Alternative approaches to sintering for the control of the interconnections may be used. For example, Zhao et al. 39 used Poly(styrene-co-divinylbenzene) spheres as porogens, creating the neck (i.e. interconnections) by addition of an adhesive which accumulated in the contact zone of the spheres. In another case, Descamps et al. 40 created the neck between PMMA spheres using acetone which slightly dissolved the spheres creating a neck at their contact. Even if sphere templating via sintering is regularly used to obtain macro-porous polymers, the detailed mechanism of the formation of the interconnections has not been studied in a sufficiently systematic manner to profit from the predictive power of an accompanying model. As discussed in more detail in Section 2, for systems where the relaxation dynamics can be described by a viscous flow, the first theoretical description of the phenomenon dates back to a theory proposed by Frenkel 41 . His model has been used quite successfully to account for a wide range of experimental results 36,42,43 . Some modifications of the existing theories were made by Bellehumeur et al. 44 and Mazur et al. 45 to take into account the viscoelastic properties of polymeric spheres which were compared with experiments. Recently, Frenkel's model has been challenged by Milner 46,47 , who established an analogy between elastic Hertzian contacts and the early stages of viscous sintering. Milner showed that the characteristic size of the flow in the neck region is incorrectly estimated in Frenkel's model and that the resulting correction leads to a slightly different scaling for the neck evolution (Section 2). This improved model was shown to describe successfully the flattening of polystyrene latex microsphere on glass surfaces 47 . Moreover, Milner extended his analysis to include predictions for the sintering of visco-elastic spheres, for which the flow cannot be considered as purely viscous. Goal of this article is to compare for the first time the scalings proposed by Frenkel and Milner by providing a detailed quantitative investigations of viscous sintering to generate porous polymeric structures with explicitly controlled pore and interconnection diameters. As template system we use spheres of paraffin with radii R in the range of 33 -125 µm radius (top row of Figure 1) which we generate via an emulsification process (Section 3). We show that for the entire range of the used control parameters, the associated mechanisms are captured slightly better by Milner's "viscous model" (Section 2) than by Frenkel's. Nevertheless, the experimental data remains unsatisfyingly scattered with respect to the theoretical prediction. We finally show that this scatter can be drastically reduced to a very satisfying comparison using Milner's "viscoelastic model" which takes correctly into account the Non-Newtonian response of the paraffin. By replicating the sintered sphere templates, we obtain macroporous polyurethanes (bottom row of Figure 1) with finely controlled pore and interconnection sizes. Since Milner's model can be readily transferred to a wide range of polymers as long as their temperature-dependent surface and visco-elastic properties are known, it will provide the reader with a powerful tool for the predictive design of macro-porous polymers.

2.1) Frenkel's and Milner's Newtonian scalings
The sintering process of spheres in contact is driven by the surface energy of the sphere packing which is proportional to its surface area. An important reduction of this surface area is Please do not adjust margins Please do not adjust margins achieved by progressively filling the contact zones between the spheres. This filling process can arise in many different ways 32 , ranging from solid-state diffusion of the sphere surface to viscous flow of the contact zone. The sintering process of spheres of paraffin (mixture of alkanes, iso-alkanes or cycloalkanes 48 ) at work in our experiments is well described by a viscous flow in the vicinity of the contact area of the spheres induced by surface forces. A comprehensive theoretical treatment of the complete dynamics is highly complex due to the non-trivial intermediate states of the system. A thorough and illuminating treatment of the early stages of sintering has been given recently by Milner 46,47 . Here we shall only provide a short summary of the main arguments. An initial dimensional analysis of the problem shows that for a Newtonian fluid of viscosity , the evolution of the radius ( ) of the circular, symmetric contact between two spheres of radius R (see top of Figure 1) must fulfil the general scaling Here, t is the time and the surface tension of the spheres. The mass density does not enter since the sintering processes are sufficiently slow to neglect inertial effects. Additionally, the considered sample heights are sufficiently small so that the pressure force induced by gravitation can be neglected with respect to surface forces. Note that provides a characteristic sintering time which highlights the important influence of the sphere radius: the smaller the spheres, the faster the sintering process.
In the regime of small deformations ( ( )⁄ <<1), Equ. (1) is expected to be well approximated by a power law 46 where C is a numerical constant. First quantitative expressions for Equ. (3) were provided by Frenkel in 1945 41 . Using a scaling argument, he predicted that = 1 2 ⁄ with = 4. The resulting expression has been used quite successfully in the past to describe experimental observations 36,42,43 . However, as Milner pointed out, his scaling does not capture properly the underlying mechanism. In the following we shall recall the main arguments. Sintering dynamics is controlled by the balance between the rate of change of the surface energy 49 , and the dissipation rate Where is the velocity field of the sphere with the indices , indicating the Cartesian coordinates x,y,z, implying implicit summation over and and is the elemental volume. In the early stage of sintering, one can approximate the geometry by overlapping spheres, as sketched in the top row of Figure 1. Milner noticed that the typical strain rate ˙≈ scales as ḣ ⁄ , where [2 − ℎ( )] is the distance between the centres of the spheres at time t. This gives simply ℎ ≈ 2 ⁄ . The key difference between Milner's and Frenkel's approach is to realise that the flow extends along the symmetry axis only up to a characteristic length (and not ), as shown in Figure  2. This is justified by the fact that at early times, the only relevant length scale of the flow, described by 2 , is given by the boundary extension, namely by . The typical volume in which the flow develops is therefore proportional to 3 , as in the Hertz contact problem 50 , and not proportional to 3 , as assumed by Frenkel. One can therefore approximate Equ.
After integration this gives the scaling ( ) 3 ∝ .
This leads to the dimensional form provided in Equ.

2.2) « Crushing » vs « smoothing » flow
It is important to keep in mind that the growth regime described in Section 2.1, which we may call the "crushing flow", is expected to dominate only for small deformations. During this "crushing flow" the velocity field is purely normal to the contact surface between the sphere and the solid, as sketched in Figure  2 (left). For larger deformations, however, a secondary "smoothing flow" takes over, which is sketched in Figure 2 (right). The smoothing flow is driven by the fact that the contact zone between the two spheres is cusp-like, with very small radii of curvature leading to an unbalanced Laplace pressure drop. As a result, the surface tension also drives a radial flow in the neck region, which becomes increasingly important. Since both processes are described by different scalings, it is important to know at what point of the sintering process the change of regime arises. For the derivation 46 , let us consider the geometry sketched in Figure 3. Two spheres of radius are in contact with a cusp of radius centred around . Let us assume early stages of sintering so that << . Then, the external boundary of the contact area between two sintering spheres is a circle whose radius can be approximated as ( ) + ( )  ( ). In the smoothing flow regime, the fluid is flowing radially, with a component tangential to the contact line. In the cross section, this flow is approximately a two-dimentional sink flow (red dashed arrows in Figure 3), with material being conveyed to fill the gap between the spheres. This approximation is justified in the early times when the angle between the two spheres at the boundary is small. The two-dimensional incompressible flow has a well-known velocity field, ( ′) = − 0 ( ) ( ) 2 ′ ( ′) 2 where 0 ( ) is the velocity of the smoothed cusp and is taken from the sink centre. Notice in his formula, that the radius of the smoothed cusp has been approximated by ( )/2 where ( ) is defined in Figure 3. This parameter is geometrically related to ( ) and provided that ( )/ << 1, one has δ(t) ~ r (t)²/R.
I.e. the rate of change of r is directly proportional to 0 . A second dynamical relation is obtained by equating, as in Section 2.1, the rate of change in surface energy ˙s urf (Equ. (4)) with the adapted dissipation rate ˙ (Equ. (5)). For the radial flow near the circular contour line one obtains Equating Equ. (10) with Equ. (4) using Equ. (9) one therefore obtains an expression for the growth rate of in the smoothing regime Thus, unlike in the crushing regime (Equ. (6)), ( ) becomes independent of time in the smoothing regime and depends only on the ratio / . Consequently, after integration, the interconnection radius r is predicted to depend linearly on time, i.e.
To obtain a quantitative estimation of the crossover time Please do not adjust margins Please do not adjust margins Using in Equ. (8), one can predict that at the transition between the two regimes, the crushing flow will have grown the contact surface to a radius = 1 4 = 0.25.
It is interesting to note that at the crossover time , the normalised radius ( )/ predicted by the smoothing flow using Equ. (12) flow is only ~ 0,08. This means that the effective validity of the approximation where only the crushing flow is taken into account, goes probably farther than the value given in Equ. (14).

2.3) Milner's visco-elastic model for paraffin spheres
In his paper 46 , Milner extended the viscous scaling of Equ. (8) to the sintering of viscoelastic liquids, leading to the prediction .
Here ( ) is the creep compliance of the visco-elastic matrix which constitutes the spheres 51 . Note that the asymptotic regime of ( ) for very long times is ( ) ~ / . Therefore, Equ. (15) is in full accordance with its viscous analogue given in Equ. (8) in the crushing time regime. Interestingly, the dependency of ( ) on ( ) 1/3 was already found by Mazur et al. 52 in the 90s. In order to compare the visco-elastic model with experimental data, it is important to be able to express the creep compliance Defining the Laplace transform of ( ) as one notices that ( ) = ℜ[̂]( = − ). As shown in Section 3.7 (Figure 8), in the experimentally accessible frequency range (10 -2 Hz < ω < 10 Hz), ( ) of our paraffin is convincingly described by a simple powerlaw ( ) ≈ − with and being temperature-dependent phenomenological coefficients. It is worth stressing that ∈ ]0,1[. Since ( ) does not show a tendency to level off in the low frequency range of the measurements (10 -2 Hz), we assume here that we can extrapolate the simple power law all the way into the frequency range relevant for our experiments, which ranges from minutes to a day, i.e. 10 -5 Hz < ω < 10 -2 Hz. In this case, an analytical expression of ( ) can be obtained because the inverse cosine-Fourier transform of − is exactly known to be 53 ( ) then follows from the fact that ( ) and ( ) are exactly related via ̂( )̂( ) = 1/ ².
One therefore finds ̂( ) = sin ( 2 ) /( 2− ), whose inverse Laplace transform is known 53 Hence using the fit parameters ν and ξ from the rheology experiments, we avail of a direct expression for the creep compliance ( ).
The visco-elastic model of Equ. (15) for our experiments therefore becomes

3.1) Generation of paraffin spheres
We generate the paraffin spheres by a dispersion method inspired by Ma et al. 30 . We use melting paraffin with a molecular weight of 341.451 g/mol from Fischer Scientific (CAS 8002-74-2). Its melting point is given at 58-62°C, which is confirmed by our DSC measurements shown in Section 3.6. To generate the spheres, we add 10 g of paraffin to a 400 mL solution of 70 °C osmosed water with 3 g of Polyvinyl alcohol (PVA) (Sigma, CAS 9002-89-5, MW 31 000-50 000 89% hydrolysate). The mixture is kept in an 800 mL beaker and is vigorously stirred with a magnet (35 mm length) for 10 min at different rotation velocities (400, 700, 950 and 1200 rpm).
Higher rotation speeds of 3200 and 4200 rpm are obtained using an ultra-Turrax (IKA ®T25). Since the paraffin is liquid at the chosen temperature, the vigorous stirring creates a paraffin-in-oil emulsion whose drops are stabilised against coalescence by the PVA. With increasing rotation speed, the paraffin is broken into increasingly smaller drops, leading to smaller spheres after solidification. We quench the emulsion by the addition of 400 mL of 4 °C osmosed water, leading to nearly instantaneous solidification of the liquid drops into spherical beads, as shown in the central row of Figure 1 and Figure 4a. We noticed that with decreasing bead size, the surface of the spheres roughens. This may be due to the increasing surface-to-volume ratio and potentially associated shrinkage phenomena since for smaller spheres, the surface-to-volume ratio increases dramatically. If the surface shrinks differently at the core of the particle during cooling/solidification, this is greatly amplified and can lead to buckling stresses.
Using numerical microscopy combined with image analysis (Section 3.2) we determine the size distributions of the obtained paraffin spheres. Figure 4d shows examples of size distributions obtained for two different rotation speeds (400, 950 rpm), putting clearly in evidence the systematic decrease of the sphere radius with rotation speed. Since the distributions are quite polydisperse, we reduce the polydispersity using

3.2) Determination of the sphere size distribution
We use the numerical microscope Keyence VHX-5000 to obtain the sphere size distributions. For the imaging, we disperse about 1 mg of dried paraffin spheres on a microscope slide and place it on the transparent tray of the microscope. Images are recorded with the transmitted lighting mode. For the image analysis, we use the software supplied with the microscope. First we apply the "brightness" mode to select the desired image area and choose the most appropriate threshold level which allows to identify (Figure 4a) and isolate spheres which are stuck together (Figure 4b). We further apply the "separation tool" to separate spheres contained in clusters, as shown in Figure 4c. After detecting at least 800 spheres, the size distributions are expressed as the probability density function defined by Here, is the measured sphere radius and ( < < + ∆ ) is the number of spheres having a radius between and + ∆ , where ∆ is the bin width of the frequency count histogram. N stands for the total number of particles. We use OriginPro 8 software for the data treatment. The obtained PDFs are then fitted to Gaussian distributions, as shown in Figure 4df.

3.3) Sintering procedure
We add 3 g of paraffin spheres with given average radius 〈 〉 (Section 3.1) to a petri dish (34 mm in diameter, VWR) and pack them manually by tapping. Small holes are drilled into the bottom and on the sides of the Petri dish beforehand to help air removal upon polymer filling. We initially close the holes from the outside with parafilm to prevent the loss of paraffin spheres. The Petri dishes containing the paraffin spheres are then put in an oven (Memmert) at temperatures ranging from 35°C to 42°C for sintering times up to 18 hours. We pre-heated the oven for 24 hours to ensure stable temperatures. Once the required sintering time is reached, the mold is cooled down at room temperature and we removed the parafilm. Using a microprocessor thermometer (HANNA instruments, HI 8757), we characterised the evolution of the temperature in the oven and at the center of the sphere template. As shown in the example of Figure 7, the fluctuations of the oven temperature are of the order of 0.5 °C and it takes about 60 min for the sphere packing to reach the temperature of the oven. In order to avoid misinterpretation of the data due to this initial equilibration time, we characterised the template morphology for samples which spent at least 2 hours in the oven.

3.4) Generation of macro-porous polyurethane
All polyurethane precursors are kindly supplied by FoamPartner (Gontenschwil, Swiss). We use the polyether-based triol Voranol 6150 and the diisocyanate M220 (both from DOW chemical) at a weight ratio of 13.4 to 1 to ensure an isocyanate index of 100. The polyurethane is obtained by simply mixing both components (Ultra Turrax, IKA ®T25 at 20K rpm, for 2 minutes). No catalyst is added to let enough time for the liquid mixture to infiltrate the sphere packing before solidification and to maintain a low reaction temperature to avoid additional sphere deformation. The mixture is poured directly onto the sintered sphere packing where it spontaneously fills the interstitial spaces between the beads. This filling process takes about 1 h. The sample is then left at room temperature for 72 h to ensure full solidification of the polyurethane. Finally, we use Soxhlet extraction by n-hexane for 6 h, at 100°C in order to remove the paraffin spheres. The resulting porous polyurethanes are let to dry overnight under a fume hood at room temperature.

3.5) Scanning electron microscopy (SEM)
We image the sphere templates and the porous polyurethanes with an SEM Hitachi SU8010 (FEG-SEM) under 1 KeV voltage acceleration to investigate the pore and interconnection morphologies. The samples are imaged directly without metallisation. Before observation of the paraffin sphere templates, they are quenched in liquid nitrogen at -196 °C and then directly broken. This helps to conserve the template structure and, in particular, the geometry of the necks. The PU sponges are cut at room temperature using a cutter knife and imaged directly (Figure 6a). The obtained images are treated with ImageJ software. In the case of the sphere templates, the pore radius and interconnection radius were measured by hand, taking into account at least 130 pores and interconnections for each data point to ensure good statistic.
In the case of the PU sponges, the interconnection sizes are calculated semi-automatically by employing the « particle analysis » tool (Figure 6c) after proper thresholding (Figure 6b) of the initial image (Figure 6a). In order to exclude zones which are not interconnections, only objects with a circularity higher than 0.5 are selected. We also only considered "particles" whose area ranges between 60 and 6000 μm². Since the circular interconnections are imaged at different angles, they commonly appear as ellipses. We use the measurement of half of the major axis of the ellipse (Frenet diametre), which corresponds to the real radius of the interconnections. For each data point we measure at least 200 interconnections to ensure good statistics.
The average values of the sphere radii   and the interconnection radii   are obtained from the arithmetic mean of the distributions, while we use the standard deviation as the error of one sample. The error bars given in the unscaled graphs of Section 4.2 are obtained by repeating each experiment three times. The error bar corresponds to the standard deviation calculated for the overall set of values. The error bars in scaled graphs (such as Figure 11) are obtained by linear error propagation. If we plot the function = / 1− its relative error is given by where ∆ and ∆ are the errors of and , respectively. Please do not adjust margins Please do not adjust margins

3.6) DSC analysis
We place 8.2 mg of the paraffin in a stainless steel capsule and put the capsule inside the instrument (PerkinElmer DSC 8500). Experiments were done under nitrogen atmosphere which allow to ensure a chemically inert environment. Figure 7 shows the curves obtained by heating/cooling between 25 °C and 80°C with a ramp of 10 °C/min. A first cycle has also been performed to erase the thermal history of the paraffin (not shown here).
The endothermic peak at 55-67°C corresponds to the melting point which is in accordance with the range given by the supplier. However, one can also see a slight shoulder in the curve around 40°C which indicates that some paraffin chains may have started to melt already before the melting point. All our experiments were typically conducted between 35 and 42°C, thus, even if our working temperatures were below the melting transition, due to the polydisperse nature of paraffin, we expected that some hydrocarbon chains may have melted at these temperatures.

3.7) Rheology of paraffin
Experiments were conducted with a rheometer MCR 702 MultiDrive® from Anton Paar at various temperatures (30, 32, 34, 36, 38, 40, 42 °C). The cone plate configuration with an angle of 5° was used. Before the tests, platelets of paraffin wax were deposited onto the peltier plate and heated up to 60 °C to transform the paraffin into its liquid state. The rheometre cone was then approached until having full contact with the wax and down to the desired gap size. The paraffin was then cooled down at 20 °C/min until 25 °C to form the solid sample between the plateaux. Afterwards, the solidified sample is heated at 5 °C/min to the desired temperature to simulate the conditions in the oven. The rheological measurements are started once the sample has reached thermal equilibrium. The dynamic viscosity of the paraffin was obtained with a frequency sweep from 0.01 up to 10 Hz (50 points log spaced) at a constant strain of 0,1 %.
he resulting  ( ) curves are shown in Figure 8. They are well fitted by a simple power law   Table 1. using ( ) and ( ) as temperature-dependent fit parameters with ∈ ]0,1[ . The obtained values are listed in Table 1.  Figure 9 shows a series of SEM images obtained for different sphere sizes and sintering temperatures after 18 hours of sintering time. Images after shorter sintering times show the same overall pattern but with less pronounced evolutions of the interconnections. Following the indications of Figure 1, we denote in the following   the average radius of the spheres and   the average radius of their interconnections. From Figure 9 we can see very clearly the influence of the sphere radius and the sintering temperature T on the kinetics of formation of the interconnections. Globally, one observes that the radius of the interconnections increases with temperature. For example, at the lowest temperature and for the largest sphere sizes, the interconnections are almost singularities as it is the case for two hard spheres in contact. With increasing temperature or decreasing sphere size, the deformation is accelerated and the spheres start to flatten earlier, forming a neck which can be observed from the side or through the footprint left after breaking the contact (i.e. the flat circles visible on the spheres). This flattening is a consequence of the "crushing flow" (Figure 2a Section 2) as described by Milner. It is clearly seen at 40°C for the largest spheres (  = 125 µm). For smaller spheres at higher temperatures, one observes that the radius of curvature in the cusp of the interconnection increases, i.e. the interconnection does not only grow in area but its boundary also becomes increasingly rounded. This is a signature of the onset of the so-called "smoothing flow" shown Figure 2b and discussed Sections 2.2. For the smallest spheres (  = 33 µm), the spheres coalesce completely at the highest temperature. This is why the item is represented by a cross in Figure 9. Figure 10 shows the macro-porous polyurethanes obtained from the sphere templates shown in Figure 9. One notices that these are accurate negative copies of the initial template, containing spherical pores with circular interconnections. Accordingly, the radii of the interconnections increase in the same manner with increasing temperature and decreasing sphere radius. For the structures with the smallest pore radii, the interstices between the spheres are too small to allow for the polymer to spread homogeneously through the template. This is why these structures are represented by a cross in Figure  10. One notices that when the radius of the interconnections becomes of the order of the sphere radius, the structure loses its spherical morphology and evolves into a very open "bicontinous like" shape, which may be interesting for applications looking to optimise permeability while maintaining a characteristic pore size and a certain mechanical stability of the material. This transition happens between 35 and 38°C for the smallest pore radius (  = 33 µm) and between 40 and 42°C for the structure with the intermediate pore radius (  = 62 µm). For the largest pore diameters, the transition happens above 42°C.

Please do not adjust margins
Please do not adjust margins Figure 11a shows an example of the evolution of the mean radius   of the interconnections with time for three sphere radii (  = 33, 62 and 125 m) at 38°C. The error bars are obtained by the procedure described in Section 3.5. Despite some scatter, all curves follow a clear trend: a rapid increase of   over a few hours is followed by a much slower evolution. Moreover, the evolution of   depends clearly on the radius   of the spheres. Figure 11b shows for the same sphere radii, how the mean radius   of the interconnections after 18 h sintering time depends on the sintering temperature T in a temperature range of 35-42°C. As expected,   increases significantly with the sintering temperature and the sphere radius  . scalings collapse the data fairly well onto a mastercurve within the experimental scatter, implying that for the chosen data both scalings capture reasonably well the main underlying mechanisms. It also confirms that the ratio ( /) evolves with temperature in the same manner for all sphere sizes, implying the absence of finite size effects. The interpretation of the overall shape of the master curves depends on how exactly ( /) depends on temperature, which is discussed in more detail later. For illustrative purposes, Figure 11 contains only a subset of the full data set. The full data set is shown in Figure 12a for different sphere sizes (R  = 33, 62, 125 m), sintering temperatures (T = 35, 38, 40, 42°C) and sintering times (2, 4, 7, 11, 18 hours). To scale all data within the same plot and to allow comparison between the models, we rewrite Equ. (3) to obtain the scaling

4.2) Quantitative analysis and interpretation
where ( ) is a function of the temperature T only. The resulting mastercurve for Frenkel's scaling is shown in Figure  12b (R 2 = 0.48), while Milner's scaling is shown in Figure 12c (R 2 = 0.61). A first observation is that while the data remain fairly Please do not adjust margins Please do not adjust margins scattered after application of Frenkel's scaling, Milner's scaling collapses the data somewhat more convincingly, especially for the lower temperatures. To quantify the quality of the scaling, we fit the scaled data to a second order polynomial with a,b and c being the fitting parameters. This captures in particular the temperature variation of the viscosity  of the paraffin, which is quite important since the experiments are done in the first endothermic shoulder in the DSC curve ( Figure  7) close to the melting transition. The temperature variation of  can be assumed negligible in comparison to that of  in our temperature range. The resulting fits are shown by the red lines in Figure 12b and c, while the resulting fit parameters and coefficient of determination R 2 are given in Table 2. Values of for Frenkel's and Milner's viscous model can therefore be approximated by using Equ. (26) and the associated fit parameters.  Please do not adjust margins Please do not adjust margins It is particularly evident by plotting the experimental results of  /  versus the theoretically predicted in Figure 13 Milner's viscous scaling collapses slightly better the data than Frenkel's scaling. Nevertheless, the spread in the collapsed data remains unsatisfyingly large, indicating the lack an important ingredient in either model. The origin of this scatter becomes readily evident by considering the rheology data of the paraffin shown in Figure 8. While the scalings used until now can take into account the dependence of the viscosity  of the paraffin on temperature, they cannot account for its dependence on the characteristic time scale of solicitation i.e. the strongly viscoelastic behaviour displayed by the paraffin in Figure 8 is not reflected in any of the models used up to this point. We therefore use the visco-elastic model provided by Milner (Equ. (15) in Section 2.3) with the appropriate modification that takes into account the rheological response of our paraffin, given by Equ. (21). The only unknown (temperature-dependent) parameter is the surface tension of the paraffin, while and are known from fitting data of the rheological measurements ( Figure 8). We therefore performed a fit to all data assuming for where ≡ 3 sin ( 2  ) 32 (1 − ) (1 − ) .
Here represents a measure, namely a set of time or temperature (determining the viscoelastic parameters = ( )and = ( ), average radii of the spheres 〈 〉 and the corresponding average radii of the interconnections 〈 〉).  Table 2. temperature dependence of ( ) we found is not physical, since ( ) has to be a decreasing function of temperature. However, a significant increase in temperature is observed only for the highest temperature T= 42 °C. For the lower temperatures, is essentially independent of temperature. At the highest temperature, the problem may presumably arise from a failure of our assumption according to which the viscoelastic behaviour of ( ) = − could be extended to the extremely low frequency range appropriated for our experiment. At this temperature, the longest relaxation time of the paraffin may be shorter than the 18 hours of our longest experiment. Another possible cause of discrepancy could be sought in the fact that for the highest temperature, the sintering dynamics is the fastest, and therefore, the largest values of  /  are measured. As a result, some quantitative discrepancies with the theory which is intended to describe the experimental data is slightly above the theoretically predicted values. This deviation may be explained by three arguments. On the one hand it is coherent with the onset of the smoothing flow, which is expected to start playing a role for / > 0.25 and which follows a different scaling, as discussed in Section 2.2. The smoothing flow would lead to larger 〈 ( )〉/〈 〉-values, as observed. On the other hand, our extrapolation of the power laws describing the rheology data (Section 3.7) may break down in this limit which concerns either high temperatures or long time scales. Besides, it is worth to mention that when neck radii approach the order of magnitude of sphere radii, cross interactions of close necks may not be negligible anymore and thus, account for the gap between experimental and theoretical values (see Figure 10). And last but not least, at these large 〈 ( )〉/〈 〉 values, interconnections start to interact on the surfaces of the spheres, leading to a shift away from the scalings which are derived for isolated interconnections. In summary we can state that both of Milner's scalings capture our experiments better than Frenkel's scaling, while the quality of the visco-elastic scaling clearly indicates the importance of taking into account the Non-Newtonian nature of the paraffin spheres.

Conclusions
In this article we present the results of thoroughly conducted sintering experiments with paraffin spheres. All experimental results are compared to the historic scaling proposed by Frenkel and to the recent models proposed by S. Milner 46,47 . While both models are similar in their approach, we believe that Milner's approximation of the flow field being confined to the contact zone (Section 2.1) is physically better justified. The resulting viscous model (Equ. (8)) leads to slightly different exponents (1/3 instead of 1/2) and pre-factors (3/32 instead of 4) with respect to Frenkel's model, with Milner's model fitting slightly better our experimental data. However, a real breakthrough in the quality of the description of the data is made using Milner's visco-elastic model (Equ. (15)), which takes into account the visco-elastic nature of the paraffin spheres. The general applicability of this model will allow a direct transfer of the predictions to other sphere systems, provided their visco-elastic properties can be quantified. Milner's visco-elastic model captures very well our data up to deformations of 〈 ( )〉/〈 〉  0.25. Beyond this point, the onset of the smoothing flow (Section 2.2), interactions between contact areas or the breakdown of our rheological approximations may play a role. In order to fully validate Milner's model and its limits, it will be important to conduct future experiments with different types of model-spheres in which both, visco-elasticity and surface energy are known and controlled in a quantitative manner. It would be important to validate not only the Non-Newtonian model, but also to access the transition to the regime in which the smoothing flow becomes dominant. Last but not least, it would be interesting to analyse the influence of the multiple contacts on each sphere 54 . Beyond the analysis of the sintering process, we successfully used the obtained sphere templates to generate macro-porous polyurethanes whose morphology can be finely tuned in a predictive manner thanks to Milner's model. Since paraffin and polyurethane can be replaced by a wide range of alternative formulations, this approach can be transferred easily to a multitude of systems which require a fine control over the pore morphology. Indeed, by the knowledge of only a few parameters: surface tension, sphere size and the temperaturedependent rheological data ( and ), the sintering behaviour of the system can be predicted