The Relative Contribution of Solutal Marangoni Convection to Thermal Marangoni Flow Instabilities in a Liquid Bridge of Smaller Aspect Ratios under Zero Gravity

The effect of solutal Marangoni convection on flow instabilities in the presence of thermal Marangoni convection in a Si-Ge liquid bridge with different aspect ratios As has been investigated by three-dimensional (3D) numerical simulations under zero gravity. We consider a half-zone model of a liquid bridge between a cold (top plane) and a hot (bottom plane) disks. The highest Si concentration is on the top of the liquid bridge. The aspect ratio (As) drastically affects the critical Marangoni numbers: the critical solutal Marangoni number (under small thermal Marangoni numbers (MaT As . 1800)) has the same dependence on As as the critical thermal Marangoni number (under small solutal Marangoni numbers (400 . MaC As . 800)), i.e., it decreases with increasing As. The azimuthal wavenumber of the traveling wave mode increases as decreasing As, i.e., larger azimuthal wavenumbers (m = 6, 7, 11, 12, and 13) appear for As = 0.25, and only m = 2 appears when As is one and larger. The oscillatory modes of the hydro waves have been extracted as the spatiotemporal structures by using dynamic mode decomposition (DMD). The present study suggests a proper parameter region of quiescent steady flow suitable for crystal growth for smaller aspect ratios of the liquid bridge.


