\section{Introduction}
Percolation of water through soil is one of the oldest problems in hydrodynamics. The fluid passes through a network of irregularly arranged interstices between solid objects. Typically, each individual thread of water passes through a narrow channel in which the equations of viscous flow can be applied. The difficulty arises from the fact that the detailed knowledge of the channels is neither available nor useful due to their complexity. We therefore typically take coarse-grained approaches, averaging over a length-scale much larger than the individual channels but still small compared to the scale of the medium. We thus define a relation between the flux, $q$, and the pressure-difference, $\Delta P$. In the simplest case of  a rigid isotropic medium, this gives rise to Darcy's law~\cite{Bat67}
\begin{equation}
q = -\frac{k}{\mu}\frac{\Delta P}{L} \/,
\label{eq:darcy}
\end{equation}
where $\mu$ is the dynamic viscosity of the fluid, $k$ is the permeability of the porous medium, and $L$ its length in the flow direction. Henceforth we shall call $q$ the Darcy flux. The permeability has the same status as all transport coefficients in hydrodynamics -- for a real system it is very difficult to calculate from first principles, but can be calculated in a model system by first solving a problem at the pore-scale and then by either analytical or numerical coarse-graining. This problem develops an additional degree of complexity if we consider that under the fluid stress the solid obstacles can move, i.e.,\  the flow itself can form channels. We treat the complication wherein the solid skeleton is deformable, i.e.,\  poroelasticity, the simplest example of which is the kitchen sponge.

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{fig1.pdf}
\caption{From left to right, a snapshot of the cross-section of one of our simulation under different levels of magnification. We magnify twice into a part of the domain to show first a sub-domain and then a channel. We apply the lubrication theory to study the flow thorough this channel, which implies that the flow inside the channel is assumed to be parabolic.}
\label{fig:schematic}
\end{figure*}

Poroelasticity play an important role in understanding the transport through a wide range of materials ranging from individual cells~\cite{charras2005non,moeendarbary2013cytoplasm}, to biological tissues, e.g.,\  soft-tissues~\cite{lai1991triphasic,yang1991possible,Auton2017}, bones~\cite{Cowin1999}, even to hydraulic fracture~\cite{detournay1988poroelastic,yarushina2013rock}. The simplest poroelastic problem is that of linear poroelasticity where we assume that the flow of the liquid is governed by the Darcy's law and the solid skeleton not only has linear constitutive relation but also undergoes small deformation. In reality, often the deformation of the solid matrix is large consequently nonlinear elastic effects have to be taken into account even if the constitutive relation is linear. Such systems are notoriously difficult to study both experimentally and numerically~\cite{MacMinn2015}.

The central question in this paper is how a coarse-grained description of the Darcy type emerges from a pore-scale model. As our model we choose a bed, a two-dimensional hexagonal lattice with defects, of soft elastic cylinders immersed in a liquid. Using both direct numerical simulations -- a set of fully coupled equations for a viscoelastic solid in contact with a Newtonian fluid -- and theory, we show that at scales that are large compared to the diameter of a cylinder the flux versus pressure-difference relationship in the system is a Darcy-like equation. When the deformability is small, as measured by the shear modulus of the solid, we obtain the Darcy equation exactly: the permeability $k$ is a constant, independent of the pressure-difference. However, as the solid skeleton becomes more deformable, the permeability becomes a nonlinear function of the pressure difference, namely for the same pressure drop we get a larger flux. Our theoretical calculations suggest that this result is largely model independent. This behavior has been already predicted from theoretical modelling at a coarse-grained level~\cite{macminn2016large} but has never been observed before in simulations or experiments.

\section{Numerical method}
We first describe our Direct Numerical Simulations (DNS). The deformable cylinders in the Newtonian fluid are modeled with a two-phase approach, defined by a variable $\phi = 0$ inside the viscoelastic solid phase and $\phi = 1$ in the fluid phase, with an evolving interface. The cylinders are organized on a hexagonal lattice, see \fig{fig:schematic}. If all the lattice sites are filled we reach the maximum solid volume fraction~\footnote{Actually we reach a packing fraction slightly lower than the maximum possible one by reducing the diameter of the cylinders such that they are initially placed in a manner that the surfaces of neighboring cylinders do not touch -- i.e.\ there is a small gap between their surfaces.}. In the rest of this paper we use a porosity $\Phi$, which is the fraction of the total volume occupied by the fluid, equal to $0.42$ by removing a certain number of randomly selected cylinders. The cylinders are made of a hyper-elastic Mooney-Rivlin \cite{bonet_wood_1997a} material characterized by a shear elastic modulus, $G$. We emphasize that the full non-linear structural response of the elastic solid is included in the simulations. The theoretical model considered later in this paper, however, is simpler. The motion of the fluid and of the viscoelastic material are governed by the conservation of momentum and the incompressibility constraint :
\begin{align}
\label{eq:NS2p}
\frac{\partial u_i}{\partial t} + \frac{\partial u_i u_j}{\partial x_j} = \frac{1}{\rho} \frac{\partial \sigma_{ij}}{\partial x_j} \;\;\; \textrm{and} \;\;\; \frac{\partial u_i}{\partial x_i} = 0,
\end{align}
where the Cauchy stress tensor $\sigma_{ij}$ has contributions from both solid and fluid stresses with a weight set by the phase variable $\phi$, i.e.,\
\begin{equation}
\label{eq:phi-stress}
\sigma_{ij} = \phi \sigma_{ij}^{\rm f} + \left( 1 - \phi \right) \sigma_{ij}^{\rm s}\/,
\end{equation}
with suffixes $^{\rm f,s}$ used to distinguish the two phases, fluid $\rm ^f$ and elastic solid $\rm ^s$. The fluid is assumed to be Newtonian and the solid is an incompressible viscous hyper-elastic material with constitutive equations:
\begin{subequations}
\begin{align}
\sigma_{ij}^{\rm f} &= -p \delta_{ij} + 2 \mu \mathcal{D}_{ij} \/, \label{eq:sf}\\
\sigma_{ij}^{\rm s} &= -p \delta_{ij} + 2 \mu \mathcal{D}_{ij} + \tau^{\rm e}_{ij}\/. \label{eq:ss}
\end{align}
\label{eq:sigma}
\end{subequations}
Here $p$ is the pressure, $\rho$ and $\mu$ are respectively the density and the dynamic viscosity both of which are assumed to be the same in the two phases\footnote{Note that, deformation of the solid becomes stationary and the value of its viscosity does not matter anymore since the solid viscosity only contributes when the solid is undergoing deformation prior to reaching the steady state.}, $\mathcal{D}_{ij}$ the rate-of-strain tensor and $\delta_{ij}$ is the Kronecker delta. The last term in $\sigma_{ij}^{\rm s}$, Eq.~(\ref{eq:ss}),  is the hyper-elastic contribution $\tau^{\rm e}_{ij}$, here modeled as a neo-Hookean Mooney-Rivlin material with the constitutive relation $\tau_{ij}^{\rm e}=G \mathcal{B}_{ij}$, where $\mathcal{B}_{ij}$ is the left Cauchy-Green tensor sometimes also called the Finger tensor. The full set of equations can be closed in a purely Eulerian manner by updating $\mathcal{B}_{ij}$ and $\phi$ with the following transport equations \citep{astarita_marrucci_1974a, larson_1988a, bonet_wood_1997a, sugiyama_ii_takeuchi_takagi_matsumoto_2011a}
\begin{align}
\label{eq:advB}
\frac{\partial \mathcal{B}_{ij}}{\partial t} + \frac{\partial u_k \mathcal{B}_{ij}}{\partial x_k} &= \mathcal{B}_{kj}\frac{\partial u_i}{\partial x_k} + \mathcal{B}_{ik}\frac{\partial u_j}{\partial x_k} \\
\frac{\partial \phi}{\partial t} + \frac{\partial u_k \phi}{\partial x_k} &= 0\/.
\label{eq:advP}
\end{align}
The algorithms have been described, used and validated against standard test cases in several earlier publications~\cite{rosti_brandt_2017a, rosti_brandt_mitra_2018a, rosti_brandt_2018a}, and more details can be found in these references.

We use a rectangular domain of size $MD \times ND$, where $D$ is the diameter of an undeformed cylinder, for three different sets of values for $M$ and $N$. We apply periodic boundary conditions in the stream-wise $x$-direction, and no-slip/no-penetration boundary conditions on the two rigid walls bounding the domain in the $y$-direction. We consider $6$ different values for the shear elastic modulus $G$, and for each one of them we impose $6$ different pressure differences from $0.5$ to $50$ to drive the flow and measure the resulting flux, resulting in a Reynolds number varying in the range $Re=\rho q D/\mu \in \left[ 10^{-5} : 10^{-4} \right]$ (in particular, we used $\rho=1$, $D=0.22$ and $\mu=1$ in our simulations). In all the cases, the numerical domain is discretised with $68$ grid points per diameter $D$. After a short transient time, the flow and the deformation of the solid skeleton reach a stationary state (see \fig{fig:schematic}). As $G$ decreases -- the solid skeleton is more deformable -- the flow changes the width and nature of the channels through which the liquid flows. We observe a general tendency of the flow to exploit the defects of the underlying lattice, by generating preferential channels which transport most of the fluid.

\section{Numerical Results}
\begin{figure}
\centering
\includegraphics[width=0.45\textwidth]{fig2.pdf}
\caption{The Darcy flux $q$ as a function of the overall pressure difference $\Delta P/L$ for different values of the deformability of the medium, $G \in \left[  0.5, 50 \right]$. In particular, the red, brown, orange, green, blue and black colors are used for $G=0.5$, $1$, $1.5$, $3$, $6$ and $50$, respectively. We also plot the results for different domain sizes (grey squares and cross). In the inset we plot the permeability $k$ normalized with $\mathcal{K}_t=C D^2 \Phi^3$ versus the dimensionless variable $\beta = \Delta P/G$. $C$ is a constant coming from the theory, equal to $0.1808$ for the present cases.}
\label{fig:darcy}
\end{figure}

In \fig{fig:darcy} we show how the Darcy flux, $q$, depends on the pressure-difference, $\Delta P$, across the domain: the most rigid case (black line) shows a linear increase of the flow rate with the pressure difference -- the standard Darcy's law for rigid porous materials. By contrast, as the material becomes more elastic, we observe a non-linear growth of the Darcy flux steeper than for its rigid counterpart, i.e.,\  a super-linear dependence of the Darcy flux on the pressure-difference. A different way to interpret the same result is to say that the permeability of the porous medium, $k$, is itself a nonlinear function of the pressure difference (or the flow rate).

We check how robust this result is in the following ways. We run simulations in three domain sizes, the smallest one being, $M=8$ and $N= 9$, in the next one we double the size in each direction, $M=16$ and $N=18$ and then obtain the largest one by again doubling the size, $M=32$ and $N=36$. We obtain the same result in all the three cases. Next, in the largest domain, we select sub-domains of the same size as the smallest domain and the relationship shown in \fig{fig:darcy} remains the same for these sub-domains too. Note that, we can reach the same porosity in different ways, depending on the position of the defects in the material; we have checked that our results remain unchanged in two different realisations of the random defects.

\section{Theoretical Model}
Next we construct a theory for this behavior. Our specific numerical simulations act as a motivation for the theory but the theory is not necessarily limited to our numerical model. Typically, the theory of porous media involves multiple scales,~\cite{Whitaker1986a,  Whitaker1986c, Mei1989, Mei1991, Auriault1992,macminn2016large, Collis2017}, ours is no exception.

In \fig{fig:schematic} we show a two-dimensional cross-section of our domain on three different spatial scales, from left to right we go from the full domain (size $L$) down to the scale of sub-domains ($\mathcal{L}$), smaller than the full domain but still much larger than the size of a single deformable circle, down to the scale of a single channel ($\ell$) between the elastic cylinders, such that $\ell \ll \mathcal{L} \ll L$. Our first step is to derive a relation between the flux and pressure difference across a deformable two-dimensional channel, e.g.,\  we solve the microscale problem sketched in the rightmost panel of \fig{fig:schematic}.

