Numerical Modeling of Heat Transfer and Thermal Stress at the Czochralski Growth of Neodymium Scandate Single Crystals

The Czochralski growth of NdScO3 single crystals along the [110]‐direction is numerically analyzed with the focus on the influence of the optical thickness on the shape of the crystal–melt interface and on the generation of thermal stresses. Due to lack of data, the optical thickness (i.e., the absorption coefficient) is varied over the entire interval between optically thin and thick. While the thermal calculation in the entire furnace is treated as axisymmetric, the stress calculation of the crystal is done three‐dimensionally in order to meet the spatial anisotropy of thermal expansion and elastic coefficients. The numerically obtained values of the deflection of the crystal/melt interface meet the experimental ones for absorption coefficients in the range between 40 and 200 m−1. The maximum values of the von Mises stress appear for the case of absorption coefficient between 20 and 40 m−1. Applying absorption coefficients in the range between 3 and 100 m−1 leads to local peaks of high temperature in the shoulder region and the tail region near the end of the cylindrical part.


Introduction
NdScO 3 belongs to the group of the rare-earth scandates with distorted perovskite structure. Thanks to structural, chemical and thermal compatibility, rare earth scandates are suitable substrates for the epitaxial growth of perovskite thin films. [1] Hence, their main application is as substrates for coherent epitaxial growth of perovskite functional films and heterostructures. Substrates can be selected to provoke a certain misfit allowing an engineering of the strain in the growing film. This strain is able to significantly affect the physical properties of the film. It has been shown that biaxial compressive strain can increase the ferroelectric phase transition temperature of BaTiO 3 thin films by nearly 500 K compared with bulk material. [2] SrTiO 3 in unstrained form is not ferroelectric at any temperature but thin films grown on DyScO 3

DOI: 10.1002/crat.202000106
show ferroelectric behavior at room temperature. [3] Among the rare-earth scandates that can be grown routinely as single crystal in adequate size, NdScO 3 is the one with largest pseudocubic lattice parameter. [1] To date, rare-earth scandates with larger lattice parameter, i.e., PrScO 3 , CeScO 3 , and LaScO 3 , cannot be grown by the conventional Czochralski technique because of their high melting point temperatures. [4] As most other perovskite-type rareearth scandates, NdScO 3 melts congruently. NdScO 3 single crystals can be grown by the conventional Czochralski method with RF-heating and automatic diameter control. The most challenging property of NdScO 3 crystals is probably the semitransparency concerning internal radiative heat transfer. As pointed out in papers as, e.g., ref. [5] there is a strong influence of the internal radiation on the temperature distribution of growing bulk crystals, e.g., from experiments with BGO (Bi 4 Ge 3 O 12 ) it is known that due to the heat losses by internal radiation via the crystal the shape of the crystal-melt interface turns out deeply convex at the initial stage of the growth and after a certain crystal length the convexity is strongly reduced.
Currently, NdScO 3 Czochralski crystals can be successfully grown at our institute up to a size of 45 mm length and 18 mm diameter. During the process of wafer manufacturing sometimes cracking occurs. We suppose that this is due to the effect of the internal stress field which was developed during the thermal history of the growth process. In the vicinity of the interface there is a thermal stress caused by the inhomogeneous temperature field (second spatial derivatives are finite). However, because of the high temperature in this region the stress can be mainly relaxed by plastic deformations such as sliding or multiplication of existing dislocations (the critical shear stress is low). [6] At later growth stages, the formerly newly grown crystal area moves away from the hot region. Toward lower temperatures the critical shear stress for plastic deformations strongly increases and eventually exceeds the local thermal stress. This causes a consolidation (freezing) of the formed dislocation structures, and not only when room temperature is reached the dislocations themselves cause an uneven thermal expansion and thus the so-called internal thermal stresses. These roughly correspond in size to the original tensions in the hot crystal, which triggered the plastic processes but with a reversed sign. [7]  Thus, even the computation of thermal and thermoelastic stress fields for steady-state situations during the growth process will give a first knowledge on the residual stress in the final crystal. Such computations are the focus of this paper. What we do not cover by these calculations is the effect of time: the longer a crystal region remains at high temperatures, the smaller will be the remaining level of thermally induced strain and stress, which finally enters the stage of freezing. [6] A temporal analysis is also the prerequisite to perform computations including the dislocation dynamics, as, e.g., by the Alexander-Haasen model. [8] Nevertheless, the steady-state equilibrium analysis of thermal stresses during the growth of oxide crystals has recently received much attention. In particular for the semitransparent sapphire have recently been carried out 3D stress analyses, both for growing according to the Kyropoulos process [9] as well as according to Czochralski and the edge-defined film-fed (EFG) process. [10] In both articles the focus is on the use of anisotropic thermoelastic material properties, in ref. [9] the relevant results are also compared with the results from isotropic data: significant differences were found in their respective distributions pattern.
With our investigations we focus on the influence of the optical absorption coefficient and of the presence or absence of a shoulder before the cylindrical part of the crystal starts. The structure in the present paper is as follows: in Section 2 the crystal growth model is presented together with the computational procedure, afterward in Section 3 the used material properties of NdScO 3 are compiled. The results for temperature field and thermal stresses will be explained and discussed in Section 4.

