1 Introduction

X-ray polarization offers a new window in high-energy astrophysics. Sensitive X-ray polarimetry is a promising method for revealing crucial and unique information about the structures and physical processes of black holes, neutron stars, and all classes of X-ray sources in the universe. The measurement of X-ray polarization was first successfully achieved in an X-ray polarimeter aboard the OSO-8 satellite for Bragg peak measurement [1] although only two separate energies (2.6 keV and 5.2 keV)were measured. In 2001, Costa used a high-resolution gas detector to measure X-ray polarization in the laboratory[2]. Since then, X-ray polarization measurements have become more widespread. Various particle physics and astrophysics missions for detecting the X-ray polarization, including IXPE [3] and eXTP [4], are being prepared. Other missions, such as XIPE [5]and PRAXyS [6], have been planned but subsequently canceled. Recently, CubeSat performed new X-ray polarization measurements of the Crab Nebula[7] using gas detector technology. In these missions, the polarization of energetic astrophysical sources was mainly measured through a narrow field of view using a high-sensitivity X-ray polarization detector, which cannot respond rapidly to transient sources such as gamma-ray bursts (GRBs). Based on the radiation mechanism, it is known that the magnetic field configuration of the GRB is reflected in the radiation polarization. With this in mind, it is hoped that the magnetic field configuration of the jet can be studied by measuring the polarization of the GRB radiation to understand this extreme explosion phenomenon[8]. The SXP is currently in the conceptual design stage in which a detector with a large area and a field of view is being developed. The detector design directly affects the detector sensitivity and needs to be optimized through simulations. The obtained data must be processed by a track reconstruction algorithm to obtain the polarization of the X-ray source. A good track-reconstruction algorithm can improve the accuracy of the results.

Algorithms for track reconstruction are also constantly being developed. Bellazzini et al. first proposed an algorithm for track reconstruction[9, 10]. Analysis of the track moments allows some properties of their distribution, such as the barycenter, the major axis, and most importantly, the projection on the pixel plane of the point of absorption of the X-ray photons, to be analyzed. Differing from the traditional moment analysis method based on the track image as well as an algorithm based on the shortest path problem in graph theory [11] and adaptive cut methods [12], the algorithm proposed here is based on the detector condition and is called the iteration cut method.

A brief introduction to SXP and the structure of the SXP detector unit is given in Sect. 2. The SXP Monte Carlo simulations are presented in Sect. 3. The algorithm for track reconstruction is proposed in Sect. 4, and the conclusions are presented in Sect. 5.

2 Detector structure

The SXP is composed of multiple detector arrays. They give an effective field of view of approximately 90° and a detection area of 300 cm² for detecting the polarization of transient X-ray transient sources. The SXP contains 96 detection modules. We use a single detector module to illustrate the algorithm in this study. A single detector module is shown inFig. 1. The SXP will be mounted on the Chinese space station to detect X-ray polarization within the energy range of 2 keV to 10 keV. The photoelectric effect has the largest reaction cross-section in the soft X-ray energy range, and results in the generation of photoelectrons with lower energy. Sufficiently distinct tracks can be obtained using thick-gas electron multiplier(THGEM) [13, 14] detection. THGEM has the advantages of high detection efficiency, simplicity, and stable operation. Track detection is the key to extracting the polarization of the X-ray source. The Topmetal-2^- chip is chosen for charge collection and reading out the electronic signals [15]. A field-programmable gate array (FPGA) is used to store information.

Fig. 1.Full-screen View

(Color online) Schematic diagram of a single detector module of the SXP, including the gas detector chamber, the THGEM, and the signal collection system.

The internal structure of the detector is divided into three parts, namely, the drift region, multiplier region, and induction region. As shown inFig. 2, a photon incident in a direction perpendicular to the drift region of the detector is converted into a photoelectron through the photoelectric effect in the working gas. When the photoelectron passes through the working gas, it ionizes the gas, giving rise to electron-ion pairs. The electric field causes the electrons to drift toward the multiplier region, where they are amplified through avalanche generation. The photoelectron track is finally detected via two-dimensional imaging. The initial direction of the photoelectron is related to the polarization of the incident photons. The ejection probability is the highest when the initial direction has the same angle as the polarization. The polarization probability is given by the distribution of cos²θ [16], where θ is the polar angle with respect to the incident X-ray direction.

Fig. 2.Full-screen View

(Color online) Photoelectric effect detection principle.

