#  An Ensemble Machine Learning Approach for Screening Covid-19 based on Urine
Parameters

Behzad Moayedi<EMAIL_ADDRESS>Abdalsamad Keramatfar<EMAIL_ADDRESS>Mohammad Hadi Goldani<EMAIL_ADDRESS>Mohammad Javad Fallahi
<EMAIL_ADDRESS>Alborz Jahangirisisakht<EMAIL_ADDRESS>Mohammad Saboori<EMAIL_ADDRESS>Leyla badiei<EMAIL_ADDRESS>Computer
Engineering Department, Amirkabir University of Technology, Tehran,Iran
Research Center for Development of Advanced Technologies, Tehran, iran  Shiraz
University of Medical Sciences, Shiraz,Iran  Yasuj University of Medical
Sciences, Yasuj, Iran

###### Abstract

The rapid spread of COVID-19 and the emergence of new variants underscore the
importance of effective screening measures. Rapid diagnosis and subsequent
quarantine of infected individuals can prevent further spread of the virus in
society. While PCR tests are the gold standard for COVID-19 diagnosis, they
are costly and time-consuming. In contrast, urine test strips are an
inexpensive, non-invasive, and rapidly obtainable screening method that can
provide important information about a patient’s health status. In this study,
we collected a new dataset and used the RGB (Red Green Blue) color space of
urine test strips parameters to detect the health status of individuals. To
improve the accuracy of our model, we converted the RGB space to 10 additional
color spaces. After evaluating four different machine learning models, we
proposed a new ensemble model based on a multi-layer perceptron neural
network. Although the initial results were not strong, we were able to improve
the model’s screening performance for COVID-19 by removing uncertain regions
of the model space. Ultimately, our model achieved a screening accuracy of 80%
based on urine parameters. Our results suggest that urine test strips can be a
useful tool for COVID-19 screening, particularly in resource-constrained
settings where PCR testing may not be feasible. Further research is needed to
validate our findings and explore the potential role of urine test strips in
COVID-19 diagnosis and management.

###### keywords:

COVID-19 , Machine Learning , Urine test , multi-layer perceptron neural
network

††journal: journal

##  1 Introduction

Severe acute respiratory syndrome coronavirus type 2 (SARS-CoV-2), which
causes coronavirus disease 2019 (COVID-19), has infected more than 550 million
people and caused more than 6.3 million deaths worldwide (as of July 6, 2022,
according to WHO data) [Announcement, 2022]. Identification of risk factors
and prognostic markers for the development of COVID-19 disease and its
consequences is urgently needed to enable early identification and follow-up
of high-risk patients. In this regard, researchers around the world tried to
find diagnosis methods based on different indicators of health including chest
scan and voice quality[Roberts et al., 2021, Kovács et al., 2021, Asiaee et
al., 2020, Despotovic et al., 2021]

Recently, biomarker studies have been encouraged around the world to diagnose
COVID-19 at an early stage, aiming to target treatment more quickly and reduce
health problems associated with the disease[Banerjee et al., 2020, Wu et al.,
2020].

A dipstick or urine test strip is an indispensable tool in medical diagnostics
to gain a rapid survey of the state of health of the patient[Yetisen et al.,
2013]. A standard urine test strip can consist of at least ten different
buffers or chemical reagents that react (with color change) when immersed in
and then removed from a urine sample. After immersion, the test can often be
read in as little as 60 to 120 seconds, although specific tests require a
longer time. Routine urinalysis with multiparameter strips is the first step
in diagnosing a wide range of diseases. The analysis includes tests for the
presence of proteins, glucose, ketones, hemoglobin, bilirubin, urobilinogen,
acetone, ascorbic acid, nitrite, and leukocytes and tests for pH and specific
gravity or tests for infection with different pathogens[Strasinger & Di
Lorenzo, 2010]. Due to the ease of performing the urine test and its non-
invasiveness, anyone without the need for special instructions can use it.
Bonetti et al. [2020], Chen et al. [2020], de Fátima Cobre et al. [2021], Liu
et al. [2020], Erdogan et al. [2021] showed that urinalysis parameters can be
used for predicting severity and detection of coronavirus disease 2019
(COVID-19). Due to the availability, simplicity, and the low cost of preparing
urine strips ($55 for a pack contains 100 strips) and high price (it costs
between $100 and $200)111https://driphydration.com/blog/rapid-pcr-test-cost/
and discomfort of PCR, carrying out this plan can help to quickly diagnose the
disease and prevent the spread of the disease.

Machine learning (ML) is a powerful and innovative tool that can help
healthcare professionals, policy makers, and other stakeholders during
decision-making processes. These computer system models learn and adapt
information using algorithms and statistical networks to analyze and draw
conclusions from patterns in data[de Fátima Cobre et al., 2021].

This paper aims to test the ability to use urine parameters to detect
Covid-19. Then we try to propose and develop a system based on machine
learning that can detect patients with corona disease based on their urine
tests. For this purpose, a new dataset has been collected, and several machine
learning models have been tested after studying the statistical relation of
parameters. Finally, based on the best model, a new ensemble model proposed.

The rest of the paper is organized as follows: Section 2 introduces the
material and methods of this paper. Section 3 presents the experimental
results in this paper. Section 4 concludes this manuscript and indicates
future works.

##  2 Material and methods

###  2.1 Data set

In order to provide data for the training of the models, samples of suspected
corona patients were collected from 4 Covid-19 test centers in the cities of
Yasuj, Shiraz, and Tehran. The method of collection was that a urine sample
and clinical information were received from each of the clients for PCR
testing, and in a short period after receiving the urine sample, an image was
captured from the urine test strip (the Combi 11 made by Machey Nagel). Figure
1 shows a urine test strip. The process of data collection includes two
stages: gathering and analyzing. Information of this dataset shows in Table 1.

Figure 1: A sample of urine test strip and its buffers Statistics | Value
---|---
Number of Samples | 5293
Average of age | 38.41
Gender distribution | Male: 3023 - Female: 2270
Distribution of healthy and sick people | Healthy: 3304 - Sick: 1989
Table 1: Statistics of collected data sets

####  2.1.1 Data Gathering

The data relating to people who were referred to 4 corona test centers in
Yasuj, Shiraz, and Tehran have been collected. The gender information of these
clients by the laboratory is shown in Table 2. Also, Table3 provides clinical
information on patients with respect to three health indicators: diabetes,
blood pressure, and smoking. The table presents the number of positive and
negative cases of each indicator for each of the four data centers (A, B, C,
and D).

For diabetes, data center B has the highest number of positive cases (112),
followed by data center A (39), data center D (26), and data center C (9).
Similarly, for blood pressure, data center B has the highest number of
positive cases (153), followed by data center A (44), data center D (39), and
data center C (15).

Regarding smoking, data center B has the highest number of positive cases
(263), followed by data center D (42), data center A (13), and data center C
(0). It is worth noting that data center C has no positive cases of smoking.

These results provide valuable insights into the prevalence of these health
indicators across different data centers. They can be used to identify areas
where more attention is needed for disease prevention and management. Further
analysis could be conducted to investigate the potential factors that
contribute to the differences observed among the data centers.

Data Center | Male | Female
---|---|---
A | 427 | 327
B | 1694 | 1230
C | 244 | 215
D | 254 | 311
Table 2: Gender of clients by different laboratories Statistics | Diabets | Blood pressure | Smoking
---|---|---|---
Data center | Positive | Negative | Positive | Negative | Positive | Negative
A | 39 | 715 | 44 | 710 | 13 | 741
B | 112 | 2812 | 153 | 2771 | 263 | 2661
C | 9 | 450 | 15 | 444 | 0 | 459
D | 26 | 538 | 39 | 525 | 42 | 522
Table 3: The clinical information of the patients according to different
laboratories

####  2.1.2 Data analysis

For analyzing the images that were captured from the urine test strip, all of
the images of a sample are converted to red, green, and blue. In addition to
this information, Figure2 Shows list of parameters of each sample.

Figure 2: List of parameters of each sample

For more analysis, Figure 3 shows a heat-map of the correlation between RGB
space of urine analysis features and the PCR test. In this image, dark colors
indicate weaker relationships, and lighter colors indicate stronger
relationships. Among the available urine analysis parameters, as a new result
of this paper, ascorbic acid is more correlated with PCR test. It is
considered a key parameter in the diagnosis of the disease. Also, there is a
weak relationship between the presence of protein in urine and the disease of
Covid-19, which was previously shown in the other studies[Bonetti et al.,
2020, Liu et al., 2020].

Figure 3: Correlation of urine analysis parameters

Figure 4 shows the correlation between the clinical and PCR test variables.
The first thing that can be examined in Figure 4 is the weak relationship
between the clinical variables and the PCR test. So results show that clinical
variables are less correlated than urine parameters. This can be seen from the
darkness of the last column of the image. Also, checking the image’s first
column shows a weak correlation between gender and smoking (0.16)[Abuse,
2020]. Examining the second column of this image shows that age has a stronger
correlation with blood pressure and diabetes variables (0.295 and
0.237)[Gurven et al., 2012, Suastika et al., 2012]. The third column shows
that diabetes and blood pressure variables are also correlated (0.272)[Bakris,
2004].

Figure 4: Correlation of clinical and PCR test values

###  2.2 Machine learning models

After going through the data collection process to develop the main part of
the system, algorithms of Sickit-learn has were used. Specifically, the neural
network class of this package called MLPClassifier and three different
classifier namly Random Forest, Logestic Regression, and Gradient Boosting
have been chosen for training. Finally, a greedy search has been engaged for
setting hyperparameters of the models222Hyperparameters space provided in the
supplementary materials.. Also, for stabilization And increasing the validity
of results, the mentioned process has been repeated 30 times, and the results
are reported as the average of all the executions.

The importance of color spaces for image classification has been demonstrated
in literature[Velastegui et al., 2021]. So in order to enrich feature
representation and possibly improve the results, we have used an ensemble
model on different color spaces. First, the RGB representation is converted to
ten spaces including hed, hsv, lab, xyz, ycbcr, ypbpr, yuv, cie, ydbdr, yiq.
Therefore, in addition to RGB, we have different color spaces for training.
Then the mentioned method was used to train the model for every 10 different
color spaces. The ensemble model is a selection model based on maximum votes
from the trained models.

Figure5 illustrates the steps from data collection to model training for
Covid-19 screening. RGB is a color model used in digital imaging and computer
graphics, and it is represented by three channels: red, green, and blue. Each
channel represents the intensity of a color component, and when combined, they
can create a wide range of colors. To extract RGB features, we use image
processing techniques to separate the three channels of an image and represent
them as individual matrices. These matrices can be used as input to a machine
learning algorithm to train the model. One common approach in using RGB
features for machine learning is to flatten the matrices and concatenate them
into a single vector, which can then be used as input to a classifier.

Figure 5: Steps from data collection to model training for Covid-19 screening

###  2.3 Experimental setup

The dataset was split into training and test sets, with 90% of the data used
for training and 10% for testing. To optimize the performance of the models,
we utilized different types of classifiers, including logistic regression,
support vector machines, random forest, and neural networks.

In particular, we used the MLPClassifier algorithm, which is a multi-layer
perceptron algorithm that trains using backpropagation. This algorithm is
available in the Scikit-learn (Sklearn) library, which is a popular machine
learning toolkit in Python. Additionally, we used the GridSearchCV class in
Scikit-learn for hyperparameter tuning. The hyperparameter settings for the
MLP model included two hidden layers with (3,2) as hidden layer sizes. We
repeated the process 30 times and reported the average of the measures to
ensure the stability of the results.

Overall, this process of training and tuning the models plays a critical role
in determining the accuracy of the model in predicting the outcome variable.
Using different classifiers and optimizing the hyperparameters, we can find
the best model that fits training dataset and achieves the highest
performance. The use of Scikit-learn and GridSearchCV makes the process more
efficient and systematic, allowing for a more comprehensive search of the
hyperparameter space.

###  2.4 Evaluation metrics

We use 4 measures for evaluating the models, including accuracy, precision,
sensitivity and recall. Accuracy shows the ratio of correct predictions of
Covid-19 to the total number of samples and is computed as:

