\section*{Acknowledgments}
\noindent
L.E.Y.-S.  acknowledges the financial support
by the Vicerrectoria de Investigaciones, Universidad de Cartagena, Colombia  through Project No. 019-2021.
 P.M. acknowledges the financial support
 by the Council of Scientific and Industrial Research (CSIR), India under Grant No. 03(1422)/18/EMR-II, and Science and Engineering Research Board (SERB), India under Grant No. CRG/2019/004059. A.B. acknowledges funding provided by the Institute of Physics Belgrade, through a grant by the Ministry of Education, Science, and Technological Development of the Republic of Serbia.
S.K.A. acknowledges support by the CNPq, Brazil through  grant 301324/2019-0.


\section*{Distinguishing simple and complex contagion processes on networks: Supplementary Material}

\section{Data description and generation of interaction networks and hypergraphs}
\label{SMsec_data}

The data sets we consider describe interactions between individuals, either in terms of physical proximity or online, and are temporally resolved.
We thus preprocess these data in order to build the  three types of structures needed for the numerical simulations: the temporally aggregated pairwise network, the aggregated hypergraph with all orders of interactions, and the aggregated triangle network.

First we build the aggregated network: we consider each interaction between two individuals as a pairwise link of a static network with weight given by the number of times (consecutive or not) that the connection has appeared in the temporal data. The aggregated network is used to simulate the simple and threshold contagion models, and the contagion events of the simplicial contagion model involving pairwise links.

We then generate the hypergraph by considering as hyperedges all the groups of simultaneous all-to-all interactions in the temporal data. Each hyperedge is stored with its weight given by the number of temporal repetitions. This is used to simulate the NL-hyper contagion process, where each hyperedge is characterized by its own infection probability.

Finally, we consider each hyperedge of order higher than 2 (involving more than 3 nodes) and we decompose it into multiple groups of three nodes (in other words we project these interactions on hyperedges of size 3). The set of the obtained three-nodes interactions, together with the hyperedges of order 2, forms the aggregated triangle network. This is used to simulate the infections taking place on triangles in the simplicial contagion model (those occurring with probability $\beta^{\Delta}$). The union of the aggregated triangle network and the pairwise network forms indeed a simplicial complex.

The unweighted networks are obtained by setting all the weights equal to 1.

\begin{table}[hb]
\centering
\begin{tabular}{|c|c|c|c|c|}
	\hline
	data set & Nb. of nodes & Nb. of pairwise & Nb. of 2nd order & Max. hyperedge \\
    &  & connections & connections
    (projected) & size\\
	\hline
	Workplace & 217 & 4274 & 772 & 4\\
	Hospital  & 75 & 1139 & 703 & 5\\
	Conference & 403 & 9565 & 2634 & 9\\
	Primary school France & 242 & 8317 & 5138 & 5\\
	High school France & 327 & 5818 & 2370 & 5\\
	University campus (48h) & 495 & 3486 & 2981 & 6\\
	Elementary school Utah (1d) & 339 & 8763 & 6868 & 8\\
	Middle school Utah (1d) & 591 & 13550 & 4453 & 7\\
	Email & 143 & 1800 & 6578 & 18\\
	\hline
\end{tabular}
\caption{\textbf{Size details about the used data sets}. The third column reports the number of pairwise connections in the aggregated network, i.e. all the connections of the networks reduced to pairwise interactions. The fourth column reports the number of 2nd order connections in the triangle network, i.e. the simplicial complex obtained by projecting higher orders on simplices of 3 nodes. The fifth column reports the number of nodes in the largest hyperedge.}
\label{tab_datasets}
\end{table}

We use the following publicly available data sets of social interactions (see also Table \ref{tab_datasets}):
\begin{itemize}
	\item \textbf{Sociopatterns}: data sets describing pairwise face-to-face interactions among individuals, collected using RFID devices \cite{SP}. These data have a temporal resolution of 20 seconds and describe contacts between individuals collected during several successive days. We consider five of these data sets: a workplace~\cite{Genois2018}, a conference~\cite{Genois2018}, a hospital~\cite{vanhems2013estimating}, a high school~\cite{mastrandrea2015contact}, and a primary school~\cite{Gemmetto2014,stehle2011high}.
	%
	\item \textbf{University campus}~\cite{sapiezynski2019interaction}: data set describing pairwise co-location interactions among  students at the campus of the Technical University of Denmark for one month, recorded by the exchange of Bluetooth radio packets between smartphones every 5 minutes for 1 month. We consider a reduced version of the original data set, involving only the first 48 hours, in order to obtain an aggregated network with a moderate density.
	%
	\item \textbf{Contacts among Utah's School-age Population}~\cite{toth2015role}: data set describing pairwise interactions among students in an elementary and a middle school in Utah. The interactions were collected with a temporal resolution of 20 seconds for 2 days via Wireless Ranging Enabled Nodes. We consider a limited version of the original data set, involving only the first day.
	%
	\item  \textbf{Email Enron}~\cite{klimt2004enron}: data set describing group interactions among a company employees via email messages: each node corresponds to an email address and a hyperedge includes the sender and all receivers of an email.
\end{itemize}

\section{Models of contagion}
\label{secSM_models}

We describe in details the four different contagion models used in the main text. For the sake of generality we define the models on weighted networks, making use of the weighted adjacency matrix $W$, which reduces to the adjacency matrix for non-weighted networks.
In order to make clear the distinction between models we give for each of them the node infection probability.

In simple contagion processes, infections take place independently on each link connecting infected to susceptible nodes at rate $\beta$ multiplied by the weight of the link (normalized such that $\max_{(i,j)}W_{ij} = 1$). The probability that a susceptible node $i$ is infected at a specific time $t$ is hence given by:
\begin{equation}
	p_i(t) = 1 - \prod_{j\in I(t)}(1-\beta W_{ij}),
	\label{eq_simple_W}
\end{equation}
where $I(t)$ is the set of infected nodes at time $t$, so the second term is the product of the independent probabilities that node $i$ is not infected by each of its infected neighbors.

\smallskip

The threshold model is a deterministic model where a susceptible node is infected when the fraction of its neighbors that are infected exceeds a  threshold $\theta$. More generally, if the network is weighted, a susceptible node is infected when the weight of its connections with infected nodes divided by the total weight of its connections exceeds the  threshold. The general probability can be described by a Heaviside function:
\begin{equation}
	p(t) = H\biggl(\frac{1}{s_i}\sum_{j\in I(t)}W_{ij}-\theta \biggr) =
	\begin{cases}
		1 \text{ if } \frac{1}{s_i}\sum_{j\in I(t)}W_{ij} \geq \theta\\
		0 \text{ if } \frac{1}{s_i}\sum_{j\in I(t)}W_{ij} < \theta
	\end{cases}\,.
	\label{eq_thresh_W}
\end{equation}
where $s_i=\sum_jW_{ij}$ is the total strength of node $i$.

\smallskip

The simplicial contagion model can be considered as a generalization of the simple model on simplicial complexes, where nodes are connected via simplices instead of only pairwise links. We only consider simplices up to the second order, i.e. three nodes interactions, neglecting structures of higher orders (which will only appear decomposed into edges of first and second order). Each susceptible node can be infected via its first-order connections, at rate $\beta$ (like in simple contagion), or via its second-order connections, if present, at rate $\beta^{\Delta}$. Second order contagions can only happen if the susceptible node is part of one (or more) closed triangles where  the other two nodes are both infected. Notice that in this case, as the model is defined on a simplicial complex, the two nodes additionally contribute independently to the first-order contagion. All infection probabilities are multiplied by the weight of the involved simplex, via the matrix $W$ at first order and tensor $W^{\Delta}$ at second order, where the element $W^{\Delta}_{ikl}$ is the weight of the closed triangle $(i,k,l)$, normalized by the same factor as the link weights. The probability per unit time that a susceptible node is infected is hence given by the combination of its first and second order interactions with infected nodes:
\begin{equation}
	p_i(t) = 1 -
 \prod_{j\in I(t)}(1-\beta W_{ij})
 \prod_{k,l\in I(t)}(1-\beta^{\Delta} W^{\Delta}_{ikl}) \ .
	\label{eq_simplicial_W}
\end{equation}

\smallskip

The non-linear hypergraph contagion model is defined on hypergraphs rather than on simplicial complexes, and  multi-body interactions are modelled through a non-linear infection term. In this case, we consider the SI epidemic model on a hypergraph with hyperedges of arbitrary size (i.e. links, triangles and also higher-order structures up to the maximum hyperedge size $M$): if in a hyperedge $h=(i_1,i_2,...i_m)$ of size $m$ there are $n_h$ infected individuals, each of the remaining $(m-n_h)$ susceptible nodes is independently infected with probability $\lambda n_h^{\nu} W^m_h$. Note that $\nu$ sets the non-linearity, $\lambda \in [0,1]$ and $m \in [2,M]$. All the infection terms are multiplied by the weight of the involved hyperedge, via the tensor $W^m$ at the $(m-1)$th order, whose element $W^m_h=W^m_{i_1,i_2,...,i_m}$ is the weight of the hyperedge $h=(i_1,i_2,...i_m)$. By assuming that each hyperedge produces an independent contagion process, the probability per unit time that a susceptible node $i$ is infected is hence given by:
\begin{equation}
	p_i(t)=1-\prod_{m=2}^M\prod_{h \in \mathcal{E}^m_i} (1-\lambda n_h^{\nu} W^m_h)
\end{equation}
where $\mathcal{E}^m_i$ is the set of all hyperedges of size $m$ containing $i$ and $n_h$ is the number of infected individuals in $h$.

\section{Classification tasks for processes on all the data sets considered}

\subsection{SI models on unweighted networks}

The analyses reported in the main text using the workplace network can be repeated with the networks described in section~\ref{SMsec_data}. The results for SI model on unweighted networks are shown in Figs.~\ref{figSM_distr_unweighted1} and \ref{figSM_distr_unweighted2}, which report the distributions of the correlations $C_1$ and $C_2$, the resulting ROC curves and the confusion matrices,
similarly to Fig.~\ref{fig2} and Fig.~\ref{fig3}(a,b) in the main text. It is evident that the results obtained in the main text can be generalized to other networks, showing the versatility of the proposed classifier (with a classification accuracy always above 0.7).
We however notice that the classification performances are worse for the networks of the two schools in France and the elementary school in Utah, with in particular
a lower ability to distinguish between simple and simplicial model when using $C_2$ (even if the AUC is still above 0.7). This could be due to the existence of a strong community structure in these networks, but future work is needed to investigate how different network topologies affect the classification based on the measure of correlations.

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_conf.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_hosp.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_highSc.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_PrSc.pdf}
\caption{\textbf{Unweighted SI model.}
Analogous of Fig.~\ref{fig2} and Fig.~\ref{fig3}(a) of the main text for four other data sets: conference, hospital, high school, and primary school. For each data set, the top panels report the distributions of $C_1$ and $C_2$. The three panels below give: in orange, the ROC curve when using $C_1$ to classify
threshold model processes against the other three; in light green, the ROC curve when using $C_2$ to classify the simplicial model against simple and threshold; in green, the ROC curve when using $C_2$ to classify the NL-hyper model against simple and threshold. The bottom panels report the confusion matrices with 4 classes (simple, complex threshold, complex simplicial and complex NL-hyper) and with 3 classes (simple, complex threshold and complex higher-order).
For the stochastic models (simple, simplicial and NL-hyper) 1000 realizations are implemented for each parameter setting, always starting with one random infected node. For the threshold model, due to its deterministic nature, we use only one realization for each different initial condition, i.e. one for each network node, for each parameter value. Confusion matrices are obtained using the same number of instances for each model.
Parameter values: $\beta \in \{0.005, 0.008, 0.014, 0.023, 0.039\}$ for both simple and simplicial models, $\beta^{\Delta} = 0.8$, $\lambda \in \{0.0001, 0.001, 0.006, 0.011, 0.015\}$, $\nu=4$,
$\theta \in \{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10\}$. Additional data sets are shown in Fig.~\ref{figSM_distr_unweighted2}.
\label{figSM_distr_unweighted1}}
\end{figure*}

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_DTU.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_Utahelem.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_Utahmid.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_email.pdf}
\caption{\textbf{Unweighted SI model.}
Same as Fig.~\ref{figSM_distr_unweighted1}
for the data sets: University campus, elementary school in Utah, middle school in Utah, and email.
\label{figSM_distr_unweighted2}}
\end{figure*}

