{'fraction_non_alphanumeric': 0.04240021020949287, 'fraction_numerical': 0.01832166829897523, 'mean_word_length': 4.651187904967602, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Boltzmann Simulations of Droplet formation in confined Channels with Thermocapillary flows

Based on mesoscale lattice Boltzmann simulations with the"Shan-Chen"model, we explore the influence of thermocapillarity on the break-up properties of fluid threads in a microfluidic T-junction, where a dispersed phase is injected perpendicularly into a main channel containing a continuous phase, and the latter induces periodic break-up of droplets due to the cross-flowing. Temperature effects are investigated by switching on/off both positive/negative temperature gradients along the main channel direction, thus promoting a different thread dynamics with anticipated/delayed break-up. Numerical simulations are performed at changing the flow-rates of both the continuous and dispersed phases, as well as the relative importance of viscous forces, surface tension forces and thermocapillary stresses. The range of parameters is broad enough to characterize the effects of thermocapillarity on different mechanisms of break-up in the confined T-junction, including the so-called"squeezing"and"dripping"regimes, previously identified in the literature. Some simple scaling arguments are proposed to rationalize the observed behaviour, and to provide quantitative guidelines on how to predict the droplet size after break-up.


I. INTRODUCTION
Droplet control and production in microfluidic devices has gained a considerable importance, due to the ample variety of applications involved. Droplets can indeed transport very small volumes of fluids, where reagents can be stored and analyzed for a variety of microfluidic applications [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. In this context, droplet manipulation based on thermocapillary stresses has also been established [13,14,[16][17][18][19][20][21][22][23]. Thermocapillary stresses arise because of a surface tension gradient induced by temperature differences, with the surface tension at the interface between two fluids being a decreasing function of the temperature in most of the practical situations encountered [13,14,16,22,24]. The resulting "Marangoni flows" induce droplet motion from the regions of high surface tension to the regions of low surface tension. The focus of this paper is on mesoscopic numerical simulations to quantify the effects of thermocapillary forces on geometry mediated droplet production in confined channels. For the confined channel, we use a T-junction geometry [25,26], where a dispersed (d) phase is injected from a side channel perpendicularly into a main channel containing a continuous (c) phase. The (periodic) production of droplets is then induced by cross-flowing and is affected by the Capillary number Ca -which quantifies the importance of the viscous forces with respect to the surface tension forces -and by the flow-rate ratio φ = Q d /Q c of the two immiscible fluids [11,25,[27][28][29][30]. Distinct mechanisms of droplet formation have been identified [11,25]: these include both the so-called "squeezing" and "dripping" regimes [11,25,26,31]. Specifically, in the squeezing mechanism the break-up process is determined by the build-up of pressure in the main channel caused by the obstruction of the dispersed phase entering from the side channel [11,25]. In the dripping regime, instead, shear forces start do be relevant and promote a smaller droplet size if compared to the squeezing regime [11,25]. The Newtonian problem has been very well characterized in a series of papers in the literature, both experimentally and numerically (see [11, 25-27, 31, 32] and references therein). In this paper, we aim to take a step forward, by exploring the influence of thermocapillarity on the dynamics and break-up properties of fluid threads in the confined T-junction geometry. Our numerical approach is based on lattice Boltzmann (LB) models [33][34][35][36], which proved to be very flexible numerical tools for the simulation of a wide variety of problems involving droplets dynamics in microfluidic geometries [28-30, 32, 37-40, 40-44]. LB allows to solve the diffuseinterface hydrodynamic equations of a binary mixture of two components [45][46][47][48][49][50] in presence of phase segregation, which confines "pure" fluids to different regions, separated by a non ideal interface with a positive surface tension. Multiphase flows with thermocapillary effects have already been investigated in the context of the LB methodology. In particular, Liu et al. [17-19, 23, 51, 52] developed LB multiphase models to account for thermocapillary flows. The LB model proposed by Liu et al. [17][18][19]23] relies on an "improved color-fluid" LB scheme [53], in which the non-ideal forces and the thermocapillary stresses are introduced using the idea of "continuum surface force" [17-19, 23, 51]. The authors have successfully described droplets dynamics and migration, and the effects of localized heating sources in confined environments. Our study differs from the works by Liu et al. [17-19, 23, 51, 52] both in terms of the LB methodology and in terms of the physical problem considered. From the methodological side, the corresponding modelling for the thermocapillary flows is developed on the top of the "Shan-Chen" model for non-ideal binary mixtures [54][55][56][57][58]. This may be useful, for example, in those situations where structured fluids and emulsion droplets [57,59] need to be simulated. The paper is organized as follows. In section II we will present the necessary mathematical background for the problem studied, showing the equations that we consider in both the continuous and dispersed phases. In subsections II A-II C we will recall the basic details of the LB scheme that we use to solve the fluid dynamical equations. In subsection II B we will illustrate the construction of the thermocapillary coupling on the top of the "Shan-Chen" model. In section III we will apply the numerical methodology to the confined T-junction geometry, and we will provide quantitative details on how the thermocapillary effects influence the transition from the squeezing dominated regime to the dripping dominated regime. We will finally rationalize the observed behaviour and provide quantitative guidelines on how to predict the droplet size after break-up. Conclusions will follow in section V. Useful benchmark tests and validation studies will be discussed in Appendix A. Some technical details on the problem of the motion of a cylindrical droplet subject to a uniform temperature gradient are reported in Appendix B.

II. CONTINUUM EQUATIONS AND NUMERICAL METHOD
The equations we consider in the continuous and dispersed phase are the Navier-Stokes (NS) equations in the form (subscripts c/d refer to the continuous/dispersed phase) Here, u and η are the velocity and the dynamic viscosity, respectively. ρ is the total density, P b the bulk pressure, and (∇u) T the transpose of (∇u). Immiscibility between the dispersed phase and the continuous phase is introduced using the so-called "Shan-Chen" model [54][55][56]60], which ensures phase separation with the formation of stable (diffuse) interfaces between the two phases with a positive surface tension σ. An additional interface force is introduced due to the temperature-dependent surface tension σ = σ(T ). The associated interfacial force is most conveniently written in the limit of a sharp interface, using the differential operator ∇ S = (δ −nn) · ∇, providing projection tangentially to the interface [17][18][19]23] where δ S is the Dirac delta function and C is the (local) curvature at the interface [17][18][19]23]. Equivalent expressions of (2) for the general case of diffuse interfaces are provided in subsection II B. In presence of thermocapillarity, we need to relate the surface tension to the temperature. For the sake of simplicity, we only consider a linear relation between the surface tension and the temperature σ(T ) = σ 0 + σ T (T − T 0 ), where T 0 is a reference temperature, σ 0 is the surface tension at T 0 , σ T = dσ/dT is the rate of change of surface tension with temperature and is an input parameter in our numerical simulations. Together with the macroscopic NS equations (1), we have considered the following advection-diffusion equation (i.e. passive-scalar equation) for the temperature field in both the continuous and the dispersed phase [61] where κ (T ) is the thermal diffusivity. Equations (1)-(3) are solved by means of an algorithm combining a multicomponent LB model coupled to an LB solver for the advection-diffusion equation (3). Technical details on the multicomponent LB and the advection-diffusion LB scheme are not new, hence they are just briefly recalled in subsections II A and II C. In subsection II B, instead, we will discuss more extensively on the modelling of thermocapillary effects, since it constitutes a novelty from the methodological point of view for the paper.

