Data classification and parameter estimations with deep learning to the simulated time-domain high-frequency gravitational waves detections

Physically, the applied Gaussian beam, which is along the z-direction with the electric field being $E_{x}^{(0)} = \psi(r, z, t, \omega_e)$, is also a strong background noise source. Therefore, to avoid such a strong background noise, we perfectly detect the transverse EMRs (e.g. along the x-direction) of the passing HFGWs (whose transports are assumed to be along the z-direction), given the z-direction electric field of the applied Gaussian beam, i.e. $E_{z}^{(0)} = 0$. This implies that if the detectors are set at the $x,y$ directions, whose photon receiving faces are parallel to the yz places, then the z-direction applied Gaussian beam itself (except its scattering along the other directions) could not be coupled to the detector. Therefore, the transverse EMRs (along the x-direction) of the HFGWs are perfectly detected by the detectors located near the waist of the applied Gaussian beam, and the detectable signals could be calculated by solving the Einstein–Maxwell equation for the proposed configuration shown in figure 1.

In the usual TT gauge/frame, the non-zero elements of the metric of the monochromatic gravitational wave passing through the detection system reads:

$$\begin{eqnarray} &&h_{11} = h_{22} = A_{\small \bigoplus} e^{ik_g\left(z-ct\right)},\, h_{12} = h_{21} = iA_{\small \bigotimes} e^{ik_g\left(z-ct\right)}, \end{eqnarray} \tag{ 2 }$$

with $A_{\small \bigoplus,\bigotimes} = A_g$ and k_g , respectively, are the dimensionless amplitude and wave number of the HFGW, and c is the light velocity. By solving the relevant Einstein–Maxwell equation, the first-order perturbation electric field along the y direction can be obtained as [45]

$$\begin{align} E_{y}^{\left(1\right)}\left(z,t,\omega_{g}\right) & = -\frac{1}{2} A_{g} B_{y}^{\left(0\right)} \kappa_{g} c\left(z+l_{1}\right) \exp \left[\mathrm{i}\left(\kappa_{g} z-\omega_{g} t\right)\right]\nonumber\\ & \quad +\frac{\mathrm{i}}{4} A_{g} B_{y}^{\left(0\right)} c \exp \left[\mathrm{i}\left(\kappa_{g} z+\omega_{g} t\right)\right], \end{align} \tag{ 3 }$$