\subsection{SI models on weighted networks}
\label{secSM_weight}

The classifier can be applied to spreading processes occurring on weighted networks: the processes spread according to the links' and hyperedges' weights,
as described in section~\ref{secSM_models}.
The classifier takes as input the correlations between the order of infection of the nodes and their strengths instead of their degrees.

The results are reported in Figs.~\ref{figSM_distr_weighted0}, \ref{figSM_distr_weighted1} and \ref{figSM_distr_weighted2}.
Overall, the classifier performs well in these cases too (accuracy always above 0.7). For some data sets, we observe however a reduced performance with respect to the unweighted case, especially for the ability to distinguish between higher-order complex contagion and simple contagion. We have checked that the weights of hyperedges of size larger than 2 are typically much smaller than the weights of the links, so that the spreading processes along links dominate even for higher order processes: as a result, higher-order complex and simple contagion become more similar processes.

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_workplace.pdf}
\caption{\textbf{Weighted SI model.} Analogous to Figs.~\ref{fig2} and \ref{fig3} for the workplace data set, using the weighted aggregated networks for the spreading processes. The parameters are chosen as in Fig.~\ref{figSM_distr_unweighted1}, except for $\theta$, which becomes a threshold for nodes strength instead of nodes degree and is consequently adapted, taking the values $\{0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22\}$.
\label{figSM_distr_weighted0}}
\end{figure*}

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_conf.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_hosp.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_highSc.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_PrSc.pdf}
\caption{\textbf{Weighted SI model.} Analogous of Fig.~\ref{figSM_distr_weighted0} for other networks: conference, hospital, high school in France, and primary school in France. Additional data sets are shown in Fig.~\ref{figSM_distr_weighted2}.
\label{figSM_distr_weighted1}}
\end{figure*}

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_DTU.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_Utahelem.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_Utahmid.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_w_email.pdf}
\caption{\textbf{Weighted SI model.} Same as Fig.~\ref{figSM_distr_weighted1} for the data sets: University campus, elementary school in Utah, middle school in Utah, and email.
\label{figSM_distr_weighted2}}
\end{figure*}

\subsection{SIR models on unweighted networks}
\label{secSM_SIR}

The above analyses can be repeated using an SIR instead of an SI model, i.e. including the possibility for nodes to recover (with a fixed probability per unit time $\mu$) after they have been infected. A recovered node cannot infect other nodes or be infected again. The process ends when none of the nodes is infected (they are all susceptible or recovered).

The results of the classifier in this case are reported for unweighted networks in Figs.~\ref{figSM_distr_SIR0}, \ref{figSM_distr_SIR1} and \ref{figSM_distr_SIR2}.
The classifier performances are slightly reduced with respect to the SI case. This can be understood by the fact that the nodes remain infected only a limited amount of time (independent of their topological properties and of their neighbors' properties). When a node recovers, it does not anymore contribute to the risk of infection of its neighbours. When reinforcement effects are needed (complex contagion), this can thus alter substantially the relationship between the structure of a node's neighborhood and its risk of becoming infected, especially if the recovery is fast.

Here, as a preliminary analysis, we have considered
simulations with $\mu=0.1$, and
the classifier performances remain good, with accuracy values always above 0.7.
A future deeper analysis will investigate in more details the impact of the recovery parameter on the classifier performance.

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_workplace.pdf}
\caption{\textbf{Unweighted SIR model.} Analogous of Figs.~\ref{fig2} and \ref{fig3} for the SIR model with $\mu=0.1$ with the workplace network. Additional data sets are shown in Figs.~\ref{figSM_distr_SIR1} and \ref{figSM_distr_SIR2}.
\label{figSM_distr_SIR0}}
\end{figure*}

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_conf.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_hosp.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_highSc.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_PrSc.pdf}
\caption{\textbf{Unweighted SIR model.} Analogous of Fig.~\ref{figSM_distr_unweighted1} for the SIR model with $\mu=0.1$. Results for additional data sets are shown in Fig.~\ref{figSM_distr_SIR2}.
\label{figSM_distr_SIR1}}
\end{figure*}

\begin{figure*}[thb]
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_DTU.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_Utahelem.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_Utahmid.pdf}
\includegraphics[width=0.41\textwidth]{figs/SM/distr+conf_SIR_email.pdf}
\caption{\textbf{Unweighted SIR model.}
Same as Fig.~\ref{figSM_distr_SIR1} for
additional data sets: University campus, elementary school in Utah, middle school in Utah, and email.
\label{figSM_distr_SIR2}}
\end{figure*}

\clearpage
\newpage

\section{Classification across data sets}
\label{sec_train1_test2}

We report here the results of classification experiments in which we train and test the classifier
on instances of the spreading processes taking place on distinct networks.

In a first experiment, we train the classifier using instances of spreading processes run on one data set, and we test on instances simulated on a second data set. The classification accuracy depends on the  networks chosen for training and testing, as reported in the table in Fig.~\ref{figSM_net_of_net}(a), where the results above 0.7 are highlighted (diagonal elements represent the cases where training and testing are implemented using the same data set).
The same results are also depicted on the right as a network of data sets: arrows start from the data set used for training and end at the data set used for testing. Only arrows associated to an accuracy above 0.6 are drawn, in purple for accuracy between 0.6 and 0.7, in turquoise above 0.7.

The second experiment instead involves training on eight different data sets and testing on the ninth one. The table in Fig.~\ref{figSM_net_of_net}(b) reports the accuracy obtained for testing on each data set. The highest accuracy is obtained when the testing data set is the workplace (we report the corresponding confusion matrix on the right of the figure).

In all these experiments we consider the same model parameter values as in Fig.~\ref{fig2}:

$\beta \in \{0.005, 0.008, 0.014, 0.023, 0.039\}$ for both simple and simplicial models,

$\beta^{\Delta} = 0.8$,

$\lambda \in \{0.0001, 0.001, 0.006, 0.011, 0.015\}$,

$\nu=4$,

$\theta \in \{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10\}$.

\begin{figure*}[thb]
\includegraphics[width=1\textwidth]{figs/SM/Train_net1_test_net2.pdf}
\caption{(a): Train on one data set, test on another one. The table reports the accuracy for each pair of data sets, while the network only reports the pairs with accuracy above 0.6, where arrows go from the training to the testing data set. (b): Train on eight data sets, test on the ninth. The table reports the accuracies and the confusion matrix is shown for the best case: testing on workplace with training implemented on the other data sets.}
\label{figSM_net_of_net}
\end{figure*}

\clearpage
\newpage

\section{Robustness with respect to incomplete observations}

The results presented in the main text have been obtained by observing the entire infection chains, i.e. all the infected nodes with their infection order. Here we explore a more realistic case where the access to infection data is limited and only a partial observation of infection chains is possible. In the two following sections we show how the proposed classification changes with respect to a reduction in the number of observed infections: first observing just a random percentage of nodes, secondly observing just the first infected.

\subsection{Observing only a fraction of the nodes}

Let us assume that only a random subset of nodes can be observed. We rely on the same numerical simulations that are presented in the main text, with the difference that the correlation between infection order and topological features is computed only for an a priori chosen subset of nodes, ignoring all the others. We repeated the experiment for 10 realizations of the random subsets, and for each
of the following subset sizes: $20\%$, $40\%$, $60\%$, $80\%$ and $100\%$ of the network nodes. For each experiment we computed the Random Forest accuracy of the classification in three classes $\{$Simple/Threshold/Higher-order$\}$, the area under the ROC curve when using $C_1$ to single out the threshold model vs. the others and the AUC when using $C_2$ to single out the higher-order models vs.
the simple and threshold ones. The results are reported in Fig.~\ref{fig_perc} for the nine data sets considered: the classification remains very efficient
in all cases.

\begin{figure*}[thb]
\includegraphics[width=0.9\textwidth]{figs/SM/cross_valid_random_h_unweighted.pdf}
\caption{Random Forest accuracy and AUC when using $C_1$ and $C_2$, as a function of the number of observed nodes in the different data sets. The $x$-axes report the number of observed nodes, corresponding to $20\%$, $40\%$, $60\%$, $80\%$ and $100\%$ of the total. The dots and bars show average and standard deviation for results collected on 10 realizations of the random subsets for each size.
\label{fig_perc}}
\end{figure*}

\clearpage
\newpage
\subsection{Observing only the first infected nodes}

As a second experiment, we try to classify contagion models by observing only the first infected nodes. The number of nodes to observe is
chosen as in the previous experiment, corresponding to  $20\%$, $40\%$, $60\%$, $80\%$ and $100\%$ of the total size.
In this case, Figure \ref{fig_first_h} shows that the performance depends more strongly on the fraction of observed nodes,
and is smaller with respect to the case where the same number of nodes have been chosen at random.
Combined with the results of Fig.~\ref{fig_perc}, this suggests that a reduction of the number of observed nodes does not have a drastic
effect if they are distributed along the entire infection chain, but the method is less robust if only the first phase of the infection pattern is observed.

\begin{figure*}[thb]
\includegraphics[width=1\textwidth]{figs/SM/cross_valid_first_h_unweighted.pdf}
\caption{Random Forest accuracy and
AUC when using $C_1$ and $C_2$, obtained by observing only the first infected nodes in each infection chain. The number of observed nodes is reported on the $x$-axes, corresponding to $20\%$, $40\%$, $60\%$, $80\%$ and $100\%$ of the total number of nodes.
}
\label{fig_first_h}
\end{figure*}

\clearpage
\newpage

\section{Influence of model parameters}
\label{secSM_param}

The classification experiments shown in the main text take into account the four different models of contagion with several choices of each model's parameters. In some cases the classifier performances depend on these parameters, while in others the parameters choice is less relevant, as discussed in the following.

\subsection{Parameters of the higher-order models}

The parameters' choice is particularly relevant when
the variation of these parameters can make two models' contagion mechanisms more similar.
For instance the simplicial model is an extension of the simple one and reduces to a simple propagation for $\beta^{\Delta} \to 0$.
The dynamics -- and hence the order at which nodes are infected -- in the simplicial model appears thus more similar to that of
simple model for smaller values of $\beta^{\Delta}$ and larger values of $\beta$: in such cases, contagion events take place more via pairwise
links than via triangles, weakening the importance of the higher-order mechanism.
This naturally affects the classification performances, as shown in Fig.~\ref{fig_paramHO}(a),
which shows the area under the ROC curve when using $C_2$ to distinguish instances of the simplicial model from instances of the simple propagation, for
varying values of $\beta$ and $\beta^{\Delta}$
(the instances of the simple model correspond to values of $\beta$ ranging between 0.005 and 0.1).
The AUC decreases when $\beta$ increases and when $\beta^{\Delta}$ decreases.

Analogously, in Fig.~\ref{fig_paramHO}(b) we show the AUC obtained when using $C_2$ to classify instances of the NL-hyper model and instances of the simple model, at varying $\nu$ and $\lambda$.
In this case, the AUC decreases slightly when $\lambda$ increases and $\nu$ decreases, but remains very high for all explored parameter values.

\begin{figure*}[thb]
\subfigure[]{\includegraphics[width=0.4\textwidth]{figs/SM/var_beta_c5_SI_simplicial_simple_all_beta_unweighted_workplace_lar.pdf}}
\subfigure[]{\includegraphics[width=0.4\textwidth]{figs/SM/var_lambda_c5_SI_NLHO_simple_all_beta_unweighted_workplace_lar.pdf}}
\caption{AUC of the ROC curve when using $C_2$
to distinguish simplicial vs. simple contagion models
(a) and NL-hyper vs. simple contagion (b), for
 varying parameters. Each colored pixel is the result of 1000 experiments. The spreading processes are simulated on the workplace data set.}
\label{fig_paramHO}
\end{figure*}

In Fig.~\ref{fig3}(d) of the main text we showed the confusion matrix obtained when classifying the four contagion models using all the nine considered networks. We depict here the confusion matrix obtained when considering the two higher-order models as a single class,
and we show the related accuracy obtained at varying the parameter $\beta$ in simple and simplicial contagion (Fig.~\ref{fig_accuracy}(a)).
Analogously we show the confusion matrix and the accuracy obtained at varying $\beta$ relative to the experiment reported in Fig.~\ref{fig3}(e), where the classifier is trained on the High School France data set and tested on the conference data set (Fig.~\ref{fig_accuracy}(b)).