A. Modelling of Multicomponent Fluid
The multicomponent LB takes into account the evolution dynamics of the discretized probability density function f s (x, t) to find at position x and time t a fluid particle of component = A, B with velocity c s . The numerical simulations described in the paper make use of the D3Q19 model [81] whose discrete velocity set is given by A majority of one of the two components then identifies the dispersed and the continuous phases in the hydrodynamic equations (1). The LB dynamics accounts (over a unitary time step) for both the relaxation and the streaming of f s (x, t) [34]: where the matrix terms A sj (the same for both components) drive the relaxation towards the local equilibrium [34] f (eq) s where c S = 1/ √ 3 is the isothermal speed of sound, u is the fluid velocity, and the weights w s are given by w s = 1/3 for s = 0, w s = 1/18 for s = 1 . . . 6 and w s = 1/36 for s = 7 . . . 18. We work in the framework of the multiple relaxation time (MRT) models [62][63][64][65]. This offers the particular advantage to control the various processes (mass diffusion, momentum diffusion, etc) independently. The MRT models make use of a "mode space" identified with a suitable basis vectors b k (k = 0, ..., 18), hence the "modes" calculated as M k = s b ks f s (k = 0, ..., 18) [62][63][64][65]. The collisional term then describes a linear relaxation of the modes, M k → (1 + a k )M k + M g k , where the relaxation frequencies −a k are related to the transport coefficients in the hydrodynamic equations (see also below). The lowest order modes are directly linked to the hydrodynamic variables as follows. The density of both components (ρ ) and the total density (ρ) coincide with the zero-th order modes The higher order modes refer to the viscous stress tensor, and also other modes which do not show-up in the hydrodynamic limit [62][63][64][65]. The term S g s is a forcing source embedding the effects of a forcing term with density g [62,64,65]. The term g = g in Eq. (7) is the total force. Given the relaxation frequencies of the momentum (−a M ), bulk (−a b ) and shear (−a s ) modes, the forcing source that ensures the correct recovery of the hydrodynamic equations can be computed exactly [65]. The hydrodynamic evolution is obtained by considering the long-wavelength limit [34,66] of the equations for the coarse-grained density and momentum. One can indeed show that the LB scheme reproduces the continuity equations and the NS equations [62,64,65]: The viscosity coefficients are determined by the relaxation frequencies of the bulk and shear modes in (5)  In equations (9), p = p = c 2 S ρ represents the internal (ideal) pressure of the mixture, while the quantity D is the diffusion flux of one component into the other (see [65] for detailed expressions). Regarding the internal forces, the "Shan-Chen" interaction model [54][55][56]65] for multicomponent LB schemes is used, where the force experienced by the particles of the -th species at x, is due to the particles of the other species at the neighbouring locations compatible with the lattice links where g AB is a parameter that regulates the intensity of phase segregating interactions. For simplicity, the sum in equation (10) extends over a set of interaction links c s coinciding with those of the LB dynamic (5). For the purposes of this paper, the effect of the internal forces is more conveniently analyzed with the help of the interaction pressure tensor P (int) [58], whose gradient is minus the internal forces, i.e., g = −∇ · P (int) . The interaction pressure tensor modifies the total pressure tensor of the model, i.e., An explicit expression for P (int) can be derived directly on the lattice [56,58]. For the case of multicomponent fluids, it reads [58] Upon Taylor expanding the expression (12), we get (we drop the x-dependency for simplicity) where we recognize a bulk pressure contribution, P b = p+c 2 S g AB ρ A ρ B , and other contributions which are proportional to the derivatives of both densities. These gradient terms establish a diffuse interface with a positive surface tension (see also section II B) whenever phase separation occurs in the model [57]. Wetting properties are introduced in the model by imposing specific values of the density in contact with the wall as described in [67,68]. The numerical simulations presented are carried out with g AB = 1.5 lbu (lattice Boltzmann units, henceforth) in (10), corresponding to a surface tension σ = 0.08 lbu and associated bulk densities ρ A = 2.1 lbu and ρ B = 0.12 lbu in the A-rich phase. The relaxation frequencies are set to a s = a b = −1.0 lbu, thus reproducing the viscous stress tensor given in Eqs. (1) with η c,d = 0.28 lbu.

