The simulation of light based on physical principles as done in path tracing algorithms combined with the performance of modern computers make it an interesting solution to light propagation problems, which are difficult if not impossible to solve numerically. Calculation of the expected image produced by the X-ray telescope for a realistic X-ray telescope and solar axion emission is one such problem. In (Christoph Krieger 2018) Christoph Krieger already wrote a basic hardcoded ⁴ raytracer for the 2014/15 data taking behind the ABRIXAS telescope. He approximated the ABRIXAS telescope as a simple lens.

For the data taking with the detector of this thesis behind the LLNL telescope, a more precise solution was desired. As such Johanna von Oy wrote a new raytracer for this purpose as part of her master thesis (Von Oy 2020). However, although this code (Von Oy 2023) correctly simulates the LLNL (and as a matter of fact also other telescopes) based on their correct layer descriptions, it still uses a similar hardcoded ⁴ approach, making it inflexible for different scenes. More importantly it makes understanding and verifying the code for correctness very difficult.

What started mainly as a toy project I wrote out of curiosity following the popular 'Ray Tracing in One Weekend' (RTOW) (Shirley 2020) book ⁵, anyhow made me appreciate the advantage of a generic raytracer (see fig. 1(a) for the final result of that book, created by this project). Generic meaning that the scene to be rendered is described by objects placed in simulated 3D world space. This reduces the amount of code that needs to be checked for correctness by orders of magnitude (reflection for one material type is a few line function, instead of having lines of code for every single reflection in a hardcoded raytracer). It also makes the program much more flexible to support different setups. In particular though, the concept of a camera independent of the scene geometry allows to visualize the scene (obviously the entire purpose of a normal raytracer!).

As a personal addition to the RTOW project, I implemented rendering to an SDL2 (Schmidt 2022b) window and handling mouse and keyboard inputs to produce an interactive real-time raytracer (this was early 2021). This sparked the idea to use the same approach for our axion raytracer. I pitched the idea to Johanna, but did not manage to convince her ⁶. While working on finalizing the axion image and X-ray finger simulations for this thesis, several big discrepancies between our raytracing results and those used for the CAST Nature paper in 2017 (Collaboration and others 2017) were noticed. The latter were done by Michael Pivovaroff at LLNL using a private in-house raytracer ⁷. The need to better verify the correctness of our own raytracing simulation lead me down the path of finishing the project of turning the RTOW based raytracer into an axion X-ray raytracer. This way verifying correctness was easier. Some additional features of the second book of the ROTW series were added (light sources among others) and other aspects were inspired by Physically Based Rendering (Pharr, Jakob, and Humphreys 2016) (in particular the idea of propagating multiple energies in each ray).

The result is TrAXer (Schmidt 2023), an interactive real-time visible light and X-ray raytracer. In essence it is two raytracers in one. First, a 'normal' visible light raytracer, which the camera uses. Secondly, an X-ray raytracer. Both use the same raytracing logic, but they differ in sampling. For visible light we emit rays from the camera into the scene, while for X-rays we emit from an X-ray source (an object placed in the scene, for example the Sun). Both of these run in parallel. Scenes are defined by placing geometric primitives with different materials (glass, metal etc.) into the scene.

The X-ray raytracer is used by choosing X-ray emitting materials and adding another object as an X-ray target (another material), which the X-ray sampling will be biased towards. There is further a specific X-ray material an object can be made of, which behaves as expected when hit by an X-ray, calculated using (Schmidt and SciNim contributors 2022). Behavior is described by the reflectivity computed via the Fresnel equations as discussed in sec. 6.1.2. To simulate reflectivity of the LLNL telescope, the depth graded multilayers of the different layers are taken into account. X-ray transmission is currently not implemented, but only because there is no need for us. The addition would be easy based on mass attenuation coefficients and the Fresnel equations for transmission.

In order to detect and read out results from the X-ray raytracing, we place "image sensors" into the scene. They use a special material, which accumulates X-rays hitting them, split into spatial pixels. For visible light they simply emit the current values stored in each pixel of the image sensor, mapped to the Viridis color scale.

