[bug?] unusual sputtering yield curves when Z1 >> Z2 with weak collisions enabled

[drobnyjt]

Description A user has pointed out an unusual result - the sputtering yield curves for normal incidence when M1 >> M2 have an unphysical appearance:

It is not yet clear if this is only happening using the default settings, implying a configuration error, or if there is a deeper, more subtle bug. Further investigation is warranted.

To Reproduce

from libRustBCA import *
import numpy as np
import matplotlib.pyplot as plt
from scripts.materials import *

carbon = {
    'symbol': 'C',
    'name': 'graphite',
    'Z': 6.0,
    'm': 12.011,
    'n': 1.136e29,
    'Es': 7.41,
    'Eb': 0.0,
    'Ec': 3.0
}

num_samples = 10000
angle = 0.0
energies = np.linspace(1, 500, 25)
ycsb = [sputtering_yield(cesium, boron, energy, angle, num_samples) for energy in energies]
yxec = [sputtering_yield(xenon, carbon, energy, angle, num_samples) for energy in energies]
plt.plot(energies, ycsb, label='Cs on B')
plt.plot(energies, yxec, label='Xe on C')
plt.title('Sputtering with libRustBCA')
plt.xlabel('E [eV]')
plt.ylabel('Y [at/ion]')
plt.legend
plt.gca().set_yscale('log')
plt.show()

Expected behavior The sputtering yield curve should monotonically increase as a function of energy until the sputtering yield peak.

Error messages, output files, or figures N/A

[drobnyjt]

Possible causes:

• electronic stopping (when Z1 >> Z2) turning electronic stopping off did not change results.
• initial particle position and boundary conditions
• weak collisions with Z1 >> Z2 (and/or M1 >> M2)
• zero bulk binding energy issue Changing Eb to 3.0 does nothing
• surface binding energy issue

[drobnyjt]

The issue has been identified - when weak collisions are on for cases where Z1 >> Z2, those weak collisions are artificially sputtering near-surface atoms.

If you make the following addition to the example above:

xenon['Z'] = 1.0

You get the following result instead:

It shows that this is indeed caused by Z1 >> Z2 and not M1 >> M2. Since Kr-C scales linearly with Z1*Z2, the "weak" collisions are actually strong enough to cause unphysical sputtering. Weak collisions are simulated at distances corresponding to N average atomic spacings away (see figure below).

Figure 1: a particle in the BCA has one strong collision with p1, one weak collision with p2, and ignore p3 because it is outside of the surface.

sputtering_yield originally defaulted to 3 weak collisions, which was default in TRIDYN - this was because weak collisions better capture the reflection of and sputtering by light ions. It is apparent that, in the case of Z1 >> Z2, weak collisions should not be used. The solution to unphysical sputtering yield curves is to set weak collisions to zero when Z1 >> Z2.

The development questions are the following:

• What should the default number of weak collisions be for the ergonomic python functions?
• If the answer remains 3, should the code emit a warning about Z1 >> Z2?
• Should it be user-specified per-material like Eb?

[dcurreli]

This is effectively a free parameter of the problem like Eb, so it should be defined per each material.It might be possible estimating it using ab-initio methods for a given solid.

[drobnyjt]

I agree, but after thinking about it I am going to investigate the weak collision loop to make sure that it's conserving energy correctly - if each weak collision isn't being initiated with the correct energy (E - T[n-1], where E0 is the incident particle and T[i] is the energy transferred in the ith collision) then that could also cause the unphysical sputtering yield curves beyond the weak collision approximation breaking down, since it would effectively generating fictional energy at every step of the cascade.

[drobnyjt]

NOTE: the x-axis should be micron, not angstrom.

Changing k (weak collision order) dramatically increases the nuclear stopping power for Z1>>Z2; this is one obvious cause of the unphysical sputtering yields - the Xenon is depositing much more energy near the surface than normal. There is still an unusual very low energy (<0.1eV) spike in the sputtered energies produced by this process - I'm continuing to investigate the origins of those sputtered atoms.

[drobnyjt]

Reopening issue for now - I still want to address the collision loop in a more complete way.

[drobnyjt]

As far as I can tell, there is no underlying bug - it looks like the problem is simply that weak collisions cannot be used when Z1>>Z2. At issue is the fact that I have, as of yet, been unable to figure out why exactly this is so - my best guess is that in reality, when Z1>>Z2, the "weak" collisions using Kr-C do not adequately account for the screening of the interaction potential by the material.

For this reason, I've removed the bug tag, as the default configuration issue has been resolved.

[drobnyjt]

I'm going to close this; although the theoretical reasons behind weak collisions not working for Z1 >> Z2 are interesting, investigation into them should be a separate issue. I consider this issue resolved by changing the default options.