The THGEMs are fabricated via standard PCB mechanical drilling and global etching processes, which allow for large-scale economic production of large-size detectors [17]. When two different negative voltages are applied to the extremities of the metal layer, electrons can drift through the etched hole and avalanche multiplication can occur, thus amplifying the charge for detection. The THGEM has several advantages, including discharge prevention, sub-millimeter position resolution, large gain, and ease of manufacture and maintenance. A thin THGEM [18-20] with 200 μm thickness, 200 μm hole diameter, 500 μm hole pitch, 30 μm rim width, and 10 μm thickness was produced.

The Topmetal-2^- is a direct charge sensor that collects drift charges directly without post-processing and reads them out into a 2D projection image. The pixel array is patterned in a matrix of 72×72 pixels with an 83 μm pitch, resulting in a charge-sensitive area of 6×6 mm². Each pixel contains a low-noise charge-sensitive preamplifier (CSA) to establish an analog signal and a discriminator with a tunable threshold to generate hits. The Topmetal-2^- is characterized by its small volume with the overall geometric size of 8 mm×9 mm×0.7 mm. The average equivalent noise charge measured is less than 15 e^-. The low noise allows for improved trail detection. Furthermore, the sensor is resistant to radiation and remains functional even after 28 h of irradiation with a ⁶⁰Co source at the dose rate of 3.6 .

3 Monte Carlo simulations

The expected track was simulated using the Geant4 [21], Garfield++ [22], and COMSOL Multiphysics [23] software packages. Geant4 is a toolkit for simulating particle-matter interactions. The detector model was constructed in Geant4 and used to simulate the initial track. The structure of the SXP is shown inFig. 3. Here, photons enter the detector through a collimated aperture and are detected by the photoelectric effect in the gas. The effect of the THGEM is modeled in the following simulations.

Fig. 3.Full-screen View

(Color online) Detector structure used in Geant4. The labels in the figure denote (1) the 100 μm thick beryllium window that transmits photons and is subjected to the pressure difference between the inside and outside of the detector chamber; (2) the vacuum collimation aperture in the kovar lid, which acts as a directional guide and has an area of 1.6 mm×1.6 mm and height of 55 mm; (3) the aluminum frame layer, which is mainly used to shield low-energy electrons and background photons; (4) the ceramic, which encases the detection gas; (5) the gas required for detection; and (6) the kovar lid and base, which represent the entire shell of the detector.

The main criteria for selecting the gas are high detection efficiency and long track length. High detection efficiency is beneficial for the detection of weak sources, whereas long tracks are useful for reconstruction. As shown inFig. 4, the detection efficiency and track length of different gases with the same thickness are simulated (panel(a) and panel (b)). It can be seen that the cross section of the photoelectric effect decreases with increasing energy of the incident photons, thus reducing the detection efficiency. The efficiency at 2 keV is lower than that at 3 keV because of the Be window, which acts as a shield against low-energy photons. Another factor worth considering is the gas diffusion coefficient. The diffusion coefficient of the gas should be as small as possible. Through experiments, it was found that the diffusion coefficients of the candidate gases were all within a certain range which did not affect the reconstruction algorithm.

Fig. 4.Full-screen View

(Color online) (a) Photon detection efficiency for different photon energies. (b) Track length distribution for different energies. The simulations were carried out in Geant4 using different detection gases with the gas thickness of 1 cm. The other conditions were the same.

For the above reasons, a mix of Neon (80%) and dimethyl ether (DME 20%) was chosen as the detection gas. Ne-based gas mixtures have longer photoelectron tracks in the lower energy bands while retaining reasonable efficiency, while the DME absorbs UV photons produced by charge propagation in the gas with small transversal diffusion and relatively high gain [24]. After completing the track simulation in Geant4, Garfield++ and COMSOL Multiphysics were used to perform the track simulations in more detail. Garfield++ is a particle detector tool that is useful for detailed simulations in which a gas is used as a sensitive medium, while COMSOL Multiphysics is a multi-physical field simulation software based on finite element calculations. The gas components and the THGEM model of the detector were input into COMSOL Multiphysics to calculate the electric field distribution. The potential distribution calculated via COMSOL Multiphysics is shown inFig. 5. The low-energy physics model in Geant4 was used to generate the secondary electrons and then reconstruct them into delta electrons in Garfield++ for simulations [25]. The Topmetal-2^- chip senses the charge signal through the exposed copper on its surface; it not only generates a signal at the incident location of each electron but also in the area before impact. To simulate the actual electronic signal collection, the charge value at each pixel position in the Garfield++ simulation was redistributed between the pixel and the surrounding pixels with a Gaussian distribution, and the appropriate σ value was selected for comparison with the experiment to obtain more accurate track simulations.

