[question] Sputtering yield comparison between RustBCA, SDTrimSP, and Eckstein

[HexSygon]

Hello Jon Drobny,

I have noticed a difference between computed Yc-c and Yd-c sputtering yields RustBCA, SDTrimSP, and Eckstein. Overall, the sputtering yields are lower for RustBCA. Please see the attached images for reference.

My colleagues have raised concerns about the accuracy of RustBCA. The sputtering yield code is relatively simple, and I've already verified it with you, so I’m confident it’s correct. I’ve also rerun the calculations and obtained the same results.

I’ve reviewed your manuscript, RustBCA: A High-Performance Binary-Collision-Approximation Code for Ion-Material Interactions, and noticed the differences between the sputtering yield methods. It seems that while RustBCA incorporates some simplifying assumptions, it also includes improvements that enhance accuracy when compared to SDTrimSP. Unfortunately, I don’t have access to TRIDYN for further comparison.

Sputtering and reflection yields for Tungsten, Carbon, and Deuterium combinations are crucial for my research. Therefore, any insight you can provide regarding the accuracy of RustBCA’s computational method, especially in comparison to SDTrimSP and Eckstein, would be greatly appreciated. Ultimately, I need to justify my decision on whether to use sputtering yields from RustBCA or SDTrimSP, and your expertise would be invaluable in guiding this choice.

Thank you for your continued support,
Alec Cacheris

[drobnyjt]

Thanks @HexSygon for writing and including the figures.

A weakness of the binary collision approximation method is that sputtering yields can be sensitive to a handful of input parameters - particularly, the binding energies, choice of electronic stopping mode, etc.

As a first step, can you confirm that all of the relevant inputs to RustBCA match the inputs used in the cited papers? If you can also provide the specific inputs to RustBCA that would be helpful.

[drobnyjt]

@HexSygon To set expectations, since the codebases are quite different I do not expect the results to be exactly the same even if the inputs are exactly the same, but based on my comparisons to (F-)TRIDYN, which should have a very similar codebase to SDTrimSP, using the same input parameters and analogous options should result in minimal differences in reflection coefficients, sputtering yields, etc. Here's some of the benchmarks in my dissertation, which can be found here:

Note: there is a typo in the above caption; these are sputtered energy-angle distributions.

The Layered Geometry example also includes a comparison to F-TRIDYN.

Additionally, there are some comparisons to experimental data in the Benchmarks page on the RustBCA Wiki, but these are not as comprehensive as I would prefer.

[HexSygon]

The only inputs I specify are the ion, target, energy, angle, and number of samples. For example, I run: sputtering_yield(deuterium, carbon, 100, 45, 1000000).

It appears that the binding energy values in the materials script for deuterium, carbon, and tungsten do NOT match the "Benchmark" page you included. I shall change these values to the "Benchmark" values and double check that they match the values used in SRIM and Eckstein (and that those values are appropriate).

I will not run this immediately, so please let me know if there is anything else to look into. I'm not sure how to specify the type of electronic stopping power or how crucial this is for instance.

[drobnyjt]

@HexSygon Thanks for getting back so quickly; that's right, I forgot you were using the python interface.

For better or for worse, I made some fairly aggressive assumptions designing the python interface because, admittedly, there are a lot of opaque options that can affect absolute quantities like the sputtering yield but don't usually change trends w.r.t. changing energy, angle, etc.

There are, however, some options that can change results significantly, but only in certain edge cases. Particularly, weak collisions, which I believe were originally developed for TRIM.SP, can improve the results for light ions but break the results for very heavy ions by artificially increasing the nuclear stopping power, as can be seen in this closed issue, which I resolved by defaulting weak collisions to zero. I chose the default options in an attempt to cover the widest range of ions, targets, and energies, but they will not be universally applicable.

These options include:
weak collision order = 0
high energy free flight paths = False
mean free path model = amorphous
electronic stopping mode = Lindhard-Scharff
interaction potential = Kr-C
scattering integral = Mendenhall-Weller

Edited to add: I believe all of these options are analogous to options available in SDTrimSP, but are hard-coded in TRIDYN and SRIM.

As of now, there's no way to change these options in the Python interface, although as you have seen, the material and particle parameters are just saved as dictionaries in materials.py and can be manually adjusted. In cases where a direct comparison to another code is warranted, I would recommend running RustBCA as a standalone code with a fully completed input file.

If it does turn out that the results are still significantly different even with all the inputs and options being the same, that suggests to me there's something unexpected going on, and I would definitely work to resolve that.

[HexSygon]

So it sounds like I should rerun the codes with the changes in binging energy as a first step?

What options do you think I should manually include in materials.py and how do you write this into the code?

[drobnyjt]

@HexSygon Yeah, I think so - ideally you would match every parameter to the source you're comparing to. Most important to change in materials.py to affect the sputtering yields are the surface binding energy Es and bulk binding energy Eb of the target.

Unfortunately for repeatability, people sometimes use those energies as free parameters in BCA codes to fit experimental data, so if the sources you're comparing to don't report those values it will be difficult to match their results.

[HexSygon]

Hey @drobnyjt, it looks like this approach worked! It seems that changing the C 'Eb' from 3 to 0 (to represent the amorphous material) made a big difference.

I do have another question, though, and I hope it’s not too basic. It’s about the atomic density of each ion/material listed in the material script. I know that the atom density (n) is calculated using the formula:

n = (ρ * N_A) / M;

where "n" is the atom density (atoms per unit volume), "ρ" is the mass density of the substance (mass per unit volume), "N_A" is Avogadro's number, and "M" is the molar mass of the substance.

