Resistivity as a Witness of Local Crystal Growth Conditions

. A detailed knowledge of the growth front geometry during physical vapor transport (PVT) growth of SiC single crystals is beneficial to achieve high quality n+ SiC substrates for power device applications. In this report we show how mapping of resistivity in SiC wafers can shed light on local growth conditions, which are very difficult-to-study in situ . We consider both thermodynamic quantities (absolute temperature T and partial pressure p N₂ ) and geometric characteristics of the growth surface relevant to growth kinetic parameters, namely atomic terrace width and atomic step velocity. Specifically, we show how an elevation map of the growth surface can be reconstructed from a spatially-resolved measurement of resistivity in a SiC wafer by integrating the spatial derivative of elevation with respect to the basal plane, which is assumed to be related to local resistivity through dependence of non-equilibrium nitrogen incorporation on the atomic step velocity.


Introduction
Silicon Carbide (SiC) semiconductor technology can help maximize the efficiency of power systems and simultaneously reduce their size, weight, and cost as compared to legacy silicon-based systems.Single crystal wafers of 4H SiC are the top choice for power electronic device substrates due to their excellent thermal and electrical properties.Demand and market size projections for SiC based systems continue to grow significantly.Recent forecasts predict that the SiC device market will grow beyond $6 billion by 2027 from $1 billion in 2021 [1].Many aspects of the SiC substrates' quality are affected by the growth front geometry during physical vapor transport (PVT) growth of SiC single crystals.Thus knowledge of this geometry throughout the growth of the crystal (boule) is beneficial to optimization of the growth process and growth equipment.Due to extremely high temperatures required for PVT growth of SiC the growth front geometry is very challenging to measure in situ with sufficient (sub-mm) resolution.An alternative of interrupting the growth process to obtain a snapshot of the growth front from the surface of an undergrown boule is commercially unappealing, especially to obtain large datasets sufficient for statistical analysis.The purpose of the presented work is to explore if spatially resolved measurements of substrate resistivity in SiC wafers can be used to reconstruct the growth front geometry evolution throughout the growth process.

Experimental
Single crystal 4H n-type SiC boules were grown by a seeded PVT method.The crystals were then fabricated into 150 mm wafers via grinding, slicing, and polishing [2].The resistivity of the wafers was mapped with 4 mm spatial resolution within a 60 mm radius circle concentric with the wafer using a tool RT-2200E by Semilab.The luminosity of an entire wafer was mapped in an Epson Perfection V850 Pro Photo scanner operated in the transmission mode.The luminosity profile extraction and data analysis were performed using ImageJ and JMP 17 software respectively.

Results and Discussion
Thermodynamics empowers us to form an equation describing equilibrium value of resistivity ρ∞ dependence on absolute temperature T and nitrogen partial pressure pN₂.From the equation of nitrogen incorporation into SiC: and the Van 't Hoff equation (ln K = ΔS/R -ΔH/RT ) we built the following ρ∞(T, pN₂) formula with two adjustable parameters A and B related to ΔS and ΔH of the reaction shown in Eqn. 1 respectively.Here we assume that resistivity is solely determined by the concentration of dissolved nitrogen atoms, i.e. ρ∞ ∝ 1/[NC]: ln ρ∞ = ln A -0.5 ln pN₂ -B/T. ( If only thermodynamics were needed to predict the local observed resistivity, we would not see a sharply defined zone of substantially lower resistivity in the basal plane (b.p.) facet area.This area is visible as an oval on the right side of the resistivity maps and transmitted light images of wafers made from boules grown off-axis, as shown in Fig. 1.
Apparently, the local geometry of the growth front is affecting nitrogen incorporation in the growing crystal.If we assume that a faster moving atomic step is less careful about rejecting impurity atoms of nitrogen adsorbed on the growth surface from being incorporated into the growing SiC crystal, we can infer the following.For a given normal growth front velocity, i.e. the rate of growth of SiC crystal perpendicular to the growth surface, more nitrogen gets incorporated in areas with lower angle between the local normal to the macroscopic (after smoothing the atomic steps) growth surface and the c axis of the crystal, i.e. the inclination angle α, and thus lower magnitude of gradient of elevation z , i.e. |∇z|, where z is measured along the c axis of the crystal, i.e. perpendicular to the basal plane.In this report, the elevation axis z pointing along the c axis of the crystal is tilted by angle θ toward the [1 � 1 � 20] direction with respect to the boule growth axis h perpendicular to the crystal seed plane and pointing along the boule growth direction.Areas with lower |∇z| have higher atomic terrace width t, which requires the atomic steps to move perpendicular to c axis of the crystal with a higher velocity v in order to yield the same normal growth front velocity along axis h [3].
To gain support for the above assumed relation between |∇z| and local resistivity ρ, we performed the forward computational task first.Based on numerous observations of grown boule surfaces, we assumed equation of a body of revolution around axis h, i.e. h = h(x,y) = h(r), where r is the radius in cylindrical coordinates (r, φ, h) given by r = √(x² + y²) and h(r) is an even-powers-only polynomial with respect to r that yields a realistic growth front surface.Due to a low value of angle θ we can approximate rotation of a surface in space by a mathematically simple "shear" distortion formula stating that the two surface equations z(x,y) and h(x,y) are related by a linear with respect to x term responsible for the tilt by the angle θ, namely