Introduction
Semiconductor single crystals with high purity are a vital requirement for novel optoelectronic devices in various fields of information technology. At present, the floating zone (FZ) method is one of the most promising techniques to obtain silicon germanium (SiGe) single crystals with superior purity. During this process, a small free melt zone stabilized by surface tension is suspended between the polycrystalline feed rod and the grown single crystal. The crystal growth is achieved by a relatively slow movement of the crystal and feed rod through the molten zone. This container less configuration is beneficial since it eliminates possible crucible contamination. However, the temperature and concentration gradients along the free surface of the liquid bridge give rise to surface tension gradients that can induce strong thermo-solutal Marangoni convective flows. The compositional variations (striations) due to time-dependent Marangoni flows may lead to inhomogeneities in the electronic properties of grown crystals and be responsible for dislocations in grown crystals. A better understanding of the growth conditions under thermo-solutal Marangoni convection is needed to prevent or minimize such undesirable structural flaws.
As mentioned earlier, the FZ technique, being a container-less growth method, has the advantage of growing high purity crystals over other crystal growth techniques. However, it also has challenges. First, the materials used have high melting points and are opac for situ observations. This presents difficulties in experiments. Second, the size of the liquid bridge is limited due to the gravity on Earth. This prevents the growth of larger crystals. The growth of larger crystals is only possible under microgravity conditions. However, the possibility of microgravity experiments is very rare and in addition such experiments are very expensive. Thus, numerical simulation studies under either micro-or zero-gravity conditions present themselves as very viable and inexpensive alternatives to study various aspects of the FZ growth technique [1]. In these simulations, mostly the half-zone model is adopted for computational efficiency by considering the half of the full-zone liquid bridge. The half zone configuration is characterized by a single toroidal vortex developed between hot and cold disks of the zone.
Literature in this area is relatively well developed. For instance, experiments and numerical simulations (see, for instance, in [2][3][4]) performed with high-Prandtl-number fluids considering pure thermal Marangoni convection showed that the flow in the liquid bridge is laminar, steady, and axisymmetric for sufficiently small thermal Marangoni numbers, but when the thermal Marangoni number exceeds a particular threshold, flow show oscillatory three-dimensional (3D) flow patterns. For low Prandtl number fluids (semiconductor melts), the steady flow breaks the axisymmetry and becomes a 3D-steady flow regime (stationary bifurcation) before the onset of a 3D oscillatory flow (Hopf bifurcation) (see, for instance, in [5][6][7]). The numerical simulations [8][9][10] and linear stability analysis [11] on thermo-solutal Marangoni convection in literature have shown that not only thermal Marangoni convection, but also solutal Marangoni convection significantly contributes to the flow structures in the liquid bridge of the FZ method. When considering the combined effects of thermal and solutal Marangoni convections for low Prandtl number fluids, the axisymmetric steady flow bifurcates to 3D oscillatory with the absence of a stationary bifurcation depending on the thermal and solutal Marangoni numbers values. It is noteworthy that a recent study on the characterization of thermal and solutal Marangoni flows of opposite directions in a SiGe liquid bridge by Jin et al. [12] has found the coexistence of both stationary and Hopf bifurcations depending on the thermal and solutal Marangoni numbers.
Numerical results of Lappa et al. [13] have confirmed the existence of different behaviors of Marangoni flows exist depending on the Prandtl number and the geometrical aspect ratio of the liquid bridge. A detailed study of the linear stability analysis for low Prandtl number fluids by Wanschura et al. [7] has found the modes with higher azimuthal wave numbers critical as the aspect ratio decreases. Similar relationships between the aspect ratio of the liquid bridge and wavenumbers for high-Prandtl numbers also been observed experimentally (see, for instance, in [14]). A parametric analysis of the influence of the aspect ratio of the full zone can also be found in the literature [15,16]. The assumption of the full zone model is the only way to capture in detail the mechanism of the instability of Marangoni flow generating in the real floating zone is one of the advantages of full zone model. They have concluded that the critical wavenumber increases if the geometrical aspect ratio decreases but the critical Marangoni number does not change much with the aspect ratio. However, the full zone model has a disadvantage that the temperature at the equatorial cross section cannot be fixed due to convective and heat transfer mechanisms establishing the temperature difference which is driving the surface flow. In the half zone model, on the other hand, the temperature difference can be fixed by holding maximum and minimum temperatures at the end disks. The no-slip condition imposed on the symmetric plane at the mid-height between the rods of the real floating zone makes the flow dynamics of the half zone model much closer to the real floating zone. Thus, the half zone model has become a very important and computationally inexpensive paradigm model for fundamental research.
These studies show that the geometrical aspect ratio of the liquid bridge is as critical as the Marangoni numbers for the thermo-solutal Marangoni flow developing in the liquid bridge. To the best of our knowledge, the dependency of flow regimes on the aspect ratio of the liquid bridge when the thermal and solutal Marangoni convections are coexist is still unknown in the literature. Thus, this was the objective of the present study.
In addition to numerical simulations, the linear stability analysis (LSA), proper orthogonal decomposition (POD), and dynamic mode decomposition (DMD) are some useful techniques to investigate the hydrodynamic processes. DMD is a leading technique in analyzing the spatio-temporal flow patterns in transitional regimes. The POD extracts a hierarchical set of orthogonal spatial modes ranked according to their contribution to the total energy of flow modes. The POD is much suitable, therefore, for extracting the statistical structures from a complex flow. In comparison, DMD assumes physically more relevant modes to the oscillatory growth and decay in time. DMD is an emerging data-driven technique to obtain linear reduced-order models for high-dimensional complex systems. In addition to linear reduced-order models, the spatial-temporal coherent structures or patterns that dominate the observed measurement data from that dynamical system can also be extracted by DMD. This is initially introduced in the fluid dynamics community by Schmid [17], and it has later been connected to this nonlinear dynamical systems theory called Koopman theory by Rowly et al. [18]. This technique comes out of the fluid mechanics community, but since then, DMD has been applied to a broad range of systems, including disease modeling, neuroscience, robotics, and finance. Part of the reason for the broad application of DMD is that it is pure data-driven and does not require any description of underlying governing equations. It can be used as a postprocessing tool for analyzing the spatio-temporal structures of experimental data.
3D flow structures are developing inside the liquid bridge floating zone system due to Marangoni convection. Investigation of the flow dynamics of those structures is a vital requirement for producing defect-free crystals. As the 3D time snapshots by numerical simulations alone do not help to understand the hydrodynamic process, the DMD is giving a proper understanding of the dynamical system.
The present study aims to investigate the influence of aspect ratio on thermo-solutal Marangoni instabilities of the liquid bridge by three-dimensional numerical simulations and also study the spatio-temporal coherent structures directly from numerical data by DMD.

