Reproducing Song et al. 2004: dry oxidation of three faces of 4H-SiC

One of the [studies we reproduced](/writing/interface-studies-reproduced) for [Questions we need Emily to answer](/writing/questions-for-emily-testing-a-theory-of-the-interface), the companion to [the essay](/writing/testing-a-theory-of-the-interface-with-an-ontology). Generated by `scripts/reproduce_studies.py` from `reproductions/song_2004.py`. Every value from the paper was read by an AI agent and has not been checked by a person; the sample and history were re-read from the PDF alone by a second agent that had not seen our earlier transcription.

The study and why we depend on it

Y. Song, S. Dhar, L. C. Feldman, G. Chung and J. R. Williams, “Modified Deal Grove model for the thermal oxidation of silicon carbide”, J. Appl. Phys. 95, 4953 (2004).

The challenge's reference for SiC. The graph holds 43 thicknesses read from its Figures 3 to 5, and the SiC growth program is the silicon growth operator unchanged with this paper's tabulated rate constants. The 1050 °C curves are held out: their constants are interpolated from the neighbouring temperatures, not read from the paper.

How well can this study be reproduced?

High for growth. One step, fully stated: dry oxygen at one atmosphere at a stated temperature. The weakness is in the data, not the history: thickness appears only on figures, and the paper does not say whether the points of a curve are separate samples.

The grade is ours and is about the paper as a source for simulation, not about the quality of the experiment: high means the sample, every step and its timing are stated well enough that a simulation needs no choices of ours; moderate that some choices are ours and are listed; low that the result depends on choices the paper does not constrain.

The sample and its history, played back

The sample is drawn from the simulation: layer heights and the read-outs are what the ISL programs compute at each moment of the history. Steps the language cannot state (a clean, a metallization, a measurement) are on the clock and marked as recorded, not simulated.

Does our sample and history match the paper's?

Left: the paper's own words, with the page of the PDF. Right: what we simulate. 5 same, 1 approximated, 2 assumed, 3 not in ours.

The timing of the history

Time is counted from the start of the first simulated step. A step the language cannot state takes no simulated time; where the paper gives its duration, the duration is in the last column. Stated means the paper gives the duration; partly stated that it gives a bound or a range; not stated; assumed that the duration in our history is ours.

C face (000-1), 1150 °C

Simulated time in this history: 6 h 00 min 01 s.

Across the samples

Timing the paper does not give

• Whether each point is a separate sample or one sample oxidized again and again. If cumulative, every point includes the earlier unloads and reloads.
• The moment time zero is counted from: loading, reaching temperature, or admitting oxygen. For the shortest times (15 and 30 min) this matters.
• Loading temperature, push rate and ramp.
• Cooling and unloading: nothing is said.
• The oxidation times themselves are only on the figure axes.

What the second reading changed

• Our stored recipe (one hold in dry oxygen at the stated temperature) is what the paper states. The second reading found nothing our encoding contradicts.
• The paper never says whether the points on a curve are separate samples or one sample re-oxidized. Our simulation treats each time as an uninterrupted oxidation from a clean surface.
• The Si-face activation energies rest on three temperatures: Figure 5 plots 950 and 1000 °C points but the tables have no constants there. We kept only the three tabulated temperatures.
• Figures 3 to 5 do not say whether the plotted thickness is from RBS or ellipsometry; in Figure 2 the two differ visibly at 6 h on the C face.

Their terms and ours

The simulation in our language

Each sample is one history file: a SAMPLE and a CONTROL_TRAJECTORY. This is the whole input; run it with PYTHONPATH=. python -m isl_harness.predict <file>.

C face (000-1), 1150 °C

# C face (000-1), 1150 °C
SAMPLE SiC_C_1150 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_C_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1150->1150 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1150->1150 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

C face (000-1), 1050 °C

# C face (000-1), 1050 °C
SAMPLE SiC_C_1050 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_C_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1050->1050 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1050->1050 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

C face (000-1), 950 °C

# C face (000-1), 950 °C
SAMPLE SiC_C_950 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_C_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 950->950 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 950->950 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

a face (11-20), 1150 °C

# a face (11-20), 1150 °C
SAMPLE SiC_a_1150 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_a_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1150->1150 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1150->1150 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

a face (11-20), 1050 °C

# a face (11-20), 1050 °C
SAMPLE SiC_a_1050 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_a_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1050->1050 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1050->1050 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

Si face (0001), 1150 °C

# Si face (0001), 1150 °C
SAMPLE SiC_Si_1150 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_Si_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1150->1150 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1150->1150 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

Si face (0001), 1050 °C

