\section{Introduction}
Nearly all surfaces in fluid-related industries are rough at their operating Reynolds numbers. These rough surfaces easily alter the aero- or hydrodynamic properties of a fluid system and induce unfavorable consequences; for instance, reduced stability of aircraft, increased fuel consumption of cargo ships, and reduced energy harvesting capacity of wind turbines~\citep{gent2000aircraft,schultz2011economic,dalili2009review}. Therefore, it is crucial to restore or replace surfaces that have degraded past a critical level. This requires monitoring roughness on industrial surfaces and efficiently predicting their effect on the flow.

One of the most important effects of surface roughness on flow is the increased drag. This additional drag is generally expressed by the Hama roughness function~\citep{hama1954boundary} or the equivalent sand-grain height. The Hama roughness function is the downward shift of the mean velocity profile  induced by roughness. The equivalent sand-grain height is the roughness height from the Moody diagram~\citep{moody1944friction} that would cause the same drag as the rough surface of interest. It is worth noting that to calculate the equivalent sand-grain height, the Hama roughness function in the fully-rough regime -- where the skin-friction coefficient is independent of Reynolds number -- needs to be determined.
To accurately determine the Hama roughness function of irregular surfaces, numerical simulations with fully resolved rough structures or experiments are needed. However, neither simulations nor experiments are cost-efficient, hence not feasible for practical purposes. To alleviate the high costs, researchers have found empirical correlations between rough surface statistics and their contribution to drag~\citep{chan2015systematic,forooghi2017toward,flack2010review,flack2020skin}. However, no universal model has emerged that is able to accurately predict the drag from a sufficiently wide range of roughness. This is mainly attributed to the high-dimensional feature space of irregular rough structures~\citep {chung2021predicting}.

Neural networks are known to be capable of extracting patterns in high-dimensional spaces, and therefore have been successfully used in various studies using fluid flow data~\citep{lee2019data,kim2020prediction,fukami2021machine}. \citet{jouybari2021data} developed a multilayer perceptron (MLP) type neural network to find a mapping of 17 different rough surface statistics to equivalent sand-grain height. They reported a state-of-the-art performance in predicting equivalent sand-grain heights from rough surface statistics. A total of 45 data samples, obtained from direct numerical simulations (DNS) and experiments, were used for training and validating their neural network. Having a large number of data samples to train a neural network is invariably advantageous; nonetheless, constructing a fluid database that is sufficiently large to train a network from scratch imposes significant computational or experimental cost and is oftentimes considered impractical. Therefore, for the practical usage of neural networks in fluid-related applications, a framework that allows neural networks to learn a generalized mapping from a limited number of data samples is needed.

In this study, we propose a transfer learning framework to improve the generalization ability of neural networks for predicting the drag on rough surfaces when only a limited number of high-fidelity data samples are available. Transfer learning is a method that improves learning performance by adapting knowledge learned from one domain to another~\citep{pan2009survey,chakraborty2021transfer}.
The method has been used, for example, for self-driving cars, which can be pre-trained with data from virtual car simulators~\citep{pan2017virtual,martinez2017beyond,akhauri2020enhanced}. Obviously, the driving environment in virtual car simulators is inherently different from the real world, however, pre-training a neural network with `approximate knowledge' can provide better initial weights for training real self-driving cars.

Similarly, for flow over rough surfaces, empirical correlations between drag and surface roughness provide an approximate knowledge of real flow physics~\citep{colebrook1937experiments,moody1944friction,chan2015systematic,forooghi2017toward,flack2010review,flack2020skin}. These correlations were developed by a fitting procedure of large experimental and numerical datasets. However, it is not straightforward to make direct use of these datasets, in particular because the information of surface topographies and flow statistics is not always accessible. Indeed, in many cases we have advanced our understanding of physics from the valuable information embedded in empirical correlations, such as those developed by~\cite{colebrook1937experiments} and~\citet{moody1944friction}.

The aim of this study is to show that transfer learning of empirical correlations can significantly improve the performance of neural networks for modeling drag on rough surfaces. Our aim is also to analyze a simple neural network with and without transfer learning to gain insight of the learning behavior, which is strongly connected to the physics of the problem. The objective of the developed neural network model is to predict drag of a particular class of roughness contained in one high-fidelty dataset, which is composed of irregular homogeneous rough surfaces. However, we foresee that transfer learning will play a central role in combining empirical relations with larger datasets from a range of sources to  develop  predictive models of a significantly larger class of rough walls than the ones considered here. The paper is organized as follows: the details of the developed transfer learning framework including the neural networks and datasets are explained in \S\ref{sec:TF}. The results of learning a mapping of surface statistics to the Hama roughness function using the transfer learning framework are shown in~\S\ref{sec:Result}. The effects of transfer learning are analyzed and discussed in~\S\ref{sec:discussion}, followed by the concluding remarks in~\S\ref{sec:conclusion}.
\begin{figure}
  \centerline{\includegraphics[width = 1.0 \linewidth]{NNTransferlearning.eps}}
  \caption{Overview of the transfer learning framework.}
\label{fig:tf_overview}
\end{figure}

\section{Transfer learning framework}\label{sec:TF}
The developed transfer learning framework is composed of two parts: (1) pre-training step and (2) fine-tuning step (figure~\ref{fig:tf_overview}). In the pre-training step, neural networks were trained to learn an `approximate’ knowledge of the mapping of surface statistics to the Hama roughness function. The training data is created by evaluating empirical correlation functions given the surface statistics of different surface topographies. In the fine-tuning step, neural networks were tuned to learn high-fidelity physics from a small DNS dataset by adapting their pre-trained `approximate’ knowledge to the domain of real physics. In \S\ref{sec:Result}, we show that transfer learning of empirical correlations improves the generalization ability of neural networks when learning the drag of rough surfaces from a small amount of data. First, however, we explain the details of the neural networks and pre-training and fine-tuning steps.

\subsection{Neural network}
We used a variation of the MLP type neural network architecture developed by~\citet{jouybari2021data}. The input of the employed network is a set of 17 different rough surface statistics calculated from the surface topography of interest. The 17 input parameters are composed of 8 primary surface statistics and 9 products of the 8 primary statistics. These parameters are selected based on their perceived importance on affecting drag on rough surfaces~\citep{jouybari2021data}.

Let $x$, $y$, and $z$ be the streamwise, wall-normal and spanwise directions, $k(x,z)$ the roughness height distribution and $A_t$ the total roughness plan area. Then, one may define the following measures of an irregular rough surface:
\begin{equation}
    \begin{array}{llcl}
   \textrm{average roughness height: }&k_{avg}&=&\int_{x,z}k dA / A_{t},\\
\textrm{average height deviation: }&Ra  &=&\int_{x,z}(|k-k_{avg}|) dA / A_{t},\\
\textrm{root-mean-squared roughness height: } & k_{rms} & =&\left (\int_{x,z}(k-k_{avg})^{2}dA/A_{t}\right )^{1/2}.
    \end{array}
\end{equation}
In addition, we define the crest height $k_{c}$, and finally, the fluid area at each wall-normal plane $A_{f}(y)$. From these quantities, the 17 statistical input parameters ($\{I_{1}, I_{2}, ..., I_{17}\}$) can be calculated as listed in table~\ref{tab:I}.
\begin{table}
  \begin{center}
\def~{\hphantom{0}}
  \begin{tabular}{ll}
      $I_{1} = k_{rms}/Ra$  & $I_{2} = Sk = \frac{1}{A_{t}}\int_{x,z}(k-k_{avg})^{3} dA / k^{3}_{rms}$ \\
      $I_{3} = kur = \frac{1}{A_{t}}\int_{x,z}(k-k_{avg})^{4} dA / k^{4}_{rms}$ & $I_{4} = ES_{x} = \frac{1}{A_{t}}\int_{x,z}|\frac{\partial k}{\partial x}| dA$ \\
      $I_{5} = ES_{z} = \frac{1}{A_{t}}\int_{x,z}|\frac{\partial k}{\partial z}| dA$  & $I_{6} = por = \frac{1}{A_{t}k_{c}}\int_{0}^{k_{c}} A_{f} dy$ \\
      $I_{7} = inc_{x} = \tan^{-1}\{\frac{1}{2}Sk(\frac{\partial k}{\partial x})\}$   & $I_{8} = |inc_{z}| = |\tan^{-1}\{\frac{1}{2}Sk(\frac{\partial k}{\partial z})\}|$ \\
      $I_{9} = I_{4}I_{5}$  & $I_{10} = I_{4}I_{2}$ \\
      $I_{11} = I_{4}I_{3}$ & $I_{12} = I_{5}I_{2}$ \\
      $I_{13} = I_{5}I_{3}$ & $I_{14} = I_{2}I_{3}$ \\
      $I_{15} = I_{4}I_{4}$ & $I_{16} = I_{5}I_{5}$ \\
      $I_{17} = I_{2}I_{2}$ &
  \end{tabular}
  \caption{Statistical input parameters of the employed neural network.}
  \label{tab:I}
  \end{center}
\end{table}
Here, $Sk$, $kur$, $por$, $ES_{x}$, $ES_{z}$, $inc_{x}$, and $inc_{z}$ indicate skewness, kurtosis, porosity, and effective slopes and average inclinations in $x$ and $z$ directions, respectively. Before training, the parameters $I_{1}$ to $I_{17}$ were normalized with the mean and standard deviation values calculated from the dataset of empirical correlations (see section~\ref{sec:pre-train} for the details of the empirical dataset).

\begin{figure}
  \centerline{\includegraphics[width = 0.9 \linewidth]{NN.eps}}
  \caption{Architecture of the used neural network.}
\label{fig:NN}
\end{figure}

Information of input parameters travels through three hidden layers with 18, 7, and 7 neurons, respectively (see figure~\ref{fig:NN}). Leaky rectified linear unit (ReLU) activations ($\max(0.01x,x)$) were applied after the hidden layers. The output of the employed neural network is a scalar value of the Hama roughness function $\Delta U^{+}$. The superscript $+$ indicates the normalization by the viscous scale $\delta_{\nu}=\nu/u_{\tau}$, where $\nu$ and $u_{\tau}$ are the kinematic viscosity and the friction velocity, respectively. In this study, the friction Reynolds number $Re_{\tau}=u_{\tau} H/ \nu=500$ was used, where $H$ is the channel half-height. At the output layer, a Sigmoid activation function, $a/(1+e^{x})$,  was applied to bound the prediction value ($a=20$ in this study).

\subsection{Pre-training step}\label{sec:pre-train}
In the pre-training step, neural networks were trained with a large dataset obtained from empirical correlation functions between rough surface statistics and the Hama roughness function $\Delta U^{+}$. We used three correlations;
\begin{equation}
  \Delta U^{+} = 2.5 \log Ra^{+} + 1.12 \log ES_{x} + 1.47;\label{EMP1}\tag{EMP1}
\end{equation}
\begin{equation}
    \Delta U^{+}= 2.5 \log \left [k^{+}_{c}\left (0.67 Sk^{2}+0.93 Sk + 1.3\right)\left (1.07\left (1-e^{-3.5 ES_{x}}\right )\right)\right ]- 3.5; \label{EMP2}\tag{EMP2}
\end{equation}
and
\begin{equation}
  \Delta U^{+} = \left\{
    \begin{array}{ll}
      2.5 \log \left [2.48 k^{+}_{rms}(1+Sk)^{2.24}\right ] - 3.5, & Sk > 0 \\
            2.5 \log \left [2.11 k^{+}_{rms}\right ] - 3.5, & Sk \approx 0, \\
      2.5 \log \left [2.73 k^{+}_{rms}(2+Sk)^{-0.45}\right ] - 3.5, & Sk < 0.
    \end{array} \right. \label{EMP3}\tag{EMP3}
\end{equation}
Here, \ref{EMP1} was proposed by \citet{chan2015systematic}, \ref{EMP2} by \citet{forooghi2017toward} and \ref{EMP3} by \citet{flack2010review,flack2020skin}.

As mentioned earlier, we will develop neural networks to predict drag of rough surfaces contained in a small high-fidelity dataset (see \S\ref{sec:finetune} for the details of the DNS dataset). The range of the roughness function in this dataset is $5.3<\Delta U^{+}<8.3$. The lower bound lies at the boundary between transitionally and fully rough regimes. \ref{EMP1} was derived using surfaces in the transitionally and fully rough regimes with $0.4<\Delta U^{+}<11.4$. \ref{EMP2} was modeled from surfaces with roughness function in the range of $5.4<\Delta U^{+}<10.1$. Finally, \ref{EMP3} was developed with rough surfaces in the fully rough regime, where most of the surfaces become fully rough when $\Delta U^{+}$ was over $5.5-6.0$.  Therefore, all three empirical models contain `approximate' knowledge of the drag of rough surfaces that can help the neural networks to learn the DNS dataset in the fine-tuning step.

A large dataset was constructed by randomly generating 10,000 irregular rough surfaces and calculating the associated Hama roughness function using \ref{EMP1}, \ref{EMP2} and \ref{EMP3}. The surfaces were constructed using the Fourier-filtering algorithm suggested by \citet{jacobs2017quantitative}. The algorithm generates self-affine surfaces with a pre-defined power spectrum $C$,
\begin{equation}
 C(q) \propto q^{-2-2h}