\begin{figure*}[h]
\includegraphics[width=0.6\textwidth]{figs/SM/accuracy_fig_3.pdf}
\caption{Confusion matrix and classification accuracy vs. $\beta$ in simple and simplicial model in the case of classification with all data sets (a)
and in the case of the classification on the conference data set after training on the High School France network (b).
For (a) five different splits of the correlation measures in training and testing sets have been considered, the plot reports mean and standard deviation.}
\label{fig_accuracy}
\end{figure*}

\clearpage
\newpage

In Figs.~\ref{fig_params_simple}, \ref{fig_params_simplical}, and \ref{fig_params_NLhyper} the distributions of $C_1$ and $C_2$ are shown
at varying parameters $\beta$ and $\lambda$. While the correlation distributions obtained with the simple contagion model
are not particularly sensitive to the value of its unique parameter $\beta$,
the effect of $\beta$ is instead important in the simplicial model, where it plays a role in the competition between contagion
via links and via triangles. Smaller values of $\beta$ (i.e. more contagions on triangles) lead to smaller values of the
correlation $C_2$ (larger in absolute value), and thus to larger differences with respect to the simple contagion model.
A similar effect is observed in the NL-hyper model at decreasing $\lambda$.

\begin{figure*}
{\includegraphics[width=1\textwidth]{figs/SM/distr_param_simple_unweighted.pdf}}
\caption{\textbf{Simple contagion model.} Distributions of correlations $C_1$ and $C_2$ at varying parameter $\beta$ in unweighted SI simulations on the nine considered networks.}
\label{fig_params_simple}
\end{figure*}

\begin{figure*}
{\includegraphics[width=1\textwidth]{figs/SM/distr_param_simplicial_unweighted.pdf}}
\caption{\textbf{Simplicial contagion model.} Distributions of correlations $C_1$ and $C_2$ at varying parameter $\beta$ (while $\beta^{\Delta}$ is fixed at 0.8) in unweighted SI simulations on the nine considered networks.}
\label{fig_params_simplical}
\end{figure*}

\begin{figure*}
{\includegraphics[width=1\textwidth]{figs/SM/distr_param_NLhyper_unweighted.pdf}}
\caption{\textbf{NL-hyper model.} Distributions of correlations $C_1$ and $C_2$ at varying parameter $\lambda$ (while $\nu$ is fixed at 4) in unweighted SI simulations on the nine considered networks.}
\label{fig_params_NLhyper}
\end{figure*}

\clearpage
\newpage

\subsection{Threshold parameter in threshold model}

The correlations measured in the threshold model are almost always around zero, since the nodes susceptibility
usually does not depend on their local neighborhood. This is hence almost independent on the threshold parameter $\theta$,
as shown in Fig.~\ref{fig_params_thresh}. We however notice that for certain networks smaller values of $\theta$ correspond to
smaller values of correlations $C_1$ and $C_2$. This may be due to the fact that the smaller is $\theta$ the less important is the
threshold effect, and the more similar the threshold model becomes to the simple contagion one.
The effect is however small so overall $\theta$ does not affect significantly the classifier performance.

\begin{figure*}[th]
{\includegraphics[width=1\textwidth]{figs/SM/distr_param_thresh_unweighted.pdf}}
\caption{\textbf{Threshold model.} Distributions of correlations $C_1$ and $C_2$ at varying parameter $\theta$ in unweighted SI simulations on the nine considered networks.}
\label{fig_params_thresh}
\end{figure*}


\section{\label{sec:level1}First-level heading}

In wake of tremendous scientific and technological advancements in the last decade,  photonic time crystals (PTCs) have been introduced as a distinct tool for light manipulation with unique features opening up new avenues in the domain of photonics.
PTCs are systems where the refractive index ($n$) varies periodically with time \cite{shaltout2016photonic, zeng2017photonic}, the study of which has brought light to many underlying features like temporal localization in between two PTCs with dissimilar topologies \cite{lustig2018topological}, localization of pulses and exponential decay in group velocity of a pulse in disordered PTCs \cite{sharabi2021disordered}, exceptional points in spatio-temporal crystals \cite{sharabi2022spatiotemporal}, and more with their intriguing physics.
The concept of PTCs, so far, has been demonstrated in transmission lines for radio frequencies \cite{reyes2015observation} and tremendous efforts are being put forward towards their realisation in optical regime using transparent conducting oxides (TCOs) \cite{lustig2021towards}.
Moreover, in PTCs, the periodic temporal variation of $n$ yields an exclusive dispersion diagram which is gaped in momentum ($k$) \cite{shaltout2016photonic,lustig2018topological, ma2019band} leading to a range of unconventional signatures like pulse amplification, pulse broadening etc. The amplification as of the currently available literature is an integral part of the gap propagation (when $k$ of the propagating pulse lies within the gap) in PTCs.
\\
The pulse amplification in PTCs is a consequence of broken continuous time translation symmetry \cite{zeng2017photonic}.
Such kind of pulse amplification in PTCs is fundamentally different from optical parametric amplification (OPA) \cite{sharabi2022spatiotemporal}.
When a pulse enters a PTC along with exponential growth in intensity, the group velocity of the propagating pulse decreases simultaneously and eventually becomes zero \cite{lustig2018topological, sharabi2021disordered}. Moreover, the pulse suffers from spatial and temporal broadening due to dispersion in the medium. An amplification of order four of magnitude is accompanied by a broadening of 2.5 times can be observed under such a scenario \cite{lustig2018topological}.
When the PTC is over, the pulse emerges as two pulses propagating in opposite directions. \\
Time reflection and refraction of the propagating pulse in the optical regime require a sudden and abrupt periodic change in $n$ that leads to
formulation of significant $k$-gaps.
While achieving a gap propagation, the strong interference among these periodically reflected and transmitted pulses give rise to eigenmodes with imaginary Floquet frequencies which result in an amplification of the pulse \cite{sharabi2021disordered}. \\
The amplification in PTCs has huge potential to fuel state-of-the-art lasing technologies and other active photonic devices.
To exercise control over such amplification it is important to characterize the amplification and its dependence.
However, the present literature lags behind the theoretical treatment that deals with pulse amplification within the k-gap of PTCs.

In this work, we report for the first time, the $k$ dependent amplification of a propagating pulse within the $k$-gap of a linear binary PTC system. We investigate the pulse dynamics in side and out side of the bandgap and study the amplification in time-reflected and transmitted pulses. Further, we quantify the amount of amplification using our proposed analytical modelling and show that the growth in the intensity of the time-reflected and transmitted waves for a pulse entering a PTC having $k$ within the BG are of the same order.
Beyond demonstrating  analytically, we establish the findings using finite difference time domain (FDTD) simulations that the amplification achieved are different at various ranges of $k$-gaps. Additionally, the study reveals that the amplification is non-uniform throughout a $k$-gap and attains a maximum at the centre of the gap.
\begin{figure}[b!]
		\includegraphics[width=8.5cm]{Figure1.eps}
		\caption{(a) Schematic of a binary photonic time crystal where permittivity varies periodically in a step-like manner with time. For $0<t<t_2$ seconds $n (t) =n_2$, and $t_2<t<T-t_2$ seconds $n (t) =n_1$ with a time period of $T$ (= $t_1+t_2$). (b) During time reflection at temporal discontinuity, reflected and transmitted waves in the same medium due to causality.}
		\label{fig1}
\end{figure}

To unveil the underlining features of the amplification in a PTC, we first analytically calculate the band diagram of a PTC using the transfer matrix method that results in ranges of $k$-gaps where the PTC manifests profound amplification.
We begin by considering a plane wave propagating in a spatially homogeneous linear binary PTC [Fig. \ref{fig1}(a)] and describe the dynamics using Maxwell's equations for the electric displacement field $D(z,t)$ coupled to the Floquet-Bloch theorem.
The plane-polarised wave propagating in $z$ direction through such time-dynamic media undergoes time reflections at the temporal interfaces, and after the reflection, both time-reflected and time-refracted waves coexist within the same medium [Fig. \ref{fig1}(b)].
These time-reflected and transmitted waves interfere to produce the Floquet-Bloch modes in PTCs.
For simplicity, we consider a non-magnetic media.
To obtain significant $k$-gaps to influence the plane wave propagation through PTC, a refractive index difference of the order of unity is required and we choose the permittivity ($\epsilon$) in general  varying from $\epsilon_2$ ($=1)$ (for $t_2$ fs) to $\epsilon_1$ ($=3$) (for $t_1$ fs) in a periodic-step like fashion in time that results in a binary PTC [see Fig. \ref{fig1} (a)].
Since in a practical photonic system the time-reflection vanishes for changes slower than the $4.3$ fs time period \cite{lustig2021towards}, the binary variation is considered to observe the time reflection in  the analytical model.
Thus we choose a time period ($T$) of 2 fs and considered $t_1=t_2=1$ fs.
During the time reflections at temporal interfaces, the momentum is conserved but temporal frequency gets modified and the instantaneous $D_x(z,t)$ within the PTC is written as a sum of forward propagating and time-reflected waves as
\begin{widetext}
\begin{equation}
D_x(z,t)=
\begin{cases}
      (a_p e^{i\omega_1 (t-pT)} + b_p e^{-i\omega_1 (t-pT)} )e^{-ik_zz}, & pT-t_1<t<pT\\
      (c_p e^{i\omega_2 (t-pT+t_1)} + d_p e^{-i\omega_2 (t-pT+t_1)} )e^{-ik_zz}, & (p-1)T<t<pT-t_1 \\
    \end{cases}
    \end{equation}
\end{widetext}
where, $p$ is the $p^{th}$ unit cell in time and the frequency of the pulse getting modified as $\omega_{1,2}=kc/n_{1,2}$ ($c$= velocity of light in free space) keeping the $k$ constant at the interfaces, and $a_p$, $b_p$, $c_p$ and $d_p$ are constants which are related by the continuity conditions at the temporal interfaces involving $D_x$ and $B_y$.
Utilizing the continuity of $D_x$ and $B_y$ at $t= pT-t_1$ and $t= pT$ interfaces we obtain the matrix equation

\begin{figure}[b!]
		\includegraphics[width=8.5cm]{Figure2.eps}
		\caption{(a) Imaginary part of $\Omega$ as a function of normalized $k$ (where, $k_0=2\pi /Tc$). (b) Band diagram of a PTC with infinite number of bilayers in time, showing allowable bands (in blue) representing ranges of $k$ where $\Omega$ is real and bandgaps (BGs) (in green) where $\Omega$ is imaginary.}
		\label{fig2}
\end{figure}
\begin{gather}
 \begin{bmatrix} c_{p}  \\ d_{p}  \end{bmatrix}
 =  \begin{bmatrix}
   P  &
   Q \\
   R &
   S
   \end{bmatrix}  \begin{bmatrix} c_{p+1}  \\ d_{p+1}  \end{bmatrix}
   \label{Transfer_matrix}