# Si face (0001), 1050 °C
SAMPLE SiC_Si_1050 {
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_Si_face
    CLEAN: not_stated
    DOPING: n_type
}
CONTROL_TRAJECTORY History {
    EVENT Clean KIND clean TIMING assumed                        # Acetone, methanol, 5% HF
    EVENT Load KIND load TIMING assumed                          # Loaded at once into the furnace
    RUN BeforeClock DURATION 1 TEMP 1050->1050 AMBIENT dry_O2 TIMING assumed VARY duration 1..1800
    RUN Oxidation DURATION 21600 TEMP 1050->1050 AMBIENT dry_O2 TIMING bound
    EVENT Unload KIND other TIMING assumed                       # Unloaded and cooled
    EVENT Measure KIND measurement TIMING assumed                # Thickness by RBS and ellipsometry
}

The programs the predictor ran on these histories: experiments/sic_interface_states/one_way_anneal.isl.template, experiments/song_sic/oxide_growth_sic.isl.template.

The model itself, experiments/song_sic/oxide_growth_sic.isl.template, as it is written. Every equation the simulation uses is in this text; the runtime holds no physics.

# ==============================================================================
# ISL MODEL: DRY THERMAL OXIDATION OF 4H-SiC (growth sector on a second material)
#
# The OPERATOR and OBSERVATION_MODEL blocks below are copied character for character
# from experiments/vitkavage_dit/oxide_growth.isl.template, the silicon program; a test
# checks that they are identical. Only the declared parameters differ.
#
# Song et al. (2004) fit their SiC data with the same linear-parabolic law, reading its
# two constants as including out-diffusion of the carbon product. The constants here are
# Arrhenius fits to that paper's Tables I and II, per crystal face
# (scripts/build_song_dataset.py). Nothing is fitted by us.
# Fields in {{braces}} are filled by scripts/run_song_sic.py.
# ==============================================================================
SAMPLE {sample_id} {{
    FAMILY: SiC_SiO2
    ORIENTATION: 4H_SiC_{face}_face
}}
STATE x_ox {{
    UNIT: cm
    INITIAL: x_i
    MIN: 0
    MAX: 1.0e-3
}}
CONTROL T {{ UNIT: K }}
CONTROL p_O2 {{ UNIT: atm }}
PARAMETER B0 {{ VALUE: {B0}; UNIT: cm^2 s^-1; STATUS: sourced; SOURCE: {source} }}
PARAMETER E_B {{ VALUE: {E_B}; UNIT: eV; STATUS: sourced; SOURCE: {source} }}
PARAMETER BA0 {{ VALUE: {BA0}; UNIT: cm s^-1; STATUS: sourced; SOURCE: {source} }}
PARAMETER E_BA {{ VALUE: {E_BA}; UNIT: eV; STATUS: sourced; SOURCE: {source} }}
# The constants above are already specific to the crystal face, so no further factor is applied.
PARAMETER r_or {{ VALUE: 1.0; UNIT: 1; STATUS: constant }}
# The paper sets the time offset to zero and reports a native oxide of about 2 angstrom.
PARAMETER x_i {{ VALUE: 2.0e-8; UNIT: cm; STATUS: sourced; SOURCE: Song et al., J. Appl. Phys. 95, 4953 (2004), native oxide about 2 angstrom, agent-read }}
PARAMETER p_ref {{ VALUE: 1.0; UNIT: atm; STATUS: constant }}
PARAMETER k_B {{ VALUE: 8.617333262e-5; UNIT: eV K^-1; STATUS: constant }}
CONTROL_TRAJECTORY {trajectory_id} {{
    RUN Oxidation DURATION {duration_s} TEMP {T_C}->{T_C} AMBIENT dry_O2
}}
# dx/dt = B / (2x + A), with A = B / (B/A).
OPERATOR F_LinearParabolicGrowth : OxideGrowth {{
    EQUATION: d(x_ox)/dt = (p_O2 / p_ref) * B0 * exp(-E_B / (k_B * T)) / (2 * x_ox + B0 * exp(-E_B / (k_B * T)) / (r_or * BA0 * exp(-E_BA / (k_B * T))))
}}
OBSERVATION_MODEL Obs_Oxide_Thickness {{
    TARGET: Thickness
    UNIT: cm
    EQUATION: Thickness = x_ox
}}
EVALUATION_CONTRACT {contract_id} {{
    FROZEN: true
    DATASET: data/derived/song_2004_sic_thickness.csv
    DATASET_DIGEST: sha256:{dataset_sha256}
    SPLIT: {split}
{targets}
}}

Our simulation against the paper

21 within, 3 near of 24 numbers. Tolerance is stated per row; nothing was fitted for this page.

Aligning the simulation: what the unstated things would have to be

Every step whose timing or temperature the paper leaves open carries a range in the history above (VARY). scripts/align_study.py searched those ranges against the paper's numbers through the same predictor, first one unknown at a time, then all together. No parameter a program declares fitted was touched. The unknowns were chosen on C_1150, C_950, Si_1150, a_1150 only; every other number is a test they played no part in.