\label{eq:self-affine}
\end{equation}
where $q$ is the wavenumber and $h$ the Hurst exponent. The Hurst exponent $h=0.8$ was used as in \citet{jacobs2017quantitative}. This power-law dependence is a well-known attribute of self-affine realistic surfaces~\citep{mandelbrot1982fractal,persson2004nature}.
Figure~\ref{fig:emp}a shows a few examples of the generated random surfaces. The randomness is imposed by choosing random amplitudes of the power spectrum and random phase shifts of Fourier modes. Accordingly, the surface statistics, such as the skewness and the effective slope, depend on the particular combination of the random amplitudes and phases of Fourier modes.

For training and validating the neural network, 7,000 and 3,000 data samples were used, respectively. Figure~\ref{fig:emp}b shows the scatter plots of the Hama roughness functions, calculated by \ref{EMP1}, \ref{EMP2} and \ref{EMP3}, against $k_{rms}^{+}$, $Sk$, and $ES_{x}$ from the 3,000 validation data samples.
The surface generation method based on the spectral density \eqref{eq:self-affine}, has generated surfaces with skewness bounded between $-1$ to $+1$ and effective slopes limited to $\lesssim 0.35$.
The bounded skewness and limited effective slopes are due to the wavy surfaces generated by Fourier modes and the self-affine power spectrum defined for realistic surfaces, as in~\citet{jacobs2017quantitative}.
The random surfaces can thus be considered to be in the low-slope regime, where the drag tends to increase with increasing effective slope~\citep{flack2014roughness}.
From figure~\ref{fig:emp}b, we observe different distributions, or knowledge, of the Hama roughness function from the different empirical models. The aim is to seed this diverse knowledge of empirical correlations in neural networks during the pre-training step.

\begin{figure}
  \centerline{\includegraphics[width = 1.0\linewidth]{random_surface_merged.eps}}%
  \caption{(a) Examples of patches of the generated surfaces in the empirical dataset. The size of the patches is $500 \times 500$ (viscous unit) in the streamwise and spanwise directions. The visualized heights of the surfaces are exaggerated five times in order to show the topography explicitly. (b) Scatter plots of the Hama roughness functions ($\Delta U^{+}$), calculated by \ref{EMP1}, \ref{EMP2}, and \ref{EMP3}, against $k_{rms}^{+}$, $Sk$, and $ES_{x}$ from the 3,000 validation data samples.}
\label{fig:emp}
\end{figure}

Three neural networks with the same architecture and the same initial weights and biases were trained simultaneously to learn the different empirical models, \ref{EMP1} to \ref{EMP3}. Let $NN^{i}$ be the neural network trained with the $i$-th empirical model, $W^{i}_{j,k}$ the weight matrix connecting the $j$-th and $k$-th layers of $NN^{i}$, $N_{w}$ the total number of weights in $NN^{i}$, $b^{i}_{j}$ the bias vector added in the $j$-th layer, $\Delta U^{+}_{i}$ the Hama roughness function predicted by the $i$-th empirical model, and $\Delta \widetilde{U}^{+}_{i}$ the Hama roughness function predicted by $NN^{i}$. The neural networks were trained to minimize a loss for pre-training:
\begin{equation}
    L_{p} = \frac{1}{3}\sum_{i=1}^{3} \left[ \left(\Delta U^{+}_{i} - \Delta \widetilde{U}^{+}_{i} \right)^{2} + \sum_{j,k} \mathrm{tr}\left\{(W^{i}_{j,k})^{T}(W^{i}_{j,k})\right\} \right].
\end{equation}
The first term in the right-hand side is the mean-squared-error loss and the second is the weight regularization loss, also used in \citet{jouybari2021data}. The sizes of the weight matrix $W^{i}_{j,k}$ and bias vector $b^{i}_{k}$ are $n_{j} \times n_{k}$ and $n_{k}$, respectively, where $n_{j}$ and $n_{k}$ are the number of neurons on the $j$-th and $k$-th layers.

The weights and biases of the neural networks were updated in the direction of minimizing the pre-training loss $L_{p}$ on the training data samples through the Adam optimizer~\citep{kingma2014adam}. The pre-training loss $L_{p}$ on the validation data samples was also evaluated each iteration of updates. Batch sizes of 16 and 32 were used for training and validating, respectively. To provide a wide distribution of data samples in the batches, i.e., preventing neural networks from learning biased data focused on average characteristics, we imposed batches to contain 50\% of the data samples with $\Delta U^{+}$ in the range 6 to 8, 25\% of the data samples with $\Delta U^{+}$ in the range 4 to 6, and the remaining 25\% of the data samples with $\Delta U^{+}$ in the range 8 to 10.
After a sufficient amount of iterations for training and validating, the set of networks ($\{NN^{1},NN^{2},NN^{3}\}$) that produced the lowest pre-training loss on the validation data samples was chosen to be merged into a single pre-trained network $NN^{pre}$. The learned sets of weights and biases in the three neural networks, $\{W^{1}_{j,k},W^{2}_{j,k},W^{3}_{j,k}\}$ and $\{b^{1}_{k},b^{2}_{k},b^{3}_{k}\}$, were averaged into the weights and biases of $W^{pre}_{j,k}$ and $b^{pre}_{k}$ as $W^{pre}_{j,k} = \sum_{i=1}^{3} W^{i}_{j,k}/3$ and $b^{pre}_{k} = \sum_{i=1}^{3} b^{i}_{k}/3$.
These pre-trained weights and biases were employed as the initial weights for the fine-tuning step (section~\ref{sec:finetune}) to transfer the knowledge of empirical correlations to the domain of high-fidelity physics.

A key factor for the performance of the fine-tuned neural network is the contents of database used in the pre-training step. If the bound of the $\Delta U^{+}$ in one's high fidelity dataset is known \textit{a priori}, then by using a narrow range of $\Delta U^{+}$ for pre-training increases the transfer learning performance significantly. We provide an example of this by training the neural networks with a bound of $\Delta U^{+}$ in Appendix~\ref{appBU}. In the main text, we assume that the range of $\Delta U^{+}$ of the high-fidelity dataset is \textit{a priori} unknown.

\subsection{Fine-tuning step}\label{sec:finetune}
A total of 35 DNS of channel flow over rough surfaces at $Re_{\tau}\approx500$ were conducted to construct a high-fidelity dataset. The 35 rough surfaces were generated following the method proposed by \citet{perez2019generating}. Figure~\ref{fig:dns}a shows a few examples of the generated surfaces. Note that different surface generation methods were used to construct the DNS and empirical datasets (section~\ref{sec:pre-train}). The method used for DNS imposes an additional constraint on the roughness height probability distribution. This constraint enables the generation of rough surfaces with distinctive characteristics. For the probability distribution, we used Weibull distributions of the form
\begin{equation}
 f(x;\lambda,s) = \frac{s}{\lambda}\left(\frac{x}{\lambda}\right)^{s-1}e^{-(x/\lambda)^{s}}.
 \label{eq:Weibull}
\end{equation}
By using different shape ($s$) and scale ($\lambda$) parameters, rough surfaces with non-zero skewness are generated. Gaussian distributions were used to generate rough surfaces with zero skewness. Rough surfaces generated with \eqref{eq:Weibull} imitate the roughness caused by wear~\citep{perez2019generating}. By using different surface generation methods in the pre-training and fine-tuning steps, we can verify that pre-training on one class of rough surfaces can help learn the drag on a different class of rough surfaces. Such verification is important, as in practice, artificially generated rough surfaces for pre-training would not perfectly imitate real-world rough surfaces.
\begin{figure}
  \centerline{\includegraphics[width = 1.0 \linewidth]{KIT_surface_merged.eps}}
  \caption{(a) Examples of patches of the generated surfaces in the DNS dataset. The size of the patches is $500 \times 500$ (viscous unit) in the streamwise and spanwise directions. The visualized heights of the surfaces are exaggerated five times in order to show the topography explicitly. (b) Scatter plots of the calculated Hama roughness function ($\Delta U^{+}$) against various surface statistics of the 35 DNS data samples. The markers {\color{red}$\circ$}, {\color{green}$\lhd$}, and $\times$ indicate the training, validation, and test data samples, respectively. The plots were drawn from an example of a combination of training and validation data samples.}
\label{fig:dns}
\end{figure}

\begin{table}
  \begin{center}
\def~{\hphantom{0}}
\begin{tabular}{lllllllllll}
Fine-tuning data    & $\Delta U^{+}$ & $k_{rms}^{+}$   & $Sk$    & $ES_{x}$    & $ES_{z}$    & $inc_{x}$     & $inc_{z}$     & $Ra^{+}$    & $kur$   & $por$   \\
F1  & 5.46  & 8.712  & -0.721 & 0.178 & 0.180 & 0.092  & 0.008  & 6.681 & 4.186 & 0.386 \\
F2  & 5.51  & 10.408 & -0.484 & 0.371 & 0.370 & -0.019 & -0.006 & 8.406 & 2.954 & 0.373 \\
F3  & 6.05  & 11.385 & -0.568 & 0.549 & 0.533 & 0.006  & 0.048  & 9.411 & 2.734 & 0.347 \\
F4  & 6.09  & 10.407 & -0.483 & 0.543 & 0.542 & 0.003  & 0.000  & 8.405 & 2.951 & 0.373 \\
F5  & 6.14  & 10.411 & -0.483 & 0.557 & 0.558 & 0.010  & 0.002  & 8.407 & 2.952 & 0.373 \\
F6  & 6.60  & 8.578  & 0.868  & 0.248 & 0.247 & -0.020 & -0.049 & 6.832 & 3.664 & 0.711 \\
F7  & 6.84  & 9.453  & -0.001 & 0.363 & 0.361 & -0.012 & -0.005 & 7.547 & 2.943 & 0.498 \\
F8  & 7.36  & 9.417  & 0.755  & 0.456 & 0.435 & 0.020  & 0.057  & 7.394 & 3.759 & 0.638 \\
F9  & 7.70  & 11.546 & 0.433  & 0.260 & 0.261 & -0.071 & -0.072 & 9.328 & 2.890 & 0.612 \\
F10 & 8.29  & 10.956 & 0.718  & 0.476 & 0.476 & 0.009  & 0.061  & 8.786 & 3.336 & 0.681 \\ \hline \hline
Test data & $\Delta U^{+}$ & $k_{rms}^{+}$   & $Sk$    & $ES_{x}$    & $ES_{z}$    & $inc_{x}$     & $inc_{z}$     & $Ra^{+}$    & $kur$   & $por$   \\
T1        & 5.32 & 9.035  & -0.770 & 0.187 & 0.176 & 0.022  & -0.062 & 6.902 & 4.306 & 0.384 \\
T2        & 5.63 & 10.406 & -0.483 & 0.399 & 0.392 & 0.030  & 0.033  & 8.397 & 2.956 & 0.372 \\
T3        & 5.88 & 9.986  & 0.003  & 0.372 & 0.373 & 0.012  & -0.019 & 7.978 & 2.936 & 0.504 \\
T4        & 5.97 & 8.958  & 0.274  & 0.242 & 0.240 & 0.022  & -0.047 & 7.256 & 2.745 & 0.567 \\
T5        & 6.17 & 8.667  & 0.132  & 0.289 & 0.289 & -0.033 & -0.013 & 6.897 & 3.030 & 0.521 \\
T6        & 6.31 & 10.439 & -0.730 & 0.230 & 0.223 & 0.091  & 0.080  & 7.993 & 4.243 & 0.387 \\
T7        & 6.39 & 9.978  & 0.001  & 0.387 & 0.382 & -0.037 & 0.096  & 7.972 & 2.934 & 0.505 \\
T8        & 6.46 & 8.889  & 0.001  & 0.271 & 0.274 & -0.012 & -0.005 & 7.100 & 2.937 & 0.498 \\
T9        & 6.46 & 9.481  & -0.133 & 0.362 & 0.363 & -0.005 & 0.015  & 7.550 & 3.015 & 0.479 \\
T10       & 6.57 & 10.408 & 0.484  & 0.371 & 0.370 & 0.019  & 0.006  & 8.406 & 2.954 & 0.627 \\
T11       & 6.59 & 9.986  & 0.003  & 0.529 & 0.517 & -0.020 & 0.028  & 7.977 & 2.937 & 0.505 \\
T12       & 6.60 & 9.512  & -0.004 & 0.194 & 0.193 & -0.068 & 0.052  & 7.597 & 2.943 & 0.495 \\
T13       & 6.65 & 9.976  & 0.001  & 0.489 & 0.472 & 0.011  & -0.007 & 7.972 & 2.926 & 0.505 \\
T14       & 6.67 & 9.980  & 0.002  & 0.537 & 0.535 & -0.006 & 0.002  & 7.973 & 2.936 & 0.502 \\
T15       & 6.80 & 9.992  & 0.003  & 0.574 & 0.551 & 0.058  & 0.025  & 7.983 & 2.933 & 0.504 \\
T16       & 6.92 & 10.429 & 0.491  & 0.444 & 0.441 & -0.116 & 0.027  & 8.416 & 2.972 & 0.629 \\
T17       & 7.00 & 9.208  & 0.864  & 0.266 & 0.267 & 0.027  & -0.012 & 7.337 & 3.633 & 0.709 \\
T18       & 7.23 & 10.407 & 0.483  & 0.543 & 0.542 & -0.003 & 0.000  & 8.405 & 2.951 & 0.627 \\
T19       & 7.27 & 10.706 & 0.120  & 0.355 & 0.355 & 0.051  & 0.007  & 8.632 & 2.751 & 0.522 \\
T20       & 7.28 & 10.400 & 0.480  & 0.554 & 0.547 & 0.004  & 0.017  & 8.399 & 2.946 & 0.627 \\
T21       & 7.34 & 9.419  & 0.754  & 0.398 & 0.380 & 0.086  & 0.029  & 7.397 & 3.755 & 0.638 \\
T22       & 7.35 & 10.886 & 0.491  & 0.284 & 0.283 & -0.042 & -0.018 & 8.792 & 2.957 & 0.627 \\
T23       & 7.59 & 10.262 & -0.418 & 0.815 & 0.808 & 0.015  & 0.007  & 8.023 & 3.573 & 0.429 \\
T24       & 7.70 & 10.688 & 0.001  & 0.557 & 0.552 & 0.000  & 0.000  & 8.532 & 2.944 & 0.496 \\
T25       & 7.77 & 11.233 & -0.001 & 0.482 & 0.483 & -0.007 & 0.041  & 8.970 & 2.946 & 0.495
\end{tabular}
  \caption{The data of roughness function ($\Delta U^{+}$) and surface statistics used for fine-tuning and testing the neural networks.}
  \label{tab:data}
  \end{center}