Let $\ell$ be the length of this channel, ${\Delta p}$ the pressure-difference across the channel and $h(\xi)$ the width of the channel as a function of the stream-wise coordinate, $\xi$. Within the range of parameters used in our simulations, we can safely assume that within this channel the Reynolds number is so low that the flow can be described by the Stokes equations. In fact, we shall go one step further and assume that it is safe to use the lubrication approximation~\cite{Bat67}. In particular, we assume that within a channel the pressure is a function of the stream-wise direction alone, i.e.,\  $p = p(\xi)$, the wall-normal component of the velocity is zero, and the velocity gradient along the stream-wise direction is much smaller than that in the wall-normal direction. In addition, the flow velocity must go to zero at the boundaries of the channel, hence the flow-rate through the channel is given by
\begin{equation}
  q = - [h^3/(12 \mu)] (dp/d\xi)\/.
\end{equation}
The flow-rate must be a constant, independent of $\xi$. The width of the channel, $h(\xi)$, is determined by the mutual interaction between the flow and the elastic property of the walls of the channel. To make further progress, we assume a Hookean response of the boundary of the channel, often called Winkler foundation~\cite{wang2005beams,dillard2018review} in other context. In this framework, the undeformed width of the channel is $h_{\rm 0}$, which together with the deformation $w(\xi)$ sets the total width of the channel, $h(\xi) = h_{\rm 0} + w(\xi)$; the elastic property of the channel walls is parameterized by a Hookean spring with a spring constant $\kappa$, such that the force-per-unit-area necessary to generate a deformation $w$ is given by $\kappa w$. Hence the pressure and the deformation are related by $w(\xi) = p(\xi)/\kappa $. This allows us to write a differential equation for $p(\xi)$ where the flow-rate, $q$, appears as a parameter -- the equation has the general form of $q = -\sigma(p)(dp/d\xi)$~\cite{rubinow1972flow}. We integrate it and enforce the result to conform to the form of Darcy's law, $q= (\mathcal{K}/\mu)({\Delta p}/\ell)$ with
\begin{eqnarray}
\quad
\label{eq:mK}
& &  \mathcal{K} = h_{\rm 0}^3 f(\beta), \quad \mbox{where} \\
\label{eq:qofp}
& &  f(\beta)=\frac{1}{12}\left[1+ \frac{3}{2}\beta + \beta^2 + \frac{1}{4}\beta^3 \right] \/.
\end{eqnarray}
Here, $\beta \equiv {\Delta p}/(\kappah_{\rm 0})$ and we have used the boundary conditions $p(0) = {\Delta p}$ and $p(\ell)=0$. To build a connection to our simulations it is appropriate to choose $\kappa$ such that $\kappah_{\rm 0}= G$. The use of the lubrication approximation coupled with the elastic properties of solid is quite commonly used to analyze flows in deformable channels, see e.g.,\  \citet{davis1986elastohydrodynamic,gro+jen04,gomez2017passive, christov2018flow}. \citet{christov2018flow} contains derivation of a similar relationship using a systematic application of asymptotics for a three dimensional channel. The only difference is that in \citet{christov2018flow} a different model for the elastic wall -- isotropic quasi-static bending of a plate under a transverse load due to the fluid pressure -- is used.

In the next step we consider the mesoscopic scale, larger than the size of single cylinders but still smaller than the scale of the whole bed, see the middle panel in \fig{fig:schematic}. The large domain contains many such mesoscopic domains of the same size. In the event of no defects, each of these  subdomain contains $m  \times n $ cylinders organized on a regular hexagonal lattice and the channels between the cylinders form a regular honeycomb lattice. In this case, each subdomain has exactly the same porosity and the same permeability.

Recall that we have randomly removed few cylinders from a regular hexagonal lattice to create our porous medium. Thus, the porosity at the scale $\mathcal{L}$, $\Phi_{\mathcal L}$, is different in different sub-domains. In the inset of \fig{fig:coarsePor} we show a representative plot of $\Phi_{\mathcal L}$, the porosity averaged over a domain of size $\mathcal{L}$, extracted from the DNS with the largest domain. In particular, we perform a volume average of the local fluid fraction $\phi$ on a domain of size $\mathcal{L}$ resulting in the porosity $\Phi_\mathcal{L}$. We incorporate this randomness into our model by choosing different values of $h_{\rm 0}$ -- undeformed width of the channel -- in different sub-domains. Thus at this mesoscale our model is a honeycomb network of channels. The length of the channels is same in all subdomains and within each subdomain all the channels have the same width set by the subdomain porosity $\Phi_{\mathcal L}$.  Consequently, the width of the channels is different in different subdomains.

For a single subdomain, our task is to calculate the effective permeability of a network of channels, $\mathcal{K}_{\mathcal L}$, where the flow rate in each channel is given as function of the pressure drop by $q = [\mathcal{K}({\Delta p})/\mu]({\Delta p}/\ell)$, where the permeability of each channel $\mathcal{K}({\Delta p})$ is the nonlinear function, $f(\beta)$ in \eq{eq:qofp}. Given a nonlinear function $f(\beta)$ there is no general method of attack known to us. We proceed therefore by assuming that $\mathcal{K}$ is independent of ${\Delta p}$. In this case, the problem corresponds to that of the effective conductivity of a honeycomb network of resistors by mapping the $q$ to current, ${\Delta p}/\ell$ to the voltage drops across the bonds of the network, and $\mathcal{K}/\mu$ to the conductivity of each of the bond, in bus-bar geometry~\cite{redner2009fractal} -- the network is connected to two parallel lines and the battery is connected across the two lines. We solve this linear problem by matrix inversion to obtain $\mathcal{K}_{\mathcal L} = \gamma\mathcal{K} (\mathcal{L}/W)$ where $W$ is the width of the sub-domain and $\gamma$ is a constant that depends on $m$ and $n$. In particular, we use $m = 8, \; n = 9$ and obtain $\gamma = 3.047$.

Thus, we obtain $\mathcal{K}_{\mathcal L} = C D^2 \Phi_{\mathcal L}^3 f(\beta)$, with a constant $C = \gamma (3/8)^3 (\mathcal{L}/W)$. Here we have used $\mathcal{L} = 8D$, and $\mathcal{L}/W = 1.125$. The form of the function in \eq{eq:qofp} suggests that the effective permeability, $\mathcal{K}_{\mathcal L}$, is solely a function of the dimensionless parameter, $\beta = \Delta P/G$, where $\Delta P$ is the pressure-difference across the sub-domain. In \fig{fig:coarsePor} we show a scatter plot of $\mathcal{K}_{\mathcal L}$ versus $\Phi_{\mathcal L}$ for different sub-domains and observe that $\mathcal{K}_{\mathcal L}$ grows as $\Phi_{\mathcal L}^3$, although with some scatter of the data. Notwithstanding this, when we plot the values of $\mathcal{K}_{\mathcal L}$ normalized by $CD^2\Phi_{\mathcal L}^3$ as a function of $\beta$ for the different sub-domains, we obtain a reasonable data-collapse in agreement with our theory.

There is another, equivalent, way to calculate the effective permeability of a mesoscale subdomain, which also requires the assumption of linearity. In particular, we incorporate the random removal of cylinders by mapping to the problem of a honeycomb network of resistors such that  each bond in the network has a conductance of $\mathcal{K}$ with probability $\mathcal{P}$ or infinite conductance (zero resistance) with probability $1-\mathcal{P}$. This probability, $\mathcal{P}$,  is different in each sub-domain. Using the expression for effective resistance of infinite but random lattices~\cite{redner2009fractal}, we obtain the effective permeability to be $\mathcal{K}_{\mathcal L} = \frac{\alpha}{\bra{\alpha + 1} \mathcal{P} - 1} \mathcal{K} $, where $\alpha$ is a geometric parameters that depends on the lattice.

In the last and final step, we average over different sub-domains to obtain an effective permeability for the whole domain. Taking the divergence of the Darcy's flux in a sub-domain, we obtain
\begin{equation}
\divg{ \sqbra{ \mathcal{K}_{\mathcal L}(\Phi_{\mathcal L}) \grad{p} } } = 0,
\end{equation}
which is a steady-state heat equation with variable diffusivity $\mathcal{K}_{\mathcal L}(\Phi_{\mathcal L})$, function of a fast variable $\Phi_{\mathcal L}$. Straightforward application of the method of multiple scales~\cite[see, e.g., ][section 9.6.2]{Fri96} shows that the effective permeability $k$ is the harmonic mean of the effective permeability of each sub-domain:
\begin{equation}
 \frac{1}{ k} = \int d\Phi_{\mathcal L} \frac{P(\Phi_{\mathcal L})}{ \mathcal{K}_{\mathcal L}(\Phi_{\mathcal L})}\/,
\label{eq:Kbar}
\end{equation}
where $P(\Phi_{\mathcal L})$ is the probability density function of the porosity of a sub-domain, $\Phi_{\mathcal L}$. The data in \fig{fig:coarsePor} justify treating $P(\Phi_{\mathcal L})$ as a Gaussian with mean value equal to the mean porosity calculated over the full domain so that \eq{eq:Kbar} is integrable. As the Gaussian is sharply peaked, the integral is well-approximated by its leading order contribution, i.e.,\  $k \approx \mathcal{K}_{\mathcal L}(\Phi)$ where $\Phi = \abra{\Phi_{\mathcal L}}$ is the mean porosity of the whole domain. This implies that the collapse we have observed for each sub-domain should also work if we plot the permeability $k$ of the full domain as a function of $\beta=\Delta P/G$ where the pressure-difference, $\Delta P$, now is across the whole domain. This is confirmed by the results in the inset of \fig{fig:darcy}.

\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{fig3.pdf}
\caption{A scatter plot of the local permeability $\mathcal{K}_{\mathcal L}$ in each sub-domain normalized with $CD^2$  as a function of the corresponding local porosity $\Phi_{\mathcal L}$. The inset shows a snapshot of porosity (averaged over the scale $\mathcal{L}$), $\Phi_{\mathcal L}$. The color scale goes from $0.3$ (black) to $0.6$ (white).}
\label{fig:coarsePor}
\end{figure}

\section{Discussion and Conclusions}
Several comments are now in order.  For the analytical calculations we have used a simple model, the Winkler foundation, whereas we have used the hyper-elastic Mooney-Rivlin model for the cylinders in our DNS. The qualitative agreement between the two shows the robustness of our results. Different elastic models will result in different expression for the function $f(\beta)$ in \eq{eq:qofp}. A crucial result of our work is to show that such a function exists, which implies that data on permeability collapse  to a single function when plotted as a function of $\beta$.

The enhanced flux, which is the most striking result of our work, has not been observed in experiments but is found analytically for certain classes of models. In particular, two among the five models discussed in \citet{macminn2016large} -- these models start from the intermediate scale denoted as $\mathcal{L}$ here-- show the possibility of super-linear response because of the nonlinear elastic behavior of the solid skeleton. These two models further assume that $\mathcal{K}_{\mathcal L}$ is independent of $\Phi_{\mathcal L}$.

In experiments~\cite{Hewitt2016,song_stone_jensen_lee_2019a}, when a fluid is forced through a deformable porous medium, the boundary between the porous material and the fluid on the inlet is normally left unconstrained. Hence under fluid pressure the boundary moves and squeezes the porous material. This decreases the permeability -- often modeled by the empirical Kozeny-Carman formula~\cite{Carman1997} -- of the porous material. Two of the models in \citet{macminn2016large}, which include the Kozeny-Carman formula, show nonlinear but sub-linear behavior.  Recent experimental measurements by \citet{song_stone_jensen_lee_2019a} also exhibit such behavior. In our case, both theory and the DNS approximate the behavior of the Kozeny-Carman function for small $\Phi_{\mathcal L}$, $\mathcal{K}_{\mathcal L} \sim \Phi_{\mathcal L}^3$. Indeed, in our simulations the boundaries are held fixed, hence by construction fluid-driven compaction is missing from our simulations. Hence, we also expect that it is possible to observe the enhanced Darcy flux in experiments, but not for very large pressure-differences where fluid-driven compaction dominates. We hope our work will encourage further experimental and numerical explorations.