B. Modelling of Thermocapillary Effects
In order to model the thermocapillary effects with a temperature-dependent surface tension, suitable terms need to be introduced on the top the bare pressure tensor of the model (12). To this aim, we focus the attention on the mechanical balance equations at a curved two-dimensional (for simplicity) interface. In what follows, the local radius of curvature will be denoted with R. First of all, we notice that there exists a "gauge invariance" in the definition of the pressure tensor [69], in that we can define P (int) modulo a tensor with zero divergence. Specifically, if we introduce the pressure tensor one can easily verify that (14) has zero divergence, meaning that the pressure tensor brings the same physics of the pressure tensor P (int) . For our purposes, we then use the pressure tensorP (TOT) = P (TOT) − P (0) . Some straightforward algebraic calculations reveal the following relation betweenP (TOT) and the normaln at the non-ideal interfacẽ where we have used ∇ρ A · ∇ρ B = −|∇ρ A ||∇ρ B | and also the property of the normaln which is taken to point from the A-rich into the B-rich phase. We then introduce the tensor with the property to project along the tangential directiont, perpendicularly ton, i.e. T ·n = 0. For our curved interface, we can then compute the pressure tensor in the directions normal (N) and tangential (T) to the interface and obtaiñ so thatP and the pressure tensor can be rewritten as [70] When analyzing the equilibrium properties (∇ ·P (TOT) = 0) with this decomposition, the mechanical balance at the interface is better understood. Indeed we get where we have introduced the gradients along the tangential and normal direction, i.e. ∇ n =n · ∇, ∇ t =t · ∇. Further, we have used that ∇ · T = − 1 Rn and also assumed that the gradient of the density fields is perpendicular to the interface (i.e. directed alongn). When we deal with cases with homogeneous and isotropic surface tension, Eq. (19) leads to the well known Laplace equation. To see this point, let us specialize our equations to the case of a cylindrical droplet described in polar coordinates, with the center of coordinates coinciding with the center of the droplet (r = 0). Because of the cylindrical symmetry, we can neglect the termt ∇ t (P N −P T ) to obtain Integrating in the normal (i.e. radial) direction, we get that is the usual Laplace law written for the bulk pressure P b evaluated in the outer (out, r = ∞) and inner (in, r = 0) regions. The surface tension σ follows In order to construct a force driven by gradients in the surface tension, we introduce "on the top" of the total pressure tensorP (TOT) a temperature-dependent term. In particular, following [60], and relying on the MRT schemes, we have added the new term to the stress modes of the multicomponent model, in such a way that the resulting hydrodynamic Eqs (1) with the coupling (2) are recovered in the limit of a thin interface. The choice for the temperature-dependent term is quite phenomenological [17][18][19]23]: it is chosen to be proportional to the squared gradient of the "order parameter", the latter taken as the density difference ρ d = ρ A − ρ B , and is also proportional to a temperaturedependent prefactorP From (16)- (17), it follows that the difference between the normal and tangential component of the pressure tensor changes as ( Hence, the balance equation at a curved interface is modified as that is producing a temperature dependent surface tension is tangential at the interface and is responsible for the Marangoni effects triggered by the tangential variation of the surface tension σ(T ) [17][18][19]23]. The simplest choice for A(T ) reproducing the motion of a curved interface from cold to hot regions is [17][18][19]23] where G 0 is a constant and G T > 0. With this choice, the surface tension is larger in regions where the temperature is lower.

