\section{Introduction}\label{sec:introduction}
}

\IEEEPARstart{M}{edical} imaging plays an increasingly prominent role in clinical diagnosis. Multimodal medical images can provide information from various aspects, thus helping physicians confirm the diagnosis and make decisions regarding future treatment. For example, Magnetic Resonance Imaging (MRI) images provide the soft-tissue structures of the body; computerized tomography (CT) scans provide bone structure and high-density tissue information; Positron Emission Tomography (PET), and Single-Photon Emission Computed Tomography (SPECT) images can show the metabolic activity of the cells of tissues. To acquire adequate information, physicians must analyze many different medical images, which is time-consuming and laborious. Multimodal medical image fusion (MMIF) can merge the complementary information of original images and present the required information in one fused image. This is clinically significant because physicians can now access more detailed information about disease-related changes, thus providing patients with more comprehensive and precise medical treatment and support.

In recent years, many deep learning-based approaches have been proven successful in image fusion. These methods can extract features from input images and construct a fused image with the information needed. For instance, Wang et al. \cite{wang2020multi} introduced a CNN-based medical image fusion algorithm. This method employed the trained Siamese convolutional network to fuse the pixel activity information and implemented a contrast pyramid to decompose the source images. For more generalized image fusion, Xu et al. \cite{xu2020fusiondn} proposed an unsupervised and unified densely connected network, FusionDN, which is capable of various kinds of image fusion tasks.

In this work, we propose a novel end-to-end feature learning framework for multimodal medical image fusion with clear edge information and detailed textures. To be more precise, we focus on anatomical (CT) and functional (MRI) image fusion. The framework contains a feature extractor, a fixed fusion strategy that does not involve any adjustable parameters, and image reconstruction. Our approach combines input images and generates a fused image with more detailed textures and less information loss. Our fusion results exceed the state-of-the-art performance both qualitatively in subjective vision and quantitatively in multiple fusion metrics.

The main contributions of this paper are summarized as follows:

\begin{enumerate}
    \item We propose a novel Dilated Residual Attention Network (DILRAN) for the feature extraction module. DILRAN coalesces the advantages of the residual attention network, the pyramid attention network, and the dilated convolutions. The proposed network has a faster convergence speed and can extract multi-scale deep semantic features.
    \item We introduce a novel fusion strategy, Softmax Feature Weighted Strategy, which achieves a good result and outperforms other fusion strategies.
    \item Extensive experiments show our proposed framework and fusion strategy exceed the state-of-the-art performance based on objective fusion metrics and subjective image quality.
\end{enumerate}

\section{Related Work}

There are several existing approaches for imaging fusion. Most traditional image fusions are based on the transform domain and are at pixel, feature and decision-level \cite{HERMESSI2021multimodal}. The usage of fuzzy logic and neural networks with multi-scale decomposition can handle uncertainties and improve efficiency for fusion images. Other iconic works for utilizing traditional image processing algorithms for fusion tasks include Possion Image Editing \cite{perez2003poisson}, cross bilateral filter \cite{shreyamsha2015image}, and non-subsampled contourlet transform \cite{da2006nonsubsampled}. Tian et al.\cite{tian2016multi} proposed an improved version of pulse coupled neural network (PCNN) to manage NSCT sub-images using a shallow learning approach. This new PCNN calculation determines the linking strength parameter using the local area singular value decomposition of the structural information factor.

Recent advances in deep learning (DL) have led to the successful use of convolutional neural networks (CNN) in various imaging tasks such as classification \cite{ilyas2018noise, wang2017residual} and image super-resolution \cite{zhang2017sharp}. CNN has the ability to capture and extract features from images, which can then be used for image reconstruction. For example, Liu et al. \cite{liu2017multi} applied a CNN to multi-focus image fusion, using the network to generate a weight map of pixel activity during fusion. Hermessi et al. \cite{hermessi2018cnn} introduced a CNN-based method in the shearlet domain for extracting feature maps of high-frequency fusion. Li et al. \cite{li2019dense} proposed DenseFuse, a deep learning architecture for infrared and visible image fusion that is trained using dense blocks. Due to its limitation for only working at a single scale, Song et al.\cite{song2019msdnet} proposed a multi-scale DenseNet (MSDNet) to overcome it. They encode the multi-scale mechanism with three filters of varying sizes for effectively capturing features at different scales. Increasing the width of the encoder network can also improve the amount of detail in the fused image. Despite this, DL has great potential in the field of medical artificial intelligence due to its ability to fit complex data, its ability to learn automatically, and its multitask adaptability.

 A multi-generator multi-discriminator conditional generative adversarial network is presented by Huang et al. \cite{huang2020gan} to fuse functional information and structural information including texture details and dense structure information. This network presents better visual effects and also preserves the approximate maximum amount of information in several MMIF datasets. Fu et al. \cite{fu2021multiscale} introduced the residual pyramid attention structure: MSPRAN, which combines the advantages of residual attention and pyramid attention mechanisms in the fusion task. The framework extracts and keeps more information than a single residual attention or pyramid attention mechanism as the number of layers increases and maintains better deep features and expression capabilities. Another network structure (MSDRA) on double residual attention \cite{li2022multiscale} combines a residual network and attention to acquire important detailed features while avoiding network gradient vanish or explosion.

\section{Proposed Method}

In this section, we provide a thorough discussion of the fusion framework we proposed and the corresponding loss function.
\subsection{Overview}

Figure \ref{overall} summarizes the overall pipeline. We introduce a novel end-to-end framework for medical image fusion. Our proposed fusion algorithm consists of a feature extractor, a fusion module, and a reconstruction module. The feature extractor aims to formulate deep semantic features of input images ($I_1$, $I_2$), and then use them as inputs to the fusion module. We introduce Dilated Residual Attention Network (DILRAN) for the feature extraction module. Next, the fusion module aims to fuse the two extracted feature maps into one map containing features from both original feature maps. The module determines for a certain pixel in the fused feature map, whether it comes from $I_1$, $I_2$, or both. The fusion rule can have many possible solutions, i.e., $max(I_{1}^i, I_{2}^i)$ for the $i$th pixel or the weighted average of both. Finally, the reconstruction module aims to reconstruct the fused image from the output of the fusion module. A sequence of convolutional layers is used when reconstructing the final fused image. Section \ref{fe}, \ref{fs}, and \ref{recon} provide an in-depth explanation of each module we mentioned above.

\begin{figure*}[h]
	\begin{center}
		\includegraphics[width=.80\textwidth]{figures/overall.png}
	\end{center}

	\caption{The overall pipeline of the proposed method, which consists of a feature extractor, a fuser, and an image reconstructor}
	\label{overall}

\end{figure*}

\subsection{Dilated Residual Attention Network}

The design principle of DILRAN is inspired by two state-of-the-art mechanisms: residual attention \cite{wang2017residual} and pyramid attention \cite{li2018pyramid}, as well as a recent paper \cite{fu2021multiscale} that relies on the residual pyramid attention mechanism. The classical residual attention network for image classification tasks \cite{wang2017residual}, which is realized by imposing the attention mechanism in the residual network architecture. Figure \ref{resid_attn} shows an example of the residual attention network. The top network path is called the trunk branch, which consists of a sequence of convolutional layers for feature processing. The bottom path is called the soft mask branch, which is formed by a downsample-upsample block with an activation function. The combination of the trunk and the soft mask branch build up the residual attention network following Equation \eqref{resid_eq}.

\begin{figure}[H]
	\begin{center}
		\includegraphics[width=.5\textwidth]{figures/resid_attn.png}
	\end{center}

	\caption{Residual attention network architecture.}
	\label{resid_attn}

\end{figure}

\begin{equation} \label{resid_eq}
    F(x) = (1 + S(x)) * T(x)
\end{equation}

Where $x$ is the feature map from the previous layer, $F(\cdot)$ is the output feature map, $S(\cdot)$ is the function for the soft mask branch, and $T(\cdot)$ is the convolution operation.

The advantage of the residual attention network is that it has a faster convergence speed without gradient vanishing or explosion, however, it may not be able to extract and learn deep features. Hence, the other mechanism we adopt is the pyramid attention network \cite{li2018pyramid}. To be more specific, we focus on the feature pyramid attention that could provide deeper semantic features and attention advantages to the output feature map. A sample pyramid attention network architecture is shown in Figure \ref{pyram_attn}. The convolution block usually contains a sequence of convolutional layers that capture features at different scales and receptive fields. We adapt the idea to replace convolutions with larger kernel filters with a sequence of convolutions with smaller kernel filters \cite{szegedy2016rethinking}. Hence, the $CB1$ in our method is a single $3 \times 3$ convolutions; the $CB2$ consists of \textit{two} $3 \times 3$ convolutions (represents a single $5 \times 5$ convolutions), and the $CB3$ consists of \textit{three} $3 \times 3$ convolutions (represents a single $7 \times 7$ convolutions). The output features can be computed by the following Equation \eqref{pyr_attn}.

\begin{figure}[h]
	\begin{center}
		\includegraphics[width=.5\textwidth]{figures/pyramid_attn.png}
	\end{center}

	\caption{Pyramid attention network architecture.}
	\label{pyram_attn}

\end{figure}

\begin{equation} \label{pyr_attn}
    F(x) = (1 + P_1(P_2(P_3(x)))) * C(x)
\end{equation}

Where $x$ is the feature map from the previous layer, $F(\cdot)$ is the output feature map, $P_i(\cdot), i \in [1,2,3]$ is the parameters of the corresponding $CB_i, i\in [1,2,3]$ in the feature pyramid attention network, $C(\cdot)$ is the convolution operation.

In \cite{fu2021multiscale}, they sequentially downsample the input feature maps to formulate the pyramid structure, and upsample again when performing the element-wise summation as shown in Figure \ref{pyram_attn}. However, this may lose the local information and fine details in the image. To solve this problem, we leverage the ${\{1,3,5\}}$-dilated convolution \cite{yu2015multi} on shallow features of the original input image to extract the multi-scale information instead of downsampling the feature map. The receptive field is expanded using three different dilated convolutions to improve the discriminative multi-scale feature extraction ability of the model. Once the multi-scale features are extracted, we concatenate those features channel-wise, and then the residual-pyramid attention paradigm is used to further extract deep features. The deep features are the output of the feature extraction module and are used in both the fusion and reconstruction modules. Our designed DILRAN architecture is shown in Figure \ref{dilran}.

\begin{figure}[h]
	\begin{center}
		\includegraphics[width=.5\textwidth]{figures/dilran_new.png}
	\end{center}

	\caption{Proposed Dilated Residual Attention Network architecture. Three pyramid convolution blocks contain one, two, and three consecutive $3 \times 3$ convolutional layers, respectively.}
	\label{dilran}

\end{figure}

\subsection{Feature Extraction Module} \label{fe}