$\displaystyle Accuracy=\frac{TP+TN}{TP+TN+FN+FP}$ (1)

Precision shows the percentages of the reported positive Covid-19 that are
correctly detected:

$\displaystyle Precision=\frac{TP}{TP+FP}$ (2)

Sensitivity is the probability of a positive Covid-19 test, conditioned on
truly being positive. In other words, it measures the percentage of positive
cases which the model identified:

$\displaystyle sensitivity=\frac{TP}{TP+FN}$ (3)

Specificity is the probability of a negative Covid-19 test, conditioned on
truly being negative. In other words, it measures the percentage of negative
cases which the model identified or the recall of negative class:

$\displaystyle specificity=\frac{TN}{TN+FP}$ (4)

Where TP represents the number of True Positive results, TN represents the
number of True Negative results, FP represents the number of False Positive
results, and FN represents the number of False Negative results.

##  3 Results

This section compares the average of two healthy and sick communities based on
different parameters. It should be noted that for each parameter of the urine
test strip, three channels of color have been formed, and a total of 33
independent T-tests have been conducted. the results are presented in Table 4.
It can be observed that the different values of blue, green, and red colors
have a significant difference (p-value $<$ 0.05) among most cases in the
studied groups. Also it can be observed that the red channel has the most
difference between groups, Followed by blue, and green.

| Channel/Parameter | Blue | Green | Red
---|---|---|---|---
P-Value | Blood | 0.38807 | 0.0091 | 0.000
Urobilinogen | 0.00001 | 0.000 | 0.000
Bilirubin | 0.00019 | 0.62763 | 0.00498
Protein | 0.000 | 0.000 | 0.06223
Nitrite | 0.00289 | 0.25097 | 0.00146
Ketone | 0.00001 | 0.000 | 0.00116
Ascorbic Acid | 0.000 | 0.000 | 0.00004
Glucose | 0.42476 | 0.45047 | 0.00012
PH | 0.16901 | 0.37185 | 0.4271
Specific Gravity | 0.05597 | 0.60356 | 0.02967
Leukocytes | 0.99463 | 0.04294 | 0.61554
Table 4: P-value of T-test for healthy and sick groups

Table 5 shows the evaluation results of the trained model using RGB color
space features with four different models. Looking at the results in the
table, we can see that the MLP model achieved the highest precision (50.97%)
and recall (61.45%) scores, indicating that it had the highest proportion of
true positives and correctly identified the most actual positive cases out of
all the models. Additionally, the MLP model achieved the highest accuracy
score (63.55%), indicating that it had the highest proportion of correct
predictions overall.

However, it’s worth noting that the Logistic Regression model achieved the
highest specificity score (65.56%), indicating that it had the highest
proportion of true negatives out of all the models. This means that the
Logistic Regression model was the most effective at correctly identifying
negative cases. Overall, the table suggests that the MLP model may be the most
effective at predicting positive cases, while the Logistic Regression model
may be the most effective at predicting negative cases. However, the choice of
model ultimately depends on the specific goals and requirements of the
analysis.

Model | Precision | Recall | Specificity | Accuracy
---|---|---|---|---
Random Forest | 49.45 | 59.07 | 63.99 | 62.11
Logistic Regression | 50.91 | 60.09 | 65.56 | 63.51
Gradient Boosting | 48.06 | 60.72 | 60.83 | 60.79
MLP | 50.97 | 61.45 | 64.81 | 63.55
Table 5: Evaluation of models

Table 6 shows the evaluation of the ensemble model. As can be observed, the
overall performance of the ensemble model is better than the model based on
the RGB color space, and based on all criteria, except recall, it has improved
the results.

Based on the results of the evaluation of the ensemble model, it seems that
the ensemble model is performing better than the individual models that were
evaluated in the table 5.

Specifically, the ensemble model has an accuracy score of 64.69%, which is
higher than the accuracy scores of all four individual models evaluated in the
previous table. Additionally, the ensemble model has a precision score of
52.63%, which is higher than the recall scores of all the individual models.
The ensemble model also has a specificity score of 66.83%, which is higher
than the specificity scores of other models. However, the recall score of the
ensemble model is 61.13%, which is lower than the precision scores of two of
the MLP model.

Overall, the results of the evaluation suggest that the ensemble model is a
promising approach for this particular task, as it achieves higher accuracy,
precision, and specificity scores than the individual models. However, the
relatively low recall score is a potential concern, and further analysis would
be needed to determine whether the ensemble model is the best approach for
this particular problem.

Model | Precision | Recall | Specificity | Accuracy
---|---|---|---|---
Random Forest | 49.45 | 59.07 | 63.99 | 62.11
Logistic Regression | 50.91 | 60.09 | 65.56 | 63.51
Gradient Boosting | 48.06 | 60.72 | 60.83 | 60.79
MLP | 50.97 | 61.45 | 64.81 | 63.55
Ensemble model | 52.63 | 61.13 | 66.83 | 64.69
Table 6: Evaluation of ensemble model trained on ten color spaces

For determining an appropriate cut-off point of sensitivity and specificity of
the test, the ROC curve can be used. The selection of the test threshold
depends on the test’s purpose and does not necessarily provide the same
sensitivity and specificity to achieve higher accuracy [Hoo et al., 2017]. A
ROC curve following the diagonal line y=x produces false positives at the same
rate as true positives. Therefore, we expect the average accuracy ROC curve
from the diagnostic test to lie above the y=x line (the reference line) in the
upper left triangle, as shown in figure6.

Figure 6: ROC curve of the model

Due to the lack of sufficient strength of the results, a part of the space in
which the model has an uncertain performance has been specified as the
uncertain region, and at the time of service, if a sample from that part of
the space with uncertainty is observed, it will be answered with the result of
”lack of sufficient information”. In order to identify uncertain parts, the
number of classifiers with positive votes between all eleven color spaces can
be used. Based on this, Table 7 shows the results of removing the uncertain
regions of space with different values. It can be observed that reducing the
response space increases specificity and decreases sensitivity. Also, the
growth rate of specificity is more than the sensitivity decrease. The purpose
of using a screening test is to extract healthy people and send the rest to
tests with higher sensitivity; for example, in table 7 when we have ten
positive classifiers, the model has a relatively good performance Because it
has a good screening and extracts about 80% of the healthy population and only
20% need to do other tests like the PCR.

| The minimum number
---
of classifiers
with positive vote
| The ratio of
---
response space
to total space
Accuracy | Recall | Specificity
6 | 1 | 0.646918 | 0.611371 | 0.668347
7 | 0.956226 | 0.659432 | 0.593706 | 0.699072
8 | 0.863899 | 0.671666 | 0.595458 | 0.718033
9 | 0.768302 | 0.687787 | 0.595977 | 0.742738
10 | 0.652327 | 0.708446 | 0.574064 | 0.782531
11 | 0.375094 | 0.780684 | 0 | 1
Table 7: The results of removing uncertain parts of the ensemble model

##  4 Conclusion

One of the essential stages of machine learning projects is the feasibility
stage, in which efforts are made to solve exploratory problems where the
success or failure of the model in performing properly in that field cannot be
determined with certainty. As the results of this paper and other
studies[Morell-Garcia et al., 2021, Zolotov et al., 2022], due to the absence
of serious distinguishing complications and the common complications caused by
Covid-19 with other diseases, it is not possible to detect Covid-19 using
urinary parameters with high accuracy. Nevertheless, we indicate the ability
of urinary parameters like ascorbic acid to screen Covid-19 patients. Based on
this, we propose an ensemble model with eleven classifiers, the model has a
relatively good performance Because it has a good screening and extracts about
80% of the healthy population and only 20% of response space need to do other
tests like the PCR. For future work, efforts will be made to conduct tests
after one-handed sampling (currently, samples are not taken on an empty
stomach at specific times and under similar conditions), and label the samples
based on the final condition of the patient to early detect whether COVID-19
will lead to severe conditions in a particular individual or not. Due to the
development of the necessary interfaces for the extraction and analysis of
urinary parameters, it is possible in future to find the items developed in
this project to evaluate covid-19 with higher accuracy and use these
parameters for screening other diseases.

##  5 Declarations

###  5.1 Ethics approval and consent to participate

The research study conducted adhered to the ethical guidelines and regulations
set forth by the Shiraz University of Medical Sciences. All methods employed
in this study were performed in strict accordance with these guidelines to
ensure the ethical treatment of participants and compliance with applicable
regulations. Prior to the commencement of the study, appropriate approvals and
permissions were obtained from the relevant institutional review board or
ethics committee. Informed consent was obtained from all participants, and
their confidentiality and privacy rights were strictly upheld throughout the
study. The data collection procedures, data analysis, and reporting followed
established protocols and guidelines, ensuring the integrity and validity of
the findings. The researchers involved in this study were experienced
professionals who were well-versed in the ethical considerations and
guidelines specific to the field of study. The study findings are thus
reliable and can be interpreted within the context of the ethical and
regulatory framework governing research in this domain. All participants were
informed about the purpose and procedures of the study and signed written
informed consent forms before data collection. The participants were assured
that their participation was voluntary and that they could withdraw from the
study at any time without any consequences. The participants were also assured
that their personal information and responses would be kept confidential and
used only for research purposes.

###  5.2 Consent for publication

Consent for publication was obtained from all participants involved in this
study. Participants were provided with a clear explanation of the purpose and
nature of the publication, as well as any potential risks or benefits
associated with it. They were assured that their personal identities and
sensitive information would remain confidential and protected. Participants
were given the opportunity to ask questions and were informed that their
participation was voluntary, and they could withdraw their consent at any time
without any negative consequences. Consent forms were signed and collected
from each participant prior to their inclusion in this publication. To ensure
anonymity, pseudonyms or identifiers have been used throughout the manuscript
to protect the privacy and confidentiality of the participants. This study
adheres to the principles outlined in the Declaration of the Shiraz University
of Medical Sciences regarding the ethical conduct of research involving human
subjects. The publication of this research has been carried out with the full
understanding and consent of all participants involved.

###  5.3 Availability of data and materials

The current study generated and/or analysed some data that are not publicly
available due to data protection and privacy reasons. These data contain
personal or sensitive information that require the consent of the individuals
or groups involved. The corresponding author can provide these data upon
reasonable request.

###  5.4 Funding and/or Conflicts of Interests/Competing Interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

###  5.5 Authors’ contributions

This paper was written by Mohammad Hadi Goldani as the corresponding author.
Behzad Moayedi, Abdalsamad Keramatfar, Mohammad Javad Fallahi, Alborz
Jahangirisisakht, Mohammad Saboori, and Leyla Badiei contributed to the
research, data analysis, and manuscript preparation. Specifically, Behzad
Moayedi, Abdalsamad Keramatfar and Mohammad Hadi Goldani were responsible for
the implementation of the machine learning algorithms, preprocessing and
evaluation of the model performance. Mohammad Javad Fallahi, Mohammad Saboori
and Leyla Badiei contributed to the data collection, analysis of the
experimental results and Alborz Jahangirisisakht contributed to data
collection. All authors have read and approved the final manuscript and agree
to be accountable for all aspects of the research, including the accuracy and
integrity of the data and its analysis. Signed by, Mohammad Hadi Goldani
(corresponding author) Behzad Moayedi Abdalsamad Keramatfar Mohammad Javad
Fallahi Alborz Jahangirisisakht Mohammad Saboori Leyla Badiei

###  5.6 Acknowledgment

The authors would like to acknowledge the department of health affairs, Shiraz
University of Medical Sciences (SUMS), Shiraz, Iran for their close
collaboration and strong support in facilitating a reliable data gathering
phase in this project. Also, The authors express their gratitude to the RCDAT
(Research Center for Development of Advanced Technologies) for research
funding.