\end{gather}
where,
\begin{subequations}
\begin{equation}
  P= e^{i\omega_2 t_2} (\cos(\omega_1 t_1) + \frac{i}{2} (\frac{\omega_1}{\omega_2}+\frac{\omega_2}{\omega_1}) \sin(\omega_1 t_1)
\end{equation}
\begin{equation}
  Q=  \frac{i}{2} e^{i\omega_2 t_2} (\frac{\omega_1}{\omega_2}+\frac{\omega_2}{\omega_1}) \sin(\omega_1 t_1)
\end{equation}
\end{subequations}
and $S=P^*$, $R=Q^*$.
To get the dispersion relation using this matrix relation we recall the $D_x(z,t)$ for the time-periodic system exhibiting $n(t)=n(t+T)$ follow the Floquet-Bloch theorem and modifies as $D_{\Omega}(t)e^{-i\Omega t} e^{ikz}$. Thus we write
\begin{equation}
    \Omega(k_z)= \frac{1}{T} \cos^{-1} [\frac{P+S}{2} ].
\end{equation}
By plotting $\Omega$ as a function of $k_z$ the band diagram for the given PTC [Fig. \ref{fig2}] with infinite number of temporal unit cells is obtained and the same is used to
recognize the BG ranges in $k$.
We obtain bands where $\Omega$ is real and bandgaps for $\Omega$ being complex [see Fig. \ref{fig2}(a)].
The first few (five in this case) of  such BGs ranges in $k$ are as follows $0.514$ to $0.731$ (which we term as BG1), $1.191$ to $1.353$ (BG2), $1.867$ to $1.933$ (BG3), $2.426$ to $2.644$ (BG4), $3.118$ to $3.229$ (BG5).
Moreover, we note the corresponding width of the BGs $0.217$, $0.162$, $0.066$, $0.218$, $0.111$ respectively (in terms of $k/k_0$).
From Fig. \ref{fig2} we observe that the wider BGs produce higher values of Im($\Omega$) within the BG.
Moreover, the PTC amplifies  plane waves with $k$ lying within these BGs.

To establish and explore the fundamental characteristics of wave propagation through PTC we simulate the pulse propagation and analyze the effect of the $k$ gaps.
We simulate the PTC by solving Maxwell's equation using FDTD method by judiciously implementing Yee's grid method and following the stability criteria for PTCs \cite{zeng2017photonic, yee1966numerical}.
We consider a gaussian pulse with full width at half maxima (FWHM) of $50$ fs.
The PTC starts at $200$ fs, and ends at $270$ fs.
At first, the pulse propagates through the vacuum and it enters the system at $200$ fs.
We simulate pulse propagation for two particular cases one for the pulse with $k$ lying within the band $k= 0.4k_0$ (which is the band propagation) [Fig. \ref{fig3}(a)] and other one for the $k$ lying within the bandgap $k= 0.62k_0$ (gap propagation) [Fig. \ref{fig3}(b)].
As we can see from Fig. \ref{fig3}(a) for the first case the propagating pulse does not get amplified and we obtain two transmitted and two time-reflected pulses at the end of the PTC. Moreover, the maximum amount of energy is transferred to the transmitted pulse leaving a relatively small amount of energy into the time-reflected pulses.
Whereas, for the gap propagation an exponential growth in the pulse intensity of $10^{14}$ is observed [Fig. \ref{fig3}(b)] and at the end of the PTC we are left with one transmitted and one time-reflected pulse, where both get equally amplified.
Here, being equipped with a $k$  not allowed to propagate through the PTC, pulse localization in space with spatial and temporal pulse broadening along with the amplification are also observed as a characteristic feature of gap propagation.

In view to predict and quantify the amount of amplification theoretically and
have a better insight into the amplification phenomenon in PTC, the transmission and reflection of a plane wave are studied.
In a practical senario, the number of temporal unit cells ($N$) is always limited by the system definition and requires to be of a finite number.
As a result, for a plane wave propagating through $N$ number of temporal unit cells, the $D$ consists of two components one forward propagating ($a_p$ or $c_p$) and time reflected ($b_p$ or $d_p$) at the $p^{th}$ unit cell depending on the portion of temporal unit cell where the wave is in.
\begin{figure}[t!]
		\includegraphics[width=8.5cm]{Figure3.eps}
		\caption{Pulse intensity profile for propagation through a PTC system depicting
(a) pulse splitting at the beginning and end of PTC due to time reflections with no amplification in case of the band propagation, (b) pulse amplification in addition to localization (in space) and pulse broadening for gap propagation. Colour bars in log scale.
}
		\label{fig3}
\end{figure}
It is evident from the PTC simulation [see Fig. \ref{fig3}] that before entering the PTC irrespective of its $k$,  the $d$ component (time-reflected part) of the field just before entering the system ($d_0$) is zero for a PTC.
In this scenario, the effect of $N$ temporal unit cells can be taken into account by simply multiplying
Eq. \ref{Transfer_matrix} $N$ times once for each  cell in time. Thus, we obtain
\begin{gather}
 \begin{bmatrix} c_{0}  \\ 0  \end{bmatrix}
 =  \begin{bmatrix}
   P  &
   Q \\
   R &
   S
   \end{bmatrix}^N   \begin{bmatrix} c_{N}  \\ d_{N}  \end{bmatrix}
   =\begin{bmatrix}
   A  &
   B \\
   C &
   D
   \end{bmatrix}   \begin{bmatrix} c_{N}  \\ d_{N}  \end{bmatrix}
   \label{TM_2}
\end{gather}
Eq. \ref{TM_2} relates the coefficients of $D$ of the wave before and after the PTC completing $N$ temporal periods.
Thus we get
\begin{subequations}
\begin{align}
 c_0 &= A c_N + B  d_N  \\
0 & = C c_N +D d_N
\end{align}
\label{subeq1}
\end{subequations}
Here we define transmission coefficient for the displacement field amplitude for PTC as
\begin{equation}
    t_N=(\frac{c_N}{c_0})_{d_0 =0}
    \label{t_n}
\end{equation}
Since the time-reflected pulse is observed at the end of the PTC, the reflection coefficient for the displacement field amplitude turns out to be
\begin{equation}
    r_N=(\frac{d_N}{c_0})_{d_0 =0}
    \label{r_n}
\end{equation}
Using Eq. \ref{TM_2}- \ref{r_n}, the amplitude transmission and reflection coefficients can be rewritten as
\begin{subequations}
\begin{align}
  t_N &=\frac{1}{A-\frac{CB}{D}}  \\
r_N &=\frac{1}{B-\frac{AD}{C}}
\end{align}
\label{t_n2}
\end{subequations}
The transmitivity and reflectivity of the intensity i.e. the transmittance and reflectance is calculated by calculating the absolute value from Eq. \ref{t_n2} as
\begin{subequations}
\begin{align}
|t_N|^2 &=\frac{1}{|A-\frac{CB}{D}|^2}  \\
|r_N|^2 &=\frac{1}{|B-\frac{AD}{C}|^2}
\end{align}
\label{abst_n2}
\end{subequations}
We plot the $|t_N|^2$ and $|r_N|^2$ as a function of momentum $k$ in Fig. \ref{fig4} (a) and (b) respectively and obtain the output intensity for both the transmitted and time reflected waves in PTC, analytically.
\begin{figure}[t]
		\includegraphics[width=8.5cm]{Figure4.eps}
		\caption{ (a) Transmittance and (b) reflectance plot of PTC showing amplification of the forward propagating and the time-reversed pulses within the BGs. The amplification is maximum at the middle of the BG and is different at different ranges of BG.
}
		\label{fig4}
\end{figure}
In Fig. \ref{fig4} (a), we observe the first band (i.e. $0<k/k_0<0.514  $) and it is evident that the transmission coefficient remains almost one for $k$ values ranging within the band which indicates that for plane waves with $k$ within the band there is no amplification which is the case of band propagation and the same scenario is true for all other bands.
A similar analysis of the reflection plot [Fig. \ref{fig4}(b)] for the bands  clearly shows that only a certain fraction of the total intensity is reflected back within the range.
As we move to the BG1 (in case of the gap propagation) the value of $t_N$ starts rising, indicating the pulse amplification within the $k$ gap.
The $t_N$ attains a maximum value at the center of the BG1.
Moreover, as we start moving away from the BG centre, the $t_N$ value starts falling and again it reaches to unity at the end of the BG1.
The same observation is true for all the other BG ranges.
Thus, for a PTC, amplification within a single $k$ gap is not uniform.
The amplification takes a certain value depending on the choice of value of the momentum and it is the maximum for a $k$ value at the centre of the BG.
It is evident from the reflectance plot in Fig. \ref{fig4}(b) for the gap propagation that the time-reflected wave generated at the end of the PTC is also equally amplified along with the transmitted wave and show similar characteristics of attending a maxima at the gap center.
Moreover, both the plots Fig. \ref{fig4} (a) and (b) show that the amplification obtained (in transmitted and time-reflected wave) at different the BG ranges is not the same and it is different for different ranges of BGs.
Additionally, we observe different BG centers are producing different amount of amplifications.
As from our analytical results, the estimated amplification at the BG centre for BGs are $\sim 10^{14}$ (BG1), $\sim 10^{11}$ (BG2), $\sim 10^{4}$ (BG3), $\sim 10^{16}$ (BG4), $\sim 10^{8}$ (BG5) respectively.

\begin{figure}[b]
		\includegraphics[width=7.5cm]{Figure5.eps}
		\caption{Comparison of the amplification of the (a, c) transmitted and (b, d) time-reflected pulse within the BG1 (a, b) and BG2 (c, d) (red dash: simulated data, blue line: analytical results).
}
		\label{fig5}
\end{figure}

To verify these properties of PTCs and validate our theoretical findings, we simulated a gaussian pulse propagation through a PTC.
This time we keep all the parameters same as in the previous case [Fig. \ref{fig3}] but select the $k$ of the pulse.
We observe the amplification of the transmitted and time-reflected pulse for BG1 and BG2 and plot them against $k$.
Observing the amplification from simulation data lead us to Fig. \ref{fig5}, where the amplification of the time reflected and transmitted pulse has been plotted and compared with the same derived from the analytical model.
For the BG1, the transmitted and time-reflected pluses are amplified in the order of $\sim 10^{14} $ in both cases [See Fig. \ref{fig5}]. Whereas, for the BG2, both pluses are amplified in the order of $\sim 10^{11} $ [See Fig. \ref{fig5}].
Thus, the simulated results are in complete agreement with the proposed analytical model and both of them show a growth in intensity of the same order.
Moreover, Fig. \ref{fig5} also verifies the fact that the amplification attains a maxima at the center of the BGs.
A careful observation of the Fig. \ref{fig5} reveals that the simulated curves are a bit widen than the analytical one and there exists a slight difference in amplifications between them.
Such difference attributes to two different schemes.
In case of analytical calculations, plane wave has been considered whereas in simulations pulse propagation is utilized.
With judiciously chosen grid resolutions and pulse parameters, the difference can be minimised to a value as low as $<5\%$.
Further to verify that the different BGs produce different amount of amplifications, we simulate the pulse propagation with $k$ values at the gap centers i.e. $0.63$ (BG1), $1.27$ (BG2), $1.91$ (BG3), $2.55$ (BG4), $3.17$ (BG5).
We compare the simulated results with the analytical one [in Fig. \ref{fig6}].
As from our simulated results, the amplification at the BG centre for BGs are $\sim 10^{14}$ (BG1), $\sim 10^{11}$ (BG2), $\sim 10^{4}$ (BG3), $\sim 10^{17}$ (BG4), $\sim 10^{9}$ (BG5), respectively which are in agreement with the analytically calculated values and verifies our observations from analytical calculations.
\begin{figure}[t]
		\includegraphics[width=7.5cm]{Figure6.eps}
		\caption{Amplification of the optical pulse at several BG centers (in red) and comparison with simulated data (in green).
}
		\label{fig6}
\end{figure}

Amplification in the gap propagation in PTCs is observed due to their periodic variation in time that break continuous time translation symmetry,
as a result does not conserve energy.
In gap propagation, the $\Omega$ becomes imaginary and as a result the forward propagating and time-reflected waves generated periodically within the PTC interfere to generate an envelope function of the form $\sim \sinh \Omega t + \cosh \Omega t$.
By comparing Figs. \ref{fig2}(a) and (b), we note that wider bandgaps result in larger values of Im($\Omega$) within the BG as well as at the gap center and larger values of Im($\Omega$) generate higher amplification within the BGs.
The BG in Fig. \ref{fig2} has been derived by considering infinite number of temporal unit cells ($N$).
By limiting the $N$ for practical applications, we limit the amount of amplification with increasing number of $N$ the amplification asymptotically approaches the theoretically predicted value in as case of infinite $N$.
In PTCs the $n$ varies periodically in time and this variation is generated by a external agency which is the driver system.
Now $\Omega$ is a function of the driving period ($T$) which controls the BGs.
Thus, the extra energy is contributed by the driver system.

In conclusion, we have shown that the amplification of the propagating optical pulse depends on the momentum within the BG of PTC, and the amplification is maximum at the center of a specific BG.
Moreover, the analytical modeling successfully quantifies the amount of amplification of the transmitted and time-reflected pulses propagating in opposite directions and establish the fact that different BGs produce different amount of amplifications.
Such new observations, involving engineering of the optical properties in PTCs, introduce a new degree of freedom to control the light and may open up a new avenue for future unconventional lasing and other active device applications in photonics.


\section{Introduction}

At high collision energy or small $x$, Deep Inelastic Scattering can be conveniently described in the dipole picture, combined with an eikonal scattering approximation for the scattering of the dipole off the target.
When gluon saturation is important, the number of gluons given by a gluon distribution is not the best way to quantify the strength of the gluonic field of a hadronic target.
In the high collision energy limit, the interaction of a high-energy probe  should be eikonal, meaning that the transverse coordinate of the probe does not change during the scattering. The basic degrees of freedom in our description are thus eikonal scattering amplitudes.

The amplitude for the simplest colored probe, a single quark, has a microscopical interpretation in terms of a path-ordered exponential, Wilson line $V({\mathbf{x}})$.
The Wilson line also provides a direct connection to the classical field in the initial stage of a heavy ion collision.
In the case of DIS, the relevant dilute probe is a color-neutral quark-antiquark dipole. Its scattering amplitude is given by the dipole amplitude
\begin{equation}
\mathcal{N}_{01} = 1-S_{01}= 1-\left<
 \frac{1}{N_\text{c}}
\mathrm{Tr} V({\mathbf{x}}_0)
V^\dag({\mathbf{x}}_1)\right>,
\end{equation}
which automatically interpolates between color transparency at $|{\mathbf{x}}_{01}| \equiv |{\mathbf{x}}_0-{\mathbf{x}}_1|=0$ and saturation at $|{\mathbf{x}}_{01}|\gtrsim 1/Q_{\text{s}}$.

This picture leads directly to the \emph{dipole picture of DIS}. Here, at leading order, the DIS scattering is factorized into the $\gamma^*$ fluctuating into a $q\bar{q}$ dipole, described by the photon light-cone wave function,  and the dipole amplitude. The total cross section is given, through the optical theorem, with the dipole then transforming back to a virtual photon. Diffractive DIS, at the focus of our attention  in Ref.~\cite{Beuf:2022kyp}, requires a specific color neutral final state after the target color field. At leading order the final state is another $q\bar{q}$ dipole, while  at NLO it includes an additional $q\bar{q}g$ state.

\section{Diffractive structure function at leading order}

The diffractive structure function $ F_2^{D(3)}(\beta, Q^2,{x_{\mathbb{P}}})$ (or equivalently the diffractive $\gamma^*$-target cross  section     $\, \textrm{d} \sigma^{\text{D}}_{\lambda,  \, q \Bar q}/\, \textrm{d} M_X^{2} \, \textrm{d} \abs{t}$, see \cite{Marquet:2007nf} for more details) is   a function of three Lorentz-invariant kinematical variables. In addition to the conventional $Q^2$ and ${x_\text{Bj}}=\beta{x_{\mathbb{P}}}$, it also depends on how the  $\gamma^*$-target energy is divided between the diffractive system $X$ and the rapidity gap, with $\beta = Q^2/(Q^2+M_X^2),$ where $M_X$ is the invariant mass of the diffractive system. In this paper we are interested in the regime where $\beta$ is not parametrically small (i.e. $M_X$ not parametrically large), so that powers of $\alpha_{\text{s}} \ln 1/\beta$ do not need to be resummed. In different regimes of $0<\beta<1$, the cross section is dominated by different kinds of partonic configurations of the virtual photon~\cite{Golec-Biernat:1999qor}. At $\beta\to 1$, one predominantly has a $q\bar{q}$ state in a longitudinal total helicity configuration, and at $\beta\approx 1/2$ in a transverse one. Although these LO $q\bar{q}$ contributions have been known already, in  Ref.~\cite{Beuf:2022kyp} we derived new expressions allowing for a completely general impact parameter dependence of the cross section:
\begin{multline}
    \frac{\, \textrm{d} \sigma^{\text{D}}_{\lambda,  \, q \Bar q}}{\, \textrm{d} M_X^{2} \, \textrm{d} \abs{t}}
    =
    \frac{N_\text{c}}{4\pi} \int_{0}^1 \!\! \, \textrm{d} z
 \!\!\! \!\!\!\!\!\!     \int\limits_{ {\mathbf{x}}_0 {\mathbf{x}}_1  \Bar {\mathbf{x}}_0  \Bar{\mathbf{x}}_1}
 \!\!\!\!\!\!   {\cal I}_{{\boldsymbol{\Delta}}}^{(2)} {\cal I}_{M_X}^{(2)}
   \widetilde{\psi}_{\gamma^{*}_\lambda \rightarrow q_{\Bar 0} \Bar{q}_{\Bar 1}}^\dagger
    \widetilde{\psi}_{\gamma^{*}_\lambda \rightarrow q_0 \Bar{q}_1}
    \big[S_{\conj{0} \conj{1}}^\dagger -1 \big]\big[S_{01} -1 \big] ,
\label{eq:loresult}
\end{multline}
where $z$ is a longitudinal momentum fraction.
The result is expressed in terms of ``transfer functions'' ${\cal I}_{\boldsymbol{\Delta}}^{(2)}$,
${\cal I}_{M_X}^{(2)}$, relating the $(q,\bar{q})$ coordinates  in the amplitude (${\mathbf{x}}_0,{\mathbf{x}}_1$) and conjugate (${\mathbf{x}}_{\conj{0}}, {\mathbf{x}}_{\conj{1}}$) to $t$ and $M_X$ as
\begin{eqnarray}
    {\cal I}_{\boldsymbol{\Delta}}^{(2)}
    &=&
=
\frac{1}{4\pi} J_0\left( \sqrt{\abs{t}}\: \norm{
z {\mathbf{x}}_{\Bar 0 0} - (1-z){\mathbf{x}}_{\Bar 1 1}
} \right),
\\
 {\cal I}_{M_X}^{(2)}
   & =&
=
\frac{1}{4\pi} J_0\left( \sqrt{z (1-z)}M_X \norm{\Bar {\mathbf{r}}-{\mathbf{r}}} \right).
\end{eqnarray}
At $\beta \ll 1$, one starts to become sensitive to higher invariant mass states in the photon, meaning Fock states with a higher number of partons. The first one of these is the $q\bar{q}g$ state. While the $q\bar{q}g$ Fock component is in general a NLO correction, in the small $\beta$ regime it is in fact the leading contribution. Thus contributions where a gluon is emitted before the target, and survives until the final state, are a special class of NLO corrections that are meaningful to consider separately.

\section{Contribution of the $q\bar{q}g$ Fock state}

\begin{figure}
\centerline{\includegraphics[width=0.2\textwidth]{ddis_gwavef1}
\includegraphics[width=0.2\textwidth]{ddis_gwavef2}
\includegraphics[width=0.2\textwidth]{ddis_gwavefinst1}
\includegraphics[width=0.2\textwidth]{ddis_gwavefinst2}
}
\caption{Contributions to the amplitude where a gluon is emitted before the target, and produced and measured in the final state.}
\label{fig:nloraddiags}
\end{figure}

In Ref.~\cite{Beuf:2022kyp} (see also \cite{Hanninen:2021byo}), we calculated the full contribution of the diagrams shown in Fig.~\ref{fig:nloraddiags} to the diffractive structure function. The calculation is performed in what we refer to as ``exact eikonal kinematics,''
where the only kinematical approximation is the eikonal interaction with the target. The results are expressed, analoguously to Eq.~\nr{eq:loresult}, in terms of ``transfer functions'' ${\cal I}_{M_X}^{(3)}$ and ${\cal I}_{{\boldsymbol{\Delta}}}^{(3)}$.
They are then multiplied by the squared gluon emission wavefunction:
\begin{multline}
{x_{\mathbb{P}}} F_{T,L, \, q \Bar q  g}^{\textrm{D}(4) \, \textrm{NLO}} ({x_\text{Bj}}, Q^2, \beta, t)
=
\!\!	\int_0^1 \!\! \frac{\, \textrm{d} z_0}{z_0}
	\frac{\, \textrm{d} z_1}{z_1}
	\frac{\, \textrm{d} z_2}{z_2}
	\delta(z_0 \! + \! z_1 \! + \! z_2 \! - \! 1)
\int_{{\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2,{\conj{\xt}}_0,{\conj{\xt}}_1,{\conj{\xt}}_2}
\\
\times		{\cal I}_{M_X}^{(3)} {\cal I}_{{\boldsymbol{\Delta}}}^{(3)}\
  \widetilde{\psi}_{\gamma^{*}_\lambda \rightarrow q_{\Bar 0} \Bar q_{\Bar 1} g_{\Bar 2}}^\dagger
    \widetilde{\psi}_{\gamma^{*}_\lambda \rightarrow q_0 \Bar{q}_1 g_2}
	\left[ 1 - S_{\conj{0} \conj{1} \conj{2}}^{(3)\dagger} \right] \left[ 1 - S_{012}^{(3)} \right] ,
\end{multline}
where the final factors represent the ``tripole'' Wilson line operators for the $q\bar{q}g$ states  to interact with the target~\cite{Beuf:2017bpd}.

The squared wavefunctions are lengthy, but straightforward expressions in terms of the transverse coordinates and momentum fractions; they can be found in Ref.~\cite{Beuf:2022kyp}. It is interesting to discuss the different scales appearing in the expressions. The photon wavefunctions involve, just as in the leading order case, Bessel $K$ functions that exponentially cuf off configurations that are too large compared to the photon virtuality $Q$. The argument of the Bessel function is the square root of
\begin{equation}
			Q^2X_{012}^2
			 =
			Q^2\left[z_{0} z_{1} {\mathbf{x}}_{01}^2 + z_{0} z_{2} {\mathbf{x}}_{02}^2 + z_{1} z_{2} {\mathbf{x}}_{12}^2\right],
\end{equation}
interpreted as the ratio of the formation time of the $q\bar{q}g$ system with coordinates ${\mathbf{x}}_{0},{\mathbf{x}}_{1},{\mathbf{x}}_{2}$ (with $ {\mathbf{x}}_{01} \equiv{\mathbf{x}}_{0}-{\mathbf{x}}_{1}$ etc.) to the lifetime of the virtual photon.  The 3-particle transfer function to $M_X$ is
\begin{equation}
  {\cal I}_{M_X}^{(3)}
        =
        2 \frac{z_0 z_1 z_2}{(4 \pi)^2} \frac{M_X}{Y_{012}} \ \mathrm{J}_1 \!\left( M_X Y_{012} \right),
\end{equation}
with  the coordinate combination conjugate to $M_X$
\begin{equation}
    {\mathbf{Y}}_{012}^2 =
        z_0 z_1 \left( {\mathbf{x}}_{\Bar 0 0} - {\mathbf{x}}_{\Bar 1 1} \right)^2 +
        z_1 z_2 \left( {\mathbf{x}}_{\Bar 2 2} - {\mathbf{x}}_{\Bar 1 1} \right)^2 +
        z_0 z_2 \left( {\mathbf{x}}_{\Bar 2 2} - {\mathbf{x}}_{\Bar 0 0} \right)^2.
\end{equation}
Finally, $t$ is conjugate to the  3-particle center-of-mass coordinate
\begin{equation}
    {\cal I}_{\boldsymbol{\Delta}}^{(3)} = \frac{1}{4 \pi} \  \mathrm{J}_0 \! \left( \sqrt{-t} \norm{z_0 {\mathbf{x}}_{\Bar 0 0} + z_1 {\mathbf{x}}_{\Bar 1 1} + z_2 {\mathbf{x}}_{\Bar 2 2}} \right).
\end{equation}

\section{Towards a full NLO computation}

There are several additional contributions that will still need to be included for a full NLO result for the diffractive structure function. From these additional contributions it is not easy to define any further subsets that would be finite by themselves, but they will all have to be considered together in future work. Firstly, in addition to the ``emission before target'' radiative corrections in Fig.~\ref{fig:nloraddiags}, there are the corresponding gluon emissions from quarks after the target. Unlike the contributions that we have included, however, they will have a collinear divergence from the gluons being emitted at a small transverse momentum with respect to the emitting (anti)quark. These collinear divergences are cancelled by wavefunction renormalization of the outgoing quarks, i.e. by ``propagator correction'' diagrams on the outgoing quark lines. These in turn also include a UV divergence, and thus they have to be considered together with virtual contributions. There are several such virtual corrections. The most straightforward of these are the loop corrections to the $\gamma^* \to q\bar{q}$ wavefunction, which have already been calculated in Refs.~\cite{Beuf:2016wdz,Hanninen:2017ddy}. There are also contributions where the gluon is emitted before the target, and then crosses the target but is reabsorbed without being measured, analoguously to the $q\bar{q}g$ contributions in~\cite{Beuf:2017bpd,Hanninen:2017ddy}. They  involve both UV divergences that must be subtracted and canceled with the ones from before and after the target, and large logarithms of longitudinal momenta that must be factorized into BK/JIMWLK evolution of the target. Finally, there are final state interactions with gluon exchanges between the quark and antiquark after the target. These are a rather novel kind of contribution in this context and there is still some discussion about the correct way of treating them.

\section{Known limiting cases}

While our calculation of the diffractive structure function in the ``exact eikonal kinematics'' is new, we have also checked that it reduces in specific limiting cases to results that are already available in the literature. In the limit of large $M_X$ we obtain the result (see \cite{Munier:2003zb} and references therein)
\begin{equation}
    {x_{\mathbb{P}}} F_{T,q\Bar q g}^{\textrm{D},(\text{MS})}
    =
    \frac{\alpha_{\text{s}} N_\text{c} C_\text{F} Q^2}{16\pi^5\alpha_{\mathrm{em}}}
\!\!\!
\!\!\!
    \int\limits_{ {\mathbf{x}}_0 {\mathbf{x}}_1 {\mathbf{x}}_2}
\!\!\!\!\!\!\!\!\!
\int_0^1
\!\!\!
\, \textrm{d} z \frac{    |\widetilde{\psi}^{\rm LO} |^2 }{z(1-z)}
 \frac{{\mathbf{x}}_{01}^2}{{\mathbf{x}}_{02}^2{\mathbf{x}}_{12}^2}
  \bigg[ \mathcal{N}_{02} + \mathcal{N}_{12} - \mathcal{N}_{01} - \mathcal{N}_{02}\mathcal{N}_{12} \bigg]^2 ,
\end{equation}
which factorizes into the LO $\gamma^*\to q\bar{q}$ wavefunction, the BK kernel for the emission of a soft gluon, and a squared Wilson line operator. This result is straightforwardly obtained from ours by first approximating the gluon as being soft $z_2\to 0$, in which case $M_X$ becomes dominated by the light cone energy of the gluon. One then removes the constraint on $M_X$ by integrating over the momentum fraction of the gluon $z_2$. This unconstrained integration introduces a divergence, which is however cured if final state emissions are included in the same kinematical approximation.

A rather more nontrivial task is to rederive the ``Wüsthoff result'' \cite{Wusthoff:1997fz,Golec-Biernat:1999qor} in the limit of large $Q^2$:
\begin{multline}
{x_{\mathbb{P}}}  F_{T,q\Bar q g}^{\textrm{D},(\text{GBW})}  = \frac{\alpha_{\text{s}}\beta}{8\pi^4}\sum_{f}e_f^2
\int_{{\mathbf{b}}}
\int_{\beta}^{1} \!\! \, \textrm{d} z
 \int_0^{Q^2} \, \textrm{d} k^2 k^4 \ln \frac{Q^2}{k^2}
 \\  \times
 \left[\int_0^\infty \, \textrm{d} r
 r  K_2(\sqrt{z}k r) J_2(\sqrt{1-z}kr)  \mathcal{N}_\text{adj}({\mathbf{b}},{\mathbf{r}},{x_{\mathbb{P}}}) \right]^2.
\label{eq:whoff}
\end{multline}
This result is characterized by an explicit logarithmic dependence on $Q^2$, and a $g\to q\bar{q}$ splitting function that can be associated with \emph{target DGLAP evolution}. It also depends on an adjoint representation dipole operator, with the small-size octet $q\bar{q}$ pair after the gluon emission acting as an effective gluon. The wavefunction of the photon splitting into this effective gluon dipole has a rank-2 traceless structure, which results in Bessel functions $J_2$ and $K_2$ in the final result. The key point in reaching this result starting from the dipole picture is to recognize that $z$ and $\beta$ in Eq.~\nr{eq:whoff} are to be understood as target ($k^-$) momentum fractions. To arrive at this form it is convenient to identify the target momentum fraction variables using invariant masses of the $q\bar{q}$ and $q\bar{q}g$ states before and after the scattering. One then  finds the large $Q^2$ limit by looking at the aligned jet kinematical limit for the dipole picture, with $z_0\ll z_1\ll z_2$. The derivation of the Wüsthoff result in the literature has not been very clearly documented, and certainly the original approach used to derive it has been very different than the one we use here. It is encouraging that as a side result our work has resulted in an independent rederivation of this widely and successfully used expression.

\subsubsection*{Acknowledgements}
{

\footnotesize
T.L and H.M are supported by the Academy of Finland, the Centre of Excellence in Quark Matter (project 346324) and projects 338263, 346567 and 321840. G.B is supported in part by the National Science Centre (Poland) under the research grant no. 2020/38/E/ST2/00122 (SONATA BIS 10). Y.M acknowledges financial support from Xunta de Galicia (Centro singular de investigación de Galicia accreditation 2019-2022); the ``María de Maeztu'' Units of
Excellence program MDM2016-0692 and the Spanish Research State Agency under project PID2020-119632GB-I00; European Union ERDF. G.B and Y.M acknowledge financial support from MSCA RISE 823947 ``Heavy ion collisions: collectivity and precision in saturation physics'' (HIEIC). This work was also supported under the European Union’s Horizon 2020 research and innovation programme by the European Research Council (ERC, grant agreement No. ERC-2018-ADG-835105 YoctoLHC) and by the STRONG-2020 project (grant agreement No. 824093). The content of this article does not reflect the official opinion of the European Union and responsibility for the information and
views expressed therein lies entirely with the authors.

}

\bibliographystyle{h-physrev4mod2}


\section{Introduction}

In the paper \cite{RGM} we have shown that the field of a magnetic
charge is described correctly by the $t$-component of the dual
electromagnetic potential $B_\mu$, so that the semi-infinite
singularities accompanying the $\varphi$-component of the usual
potential $A_\mu$, that must be taken into account during some
mathematical calculations, can not be considered as representing
real physical characteristics of the magnetic charge. In \cite{RGM}
our consideration was restricted to exclusively the static
spherically symmetric dyonic case that ideally suited our objective
of giving simple and clear arguments in favor of our novel approach
to the description of magnetic charges without Dirac strings. In the
present paper we shall expand our analysis to the stationary axially
symmetric dyonic black holes for which the effect of rotation
introduces additional technical difficulties; however, these
difficulties will be circumvented in an elegant way, clearly
confirming the physical conclusions of the previous paper in a more
general situation.

In the next section we shall introduce the nonzero components of the
dual electromagnetic potential $B_\mu$ within the framework of the
well-known Ernst formulation of the stationary axially symmetric
problem \cite{Ern2} and define explicitly the magnetic potential
$\Phi_m$ which, similar to the electric potential $\Phi_e$
introduced long ago by Carter \cite{Car}, takes constant value on
the black-hole horizon. The advantages of the enhanced Ernst
formalism are illustrated here by the example of the dyonic
Kerr-Newman black hole \cite{Car,MGa} for which a complete set of
the corresponding potentials will be constructed. The validity of
the original Tomimatsu mass integral \cite{Tom} will be proven in
Sec.~III with the aid of the symmetrical representation of the
electromagnetic energy-momentum tensor. Discussion of the results
obtained and conclusions can be found in Sec.~IV.

\section{The enhanced Ernst formalism}

In the theory of exact solutions of the Einstein-Maxwell equations,
the Ernst formalism, developed in two papers \cite{Ern1,Ern2} in
1968, occupies an outstanding place as constituting the basis for
various solution generating techniques and different approaches to
the multipole analysis of vacuum and electrovac spacetimes. In
particular, Ernst trivialized the derivation of the Kerr \cite{Ker}
and Kerr-Newman \cite{KNe} black hole solutions that were originally
obtained by means of hardly reproducible procedures.

The main idea of Ernst's formalism is to use the Papapetrou line
element \cite{Pap},
\begin{equation} d s^2=f^{-1}[e^{2\gamma}(d\rho^2+d z^2)+\rho^2 d\varphi^2]-f(d
t-\omega d\varphi)^2, \label{Papa} \end{equation}
describing a generic stationary axisymmetric electrovac field, with
three arbitrary functions $f$, $\gamma$ and $\omega$ depending on
the cylindrical coordinates ($\rho,z$), for reducing the
corresponding set of the Einstein-Maxwell equations to a fundamental
system of two differential equations for the complex potentials ${\cal{E}}$
and $\Phi$ of the following elegant form:
\begin{eqnarray} ({\rm Re}\,{{\cal{E}}}+\Phi\bar\Phi)\Delta{{\cal{E}}}=
(\mbox{\boldmath $\nabla$}{{\cal{E}}}+2\bar\Phi\mbox{\boldmath $\nabla$}\Phi)\mbox{\boldmath $\nabla$}{{\cal{E}}}, \nonumber\\ ({\rm
Re}\,{{\cal{E}}}+\Phi\bar\Phi)\Delta{\Phi}=
(\mbox{\boldmath $\nabla$}{{\cal{E}}}+2\bar\Phi\mbox{\boldmath $\nabla$}\Phi)\mbox{\boldmath $\nabla$}{\Phi}, \label{EE} \end{eqnarray}
where $\Delta$ and $\mbox{\boldmath $\nabla$}$ are the usual three-dimentional
Laplacian and gradient operators, respectively, and a bar over a
symbol means complex conjugation.

The potentials ${\cal{E}}$ and $\Phi$ are related to the metric functions
$f$, $\omega$ and to the $\varphi$ and $t$ components of the
electromagnetic 4-potential $A_\mu=(0,0,A_\varphi,A_t)$ by the
equations
\begin{equation} {\cal{E}}=f-\Phi\bar\Phi+i\chi, \quad \Phi=-A_t+iA^{'}_\varphi,
\label{EF} \end{equation}
and by the systems of the first-order differential equations
\begin{eqnarray} &&\partial_\rho\omega=-\rho f^{-2}[\partial_z\chi+2{\rm Im}
(\bar\Phi\partial_z\Phi)], \nonumber\\
&&\partial_z\omega=\rho f^{-2}[\partial_\rho\chi+2{\rm Im}
(\bar\Phi\partial_\rho\Phi)], \label{wx} \end{eqnarray}
and
\begin{eqnarray} &&\partial_\rho A^{'}_\varphi=\rho^{-1}f(\partial_z A_\varphi
+\omega\partial_z A_t), \nonumber\\ &&\partial_z
A^{'}_\varphi=-\rho^{-1}f(\partial_\rho A_\varphi
+\omega\partial_\rho A_t), \label{Afp} \end{eqnarray}
so that the knowledge of ${\cal{E}}$ and $\Phi$ permits one to find the
functions $f$, $\omega$, $A_t$ and $A_\varphi$ from
(\ref{EF})-(\ref{Afp}), while for the determination of the remaining
metric function $\gamma$ one has to solve the system
\begin{eqnarray} \partial_{\rho}\gamma&=&\frac{1}{4}\rho
f^{-2}[(\partial_{\rho}{\cal
E}+2\bar\Phi\partial_{\rho}\Phi)(\partial_{\rho}\bar{\cal
E}+2\Phi\partial_{\rho}\bar\Phi)- (\partial_{z}{\cal
E}+2\bar\Phi\partial_{z}\Phi)(\partial_{z}\bar{\cal
E}+2\Phi\partial_{z}\bar\Phi)] \nonumber\\ &&- \rho
f^{-1}(\partial_{\rho}\Phi\partial_{\rho}\bar\Phi
-\partial_{z}\Phi\partial_{z}\bar\Phi), \nonumber\\
\partial_{z}\gamma&=&\frac{1}{2}\rho f^{-2}{\rm Re}[(\partial_{\rho}{\cal
E}+2\bar\Phi\partial_{\rho}\Phi) (\partial_{z}\bar{\cal
E}+2\Phi\partial_{z}\bar\Phi)] - 2\rho f^{-1}{\rm
Re}(\partial_{\rho}\bar\Phi\partial_{z}\Phi), \label{ge} \end{eqnarray}
the integrability condition of which are equations (\ref{EE}).

Note that the potential $A^{'}_\varphi={\rm Im}\,\Phi$ is regarded
in the Ernst formalism as an auxiliary function, the knowledge of
which makes it possible the calculation of the corresponding
magnetic component $A_\varphi$ of the 4-potential $A_\mu$. However,
in our preceding paper \cite{RGM} we have already shown that
$A_\varphi$ does not describe correctly the field of the magnetic
charge, so it seems desirable to supplement the above formalism with
the nonzero components of the dual electromagnetic 4-potential
$B_\mu=(0,0,B_\varphi,B_t)$ that are related to the components $A_t$
and $A_\varphi$ by the first-order differential equations. Indeed,
using the one-form $B=B_t  dt+B_\varphi d\varphi$, we obtain the
desired relations by means of the formula
\begin{equation} dB=\star F, \label{BF} \end{equation}
where the star denotes the Hodge dual, and $F$ is the usual
electromagnetic 2-form. Taking into account that, on the one hand,
\begin{equation} dB=d(B_\nu dx^\nu)=\partial_a B_\nu dx^a\wedge dx^\nu,
\label{lBF} \end{equation}
and, on the other hand,
\begin{eqnarray} \star F&=&\star d(A_t  dt+A_\varphi d\varphi) \nonumber\\
&=&(g^{t\beta}\partial_a A_t+g^{\varphi\beta}\partial_a A_\varphi)
g^{ab}\sqrt{-g}\,\varepsilon_{b\beta\gamma\delta}dx^\gamma\wedge
dx^\delta, \label{rBF} \end{eqnarray}
($a,b\in\{\rho,z\}$), we get from (\ref{lBF}) and (\ref{rBF}), by
first equating the coefficients at $d\rho\wedge dt$ and $dz\wedge
dt$, the system of differential equations for $B_t$ in terms of
$A_t$ and $A_\varphi$, namely,
\begin{eqnarray} &&\partial_\rho B_t=\rho^{-1}f(\partial_z A_\varphi
+\omega\partial_z A_t), \nonumber\\ &&\partial_z
B_t=-\rho^{-1}f(\partial_\rho A_\varphi +\omega\partial_\rho A_t),
\label{BtA} \end{eqnarray}
and then, by equating the coefficients at $d\rho\wedge d\varphi$ and
$dz\wedge d\varphi$, the analogous system for the determination of
$B_\varphi$:
\begin{eqnarray} &&\partial_\rho B_\varphi=\rho^{-1}f[(\rho^2
f^{-2}-\omega^2)\partial_z A_t -\omega\partial_z A_\varphi)], \nonumber\\
&&\partial_z B_\varphi=-\rho^{-1}f[(\rho^2
f^{-2}-\omega^2)\partial_\rho A_t -\omega\partial_\rho A_\varphi)].
\label{BfA} \end{eqnarray}

A simple inspection of formulas (\ref{Afp}) and (\ref{BtA}) shows
that the $t$-component of the dual potential $B_\mu$ coincides with
the auxiliary potential $A^{'}_\varphi$ of the Ernst formalism,
i.e.,
\begin{equation} B_t=A^{'}_\varphi. \label{BAp} \end{equation}
Curiously, it is also possible to identify the component $B_\varphi$
(up to a sign) as the potential $B_2$ introduced in the paper
\cite{Kin} by Kinnersley as part of various matrix potentials of his
solution generating method; in particular, it arises as the
imaginary part of Kinnersley's potential $\Phi_2$.

The knowledge of the full set of the components $A_t$, $A_\varphi$,
$B_t$ and $B_\varphi$ of the 4-potentials $A_\mu$ and $B_\mu$ allows
one to analyze the electric and magnetic fields of the dyonic black
hole solutions in a symmetrical way advocated long ago by Schwinger
\cite{Sch}. For example, as was shown by Carter \cite{Car}, the
electric potential $\Phi_e$ determined as the combination
\begin{equation} \Phi_e=-A_t-\omega^{-1}A_\varphi, \label{Fe} \end{equation}
assumes constant value on the black-hole horizon. Having introduced
explicitly the dual potential $B_\mu$ into the Ernst formalism, we
can now define ``symmetrically'' the magnetic counterpart of
$\Phi_e$ as
\begin{equation} \Phi_m=B_t+\omega^{-1}B_\varphi, \label{Fm} \end{equation}
the magnetic potential $\Phi_m$ also taking constant value on the
horizon, which may be considered an important result following from
our approach.

The above said can be well illustrated by the dyonic Kerr-Newman
black hole solution. Following \cite{MGa}, we write its defining
Ernst potentials ${\cal{E}}$ and $\Phi$ in the form
\begin{eqnarray} {\cal{E}}&=&\frac{\sigma x-m-iay}{\sigma x+m-iay}, \quad \Phi=\frac{q+ip}{\sigma
x+m-iay}, \nonumber\\ x&=&\frac{1}{2\sigma}(r_++r_-), \quad
y=\frac{1}{2\sigma}(r_+-r_-), \quad r_\pm=\sqrt{\rho^2+(z\pm\sigma)^2},
\nonumber\\ \sigma&=&\sqrt{m^2-a^2-q^2-p^2}, \label{EFKN} \end{eqnarray}
where the parameters $m$, $a$, $q$ and $p$ stand, respectively, for
the mass, angular momentum per unit mass, electric and magnetic
charges of the black hole (we restrict our consideration to the
real-valued $\sigma$ only).

The corresponding metric functions $f$, $\gamma$ and $\omega$ have
the form
\begin{eqnarray} f&=&\frac{\sigma^2(x^2-1)-a^2(1-y^2)}{(\sigma x+m)^2+a^2y^2}, \quad
e^{2\gamma}=\frac{\sigma^2(x^2-1)-a^2(1-y^2)}{\sigma^2(x^2-y^2)},
\nonumber\\ \omega&=&-\frac{a(1-y^2)[2m(\sigma
x+m)-q^2-p^2]}{\sigma^2(x^2-1)-a^2(1-y^2)}, \label{mfKN} \end{eqnarray}
while for the electric and magnetic components $A_t$ and $A_\varphi$
of the 4-potential $A_\mu$ we have the expressions
\begin{eqnarray} A_t&=&-\frac{q(\sigma x+m)-apy}{(\sigma x+m)^2+a^2y^2}, \nonumber\\
A_\varphi&=&-py+\frac{a(1-y^2)[q(\sigma x+m)-apy]}{(\sigma x+m)^2+a^2y^2},
\label{Atf} \end{eqnarray}
where the integration constant on the right-hand side of $A_\varphi$
has been chosen equal to zero, thus determining the case with two
magnetic ``strings''.

Turning now to the components of the dual 4-potential $B_\mu$, we
see that $B_t$ is obtainable as just the imaginary part of the Ernst
potential $\Phi$, while $B_\varphi$ must be found by solving the
system (\ref{BfA}). The resulting expressions are
\begin{eqnarray} B_t&=&\frac{p(\sigma x+m)+aqy}{(\sigma x+m)^2+a^2y^2}, \nonumber\\
B_\varphi&=&-qy-\frac{a(1-y^2)[p(\sigma x+m)+aqy]}{(\sigma x+m)^2+a^2y^2},
\label{Btf} \end{eqnarray}
where the choice of the integration constant in $B_\varphi$ is the
same as for $A_\varphi$ and defines a pair of electric ``Dirac
strings''.

The only plausible conclusion that can be drawn from the structure
of the components (\ref{Atf}) and (\ref{Btf}) is that the field of
the electric charge $q$ in the dyonic Kerr-Newman solution is
described by $A_t$, and the field of the magnetic charge $p$ is
determined by $B_t$, both components $A_t$ and $B_t$ being well
behaved and asymptotically flat. In turn, the components $A_\varphi$
and $B_\varphi$ possessing the string singularities do not define
the singularity structure of this dyonic black hole solution,
playing exclusively auxiliary mathematical role in some
calculations. For instance, the components $A_\varphi$ and
$B_\varphi$ are needed for the evaluation of the electric  and
magnetic potentials $\Phi_e$ and $\Phi_m$ on the horizon ($\rho=0$,
$-\sigma<z<\sigma$, or $x=1$):
\begin{eqnarray} \Phi^H_e&=&\left.-A_t-\omega^{-1}A_\varphi\right|_{x=1}=
\frac{q(m+\sigma)}{(m+\sigma)^2+a^2}, \nonumber\\
\Phi^H_m&=&\left.B_t+\omega^{-1}B_\varphi\right|_{x=1}=
\frac{p(m+\sigma)}{(m+\sigma)^2+a^2}. \label{Femh} \end{eqnarray}

After introducing the angular momentum $J=ma$, and also recalling
that $\omega$ takes constant value on the horizon, so that
\begin{equation} \omega^{-1}(x=1)\equiv\Omega^H=\frac{a}{(m+\sigma)^2+a^2},
\label{Om} \end{equation}
one can see that the above formulas verify the Smarr mass relation
\cite{Sma}
\begin{equation} m=\sigma+2J\Omega^H+q\Phi^H_e+p\Phi^H_m. \label{Sma} \end{equation}

We now turn to the discussion of the evaluation of the Komar mass
\cite{Kom} in the dyonic Kerr-Newman solution, the issue that also
addresses the question of the distribution of that mass.

\section{Validating Tomimatsu's mass integral formula}

To calculate the Komar \cite{Kom} mass $M$ of a rotating charged
black hole, Tomimatsu \cite{Tom} derived a simple formula
\begin{equation} M=-\frac{1}{8\pi}\int_H \omega\partial_z\chi\, d\varphi dz,
\label{TMI} \end{equation}
where the integral is taken over the horizon of the black hole.
Formula (\ref{TMI}) was widely used for years in application to
non-isolated black holes in the presence of other black holes or
exterior gravitational fields. In the case of the dyonic Kerr-Newman
black hole, (\ref{TMI}) assumes the form
\begin{equation} M=-\frac{1}{4}\omega^H[\chi(y=1)-\chi(y=-1)], \label{TKN} \end{equation}
where both $\omega$ and $\chi$ must be taken on the horizon ($x=1$).
It is not difficult to verify that the corresponding $M$ calculated
with the help of (\ref{TKN}) coincides with the mass parameter $m$
in (\ref{EFKN}).

However, the validity of the mass formula (\ref{TMI}) in the
presence of magnetic charge was questioned in the paper \cite{CGa}.
The authors of \cite{CGa} used during their calculations the
conventional representation of the electromagnetic energy-momentum
tensor
\begin{equation}
T{^\mu}{_\nu}=\frac{1}{4\pi}\left(F^{\mu\alpha}F_{\nu\alpha}-\frac{1}{4}\delta^\mu_\nu
F^{\alpha\beta}F_{\alpha\beta}\right), \label{TE1} \end{equation}
which led them to a specific dyonic configuration with two magnetic
Dirac strings and an additional electromagnetic term in the
integrand of (\ref{TMI}), both strings carrying portions of nonzero
mass, so that the mass parameter $m$ becomes the sum of three
different contributions -- one coming from the surface integral
evaluated on the horizon, and two others arising from the singular
``massive'' Dirac strings. Although the version of the dyonic
Kerr-Newman black hole presented in \cite{CGa} is manifestly
physically inconsistent (see \cite{GMR} for the discussion of
unphysical features of that model), the mathematical computation of
the Komar integral performed in \cite{CGa} looks correct (albeit
some misprints). The explanation of such a seemingly puzzling
situation is quite simple in the framework of the ideas developed in
\cite{RGM} and in the present paper: the pathologies of a specific
representation of the electromagnetic energy-momentum tensor
formally taking part in the calculation of the Komar mass integral
should not be ascribed to the dyonic model itself since the
singularity structure of the magnetic charge is determined by the
well-behaved component $B_t$, and not by the function $A_\varphi$.
In this respect, the desire to automatically associate the
singularities of the auxiliary potentials with the intrinsic
properties of the dyonic black hole would have forced the authors of
\cite{CGa}, after using a different representation of
$T{^\mu}{_\nu}$ involving say the dual electromagnetic tensor
$\tilde F^{\mu\nu}$ only, to draw a new conclusion that it is the
{\it electric} string singularities of the component $B_\varphi$
that contribute to the expression of the Komar mass, with zero
contribution coming from magnetic charge.

As has already been shown in \cite{RGM}, the choice of the
energy-momentum tensor $T{^\mu}{_\nu}$ in the symmetrical
representation
\begin{equation} T{^\mu}{_\nu}=\frac{1}{8\pi}(F^{\mu\alpha}F_{\nu\alpha}
+\tilde{F}^{\mu\alpha}\tilde{F}_{\nu\alpha}), \label{Tsym} \end{equation}
where
\begin{equation} F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu, \quad \tilde
F_{\mu\nu}=\partial_\mu B_\nu-\partial_\nu B_\mu,\label{FAB} \end{equation}
permits one to avoid singular sources during the calculation of the
Komar mass integral, reducing the calculational procedure to
exclusively the integrals over the black hole horizon. Although the
paper \cite{RGM} treated the static case, the rotation of the black
hole does not really change the qualitative picture of the
non-rotating model, and below we shall demonstrate that the Komar
mass of the dyonic black hole is obtainable straightforwardly by
means of the original Tomimatsu's mass integral formula (\ref{TMI}),
without the need to consider any singular terms outside the horizon.

In his article \cite{Tom}, Tomimatsu started with the same standard
integral for the calculation of the Komar mass that has been
recently used in the papers \cite{RGM} and \cite{CGa}:
\begin{equation} M_K=\frac{1}{4\pi}\int_{\infty} D^\nu k^\mu
d\Sigma_{\mu\nu}=\frac{1}{4\pi}\int_{\partial \mathcal{M}} D^\nu
k^\mu d\Sigma_{\mu\nu} +\frac{1}{4\pi}\int_{\mathcal{M}}D_\nu D^\nu
k^\mu dS_\mu, \label{KI} \end{equation}
with the same decomposition into the surface and bulk integrals.

By choosing the horizon of the black hole as $\partial \mathcal{M}$,
Tomimatsu computed the first integral on the right-hand side of
(\ref{KI}) and obtained
\begin{eqnarray} \frac{1}{4\pi}\int_{H} D^\nu k^\mu d\Sigma_{\mu\nu}&=&
\frac{1}{8\pi}\int_{H}[-\omega\partial_z\chi+
2\omega{\rm Im}(\Phi\partial_z\bar\Phi)]d\varphi dz \nonumber\\
&=&\frac{1}{8\pi}\int_{H}[-\omega\partial_z\chi+
2\omega(A_t\partial_z B_t-B_t\partial_z A_t)]d\varphi dz,
\label{TI1t} \end{eqnarray}
and he also rewrote the bulk integral on the right-hand side of
(\ref{KI}) in the form
\begin{equation} \frac{1}{4\pi}\int_\mathcal{M}D_\nu D^\nu k^\mu dS_\mu
=-2\int_\mathcal{M} T{^t}{_t}\sqrt{-g}\,d^3x, \label{TI2t} \end{equation}
and the correctness of formulas (\ref{TI1t}) and (\ref{TI2t}) was
not objected in \cite{CGa}. The authors of \cite{CGa}, however,
questioned Tomimatsu's result of computing the integral on the
right-hand side of (\ref{TI2t}), namely,
\begin{equation} -2\int_\mathcal{M} T{^t}{_t}\sqrt{-g}\,d^3x=
-\frac{1}{4\pi}\int_{H}\omega{\rm
Im}(\Phi\partial_z\bar\Phi)d\varphi dz, \label{TIr} \end{equation}
which, together with (\ref{TI1t}), gives formula (\ref{TMI}). Though
they rightly pointed out that the representation (\ref{TE1}) of the
energy-momentum tensor used by Tomimatsu requires additionally
taking account of two singular string sources, which modifies the
horizon contribution (\ref{TMI}) of the Komar mass, they still
erroneously ascribed the {\it formal} mass distribution due to
singularities of the auxiliary function to the genuine dyonic
Kerr-Newman space. Actually, we have a strong impression that
Tomimatsu obtained his formula (\ref{TMI}) after deliberately
suppressing the additional electromagnetic term discussed in
\cite{CGa}, with the idea to get a physically consistent expression
for the Komar mass of a black hole. On the other hand, the authors
of \cite{CGa} have restored the additional electromagnetic term in
Tomimatsu's formula (\ref{TMI}) for mathematical consistency, but
this has led them to the physically incorrect result for the mass
distribution in a dyonic black hole.

Remarkably, the validity of Tomimatsu's mass integral (\ref{TMI})
can be readily demonstrated by employing the symmetrical
representation of the electromagnetic energy-momentum tensor
(\ref{Tsym}) for the evaluation of the integral on the right-hand
side of (\ref{TI2t}). Then, following the steps outlined in the
paper \cite{RGM} for that representation, the bulk integral on the
left-hand side of (\ref{TIr}) reduces to the surface integral over
the horizon, yielding
\begin{equation} -2\int_\mathcal{M} T{^t}{_t}\sqrt{-g}\,d^3x=
\frac{1}{4\pi}\int_{H}(A_t\partial_z B_\varphi-B_t\partial_z
A_\varphi)d\varphi dz, \label{TIc} \end{equation}
where, at the last stage of the computation, we have used the
substitutions
\begin{equation} \rho^{-1}f[(\rho^2 f^{-2}-\omega^2)\partial_\rho A_t
-\omega\partial_\rho A_\varphi)] =-\partial_z B_\varphi \label{dzB}
\end{equation}
and
\begin{equation} \rho^{-1}f[(\rho^2 f^{-2}-\omega^2)\partial_\rho B_t
-\omega\partial_\rho B_\varphi)] =\partial_z A_\varphi, \label{dzA}
\end{equation}
the latter relation being the corollary of the first equations of
the systems (\ref{BtA}) and (\ref{BfA}).

Now, combining formulas (\ref{TI1t}) and (\ref{TIc}) in one, and
also taking into account that $\omega$ assumes constant value on the
horizon, we get for the Komar integral (\ref{KI}) the expression
\begin{equation} M_K= \frac{1}{8\pi}\int_{H}[-\omega\partial_z\chi+ 2\omega
A_t\partial_z(B_t+\omega^{-1}B_\varphi) + 2\omega
B_t\partial_z(-A_t-\omega^{-1}A_\varphi)]d\varphi dz, \label{KMp}
\end{equation}
and lastly, after noting that the second and third terms in the
integrand of (\ref{KMp}) vanish because these contain the
derivatives of the potentials $\Phi_e$ and $\Phi_m$, both potentials
taking constant values on the horizon, we obtain the final
expression for the Komar mass
\begin{equation} M_K=-\frac{1}{8\pi}\int_H \omega\partial_z\chi\, d\varphi dz,
\label{KMn} \end{equation}
which fully coincides with Tomimatsu's formula (\ref{TMI}).

Therefore, the use of the symmetrical representation (\ref{Tsym}) of
the electromagnetic energy-momentum tensor during the calculation of
the Komar mass integral leads straightforwardly to the original
formula obtained by Tomimatsu in the paper \cite{Tom}. We think this
gives us a nice example of a brilliant physical intuition prevailing
over scholastic mathematical estimates.

\section{Discussion and conclusions}

The derivation of formula (\ref{KMn}) involving exclusively the
integrals over the event horizon unequivocally suggests that the
whole Komar mass evaluated in this way is located inside the horizon
of the black hole. In this respect, it seems remarkable that in the
generic expression (\ref{KI}) for the Komar mass the integration is
set to be performed over a sphere of infinite radius, thus giving an
opportunity to use, if necessary, singular functions during the
computational process. The presence of the electromagnetic field
obviously complicates the evaluation of the Komar mass, both
technically and conceptually, compared to the pure vacuum case
since, as we have seen in our previous paper and in the present one,
the correct choice of the representation of the electromagnetic
energy-momentum tensor is required to avoid the presence of
artificial singularities in the dyonic black holes; consequently, in
the case when an unsymmetrical representation of the energy-momentum
tensor is employed, a very accurate physical interpretation of the
results obtained is needed. Thus, the use of the representation
(\ref{TE1}) in the paper \cite{CGa} urged the authors of that paper
to evaluate the mass integral (\ref{KI}) with the help of the
pathological $\varphi$-component of the potential $A_\mu$. So, it is
not a surprise that they could only arrive, within the framework of
their approach, at the mass distribution spreading along the whole
symmetry axis, and this purely technical result was erroneously
claimed by them an intrinsic property of the dyonic Kerr-Newman
black hole. At the same time, what those authors really did was
simply calculating in a not rational way the same value of the Komar
mass (located entirely inside the black hole horizon) that otherwise
follows directly from Tomimatsu's formula (\ref{TMI}) when the
symmetrical representation of the energy-momentum tensor is used. It
is also clear that since a certain part of the total Komar mass $m$
calculated in the paper \cite{CGa} for the Kerr-Newman dyon comes
from the string singularities, the horizon contribution there
differs from the value obtainable by means of Tomimatsu's formula in
the absence of Dirac strings, which explains the appearance of the
additional electromagnetic term in the mass formula of the paper
\cite{CGa}.

Summarizing the results obtained in our short series of two papers,
it should be first of all pointed out that the knowledge of only the
4-potential $A_\mu$ is generically not sufficient for a correct
description of the electromagnetic field which also requires the
knowledge of the dual 4-potential $B_\mu$. In the case of stationary
axially symmetric fields, these potentials $A_\mu$ and $B_\mu$ have
the nonzero $t$- and $\varphi$-components, namely $A_t$,
$A_\varphi$, $B_t$ and $B_\varphi$, among which it is the
$t$-components $A_t$ and $B_t$ that are the basic key functions
defining the physical properties of the electric and magnetic field,
respectively, in particular their singularity structure, while the
$\varphi$-components $A_\varphi$ and $B_\varphi$ play an auxiliary
role in the description of the electromagnetic field, and the
singularities of the functions $A_\varphi$ and $B_\varphi$ are not
characteristic of the proper electric or magnetic field.

We have shown that the use of a specific representation of the
electromagnetic energy-momentum tensor is able to provoke erroneous
interpretations of the physical properties of dyonic black holes:
thus, the choice of the canonical representation (\ref{TE1}) for
$T{^\mu}{_\nu}$ in the Komar mass integral leads to appearance of
magnetic Dirac strings \cite{Dir} as the sources of mass, while the
representation of $T{^\mu}{_\nu}$ involving only the dual
electromagnetic tensor $\tilde F^{\mu\nu}$ (see formula (8) of
\cite{RGM}) gives rise to massive Dirac strings generated by the
electric charge. This naturally singles out the symmetrical
representation (\ref{Tsym}) of $T{^\mu}{_\nu}$ as the most
appropriate one for the dyonic solutions because no contributions
due to string singularities emerge during the evaluation of the mass
integral with the help of (\ref{Tsym}).

It follows directly from our analysis that Dirac strings (magnetic
and electric ones) must be excluded from the physical picture of
dyonic spacetimes. Nevertheless, the semi-infinite singularities
that are characteristic mathematical attributes of the components
$A_\varphi$ and $B_\varphi$ in the presence of nonzero magnetic and
electric net charges still remain a legitimate part of the general
mathematical toolkit and are expected to be taken into account as
purely mathematical objects in some calculations involving the
functions $A_\varphi$ and $B_\varphi$.

Bearing in mind our basic idea that electromagnetism is necessarily
a theory of two electromagnetic potentials, we have slightly
enlarged the well-known Ernst formalism by explicitly introducing
into it the components $B_t$ and $B_\varphi$ of the dual
electromagnetic potential $B_\mu$. This improves the formalism in
two ways. First, it permits now a unified symmetrical treatment of
the electric and magnetic fields, in particular the introduction for
the first time of the magnetic potential $\Phi_m$ which takes
constant value on the horizon, half a century later than Carter's
electric potential $\Phi_e$ \cite{Car}. Second, after our amendment,
the Ernst formalism looks not only more complete but also
logistically refined: the Ernst auxiliary magnetic function
$A^{'}_\varphi$, which was needed before just for computing the
``genuine'' component $A_\varphi$ of $A_\mu$, and which we
identified as the $t$-component of the dual potential $B_\mu$, now
plays, alongside $A_t$, the leading role in the description of the
electromagnetic field, while $A_\varphi$ plays the role of an
auxiliary function. This, in our opinion, enriches the Ernst
formalism conceptually, as the knowledge of the electromagnetic
Ernst potential $\Phi=-A_t+iB_t$ supplies us directly with the
explicit expressions of the physical components of the
electromagnetic 4-potentials determining the intrinsic properties of
the electromagnetic field, without the need of finding $A_\varphi$.
We notice in this respect that it is the component $B_t$, and not
$A_\varphi$, that takes part for instance in the definition of the
relativistic multipole moments of the electromagnetic field
\cite{Sim,HPe,SAp,FCH,MMeR}, which gives us another good
illustration of a generic secondary role of the component
$A_\varphi$ in the physical analysis.

We hope that our present paper, as well as the paper \cite{RGM},
presenting some new ideas about the description of magnetic charges,
could be also helpful in the search and experimental detection of
the dyonic sources. Of course, a natural expectation would be that
some known elementary particles, in addition to electric charges
they have, might also carry magnetic charges, such particles thus
being the dyonic objects. Taking the dyonic Kerr-Newman solution
considered in Sec.~II as the simplest model for a stationary dyon,
we observe that the corresponding magnetic dipole moment of the
source is $aq$, while the electric dipole moment is equal to $-ap$,
the latter moment arising due to rotation of the magnetic charge.
Therefore, the presence of the electric dipole moment in elementary
particles might be considered in principle as indirect indication
that the particles are endowed with nonzero magnetic charges.

\section*{Acknowledgments}

This work was partially supported by CONACyT of Mexico and by
``Secretar\'ia de Educaci\'on, Ciencia, Tecnolog\'ia e Innovaci\'on
de la Ciudad de M\'exico (SECTEI)'' of the Mexico City.