Does this 'n' value reflect the incoming ion density or should I set an ion 'n' value to what I see in my experiement? I use RustBCA for DIII-D ion-divertor material purposes and I know my plasma ion density for each species, is this the value I should be setting as 'n' for the ion while the target can remain to be the RustBCA preset value?

Sorry if this is a bit much, just trying to wrap my head around everything.

Please let me know what you think and thanks again.

[drobnyjt]

That's great to hear, and makes a lot of sense. In previous work I've done on divertor materials we assumed it was zero for a similar reason. Since the bulk binding energy is subtracted from every atom that starts in the target material, it primarily increases the sputtering threshold energy and decreases sputtering yields.

No worries; your question is a good one that touches on a conscious design decision I made when writing RustBCA. First things first, the incident ion species defined in the [particle_parameters] block in the full input file do not have a prescribed density or even fluence, just a number of samples, also known as "computational ions", denoted N, or num_samples in the Python interface. A quick note here, just in case - RustBCA will not emit a warning or error if you include a non-existent option for a block in an input file, it will just ignore it, so if you put a line like n = [ 1e19 ] in the [particle_parameters] block, it will be, for better or worse, silently ignored. Only the [material_parameters] block that describes the target has number densities, and those are as you have described.

Unlike TRIDYN, which asks for an ion fluence and includes an ad hoc diffusion and erosion model to try to capture the evolution of the target over time, RustBCA is intended to be used as a Monte Carlo solver for only one instant in time where the material remains completely static. This choice was made to make RustBCA more useful for integrated modeling, where other physics - diffusion, erosion, etc. - can be handled by distinct high fidelity, time-dependent models as opposed to being hard-coded into RustBCA.

If you are looking to calculate flux- or fluence-dependent quantities, or time-dependent quantities like erosion in nm/s, the best way to do that is to run RustBCA with a sufficiently high number of computational ions to get the precision or smoothness you need, and then multiply the appropriate yield or distribution function by the flux or fluence to get the instantaneous value of that quantity. To give an explicit example, the sputtering yield Y from RustBCA can be used to calculate the instantaneous erosion rate like so:

dx/dt = Y [atoms / ion] * incident ion flux [ions / m2 / s] / target atomic density [1/m3]

Hopefully that makes sense and paints a better picture of where RustBCA is intended to fit in when modeling plasma-material interactions; definitely feel free to ask any more questions you have if what I've written isn't clear.

[HexSygon]

Hey jon,

I just realized that you stated: " (For) divertor materials we assumed it (Eb) was zero for a similar reason. Since the bulk binding energy is subtracted from every atom that starts in the target material, it primarily increases the sputtering threshold energy and decreases sputtering yields."

Changing the Carbon Eb from 3 to 0 resulted in a decrease sputtering threshold and icnreased C-C sputtering.

Did you mean to say something else or is their another point I am missing? Perhaps you mean that changing C Eb to zero reduces C-W sputtering? I'm going to look into this if so.

[HexSygon]

This also makes me wonder if the W Eb preset of 0 is appropriate or not...

[drobnyjt]

Sorry for the confusion, I meant that a non-zero Eb increases the sputtering threshold. Your observations are correct.

Eb=0 for tungsten in materials.py because that was what was used in F-TRIDYN for the PSI-SciDAC project predating RustBCA (and for which RustBCA was originall developed). I believe the original default value, dating back to TRIM, was 3.0 eV, which if I recall correctly, was used for all metals for which a theoretical value was unknown.

[HexSygon]

Intersting, so there could be a large uncertainty in the tungsten Eb. Using W Eb =3 rather than zero would reduce my sputtering yields and align my results better.

Do you have any more information or references regarding W Eb? I found some talking about W Es but not Eb.

[drobnyjt]

I know of a couple sources that discuss it off the top of my head, although they are not specific to tungsten. Eckstein's 1991 text Computer Simulation of Ion-Solid Interactions has the following to say:

"
The bulk binding energy Eb is used because it is assumed that a target atom is
bound at its lattice site with a certain energy. The energy needed to remove a
target atom from its lattice site by an elastic energy transfer is Eb so that the
recoiling atom has an energy E = T - Eb after leaving its lattice site. There
is disagreement about the energy that should be used for Eb. Some argue that
Eb should be equal to the vacancy formation energy, which is of the order of
one to a few eV (see Table 6.1 for a few data taken from [6.10]), while in some
MARLOWE [6.22] and most TRIM.SP [6.23] calculations Eb = 0 is chosen.
"

Hofsäss and Stegmaier recently proposed a formulation of the binary collision approximation that replaces the concept of surface binding energies with only bulk binding energies (that are larger in value). It includes the following passage about previous BCA codes:
"
The nuclear collisions within the collision cascade are either treated as inelastic (without obeying momentum conservation) by using a bulk binding energy of typically Eb = 3 eV or (as default) elastic with zero bulk binding energy. TRIDYN as well as SDTrimSP both use Eb = 0 eV as default value.
"

So either a few eV or 0 eV are common assumptions, and Hofsäss and Stegmaier go further and assume Es=0 and increase the bulk binding energy to compensate, but as far as I know there is little persuasive evidence that one option or the other is necessarily more correct. In principle it should be possible to check using molecular dynamics, but I do not know anyone who has done so in detail.

[HexSygon]

Very interesting Jon. I’m forwarding this message along to my colleagues. Thank you very much for this. It’s probably reasonable to do a parameter scan and simply discuss the differences.

I will also compare sputtering and reflection yields between single and compound materials and keep you in the loop. I will certainty reference your RustBCA work.