## References

  * Abuse [2020] Abuse, N. (2020).  Are there gender differences in tobacco smoking.
  * Announcement [2022] Announcement, W. (2022).  Who announcement.  https://www.who.int/publications/m/item/weekly-epidemiological-update-on-covid-19—6-july-2022.
  * Asiaee et al. [2020] Asiaee, M., Vahedian-Azimi, A., Atashi, S. S., Keramatfar, A., & Nourbakhsh, M. (2020).  Voice quality evaluation in patients with covid-19: An acoustic analysis.  Journal of Voice, .
  * Bakris [2004] Bakris, G. L. (2004).  The importance of blood pressure control in the patient with diabetes.  The American journal of medicine, 116, 30–38.
  * Banerjee et al. [2020] Banerjee, A., Ray, S., Vorselaars, B., Kitson, J., Mamalakis, M., Weeks, S., Baker, M., & Mackenzie, L. S. (2020).  Use of machine learning and artificial intelligence to predict sars-cov-2 infection from full blood counts in a population.  International immunopharmacology, 86, 106705.
  * Bonetti et al. [2020] Bonetti, G., Manelli, F., Bettinardi, A., Borrelli, G., Fiordalisi, G., Marino, A., Menolfi, A., Saggini, S., Volpi, R., Adamini, R. et al. (2020).  Urinalysis parameters for predicting severity in coronavirus disease 2019 (covid-19).  Clinical Chemistry and Laboratory Medicine (CCLM), 58, e163–e165.
  * Chen et al. [2020] Chen, J., Wu, L., Zhang, J., Zhang, L., Gong, D., Zhao, Y., Chen, Q., Huang, S., Yang, M., Yang, X. et al. (2020).  Deep learning-based model for detecting 2019 novel coronavirus pneumonia on high-resolution computed tomography.  Scientific reports, 10, 1–11.
  * Despotovic et al. [2021] Despotovic, V., Ismael, M., Cornil, M., Mc Call, R., & Fagherazzi, G. (2021).  Detection of covid-19 from voice, cough and breathing patterns: Dataset and preliminary results.  Computers in Biology and Medicine, 138, 104944.
  * Erdogan et al. [2021] Erdogan, O., Ok, F., Carkci, S., & Durmus, E. (2021).  Is there an association between urine biochemical parameters on admission and the severity of covid-19?  International Journal of Clinical Practice, 75, e14809.
  * de Fátima Cobre et al. [2021] de Fátima Cobre, A., Stremel, D. P., Noleto, G. R., Fachi, M. M., Surek, M., Wiens, A., Tonin, F. S., & Pontarolo, R. (2021).  Diagnosis and prediction of covid-19 severity: can biochemical tests and machine learning be used as prognostic indicators?  Computers in biology and medicine, 134, 104531.
  * Gurven et al. [2012] Gurven, M., Blackwell, A. D., Rodríguez, D. E., Stieglitz, J., & Kaplan, H. (2012).  Does blood pressure inevitably rise with age? longitudinal evidence among forager-horticulturalists.  Hypertension, 60, 25–33.
  * Hoo et al. [2017] Hoo, Z. H., Candlish, J., & Teare, D. (2017).  What is an roc curve?
  * Kovács et al. [2021] Kovács, A., Palásti, P., Veréb, D., Bozsik, B., Palkó, A., & Kincses, Z. T. (2021).  The sensitivity and specificity of chest ct in the diagnosis of covid-19.  European radiology, 31, 2819–2824.
  * Liu et al. [2020] Liu, R., Ma, Q., Han, H., Su, H., Liu, F., Wu, K., Wang, W., & Zhu, C. (2020).  The value of urine biochemical parameters in the prediction of the severity of coronavirus disease 2019.  Clinical Chemistry and Laboratory Medicine (CCLM), 58, 1121–1124.
  * Morell-Garcia et al. [2021] Morell-Garcia, D., Ramos-Chavarino, D., Bauça, J. M., Argente del Castillo, P., Ballesteros-Vizoso, M. A., Garcia de Guadiana-Romualdo, L., Gómez-Cobo, C., Pou, J. A., Amezaga-Menéndez, R., Alonso-Fernández, A. et al. (2021).  Urine biomarkers for the prediction of mortality in covid-19 hospitalized patients.  Scientific Reports, 11, 1–13.
  * Roberts et al. [2021] Roberts, M., Driggs, D., Thorpe, M., Gilbey, J., Yeung, M., Ursprung, S., Aviles-Rivero, A. I., Etmann, C., McCague, C., Beer, L. et al. (2021).  Common pitfalls and recommendations for using machine learning to detect and prognosticate for covid-19 using chest radiographs and ct scans.  Nature Machine Intelligence, 3, 199–217.
  * Strasinger & Di Lorenzo [2010] Strasinger, S. K., & Di Lorenzo, M. S. (2010).  Anlisis de orina y de los lquidos corporales/Urinalysis and body fluids.  Ed. Médica Panamericana.
  * Suastika et al. [2012] Suastika, K., Dwipayana, P., Semadi, M. S., & Kuswardhani, R. T. (2012).  Age is an important risk factor for type 2 diabetes mellitus and cardiovascular diseases.  Glucose Tolerance, (pp. 67–80).
  * Velastegui et al. [2021] Velastegui, R., Yang, L., & Han, D. (2021).  The importance of color spaces for image classification using artificial neural networks: a review.  In Computational Science and Its Applications–ICCSA 2021: 21st International Conference, Cagliari, Italy, September 13–16, 2021, Proceedings, Part II 21 (pp. 70–83).  Springer.
  * Wu et al. [2020] Wu, J., Zhang, P., Zhang, L., Meng, W., Li, J., Tong, C., Li, Y., Cai, J., Yang, Z., Zhu, J. et al. (2020).  Rapid and accurate identification of covid-19 infection through machine learning based on clinical available blood test results.  MedRxiv, .
  * Yetisen et al. [2013] Yetisen, A. K., Akram, M. S., & Lowe, C. R. (2013).  based microfluidic point-of-care diagnostic devices.  Lab on a Chip, 13, 2210–2251.
  * Zolotov et al. [2022] Zolotov, E., Sigal, A., Havrda, M., Jeřábková, K., Krátká, K., Uzlová, N., & Rychlík, I. (2022).  Can renal parameters predict the mortality of hospitalized covid-19 patients?  Kidney and Blood Pressure Research, 47, 309–319.


#  Development of the fully Geant4 compatible package for the simulation of
Dark Matter in fixed target experiments

B. Banto Oberhauser  0009-0006-4795-1008 ETH Zürich, Institute for Particle
Physics and Astrophysics, CH-8093 Zürich, Switzerland  P. Bisio
0009-0006-8677-7495 INFN, Sezione di Genova, 16147 Genova, Italia  Universitá
degli Studi di Genova, 16126 Genova, Italia  A. Celentano  0000-0002-7104-2983
INFN, Sezione di Genova, 16147 Genova, Italia  E. Depero  0000-0003-2239-1746
ETH Zürich, Institute for Particle Physics and Astrophysics, CH-8093 Zürich,
Switzerland  R. R. Dusaev  0000-0002-6147-8038 Tomsk Polytechnic University,
634050 Tomsk, Russia  D. V. Kirpichnikov  0000-0002-7177-077X Institute for
Nuclear Research, 117312 Moscow, Russia  M. M. Kirsanov  0000-0002-8879-6538
<EMAIL_ADDRESS>Institute for Nuclear Research, 117312 Moscow, Russia
N. V. Krasnikov  0000-0002-8717-6492 Institute for Nuclear Research, 117312
Moscow, Russia  Joint Institute for Nuclear Research, 141980 Dubna, Russia  A.
Marini  0000-0002-6778-2161 INFN, Sezione di Genova, 16147 Genova, Italia  L.
Marsicano  0000-0002-8931-7498 INFN, Sezione di Genova, 16147 Genova, Italia
L. Molina-Bueno  0000-0001-9720-9764 Instituto de Física Corpuscular,
Universidad de Valencia and CSIC, Carrer del Catedrátic José Beltrán Martinez,
2, 46980 Paterna, Valencia, Spain  M. Mongillo  0009-0000-7331-4076 ETH
Zürich, Institute for Particle Physics and Astrophysics, CH-8093 Zürich,
Switzerland  D. Shchukin  0009-0007-5508-3615 P.N. Lebedev Physical Institute,
Moscow, Russia, 119991 Moscow, Russia  H. Sieber  0000-0003-1476-4258 ETH
Zürich, Institute for Particle Physics and Astrophysics, CH-8093 Zürich,
Switzerland  I. V. Voronchikhin  0000-0003-3037-636X Tomsk Polytechnic
University, 634050 Tomsk, Russia

###### Abstract

The search for new comparably light (well below the electroweak scale) feebly
interacting particles is an exciting possibility to explain some mysterious
phenomena in physics, among them the origin of Dark Matter. The sensitivity
study through detailed simulation of projected experiments is a key point in
estimating their potential for discovery.

Several years ago we created the DMG4 package for the simulation of DM (Dark
Matter) particles in fixed target experiments. The natural approach is to
integrate this simulation into the same program that performs the full
simulation of particles in the experiment setup. The Geant4 toolkit framework
was chosen as the most popular and versatile solution nowadays.

The simulation of DM particles production by this package accommodates several
possible scenarios, employing electron, muon or photon beams and involving
various mediators, such as vector, axial vector, scalar, pseudoscalar, or spin
2 particles. The bremsstrahlung, annihilation or Primakoff processes can be
simulated.

The package DMG4 contains a subpackage DarkMatter with cross section methods
weakly connected to Geant4. It can be used in different frameworks.

In this paper, we present the latest developments of the package, such as
extending the list of possible mediator particle types, refining formulas for
the simulation and extending the mediator mass range. The user interface is
also made more flexible and convenient.

In this work, we also demonstrate the usage of the package, the improvements
in the simulation accuracy and some cross check validations.

## New version program summary

_Program title:_ DMG4
_CPC Library link to program files:_
_Code Ocean capsule:_
_Licensing provisions:_ GNU General Public License 3 (GPL)
_Programming language:_ c++
_Nature of problem:_ For the simulation of Dark Matter production processes in
fixed target experiments a code that can be easily integrated in programs for
the full simulation of experimental setup is needed.
_Solution method:_ A fully Geant4 compatible DM simulation package DMG4 was
presented in 2020. We present numerous further developments of this package.

##  I Introduction

The _DMG4_ package is designed for simulating the production of feebly
interacting particles beyond the Standard Model (SM) present in many models
trying to explain Dark Matter and various anomalies observed in particle
physics. It deals with particles of masses well below the electroweak scale,
extending down to $\sim 1$ keV. The experimental technique involving a fixed
thick target is very popular in the searches for these particles [1, 2, 3, 4].
It is convenient to simulate their production in the same program that is used
for the simulation of the whole setup. For this reason, the DMG4 package by
construction is made fully compatible with the Geant4 toolkit [5], which is a
foundation for the majority of such programs. The initial description of the
DMG4 package is presented in ref. [6]. In this paper we describe its further
developments.

While the most popular DM model introduces a vector particle, “dark photon”
[7] as a mediator between SM and DM particles, different quantum numbers and
couplings of such mediators are possible [8, 9]. We tried to provide the
possibility of simulation for as wide as possible class of models that can be
relevant for the searches for BSM (beyond SM) particles at accelerators with
moderate beam energies $\sim\mathcal{O}(10-100)$ GeV. The recent extensions we
have implemented are presented in this article.

Several relevant production processes involve $2\rightarrow 3$ interactions on
a heavy nucleus, for which is difficult to derive exactly all formulas needed
for efficient simulation. Various approximations are often used, introducing
deviations from results obtained at exact tree level (ETL). The latter results
serve as reference, but cannot be used in the codes that perform event-by-
event sampling. We performed a number of studies and developments in order to
improve the accuracy of simulation, which we present in this article.

The DM mediators can be produced through annihilation processes, either of
electromagnetic shower positrons or beam positrons with target electrons. The
corresponding exact formulas were derived and used in the package. However,
for such processes, there were other implementation problems that we have
recently addressed. We present the corresponding solutions here.

##  II DMG4 package structure

Structurally, the DMG4 package consists of two main components: a set of
lower-level routines responsible for DM-related calculus, and a layer of
Geant4 interfaces realization. This way various numerical tests and model
benchmarking can be implemented using only the numerical routines part, before
it gets embedded into Geant4 pipeline.

As for numerical routines, every process subclasses common parent root
DarkMatter (being an _abstract product_ in terms of _factory pattern_)
facilitating common requirements for dark matter models in terms of mean free
path, spatial dependencies on cross sections, decay width, etc. Subclasses
implementing a particular model (_concrete products_) implement its interface
(_DM production process_) by overriding certain virtual methods of DarkMatter
base class.