Fig. 5.Full-screen View

(Color online) Distribution of cross-sectional potential in the center of the THGEM. The vertical and horizontal coordinates represent the distances. The structure of the hole is clearly visible.

Using Garfield++ simulations, which are intrinsically slow, a Gaussian distribution for the diffusion can be added to the slender electron track information produced by Geant4 to make the tracks thicker to simulate the signal sensed by the electrons on the Topmetal-2^- chip. The width of the Garfield++ track is first fitted with a Gaussian distribution to obtain the corresponding σ value. The fitted Gaussian distribution is then added to the Geant4 track to obtain a track width similar to that simulated with Garfield++. At the same time, different apertures in the THGEM directly affect the transverse diffusion of electrons, leading to different track thicknesses. The effect of different pore sizes on the track can be easily simulated by varying the σ value of the Gaussian distribution.Fig. 2 shows an example of a track obtained through simulations. The size of each histogram bin is equivalent to the size of each pixel in the Topmetal-2^-. Panel (a) shows the track used in the reconstruction. The size of the square represents the charge value at the coordinate position. The label numbers at the horizontal and vertical axes represent the detector readout positions, which coincide with the Topmetal-2^- pixels.

4 Track reconstruction algorithm

4.2 Event reconstruction

Electrons may be scattered multiple times as they pass through the gas, thus deviating from the initial ejection direction. When electrons are close to stopping, they release most of their energy in the form of a Bragg peak. Thus, only a modest charge is accumulated at the initial interaction point. In this case, monitoring only the start of the track is sufficient to achieve full reconstruction. The repeated iteration method is divided into the following four main steps, of which the first two are similar to those in the traditional method.

1) The barycenter and gravity line are calculated from the location and charge (as the weight) of the track inFig. 2(a). This line roughly marks the direction of the track trend.

2) A line represented by the green line inFig. 2(a) is traced at the barycenter perpendicular to the gravity line. Because Bragg peaks are accumulated at the tail of the track, the head and tail can be separated using third-order moment calculations to determine the half-zone where the photoelectron interaction point is located. The third-order moment M_d is calculated as

$M_{\text{d}}={\displaystyle {\sum }_{\text{pixels}}}(d^3\times q)\mathrm{,}$ (1)

where d is the distance between the pixel and the line perpendicular to the gravity line (d may be negative), q is the charge of the pixel, and ∑_pixels indicates summation over all the pixels. Thus, by comparing the d and M_d of each pixel with similar or dissimilar signs, it can be determined whether the pixel is within the region where the photoelectron Bragg peak is located, or in the region where the photoelectrons are produced, which is indicated by the semicircle in the figure. Tracks with M_d = 0 are symmetrical in shape and difficult to reconstruct; thus, such tracks are rejected.

3) Multiple lines perpendicular to the gravity line are traced to re-divide the track, as shown by the multiple green lines inFig. 2(a). The length of the track determines the direction of the subdivision, and for excessively long tracks, the track is subdivided along the direction that shortens the starting half of the track. The track length is calculated as $M_2^{\text{min}}/M_2^{\text{max}}$. Each point on the track has the distance of $M_2^{\text{min}}$ from the gravity line and the distance of $M_2^{\text{max}}$ from the line perpendicular to the gravity line. The sum of the squares of $M_2^{\text{min}}$ and $M_2^{\text{max}}$ represents the square of the distance of the point from the barycenter M₂. M₂ is given by

$M_2={\displaystyle {\sum }_n^{n=\text{pixels}}}{(x_n-x_0)}^2+{(y_n-y_0)}^2\mathrm{,}$ (2)