Together: BeforeClock duration 1 s (was 1).

As reproduced: 21 of 24 within tolerance, median error 3.0%. Aligned: 21 of 24 within tolerance, median error 3.0%. On the held-out samples only: 8 of 10 within tolerance, median error 6.5% as reproduced, 8 of 10 within tolerance, median error 6.5% aligned.

Choices of program or material property, each scored on its own and on top of the aligned history:

What was searched. The paper states the gas, the temperature and the times on the figure axes. It does not say when time zero is, or whether the points of a curve are separate pieces. Time zero was written as one unknown step before the oxidation, time in oxygen at temperature before the clock starts (1 s to 30 min; the language does not allow a zero-length step, so 1 s stands for none). The opposite sign, a clock that starts before the oxygen arrives, and the cumulative reading were declared as variants. The unknown was chosen on the 1150 and 950 °C curves (14 points); the 1050 °C curves (10 points), whose rate constants we interpolate, are the test.

A long offset is ruled out; a short one cannot be told from none. On the 14 chosen points, 13 are within tolerance for any offset from 1 s to 148 s (median error 1.1% at 1 s, 2.4% at 148 s), 12 at 8.6 min (3.0%) and 8 at 30 min (7.3%, every point high). The search keeps 1 s. This is not independent of the paper: Song et al. fitted their rate constants with the time offset set to zero, and our program uses those constants, so a zero offset is built into the numbers being compared.

The held-out 1050 °C curves do not change: 8 of 10 within tolerance, median error 6.5%, before and after. The two that fail are the Si face at 2 h and 6 h, 13% and 14% high; the C face at 1050 °C is 5 to 8% low at every time. Those are errors of the interpolated constants, not of the timing.

The variants. A clock that starts 5 min before the oxygen reaches the sample: 21 of 24 (median 3.4%), and on the 1050 °C points alone 8 of 10 (median 8.2%, computed from the same rows); 10 min: 20 of 24 (median 6.2%), 7 of 10 at 1050 °C. One sample per curve, oxidized cumulatively, with 5 min in oxygen before the clock at each reload: 21 of 24 (median 5.3%); on the chosen points the median rises from 1.1% to 4.1%, on the 1050 °C points it falls from 6.5% to 5.4% (signed +0.2%). With 15 min per reload: 12 of 24 (median 10.1%, every point high). So the data allow at most a few minutes of extra oxygen per reload if the curves are cumulative, and they do not choose between separate pieces and one piece re-oxidized.

8 within, 2 near of 10 numbers with the aligned settings (the samples used to choose the unknowns are left out).

The thickness curves cannot test their own time zero, because the program's rate constants were fitted to these curves under the paper's assumption of no offset. A test needs either the paper's 0.25 h and 0.5 h points (not in our dataset) or another laboratory's curve at one of these temperatures.

The remaining error at 1050 °C is in the constants: at 1050 °C the program's B and B/A come from an Arrhenius fit to the paper's tables, interpolated from the neighbouring temperatures rather than read from the paper. No history change reaches it.

What our language and ontology cannot yet say

• A clock reference for a RUN: TIME_ZERO at_load | at_temperature | at_gas_on (or a signed OFFSET <s>), so that “time” on a figure axis can be stated as what it is. Here the offset had to be an extra 1 s step, because a RUN cannot have zero duration and VARY cannot go negative.
• Measurements part-way through a history: MEASURE <observable> AFTER <step> followed by more steps, so a cumulatively oxidized sample (oxidize, unload, measure, reload, oxidize) is one history with several measured points. Today each point has to be a separate sample built by a variant function.
• An optional or repeated step (REPEAT <n> { … }) for the reload of a cumulative series, with its own VARY on the per-reload time.
• OFFCUT as a field a program reads, so that this on-axis data and Kobayashi's 4° off-axis epilayers are not the same sample to the growth program; and a word for the measurement technique (RBS or ellipsometry) on an observable.

With heating and cooling transitions

No paper we hold states how its wafers heated up or cooled down. The project owner specified scenarios for one Vitkavage wafer (isl_harness/vitkavage_transitions.py); isl_harness/thermal_transitions.py applies the same numbers to every sample here: a heat-up in nitrogen before the first furnace step (a hot push of 1 to 5 min, or a cold-furnace ramp of 60 to 120 min), the owner's radiative cooling knots after the last furnace step or after the paper's own pull (fast, or every time multiplied by 5), and ramps around an anneal at 300 to 500 °C. The paper's own steps are kept as written. All transitions are assumptions, written in the history with TIMING assumed.

The paper's own conclusions

What we would like Emily to check

1. In work of this kind, would each oxidation time be a fresh 5 mm piece, or the same piece put back in the furnace?
2. Is the thickness in Figures 3 to 5 the ellipsometric one?
3. For the Si face, are three temperatures enough for you to accept an activation energy of 3.12 eV for B?