The detailed architecture for the feature extractor is shown in Figure \ref{overall}. The input image is first passed to a residual attention network with a $1 \times 1$ convolutional layer to obtain a 64-dimensional shallow feature map. Then, the shallow feature is passed to the DILRAN module to obtain deep semantic features. To ensure the final output feature map contains both the local and global information of the original input image, and to stabilize the training process, another residual attention network is utilized between shallow and deep semantic features. Finally, the output feature map is obtained and used in the following modules.

\subsection{Fusion Strategy} \label{fs}
The fusion strategy in the fusion module is used to fuse the extracted features of input images into a single feature map. In this section, We introduce a novel fusion strategy termed ``Softmax Feature Weighted Strategy'' that exceeds the state-of-the-art performance. The advantage of the proposed fusion strategy is validated in Section \ref{exp_res}.

\subsubsection{Softmax-based weighted strategy}
We obtain two output feature maps $f_1, f_2$ from the extraction module for input images $I_1, I_2$, respectively. The output feature map from the extraction module can be used to generate the corresponding weight map that indicates the amount of contribution of each pixel to the fused feature map \cite{lahoud2019zero}. First, to get the weight map, we take the Softmax \cite{bridle1989training} operation to the feature map which can be realized by Equation \eqref{softmax_op}, where $x_i$ is the $i$th channel of the output feature map $x$.

\begin{equation} \label{softmax_op}
    S(x_i) = \frac{\exp(x_i)}{\sum_{i} \exp(x_i)}
\end{equation}

After we obtained the Softmax output, we compute the matrix nuclear norm ($\|\cdot\|_{*}$), which is the summation of its singular values. Finally, we obtain the weights for the output feature map by taking the weighted average of the maximum value of the nuclear norm. The formula is given in Equation \eqref{sfnn}.

\begin{align} \label{sfnn}
    &W_1 = \frac{\phi(\|S(x_{i})^1\|_{*})}{\phi(\|S(x_{i})^1\|_{*}) + \phi(\|S(x_{i})^2\|_{*})} \\
    &W_2 = \frac{\phi(\|S(x_{i})^2\|_{*})}{\phi(\|S(x_{i})^1\|_{*}) + \phi(\|S(x_{i})^2\|_{*})}
\end{align}

Where $S(x_{i})^j, j \in [1,2]$ is the weight map after Softmax operation for the feature map $f_j, j \in [1,2]$, $\phi(\cdot)$ can be any arbitrary functions. The final fused feature map is then given by $f = W_1 * f_1 + W_2 * f_2$.

\subsection{Reconstruction Module} \label{recon}

The input of the reconstruction module is the fused feature map from the fusion module. The reconstruction module is used to generate the fused human-visible image from the fused feature map while retaining as many details as possible. Our reconstruction module utilizes three consecutive $3 \times 3$ convolutional layers with output channels 64, 32, and 1 to reduce the channel from 64 to 1, indicating the output image is gray-scale. An overview of the reconstruction module is shown in Figure \ref{img_recon}.

\begin{figure}[h]
	\begin{center}
		\includegraphics[width=.5\textwidth]{figures/recon.png}
	\end{center}

	\caption{Image reconstruction module architecture. There are three consecutive $3 \times 3$ convolutional layers, each with the output channel 64, 32, and 1, respectively.}
	\label{img_recon}

\end{figure}

\subsection{Loss function} \label{loss_func}

We hypothesize that in order to make the fused image as close as two input images, we need to minimize the distance between the image reconstructed from deep semantic features and the original input image. Then, after we apply the fusion strategy and obtain the fused deep features, the reconstructor will produce a more realistic and reasonable fused image. Our goal is to coalesce input images to form an output fused image that retains as many details as the input. In other words, the fused image should have a strong correlation with both input images. We argue that the gap between fused images and input images can be mainly measured by the common regression loss function ``mean squared error'' (MSE). We also add gradient loss \cite{ma2020structure, fabbri2018enhancing} to model the fine details of textures in the reconstructed image and perceptual loss \cite{johnson2016perceptual} to model the high-level semantic similarity between reconstructed and input images. In detail, our loss function is defined as follows:

\begin{equation} \label{mse}
    MSE = \sum_{j=1}^{M} \|I_o - I_j\|_{F}^{2}
\end{equation}

\begin{equation} \label{grad_img}
    \nabla_I = \sum_{j=1}^{M} \|\nabla I_o - \nabla I_j\|_{2}^{2}
\end{equation}

\begin{equation} \label{perp_loss}
    Percep = \frac{1}{C*W_i*H_i} \sum_{j=1}^{M} \sum_{k=1}^{C} \|f_{i}^{k}(I_o) - f_{i}^{k}(I_j)\|_{2}^{2}
\end{equation}

\begin{equation} \label{total_loss}
    \mathcal{L}(\theta) = \frac{1}{W*H} (MSE + \lambda_1 * \nabla_I) + \lambda_2 * Percep
\end{equation}

Equation \eqref{mse} is a variant of the MSE loss, where $M$ is the number of input images, $I_o$ is the output image, $\|\cdot\|_F$ is the matrix Frobenius norm. Image gradient loss is given in Equation \eqref{grad_img}, which is realized by the $\mathcal{L}_2$ norm of the image gradient in $x$ and $y$-direction. Equation \eqref{perp_loss} is the perceptual loss \cite{johnson2016perceptual}, where $f_{i}^k(x)$ is the the $k$th channel in $i$th layer (with size $W_i \times H_i$) from the pre-trained VGG16 network \cite{simonyan2014very} with input image $x$, and $C$ is the number of channels. We prefer $i$ to be large, i.e., the deeper layer of the VGG network. Finally, the total loss function is given in Equation \eqref{total_loss}, $\theta$ is the set of network weights to be optimized, $\lambda_1, \lambda_2$ are weight balancing factors of the gradient and perceptual loss, respectively. The total loss has been normalized by the total pixels in the image with width $W$ and height $H$.

\section{Experimental Results} \label{exp_res}

In this section, we will discuss the dataset we used, the experimental setup, evaluation metrics, and the experimental results we obtained.
\subsection{Data}

We adapt the commonly used medical image fusion dataset ``The Harvard Whole Brain Atlas'' \footnote{Available at \url{https://www.med.harvard.edu/aanlib/}} \cite{summers2003harvard} for this work. The dataset we obtained contains three sets of co-registered multimodal medical image pairs: MRI-CT (184 pairs), MRI-PET (269 pairs), and MRI-SPECT (357 pairs). All images are with size $256 \times 256$, MRI and CT are single-channel images with pixel intensities in the range of $[0, 255]$, and PET and SPECT are three-channel images with pixel intensities in the range of $[0, 255]$. In terms of the computational complexity of our model and the time given for this work, we only focus on the MRI-CT fusion task.

\subsection{Experimental Setup}

For the MRI-CT fusion task, we have 184 image pairs. We randomly select 20 image pairs for testing, so that these data are not included in the training process. For the rest of the data, we further split to training and validation data with the ratio of $0.8:0.2$, i.e., 80\% of the data for training, 20\% for validation. All images are normalized in the range of $[0, 1]$ before training. During the training phase, we take out the fusion module, and only train the feature extractor and reconstructor as discussed in Section \ref{loss_func}. The stochastic gradient descent with Adam update rule \cite{kingma2014adam} is utilized as the optimizer, the learning rate is set to 0.0001 and the model is trained for 100 epochs. The weight balancing factors $\lambda_1$ and $\lambda_2$ for image gradient loss and perceptual loss are all set to 0.2, and we use the third VGG block \cite{simonyan2014very} output to calculate the perceptual loss. To prove the effectiveness of the proposed method, we compare it with the following methods: zero-shot learning for medical image fusion \cite{lahoud2019zero}, MSRPAN \cite{fu2021multiscale}, and MSDRA \cite{li2022multiscale}. The parameters of these methods are the default setting suggested by the authors. The experimental environment is a Windows 10 (64 bit) system computer with the Intel Core I7 CPU, NVIDIA GeForce RTX 3050 Ti GPU, and 16 GB of RAM.

\subsection{Evaluation Metrics}

There are many fusion metrics have been proposed in the past few years. However, different metrics reflect the different perspectives of fusion performance. Hence, we select six different, but commonly used fusion metrics in this work: Peak signal-to-noise ratio (PSNR), Structural Similarity (SSIM) \cite{wang2004image}, Feature SSIM (FSIM) \cite{zhang2011fsim}, Mutual information (MI) \cite{qu2002information}, pixel-wise feature MI ($FMI_{pixel}$) \cite{haghighat2011non}, and Information Entropy (EN). PSNR is the ratio between the maximum signal power and the signal noise power. It is used to measure the quality of the reconstructed image. The larger the PSNR value is, the better the output image quality. Since PSNR is a reference-based metric, for input image $I_1, I_2$ and fused image $I_o$, we first compute the $PSNR(I_1, I_o)$ and $PSNR(I_2, I_o)$, and then the average of these two numbers to obtain the final PSNR value. SSIM is a metric to measure the structural similarity of two images from brightness, contrast, and structure. The SSIM is calculated in the same way as the PSNR. FSIM is a full reference fusion metric that measures how close the phase congruency and magnitude of the gradient are between the fused image and source image(s). MI measures the similarity of the image intensity distributions between two or more images, which is calculated in the same way like the PSNR as well. Pixel-wise feature MI is a non-reference image fusion metric that measures the mutual information of the image features between fused image and source image(s). Finally, information entropy measures the amount of information in the fused image. For all metrics, the larger the value, the better the fusion quality.

\subsection{Selection of fusion strategy}

As we discussed in Section \ref{fs} and Equation \eqref{sfnn}, our candidate functions for $\phi(\cdot)$ are $max()$, $mean()$, and $sum()$ functions to determine the weights from the matrix nuclear norm. We also investigate the squared version of the three functions above, e.g., $\phi(\cdot) = max()^2$. Next, we provide detailed quantitative results when different fusion strategies are used in Table \ref{dif_fs}. FER \cite{fu2021multiscale} and FL1N \cite{li2022multiscale} are two fusion strategies proposed previously for the same task. The SFNN is what we proposed in this work, and selection methods are described above. Our proposed fusion strategy performed well on four metrics (PSNR, FMI, FSIM, EN). Figure \ref{fs_qualres} shows the qualitative results of different fusion strategies. For SFNN, we select two strategies that produce the best image quality visually, and also notice that different choices of $\phi$ in our proposed strategy do not affect the metrics very much. Although our SFNN-$max^2$ achieves the best result in terms of the objective metrics, it has a low visual fidelity, hence we select the second optimal SFNN-$max$ as our ultimate strategy in this work. Compared with the FER strategy \cite{fu2021multiscale} in Figure \ref{fig:fer}, our results have better fidelity, and the inner tissue boundary is more clear. Similarly, our results produce a more brightness edge than FL1N strategy \cite{li2022multiscale} in Figure \ref{fig:fl1n}, which is better when separating between the edge and tissues.

\begin{table}[H]
\renewcommand{\arraystretch}{1.3}
\caption{Comparison between different fusion strategies, bold values represent the best results, and values in blue represent the second-best results.}
\centering
\resizebox{\columnwidth}{!}{
\begin{tabular}{ccccccc}
\hline & PSNR & SSIM & MI & FMI-Pixel & FSIM & Entropy \\
\hline FER \cite{fu2021multiscale} & $13.944$ & $0.736$ & $4.551$ & $0.876$ & $0.806$ & $8.720$ \\
FL1N \cite{li2022multiscale} & $15.979$ & $0.739$ & $4.569$ & $0.878$ & $0.813$ & $9.782$ \\
SFNN ($mean$) & $15.876$ & ${\color{blue}0.740}$ & $\textbf{4.578}$ & $0.876$ & $0.812$ & $9.772$ \\
SFNN ($mean^2$) & $16.311$ & $\textbf{0.741}$ & ${\color{blue}4.567}$ & ${\color{blue}0.883}$ & $0.819$ & $9.812$ \\
SFNN ($max$) & ${\color{blue}16.413}$ & ${\color{blue}0.740}$ & $4.558$ & $\textbf{0.891}$ & ${\color{blue}0.820}$ & ${\color{blue}9.816}$ \\
SFNN ($max^2$) & $\textbf{18.258}$ & ${\color{blue}0.740}$ & $4.451$ & $\textbf{0.891}$ & $\textbf{0.829}$ & $\textbf{9.887}$ \\
SFNN ($sum$) & $15.876$ & ${\color{blue}0.740}$ & $\textbf{4.578}$ & $0.876$ & $0.812$ & $9.772$ \\
SFNN ($sum^2$) & $16.311$ & ${\color{blue}0.740}$ & $\textbf{4.578}$ & $0.876$ & $0.819$ & $9.812$ \\
\hline
\end{tabular}
}
\label{dif_fs}
\end{table}

\begin{figure}[h!]
     \centering
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge1.png}
         \caption{CT}
         \label{fig:ct_orig}
     \end{subfigure}
     \hfill
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge2.png}
         \caption{MRI}
         \label{fig:mri_orig}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge4.png}
         \caption{FER \cite{fu2021multiscale}}
         \label{fig:fer}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge5.png}
         \caption{FL1N \cite{li2022multiscale}}
         \label{fig:fl1n}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge6.png}
         \caption{SFNN-Max}
         \label{fig:sfnnMa}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/dif_fs/enlarge7.png}
         \caption{SFNN-Mean}
         \label{fig:sfnnMe}
     \end{subfigure}
     \caption{Qualitative comparison of different fusion strategies between CT and MRI. (a) and (b) are source images, (c) is the FER fusion strategy \cite{fu2021multiscale}, (d) is the FL1N fusion strategy \cite{li2022multiscale}, and (e), (f) are the strategies we proposed in this work.}
     \label{fs_qualres}