These things are visible in the example screenshot showing the "CAST" scene in fig. 1(b). The view is seen from behind and above the setup with a narrow field of view (like a telephoto lens). Only relevant pieces of the setup are placed into the scene (partially for performance, although bounding volume hierarchies help there, mostly for simplicity). The Sun is visible on the right side, emitting a yellow hue for the visible rays and an invisible X-ray spectrum following the expected solar axion flux (with correct radial emission). The long blueish tube is the inner bore of the CAST magnet. Then, with red hues (in visible light) we have the LLNL telescope. In the bottom left is the image sensor, which has the same physical size as the center GridPix of the Septemboard detector and is placed in the focal point of the telescope. As every object of the scene is floating above the earth, combined with the telephoto like view, the sensor looks a bit odd. ⁸

The Lawrence Livermore National Laboratory, introduced in sec. 5.1.4, is the X-ray telescope of interest for us. Let us now first get an overview of the telescope and its construction itself, before we move on to comparing the TrAXer simulations with PANTER measurements and raytracing simulations from Michael Pivovaroff from LLNL.

Its design follows the technology developed for the NuSTAR mission (Harrison et al. 2013, 2006; Koglin et al. 2009; Craig et al. 2011; Hailey et al. 2010). This means it is not a true Wolter optic (meaning parabolic and hyperbolic mirrors), but instead uses a cone approximation (Petre and Serlemitsos 1985), which are easier to produce and in particular for CAST provide more than enough angular resolution.

It consists of two sets of 13 layers. Each mirror has a physical length of \(\SI{225}{mm}\) and thickness of \(d_{\text{glass}} = \SI{0.21}{mm}\). These mirrors are made of glass, which are thermally-formed ("slumped") into the conical shapes. Production of the mirrors happened at the National Space Institute at the Technical University of Denmark (DTU Space) as part of the PhD thesis of Anders Jakobsen (Jakobsen 2015). The first sets of layers \(i\) describe truncated coni with a remaining height of \(\SI{225}{mm}\), an arc angle of \(\SI{30}{°}\), a radius on the opening end ⁹ of \(R_{1,i}\) and a cone angle of \(α_i\). That is for each layer a non-truncated cone would have a total height of \(h = R_{1,i} / \tan{α_i}\). The second set of 13 layers only differs by their radii, \(R_{4,i}\) on the opening end and an angle of \(3 α_i\). In the horizontal direction along the optical axis of the telescope, there is a spacing of \(x_{\text{sep}} = \SI{4}{mm}\) between each set of layers. Due to the tilted angle the physical distance between each set is minutely larger. Combined the telescope thus has a length of roughly \(\SI{454}{mm}\).

One can define additional helpful radii, \(R_{2,i}\) the minimum radius of the first set of cones, \(R_{4,i}\) the maximum radii of the opening end of the second set of mirrors and finally \(R_{3,i}\), the radius at the mid point where the hypothetical extension of the first layer and last layer would touch at \(x = \SI{2}{mm}\) from each set of mirrors. Equations \eqref{eq:llnl_rest:radii_eqs} specify how these numbers are related with \(l\) the length of the mirrors.

which we can rewrite to compute everything from \(R_{1,i}\) by:

I emphasize this, because this conical nature describes the entire telescope as long as the radii \(R_{1,i}\), \(R_{4,i}\) and angle \(α_i\) are known. ¹⁰

Alternatively, this can also be used to construct an optimal conical telescope from a starting radius \(R_{1,0}\) (typically called the "mandrel"), if the iterative condition \(R_{3,i+1} = R_{1,i} + d_{\text{glass}}\) is employed (in this case some target parameters are needed of course, like focal length, \(x_{\text{sep}}\) and a few more).

To calculate the focal length \(f\) of an X-ray telescope of the Wolter type, we can use the Wolter equation,

\[ \tan(4α_i) = \frac{R_{3,i,\text{virtual}}}{f} \]