Figure 1: Classes of the main DMG4 interfaces shown as UML diagram [10]. For
brevity, in bremsstrahlung part only dark photons and dark axials classes are
shown.

The base DarkMatter class is extended by the _DM annihilation process_
interface within DarkMatterAnnihilation. The _DM annihilation process_ appends
contract with routines yielding certain resonant-specific quantities
(additional yield factor and angular distribution PDF). The diagram 1
illustrates relationships between interfaces and corresponding implementation
classes with dark axial and dark photon modules as examples of concrete
products. A particular sampler instance can be associated with physics model
at the most abstract level.

Component diagram 2 provides a wider view on how particular numerical routines
are exposed in terms of Geant4 API (namely, G4ParticleDefinition and
G4VDiscreteProcess). Corresponding interfaces are implemented within the
interim layer composed of a set of processes and related particles while the
cross-section calculus is expressed by means of internal classes.

Figure 2: UML Components diagram [10] of DMG4 package exposing
Geant4-compatible facade.

The DM sector particles introduced in the package are listed in Table 1. The
PDG codes are ascribed according to the slightly extended rules in [11]. As
compared to early versions, we reduced the number of mediator particle
classes, which was possible due to usage of a new class
_DarkMatterParametersFactory_. Using this class it is more convenient to
specify the decay modes and some other parameters of particles. On the other
hand, we added several DM particles of models with semivisible decay modes.
Note that in many configurations we assume the mediators to be stable,
although in full models they decay into stable DM particles. But for the
missing energy signature simulation it is unimportant as the energy carried
away is the same.

Table 1: DM particles defined in the package DMG4 Name | PDG ID | emitted by | spin | parity | Model | stable? | decay
---|---|---|---|---|---|---|---
DMParticleAPrime | 5500022 | $e^{+},e^{-}$ | 1 | 1 |  | true | -
| 5500122 |  |  |  | Visible X | false | $e^{+}e^{-}$
| 5500222 |  |  |  | B - L | false |  $e^{+}e^{-}$, $\mu^{+}\mu^{-}$ etc.
| 5500322 |  |  |  | Inelastic DM | false | $\chi_{1},\chi_{2}$
DMParticleScalar | 5400022 | $e^{+},e^{-}$ | 0 | 1 |  | true | -
| 5400022 |  |  |  | Visible X | false | $e^{+}e^{-}$
DMParticlePseudoScalar | 5410022 | $e^{+},e^{-}$ | 0 | -1 |  | true | -
| 5410122 |  |  |  | Visible X | false | $e^{+}e^{-}$
DMParticleAxial | 5510022 | $e^{+},e^{-}$ | 1 | -1 |  | true | -
| 5510122 |  |  |  |  |  | $e^{+}e^{-}$
DMParticleZPrime | 5500023 | $\mu$ | 1 | 1 | $L_{\mu}-L_{\tau}$ | true | -
| 5500023 | $\mu$ |  |  | $L_{\mu}-L_{\tau}$ | false |  $\nu\nu$, $\mu^{+}\mu^{-}$
DMParticleALP | 5300122 | $\gamma$ | 0 | -1 |  | false | $\gamma\gamma$
DMParticleChi1 | 5200014 | mediator | 1/2 | 1 | Inelastic DM | true | -
DMParticleChi2 | 5200013 | mediator | 1/2 | 1 | Inelastic DM | false | $\chi_{1}e^{+}e^{-}$
DMParticleChi | 5200012 | mediator | 1/2 | 1 |  | true | -

The dark sector particles are assumed to be stable as they are anyway
undetectable. This does not apply to the $\chi_{2}$ appearing in models of
inelastic Dark Matter[12]. For this particle, the 3-body decay is simulated
using the matrix element sampling implemented in the new
G4iDM3bodyDecayChannel class, as explained in Sec.IV.

The current version of DMG4 package contains the following processes of DM
production:

  * •

Bremsstrahlung-like process of the type $bN\to bNX$, where $b$ is a projectile
(can be $e^{-},e^{+},\mu^{-},\mu^{+}$), and $X$ is a DM mediator. Apart from
the mediators listed in the table, we added a process with a spin 2 particle
without adding a corresponding particle class. It means it is possible to
simulate the invisibly decaying spin 2 mediator.

  * •

Primakoff process of photon conversion $\gamma N\to aN$, where $a$ is an
axion-like particle (ALP) [13]

  * •

Resonant in-flight positron annihilation on atomic electrons
$e^{+}e^{-}\rightarrow X\rightarrow\chi\chi$, where $\chi$ is a dark matter
mediator decay product [14].

The physics for a simulation run is configured in the function
DarkMatterPhysicsConfigure called from the constructor of the steering class
DarkMatterPhysics. The configuration is done using the singleton class
_DarkMatterParametersFactory_. A user specifies a number of parameters, then
the steering class performs the needed actions:

  * •

Creates an instance of one of the concrete classes corresponding to the needed
process and derived from the base class DarkMatter, for example DarkPhotons

  * •

Instantiates and registers the needed particles and processes provided by the
DMG4 package in terms of the native Geant4 API.

The following main parameters are to be specified: 1. The process code; 2. The
mixing (or coupling) parameter $\epsilon$; 3. The mediator mass; 4\. The decay
mode; 5. The cut-off minimal energy of particles that can initiate the
processes of DM production. The latter is needed to avoid the simulation of
very soft DM particles that are anyway undetectable. In the new versions, this
parameter also specifies, for some processes or for a part of the mediator
mass range, a minimal mediator energy that is allowed to be generated.

##  III Package DarkMatter

This class, which we try to keep as weakly connected to Geant4 as possible,
contains the methods that return total and differential cross sections, first
of all for the bremsstrahlung processes. For these processes the cross
sections can be calculated numerically at Exact Tree Level (ETL) or using the
formulas of Weizsaker-Williams (WW) and improved Weizsaker-Williams (IWW)
approximations [1].

As explained in [6], for the electron beam it contains the tabulated K-factors
that correct the total cross sections obtained in IWW approximations for the
electron beam to the values calculated at ETL. As compared to the previous
version, we extended the corresponding tables down to the mediator masses of 1
keV and, for some classes, up to 3 GeV.

In addition to this, we also tabulated the differential cross sections
$d\sigma/dx$ for the electron beam and mediator masses below 1 MeV, where it
can have non-trivial shape, significantly different from the sharply peaked
one at around $x=1$ for heavier mediators.

The Lagrangians and initial cross section formulas used in the package can be
found in the previous paper [6]. The Lagrangian for the new model with spin 2
mediator looks like this [15, 16, 17]:

$\displaystyle\mathcal{L}^{\rm G}_{\rm eff}\supset-\frac{ic^{\rm
G}_{ee}}{2\Lambda}G^{\mu\nu}\left(\overline{e}\gamma_{\mu}\overleftrightarrow{D}_{\nu}e-\eta_{\mu\nu}\overline{e}\gamma_{\rho}\overleftrightarrow{D}^{\rho}e\right)$
(1) $\displaystyle+\frac{c^{\rm
G}_{\gamma\gamma}}{\Lambda}G^{\mu\nu}\left(\frac{1}{4}\eta_{\mu\nu}F_{\lambda\rho}F^{\lambda\rho}+F_{\mu\lambda}F^{\lambda}_{\nu}\right)+\frac{c^{\rm
G}_{\rm DM}}{\Lambda}G^{\mu\nu}T^{\rm DM}_{\mu\nu},$

where $e$ is the label of the SM electron,
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is a stress tensor of
the SM photon field $A_{\mu}$, $D_{\mu}=\partial_{\mu}-ieA_{\mu}$ is a
covariant derivative of the $U(1)$ gauge field, and $\Lambda$ is the
dimensional parameter for spin-2 interactions, that is associated with the
scale of new physics; $c_{ee}^{\rm G}$ and $c_{\gamma\gamma}^{\rm G}$ are
dimensionless couplings for the electron and photon respectively. We choose
the universal coupling $c^{\rm G}_{ee}=c^{\rm G}_{\gamma\gamma}$ throughout
the paper.

Finally, since the early versions of the package, significant progress has
been made in the cross sections for the case of a muon beam and in the
annihilation processes. These, alongside other developments, are described in
the following sections.

##  IV Inclusion of inelastic Dark Matter models

In multi-generational dark sector models, such as scenarios involving
inelastic Dark Matter, the dark photon exhibits off-diagonal couplings to a
pair of particles, denoted as $\chi_{1}$ and $\chi_{2}$, with different
masses. The heavier particle, $\chi_{2}$, is unstable and decays via the
channel $\chi_{2}\rightarrow\chi_{1}e^{-}e^{+}$. In such cases, the dark
photon decay cascade leads to both visible and invisible final states and is
consequently referred to as semivisible $A^{\prime}$.

To enable the simulation of these next-to-minimal models, two new particles,
DMParticleChi1 and DMParticleChi2, have been added to the DMG4 package.
Notably, two distinct semivisible models are now included: inelastic Dark
Matter (iDM) [12] and inelastic Dirac Dark Matter (i2DM) [18]. The decay
widths for these models are taken from Refs. [19, 20] and Ref.[18]
respectively.

To account for the specific 3-body decay process
$\chi_{2}\rightarrow\chi_{1}e^{-}e^{+}$, the DMG4 class G4iDM3bodyDecayChannel
has been implemented. This class, derived from G4VDecayChannel and based on
G4DalitzDecayChannel, performs the sampling of the spin-averaged matrix
element, as illustrated in Fig. 3. The kinematics of the final state particles
is subsequently calculated based on the obtained squared invariant masses of
the $e^{+}e^{-}$ and $\chi_{1}e^{+}$ systems.

Figure 3: Dalitz plot for $\chi_{2}\rightarrow\chi_{1}e^{+}e^{-}$ obtained
through the matrix element sampling from the method G4iDM3bodyDecayChannel for
$m_{A^{\prime}}=0.25$ GeV, $m_{A^{\prime}}=3m_{\chi_{1}}$ and $f=0.4$.

##  V Angular sampling of ALP production via Primakoff effect

In order to simulate the final state kinematics of the ALP produced in the
Primakoff process $\gamma N\rightarrow aN$, the sampling of the emission angle
was implemented in DMG4. The differential cross section with respect to the
angle of ALP emission is obtained in the limit when $m_{a}\ll E_{a}$ and
$\theta_{a}\ll 1$ and is used in the sampling scheme, as presented in our
previous work [13]. In this approximation the energy of the emitted ALP
$E_{a}$ can be calculated based on the angle $\theta_{a}$. The results of the
angular sampling are displayed in Fig. 4 and compared to the reference
distribution function, which is a single differential cross section.

Figure 4: Comparison between results obtained using the sampling from the
method SimulateEmissionWithAngle3 and the reference function for the
differential cross-section $d\sigma/d\theta_{a}$ of the Primakoff process
$\gamma N\rightarrow aN$ [13]. The cross-section was calculated for an
incoming photon energy $E_{\gamma}=50$ GeV.

##  VI Developments in muon-philic processes and cross sections
implementations

In this section, we describe the recent developments in the simulation of
light dark matter production in bremsstrahlung-like reactions of muons on
nuclei, $\mu N\rightarrow\mu NX$. In particular, we describe the extension of
muon-philic models to different mediator types and the improvements in the
underlying production cross sections implementation.

###  VI.1 Broadening the class of mediators

In our previous work [6], the emission of a light mediator from muons was
associated with the $L_{\mu}-L_{\tau}$ gauge vector boson, $X=Z^{\prime}$,
[21, 22, 23, 24, 25, 26]. Within this model, the illustrative Lagrangian of
the theory is given by

$\mathcal{L}\supset\mathcal{L}_{SM}-\frac{1}{4}F_{\mu\nu}^{\prime}F^{\mu\nu\prime}+\frac{m_{Z^{\prime}}^{2}}{2}Z_{\mu}^{\prime}Z^{\mu\prime}-g_{Z^{\prime}}Z_{\alpha}^{\prime}\big{(}\bar{\mu}\gamma^{\alpha}\mu-\bar{\tau}\gamma^{\alpha}\tau+\bar{\nu}_{\mu}\gamma^{\alpha}P_{L}\nu_{\mu}-\bar{\nu}_{\tau}\gamma^{\alpha}P_{L}\nu_{\tau}\big{)},$
(2)