C. Modelling of Advection-Diffusion Equation for the Temperature Field
For the temperature dynamics, following [61], we add another set of populations g s (x, t). At difference with the multicomponent case, instead of using the MRT approach, the collisional part of the LB dynamics is now taken to be a BGK relaxation dynamics for simplicity where τ g is a characteristic time governing the relaxation dynamics towards the local equilibrium distribution As in the case of the multicomponent part of the model, the hydrodynamic evolution for the temperature results from the long-wavelength limit [34,66] of the equations for the coarse-grained temperature fields, the latter computed as the "density" of the g s population, i.e. T = s g s . One can show that Eqs (3) are recovered with κ (T ) c,d = c 2 S (τ g − 1/2). Let us also note that, in order to get the right hydrodynamic limit of the temperature evolution, an extra body-force term is in principle needed in order to avoid spurious terms in the continuum limit [71]. The importance of this term may depend on the applications. In all our simulations, we have checked that it is negligible. Notice that for more complex cases where thermal effects may enter the equation of state, one can also devise different strategies [72,73].

III. T-JUNCTION SET-UP
In this section we report the main details for the T-junction set-up used for the numerical simulations. The Tjunction is embedded in a rectangular parallelepiped with size L x × L y × H, and both the main and side channels have a square cross-section resolved with H × H = 32 × 32 lattice cells. The resolution used for H is chosen to be compatible with those used for the successful benchmark tests described in Appendix A. The main channel and the side channel lengths are resolved with a L x = 640 lattice cells and L y = 160 lattice cells, respectively. All the simulation details are summarized in table I. In order to highlight the thermocapillary effects in the transition from squeezing to dripping, at changing a variety of flow-rates in both the continuous and dispersed phases, we decided to keep the viscous ratio λ = η d /η c and thermal diffusivity ratio χ = κ (T ) c /κ (T ) d fixed to unity. Moreover, the Reynolds number, Re = ρv c H/η c , is small, hence it is expected a minor influence on the thread deformation and break-up [74]. We then choose the following dimensionless numbers [25,27,31] to fully characterize the system. First, the Capillary number Ca that is calculated in the continuous phase and the flow-rate ratio where Q d = v d H 2 and Q c = v c H 2 are the flow-rates at the two inlets of the side and main channel. We emphasize that the flow-rates are input parameters for the T-Junction set-up. In addition to the two control parameters Ca and φ, to quantify the effects of thermocapillarity, we use the Marangoni number where ∆σ H represents the variation in the surface tension over a distance H (see also Appendix B). In the outlet, we impose a fixed pressure and also use Neumann-type boundary conditions for the velocity fields. A Dirichlet-type boundary condition is set at the inlets of the main and side channel, by imposing the pressure gradient computed from the exact solution of a Stokes flow in a square duct [38]. As for the temperature boundary conditions, we set a constant temperature at the inlet (T = T inlet ) and outlet (T = T outlet ) of the main channel. In contact with the walls of the main and side channels, we set a Dirichlet-type boundary condition for the temperature field that is a space-dependent one: it coincides with the linear profile going from T inlet to T outlet . In the inlet of the side channel we set the temperature to be consistent with such linear profile. By suitably changing the parameter G T > 0 in (22), while setting T inlet − T outlet = ±1.4 lbu, we are therefore able to explore situations where the temperature gradient promotes migration in the positive/negative main channel direction, dependently on the sign of the temperature gradient ∇T = (T inlet − T outlet )/L x (see also Fig. 1). Notice that we use the same temperature gradient that we tested in the benchmark tests in Appendix A, hence we already know the characteristic thermocapillary velocities developed in such temperature gradients (see Fig. 8). Moreover, in presence of a non zero temperature gradient, the surface tension σ = σ(T ) used to compute the Capillary number (25) is the one corresponding to the temperature at the T-junction [82], i.e., where the break-up process takes place. In addition, wetting properties need to be introduced at the boundaries. To this aim, we impose perfect wetting for the continuous phase.