\end{table}

Details of the methodology and validation of the conducted DNS are reported in the study of~\citet{jiasheng2021dns}. Here, we provide a brief overview.
A pseudo-spectral incompressible Navier–Stokes solver, SIMSON~\citep{chevalier2007simson} with an immersed boundary method~\citep{goldstein1993modeling}, was employed for the channel simulations. Periodic boundary conditions were used in the streamwise and spanwise directions. No-slip conditions were applied at the top and bottom planes of the channel. Roughness effects were added to both the top and bottom planes. A minimal channel approach was adopted to find a small computational domain size that can accurately reproduce the Hama roughness function calculated from a full channel simulation~\citep{chung2015fast}.
Further details on the domain size and grid resolution are provided in Appendix~\ref{appDNS}.

The scatter plots of the calculated Hama roughness function on the 35 data samples against the surface statistics involved in calculating the input of neural networks are shown in figure~\ref{fig:dns}b. The generated surfaces in the fine-tuning and pre-training steps show moderately different characteristics. For instance, the surfaces generated in the fine-tuning step show higher limits of the effective slopes (up to $\sim 0.8$) compared to those in the pre-training step (up to $\sim 0.35$). However, note that despite the larger effective slope, the surfaces are not in the regime where drag tend to decrease with increasing effective slope~\citep{flack2014roughness}. The 35 data samples were split into 6, 4, and 25 data samples for training, validation, and test, respectively. More than 70\% of the data samples were used for testing, i.e.,~not used during the fine-tuning step, to fairly evaluate the generalization ability of the developed neural networks. Note that a total of 10 data samples were used in the fine-tuning step, and these fine-tuning data samples are completely separated from the test data samples. The data used for fine-tuning and testing the neural networks is provided in table~\ref{tab:data}. The surface statistics and the Hama roughness function in the test dataset are distributed widely in the total dataset as shown by the $\times$ markers in figure~\ref{fig:dns}b.

To reduce the uncertainty arising from splitting a very small dataset into training and validation sets, we adopted an approach based on averaging many neural networks.
More specifically, we employed 210 neural networks, $\{NN^{1}, NN^{2}, ... , NN^{210}\}$, to learn from the 210 different data combinations derived by selecting 6 training and 4 validation data samples from the 10 fine-tuning data samples
(${10 \choose 6} ={10 \choose 4} =10!/(6!\,4!)=210$).
The weights and biases of these 210 neural networks ($\{W^{1}_{j,k},W^{2}_{j,k},...,W^{210}_{j,k}\}$ and $\{b^{1}_{k},b^{2}_{k},.., b^{210}_{k}\}$) were initialized by the pre-trained weights and biases, $W^{pre}_{j,k}$ and $b^{pre}_{k}$. Similar to the pre-training step, the weights and biases of each neural network $NN^{i}$ were updated by the Adam optimizer that minimizes a loss for fine tuning $L^i_{f}$,
\begin{equation}
    L^i_{f} = \left (\Delta U^{+} - \Delta \widetilde{U}_{i}^{+}\right )^{2} + \frac{1}{N_w} \sum_{j,k} \mathrm{tr}\left\{(W^{i}_{j,k})^{T}(W^{i}_{j,k})\right\}.
\end{equation}
Here, $\Delta U^{+}$ is the ground truth, and $\Delta \widetilde{U}_{i}^{+}$ is the roughness function predicted by $NN^{i}$.
After a sufficient amount of updates, the networks with the smallest fine-tuning loss on the validation data samples were chosen to predict the roughness function $\Delta U^{+}$ on test data (section~\ref{sec:Result}).

\section{Prediction of roughness functions}\label{sec:Result}
We compare the predicted Hama roughness functions on the 25 test data samples that were not included in the fine-tuning step. The roughness functions are calculated by the following prediction methods: (1) neural networks trained with transfer learning (NNTF), (2) neural networks trained without transfer learning (NN), and (3) empirical correlations (EMP). Neural networks in NN were trained without the pre-training step, i.e., the weights and biases of the networks were randomly initialized instead of being initialized by $W^{pre}_{j,k}$ and $b^{pre}_{k}$. Our focus here is on the comparison between NNTF and NN to demonstrate the effects of transfer learning. However, we also include EMP predictions in the comparison since these methods are widespread in the literature.

The performance of neural networks can vary greatly depending on the selection of data for training and validation. This is because it is generally hard to form a distribution covering the entire feature space from a limited number of data samples. By ensemble averaging the predictions from neural networks learned from different data combinations, we can partly alleviate this problem.
The ensemble predictions from NNTF and NN are thus obtained by averaging the prediction values of the 210 neural networks as,
\begin{equation}
    \Delta \widetilde{U}^{+} = \sum_{i=1}^{210} \Delta \widetilde{U}_{i}^{+}/210.
    \label{eq:u_ensemble}
\end{equation}
Similarly, we use ensemble predictions of the empirical models $\{$\ref{EMP1}, \ref{EMP2}, \ref{EMP3}$\}$ in EMP as,
\begin{equation}
    \Delta \widetilde{U}_{emp}^{+} = \sum_{i=1}^{3} \Delta \widetilde{U}_{i,emp}^{+}/3,
\end{equation}
where $\Delta \widetilde{U}_{i, emp}^{+}$ is the Hama roughness function predicted by the $i$-th empirical model.

\begin{figure}
  \centering
  \includegraphics[width = 0.5
  \linewidth]{Ensemble_err.eps}
  \caption{Average and maximum errors from the ensemble predictions performed by NNTF, NN, and EMP.}
\label{fig:ensemble}
\end{figure}

Figure~\ref{fig:ensemble} shows the maximum and average errors calculated from the ensemble predictions by NNTF, NN, and EMP. The predictions from NN exhibit 7.62\% and 23.05\% of average and maximum errors, respectively, while the corresponding error associated with the predictions from NNTF are 6.25\% and 14.15\%. Note that both NN and NNTF are performing better than EMP here.
The most significant advantage of using transfer learning is in the reduction of the maximum error; nearly ten percent of decrease is achieved by using NNTF instead of NN. As we will discuss more quantitatively in the next section, this error reduction demonstrates the capability of NNTF in learning a better generalized mapping from a small amount of DNS data.

To further investigate transfer learning, the best and worst performing neural networks out of the 210 networks in NNTF and NN were extracted and compared in figure~\ref{fig:bestandworst}. The errors from the best performing networks in NNTF and NN (figure~\ref{fig:bestandworst}a) show a similar trend as the ensemble predictions (figure~\ref{fig:ensemble}); both average and maximum errors are clearly reduced by using transfer learning. We also note that the ensemble predictions from NNTF (figure~\ref{fig:ensemble}) marginally perform better compared to both the best performing neural network in NN and empirical correlation in EMP (figure~\ref{fig:bestandworst}a).

The advantage of using transfer learning is more clearly demonstrated in the worst performing case (figure~\ref{fig:bestandworst}b). The maximum error of the worst performing neural network in NN is nearly 60\%, showing that this network has failed to provide a reliable prediction. On the other hand, the network in NNTF exhibits significantly smaller errors, indicating that transfer learning is helping the network to provide a safer prediction even in the worst scenario. Therefore, these results show that transfer learning enables the networks to learn a more generalized mapping from limited DNS data by adapting an `approximate’ knowledge from empirical correlations. It is also important to emphasize that ensemble predictions should be used in practice to reduce the uncertainty arising from a limited pool of data as can be seen from the wide range of errors that occur with the best and worst performing neural networks in figure~\ref{fig:bestandworst}.

\begin{figure}
  \centerline{\includegraphics[width = 1.0 \linewidth]{Bestandworst_err.eps}}
  \caption{Average and maximum errors from the (a) best and (b) worst performing neural networks and empirical correlations.}
\label{fig:bestandworst}
\end{figure}

\section{Analysis of transfer learning and approximate knowledge}\label{sec:discussion}
In the previous section, we demonstrated that transfer learning of empirical correlations improves the learning performances of neural networks. In this section, we analyze the effects of transfer learning and why `approximate knowledge' in the empirical correlations can help neural networks.

\subsection{Effects of transfer learning on weight generalization}\label{sec:weights}
In this section, we investigate how transfer learning improves the generalization ability by characterizing the learning weights of the networks. Let $w^{i}_{l}$, where $l=1,2,...,N_{w}$, be the vectorized elements of the weight matrices $W^{i}_{j,k}$ in $NN^{i}$ ($N_{w}=488$ in this study). Then, we define the deviation of weights $w^{i}_{l}$ from different $NN^{i}$ as
\begin{equation}
  \sigma_{l} = \frac{\sqrt{\sum_{i=1}^{210}\left(w^{i}_{l} - \mu_{l}\right)^{2}/210}}{\sum_{i=1}^{210}\left|w^{i}_{l}\right|/210},
\end{equation}
where $\mu_{l}= \sum_{i=1}^{210}w^{i}_{l}/210$. Note that the networks in NNTF were initialized with the pre-trained weights, and the networks in NN were initialized with the same random weights. Therefore, the deviation $\sigma_{l}$ indicates how a weight in a neural network is updated differently for different combinations of training and validation data samples (i.e., different $NN^{i}$). In the ideal case with an infinite amount of data samples, the deviation will approach zero as the network would converge to specific weights representing a unique generalized mapping of the real physics.