Numerical Simulation
The model of a half-zone liquid bridge between two parallel concentric cold and hot disks of radius a with a distance L apart ( Figure 1) is adopted. The aspect ratio of the liquid bridge is defined as A s = L/a. We assumed (i) the liquid surface is nondeformable, cylindrical, and adiabatic; (ii) the melt (a mixture of Silicon and Germanium) is an incompressible Newtonian fluid; and (iii) the system is under zero gravity, thus the liquid bridge keeps its cylindrical shape as we assume zero gravity.
We consider a cylindrical coordinates system (r, θ, z), and the representative length a, velocity ν/a, temperature difference ∆T, concentration difference ∆C (between top and bottom walls), and time a 2 /ν, where ν(= µ/ρ) denotes the kinematic viscosity given in terms of melt density ρ and melt viscosity µ. Then, we obtain the non-dimensional governing equations of the melt from the overall conservation of mass, the balance of momentum, the balance of energy, and the conservation of species mass, ∂T ∂t where the coordinates (r, θ, z) and time t are, hereafter, nondimensional, and u = (u r , u θ , u z ), p, T, and C denote the dimensionless flow velocity, pressure, temperature, and silicon concentration, respectively. The associated Schmidt and Prandtl numbers are defined by Sc = ν/D and Pr = ν/α, where α and D represent, respectively, the melt thermal diffusivity and the diffusion coefficient of Si in the melt. The boundary conditions on the walls are non-slip (u = 0) for flow velocity components, the top cold disk surface (z = A s ) has T = 0 and C = C Si = 1, the bottom hot disk surface (z = 0) has T = 1 and C = 0. The conditions for the free surface at r = 1 are Here, we choose the same definition of the thermal and solutal Marangoni numbers as in the literature [6,19] Ma T and Ma C which, respectively, represent the magnitude of the surface tension due to the temperature and concentration gradients relative to the viscous force, where σ T and σ C denote the surface tension coefficients of the temperature and concentration fields, respectively. As shown in Figure 1, the directions of Marangoni forces along the free surface are the same direction in the present study. The grid resolution used in the present study is N r × N θ × N z = 40 × 40 × 60 and the details of the nonuniform grid refinement, numerical schemes, and the validation of transient OpenFOAM solver can be found in our previous study [20]. The physical properties of Si x Ge 1−x liquid bridge are from Abbasog et al. [21]: the kinematic viscosity ν = 1.4 × 10 −7 m 2 /s, the thermal diffusivity α = 2.2 × 10 −5 m 2 /s, the diffusion coefficient D = 1.0 × 10 −8 m 2 /s, which corresponds to Pr = 6.37 × 10 −3 and Sc = 14.

Dynamic Mode Decomposition (DMD)
Dynamic mode decomposition essentially gives a coupled system of spatial temporal modes, and it was recently used to extract an approximate periodic system from unknown data samples. The Marangoni convection along the free surface of the liquid bridge gives rise to the development of periodic flow structures and DMD is useful to analyze the dynamical behavior of such structures.
A temporal sequence of data snapshots from numerical simulation are given by a matrix X N 1 , where x i denotes the i-th flow field (N = 1, 2, 3, ..., N). A constant sampling time ∆t of the data sequence is chosen according to the Nyquist criterion (see in [17] for details). A linear mapping A between a flow field x i and the subsequent flow field x i+1 is which is considered as a linear approximation of solving the governing equations. The dynamic characteristics of the system is given by the eigenvalues and eigenvectors of the matrix A. Furthermore, Equation (11) can be written in the matrix form where r is the residual vector and e N−1 ∈ R N−1 . The singular value decomposition (SVD) is applied since SVD is robust to numerical errors and noises in the data. Substituting the SVD of X N−1 1 (= UΣW H ) into Equation (12), then with U contains the proper orthogonal modes of X N−1

1
. As the A related toÃ via a similarity transformation, the dynamic modes Φ i can be expressed as where y i is the i-th eigenvector ofÃ, i.e.,Ãy i = µ i y i . The logarithmic transform of µ i , λ i = log(µ i )/∆t (15) represents the stability of extracted dynamics oscillatory mode, where, the real part of the eigenvalues, [λ], represents the growth/decay rate, and the imaginary part, [λ], represents the frequency (wavelength) of each mode. For each case, 250 snapshots of the instantaneous concentration fields with a sampling rate ∆t = 2 were processed in the present study using the DMD algorithm.

