Unexpected anisotropy in elastic scattering of low-energy neutrons at finite temperature

Dear Geant4 experts,

I am currently simulating elastic scattering (hadElastic) of low-energy neutrons (in the meV range), and I have encountered unexpected behavior related to thermal vibrations.
I would appreciate your insights on this issue.

The simulation is based on a modified version of the B4 example.
I am using Geant4 version 11.3.2.
In the DetectorConstruction class, I define a 1 cm × 1 cm × 1 cm hydrogen sample using G4Material.
The world volume is set to Galactic, and the hydrogen sample is placed at the origin.
The physics list is set to FTFP_BERT_HP.
In the PrimaryGeneratorAction, I define an incident neutron with a kinetic energy of 25 meV, originating from the position (0, 0, -100 mm) and directed along the momentum vector (0, 0, 1).
In the SteppingAction, I record the final momentum vector of the primary neutron (trackID = 1) after interactions.
I store both the raw (non-normalized) and normalized vectors.
The total number of scattering events** that occurred within the sample is also recorded.
These data are passed through EventAction and RunAction and written into a ROOT file.

I ran the simulation at three different sample temperatures: 0 K, 300 K, and 10000 K, to study the temperature dependence of neutron energy transfer and scattering angles.
For simplicity, I analyzed only events in which the neutron scattered exactly once.
As shown in the attached plots, I observed unexpected behavior.
The figures show the x (left), y (center), and z (right) components of the final momentum vector (non-nomarized) of neutrons that scattered once.
Black lines indicate 0 K, red for 300 K, and green for 10000 K.
At 0 K, where the sample has no thermal motion, neutrons do not gain energy from the sample.
As expected, the distributions of momentum components are nearly flat around 0, with the z-component boosted due to the neutron’s initial direction along +z.
However, at finite temperatures (300 K and 10000 K), I found that the x- and y-components of the momentum vector exhibit different distributions.
This is surprising because the simulation should be symmetric around the z-axis, and such anisotropy is not expected.
Could anyone advise me on what might be causing this asymmetry and how to resolve it?

Thank you very much in advance.

Geant4 Version: 11.3.2
Operating System: macOS 12.7.4
Compiler/Version: Apple Clang 14.0.0
CMake Version: 3.26.4