Change alpha particle momentum initialization

[PhysicsDan]

Currently, the alpha particle energies are initialized based on the two-step process:
$\alpha_0$ channel
p + $^{11}$ B $\rightarrow$ $^8$ Be + $\alpha$ + 8.59 MeV
$^8$ Be $\rightarrow$ 2 $\alpha$ + 0.09 MeV

The dominant reaction channel is the $\alpha_1$ channel:
p + $^{11}$ B $\rightarrow$ $^8$ Be* + $\alpha$ + 5.557 MeV
$^8$ Be* $\rightarrow$ 2 $\alpha$ + 3.124 MeV

I have modified the two lines where the fusion and decay energies are defined to use the $\alpha_1$ channel.
In the future, this could be modified to calculate the probability of each channel but this is a simple update that should match the expected alpha distribution better.

Paper with plots of the two channels and nuclear energy levels: https://doi.org/10.1016/j.radphyschem.2022.110727

@roelof-groenewald - this may be of interest to you

ToDo

• Check tests pass correctly
• Update the checksum for the tests as alpha particle energies have changed

[PhysicsDan]

It looks like there is an issue with the pB fusion test.
It is failing on line 497 of the analysis
assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, rtol=1.e-2))
I will have a look and see if I can find the issue

Update:
Looking at 2D only for now.
PASS: Conservation of energy check. The current threshold is 1e-8, at 1e-6 the check passes for all cases.
PASS: The energy of the alpha particle emitted in step one p(11B, 8Be*)alpha is within tolerance of analytical calculation.

FAIL: The energies of the alpha particles produced in the second stage (Be8* -> 2 alpha) are not within the relative tolerance of the analytical calculation. In particular, there is a larger variation in the function check_initial_energy2.

[PhysicsDan]

Now the code reproduces the expected alpha particle energy distribution with this simple model.

This Figure is a plot retrieved from the test simulation and matches the distribution in Spraker et al (2012) quite well.

Why the test is failing?
The minimum alpha particle energy in the simulation does not match the value calculated by the analysis script. The discrepancy is in the range of 10 - 120 keV. Once this issue is fixed the test will pass.

I can't see where the issue comes from - if anyone has any suggestions let me know.

[roelof-groenewald]

Thanks for the update @PhysicsDan! I'll take a look at the changes and also see if I can understand where the test failure is coming from later this week.

[PhysicsDan]

Thanks, I had another look and can't find the issue.

The values that the simulation outputs are compared against are calculated using:

max_p_alpha23 = 0.5*(np.sqrt(p_sq_after) + \
                np.sqrt(4*m_a*energy_alpha2_plus_3_theory - p_sq_after))
min_p_alpha23 = 0.5*(np.sqrt(p_sq_after) - \
                np.sqrt(4*m_a*energy_alpha2_plus_3_theory - p_sq_after))
max_energy_alpha23 = max_p_alpha23**2/(2.*m_a)
min_energy_alpha23 = min_p_alpha23**2/(2.*m_a)

The max and min energy values for the momentum are calculated using the same variables but only the calculation of the maximum alpha energy passes.
e.g.
alpha3 min (FAIL)
Simulation: 0.04063286823092652 MeV
Expectation: 0.03910416430516294 MeV
alpha3 max (PASS)
Simulation: 5.26719008142914 MeV
Expectation: 5.275499859832506 MeV

This most recent commit modifies the energy released in the fusion step and decay step to make sure it matches the energy calculated based on the $\Delta mc^2$ value.

E_fusion: 5.55610759 MeV
E_decay: 3.12600414 MeV
E_total: 8.68211173 MeV # From Janis text book (8.68212 MeV)

Another note: I noticed all the masses include the electron mass.
I have continued with this convention although I don't understand the reason.

[roelof-groenewald]

I haven't had a chance yet to investigate the test failure but I wanted to comment so long on your note

Another note: I noticed all the masses include the electron mass.
I have continued with this convention although I don't understand the reason.

We are aware of this, but opted to consistently use masses with electrons included in the fusion module (see #4480). The reasoning was that in other parts of WarpX we allow ionization of a species but, for simplicity, that process does not change the particle mass. Therefore, WarpX already has other modules where we assume keeping the electron mass is an okay approximation and we decided to stick with that approach in the fusion module.

[roelof-groenewald]

@PhysicsDan I believe I understand the reason for the error you are seeing and it is not actually due to a bug or an error, just statistical variation. The failing test is meant to check whether the range of alpha energies fall within the expected range, but is doing this by checking that the maximum and minimum values are close to the theoretical min and max values. How close the observed energies will get to those bounds is simply a matter of statistical sampling. I ran the test using different random seeds and recorded the difference between the observed energy bounds and the theoretical ones. The results were as follows (I show the minimum bound differences as negative values to make the symmetry between max and min "errors" clear):

We see that the observed "error" varies strongly with the random seed but more importantly the magnitude of the difference between observed and theoretical energy bounds is the same for the min and max bounds. Therefore, I believe that the only reason the max-value checks are passing while the min-value is failing is because a relative error is used so the min values (being close to 0) has a lower tolerance than the max values. I suggest we change the test to check that indeed the minimum observed energy is >= to the the theoretical min (and vice versa for the max) and perhaps use a second check for values close to the bounds with absolute tolerance set to ~0.01 MeV.

[PhysicsDan]

Thanks, good spot!

Here is just changed the test tolerance to:
20 keV for test 1
and
200 keV for test 2 (The statistics are bad for this test hence the tolerance had to be higher).

I left the tolerance of the maximum values unchanged.

[roelof-groenewald]

Could we change this and the other similar calls to the following instead?

Suggested change

	assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, atol=3.218e-15)) # 0.02 MeV abs tol
	assert(0 <= np.amin(energy_alpha2_simulation) - min_energy_alpha23 < 3.218e-15) # 0.02 MeV abs tol