where \(α_i\) is the angle of the \(i^{\text{th}}\) layer and \(R_{3,i,\text{virtual}}\) is the radius corresponding to the virtual reflection point. Due to the double reflection inherent to a Wolter telescope, the incoming ray picks up a total angle of \(4α\) (\(0\) to \(2α\) after first reflection of the shell with angle \(α\), then \(2α\) to \(4α\) after reflecting from second shell of angle \(3α\)). The virtual radius is the radius at the point if one extends the incoming ray through the first set of mirrors and reflects at a single virtual mirror of an angle \(4α\) precisely at the midpoint of the telescope. This is illustrated in fig. 3, which is an adapted version of fig. 4.6 from A. Jakobsen's PhD thesis (Jakobsen 2015) to better illustrate this.

In addition to the PhD thesis mentioned, there is a paper (Aznar et al. 2015) about the telescope for CAST and initial results. Note though, that both the PhD thesis and the paper contain contradicting and wrong information about the details of the telescope design. Neither (nor combined do they) describe the real telescope. ^{11 , 12} Thanks to the help of and personal communication with Jaime Ruz and Julia Vogel we were able to both clear up confusions and find the original numbers of the telescope that were actually built. Table 37 gives an overview of the telescope design. It is adapted from (Aznar et al. 2015), but modified to use the correct numbers. Then in tab. 38 is the list of all layers and their corresponding radii in detail as given to me by Jaime Ruz.

Having discussed the physical layout of the telescope, let us quickly talk about the telescope coatings. In contrast to telescope like the ABRIXAS or XMM-Newton optics, which use a single layer gold coating, the NuSTAR design uses a depth graded multilayer coating as introduced in sec. 6.1.2. There are four different 'recipes', used for different layers. All recipes are a depth graded multilayer of different number of Pt/C layers.

Note that in this terminology the high \(Z\) material (platinum) is actually below the low \(Z\) material (carbon) contrary to what Pt/C might imply!

Table 39 gives an overview of which recipes are used for which layer and how they differ. \(d_{\text{min}}\) is the minimum thickness of one multilayer and \(d_{\text{max}}\) the maximum thickness. The top most multilayer is the thickest. \(Γ\) is the ratio of the top to bottom material in each multilayer.

For example recipe 1 has a top multilayer of carbon of a thickness \(Γ · d_{\text{max}} = \SI{10.125}{nm}\) on top of \((1 - Γ) · d_{\text{max}} = \SI{12.375}{nm}\) of platinum. Because recipe 1 only has 2 such multilayers, the bottom Pt/C layer has a combined thickness of \(d_{\text{min}}\).

To reiterate the equations that describe the layer thicknesses for \(N\) layers as given in sec. 6.1.2, a depth-graded multilayer is described by the equation:

where \(d_i\) is the depth of layer \(i\) (out of \(N\) layers),

\[ a = d_{\text{min}} (b + N)^c \]

and

\[ b = \frac{1 - N k}{k - 1} \]

with

\[ k = \left(\frac{d_{\text{min}}}{d_{\text{max}}}\right)^{\frac{1}{c}}. \]

In all four recipes used for the LLNL telescope the parameter \(c\) is set to 1.

Calculation of the reflectivity follows the equations also mentioned in the theory section. The reflectivities of each layer are calculated with (Schmidt and SciNim contributors 2022), which is used as part of the material description in TrAXer.

The LLNL telescope was tested and characterized at the PANTER X-ray test facility in Munich in July 2016. It was reported on these tests during the \(62^{\text{nd}}\) CAST collaboration meeting (CCM) on <2016-09-26 Mon> and again in the CCM on <2017-01-23 Mon>. ¹³ These results are valuable to us, because they can be used as a cross-check to verify our TrAXer based raytracing simulations.

At the PANTER facility the X-ray source sits \(\SI{130.297}{m}\) away from the center of the telescope (the plane of the virtual reflection). It can be approximated as a source of radius \(\SI{0.42}{mm}\). The detector used to measure at PANTER is installed on the optical axis (i.e. directly 'in front of' the source). The telescope itself is offset by the radius of its shells, such that its optical axis aligns with the optical axis defined by the X-ray beam (instead of aligning the telescope entrance with the beam axis and the detector offset from it!). In this way the light entering the telescope does not enter perpendicular to the opening, but under an angle.

