Correct Compton shell selection and Doppler broadening

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.

[paulromano]

Description

The Compton Doppler-broadening algorithm on the current develop branch does not sample the bound-electron impulse approximation correctly. In particular, it selects a subshell using only the number of electrons in that shell. If the incident photon cannot ionize the selected shell, or if the kinematic upper bound on longitudinal electron momentum is negative, the code falls back to the free-electron Compton energy. However, a negative upper bound is not an invalid condition: it means that only part of the negative-momentum tail of the shell's Compton profile is accessible.

For shell $i$, the physically allowed longitudinal-momentum interval is

$$ -\frac{1}{\alpha_\mathrm{fs}} \le p_z \le p_{z,\max,i}, $$

where $\alpha_\mathrm{fs}$ is the fine-structure constant. The probability of selecting the shell is therefore proportional to

$$ P_i \propto f_i H(E-E_{b,i}) \int_{-1/\alpha_\mathrm{fs}}^{p_{z,\max,i}} J_i(p_z),dp_z, $$

where $f_i$ is the number of electrons in the shell, $E_{b,i}$ is its binding energy, and $J_i$ is its symmetric Compton profile. Thus, shell selection depends on both electron occupancy and the fraction of the profile that is kinematically accessible.

The develop implementation also samples only nonnegative $p_z$, chooses randomly between the two positive solutions of the scattered-energy quadratic, and does not apply the outgoing-energy rejection factor from the approximate relativistic impulse approximation. These choices distort the broadened energy distribution, particularly near shell thresholds and for forward scattering.

Relationship to #3874

#3874 correctly identifies that inaccessible shells should not be accepted and that negative values of $p_{z,\max}$ need to be handled. However, its proposed implementation has three blocking problems:

1. The negative $p_{z,\max}$ branch samples outside the tabulated CDF. OpenMC stores only the nonnegative half of each symmetric Compton profile, so its integrated value approaches approximately $1/2$, not 1. For $p_{z,\max}<0$, Fix a bug in compton scattering electron shell selection #3874samples a CDF value between the integral at $|p_{z,\max}|$ and 1. Values above the end of the half-profile CDF can make the subsequent lookup run past the profile table.

2. The sign of $p_z$ is lost and the wrong energy root can be selected.Fix a bug in compton scattering electron shell selection #3874 inverts the profile using $|p_{z,\max}|$ but leaves the sampled momentum nonnegative. It then always selects the lower positive root of the scattered-energy quadratic. The physical branches are instead determined by the sign of the sampled momentum: $p_z<0$ corresponds to an energy below the free-electron Compton energy and requires the lower root, whereas $p_z>0$ corresponds to an energy above it and requires the upper root.

3. Shell selection is not weighted by the accessible profile mass.Fix a bug in compton scattering electron shell selection #3874samples shells by electron occupancy and rejects only on the binding threshold. That threshold check is necessary but insufficient. At fixed energy and angle, different shells have different values of $p_{z,\max,i}$ and therefore different accessible fractions of their momentum distributions. Selecting among all energetically open shells using electron occupancy alone produces the wrong shell probabilities.

Proposed approach

This PR replaces the existing sampler with the approximate relativistic impulse approximation procedure described in Sec. 3.4.8 of Kaltiaisenaho's Master's thesis:

• Build an integrated half-profile for each shell, validate its tail, normalize the full symmetric profile, and precompute the negative-branch mass.
• Extrapolate the finite tabulated profile using the log-linear tail prescribed in the thesis, $J_i(p_z) = J_i(p_{z,N}) \exp\left[a_i(p_z-p_{z,N})\right], p_z > p_{z,N},$ and integrate and invert that tail analytically.
• Reject shells below their binding threshold or with no kinematically accessible momentum, and calculate their accessible profile
  mass.
• Sample a signed $p_z$ conditional on $-1/\alpha_\mathrm{fs}\le p_z\le p_{z,\max,i}$. When $p_{z,\max,i}>0$, select the negative and positive branches in proportion to their integrated profile masses.
• Solve the scattered-energy quadratic using a numerically stable formulation and choose the root associated with the sign of $p_z$.
• Enforce the shell binding-energy limit and apply the $E'/E$ rejection factor from the approximate relativistic impulse approximation.

The implementation uses a two-stage strategy to retain the accuracy of this procedure without paying the full cost of constructing the conditional shell distribution for ordinary collisions. It first makes two inexpensive occupancy-based shell proposals and accepts each according to its accessible profile mass. If neither proposal succeeds, it evaluates the kinematics for every shell and samples directly from the conditional mass function $f_i M_i$, where $M_i$ is the accessible profile mass. This is statistically equivalent to repeated shell rejection but avoids pathological runtimes when every shell has very little accessible mass, such as near forward scattering.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

A great PR!
I have a few concerns and suggestions.

[paulromano]

Thanks for the review @GuySten. Your comments/suggestions have been addressed.

[GuySten]

Looks good to me.
Thanks @paulromano for addressing the problem.

[amandalund]

Very nice correction @paulromano! The physics looks sound to me.