The distributions of the deviation $\sigma_{l}$ calculated from NNTF and NN are compared in figure~\ref{fig:histo}. A number of weights with large $\sigma_{l}$ are observed in NN, indicating that the updates of weights largely deviate depending on the data. This implies that the current data is insufficient for NN to learn a generalized mapping. On the other hand, the deviations are significantly reduced in NNTF, implying that the networks in NNTF are learning a better-generalized mapping compared to NN. This is because the networks in NNTF were initialized with the pre-trained weights containing the `approximate’ knowledge. These weights provide a better initial state for the weight optimization problem in the fine-tuning step.

\begin{figure}
  \centerline{\includegraphics[width = 0.5 \linewidth]{PDF_weights.eps}}
  \caption{Distributions of the deviations of weights $\sigma_{l}$ in NNTF (red) and NN (black).}
\label{fig:histo}
\end{figure}

\subsection{Effects of transfer learning on sensitivities of input surface statistics}\label{sec:sensitivity}
This section investigates how transfer learning affects the neural networks in learning the sensitivities of the input surface statistics.
The sensitivities of input surface statistics in predicting drag are obtained by calculating the derivatives of drag with respect to each input (surface statistical measure).
A similar analysis was employed by ~\citet{kim2020prediction}, for evaluating the influences of flow variables in predicting heat transfer characteristics.
We define the sensitivity of the input $I_{i}$ in predicting drag $\Delta \widetilde U^{+}$ (defined in equation~\ref{eq:u_ensemble}) as:
\begin{equation}
    S_{i} = \left\langle\left|\frac{\partial \Delta\widetilde U^{+}}{\partial I_{i}}\right|\right\rangle . \label{eq:sensitivity}
\end{equation}
Here, $\langle\cdot\rangle $ is the averaging operation along with the calculated derivatives from the 25 test DNS dataset. The derivative $\frac{\partial \Delta\widetilde U^{+}}{\partial I_{i}}$ can be analytically calculated as neural networks are composed of differentiable function compositions and matrix multiplications. However, more efficiently, it can be done by using an automatic differentiation algorithm. In this study, we used the automatic differentiation algorithm implemented in PyTorch~\citep{paszke2017automatic}.

\begin{figure}
  \centerline{\includegraphics[width = 0.6 \linewidth]{Sensitivity.eps}}
  \caption{Sensitivities of the input surface statistics in NNTF (red) and NN (black). The sensitivities are normalized with their total sum ($S_{i}/\sum_{i}S_{i}$). The input surface statistics are sorted in descending order from the largest to smallest sensitive surface statistics in NNTF.}
\label{fig:sen}
\end{figure}

Figure~\ref{fig:sen} shows the sensitivities of the input surface statistics when trained with NNTF and NN. It is observed that for both  neural networks, $I_{12}=ES_{z}\times Sk$,  $I_{10}=ES_{x}\times Sk$ and $I_{6}=por$ have the highest sensitivities. In other words, these are the most important surface measures when it comes to predicting the friction drag for the particular class of rough surfaces contained in the high-fidelity dataset.
These statistics are not explicitly included in the empirical models \eqref{EMP1}, \eqref{EMP2} and \eqref{EMP3}. Therefore, this indicates that the neural networks are learning the mapping of surface statistics to drag beyond the empirical models. However, the importance of $I_{12}, I_{10}$ and $I_{16}$ is partially  implied by previous studies.
The importance of the product of effective slope and skewness is partly implied in the study of \citet{forooghi2017toward}, as their empirical model (\ref{EMP2}) is nonlinear with respect to skewness and effective slope. The importance of porosity in predicting drag was recently emphasized by \cite{jouybari2021data}.
Although the result shown here is for a particular class of rough surfaces, it demonstrates the necessity of expanding the feature space of surface statistics for drag prediction to more complex nonlinear combinations of the most basic statistical moments.

To understand what influence transfer learning has on sensitivities, we can identify from figure~\ref{fig:sen}  the input statistics that are ranked differently by NN and NNTF in terms of importance. The surface statistics related to effective slopes and skewness, e.g., $I_{12}=ES_{z}\times Sk$, $I_{10}=ES_{x}\times Sk$, $I_{2}=Sk$, $I_{4}=ES_{x}$, and $I_{5}=ES_{z}$, are found to be more sensitive in NNTF. These higher sensitivities in NNTF are mainly due to the inclusion of the approximate knowledge about the drag dependencies on effective slope and skewness in the empirical correlations. Conversely, the sensitivity of the input defined as the square of skewness ($I_{17}=Sk \times Sk$) is smaller in NNTF compared to NN. This indicates that NNTF learns that the square of skewness compromises the information about the sign of skewness, which is known to be important for drag prediction~\citep{jelly2018reynolds}. In addition, NNTF reduces the sensitivities of statistics related to kurtosis ($I_{3}=kur$, $I_{13}=ES_{x}\times kur$, and $I_{11}=ES_{z}\times kur$). As kurtosis measures the outliers~\citep{decarlo1997meaning}, the reduction indicates that the effects of the outlier roughness heights are relaxed in NNTF.

\subsection{Approximate knowledge in the empirical correlations}\label{sec:approximate}
As briefly introduced in $\S$~\ref{sec:pre-train}, the empirical correlations contain different types of approximate knowledge. This is because the empirical correlations were fitted from different types of surfaces~\citep{chung2021predicting}; \ref{EMP1} was derived from data of a pipe with sinusoidal irregular roughness; \ref{EMP2} was developed from data of rough surfaces with different arrangements and size distributions of roughness elements and finally; \ref{EMP3} was constructed from a wide range of realistic rough surfaces, such as gravels and commercial steel pipes. As a consequence, each of the relations \ref{EMP1}, \ref{EMP2} and \ref{EMP3} depends on a different measure of roughness height ($Ra^{+}$, $k_{c}^{+}$, or $k_{rms}^{+}$) and different combinations of surface statistics ($Sk$ and/or $ES_{x}$).
Figure~\ref{fig:contour} shows contour levels of the drag predicted by  \eqref{EMP1}, \eqref{EMP2} and \eqref{EMP3} in a map spanned by skewness $Sk$ and effective slope $ES_x$. Note that the roughness heights for \ref{EMP1} to \ref{EMP3} ($Ra^{+}$, $k_{c}^{+}$, and $k_{rms}^{+}$) are chosen by their respective values that produce $\Delta U^{+}=-3.5$ at a common location ($ES_{x} = 0.3$ and $Sk = 0$) as in \cite{chung2021predicting}.

\begin{figure}
  \centerline{\includegraphics[width = 1.0 \linewidth]{contours.eps}}
  \caption{Contours of the roughness function calculated from \ref{EMP1}, \ref{EMP2}, \ref{EMP3}. The roughness heights $(Ra^{+},k_{c}^{+},k_{rms}^{+})=(0.234,0.474,1.106)$) are used for calculating the empirical correlations. The red {\color{red}$\times$} and {\color{red}$\circ$} markers indicate the test and fine-tuning DNS data samples, respectively. The black $\circ$ markers indicate the fitting data samples used for deriving the empirical datasets (data extracted from~\cite{chung2021predicting}).}
\label{fig:contour}
\end{figure}

The correlation \ref{EMP1} contains the approximate knowledge that drag tends to increase with increasing effective slope, while it does not contain any dependency of drag with respect to skewness. This is because \ref{EMP1} was derived from fitting data with a narrow regime of skewness as shown in the black $\circ$ markers in EMP1 of figure~\ref{fig:contour}. Oppositely, the correlation \ref{EMP3} contains the knowledge that drag tends to increase with increasing skewness, without the dependency on effective slope. This is because the majority of the fitting data for \ref{EMP3} (black $\circ$ markers in EMP3 of figure~\ref{fig:contour}) lies on an insensitive regime of effective slope -- between the sparse regime ($ES<0.3-0.6$) and dense regime ($ES>0.4-3.0$) -- with relatively small effects on drag~\citep{jimenez2004turbulent}. \ref{EMP2} shows a nonlinear dependence of drag on the skewness and the effective slope.
As a result, this model includes approximate knowledge about the drag dependency on the product of skewness and effective slope as also discussed in \S\ref{sec:sensitivity}.

As the empirical correlations are derived by fitting data from a particular set of surfaces, a naive extrapolation of them for predicting drag of other set of rough surfaces may lead to large errors. For example, if the correlations are directly applied to predict drag on the current test surfaces of the DNS dataset (red {\color{red}$\times$} markers in figure~\ref{fig:contour}), the average and maximum errors of \ref{EMP1}, \ref{EMP2}, and \ref{EMP3} are ($17\%$, $29\%$), ($8\%$, $17\%$) and ($25\%$, $44\%$), respectively. Note that the errors from \ref{EMP2} are the smallest due to the similar surface spaces between the current DNS data and the fitting data.
Interestingly, large errors of the empirical models do not necessarily lead to inaccurate fine-tuned neural networks. To demonstrate this, we  fine-tuned neural networks pre-trained only with one empirical correlation. We found that the average and maximum errors of the ensemble networks fine-tuned from \ref{EMP1}, \ref{EMP2}, and \ref{EMP3}, were ($6\%$, $18\%$), ($6\%$, $17\%$), and ($7\%$, $14\%$), respectively. Despite the high errors of the empirical correlations of \ref{EMP1} and \ref{EMP3}, their resulting fine-tuned networks are as effective as the fine-tuned network from \ref{EMP2}. This indicates that an empirical model for pre-training does not necessarily need to be accurate, it merely needs to provide an approximate knowledge of the dependencies of on surface statistics.
Moreover, NNTF, which is pre-trained with all of the approximate knowledge provided by the three empirical models, achieves a more favorable overall performance (figure~\ref{fig:ensemble}), compared to the fine-tuned networks pre-trained with any single EMP. This is because the knowledge contained in each empirical model contributes with information of drag dependencies.
This can be qualitatively shown by visualizing knowledge domains.

Figure~\ref{fig:knowledge} compares the knowledge domains of high-fidelity physics (DNS domain), approximate knowledge in the empirical models (EMP domain), and non-physical knowledge from the randomly initialized neural networks (Non-physics domain).
To visualize the knowledge domains, we extracted principal axes of the data samples composed of $(ES_{x}, Sk, \Delta U^{+})$ in each domain. Note that $\Delta U^{+}$ in the DNS, EMP, and Non-physics domains are obtained from simulations, empirical models, and randomly initialized neural networks, respectively. The DNS domain is composed of rough surfaces in the DNS dataset, while the EMP and Non-physics domains are composed of rough surfaces in the empirical dataset.
After computing the principal axes in each domain, ellipsoids that approximately bound the DNS, EMP, and Non-physics domains along their principal axes are visualized. As shown in figure~\ref{fig:knowledge}a, no particular correlation (or rather alignment) can be found between the Non-physics domain and the DNS domain. On the other hand, an alignment between the EMP domain and the DNS domain is observed in figure~\ref{fig:knowledge}b. Therefore, it is expected that neural networks with approximate knowledge from the EMP domain can more easily adapt to the DNS domain compared to those without any physical knowledge.

\begin{figure}
  \centerline{\includegraphics[width = 1.0 \linewidth]{knowledge_domains.eps}}
  \caption{Visualization of the knowledge domains of the DNS domain (red), EMP domain (blue), and Non-physics domain (green). Comparisons between the (a) DNS and Non-physics domains and the (b) DNS and EMP domains are shown.}
\label{fig:knowledge}
\end{figure}

We also quantified the alignment between the knowledge domains by calculating the angles between the principal axes. The angles between the principal axes in the EMP and DNS domains are computed as ($16.5^{\circ}$, $18.1^{\circ}$, $8.9^{\circ}$), while the angles between the principal axes in the Non-physics and DNS domains are computed as ($120.9^{\circ}$, $119.9^{\circ}$, $12.3^{\circ}$).
Accordingly, the aligned `approximate knowledge' from EMPs assists the networks to adapt to the DNS domain; thus, a better performance of neural networks can be achieved by transfer learning of empirical correlations.

\section{Conclusions}\label{sec:conclusion}
We have developed a transfer learning framework to learn the Hama roughness function from a limited number of DNS data samples. The framework is composed of two steps: (1) pre-training step and (2) fine-tuning step. In the pre-training step, neural networks learn an `approximate' knowledge from empirical correlations. In the fine-tuning step, the networks are fine-tuned using a small DNS dataset. Neural networks trained with transfer learning have shown a significant improvement in predicting the roughness function on test data samples not included in the fine-tuning step. This is because the `approximate' knowledge in empirical correlations is aligned to high fidelity physics, which assists neural networks to learn better-generalized weights. In addition, a sensitivity analysis showed that the neural networks with transfer-learning were clearly emphasizing the importance of certain input surface statistics (effective slopes, skewness and porosity) in predicting drag. These extracted statistics are in general  good agreement with what has been reported by other investigations, but also highlight that nonlinear functions of several surfaces statistics could provide high correlation with drag.