Crystal Growth Model, Computational Procedure, and Exemplary Results
A typical NdScO 3 crystal grown at IKZ is shown in Figure 1. Such crystals are grown with a pulling rate of 1 mm h −1 and a crystal rotation rate of 10 rpm. As seen in Figure 1, the deflection of the crystal/melt interface is about 3-4 mm. The frequency of the RF induction heating is fixed to 80 kHz.

Model of Crystal Growth Arrangement
The crystal growth furnace is a conventional Czochralski furnace with RF-induction heating and automatic diameter control. A schematic diagram of the setup is shown in Figure 2. The geometry includes crystal and melt, and around the melt the crucible. Between crucible and afterheater is fixed an inclined baffle, which is to prevent spiral formation. [4,11] Due to the high melting point temperature of NdScO 3 all mentioned parts are made from iridium. The space between the crucible and the ceramic susceptor is nearly completely filled up by an isolation bulk material. The water-cooled coil made from copper is arranged around the susceptor. The most outer part is the water-cooled enclosure wall. The space inside the furnace is filled by argon.
The input weight of NdScO 3 material is 200 g. The crystal shoulder has a growth length of 5 mm, along which the initial radius of 2 mm is increased to the actual crystal radius of 9 mm. After the length of the cylinder part of 40 mm, on a growth length of 5 mm, the radius will be tailed to 7 mm. This results in a crystal weight of around 74 g. From the observation a height of the melt meniscus of 3 mm is used as a number. For a crucible radius of 19 mm, the final level of the melt is roughly 18 mm.
In addition to the full crystal length of 40 mm, calculations were performed also using half crystal length. For the whole series of absorption coefficients (see Section 3: material data) the convex deflection of the crystal-melt interface was larger than for the full crystal length case, especially clearly in shoulder growth mode. However, smaller values of thermal stresses were found, so that in this article only calculations at full crystal length will be presented.

