Here are my 2cts, after having spent quite some time trying to understand the role of D50 in ICC profiles and applying the math to see if I got the concepts right.
I think that one of the main points is that one should not look at the intermediate XYZ values in an RGB > RGB ICC conversion. Independently from ICC profiles, a conversion between two matrixbased colorspaces (i.e. defined by some primaries, like sRGB, AdobeRGB, Rec.2020, etc…) processes in the following way:

RGBin > XYZ(WPin) using the source colorspace primaries, therefore giving XYZ values relative to the WP of the source colorspace (for ex. D65 in the sRGB case)

chromatic adaptation from source to destination WP. If the Bradford method is used, this bold down to multiplying a 3x3 matrix that I will call WPin_to_WPout

XYZ(WPout) > RGBout using the primaries of the destination colorspace
More schematically:
XYZ(WPin) = RGBin_to_XYZ * RGBin
XYZ(WPout) = WPin_to_WPout * XYZ(WPin)
RGBout = XYZ_to_RGBout * XYZ(WPout)
Since the Bradford chromatic adaptation is a linear operation, an arbitrary number of intermediate adaptations can be inserted without changing the result (apart maybe from numerical precision considerations). Therefore, the above is equivalent to
XYZ(WPin) = RGBin_to_XYZ * RGBin
XYZ(D50) = WPin_to_D50 * XYZ(WPin)
XYZ(WPout) = D50_to_WPout * XYZ(D50)
RGBout = XYZ_to_RGBout * XYZ(WPout)
This is what an ICC conversion actually does. However, instead of explicitly doing the chromatic adaptations at each conversion, the ICC profiles store the precalculated D50adapted primaries.
This becomes more clear if we rewrite the above like this:
XYZ(D50) = WPin_to_D50 * RGBin_to_XYZ * RGBin
RGBout = XYZ_to_RGBout * D50_to_WPout * XYZ(D50)
The two matrix products WPin_to_D50 * RGBin_to_XYZ
and XYZ_to_RGBout * D50_to_WPout
are precomputed when the ICC profile is generated, thus providing the D50adapted primaries that are stored in the ICC profile instead of the original primaries found in the colorspace definition.
To be more precise, the ICC profile stores WPout_to_D50 * RGBin_to_XYZ
for the destination profile, and the resulting matrix is then inverted to provide the XYZ>RGB conversion…
Why D50? So far I could not grab wether the choice of D50 is totally arbitrary or if it has some practical consequences. Someone told me that it has to do with the fact that D50 is the reference illuminant when viewing printed images, and that ICCs where originally developed for the print industry, but I still don’t know if replacing D50 with dome other reference illuminant would change anything…
One practical consequence is that the luminance Y
that one computes from the ICC primaries is relative to D50, and so rather far from daylight conditions. I have the idea that for the B&W conversion of a daylight scene one would rather prefer to use the D65adapted luminance, but I do not have strong arguments to support my intuition…
Hope this helps!