where $F_{\mu\nu}^{\prime}$ is the field strength tensor associated with the
$Z_{\mu}^{\prime}$ field, $m_{Z^{\prime}}$ the mass of the gauge boson,
$g_{Z^{\prime}}=\epsilon_{Z^{\prime}}e$ the coupling to SM leptons and $P_{L}$
the left-handed chiral projection operator. The decay width associated with
the _vanilla_ model is given by the purely invisible channel to SM neutrinos,
$Z^{\prime}\rightarrow\bar{\nu}\nu$, such that

$\Gamma_{Z^{\prime}\rightarrow\bar{\nu}\nu}=\frac{\alpha_{Z^{\prime}}m_{Z^{\prime}}}{3},$
(3)

where $\alpha_{Z^{\prime}}=g_{Z^{\prime}}^{2}/(4\pi)$.
The addition to the aforementioned model, the formulas for other types of
light muon-philic mediators is provided, namely for scalar and pseudoscalar
particles, respectively $X=S$ and $X=P$, with corresponding simplified
Lagrangians

$\mathcal{L}\supset\mathcal{L}_{SM}+\frac{1}{2}(\partial_{\mu}S)^{2}-\frac{1}{2}m_{S}^{2}S^{2}+g_{S}S\overline{\mu}\mu,$
(4)
$\mathcal{L}\supset\mathcal{L}_{SM}+\frac{1}{2}(\partial_{\mu}P)^{2}-\frac{1}{2}m_{P}^{2}P^{2}+ig_{P}P\overline{\mu}\gamma_{5}\mu,$
(5)

where it is assumed that the mediators do not need to couple to neutrinos,
have masses $m_{S}$ and $m_{P}$ and muon-specific couplings $g_{S}$ and
$g_{P}$.

###  VI.2 Calculations of the production cross sections

As introduced in Sec. III, the differential and total production cross
sections for the reaction $\mu N\rightarrow\mu NX$ can be computed within the
WW approximation. Within this phase-space approximation, the $2\rightarrow 3$
process depicted above is factorized through the use of the equivalent photon
flux approximation [27, 28] into a $2\rightarrow 2$ process, provided the
virtual photon flux

$\chi^{\text{WW}}=\int_{t_{\text{min}}}^{t_{\text{max}}}dt\
\frac{t-t_{\text{min}}}{t^{2}}F^{2}(t),$ (6)

with $t_{\text{min}}$ and $t_{\text{max}}$ being the minimum and maximum
momentum transfer to the nucleus, and $F^{2}(t)$ the elastic form factor
associated with the nucleus (see e.g. [27]). As such, the double-differential
cross section reads [29]

$\frac{d^{2}\sigma^{X}_{2\to
3}}{dxd\cos\theta_{X}}\Bigr{|}_{\text{WW}}\simeq\frac{\alpha\chi^{\text{WW}}}{\pi(1-x)}E_{0}^{2}x\beta_{X}\frac{d\sigma^{X}_{2\to
2}}{d(pk)}\Bigr{|}_{t=t_{\text{min}}},$ (7)

with $X=Z^{\prime},S,P$, $x$ and $\theta_{X}$ being respectively the
fractional energy of the $X$ boson and its emission angle, $\alpha$ the fine
structure constant, $\beta_{X}$ the corresponding Lorentz $\beta-$factor, $p$
and $k$ the four-momenta associated with the initial-state muon and final-
state $X$ boson, and $E_{0}$ the muon energy before the interaction. Together
with the double differential cross section of Eq. (7), we implemented within
the DarkMatter package the single-differential cross section for $x$, obtained
by performing a complicated analytical integration of Eq. (7). This allowed us
to obtain much better run-time performance of the event sampling code. The
total cross section then reads [30]

$\sigma_{2\rightarrow
3}^{X}\Bigr{|}_{\text{WW}}=\int_{x_{\text{min}}}^{x_{\text{max}}}dx\
\bigg{(}\epsilon_{X}^{2}\alpha^{3}\sqrt{1-\frac{m_{X}^{2}}{E_{0}^{2}}}\frac{1-x}{x}\sum_{i=1}^{6}I_{i}^{X}(x,\tilde{u})\bigg{\rvert}_{\tilde{u}=\tilde{u}_{\text{min}}}^{\tilde{u}=\tilde{u}_{\text{max}}}\bigg{)},$
(8)

where the six special functions $I_{i}$, $i=1,2,...6$ and the definition of
the variable $\tilde{u}_{\text{min}}$ and $\tilde{u}_{\text{max}}$ can be
found in [30]. This provides one with a more precise expression than within
the IWW approximation, where the dependence on $(x,\ \theta_{X})$ is neglected
in Eq. (6). As a result, we don’t need K-factors and tables for muon-beam-
related simulations. For an appreciable estimate of the precision of the
analytical results of the WW approximation integration, a comparison of the
single-differential cross section $d\sigma_{2\rightarrow
3}^{X}/dx\rvert_{\text{WW}}$ with both ETL and IWW results is shown in Fig. 5
for different $m_{X}$ in the $X=Z^{\prime}$ scenario. The relative error
between the computations at ETL and within the WW phase-space approximation is
$\leq 2\%$, while both the WW numerical and analytical methods agree well
within $<1\%$. The reader is referred to [29] and [30] for a more in-depth
discussion of the errors.

Figure 5: (Top) Single-differential cross-section $d\sigma_{2\rightarrow
3}^{X}/dx\rvert_{\text{approx.}}$ for the process $\mu N\rightarrow\mu
NZ^{\prime}$ in both ETL, WW and IWW approximation. In the WW approximation,
the results from the analytical expression of Eq. (8) and from the numerical
integration of Eq. (7) are shown. (Bottom) Relative error with respect to the
ETL computation, defined as
$(\mathcal{O}_{\text{approx.}}-\mathcal{O}_{\text{ETL}})/\mathcal{O}_{\text{ETL}}$.

As a final test, the implementation of the $Z^{\prime}$ vector boson emission
process in the mass limit $m_{Z^{\prime}}\rightarrow 0$ is compared to the
Geant4 treatment of SM muon bremsstrahlung, $\mu N\rightarrow\mu N\gamma$. The
comparison is driven against the absolute yields of both processes. As such,
those are extracted from a minimal Geant4 simulation of muons impinging on a
thick block of lead (Pb). The parameters and production cuts of the
G4MuBremsstrahlungModel are chosen to be comparable with the one of DarkZ,
namely similar values of $t_{\text{min}}$ and $t_{\text{max}}$, and equal
parameters $x_{\text{min}}$ and $x_{\text{max}}$ (respectively
$v_{\text{cut}}$ and $v_{\text{max}}$ within Geant4, see [31]), as well as a 1
GeV production cut (emission threshold). This later choice implies that the
cross sections do not depend on the $Z^{\prime}$ mass for values below 0.1
MeV. The results for mono-energetic 160 GeV muons with
$g_{Z^{\prime}}=\epsilon_{Z^{\prime}}e=1$ are shown in Fig. 6. Within the
statistical uncertainties, it is found that the absolute yields agree within
$1\%$.

Figure 6: Comparison of the absolute yield for both SM bremsstrahlung, $\mu
N\rightarrow\mu N\gamma$ (G4MuBremsstrahlung), and $Z^{\prime}-$strahlung,
$\mu N\rightarrow\mu NZ^{\prime}$ (DMProcessDMBrem), in the limit
$m_{Z^{\prime}}\rightarrow 0$ as simulated with mono-energetic muons impinging
on a thick block of lead. Both simulations are performed with Geant4, and
light dark matter is generated through DMG4. Similar parameters and production
cuts are chosen in both simulations.

###  VI.3 Final-state muon kinematics

As compared to the experiments in electron/positron beams, for the muon beam
the sampling of the recoil muon kinematic variables is much more important.
This is due to the fact that the missing energy in this case is determined by
measuring the recoil muon momentum by tracker. Similarly to Eq. (7), for the
recoil muon variables we have

$\frac{d^{2}\sigma^{X}_{2\to
3}}{dyd\cos\psi_{\mu}^{\prime}}\Bigr{|}_{\text{WW}}\simeq\frac{\alpha\chi^{\text{WW}}}{\pi(1-y)}E_{0}^{2}y\beta_{\mu}^{\prime}\frac{d\sigma^{X}_{2\to
2}}{d(pp^{\prime})}\Bigr{|}_{t=t_{\text{min}}},$ (9)

where $y$ is the final-state muon fractional energy, $\psi_{\mu}^{\prime}$ is
its emission angle, $\beta_{\mu}^{\prime}$ is its Lorentz $\beta-$factor and
$p^{\prime}$ is its four-momentum. The corresponding double differential cross
section is implemented in the package DarkMatter. For illustrative purposes,
the results of events sampling for both the $y$ and $\psi_{\mu}^{\prime}$
variables are shown Fig. 7, as extracted from a DMG4-based simulation of light
$Z^{\prime}$ production with fixed mass $m_{Z^{\prime}}=100$ MeV. The validity
of the implementation is inferred through a Kolmogorov-Smirnov (KS) test [32,
33] for the goodness of the fit, comparing both the target partial
distribution function (PDF) and the corresponding sampled histogram, and
indicating for the test mass range $m_{Z^{\prime}}=10-1000$ MeV that the
identity null hypothesis is accepted at a 0.05 significance level.

Figure 7: Comparison between the target partial distribution function and the
sampled events from the DarkMatter function SimulateEmissionByMuon2 for final-
state muons within a DMG4-based simulation of the process $\mu N\rightarrow\mu
NZ^{\prime}$. (Left) The fractional muon energy distribution, $y$. (Right) The
final-state muon emission angle distribution, $\psi_{\mu}^{\prime}$.

##  VII Developments in the annihilation processes

In this section, we describe the main DMG4 developments and improvements
concerning light dark matter production via resonant $e^{+}e^{-}$
annihilation. We refer the reader to our previous work [6] for a comprehensive
discussion regarding the general implementation of this process in the
package.

###  VII.1 Extension to different models

For the four mediator cases already implemented in the DMG4 package,
$e^{+}e^{-}\rightarrow X\rightarrow\chi\chi$, where $X$ can be a vector,
axial-vector, scalar, and pseudo-scalar particle, we introduced the
possibility for the final state DM particles $\chi$ to be (complex) scalars.
Specifically, we considered the following Lagrangians, where $\Phi$ is the LDM
complex scalar field, $\psi$ is a SM charged field. $S$, $P$, $A_{\mu}$,
$V_{\mu}$ are the mediators field, while
$\epsilon_{V},\epsilon_{A},\epsilon_{S},\epsilon_{P}$ are the mixing (or
coupling) parameters. Finally, $m_{V},m_{A},m_{S},m_{P}$ are the masses of the
mediators.
Vector case:

$\mathcal{L}\supset\mathcal{L}_{SM}-\frac{1}{4}V_{\mu\nu}^{2}+\frac{1}{2}m_{V}^{2}V_{\mu}^{2}-\sum_{\psi}e\epsilon_{V}V_{\mu}\bar{\psi}\gamma^{\mu}\psi+g^{D}_{V}V_{\mu}(\Phi^{*}\partial^{\mu}\Phi-(\partial^{\mu}\Phi)\Phi^{*})+\frac{1}{2}(\partial_{\mu}\Phi)^{2}-\frac{1}{2}m^{2}_{\Phi}\Phi^{*}\Phi$
(10)

Axial vector case:

$\mathcal{L}\supset\mathcal{L}_{SM}-\frac{1}{4}A_{\mu\nu}^{2}+\frac{1}{2}m_{A}^{2}A_{\mu}^{2}-\sum_{\psi}e\epsilon_{A}A_{\mu}\bar{\psi}\gamma_{5}\gamma^{\mu}\psi+g^{D}_{A}A_{\mu}(\Phi^{*}\partial^{\mu}\Phi-(\partial^{\mu}\Phi)\Phi^{*})+\frac{1}{2}(\partial_{\mu}\Phi)^{2}-\frac{1}{2}m^{2}_{\Phi}\Phi^{*}\Phi$
(11)