\end{figure}

\subsection{Comparison with other methods}

Table \ref{main_res} presents the quantitative results of the MRI-CT fusion task between our proposed method and other baseline methods. Our method performs well on multiple fusion metrics: PSNR, pixel-wise FMI, FSIM, and Entropy. The largest PSNR value indicates our fused image contains less noise and achieves the best image quality. The largest FMI value demonstrates our fused image contains as many source image features as possible. The largest FSIM value suggests our fused image has less information loss at the feature level. Finally, the largest entropy shows our fused image contains more information and details.

\begin{table}[h]
\renewcommand{\arraystretch}{1.3}
\caption{Comparison between different methods, bold numbers represent optimal values}
\centering
\resizebox{\columnwidth}{!}{
\begin{tabular}{ccccccc}
\hline & PSNR & SSIM & MI & FMI-Pixel & FSIM & Entropy \\
\hline Zero-shot \cite{lahoud2019zero} & $13.525$ & $0.681$ & $4.633$ & $0.836$ & $0.738$ & $4.279$ \\
MSDRA \cite{li2022multiscale} & $15.693$ & $0.697$ & $4.586$ & $0.867$ & $0.797$ & $8.167$ \\
MSRPAN \cite{fu2021multiscale} & $14.528$ & $\textbf{0.741}$ & $\textbf{4.652}$ & $0.874$ & $0.808$ & $8.969$ \\
Ours & $\textbf{16.413}$ & $0.740$ & $4.558$ & $\textbf{0.891}$ & $\textbf{0.820}$ & $\textbf{9.816}$ \\
\hline
\end{tabular}
}
\label{main_res}
\end{table}

From Figure \ref{ca_qualres}, it is also obvious from the enlarged red box that the proposed method preserves the edge and detail information well. Compared to other methods, Figure \ref{fig:ca_zs} does not retain the edge information from CT, and details of inner tissues are also distorted; Figure \ref{fig:ca_fer} has a favorable visual appearance, but the fidelity is low and lost the tissue boundary information from the MRI source image. Figure \ref{fig:ca_fl1n} has an uncleared boundary so it is hard to differentiate between boundaries and tissues. Finally, our proposed method has a better intensity contrast between edges and tissues, retains important information from both source images, and results in better fidelity.

\begin{figure}[h!]
     \centering
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_1.png}
         \caption{CT}
         \label{fig:ct_ca_orig}
     \end{subfigure}
     \hfill
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_2.png}
         \caption{MRI}
         \label{fig:mri_ca_rig}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_3.png}
         \caption{Zero-learning \cite{lahoud2019zero}}
         \label{fig:ca_zs}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_4.png}
         \caption{MSPRAN \cite{fu2021multiscale}}
         \label{fig:ca_fer}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_5.png}
         \caption{MSDRA \cite{li2022multiscale}}
         \label{fig:ca_fl1n}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/comp_algo/coalg_6.png}
         \caption{Ours}
         \label{fig:ca_our}
     \end{subfigure}
     \caption{Qualitative comparison of different methods between CT and MRI. (a) and (b) are source images, (c) is the zero-learning method \cite{lahoud2019zero}, (d) is the MSPRAN method \cite{fu2021multiscale}, (e) is the MSDRA method \cite{li2022multiscale}, and (f) is the method we proposed in this work.}
     \label{ca_qualres}
\end{figure}

\subsection{Ablation study}

To test the effectiveness of our proposed loss function, we ablate the image gradient loss and the perceptual loss as described in Equation \eqref{grad_img} and \eqref{perp_loss}. Hence, the loss function becomes only the MSE loss in Equation \eqref{mse}. Table \ref{abl} shows the quantitative results for the two loss functions on six different fusion metrics.

\begin{table}[h]
\renewcommand{\arraystretch}{1.3}
\caption{Ablation study on the loss function, bold numbers represent optimal values.}
\centering
\resizebox{\columnwidth}{!}{
\begin{tabular}{ccccccc}
\hline & PSNR & SSIM & MI & FMI-Pixel & FSIM & Entropy \\
\hline Only MSE & $\textbf{16.519}$ & $0.739$ & $4.428$ & $0.890$ & $\textbf{0.820}$ & $8.72$ \\
All loss & $16.413$ & $\textbf{0.740}$ & $\textbf{4.558}$ & $\textbf{0.891}$ & $\textbf{0.820}$ & $\textbf{9.816}$ \\
\hline
\end{tabular}
}
\label{abl}
\end{table}

We observe that using the ablated loss function (MSE loss) only achieves a higher PSNR value by a small margin. The rest of the metrics except FSIM, are all worse than using all three losses (Equation \eqref{total_loss}). FSIM values are equal across both loss functions. It is evident that gradient and perceptual loss will force the network to learn the high-level semantic features and focus more details on the features that would potentially contain important information, rather than purely minimizing pixel intensity differences using the traditional MSE loss.

\section{Discussion and Limitations}

Through the comparison of different fusion results visually, we can see that the fused image using our proposed framework have a clearer representation of boundaries and details of inner tissues. The image produced by our framework is more natural than other baseline methods, where other images are either a little bit sharp, contain artificial noise, or the boundary information is not clear enough. The quantitative results validate our findings by performing well on several fusion metrics.

We do not test our proposed framework on the other two fusion tasks (MRI-SPECT and MRI-PET). The SPECT and PET images are 3-channel and MRI images are 1-channel gray-scale. To achieve the fusion task, we could transfer the 3-channel images into YCbCr color space and only focus on the Y channel for structural details and brightness information. The Y channel and the single channel MRI can be used as inputs to our framework. Then, the output can be interpreted as the fused Y-channel. We combine it with Cb and Cr channels and convert it back to RGB to get the final fused 3-channel image. We expect our proposed method also performs well on the 3-channel image fusion task.

Our proposed fusion strategy does not perform well in some cases (see Figure \ref{failed_qualres}) and may introduce fusion noise to the result. In future work, we will focus on reducing the noise while maintaining the overall fusion quality.

\begin{figure}[h!]
     \centering
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/fail/fail_1.png}
         \caption{CT}
         \label{fig:ct_forig}
     \end{subfigure}
     \hfill
          \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/fail/fail_2.png}
         \caption{MRI}
         \label{fig:mri_forig}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/fail/fail_4.png}
         \caption{MSDRA \cite{li2022multiscale}}
         \label{fig:f_msdra}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figures/fail/fail_3.png}
         \caption{Ours}
         \label{fig:fours}
     \end{subfigure}
     \caption{Failure case of our proposed method, the features from CT do not handle properly and the details from MRI can barely be distinguished. MSDRA \cite{li2022multiscale} produces a slightly better visual appearance but it suffers from the artificial noise in the image.}
     \label{failed_qualres}
\end{figure}

\section{Conclusion}

In this paper, we propose a novel network architecture based on multi-scale feature extraction and fusion strategy for multimodal medical image fusion. The original multimodal images are passed into the feature extractor to multi-scale deep semantic features. The feature maps obtained from the feature extractor are fused based on the fixed Softmax-based weighted feature fusion strategy we proposed in this work. The fused feature map is used as the input to the reconstruction module to obtain the final fused image. Our fusion strategy is fixed, which means that there is no parameter that needs to be updated in both the training and inference phase. Thus, it can achieve real-time image fusion. Extensive experiments show our proposed method is superior to several baseline methods in both subjective visual appearance and objective fusion metrics.

Our network architecture is based on the MSRPAN \cite{fu2021multiscale}, which improves on their pyramid structure. Experiments also validate the effectiveness of our method compared to other reference methods. However, our method performs not well on some corner cases and might introduce artificial noise in the fused image as we discussed previously. Despite the limitations, our method could provide a reliable reference for disease diagnosis in the real-life clinical routine.