The setup of having such an effective point source implies that the best focal point will not be in the physical focal spot defined by the focal length, but slightly behind. At PANTER this was measured to be at \(\SI{1519}{mm}\) instead of \(\SI{1500}{mm}\) and provides a first check of our raytracer. When reproducing the same setup in TrAXer we find the smallest, most symmetric spot around the \(\SIrange{1519}{1520}{mm}\) mark as well. ¹⁴

In this setup three different X-ray fluorescence lines were measured. \(\ce{Al}\) Kα (\(\SI{1.49}{keV}\)), \(\ce{Ti}\) Kα (\(\SI{4.51}{keV}\)) and \(\ce{Fe}\) Kα (\(\SI{6.41}{keV}\)). We can compare the images between measurements and simulations as well as compute the 'half power diameter' (HPD). The HPD is the radius around the image center (based on the weighted mean center) in which \(\SI{50}{\%}\) of the flux is contained. This is computed based on the encircled energy function (EEF): integrate the flux in radial slices around the center, compute a radial cumulative distribution function and find the radius corresponding to \(\SI{50}{\%}\).

Let's consider the aluminum Kα lines. Fig. 4 compares the PANTER data, with the LLNL raytracing result and our raytracing result. ¹⁵ All three plots show the inner \(\num{3}·\SI{3}{mm²}\) of the image. Generally, the agreement is very good. The "bow tie" shape can be seen in all three figures.

Computing the encircled energy function yields a figure as shown in fig. 5. We see the HPD shown as the red line at a radius of \(r_{\text{HPD}} \approx \SI{0.79}{mm}\). This is converted to an HPD in arc seconds via

\[ α_{\text{HPD}} = \arctan\left( \frac{r_{\text{HPD}}}{f} \right) \]

where \(f\) is the focal length of the optic. This leads to an HPD of \(\SI{216.38}{\arcsecond}\), which is slightly above both the LLNL raytracing simulation (PANTER model) and the PANTER data. Doing this for all three fluorescence lines, we get the numbers shown in tab. 40, where the differences are even smaller. Therefore, our implementation of the raytracer including figure errors yields very compatible results to the LLNL model developed for PANTER. The usage of TrAXer as the tool of choice for the calculation of the axion image seems justified.

See sec. 37.3.1 below for more details on how the figure error is introduced and its parameters determined.