20
Growing and Forming of Semiconductor Layers where const is an immaterial constant that does not affect the shape, but only offsets the surface vertically.Here, using approximation cos θ ≈ 1 when θ is close to zero, we neglect the differences between the [1 � 1 � 20] direction and axis x formed by projecting this direction onto the plane perpendicular to h.The y axis perpendicular to h, z and x remains unchanged, because the tilt by θ occurs around this axis.Thus calculated z(x,y) surface is shown in the top of Fig. 2 viewed along the y axis.Right below it in black we show h(x) cross-section at y = 0 with highly exaggerated (you can think of it as step bunching with a constant macro-step height) atomic steps to show the variation of terrace width t.Note multiple local minima of t along the x axis.Since we formed an analytical expression for z(x,y) surface, we differentiated it analytically and thus arrived at the expression for |∇z| as a function of x and y, which is graphed as a filled contour map in the bottom right of Fig. 2. Note that |∇z| map in Fig. 2 qualitatively matches Fig. 1, which lends support to the main assumption of this work: the larger the spacing between the atomic steps or contour lines of z, i.e. the lower the |∇z|, the lower the local resistivity of the growth crystal.The b.p. oval is not sharp and not as intense in our simulation, because we ignored spontanteous formation of an almost-perfectly flat (|∇z| ≈ 0) b.p. facet during growth and simplified the boule surface geometry by using a smooth (edgefree) curved surface formula with a polynomial profile.Now we have to find a way to convert experimentally observed maps of resistivity (tables of ρ for varying x and y) into h(x,y) maps by performing the reverse computational task, namely integrate spatially varying |∇z| calculated from ρ (see the next two paragraphs) over the measured area to arrive at z, and finally, using Eqn.3, produce the spatial distribution of h.This task is not trivial, because to integrate z over an area in the x-y plane we need to know the values of both partial derivatives ∂z/∂x and ∂z/∂y, and not just their mathematical combination into a scalar, namely the gradient vector magnitude |∇z| = √[(∂z/∂x)² + (∂z/∂y)²].We are leaving for future work a nontrivial approach we envisioned to perform this integration by taking advantage of an often accurate boundary condition Solid State Phenomena Vol.362 h(rmax) = const, where rmax is the radius of the edge of a boule.In this report we narrow down our ambition to producing a 2D h(x) profile at y = 0, which is still valuable at describing the growth front geometry, since it can (1) capture asymmetry along [1 � 1 � 20] direction, (2) reveal local maxima and minima of h, and (3) describe how curved a growth front is near the edge.We can easily perform integration of z with respect to x at y = 0, because the enabling assumption of ∂z/∂y = 0 at y = 0 is supported by almost-always present top-bottom mirror symmetry of dopant distribution, e.g. two wafers shown in Fig. 1.One necessary step before integration is to convert the always-non-negative magnitude of gradient vector into its sole component at y = 0 is to determine the sign of ∂z/∂x.To do this we compare x to the x of the b.p. facet center xbp.If x > xbp, ∂z/∂x is negative (steps go down while x is increasing, e.g.black steps in Fig. 2 near the right edge of the drawing); otherwise the sign of ∂z/∂x is positive.This logic works fine for for ∂z/∂x = 0 observed in the b.p. facet.The xbp value is found as the x value for the data row with the lowest ρ in the data table for a given wafer.
Deriving mathematical transformation of ρ into |∇z| is difficult from first principles.The exact mechanism of incorporating excess nitrogen into the growing crystal and thus pushing ρ below the amount prescibed by the equilibrium Eqn. 2 as a function of atomic step velocity is not trivial.If we assume local equilibrium between nitrogen adsorbed on the growth surface and nitrogen dissolved in the bulk of SiC, the concentration of nitrogen adsorbed on a terrace of the growth surface must be increased near an atomic step invading a terrace at speed comparable to the speed of diffusion of adsorbed nitrogen.The situation with the local adsorbed nitrogen concentration is somewhat analogous to the phenomenon known as the "bow wave", when water level is elevated in front of a moving ship due to the pushing action of the ship's hull: the higher the ship's speed, the higher the bow wave.Just like water in the bow wave does not have time to flow away from the moving hull, nitrogen adsorbed onto a terrace pushed by the invading atomic step does not have time to redistribute and ends up piling up in front of a moving atomic step.This local pile-up of nitrogen gets "recorded" by an increased concentration of the dopant incorporated in the bulk.Even if we have a first-principles equation relating the atomic step velocity to excess nitrogen incorporation, the macroscopic growth rate determining this velocity varies throughout the growth process and from boule to boule, so we simply do not know the key input into this theoretical equation for a given wafer.
What we can do to deal with the deficiency of the theoretical prediction is to form an asymptotically correct empirical equation with minimal mathematical complexity and minimal set of parameters that adequately reflect the dependence.We define deficit of resistivity as where ρ is the observed local value of resistivity.Apparently, Δρ vanishes with increasing |∇z| as the velocity of the atomic steps vanishes, and the grown crystal approaches equilibrium.At the other extreme, within the b.p. facet where |∇z| ≈ 0 the deficit of resistivity exhibits its maximum (finite) value Δρbp.To accommodate both extremes within the same formula we can assume the trend is described by an exponential decay equation where τ is an empirical dimensionless decay constant with a graphical meaning of a value of |∇z| at which the initial (at |∇z| = 0) slope to the curve defined by Eqn. 5 intersects the horizontal axis, i.e. zero level of Δρ.It should be noted that both Δρbp and τ parameters are not material properties and should depend on the macroscopic growth rate and local temperature, which theoretically should affect the nitrogen pile-up in front of a moving atomic step, since we are dealing here with adsorptiondiffusion problem, where diffusion of nitrogen occurs both (1) between the bulk of the gas phase and the growth surface and (2) within the adsorbed layer, away from atomic steps invading terraces and pushing adsorbed nitrogen resulting in the near-step nitrogen pile-up we described in the previous paragraph.