IV. RESULTS FOR T-JUNCTION GEOMETRY
We start our numerical investigation by getting a qualitative idea of the thermocapillary effects in the confined T-junction geometry. Specifically, in Fig. 1 we present a time series of 3D snapshots of density contours overlaid on the viscous stress magnitude at fixed Ca for different scenarios, depending on the sign of the temperature gradient. The Marangoni number Ma = 1.26 is kept fixed. Moreover, we have chosen a relatively small Ca = 0.001, which falls into the squeezing regime [25,27,28,31], where the incoming thread tends to occupy and obstruct the cross-section of the main channel before breaking at the junction. In this regime, the interfacial forces are dominant, and therefore we expect the thermocapillary effects to be more pronounced. Cases with zero temperature gradient -hence no thermocapillary effects -are reported in the left column (Panels (a)-(e)), and they are useful to compare with the cases where the temperature gradient is different from zero. In the middle column we report cases with a positive temperature gradient and Ma = 1.26 (Panels (f)-(j)), while the right panel reports cases with negative temperature gradient and Ma = 1.26 (Panels (k)-(n)). Notice that on the top row we have reported the temperature profile in a slice taken at z = H/2 for the three cases analyzed. Without thermocapillary effect, the thread starts to enter the main channel (Panel (b)), obstructs the main channel (Panel (c)), the continuous phase squeezes the dispersed phase (Panel (d)), and finally the droplet generation takes place at the junction (Panel (e)). The droplet formation process is visibly influenced by the switch-on of a non-zero temperature gradient: when Ma = 1.26, ∇T > 0 we observe that the droplet size is smaller than the corresponding case with ∇T = 0; additionally, the break-up process is faster, leading to the production of two droplets within the same time it takes to generate one droplet with ∇T = 0. On the contrary, when Ma = 1.26, ∇T < 0, the break-up process is delayed, and no droplet formation is observed within the time lapse of a few shear times τ shear = H/v c , which is taken as the characteristic unit of time for the break-up process. The case with negative temperature gradient is also visibly characterized by a reduction of the viscous stress at the junction, which may be taken as a further indication that the whole process is evidently affected by thermocapillary stresses.
To go deeper and characterize the effects of thermocapillarity on the droplet size at small Ca, we have further performed different numerical simulations at changing the flow-rate ratio, the Marangoni number, and the sign of the temperature gradient. This is necessary to compare our results for the droplet size with the scaling laws which have been proposed in the literature [25,31]. As already highlighted in Fig. 1, in the regime of small Ca, the dispersed phase enters the main channel until it obstructs the flow. Afterwards, due to the squeezing of the continuous phase, the dispersed thread breaks at the T-junction [25,27]. Assuming that the thread squeezes at a rate proportional to Q c , a scaling law for the droplet volume has been established in the literature [11,25] V For a fixed Capillary number, the constants α 1 and α 2 (which are of the order one) are determined by the geometrical details of the channels [25,31]. As already put forward by other authors [11,25], the relation (28) allows to match together two scaling laws expected in the regimes where Q d is smaller or larger than Q c . When Q d is larger than Q c , one expects that the time to squeeze the droplet is larger than the time it takes to block the channel, which implies the droplet size to scale linearly with the flow-rate ratio. On the other hand, when Q d is much smaller than Q c , one expects the dimensionless droplet size to be independent of the flow-rate ratio [11,25]. Equation (28) has already been verified in experiments [3,7,11] and also in numerical simulations [25,[27][28][29][30][31]. It is therefore important to check both the prediction of a linear scaling in φ, as well as a possible dependency of the scaling constants α 1 and α 2 on the Marangoni number and the sign of the temperature gradient. To this aim, in the left Panel of Fig. 2, we plot the dimensionless droplet volume versus the flow-rate ratio at Ca = 0.0014, for different values of Ma and different signs of the temperature gradient. For all the cases investigated, we confirm the validity of the linear scaling relation (28), with the curves that are steeper for negative temperature gradients. Notice that the effect of thermocapillarity is larger (smaller) at the larger (smaller) φ, which may be taken as a clear indication that one needs a sufficiently large squeezing time in order for the Marangoni stresses to show-up with a quantitative effect. To make further progress, from the data reported in the left Panel of Fig. 2, we have obtained the constants α 1 and α 2 from the linear fit based on Eq. (28). These two coefficients are displayed against the Marangoni number in the right Panel of Fig. 2. Notice that, in order to visualize together the effects of both positive and negative temperature gradients, we use The effect is actually more visible on the slope coefficient α 2 . The plot of log α 2 indeed shows two distinguished slopes, dependently on the sign of the temperature gradient. This change in slope is somehow indicative that the geometry mediated break-up process is affected in a non trivial way by the thermocapillary effects: for ∇T < 0 the thermocapillary effects oppose to the squeezing process, producing larger droplet sizes and somehow "pushing" the system even more towards the squeezing regime, favoring the inhibition of the main channel by the dispersed phase (see also right column in Fig. 1). When ∇T > 0, instead, the thermocapillary stresses promote droplet break-up and somehow spoil distinctive features of the squeezing regime, as also evident from the flatter curves displayed in the left Panel of Fig. 2. The latter observation stimulates another interesting analysis, concerning the importance of the thermocapillary effects in the transition from the squeezing regime to the dripping regime [25,[27][28][29][30][31], where the increase of the Capillary number induces a reduced droplet size. For the specific case of our numerical simulations, this analysis is displayed in Fig. 3, where we monitor the droplet size by changing the Capillary number for fixed flow-rate ratio. We considered two cases with flow-rate ratio φ = 1.0 (left Panel) and φ = 0.5 (right Panel). Furthermore, the effects of thermocapillarity are also shown, in that we have conducted simulations by switching-on a positive (blue triangles) and a negative (black squares) temperature gradient for the same Marangoni number Ma = 1.26. For a given Ma and ∇T , the droplet size decreases as a function of Ca, due to the increased importance of the viscous  stresses over the surface tension forces. We actually see that for the largest Ca analyzed the droplet sizes come out to be comparable for all the cases, while at smaller Ca data with a positive (negative) temperature gradient underestimate (overestimate) the "bare" data, i.e. those with zero temperature gradient. These facts said, we have tested the possibility to design a quantitative guideline to predict the droplet size, once the bare Newtonian data and the input thermocapillary parameters are known. Specifically, we checked whether a suitable redefinition of the Capillary number in terms of the input thermocapillary parameters allows to rescale the results with positive/negative temperature gradients on the bare Newtonian data. Notice that, once the temperature gradient and the Marangoni number are assigned, a characteristic Marangoni velocity v m naturally follows, which scales linearly with Ma. With such velocity at hands, one can define a rescaled Capillary number as Equation (29) takes into account the balance between the viscous forces and surface tension forces, including also a modification coming from the Marangoni velocity v m , whose sign is positive/negative for positive/negative temperature gradients. We have then considered the data in the main Panels of Fig. 3 and tested this rescaling. Results are reported in the insets of Fig. 3. Remarkably, by using a unique value of the Marangoni velocity v m = ±0.00035 lbu (for Ma = 1.26), we are able to rescale the thermocapillary data on the top of the bare Newtonian data. This holds true for both the flow-rate ratios φ = 1.0 and φ = 0.5. We emphasize that we have not predicted the value of the thermocapillary velocity v m , but just found it at a fixed Ma on the basis of the best collapse of the data. Nevertheless, we notice that the value falls well in the ballpark of the values of steady migration velocities that we have measured and reported in the benchmark tests described in Appendix A. We further checked if the rescaling (29) performs well at changing the Marangoni number. In Fig. 4 we report the droplet size as a , hence no extra fitting parameters have been introduced. As we can see from the figure, all the data nicely collapse on a single master curve, which is a clear indication that the proposed argument is robust at changing the Marangoni number, at least up to Ma ≈ 1. Finally, some points of discussion emerge, to connect the results of the present numerical simulations to experimental data, especially with the aim of proposing a set of experimental suggestions. Experimental literature indeed exists (see [13][14][15][16] and references therein), reporting on the motion of droplets and threads in confined channels with thermocapillary forces. For example, the authors in [15] report on the experimental study of temperature-induced migration of a bubble in a Hele-Shaw cell. In the experiments, the bubble velocity hits a maximum value 800µ m/s, and it depends on the temperature gradient ( power supply of 65 − 300 mW), bubble radius (20 − 400µm) and the wall thickness (22 and 37µm). Thus, in an experiment involving the motion of threads in microfluidic T-junctions with characteristic flow-rates of the order of Q c,d ≈ 0.1µL/min [11,75] and channel lengths of the order of a few tens of micrometers (say H ≈ W ≈ 100µm), one expects velocities of the order of v c,d 100µm/s, thus the thermocapillary effects are expected to show-up in sizeable effects.  (29). Notice that the value of the thermocapillary velocity vm is changing with Ma, and is determined by simply rescaling the value of vm found for Ma = 1.26 in Fig. 3.

V. CONCLUSIONS
One of the main features that drives the interest in developing and exploiting microfluidic devices is the possibility to control the process of droplet formation as a function of the flow and material properties [1][2][3][4][5]. The complex flows in these devices are naturally accompanied by confinement, which has an important effect on the droplet break-up process. In this context, droplet manipulation based on thermocapillary stresses has been demonstrated both experimentally and numerically (see [13,14,[16][17][18][19]23] and references therein). Thermocapillary stresses arise along an interface because of a temperature dependency of the surface tension. Based on lattice Boltzmann simulations with the "Shan-Chen" model for non-ideal binary mixtures [54][55][56][57][58], we have investigated the combined effects of thermocapillarity and confinement on the dynamics and break-up of fluid threads in a microfluidic T-junction geometry, where a dispersed phase is injected perpendicularly into a main channel containing a continuous flowing phase which induces periodic formation of droplets. The numerical model for multicomponent fluids [55,58] has been first extended to include the effects of thermocapillarity, and then has been applied to the simulations in the confined T-junction geometry, by focusing the attention on the operational regimes of "squeezing" and "dripping" [25] with moderate flow-rate ratios. We have changed various parameters, including the temperature gradient ∇T , the Capillary number Ca, and the Marangoni number Ma, to get a quantitative understanding on how the thread dynamics and subsequent break-up are affected by the thermocapillary forces. To this aim, we wish to stress the crucial role played by numerical simulations, in allowing a systematic exploration of the importance of the relevant free parameters in the model. With the data at Ma = 0 at hand (the "bare" data), we have indeed tested the possibility to design guidelines to predict the droplet size for the cases where thermocapillary effects are present (Ma = 0). Numerical results indicate that a suitable redefinition of the Capillary number based on a characteristic thermocapillary velocity well accounts for the observed changes and allows to well rescale the thermocapillary data on the top of the bare data. Remarkably, this simple argument is robust at changing both the flow-rate ratio and the Marangoni number, up to Ma ≈ 1. Quantifying such rescaling properties at higher Ma sets interesting challenges for the future: for the description of the thermocapillary migration of a droplet at very large Ma, previous studies suggest non trivial dependencies on Ma (see [17] and references therein); our model is able to simulate flows at those Marangoni numbers (see Appendix A), and therefore it is definitively warranted a systematic investigation, to verify if the aforementioned rescaling holds and if the rescaling velocity is still linearly dependent on Ma or not [17]. Along these lines, it is also definitively warranted a systematic exploration with fluids with dissimilar properties (viscosity ratio and/or thermal conductivities) as well as comparisons with other existing LB methodologies for thermocapillary flows [17][18][19]23]. We actually expect our numerical work can stimulate further investigation and can also bring forth experimental suggestions for droplet control in microfluidic geometries. To perform benchmark numerical simulations for the thermocapillary migration, we choose a two dimensional (2D) set-up. Results are available from the literature [76][77][78][79][80], concerning the motion of a droplet with radius R in steady flows driven by thermocapillary effects in three dimensions (3D). As for the 2D case, the corresponding solution is straightforward to obtain, and details are reported in the Appendix B. The choice of a 2D set-up was initially done to easily approach data at larger resolution, to appreciate the convergence of the numerical results to the continuum limit predictions. Specifically, one can solve the problem of thermocapillary motion of a cylindrical droplet subject to a uniform temperature gradient in the absence of gravity. Assuming negligible advective effects on the momentum and heat transports, one can determine the constant speed v s of translation in the direction of the temperature gradient In the above equation, ∆σ 2R = ∆σ 2R (κ comes from the fact that the temperature at the interface is fixed by the ratio of the thermal conductivities (see Appendix B for details). Notice that for the case of a droplet, the Marangoni number is defined as We performed numerical simulations in a 2D box with size L x × L y = 200 × 640 lattice cells, with two walls located at y = 0 and y = L y . Periodic boundary conditions are applied in the x direction. We performed different numerical tests, at varying the droplet size and the temperature gradient. To start with, in figure 5 we report the evolution in time of the center of mass of a droplet (relative to its initial value) steadily translating in a constant temperature gradient. The droplet radius is resolved with a variable number of lattice cells, and the temperature gradient is consequently adjusted to match the same Marangoni number, Ma = 1.25. As we can see, above R ≈ 15 lattice cells, results show convergence. These converged results are then compared against their theoetical predictions. Specifically, from now on, we report details on a case with droplet radius R = 16 lattice cells placed in the middle of the box, while the temperature of the upper and lower walls are fixed to be T up = −0.7 lbu and T down = 0.7 lbu, respectively. We first consider a case with equal thermal conductivities inside and outside the droplet, κ .166 lbu and same viscosity η c = η d = 0.28 lbu. The tipical velocity vector profiles in the reference frame of the moving droplet are shown in Fig. 6. The velocity vectors reveal recirculations, and the flow pattern is also compared with the analytical solution (right Panel) described in Appendix B. Although not used in the numerical simulations with the T-Junctions, the model is capable of handling also higher Marangoni numbers. Some velocity vectors profiles are reported in Fig. 7. Further quantitative check is offered in Fig. 8, where we analyze the steady velocity v s as a function of the parameter G T in (22). Specifically, the constant G T in (22) has been changed from 0.0 to 1.0. The inset of the left Panel reports the steady velocity (in lbu) that is developed due to the thermocapillary migration at changing G T . Notice that for the set of parameters chosen, the Marangoni number -hence and the steady velocity -are well fitted with a linear behaviour, Ma ≈ 2.52 G T (v s ≈ 0.0008 G T ). The right Panel reports a direct check of Eq. (A1): the black line stands for the rhs of Eq. (A1), with the jump in surface tension ∆σ 2R estimated from the modified pressure tensor described in Section II A. The points represent the steady velocity of migration obtained from the numerical simulations. To complete our benchmarking tests, the analysis reported in the right Panel of Fig. 8 is extended to fluids with dissimilar thermal conductivities. Results are reported in Fig. 9. Overall, Figs. 8.-9 reveal that the agreement between the measured velocity and the expected one is satisfactory, although the measured velocity tends to slightly underestimate the rhs of Eq. (A1) for the largest Ma. This echoes the results of Liu et al. [17], in that the thermocapillary velocity normalized by the corresponding hydrodynamic prediction (valid at small Ma) shows a decreasing behaviour in terms of Ma, especially pronounced for Ma above a few units. To conclude, we remark that in all the tests that we run, we noticed that a characteristic droplet diameter of the order of a few tens of grid points is enough to reproduce well the theoretical prediction of sharp interface hydrodynamics. Since 3D numerical simulations are affordable with such resolutions, this fact is instrumental to up-scale the problem and perform full 3D simulations in a confined T-junction in presence of thermocapillarity, as described in section III.  the temperature gradient. The droplet is centered at the origin, and the temperature profile far away from the droplet is set by T ∞ = T 0 + x(∆T )/R. We also assume a negligible deformation for the droplet and consider an unperturbed cylindrical shape. We then use a cylindrical coordinates system (r, θ) taken from the center of the droplet. The equations of continuity are expressed as where P ∞ denotes the pressure far from the drop. At the drop surface r = R, we need to impose the condition of zero normal velocity and the continuity of the tangential velocity component and the tangential stress u r,c = u r,d = 0,