%% %% =================================== %% OModel.m %% by Bart Trzynadlowski, 2011-2013 %% =================================== %% %% Model for oxygen precipitation in silicon. Includes both the full and reduced %% kinetic precipitation models described in my doctoral dissertation ("Reduced %% Moment-Based Models for Oxygen Precipitates and Dislocation Loops in %% Silicon," University of Washington, 2013). %% %% For usage instructions, please read the description of the OModel() function, %% below. Further questions may be directed to me, Bart Trzynadlowski %% (bart.trzy@gmail.com), or Prof. Scott Dunham (dunham@ee.washington.edu). %% %% %% Acknowledgments %% --------------- %% This model was developed at the Nanotechnology Modeling Laboratory %% (University of Washington, Dept. of Electrical Engineering) under the %% supervision of Prof. Scott Dunham. VASP calculations of small oxygen clusters %% and point defects were performed by Renyu Chen. Funding was generously %% provided by SiWEDS member companies (namely Texas Instruments Inc.), Sony %% Corp., and the Semiconductor Research Corp. %% %% %% Version History %% --------------- %% %% 2013.03.20: %% - Initial release included in doctoral dissertation. %% % % [n, t, temp, CO, m0, m1, ncrit, SI, m0DL, m1DL, ncritDL, f, idxK, r] = % OModel(model, r, CO, m0, nAvg, netCI, waferThickness, 'steps', % steps) % OModel(model, r, CO, m0, nAvg, netCI, waferThickness, 'steps', % steps, params) % OModel(model, r, CO, m0, nAvg, netCI, waferThickness, 'points', % temps, times) % OModel(model, r, CO, m0, nAvg, netCI, waferThickness, 'points', % temps, times, params) % % Runs the oxygen model. Outputs are optional; none of them, all of them, or % any number in between may be specified but their order cannot be changed. % Please pay careful attention to the use of units for input parameters, which % is not always consistent. % % The default model parameters were fitted to experiments conducted with FTIR % spectroscopy. All interstitial oxygen concentrations were interpreted using % the "new" (1983) ASTM calibration factor (ASTM F121-83): 4.90 ppma-cm or % 2.45e17 cm^-2. Oxygen concentrations should be normalized to this standard % for best results. % % % Inputs: % % model Selects the model type: 'rkpm' for the reduced model and % 'fkpm' for the full model. If you are unsure, use % 'rkpm'. The full model is not necessarily more accurate % and some initial conditions (m0 and nAvg) cannot be % represented in it, resulting in undefined behavior for % anything other than values near 0 for those quantities. % % r Rediscretization factor. This is called the "sample % discretization factor", S, in the dissertation. Controls % the spacing between samples in size space. A value of 1 % corresponds to unit spacing, which in the full model % generates an impossibly large number of equations but is % valid in the RKPM model. Default physical parameters % were fitted assuming a value of 1.1, which is % recommended. Values in the range of 1.05-1.3 are most % reasonable. % % The value passed will be adjusted to ensure that the % sample point corresponding to size k = 72. Sometimes % this fails in which case a slightly different r should % be tried. % % CO Initial interstitial oxygen concentration [cm^-3]. % % m0 In the RKPM model, the initial concentration of oxygen % precipitates [cm^-3]. Avoid values less than 1. % % nAvg In the RKPM model, the average size of initial oxygen % precipitates [oxygen atoms]. Avoid values less than 73. % % netCI Initial point defect concentrations specified as the net % silicon interstitial concentration: CI - CV [cm^-3]. % This is the only quantity that may be negative. % % waferThickness Wafer thickness [um]. The model is zero-dimensional and % solves for a single spatial sample point that can be % interpreted as either existing in the bulk or as the % average value over the entire wafer. Oxygen out- % diffusion is assumed negligible and not simulated but % the diffusion of point defects to the wafer surface and % recombination there is approximated. If unknown, % reasonable values are 400 or 500 [um]. % % steps When the temperature schedule format is 'steps', a % vector of process steps and ramp rates defines the % schedule. Its format is: % % [temp1, time1] % [temp1, time1, ramp12, temp2, time2] % [temp1, time1, ramp12, temp2, time2, ramp23, temp3, % time3, ...] % % At least a single temperature [C] and time [hr] must be % specified. Subsequent steps can be added in triplets % consisting of the ramp rate [C/min] from the previous % step, temperature, and time. % % temps When the temperature schedule format is 'points', a % sequence of temperatures and corresponding absolute % times define the schedule. The temperature vector is % specified in units of [C]. Temperatures are interpolated % linearly between adjacent points. % % times When the temperature schedule format is 'points', this % defines the absolute time [sec] of each temperature. % The temperature and time vectors must be of equal % lengths. The first time must be 0, subsequent times must % appear in ascending order, and a minimum of two points % are required. Duplicate time values (i.e., zero-length % steps) are not allowed. % % params Allows the physical parameters to be specified as a % vector: % % [Alpha750, Alpha1050, Css750, Css1050] % % Alpha750 is the precipitate surface energy [J/m^2] at % T = 750 C and Alpha1050 is the value at T = 1050 C. The % two points define a linear relationship between the % surface energy and temperature. Similarly, Css750 and % Css1050 are the solubility concentrations [cm^-3] of % oxygen in silicon at T = 750 C and 1050 C, respectively. % They are used to construct an Arrhenius function: % % Css = Css0 * exp(-Ea / kBT) % % where Css0 is the prefactor and Ea is the activation % energy. % % The default parameters are the best-fit values from my % dissertation: % % Alpha750 = 0.1915 [J/m^2] % Alpha1050 = 0.2565 [J/m^2] % Css750 = 4.85551e15 [cm^-3] % Css1050 = 2.2888e17 [cm^-3] % % % Outputs: % % n Sizes of each solution term. Values of 0 at the front of % the vector correspond to non-precipitate solution terms % (e.g., point defects, dislocation loops). CO is size 1. % % t Time points at which solutions are available [sec]. % % temp Temperature profile [C]. % % CO Interstitial oxygen concentration [cm^-3]. % % m0 Concentration of all oxygen precipitates of size 72 or % larger [cm^-3]. The RKPM model solves this directly. The % FKPM model computes this as a post-processing step by % integrating the precipitate size distribution at each % available time step. % % m1 Concentration of oxygen atoms in all precipitates of % size 72 or larger [cm^-3]. Solved directly by the RKPM % model and integrated during post-processing by the FKPM % model. % % ncrit Critical size of oxygen precipitates [oxygen atoms]. % This is a valuable diagnostic aid. Precipitates larger % than this will tend to grow while smaller ones will % shrink. Negative values are equivalent to an infinitely % large critical size (all precipitates dissolve). % % SI Silicon interstitial supersaturation: CI/CI*, where CI* % is the thermal equilibrium concentration of % interstitials at the current temperature. Because fast % recombination is assumed, CV/CV* = 1/SI. % % m0DL Concentration of all dislocation loops [cm^-3]. Faulted % dislocation loops are nucleated at the surface of oxygen % precipitates due to high interstitial supersaturations % and local precipitate-induced strain. % % m1DL Concentration of silicon interstitials in all % dislocation loops [cm^-3]. % % ncritDL Critical size of dislocation loops [interstitial atoms]. % % f Complete solution vector [cm^-3]. Non-zero values of % n(i) are the size, in oxygen atoms, of the precipitate % given by f(:,i). % % idxK Index of size k (72) in the solution vector. f(1,idxK) % would therefore be f(k) (using my dissertation's % notation) at the very first time step. % % r Actual rediscretization factor used in the simulation. % This is not in general the same as the value specified % but is instead adjusted to ensure k = 72. % % % Examples: % % 1. A single-step process, 750 C for 10 hr. % % OModel('rkpm', 1.1, 9e17, 1, 73, -4e12, 400, 'steps', [750 10]) % % 2. As above but specified with 'points'. % % OModel('rkpm', 1.1, 9e17, 1, 73, -4e12, 400, 'points', ... % [750 750], [0 10*3600]) % % 3. A two-step process, 800 C for 4 hr, 1050 C for 16 hr, with a 10 C/min % ramp rate between steps. % % OModel('rkpm', 1.1, 9e17, 1e5, 1e4, -4e12, 400, 'steps', ... % [800 4 10 1050 16]) % % 4. As above but specified with 'points'. % % OModel('rkpm', 1.1, 9e17, 1e5, 1e4, -4e12, 400, 'points', ... % [800 800 1050 1050], ... % [0, 4*3600, 4*3600+1500, 4*3600+1500+16*3600]) % % 5. Obtaining the first 4 output vectors. % % [n, t, temp, CO] = OModel(...); % % 6. Obtaining the first 7 output vectors. % % [n, t, temp, CO, m0, m1, ncrit] = OModel(...); % function varargout = OModel(model, r, CO, m0, nAvg, netCI, waferThickness, format, varargin) % Default parameters (Alpha750, Alpha1050, Css750, Css1050) parameters = [ 0.1915 0.2565 4.85551e15 2.2888e17 ]; % Clear outputs so MATLAB doesn't complain if we abort early varargout = cell(1, nargout); % % Validate model type % global modelType; % 'f' for full model, 'r' for RKPM if (model(1) == 'f') || (model(1) == 'F') modelType = 'f'; elseif (model(1) == 'r') || (model(1) == 'R') modelType = 'r'; else fprintf('Error: Invalid model selection. Use either ''full'' or ''rkpm''.\n'); return; end % % Validate schedule format, fetch parameters, and construct temperature % schedule. The temperature schedule is described by two vectors, % constructed here: % % 1. tempProfile: a vector of temperatures [C] % 2. timeProfile: a vector of corresponding times [sec] to interpolate % temperatures between. % if strcmpi(format, 'steps') % Get thermal steps if nargin < 9 fprintf('Error: Step vector argument is missing.\n'); return; end steps = varargin{1}; if length(steps) < 2 fprintf('Error: Step vector must contain at least two elements: temperature [C], time [hr].\n'); return; end % Get optional model parameters if nargin == 10 parameters = varargin{2}; elseif nargin > 10 fprintf('Error: Too many arguments for ''steps'' format.\n'); return; end % % The step vector specifies the process as a series of temperature % steps consisting of a temperature and time, connected with a ramp % rate. The format is: % % temp1, time1, [ramp1, temp2, time2, ...] % % where [] indicates optional elements. Temperatures are in [C], % times in [hr], and ramp rates in [C/min]. Each additional step % requires a ramp rate, temperature, and time. The ramp rate is % used to raise/lower the temperature linearly over time until the % next temperature is reached. Then, that temperature will be held % constant for the specified amount of time. % if length(steps) < 2 fprintf('Error: The minimum step vector size is 2 elements: temperature [C], time [hr].\n'); return; end % Begin building temperature profile tempProfile = [ steps(1) steps(1) ]; % temperature 1 timeProfile = [ 0 steps(2) ] * 3600; % from t = 0 to time1 % Remaining elements must be triplets of: ramp rate, temp, time num = length(steps) - 2; % elements remaining in vector idx = 3; if mod(num, 3) ~= 0 fprintf('Error: Additional temperature steps require 3 elements: ramp rate [C/min], temperature [C], time [hr].\n'); return; end % Construct the rest of the temperature profile for i = 1:(num/3) rampRate = steps(idx); temp = steps(idx+1); time = steps(idx+2); idx = idx + 3; % Ramp up to next temperature rampTime = abs(60*((temp-tempProfile(end))/rampRate)); % [sec] tempProfile(end+1) = temp; timeProfile(end+1) = timeProfile(end) + rampTime; % Hold for specified time tempProfile(end+1) = temp; timeProfile(end+1) = timeProfile(end) + time*3600; end elseif strcmpi(format, 'points') % % This mode takes a temperature list and a corresponding time list % (in seconds) to define the temperature schedule. % % Get temperatures if nargin < 9 fprintf('Error: Temperature vector argument is missing.\n'); return; end tempProfile = varargin{1}; if length(tempProfile) < 2 fprintf('Error: Temperature vector must specify at least two temperatures.\n'); return; end % Get times if nargin < 10 fprintf('Error: Time vector argument is missing.\n'); return; end timeProfile = varargin{2}; if length(timeProfile) ~= length(tempProfile) fprintf('Error: Temperature and time vectors must have equal lengths.\n'); return; end if timeProfile(1) ~= 0 fprintf('Error: First time must be 0.\n'); return; end % Get optional model parameters if nargin == 11 parameters = varargin{3}; elseif nargin > 11 fprintf('Error: Too many arguments for ''points'' format.\n'); return; end else fprintf('Error: Invalid temperature schedule format. Use either ''steps'' or ''points''.\n'); return; end % % Physical parameters (temperature-independent) % InitConstants(parameters, waferThickness); % % Global temperature callback. Interpolates between steps in the time % profile. % global Temperature; % function callback to compute temperature [K] as a function of time [sec] Temperature = @(t) interp1(timeProfile, tempProfile, t) + 273.15; % % Solution vector layout and sample points % global n idxNetCI idxM0DL idxM1DL idxM0Hi idxM1Hi idxSizeOne idxK; [error, r] = InitSolutionLayout(r); if error == true return; end % % Initialize solutions % UpdateTemperature(0); % set initial temperature fInit = zeros(1,length(n)); % row vector of solutions fInit(idxSizeOne) = CO; fInit((idxSizeOne+1):idxK) = ones(1,idxK-idxSizeOne); % precipitates up to size k initialized to 1 fInit(idxM0Hi) = m0; fInit(idxM1Hi) = nAvg*fInit(idxM0Hi); fInit(idxNetCI) = netCI; % % Run the simulation % PrintSimulationSettings(r, tempProfile, timeProfile, fInit, waferThickness); OutputCallback(0, fInit, 'init'); % initialize progress indicator [t, f] = ode15s(@ODECallback, [0 timeProfile(end)], fInit); OutputCallback(0, squeeze(f(end,:)), 'done'); % % Return all results that were requested % % n: sample point vector if nargout >= 1 varargout{1} = n; end % t: time points vector if nargout >= 2 varargout{2} = t; end % temp: temperature vector [C] if nargout >= 3 varargout{3} = t; for i = 1:length(t) varargout{3}(i) = Temperature(varargout{3}(i)) - 273.15; end end % CO: vector of interstitial oxygen if nargout >= 4 varargout{4} = squeeze(f(:,idxSizeOne)); end % m0: vector of m0 (precipitates w/ n >= k) m0 = []; m1 = []; m2 = []; % unused if nargout >= 5 if modelType == 'f' % Full model, compute all moments at each time step for i = 1:length(t) [m0(end+1), m1(end+1), m2(end+1)] = ComputeMoments(squeeze(f(i,:)), idxK, length(n)); end else % RKPM model, copy moments from solution vector m0 = squeeze(f(:,idxM0Hi)); m1 = squeeze(f(:,idxM1Hi)); end varargout{5} = m0; end % m1: vector of m1 (all oxygen atoms in precipitates w/ n >= k) if nargout >= 6 varargout{6} = m1; end % ncrit: vector of critical sizes of oxygen precipitates SI = []; % also compute SI (temperature-dependent) here ncrit = []; ncritDL = []; if nargout >= 7 global CIstar; for i = 1:length(t) UpdateTemperature(t(i)); % recomputes CI* for this temperature SI(end+1) = ComputeCI(squeeze(f(i,:))) / CIstar; ncrit(end+1) = CriticalSizeEstimated(f(i,idxSizeOne), 1/SI(end)); ncritDL(end+1) = LoopCriticalSize(m1(i)/m0(i), 1/SI(end)); end varargout{7} = ncrit; end % SI: vector of CI/CI* if nargout >= 8 varargout{8} = SI; end % m0DL: vector of dislocation loop m0 if nargout >= 9 varargout{9} = squeeze(f(:,idxM0DL)); end % m1DL: vector of dislocation loop m1 if nargout >= 10 varargout{10} = squeeze(f(:,idxM1DL)); end % ncritDL: vector of critical sizes of dislocation loops if nargout >= 11 varargout{11} = ncritDL; end % f: solution vector (2D -- time, solution) if nargout >= 12 % If RKPM model, fill in the estimator if modelType == 'r' for i = 1:length(t) k = n(idxK); nAvg = max(m1(i)/m0(i),k); f(i,idxK) = FkEstimator(m0(i), nAvg); end end varargout{12} = f; end % idxK: index in solution vector of size k (note: not solved in RKPM model) if nargout >= 13 varargout{13} = idxK; end % r: sample discretization factor (may have been recomputed) if nargout >= 14 varargout{14} = r; end end % % status = OutputCallback(t, f, flag) % % Progress indicator. Called by the solver for each time step to display % information. % function status = OutputCallback(t, f, flag) persistent step; global modelType Temperature n idxNetCI idxM0DL idxM1DL idxM0Hi idxM1Hi idxSizeOne idxK Css CIstar CVstar; if isempty(flag) % Only print every 100 steps step = step + 1; if mod(step, 100) == 0 tempC = Temperature(t) - 273.15; CO = f(idxSizeOne); SI = ComputeCI(f)/CIstar; if modelType == 'f' % compute moments for full model (time consuming!) [m0, m1, m2] = ComputeMoments(f, idxK, length(n)); else % RKPM model solves moments m0 = f(idxM0Hi); m1 = f(idxM1Hi); end nAvg = m1/m0; ncrit = CriticalSizeEstimated(CO, 1/SI); m0DL = f(idxM0DL); m1DL = f(idxM1DL); fprintf('Step: %d T=%g C t=%1.2f sec (%1.1f hr) CO=%1.2e (CO/Css=%1.1f) CI/CI*=%1.1f m0=%1.2e m1=%1.2e nAvg=%1.2e (ncrit=%1.2e) m0DL=%1.2e m1DL=%1.2e\n', step, tempC, t, t/3600, CO, CO/Css, SI, m0, m1, nAvg, ncrit, m0DL, m1DL); end elseif flag == 'init' step = 0; fprintf('Solver Progress\n'); fprintf('---------------\n'); else % flag == 'done' fprintf('\nResults\n'); fprintf('-------\n'); CO = f(idxSizeOne); CI = ComputeCI(f); CV = ComputeCV(f); m0 = f(idxM0Hi); m1 = f(idxM1Hi); nAvg = m1 / m0; m0DL = f(idxM0DL); m1DL = f(idxM1DL); nAvgDL = m1DL / m0DL; fprintf('CO = %g cm^-3 (CO/Css = %1.2f)\n', CO, CO/Css); fprintf('m0 = %g cm^-3\n', m0); fprintf('m1 = %g cm^-3\n', m1); fprintf('nAvg = %g (Avg. radius: %1.1f nm)\n', nAvg, Radius(nAvg)/1e-7); fprintf('CI = %g cm^-3 (CI/CI* = %1.2f)\n', CI, CI/CIstar); fprintf('CV = %g cm^-3 (CV/CV* = %1.2f)\n', CV, CV/CVstar); fprintf('m0DL = %g cm^-3\n', m0DL); fprintf('m1DL = %g cm^-3\n', m1DL); fprintf('nAvgDL = %g (Avg. radius: %1.1f nm)\n', nAvgDL, LoopRadius(nAvgDL)/1e-7); fprintf('\nFinished.\n\n'); end status = 0; end % % PrintSimulationSettings(): % % Prints a descriptive summary of the simulation settings. % function PrintSimulationSettings(r, tempProfile, timeProfile, fInit, waferThickness) global kB n idxSizeOne idxM0Hi idxM1Hi CIstar CVstar Css Css0 CssEa Alpha750 Alpha1050 modelType; % Basic settings fprintf('Simulation Settings\n'); fprintf('-------------------\n'); fprintf('Model: '); if modelType == 'f' fprintf('Full\n'); else fprintf('RKPM\n'); end fprintf('Sampling Factor: %g (%d equations)\n', r, length(n)-4); fprintf('Wafer Thickness: %g um\n', waferThickness); % Model parameters fprintf('Surface Energy: %g J/m^2 (T = 750 C), %g (T = 1050 C)\n', Alpha750, Alpha1050); fprintf('Css: %g cm^-3 (T = 750 C), %g (T = 1050 C)\n', Css0*exp(-CssEa/(kB*(750+273.15))), Css0*exp(-CssEa/(kB*(1050+273.15)))); % Temperature profile fprintf('Temperature Schedule: \n'); for i = 2:length(tempProfile) deltaTime = timeProfile(i) - timeProfile(i-1); temp1 = tempProfile(i-1); temp2 = tempProfile(i); if temp1 == temp2 UpdateTemperature(timeProfile(i)); fprintf('\t%g C (%1.2f hr)\n', temp1, deltaTime/3600); fprintf('\t\tCss = %g cm^-3\n', Css); fprintf('\t\tCI* = %g cm^-3\n', CIstar); fprintf('\t\tCV* = %g cm^-3\n', CVstar); else fprintf('\t%g -> %g C (%g C/min)\n', temp1, temp2, abs(temp2-temp1)/(deltaTime/60)); end end % Initial conditions UpdateTemperature(0); % reset to initial temperature CO = fInit(idxSizeOne); CI = ComputeCI(fInit); CV = ComputeCV(fInit); m0 = fInit(idxM0Hi); m1 = fInit(idxM1Hi); nAvg = m1 / m0; fprintf('CO = %g cm^-3 (CO/Css = %1.2f)\n', CO, CO/Css); fprintf('m0 = %g cm^-3\n', m0); fprintf('m1 = %g cm^-3\n', m1); fprintf('nAvg = %g (Avg. radius: %1.1f nm)\n', nAvg, Radius(nAvg)/1e-7); fprintf('CI = %g cm^-3 (CI/CI* = %1.2f)\n', CI, CI/CIstar); fprintf('CV = %g cm^-3 (CV/CV* = %1.2f)\n', CV, CV/CVstar); fprintf('\n'); end % % InitSolutionLayout(r): % % Given the sample discretization factor, r, defines all components of the % solution vector, generates the sample points (n), and defines size k. % % Sets error to true if an unrecoverable error occurred and the simulation % should be aborted. Returns the new, adjusted discretization factor in newR. % function [error, newR] = InitSolutionLayout(r) % Pointers: indices of solutions in the overall solution vector, f global idxNetCI idxM0DL idxM1DL; % point defects and dislocation loops global idxM0Hi idxM1Hi; % oxygen precipitate moments (n >= k) global idxSizeOne; % interstitial oxygen global idxChangeOver; % sample spacing becomes > 1 after this global idxK; % size k precipitate (not solved in RKPM) % Samples global n; % size n(i) of each solution vector element i global maxSmallClusterSize; % small cluster/macroscopic transition (between 3 and nChangeOver) % Imported... global modelType; newR = r; % % Size distribution and layout of solution variables is described below. % % Index Description Pointer % 1 NetCI (silicon interstitials minus vacancies) idxNetCI % 2 m0DL (m0 for dislocation loops) idxM0DL % 3 m1DL (m1 for dislocation loops) idxM1DL % 4 -- % 5 -- % 6 -- % 7 -- % 8 m1 (DFA approximation of m1 in RKPM case) idxM1Hi % 9 m0 (0th moment) idxM0Hi % 10 CO (oxygen interstitials; i.e., the solute) idxSizeOne % 11 CO2 (O2, all point defect states in relative equilibrium) % 12 f3 (size 3 precipitate) % ... % % The small cluster/macroscopic transition point, maxSmallClusterSize, % indicates where macroscopic energies begin. The transition from % maxSmallClusterSize-1 -> maxSmallClusterSize uses small cluster energies % for both, maxClusterSize -> maxClusterSize + 1 uses macro energies for % both. % % The change-over from discrete (single spacing) to interpolated equations, % at index idxChangeOver, must be between size maxSmallClusterSize and index % idxK. % % For the moment-based model, the maximum size is k, beyond which a moment- % based approximation is used. For the full model, solutions are tracked to % a much larger size. % idxNetCI = 1; idxM0DL = 2; idxM1DL = 3; idxM0Hi = 9; idxM1Hi = 8; idxSizeOne = 10; % % Set up boundaries of different regions (set these carefully) % nChangeOver = 10; % size beyond which sample spacing can be > 1 maxSmallClusterSize = 3; % sizes 2 to here are small clusters k = 72; % size at which moments begin in moment model if modelType == 'f' maxSampleSize = 1e9; % maximum sample size for full model else maxSampleSize = k; % in RKPM mode, moment model begins at k (the sample at k won't actually be solved) end if (nChangeOver <= maxSmallClusterSize) || (nChangeOver >= k) fprintf('Internal error: Discrete change-over size is invalid. Please change it.\n'); error = true; return; end % % Generate samples with logarithmic spacing such that a sample at size k % exists. % idxChangeOver = idxSizeOne-1+nChangeOver; if r == 1.0 idxK = k-nChangeOver; else idxK = log((k-nChangeOver)*(r-1)+1)/log(r); % index at which size k currently exists (not likely to be integral) end idxK = floor(idxK); % make it integral newR = fzero(@(x) x^idxK - (k-nChangeOver)*(x - 1) - 1, r); % compute ratio necessary for n(idxK) = k idxK = idxK + idxChangeOver; % Leading zeros are for non-oxygen solutions that have no 'size' n = [ zeros(1,idxSizeOne-1) Samples(maxSampleSize, nChangeOver, newR) ]; % Sometimes the math above fails... if abs(n(idxK)/k-1) > 0.0001 fprintf('Internal error: Unable to generate sample point for k = %d. Please choose a different rediscretization factor, r.\n', k); error = true; return; end % Success! error = false; end % % n = Samples(maxSize, changeOverSize, ratio) % % Generate a size space sample vector, n(i). Sample spacing is 1 up to the % change-over size and then increases by a factor of "ratio" for each sample % point. The spacing between the change-over size and the following sample is % guaranteed to be 1. % % I didn't know about MATLAB's logspace() when I wrote this :) This could be % made much simpler... % function n = Samples(maxSize, changeOverSize, ratio) % Spacing is 1 for [1:changeOverSize] n = 1:1:changeOverSize; % Logarithmic samples delta = 1; i = changeOverSize + 1; while n(i-1) < maxSize n(i) = n(i-1) + delta; delta = delta * ratio; i = i + 1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Parameters and Constants %% %% Physical constants and parameters are defined here. A global callback, %% Temperature(t), which returns the temperature in K at time t, must be %% defined. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % InitConstants(params, waferThickness) % % Initializes all temperature-independent global physical constants and % parameters. Wafer thickness is supplied in microns. % function InitConstants(params, waferThickness) global kB CSi VolSiO2 VolSi uSi KSi KSiO2; global reactDist Alpha750 Alpha1050 Css0 CssEa; global HfOVRelaxed HfO2Relaxed HfO2VRelaxed HfV; global EbVV EbOV EbOO EbOOnVm; global aSi StrainOi StrainOn StrainV; global b coreRadius HfSF HfCore KDL indEpsSF; global surfReactDist; % Conversion factor for [Pa] -> [eV/cm^3] PaToLocal = 6.24151e12; % 1 [Pa] = 6.24150974e12 [eV/cm^3] % Physical constants kB = 8.62e-5; % Boltzmann constant [eV/K] CSi = 5e22; % silicon lattice site density [cm^-3] VolSiO2 = 4.35e-23; % volume of SiO2 molecule [cm^3] (Ref: Hoessinger, http://www.iue.tuwien.ac.at/phd/hoessinger/node28.html) VolSi = 2.00e-23; % volume of Si atom [cm^3] uSi = 6.49e10*PaToLocal;% shear modulus of Si [Pa] -> [eV/cm^3] (Ref: Hull, Properties of Crystalline Silicon) KSi = 9.78e10*PaToLocal;% bulk modulus of Si [Pa] -> [eV/cm^3] (Ref: Hull) KSiO2 = 3.69e10*PaToLocal;% bulk modulus of SiO2 [Pa] -> [eV/cm^3] (Ref: Europhys. Lett., 57 (3), pp. 375-381 2002) % Precipitates reactDist = 5e-8; % interface reaction distance [cm] Alpha750 = params(1); % surface energy at 750C [J/m^2] Alpha1050 = params(2); % surface energy at 1050C [J/m^2] Css750 = params(3); % solid solubility at 750C [cm^-3] Css1050 = params(4); % solid solubility at 1050C [cm^-3] T750 = 750 + 273.15; T1050 = 1050 + 273.15; CssEa = log(Css750/Css1050)*kB/(1/T1050-1/T750); % solubility activation energy [eV] Css0 = Css750*exp(log(Css750/Css1050)/(T750/T1050-1)); % solubility prefactor [cm^-3] % Small cluster formation enthalpies (relative to perfect Si and %interstitial O) without strain (fully relaxed, no strain) HfOVRelaxed = 2.02; % OV [eV] HfO2Relaxed = -0.39; % O2 [eV] HfO2VRelaxed = 0.51; % O2V [eV] % Vacancy formation enthalpy HfV = 3.5; % [eV] % Binding energies, estimated from ab initio formation energies EbVV = -1.7; % [eV], literature estimates are from 1.5-2.0 EbOV = HfOVRelaxed-HfV; % [eV] EbOO = HfO2Relaxed; % [eV] EbOOnVm = EbOO; % assume same % Estimates of small cluster linear (transformational, in radial direction) % strain components. From ab initio and valid only for a volume of % (2*aSi)^3. aSi = 5.431e-8; % silicon lattice constant [cm] StrainOi = 0.002325; % first Oi StrainOn = 0.001146; % each additional Oi StrainV = -0.004344; % each V % Dislocation loops b = sqrt(3)/3; % Burgers vector magnitude [aSi] coreRadius = b; % dislocation core radius [aSi] HfSF = 0.0152; % stacking fault energy [eV/atom] HfCore = 0; % core energy [eV/aSi] KDL = 72/sqrt(2); % dislocation energy prefactor (sqrt(2) was for semi-circular geometry) [GPa] indEpsSF = 0.996; % normalized induced strain of a stacking fault % 0D surface boundary condition distance factor surfReactDist = ((waferThickness * 1e-4)^2)/12; end % % UpdateTemperature(t): % % This should be called at each time step of the solver. Computes the % temperature at time t and recalculates all global temperature-dependent % parameters. % function UpdateTemperature(t) global kBT DO Css Alpha HAtomic DI DV CIstar CVstar; % exported global Temperature kB CSi Css0 CssEa Alpha750 Alpha1050; % imported % Get current temperature T = Temperature(t); % Physical constants kBT = kB*T; % Oxygen parameters DO = 0.13*exp(-2.53/kBT); % oxygen diffusivity [cm^2/sec] Css = Css0 * exp(-CssEa/kBT); % oxygen solid solubility [cm^-3] HAtomic = kBT*log(Css/CSi); % per-atom formation energy (i.e., Gp) [eV] Alpha = 6.24151e14 * interp1([750 1050]+273.15, [Alpha750 Alpha1050], T, 'linear', 'extrap'); % surface energy [eV/cm^2] % Silicon interstitials and vacancies DI = 51.4*exp(-1.77/kBT); % interstitial diffusivity [cm^2/sec] DV = 3.07*exp(-2.12/kBT); % vacancy diffusivity [cm^2/sec] CIstar = (2980/51.4)*5e22*exp((-4.95+1.77)/kBT); % interstitial equilibrium [cm^-3] CVstar = (86/3.07)*5e22*exp((-4.56+2.12)/kBT); % vacancy equilibrium [cm^-3] end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Dislocation Model %% %% Functions for the dislocation model: energetics, etc. The formation energy %% is repeatedly used here in different forms. Care must be taken that all of %% these forms are consistent with each other. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = LoopStrainEnergy(stress) % % Given a stress in [Pa], returns the stress energy per atom [eV] based on the % induced strain of the stacking fault. % function z = LoopStrainEnergy(stress) global VolSi indEpsSF; u = 1.602e-19; % unit conversion factor [J/eV] v = VolSi * 1e-6; % atomic volume [cm^3] -> [m^3] because strain is [Pa] = [N/m^2] z = -v*indEpsSF*stress; % strain energy per SF atom [J] z = z / u; % [eV] end % % z = PrecipitateStress(n, SV) % % Computes the tangential stress (sigma_theta,theta; sigma_phi,phi [Pa]) at the % very edge of an oxygen precipitate, in the silicon matrix. Equation B.7 from % the Senkader thesis is used: S. Senkader, "Physical Modeling and Simulation of % Oxygen Precipitation in Silicon", Vienna University of Technology, 1996. The % sign convention used by Senkader seems to be opposite of what is expected, so % the equation is inverted here. % function z = PrecipitateStress(n, SV) global uSi KSiO2; mOpt = MOptimal(n,SV); u = uSi / 6.24151e12; % [Pa = N/m^2] K = KSiO2 / 6.24151e12; % [Pa = N/m^2] A = -3*K*eT(n,mOpt)/(3*K+4*u); z = -2*u*A; % stress at r = rp (w/ inverted sign) end % % ncrit = LoopCriticalSize(nOx, SI) % % Dislocation loop critical size as a function of CI/CI* near an oxygen % precipitate of size nOx (e.g., m1/m0). Dislocation loop formation energy takes % the form: % % GtotDL(n) = A*n + B*sqrt(n) + C*sqrt(n)*[log(D*sqrt(n))-1] % % Differentiating with respect to n and setting to 0 allows ncrit to be solved: % % 0 = A + 0.5*B/sqrt(n) + 0.5*C/sqrt(n)*[log(D*sqrt(n))-1] % 0 = A*sqrt(n) + 0.5*B + 0.5*C*[log(D*sqrt(n))-1] % % Setting x = sqrt(n), we get: % % x = -(0.5*B + 0.5*C*[log(D*x)-1])/A % % Iterating this a few times provides a solution. % function ncrit = LoopCriticalSize(nOx, SI) global kB kBT HfSF HfCore KDL b coreRadius; SV = 1/SI; rf = sqrt(sqrt(3)/(8*pi)); % radius factor [aSi/sqrt(atoms)] u = .9998; % unit conversion factor0 [GPa*aSi^3] -> [eV] % Energy equation: Gtot(n) = A*n + B*sqrt(n) + C*sqrt(n)*[log(D*sqrt(n))-1] A = -kBT*log(SI) + HfSF + LoopStrainEnergy(PrecipitateStress(nOx, SV)); B = 2*pi*rf*HfCore; C = u*0.5*KDL*b*b*rf; D = rf*(8/coreRadius); % Iterate x = sqrt(1e4); % initial guess for sqrt(n) for i = 1:6 x = -0.5*(B+C*(log(D*x)))/A; end ncrit = x*x; end % % z = LoopCnStar(n, nOx, SV) % % Equilibrium concentration of interstitials with a size n dislocation loop near % an oxygen precipitate of size nOx. % function z = LoopCnStar(n, nOx, SV) global kB kBT CIstar HfSF HfCore KDL b coreRadius; % dGtotDL/dn = A + 0.5*B/sqrt(n) + 0.5*C/sqrt(n)*[log(D*sqrt(n))-1] % Here we compute the enthalpy portion of this (log(SI) term not included). rf = sqrt(sqrt(3)/(8*pi)); % radius factor [aSi/sqrt(atoms)] u = .9998; % unit conversion factor0 [GPa*aSi^3] -> [eV] A = HfSF + LoopStrainEnergy(PrecipitateStress(nOx, SV)); B = 2*pi*rf*HfCore; C = u*0.5*KDL*b*b*rf; D = rf*(8/coreRadius); dHdn = A + 0.5*B/sqrt(n) + 0.5*C/sqrt(n)*(log(D*sqrt(n))-1); % Cn* = CI* * exp{dHtotDL(n)/dn) / kBT} z = CIstar * exp(dHdn/kBT); end % % z = LoopFnStar(n, m0DL, m0Ox, nAvgOx, SI) % % Equilibrium concentration of size n dislocation loops assuming % heterogeneous nucleation at oxygen precipitate sites. % function z = LoopFnStar(n, m0DL, m0Ox, nAvgOx, SI) global kB kBT b HfSF HfCore KDL coreRadius reactDist CSi; if n == 0 z = 0; return; end SV = 1/SI; % Compute energy terms u = .9998; % unit conversion factor [GPa*aSi^3] -> [eV] r = sqrt(n*sqrt(3)/(8*pi)); % radius [aSi] Gperim = 2*pi*r*HfCore; Gself = u*0.5*KDL*b*b*r*(log(8*r/coreRadius)-1); % elastic self-energy G = -n*kBT*log(SI) + n*(HfSF + LoopStrainEnergy(PrecipitateStress(nAvgOx, SV))) + Gperim + Gself; % fn*, the 4*2*pi*r/a0 factor is because there are 4 <111> planes z = (8*pi*Radius(n)/reactDist)*(m0Ox-m0DL)*exp(-G/kBT); end % % z = LoopRadius(n) % % Dislocation loop radius [cm]. % function z = LoopRadius(n) global aSi; z = aSi * sqrt(n*sqrt(3)/(8*pi)); end % % z = LoopLambda(n) % % Kinetic growth factor [cm] for dislocation loops. % function z = LoopLambda(n) global aSi; r = sqrt(n*sqrt(3)/(8*pi)); c = sqrt(3)/4; z = 4*pi*pi*aSi*r/log(8*r/c); end % % z = LoopGRate(n, CI) % % Dislocation loop growth rate. % function z = LoopGRate(n, CI) global DI; if n == 0 z = 0; else z = DI*LoopLambda(n)*CI; end end % % z = LoopDRate(n, nOx, SV) % % Dislocation loop dissolution rate near a size nOx oxygen precipitate. % function z = LoopDRate(n, nOx, SV) global DI; if n == 0; z = 0; else z = DI*LoopLambda(n)*LoopCnStar(n, nOx, SV); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Energies of Small Oxygen Clusters %% %% Energies and optimal number of point defects for discrete, small oxygen %% clusters, including effective quantities based only on the number of oxygen %% atoms, n, assuming all possible point defect states, m, are in relative %% equilibrium. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = HStrainDiscrete(n,m) % % Computes an estimate of the strain energy for discrete oxygen clusters based % on ab initio calculations. Only valid for sizes n << 64. % function z = HStrainDiscrete(n,m) global uSi KSiO2 aSi StrainOi StrainOn StrainV; % eT(n,m): misfit strain, computed from ab initio estimates eT = StrainOi + (n-1)*StrainOn + m*StrainV; eC = eT / (1+4*uSi/(3*KSiO2)); % Strain energy. Volume is based on 2x2x2 cell used in VASP for strain % calculations: (8*aSi)^3. Because of this fixed volume, these discrete % strain results become invalid as the precipitate size approaches the % volume of a 64 silicon atom cell. z = ((8*aSi)^3)*6*uSi*eT*eC; end % % z = HfDiscrete(n,m) % % Computes the formation enthalpy of a small n, m cluster. A very simple % heuristic based on limited ab initio data is used. % % OV is included here, unlike in the solved model. Not intended for use at % run-time but can be used to pre-compute discrete cluster values for sizes % 3 and above. % function z = HfDiscrete(n,m) global HfOVRelaxed HfO2Relaxed HfO2VRelaxed HfV EbVV EbOO EbOV EbOOnVm; % Strain energy term StrainEnergy = HStrainDiscrete(n,m) - n*HStrainDiscrete(1,0); % Handle cases for which we have ab initio data if (n == 1) if (m == 0) z = 0; % Oi is 0 by definition return; elseif (m == 1) z = HfOVRelaxed+StrainEnergy; return; end elseif (n == 2) if (m == 0) z = HfO2Relaxed+StrainEnergy; return; elseif (m == 1) z = HfO2VRelaxed+StrainEnergy; return; end end % % Compute energy by first binding Oi's to dangling bonds caused by V's. % Then, any remaining Oi should be paired to the precipitate. % % Formation energy of all vacancies Hf = m*HfV; % Reduce energy by number of paired V's freeO = n; numVV = floor(m/2); numDanglingBonds = 4*m-2*numVV; Hf = Hf + EbVV*numVV; % Each Oi can satisfy two dangling bonds if (freeO >= (numDanglingBonds/2)) Hf = Hf + EbOV*numDanglingBonds/2; freeO = freeO - numDanglingBonds/2; else Hf = Hf + EbOV*freeO; freeO = 0; end % Remaining Oi bind to each other and then precipitate if (freeO >= 2) Hf = Hf + EbOO; freeO = freeO - 2; end Hf = Hf + EbOOnVm*freeO; % Add in discrete strain energy, less n*HStrain(Oi) (which is our reference) z = Hf + StrainEnergy; end % % z = DeltaGfvSmallCluster(n, SV) % % Change in energy of small cluster formation due to enthalpy of formation and % vacancy entropy of mixing terms. All possible m states (i.e., number of point % defects incorporated) are assumed to be in relative thermal equilibrium: % % fn* = f*n,0 + f*n,1 + ... + f*n,n % % From this, we can compute a single DeltaGTotal for a size n small cluster: % % DeltaGTotal(n) = -n*kBT*log(CO/CSi) % - kBT*log(exp{-Hf_n,0/kBT} + SV*exp{-Hf_n,1/kBT} + ...) % % And finally, the definition of DeltaGfv is just: % % DeltaGfv = -kBT*log(exp{-Hf_n,0/kBT} + SV*exp{-Hf_n,1/kBT} + ...) % % This term is defined because it is not convenient to break down the effective, % aggregate energy for small clusters of size n any further. It is this term % that must be used to enforce continuity between small clusters and % macroscopic precipitates. % function z = DeltaGfvSmallCluster(n, SV) global kBT; % Compute term inside log first z = 0; for m = 0:n z = z + (SV^m)*exp(-HfDiscrete(n,m)/kBT); end % Return the final expression z = -kBT*log(z); end % % z = DeltaGTotalSmallCluster(n, CO, SV) % % Total change in formation energy upon small cluster formation. This energy is % continuous with the small cluster form at the transition point. % function z = DeltaGTotalSmallCluster(n, CO, SV) global kBT CSi; z = -n*kBT*log(CO/CSi) + DeltaGfvSmallCluster(n,SV); end % % UpdateSmallClusterEnergies(SV) % % As an optimization, computes DeltaGfvSmallCluster/kBT for all sizes up to % the transition point, storing them in a global array. % % This must be called each time step before computing any quantities that depend % on precipitate energies (including dissolution rates)! % function UpdateSmallClusterEnergies(SV) global maxSmallClusterSize kBT; % imported global DeltaGfvSmallClusterNoKBT; % exported % Compute unitless formation energy factors DeltaGfvSmallClusterNoKBT = zeros(maxSmallClusterSize,1); for n = 1:maxSmallClusterSize DeltaGfvSmallClusterNoKBT(n) = DeltaGfvSmallCluster(n, SV)/kBT; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Energies of Macroscopic Oxygen Precipitates %% %% Energies and optimal number of point defects for macroscopic (large size %% limit) precipitates. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = eT(n, m) % % Translational strain (linear strain along radial direction) of an oxygen % precipitate. % function z = eT(n, m) global VolSiO2 VolSi; z = (VolSiO2*n/(2*VolSi*(m+n/2)))^(1/3) - 1; end % % z = MOptimal(n, SV) % % Optimal number of point defects (vacancies/interstitials) to eject. That % is, the number that minimizes the free energy of a size n precipitate. Full % derivation appears in dissertation. SV is the vacancy supersaturation. % function z = MOptimal(n, SV) global kBT VolSiO2 VolSi uSi KSiO2; % Value of m for 0 strain (point about which eT is linearized) m0 = n*((VolSiO2/(2*VolSi)) - 1/2); % Optimal value of m based on linear approximation of eT(n,m) z = 0.5 * n * ((3/4)*kBT*log(SV)*(VolSiO2/(VolSi*VolSi*uSi))*(1+4*uSi/(3*KSiO2)) + (VolSiO2/VolSi) - 1); end % % z = HStrain(n,m) % % Precipitate strain energy. % % NOTE: Minimum strain energy, HStrain(n,mOpt), should probably be computed in a % more optimal fashion than calling this directly. % function z = HStrain(n,m) global VolSi VolSiO2 uSi KSiO2; % eT(n,m): misfit strain eT = (VolSiO2*n/(2*VolSi*(m+n/2)))^(1/3) - 1; % eC eC = eT / (1+4*uSi/(3*KSiO2)); % Strain energy z = (6*4*pi/3)*(Radius(n)^3)*uSi*eT*eC; end % % z = GSurface(n) % % Surface energy component of a precipitate (temperature-dependent). % function z = GSurface(n) global Alpha; z = 4*pi*Alpha*Radius(n)^2; end % % z = GfMacro(n, SV) % % Returns the formation energy of a size n precipitate using the macroscopic % form of the equations. This is the formation energy with entropy of mixing % terms for oxygen interstitials and vacancies excluded. % % Gf(n) = n*HAtomic + GSurface(n) + HStrain(n,mOpt(n)) % % Note that this term alone is not continuous at the transition point between % macroscopic precipitates and small clusters. % function z = GfMacro(n,SV) global HAtomic; if (n <= 1) z = 0; % by definition else mOpt = MOptimal(n,SV); z = n*HAtomic+GSurface(n)+HStrain(n,mOpt); end end % % z = DeltaGfvMacro(n, SV) % % Change in energy upon precipitate formation due to formation energy and % vacancy entropy of mixing terms. Continuity at the transition size between % small clusters and macroscopic precipitates is enforced by applying an % offset. % function z = DeltaGfvMacro(n,SV) global kBT; DeltaGfvMacro3 = -MOptimal(3,SV)*kBT*log(SV) + GfMacro(3,SV); DeltaGOffset = DeltaGfvSmallCluster(3,SV) - DeltaGfvMacro3; z = -MOptimal(n,SV)*kBT*log(SV) + GfMacro(n,SV) + DeltaGOffset; end % % z = DeltaGTotalMacro(n, CO, SV) % % Total change in formation energy upon macroscopic precipitate formation. This % energy is continuous with the small cluster form at the transition point. % function z = DeltaGTotalMacro(n, CO, SV) global kBT CSi; z = -n*kBT*log(CO/CSi) + DeltaGfvMacro(n,SV); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Growth and Dissolution Rates for Oxygen Precipitates %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = Radius(n) % % Radius of a large ("macroscopic") oxygen precipitate in [cm]. % function z = Radius(n) global VolSiO2; z = ( 3 * n * ((VolSiO2/2) / (4*pi)) )^(1/3); end % % z = RatePrefactor(n) % % Common prefactor for growth/dissolution rates. % function z = RatePrefactor(n) global reactDist DO; r = Radius(n); z = 4 * pi * (r^2) * (DO/(reactDist+r)); end % % z = CnStar(n, SV) % % Concentration of oxygen interstitials at equilibrium with a size n macroscopic % precipitate (near the interface). % % Cn* is obtained by differentiating the total formation energy change with % respect to n, setting it to 0, and solving for CO = Cn*: % % Cn* = Css * (CI/CI*)^gammaI * exp{(Gf(n)-Gf(n-1))/kBT} % % Where gammaI = mOpt(n) - mOpt(n-1). % % Note that this is intended to be used with DRate(n), hence it looks backwards % when performing the differentiation. When using this in the growth rate, we % must be consistent, so passing in n+1 is a good idea. % function z = CnStar(n, SV) global kBT Css; mOpt_n = MOptimal(n,SV); mOpt_p = MOptimal(n-1,SV); deltaM = mOpt_n-mOpt_p; deltaGfExc = GSurface(n)+HStrain(n,mOpt_n)-GSurface(n-1)-HStrain(n-1,mOpt_p); SI = 1/SV; z = Css*(SI^deltaM)*exp(deltaGfExc/kBT); end % % z = GRate(n, CO) % % Growth rate from size n to n+1. % function z = GRate(n, CO) z = RatePrefactor(n)*CO; end % % z = DRate(n, CO, SV) % % Dissolution rate of a precipitate from size n to n-1. SV = CV/CV*. % % NOTE: Very important to keep this code in sync with the solver callback's % code for small precipitates. % function z = DRate(n, SV) global maxSmallClusterSize DeltaGfvSmallClusterNoKBT CSi; % % There are a few different size regimes to handle (and the transitions % between them). Sizes 1 and 2 are special and handled directly in the % solver callback. % if (n < 3) z = 0; % error, DRate() cannot be used with n < 3 % % Dissolution of Small Cluster -> Small Cluster % --------------------------------------------- % % Dissolution of size 3 -> size 2: % % d(3) = g(2)*(f2*)/(f3*) = g(2)*exp{(Gtot_3-Gtot_2)/kBT} % % Gtot_n,m = -n*kBT*log(CO/CSi) - m*kT*log(CV/CV*) + Hf_n,m % fn,m* = CSi*exp(-Gtot_n,m/kBT) = CSi*(CO/CSi)^n*(CV/CV*)^m * % exp{-Hf_n,m/kBT} % fn* = fn,0*+fn,1*+...+fn,mMax* = CSi*exp(-Gtot_n/kBT) % % The overall energy term for size n, Gtot_n, can be computed from the % above relations. % % Gtot_n = -kT*log((fn,0*+...+fn,mMax*)/CSi) % = -kT*log((CO/CSi)^n*((CV/CV*)^0*exp{-Hf_n,0/kBT}+ % (CV/CV*)^1*exp{-Hf_n,1/kBT}+ ... % (CV/CV*)^mMax*exp{-Hf_n,mMax/kBT})) % % Size 2 only has two m-states and size 3 has three, resulting in: % % d(3) = g(2)*(CSi/CO)^3*(CO/CSi)^2 * % exp{-log(exp{-Hf_3,0/kBT}+(CV/CV*)*exp{-Hf_3,1/kBT}+ % (CV/CV*)^2*exp{-Hf_3,2/kBT}) % - -log(exp{-Hf_2,0/kBT}+(CV/CV*)*exp{-Hf_2,1/kBT})} % % Easier and equivalent way is to factor out CO/CSi terms leaving the % DeltaGfv terms (vacancy entropy of mixing and formation energy) and take % the difference of those. % elseif (n <= maxSmallClusterSize) % Note: this deltaG is unitless (kBT is factored out and cancels with % the denominator of the exponential, which is why it is omitted). deltaGf = DeltaGfvSmallClusterNoKBT(n) - DeltaGfvSmallClusterNoKBT(n-1); z = RatePrefactor(n-1)*CSi*exp(deltaGf); % % Dissolution of Macroscopic Precipitates % --------------------------------------- % % Dissolution of a size n (macro) -> n-1 (macro): % % d(n) = g(n-1)*(fn*)/(f(n-1)*) = g(n-1)*exp{(Gtot_n-Gtot_n-1)/kBT} % % Where % % Gtot_n = -n*kT*log(CO/CSi) - mOpt(n)*kT*log(CV/CV*) + Gf(n) % Gf(n) = n*Gp + Gsurf(n) + Hstrain(n,mOpt(n)) % % Resulting in: % % d(n) = g(n-1)*(CSi/CO)*(CV*/CV)^(DeltaM)*exp{(Gf(n)-Gf(n-1))/kBT} % DeltaM = mOpt(n)-mOpt(n-1) % % Note that CO will cancel with CO in g(n-1). % else z = RatePrefactor(n-1) * CnStar(n, SV); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Diagnostics %% %% Functions that compute quantities useful for diagnostic and model development %% purposes. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = FnStar(n, CO, SV) % % Equilibrium concentration of size n precipitate, for debugging purposes. % function z = FnStar(n, CO, SV) global maxSmallClusterSize kBT CSi; if (n <= maxSmallClusterSize) Gtot = DeltaGTotalSmallCluster(n, CO, SV); else Gtot = DeltaGTotalMacro(n, CO, SV); end z = CSi*exp(-Gtot/kBT); end % % z = CriticalSizeEstimated % % This is a new expression derived on 14 Nov 2012. It relies on the fact % that the minimum strain, eT(n,mopt), is independent of n. % % A negative result means precipitates cannot form. % function z = CriticalSizeEstimated(CO, SV) global Alpha KSiO2 VolSi VolSiO2 uSi kBT Css; R = 3*(VolSiO2/2)/(4*pi); X = (3/2) * kBT * log(SV) * (VolSiO2/(4*VolSi*VolSi*uSi)) * (1+(4*uSi)/(3*KSiO2)) + (VolSiO2/(2*VolSi)) - (1/2); eTmin = (VolSiO2/(2*VolSi*((1/2)+X)))^(1/3) - 1; HsPrime = 8*pi*R*uSi*eTmin*eTmin/(1+(4*uSi)/(3*KSiO2)); z = (4*pi*Alpha*(2/3)*(R^(2/3)) / (kBT*log((CO/Css)*(SV^X))-HsPrime))^3; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Differential Equations %% %% The complete set of kinetic rate equations that comprise the core of the %% model are defined here. The ODE solver callback evaluates the right-hand side %% of each equation at a particular time step. %% %% NOTE: Because of small oxygen clusters, which are a special case, the ODE %% callback and supporting function must be kept in careful sync with the code %% which computes dissolution rates (based on precipitate energies). %% %% O2 and O2V are both assumed to be present in f2. O3, O3V, and O3V2 are %% assumed for f3. f1 only includes Oi (i.e., CO), not OV. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % z = ComputeCI(f) % % CI and CV are assumed to be in equilibrium: CICV=CI*CV*. Only interstitial % emission is tracked for oxygen precipitation, so it makes sense to only % consider "net CI". This results in two simultaneous equations, allowing CI % and CV to be computed on demand: % % CI*CV = (CI*)*(CV*) % NetCI = CI - CV % function z = ComputeCI(f) global idxNetCI CIstar CVstar; NetCI = f(idxNetCI); z = 0.5*(NetCI+sqrt(NetCI^2+4*CIstar*CVstar)); end % % z = ComputeCV(f) % % Computes CV based on the current NetCI, CI*, and CV*. % function z = ComputeCV(f) global idxNetCI; NetCI = f(idxNetCI); z = ComputeCI(f) - NetCI; end % % z = ComputeIEjectionRate(reactTerm, df, maxIdx) % % Computes the rate of change of silicon interstitials due to precipitate growth % for this time step. % % maxIdx is a parameter to allow us to avoid counting anything that is already % accounted for by the moments. This function sums from size 2 up until maxIdx. % function z = ComputeIEjectionRate(df, maxIdx) global n idxSizeOne; % Small clusters % ... we assume no interstitial ejection for these ... z = 0; % % Macroscopic precipitates: assume fixed ratio of interstitials. This is % justified because mOpt/n typically stays within a small range about 0.5. % Given the inaccuracies in the dislocation loop model, an error in the % exact number of ejected interstitials ejected is presumed to be less % important than an error in the strain energy, which effects precipitate % growth rates. Therefore, mOpt is used to compute energies but 0.5*n is % used to emit interstitials. % for i=(idxSizeOne+3):maxIdx if (i == length(n)) sampFactor = (2*(n(i)-n(i-1))/2); else sampFactor = (n(i+1)-n(i-1))/2; end ejected = 0.5*n(i); z = z + sampFactor*ejected*df(i); end end % % [m0, m1, m2] = ComputeMoments(f, startIdx, endIdx) % % Computes the moments of the distribution between the start and ending indices. % The input solution vector, f, must be one-dimensional (i.e., no time % components). % function [m0, m1, m2] = ComputeMoments(f, startIdx, endIdx) global n; m0 = 0; m1 = 0; m2 = 0; for i = startIdx:endIdx if (i == length(n)) sampFactor = (2*(n(i)-n(i-1))/2); else sampFactor = (n(i+1)-n(i-1))/2; end m0 = m0 + sampFactor*f(i); m1 = m1 + n(i)*sampFactor*f(i); m2 = m2 + n(i)*n(i)*sampFactor*f(i); end end % % [a, b, c, d] = Discretize(np, nc, nn, CO, SV) % % Computes the Kobayashi discretization coefficients for the solution variable % at index i (size nn, with size np being the previous sample point and nn being % the next). The solution will have the form: % % df(i)/dt = (a*f(i-1)-b*f(i)) - (c*f(i)-d*f(i+1)) % % The first term, a*f(i-1)-b*f(i) corresponds to the flux from n(i-1) -> n(i). % The second corresponds to the flux from n(i) -> n(i+1). When i is the end % point, c and d are invalid and should be ignored. The form of the equation is % then: % % df(i)/dt = a*f(i-1) - b*f(i) % % Alternatively, it may be possible to just assume df(i)/dt = 0 because care % should have been taken to make the maximum end point large enough that the % solution drops to 0 there. % % This function should only be used from i = idxChangeOver onwards. % % The derivation of the coefficients appears in Appendix B of: S. Kobayashi, % Journal of Crystal Growth, vol. 174, pp. 163-169 (1997). The same notation is % used here. % function [a, b, c, d] = Discretize(np, nc, nn, CO, SV) global n; mW = np; % previous sample point mP = nc; % current (center) mE = nn; % next delta_mP = (mE-mW)/2; delta_mw = mP-mW; delta_me = mE-mP; % PP and QP (roughly the inflow and outflow, relative to current point, nc) pP = GRate(mP,CO); qE = DRate(mE,SV); PP = (pP-qE)/(1-(qE/pP)^delta_me); %PP = (pP^delta_me)*(pP-qE)/(pP^delta_me-qE^delta_me); pW = GRate(mW,CO); qP = DRate(mP,SV); QP = (pW-qP)/((pW/qP)^delta_mw-1); %QP = (qP^delta_mw)*(pW-qP)/(pW^delta_mw-qP^delta_mw); % PW (analogous to PP above, w/ changes: P -> W, E -> P, delta_me -> delta_mw PW = (pW-qP)/(1-(qP/pW)^delta_mw); %PW = (pW^delta_mw)*(pW-qP)/(pW^delta_mw-qP^delta_mw); % QE (analogous to QP above, w/ changes: W -> P, P-> E, delta_mw -> delta_me QE = (pP-qE)/((pP/qE)^delta_me-1); %QE = (qE^delta_me)*(pP-qE)/(pP^delta_me-qE^delta_me); % Return coefficients a = PW / delta_mP; b = QP / delta_mP; c = PP / delta_mP; d = QE / delta_mP; end % % z = FkEstimator(m0, nAvg) % % Empirical f(k) estimator. % function z = FkEstimator(m0, nAvg) p = 2; p0 = 5e-6; p1 = 0.1; z = m0/(((p0*(nAvg-72))^p+p1)*(nAvg-72)+1); end % % df = ODECallback(t, f) % % ODE solver callback. Computes df/dt at time t for all solution variables. % function df = ODECallback(t, f) global modelType kBT n idxChangeOver idxNetCI idxSizeOne idxM0DL idxM1DL idxM0Hi idxM1Hi idxK surfReactDist reactDist CSi DO Css DI CIstar DV CVstar; % Update temperature and parameters UpdateTemperature(t); % Clear solution df = zeros(length(n),1); % % Silicon interstitials and vacancies: % % There is no explicit I/V recombination term and we assume that % CICV = CI*CV*. % % What is actually solved for is NetCI = CI - CV. % CI = ComputeCI(f); CV = ComputeCV(f); % Point defect supersaturations SI = CI/CIstar; SV = CV/CVstar; % Precompute small cluster formation energies for this time step UpdateSmallClusterEnergies(SV); % % Size 2 is special. O2 is mobile and an O2V state (one point defect % incorporated) exists. It is treated the same as other, larger small % clusters: the O2 and O2V states are assumed to be present in % proportion to their relative equilibrium concentrations. % % Relevant reactions are: % % O + O <-> O2 (f2 contains O2 and O2V) % f2 + O <-> f3 (all O3Vm precipitates) % f3 + O <-> f4 % % Note that O+O only forms the O2 (not O2V) state. The dissolution rate % for size 2 is adjusted to only take into account the clusters in the % O2 state: % % f2* = f2,0* + f2,1* % % Where f2,m* = CSi*exp{-Gtot_2,m/kBT} = CSi*(CO/CSi)^2 * % (CV/CV*)^m*exp{-Gf_2,m/kBT} % % f2,0*/f2* = exp{-Gf_2,0/kBT}/(exp{-Gf_2,0)/kBT}+ % (CV/CV*)*exp{-Gf_2,1/kBT}) % % Remaining solute CO = f(idxSizeOne); % Size 2 energies HfO2 = HfDiscrete(2,0); HfO2V = HfDiscrete(2,1); % O + O <-> O2 FractionO2 = 1/(1+SV*exp((HfO2-HfO2V)/kBT)); Rate_O_O_To_f2 = 4*pi*reactDist*DO*(CO*CO - CSi*exp(HfO2/kBT)*FractionO2*f(idxSizeOne+1)); % f2 + O <-> f3 Rate_f2_O_To_f3 = GRate(2,CO)*f(idxSizeOne+1) - DRate(3,SV)*f(idxSizeOne+2); % Sizes 2 and 3 df(idxSizeOne+1) = Rate_O_O_To_f2 - Rate_f2_O_To_f3; df(idxSizeOne+2) = Rate_f2_O_To_f3 - (GRate(n(idxSizeOne+2),CO)*f(idxSizeOne+2) - DRate(n(idxSizeOne+3),SV)*f(idxSizeOne+3)); % Conservation of solute df(idxSizeOne) = - 2*df(idxSizeOne+1) - 3*df(idxSizeOne+2); % % Sizes 4 to nChangeOver-1: Discrete chemical rate equations (sample % spacing is 1) % for i = (idxSizeOne+3):(idxChangeOver-1) df(i) = GRate(n(i)-1,CO)*f(i-1)-DRate(n(i),SV)*f(i) - (GRate(n(i),CO)*f(i)-DRate(n(i)+1,SV)*f(i+1)); df(idxSizeOne) = df(idxSizeOne) - n(i)*df(i); end % % Sizes nChangeOver and above: Interpolated rate equations. % % For full model, generate equations up to maximum size. For RKPM model, % only up to size k-1 (k is special and has to be estimated). % if modelType == 'f' maxIdx = length(n); else maxIdx = idxK - 1; end for i = idxChangeOver:maxIdx % Compute reaction rate and change in solute if (i == length(n)) % Compute interpolation coefficients assuming a fictitious sample % point same distance out from previous one [a, b, c, d] = Discretize(n(i-1), n(i), n(i)+(n(i)-n(i-1)), CO, SV); %df(i) = 0; df(i) = a*f(i-1) - b*f(i); df(idxSizeOne) = df(idxSizeOne) - (n(i)*2*(n(i)-n(i-1))/2) * df(i); else % Compute interpolation coefficients [a, b, c, d] = Discretize(n(i-1), n(i), n(i+1), CO, SV); % Equations if (i == (idxK-1)) && (i == maxIdx) % This occurs only in RKPM. The expression computed here is % incomplete and lacks the reverse term (outbound flux to f(k)). % It will be filled in later when dm0/dt is computed. df(i) = a*f(i-1) - b*f(i); % When df(k) is filled in later, df(idxSizeOne) will be updated else df(i) = a*f(i-1) - b*f(i) - (c*f(i) - d*f(i+1)); df(idxSizeOne) = df(idxSizeOne) - (n(i)*(n(i+1)-n(i-1))/2) * df(i); end end end % % Moments of high distribution % m0 = max(0, f(idxM0Hi)); m1 = f(idxM1Hi); k = n(idxK); nAvg = max(m1/m0, k); % % Full model case % if modelType == 'f' % % m1 (high part of distribution) % for i = idxK:length(n) if (i == length(n)) sampFactor = (2*(n(i)-n(i-1))/2); else sampFactor = (n(i+1)-n(i-1))/2; end df(idxM1Hi) = df(idxM1Hi) + n(i)*sampFactor*df(i); end % % m0 % [a, b, c, d] = Discretize(n(idxK-1), n(idxK), n(idxK)+1, CO, SV); Ik = a*f(idxK-1) - b*f(idxK); sampFactorK = (n(idxK)+1-n(idxK-1))/2; df(idxM0Hi) = sampFactorK*Ik; if (f(idxM0Hi) < 1e-10) && (df(idxM0Hi) < 0) df(idxM0Hi) = 0; end % % RKPM model case % else % % Size k is special. There is no solution at idxK+1 so we must estimate % f(k). % [a, b, c, d] = Discretize(n(idxK-1), n(idxK), n(idxK)+1, CO, SV); sampFactorK = (n(idxK)+1-n(idxK-1))/2; f(idxK) = FkEstimator(f(idxM0Hi), nAvg); % % Delta function approximation % % Compute gammas. gamma3 includes adjustment factor for when nAvg -> k. gamma2 = RatePrefactor(nAvg)/DO; gamma3 = DRate(nAvg,SV)/DO - DRate(k,SV)*f(idxK)/(DO*m0); % Compute dm1/dt % To-do: why is n(idxK)+n(idxK)-n(idxK-1) unstable for discretization? [a, b, c, d] = Discretize(n(idxK-1), n(idxK), n(idxK)+1, CO, SV); sampFactorK = (n(idxK)+1-n(idxK-1))/2; Ik = a*f(idxK-1) - b*f(idxK); % flux from k-1 -> k df(idxM1Hi) = sampFactorK*k*Ik + DO*m0*(CO*gamma2-gamma3); % % m0 and connection back to f(k-1) % % The use of sampFactorK here is in order to be consistent with % ComputeMoments(), otherwise, plugging in the full f(k) here will % not result in the same moments as those obtained from a full % simulation where ComputeMoments() is invoked at each time step. % df(idxM0Hi) = sampFactorK*Ik; % Recompute Ik interpolation as for outgoing flux at k-1 -> k [a, b, c, d] = Discretize(n(idxK-2),n(idxK-1),n(idxK),CO,SV); Ik = c*f(idxK-1)-d*f(idxK); % Connect it back to df(k-1) equation partially computed earlier df(idxK-1) = df(idxK-1) - 1*Ik; % % Conservation of solute (m1 and f(k-1)) % df(idxSizeOne) = df(idxSizeOne) - (n(idxK-1)*(n(idxK)-n(idxK-2))/2) * df(idxK-1) - df(idxM1Hi); end % % Dislocation loop nucleation and growth % nc = LoopCriticalSize(nAvg, SI); if (imag(nc) ~= 0) || (nc == NaN) nc = 0; end m0DL = max(0, f(idxM0DL)); m1DL = max(0, f(idxM1DL)); if m0DL ~= 0 m1DLhat = m1DL/m0DL; else m1DLhat = 0; % loop growth rate for n = 0 will return 0 end if (nc > 1e3) df(idxM0DL) = 0; df(idxM1DL) = (LoopGRate(m1DLhat, CI)*m0DL - LoopDRate(m1DLhat, nAvg, SV)*m0DL); else df(idxM0DL) = LoopGRate(nc, CI)*LoopFnStar(nc, m0DL, m0, nAvg, SI); df(idxM1DL) = nc*df(idxM0DL) + LoopGRate(m1DLhat, CI)*m0DL - LoopDRate(m1DLhat, nAvg, SV)*m0DL; end % % Point defects % % 0D boundary condition for point defects (approximation of diffusion to and % absorption at wafer surfaces. ISurfaceTerm = -((DI/surfReactDist)*(CI-CIstar)-(DV/surfReactDist)*(CV-CVstar)); % Oxygen precipitates eject interstitials, loops consume them df(idxNetCI) = ISurfaceTerm + ComputeIEjectionRate(df, maxIdx) - df(idxM1DL); if modelType == 'r' df(idxNetCI) = df(idxNetCI) + 0.5*df(idxM1Hi); end % For some reason, calling the progress indicator here manually is faster % than using a solver callback passed to ode15s! OutputCallback(t, f, []); end