Most studies in this field using homogenization to understand the fluid flow through rigid/deformable/active porous media~\cite{Whitaker1986a, Whitaker1986c, Mei1989, Mei1991, Auriault1992, Collis2017}, adopt a continuum description for both the solid and fluid phase and couple them through the kinematic interface conditions. The constitutive equations for the pore-scale description of the problem are coarse-grained to obtain a description of the equivalent fluid-solid interaction. The transport coefficients of the resultant equations depend on the solvability conditions (closure problem) of the homogenization techniques. However, the closure problem remained unsolved in all these models, thus no explicit Darcy-like relation was obtained. Our theory stands apart from such models. We write down a Darcy-like relation between flow-rate and pressure-difference with an explicit expression for the permeability that depends on the shear modulus of the solid skeleton. Furthermore, we show that the nonlinear flow-rate versus pressure-difference relation is an intrinsic property of the medium rooted in the pore scale rearrangement induced by fluid flow. The weakest link in our theory is the assumption of linearity to calculate the effective permeability of a network of channels.

\section*{Acknowledgements}
Our work is inspired by experiments being done by John Wettlaufer and his group. We thank John for helping us at every stage of this work. We thank Dominic Vella for introducing us to the Winkler foundations. MR and LB are supported by the European Research Council Grant no. ERC-2013-CoG-616186, TRITOS and by the Swedish Research Council Grant no. VR 2014-5001. SP acknowledges the support of the Swedish Research Council Grant no. 638-2013-9243. DM acknowledges the support of the Swedish Research Council Grant no. 638-2013-9243 as well as 2016-05225. We gratefully acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing).


\section{Introduction}

Understanding nucleon structure in terms of the fundamental degrees of freedom of Quantum Chromodynamics (QCD) is one of the main goals in the theory of strong interactions.   In recent years it became clear that  exclusive reactions may provide information
about hadron structure
encoded in so-called Generalized Parton Distributions  \cite{Ji,Radyushkin} (GPDs).
For each quark flavor $q$
there are eight GPDs. Four correspond to  parton helicity-conserving (chiral-even) processes,  denoted
by $H^q$,  $\tilde H^q$,  $E^q$ and  $\tilde E^q$, and
four correspond to parton helicity-flip (chiral-odd) processes  \cite{ji,diehl},  $H^q_T$,  $\tilde H^q_T$,  $E^q_T$ and  $\tilde E^q_T$.
The GPDs depend on three kinematic variables: $x$, $\xi$ and $t$. In  a symmetric frame,  $x$ is the   average longitudinal momentum fraction of the struck parton before and after the hard interaction and $\xi$ (skewness) is half of the longitudinal momentum fraction transferred to  the struck parton. The skewness can be expressed in terms of the   Bjorken variable $x_B$  as
$\xi\simeq x_B/(2-x_B)$. Here $x_B=Q^2/(2p\cdot q)$ and $t=(p-p^\prime)^2$, where $p$ and $p^\prime$ are the initial and final four-momenta of the nucleon.

When the theoretical calculations for longitudinal virtual photons were compared with the JLab  $\pi^0$ and $\eta$ data\cite{clas1,clas2,clas3}
they were  found  to  underestimate the measured cross sections by more than an order of magnitude in their accessible kinematic regions.
The failure to describe the experimental results with quark helicity-conserving operators  stimulated a consideration of the role of the  chiral-odd quark helicity-flip processes. Deeply virtual meson  electroproduction (DVMP), and in particular $\pi^0$ production in the reaction $ep\to e^\prime p^\prime \pi^0$, was
identified~\cite{Ahmad:2008hp,G-K-09,Goloskokov:2011rd}
as especially sensitive to the quark helicity-flip subprocesses.
During the past few years, two parallel theoretical approaches - \cite{Ahmad:2008hp,GL}~(GL) and \cite{G-K-09,Goloskokov:2011rd}~(GK) have been developed utilizing the  chiral-odd GPDs in the calculation of  pseudoscalar meson electroproduction. The GL and GK approaches, though employing different models of  GPDs, lead to {\it transverse} photon amplitudes that are much larger than the longitudinal amplitudes.

\section{Definition of structure functions}

The unpolarized reduced meson cross section is described
 by 4 structure functions $\sigma_T$, $\sigma_L$, $\sigma_{TT}$ and $\sigma_{LT}$:
\begin{equation}
2\pi\frac{d^2\sigma(\gamma^*p\to p\pi^0)}{dtd\phi_\pi} =
\frac{d\sigma_T}{dt} +
\epsilon  \frac{d\sigma_L} {dt}+
\epsilon  \frac{d\sigma_{TT}} {dt}  \cos 2\phi+
\sqrt{2\epsilon (1+\epsilon)}   \frac{d\sigma_{LT}} {dt}  \cos \phi.
\end{equation}

\noindent
References~\cite{Goloskokov:2011rd,GL} obtain the following relations for unpolarized structure functions:

\begin{equation}
\label{SL}
\frac{d\sigma_{L} }{dt}= \frac{4\pi\alpha}{k^\prime}\frac{1}{Q^4}\left\{ \left( 1-\xi^2 \right) \left|\langle \tilde{H} \rangle\right|^2 -2\xi^2 {Re}\left[ \langle \tilde{H} \rangle^* \langle \tilde{E} \rangle \right] - \frac{t^\prime}{4m^2} \xi^2 \left| \langle \tilde{E} \rangle \right|^2 \right\}
\end{equation}