Growing and Forming of Semiconductor Layers
It should be noted that far from b.p. facet this equation is very sensitive to error in Δρ calculation, because |∇z| diverges as Δρ approaches 0. To avoid this divergence at the expense of losing asymptotic accuracy one can use a linear approximation of Eqn. 6 formed by a partial sum of only up to the 1 st order term of its Taylor expansion around Δρbp : Next, we use top-down mirror symmetry and assumption of smooth contour lines of z(x, y) to put ∂z/∂y = 0 at y = 0 and give the ∂z/∂x partial derivative + or -sign for positions left or right of the b.p. facet's center respectively.Then we numerically integrate a table of ∂z/∂x|y=0 with respect to x to produce a table of z(x), which lets us generatre a table of h(x) using Eqn. 3. The decay constant τ can be iteratively refined by targeting an h(x) profile that is not tilting either left or right.It turns out that the tilt of reconstructed h(x) is quite sensitive to τ.In cases where temperature and/or growth rate varies from wafer to wafer for a given boule, assuming a constant τ produces a variation of tilt of the h(x) profile across a boule.
In order for our approach to work we need to estimate ρ∞ in Eqn. 4. If we know all four parameters in Eqn. 2 for a given wafer, we are all set.If we do not (as is the case for our test wafers shown below), we have to make the best of the situation by assuming that the highest observed for a given wafer value of ρ approximates the equilibrium value ρ∞ ≈ max(ρ).Caveat: for the point with maximum value of ρ Eqn.6 will diverge, so the table of results will have a corresponding empty cell.The value of constant Δρbp for a given wafer is then found as Δρbp ≈ max(Δρ).An example of the implementation of the above approach to reconstruction of the growth front profile for one wafer is given in Fig. 3.The value of τ was tuned to equal 0.039 in order to yield a near-zero slope of h(x) at the wafer center.Note how the presence of a large b.p. facet on the right makes the h profile highly asymmetrical.