where x_n and y_n represent the coordinates of each pixel, and x₀ and y₀ are the coordinates of the barycenter. The track length value is calculated at the best applicable energy point of the traditional algorithm, which is approximately a Gaussian distribution, and the mean value of the Gaussian distribution is used as the initial threshold for determining if the track length is long or short. If the track length is smaller than this value, the track is a short track. Otherwise, it is a long track. At the same time, the mean value of the Gaussian distribution is used as one of the parameters that can be adjusted later. Furthermore, for tracks that are too short, a re-division needs to be performed along a direction that makes the half-zone of the starting track longer. The spacing between the green lines is set as a modifiable parameter. The degree of curvature of the track can be determined by calculating ${\displaystyle {\sum }_{\text{pixels}}}{M}_2^{\text{min}}$. Unlike the third-order moment calculation, the influence of the charge value is removed in the second-order moment calculation and only the shape of the track is considered, thus preventing the uneven charge distribution of the track from affecting the results. The variation in the second-order moment is expressed by

$\left|{\displaystyle \sum _{\text{out}}{M}_2^{\text{min}}-}{\displaystyle \sum _{\text{last}}{M}_2^{\text{min}}}\right|/{\displaystyle \sum _{\text{pixels}}{M}_2^{\text{min}}\mathrm{,}}$ (3)

where ${\displaystyle {\sum }_{\text{cut}}}{M}_2^{\text{min}}$ is the summation of the second-order moments corresponding to the new green line after the move, and ${\displaystyle {\sum }_{\text{last}}}{M}_2^{\text{min}}$ is the second-order moment summation obtained after the last move. The reason for subtracting the two summations and taking the absolute value is because the result represents the amount of change in the second-order moment, and only the change in the trail needs to be detected. Each time the green line is moved, the point of charge at the other end of the green line is not involved in the calculation of the new second-order moment summation. The initial second-order moment sum of the first move is therefore calculated as (${\displaystyle {\sum }_{\text{pixels}}}{M}_2^{\text{min}}$). The variation is summed over the initial second-order moment to normalize the variation of the second-order moments for different tracks. The second-order moments of the green lines are computed in sequence, and the final split position chosen is the position with the most dramatic change. The above equation provides the maximum value at which a split will occur. To prevent the green line from dividing the outside of the track, a minimum track head length was set to preserve the first part of the track and the green line was forbidden from crossing this position. The minimum track head length is also a variable that can be adjusted. The actual parameters used were obtained after several rounds of empirical testing.

4) The barycenter and center of gravity line are calculated in the remaining half of the partition. The center of gravity is identified in the diagram by the interaction point (marked by an asterisk inFig. 6(a))generated by the algorithm. The remaining track gravity line is illustrated in red in the diagram, while the arrow indicates the polarization direction obtained via the algorithm because the initial direction of photoelectron emission is also the polarization angle.

Fig. 6.Full-screen View

(Color online) A detailed diagram of the reconstruction process (a) and a typical track image (b). This track was produced by a 8 photon in neon (80%) mixed with dimethyl ether (DME 20%) at 0.8 atmosphere. The horizontal and vertical coordinates correspond to pixels with the same size as the actual Topmetal-2^- pixels. In b, the photoelectron hits in the figure correspond to the actual positions of the chip pixels. The details of Fig. 6(a) are described in Section 4.

Figure 6 shows the tracks obtained from the simulation and the reconstruction process. Some of these parameters can be modified at the beginning of the reconstruction, including the pixel blocks used for filtering circular tracks, and the value of the minimum preserved track head length. These parameters have a direct impact on all aspects of the reconstruction, including the selection of pre-reconstruction events. Monte Carlo simulations were used to generate a certain number of photoelectron tracks. The above algorithm was then used to analyze the tracks line by line. The initial emission angle of each track was obtained through parameter optimization. The overall angle distribution is shown inFig. 7, where the horizontal coordinate represents the angle and the vertical coordinate represents the number of instances.

Fig. 7.Full-screen View

(Color online) Superposition distribution of track reconstruction angle. The selected photons were 100% polarized with the energy of 8 keV. The photons entered the detector vertically with the angles θ and φ fixed to π/4 and π/2, respectively. θ is the azimuth angle of the X-ray electric vector, and θ is the polar angle with respect to the incident X-ray direction.

4.3 Event fitting analysis

The angular distribution of the K-shell photoelectron emission in the non-relativistic region is theoretically given by [26]:

$\frac{\text{d}{\sigma }_{\text{ph}}^K}{\text{d}\Omega }\propto \frac{{{\displaystyle \sin }}^2\varphi {{\displaystyle \cos }}^2\theta }{{(1-\beta \cos \varphi )}^4}\mathrm{,}$ (4)

where β is the ratio of the photoelectron velocity to the speed of light, θ is the azimuth angle of the X-ray electric vector, and φ is the polar angle with respect to the incident X-ray direction.