\ifpeerreview \else
\section*{Acknowledgments}

The authors would like to thank Dr. David B. Lindell for his advice in this work and Jingtao Zhou for his suggestions on the figures.

\fi

\bibliographystyle{IEEEtran}


\section{Introduction} \label{section1}
As benchmark systems to accurately measure the masses and radii for a variety of stars with negligible model dependence,
binary stars play a fundamental role in our understanding of stars and even the Universe \citep{Torres2010, Chen2020, Lampens2021}.
Eclipsing binary (EB) systems hosting at least one pulsating star are much more valuable,
since they will provide strong constraints on the input physics of asteroseismic models
for the pulsating component and also offer the possibility to carry out mode identification
through the technique of eclipse mapping \citep{Nather1974, Mkrtichian2018, Chen2022}.
Meanwhile, the study of stellar oscillations will unravel the structure and dynamics of stellar interiors by means of asteroseismology \citep{Aerts2021},
which in turn helps us to probe the influences of tidal forces, mass transfer, and angular momentum transfer between the components \citep{Murphy2018, Guo2021, Kovalev2022}.

Owing to the remarkable success of the space telescopes such as $\it Kepler$ \citep{Borucki2010} and $\it TESS$ \citep{Ricker2015},
almost all types of pulsating stars across the whole Hertzsprung-Russell Diagram,
ranging from extremely short-period pulsating white dwarfs to massive $\beta$-Cep pulsators to long-period red giant variable stars,
have been discovered in EB systems.
\cite{Gaulme2019} performed a systematic search for pulsating stars in the $\it Kepler$ EB catalog\footnote{\url{http://keplerebs.villanova.edu/}}
and identified 303 systems with stellar pulsators.
Combining the high-precision space-based photometry and ground-based spectroscopic observations,
researchers have made significant advances in the study of pulsating EB.
A total of 14 double-lined spectroscopic EB with an oscillating red giant component have been found by the NASA $\it Kepler$ mission \citep{Benbakoura2021}.
After comparing the masses and radii of the red giants given by dynamic modelling of EB systems with the results from the asteroseismic scaling relations,
\cite{Gaulme2016} and \cite{Benbakoura2021} concluded that the asteroseismic radii and masses are systematically overpredicted by 5\% and 15\%, respectively.
A number of detailed asteroseismic studies on individual EB systems with classical pulsators (e.g., $\delta$ Scuti and $\gamma$ Doradus) have been published
(e.g. \citealt{Hambleton2013, Schmid2016, Guo2017, Zhang2018, Chen2021}).
\cite{Li2020} performed an ensemble asteroseismic study of 35 $\it Kepler$ $\gamma$ Dor stars in EB systems
and discovered that $\gamma$ Dor stars in binaries tend to show slower near-core rotation rates compared with that of single stars.
Space missions have also contributed to the understanding of EB with compact pulsators,
including pulsating extremely low mass white dwarfs (e.g. \citealt{Zhang2016, Wang2020, Kim2021}), hot subdwarfs (e.g. \citealt{Baran2018, Baran2021, Dai2022}),
and canonical white dwarfs (e.g. \citealt{Parsons2017, Garza2021}).
Another highlight from $\it Kepler$ and $\it TESS$ is the discovery and characterisation of heartbeat stars, triggering the beginnings of tidal asteroseismology
(e.g. \citealt{Hambleton2013, Bowman2019, Handler2020}).
Thus far, about 180 heartbeat stars have been identified in the entire $\it Kepler$ data \citep{Kirk2016}
and  \cite{Prsa2021} discovered 22 heartbeat binaries in the 2 minute TESS data up to Sector 26, but a handful of them have been closely studied.

In this work, we use high-precision $\it TESS$ photometry and high-resolution spectroscopic measurements
 to characterize a neglected eccentric binary system FX UMa (TIC 219707463).
 In Section \ref{section2}, we describe the observations, radial velocity extraction,
 spectroscopic orbital elements,
 spectral disentangling, and determination of the atmospheric parameters.
The final orbital solution and physical parameters for FX UMa are given in Section \ref{section3}.
 In Section \ref{section4}, we perform a Fourier analysis of the residual light curve and discuss the pulsations.
 Finally, we summarize our findings in Section \ref{section5}.

\section{Observations and data reductions} \label{section2}
\subsection{TESS photometry}
Until now, FX UMa was observed by $\it TESS$ in a 2-min cadence mode during three non-contiguous sectors: 14, 20, and 40.
The $\it TESS$ light-curve (LC) files,  produced by the $\it TESS$ Science Processing Operations Center \citep{Jenkins2016},
were downloaded from the Mikulski Archive for Space Telescopes (MAST\footnote{\url{https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html}}).
For the EB system FX UMa, we made use of the Simple Aperture Photometry data labelled SAP\underline{~~}FLUX,
since PDCSAP\underline{~~}FLUX (Pre-Search Data Conditioning; \citealt{Smith2012}) does not detrend the LC of eclipsing binaries perfectly.
Visual examination of the raw LC reveals negligible signatures of unphysical trends.
So we didn't do anything except discard five outlying data points.
At the end, a Python package called $\it Lightkurve$\footnote{\url{https://docs.lightkurve.org/}} \citep{Lightkurve2018} was applied  to remove all nan values
and stitch together multiple sectors of $\it TESS$ observations after normalizing them.
The final $\it TESS$ light curve of FX UMa  is plotted in Figure \ref{fig1},
which shows a periodic brightening near the eclipse, the typical feature of heartbeat stars.
The zoom panel presents a short segment of the LC, where the oscillations can be clearly seen.

\begin{figure}
\center
\includegraphics[scale=0.7]{fig1.pdf}
\caption{The full 2-min $\it TESS$ light curve of FX UMa from sectors 14, 20, and 40.
The zoom panel is a closer view of a selected region of the light curve.}
\label{fig1}
\end{figure}

\subsection{SONG high-resolution spectra}
High-resolution spectroscopic observations of FX UMa were performed by using the 1m automated Hertzsprung $\it SONG$ telescope
at the Teide Observatory on the island of Tenerife, Spain \citep{Andersen2014, Fredslund2019}.
The $\it SONG$ spectrograph consist of 51 spectral orders in the wavelength range of 4400 -- 6900 \AA.
The observations were carried out with slit no. 5, corresponding a spectral resolution of 77000.
We obtained a total of 34 spectra from 2019 October 28 to 2020 April 29,
one of which was discarded in this study due to low signal-to-noise ratio (S/N).
The 1D spectra, including the blaze function from the master (summed) flat field, were extracted with
the $\it SONG$ spectral-reduction pipeline called $\it songwriter$ \citep{Ritter2014, Grundahl2017}.
The detailed description of the Hertzsprung $\it SONG$ telescope is well-documented in https://soda.phys.au.dk/.

\subsubsection{Radial velocity measurements}
The broadening-function technique\footnote{\url{http://www.astro.utoronto.ca/\~rucinski/SVDcookbook.html}}
(BF; \citealt{Rucinski1992,Rucinski2002}) was used to
determine radial velocities (RVs) from the 1D extracted spectra of FX UMa.
Compared with more familiar cross-correlation function,
 the BF method generally improves the ability to measure the Doppler shifts
 from the complex spectra of double-lined spectroscopic binaries
 showing substantial rotational broadening and overlapping spectral lines \citep{Rucinski2002, Bavarsad2016, Clark2019, chen101975}.
 This is essential for a short-period binary with at least one fast-rotating component star like FX UMa.

 We first removed the spectrograph blaze function from the 1D spectra and normalized each spectrum by its continuum order by order
using the open-source spectroscopic tool iSpec \citep{Blanco2014, Blanco2019}.
The orders of each spectrum were then merged,
since single spectral order covers a narrow spectral range of about 4 nm,
which makes it hard to give strong BF peaks.
Afterwards, we computed the BFs of our target spectra by employing a modified version of the BF software suite provided publicly
by \cite{Rawls2016}\footnote{\url{https://github.com/mrawls/BF-rvplotter}}.
A high-resolution PHOENIX synthetic spectrum \citep{Husser2013} was selected as the BF template.
The portion of spectra we computed the BFs are in the wavelength region 4900 -- 5500 \AA.
This is because this region not only staves off strong hydrogen Balmer lines,
but also contains most of the information on the velocities.
Following the same procedure of \cite{Rawls2016}, we adopted gauss filter to process original BF smoothly to eliminate un-correlated, small-scale noise.
The normalized smoothed BF for FX UMa is shown in Figure \ref{fig2}.
We can see clearly from it the BF in velocity space displays two peaks,
whose positions are equal to the RVs of both components of FX UMa.
The geocentric (uncorrected) RVs, marked by the blue vertical lines in Figure \ref{fig2},
were obtained by rotational profile fitting to the smoothed BFs.
The barycentric velocity corrections provided by the SONG pipeline were then applied to them
to yield the final RV measurements, which are given in Table \ref{tab1}.
The uncertainty of our measurements comes from the error in fitting a rotational profile to each BF profile
with an open-source software package called Non-Linear Least-Squares Minimization and Curve-Fitting for Python (LMFIT)\footnote{\url{https://lmfit.github.io/lmfit-py/}}.

\subsubsection{Spectroscopic orbital elements}
In order to solve for the spectroscopic orbital parameters of FX UMa,
we used the $\it rvfit$ code \citep{Iglesias2015} to fit our RV measurements.
Using an adaptive simulated annealing algorithm,
the $\it rvfit$ code can automatically fit the RVs of stellar binaries and exoplanets.
It is also a user-friendly code that converges to a global solution minimum without the need to provide preliminary parameter values (see \citealt{Iglesias2015}, for details).
In the analysis of FX UMa, we fixed the orbital period to 4.50725 days,
which was taken from the Data Validation Report Summary, provided by the $\it TESS$ Science Processing Operations Center Pipeline.
The other six orbital parameters, i.e., epoch of periastron passage ($T_{p}$),
argument of the periastron ($\omega$), eccentricity ($e$),
systematic velocity ($v_{\gamma}$), and semi-amplitudes of RVs for both components ($K_{1}$, $K_{2}$), were kept free during the analysis.
The initial values for $v_{\gamma}$, $K_{1}$ and $K_{2}$ were estimated by visually examining the phased RV curves of FX UMa.
The initial values of the remaining three parameters ($T_{p}$, $\omega$, $e$) were set to 2458684.2139, 0, and 0, respectively.
To cover reasonable models, we allowed the six adjustable parameters to vary in wide ranges.
The fitted parameters and other derived quantities for the best-fit model are given in Table \ref{tab2}.
Figure \ref{fig3} displays the theoretical RV curve fits to the observed ones and the differences between them.

\begin{figure}
\center
\includegraphics[scale=0.48]{fig2.pdf}
\caption{BF plots for FX UMa. Each panel stands for one spectroscopic observation,
for which the double-rotational profile fit to the normalized smoothed BF are plotted.
The central position of each BF, corresponding to the radial velocity of each component of FX UMa, is marked by the blue vertical line.}
\label{fig2}
\end{figure}

\begin{deluxetable}{lrrr}\label{tab1}
\tablecolumns{4}
\tablewidth{0pc}
\tabletypesize{\scriptsize}
\tablecaption{Radial velocities for FX UMa extracted from SONG spectra}
\tablehead{\colhead{BJD}    &\colhead{Phase}		&\colhead{$RV_{1}$}             &\colhead{$RV_{2}$}              \\
\colhead{(2400000+)}      		&\colhead{}		&\colhead{(km s$^{-1}$)}       &\colhead{(km s$^{-1}$)}        }
\startdata
  58784.616338  &0.276     &-70.94$\pm$0.19     	&53.47$\pm$0.18\\
  58784.704294  &0.295     &-68.58$\pm$0.24      	&50.57$\pm$0.23\\
  58785.581477  &0.490     &-39.14$\pm$0.24     	&19.96$\pm$0.27\\
  58787.587152  &0.935     &140.26$\pm$0.20     	&-157.86$\pm$0.19\\
  58787.760165  &0.973     &99.67$\pm$0.34       	&-121.47$\pm$0.40\\
  58804.530779  &0.694     &13.61$\pm$0.32    	&-31.63$\pm$0.30\\
  58804.694505  &0.730     &26.77$\pm$0.21    	&-45.09$\pm$0.22\\
  58805.731430  &0.960     &119.89$\pm$0.17   	&-138.42$\pm$0.16\\
  58806.746889  &0.186     &-80.72$\pm$0.17     	&60.49$\pm$0.15\\
  58811.552382  &0.252     &-74.31$\pm$0.15     	&55.93$\pm$0.16\\
  58813.720354  &0.733     &26.80$\pm$0.18    	&-46.66$\pm$0.19\\
  58814.512091  &0.909     &130.09$\pm$0.16   	&-148.96$\pm$0.17\\
  58814.780540  &0.968     &105.00$\pm$0.15   	&-125.25$\pm$0.14\\
  58817.532391  &0.579     &-18.67$\pm$0.31      	&1.66$\pm$0.40\\
  58818.541681  &0.803     &57.62$\pm$0.14    	&-77.23$\pm$0.16\\
  58819.605796  &0.039     &-39.30$\pm$0.27     	&16.94$\pm$0.27\\
  58821.678834  &0.499     &-35.34$\pm$0.19     	&18.55$\pm$0.21\\
  58827.511274  &0.793     &53.70$\pm$0.15    		&-71.89$\pm$0.16\\
  58828.575788  &0.029     &-28.72$\pm$0.28      	&7.28$\pm$0.26\\
  58829.736020  &0.286     &-70.76$\pm$0.21     	&51.66$\pm$0.22\\
  58830.536858  &0.464     &-42.61$\pm$0.19     	&23.21$\pm$0.20\\
  58910.612349  &0.230     &-76.92$\pm$0.15     	&57.71$\pm$0.15\\
  58911.533949  &0.434     &-48.47$\pm$0.18     	&28.70$\pm$0.19\\
  58912.439972  &0.635     &-3.72$\pm$0.34    		&-16.09$\pm$0.33\\
  58914.423942  &0.075     &-66.84$\pm$0.19     	&49.30$\pm$0.19\\
  58915.544257  &0.324     &-63.94$\pm$0.19     	&47.22$\pm$0.20\\
  58916.524161  &0.541     &-28.26$\pm$0.23      	&8.56$\pm$0.25\\
  58929.502502  &0.421     &-50.92$\pm$0.21     	&30.80$\pm$0.21\\
  58956.422493  &0.393     &-55.58$\pm$0.27     	&35.18$\pm$0.27\\
  58960.427192  &0.282     &-71.02$\pm$0.17     	&51.35$\pm$0.16\\
  58962.465837  &0.734     &27.77$\pm$0.24    	&-46.34$\pm$0.22\\
  58965.435556  &0.393     &-55.04$\pm$0.22     	&35.32$\pm$0.21\\
  58969.422138  &0.278     &-70.00$\pm$0.21     	&53.36$\pm$0.19\\
\enddata
\end{deluxetable}

\begin{deluxetable}{lr}\label{tab2}
\tablecolumns{2}
\tabletypesize{\scriptsize}
\tablecaption{Radial-velocity fitting results for the SB2 system FX UMa}
\tablehead{\colhead{Parameter}   		&\colhead{Value} }
\startdata
\multicolumn{2}{c}{Adjusted Quantities}\\
\hline
$P_{orb}$ (d)			&4.50725$^{\star}$\\
$T_p$ (BJD)		&2458688.59179 $\pm$ 0.00034\\
$e$			&0.54765 $\pm$ 0.00026\\
$\omega$ (deg)		&50.269 $\pm$ 0.049\\
$v_{\gamma}$ (km/s)	&-9.475 $\pm$ 0.024\\
$K_1$ (km/s)		&110.197 $\pm$ 0.057\\
$K_2$ (km/s)		&110.251 $\pm$ 0.056\\
\hline
\multicolumn{2}{c}{Derived Quantities}\\
\hline
$M_1\sin ^3i$ ($M_\odot$)	&1.4656 $\pm$ 0.0019\\
$M_2\sin ^3i$ ($M_\odot$)	&1.4649 $\pm$ 0.0019\\
$q = M_2/M_1$		&0.99951 $\pm$ 0.00072\\
$a_1\sin i$ ($10^6$ km)	&5.7146 $\pm$ 0.0032\\
$a_2\sin i$ ($10^6$ km)	&5.7175 $\pm$ 0.0031\\
$a  \sin i$ ($10^6$ km)	&11.4321 $\pm$ 0.0045\\
\hline
\multicolumn{2}{c}{Other Quantities}\\
\hline
$\chi^2$			&1522.83\\
$N_{obs}$ (star 1)		&33\\
$N_{obs}$ (star 2)	&33\\
Time span (days)		&184.81\\
$rms_1$ (km/s)		&1.22\\
$rms_2$ (km/s)		&1.63\\
\enddata
\tablecomments{\footnotesize $^{\star}$ Parameter fixed beforehand.}
\end{deluxetable}

\begin{figure}
\center
\includegraphics[scale=0.7]{fig3.pdf}
\caption{RVs of both component stars of FX UMa as a function of phase.
The blue solid lines represent the theoretical RV curves for the both components,
derived by using the $\it rvfit$ code \citep{Iglesias2015}.
The red dot-dashed lines represent the best-fitting models from the PHOEBE (Model 1) that were constrained by both LC and RV observations,
which will be described in Section \ref{section3}.
The yellow dashed line in the top panel shows the systemic velocity of $v_{\gamma}$ = -9.475 km/s.
The residuals from the best-fit model are presented in the bottom two panels.}
\label{fig3}
\end{figure}

\subsubsection{Atmospheric parameters from disentangled spectra}
As a double-lined spectroscopic binary system (SB2), each observed spectrum of FX UMa is a composite of individual spectra of component stars.
We attempted to reconstruct the spectra of the individual components
using the spectral disentangling technique (for a summary of different methods see e.g., \citealt{Pavlovski2010}).
The spectral disentangling tool FDBinary\footnote{\url{http://sail.zpf.fer.hr/fdbinary/}} \citep{Ilijic2004}
was employed in this study for performing spectral decomposition.
Without the use of template spectra, FDBinary requires six orbital parameters ($P$, $T_{p}$, $\omega$, $e$, $K_{1}$, $K_{2}$)
to outline the shape of the Keplerian RV curve.
In the runs, we focused on the spectral interval of 4900 -- 5500 \AA\ and the orbital period ($P$) was kept at 4.50725 days.
The initial values for $T_{p}$, $\omega$, $e$, $K_{1}$, and $K_{2}$ were taken from the previous RV fitting results
and we let them free in the analysis.
Based on the amplitude of each star's BF in Figure \ref{fig2}, the light contributions of the two component stars of FX UMa
were estimated to be 0.49 and 0.51, respectively.
One of 33 $\it SONG$ spectra that were used in RV analysis was discarded due to a relatively low S/N.
The remaining 32 $\it SONG$ spectra of FX UMa were processed together in FDBinary and
at last we obtained the disentangled spectra of each component star, as shown in Figure \ref{fig4}.

The atmospheric parameters for both component stars were then obtained from their disentangled spectra.
All this was done using the synthetic spectra fitting technique implemented in the code iSpec \citep{Blanco2014, Blanco2019}.
iSpec integrates  a broad variety of radiative transfer codes, model atmospheres, solar abundances, and atomic line lists.
We fit the disentangled spectrum of each star separately.
In this analysis, we employed the SPECTRUM radiative code \citep{Gray1994},
the ATLAS9 Castelli model atmospheres \citep{Kurucz2005},
the Grevesse 2007 solar abundances \citep{Grevesse2007},
and the third version of the Vienna Atomic Line Database (VALD3; \citealt{Ryabchikova2015})
to produce synthetic spectra on the fly.
Since the surface gravity log $g$ is usually not well constrained with spectroscopy, we
fixed log $g$ to the values derived from the stellar mass and radius (log $g_{1}$ =  log $g_{2}$ = 4.3, see the section \ref{section3}).
We iterated between the spectroscopic analysis and binary models to obtain log g.
We adopted a resolution of $R$ = 77000 appropriate to the spectral resolution of $\it SONG$ spectra.
The radial velocity for each component was set to 0, since FDBinary provides a pair of disentangled spectra with zero RV.
In this work, we didn't analyze abundances of specific elements.
Following \cite{Blanco2014}, the limb-darkening coefficient was fixed to 0.6.
The adjustable parameters consist of effective temperature $T_{eff}$, metallicity $[M/H]$,
alpha enhancement $[\alpha/Fe]$, microturbulence velocity $v_{mic}$,
and projected rotational velocity $v_{rot}$sin$i$.
The application can automatically calculate the macroturbulence velocity
based on an empirical relation established by Gaia-ESO Survey.
The resulting atmospheric parameters for both stars in FX UMa are presented in Table \ref{tab3}.
Figure \ref{fig4} displays the observed composite, disentangled, and fitted synthetic spectra
of both components of FX UMa.

\begin{figure}
\center
\includegraphics[scale=0.55]{fig4.pdf}
\caption{One observed SONG spectrum (middle, BJD = 2458814.512091, $\phi$ = 0.909) and
the disentangled spectra from FDBinary for the two components in FX UMa (upper and lower).
The zoom panel shows a small portion of these spectra to better see the details.
The besting-fitting theoretical spectra for the two stars, obtained by the iSpec code, are shown as the red solid lines.
The spectra of the stars 1 and 2 have been shifted vertically by +0.3 and -0.3, respectively, for comparison purposes.}
\label{fig4}
\end{figure}

\section{Binary modelling} \label{section3}
As seen in Figure \ref{fig1}, the LCs of FX UMa show a periodic, broad brightening near the eclipse,
which look distinctly different from intrinsic star variability and instead are characteristic of heartbeat stars.
Apart from the eclipse-like light changes and the ``heartbeat-like" profile, the LCs outside of eclipses present
multi-periodic light variations with characteristic of hybrid $\delta$ Sct - $\gamma$ Dor oscillations.
To obtain the physical parameters of FX UMa and probe its pulsational properties in detail,
we carried out a simultaneous fit to our double-lined RVs and $\it TESS$ LCs
using the PHysics Of Eclipsing BinariEs (PHOEBE, \citealt{Andrej2018, Conroy2020})
code\footnote{\url{http://phoebe-project.org/}} in the detached mode.

\subsection{PHOEBE setup}
Our analyses were separately carried out in two cases: LCs and RVs were simultaneously solved with basic parameters (hereafter Model 1), and
were analyzed by including several nuisance parameters in basic parameters (hereafter Model 2).
First, we fixed the gravity-brightening coefficient and the bolometric albedo to their typical values
suitable for stars with radiative envelopes, i.e. $A_{1,2}$ = 1 and $g_{1,2}$ = 1.
The limb-darkening coefficients for each star were automatically interpolated from the PHOEBE's built-in lookup table,
with the atmosphere models set as PHOENIX.
We assumed a pseudo-synchronous rotation for the eccentric binary system FX UMa.
The synchronicity parameter was computed by using the formalism of $F_{1, 2}$ = $\sqrt{(1~+~e)/(1~ - ~e)^3}$ \citep{Andrej2018},
where $e$ is the orbital eccentricity from the above RV fit.
The initial values for the mass ratio ($q$ = $M_{2} / M_{1}$), orbital eccentricity ($e$), argument of periastron ($\omega$),
projected orbital semi-major axis ($a_{orb}$ sin$i$), and systemic velocity ($v_{\gamma}$), were taken from Table \ref{tab2}.
The input effective temperatures for both components ($T_{\rm eff, 1}$,  $T_{\rm eff, 2}$) were taken from our spectral analysis.
Since iSpec estimates the errors based on the co-variance matrix \citep{Blanco2014} and
in Table \ref{tab3} the uncertainties of the effective temperatures were probably underestimated,
we assumed the errors of $T_{\rm eff, 1}$ and $T_{\rm eff, 2}$ to be 100 K.
The orbital ephemeris, i.e. reference time of superior conjunction ($t_{0}$) and orbital period ($P_{orb}$),
were obtained by the TESS light curves of FX UMa.
We utilized the emcee sampler \citep{Foreman2019, Foreman2013} built into PHOEBE to
explore the parameter spaces, find the optimal solution, and determine the uncertainties.
In a word, Model 1 has a total of 13 parameters that were fitted:
$t_{0}$, $P_{orb}$, $q$, $e$, $\omega$, $a_{orb}$ sin$i$, $v_{\gamma}$, $T_{\rm eff, 1}$,  the temperature ratio $T_{\rm eff, 2}$ / $T_{\rm eff, 1}$,
orbital inclination ($i$), the sum and ratio of fractional radii (($R_{\rm equiv, 1}$ + $R_{\rm equiv, 2}$) / $a_{orb}$, $R_{\rm equiv, 2}$ / $R_{\rm equiv, 1}$), and the passband luminosity of star 1 ($L_{pb}$).
Besides, we also marginalized over three noise nuisance parameters ($\sigma$$_{lnf,lc01}$, $\sigma$$_{lnf,rv01}$, $\sigma$$_{lnf, rv02}$),
since the observational uncertainties may be underestimated by a constant factor (see section 3.5 of the work of \cite{Conroy2020}).
Based on our spectroscopic results and light curve properties of FX UMa,
we initialized chains with Uniform or Gaussian distributions on each of the 13 parameters.
We then ran the emcee sampler using 50 walkers for 8000 iterations.
Based on the autocorrelation times of the chains, PHOEBE automatically determines defaults for burn-in and thinning.
The ``burn-in" was set to 1888 (twice the default value), where most of the walkers have converged in examining the ln(probability) plot,
and walkers with ln(probability) $<$ 3.318e5 were rejected.

\subsection{Uncertainties and stellar parameters}
The median value of each parameter's posterior distribution is adopted as the best-fitting value and given in Table \ref{tab3},
in which the upper and lower uncertainties are obtained at the 16th and 84th percentiles, respectively.
As we can see from Table \ref{tab2} and Table \ref{tab3},
there is a negligible difference in $q$, $e$, $\omega$, $a_{orb}$ sin$i$, and $v_{\gamma}$
between the independent RV fit by the code $\it rvfit$ and the PHOEBE binary model that were constrained by both LC and RV observations.
For comparison, the synthetic RV curves for the two stars, produced by the final PHOEBE binary model (Model 1), are overplotted in Figure \ref{fig3}.
As can be seen from this figure, both fits are consistent with each other.
Based on the RV residuals from the PHOEBE fit, the root mean square errors for the two components are estimated to be 1.23 and 1.40, respectively,
which are also in good agreement with the calculated scatter from the $\it rvfit$.
However, the formal PHOEBE errors are 10 times smaller than those of the RV-only fit.
If the uncertainties from PHOEBE were reasonable,
they would result in unrealistic mass errors for the two component stars ($\sim$0.0003$M_{\sun}$), as shown in Table \ref{tab3}.
In order to estimate the actual uncertainties in the model, we included $A_{1,2}$, $g_{1,2}$, $F_{1, 2}$, and third light ($l_{3}$)
as additional nuisance parameters in our fits (hereafter Model 2).
We initialized emcee with 90 walkers, each of which was computed for 40,000 steps (3.6 million models).
The ``burn-in" was set to 20,000 and walkers with ln(probability) $<$ 3.305e5 were rejected.
The fitted parameters are given in the Model 2 column of Table \ref{tab3}.

The synthetic light curves, computed with the Model 1 and Model 2 parameters of Table \ref{tab3},
are shown as solid lines in the top panel of Figure \ref{fig5}.
The unbinned and binned light residuals corresponding to the two models are presented in the middle and bottom panels, respectively.
As can be seen in Figure \ref{fig5}, there are almost no systematic trends in the LC residuals.
Overall, the two models describe the observed light curve satisfactorily.
However, the errors of the fitted parameters were unreasonably small.
This may be because FX UMa is a partially eclipsing binary with very similar components and only one sharp eclipse can be seen in its light curve.
It is difficult to constrain the inclination, relative radii, and temperature ratio for such system.
 In Appendix, we show the parameter posterior distributions and their interdependencies for the two PHOEBE models.
Based on the full corner plots, we can see that most parameters do not obey Gaussian distributions
and there are correlations between the model fitting parameters.
The formal PHOEBE uncertainties are based solely on the posterior distributions of the fitted parameters and
do not include correlations among parameters or systematic effects \citep{Conroy2020, Jones2022}.
For these reasons, we suggest that the uncertainties in the two models are seriously underestimated and should be taken just as a proof of concept.

As can be seen from the posterior distributions of the Model 2, it seems that an extreme number of iterations may yield near Gaussian-like distributions.
But when we tried to restart additional 10,000 iterations from the emcee solution of Model 2, the computer automatically killed the PHOEBE program for a short time due to insufficient memory.

\begin{deluxetable}{lrr}\label{tab3}
\tablecolumns{3}
\tablewidth{0pc}
\tablecaption{Atmospheric parameters and binary model parameters for FX UMa}
\tablehead{\multicolumn{3}{c}{iSpec analysis}}
\startdata
Parameters    &Star 1    &Star 2\\
\hline
$T_{eff}$ (K)		         	&7392$\pm$44			&7390$\pm$44	\\
log $g$ (dex)				&4.3$^{\star}$			&4.3$^{\star}$		\\
$[M/H]$ (dex)				&-0.41$\pm$0.02		&-0.43$\pm$0.02	\\
$[\alpha/Fe]$ (dex)			&0.23$\pm$0.03		&0.23$\pm$0.03	\\
$v_{mic}$ (km s$^{-1}$)		&4.34$\pm$0.12		&4.29$\pm$0.12	\\
$v_{rot}$sin$i$ (km s$^{-1}$)	&65.4$\pm$0.8			&66.8$\pm$0.8	\\
\hline
\multicolumn{3}{c}{PHOEBE analysis$^{\dagger}$}\\
\hline
Parameters    &Model 1    &Model 2\\
\hline
$t_{0}$			  	  	&1684.21629$_{-1e-5}^{+1e-5}$		 	 &1684.21601$_{-3e-5}^{+4e-5}$\\
$P_{orb}$ (days)                	&4.5072524$_{-4e-7}^{+2e-7}$				 &4.5072537$_{-5e-7}^{+4e-7}$ \\
$i$ (deg)                         		&78.59$_{-0.02}^{+0.02}$					 &78.59$_{-0.02}^{+0.05}$\\
$e$						&0.54756$_{-3e-5}^{+5e-5}$				 &0.54809$_{-7e-5}^{+6e-5}$\\
$\omega$ (deg)			&50.147$_{-0.003}^{+0.004}$				 &50.246$_{-0.030}^{+0.011}$\\
$v_{\gamma}$ (km s$^{-1}$)	&-9.485$_{-0.002}^{+0.001}$				 &-9.446$_{-0.003}^{+0.003}$\\
$q=M_{2}/M_{1}$                 	&0.9961$_{-1e-4}^{+1e-4}$				 &1.0002$_{-1e-4}^{+5e-4}$\\
$a_{orb}$ sin$i$ ($R_{\sun}$) 	&16.4385$_{-1e-3}^{+8e-4}$               		 &16.4412$_{-4e-4}^{+2e-3}$ \\
$T_{\rm eff,1}$ (K)			&7356$_{-10}^{+8}$						 &7451$_{-8}^{+10}$\\
$T_{\rm eff,2}/T_{\rm eff,1}$	&0.984$_{-0.002}^{+0.003}$				&0.981$_{-0.001}^{+0.005}$\\
($R_{\rm equiv, 1}$ + $R_{\rm equiv, 2}$) / $a_{orb}$		&0.1797$_{-1e-4}^{+2e-4}$				&0.1797$_{-2e-4}^{+1e-4}$\\
$R_{\rm equiv, 2}$ / $R_{\rm equiv, 1}$				&0.907$_{-0.005}^{+0.003}$				&0.960$_{-0.004}^{+0.004}$\\
$L_{pb}$ 					&7.09$_{-0.04}^{+0.03}$					&6.71$_{-0.05}^{+0.02}$\\
$l_{3}$					&0$^{\star}$							&0.0091$_{-7e-4}^{+6e-4}$\\
$F_{1}$					&4.09$^{\star}$							&5.68$_{-0.09}^{+0.04}$\\
$F_{2}$					&4.09$^{\star}$							&3.24$_{-0.13}^{+0.03}$\\
$A_{1}$					&1.0$^{\star}$							&0.916$_{-0.007}^{+0.006}$\\
$A_{2}$					&1.0$^{\star}$							&0.937$_{-0.005}^{+0.013}$\\
$g_{1}$					&1.0$^{\star}$							&0.919$_{-0.013}^{+0.044}$\\
$g_{2}$					&1.0$^{\star}$							&0.919$_{-0.015}^{+0.004}$\\
\hline
\multicolumn{3}{c}{Stellar parameters derived by the PHOEBE code}\\
\hline
$M_{1}$ ($M_{\sun}$) 		&1.5605$_{-3e-4}^{+4e-4}$ 				&1.5581$_{-8e-4}^{+2e-4}$\\
$M_{2}$ ($M_{\sun}$) 		&1.5543$_{-4e-4}^{+5e-4}$ 				&1.5585$_{-2e-4}^{+4e-4}$\\
$R_{1}$ ($R_{\sun}$)                &1.579$_{-0.001}^{+0.003}$               		&1.538$_{-0.003}^{+0.004}$\\
$R_{2}$ ($R_{\sun}$)                &1.433$_{-0.004}^{+0.005}$               		&1.477$_{-0.002}^{+0.003}$\\
$log$ $g_{1}$ (dex)			&4.234$_{-0.001}^{+0.001}$				&4.256$_{-0.002}^{+0.002}$\\
$log$ $g_{2}$ (dex)			&4.317$_{-0.003}^{+0.002}$				&4.291$_{-0.002}^{+0.002}$\\
\enddata
\tablecomments{\footnotesize $^{\star}$ Fixed during the fitting procedure.
$^{\dagger}$ The uncertainties are the formal ones from the PHOEBE fit, and they are probably underestimated.}
\end{deluxetable}

\begin{figure}
\center
\includegraphics[scale=0.65]{fig5.pdf}
\caption{The top panel shows the phased TESS light curve of FX UMa with the fitted models.
The red and cyan solid lines were computed with the Model 1 and Model 2 parameters of Table \ref{tab3}, respectively.
The unbinned and binned light residuals corresponding to the two models are presented in the middle and bottom panels, respectively.}
\label{fig5}
\end{figure}

\section{Pulsation characteristics} \label{section4}
The residual LC of FX UMa was obtained by removing the modeled EB light curve from the original TESS observations.
We utilized the Period04 software \citep{Lenz2005} to search for frequencies of pulsation in the residual LC of FX UMa.
This was done through an iterative pre-whitening process.
We stopped the frequency search when in the Fourier amplitude spectrum none of peaks satisfied S/N ratio $\geq$ 4 rule.
The search was restricted in the frequency range from 0 to 80 day$^{-1}$.
We have also checked for frequencies up to the Nyquist frequency of two-minute cadence TESS data ($\simeq$ 359 day$^{-1}$),
but did not find any peaks beyond 80 day$^{-1}$.
Finally, we detected a total of 103 frequencies with S/N $\geq$ 4.
With the frequency resolution of 0.002 day$^{-1}$, we performed a search for potential
orbital harmonics ($f_{i}$ = $Nf_{orb}$, $f_{orb}$ = 0.22186467$\pm$0.00000001 day$^{-1}$) and combination frequencies
using our own codes and Period04, respectively.
We identified 12 independent frequencies and 17 multiples of the orbital frequency, both of which are listed in Table \ref{tb4}.
The remaining 74 frequencies presented in Table \ref{tb5} are found as probable combination frequencies.
Figure \ref{fig6} shows the Fourier amplitude spectrum for the residual LC of FX UMa.

\subsection{Tidally Excited Modes}
More than 20\% heartbeat stars have been observed to show tidally excited oscillations (TEOs, \citealt{Kurtz2022}), driven by dynamical tides.
Most of the observed TEOs occur at harmonics of the orbital frequency ($f_{i}$ = $Nf_{orb}$),
which are likely triggered by the linear dynamical tide.
That is to say, the detected $Nf_{orb}$ peaks are thought to the signature of tidally excited modes.
At the same time we recognize that the imperfect removal of the binarity-induced light variations
can result in alias peaks of the form $Nf_{orb}$ with low amplitudes.
As seen in Table \ref{tb4}, there are a total of 17 orbital frequency harmonics
in the range of 0.6 to 28.4 day$^{-1}$ ($N$ = 3 to 128).
Thereinto, the most prominent peak is the $f_{4}$ = 3.549841$\pm$0.000003 day$^{-1}$ $\simeq$ 16$f_{orb}$
with an amplitude of 0.562 mmag and a S/N of 80.3, which can not be attributed to imperfect light curve modeling.
Following \cite{Guo2019}, we consider the 10 $Nf_{orb}$ peaks with S/N $>$ 10 as high probability TEOs.
However, we cannot rule out the possibility that other $Nf_{orb}$ peaks (4 $\leq$ S/N $\leq$ 10) are actually real TEOs.

In the lower frequency $g$-mode regime, the four significant frequencies of $f_{2}$, $f_{3}$, $f_{5}$, and $f_{9}$
are not multiples of the $f_{orb}$ and may be self-excited $\gamma$ Dor-type $g$ modes.
After checking whether there is any connection between these anharmonic frequencies and the orbital frequency,
we found that they can pair up and sum to give exact harmonics of the orbital frequency:
$f_{5}$ + $f_{9}$ $\simeq$ $f_{20}$ $\simeq$ 20$f_{orb}$  and $f_{2}$ + $f_{3}$ $\simeq$ 14$f_{orb}$.
This implies the existence of non-linear tidal processes in the eccentric binary system FX UMa.
When the amplitude of a linear TEO mode exceeds the parametric instability threshold,
then it may experience non-linear resonance mode coupling and decay into (or more) daughter modes \citep{Weinberg2012, Yu2020}.
Observationally, the sum of the daughters's frequencies is equal to the frequency of the parent mode.
Non-linear tidal oscillations have been found in some eccentric binary systems,
such as KOI-54 \citep{Burkart2012,Guo2022}, KIC 4544587 \citep{Hambleton2013},
KIC 3858884 \citep{Manzoori2020}, and KIC 3230227 \citep{Guo2020}.
Therefore, we argue that $f_{2}$, $f_{3}$, $f_{5}$, and $f_{9}$
can be considered as the non-linearly excited daughter modes of different parent modes.
The parent mode of $f_{5}$ and $f_{9}$ is resonantly driven by a linear dynamical tide at $f_{20}$ $\simeq$ 20$f_{orb}$.
Interestingly, we did not detect the parent mode of $f_{2}$ and $f_{3}$, which was supposed to be at 14 times the orbital frequency.
The same situation exists for the eccentric binary system KIC 4544587 \citep{Hambleton2013}.
The authors concluded that the two daughter modes probably come from non-linear driving by the equilibrium tide.
Following \citep{Hambleton2013}, the parent mode of $f_{2}$ and $f_{3}$ is the component of the equilibrium tide
that pulsates at an orbital harmonic of 14$f_{orb}$.

\subsection{$\delta$ Scuti-type $p$ Modes}
In the high-frequency region, we obtained 8 independent frequencies ranging from 20 to 32 day$^{-1}$.
These frequencies are typical $\delta$ Sct-type pressure modes.
The two component stars of FX UMa are almost identical
and their physical properties agree well with those of a typical $\delta$ Scuti pulsator.
So it is very hard to ascertain which star the observed $\delta$ Sct-type pulsations originated from.
We got the pulsation constants ($Q$) of all these modes using the physical parameters of the star 1 and the equation of $Q=P_{pul}(\bar{\rho}_{1}/\bar{\rho}_{\sun})^{1/2}$,
where $P_{pul}$ is the pulsation period and $\bar{\rho}_{1}$ is the mean density, $\bar{\rho}_{1}$ = $M_{1}$/(4$\pi$$R_{1}$$^3$/3).
The $Q$ values are in the range of 0.021 to 0.033 days, suggesting low-order $p$ mode oscillations of $\delta$ Sct stars \citep{Breger2000}.

\begin{deluxetable*}{lrrrrcc}\label{tb4}
\tablecolumns{7}
\tablewidth{0pc}
\tablecaption{Independent oscillation frequencies and orbital harmonic frequencies for FX UMa}
\tablehead{\colhead{ID}   &\colhead{Frequency}   &\colhead{Amplitude}  &\colhead{Phase}   &\colhead{S/N}  &\colhead{Remark}  &\colhead{Comment}\\
\colhead{}    &\colhead{(day$^{-1}$)}     &\colhead{(mmag)}   &\colhead{(rad/2$\pi$)}   &\colhead{}    &\colhead{} &\colhead{}}
\startdata
$f_{1}$      &22.176698$\pm$0.000045     	&0.913$\pm$0.039    	&0.658$\pm$0.006    	&89.1        &\nodata		&$\delta$ Sct  \\
$f_{2}$	 &1.152015$\pm$0.000004	&0.629$\pm$0.007		&0.117$\pm$0.001		&55.1	&\nodata		&Non-linear TEO\\
$f_{3}$	 &1.954066$\pm$0.000003	&0.625$\pm$0.004		&0.508$\pm$0.001		&59.4	&\nodata		&Non-linear TEO\\
$f_{4}$	 &3.549841$\pm$0.000003	&0.562$\pm$0.003		&0.722$\pm$0.001		&80.3	&16$f_{orb}$	&Linear TEO\\
$f_{5}$	 &1.651120$\pm$0.000005	&0.528$\pm$0.007		&0.713$\pm$0.001		&47.9	&\nodata		&Non-linear TEO\\
$f_{6}$	 &20.903592$\pm$0.000011	&0.458$\pm$0.014		&0.775$\pm$0.003		&52.2	&\nodata		&$\delta$ Sct	\\
$f_{7}$	 &3.327977$\pm$0.000005	&0.454$\pm$0.003		&0.368$\pm$0.001		&60.2	&15$f_{orb}$	&Linear TEO\\
$f_{8}$	 &20.439244$\pm$0.000036	&0.378$\pm$0.017		&0.511$\pm$0.008		&47.1	&\nodata		&$\delta$ Sct	\\
$f_{9}$	 &2.786152$\pm$0.012667	&0.394$\pm$0.103		&0.255$\pm$0.264		&45.4	&\nodata		&Non-linear TEO	\\
$f_{10}$	 &26.708523$\pm$0.000022	&0.324$\pm$0.011		&0.418$\pm$0.006		&38.3	&\nodata		&$\delta$ Sct	\\
$f_{11}$	 &23.263153$\pm$0.000005	&0.327$\pm$0.004		&0.764$\pm$0.001		&33.6	&\nodata		&$\delta$ Sct	\\
$f_{14}$ 	 &2.884247$\pm$0.000012	&0.228$\pm$0.006		&0.681$\pm$0.003		&27.0	&13$f_{orb}$	&Linear TEO\\
$f_{16}$	 &1.109313$\pm$0.002787	&0.212$\pm$0.027		&0.219$\pm$0.175		&18.5	&5$f_{orb}$	&Linear TEO\\
$f_{17}$	 &22.161595$\pm$0.008500	&0.194$\pm$0.044		&0.978$\pm$0.169		&19.0	&\nodata		&$\delta$ Sct	\\
$f_{18}$	 &22.352218$\pm$0.000016	&0.194$\pm$0.007		&0.789$\pm$0.005		&18.5	&\nodata		&$\delta$ Sct	\\
$f_{19}$	 &31.599018$\pm$0.019344	&0.190$\pm$0.039		&0.144$\pm$0.107		&28.2	&\nodata		&$\delta$ Sct	\\
$f_{20}$	 &4.437298$\pm$0.000011	&0.153$\pm$0.004		&0.042$\pm$0.003		&24.7	&20$f_{orb}$ $\simeq$ $f_{5}$ + $f_{9}$ &Linear TEO \\
$f_{24}$	 &1.553085$\pm$0.002818	&0.141$\pm$0.026		&0.150$\pm$0.234		&12.7	&7$f_{orb}$	&Linear TEO\\
$f_{32}$	 &2.218628$\pm$0.000017	&0.121$\pm$0.004		&0.258$\pm$0.004		&11.6	&10$f_{orb}$	&Linear TEO\\
$f_{36}$	 &0.665672$\pm$0.000021	&0.098$\pm$0.004		&0.592$\pm$0.007		&8.5		&$3f_{orb}$	&Probable TEO\\
$f_{40}$	 &1.331203$\pm$0.016790	&0.097$\pm$0.022		&0.655$\pm$0.173		&8.5		&6$f_{orb}$	&Probable TEO\\
$f_{41}$	 &1.996715$\pm$0.000021	&0.093$\pm$0.004		&0.048$\pm$0.006		&8.9		&$9f_{orb}$	&Probable TEO\\
$f_{42}$	 &28.399938$\pm$0.004343	&0.087$\pm$0.013		&0.816$\pm$0.116		&9.5		&128$f_{orb}$ $\simeq$ $f_{10}$ + $f_{11}$ - $f_{13}$	&Probable TEO\\
$f_{43}$	 &4.659156$\pm$0.000016	&0.093$\pm$0.003		&0.928$\pm$0.006		&14.2	&$21f_{orb}$	&Linear TEO\\
$f_{64}$	 &4.881017$\pm$0.000022	&0.067$\pm$0.003		&0.249$\pm$0.008		&10.6	&22$f_{orb}$	&Linear TEO\\
$f_{65}$	 &5.102835$\pm$0.000025	&0.067$\pm$0.004		&0.676$\pm$0.006		&10.7	&23$f_{orb}$	&Linear TEO\\
$f_{74}$	 &4.215388$\pm$0.000023	&0.060$\pm$0.004		&0.347$\pm$0.008		&9.9		&$19f_{orb}$	&Probable TEO\\
$f_{78}$	 &0.887649$\pm$0.013698	&0.061$\pm$0.014		&0.798$\pm$0.175		&5.2		&$4f_{orb}$	&Probable TEO\\
$f_{97}$	 &20.411762$\pm$0.077050	&0.061$\pm$0.059		&0.733$\pm$0.280		&7.7		&92$f_{orb}$ $\simeq$ $f_{11}$ + $f_{37}$ - $f_{14}$	&Probable TEO\\
\enddata
\end{deluxetable*}

\begin{figure}
\center
\includegraphics[scale=0.8]{fig6.pdf}
\caption{The amplitude spectrum of the residual LC of FX UMa after subtracting the binarity-induced light variations.
The red, green, and blue upside-down triangles mark the frequencies of the detected non-linear TEOs, linear TEOs, and
$\delta$ Sct pulsations, respectively.
The remaining significant peaks are the combination frequencies, as summarized in Table \ref{tb5}.}
\label{fig6}
\end{figure}

\section{Summary}\label{section5}
 We report the first result of the study of pulsating eclipsing binaries
 combining high-precision $\it TESS$ photometry and high-resolution $\it SONG$ spectroscopic observations.
 In this work, we have carried out a detailed analysis of a bright double-lined spectroscopic binary FX UMa.
The major results can be summarized as follows:
\begin{enumerate}
\item FX UMa was observed by $\it TESS$ in 2 minutes cadence mode during three non-contiguous sectors: 14, 20, and 40.
The $\it TESS$ light curve shows a periodic, broad brightening near the eclipse, the typical feature of heartbeat stars.
In addition to the eclipse-like light changes and the “heartbeat-like” profile, the LC in the outside eclipses clearly display multi-periodic light variations.

\item We obtained a total of 33 $\it SONG$ high-resolution spectra from 2019 October 28 to 2020 April 29.
$\it SONG$ spectra confirm that FX UMa is a double-lined spectroscopic binary system.
The radial velocities of the two component stars were extracted by using the broadening-function technique.
The independent RV fitting yields a mass ratio of $q$ = 0.99951 $\pm$ 0.00072
and a high eccentricity of $e$ = 0.54765 $\pm$ 0.00026 for this binary system.

\item We reconstructed the individual spectra of each component star from the observed composite spectra
with the spectral disentangling tool FDBinary.
The atmospheric parameters for both component stars were then obtained through fitting their disentangled spectra.
The two components of FX UMa have almost exactly the same atmospheric parameters,
 with $T_{eff}$ = 7391.7$\pm$43.7 K, $[M/H]$ = -0.41$\pm$0.02 dex, $v_{rot}$sin$i$ = 65.40$\pm$0.80 km s$^{-1}$ for star 1
 and $T_{eff}$ = 7389.8$\pm$43.9 K, $[M/H]$ = -0.43$\pm$0.02 dex, $v_{rot}$sin$i$ = 66.79$\pm$0.83 km s$^{-1}$ for star 2.

 \item We performed a simultaneous fit to our double-lined RVs and $\it TESS$ light curves with the PHOEBE code.
 The fitting results indicate that FX UMa is a detached, eccentric binary system with
 an inclination of about 78.59 degrees.
 Combined with the radial-velocity measurements, our light-curve modeling
yields absolute masses and radii of the two components, as given in Table \ref{tab3}.
It is worth re-mentioning here that the uncertainties of the fitted parameters were seriously underestimated and should be taken just as a proof of concept.

 \item We utilized the Period04 software to extract significant frequencies from the residual LC of FX UMa,
 obtained by removing the modeled EB light curve from the original TESS observations.
 We detected a total of 103 frequencies with S/N $\geq$ 4, including 12 independent frequencies,
 17 multiples of the orbital frequency, and 74 combination frequencies.
 The eight independent frequencies in the range of 20 to 32 day$^{-1}$
 are typical low-order pressure modes of $\delta$ Scuti pulsators.
 At present it is hard to find out which star the observed $\delta$ Sct-type pulsations originated from, since the two components of FX UMa are almost identical.
 Most of the observed TEOs oscillate at harmonics of the orbital frequency ($f_{i}$ = $Nf_{orb}$).
 Consequently, $Nf_{orb}$ peaks have been viewed as the signature of tidally excited modes triggered by the linear dynamical tide.
The ten $Nf_{orb}$ peaks with S/N $>$ 10 have very high amplitudes and are considered as high-probability TEOs.
The remaining $Nf_{orb}$ peaks (4 $\leq$ S/N $\leq$ 10)  may be originated from the imperfect removal,
or they are actually real TEOs. We found that the four anharmonic frequencies ($f_{2}$, $f_{3}$, $f_{5}$, and $f_{9}$)
can pair up and sum to give exact harmonics of the $f_{orb}$:
$f_{5}$ + $f_{9}$ $\simeq$ $f_{20}$ $\simeq$ 20$f_{orb}$  and $f_{2}$ + $f_{3}$ $\simeq$ 14$f_{orb}$.
They are probably attributed to the non-linearly excited daughter modes of different parent modes
that are resonantly driven by the linear dynamical tide.
\end{enumerate}

Heartbeat stars with TEOs provide unique opportunities to test theories of stellar tides and their interaction with pulsation,
and with angular momentum \citep{Kurtz2022}.
As summarized in Table 1 of \cite{Guo2021}, 22 heartbeat binaries have been observed to show tidally excited oscillations,
but only a handful of them have been studied in detail.
The discovery of linear and non-linear tidal oscillations in the SB2 system FX UMa
presents us with a new opportunity.
During our analysis, additional $\it TESS$ observations for this object were released for Sector 47.
However, the two-minute cadence data of FX UMa in Sector 47 were not available at the time of publication.
Based on the Web $\it TESS$ Viewing Tool \footnote{\url{https://heasarc.gsfc.nasa.gov/cgi-bin/tess/webtess/wtv.py}},
we observe further that FX UMa will be observed by $\it TESS$ during Sector 60.
Long time series of $\it TESS$ photometry help to resolve individual pulsations.
We expect more TEOs will be reported in the future.

\begin{acknowledgments}
This work is supported by the Sichuan Science and Technology Program (Grant No. 2020YFSY0034),
the National Natural Science Foundation of China (Grant No. 12003022, 12173028),
the National Key Basic R\&D Program of China (Grant No. 2021YFA1600401),
the Sichuan Youth Science and Technology Innovation Research Team (Grant No. 21CXTD0038),
and the Innovation Research Team funds of China West Normal University (Grant No. KCXTD2022-6).
R.A.B acknowledges the support from the Major Science and Technology Project of Qinghai Province (Grant No. 2019-ZJ-A10).
This research is based on observations made at the Hertzsprung SONG telescope
operated at the Spanish Observatorio del Teide on the island of Tenerife
by the Aarhus and Copenhagen Universities and by the Instituto de Astrofísica de Canarias.
Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement no.: DNRF106)

\end{acknowledgments}

\software{$\it Lightkurve$ \citep{Lightkurve2018},
iSpec \citep{Blanco2014, Blanco2019},
BF-rvplotter (\url{https://github.com/mrawls/BF-rvplotter}),
LMFIT \citep{Newville2021},
$\it rvfit$ \citep{Iglesias2015},
FDBinary \citep{Ilijic2004},
PHOEBE, \citep{Andrej2018, Conroy2020},
emcee \citep{Foreman2019, Foreman2013},
Period04 \citep{Lenz2005}.
}