Results and Discussion
The flow diagram (Ma C , Ma T ) for three different aspect ratios A s = 0.25, A s = 0.5, and A s = 1 are summarized in Figure 2, including the previous DNS results of Minakuchi et al. [22] for A s = 0.5. The results indicate that, when Ma T and Ma C are sufficiently small, the flow in the liquid bridge is steady and axisymmetric (m = 0). When Ma T as well as Ma C are larger, the flow becomes 3D oscillatory or 3D steady depending on the aspect ratio of the liquid bridge.  Figure 2. This direct transition from an axisymmetric steady to a 3D time-dependent flow has been reported in many numerical works [3,4] for high-Pr number fluids. In this study, due to high-Sc number, the primary bifurcation corresponds to the one of high-Pr number fluids. On the other hand, the Marangoni flow of low-Pr number fluids becomes oscillatory via a 3D steady flow when considering pure thermal Marangoni convection [5,19].
For a shorter liquid bridge, A s = 0.25, when Ma C increases, m-fold symmetric structures for concentration distributions in the middle plane as in Figure 3a,c, and the azimuthal wavenumbers are m = 11, 12, and 13 depending on Ma T . The azimuthal wavenumbers m = 4, 5, 6, and 7 dependent on Ma T have been reported by Minakuchi et al. [22] at A s = 0.5. Figure 4 depicted that the azimuthal wavenumber of oscillatory modes drastically decreases as the aspect ratio increased. The high-frequency wave mode appears, because there is not enough space for the large circulation in the liquid bridge for a flatter aspect ratio. According to Lappa et al. [15], in the case of a short liquid zone, Marangoni convection (known as surface phenomena) is confined to a region near the free surface which prevents the coalescence of the opposite convection cell. Thus, several layers of vortices emergence by continuity in the interior of the liquid bridge which leads to the occurrence of the traveling waves with high frequency. A similar dependence of azimuthal wave numbers of oscillatory modes (m) on A s has been observed in previous studies in the case of pure thermal Marangoni convection [6]. The relationship between m and A s obtained from the present study is of a roughly constant value, mA s = 2.0-2.75 ( Figure 4). This range is consistent with those obtained by Chen et al. [6] for the Pr = 0.02 fluid with 0.4 ≤ A s ≤ 2.0, and is also in good agreement with those obtained by Wanschura et al. [7] especially in the present range, 0.5 ≤ A s ≤ 1.3. It is evident that mA s becomes nearly constant when the aspect ratios are relatively smaller than A s = 1. The present results also suggest that solutal Marangoni convection does not affect the azimuthal wavenumber of oscillatory modes at small Ma T ≤ 1800A −1 s . In the region, Ma T ≤ 1800A −1 s , the critical values, Ma C,Cri , are dependent on the aspect ratio of the liquid bridge. Figure 5a shows an detailed example at Ma T = 0 and Ma T = 1400. The Ma C,Cri decreases with A s , which dependency of Ma C,Cri on A s corresponds to the results of pure thermal Marangoni convection for a high-Pr number (Pr = 1) fluids [3] shown in Figure 5b as dashed line.
According to Chun et al. [23], instability mechanism for low-Pr number fluids is not based on the amplification of temperature disturbance because of high thermal diffusivity leads to spread out any possible disturbance over the liquid bridge. In the present study, due to high-Sc number, a concentration disturbance leads to a disturbance of the concentration gradient and thus to a surface tension gradient that triggers a distortion on the velocity field. The distortion of the velocity induces distortion of the concentration field. This coupling mechanism is responsible for the growth or decay of the initial concentration disturbance. Thus, the instability mechanism emerged as hydrosolutal in nature.
The frequency spectrum of the time-dependent concentration at a sampling point of (r, θ, z) = (0.75, 0, 0.5A s ) when (Ma T , Ma C ) = (9100, 3572) in Figure 6 shows that the azimuthal traveling waves are associated with harmonics of their fundamental frequencies.
The dynamic mode decomposition results are used to determine fundamental frequencies of azimuthal traveling waves developed in the liquid bridge. Representative dynamic modes for A s = 0.25 are displayed in Figure 7 using the concentration field in the middle plane (z = 0.5A s ); their respective eigenvalues and calculated frequencies are given in Case (i) of Table 1