Scalar case:

$\mathcal{L}\supset\mathcal{L}_{SM}+\frac{1}{2}(\partial_{\mu}S)^{2}-\frac{1}{2}m^{2}_{S}S^{2}-\sum_{\psi}e\epsilon_{S}S\bar{\psi}{\psi}+g^{D}_{S}m_{S}S\Phi^{*}\Phi+\frac{1}{2}(\partial_{\mu}\Phi)^{2}-\frac{1}{2}m^{2}_{\Phi}\Phi^{*}\Phi$
(12)

Pseudo-scalar case:

$\mathcal{L}\supset\mathcal{L}_{SM}+\frac{1}{2}(\partial_{\mu}P)^{2}-\frac{1}{2}m^{2}_{P}P^{2}-\sum_{\psi}e\epsilon_{P}P\bar{\psi}\gamma^{5}{\psi}+g^{D}_{P}m_{P}P\Phi^{*}\Phi+\frac{1}{2}(\partial_{\mu}\Phi)^{2}-\frac{1}{2}m^{2}_{\Phi}\Phi^{*}\Phi\;\;.$
(13)

The corresponding cross sections for the $e^{+}e^{-}\rightarrow
X\rightarrow\Phi\Phi$ processes implemented in the code reads

$\sigma_{e^{+}e^{-}}=\frac{4\pi\alpha_{EM}\alpha_{D}\varepsilon^{2}}{\sqrt{s}}q\frac{\mathcal{K}}{(s-m_{X}^{2})^{2}+\Gamma_{X}^{2}m_{X}^{2}}\;\;$
(14)

where $s$ is the invariant mass of the $e^{+}e^{-}$ system, $m_{X}$ the
mediator mass, $q=\frac{\sqrt{s}}{2}\sqrt{1-\frac{4m_{\Phi}^{2}}{s}}$,
$\Gamma_{X}$ is the intermediate DM particle decay width to dark particles
$\Phi$, $\alpha_{EM}$ is the electromagnetic fine structure constant, and
$\alpha_{D}\equiv\frac{\left(g^{D}_{X}\right)^{2}}{4\pi}$ is the $\Phi$
coupling squared to the dark particles $\chi$. Finally, $\mathcal{K}$ is a
kinematic factor that reads, respectively, $2/3q^{2}$ for the vector and
axial-vector mediator, and $m^{2}_{X}/4$ for the scalar and pseudo-scalar
cases. Finally, the $X$ decay widths read:

$\displaystyle\Gamma_{V\rightarrow\Phi^{*}\Phi}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{D}}{12}m_{V}\bigl{(}1-4r^{2}\bigr{)}^{\frac{3}{2}}$
(15) $\displaystyle\Gamma_{A\rightarrow\Phi^{*}\Phi}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{D}}{12}m_{A}\bigl{(}1-4r^{2}\bigr{)}^{\frac{3}{2}}$
(16) $\displaystyle\Gamma_{S\rightarrow\Phi^{*}\Phi}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{D}}{4}m_{S}\bigl{(}1-4r^{2}\bigr{)}^{1/2}$ (17)
$\displaystyle\Gamma_{P\rightarrow\Phi^{*}\Phi}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{D}}{4}m_{S}\bigl{(}1-4r^{2}\bigr{)}^{1/2},$ (18)

where $r=m_{\Phi}/m_{X}$.

For the vector-mediator model, we also considered an inelastic (i.e., non-
diagonal) model foreseeing two fermionic LDM states $\chi_{2}$ and $\chi_{1}$
with different masses $m_{2}$ and $m_{1}$, with $\Delta=m_{2}-m_{1}>0$ [20].
In this model, at tree level only the $A^{\prime}-\chi_{2}-\chi_{1}$ vertex
exists. The corresponding Lagrangian reads:

$\mathcal{L}\supset\mathcal{L}_{SM}-\frac{1}{4}V_{\mu\nu}V^{\mu\nu}+\frac{1}{2}m^{2}_{V}V^{2}-\sum_{\psi}e\varepsilon_{V}V^{\mu}\bar{\psi}\gamma_{\mu}{\psi}+\sum_{i}\bar{\chi}_{i}(\not{\partial}-m_{i})\chi_{i}+g_{D}A^{\mu}(\bar{\chi}_{2}\gamma_{\mu}\chi_{1}+\bar{\chi_{1}}\gamma_{\mu}\chi_{2})$
(19)

The corresponding annihilation cross section reads:

$\sigma_{e^{+}e^{-}}=\frac{8\pi\varepsilon_{V}^{2}\alpha_{D}\alpha_{EM}}{\sqrt{s}}\frac{q}{(s-m^{2}_{V})^{2}+\Gamma^{2}m^{2}_{V}}(m_{1}m_{2}+E_{1}E_{2}-\bar{p}^{2}/3)\;\;,$
(20)

where $q$ is the $\chi_{2}$ and $\chi_{1}$ three-momentum in the $V$ rest
frame, while $E_{2}$ and $E_{1}$ are the final state particles total energies
in the same frame. The decay width $\Gamma$ for the
$V\rightarrow\chi_{2}\chi_{1}$ channel is:

$\Gamma_{V\rightarrow\chi_{2}\chi_{1}}=\frac{\alpha_{D}}{3m^{2}_{V}}2q\left(3m_{1}m_{2}+m^{2}_{V}-\frac{m^{2}_{1}+m^{2}_{2}}{2}-\frac{(m^{2}_{2}-m^{2}_{1})^{2}}{2m^{2}_{V}}\right)$
(21)

In this DMG4 software version, we also included the possibility to simulate
the $e^{+}e^{-}\rightarrow Z^{\prime}\rightarrow\bar{\nu}\nu$ process, where
the intermediate $Z^{\prime}$ particle is the massive force mediatior
associated with a new SM gauge symmetry involving the neutrino sector.
Specifically, we considered the two cases of a $B-L$ gauge symmetry and a
$L_{\mu}-L_{\tau}$ one [34, 35]. While the first case resembles that of a
vector dark photon by making the substitution $\varepsilon e\leftrightarrow
g$, where $g_{Z^{\prime}}$ is the gauge coupling constant, for the second
scenario at tree level the $Z^{\prime}$ couples only to second and third
generation leptons. Nevertheless, thanks to the photon-$Z^{\prime}$ coupling
introduced at next-to-leading order by loops involving muon and tau leptons,
an effecting $Z^{\prime}-e^{+}-e^{-}$ vertex appears, with coupling
$e\Pi(q^{2})$, where $q^{2}$ is the $Z^{\prime}$ four-momentum squared. The
function $\Pi$ reads:

$\Pi(q^{2})=\frac{e\,g_{Z}^{\prime}}{2\pi^{2}}\int_{0}^{1}dx\,x(1-x)\ln\frac{m^{2}_{\tau}-x(1-x)q^{2}}{m^{2}_{\mu}-x(1-x)q^{2}}\;\;.$
(22)

In the code, we were able to introduce the full analytical expression for the
$\Pi$ function, allowing to optimize the corresponding computation time. This
is given by the formula:

$\begin{split}\frac{2\pi^{2}}{eg_{Z^{\prime}}}\Pi\left(q^{2}\right)=\frac{1}{3}\left[\frac{1}{2}\log\left(\frac{r_{\tau}}{r_{\mu}}\right)+2\left(r_{\mu}-r_{\tau}\right)+\right.\\\
-\left(1+2r_{\mu}\right)\sqrt{1-4r_{\mu}}\operatorname{coth}^{-1}\left(\sqrt{1-4r_{\mu}}\right)+\\\
\left.+\left(1+2r_{\tau}\right)\sqrt{1-4r_{\tau}}\operatorname{coth}^{-1}\left(\sqrt{1-4r_{\tau}}\right)\right],\quad\end{split}$
(23)

with $r_{\mu}=m_{\mu}^{2}/q^{2}$ and $r_{\tau}=m_{\tau}^{2}/q^{2}$.

###  VII.2 Narrow-width resonance production

The introduction of all new DMG4 processes in Geant4 is based on the
G4VDiscreteProcess class - each process is described by a new class inheriting
from the latter. In particular, each class has to provide a concrete
implementation of the two pure virtual methods GetMeanFreePath and
PostStepDoIt; in particular, the first method is responsible of computing the
positron mean free path for the annihilation process at the beginning of each
Geant4 step, when the step-length is not yet determined, from the knowledge of
the corresponding G4Track. It follows that, if the cross section for the
annihilation process is computed solely from the value of the positron energy
at the beginning of the step $E_{i}$, without accounting for the full cross-
section energy dependence, the obtained result would under-estimate (over-
estimate) the average value along the step if $E_{i}>E_{\text{MAX}}$
($E_{i}<E_{\text{MAX}}$), where $E_{\text{MAX}}$ is the positron resonant
energy in the laboratory frame. While in the case of a broad resonance this
effect is not critical, if the typical energy loss across a Geant4 step is
significantly larger than the mediator width $\Gamma$, it could induce a
significant distortion of the simulation results in case of a narrow resonance
- this is the case, for example, of the $Z^{\prime}$ model described
previously in the text. Two ad-hoc modifications were introduced in the code
to account for this effect.

####  VII.2.1 Energy dependence of the cross section

As discussed previously, the G4VDiscreteProcess-based implementation requires
to compute the cross section for the annihilation process at the beginning of
each Geant4 step, before the latter is determined. To account for the energy
dependency of the cross section, we followed the strategy implemented in
Geant4 for other electromagnetic-like processes and described in the manual,
based on the sole assumption that the cross section $\sigma(E)$ has only one
maximum $\sigma_{\text{MAX}}$, for $E=E_{\text{MAX}}$. Specifically, our
implementation assumes that, during the step still to be computed, the cross
section $\sigma(E)$ will be lower than a certain value $\sigma_{m}$. The mean
free path for the process is then computed via $\sigma_{m}$. If this process
is selected by Geant4 as that taking place and limiting the next step length,
the final state sampling, implemented via the PostStepDoIt method, is invoked
with reduced probability $p_{red}=\sigma(E_{f})/\sigma_{m}$, where $E_{f}$ is
the positron energy at the end of the step; alternatively, nothing happens and
the positron tracking continues. In the following we call this procedure
“Maximum-Rejection”. We underline that it is not obvious that in this
procedure $p_{red}$ should be computed considering the cross section at the
end of the step and not the average value along the step111Incidentally, we
note that in the Geant4 physics reference manual for release 11.2, rev 8.0, in
Sec. 7.4 “Correcting the Cross Section for Energy Variation” the prescription
to use of $\sigma(E_{f})$ in the computation of $p_{red}$ is supported by a
reference to [39]. However, in that work no justification of this method is
provided.. To empirically verify this, we made a dedicated toy-MC study,
presented in Appendix A.

The value $\sigma_{m}$ is selected as follows, dividing the positron energy
range in three dinstinct intervals, also shown in Fig. 8. The first case
corresponds to the condition $E_{i}<E_{MAX}$: in this case, the positron
energy is already lower than the cross section maximum, and $\sigma(E)$ will
further decrease along the positron trajectory; therefore,
$\sigma_{m}=\sigma(E_{i})$. In the second case, instead,
$E_{MAX}<E_{i}<E_{MAX}/\xi$, where $0<\xi<1$ is a dimensionless parameter to
be appropriately tuned; this case corresponds to a positron with initial
energy close to that of the resonance, and assumes $\sigma_{m}=\sigma_{MAX}$.
Finally, if $E_{i}>E_{MAX}/\xi$, then $\sigma_{m}=\sigma(\xi E_{i})$.

Figure 8: Cross section $\sigma(E_{+})$ for the annihilation process
$e^{+}e^{-}\rightarrow X\rightarrow\chi\chi$ for a narrow-width mediator. The
figure highlights the three regions (I), (II), and (III) in which the energy
interval is divided into to handle the energy dependence of the cross section.