near the waist of the applied Gaussian beam. Above, $E_{y}^{(0)}(x, z, t, \omega_{e}) = 2 x(l/W^{2}-i\kappa_{e}/(2R)\int E_{x}^{(0)}(x, y, z, t, \omega_e){\mathrm{d}}y$ and $B_{z}^{(0)}(x,y,z,t,\omega_{e}) = i\omega_{e}^{-1}(\partial E_{x}^{(0)}(x, y, z, t, \omega_{e})/\partial y-\partial E_{y}^{(0)}(x, z, t, \omega_{e})/\partial x),$ are the y-direction electric field and the z-direction magnetic field of the applied Gaussian beam, respectively.

Strictly speaking, as the applied Gaussian beam is naturally defined in the detector proper frame (PDF), the detectable signals should be calculated by solving the Einstein–Maxwell equation in the PDF [47–49], wherein the non-zero elements of the metric of the passing HFGW alternatively reads

$$\begin{eqnarray} &&\left\{ \begin{array}{ll} h_{00} = 2R_{0i0j}x^ix^j\left(\frac{-1+e^{-i\omega_gz}}{\omega_g^2 z^2}+\frac{i}{\omega z}\right)\\ \\ h_{ij} = 2R_{ikjl}x^kx^l\left(\frac{2i\left(1-e^{-i\omega_gz}\right)}{\omega_g^3 z^3}+\frac{i\left(1+e^{-i\omega_gz}\right)}{\omega^2 z^2}\right)\\ \\ h_{0i} = 2R_{0jik}x^jx^k\left(\frac{i\left(1-e^{-i\omega_gz}\right)}{\omega_g^3 z^3}+\frac{\left(e^{-i\omega_gz}-i/2\right)}{\omega^2 z^2}\right), \end{array} \right. \end{eqnarray} \tag{ 4 }$$

with $R_{ijkl} = [h_{il,jk}+h_{jk,il}-h_{jl,ik}-h_{ik,jl}],\,i,j,k,l = 0,1,2,3$ is the linearized Riemann tensor. Again, we solve numerically the Einstein–Maxwell equation with the above metric for the same detection configuration and then get the time-domain PPF signals in the PDF.

Obviously, if the frequency of the Gaussian beam equals to that of the passing HFGW, i.e. $\omega_{g} = \omega_{e}$, the density of the first-order transverse PPF along the x direction can be calculated by:

$$\begin{equation} \begin{aligned} &n_{x}^{\left(1\right)}\left(x,y,z,t\right) = \frac{1}{\hbar\omega_{e}\mu_{0}}Re\left[E_{y}^{\left(1\right)}\left(\omega_e,\omega_g,t\right) B_{z}^{\left(0\right)}\left(\omega_e,t\right)\right]. \end{aligned} \end{equation} \tag{ 5 }$$

Here, $E_{y}^{(1)}(\omega_e,\omega_g,t)$ and $B_{z}^{(0)}(\omega_e,t)$ are the electric fields of the first-order PPF signal and the magnetic field of the applied Gaussian beam. Correspondingly, the first-order energy flux density of the first-order PPF is along the x-direction, i.e. the transverse of the z-direction transporting HFGW. Thus, $n_{x}^{(1)}$ refers to the transverse first-order PPF, which is certainly related to the physical parameters such as the $\omega_{e},\omega_{g}$, and A_g , etc. Specifically, for a photon detector located at x = 0 (i.e. at the waist of the applied Gaussian beam), figure 3 shows the time-domain photon number distributions of the transverse first-order PPF signals by numerically solving the Einstein–Maxwell equation, either in the usual TT gauge/frame (blue line) or in the PDF (red line). Obviously, the signal behaviors calculated either in the usual TT gauge/frame or in PDF are almost the same (although their strengths are different), and the direct response (without any scattering), shown as the black line in figure 3, of the applied Gaussian beam (served as background noise) is negligible (as it transports along the z direction and thus its polarization is along the x-direction). Noted that the ±-sign in the vertical axis of the figure 3 refers to the positive- and negative of the x directions, respectively. As the detector is located at x = 0, the detected number for the positive x-direction photons and that for negative x-direction are the same, i.e. $|n_{x = 0}^{(1)}|$. Given the account rate of the current single-photon detector, which can reach GHz order, the signals are really sufficiently strong to be detected experimentally.

Zoom In Zoom Out Reset image size

Consequently, at the detection time $t = t_s$, the photon number of the first-order PPF signals detected by the detector (located at x = 0) can be calculated as can be calculated as [45], i.e.

$$\begin{equation} \begin{aligned} &N_{x = 0}^{\left(1\right)}\left(x = 0,t_s\right) = \int_{t_s}^{t_s+t_{0}}\iint_{\Delta s} |n_{x}^{\left(1\right)}\left(x = 0,t\right)|{\mathrm{d}}y{\mathrm{d}}z{\mathrm{d}}t, \end{aligned} \end{equation} \tag{ 6 }$$

with t₀ being the detection duration, Δs the receiving area of the detector, and $t_s = kt_{0},\,k = [0,1,2, \ldots, 100]$ are the time series of the detections.

As mentioned above, the behaviors of the first-order transverse PPFs calculated in either the usual TT gauge/frame or in the PDF are almost the same, although their strengths are different. Basically, the different signal strengths refer only to the choice of the signal-noise ratio (SNR) during the detected data processing. Specifically, let us first assume that the time domain signals $|n^{(1)}_{x = 0}(t_s)$, generated by using a photon detector with the photon receiving surface $0\lt y\lt0.15m$, $0\lt z\lt0.1m$ (thus $\Delta s = 0.015$ m²) [33], has been detected by the single-photon detector with the duration t₀ (e.g. $t_{0} = 0.08/f_{e}$ here). Next, the noise data should be generated, as noise is unavoidable in any signal detection. For simplicity, in the present work we only consider the Gaussian white noises (GWNs). Therefore, the simulatively detected data for the processing can be generated by superposing the signals, shown in figure 3 with the simulated GWNs [57], with the controllable SNR [58]:

$$\begin{eqnarray} &&{\mathrm{SNR}} = \frac{\sum_{i = 1}^Ns_i^2\left(t_s\right)}{\sum_{i = 1}^Nn_i^2\left(t_s\right)}. \end{eqnarray} \tag{ 7 }$$

Here, N is the data quantity, $s_i(t_s)$ and $n_i(t_s)$ are the amplitudes (i.e. photon number) of the first-order PPF signals and the added GWNs, respectively. In practice, signal detection is performed by multiple repeats, and thus these amplitudes are the accumulation of the many detection events typically such as ${\sim}10^4$ ones. In the present work, we assume that the above data amplitudes are generated by the accumulation of 10⁴ detected events. Each of the time-domain detected data is obtained by using the single-photon detector (located at x = 0) to detect the transverse first-order PPF signal, with the account rate being at GHz order. This is why the data amplitude in figure 4 is now sufficiently large for processing.

Zoom In Zoom Out Reset image size

The simulated time domain detection data can be generated as follows: First, we extract 101 dots from the blue line in figure 4, with $t_s = [0, t_0, 2\times t_0, \ldots, 100\times t_0]\times 10^{-9}$s, to get a $N = 1\times101$-dimension data as the pure time-domain first-order PPF signals. Here, the frequency and dimensionless amplitude of the simulatively-detected GW are set as: $\omega_{g} = 2\pi\times10^{9}$ Hz and $A_{g} = 1.0\times10^{-31}$, respectively. Then, we add a GWN to these data with the SNR = $[0,1]$ and get a simulated 101-dimension time-domain detection data.

In our study, we aim to generate $10\,000$ sets of simulated time-domain detection data. This total is divided into two subsets. The first subset, encompassing 8000 sets, is generated by randomly incorporating various high-frequency gravitational wave (HFGW) signals. These signals have amplitudes within the range: $A_g = 1.0\times[10^{-31}, 10^{-29}]$ and frequencies within the range: $\omega_g = 2\pi\times[10^{8}, 10^{10}]$ Hz, this data is then combined with GWNs, featuring an SNR of $(0,1]$, and designated as positive samples. Each data set within these positive samples is labeled as 1. The remaining 2000 sets form the second subset, composed of GWNs with varying means and variances. These sets are generated to function as negative samples, with each data set marked as 0. After normalizing both the positive- and negative samples, we disrupt them to generate the simulated time-domain detection data of $10\,000$ sets, shown in figure 5(a). As seen in figure 5(b), the generated time-domain data is split into the training-, validation-, and the test datasets, with the ratio is set as $6:2:2$.

Zoom In Zoom Out Reset image size

We use a two-class CNN model, whose structure and the relevant parameters are shown in table 1, to calibrate the positive sample data. Generally, a CNN model mainly contains the convolution-, ReLU-, maxpooling-, flatten-, linear-, dropout- and softmax layers. Various layers are incorporated within a CNN model, each serving a unique purpose. Convolutional layers are responsible for feature extraction, employing filters that discern patterns within the data. Rectified Linear Units or ReLU layers infuse non-linearity into the model, preserving positive values and eliminating negative values from the convolution output. The complexity of the model is reduced by Max pooling layers, which downscale the dimensions of the input data, retaining only the maximum values from segments of the input. Finally, the output of the trained model can be used to identify if the input $1\times101$ dimension data is the positive sample one or not.

During the training, we add a L2-regularization proceeding and the dropout layer to prevent overfitting. Here, L2-regularization is realized by adding a L2 norm, with the weight λ into the loss function, as its penalty term; the dropout layer is used to temporarily discard certain units from the neural network with certain probabilities. With the CNN architecture shown in table 1, in what follows, we use Tensorflow [59] (see, e.g. www.tensorflow.org) to build the CNN model and then implement the train on the above-mentioned spatial domain training set. By changing the hyperparameters to compare the performance of the trained models, one can optimize the relevant hyperparameters. In the present work, the optimized hyperparameter of the L2-regularization can be obtained after the tuning as λ = 0.01, which implies that the probability of temporarily discarding units is optimized to 0.3. Similarly, the training times (i.e. epochs), the learning rate of Adam optimizer (called lr), and the BatchSize (i.e. the number of training data bars) can also be optimized to the desired values; typically, e.g. epochs = 100, lr = 0.001, and BatchSize = 128. After these hyperparameter optimizations, highly accurate training can be realized. The CNN architecture can be used to train a new CNN model with the datasets in the simulated time-domain detection data generated above. It is shown in figure 6(a) that the accuracies of both the training- and validation sets can reach to 95%, and the loss is decreased to 0.18. This represents that the CNN training model used has sufficiently high accuracies on both the training- and validating sets. Using such a trained CNN model to the test dataset, we get the ROC curve shown in figure 6(b) with the AUC being 0.986, which is very close to the optimal value 1.

Zoom In Zoom Out Reset image size

With the above CNN model, we now develop the data; the frequency and dimensionless amplitude of the inputting HFGWs are assumed as $\omega_g = 2\pi\times10^9$ Hz and $2.0\times10^{-31}$, respectively. The original positive samples with 3000 data are obtained by the pure first-order PPF signal superposed by the GWNs with SNR = 1. The GWNs also generate the negative samples with 2000 data. Then, the simulated detection databases (with a total of 5000 samples) can be obtained by mixing the above positive- and negative samples, which are typically shown in figure 7(a), wherein the dimensions of each data are $1\times101$. Following this, we apply the previously validated, trained CNN model to analyze a 5000-group sample dataset. From this analysis, we successfully identify 3517 instances of effective positive samples. In this case, 2000 data are stochastically selected as the positive sample signals, as shown in figure 7(b). They will be further processed to extract the information of the simulatively detected HFGWs.

Zoom In Zoom Out Reset image size

Certainly, many algorithms can be used to implement the expected data classification. We now compared the performance of the CNN model demonstrated above with the other classifiers, typically such as the K-nearest neighbor (KNN) [60], random forest (RF) [61], support vector machines (SVM) [62], naive Bayes classification (NB) [63] and logistic regression (LR) [64] algorithms, etc. For the comparison, we use the same test and training sets generated in section 3.1, and then run the different algorithm for the classification. The accuracies using different algorithms are listed in figure 8. It is seen that, compared with the other machine learning algorithms, the accuracy of the used CNN model is higher than that of the other algorithm by about 10%.

Zoom In Zoom Out Reset image size

We now discuss the possible improvements by using the trained model, compared with the usual matched filtering algorithm without the CNN pre-processing. Given that the matched filtering algorithm requires significantly powerful computing power and relatively long processing times, we simply use the Pearson correlation coefficient, defined by equation (8), to quantify the correlation between the same two vectors. It is shown in figure 9 that, compared with the usual matched filtering algorithm, the execution time of the CNN is shorted by almost two orders. However, the classification accuracy is almost the same.

Zoom In Zoom Out Reset image size

Firstly, we discuss how to extract the frequency information of the simulatively detected HFGWs from the classified time-domain detection data shown in figure 10(b) by using the frequency matchings. Suppose that the frequency of the simulatively detected HFGW is in the range of $\omega_{g} = 2\pi\times[0.84, 1.16]\times10^{9}$ Hz and the dimensionless amplitude is $1.0\times10^{31}$, we get the template library data by using equation (6) to generate 2000 data (their frequency difference is the same) of the pure first-order transverse PPF signals (without any noise). Figure 10 shows the established frequency template library, wherein each data set has a defined frequency label and each point represents template data, whose dimension is also set as $1\times101$.

Zoom In Zoom Out Reset image size

It is known that, the correlation between two vectors: $X = (X_{1}, X_{2}, \ldots, X_{n})$ and $Y = (Y_{1}, Y_{2}, \ldots, Y_{n})$, can be described by the Pearson correlation coefficient [65]:

$$\begin{eqnarray} &&r = \frac{\sum_{i = 1}^{n}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sqrt{\sum_{i = 1}^{n}\left(X_{i}-\bar{X}\right)^{2}} \sqrt{\sum_{i = 1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}}. \end{eqnarray} \tag{ 8 }$$

Here, $\bar{X}$ and $\bar{Y}$ are the mean values of X- and Y vectors, respectively. Obviously, the value range of this coefficient is: $r\epsilon [-1,1]$; the larger the value of r corresponds, the higher the matching degree, and thus r = 1 means the perfect matching. On the contrary, $r\unicode{x2A7D} 0$ represents the mismatch of the two vectors, and $r = -1$ means that the two vectors are negatively correlated (i.e. their parameters are completely independent).

By matching each data in the positive sample database with the data in the frequency template library shown in figure 10 (with also $N_t = 2000$ data) one by one, the frequencies corresponding to the largest Pearson correlation coefficients are output as the estimated frequencies of the simulatively detected HFGWs. It is seen from figure 11 that the estimated frequencies reveal a statistical distribution, which can be fitted as a Gaussian curve (i.e. the red line in figure 11). The vertical axis in figure 11 represents the number of positive sample data.

Zoom In Zoom Out Reset image size

The accuracy of the above frequency parameter estimations is described below. According to the central limit theorem [66], the errors deviate stochastically from their average value. Therefore, the confidence level of the estimation is defined as the probability that the estimated value of the parameter falls in the corresponding area [67]. The confidence interval [67] refers to the error range between the sample statistical value and the overall parameter value at a certain confidence level. Of course, a larger confidence interval means a higher confidence level. From the statistical distribution shown in figure 11, we argue that the mean value of the estimated frequency is $\bar{\omega}_g = 2\pi\times1.0\times10^9$ Hz, which is the same as the input frequency: $\omega_g = 2\pi\times1.0\times10^9$ Hz of the simulatively detected HFGW. The standard deviation is calculated as: $\delta \omega_g = 2\pi\times3.046\times 10^{6}$ Hz. With the present 2000 sample amounts, the confidence interval in: ($\bar{\omega}_g-z\delta\omega_g$, $\mu+z\delta\omega_g$) = $2\pi\times(0.9958, 1.0042)\times 10^{9}$ Hz, with the variance multiple z = 1.96 [67] (describing the degree of deviation from the variance) at 95% confidence level can be obtained.

Continuously, we discuss how to extract the dimensionless amplitude of the simulatively detected HFGWs from the positive samples shown in figure 7(b). Similarly, we first generate the amplitude template library of the HFGWs wanted to be detected. It is seen from equation (7) that the amplitude of the input GW is linearly related to the amplitude of the first-order PPF signal without noise. Given that the frequency information of the simulatively detected HFGWs had been extracted by the above frequency matchings, the amplitude template library can be generated simply as follows. Inputting a series of HFGWs with different dimensionless amplitudes but the same frequency $\omega_g = \bar{\omega}_g$ (i.e. the mean value of the estimated frequency by the above frequency matchings), a series of the noise-free first-order transverse PPF signals generated simulatively are served as the amplitude templates library. Here, the amplitudes of the input HFGWs are selected as in the range: $A_g = [10^{-32}, 10^{-28}]$ with the step being $\delta A_g = 10^{-32}$. Figure 12 shows specifically the template data of the simulatively detected HFGW (with the estimated frequency: $\omega_a = \bar{\omega}_g$), whose dimensionless amplitudes are assumed to be in the range: $[10^{-32}, 10^{-28}]$. In the figure, each point represents a template data with the dimensions being $1\times101$.

Zoom In Zoom Out Reset image size

Different from the Pearson correlation coefficients used above for the frequency matching, the amplitude information can be extracted simply by calculating the amplitude data correlation between the positive sample data and the template data. This is measured by the Euclidean distance [68]:

$$\begin{eqnarray} &&d = \sqrt{\sum_{i = 1}^{n}\left(X_{i}-Y_{i}\right)^{2}}, \end{eqnarray} \tag{ 9 }$$

between the vector $X = (X_{1},X_{2},\ldots,X_{n})$ in the positive sample and the vector $Y = (Y_{1},Y_{2},\ldots,Y_{n})$ in the template. Certainly, the smaller the Euclidean distance corresponds, the greater the magnitude correlation. To this end, we match each data in the positive sample database with the amplitude template library one by one. Alternatively, a set of template data with the smallest Euclidean distances is treated as the possible amplitude values of the simulatively detected HFGWs. The matched results are also shown as a statistical distribution, shown schematically in figure 13. Here, the horizontal axis represents the output matching amplitudes, and the value on the vertical axis is the number of the positive sample data whose amplitudes fall in the amplitude range of the corresponding histograms. Again, the estimated amplitudes can also be fitted as a Gaussian curve, shown as the red line in figure 13. It shows that the mean value of the estimated amplitudes is: $\bar{A}_g = 2.733\times 10^{-31}$, which is very close to the input value: $A_g = 2.0\times10^{-31}$, with the standard deviation being $\delta A_g = 1.572\times 10^{-32}$. The confidence level of such an estimation is at 95% in the confidence interval: $(2.4474, 3.4546)\times 10^{-31}$.

Zoom In Zoom Out Reset image size

The reliability of the data processing method presented above could be tested by using multiple datasets for the different simulative detections.

In fact, much detection data can be simulatively generated for the different parameters of the HFGWs, as well as the background noises. For example, for the simulative detecting of the HFGWs, whose frequencies are in the range $\omega_g = 2\pi\times [10^8, 10^{12}]$ Hz and the dimensionless amplitudes are in the range $A_g = 1.0\times[10^{-33}, 10^{-29}]$, we can specifically generate the simulative test dataset shown in figure 14(a). The ROC curve of the trained CNN model to identify the positive samples with the test dataset is shown in figure 14(b), with the AUC being 0.978. One can see again that the trained CNN model possesses a sufficiently high resolution to identify whether the test dataset contains the HFGW signals accurately. Thus, the parameter estimation for the simulatively detected HFGWs is significatively accurate.

Zoom In Zoom Out Reset image size

Certainly, for the input of any HFGW signal, the corresponding positive samples (containing the HFGW-induced first-order PPF signals) with 2000 or more groups of data can be obtained. By matching them with the data in the corresponding frequency- and amplitude template libraries, the confidence interval of amplitude and frequency in the HFGW signals can be effectively estimated. In table 2, we also list a series of estimated results at $95\%$ confidence level. This indicates that the proposed data processing method, based on deep learning, is feasible for the practice EMR detections of the HFGWs.

Table 2. Multi-data testings to extract the parameters of the simulatively detected HFGWs.

Frequency estimations		Amplitude estimations
ωg ( Hz)	The confidence interval at 95% confidence level	Ag	The confidence interval at $95\%$ confidence level
$2\pi\times3.0\times10^{8}$	$2\pi\times(2.988, 3.012)\times 10^{8}$	$2.0\times10^{-31}$	$(2.468, 2.958)\times10^{-31}$
$2\pi\times5.0\times10^{8}$	$2\pi\times(4.980, 5.021)\times 10^{8}$	$5.0\times10^{-31}$	$(6.234, 7.381)\times 10^{-31}$
$2\pi\times7.0\times10^{8}$	$2\pi\times(6.974, 7.026)\times 10^{8}$	$7.0\times10^{-31}$	$(8.726, 10.338))\times 10^{-31}$
$2\pi\times3.0\times 10^{9}$	$2\pi\times(2.989, 3.011)\times 10^{9}$	$2.0\times10^{-30}$	$( 2.489 ,2.973)\times 10^{-30}$
$2\pi\times5.0\times10^{9}$	$2\pi\times(4.953, 5.051)\times 10^{9}$	$5.0\times10^{-30}$	$(6.237, 7.383)\times 10^{-30}$
$2\pi\times7.0\times10^{9}$	$2\pi\times(6.972, 7.078)\times 10^{9}$	$7.0\times10^{-30}$	$(1.05, 1.25)\times 10^{-29}$