Fig. 3.
In the upward order: reconstruction of growth front profile from raw resistivity ρ data to growth front profile (all profiles are at y = 0) for assumed τ = 0.039 in Eqn. 6.Note large b.p. facet manifested as a flat shelf at maximum value of z.Full heat map of resistivity is shown in the bottom right.
To scale up the implementation of our approach we reconstructed growth front profiles for 99 wafers cut from three boules (see Fig. 4).The h(x) profiles were offset vertically by a distance between a given wafer from the seed end of the boule to show to-scale evolution of growth front during growth.A global constant τ = 0.03 was assumed.In order to avoid the above divergence of Eqn. 6 and associated unrealistic profiles for wafers with near-center ρ values close to the wafer-level maximum ρ, Eqn.7 was used instead of Eqn. 6. Due to variation of growth front temperature over the course of growth of each boule, a uniform approach with a constant τ to reconstruction of h(x) profiles for multiple wafers cut from multiple boules resulted in too unrealistic variation of average slope of h(x) profiles.To suppress this variation without a tedious manual refinement of τ for every single wafer we subtracted the wafer-center value of ∂z/∂x|y=0, x=0 from every raw ∂z/∂x value (see middle row of graphs in Fig. 4) before z integration.This correction essentially individually shears in the vertical direction every raw h(x) profile until it reaches zero slope at wafer center, i.e. until ∂h/∂x|y=0, x=0 = 0.The downside of this easy-to-implement correction is deviation of slope of the b.p. facet (shelf of z(x) profile) from ideal zero value in z(x) profiles.