Critical Ma T under a Weak Solutal Marangoni Convection
When Ma C ≤ 700A −1 s as Ma T increases, the axisymmetric steady flow becomes 3D oscillatory for shorter liquid bridges (A s = 0.25 and 0.5). We note that the flow transition for a long liquid bridge (A s = 1) is from axisymmetric steady (Figure 9a) to 3D steady (Figure 9b). The critical thermal Marangoni number Ma T,Cri decreases with the aspect ratio of the liquid bridge and a similar trend can be observed for the results of pure thermal Marangoni convection by Imaishi et al. [19] as shown in Figure 5b. Under the presence of weak solutal Marangoni convection, the critical thermal Marangoni number is roughly Ma T,Cri A s ≈ 3200 at the small solutal Marangoni numbers, 400 Ma C A s 800, as shown by the dashed lines in Figure 2a,c.
The higher azimuthal wave numbers (m = 6 in Figure 3d and m = 7 in Figure 3e) are observed for a short liquid bridge (A s = 0.25) as Ma T increases due to the limited space inside the liquid bridge, therefore the characteristic length for A s 1 is the height of the liquid bridge L. Figure 3 shows that the flow patterns with m = 6 and 7 do not distinctly appear compared to m = 11 and 12. This observation is confirmed by Figure 8 as the concentration values corresponding to m = 11 and 12 are relatively higher than those corresponding to m = 6 and 7. Figure 10 shows the representative dynamic modes (m = 6) corresponding to second leading eigenvalues (λ 2 ) for (Ma T , Ma C ) = (13,300, 1786), (13,300,2144), and (13,300, 2500) at A s = 0.25. For these three cases, dynamic modes corresponding to dominant eigenvalues (λ 1 ) are shown m = 7 (see Figure 7e). It is concluded that the azimuthal traveling waves of m = 6 (passively) and m = 7 are superimposed when (Ma T , Ma C ) = (13,300, 1786), (13,300,2144), and (13,300, 2500).
For A s = 1, when Ma C ≥ 1250 as Ma T increases, after a certain threshold, the steady axisymmetric flow departs to rotating oscillatory flow with m = 2 (Figure 9c), and eventually, the flow becomes 3D steady (Figure 9b) while the azimuthal wavenumber remains unchanged if Ma T further increases. The extracted spatio-temporal coherent structures ( Figure 11) and respective frequencies (case(ii) of Table 1) of traveling waves associated with rotating oscillatory flow regime are computed by DMD. The non-dimensional concentration with respect time calculated by DNS for (Ma T , Ma C ) = (3150, 1490) and (1350, 1786) in rotating oscillatory flow regime at the sampling point of (r, θ, z) = (0.5, 0, 0.5A s ) and A s = 1 is presented in Figure 12

Conclusions
The dependency of flow regimes of thermo-solutal Marangoni convection in the halfzone liquid bridge on aspect ratios (A s ) has been investigated using three-dimensional numerical simulations under zero-gravity. The dependency of Ma C,cri on A s for weak solutal Marangoni convection is in good agreement with the previous results of the case of pure thermal Marangoni convection for high-Pr number fluids. Ma T,Cri values are decreased as A s increases, which is similar to the previous results for the case of pure thermal Marangoni convection. From our results, for a short liquid bridge (A s ≤ 1), the critical Ma T,Cri A s in the case of a weak solutal Marangoni convection is roughly constant 3200. The Marangoni convection is limited by the height of the liquid bridge which determines the characteristic length of the Marangoni convection that becomes comparable to the radius with increasing A s (the flow behavior changes at A s ≈ 1), consistent with the previous analysis of the pure thermal Marangoni convection in tall liquid bridges with A s 1.
The flow pattern in the flat liquid bridge is analyzed by dynamic mode decomposition (DMD). The azimuthal wavenumbers of oscillatory modes are increased when the aspect ratio of the liquid bridge gets smaller. From our results, mA s ≈ 2.0-2.75 is obtained, which is a similar relationship obtained in previous studies in the case of pure thermal Marangoni convection. The azimuthal wavenumbers of oscillatory modes are insusceptible of solutal Marangoni convection.