The previous section provides us with reasonable guarantees that our raytracer produces compatible results to PANTER measurements and the LLNL raytracing results. In particular for the axion-electron coupling we do not wish to reuse the raytracing prediction for the solar axion image computed by Michael Pivovaroff at LLNL, as used in the CAST Nature 2017 paper (Collaboration and others 2017) paper. ¹⁸ This is because of the different radial emission profile (compare fig. #fig:theory:solar_axion_flux:radial_dependence_ksvz_dfsz). And for the chameleon, due to its production in the solar tachocline, this is of course absolutely mandatory.

To simulate the solar axion image, we use the scene already shown in fig. 1(b). The Sun is placed as an X-ray emitter ¹⁹ at \(\SI{0.989}{AU}\) away ²⁰ from the telescope. It emits X-rays of energies and from radii corresponding to the distribution seen in fig. #fig:theory:solar_axion_flux:flux_vs_energy_and_radius calculated using the opacity calculation code (Von Oy 2023) also developed by Johanna von Oy during her master thesis (Von Oy 2020). The calculations are based on the AGSS09 solar model (Asplund et al. 2009; Serenelli et al. 2009). The telescope is rotated under the angle it was installed at CAST, as deduced from the X-ray finger taken at CAST, fig. #fig:cast:xray_finger_centers (about \(\SI{14}{°}\)). The image sensor is placed slightly in front of the focal point. \(\SI{1492.93}{mm}\) away from the center of the telescope (to the virtual reflection point) instead of \(\SI{1500}{mm}\) to account for the median conversion point of the X-rays of the expected axion electron X-ray flux (as mentioned in sec. 13.10.5 and calculated in 36).

With all these in place, the axion image comes out as seen in fig. 39(a). In comparison fig. 39(b) shows the axion image as computed by the LLNL raytracing code (albeit only for the axion-photon coupling and related emission). The similarity between the two images is implying our raytracing simulation is sensible. Small differences are expected, not only because of a different solar emission. The different rotation angle is due to a different assumed rotation of the telescope. Our raytracing image shows a larger bulge in the center at very low intensity. The absolute size of the major flux region is very comparable. The slight asymmetry in our image is due to it being effectively \(\SI{8}{mm}\) in front of the focal point.

We can of course also attempt to reproduce an X-ray finger run using TrAXer. Essentially it is not too dissimilar from the PANTER measurements, with two differences:

1. source distance and size
2. source placement. In contrast to the PANTER measurement, where the source is on the optical axis, at CAST for X-ray finger measurements the source is in the magnet bore and therefore in front of the telescope (i.e. away from the optical axis).

The main reason this is not part of the main thesis is that this is neither very important and at the same time we don't actually know the dimensions or distance of the X-ray finger used at CAST very well.

-> Decide final arguments

./raytracer \
    --width 1200 --speed 10.0 --nJobs 32 --vfov 45 --maxDepth 10 \
    --llnl --focalPoint \
    --sourceKind skXrayFinger \
    --sensorKind sSum \
    --energyMin 2.5 --energyMax 3.5 \
    --usePerfectMirror=false \
    --sourceDistance 10.m \
    --sourceRadius 2.0.mm 

[INFO] Writing buffers to binary files. [INFO] Writing file: out/buffer₂₀₂₃-11-12T17:16:36+01:00_{typeuint32len1440000width1200height1200.dat} [INFO] Writing file: out/counts₂₀₂₃-11-12T17:16:36+01:00_{typeintlen1440000width1200height1200.dat} [INFO] Writing file: out/image_sensor02023-11-12T17:16:36+01:00__{dx14.0dy14.0dz0.1typefloatlen1000000width1000height1000.dat}

Note the lack of --sourceOnOpticalAxis, which puts the source behind in at the center of the magnet bore instead of the optical axis.

./plotBinary \
    --dtype float \
    -f out/image_sensor_0_2023-11-12T17:16:36+01:00__dx_14.0_dy_14.0_dz_0.1_type_float_len_1000000_width_1000_height_1000.dat \
    --invertY \
    --out ~/phd/Figs/raytracing/xray_finger_10m_2mm_3keV.pdf \
    --inPixels=false \
    --title "X-ray finger in 10 m distance, 2 mm radius at 2.5-3.5 keV"

And now using 14.2 m and 3 mm source, which is what is mentioned in (Jakobsen 2015) under fig. 4.32. While I don't understand how these numbers are supposed to make sense, they produce a better image.

./raytracer \
    --width 1200 --speed 10.0 --nJobs 32 --vfov 45 --maxDepth 10 \
    --llnl --focalPoint \
    --sourceKind skXrayFinger \
    --sensorKind sSum \
    --energyMin 2.5 --energyMax 3.5 \
    --usePerfectMirror=false \
    --sourceDistance 14.2.m \
    --sourceRadius 3.0.mm 

[INFO] Writing buffers to binary files. [INFO] Writing file: out/buffer₂₀₂₃-11-12T17:21:37+01:00_{typeuint32len1440000width1200height1200.dat} [INFO] Writing file: out/counts₂₀₂₃-11-12T17:21:37+01:00_{typeintlen1440000width1200height1200.dat} [INFO] Writing file: out/image_sensor02023-11-12T17:21:37+01:00__{dx14.0dy14.0dz0.1typefloatlen1000000width1000height1000.dat}

./plotBinary \
    --dtype float \
    -f out/image_sensor_0_2023-11-12T17:21:37+01:00__dx_14.0_dy_14.0_dz_0.1_type_float_len_1000000_width_1000_height_1000.dat \
    --invertY \
    --out ~/phd/Figs/raytracing/xray_finger_14.2m_3mm_3keV.pdf \
    --inPixels=false \
    --title "X-ray finger in 14.2 m distance, 3 mm radius at 2.5-3.5 keV"

Yielding (10 m) and (14.2 m):

I still believe it is a bit weird as to why the 10 m case is definitely too large. Comparing our result to (from ./../org/Doc/StatusAndProgress.html)

also looks a bit smaller than ours.

But well, maybe there is still something wrong, but I think this will be for someone else to understand. :/

In any case, comparing the 14.2 m simulation to our real data: shows a near perfect match! So all in all this is a success!

Need:

• light sources (4h of work at most)
• cylinders, hyperboloids, paraboloids as objects (once we figure out one, the rest should be relatively easy)
• placing different telescope layers etc. (2h)

In theory this should be possible as an extensive weekend project! I'd say this is definitely worth it.

<2023-11-02 Thu 19:11>, damn I cannot stress how good it felt to just turn the old TODO there into a DONE! When I wrote the above lines, I didn't think I'd end up finishing that raytracer, to be honest. At the same time I also didn't think the thesis would still be worked on at the end of 2023! "Extensive weekend project" – yeah. I mean my estimates weren't that far off (within π I guess), but the additional work to get everything working correctly, implementing the actual X-ray raytracer on top of the visible light based raytracer using spectral radiance and not an RGB based approach, thinking and implementing how to essentially have 2 raytracers in one with ImageSensors as a replacement for the camera for X-rays etc. In the end it was a good 4 weeks of work I guess instead of a weekend. :)