Model Equations, Computational Procedure, and Exemplary Results
The goal of our calculations is the thermal stresses in the crystal, at the end of the actual growth process, where the crystal is still in contact with the melt. Therefore, we solve the basic set of equations for linear thermoelasticity just for the crystal (later called Step 2) but as a full 3D problem. In order to obtain the temperature field in the crystal for this computation we solve the axisymmetric equations for heat transport in the entire furnace comprising crystal, melt and gas and all equipment parts (later called Step 1): The Maxwell equations for the RF-induction process, the heat transport equation, which includes conduction and in the melt also convection, and which is coupled in the semitransparent crystal with the discrete transfer method for handling the radiation, the Navier-Stokes equation for melt and gas convection. The radiation exchange between all surfaces of the enclosure will be calculated via view angle factors.
A crucial issue is the semitransparency of the crystal. To the best of our knowledge, there is no information about the absorption coefficients of Nd-scandate available, neither as a constant nor as a function of temperature and wavelength. Therefore, we are going to vary the absorptive capacity of the crystal in the model calculations over a wide range, so that the optical thickness of the crystal is optionally optically thin, optically intermediate and optically thick. For handling this wide range of the optical thickness we need to use the radiative transfer equation (RTE) instead of typically used approximations. The RTE describes the energy transfer, i.e., the conservation of the radiation intensity, as an integrodifferential equation. This equation basically covers three items: emission of radiation of the medium itself, absorption and scattering of radiation intensity, and as the last the scattering integral over the space angle for the probability of a ray coming from one direction and being scattered into another direction. In case of the oxide crystals, scattering of thermal radiation can be neglected because the size of crystal defects is much smaller than the present wavelength of around 1 m and larger. The numerical solution of the RTE is very expensive in terms of number of arithmetic operations, i.e., in terms of time. Therefore, especially for the case of medium to large optical thickness also physically simplified equations are used instead of the RTE.
Probably the most used of this type is the so-called P 1 approximation, which is a linear equation (Helmholtz type) for the incident radiation. The P 1 equation accounts for incident and emitting radiation and also covers the scattering. Due to the physical simplification, the P 1 equation is applicable in the case of optical thickness larger than 1. An even greater physical simplification of the radiation equation represents the Rosseland approximation. This is not a separate equation, but an extension of the coefficient of thermal conductivity by a term of a highly nonlinear radiative conductivity. This term also includes the coefficients for absorption and scattering. This simple Rosseland approximation is numerically very effective, but only usable at an optical thickness about larger than 3, where the radiation propagation is more and more similar to heat conduction.
In model calculations for semitransparent crystals both the full RTE and the simplified P 1 -approximation is used. The latter is used in case that the optical properties are known with sufficient accuracy and in case that together with the characteristic length there is really a sufficiently large optical thickness: For instance, for semitransparent oxidic crystals, Stelian et al. [12,13,10] have applied a 3D numerical modelling to investigate thermal stress distribution in sapphire and langatate semitransparent crystals grown by the Czochralski (CZ) and by the edgedefined film-fed growth (EFG) technique. As method for treating the internal radiation in sapphire crystals (optically thick) they applied the P 1 -approximation, which was successfully validated by comparing the computed shape of the solid-liquid interface to the experimental results.
The both main steps of our computational procedure are structured as follows.
Step 1: Axisymmetric Computation of Heat and Mass Flow in the Entire Furnace: For solving the heat and mass transfer problem we have used the software CGSim [14] which is a tool for modelling several techniques of crystal growth from the melt. CGSim is based on the finite volume method. For the solution of the RTE, a variant of Discrete Transfer Method is implemented in CGSim. At the surface of the crystal with the transparent gas we set specular boundary condition (Fresnel type) for incident, transmitted and reflected rays.
The computational geometry is meshed by quadrilateral elements, except the argon domain and the wall of the enclosure which are meshed by triangular elements. The element size in crystal and melt is between 0.3 and 0.35 mm, a size reduction between 0.18 and 0.23 mm provides results that are only slightly changed.
When processing a data set by the CGSim software, the individual features for modelling the numerous characteristics of the Czochralski method are activated step by step (we call them runs) in order to ensure the convergence of the calculation process. In detail, the first run over, e.g., 1000 iterations includes heat conduction in all materials, radiation between all walls inside the argon filled cavity, the RF-heating in the electrically conductive components by using the initial value of the RF-current, the crystal rotation rate and the convection in melt and argon, the initial shape of the crystal-melt interface (where the melting point temperature is assigned as Dirichlet boundary condition, with the exception of the line element that touches the triple point and of the triple point itself) and the target growth rate here used as a fixed value. The turbulent viscosity both for melt and gas flow was computed by a one-equation-turbulence model; the turbulent Prandtl number was set to the typically used value of 0.9.
The second run takes the results of the first run as initialization and typically requires about 3000 iterations. In this run, the RF power is adapted to meet the melting point temperature at the triple point, where crystal, melt, and gas phase come together.
Only in the third run, which typically needs between 2 × 10 4 and 6 × 10 4 iterations, the semitransparency of the crystal and the relocation of the shape of the crystal-melt interface are switched on, too. The PID algorithm for calculating the correct RF power is still activated. The largest number of iterations is needed in case of very small absorption coefficients.
The final shape of the crystal-melt interface is found by minimizing the difference of the local crystallization rate and the pulling rate. The iterative procedure of further change of the RF current is finished when at the triple point the temperature differs from the melting point temperature by no more than an accuracy limit (10 −4 ) and when the crystallization rate differs from the target growth rate by no more than an accuracy limit (5 × 10 −4 ). Finally, this computation delivers an axisymmetric temperature field, which will be used in Step 2.
As an example, these field distributions are plotted for the core region in Figure 3, for the case of an absorption coefficient of 100 m −1 , all other material properties are as compiled in Section 3.
In the middle of the gas ring around the crystal, a temperature drop of around 300 K extends from the surface of the melt toward the roof of the inner chamber: the buoyancy flow forms here as several rolling cells, above the baffle with an absolute velocity in the order of 10 −4 m s −1 , below it with over 10 −3 m s −1 . The gas velocity is too small to have a visible influence on the temperature field. In the melt, there is a temperature range of about 210 K between the crystal-melt phase boundary and the outermost point at the bottom. This leads to an upward flow at the crucible wall and a downward flow in the axis of rotation, which here has a velocity maximum of about 0.03 m s −1 . Due to the relatively high Prandtl number of about 6, the temperature field of the melt is strongly influenced by the flow, up to the formation of temperature boundary layers.
Step 2: 3D Calculation of Thermal Stresses: After CGSim run, we take from the crystal domain the geometric and the temperature field data and expand both to 3D by rotating it by 360°. The 3D data are transferred to the COMSOL software, [15] which is based on the finite element method. The surface and volume of the crystal are meshed by triangles and tetrahedrons, respectively. The mesh is extremely fine at the edges of the geometric units and along surfaces with high curvature.
While temperature field and interface shape are still axially symmetric, the anisotropic coefficients of thermal expansion and elasticity might generate 3D displacements and thermal stresses as well as also 3D von Mises stress.
As an exemplary result, the von Mises stress for the crystal in Figure 3 is plotted in different views in Figures 4 and 5. Growth direction is [110] and the thermoelastic material data were converted with respect to this orientation (for details see Section 3).
According to the basic equations of linear thermoelasticity, thermal stress occurs when the second derivative of the temperature field is finite. This fact is obviously also reflected in the comparison of Figures 3 (right) and 4a: the largest stress values are present at the ends of the crystal cylinder, i.e., where the shoulder ends and where the tailing of the crystal begins shortly before the end of the growth process. These regions also show a more complicated temperature distribution than the central cylindrical part: this will be analyzed in Section 4.1. The maximum thermal stress at the phase boundary (Figure 4b) is about half the size of that at the remaining crystal surface (Figure 4a). In Figure 4b, the anisotropy can be clearly seen, resulting only from the anisotropic thermoelastic material data. As Figure 5a,b shows, the stress maxima are still somewhat higher on the inside than on the surface. The section along the y-axis (Figure 5a) shows the asymmetry of Figure 4b inside, while the section in x-direction (Figure 5b) is almost perfectly symmetric.