We have shown that the prediction performance is enhanced when the `approximate' knowledge of the empirical dataset is well aligned to the high-fidelity dataset. Therefore, it is advantageous to employ empirical correlations developed for classes of rough surfaces that are similar to one's high-fidelity dataset.
The current NNTF, trained with irregular homogeneous rough surfaces, would not be effective when predicting drag on regular rough surfaces (e.g.,~cuboids, bars, etc.) or inhomogeneous rough surfaces. Similarly, predicting values of roughness function beyond training is also expected to not be effective due to the general limitations of data-driven methods in extrapolating.

To further increase the generalization ability of neural networks, ultimately, the dataset has to be expanded. A large database would also help to study classes of rough surfaces with undiscovered empirical correlations, where pre-training is currently not possible. Accordingly, a collective community effort is needed to construct and train Big Data of flows over many classes of rough surfaces in different flow conditions. One step in this direction, is the online database \href{http://roughnessdatabase.org/}{Roughness database}~\citep{roughnessdatabase} currently under construction. Moreover, new primary surface statistics that strongly affect drag may be discovered in the future by extending the sensitivity analysis in \S\ref{sec:sensitivity} with a larger surface parameter space and Big Data.
Also, expansion of the network structure would enable one to fully leverage information inside large databases. The structure of the neural network in this study is very simple, as our focus has been to show that transfer learning can improve the performance of a given network. Thus, this network structure should not be considered as the optimal model for predicting drag. In addition to expanding the database, further optimizing the neural networks will also improve the performance for drag prediction.

Finally, while we have considered transport of momentum (i.e., drag coefficient), empirical relations are widespread also for characterizing transport of energy (Nusselt number) and mass (Sherwood number). Therefore, the proposed framework can also be applied to other engineering applications where only a limited amount of high-fidelity data is available, but a significant amount of knowledge has been accumulated.

\section*{Acknowledgements}{This work was supported by the Swedish Energy Agency under Grant number 51554-1, Swedish Foundation for Strategic Research (FFL15-001), and Friedrich und Elisabeth Boysen-Foundation (BOY-151). SL and SB would like to thank Dr. Oscar Beijbom for the helpful discussions regarding the neural networks.}

\section*{Declaration of interests}{The authors report no conflict of interest.}


\section{Introduction}
\input{intro.tex}

\input{section2.tex}

\input{section3.tex}

\input{section4.tex}

\section{Conclusions}
\label{sec:conclu}

In this letter, a detailed analysis of the $K \to \pi a$ decay has been presented, in view of the recent NA62 measurement and
the foreseen updates from the KOTO experiment. Assuming flavor and CP conserving ALP couplings with fermions, the dominant contribution
to the $K \to \pi a$ decays arises from the penguin diagrams, manly proportional to the $c_t$ coupling. NA62 and E949 experiments
bound $c_t$ to be smaller than $6 \times 10^{-4}$ in most of the allowed $m_a$ region, for the chosen reference value of the PQ symmetry
breaking scale $f_a=1$ TeV. Expected KOTO results can provide comparable bounds on the $c_t$ coupling. Subdominant tree--level diagrams
can, however, contribute to the $K \to \pi a$ decay process. The tree-level amplitudes contributing to the $K^+ \to \pi^+ a$ and
$K^0_L \to \pi^0 a$ decays have been derived in this letter, following the Lepage--Brodsky \cite{Lepage:1980fj,Szczepaniak:1990dt} technique.
Assuming all ALP couplings with fermions to be vanishing, but $(c_u,c_s)$, independent limits on the the light quark--ALP couplings
$(c_u,c_s) \sim 0.05$ can be obtained in the allowed $m_a$ mass range, from NA62 and E949 experiments. KOTO experiment, conversely, will
not provide competitive bounds on the $(c_d,c_s)$ couplings, due to the SM CP violating suppression of the $K^0_L$ channel.

If a universal ALP coupling to fermions, $c_{a\Phi}$, is assumed, the results presented in this letter confirm and refine previous analysis
\cite{Izaguirre:2016dfi,Dolan:2017osp,Gavela:2019wzg}, where only penguin contributions were considered. However, in general MFV ALP
frameworks it may not be unconceivable to assign different, but flavor universal, PQ charges to the up and down quark sectors. In this case
the simpler ''penguin dominated'' result may be significantly modified. The subdominant tree-level contribution proportional to $c_\downarrow$
can strongly interfere with the measurement of $c_\uparrow$, introducing an uncertainty band of roughly one order of magnitude. Similar
analysis should be extended to all meson decays in ALP\cite{newUS}.

\section{Acknowledgements}

The authors thank O. Sumensari for useful discussions in the early stage of this paper. A.G. and S.R. acknowledge support from the European
Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curiegrant agreements 690575 (RISE InvisiblesPlus) and
674896 (ITN ELUSIVES). This project has  received funding/support from the European Union’s Horizon 2020 research and innovation programme
under the Marie Sklodowska-Curie grant agreement No 860881-HIDDEN.

\bibliographystyle{unsrt}

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

Light pseudoscalar particles naturally arise in many extensions of the Standard Model (SM) of particle physics. In
particular they are a common feature of any BSM model endowed with (at least) a global $U(1)$ symmetry spontaneously
broken at some high scale $f_a \gg v$. Small breaking terms of the global symmetry are then necessary in order to provide
a mass term, $m_a \ll f_a$, to the generated pseudo Nambu-Goldstone boson (pNGB). Therefore, it may be not inconceivable
that the first hint of new physics at (or above) the TeV scale could be the discovery of a light pseudoscalar state.

Sharing a common nature with the QCD axion \cite{Peccei:1977hh,Wilczek:1977pj,Weinberg:1977ma}, these class of pNGBs
are generically referred to as Axion-Like Particles (ALPs). The key difference between the QCD axion and a generic
ALP can be summarized in the well-known constraints \cite{Weinberg:1977ma}:
\begin{eqnarray}
m_a f_a \approx m_\pi f_\pi
\label{mafarelation}
\end{eqnarray}
that bounds the QCD axion mass and the $U(1)_{PQ}$ symmetry breaking scale via QCD instanton effects. Present bounds
on the reference QCD invisible-axion models, like the DNSZ and KVSZ ones \cite{Kim:1979if,Shifman:1979if,Zhitnitsky:1980tq,
Dine:1981rt}, force the axion mass to be typically in the sub-eV range with the symmetry breaking scale $f_a \gtrsim 10^{11}$
GeV. In a generic ALP framework, instead, one can assume the ALP mass being determined by some unspecified UV physics,
besides the usual QCD anomalous contribution, and then relation in Eq.~(\ref{mafarelation}) has not to be enforced. Consequently,
the ALP mass and the $U(1)_{PQ}$ breaking scale $f_a$ can be taken as independent parameters and in the range phenomenologically
of interest at present or near future colliders. In this context, a generic ALP can be seen as the generalization of non
fine-tuned axion models, and permit to look without prejudice to all the possible light pseudo--scalar signatures in Cosmology,
Astro--Particle and Collider physics.

The ALP parameter space has been intensively explored in several terrestrial facilities, covering a wide energy range~
\cite{Mimasu:2014nea,Jaeckel:2015jla,Bauer:2017ris,Brivio:2017ije,Alonso-Alvarez:2018irt,Harland-Lang:2019zur,
Baldenegro:2018hng,MartinCamalich:2020dfe}, as well as by many astrophysical and cosmological probes~\cite{Cadamuro:2011fd,
Millea:2015qra,DiLuzio:2016sbl}. The synergy of these experimental searches allows to access several orders of magnitude in
ALP masses and couplings, cf.~e.g.~Ref.~\cite{Irastorza:2018dyq} and references therein. While astrophysics and cosmology impose
severe constraints on ALPs in the sub-KeV mass range, the most efficient probes of weakly-coupled particles in the MeV-GeV
range come from experiments acting on the precision frontier~\cite{Essig:2013lka}. Fixed-target facilities such as
E949~\cite{Artamonov:2009sz,Artamonov:2008qb,Adler:2008zza}, NA62~\cite{CortinaGil:2020fcx,CortinaGil:2021nts} and KOTO~\cite{Ahn:2018mvc}
and the proposed SHiP \cite{Alekhin:2015byh} and DUNE \cite{Kelly:2020dda} experiments can be very efficient to constrain long-lived
particles. Furthermore, the rich ongoing research program in the $B$-physics experiments at LHCb~\cite{Aaij:2015tna,Aaij:2016qsm}
and the $B$-factories~\cite{Masso:1995tw,Bevan:2014iga,Izaguirre:2016dfi,Dolan:2017osp,Kou:2018nap,CidVidal:2018blh,
deNiverville:2018hrc,Gavela:2019wzg,MartinCamalich:2020dfe} offers several possibilities to probe yet unexplored ALP couplings.

The main goal of this letter, is the detailed analysis of the $K \to \pi a$ decays, in view of the recent NA62 updated
measurement and the foreseen updates of the KOTO results. While one-loop penguin contributions are usually considered
\cite{Izaguirre:2016dfi,Dolan:2017osp,Gavela:2019wzg}, the tree-level process contributing to the decay has most often
been overlooked in literature typically for two reasons. First of all one expects the penguin diagrams to dominate,
despite the loop suppression, being proportional to the mass of the virtual top quark running in the loop. However,
in the $K\to \pi a$ case, one has to take into account that the top--loop contribution is CKM suppressed compared to the
tree-level roughly by a factor $\lambda^4$, and this can partially compensate for the top mass enhancement. The second
reason is because the tree-level contributions show a much more complicated hadronization structure, being the ALP
emitted ``inside'' the initial or final meson and so require a dedicated treatment \cite{Lepage:1980fj,Szczepaniak:1990dt}.
An alternative approach for calculating $K \to \pi a$ decays using chiral perturbation theory can be found, for example,
in \cite{Georgi:1986df,Bauer:2021wjo}.

The letter is organized as follows. In Sec.~\ref{sec:effective} the effective Lagrangian describing the flavor-conserving
interactions between the ALP and SM fermions up to dimension five is introduced. In Sec.~\ref{sec:hadronization} the
tree-level contribution to the $K\to \pi a$ decay is going to be thoroughly discussed. Then, the one-loop penguin
contributions are shortly reviewed. Finally, in Sec.~\ref{sec:pheno} the interplay between the tree-level and one-loop
contributions to the $K \to \pi a$ Branching Ratio is thoroughly discussed. Concluding remarks are deferred to
Sec.~\ref{sec:conclu}.

\section{Effective ALP-SM Fermion Lagrangian}
\label{sec:effective}

The most general effective Lagrangian describing ALP interactions with SM quarks, including operators up to dimension five, reads:
\begin{eqnarray}
\delta \mathcal{L}^a_{\mathrm{eff}}
 &=&  \frac{\partial_\mu a}{f_a} \left[\overline{U}\, \gamma^\mu \left(C^{(u)}_{L} P_L  \,+\,C^{(u)}_{R} P_R \right) \,U
                            \,+\,\overline{D}\, \gamma^\mu \left(C^{(d)}_{L} P_L \,+\,C^{(d)}_{R} P_R\right) \,D \right]

\label{eq:ALP_SMquarks}
\end{eqnarray}
where $f_a$ is the $U(1)_{PQ}$ symmetry breaking scale, $U$ and $D$ the SM up and down flavour triplets and $C_{L,R}$ are general
$3\times 3$ hermitian matrices. One can, however, heavily reduce the number of independent parameters imposing that the only source
of flavor--violation in the model arises through the SM Yukawa couplings. Therefore, once the Minimal Flavor Violation (MFV) ansatz
is assumed, the ALP--quarks Lagrangian becomes
\begin{eqnarray}
\delta \mathcal{L}^{a,MFV}_{\mathrm{eff}} &=&
        -\frac{\partial_\mu a}{ 2 f_a} \sum_{i=quarks} \hspace{-0.25cm} c_i\, \overline{\psi}_i \gamma^\mu \gamma_5 \,\psi_i
   \,=\, i \frac{a}{f_a} \sum_{i=quarks} \hspace{-0.25cm} c_i m_i \, \overline{\psi}_i \gamma_5 \,\psi_i
\label{eq:ALP_SMMFV}
\end{eqnarray}
and depends only on six independent flavor diagonal couplings, $c_i$, once the vector fermionic current conservation is implied. Therefore,
in the MFV framework, flavor--violating ALP couplings can be  generated only at loop level, and is proportional to the SM CKM mixings.
To further reduce the number of independent parameters, one can additionally assume universal couplings for the up and down sectors, in the
following denoted as $c_\uparrow$ and $c_\downarrow$, respectively. It may not be unconceivable, in fact, that an hypothetical UV--complete
model provides different universal PQ charges to the up and down quark sectors \cite{Brivio:2017sdm,Merlo:2017sun,Alonso-Gonzalez:2018vpc}.
Finally, in the most constrained scenario, one can assume a unique ALP-fermion coupling, often denoted as $c_{a\Phi}$ in the literature,
as it originates from the dimension five ALP-Higgs operator
\begin{eqnarray}
\mathcal{O}_{a\Phi} \equiv i\,\frac{\partial_\mu a}{f_a} \left( \Phi \stackrel{\leftrightarrow}{D}_\mu \Phi\right),
\end{eqnarray}
once the Higgs field is accordingly redefined \cite{Georgi:1986kr}. These additional assumptions will become useful in simplifying
the phenomenological analysis and their implication will be discussed in Sec.~\ref{sec:pheno}. A general discussion on tree--level
flavor--violating ALPs couplings to fermions can be found in \cite{MartinCamalich:2020dfe}, while for a recent analysis on
CP--violating ALP couplings to fermions one is referred, for example, to \cite{DiLuzio:2020oah}.

It might be useful, for simplifying intermediate calculations and explicitly showing the mass dependence of ALP-fermion couplings, to write
the effective Lagrangian of Eq.~(\ref{eq:ALP_SMMFV}) in the ``Yukawa'' basis instead of the ``derivative'' one. The two versions of the
effective Lagrangian are equivalent up to operators of $O(1/f_a^2)$.

\section{Meson Hadronization in Flavor Changing Processes}
\label{sec:hadronization}

Using the effective Lagrangian, implemented with the flavor conserving assumption, of Eq.~\ref{eq:ALP_SMMFV} one can calculate
the hadronic decay rates of mesons in ALPs. In the following, due to their experimental relevance, the pseudo-scalar meson
hadronic decays:
\begin{eqnarray}
K^+ \rightarrow \pi^+ \, a \qquad {\rm and} \qquad K^0_L \rightarrow \pi^0 \, a
\label{invKPidecay}
\end{eqnarray}
will be mainly considered, with the ALP sufficiently long-living to escape the detector without decaying or decaying into invisible
channels. In such a case the only possible ALP signature is its missing energy/momentum.\footnote{This study will be generalized to
a wider class of scalar meson decay in a subsequent work \cite{newUS}.} In the following, with $M_{K^{+,0}}, P_{K^{+,0}}$ and
$M_{\pi^{+,0}}, P_{\pi^{+,0}}$ the mass and 4--momentum of the $K^{+,0}$ and $\pi^{+,0}$ mesons will be denoted, respectively,
while the ALP mass and 4--momentum will be indicated with $m_a$ and $k_a$.

The processes at hand receive contributions both from tree--level and one--loop diagrams. Typically, for this kind of processes,
the one--loop penguin diagrams with a top virtual exchange give the dominant contribution, as the $m_t/m_{q,ext}$ ratio largely
compensate the loop suppression factor. The tree--level and one--loop decay channels involve, however, quite different
hadronization structures that have to be properly taken into account for correctly comparing the relative contributions to the $K$
decay rates. In the following subsections it is briefly illustrated how to handle both contributions in a simple way, applying
different hadronization methods \cite{Lepage:1980fj,Szczepaniak:1990dt,Carrasco:2016kpy}.

\subsection{The tree-level s-channel process}

Charged pseudo-scalar meson decays proceed through the s--channel tree-level diagrams of Fig.~\ref{fig:tree_charged}. Here
only the diagrams where the ALP is emitted from the $K^+$ meson are shown, the ones where the ALP is emitted from the $\pi^+$ follow
straightforwardly. The tree-level diagram with the ALP emitted from the $W^+$ internal line automatically vanishes, being the
WW-ALP coupling proportional to the fully antisymmetric 4D tensor.
The hadronization of the s--channel can be done following the Lepage--Brodsky technique \cite{Lepage:1980fj,Szczepaniak:1990dt}.
Leptonic pseudo--scalar decays, $\mathcal{P} \to \ell \,\nu_\ell \, a$, deserve a similar treatment and have been derived in
\cite{Aditya:2012ay,Hazard:2016fnc}. Let's recall briefly the notation.

\begin{figure}[!ht]
\centering
\includegraphics[scale=0.16]{PICTURES/s_emission_1.jpg}\hspace{2 cm}\includegraphics[scale=0.16]{PICTURES/u_emission_1.jpg}
\caption{Tree level contributions to the $K^+ \rightarrow \pi^+ \, a$ amplitude, with the ALP emitted from the $K^+$ meson. Diagrams
where the ALP is emitted from the $\pi^+$ quarks are straightforward.
\label{fig:tree_charged}}
\end{figure}

The parent meson $K^+$ constituent quarks annihilate into a virtual W boson that then produces the final $\pi^+$ meson partons. The hadronic
process can be factorized as $\bra{\pi^+}\bar{u}\,\Gamma_{(\pi)}\,d\ket{0}\bra{0}\bar{s}\,\Gamma_{(K)}\,u\ket{K^+}$ with the operator
insertion $\Gamma_{(\pi)}\otimes\Gamma_{(K)}$ being $\gamma^\mu P_L \otimes \Gamma_{\mu}$ or $\Gamma^{\prime}_{\mu} \otimes \gamma^\mu P_L$
depending if the ALP is emitted by the initial or final mesons, with
\begin{eqnarray}
\Gamma_\mu \hspace{-0.15cm} &=& \hspace{-0.15cm} \frac{4 G_F}{\sqrt{2}} V^*_{us}V_{ud}
          \left(\frac{c_s\,m_s}{f_a}\,\gamma_5\,\frac{\slashed{k}_a-\slashed{p}_{\bar{s}}+m_s}{m_a^2-2 k_a \cdot p_{\bar{s}}}\,\gamma_\mu P_L
   \,-\, \frac{c_u\,m_u}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{k}_a-\slashed{p}_u-m_u}{m_a^2 - 2 k_a \cdot p_u}\,\gamma_5 \right)
          \label{eq:GammaK} \\
\Gamma'_\mu \hspace{-0.15cm} &=& \hspace{-0.15cm}  \frac{4 G_F}{\sqrt{2}} V^*_{us}V_{ud}
          \left(\frac{c_u\,m_u}{f_a}\,\gamma_5 \,\frac{\slashed{k}_a+\slashed{p}'_u+m_u}{m_a^2-2 k_a\cdot p'_u}\, \gamma_\mu P_L\,
   \,-\, \frac{c_d\,m_d}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{k_a}+\slashed{p}'_{\bar{d}}-m_d}{m_a^2-2 k_a \cdot p'_{\bar{d}}}\,\gamma_5\right).
\label{eq:GammaPi}
\end{eqnarray}
In Eqs.~(\ref{eq:GammaK}) and (\ref{eq:GammaPi}) with $p_s,p_u$ and $p'_d,p'_u$ the quark momenta of the initial and final meson are denoted,
respectively. In deriving these Equations the ``Yukawa'' basis for the ALP-SM quark couplings has been explicitly used, providing a simpler
$\gamma$ structure and showing explicitly the mass dependence of the quark-ALP couplings.

The vector and axial matrix elements can be parameterized in terms of the meson decay constants $f_K$ and $f_\pi$ as:
\begin{eqnarray}
& &\bra{0}\bar{s}\,\gamma^\mu \,\gamma_5 \,u \ket{K^+} = i f_K P_K^\mu \qquad , \qquad \bra{0}\bar{s}\,\gamma^\mu\,u \ket{K^+} = 0 \\
& &\bra{0}\bar{d}\,\gamma^\mu \,\gamma_5 \,u \ket{\pi^+} = i f_\pi P_\pi^\mu \,\,\,\qquad ,\qquad \,\bra{0}\bar{d}\,\gamma^\mu\,u\ket{\pi^+} = 0 \,.
\end{eqnarray}
To compute the $\bra{0} \bar{s}\, \Gamma_\mu \,u\,\ket{K^+}$ and $\bra{\pi^+}\bar{u}\, \Gamma'_\mu\,d\,\ket{0}$ hadronic matrix elements,
one has to assume a model for describing the effective quark--antiquark distribution inside the meson emitting the ALP. Following
\cite{Lepage:1980fj,Szczepaniak:1990dt,Aditya:2012ay}, the ground state of a meson $M$ is parameterized with the wave--function
\begin{equation}
\Psi_M(x)=\frac{1}{4}\phi_M(x)\gamma^5(\s{P}_M + g_M(x) M_M),
\label{eq:wave_parameter}
\end{equation}
where ${P}_M$ and $M_M$ denote the momentum and the mass of the meson emitting the ALP\footnote{In the literature the functions $g_M(x)$
are conventionally assumed to be constants, $g_M \in [0,1]$.}. In Eq.~(\ref{eq:wave_parameter}), with $x$ one typically denotes the
fraction of the momentum carried by the heaviest quark in the meson. The function $\phi_M(x)$ describes the meson quarks momenta
distribution, that for heavy and light mesons reads, respectively:
\begin{eqnarray}
\phi_H(x) \propto  \left[\frac{\xi^2}{1-x}+\frac{1}{x}-1\right]^{-2} \qquad , \qquad \phi_L(x) \propto  x(1-x) \,,
\end{eqnarray}
with the normalization fixed such that:
\begin{equation}
\int_0^1 dx\,\phi_M(x) = 1.
\label{eq:normalization}
\end{equation}
The parameter $\xi$ in $\phi_H(x)$ is a small parameter typically of $O(m_q/m_Q)$, being $q$ and $Q$ the light and heavy quark in the
meson. The mass function $g_M(x)$ is usually taken to be a constant varying from $g_H(x) \approx 1$ and $g_L(x) \ll 1$ for a heavy or a
light meson, respectively. The hadronic matrix element can then be obtained by integrating over the momentum fraction $x$ the trace of
the $\Gamma^\mu$ amplitude over the meson wave--function $\Psi_M(x)$:
\begin{equation}
\bra{0} \bar{Q}\,\Gamma^{\mu}\,q \ket{M} \equiv i f_M \int_0^1 dx\,\mathrm{Tr}\left[\Gamma^{\mu} \Psi_M (x)\right]\,.
\label{eq:mesonic_parametrization}
\end{equation}
In Eqs.~(\ref{eq:wave_parameter}-\ref{eq:mesonic_parametrization}), a slightly different notation with respect to the cited literature
is used, in particular, as for weak decays colors are left unchanged, all the color matrices/traces have been removed from scratch.
Moreover, the functions $\Phi_M(x)$ have been normalized to one, in such a way that in Eq.~(\ref{eq:mesonic_parametrization}) the mesonic
form factor can be explicitly factorized.

By inserting Eqs.~(\ref{eq:GammaK}-\ref{eq:GammaPi}) into Eq.~(\ref{eq:mesonic_parametrization}), and making the following assignments for
the initial and final quark momenta:
\begin{eqnarray}
p_{\bar{s}}  &=& x P_K \qquad , \qquad p_u  = (1-x) P_K \nonumber \\
p'_{\bar{d}} &=& x P_\pi \qquad \,, \, \qquad p'_u = (1-x) P_\pi \nonumber
\end{eqnarray}
one obtains the following decay amplitudes for the $K^+$-ALP and $\pi^+$-ALP emission processes:
\begin{eqnarray}
\mathcal{M}_{K^+} &=& \frac{G_F}{\sqrt{2}} (V^*_{us} V_{ud})\,f_K\, f_\pi \, (k_a\cdot P_\pi) \frac{M_K}{f_a} \,\times\, \label{eq:MKALP} \\
     & & \hspace{0.75cm} \times \int^1_0 \left\{ \frac{c_s\,m_s\,\theta(x-\delta^K_a)}{m_a^2-2\,x\,k_a\cdot P_K} -
         \frac{c_u\,m_u\,\theta(1-x-\delta^K_a)}{m_a^2-2\,(1-x)\,k_a\cdot P_K} \right\} \,\phi_K(x)\,g_K(x)\, dx \nonumber \\
\mathcal{M}_{\pi^+} &=& \frac{G_F}{\sqrt{2}} (V^*_{us} V_{ud})\,f_K\, f_\pi\, (k_a\cdot P_K) \frac{M_\pi}{f_a} \,\times\, \label{eq:MPiALP}\\
     & & \hspace{0.75cm} \times \int^1_0 \left\{ \frac{c_d\,m_d\,\theta(x-\delta^\pi_a)}{m_a^2-2\,x\,k_a\cdot P_\pi} -
        \frac{c_u\,m_u\,\theta(1-x-\delta^\pi_a)}{m_a^2-2\,(1-x)\,k_a\cdot P_\pi} \right\} \,\phi_\pi(x)\,g_\pi(x)\,dx \nonumber
\end{eqnarray}
with $\delta^M_a = m_a^2/(2 k\cdot P_M)$ an explicit cutoff introduced in the fractional momentum to remove the unphysical singularities
appearing in the integrals.
The result in Eq.~(\ref{eq:MKALP}) is in agreement with the decay amplitude for $B^\pm \to \ell \,\bar{\nu}_\ell \, a$ calculated in
\cite{Aditya:2012ay} once $c_u=c_s=1$ is assumed, the pion hadronic current is replaced by the leptonic one, and K quantities replaced by
the corresponding B meson ones. Eqs.~(\ref{eq:MKALP}) and (\ref{eq:MPiALP}) represent the main result of this section.

Few comments are in order. To numerically evaluate the $K^+\to\pi^+\, a$ branching ratio one has to assume a specific form of the hadronic
functions $\phi_K(x)$, $g_K(x)$, $\phi_\pi(x)$ and $g_\pi(x)$, and assign a low energy meaning to the quark masses. This, inevitably,
introduces some model dependence in the calculation. To have an order of magnitude estimate of the $\mathcal{M}_{K^+}$ amplitude, one
can consider the two extreme cases and treat the K-meson either as a light meson (i.e. assuming an exact global $SU(3)$ symmetry) or as
an heavy meson (i.e. $m_s \gg m_u,m_d$). In the light meson approximation, substituting the light quarks with the corresponding partons,
i.e. $\hat{m}_s=\hat{m}_u=M_K/2$, one obtains, for a massless ALP:
\begin{eqnarray}
\mathcal{M}^L_{K^+} \approx - \frac{3\, G_F \,f_K\,f_\pi}{4\sqrt{2}} (V^*_{us} V_{ud}) \frac{M^2_K}{f_a}\,g_K\,\left(c_s-c_u\right) \,.
\label{eq:MKLALPapprox}
\end{eqnarray}
Conversely, in the heavy meson approximation, one can assume $\hat{m}_u = \xi M_K$ and $\hat{m}_s = (1-\xi) M_K$, with $\xi = m_u/m_s$.
Moreover, approximating $\phi_K(x) \approx \delta(1-x-\xi)$ as suggested in \cite{Bhattacharya:2018msv}, one obtains, in the $m_a=0$ limit:
\begin{eqnarray}
\mathcal{M}^H_{K^+} \approx  - \frac{G_F \,f_K\,f_\pi}{2\sqrt{2}} (V^*_{us} V_{ud}) \frac{M^2_K}{f_a}\,g_K\,\left(c_s-c_u\right) \,.
\label{eq:MKHALPapprox}
\end{eqnarray}
From the approximate formulas of Eq.~(\ref{eq:MKLALPapprox}) and (\ref{eq:MKHALPapprox}) one can estimate roughly the order of magnitude
of the uncertainties introduced in the calculation by the hadronization procedure for the K meson. To reproduce the numerical results of
the following section, an intermediate approach will be, instead, considered: the heavy meson function will be used, with the two partons
defined as $\hat{m}_u = m_u + \Lambda$ and $\hat{m}_s = m_s +\Lambda$ with $\Lambda = (M_K-m_u-m_s)/2$ a parameter of order $\Lambda_{QCD}$.
Conversely, for the estimation of the $\mathcal{M}_\pi$ amplitude, one can safely assume to parametrize the pion using the light meson
wave-function. If one take $g_\pi(x) \approx 0$, as customarily suggested in the literature the pion contribution automatically vanishes.
A conservative estimate can however be obtained by setting, for example, $g_\pi/g_K \approx M_\pi/M_K$, which predicts the following upper
bound to the ratio
\begin{eqnarray}
R_{\pi K} = \left| \frac{\mathcal{M}_{\pi^+}}{\mathcal{M}_{K^+}} \right| \lesssim \left(\frac{M_\pi}{M_K}\right)^3 \simeq 1.\times 10^{-2} \nonumber \,.
\label{MpiKapprox}
\end{eqnarray}
For this reason, even in the numerical calculation one can neglect the ALP-$\pi$ emission as expected on a general ground, once same
order ALP couplings to $u,d$ and $s$ quarks are assumed.

Finally, from Eqs.~(\ref{eq:MKLALPapprox}) and (\ref{eq:MKHALPapprox}) it appears evident the presence of an ``accidental''
cancelation if $c_s=c_u$ is assumed. This cancellation is still partially at work even when the full $\phi_K(x)$ is used and indicates a
possible underestimation of the $\mathcal{M}_K$ amplitude (and consequently on the ALP-quark coupling limits) in a ``universal'' ALP--SM
quark coupling scenario compared to the general case.

\subsection{The tree-level t-channel process}

Neutral pseudo-scalar meson decays proceed through the t--channel tree-level diagrams of Fig.~\ref{fig:tree_neutral}. Here only the
diagrams where the ALP is emitted from the $K^0$ meson are shown, the ones where the ALP is emitted from the $\pi^0$ follow
straightforwardly. The tree-level diagram with the ALP emitted from the $W^+$ internal line automatically vanishes, as for the
s--channel case. The hadronization of the t--channel can be done along the lines depicted in \cite{Lepage:1980fj,Szczepaniak:1990dt}.

\begin{figure}[!ht]
\centering
\includegraphics[scale=0.16]{PICTURES/ns_neut_2.jpg}\hspace{2 cm}\includegraphics[scale=0.16]{PICTURES/nd_neut_2.jpg}
\caption{Tree level contribution to the amplitude for the $K^0 \rightarrow \pi^0 \, a$ decay, with the ALP emitted from the $K^0$
meson. Diagrams where the ALP is emitted from the $\pi^0$ quarks are straightforward. Similar diagrams can be depicted for the
CP conjugate process $\bar{K}^0 \rightarrow \pi^0 \, a$.
\label{fig:tree_neutral}}
\end{figure}

\noindent
Using the conventions defined in Eqs.~(\ref{eq:wave_parameter}) and (\ref{eq:normalization}), the neutral hadronic matrix element for the
$K^0$ transition reads:
\begin{equation}
\bra{\pi^0}\,\Gamma_1 \otimes \Gamma_2 \ket{K^0} \equiv
   - \frac{f_K f_\pi}{\sqrt{2}} \int_0^1 dx \int_0^1 dy\,\mathrm{Tr}\left[\Psi_\pi(y)\Gamma_{(1)}\Psi_K(x)\Gamma_{(2)} \right]\,,
\label{eq:tchannel_parametrization}
\end{equation}
with the $1/\sqrt{2}$ factor taking care of the Clebsch--Gordon suppression of the decay into the neutral $\pi$ meson.

The operator insertion $\Gamma_{(1)}\otimes\Gamma_{(2)}$ is defined as $\gamma_\mu P_L \otimes \Gamma^{\mu}_{(\bar{q})}$ or
$\Gamma^{\mu}_{(q)} \otimes \gamma_\mu P_L$ depending if the ALP is emitted from the $\bar{q}$ anti--quark (i.e. $\bar{s}$ for the
$K^0$ and $\bar{u}$ for the $\pi^0$ meson) or from quark $q$ (i.e. the $d$ quark for the $K^0$ and the $u$ quark for the $\pi^0$ meson),
respectively, with
\begin{eqnarray}
\Gamma^\mu_{(\bar{q})}\hspace{-0.15cm} &=& \hspace{-0.15cm} - \frac{4 G_F}{\sqrt{2}} V^*_{us}V_{ud}
     \left(\frac{c_s\,m_s}{f_a}\,\gamma_5\,\frac{\slashed{k}_a-\slashed{p}_{\bar{s}}+m_s}{m_a^2-2 k_a \cdot p_{\bar{s}}}\,\gamma_\mu P_L \,-\,
     \frac{c_u\,m_u}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{k}_a+\slashed{p}'_{\bar{u}}-m_u}{m_a^2-2 k_a\cdot p'_{\bar{u}}}\,\gamma_5\right)
     \hspace{0.4cm}
\label{eq:Gamma_qb} \\
\Gamma^\mu_{(q)} \hspace{-0.15cm} &=& \hspace{-0.15cm} - \frac{4 G_F}{\sqrt{2}} V^*_{us}V_{ud}
          \left(\frac{c_d\,m_d}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{p}_d-\slashed{k}_a+m_d}{m_a^2-2 k_a\cdot p_d}\,\gamma_5 \,+\,
          \frac{c_u\,m_u}{f_a}\,\gamma_5\,\frac{\slashed{k}_a+\slashed{p}'_u+m_u}{m_a^2-2 k_a \cdot p'_u}\,\gamma_\mu P_L\right)
\label{eq:Gamma_q}
\end{eqnarray}
The $\bar{K}^0\to\pi^0\,a$ decay amplitude can be obtained similarly:
\begin{equation}
\bra{\pi^0}\,\bar{\Gamma}_1 \otimes \bar{\Gamma}_2 \ket{\bar{K}^0} \equiv
   -\frac{f_K f_\pi}{\sqrt{2}} \int_0^1 dx \int_0^1 dy\,\mathrm{Tr}\left[\Psi_\pi(y)\bar{\Gamma}_{(1)}\Psi_K(x)\bar{\Gamma}_{(2)} \right]\,,
\label{eq:tchannel_bparametrization}
\end{equation}
with the operator insertions $\bar{\Gamma}_{(1)}\otimes\bar{\Gamma}_{(2)}$ being $\gamma_\mu P_L \otimes \bar{\Gamma}^{\mu}_{(q)}$ or
$\bar{\Gamma}^{\mu}_{(\bar{q})} \otimes \gamma_\mu P_L$ with
\begin{eqnarray}
\bar{\Gamma}^\mu_{(q)}\hspace{-0.2cm} &=& \hspace{-0.15cm} - \frac{4 G_F}{\sqrt{2}} V_{us}V^*_{ud}
     \left(\frac{c_s\,m_s}{f_a}\,\gamma_5\,\frac{\slashed{p}_s-\slashed{k}_a+m_s}{m_a^2-2 k_a \cdot p_s}\,\gamma_\mu P_L \,+\,
     \frac{c_u\,m_u}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{k}_a+\slashed{p}'_u+m_u}{m_a^2-2 k_a\cdot p'_u}\,\gamma_5\right)
     \hspace{0.5cm}  \label{eq:bGamma_qb} \\
\bar{\Gamma}^\mu_{(\bar{q})} \hspace{-0.2cm} &=& \hspace{-0.15cm} - \frac{4 G_F}{\sqrt{2}} V_{us}V^*_{ud}
     \left(\frac{c_d\,m_d}{f_a}\,\gamma_\mu P_L\,\frac{\slashed{k}_a-\slashed{p}_{\bar{d}}+m_d}{m_a^2-2 k_a\cdot p_{\bar{d}}}\,\gamma_5\,-\,
     \frac{c_u\,m_u}{f_a}\,\gamma_5\,\frac{\slashed{k}_a+\slashed{p}'_{\bar{u}}-m_u}{m_a^2-2 k_a \cdot p'_{\bar{u}}}\,\gamma_\mu P_L\right)
                     \label{eq:bGamma_q}
\end{eqnarray}
Adopting the same phase conventions as \cite{Buras:1998raa}, one defines the neutral Kaon mass eigenstates:
\begin{eqnarray}
K^0_L &=& \frac{1}{\sqrt{2(1+|\tilde{\epsilon}|^2)}}\left( \left(1+\tilde{\epsilon}\right) K^0 + \left(1-\tilde{\epsilon}\right) \bar{K}^0\right) \\
K^0_S &=& \frac{1}{\sqrt{2(1+|\tilde{\epsilon}|^2)}}\left( \left(1+\tilde{\epsilon}\right) K^0 - \left(1-\tilde{\epsilon}\right) \bar{K}^0\right).
\end{eqnarray}
By making the following assignments for the initial and final quark momenta,
\begin{eqnarray}
p_{\bar{s}}  &=& x P_K \qquad , \qquad p'_u  = (1-x) P_K \nonumber \\
p_d &=& y P_\pi \qquad \, , \,\qquad p'_{\bar{u}} = (1-y) P_\pi \nonumber
\end{eqnarray}
the amplitude for the $K^0_L \to \pi^0\,a$ decay, when the ALP emitted by the $K^0_L$ meson, reads:
\begin{eqnarray}
\mathcal{M}_{K^0_L} &=&
  -\frac{ \tilde{\epsilon}\,G_F}{2 \sqrt{2}}\,\Re[V^*_{us} V_{ud}]\,f_{\bar{K}}\, f_{\pi} \, (k_a\cdot P_\pi) \frac{M_K}{f_a} \,\times\,
    \hspace{0.5cm} \nonumber \\
 & & \hspace{-0.75cm} \times \int^1_0 \left\{ \frac{c_s\,m_s\,\theta(x-\delta^K_a)}{m_a^2-2\,x\,k_a\cdot P_K} -
    \frac{c_d\,m_d\,\theta(1-x-\delta^K_a)}{m_a^2-2\,(1-x)\,k_a\cdot P_K} \right\} \,\phi_K(x)\,g_K(x)\, dx
\label{eq:MK0LALP}
\end{eqnarray}
once the trivial integration in $y$ is performed. As one can notice from Eq.~(\ref{eq:MK0LALP}), the amplitude for the $K^0_L$ decay is proportional
to the oscillation CP violation parameter $\tilde{\epsilon}$, as expected from general considerations on the CP properties of $K^0$ and $\bar{K}^0$
decays, and from the absence of imaginary part in the CKM for the tree--level diagram, i.e. $\Im[V^*_{us} V_{ud}]=0$. Consequently the
$K^0_L$ decay amplitude is suppressed by $O(10^{-3})$, with respect to the corresponding tree--level charged process. Conversely, the $K^0_S
\to \pi^0\, a$ decay amplitude can be obtained from the result of Eq.~\ref{eq:MK0LALP}, simply removing the $\tilde{\epsilon}$ suppression factor,
and would be of the same order as the charged process.

The amplitude contribution when the ALP is emitted from the $\pi^0$ meson is not explicitly reported here, showing the same $M_\pi/M_K$
suppression as for the charged decay case.

\subsection{The one-loop process}

Charged and neutral pseudo-scalar $K \to \pi \,a$ meson decays, assuming flavour conserving fermion-ALP interactions of Eq.~(\ref{eq:ALP_SMMFV}),
receive contributions at one--loop level\cite{Izaguirre:2016dfi,Bauer:2017ris,Gavela:2019wzg} from the diagrams shown in
Fig.~\ref{fig:figure_loop}. In the following only the contribution arising from fermion-ALP interaction will be considered.
\begin{figure}[!th]\center
\includegraphics[scale=0.15]{PICTURES/Loopcff.jpg}\hspace{1 cm}\includegraphics[scale=0.15]{PICTURES/Loopcw.jpg}
\caption{One-loop penguin contributions \label{fig:figure_loop}}
\end{figure}

In this kind of processes, only one quark line participate to the ALP emission, the other quark being a spectator. Customarily
the hadronization of a matrix element between two pseudo-scalar meson mediated by a vector current, where one of the quark does
not interact can be factorised as
\begin{eqnarray}
\bra{P}\bar{q}_1\gamma^\mu Q_2\ket{M}=f_+(q^2)(P_M+P_P)^\mu+f_0(q^2)\,q^\mu
\label{loop_formfactors}
\end{eqnarray}
with $q=P_M-P_P$. The form factors $f_{+,0}(0)=1$ in the isospin symmetric limit, while the non approximated, $q^2$ dependent, form factors
are obtained from LQCD calculation \cite{Carrasco:2016kpy}. From Eq.~(\ref{loop_formfactors}) the amplitude for the $K^+ \to \pi^+ \,a$
decay reads:
\begin{equation}
{\mathcal M}^{L}_{K^+} = \frac{G_F\, m^2_t}{4\sqrt{2}\pi^2} (V_{ts}V_{td}^*) \frac{M^2_{K^+}}{f_a}\left(1-\frac{M^2_{\pi^+}}{M^2_{K^+}}\right)
                           \! \left[f_+(m_a^2) + \frac{m_a^2}{M^2_{K^+}-M^2_{\pi^+}} f_-(m_a^2)\right]\!\sum_{q=u,c,t}\! \!c^{(q)}_{sd}
\label{loopAmplitudeK+}
\end{equation}
with the coefficient
\begin{equation}
c^{(q)}_{sd} = \frac{V_{qi}V_{qj}^*}{V_{ts}V_{td}^*} \left[3\,c_W \frac{g(x_q)}{x_t}-\frac{c_{q}\,x_q}{4\,x_t}
                   \ln\left(\frac{f_a^2}{m_q^2}\right)\right]
\label{ALPflavorViol}
\end{equation}
opportunely normalized in order to factorize out all the relevant scale dependences. The penguin with the ALP emitted from the
internal W line is included for completeness, even if in the following phenomenological analysis $c_W=0$ will be assumed.
The dominant contribution from the penguin diagram is mostly proportional to the $c_t$ coupling. For the $K$ meson decay, with the
charm contribution roughly accounting for $10\%$ of the total contribution.

One-loop diagrams, with the ALP emitted from the initial/final quarks can be safely neglected being suppressed by at least a factor
$m_s^2/m_W^2 \approx 10^{-6}$ with respect to the penguin contributions, as they arise at third order in the external momenta expansion.
Therefore, no sensitivity on the ALP--down quark couplings can emerge in the $K\to \pi \,a$ decays from one loop diagrams.

An order of magnitude of the tree vs loop amplitude ratio is obtained from comparing Eqs.~(\ref{eq:MKLALPapprox}) and (\ref{loopAmplitudeK+}),
giving:
\begin{eqnarray}
R_{T/L} = \left| \frac{\mathcal{M}^T_{K^+}}{\mathcal{M}^L_{K^+}} \right| \approx 2\,\pi^2 \frac{f_K \, f_\pi}{m_t^2}
\left|\frac{V^*_{us} \, V_{ud}}{V^*_{ts} \, V_{td}} \right| \simeq 1. \times 10^{-2} \,.
\end{eqnarray}
showing the expected level of suppression. Even if the tree vs loop ratio is at the per cent level, the tree level diagrams may have a
non negligible impact in the measurement of the $K \to \pi \,a$ decays, as in principle they depend on different and less constrained,
down quark--ALP couplings.

Finally, the loop contribution to the $K^0_L \to \pi^0 \,a$ decay can be easily obtained from Eq.~(\ref{loopAmplitudeK+}) and reads:
\begin{equation}
{\mathcal M}^{Loop}_{K^0_L} = \frac{G_F\, m^2_t}{4\,\sqrt{2}\pi^2} \Im(V_{ts}V_{td}^*) \frac{M^2_{K^+}}{f_a}\left(1-\frac{M^2_{\pi^+}}{M^2_{K^+}}\right)
                            \left[ f_+(m_a^2) + \frac{m_a^2}{M^2_{K^+}-M^2_{\pi^+}} f_-(m_a^2)\right]\, c^{(t)}_{sd} \nonumber
\label{loopAmplitudeK0L}
\end{equation}
being proportional to the non vanishing imaginary part of the CKM matrix.

\section{Bounds on ALP-fermion couplings}
\label{sec:pheno}

Armed with the tree--level and one--loop, charged and neutral, $K\to\pi\,a$ decays amplitudes obtained in the previous section,
one can bound the ALP-fermion couplings using the experimental limits provided by the NA62~\cite{CortinaGil:2020fcx,CortinaGil:2021nts},
E949~\cite{Artamonov:2009sz,Artamonov:2008qb,Adler:2008zza} and KOTO~\cite{Ahn:2018mvc} experiments. The main assumption underlying
the following phenomenological analysis is that the ALP lifetime is sufficiently long for escaping the detector (i.e. $\tau_a\gtrsim
100$ ps) or alternatively the ALP is mainly decaying in a, not better specified, invisible sector. Visible ALP decays have been studied,
for example, in \cite{DiLuzio:2020oah,Kelly:2020dda,Bauer:2017ris,Gavela:2019wzg}.

The tree--level amplitudes of Eqs.~(\ref{eq:MKALP}), (\ref{eq:MPiALP}) and (\ref{eq:MK0LALP}) depend on the ALP couplings with $s,d$ and
$u$ quarks, while the one--loop ones reported in Eqs.~(\ref{loopAmplitudeK+}), (\ref{ALPflavorViol}) and (\ref{loopAmplitudeK0L}), are
typically dominated by the ALP coupling with the heaviest quark running in the loop, the $t$ quark, being the $c,u$ contributions suppressed
by the $m_{u,c}/m_t$ mass ratio barring Cabibbo enhancements. Being the focus of this paper on ALP-fermion couplings, for the rest of the
section $c_W=0$ will be assumed. The interplay between the simultaneous presence of $c_W$ and $c_t$ has been discussed in detail in
\cite{Gavela:2019wzg}.

An analysis of the $K\to \pi\,a$ decay with completely general, but flavor conserving, ALP-quark couplings, would require to consider a
five-parameters fit, $(c_u,c_d,c_c,c_s,c_t)$ beside the ALP mass $m_a$. In order to obtain meaningful information about the ALP-fermion
couplings different simplifying assumptions have to be introduced. The phenomenological approach followed in this section will be twofold.
First of all, in Sec.~4.1, all ALP-fermion couplings, introduced in the Lagrangian of Eq.~(\ref{eq:ALP_SMMFV}), will be assumed independent.
Then, using the tree--level amplitudes for the s-- and t--channels, limits on $(c_u,c_s)$ and $(c_d,c_s)$ will be obtained respectively
from the charged and neutral $K$ meson decays, setting all the other ALP-quark couplings to $0$.
Afterwards, in Sec.~4.2, only two independent family universal ALP-fermion couplings, $c_\uparrow$ for the up quarks and $c_\downarrow$
for the down ones, will considered, for sake of simplicity. Under this assumption, the interplay between the tree--level and loop
contributions to the $K\to \pi \,a$ decay will be thoroughly discussed. Limits for the universal ALP--fermion coupling $c_{a\Phi}$ can
be then obtained straightforwardly.

\subsection{Tree--level Contributions}
\label{pheno-tree}

The tree--level amplitudes for charged and neutral $K$ decays, given by Eq.~(\ref{eq:MKALP}) and (\ref{eq:MK0LALP}) for charged and neutral
channels respectively, in the most general case depend on four parameters: the ALP couplings to the three light quarks, $c_u, c_d$ and
$c_s$ and the ALP mass, $m_a$. As derived in Eq.~(\ref{MpiKapprox}) the diagram with the ALP emitted by the pion contribution is strongly
suppressed, ranging from $10^{-2}$ (in the most conservative case) to $0$ if $g_\pi=0$ is assumed. Therefore, it seems reasonable, in the
following, to neglect the $\pi$-ALP emission diagrams, and, consequently the $K^+$ decay rate depends only on the $(c_u,c_s)$ ALP-fermion
couplings, while the $K^0$ decay rates only on $(c_d,c_s)$ ones.

\begin{figure}[t!]
\centering
\includegraphics[scale=0.28]{PICTURES/NewCharged_Final.pdf}\hspace{0.5 cm}\includegraphics[scale=0.28]{PICTURES/NewNeutral_Final.pdf}
\caption{Excluded parameter regions derived from tree-level channels for charged (left) and neutral (right) $K \to \pi a$ decays. In
the left--plot limits are derived from NA62 (pink) and E949 (cyan) experiments. In the right--plot, bounds are obtained from present (cyan)
and expected (pink) KOTO data.}
\label{fig:tree_plot}
\end{figure}

The left plot in Fig.~\ref{fig:tree_plot} shows the allowed regions of parameters as function of the ALP mass $m_a$ for the chosen
reference value $f_a=1$ TeV. The shaded gray area is excluded by present experimental data. The upper and lower contours, delimiting
the colored shaded area represent the bounds obtained setting $c_s=c_u$ and $c_s=-c_u$ respectively. The exclusive limit on $c_s\,(c_u)$
with $c_u\, (c_s) =0$ lies inside the colored shaded area. As noticed in Sec.~\ref{sec:hadronization}, for $c_s=c_u$ the $K^+$ decay
rate gets suppressed by an accidental parametric cancellation, leading to a less stringent bound on the ALP--fermion couplings. The
shaded colored area represent consequently the typical uncertainty in the bound prediction from $K^+\to \pi^+\,a$ decay rates once
letting the couplings $(c_u,c_s)$ freely varying in the $|c_s/c_u|\leq 1$ range. The pink contours and shaded region are obtained from
NA62 data \cite{CortinaGil:2020fcx,CortinaGil:2021nts} while the cyan ones refer to bounds obtained from the E949 \cite{Artamonov:2009sz}
experiment. For $m_a$ values below $0.15$ GeV the two experiments provide similar results, with a slight edge in favor of NA62, bounding
$(c_u,c_s) \lesssim 0.05$. In the $m_a > 0.15$ region, latest NA62 measurements has instead improved the sensitivity of E949 by roughly
a factor 10, bounding $(c_u,c_s) \lesssim 0.01$. For $m_a\approx m_{\pi^0}$ both experiments loose sensitivity. No significative effects
are obtained in this plot from the $\pi$-ALP emission diagrams, once the $c_d$ parameter is assumed to lie in the perturbative range.

A similar analysis, for $K^0_L\to \pi^0\,a$ decay is presented in the right plot of Fig.~\ref{fig:tree_plot}, where the cyan and pink regions
are obtained using the present\cite{Ahn:2018mvc} and expected KOTO experiment data, respectively. For this plot the upper and lower contours,
delimiting the shaded area represent the bounds obtained setting on $c_s=c_d$ and $c_s=-c_d$ respectively. KOTO experiment results much less
sensitive to the $(c_d,c_s)$ ALP-fermion couplings, as CP violation in the tree--level processes can occurs only through the CP-mixing
$\tilde{\epsilon}$ parameter, thus suppressing this channel by roughly a factor $10^{-3}$. Present KOTO data do not provide any real constraint,
with the prospect that future data could reach sensitivity to the perturbativity region\footnote{Notice that KOTO experiment provides only
mass independent limits on the $K^0_L\to \pi^0\,a$ branching ratio and consequently the loss of sensitivity in the $m_a \approx m_{\pi^0}$
region does not show up in the plot.}.

The results showed in Fig.~\ref{fig:tree_plot}, even if not looking flashy, represent, nonetheless, the most stringent model--independent
bounds on light quark couplings to ALP, for an  ALP mass in the sub--GeV range, once flavour conserving, but not flavor universal
ALP-fermion interactions, are assumed.

\subsection{Interplay between tree--level and one--loop contributions}
\label{pheno-loop}

\begin{figure}[t!]
\centering
\hspace{-3 ex}\includegraphics[scale=0.4]{PICTURES/NewLoop_Final.pdf}
\caption{Excluded parameter regions for an universal ALP--up quark coupling $c_\uparrow$ derived from NA62 (pink), E949 (cyan) and
KOTO expected (dashed line) experiments.}
\label{fig:sum_plot}
\end{figure}

To constrain, simultaneously, tree--level and one--loop contributions to the $K\to \pi\,a$ decay one has to adopt simplified  frameworks.
Following \cite{Gavela:2019wzg}, one can consider the scenario of universal ALP-quark coupling, $c_{a\Phi}$. From the analysis of
Sec.~\ref{sec:hadronization} one easily realizes that in this scenario, the top-penguin loop contribution dominates the charged and
neutral $K$ decay, once $c_W=0$ is assumed. The full cyan and pink lines in Fig.~\ref{fig:sum_plot}, represent the limits on $c_{a\Phi}$
obtained from E949 and NA62 respectively as function of the ALP mass $m_a$. The dashed gray line represents, instead, the $c_{a\Phi}$
limits from the expected KOTO upgrade. These results are in agreement with the bounds presented in \cite{Gavela:2019wzg} and show that
$K$ meson decays typically constrain $c_{a\Phi}\lesssim 10^{-3}$ in the sub-GeV ALP mass range.

\begin{figure}[t!]
\centering
\includegraphics[scale=0.28]{PICTURES/NewCorrelation2_Final.pdf}\hspace{0.5 cm}
\includegraphics[scale=0.28]{PICTURES/NewCorrelation1_Final.pdf} \\
\includegraphics[scale=0.28]{PICTURES/NewCorrelation2m_Final.pdf}\hspace{0.5 cm}
\includegraphics[scale=0.28]{PICTURES/NewCorrelation1m_Final.pdf}
\caption{Excluded parameter regions for universal ALP--up and down quark couplings $c_\uparrow$ and $c_\downarrow$ derived from NA62 (cyan)
and KOTO (pink and dashed pink line) and $Y(ns)\to \gamma \,a$ (full vertical black line) experiments. The upper plots refer to
the $\rm{sign}(c_\downarrow) = \rm{sign}(c_\uparrow)$ case, while in the lower ones $\rm{sign}(c_\downarrow) = - \rm{sign}(c_\uparrow)$ has
been chosen.}
\label{fig:corr_plot}
\end{figure}

In general MFV ALP frameworks, however, it may not be unconceivable to assign different, but flavor universal, PQ charges to the up and
down quark sectors, see for example \cite{Arias-Aragon:2017eww,Merlo:2017sun,Alonso-Gonzalez:2018vpc}, that in the following will be
denoted as $c_\uparrow$ and $c_\downarrow$, respectively. In this scenario, one--loop amplitudes only depend from $c_\uparrow$ while
the tree--level amplitudes are practically proportional to a linear combination of $c_\uparrow$ and $c_\downarrow$, as evident for example
in the simplified amplitudes of Eqs.~(\ref{eq:MKLALPapprox}) and (\ref{eq:MKHALPapprox}). From the tree--level analysis, summarized in
Fig.~\ref{fig:tree_plot}, one learns that present data limit $c_\downarrow$ to be typically below $10^{-1}$. Indeed, to study the interplay
between tree--level and one--loop the reference value $c_\downarrow=\pm 0.05$ has been chosen, somehow in the ridge of the parameters
allowed from the previous analysis on tree--level contributions. The blue and brown shaded regions showed in Fig.~\ref{fig:sum_plot}
represent the variability of NA62 and E979 bounds on $c_\uparrow$ once $c_\downarrow$ is let varying in the $[-0.05,0.05]$ range.
The presence of the tree-level contribution can modify the bounds on $c_\uparrow$ extracted from penguin diagrams of roughly one order
of magnitude, in all the $m_a$ range. The expected KOTO limits on the $K^0_L\to \pi^0\,a$ decay is reported in Fig.~\ref{fig:sum_plot}
as a black dashed line, giving a practically constant bound $c_\uparrow \lesssim 1 \times 10^{-3} $ over all the $m_a$ range of interest,
yet not competitive with the charged sector one. Notice that, however, the neutral $K$ decay sector does not suffer from any relevant
interference from the tree--level processes, largely suppressed from the CP violating parameter $\tilde{\epsilon}$. A simultaneous
measurement of the charged and neutral $K\to \pi\,a$ decays, may thus in a not too far future may give independent indications on the
relative size of the ALP--light quark couplings.

Finally, in Fig.~\ref{fig:corr_plot}, a summary on the combined bounds on $(c_\uparrow, c_\downarrow)$ is presented for two reference values
of the ALP mass $m_a=0$ GeV and $m_a=0.2$ GeV. For the two upper plots $\rm{sign}(c_\downarrow) = \rm{sign}(c_\uparrow)$ has been taken. In
the lower plots, where $\rm{sign}(c_\downarrow) = - \rm{sign}(c_\uparrow)$ has been considered, a partial cancellation between one--loop and
tree-level contributions takes place. In this second scenario, the $c_\downarrow$ constraint from the $\Upsilon(ns)$ decays at Babar
and Belle (full vertical black line) derived by \cite{Merlo:2019anv} can contribute to close this flat direction.