Fig. 4.
In the upward order: reconstruction of growth front profiles for three boules from raw ρ data (all profiles are at y = 0) for assumed τ = 0.03 in Eqn. 7. The color encodes Slice #, i.e. the distance from a given wafer to the seed end of an ingot measured in units of the multiwire saw pitch.Fig. 4 demonstrates a diversity of growth front geometries and their evolutions from typical dome (convex) shape at the early stages of the growth, which later during growth can take alternate paths of either convex shape observed for some boules or non-convex shape observed for other boules.An occasional physically impossible (due to continuous addition of SiC everywhere on the growth front) intersection of growth fronts for adjacent wafers can be explained by complexity of time profiles of partial pressure of nitrogen and growth front temperature.Moreover, growth front temperature for a given wafer can exhibit an appreciable difference between the center and the periphery.This temperature non-uniformity further distorts the results of our simplistic reconstruction shown in Fig. 4, which assumes that ρ only depends on |∇z|.The reconstruction can be improved by incorporating knowledge of growth front temperature and its distribution at different stages of growth as well as the time profile of nitrogen partial pressure to be input into Eqn. 2 for a better estimate of ρ∞.
The most unfortunate feature of the ρ maps we had at our disposal is their rather limited (60 mm) radius, which leaves unmapped the outer 15 mm wide margin of a 150 mm diameter wafer and even a larger margin of the growth surface.One plausible remedy for this limitation is replacing resistivity with optical density (absorbance).The latter can be measured with high spatial and tonal resolution using a flatbed scanner run in the transmission mode, e.g. the same wafer image shown in Figs. 1 and  5.The Beer-Lambert law relates absorbance of a transparent sample to concentration of lightabsorbing species.To calculate absorbance A one needs to know the ratio of absolute intensities of light passing through the sample of unknown concentration I and through the reference (undoped) sample I0 : A = log10(I0/I).The Beer-Lambert law states A = εℓc, where ε is the molar attenuation coefficient or absorptivity of the attenuating species, ℓ is the optical path length, and c is the concentration of the attenuating species.Ideally, one wants to put a piece of undoped polished 4H-SiC wafer of equal thickness and a piece of completely light-absorbing black cardboard together with the sample wafer into a flatbed scanner, so that one has internal standards for both complete light absorption and complete light transmission through the reference sample.The latter is better than just air, because it has reflection losses properly represented.With this setup and three values of luminosity L a scanner will output (L0 for reference, Lblack for complete absorption, and L for a selected sample area) the local absorbance is calculated as A = log10([L0 -Lblack] / [L -Lblack]).Unfortunately, our image shown in Fig. 1 did not come from the above optimal setup, so we could only do an uncalibrated proof-of-concept rough reconstruction where greyscale value (luminosity) counted from the darkest pixel found in the b.p. facet area.Since we use luminosity as an alias of resistivity here, we labeled it ρ, a.u. in Fig. 5, which shows a growth front geometry reconstruction analogous to one shown in Fig. 3, but with full coverage of the wafer's diameter (x spanning the range from -75 to 75 mm) and higher spatial resolution of raw data.Despite the above imperfections, the reconstructed profile appears realistic, based on our experience, which makes quick and inexpensive flatbed scanner transmitted light images a promising alternative to resistivity maps as the source of raw data for the growth front geometry reconstruction.

Summary
Elevation profile of the growth surface can be reconstructed from a spatially-resolved measurement of resistivity in a SiC wafer by first mathematically transforming resistivity onto the magnitude of the elevation gradient (where elevation is measured along the c axis of the crystal) and then integrating the spatial derivative of elevation with respect to the basal plane along the

Fig. 1 .
Fig. 1.Left: experimental resistivity ρ heat map (darker regions have lower ρ) that was the target of simulation shown in Fig. 2. Note a local minimum of resistivity between the b.p. facet oval and the wafer's center qualitatively reproduced by the |∇z| map in Fig. 2. Right: transmitted light image of a different wafer obtained in a flatbed scanner.Ignore imaging artifacts, namely interference fringes (Newton rings) formed due to a small air gap between the scanner glass and the polished wafer.Both wafers have 150 mm diameter.

Fig. 2 .
Fig. 2. Simulation of crystal growth surface using a body of revolution with an even-powers-only polynomial radial profile h(r) shown in black with exaggerated atomic steps converted into z(x,y) surface shown at the top using Eqn.3, contour lines of z(x,y) at the bottom left, and area-filled contour map of magnitude of gradient of elevation |∇z| at bottom right.The [1 � 1 � 20] direction points right.

Fig. 5 .
Fig. 5.In the upward order: reconstruction of growth front profile from raw flatbed scanner luminosity ρ, a.u.data to growth front profile (all profiles are at y = 0) for assumed τ = 0.036 in Eqn. 6. Full-wafer raw transmitted light image is shown in the bottom right.Graphs of |∇z| (blue, always non-negative) and ∂z/∂x (red, equal in absolute value, but switching sign from + to -) are overlaid.Note more realistic sharper edges of the b.p. facet due to better localization and finer spatial resolution of the light probe vs. the resistivity probe, which averages resistivity over a diffuse spot several mm in diameter.