Material Properties
Due to the very high melting point temperature of NdScO 3 (2491 K [4] ) up to now only some properties of solid and melt related to the heat transfer and fluid flow have been measured. NdScO 3 has an orthorhombically distorted perovskite structure. Geometrically speaking: it has a distorted cubic symmetry, i.e., in this case the symmetry is reduced to orthorhombic symmetry (space group Pbnm, no. 62). The lattice parameters are a = 5.575 Å, b = 5.776 Å, and c = 8.003 Å. [1] The coefficient of thermal expansion was measured by Badie [16] and reported as average values over the temperature range from 1000 to 2000°C. In However, the current crystal growth direction is the [110] direction as shown in Figure 6. Accordingly, the crystal lattice The elastic tensor of NdScO 3 was measured by Pestka et al. [17] using resonant ultrasound spectroscopy (RUS) in combination with ab initio calculations that were used as a means of estimating the initial elastic constants used in RUS measurements. They In ref. [18] the transmission properties of several rare earth (RE) scandates were investigated, among them NdScO 3 , through polished (110) wafers. NdScO 3 has several absorption bands both in the visible and in the infrared range. These bands correspond with the energy levels of Nd 3+ . Beyond a wavelength of 8000 nm NdScO 3 is opaque. The absorption spectrum is quite complicated and moreover, it was obtained at room temperature. Therefore, we kept things simple and did not try to feed some transparency bands into our software. Instead we have taken the following procedure: According to the characteristic length L ≈ 22 mm of a typical crystal, which was calculated by the third root of the crystal volume, we have defined different value ranges of a grey absorption coefficient in order to meet the various branches of the dimensionless optical thickness : for = 0.01-3 [m −1 ] we get < 0.1, i.e., the medium is optical thin, for = 400, 800 [m −1 ] we get > 10, i.e., the medium is optical thick, and for = 10-200 [m −1 ], the medium is modelled as moderately optical thick.
The refractive index as a constant for all temperatures was estimated to be about 2, based on detailed knowledge of the data of Sc 2 O 3 and DyScO 3 . [19,20] The heat transfer and fluid flow related properties of solid and liquid NdScO 3 are not yet fully measured. From ref. [4], we take the enthalpy of fusion ΔH f = 89.2 kJ mol −1 and the heat capacity of solid NdScO 3 c p = 630 J (kg K) −1 . The heat conductivity of solid NdScO 3 was fixed to the isotropic value of = 2.22 W (m K) −1 , as taken from DyScO 3 . [21] However, except the case of very large optical absorption, the heat transfer by radiation is by far dominating over the heat transfer by thermal conduction, as can be shown by the Stark number: [22,23] The Stark number is the ratio of the energy transport by conduction to that by radiation: The NdScO 3 melt is assumed to be radiatively opaque with a surface emissivity = 0.4. The other melt properties are taken as those from the DyScO 3 melt published in ref. [24]: heat conductivity = 4 W (m K) −1 , heat capacity c p = 500 J (kg K) −1 , dynamic viscosity = 0.04986 Pa s and thermal expansion coefficient = 4 × 10 −5 K −1 . These data result into a Prandtl number of Pr ≈ 6.2 which is a typical value of an oxide melt.