The corresponding equation for the response amplitude after projecting the track onto a two-dimensional plane for processing is

$M(\theta )=A\times {{\displaystyle \cos }}^2(\theta -B)+C\mathrm{.}$ (5)

The emission direction of polarized photoelectrons obeys such a distribution in theory. This equation was also used to fit the data shown inFig. 7. In practice, modulation factors are often used to describe the polarization of the resulting photons because of slight deviations due to limited detection capability and the precision of the algorithm itself. The modulation factor is calculated as

$f_{\text{eff}}=\frac{M_{\text{max}}-M_{\text{min}}}{M_{\text{max}}+M_{\text{min}}}\mathrm{,}$ (6)

where M_max and M_min correspond to the maximum and minimum values of the angle distribution inFig. 7. The track reconstruction efficiency and the corresponding pre-screening efficiency of the SXP are shown inFig. 8. The low-energy portion is difficult to reconstruct because the track lengths are too short. The middle energy portion is the most suitable area for the algorithm. At higher energies and longer track lengths, the modulation factor decreases because the parameters of the algorithm need to be adjusted. The circles show the efficiency of the pre-selection, and the ordinate the proportion of tracks retained for reconstruction. It can be observed that only a few cases have been screened out, and that the vast majority have been retained for reconstruction.

Fig. 8.Full-screen View

(Color online) Modulation factor as a function of photon energy. The selected photons were 100

Based on the conditions set in the simulations, a 4.5 keV X-ray polarization source was used in the experiment to reconstruct the measurements. Figure 9 shows the actual structure of the detector. The upper part of the structure is the SXP detector, while the center is the THGEM. The bottom part is the Topmetal-2^-. Both sides were equipped with gas pipes for ventilation and power supply connection wires for the application of high voltages. The signal collected by the Topmetal-2^- is shown inFig. 10.

Fig. 9.Full-screen View

(Color online) The actual detector in the laboratory. The upper part is the detector chamber and the THGEM, and the bottom part is the Topmetal-2^-.

Fig. 10.Full-screen View

(Color online) Experimentally obtained tracks. (a) The THGEM structures were used to multiply the resulting track. Some background noise was present. (b) The tracks were processed with an image erosion and dilation algorithm for reconstruction, which almost completely removed the effect of the background.

The experimental tracks were reconstructed, and their angular distributions are shown inFig. 11. The experimental tracks yielded the modulation factor of ∼ 28.60 ± 0.86%, while the event rejection rate was 46.68%. The modulation factor of the unpolarized X-ray source was 2:057%. The tracks obtained from the THGEM experiment were very rough and poorly reconstructed because of the large pore size. There were very few tracks that could be easily reconstructed, as shown inFig. 10. In addition to the event selection rules described in Section 4.1 (the range for the circular track judgment was nine pixels), there were tracks distributed at the edge of the image in the experiment that were not involved in the reconstruction. For a track to be reconstructed, the track must have a sufficient number of pixels in the central area. Both the size of the central area and the number of pixels used as the threshold are adjustable input parameters. The main reason for the inconsistency between the experimental and simulation results is the influence of the experimental background; for example, the X-ray sources were not 100% polarized, and readout noise or dark current were present. Meanwhile, the spatial resolution was insufficient, resulting in possible inaccuracies in the reconstructed photoelectron emission direction. The spatial resolution will be optimized in future studies.

Fig. 11.Full-screen View

(Color online) Angular distribution of experimentally obtained tracks after reconstruction. A 4.5 keV X-ray polarization source was used in the left figure. A 4.5 keV unpolarized X-ray source was used in the right figure.

5 Discussion and Summary

The SXP adopts the most popular soft X-ray polarization detection method at present. The appropriate experimental gas and structure were selected through experiments and simulations. The iteration cut algorithm can provide reconstruction angles that are important for measuring the polarization and relevant for detection sources in the 210 energy range. The modulation factor reaches as high as 57 % above 7 . At low energies, the modulation factor is relatively low because of the short photoelectron tracks and the effect of electron diffusion in the THGEM avalanche. Experiments have shown that reducing the hole size of the THGEM structure can yield better results. The next step will be to optimize the hardware structure of the SXP and the algorithm, for example by using deep learning. It is envisaged that high-precision X-ray polarization detection will have a significant impact in revealing the physical mechanisms of astrophysical sources in the future.