I moved the "light target" to the end of the magnet (on the telescope side) in order to get better statistics for regular X-ray finger runs. While doing this I realized that with --sourceOnOpticalAxis the target at the beginning is problematic, because the telescope may not be fully illuminated at all! But the below shows that it does not make an important difference.

UPDATE <2023-11-10 Fri 21:28>: I forgot to add --ignoreMagnet to all commands above! Need to rerun them tomorrow!!!

• [X] FOR THAT adjust the light target size, make it a parameter! -> Done.

./raytracer \
    --width 1200 --speed 10.0 --nJobs 32 --vfov 45 --maxDepth 10 \
    --llnl --focalPoint \
    --sourceKind skXrayFinger \
    --rayAt 1.013 \
    --sensorKind sSum \
    --energyMin 1.48 --energyMax 1.50 \
    --usePerfectMirror=false \
    --ignoreWindow \
    --sourceDistance 130.297.m \
    --sourceRadius 0.42.mm \
    --telescopeRotation 90.0 \
    --sourceOnOpticalAxis \
    --ignoreMagnet \
    --targetRadius 40.mm

[INFO] Writing buffers to binary files. [INFO] Writing file: out/buffer₂₀₂₃-11-11T09:26:25+01:00_{typeuint32len1440000width1200height1200.dat} [INFO] Writing file: out/counts₂₀₂₃-11-11T09:26:25+01:00_{typeintlen1440000width1200height1200.dat} [INFO] Writing file: out/image_sensor02023-11-11T09:26:25+01:00__{dx14.0dy14.0dz0.1typefloatlen1000000width1000height1000.dat}

Once you deem that enough statistics has been accumulated, press F5 on your keyboard to save all buffers to binary files. The output filenames will be printed to the terminal.

And then with the produced binary files we can plot them and compute the HPDs:

./plotBinary \
    --dtype float \
    -f out/image_sensor_0_2023-11-11T09:26:25+01:00__dx_14.0_dy_14.0_dz_0.1_type_float_len_1000000_width_1000_height_1000.dat \
    --invertY \
    --out /tmp/test_me.pdf \
    --inPixels=false \
    --title "Al Kα, 1.49 keV, 3x3 mm, figure errors, source optical axis" \
    --xrange 1.5

For my development notes taken when I fixed the implementation to what it is now in the above, see ./../org/Doc/StatusAndProgress.html.