Results and Discussion
First, we analyze the temperature fields calculated as a function of different absorption coefficients. The temperature field in the growth arrangement is imprinted by the balance of heat input and output. Feeding is from the inductively heated crucible via the boundary of the melt and by the release of the heat of solidification at the melt-crystal phase boundary.
The heat losses are caused by heat transport to the colder inert gas and to the colder housing parts in the upper part of the growth arrangement. The heat input is governed by the constraint that at the triple point the temperature should be the melting point temperature. The removal of the heat is greatly influenced by the radiative transport in the crystal and the adsorption coefficient does not only trigger the overall heat removal but also the local heat fluxes in the crystal. Radiative heat transfer is not only defined by material parameters such as the adsorption coefficient but also by the geometrical boundary of the semitransparent object, i.e., the crystal. In this sense we are interested in the influence of the geometry of the shoulder region on heat transport and stress distribution. Therefore, we consider besides the normal situation of a distinct shoulder between the seed and the crystal cylinder ("shoulder growth") also the case when after the seed dipping the diameter of the crystal cylinder is reached very quickly ("Flat top surface growth"). Figure 7 shows the axial temperature profiles (i.e., at r = 0) of melt and crystal; the profiles start at the bottom of the crucible (left) and go up to the top of the seed (right).

Case: Shoulder Growth
Due to the differently sized deflections of the convex crystal/melt interface shape, the axial temperature profiles intersect the melting point temperature at different values of the zcoordinate. The profiles are shown for various absorption coefficients ranging from 0.01 to 800 m −1 .
In the case of a large absorption coefficient such as 200 m −1 , the radiation intensity is reduced to a about 1/3 after only 5 mm of travel, i.e., in relation to the entire crystal there is only a weak radiation transport in addition to the immanent heat conduction. This leads to a relatively large temperature gradient throughout the crystal. In case of a small absorption coefficient, e.g., 10 m −1 , even after the full path of the characteristic length (about 22 mm) only 20% of the radiation intensity is absorbed. This strong radiation transport accompanies the heat conduction and results in a smaller vertical temperature gradient in the crystal.
In Figure 7 there are two more conspicuous features. The first one: the temperature profile is not monotonically decreasing for each parameter value. For some values, local temperature peaks  appear in both the shoulder and tail sections. In the tail section the peak temperature exceeds even the melting point temperature. It should be noted that first in the calculation the physical properties of the solid were used everywhere. Second for a remelting latent heat is required which would lower the temperature locally. We will discuss these temperature peaks in more detail below.
The second peculiarity: For the absorption coefficients 0.01 and 1 m −1 Figure 7 shows a strong temperature drop from a high temperature in the melt toward the melting point temperature at the crystallization front. The origin of this drop can be better seen in Figure 8a, where the temperature isolines in crystal and melt are shown.
The interface is going deep into the melt with a blunt end at the bottom. This geometrical shape causes a counter rotation to the overall dominating buoyancy convection induced by the crystal rotation, which means a thin boundary layer is formed. The boundary layer has a thickness of about half a mm and spans a temperature difference of about 70 K for = 0.01 m −1 . The boundary layer is numerically resolved with 6 finite volume mesh layers. In the crystal the temperature gradient is small, much smaller than in the melt, because the heat removal from the interface is mainly governed by radiation in the nearly transparent crystal. The temperature gradient is about −10 K cm −1 in the cylindrical part of the crystal.
For the cases with a moderate (Figure 8b) and a high absorption coefficient (Figure 8c) the situation is quite different. In the latter situation we observe a steep temperature gradient in the crystal (about −50 K cm −1 ) because the heat is transported mainly by conduction. Interestingly, the total heat input by RF heating is only slightly different from that for the case in Figure 8a (4% less). The most intriguing case is that of Figure 8b: For the crystal region near the interface, the temperature calculation results in a local maximum. This means that the amount of heat entering Figure 9. Starting with the case of Figure 8b (which is drawn here again in (a)): the change to straight cylinder geometry and to reduced refractive indices causes the local temperature peak to disappear, isoline gap is ΔT = 1 K.
at the melt-crystal interface cannot be transported to the cooler areas sufficiently quickly. This is an effect of the nonlinear behavior of the internal radiation. As mentioned earlier the radiation in the crystal is influenced by the boundary of the crystal, i.e., by its shape and also by the refraction index at the boundary. We demonstrate this influence by computing the temperature field for a geometry without slimming at the tail (Figure 9b) and for different refraction indices (Figure 9c,d).
A straight cylinder slightly reduces the effect: In case of a straight cylinder geometry ( Figure 9b) the local peak is only half as large, only 7 K instead of 15 K above the melting point temperature (Figure 9a). A further reduction of the peak can be achieved by reducing the refractive index from 2 to 1.5: the peak size is going from 7 to 1 K above melting point temperature (Figure 9c).
Finally, in the case of a refractive index of 1.33 there is no temperature peak at all (Figure 9d). The reduction of the refractive index means a decrease of reflection at the crystal boundary. In this respect, the above results correspond to the finding of ref. [25] that the Fresnel reflection leads to a distortion of the isotherms. The isoline distributions in the crystal shoulder shown in ref. [25] are very similar to those parts in Figures 9d and 3 (right). Because in the case of the shoulder the relatively low temperature spreads from the shoulder tip far into the shoulder volume (around the axis), in ref. [25] it is called a cooled domain.
However, the phenomenon immediately above the crystalmelt phase boundary is definitely a kind of overheating. This phenomenon of enclosing of the radiation energy in a bulky medium (i.e., medium with moderate aspect ratio, not in thin films or similar) is understandable in terms of basic radiation physics (see, e.g., refs. [26,27]).
The moderately transparent medium absorbs radiation internally and also emits radiation internally from its volume. Some part of that radiation reaches the semitransparent surface and passes into the environment. The rest of the radiation will again be partially absorbed, repeatedly reflected at the surfaces, and so on. According to ref. [27], internally emitted radiation is almost perfectly diffuse, hence, a considerable part of this radiation reaches a surface at angles larger than the critical angle of total reflection.
In turn, the critical angle of total reflection gets smaller with increasing refractive index. Hence, when simulating the current problem but with assumed smaller refractive indices (as in Figure 9) obviously a larger amount of radiation can escape into the environment which causes the overheating to disappear (Figures 9b-d).