That way we simultaneously check that the minimum energy is not less than the theoretical minimum and that the minimum value is close to the expected minimum. I think it would then also be good to swap the maximum energy checks to the same approach otherwise someone looking at the code in the future might be confused by the asymmetry. What do you think?

[PhysicsDan]

I made these changes. However, the test now fails.

  File "./analysis_proton_boron_fusion.py", line 587, in check_initial_energy2
    assert(0 <= np.amin(energy_alpha2_simulation) - min_energy_alpha23 < 3.218e-14)

It appears in check_initial_energy2 that the difference between the limits. I.e, the np.amin(energy_alpha2_simulation) - min_energy_alpha23 < 0 and
np.amax(energy_alpha2_simulation) - max_energy_alpha23 > 0.

check_initial_energy1 passes - we are in the center of the momentum frame

check_initial_energy2 fails - we are in the boron ion rest frame
I assume that there is some small error in the transformation between reference frames causing this.

[roelof-groenewald]

That is unfortunate! Do you know if this same issue was present before changing the reaction pathway (this PR)? If so, we can change the check_initial_energy2 function back to the way you had it when the test passed (and merge this PR) and create a new issue to look into the origin of the energy exceeding the expected limits.

[PhysicsDan]

No, I'll have a look later this week and let you know. Yes, no problem if it is an issue in the previous version I'll revert the test and create an issue.

[PhysicsDan]

Original Version (check_initial_energy2). This is what we expect.

My branch:

For some reason, the minimum value of the 2nd and 3rd alphas generated is less than calculated in the analysis.

The larger absolute variation is ok because the energy of the 2nd and 3rd alpha particles is also slightly larger via the alpha1 channel.

[PhysicsDan]

Apologies for the spam. I had a look at the check_initial_energy1 as well.

Main Branch (FAIL)

My branch (PASS)

TLDR:
The main dev branch will fail if we check if the alpha particle energy is within the limits we calculate theoretically in the check_initial_energy1 function (both max and min are wrong) but pass the check_initial_energy2.

My branch will pass check_initial_energy1 but fails check_initial_energy2 as some alpha particles are generated with a lower energy than expected. The maximum values are correct.

Could this be some floating point issue that is amplified by the Lorentz transformation that will be performed during the simulation for check_initial_energy2? Perhaps this is only apparent as the minimum alpha energies are ~x7 lower in the alpha1 channel compared to alpha0 channel.

Summary: In my opinion, this PR is still a net benefit as the alpha spectrum generated matches better what has been observed experimentally (despite the issue above).

If you agree I can revert the tests so that they pass with the updated E_fusion and E_decay values

[roelof-groenewald]

Thanks for checking @PhysicsDan. I agree with you that since the current development version also seem to fail this more rigorous check, albeit in the first test instead of the second, we can merge this to gain the benefit of the more accurate decay channel. So please do revert the check in check_initial_energy2 to the case where the test passed and let's also summarize the observed discrepancy in a new issue.

[roelof-groenewald]

So just to make sure, the problem is actually only observed with the second case (where the COM frame and labframe are not the same)? If so, the check for energy in case 1 could be as we discussed previously (checking the energy range falls strictly between the theoretical min and max energy) while the second case's check could then be updated once the underlying issue is fixed.
That change for the case 1 tests is not necessary before merging. I just wanted to point out that case 1 does not have any remaining issues.

[PhysicsDan]

Yes to confirm check_initial_energy1 does pass with

        assert(0 <= max_energy_alpha23 - np.amax(energy_alpha2_simulation) < 3.218e-15)
        assert(0 <= np.amin(energy_alpha2_simulation) - min_energy_alpha23 < 3.218e-15)
        assert(0 <= max_energy_alpha23 - np.amax(energy_alpha3_simulation) < 3.218e-15)
        assert(0 <= np.amin(energy_alpha3_simulation) - min_energy_alpha23 < 3.218e-15)

There are no issues with case 1 in this PR. However, this test would fail with the previous alpha0 channel implementation.

[roelof-groenewald]

Thanks for updating this!

Please remember to create a new issue from the observed inconsistency with the alpha particle energies and link this PR to it. (You could even create an issue from the comment in which you explain the observed errors).

[PhysicsDan]

Yes no problem, I will do that now.

Just to note I had to change the check that compares the maximum values to a relative tolerance. The percentage variation is ~constant across a range of proton energies.