\begin{equation}
\label{ST}
\frac{d\sigma_{T}}{dt} = \frac{4\pi\alpha}{2k^\prime Q^4} \left[ \left(1-\xi^2\right) \left|\langle H_T \rangle\right|^2 - \frac{t'}{8m^2} \left|\langle \bar{E}_T \rangle\right|^2\right]
\end{equation}

\begin{equation}
\label{SLT}
\frac{d\sigma_{LT}}{dt} = \frac{4\pi\alpha}{\sqrt{2}k^\prime Q^4} \xi \sqrt{1-\xi^2} \frac{\sqrt{-t'}}{2m} { Re} \left[ \langle H_T\rangle^* \langle\tilde{E}\rangle \right]
\end{equation}

\begin{equation}
\label{STT}
\frac{d\sigma_{TT}}{dt} = \frac{4\pi\alpha}{k^\prime Q^4}\frac{t'}{16m^2}\left|\langle \bar{E}_T \rangle\right|^2
\end{equation}

\noindent
Here $m$ is the mass of the proton, $t^\prime =t-t_{min}$, where $|t_{min}|$ is the minimum value of $|t|$ corresponding to $\theta_\pi =0$, $k^\prime(Q^2,x_B)$ is a phase space factor and  $\bar E_T = 2\widetilde H_T + E_T$.
The brackets $\langle  H_T \rangle$ and $\langle \bar E_T \rangle$ denote
the convolution of the elementary process
$\gamma^*q\to q\pi^0$
with the GPDs $H_T$ and $\bar E_T$. They are called them generalized form factors.

\begin{figure*}[t!]
\centering
\includegraphics {fig16_scissored.pdf}
\caption
{
(Color online)
The extracted structure functions vs. $t$ for
the $\pi^0$ (left column)
and $\eta$ (right column). The top row presents data for the kinematic point  ($Q^2=$1.38~GeV$^2$,$x_B$=0.17)
and bottom row for the kinematic point ($Q^2=$2.21~GeV$^2$,$x_B$=0.28).
The data and curves are as follows:
black circles - $d\sigma_U/dt =d\sigma_T/dt +\epsilon d\sigma_L/dt$,
blue triangles - $d\sigma_{TT}/dt$,
red  squares - $d\sigma_{LT}/dt$.
The error bars are statistical only. The gray bands are our estimates of the absolute normalization systematic uncertainties on $d\sigma_U/dt $.
The curves are theoretical predictions produced  with the GPG model of Goloskokov and Kroll~\cite{Goloskokov:2011rd}.
}
\label{fig:structure}
\end{figure*}

\section{Experimental data}
The cross section of the reaction $ep\rightarrow ep(\pi^0/\eta)$ measured by
the CLAS collaboration at Jlab in bins of $Q^2$, $x_B$, $t$ and $\phi$ were published in
 Refs.
 \cite{clas1,clas2,clas3}.
Structure functions $\sigma_U=\sigma_T+\epsilon\sigma_L$, $\sigma_{LT}$ and $\sigma_{TT}$ have been extracted from the angular distributions. These functions were compared
with the predictions of the GPD models \cite{Goloskokov:2011rd,GL}. The result confirmed that the measured unseparated cross sections are much larger than expected from
leading-twist handbag calculations which are dominated by longitudinal photons.
As an example, the comparison of the $\pi^0$ and $\eta$ structure functions is shown in Fig.~\ref{fig:structure} for two kinematical bins in $x_B$ and $Q^2$.
The structure functions
$\sigma_U$ and $\sigma_{TT}$ for $\eta$ are, respectively,  factors of 2.5 and 10 smaller than for $\pi^0$.
However, the GK GPD model \cite{Goloskokov:2011rd}  (curves)  follows the experimental data.
Taken together, the $\pi^0$ and $\eta$ results stringthen
the statement about the transversity GPD dominance in the pseudoscalar electroproduction process.

\begin{figure}[t!]
\vspace*{-10 mm}
\centerline{
\includegraphics[height=15cm]{ht_et_pi0_eta.pdf}
}
\vspace*{-27mm}
\caption{(Color online) Data points: CLAS.
Top left: $\left|\langle H_T \rangle^\pi\right|$,
top right: $\left|\langle \bar{E}_T \rangle^\pi\right|$,
bottom left: $\left|\langle H_T \rangle^\eta\right|$,
bottom right: $\left|\langle \bar{E}_T \rangle^\eta\right|$ as a function of -$t$ for different kinematics:
($Q^2$=1.2 $GeV^2$,$x_B$=0.15) black,
($Q^2$=1.8 GeV$^2$,$x_B$=0.22) red,
($Q^2$=2.2 GeV$^2$,$x_B$=0.29) blue,
($Q^2$=2.7 GeV$^2$,$x_B$=0.34) magenta.
}
\label{fig:ht_et_pi0_eta}
\end{figure}

\section{Generalized form factors}
The squared magnitudes of the
generalized form factors
$\left|\langle H_T \rangle\right|^2$ and
$\left|\langle \bar{E}_T \rangle\right|^2$
may be directly extracted from the experimental data (see Eqs. \ref{ST} and \ref{STT} ) in the framework of GPD models.
\begin{equation}\label{ET}
\left|\langle \bar{E}_T \rangle^{\pi,\eta}\right |^2=\frac{k^\prime Q^4}{4\pi\alpha}     \frac{16m^2}{t'} \frac{d\sigma^{\pi,\eta}_{TT}}{dt}
\end{equation}
\begin{equation}\label{HT}
\left| \langle H_T \rangle^{\pi,\eta}\right|^2=\frac{2k^\prime Q^4}{4\pi\alpha}       \frac{1}{1-\xi^2} \left [ \frac{d\sigma^{\pi,\eta}_{T}}{dt} + \frac{d\sigma^{\pi,\eta}_{TT}}{dt} \right ].
\end{equation}
\noindent
Figure~\ref{fig:ht_et_pi0_eta} presents the modulus of the generalized form factors
$\left|\langle H_T \rangle^\pi\right|$, $\left|\langle \bar{E}_T \rangle^\pi\right|$, $\left|\langle H_T \rangle^\eta\right|$ and $\left|\langle \bar{E}_T \rangle^\eta\right|$
 for 4 different kinematics.
  Note the dominance of the $\left|\langle \bar{E}_T \rangle\right|$ over $\left|\langle H_T \rangle\right|$ for both $\pi^0$ and $\eta$.
Generalized form factors $\langle H_T \rangle^\pi$ and $\langle \bar{E}_T \rangle^\pi$ are shown in more detail in  Fig.~\ref{fig:ht_et_pi0}.
The  $\langle \bar{E}_T \rangle^\pi$ formfactor has steeper t-dependemce than $\langle H_T \rangle^\pi$.
The t-slope parameters, obtained by an exponential fit of the form $e^{bt}$, are
 $b(\langle \bar{E}_T \rangle)$=1.27~$GeV^{-2}$ and $b(\langle H_T \rangle)$=0.98~$GeV^{-2}$.

\section{Flavor decomposition}

In  electroproduction of $\pi^0$ and $\eta$ mesons the GPDs $F_i$ appears in the following combinations
\begin{equation}
F^\pi_i=\frac{1}{\sqrt{2}}[e_uF^u_i-e_dF^d_i]
\end{equation}
\begin{equation}
F^\eta_i=\frac{1}{\sqrt{6}}[e_uF^u_i+e_dF^d_i-2e_sF^s_i]\\
\end{equation}
\noindent
The $q$ and $\bar q$ GPDs contribute in the quark combinations $F^q_i-F^{\bar q}_i$.
Hence there is no contribution from the strange quarks if we assume that
$F^s_i\simeq F^{\bar s}_i$.
For flavor decomposition we have to take into account the decay constants $f_\pi$ and $f_\eta$, the chiral  condensate
constants $\mu_{\pi^0}=$2.57~GeV, $\mu_1=$0.958~GeV and $\mu_8$=2.32~GeV, and the contribution from singlet and octet $\eta$ states \cite{Goloskokov:2011rd}.

\begin{equation}
F^\eta_i=F^{\pi}_i\left(\cos\theta_8-\sqrt{2} \frac{\mu_1}{\mu_8}\frac{f_1}{f_8}\sin\theta_1\right)\frac{f_8}{f_{\pi^0}}\frac{\mu_8}{\mu_{\pi^0}}=\frac{ F^8_i}{k_\eta},
\end{equation}
\noindent
where the mixing angles are:
$\theta_8=-21.2^o$ and $\theta_1=-9.2^o$.
The octet and singlet wave functions are very similar and the decay constants are close as well
$f_8=1.26f_\pi$ and
$f_1=1.17f_\pi$.
The overall factor for the $\eta$ meson is
$k_\eta=0.863$.
Using $e_u=\frac{2}{3}$ and $e_d=-\frac{1}{3}$ we will end up with
\begin{equation}\label{decom}
F^\pi_i=\frac{1}{3\sqrt{2}}[2F^u_i+F^d_i]
\end{equation}
\begin{equation}
{k_\eta}{F^\eta_i}=\frac{1}{3\sqrt{6}}[2F^u_i-F^d_i].
\end{equation}
\noindent
Experimentally we have  access only to the $\left|\left<F_i^\pi\right>\right|^2$ and $\left|\left<F_i^\eta\right>\right|^2$
(see Eq.~\ref{ET}-\ref{HT}).
The final equation for the $\left<H_T\right>$ convolution reads

\begin{equation}
\frac{1}{18}       \left | 2\langle H_T \rangle^u + \langle H_T \rangle^d  \right  |^2 =  \left |\langle H_T \rangle^\pi \right |^2\\
\end{equation}
\begin{equation}
 \frac{1}{54}       \left  | 2\langle H_T \rangle^u - \langle H_T \rangle^d  \right  |^2 = k_\eta^2\left |\langle H_T \rangle^\eta \right |^2.
\end{equation}
\noindent
and simular equations for $\langle \bar{E}_T \rangle$.

\begin{figure}[t!]
\vspace*{-10 mm}
\centerline{
\includegraphics[width=10cm]{ht_et_pi.pdf}
}
\vspace*{-25 mm}
\caption
{
(Color online) Generalized form factors $|\langle H_T \rangle^\pi|$   and $|\langle \bar{E}_T \rangle^\pi|$
as a function of -$t$ for $Q^2=2.2~GeV^2$  and $x_B$=0.27.
Top: $|\langle \bar{E}_T \rangle^\pi|$ in red;
Bottom: $|\langle H_T \rangle^\pi|$ in blue.
}
\label{fig:ht_et_pi0}
\end{figure}

\begin{figure}[t!]
\vspace*{-10 mm}
\centerline{
\includegraphics[width=15cm]{ud_3_s.pdf}
}
\vspace*{-35 mm}
\caption
{
(Color online) Flavor separated generalized form factors $\langle H_T \rangle$ and $\langle \bar{E}_T \rangle$
as a function of -$t$ for $Q^2=2.2~GeV^2$  and $x_B$=0.27.
Top: $\langle H_T \rangle^u$ (red) and $\langle H_T \rangle^d$ (blue);
Bottom: $\langle \bar{E}_T \rangle^u$ (red) and $\langle \bar{E}_T \rangle^d$ (blue).
}
\label{fig:ud_3}
\end{figure}

\begin{figure*}[t!]
\vspace*{-25 mm}
\centerline{
\includegraphics[width=15cm]{den_01_02.pdf}
}
\vspace*{-50 mm}
\caption
{
(Color online) Impact parameter density  of quarks which are transversely
polarized along $b_x$-axis
in an unpolarized
proton.
Top left and right panels are for u-quarks with $x$=0.1 and $x$=0.2, respectively.
Bottton left and right panels are for d-quarks with $x$=0.1 and $x$=0.2.
}
\label{fig:den_01-02}
\end{figure*}

The solution of these equations will lead to the flavor decomposition of the generalized form factors
$\langle H_T \rangle^u$ and  $\langle H_T \rangle^d$ as well as  $\langle \bar{E}_T \rangle^u$ and $\langle \bar{E}_T \rangle^d$.
However, the  convolution integrals have real and imaginary parts. So it is impossible to  solve these equations unambiguously with only two equations in hands.
So, in order to estimate the form factors, we make an ad hoc assumption that the
relative phace  $\Delta\phi$ between $\langle H_T \rangle^u$ and $\langle H_T \rangle^d$  equals 0 or 180 degrees.
Ignoring an overall phase, the form factors are then real, and we arbitrarily choose the solution with
  $\langle H_T \rangle^u$ and  $\langle \bar{E}_T \rangle^u$ positive.
Fig.~\ref{fig:ud_3}  presents
$\langle H_T \rangle^u$, $\langle H_T \rangle^d$, $\langle \bar{E}_T \rangle^u$ and $\langle \bar{E}_T \rangle^d$
for one kinematic point $(Q^2=2.2~GeV^2,x_B=0.27)$ calculated with this assumption.
Note the different signs of $\langle H_T \rangle^u$ and  $\langle H_T \rangle^d$  and the same sign
of $\langle \bar{E}_T \rangle^u$ and $\langle \bar{E}_T \rangle^d$.

\section{Quark spin densities in the transverse plane and generalized transversity distributions}
Two-dimensional Fourier transforms of GPDs $H(x,\xi=0,-{\vec \Delta^2)}$ and $\bar E_T(x,0,-{\vec \Delta^2)}$, where $\vec \Delta^2=-t$,
define the spin density of the polarized quarks in an unpolarized proton \cite{Diehl:2005jf}.
\begin{equation}
H(x,{\vec b})=\int \frac{d^2{\vec \Delta}}{(2\pi)^2} e^{-i{\vec b \vec \Delta}} H(x,0,-{\vec \Delta}^2)
\end{equation}
\begin{equation}
\bar E_T(x,{\vec b})=\int \frac{d^2{\vec \Delta}}{(2\pi)^2} e^{-i{\vec b \vec \Delta}} \bar E_T(x,0,-{\vec \Delta}^2)
\end{equation}
\noindent
The GPDs $\bar E_T(x,\xi=0,t)$ and $H(x,\xi=0,t)$ were parametrized in the form \cite{Diehl:2013}
$$
\bar F^q(x,\xi=0,t)=g^q(x) \cdot \exp[(f^q(x)t],
$$
where $g^q(x)$   and $f^q(x)$ are the GPD forward limit and profile functions respectively.
The Fourier transform for this parametrization reads

\begin{equation}
\bar F^q(x,b)= \frac{1}{4\pi} \frac{g^q(x)}{f^q(x)}\exp[\frac{-b^2}{4f^q(x)}].
\end{equation}

\noindent For quarks polarized along $b_x$-axis the impact parameter density reads \cite{Diehl:2005jf}

\begin{equation}
\delta^q(x,{\vec b})= \frac{1}{2} [ H^q(x,{\vec b})-\frac{b_y}{m} \frac{\partial}{\partial  b^2}\bar E^q_T(x,{\vec b}) ].
\end{equation}
\noindent
The GPD $H(x,{\vec b})$ describes the density of unpolarized quarks and $\bar E_T(x,{\vec b})$ is related to the distortion of the polarized quark distribution in the transverse plane.
We can map the $u$ and $d$-quark transverse spin density distributions as a function of Feynman $x$ based on the GPD model  \cite{Goloskokov:2011rd} tuned to describe the CLAS data.
For example, Fig.~\ref{fig:den_01-02} shows  the impact parameter density  of transversely
polarized quarks along the $b_x$-axis  in an unpolarized  proton for  Feynman $x$=0.1 and $x$=0.2.
Note the distortion along $b_y$-axis,  similar for  u and d-quarks.
Looking at the transverse quark density distribution,
we can say that this width is diminished as $x\to 1$.
This behavior is typical  for the GPD models.

\section{Conclusion}
Differential cross sections of exclusive $\pi^0$ and $\eta$ electroproduction have been obtained in the few-GeV region at more than 1800 kinematic points in  bins of $Q^2, x_B$, $t$ and $\phi_\pi$.
Virtual photon structure functions
$\sigma_U$, $\sigma_{TT}$ and $d\sigma_{LT}$ have been obtained. It is found that $\sigma_U$ and $\sigma_{TT}$ are comparable in magnitude with each other, while $\sigma_{LT}$ is very much smaller than either.
Generalized form factors of the transversity GPDs $\langle H_T \rangle^{\pi,\eta}$ and $\langle \bar{E}_T \rangle^{\pi,\eta}$ were directly extracted from the experimental observables
for the first time. It was found that the GPD $\bar E_T$ dominates in  pseudoscalar meson production.
The combined $\pi^0$ and $\eta$ data opens the way for the flavor decomposition of the transversity GPDs.
Within some simplifying assumptions, the  decomposition has been  demonstrated.
 The spin density of  polarized  $u$ and $d$-quarks  in the transverse plane was evaluated for different values   of $x$
 from the  GPD model tuned to described the experimental data.

\section*{Acknowledgments}
The author  thanks G. Goldstein, S. Goloskokov, P. Kroll, S. Liuti and A.~Radyushkin  for many informative discussions and making available the results of their calculations.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177.


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

The role of data today is no longer limited to the verification of first-principles-derived-models but it is also being used to learn such models. This is particularly important in non-autonomous nonlinear dynamical systems that describe a multitude of problems from science and engineering.
Recent methods leverage the fact that most dynamical equations governing physical systems contain a few terms, making them sparse in high-dimensional nonlinear function space \citep{brunton2016discovering,rudy2017data}. By constructing an appropriate feature library based on the data coordinates, one can apply sparse regression to discover the governing equations of the dynamical system. Initial work in this field has mainly focused on the behavior of the approach with regards to noise, multi-fidelity data etc., however few have tried to improve upon the sparse regression algorithm at the core of the approach.

This is exactly the first focus area of the current work. We first look at the sparse regression method most commonly employed in this field: LASSO \cite{tibshirani1996regression}. We point out the theoretical underpinnings of LASSO along with accuracy bounds.
Although LASSO works well for uncorrelated features, for highly correlated it tends to choose a feature at random from each of the correlated groups \cite{gauraha2017dual}.
When the feature library is built from data, higher order features are more correlated with each other, where LASSO potentially runs into difficulties.
To alleviate these difficulties, we propose to solve the dual of LASSO to learn the governing equations.
Even in the case of correlated features, the dual LASSO has a unique solution, which then allows us to correctly choose the present features.
The second part of this work deals with the case when the exact function blocks that describe the dynamical system are not present in the feature library. We propose and investigate an approach to handle such cases by using an appropriate family of orthogonal functional basis to span the feature library.
As the result may not be sparse in terms of the orthogonal functions, we resort to sequentially thresholded least squares (STRidge) algorithm \cite{rudy2017data} that does not impose sparsity, but chooses the dominant features.

These approaches are demonstrated on the Lorenz 63 system \cite{lorenz1963deterministic}and the quadratic Lorenz system \cite{eminauga2015modified}.
We compare the performance of the LASSO and dual LASSO, and then use the approach of data-driven orthogonal regression to learn the feature library.

\subsection{General Methodology}\label{sec:general_method}

Let us assume that we have $n$ state space parameters ($x_1, \ldots, x_n$), with measurements for $x_i$ and $\dot{x_i} = dx/dt$ at times $t = 1, \ldots, T$ (denoted by a superscript).
If only state observations are available, the rate parameters can be computed using finite difference.
This is followed by constructing a non-linear library of features using the state space parameters. The span of these features now describes the feature space.
Typically we would construct this feature space through a class of functions that are dense in the space that our dynamical system lives in.
In this work we assume a polynomial feature library.
After constructing the feature library (say $X$), we formulate the regression problem as $\dot{X} = XW + \varepsilon$, where $\dot{X}_{(t, j)} = \dot{x_j}^t$ and $W$ are the unknown weights, with $\varepsilon$ being the noise.

Often, dynamical models have a functional form where the terms live in a smaller manifold in the feature space.
We thus use sparse regression to pick the relevant features.
These features, with their corresponding coefficients describe the functional form of the governing equations. This procedure is pictorially illustrated in figure \ref{fig:schematic}. For noisy observations, we carry out data filtering and smoothing at the input step and perform multiple sparse regression solves until convergence.
However we do not focus on this, and we assume that we have noise free observations.

\begin{figure}[h]
	\centering
	 \includegraphics[width=\linewidth]{schematic_v3.png}
	\caption{Schematic of the generic methodology used.}
	\label{fig:schematic}
\end{figure}
\vspace*{-0.25in}

\section{Sparse regression over fixed feature space} \label{sec:fixed}

In this case, we assume that the feature library is fixed, and that we wish to find either the exact sparse equation form from this library or the closest approximation to the governing equation only from the terms in the library.
The highest polynomial degree in the feature space ($X$) is $p$. Then, the feature space contains terms of the form $(x^t_1)^{p_1} \ldots (x^t_i)^{p_i} \ldots (x^t_n)^{p_n}$, such that $p_1 + \ldots + p_n \leq p$. The number of solutions to this system (\ie the number of terms in the feature library) is $m = {n+p \choose n} = \frac{(n+p)!}{n!p!}$.
Whereas, empirically the number of distinct terms in the governing equations is $\mathcal{O}(n)$. Thus even for small enough $p$, the terms in the feature library are much more in number than those to be chosen, which justifies sparse regression to select the features. Let us denote the coefficient matrix, obtained from the sparse regression solve by $W$. The optimization problem with some form of penalty ($\mathcal{P}$) is:
\begin{equation} \lb{reg}
\min_{W} \mathcal{L}(W) = \left[ \( \dot{X} - XW \)^2 + \mathcal{P}(W) \right] \ \ \text{where} \ \ \dot{X} \in \RR^{T \times n} \ , \ X \in \RR^{T \times m} \text{ and, } W \in \RR^{m \times n}
\end{equation}
To further select the features appropriately, we use our knowledge of the underlying physics of the dynamical system. Often, we have knowledge about the general scales of the state space parameters (for example, in the Lorenz system, the state space parameters represent the rate of convection, horizontal temperature variation and the vertical temperature variation). Thus, we know the the typical scales of these parameters from physics. Let us denote the scale of $x_i$ by $L_i$, then $(x^t_1)^{p_1} \ldots (x^t_n)^{p_n} \approx L_1^{p_1} \ldots L_n^{p_n}$. Thus, instead of selecting features based on their coefficient values, we select features by looking at their net magnitude, by replacing the state space parameters with their respective scales.
We refer to this as `scale based thresholding'.

\subsection{Sparsity through $L_1$ regularization (LASSO)} \label{sec:lasso}

We now look at the use of LASSO for imposing sparsity. As is well-known, the LASSO penalty is $\mathcal{P}(W) = \lambda ||W||_1$, which serves as a convex counterpart to the non-convex $L_0$ norm. \citet{kakade2009complexity} showed that the Rademacher averages for linear regressors with $L_1$ penalty are bounded by $\mathcal{R}_n(\mathcal{F}_W) \leq X W_{\max} \sqrt{\frac{2 \log(m)}{n}}$, where $W_{\max} = \max ||W||_1$. We can then use the contraction property, assuming the Lipschitz constant of the LASSO problem to be $L$ (more details and methods to compute it may be found in \cite{konevcny2016mini,konecny2013semi}) to obtain the bound of Rademacher averages of the loss function class $\mathcal{G}$. Using the symmetrization lemma, this yields a bound on the expected maximum error (uniform deviation), as given by \eref{symm}.
\begin{equation} \lb{symm}
\E \max_{g \in \mathcal{G}}\left[ \E g(x_i) - \frac{1}{n} \sum_{i = 1}^{n} g(x_i) \right] \leq 2 L XW_{\max} \sqrt{\frac{2 \log(m)}{n}}
\end{equation}
Note that even though we cannot explicitly compute the right hand side, we have $m >>> n$, which renders this accuracy bound impractical. Thus, LASSO in the case of exceedingly high number of features does not provide good accuracy bounds. To alleviate this problem, we first use the SAFE bounds provided by \citet{ghaoui2010safe} to remove the irrelevant features and then apply LASSO over the remaining library. Even after weeding out the absent features, another important task that remains is to choose the hyperparameter $\lambda$. Literature widely suggests the order of $\lambda$ to be $\lambda \approx \sqrt{T\log(m)}$ \cite{bickel2009simultaneous}, however choosing this parameter independent of the degree of correlation is not recommended \cite{dalalyan2017prediction}. It is further shown that the larger the correlation amongst the features, the smaller the value of $\lambda$ \cite{hebiri2013correlations}. \citet{hebiri2013correlations} also prove that LASSO can achieve fast rate bounds (\ie error bounded by $\lambda^2$) even for correlated features. Even though reassuring, in practice this requires significant amount of tuning and cross validation, which is a major issue in the case at hand, as typically we only have sparse and limited observations of the physical states. Further, even though the LASSO fit is unique in case of high correlations, they are problematic for parameter estimation and feature selection \cite{osborne2000lasso,tibshirani2013lasso}. The pitfalls of LASSO (even after removing the irrelevant features using the SAFE rules) are that it requires significant hyperparameter tuning and it is extremely sensitive to $\lambda$ for correlated features (observed empirically). These motivate us to instead formulate a new approach to solve the sparse regression problem.

\subsection{Sparsity through dual LASSO regression} \label{sec:dual}

To overcome the difficulties in the application of LASSO (along with the SAFE rules), we formulate and solve its dual problem. In order for the LASSO solution (that is, the feature and parameter choices) to be unique, the feature matrix must satisfy the irrepresentability condition (IC) of \citet{zhao2006model} and beta-min condition of \citet{buhlmann2011statistics}. The feature library violates the IC for highly correlated columns, leading to an unstable feature selection. However, even for highly correlated features, the corresponding dual LASSO solution is always unique \cite{tibshirani2013lasso}. The dual problem is given by \eref{dual} \cite{osborne2000lasso}, which is strictly convex in $\theta$ (implying a unique solution).
\begin{equation}\lb{dual}
\max_{\theta} \mathcal{D}(\theta) = ||\dot{X}||_2^2 - ||\theta - \dot{X}||_2^2 \quad \quad \text{such that} \quad \quad ||X^T \theta||_{\infty} \leq \lambda
\end{equation}
Now, let $\hat{W}$ be a solution of \eref{reg} with LASSO penalty and $\hat{\theta}$ be the unique solution to the corresponding dual problem \eref{dual}. Then stationarity condition implies $\hat{\theta} = \dot{X} - X\hat{W} \lb{sta}$.
Even though LASSO does not have a unique $\hat{W}$, the fitted value $X\hat{W}$ is unique, as the optimization problem \eref{reg} is strongly convex in $XW$ for $\mathcal{P}(W) = \lambda ||W||_1$. We make use of this by first computing a solution to the primal LASSO problem and then computing the unique dual solution by using the primal fitted value and \eref{sta}. Once we have the unique dual solution $\hat{\theta}$, we do feature selection by using the dual active set, which is same as the primal active set with high probability under the IC \cite{gauraha2017dual}. It is also be shown that features discarded by SAFE rules are the inactive constraints over the optimal dual solution \cite{wang2013lasso}, which justifies omitting them. KKT conditions imply:
\begin{equation} \lb{kkt}
\hat{\theta}^T X_i
\begin{cases}
= \text{sign}(\hat{W}_i) \text{ if } \hat{W}_i \neq 0 \\
\in (-1, 1)  \text{ if } \hat{W}_i = 0
\end{cases}
\end{equation}
\Eref{kkt} gives us a direct way to compute the active dual set (which is the same as the active primal set with high probability) once we have $\hat{\theta}$. We discard the features for which $\hat{\theta}^T X_i \in (-1, 1)$ and retain the others. This does not give us a good fit of the solution, so to compute the coefficients accurately, we perform ridge regression ($\mathcal{P}(W) = \lambda_2 ||W||_2^2$ over these active features. This pseudocode for this is given in algorithm \ref{alg:dual}. We refer to this algorithm as `dual LASSO'.

\begin{algorithm}
	\caption{Sparse regression using dual LASSO} \label{alg:dual}
	\begin{algorithmic}
		\Require state parameters: $\bm{x} = x_i^t, \dot{\bm{x}} = \dot{x}_i^t$; LASSO penalty $\lambda$, ridge penalty $\lambda_2$
		\State Compute the primal LASSO solution: $\hat{W} = \min_{W} \left[ \( \dot{X} - XW \)^2 + \lambda ||W||_1 \right]$
		\State Compute the unique dual solution $\hat{\theta} = \dot{X} - X\hat{W} \lb{sta}$
		\State Compute dual active set (same as primal active with h.p.) $S^d = \{1 \leq j \leq m : \hat{\theta}^T X_j \notin (-1, 1) \}$
		\State Construct reduced feature matrix $\tilde{X}$ by only considering features whose indices are in $S^d$
		\State Solve ridge regression $W^* = \min_{W} \left[\( \dot{X} - XW \)^2 + \lambda_2 ||W||_2^2 \right]$
	\end{algorithmic}
\end{algorithm}
\vspace*{-0.1in}

\section{Sparse Regression Over Adaptive Feature Space} \label{sec:adap}

In this section, we look at cases where the feature library is not known. If we have no prior belief over the form of the equations, we may not be able to construct an efficient feature library. In such situations, learning this library from data might be the most advantageous choice.
One approach is to make the use of orthogonal functions of some parametric family, such as polynomials to construct this library. This ensures no correlation between the features and affords maximum freedom to the algorithm. The drawback in this case is that the regressor may not be sparse over this feature library.

\subsection{Adaptive feature space growth through orthogonal regression} \label{sec:ortho_reg}

We propose the following approach for this purpose: Starting with an empty library, we recursively add a feature to it and compute the corresponding loss function of the resulting fit by using STRidge (\sref{stridge}). If the loss function decreases by more than a certain fraction, we keep this feature. Otherwise we discard it and look at the next orthogonal feature. Once every few of the addition timesteps, we perform a removal step.
That is, we discard the feature(s) that do not result in a significant increase in the loss function (if any).
This ensures that we do not keep lower order functions that may not be required to describe the equations as higher order functions are added.
Our algorithm is inspired by previous greedy feature development algorithms by \citet{efron2004least} (LARS), \citet{zhang2009adaptive} (FoBa) etc.
However these require pre-determined full possible feature space, whereas we construct new features on the fly. Once the equations are obtained in terms of these orthogonal polynomials, we distill their sparse forms by using symbolic equation simplification \cite{bailey2014automated}.

\subsection{Sequentially thresholded ridge regression} \label{sec:stridge}

To compute regressors over the orthogonal feature space, we use sequentially thresholded ridge regression (STRidge), proposed by \citet{rudy2017data}. The idea is simple: we iteratively compute the ridge regression solution with decreasing penalty proportional to the condition number of $X$, and discard the components using scale based thresholding (\sref{fixed}).
As the feature matrix is orthonormal by construction, the analytical solution is $W = \(1 + \lambda \)^{-1}X^T\dot{X}$.
This algorithm is effective in choosing optimal features without necessitating sparsity.
The overall pseudocode for learning the governing equations through adaptively growing the feature library is given by algorithm \ref{alg:adap}, and the corresponding results are presented in \sref{res}.

\begin{algorithm}
	\caption{Learning the governing equations through adaptive growth of the feature library} \label{alg:adap}
	\begin{algorithmic}
		\Require state parameters: $\bm{x} = x_i^t, \dot{\bm{x}} = \dot{x}_i^t$; orthogonal family $F_j(\bullet)$; feature addition / removal thresholds: $r_{a} \ (\leq 1) , r_{r} \ (\geq 1), \lambda_0$; removal step frequency $k_r$
		\State Initialize: $X = \emptyset, W = \bm{0}, t = 0, \mathcal{L} = \infty$
		\While{True}
		\State $X_{t} = \text{append}(X, F_k(\bm{x}))$
		\State Solve the STRidge problem: $W_{t} = \text{STRidge}(\dot{X}, X_{t}, \lambda_0)$
		\State Compute the loss $\mathcal{L}_{t} = \( \dot{X} - X_{t} W_{t} \)^2$
		\If{$\mathcal{L}_t \leq r_a \mathcal{L}$}
		\State $X = X_t \ ; \ W = W_t$
		\EndIf
		\If{$\mod(k, k_r) == 0$}
		\For{$i = 1, \ldots, X.\text{shape}[1]$} \quad \quad \quad \quad \quad \quad \quad \quad \quad \ \ \ (number of columns of $X$)
		\State $X_t = \text{append}\(X[:, 1:i-1], X[:, i+1:\text{end}]\)$ \quad (ignore the $i^{th}$ column of $X$)
		\State Solve the STRidge problem: $W_{t} = \text{STRidge}(\dot{X}, X_{t}, \lambda_0)$
		\State Compute the loss $\mathcal{L}_{t} = \( \dot{X} - X_{t} W_{t} \)^2$
		\If{$\mathcal{L}_t \leq r_r \mathcal{L}$}
		\State $X = X_t \ ; \ W = W_t$
		\EndIf
		\EndFor
		\EndIf
		\State $k = k+1$.
		\State \textbf{break} if no change in feature space over multiple iterations.
		\EndWhile
		\State Perform symbolic simplification of $\dot{X} = XW$ to obtain the final form of the equations
	\end{algorithmic}
\end{algorithm}

\section{Results} \label{sec:res}

For the first part (\sref{fixed}), our testbed will be the Lorenz 63 system ($n = 3$), given by \eref{lorenz}.
\begin{equation}\lb{lorenz}
\dot{x} = 10(yz - x) \ ; \quad \quad \quad \dot{y} = x(28 - z) \ ; \quad \quad \quad \dot{z} = xy - 2.667 z
\end{equation}
For the first part, we consider three polynomial feature libraries with $p = 3, 10$ and $20$ ($m = 20, 286$ and $1771$). The idea behind considering larger orders ($p$) is that it highlights the poor performance of LASSO for highly correlated features.
\Fref{residual} shows the residual after the LASSO and the dual LASSO model fits. We can see that the residuals for both methods are comparable for all of the feature libraries, which empirically shows that the model fit even for a high degree of correlations is good for both the methods.
\Fref{bar} plots the number of non-zero features in the equations for different $p$ values. LASSO has a much higher number of non-zero terms, and this number increases significantly with $p$ (and $m$), indicating instability of the solution.
Dual LASSO performs very well, and the number of present features does not change for the most part with $p$.
\Fref{coeffs} plots the absolute weights for the components for the $p = 3$ case for the $\dot{y}$ equation.
Dual LASSO retrieves the correct features (with accurate weights), while LASSO detects the correct features but also detects high order features that have low weights and are highly correlated to each other. This serves as a great validation of the superiority of dual LASSO over conventional LASSO for model discovery.

Now we look at the results for \sref{adap}, on the quadratic Lorenz system (\eref{lorenz_quad}) \cite{eminauga2015modified}.
\begin{equation}\lb{lorenz_quad}
\dot{x} = 10(yz - x) \ ; \quad \quad \quad \dot{y} = x(28 - z) \ ; \quad \quad \quad \dot{z} = (xy)^2 - 2.667 z
\end{equation}
We start with an empty feature library and $W = \bm{0}$ and iteratively grow the feature space using algorithm \ref{alg:adap} using Legendre polynomials (denoted by $\leg_p(\bullet)$), $r_a = 0.75, r_r = 1.25, \lambda_0 = 1$ and removal step working every 10 addition steps ($k_r = 10$). The final obtained model is given in \eref{final_adap}. We use symbolic simplification \cite{bailey2014automated} to simply this model. Note that the model is not unique before symbolic simplification (as $\leg_0(x) = \leg_0(y) = \leg_0(z)$). We also do scale based thresholding (\sref{fixed}) once after the symbolic simplification to remove any terms that may remain due to approximate factorization. The results before and after scale based thresholding are given in \eref{fin_1}. We can see that this algorithm does a good job at identifying the correct governing equations by iteratively building the feature library, and does not include any incorrect features.

\vspace*{-0.1in}
\begin{align} \lb{final_adap}
\begin{aligned}
	\dot{x} =& 9.93 \leg_1(y) \leg_1(z) - 9.89 \leg_1(x) \\
	\dot{y} =& 27.66 \leg_1(x) - 1.04\leg_1(x)\leg_1(z)
\end{aligned}
\ \
\begin{aligned}
	\dot{z} = &0.43 \leg_2(x) \leg_2(y) + 0.22 \leg_2(x) + 0.21 \leg_2(y) \\
	 -&2.62 \leg_1(z) + 2.09 \leg_0(x) - 0.22 \leg_0(y) - 1.95 \leg_0(z)
\end{aligned}
\end{align}
\begin{align} \lb{fin_1}
\underbrace{\begin{aligned}
\dot{x} &= 9.93 yz - 9.89x \\
\dot{y} &= 27.66x - 1.04xz \\
\dot{z} &= 0.97(xy)^2 + 0.007(x^2 - y^2) \\
&- 2.62z + 0.027
\end{aligned}}_{\text{after symbolic simplification}} \quad \implies \quad
\underbrace{\begin{aligned}
\dot{x} &= 9.93 yz - 9.89x \\
\dot{y} &= 27.66x - 1.04xz \\
\dot{z} &= 0.9675(xy)^2 - 2.62z
\end{aligned}}_{\text{after scale based thresholding}}
\end{align}

\section{Conclusions and Future Work} \label{sec:concl}

In this work, we investigated LASSO, proposed dual LASSO and data driven feature learning approaches to solve the problem of discovering governing equations only from state parameter data.
Future work directions involve extending the ideas of feature library building where one can construct the function to be added through a mix of a larger family of orthogonal functions. Approaches to involve kernel compositionality \cite{duvenaud2013structure} to build these libraries can also be investigated.

\begin{figure}[thp!]
	\centering
	\vspace{0.135in}
	\subfloat[Residuals after the model fit\label{fig:residual}]{\includegraphics[trim=0 0 0 65, width=0.5\textwidth]{residual.pdf}}
	\subfloat[Number of nonzero terms in the equations for $\dot{x}, \dot{y}, \dot{z}$\label{fig:bar}]{\includegraphics[trim=0 0 0 65,width=0.5\textwidth]{nonzero_bar.pdf}} \\
	\subfloat[The absolute weights in the equation for $\dot{y}$ \label{fig:coeffs}]{\includegraphics[trim=0 0 0 40,width=0.65\textwidth]{nonzero_coeffs_with_terms.png}}
	\caption{Normalized (combined) residual after model fit, corresponding number of non-zero coefficients in the equation for $\dot{y}$, and the absolute weights for the components in the $\dot{y}$ eqation ($p = 3$).}
	\label{fig:Y_post}
\end{figure}

\bibliographystyle{plainnat}


\section{Introduction}
With flexible mobility and deployment, unmanned aerial vehicles (UAVs) or drones have found a plethora of new applications in the recent years. Typical use cases of commercial UAVs include cargo delivery, surveillance and inspection, search and rescue, aerial photography, among the others. As projected by the Federal Aviation Administration (FAA), the UAV industry will generate more than 82 billion dollars for the U.S. economy alone and create more than 100,000 new jobs in the next decade. The prosperous global market of UAVs is also envisioned to  bring new and valuable opportunities to the future wireless communication industry, such as the forthcoming fifth-generation (5G) cellular network. On one hand, to enable reliable communications for widely deployed UAVs in the future, a promising solution is to integrate UAVs into the future 5G cellular network as new aerial users served by the terrestrial base stations (BSs). In fact, recent studies by the 3rd generation partnership project (3GPP) have demonstrated the viability to support the basic communication requirements for UAVs with the existing cellular network\cite{zeng2018cellular}. On the other hand, with continuously miniaturized BSs/relays, it becomes more feasible to mount them on UAVs and make complementary aerial communication platforms in 5G cellular network to provide or enhance the communication services for the terrestrial users\cite{zeng2016wireless}.

\begin{figure}[!t]
\centering
\includegraphics[width=6in]{UAV_Legi.eps}
\caption{Eavesdropping and jamming attacks to legitimate UAV communications in  5G  wireless network.}\label{UAV_Legi}
\vspace{-12pt}
\end{figure}
However, the integration of UAVs into 5G also brings new challenges that need to be addressed. Compared to terrestrial wireless channels that in general suffer from severe path-loss, shadowing, and multi-path fading, the high altitude of UAVs generally leads to more dominant line-of-sight (LoS) channels with the ground nodes. Although the strong LoS links can be exploited to improve the communication performance, they also cause severe interference to the terrestrial communications and thus effective air-ground interference management techniques are needed\cite{zeng2018cellular}. Moreover, the unique LoS-dominant UAV-ground channel poses new security challenges to the future wireless network with  various UAV applications as shown in Fig.\,\ref{UAV_Legi}. Specifically, in Fig.\,\ref{UAV_Legi}(a) and Fig.\,\ref{UAV_Legi}(b), the UAV is deployed as an aerial BS or relay that sends/forwards private messages to a legitimate ground node, respectively. However, the strong UAV-ground LoS links also enhance the reception quality of the terrestrial eavesdroppers. As a result, even a remote eavesdropper can take advantage of the LoS link to overhear the air-to-ground (A2G) communication clearly. Furthermore, in Fig.\,\ref{UAV_Legi}(c) and Fig.\,\ref{UAV_Legi}(d), the UAV acts as a mobile hub to collect private data from terrestrial Internet-of-things (IoT) devices such as sensors, or an aerial user that receives critical control signals from its serving ground BS. In such cases, a ground jammer can also exploit the strong ground-to-air (G2A) LoS links to launch more powerful jamming attacks to interfere with the UAV, thus significantly degrading its reception quality and even resulting in UAV communication failure. Therefore, compared to terrestrial communications, A2G and G2A communications are more susceptible to terrestrial eavesdropping and jamming, respectively. As such, how to effectively safeguard the legitimate UAV communications in the future wireless network is a new and challenging problem to resolve.

\begin{figure}[!t]
\centering
\includegraphics[width=6in]{UAV_Mali.eps}
\caption{Eavesdropping and jamming attacks to terrestrial wireless network by malicious UAV nodes.}\label{UAV_Mali}
\vspace{-12pt}
\end{figure}
On the other hand, UAVs may be a potentially new security threat to the terrestrial cellular network if they are misused by unauthorized parties for malicious purposes. As shown in Fig.\,\ref{UAV_Mali}(a), if the UAV is deployed as a malicious eavesdropper, it could more clearly overhear the terrestrial private communications between the ground BSs and their users even in a wide area as compared to the traditional terrestrial eavesdroppers, due to the more LoS-dominant G2A links. Moreover, the UAV may also be used as a malicious jammer to disrupt both the uplink and downlink terrestrial communications in the cellular network, as shown in Fig.\,\ref{UAV_Mali}(b). Its strong A2G links thus impose more severe as well as more  widely spread interference to legitimate ground communications than conventional terrestrial jammers with the same jamming signal power. In addition, the malicious UAVs can leverage their high mobility to flexibly change locations so that they can keep track of their moving targets over time and thus overhear/jam their communications more effectively. In light of the above, another challenging security problem arises, i.e., how to efficiently protect the terrestrial communications against the powerful  UAV eavesdropping/jamming attacks?

Considering the above two new UAV security issues in future wireless networks, this article aims to identify their key challenges and propose  promising solutions to tackle them, mainly from a physical-layer (PHY) design perspective. Note that PHY security has been extensively studied in the literature to combat the malicious eavesdropping and/or jamming attacks in terrestrial communication networks (see, e.g. \cite{yang2015safeguarding,jameel2018comprehensive} and the references therein). Typical anti-eavesdropping techniques include artificial noise (AN), cooperative jamming, transmit beamforming, etc. \cite{mukherjee2014principles}; while direct sequence spread spectrum (DSSS) and frequency hopping spread spectrum (FHSS) are two commonly used anti-jamming techniques at PHY\cite{mpitziopoulos2009survey}. However, such techniques consider mainly terrestrial nodes as legitimate or malicious users in the network and thus they may be ineffective to deal with the new and more challenging LoS-induced security issues with UAVs. As such, it is necessary to revisit the conventional techniques for UAV security and devise new solutions to safeguard future wireless networks more effectively.

\section{Securing Legitimate UAV Communications in Wireless Network}
As compared to terrestrial communications, legitimate A2G/G2A communications are more susceptible to ground eavesdropping/jamming due to the LoS-dominant channels. In this section, we propose promising countermeasures to improve UAV communication security and reliability.

\subsection{Protecting A2G Communication from Terrestrial Eavesdropping}
As shown in Fig.\;\ref{UAV_Legi}(a) and Fig.\;\ref{UAV_Legi}(b), for A2G communications where UAVs transmit confidential messages  to their ground nodes, the received signal strength is strong over a large area on the ground due to the LoS-dominant channels, which makes the prevention of terrestrial eavesdropping highly difficult. Nonetheless, by exploiting the high altitude of UAVs and their high mobility in the three-dimensional (3D) space, secure A2G communication may still be achievable, with the use of the following techniques.

\subsubsection{\bf UAV 3D Beamforming}
The miniaturization of light-weight antennas and steadily increasing payload of UAVs have made it feasible to equip UAVs with multiple antennas such as planar arrays and thereby enable their 3D beamforming towards the ground, as shown in Fig.\;\ref{tech1}(a). Different from the traditional two-dimensional (2D) beamforming (e.g., with linear array) by ground BSs which controls the beam pattern only in the azimuth plane, 3D beamforming is able to adjust the beam pattern in both elevation and azimuth planes with more refined beam resolution, which can be leveraged to null the legitimate user's signal in the directions of the eavesdroppers more effectively.
Furthermore, the dominant LoS A2G channel also makes transmit beamforming more efficient as compared to that in terrestrial channels with more complicated multi-path effects due to the rich scattering environment. As shown in Fig.\;\ref{tech1}(a), even equipped with 3D beamforming capability, a ground BS may lack sufficiently separated elevation angles with the legitimate user and the eavesdropper due to its low altitude above the ground; as a result, it cannot avoid information leakage if the two nodes are located in the same direction from it, even though they have different horizontal distances with it. In contrast, a UAV transmitter at a much higher altitude can still exploit the ground nodes' different elevation angles ($\theta_1$ and $\theta_2$) with it to cancel the legitimate user's signal at the eavesdropper by e.g., zero-forcing (ZF) based precoding. Furthermore, 3D beamforming can be jointly employed with transmitting AN to jam or interfere with the ground eavesdroppers and thereby achieve better security \cite{wu2016secure}. Last but not least, the wiretap channels of the  ground eavesdroppers can be practically obtained with high accuracy at the UAV if their locations are known and thus their dominant LoS channel components can be resolved. This further enhances the robustness of 3D beamforming/jamming by the UAV.
\begin{figure}[!t]
\centering
\includegraphics[width=0.65\textwidth]{tech12.pdf}
\caption{Approaches for safeguarding legitimate A2G communications against terrestrial eavesdropping.} \label{tech1}
\end{figure}

\subsubsection{\bf Joint UAV Positioning/Trajectory and Communication Design}
As compared to terrestrial nodes, UAVs can be deployed more flexibly and move more freely in 3D space, which can be utilized to improve the communication throughput \cite{JR:wu2017joint,zeng2016throughput} or secrecy rate \cite{zhang2018securing} with the ground users.
For example, a UAV transmitter at a given horizontal position can avoid signal blockage with legitimate users due to a high-rise building by increasing its altitude, or incur blockage and thus more propagation loss with the ground eavesdropper by lowing its altitude (see Fig.\;\ref{tech1}(b)) so as to achieve better secrecy performance. On the other hand, for flying UAVs, a proper trajectory design of the UAV helps shorten its communication distances with legitimate users while enlarging the wiretap channel distances to alleviate the ground eavesdropping. Furthermore, communication scheduling and power control can be jointly designed with the UAV trajectory to achieve the optimum secrecy performance. For example, as shown in Fig.\;\ref{tech1}(b), when a UAV has to fly over some eavesdroppers due to its mission requirement, the legitimate channel may not be always stronger than the wiretap channel along its trajectory.
In this case, the UAV can increase its transmit power/rate when it flies closer to a legitimate user for communication while being further away from the eavesdroppers, and decrease transmit power/rate or even stop transmission otherwise. Finally, the UAV positioning/trajectory and communication co-design can be adaptive in real time to deal with moving eavesdroppers on the ground that intend to improve the eavesdropping performance and/or hide themselves from getting exposed.

\subsubsection{\bf Multi-UAV Cooperation}
In practice, a single UAV only has limited communication and maneuvering capability and thus may not achieve the desired secure communication performance in some challenging scenarios (e.g., with multiple collusive eavesdroppers over a large  area). This thus motivates the deployment of multiple collaborative UAVs to achieve more efficient secure communications  \cite{JR:wu2017joint,zhong2018secure}.
\begin{comment}
or a swarm of To maximally improve the secrecy performance, different UAVs may be assigned with different roles.by leveraging e.g., high-speed wireless fronthaul. via direct UAV-UAV communications
\end{comment}
For example, as shown in Fig.\;\ref{tech1}(c), based on the locations of eavesdroppers, the ground users can be grouped into different  clusters with each cluster served by a single UAV. As such, it may not be necessary for the UAV to fly over the eavesdroppers as in Fig.\;\ref{tech1}(b), which thus helps reduce the information leakage. Alternatively, some UAVs may act as aerial jammers  \cite{JR:Anli18WCL,zhong2018secure}, which are deployed above a group of nearby eavesdroppers to degrade their signal reception by sending AN signals \cite{JR:Anli18WCL,zhong2018secure}. This in turn provides higher flexibility for the deployment or trajectory design of other UAV transmitters to achieve better secrecy communication performance.

\subsection{Securing G2A Communication Against Terrestrial Jamming}
In contrast to legitimate ground receivers in A2G communications, legitimate UAV receivers in G2A communications (see Fig.\;\ref{UAV_Legi}(c) and Fig.\;\ref{UAV_Legi}(d)) are exposed in the sky and thus more vulnerable to jamming attacks by ground adversary nodes that send AN signals to interfere with the UAV so as to reduce its receive signal-to-interference-plus-noise ratio (SINR) for decoding.
\begin{comment}
In practice, jammers can be either proactive or reactive. While a proactive jammer does not sense for ongoing legitimate transmissions but jams the channel persistently, a reactive jammer transmits interfering signals only when it detects the existence of legitimate transmissions, which is thus more energy-efficient and less likely to be detected. To achieve reliable transmission, the counter strategies generally include using spread spectrum techniques as well as those improving the SINR.
\end{comment}
Although DSSS and FHSS are two widely applied anti-jamming techniques, they usually result in low spectrum efficiency and may not be sufficient to deal with the strong G2A jamming due to the LoS-dominant channel. Fortunately, the techniques proposed in the preceding subsection can also be applied to cope with the malicious jamming more effectively. For example, the 3D receive beamforming at a legitimate UAV can provide higher spatial resolution than 2D beamforming (with a linear array) and thus more efficiently cancel the interference from the ground jammers, as long as they are not so close to the legitimate ground transmitters so that their elevation angles with the UAV can be resolved. In the case that the legitimate ground transmitter is a multi-antenna BS as shown in Fig.\;\ref{UAV_Legi}(d), the 3D transmit and receive beamforming can be jointly employed at the BS and UAV to maximally improve the legitimate link SINR. Moreover, the 3D high mobility of UAV can be leveraged to optimize its position or trajectory by shortening/enlarging the distances from legitimate transmitters/malicious jammers. Finally, in the case that the UAV receiver needs to move away from the jammed area to avoid the strong interference, the device-to-device (D2D) communications among the legitimate ground nodes (e.g., the ground sensors shown in Fig.\;\ref{UAV_Legi}(c)) can be exploited to forward the data of the ground nodes in the jammed area via a separate (unjammed) channel to other nodes that are sufficiently far away from the jammer but close to the UAV position/trajectory for more reliable uploading to the UAV.

\section{Safeguarding Terrestrial Network Against Malicious UAV Attacks}
\begin{comment}
another effective countermeasure
With the recent standardization progress of supporting UAVs with cellular networks, it is anticipated that more and more UAVs will be connected to 5G-and-beyond wireless networks.
\end{comment}
The LoS-dominant air-ground channels also make the terrestrial communications highly exposed to malicious UAVs. For example, a single UAV eavesdropper or jammer is capable of wiretapping or contaminating the transmissions within multiple cells. In addition, the flexible mobility and deployment of UAVs may be harnessed by malicious users to launch more aggressive attacks to terrestrial networks. It is therefore of paramount importance to develop advanced countermeasures for combating such attacks by malicious UAVs.

\subsection{Anti-UAV Eavesdropping of Terrestrial Communication
\begin{comment}
Similarly, all defense strategies to improve the secrecy performance should aim to maximize SINR difference at legitimate ground receivers and malicious UAV eavesdroppers.
\end{comment}
As shown in Fig.\;\ref{UAV_Mali}(a), thanks to the different altitudes of a legitimate ground receiver and a UAV eavesdropper, their well separated elevation angles can be exploited by the 3D transmit beamforming at a ground BS to achieve efficient signal nulling and hence secure  communication in the downlink. However, for uplink transmissions where the legitimate ground users may not be equipped with a large  number of antennas, how to combat against malicious UAV eavesdropping becomes more challenging. To resolve this issue, we propose two promising techniques based on different cooperation mechanisms in the terrestrial network, namely multi-hop D2D relaying and cooperative remote jamming, elaborated as follows.
\begin{figure}[!t]
\centering
\includegraphics[width=0.7\textwidth]{UAV_Mali_tech3.pdf}
\caption{Approaches for safeguarding terrestrial communications against UAV eavesdropping.}\label{tech2}
\end{figure}

\subsubsection{\bf Multi-hop D2D Relaying}
Given the LoS advantage of G2A channels over terrestrial fading channels, the transmitted signal that needs to be received reliably at a ground receiver in terrestrial communication is more likely to be wiretapped by the UAV under the same link distance. To overcome this difficulty, one viable solution is to adopt multi-hop D2D relaying in the terrestrial network as shown in Fig.\;\ref{tech2}(a), i.e., utilizing nearby ground nodes as legitimate relays to help forward the confidential message from a source (legitimate ground node) to its destination (ground BS). By judiciously selecting the relay nodes (e.g., those with blocked LoS paths with the UAV) and designing their cooperative transmission (e.g., transmit power control and distributed space-time coding), the spatial/multi-path diversity gain can be reaped to significantly enhance the secrecy communication rate.
\begin{comment}
and even the distributed beamforming gain, which can greatly increase the communication reliability. On the other hand, the shortened transmission distance in each hop allows the relay node to reduce the transmit power,  which could drastically decrease the probability of being detected and hence eavesdropped by malicious UAVs.
\end{comment}
However, on the other hand, an excessive number of hops may entail longer end-to-end delay as well as lower spectrum efficiency, and also increase the risk of being eavesdropped due to the more  exposure. As such, the number of hops, the transmit power in each hop, and the relay selection as well as the relay protocol need to be jointly designed to reconcile the above trade-off.

\subsubsection{\bf Cooperative Remote Jamming}
Another effective means to prevent UAV eavesdropping is the active jamming. However, its major weakness is that both the wiretap link and the legitimate link are degraded by the AN signal. To mitigate such adverse effect, a promising approach is by employing terrestrial nodes (e.g., ground BSs) at favorable locations for cooperative {\it remote} jamming (e.g., those who are close to the UAV eavesdroppers but distant from the legitimate ground receivers of interest). As shown in Fig.\;\ref{tech2}(b), a set of cooperative terrestrial BSs that are far away from the cell of the legitimate ground user  but still have strong LoS channels with the UAV eavesdropper, are selected as  jammers. As such, the signal reception at the UAV eavesdropper is dramatically degraded whereas only negligible interference is perceived at the legitimate ground receiver (serving BS) due to the more severe propagation loss and multi-path fading over terrestrial channels. In practice, depending on the traffic load  of the terrestrial network, the jamming signals of the cooperative ground BSs can be either random  message signals sent to their respective ground users or AN signals for dedicated jamming.

\subsection{Anti-UAV Jamming to Terrestrial Communication}
The jamming attack by malicious UAVs is another challenging issue for safeguarding terrestrial communications against their strong LoS jamming. Fortunately, for uplink transmissions (see Fig.\;\ref{UAV_Mali}(b)), the receive 3D beamforming can be applied at ground BSs to mitigate the UAV jamming signals efficiently. Furthermore, a neighboring idle BS that receives the jamming signal can also forward it via the high-speed backhaul link to the serving BS for cooperative jamming signal cancellation. While for downlink transmissions with single-antenna legitimate receivers (e.g., sensors), the above approaches become infeasible and anti-UAV jamming becomes more challenging.
\begin{comment}
In addition, even for mobile receivers (e.g., ground vehicles), the UAV's jamming effect cannot be simply mitigated by escaping away from the jammed area, since it requires a very long duration to completely escape. In light of this, one
\end{comment}
One possible solution for this scenario  is again by exploiting the potential cooperative D2D communication of the ground nodes and their probabilistic LoS channels with the UAV jammer. Specifically, the terrestrial message is first sent to a ground node that is near to the legitimate user and also has blocked LoS link with the UAV jammer for reliable decoding. Then, the decoded message is forwarded to the legitimate user over an unjammed channel via D2D communication.

\section{Numerical Results and Discussion}
Numerical results are provided in this section to demonstrate the effectiveness of some of the proposed techniques in the previous two sections, respectively.

\subsection{UAV-assisted Jamming}
\begin{figure}[t]
\centering
\subfigure[Simulation setup.]{\includegraphics[width=0.45\textwidth]{UAV_coop_jammer_setup.pdf}}
\subfigure[Secrecy rate versus UAV transmit power.]{\includegraphics[width=0.48\textwidth]{multiUAV_coop3.eps}}
\caption{Secrecy rate performance with UAV-assisted jamming.}\label{multiUAV:coop} \vspace{-4mm}
\end{figure}
First, we show the performance of the proposed UAV positioning design with cooperative jamming against the terrestrial eavesdropping. As shown in Fig.\;\ref{multiUAV:coop}(a), we consider that a ground eavesdropper intends to wiretap the A2G communication from a UAV-BS to a ground user. To  improve the secrecy rate, a cooperative UAV-jammer is deployed right above the eavesdropper with the fixed jamming power given by $P_{J}=5$ dBm. The two UAVs are deployed at the altitude of $200$ meters (m) such that the A2G links approximately follow the free-space path loss model \cite{JR:wu2017joint}. The reference signal-to-noise ratio (SNR) at the distance of $1$ m is set as $80$ dB.
\begin{comment}
As the secrecy rate in the considered setup mainly depends on the UAV-BS's location
\end{comment}
For comparison, we consider the following three schemes: 1) optimized UAV-BS location with jammer; 2) optimized UAV-BS location without jammer; and 3) UAV-BS location fixed  at $(0, 200)$ m without jammer, i.e., right above the legitimate user. The optimal UAV-BS locations in 1) and 2) are obtained by one dimensional search. In Fig.\;\ref{multiUAV:coop}(b), we plot the achievable secrecy rates in bits per second per Hertz (bps/Hz) by different schemes versus the UAV-BS transmit power, $P_U$. First, it is observed that compared to hovering at the fixed location, optimizing the UAV-BS locations achieves significant secrecy rate gains, even without the helping UAV-jammer. This is because the UAV-BS can enlarge the difference in the receive SNRs of the legitimate user and the eavesdropper by properly positioning itself away from both of them, i.e., hovering at the left side of $(0, 200)$. Second, one can observe that for the low UAV-BS transmit power, deploying the UAV-jammer with fixed jamming power degrades the secrecy rate; while for the high UAV-BS transmit power, the secrecy rate is substantially improved as compared to the other two schemes. This is expected as in the former case, the eavesdropping rate is quite small and hence deploying a dedicated jammer has a more detrimental effect on the legitimate receiver than the eavesdropper. While in the latter case where the eavesdropping rate becomes high, deploying the UAV-jammer in the eavesdropper's vicinity can effectively disrupt its signal reception. Meanwhile, such a cooperative UAV-jammer also enables the UAV-BS to move closer to the user for improving the secrecy rate. For example, when $P_{U}=15$ dBm, the optimized locations of the UAV-BS without and with the UAV-jammer are $(-100, 200)$ m and $(-89.5, 200)$ m, respectively. It is worth pointing out that instead of fixing the UAV-jammer's power and location in this example, we can optimize them along with the UAV-BS's location  so that the secrecy rate with the UAV-jammer is always better than that without using it \cite{zhang2018securing,zhong2018secure,JR:Anli18WCL}.

\subsection{Cooperative Remote Jamming}
\begin{figure}[hbtp]
\centering
\subfigure[Cellular network setup.]{\includegraphics[width=0.5\textwidth]{sys_remote.eps}}
\subfigure[Secrecy rate versus number of cooperative BSs.]{\includegraphics[width=0.48\textwidth]{remote_jam.eps}}
\caption{Secrecy rate performance with remote jamming.}\label{remote_jam}\vspace{-12pt}
\end{figure}
Next, we evaluate the performance of the proposed cooperative remote jamming technique for protecting the terrestrial communication from the UAV eavesdropping. As shown in Fig.\,\ref{remote_jam}(a), we consider a cellular network where a UAV eavesdropper aims to overhear the downlink communication in a specific cell. To combat against the UAV's strong LoS eavesdropping link, multiple remote BSs send independent jamming signals to the UAV eavesdropper. In Fig.\,\ref{remote_jam}(b), we plot the secrecy rate of the legitimate communication versus the number of cooperative BSs with  different BS transmit power, denoted as $P_T$. The cell radius is assumed to be 800 m. The UAV altitude is set to be 200 m, and the horizontal distance between the UAV and the legitimate BS transmitter is set to be 1000 m. The location of the legitimate receiver is randomly generated in the cell of the serving BS. The ground user and the UAV are assumed to be equipped with an isotropic antenna, while each BS employs a fixed antenna pattern with 10-degree downtilt angle \cite{zeng2018cellular}. From Fig.\,\ref{remote_jam}(b), it is observed that the secrecy rate is approximately  zero when there are no cooperative jamming BSs, regardless of the BS's transmit power. This is because the G2A eavesdropping channel with dominant LoS is much better than the terrestrial fading channel of the legitimate link, thus leading to virtually zero secrecy rate. However, the secrecy rate is observed to increase rapidly as the number of cooperative BSs increases. This is expected since the UAV eavesdropper suffers from more interference from the increasing number of jamming BSs, which only cause negligible interference to the legitimate receiver due to the more severe pathloss and multi-path fading over the terrestrial channels. Furthermore, it is observed that increasing $P_T$ further improves the secrecy rate, despite that the achievable rates of both the legitimate link and the UAV eavesdropping link increase with $P_T$. The reason lies in that the former increases much faster than the latter, which is severely limited by the increasing G2A interference with higher $P_T$.

\section{Conclusions and Future Work}
\begin{comment}
we identify two new research paradigms in future wireless network with UAVs from the
\end{comment}
In this article, we focus on addressing two new and challenging security issues arising from the LoS-dominant UAV-ground channels in future wireless networks, from a PHY design perspective. We propose promising approaches and techniques to achieve secure UAV-ground communications in the presence of terrestrial eavesdroppers/jammers, as well as secure terrestrial communications against the malicious eavesdropping/jamming attacks by the UAVs. In future wireless networks, legitimate aerial and terrestrial communications both need to be protected against sophisticated attacks that may involve collusive eavesdroppers and jammers on the ground as well as over the air. As such, the proposed approaches in this article need to be further investigated to address more challenging scenarios in practice. Besides the PHY security issues discussed in this article, there are other related problems that need to be addressed, as listed below to motivate future work:
\begin{itemize}
  \item {\it UAV-assisted terrestrial adversary detection:} With the rapid advances of high-definition optical cameras, it would be an appealing solution to deploy them on UAVs for detecting, identifying, and tracking  malicious nodes on the ground. In addition, thanks to the high altitude and favorable LoS-dominant link, the UAV generally has a stronger spectrum sensing ability than conventional terrestrial nodes. Thus, the UAV can also help achieve more accurate terrestrial jamming detection.
  \item {\it Malicious UAV detection:} When there are malicious UAVs, it is imperative to detect their presence and also track their moving locations. For active UAVs such as UAV jammers, they can be detected/localized by using conventional signal sensing and ranging techniques. Whereas for passive UAVs such as UAV eavesdroppers, more sophisticated detection methods are generally required, such as radar and computer vision based methods\cite{guvencc2017detection}.
  \item {\it UAV spoofing:} In addition to the eavesdropping and jamming attacks, legitimate UAV communications may suffer from a more advanced spoofing attack such as the global position system (GPS) spoofing\cite{kerns2014unmanned}, where the ground adversary aims to deceive the UAV's navigation system by sending counterfeit GPS signals.
  \item {\it UAV-aided wireless surveillance and intervention:} By exploiting its high mobility and LoS channel with the ground, UAV can  act as a legitimate eavesdropper or jammer for realizing the surveillance of suspicious communications or intervention of malicious communications on the ground, respectively\cite{xu2017surveillance}. In contrast to the techniques proposed in this article  to defend against UAV eavesdropping and jamming, the objective is reversed in this case  to design efficient eavesdropping and jamming schemes for legitimate UAVs.
\end{itemize}
\begin{comment}
By launching this attack, the target UAV may deviate from its pre-determined trajectory and fly towards a target point designated by the spoofer for illegal purposes. Potential countermeasures include cryptographic signal authentication techniques and antenna-oriented non-cryptographic techniques, which requires antenna movement to differentiate between spoofed and authentic signal sets\cite{kerns2014unmanned}.
\end{comment}

\bibliographystyle{IEEEtran}