Case: Flat Top Surface Growth
Now we discuss the case, where the crystal is grown quickly from the seed to the cylindrical part resulting in an almost flat top of the cylinder. This configuration requires only very little more induced power. Even in the extreme case of the smallest absorption coefficient (0.01 m −1 ) the power value is higher by only 1.4%. But there is a clear effect on the axial temperature profile. As Figure 10 shows, the local temperature peak in the vicinity of the seed (in Figure 7 for shoulder growth) disappears for = 1-100 m −1 . Indeed, the local peak near the seed is not levelled out but the "valleys" along the axis of the crystal cylinder are filled up due to the slightly enlarged heat throughput. As a result, the temperature in the melt is raised and the buoyancy flow is stronger, both reducing the crystal/melt interface deflection and altering its shape (Figure 11, especially part (a)). Additionally, the magnitude of the local peak just above the crystal/melt interface and the temperature of the cylindrical part is enlarged for all values. Figure 11b,c is nearly identical to the corresponding figures of the shoulder growth because the Grashof number of the buoyancy convection is enlarged by only less than 1%.

Comparison to Experiments
As mentioned at the beginning of Section 2, crystals grown at our institute have an interface deflection in the range Δz = 3-4 mm. Compared to those experimental results, only the computations with the absorption coefficient in the green colored section in Figure 12 lead to the same interface deflections. Thus, assuming an absorption coefficient within the range between 40 and 200 m −1 gives realistic results. This range roughly meets those values which were used by other authors for simulations: 90 m −1 for langatate, [12] two wavelength bands of 20 and 254 m −1 for sapphire, 3 and 95 m −1 for Bi 4 Ge 12 O 3 , and 10 m −1 and 347 m −1 for B 12 SiO 20 (all in ref. [25]). www.advancedsciencenews.com www.crt-journal.org  . For flat top surface growth: temperature field for selected cases of optical absorption coefficients, the isoline gap is 2 K for the crystal, and 10 K for the melt. However, in order to cover the overheating in the crystal, the change of gap is in a) at 6 K above Tmp and in b) 16 K above Tmp.
As analyzed in Figure 11, the absence of a shoulder provides a smaller deflection of the crystal/melt interface. In more detail this can be seen in Figure 12a. However, in the interesting (green) region of absorption coefficients the interface deflection is almost the same (Figure 12a ) as well as the maximum von Mises stress at the crystal/melt interface ( Figure 12b) and at the ab-plane (Figure 12c). Thus, avoiding the shouldering of the crystal geometry does not bring any benefit. Similarly, in terms of a gray absorbingemitting material, if NdScO 3 has really a coefficient within the green marked range, then the effect might be rather small.

Conclusion
Based on axisymmetric temperature field and crystal/melt interface computations at the Czochralski method, we have figured out 3D thermal stresses in NdScO 3 crystals. In order to come up with the lack of true data we have modelled the optical absorption property of the crystal by the large range of coefficients between 0.01 and 800 m −1 . Applying absorption coefficients in the range between 3 and 100 m −1 leads to local temperature maxima in the shoulder region adjacent to the seed and in the tail of the cylindrical part at the end of the crystal. The occurrence of these maxima can be explained by the specific physics of thermal radiation. In some cases we even observed a temperature exceeding the melting point temperature just above the crystal/melt interface. In real crystal growth experiments, no indication of remelting inside the growing crystal has been observed yet. However, a remelting can be observed only in a transient calculation allowing a phase transition also inside the crystal. Further studies are needed to elucidate this issue.
The experimental melt/crystal interface deflection values agree with computational cases which have a grey absorption coefficient in the range between 40 and 200 m −1 . However, large values of thermal stress appear just in case of intermediate optical absorption coefficients (1-200 m −1 ). Hence, this might be part of an explanation for large residual thermal stresses that sometimes lead to mechanical cracking during the process of wafer preparation.