The parameter $\xi$ controls the transition between case (II) and case (III),
and is connected to the typical positron energy loss across a step. If the
typical value of the continuous energy loss across a step is $\Delta E$, then
the value of $\xi$ can be set from the requirement that, for a positron
initial energy $E_{i}$ in the interval $[E_{R},E_{i}-n\Delta_{E}]$, where
$n>1$ is a safety factor, case (II) should be applied, since it is possible
that, in the next step, the particle will loose enough energy to cross the
resonant energy. All values of $\xi$ in the interval
$[\frac{E_{R}}{E_{R}+n(\Delta E)_{MAX}},1]$ satisfy this requirement. In the
code, the lower limit of the interval is used.

####  VII.2.2 Dynamical step limitation

The approach presented in the previous section only considers continuous
energy losses for the positron, but not the discrete ones associated to
bremsstrahlung emission, that may significantly alter the positron energy.
Also, this approach is expected to work well for geometries with many
geometrical boundaries, forcing many positron steps, but can be inaccurate for
large homogeneous volumes. To partially compensate for this, we implemented a
step limitation for positrons. The intrinsic scale of the positron step
limitation $\delta_{0}$ for narrow resonance production is given by the ratio
of the modified resonance width $\Gamma^{*}=\Gamma M_{X}/m_{e}$, where the
factor $M_{X}/m_{e}$ is a kinematic correction induced by the Lorentz
transformation from the center-of-mass frame to the laboratory one, and the
positron energy loss per unit length $-\frac{dE}{dx}$,
$\delta_{0}=-\Gamma^{*}/\frac{dE}{dx}$.

In order to optimize the computation time, we use a dynamical step limitation
for positrons, based on the relative cross section variation with respect to
the energy: the more the cross section is expected to vary for a given
positron energy, the smaller the maximum step allowed for the positron. The
resonant annihilation cross section can indeed be written as
$\sigma=f(E_{+})\frac{1}{(E-E_{MAX})^{2}+(\Gamma^{*})^{2}}$, where $f$ is a
smooth-varying function of the positron energy. It follows that the relative
cross section variation with respect to the energy reads:

$V(E)=\frac{1}{\sigma}|\frac{d\sigma}{dE}|\simeq\frac{2|E-E_{R}|}{(E-E_{R})^{2}+(\Gamma^{*})^{2}}\;\;.$
(24)

This function has two maxima at $E=E_{R}\pm\Gamma^{*}$, where
$V=1/\Gamma^{*}$. To avoid numerical instabilities in the interval between
these two values, associated with the zero at $E=E_{R}$, we slightly modified
it to be $V^{\prime}=V$ if $|E-E_{R}|>\Gamma^{*}$, and
$V^{\prime}=1/\Gamma^{*}$ otherwise. In the code, we thus define the positron
step limit as $\delta_{MAX}=\delta_{0}\cdot 1/V^{\prime}\cdot(1/\Gamma^{*})$.

###  VII.3 Atomic effects in $e^{+}e^{-}$ annihilation

In the previous DMG4 version, all $e^{+}e^{-}$ annihilation processes were
handled by considering the atomic electron to be at rest in the laboratory
frame. Atomic motions reflect to a broadening of the annihilation cross
section $\sigma(E_{+})$, since a given positron energy can reflect to
different values of the Mandelstam variable $s$ for the electron-positron
system, defined by the equation:

$s=2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})\;\;,$ (25)

where $(E_{-})$ ($P_{-}$) is the atomic electron total energy (three-
momentum), and $z$ is the cosine of the angle between the $e^{+}$ and $e^{-}$
momenta. From this, it follows that a positron whose energy won’t match the
resonance mass if annihilating with an at-rest electron can do so for an
orbiting one and vice-versa, resulting to a broad range $\Delta_{E}$ of
energies that can contribute to the process. To illustrate this effect, we
consider the case of a $m_{X}=225$ MeV mediator, and an electron kinetic
energy of about 10 keV (typical value for lead): the positron energy in the
laboratory frame resulting to $s=m_{X}^{2}$ spans a range from 40.7 GeV to
60.3 GeV (for comparison, the at-rest value is 48.6 GeV). Similarly, for
$m_{X}=17.7$ MeV, considering an electron kinetic energy of about 290 eV
(typical value for carbon), the resonant positron energy in the laboratory
frame spans from 273 MeV to 291 MeV. From this consideration, it follows that
the at-rest approximation is justified for large-width resonances, when
$\Gamma_{X}\frac{m}{m_{e}}\gtrsim\Delta_{E}$, while it can’t be applied in the
narrow width scenario.

In this DMG4 software version, we implemented a full calculation involving
atomic motions as follows. For the simpler case of massless final state
particles, such as in the reaction $e^{+}e^{-}\rightarrow
Z^{\prime}\rightarrow\nu\bar{\nu}$, the general expression for the cross
section can be written as
$\sigma=A(E_{+})\frac{1}{(s-m^{2}_{X})^{2}+m^{2}_{X}\Gamma^{2}}$, where $A$ is
a smooth-varying function of the impinging positron energy, and $m_{X}$,
$\Gamma$ are, respectively, the mediator mass and total width. By substituting
the previous expression for $s$, this translates to:

$\sigma=f(E_{+})\frac{1}{(2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})-m^{2}_{X})^{2}+m^{2}_{X}\Gamma^{2}}\;\;.$
(26)

The code computes, for each interaction, the average value of this expression,
by considering an isotropic electron motion, i.e. a flat PDF for $z$ in the
$[-1,1]$ interval:

$\displaystyle\bar{\sigma}$
$\displaystyle=\int_{-1}^{1}dz\,p(z)\,\frac{A(E_{+})}{(2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})-m^{2}_{X})^{2}+m^{2}_{X}\Gamma^{2}}=$
(27)
$\displaystyle=\frac{A(E_{+})}{4m_{X}\Gamma}\left[\arctan{\frac{m^{2}_{X}+2P_{+}P_{-}-2m^{2}_{e}-2E_{+}E_{-}}{m_{X}\Gamma}}-\arctan{\frac{m^{2}_{X}-2P_{+}P_{-}-2m^{2}_{e}-2E_{+}E_{-}}{m_{X}\Gamma}}\right]\;\;,$
(28)

where $E_{+}$ and $E_{-}$ ($P_{+}$ and $P_{-}$) are the positron and the
electron total energies (momenta) in the laboratory frame, and
$p(z)=\frac{1}{2}$ within the integration interval.

Concerning the atomic electron energies, the code currently implements two
possible models. The first, simplified model assumes, for each atomic shell, a
fixed kinetic energy equal to the binding energy magnitude, $T_{-}=|B|$, and
computes the average cross section by summing over the atomic shells,
considering for each a weight proportional to the occupation number. The
shells atomic energies for the target material are obtained directly from
Geant4, using the G4Element::GetAtomicShell method. Similarly, the occupation
numbers are obtained from the G4Element::GetNbOfShellElectrons method. In the
second, more realistic model, used by default, an exponential PDF distribution
222We incidentally note that this model is the same adopted by the Geant4
software to handle atomic motion effects in Compton scattering. for $T_{-}$ is
considered, $p(T_{-})=\frac{1}{|B|}\exp(-T_{-}/|B|)$. The contribution for
each shell is separately computed, by integrating the expression for
$\bar{\sigma}$ over $p(T_{-})$, and summing the contributions with weights
each equal to the shell occupation number. The code computes the integral over
$p(T_{-})$ numerically through the Monte Carlo method. To speed-up the
computation time, since the expression for $p(T_{-})$ is independent from the
positron kinematics, a set of random electron kinetic energies ${T_{-}}$ is
pre-generated for each atomic shell at initialization time. The effect of the
inclusion of atomic motion description in the software is illustrated in Fig.
9, showing the annihilation cross section as a function of the positron energy
for a 225 MeV resonance, having a width of 7 keV. The black curve is the
result obtained without considering atomic motions, while the red (green)
curve are the results obtained from the first and second model described
before. In the first case, the discreteness of the atomic shells kinetic
energies is clearly visible.

Figure 9: Cross section $\sigma(E_{+})$ for the annihilation process
$e^{+}e^{-}\rightarrow Z^{\prime}\rightarrow\nu\bar{\nu}$ for a 225 MeV
$Z^{\prime}$ with $\Gamma=7$ keV, decaying into neutrinos. The black curve is
the result obtained without including atomic effects, while the red (green)
curve are the results obtained from the simplified (realistic) atomic model
currently implemented in DMG4.

This approach is slightly modified for the more general case involving equally
massive particles in the final state, such as in the LDM production process
$e^{+}e^{-}\rightarrow X\rightarrow\chi\chi$, since in this case the function
$A(E_{+})$ is proportional to the final-state momentum $q$, that intrinsically
contains a kinematic threshold333The extension of the code to handle final
state with particles having different masses is foreseen for the next DMG4
release.. Starting from the previous for $\bar{\sigma}$ in Eq. 27, and writing
explicitly the $q$ factor, the expression for $\bar{\sigma}$ is:

$\bar{\sigma}=\int_{z_{min}}^{z_{max}}dz\,p(z)\,\frac{qB(E_{+})}{(2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})-m^{2}_{X})^{2}+m^{2}_{X}\Gamma^{2}}\;\;,$
(29)

where $B(E_{+})\equiv A(E_{+})/q$. The integration range can more easily be
determined by changing the integration variable from $z$ to $s$ for fixed
$E_{+}$, and explicitly setting $p(z)=\frac{1}{2}$:

$\bar{\sigma}=\frac{1}{4E_{+}P_{-}}\int_{s(z=z_{max})}^{s(z=z_{min})}ds\frac{B(s)q(s)}{(s-M^{2})^{2}+M^{2}\Gamma^{2}}\;\;,$
(30)

and considering the threshold requirement $s>s_{min}=4m^{2}_{\chi}$. From the
identity $s=2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})$, it follows that
$z<z_{max}=\frac{2m^{2}_{e}+2E_{+}E_{-}-s_{min}}{2E_{+}P_{-}}$. Note that if
$z_{max}<-1$, then the positron energy is incompatible with $s_{min}$ for all
angles, and thus $\bar{\sigma}=0$. If, instead, $z_{max}>1$, the integral is
still performed only in the physical region up to $z_{max}=+1$. The result of
the integral, assuming a constant value for $B$ and defining
$x=\sqrt{s-4m^{2}_{\chi}}$ is:

$\bar{\sigma}=\frac{B}{8E_{+}P_{-}}\cdot\left[P(x(s(z=-1)))-P(x(s(z=z_{max})))\right]\;\;,$
(31)

where $P(x)$ is given by:

$\displaystyle P(x)=$
$\displaystyle\left[\frac{1}{2\sqrt{2(\eta^{2}+\Delta^{2})}}\log{\frac{x^{2}-\sqrt{2(\eta^{2}+\Delta^{2})}x+\eta^{2}}{x^{2}+\sqrt{2(\eta^{2}+\Delta^{2})}x+\eta^{2}}}\right.+$
(32)
$\displaystyle+\frac{1}{2\sqrt{\eta^{2}/2-\Delta^{2}/2}}\left.\left(\arctan{\frac{2x-\sqrt{2(\eta^{2}+\Delta^{2})}}{\sqrt{2\eta^{2}-2\Delta^{2}}}}+\arctan{\frac{2x+\sqrt{2(\eta^{2}+\Delta^{2})}}{\sqrt{2\eta^{2}-2\Delta^{2}}}}\right)\right]\;\;,$
(33)

where $\Delta^{2}=m_{X}^{2}-4m^{2}_{\chi}$ and
$\eta^{4}=\Delta^{4}+\Gamma^{2}m^{2}_{X}$. Finally, to handle the residual
$s-$ (and thus $z-$) dependency in $B$ for fixed values of $E_{+}$ and
$E_{-}$, this term is computed event-by-event via Monte Carlo integration,
taking 20 uniformly-spaced $z-$ values in the allowed interval and averaging
over the corresponding $B$ values. As an example, Fig. 10 shows the obtained
cross section, considering the more realistic model for atomic electrons
energies, for a 270 MeV mediator with $\Gamma=1$ keV, decaying into 90 MeV
final state particles. In this configuration, for the annihilation with at-
rest electrons, the positron threshold energy would be approximately 25 GeV,
but thanks to atomic motions also lower-energy positrons can undergo this
reaction.

Figure 10: Cross section $\sigma(E_{+})$ for the annihilation process
$e^{+}e^{-}\rightarrow X\rightarrow\chi\chi$ for a 240 MeV mediator with
$\Gamma\simeq 1$ keV, decaying into 80 MeV final state particles, after
including atomic effects. The black curve is the result obtained without
considering atomic effects, the red curve is the result obtained ignoring the
$q$ dependency of the cross section, and finally the green curve is the full
result accounting for the kinematic threshold. Electron atomic energies were
simulated using the realistic model previously described.

##  VIII Conclusion

The package DMG4 for the simulation of light dark matter production in fixed
target experiments is created. It can be used in full simulation programs
based on the Geant4 framework. The subpackage DarkMatter containing a
collection of cross sections is weakly connected to Geant4, thus can be used
in programs based on other frameworks. The package is widely used in the NA64
experiment.

In this paper, numerous developments of the the package are presented, in
particular the WW approximation cross sections for the muon beam [30], models
with semivisible $A^{\prime}$ (inelastic Dark Matter) [18] and several
important improvements for the annihilation processes.

The package is available at http://mkirsano.web.cern.ch/mkirsano/DMG4.tar.gz
and https://gitlab.cern.ch/P348/DMG4 (internal NA64 repository). It is
recommended also to contact the corresponding author Mikhail Kirsanov about
the usage.

##  IX Acknowledgements

This work was supported by the Ministry of Science and Higher Education (MSHE)
and RAS (Russia), MSHE Project FSWW-20230003, Tomsk Polytechnic University
within the assignment of MSHE (Russia), the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (Grant
agreement No. 947715 - POKER Starting Grant). This work is partially supported
by ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data
and Quantum Computing, funded by European Union – NextGenerationEU. This work
is also supported by RyC-030551-I, PID2021-123955NA-100 and CNS2022-135850
funded by MCIN/AEI/FEDER, UE (Spain) and by ETH Zurich and SNSF Grants No.
186181, 186158, 197346, 216602 (Switzerland).

##  X Appendix A

In order to check the the procedure “Maximum-Rejection” we considered a simple
model where there is a pseudo-positron particle moving forward in a material
and loosing a constant, deterministic amount of energy $\alpha$ per unit
length. Therefore, if the pseudo-positron pseudo-positron enters in this
material at $z=0$ with initial energy $E_{0}$, the energy at depth $z$ in the
material is given by $E(z)=E_{0}-\alpha z$, for $0<z<E_{0}/\alpha$. The only
other process available for this particle is the annihilation process, with a
simplified cross section formula
$\sigma(E)=\sigma_{0}\Gamma^{2}/((E-E_{R})^{2}+\Gamma^{2})$.

If a mono-chromatic beam of pseudo-positrons impinges on the material at
$z=0$, the number of particles at a given depth in the material $N(z)$ can be
computed from the following equation:

$\frac{dN}{dz}=-N(z)k\sigma_{0}\frac{\Gamma^{2}}{(E(z)-E_{R})^{2}+\Gamma^{2}}=-N(z)k\sigma_{0}\frac{\Gamma^{2}}{(E_{0}-\alpha
z-E_{R})^{2}+\Gamma^{2}}\;\;,$ (34)

where $k$ is the number density of electrons per unit volume in the material.
This equation can be solved exactly, obtaining:

$N(z)=N_{0}\cdot
e^{-\frac{k\sigma_{0}\Gamma}{\alpha}atan\frac{E_{0}-E_{R}}{\Gamma}}\cdot
e^{-\frac{k\sigma_{0}\Gamma}{\alpha}atan\frac{E_{0}-E_{R}-\alpha z}{\Gamma}}$
(35)

This gives a number of interactions per unit length, $-\frac{dN}{dz}$, equal
to:

$-\frac{dN}{dz}=N_{0}\cdot
e^{-\frac{k\sigma_{0}\Gamma}{\alpha}atan\frac{E_{0}-E_{R}}{\Gamma}}\cdot
e^{-\frac{k\sigma_{0}\Gamma}{\alpha}atan\frac{E_{0}-E_{R}-\alpha
z}{\Gamma}}\cdot k\sigma_{0}\frac{1}{1+\left(\frac{E_{0}-E_{R}-\alpha
z}{\Gamma}\right)^{2}}.$ (36)

To check the “Maximum-Rejection” procedure we made a toy Monte Carlo
simulation, accounting for the energy variation across the step as discussed
in Sec. VII.2.1. Specifically, we implemented a first version computing the
reduced probability $p_{red}$ through the cross section at the end of the
step, as recommended in the Geant4 manual, and also by using the mean cross
section value along the step. The result is shown in Fig. 11. This test
confirms the solution suggested by the Geant4 manual.

Figure 11:  The result of checking the “Maximum-Rejection” recipe. Green curve
- exact solution, blue points with error bars - using $E_{f}$ as specified in
the recipe, red points with error bars - using mean cross section along the
step instead.

## References

  * Bjorken  _et al._ [2009] J. D. Bjorken, R. Essig, P. Schuster, and N. Toro, Phys. Rev. D 80, 075018 (2009), arXiv:0906.0580 [hep-ph] .
  * Izaguirre  _et al._ [2013] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D 88, 114015 (2013), arXiv:1307.6554 [hep-ph] .
  * Batell  _et al._ [2009] B. Batell, M. Pospelov, and A. Ritz, Phys. Rev. D 80, 095024 (2009), arXiv:0906.5614 [hep-ph] .
  * Izaguirre  _et al._ [2015] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D 91, 094026 (2015), arXiv:1411.1404 [hep-ph] .
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Bondi  _et al._ [2021] M. Bondi, A. Celentano, R. R. Dusaev, D. V. Kirpichnikov, M. M. Kirsanov, N. V. Krasnikov, L. Marsicano, and D. Shchukin, Comput. Phys. Commun. 269, 108129 (2021), arXiv:2101.12192 [hep-ph] .
  * Holdom [1986] B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Battaglieri  _et al._ [2017] M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter_ (2017) arXiv:1707.04591 [hep-ph] .
  * Beacham  _et al._ [2020] J. Beacham _et al._ , J. Phys. G 47, 010501 (2020), arXiv:1901.09966 [hep-ex] .
  * ISO Central Secretary [2012] ISO Central Secretary,  _Information technology – Object Management Group Unified Modeling Language (OMG UML) – Part 1: Infrastructure_, Standard ISO/IEC 19505-1:2012 (International Organization for Standardization, Geneva, CH, 2012).
  * [11] Monte Carlo Particle Numbering Scheme, https://pdg.lbl.gov/2019/reviews/rpp2019-rev-monte-carlo-numbering.pdf, accessed: 2021-01-25.
  * Tucker-Smith and Weiner [2001] D. Tucker-Smith and N. Weiner, Phys. Rev. D 64, 043502 (2001), arXiv:hep-ph/0101138 .
  * Dusaev  _et al._ [2020] R. R. Dusaev, D. V. Kirpichnikov, and M. M. Kirsanov, Phys. Rev. D 102, 055018 (2020), arXiv:2004.04469 [hep-ph] .
  * Marsicano  _et al._ [2018] L. Marsicano, M. Battaglieri, M. Bondí, C. D. R. Carvajal, A. Celentano, M. De Napoli, R. De Vita, E. Nardi, M. Raggi, and P. Valente, Phys. Rev. Lett. 121, 041802 (2018), arXiv:1807.05884 [hep-ex] .
  * Voronchikhin and Kirpichnikov [2022] I. V. Voronchikhin and D. V. Kirpichnikov, Physical Review D 106, 10.1103/physrevd.106.115041 (2022).
  * Lee  _et al._ [2014] H. M. Lee, M. Park, and V. Sanz, Eur. Phys. J. C 74, 2715 (2014), arXiv:1306.4107 [hep-ph] .
  * Kang and Lee [2020] Y.-J. Kang and H. M. Lee, The European Physical Journal C 80, 10.1140/epjc/s10052-020-8153-x (2020), arXiv:2001.04868 .
  * Filimonova  _et al._ [2022] A. Filimonova, S. Junius, L. L. Honorez, and S. Westhoff, Journal of High Energy Physics 2022, 10.1007/jhep06(2022)048 (2022).
  * Mohlabeng [2019] G. Mohlabeng, Phys. Rev. D 99, 115001 (2019), arXiv:1902.05075 [hep-ph] .
  * Cazzaniga  _et al._ [2021] C. Cazzaniga _et al._ (NA64), Eur. Phys. J. C 81, 959 (2021), arXiv:2107.02021 [hep-ex] .
  * He  _et al._ [1991a] X. G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 43, 22 (1991a).
  * He  _et al._ [1991b] X.-G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 44, 2118 (1991b).
  * Foot  _et al._ [1994] R. Foot, X. G. He, H. Lew, and R. R. Volkas, Phys. Rev. D 50, 4571 (1994), arXiv:hep-ph/9401250 .
  * Altmannshofer  _et al._ [2016] W. Altmannshofer, S. Gori, S. Profumo, and F. S. Queiroz, JHEP 12, 106, arXiv:1609.04026 [hep-ph] .
  * Kile  _et al._ [2015] J. Kile, A. Kobach, and A. Soni, Phys. Lett. B 744, 330 (2015), arXiv:1411.1407 [hep-ph] .
  * Park  _et al._ [2016] J.-C. Park, J. Kim, and S. C. Park, Phys. Lett. B 752, 59 (2016), arXiv:1505.04620 [hep-ph] .
  * Liu  _et al._ [2017] Y.-S. Liu, D. McKeen, and G. A. Miller, Phys. Rev. D 95, 036010 (2017), arXiv:1609.06781 [hep-ph] .
  * Liu and Miller [2017] Y.-S. Liu and G. A. Miller, Phys. Rev. D 96, 016004 (2017), arXiv:1705.01633 [hep-ph] .
  * Kirpichnikov  _et al._ [2021] D. V. Kirpichnikov, H. Sieber, L. M. Bueno, P. Crivelli, and M. M. Kirsanov, Phys. Rev. D 104, 076012 (2021), arXiv:2107.13297 [hep-ph] .
  * Sieber  _et al._ [2023] H. Sieber, D. V. Kirpichnikov, I. V. Voronchikhin, P. Crivelli, S. N. Gninenko, M. M. Kirsanov, N. V. Krasnikov, L. Molina-Bueno, and S. K. Sekatskii, Phys. Rev. D 108, 056018 (2023), arXiv:2305.09015 [hep-ph] .
  * Bogdanov  _et al._ [2006] A. G. Bogdanov, H. Burkhardt, V. N. Ivanchenko, S. R. Kelner, R. P. Kokoulin, M. Maire, A. M. Rybin, and L. Urban, IEEE Trans. Nucl. Sci. 53, 513 (2006).
  * Kolmogorov [1933] A. Kolmogorov, G. Ist. Ital. Attuari. 4, 83 (1933).
  * Smirnov [1948] N. Smirnov, Ann. Math. Statist. 19, 279 (1948).
  * Andreev  _et al._ [2022a] Y. M. Andreev _et al._ (NA64), Phys. Rev. D 106, 032015 (2022a), arXiv:2206.03101 [hep-ex] .
  * Andreev  _et al._ [2022b] Y. M. Andreev _et al._ (NA64), Phys. Rev. Lett. 129, 161801 (2022b), arXiv:2207.09979 [hep-ex] .
  * Note [1] Incidentally, we note that in the Geant4 physics reference manual for release 11.2, rev 8.0, in Sec. 7.4 “Correcting the Cross Section for Energy Variation” the prescription to use of $\sigma(E_{f})$ in the computation of $p_{red}$ is supported by a reference to [39]. However, in that work no justification of this method is provided.
  * Note [2] We incidentally note that this model is the same adopted by the Geant4 software to handle atomic motion effects in Compton scattering.
  * Note [3] The extension of the code to handle final state with particles having different masses is foreseen for the next DMG4 release.
  * Ivanchenko  _et al._ [1991] V. N. Ivanchenko, A. D. Bukin, M. S. Dubrovin, S. I. Eidelman, N. A. Grozina, and V. A. Tayursky, in  _Workshop on Detector and Event Simulation in High-energy Physics (MC ’91)_ (1991) pp. 79–85.