\section{Introduction}

After the experimental realization of Bose-Einstein condensate (BEC) of
$^{52}$Cr \cite{pfau,4}, $^{164}$Dy \cite{7}, and $^{168}$Er \cite{8}
atoms with large magnetic dipolar interaction, there have been renewed interest in the study of  cold atoms,
both theoretically and experimentally. The anisotropic long-range dipolar interaction acting in all partial waves in these atoms
is basically different from isotropic short-range  S-wave interaction acting in nondipolar atoms.  Polar bosonic molecules with much larger
permanent  electric
dipole moment  are also candidates  for possible BEC
experiments \cite{bosmol}.
Dipolar BECs have novel distinct properties. Due to the anisotropic  dipolar interaction, the stability of a
dipolar BEC depends on the scattering length as well as
the trap geometry \cite{4,12,14}.
The shock and sound waves also propagate anisotropically in a dipolar BEC
\cite{17,17A}. Anisotropic collapse has been observed and studied
in a dipolar BEC of $^{52}$Cr atoms \cite{18}.
Anisotropic rotons \cite{roton} and anisotropic critical superfluid
velocity \cite{crit}
have been suggested and studied in dipolar BECs. Distinct stable checkerboard, stripe, and star configurations in dipolar BECs
have been identified in a two-dimensional (2D) optical lattice  as stable Mott insulator \cite{capo} as well as superfluid soliton \cite{star}
states. Anisotropic solitons in 2D have also been suggested in dipolar
BECs \cite{anisol}. A new possibility of studying universal properties of dipolar BECs at unitarity has been suggested \cite{uni}.

After the realization of BEC of alkali-metal atoms, degenerate
nondipolar gas of fermionic $^6$Li \cite{li}, $^{40}$K \cite{k}, and
$^{87}$Sr \cite{sr} atoms were observed. Later, superfluid states of
paired $^6$Li \cite{hulet} and $^{40}$K \cite{jin} Fermi atoms have been
investigated. More recently a degenerate dipolar gas of fermionic
$^{161}$Dy atoms with large magnetic dipole moment has been created and
studied \cite{dipf}. Realization of quantum degeneracy in $^{161}$Dy
atoms should be considered as a doorway for the study of anisotropic
superfluidity in dipolar fermions. Fermionic polar molecules, such as
$^{40}$K$^{87}$Rb, of large permanent electric dipole moment are also
considered for this purpose \cite{fermol}. { The $^{40}$K$^{87}$Rb molecule in the
singlet rovibrational ground state has an electric dipole moment of 0.6 Debye,
thus leading to a dipolar interaction larger than in the case of $^{161}$Dy
atoms by more than an order of magnitude \cite{pfau}.}
One advantage of studying the
effect of dipolar interaction in a degenerate dipolar Fermi gas over
that in a dipolar BEC is the remarkable stability of the degenerate
Fermi gas. A BEC is usually fragile or short-lived for many experiments
due to three-body loss through molecule formation. The three-body loss
is highly suppressed in a degenerate Fermi gas due to Pauli repulsion
among identical fermions. On the other hand a theoretical study of a BEC
is much simpler than that of a degenerate Fermi gas \cite{frmp} due to
the existence of a simple order parameter and a simple mean-field
Gross-Pitaevskii equation in the former case.

A microscopic description of a degenerate Fermi gas is complicated due to the difficulties with the antisymmetrization of a many-fermion system and finding an appropriate simple order parameter. Drastic approximations are often necessary to achieve this goal for fermions. There  have been several different theoretical descriptions for a degenerate nondipolar Fermi gas \cite{dfer}.
There have also been a few studies of the degenerate dipolar Fermi gas
 employing different types of approximations \cite{fermidipth}.  However, for a description of some macroscopic observables of a degenerate Fermi gas, a simple hydrodynamic description \cite{hydro}
 often seems appropriate which does not require a precise antisymmetrization of the dynamics \cite{frmp,st,lca}. The effect of antisymmetrization is included approximately via an appropriate energy or Lagrangian density. A minimization of the energy leads to a well-founded variational approximation \cite{var}. An Euler-Lagrange equation derived from such a Lagrangian provides further improvement over the variational results. Here we present a description of the degenerate spin polarized Fermi gas based on such a hydrodynamic model \cite{hydro,st}.

The present theoretical formulation for the degenerate dipolar Fermi gas starts in section \ref{IIA}
with the standard equations of classical hydrodynamics \cite{hydro} with the appropriate equation of state including the kinetic energy of fermions filling the Fermi sea and the dipolar interaction among them.  The equivalent of the Thomas-Fermi (TF) approximation for the dipolar gas is obtained by setting the velocity field equal to zero in the hydrodynamic equations. An equivalent classical energy   density is written using the local-density approximation (LDA) \cite{frmp}.  A quantum pressure term is then introduced in the LDA energy density. With such a quantum pressure term, a nonlinear Schr\"odinger-type equation is derived for the density of fermions. For a moderate number of fermions (greater than 100 or so), the quantum pressure term is negligibly small and hence has insignificant effect on the result. {  However, the inclusion of the quantum
pressure term allows one to write a dynamical equation for fermions in the form of a nonlinear Schr\"odinger
equation
to study the dynamics. The LDA or the TF
approximation, on the other hand, allows only the investigation of static properties of the system. For example, these approximations cannot be used to study the sound propagation dynamics in fermions as in section
\ref{IIIB}.}
  In section \ref{IIB},  we present a Gaussian variational approximation for the problem described by the LDA energy density.
In section \ref{IIC},  simplified models are derived in reduced dimensions, appropriate for cigar and disk shapes of the degenerate dipolar gas when there is a strong trap in radial or in axial directions, respectively. In section \ref{IID}, using the hydrodynamic model, we obtain an analytic expression for the anisotropic sound velocity in the degenerate dipolar gas.
In section \ref{IIIA}, we present numerical results for stationary properties $-$ shape, size, and chemical potential $-$ of a trapped degenerate Fermi gas of $^{161}$Dy atoms, and compare with results from appropriate models in reduced dimensions and variational approximation. In  section \ref{IIIB}, the anisotropic sound propagation in an infinite dipolar degenerate gas is demonstrated numerically  and the velocities so obtained are compared with the analytical results of section \ref{II} D.
Finally, in section IV we present a
summary of our study.

\section{Analytical Consideration}
\label{II}
\subsection{Hydrodynamic Model}
\label{IIA}

The normal one-component Fermi gas can be in the collisionless regime where collision is rare or in the collisional hydrodynamic
regime where frequent collision due to the dipolar interaction  allow the system to settle to local equilibrium, where
the system can be described by simple hydrodynamic equation rather than a detailed multi-particle description.
Here we consider the system in such a configuration. A semi-quantitative estimate for the validity of a hydrodynamic description
is given by the condition that relaxation time $\tau_R$ is small compared to the time scale $1/\omega$
defined by average trap frequency $\omega$  \cite{17,17A}.  For contact interaction this condition can be expressed
in terms of the scattering length to measure the strength of atomic interaction \cite{17B}.  The
strength of dipolar interaction is usually measured in terms of the convenient
length  $L_d\equiv 3a_{dd}= \mu_0\bar \mu^2 m/(4\pi \hbar^2)$, where $\bar \mu $ is the magnetic moment of an atom of mass $m$ and  $\mu_0$ the permeability of free space.
Using this length scale to measure dipolar interaction, the condition for the
validity of hydrodynamic description can be expressed as \cite{17A,17B}
\begin{equation}\label{cond}
(\omega  \tau_R)^{-1}
\approx   \left( N^{1/3}\frac{L_d}{a_{\mathrm{ho}}}\right)^2F\left( \frac{T}{T_F}\right)  > 1,
\end{equation}
where $N$ is the number of atoms,
$a_{\mathrm{ho}}=\sqrt{\hbar/(m\omega)}$ the oscillator length, $T$  the temperature, $T_F$ the Fermi temperature and
$F(x)$  an universal function and can be taken of the order of
unity \cite{17B} under experimental conditions. Hence for sufficently
large  $N$ and/or sufficiently small   $a_{\mathrm{ho}}$
the system should enter the hydrodynamic regime. For $^{161}$Dy atoms of magnetic moment $10\mu_B$ with $\mu_B$ the Bohr magneton, the length $L_d\equiv 3 a_{dd}\approx 3\times 130 a_0$ where $a_0$ is the Bohr radius and for
an oscillator length $a_{\mathrm{ho}}=0.25 $ $\mu$m and $N=5000$ the system is already in the hydrodynamic regime.  For polar Fermi molecule
$^{40}$K$^{87}$Rb of electric dipole moment $d= 0.566$ Debye in the singlet
rovibrational ground state \cite{fermol}, the length
 $L_d =md^2/(4\pi\hbar^2\epsilon_0)\approx 3 \times 2000 a_0$ \cite{pfau}, where $\epsilon_0$ is the permitivity of free space,
 and the system should enter
the hydrodynamic regime with larger $a_{\mathrm{ho}} = 1$  $\mu$m (weaker trap)
and much smaller $N $ of few hundreds.

At sufficiently low temperature the    normal Fermi gas enters
the degenerate phase.
Most macroscopic properties, like shape, density, chemical potential, sound propagation etc., of this system can be described
by the  Landau hydrodynamical equations \cite{hydro}.
A velocity field
is introduced as the gradient of a velocity potential $\Phi$ of flow, usually related to phase of the order parameter,  by
$
{\bf v}({\bf r},t)=   \nabla \Phi({\bf r},t),
$
subject to the irrotational condition
$\nabla \times {\bf v}({\bf r},t)=0$.
The continuity  and the flow equations are then given by
\cite{hydro}
\begin{eqnarray}\label{h1}
&&\frac{\partial n({\bf r},t)}{\partial t}+\nabla \cdot [n({\bf r},t) \bf v({\bf r},t) ]=0,\\
&& m\frac{\partial \bf v(\bf r,t)}{\partial t}+\nabla\left[
\frac{ m {{\bf v(r},t)}^2}{2}+\mu(n,{\bf r})+
V({\bf r})\right]=0,\label{den}
\end{eqnarray}
where $n({\bf r},t)$ is density at space point $\bf r$ and time $t$,
$V(\bf r)$ is an external trap, usually taken
as
\begin{equation} \label{pot}
V({\bf r})=\frac{1}{2}m\omega^2 (\nu^2 \rho^2+\lambda^2 z^2),
\end{equation}
with $\omega$ the trap frequency and $\nu$ and $\lambda$ are anisotropy parameters.
The bulk chemical potential
$\mu(n, \bf r)$ is   determined by the equation of state of the
uniform Fermi gas of dipolar atoms and is given by
\begin{eqnarray}\label{chemx}
\mu(n,{\bf r})=
\frac{\hbar^2 [6\pi^2 n({\bf r})]^{2/3} }{2m} +
\int d{\bf r'} n({\bf r'}) U_{dd}({\bf r-r'}),
\end{eqnarray}
where the first term on the right-hand side (rhs) is the Fermi energy $E_F\equiv{\hbar^2
(6\pi^2 n)^{2/3} }/({2m})
$ \cite{frmp}
and the last
term is the contribution from dipolar interaction energy \cite{st}, with $U_{dd}(\bf r -r ')$
is the dipolar potential.

The hydrodynamic description is
valid for a macroscopic
observable with its  characteristic excitation wave-length $\lambda$
much larger than the healing length \cite{frmp}. A safe condition to satisfy this criterion is to take the wave length to be much larger than
de Broglie wavelength at the
Fermi surface, i.e., \cite{frmp,am}
\begin{equation}\label{con}
\lambda  \gg 2\pi/k_F
\end{equation}
with Fermi momentum $k_F$ defined by
$E_F=\hbar^2k_F^2/2m$.
{ Truly speaking, a degenerate Fermi
gas may not be fully irrotational and allows rotational components
in the velocity field, not allowed in the present
hydrodynamic formulation. This fact should influence only the rotational
properties \cite{frmp} of the degenerate Fermi gas not
considered in this paper. We also assume the absence of any velocity dependent frictional force.}

An approximate TF profile for  density  can now be obtained by setting  velocity   ${\bf v}=0$ in   (\ref{den}), when \cite{frmp}
\begin{equation}\label{tf}
\mu(n,{\bf r})+V({\bf r})=\mu_0^{\mathrm{TF}},
\end{equation}
where $\mu_0^{\mathrm{TF}}$ is the chemical potential of the trapped gas.  When   (\ref{tf}) is solved
subject to the appropriate normalization condition, we obtain both the chemical potential $\mu_0^{\mathrm{TF}}$ and  the density $n(\bf r)$.

There is an equivalent description of the trapped degenerate Fermi gas in the  LDA, based on the assumption that, locally the dipolar Fermi gas would behave like a uniform gas, so that the energy density can be written as the energy of the uniform system times the local density \cite{frmp}.
For the degenerate dipolar Fermi system,  the classical energy
density per particle  is given by \cite{fermidipth,st}
\begin{eqnarray}
 \label{lagcl} {\cal E}_{\mathrm{cl}}({\bf r})&&= \frac{3}{5}
\frac{\hbar^2}{2m}(6\pi^2N)^{2/3} [n({\bf r})]^{5/3}+V({\bf r}) n({\bf r})\nonumber \\
&&+N\frac{1}{2}\int d{\bf r'}n({\bf r}) n({\bf r'}) U_{dd}({\bf r-r'}), \end{eqnarray}
where $n(\bf r)$ is the density of fermions normalized as $\int d{\bf r}n({\bf
r}) =1$.
The first term on the rhs of   (\ref{lagcl}) is the total
kinetic energy of the spin polarized fermions filling all levels up to the Fermi sea, the second
term is energy in the trap, and the third term describes the dipolar interaction. A minimization of the classical energy (\ref{lagcl})
leads to the TF condition (\ref{tf}).

As the degenerate Fermi gas is a quantum system, a quantum pressure term  when included in   (\ref{lagcl})
yields the following expression  for net energy density
\begin{eqnarray}\label{grad}
\label{lag} &&{\cal E}({\bf
r})= \frac{\hbar^2}{8m}|\nabla_{\bf r}\sqrt{n(\bf r)}|^2 +{\cal E}_{\mathrm{cl}}({\bf r}).
 \end{eqnarray}
This gradient correction term \cite{grad} to the TF energy density (\ref{lagcl})
takes into account the additional kinetic
energy due to spatial variation of density (near the surface).
Such a surface term was first considered by von Weisz\"acker \cite{von,toigo,kohn} in the description  a large nuclei.
Previous descriptions of a degenerate Fermi gas considered different
coefficients in this term
\cite{dfer,grad}.
The energy density (\ref{lag}) has
successfully used in many problems of Fermi gas
\cite{toigo,kohn,as,many}.

With the  Lagrangian  density ${\cal L}({\bf r})={\cal E}({\bf r})-\mu_0n(\bf r)$
the Euler-Lagrange equation
is given by  \cite{as}
\begin{eqnarray}  \label{gp3d}
 \mu_0 \sqrt{n{(\bf r)}}
&& =  \biggr[ -\frac{\hbar^2\nabla^2 }{8m} + V({\bf r})
+
\frac{\hbar^2}{2m}[6\pi^2 N n({\bf r})]^{2/3}\nonumber \\
&&+  N\int U_{dd}({\bf r -r'})n({\bf r'})d{\bf r'}
\biggr] \sqrt{n({\bf r})},
\end{eqnarray}
with $\mu_0$ the chemical potential.
The derivarive term in   (\ref{gp3d}) term contributes much less than the
dominant ``Fermi energy" term $\hbar^2[6\pi^2Nn({\bf r})]^{2/3}/(2m)$ and its neglect
 leads to the  TF relation (\ref{tf}). The   dipolar interaction in   (\ref{gp3d}) is taken as
$
 U_{dd}({\bf R}) = 3
a_{dd}\hbar^2(1-3\cos^2	\theta)$ $/(mR^3),$ $\quad {\bf R=r-r'}
$ where $	\theta$ is the angle between the $\bf R$ and the polarization
direction $z.$

The condition (\ref{cond}) refers to the validity of a
hydrodynamic  description of the system ($\bf v \ne 0$) \cite{17A,17B}. A reliable stationary description obtained
using   (\ref{gp3d}) for $\bf v =0$,
of density and other macroscopic properties, such as sound
velocity, as considered in this paper,
can be obtained for a smaller number of fermions
consistent with condition (\ref{con}) \cite{frmp,am}, provided the contribution of the kinetic energy term $\hbar^2[6\pi^2N n({\bf r})]^{2/3}
/(2m)$ in   (\ref{gp3d}) is much larger than that of  other terms.

It is convenient to write a dimensionless form of   (\ref{gp3d})
with the potential (\ref{pot})
as
\begin{eqnarray}  \label{3d1}
 \mu_0 \sqrt{n{(\bf r)}}
&& =  \biggr[ -\frac{\nabla^2 }{8} + \frac{1}{2}(\nu^2\rho^2+\lambda^2z^2)
+
\frac{1}{2}[6\pi^2 N n({\bf r})]^{2/3}\nonumber \\
&&+  3a_{dd}N\int
\frac{1-3\cos^2 	\theta }{R^3}
n({\bf r'})d{\bf r'}
\biggr] \sqrt{n({\bf r})},
\end{eqnarray}
where energy, length and density are expressed in units of
$\hbar\omega$, $l_0=\sqrt{\hbar/m\omega}$ and $l_0^{-3}$.

\subsection{Variational Approximation}
\label{IIB}

  The energy density corresponding to the dimensionless equation (\ref{3d1})
 can be written as
\begin{eqnarray}
{\cal E}({\bf r})&&= \frac{1}{8}|\nabla_{\bf r}\sqrt{n(\bf r)}|^2
+\frac{3a_{dd}N}{2}\int d{\bf r'}n({\bf r})
n({\bf r'}) \frac{1-3\cos^2 	\theta }{R^3}
\nonumber \\
&&+\frac{1}{2}(\nu^2\rho^2+\lambda^2 z^2)n({\bf r})+\frac{3}{10}
(6\pi^2N)^{2/3}
[n({\bf r})]^{5/3}
.
\end{eqnarray}
A variational approximation for the problem  can be obtained with the following Gaussian ansatz for density \cite{var}
\begin{eqnarray}
n({\bf r})= \frac{1}{\pi^{3/2} w_\rho^2 w_z}  \exp\left[-\frac{\rho^2}{w_\rho^2} -\frac{z^2}{w_z^2}  \right],
\end{eqnarray}
where $w_\rho$ and $w_z$ are the variational widths along radial $\rho$ and axial $z$
directions.
With this density,
the effective energy per particle  of the system
$E=\int d{\bf r}{\cal E}({\bf r})$
is
\begin{eqnarray}\label{energy}
E&&= \frac{1}{8}\left[\frac{1}{w_{\rho}^2}+
\frac{1}{2w_{z}^2}
\right]
+\left[\frac{\nu ^2 w_{\rho}^2}{2}+\frac{\lambda^2 w_z^2}{4}\right]
\nonumber \\
&&-\frac{Na_{dd}f(\kappa)   }{\sqrt{2\pi} w_{\rho}^2w_{ z}}
+ \sqrt{\frac{3}{5}}  \frac{9 (6\pi^2  N)^{2/3} }{50
\pi w_{\rho}^{4/3} w_{z}^{2/3}},
\end{eqnarray}
where $ \kappa = w_\rho/w_z$ and
\begin{eqnarray}
f(\kappa)=\frac{1+2\kappa^2}{1-\kappa^2} -\frac{3\kappa^2\mbox{tanh}^{-1}
\sqrt{1-\kappa^2}}{(1-\kappa^2)^{{3/2}}}. \label{eqn:fkappa}
\end{eqnarray}
The  variational equations are obtained by minimizing the energy
(\ref{energy}) by $\partial E/\partial w_z= \partial E/\partial
w_\rho=0$ \cite{4,var}:
\begin{eqnarray}\label{var1}
 {w_{\rho}}\nu^2  &&=
\frac{1}{4w_{\rho}^3} -\frac{
a_{dd}}{\sqrt{2\pi}} \frac{Ng(\kappa)}{w_{\rho}^3w_{z}}
+\sqrt{\frac{3}{5}}\frac{6(6 \pi^2 N)^{2/3}}
{25\pi w_{z}^{2/3}w_{\rho}^{7/3}}
, \\  {w}_{z} \lambda ^2 &&=
\frac{1}{4w_{z}^3}- \frac{ a_{dd}}{\sqrt{2\pi}}
\frac{2Nh(\kappa) }{w_{\rho}^2w_{z}^2}
+\sqrt{\frac{3}{5}}\frac{6(6 \pi^2 N)^{2/3}}{25\pi w_{z}^{5/3}
w_{\rho}^{4/3}},\label{var2}
\end{eqnarray}
where
\begin{eqnarray}
&& g(\kappa) = \frac{2-7\kappa^2-4\kappa^4}{(1-\kappa^2)^2} +
\frac{9\kappa^4\mbox{tanh}^{-1}\sqrt{1-\kappa^2}}{(1-\kappa^2)^{5/2}}, \\
&& h(\kappa) = \frac{1+10\kappa^2-2\kappa^4}{(1-\kappa^2)^2} -
\frac{9\kappa^2\mbox{tanh}^{-1}\sqrt{1-\kappa^2}}{(1-\kappa^2)^{5/2}}.
\end{eqnarray}
The chemical potential $\mu_0$ of the system per particle is
\begin{eqnarray} \label{muvar}
\mu_0 &&= \frac{1}{8}\left[\frac{1}{w_{\rho}^2}+
\frac{1}{2w_{z}^2}
\right]
+\left[\frac{\nu^2 w_{\rho}^2}{2}+\frac{\lambda^2 w_z^2}{4}\right]
\nonumber \\
&&-\frac{2Na_{dd}f(\kappa)   }{\sqrt{2\pi} w_{\rho}^2w_{ z}}
+ \sqrt{\frac{5}{3}}  \frac{9 (6\pi^2  N)^{2/3} }{50
\pi w_{\rho}^{4/3} w_{z}^{2/3}}.
\end{eqnarray}

\subsection{Approximate density for cigar and disk shapes}
\label{IIC}

In many situations of experimental interest, the   Fermi gas
could be subject to
a strong trap either in the polarization $z$ direction or in the
transverse radial $\rho$ plane.
The system then has a one-dimensional (1D) cigar  or 2D disk shape, respectively. In such cases simplified equations in lower dimensions could be very useful  \cite{luca}.

First we consider the reduced 1D equation for a cigar shape.
We assume that there is a strong trap
in the $x-y$ plane
and that the
density in this  plane is given by the Gaussian ground state
\cite{luca}
$n(\rho)=\exp(-\rho^2/d_\rho^2)/(\sqrt \pi d_\rho),  d_\rho = \sqrt {1/(2\nu)}$
of the trap  $\nu^2\rho^2/2$,
so that the 3D density $n({\bf r})$ satisfies
\begin{eqnarray}
\sqrt {n({\bf r})}\equiv  \sqrt{n(\rho)}\sqrt{n(z)}= \frac{ \sqrt{n(z)}
}{\sqrt \pi d_\rho}\exp\left[-
\frac{\rho^2}{2d_{\rho}^2}
\right] .
\end{eqnarray}
Substituting this density in   (\ref{3d1}),
and multiplying this equation by the corresponding Gaussian $\sqrt{n(\rho)}$
and integrating out the $\rho$ dependence
we obtain the reduced 1D equation \cite{1D,lp}
\begin{eqnarray}  \label{1D}
&&\mu_{1D} \sqrt{n(z)}
 =  \biggr[ -\frac{\partial_z^2}{8} +\frac{\lambda^2 z^2}{2}
+   \frac{3[6  N n(z)]^{2/3}}{10 d_\rho^{4/3}  }
 \nonumber \\
&&+  \frac{2a_{dd}N}{d_\rho^2}\int_{-\infty}^{\infty}
\frac{dk_z}{2\pi}e^{ik_z z}n(k_z)s_{1D}\left( \frac{k_zd_\rho}{\sqrt 2} \right)  \biggr]\sqrt{n(z)}
,\\
\label{1Dvdd}
&&s_{1D}(\zeta)=\int_0^\infty du \frac{2\zeta^2 -u}{u+\zeta^2}
   e^{-u},
\end{eqnarray}
where
$n(k_z)=\int_{-\infty}^{\infty} e^{-ik_z z}n(z) dz,$
and  $\mu_{1D}$ is the chemical potential.
An approximate variational solution of   (\ref{1D}) is possible with the following ansatz for density \cite{lp}
 $ n(z)={1}/({\sqrt \pi  w_z})\exp[ -{z^2}/{ w_z^2}  ],$
while the width $w_z$ is determined by solving   (\ref{var2}) with $w_\rho =d_\rho$ and $\kappa= d_\rho/w_z$.

Next we consider the reduced 2D equation suitable for a disk shape
 with a strong axial trap. The dipolar Fermi gas is assumed
to be in the ground state \cite{luca} $n(z)=\exp(-z^2/d_z^2)/(\sqrt \pi d_z), d_z=\sqrt{1/(2\lambda)}$,
of the axial trap $\lambda^2z^2/2$
and the 3D density can  be approximated as
\begin{eqnarray}
\sqrt {n({\bf r})}\equiv  \sqrt{n(\rho)}\sqrt{n(z)}= \frac{\sqrt{n(\rho)}
}{ \pi^{1/4} \sqrt{d_z}}\exp\left[-
\frac{z^2}{2d_z^2}
\right] .
\end{eqnarray}
Using this ansatz in   (\ref{3d1}),
and multiplying   by the corresponding Gaussian $\sqrt{n(z)}$
and integrating out the $z$ dependence
we get
 the   effective 2D equation \cite{lp,2D}
\begin{eqnarray}
&&\mu_{2D}\sqrt{n(\rho)}=
\biggr[ -\frac{\nabla^2_\rho}{8}+\frac{1}{2}\nu^2 \rho^2+
\sqrt{ \frac{3}{5}  } \frac{[6N
n(\rho)]^{2/3}\pi}{2d_z^{2/3}} \nonumber
\\
\label{2D}
&&   + \frac{4\pi a_{dd}N}{\sqrt{2\pi}d_z}
\int\frac{d^2k_\rho}{(2\pi)^2}e^{i{\bf k_\rho .
\rho}}
n({\bf k_{\rho}})  h_{2D}\left( \frac{k_\rho d_z}{\sqrt 2}  \right)
\biggr]  \sqrt{n(\rho)},
\end{eqnarray}
where
$n({\bf k_{\rho}})=\int e^{i\bf k_\rho . \rho }n(\rho) d^2\rho,
h_{2D}(\xi)=2-3\sqrt \pi \xi e^{\xi^2} [1-\mbox{erf}(\xi)].$
 An approximate variational solution of   (\ref{2D}) is possible with the following ansatz for density \cite{lp}
$
  n(\rho)={1}/{( \pi  w_\rho^2)}\exp[ -{\rho^2}{ w_\rho^2}  ],
$
while the width $w_\rho$ is determined by solving   (\ref{var1}) with $w_z =d_z$ and $\kappa= w_\rho/d_z$.

\subsection{Sound propagation in a uniform  Fermi gas}
\label{IID}

To find the sound velocity in a uniform Fermi gas, we evaluate the
Bogoliubov spectrum using the linearized hydrodynamic equations. We consider
a {\it stationary} Fermi gas in a box with periodic boundary condition with the trap $V$ fixing just the allowed plane wave solution. Then   (\ref{den}), after the inclusion of the gradient term of   (\ref{grad}),
reduces to
\begin{equation}\label{h3}
-\frac{\hbar^2\nabla^2 \sqrt{n({\bf r},t)}}{8m\sqrt{n({\bf r},t)}}+\mu(n,{\bf r})+
m \frac{\partial \Phi({\bf r},t)}{\partial t}=0.
\end{equation}
Now we allow small perturbation in $n$ and $\Phi$ around their equilibrium
values $n_0$ and $\Phi_0$ by $n({\bf r},t) \approx n_0+\bar n({\bf r},t)$
and $\Phi  ({\bf r},t)\approx \Phi_0+\bar \Phi({\bf r},t)$, then   (\ref{h1}) leads to
\begin{equation}\label{h4}
\frac{\partial \bar n({\bf r},t)}{\partial t}+n_0\nabla^2 \bar \Phi({\bf r},t) =0.
\end{equation}
In   (\ref{h3}), we need to use $\sqrt {n({\bf r},t)}  \approx \sqrt{n_0}+
\bar n({\bf r},t)/(2\sqrt{n_0})$ and $ [n({\bf r},t)]^{2/3}  \approx n_0^{2/3}+ 2\bar n({\bf r},t)/(3{n_0}^{1/3})$, while $\mu(n,{\bf r})\approx
\tilde \mu_0+\bar\mu (n,{\bf r})$ with
\begin{equation}\label{h5}
\bar \mu(n,{\bf r})=\frac{\hbar^2(6\pi^2)^{2/3}\bar n}{3mn_0^{1/3}}+\int
d{\bf r'}\bar n({\bf r'},t)U_{dd}({\bf r-r'}),
\end{equation}
where $\tilde  \mu_0$ is the stationary value of $\mu$. Then   (\ref{h3}) leads to
\begin{equation}\label{h6}
-\frac{\hbar^2\nabla^2 {\bar n({\bf r},t)}}{16mn_0}+\bar \mu(n,{\bf r})+
m \frac{\partial\bar  \Phi({\bf r},t)}{\partial t}=0.
\end{equation}
Assuming the perturbations $\bar n$ and $\bar \Phi$ in plane-wave forms
$\bar n =\bar n_0\exp[i({\bf k \cdot r})-\omega t]$
and $\bar \Phi =\bar \Phi_0\exp[i({\bf k \cdot r})-\omega t],$ $\bar \mu$ has the form $\bar\mu =\bar n_0 \bar \mu_0\exp[i({\bf k \cdot r})-\omega t]$ with
\begin{equation}\label{h10}
\bar \mu_0= \frac{\hbar^2(6\pi^2)^{2/3}}{3mn_0^{1/3}}+   \frac{4\pi  a_{dd}\hbar^2}{m} (3\cos^2\theta -1),
\end{equation}
where $\theta$ is the  angle  between  the propagation direction and  polarization direction
$z$, and the last term in   (\ref{h10}) is just the Fourier transform of the dipolar potential. Then
   (\ref{h4}) and (\ref{h6}) become:
\begin{eqnarray}
&&-i\omega \bar n_0  -n_0k^2\bar \Phi_0 =0,\\
&&\left[ \frac{\hbar^2k^2}{16mn_0} + \bar\mu_0  \right]\bar n_0 -i\omega m \bar\Phi_0 =0.
\end{eqnarray}
The condition of existence of the solution to this set of equations leads to
the Bogoliubov spectrum
\begin{equation}
\epsilon_k\equiv \hbar \omega = \sqrt{\frac{\hbar^2k^2}{4m}\left[
\frac{\hbar^2k^2}{4m}+4n_0\bar \mu_0  \right]
}.
\end{equation}
The sound velocity, defined as $c_s(\theta)=\lim_{k\to 0}(\epsilon_k/\hbar k)$
for a uniform dipolar Fermi gas can be written as
\begin{equation}\label{vel}
c_s(\theta) = \sqrt{\frac{v_F^2}{3}+ \frac{4\pi n_0 a_{dd}\hbar^2}{m^2}(3\cos^2\theta-1)},
\end{equation}
where  $v_F\equiv \hbar k_F/m =\hbar(6 \pi ^2n_0)^{1/3}/m$ is the
Fermi velocity.
 For a nondipolar Fermi gas ($a_{dd}=0$),   (\ref{vel}) leads to the well-known
 velocity of $v_F/\sqrt 3$ \cite{frmp,lca}.

{The angle-dependent second term under the square root in  (\ref{vel}) is responsible
for anisotropic sound velocity. Specifically, for $\theta > 54.73^\circ$ degrees, this term is negative
corresponding to a decrease in sound velocity. For large dipolar interaction, for $\theta$ greater than
a critical value and for large density $n_0$
the sound velocity could be imaginary corresponding to no propagation. However, in this study we shall only consider moderate values of density and dipolar interaction, that allow anisotropic sound propagation in all directions.
}
\section{Numerical Results}

\label{III}
\subsection{Stationary properties of trapped  dipolar Fermi gas}

\label{IIIA}

For a trapped 3D Fermi system we solve   (\ref{3d1}) numerically after discretization \cite{CPC}. The divergence of the dipolar term at short distances
has been handled by treating this term in momentum ($\bf k$) space \cite{12}.

For numerical calculation in section \ref{IIIA}, we consider a degenerate Fermi gas of  $^{161}$Dy atoms with $a_{dd} \approx 130a_0$ and employ
the oscillator length $l_0 = 1$ $\mu$m.

\begin{figure
\begin{center}
\includegraphics[width=.7\linewidth]{fig1a.jpg}
\includegraphics[width=.7\linewidth]{fig1b.jpg}
\end{center}

\caption{  (a)  Radial density of a disk-shaped
Fermi gas of $^{161}$Dy atoms
in $x-y$ plane
$n(x,y=0)=\int_{-\infty}^{\infty} dz n(x,y=0,z)$ from the 3D
(\ref{3d1}) for  $\nu=1, \lambda=10$
compared
with its variational (var) approximation and
the solution of the reduced 2D   (\ref{2D}) for
$N=100, 500, 10000$.
(b) Axial density of a cigar-shaped
Fermi gas of $^{161}$Dy atoms
along polarization $z$ axis
$n(z)=\int_{-\infty}^{\infty}   \int_{-\infty}^{\infty} dx dy n(x,y,z)$ from the 3D
(\ref{3d1}) for  $\nu=1, \lambda=0.1$
compared
with its variational (var) approximation and
the solution of the reduced 1D   (\ref{1D}) for
$N=100, 1000$.
}

\label{fig1}
\end{figure}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.7\linewidth]{fig2a.jpg}
\includegraphics[width=.7\linewidth]{fig2b.jpg}
\includegraphics[width=.7\linewidth]{fig2c.jpg}
\end{center}

\caption{  (a)  Reduced chemical potential $\mu/(\hbar \omega)$,
and rms sizes $\langle x \rangle $, and $\langle z \rangle $ of a cigar-shaped Fermi
  gas of $^{161}$Dy atoms with
$\nu=1,\lambda =0.3$ versus   $N$ from a numerical (n) solution of the
3D   (\ref{3d1}) and variational  (v) approximation [   (\ref{var1}), (\ref{var2}) and (\ref{muvar})]. (b) The same for a  disk-shaped gas with
$\nu =1, \lambda =10$. (c) The same for the trapped  Fermi gas of 1000 $^{161}$Dy atoms
 versus the   parameter  $\lambda$ with $\nu=1$  from a numerical solution of
    (\ref{3d1}) and variational approximation.
}

\label{fig2}
\end{figure}

The anisotropic dipolar interaction is partially attractive in certain angles and repulsive in others and contributes very little  in a spherically
symmetric trap. The dipolar interaction contributes attractively in the cigar-shaped configuration along the polarization $z$ direction and repulsively in the disk-shaped configuration perpendicular to the polarization $z$ direction.
Hence we will mostly consider the degenerate dipolar Fermi gas in cigar and
disk shapes.
The 3D model (\ref{3d1}) is very strongly nonlinear with nonlinearity ${\cal N}=(6\pi^2 N)^{2/3}/2$, leading to large nonlinearities of ${\cal N}=760$ and 3526 for
$N=1000$ and 10000, respectively.

First, we compare the results of the reduced 2D
density $n(x,y)$ of a disk-shaped degenerate Fermi gas of
$^{161}$Dy atoms  with $\lambda=10, \nu=1$  as obtained from the 3D and 2D models (\ref{3d1}) and (\ref{2D}), respectively, and from the variational
approximation to the 3D model. The reduced 2D density $n(x,y)$ is defined by $n(\rho)\equiv n(x,y)=\int_{-\infty}^{\infty} dz n(x,y,z)$. In  figure \ref{fig1} (a) we plot the reduced 2D density along $x$ axis $n(x,y=0)$ for different  $N$. Next, we compare the results of reduced 1D
density $n(z)$ of a cigar-shaped degenerate Fermi gas of
$^{161}$Dy atoms  with $\lambda=0.1, \nu=1$  as obtained from the 3D and 1D models (\ref{3d1}) and (\ref{1D}), respectively, and from the variational
approximation to the 3D model. The reduced 1D density $n(z)$ is defined by $n(z)=\int_{-\infty}^{\infty} dx \int_{-\infty}^{\infty}dy  n(x,y,z)$. In  figure \ref{fig1} (b) we plot the reduced 1D density along $z$ axis $n(z)$ for different  $N$. From  figure \ref{fig1} we find that in both cigar and disk shapes the 1D and 2D models perform fairly well, even for large nonlinearities, when compared with the results of the full 3D model.
To see if the condition (\ref{con}) is satisfied by the shapes in  figure \ref{fig1}, we can use the TF estimate for $k_F$ of a trapped degenerate Fermi gas
\cite{frmp}: $k_F \approx (48N)^{1/6}/l_0$. For $100<N<10000$, as in  figure \ref{fig1},  the de Broglie wave length at the Fermi surface is
$2\pi/k_F \approx l_0$ . The shapes and sizes in  figure \ref{fig1} are much
larger than this value, consistent with the condition (\ref{con}).

Next we consider a
cigar-shaped degenerate Fermi system of $^{161}$Dy atoms with trap
parameters $\nu=1, \lambda =0.3$. In  figure \ref{fig2} (a) we plot the
numerical and variational results for the chemical potential $\mu$ and
rms sizes $\langle x \rangle $ and $\langle z \rangle $ of this system
for  $10< N <10000$.  In  figure \ref{fig2} (b)
we plot the same for a disk-shaped
degenerate Fermi system with trap
parameters $\nu=1, \lambda =10$. Finally, in  figure \ref{fig2} (c) we plot
the same quantities for $N=1000$ versus the trap anisotropy $\lambda$
for $\nu=1$. In all cases
the variational results are in good
agreement with the  3D model. {  For medium to
small number ($N<10000$)
 of trapped $^{161}$Dy
atoms as considered here and also of experimental interest,
the effect of the dipolar term in Eq. (\ref{3d1}) (the last term in this equation)
is small compared to the Fermi energy term (the last but one term there). Hence
the effect of the dipolar term in figures \ref{fig2} is small. The plots in these
figures only change by less than about two to four percent if we set $a_{dd}=0$.
}

\subsection{Anisotropic sound propagation in uniform dipolar Fermi gas}

The hydrodynamic analytical result for sound velocity (\ref{vel}) in a
uniform dipolar Fermi gas shows a clear anisotropy through the angle
$\theta$ between the propagation direction $\bf r$ and the polarization
direction $z$. In this expression for sound velocity there are two
competing terms under the square root. The first term $v_F^2 /3$
involving the Fermi velocity is isotropic and proportional to
$n_0^{2/3}$ whereas the second term proportional to dipolar interaction
is anisotropic and proportional to density $n_0$. The anisotropy in
sound velocity $c_s(\theta)$ will manifest strongly for large {
strength
$a_{dd}$ and for large }
$n_0$ as
the anisotropic dipolar term in   (\ref{vel}) will increase more
rapidly with $n_0$ than the isotropic term. In  figure \ref{fig3} (a) we
plot the sound velocity $c_s(\theta)$ for angles $\theta=0$ and $\pi$ as
a function of density $n_0$ {for a degenerate Fermi gas with $a_{dd}=130a_0$
($^{161}$Dy atom)  and $2000a_0$ (polar $^{40}$K-$^{87}$Rb molecules in the singlet
rovibrational ground state \cite{pfau})
 using the analytical result (\ref{vel}).}
The sound velocity for
nondipolar atoms ($a_{dd}=0$) is also plotted.  From this plot we find
that the anisotropy in sound propagation as measured by the difference
$[c_s(0)-c_s(\pi/2)]$ is sizable only for a relatively large density {$n_0=10^{15}$ cm$^{-3}$ for $^{161}$Dy atoms. However, for polar $^{40}$K-$^{87}$Rb molecules
the anisotropy is appreciable for a relatively low density of $n_0=10^{13}$ cm$^{-3}$.
From figure \ref{fig3} (a) we see that for $^{40}$K-$^{87}$Rb molecules, the sound velocity is imaginary
for $\theta=0$ for a density of about $n_0=5\times10^{13}$ cm$^{-3}$.}

\begin{figure}[!b]
\begin{center}
\includegraphics[width=.7\linewidth]{fig3a.jpg}
\includegraphics[width=.7\linewidth]{fig3b.jpg}
\includegraphics[width=.7\linewidth]{fig3c.jpg}
\end{center}

\caption{ (a)  The analytical sound velocity $c_s(\theta)$
given by   (\ref{vel}) in axial $z$ ($\theta=0$) and radial $\rho$ ($\theta=\pi/2$) directions
versus density $n_0$ for $a_{dd}=130 a_0$ and 2000$a_0$
as well as the velocity for $a_{dd}=0$.
 (b) The analytical (ana) sound velocity $c_s(\theta)$ in axial $z$ ($\theta=0$) and radial $\rho$ ($\theta=\pi/2$) directions versus $a_{dd}$
compared with the results of numerical simulation (num).
(c) The analytical sound velocity $c_s(\theta)$ versus polar angle $\theta$
compared with the results of numerical simulation.
}

\label{fig3}
\end{figure}

To study sound propagation,  the present stationary (static) 3D model (\ref{3d1}) is generalized to include time variation by replacing $\mu$ by the usual time derivative $i\hbar \partial/\partial t$.   Consequently, the infinite dipolar Fermi gas satisfies the following Galilei-invariant equation
\cite{as}
\begin{eqnarray}  \label{3d2}
&& i\hbar \frac{\partial}{\partial t} \sqrt{n_0{(\bf r)}}
 =  \biggr[ -\frac{\hbar^2 \nabla^2 }{2(2m)} +2V({\bf r})
+2
\frac{\hbar^2}{2m}[6\pi^2  n_0({\bf r})]^{2/3}\nonumber \\
&&+ 12a_{dd}\frac{\hbar^2}{2m}\int \frac{1-3\cos^2 	\theta}{R^3}n_0({\bf r'})d{\bf r'}
\biggr] \sqrt{n_0({\bf r})},
\end{eqnarray}
where density  $n_0({\bf r})=Nn({\bf r})$ is not normalizable for infinite hydrodynamics.
Equation (\ref{3d2}) is consistent \cite{as}
with the time-dependent hydrodynamic
equation (\ref{den})  after  the inclusion
of the gradient correction term  \cite{grad}.

\begin{figure}[!b]
\begin{center}
\includegraphics[width=\linewidth]{fig4a.jpg}
\includegraphics[width=\linewidth]{fig4b.jpg}
\end{center}

\caption{ The contour plot of 3D density in $x-z$ plane
$n_0(x,y=0,z)$ during sound propagation on a background density
$n_0= 10^3$ $\mu$m$^{-3}$ of a degenerate dipolar $^{161}$Dy gas with $a_{dd}
=130a_0$ at times
$t=$
(a) $0.15t_0$, (b) $0.2t_0$, (c) $0.25t_0$, and (d) $0.3t_0$ ($t_0=0.005$ s).
The contour plot of 3D density in $x-y$ plane
   $n_0(x,y,z=0)$  during the same sound propagation at times
(e) $0.15t_0$, (f) $0.2t_0$, (g) $0.25t_0$, and (h) $0.3t_0$.
}

\label{fig4}
\end{figure}

{ The anisotropy of the dipolar interaction would be prominent at low to medium density for large
dipolar interaction.  To illustrate the anisotropic sound propagation we will consider two examples:
$^{161}$Dy atoms at a medium density of 10$^{15}$ cm$^{-3}$ and the polar molecules
$^{40}$K-$^{87}$Rb at the low density of 10$^{13}$ cm$^{-3}$.
First we consider
the numerical simulation of sound propagation in the infinite   $^{161}$Dy gas
($a_{dd}=130a_0$)
at a background density of
$n_0=10^3 $ $\mu$m$^{-3}=10^{15}$ cm$^{-3}$.  }
The numerical simulation is initiated with
a 3D Gaussian pulse
at the center of the uniform 3D background
density given by $n_0({\bf r})
 = (10^3 + 10^2e^{-r^2/2w^2})$ $\mu$m$^{-3}$, $w=2 $ $\mu$m,
subject to a weak expulsive Gaussian potential
$V ({\bf r}) = 0.00001e^{-r^2/2w^2}$ $\mu$m$^{-2}$.
With this initial condition   (\ref{3d2}) is solved by
 real-time propagation \cite{CPC} to study sound waves.
An ellipsoid-like sound wave front is found to emerge outwards upon time propagation. From a numerical study of this wave front the sound velocity
in different directions is calculated.

Typical anisotropic sound propagation in a $^{161}$Dy degenerate gas
of density 10$^{15}$ cm$^{-3}$ is shown in  figure \ref{fig4}.
The sound propagation is illustrated via contour plots of 3D density $n_0(x,y,0)$
in the $x-y$ plane
 and of density $n_0(x,0,z)$ in the $x-z$ plane.
In  figures \ref{fig4} (a), (b), (c), and (d), to illustrate
the anisotropic sound propagation in the $x-z$ plane,
we show the contour plot of  $n_0(x,0,z)$ at times $t= 0.15t_0, 0.2t_0,
0.25t_0$ and $0.3t_0$, respectively,  with
$t_0=2ml_0^2/\hbar^2\approx 0.005$ s the time scale and
$l_0=1$ $\mu$m the
length scale used in the numerical solution of
  (\ref{3d2}). However, the propagation in the radial $x-y$ plane
is isotropic. This isotropic  sound propagation is shown in  figures
\ref{fig4} (e), (f), (g), and (h)
via the contour plot of  $n_0(x,y,0)$
at times $t= 0.15t_0, 0.2t_0,0.25t_0$ and $0.3t_0$, respectively.
A clean wave front of high density is identified in these contour
plots encompassing a region of low density $-$ a typical panorama in sound propagation.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\linewidth]{fig5a.jpg}
\includegraphics[width=\linewidth]{fig5b.jpg}
\end{center}

{
\caption{  The contour plot of 3D density in $x-z$ plane
$n_0(x,y=0,z)$ during sound propagation on a background density
$n_0= 10$ $\mu$m$^{-3}$ of a degenerate dipolar molecular $^{40}$K-$^{87}$Rb gas with $a_{dd}
=1000a_0$ at times
$t=$
(a) $0.2t_0$, (b) $0.5t_0$, (c) $0.8t_0$, and (d) $1.1t_0$ ($t_0=0.004$ s, $l_0=1$ $\mu$m).
The contour plot of 3D density in $x-y$ plane
   $n_0(x,y,z=0)$  during the same sound propagation at times
(e) $0.2t_0$, (f) $0.5t_0$, (g) $0.8t_0$, and (h) $1.1t_0$. }}

\label{fig5}
\end{figure}

The sound propagation
along $z$ direction  (polar angle $\theta=0$)
has a velocity ($c_s=v_F/\sqrt 3$) larger than that
for a nondipolar system ($a_{dd}=0$)
of same density, whereas that in the $x-y$ plane
(polar angle $\theta =\pi/2$)
has a velocity smaller than that for a nondipolar
system as can be seen from   (\ref{vel}).
The analytical  sound velocity $c_s$
for nondipolar atoms  of density $10^{15}$ cm$^{-3}$ and atomic mass 161
is 8.88 mm/s compared to the numerically obtained  velocity of 8.5 mm/s.
For $^{161}$Dy atoms of density $n_0=10^{15}$ cm$^{-3}$
with $a_{dd} = 130a_0$,
the analytical radial sound velocity
is $c_s(\pi/2) = 8.1$ mm/s (numerical 7.0 mm/s), and
the analytical axial  sound velocity is $c_s(0) = 10.3$ mm/s
 (numerical 10.0 mm/s). In  figure \ref{fig3} (b), we present the variation
of axial ($\theta =0$) and radial ($\theta=\pi/2$) sound velocities
versus $a_{dd}$ as obtained from numerical simulation and analytical
consideration (\ref{vel}). In  figure \ref{fig3} (c), the analytical result
for velocity versus the polar angle $\theta$ is compared with the
numerical result. The agreement between
the analytical result
(\ref{vel}) for sound velocity and the result of numerical simulation is
satisfactory considering the very large nonlinearities present in the
system due to large density ($n_0$) and large dipolar interaction.

{
To have a larger anisotropy in sound propagation than in a gas of $^{161}$Dy atoms at a lower density
a stronger dipolar interaction, as in polar molecules, is needed. A convenient polar fermionic molecule
molecule  $^{40}$K-$^{87}$Rb has
a permanent electric dipole moment 0.052 Debye ($a_{dd}\approx 20a_0$) for the triplet rovibrational ground state and 0.566 Debye ($a_{dd}\approx 2000a_0$) for the singlet rovibrational ground state \cite{pfau}.  For illustration, we consider a uniform gas of $^{40}$K-$^{87}$Rb molecules of density $n_0=10^{13}$  cm$^{-3}$ with   $a_{dd} = 1000a_0$.
The numerical simulation is initiated with
a 3D Gaussian pulse
at the center of the uniform 3D background
density given by $n_0({\bf r})
 = (10 + e^{-r^2/2w^2})$ $\mu$m$^{-3}$, $w=2 $ $\mu$m,
subject to a weak expulsive Gaussian potential
$V ({\bf r}) = 0.00001e^{-r^2/2w^2}$ $\mu$m$^{-2}$.
With this initial condition   (\ref{3d2}) is solved by
 real-time propagation \cite{CPC} to study sound waves.
Typical anisotropic sound propagation in a $^{40}$K-$^{87}$Rb degenerate gas
of density 10$^{13}$ cm$^{-3}$ is shown in  figure 5.
Because of the stronger dipolar interaction, a larger anisotropy has appeared
at a lower density in
figure 5 compared to figure \ref{fig4}. For $^{40}$K-$^{87}$Rb molecules
of density $n_0=10^{13}$ cm$^{-3}$
with $a_{dd} = 1000a_0$,
the analytical radial sound velocity
is $c_s(\pi/2) = 2.05$ mm/s (numerical 2.4 mm/s), and
the analytical axial  sound velocity is $c_s(0) = 3.03$ mm/s
 (numerical 3.4 mm/s).
 }

\label{IIIB}

\section{Summary}

We developed a 3D theoretical formulation for describing certain macroscopic observables of a degenerate dipolar Fermi gas appropriate for studying stationary properties, such as, shape, size, energy, chemical potential etc. of a trapped system. The effect of dipolar interaction is negligible in the spherically-symmetric configuration and
the dipolar interaction manifests strongly in the asymmetric cigar and disk shapes. Simple reduced equations in 1D and 2D suitable for studying the trapped degenerate dipolar Fermi gas in cigar and disk shapes, respectively,
are derived. Also, a Gaussian variational approximation for studying these macroscopic properties is derived.  We apply the present formalism to study the stationary properties of a trapped degenerate Fermi gas of $^{161}$Dy atoms. The stationary properties of the 3D model under appropriate conditions are found to be in satisfactory agreement with those from the reduced 1D and 2D models as well as with variational approximation.

Using the present 3D model we also obtained analytical results for
anisotropic sound velocity
as a consequence of the anisotropic dipolar interaction
in a dipolar Fermi gas in agreement with
numerical simulation of the 3D   (\ref{3d2}).
The sound velocity is larger along the polarization direction than
in the transverse plane.

\ack

We thank FAPESP  and  CNPq (Brazil)  for  support.
We thank
Luca Salasnich for useful comments.

\section*{Reference}


\section{Introduction}

Early-Type Galaxies (ETGs) are thought to be a homogeneous class of E and S0 galaxies (\cite[Baldry \textit{et al.} 2004]{Baldry2004}). They are predominantly old and inert, have a high central surface brightness, and cover a wide range of luminosities (\cite[Driver \textit{et al.} 2006]{driver2006}). ETGs are classed as smooth, highly concentrated, and spheroidal systems with a lack of spiral arms. Their luminosity profiles tend to follow the de Vaucouleurs I(r)$\sim$r$^{1/4}$ law  or a more general Sersic law distribution (\cite[D'Onofrio \textit{et al.} 2011]{donofrio2011}).

ETGs are thought to be quiescent at zero redshift, leading to the assumption that they are devoid of both cold gas and dust (\cite[Bregman, Hogg \& Roberts 1992]{bregman1992}). Recent detections of ETGs in infrared and sub-mm regimes have revealed largely unexpected amount of cold gas and dust (eg. \cite[Leeuw \textit{et al.} (2004)]{leeuw2004}). Our interests lie in exploring this dusty presence, with the aims of gaining further understanding about the processes which form ETGs.

\begin{figure}
\begin{center}
 \includegraphics[width=5.0in]{proceedingsplot.eps}
 \caption{Optical Colour-Magnitude Diagram for our samples. The highly concentrated (C$\ge$2.86) optical GAMA sources are highlighted by green plus signs above a background of all the GAMA sources (black dots). On top of this is the sub-mm selected sources from the combined GAMA/H-ATLAS catalogues, shown as dark blue open triangles. The Red Sequence (RS) line (\cite[Bernardi \textit{et al.} 2010]{Bernardietal10}) is shown to indicate the regime above which we expect RS galaxies to sit. This plot shows that our sample of ETGs could not have been selected with a direct colour cut, as both blue and red ETGs are showing up. This is confirmed through visual classifications of the samples.}
   \label{fig1}
\end{center}
\end{figure}

\begin{figure}
 \hspace*{-1cm}
 \includegraphics[width=0.6\textwidth]{multiplot.eps}
 \hspace*{-1cm}
\includegraphics[width=0.6\textwidth]{multiplot2.eps}
\caption{[Left] The distribution of isothermal dust temperatures for the ETG sample with concentration index $\ge$ 2.86. Our mean temperature is the long dashed line and the short dashed line is the mean temperature from  \cite[Amblard \textit{et al.} (2010)]{Amblardetal10}  (for the H-ATLAS Science Demonstration Phase (SDP) sample). [Right] The distribution of dust masses fit by modified Planck functions, emphasising our mean dust mass (long dashed line) and that of the SDP ETG sample of \cite[Rowlands \textit{et al.} (2011)]{Rowlandsetal11} (short dashed line).}
   \label{fig2}

\end{figure}

\section{Selection Techniques}

With the aim of creating a highly robust sample of sub-mm selected ETGs, we experimented with multiple proxies for morphology. Proxies tested were optical colour (eg. \cite[Kaviraj \textit{et al.} 2010]{kaviraj2010}), r band Sersic index (eg. \cite[Leeuw \textit{et al.} 2008]{leeuw2008}) and r band Concentration index (eg. \cite[Blanton \etal\ 2003]{blanton2010)}. We used Sersic index values from the GAMA single Sersic fits to the objects, using GALFIT3 (\cite[Kelvin \textit{et al.} submitted]{Kelvininprep}). The Concentration values are the ratio of SDSS Petrosian 90 radius to Petrosian 50 radius for each galaxy. We combined the GAMA I multi-wavelength survey with Herschel-ATLAS Phase 1 far infrared and sub-mm data to create a catalogue of sub-mm detected ETGs.

We examined the behaviour of the three proxies for both our matched multi-wavelength sample and an optical morphologically classified SDSS sample from \cite[Nair \& Abraham (2010)]{NairAbe2010}. From this we concluded that a straight colour cut is not a viable approach to take, as ETGs tend to be red, but not all red galaxies tend to be ETGs. For example, many red spirals were present in such a cut. Also, as our sub-mm approach picks out dusty galaxies, we may expect some blue ETGs to be present. A Sersic index cut at $\ge$3.5 seemed reasonable, but a concentration index cut of $\ge$2.86 gave us a sample which matched best with the optical visually classified ETGs in the Nair $\&$ Abraham control sample.

We have applied our concentration index criterion, with additional criteria for redshift, ellipticity and apparent effective radius, to our GAMA/H-ATLAS data. With the addition of visual classifications to root out any remaining late-type galaxies, this resulted in a condensed sample of 239 lenticular and elliptical galaxies, shown as open triangles in Fig.\,\ref{fig1}. These sub-mm detected ETGs mostly have intermediate absolute magnitudes (-19.5$\sim$M$_{g}\sim$-21.5).

\section{Sub-mm Single Temperature Fits}
Modified Planck functions were fit to the sub-mm regime of our sample's observed data. The resultant cold dust temperatures and masses are displayed in Fig.\,\ref{fig2}. We find a mean dust temperature of 22.5$^{+10}_{-6}$K (for $\beta=$2.0) and a mean dust mass of 2.038$\times$10$^{8}$M$_{\odot}$. Fig.\,\ref{fig2}  shows that the sub-mm sample contains some ETGs with high dust masses and mainly low dust temperatures compared to these references. The low dust temperatures will need to be checked with additional flux points, including PACS data, covering the peak of the cool emission. Unlike the H-ATLAS galaxies which have been studied up to this point, our sample exhibits only intermediate stellar masses, not high stellar masses.


\section{Introduction}

\textbf{Introduction} During the past decades, the covariant density functional theory (CDFT) has received wide attention due to its successful descriptions of both ground-state and excited-state properties of nuclei all over the nuclear chart.
In this report, we will mainly focus on our recent progress in the localized covariant density functional constrained by the relativistic Hartree-Fock (RHF) theory \cite{Liang2012b}.

Research on quantum mechanical many-body problems is essential in many fields of modern physics.
By reducing the many-body problem formulated in terms of $N$-body wave functions to the one-body level with local density distributions $\rho(\mathbf{r})$, the density functional theory (DFT) of Kohn and Sham \cite{Kohn1965} has accomplished great success.
In particular, this theory states that there always exists a local single-particle potential $V_{\rm KS}(\mathbf{r})$ so that the exact ground-state density of the interacting system can be reproduced by non-interacting particles moving in such local potential.

The DFT has been also widely used in nuclear physics since the 1970s \cite{Bender2003}.
In particular, the Lorentz symmetry taken into account in the CDFT puts stringent restrictions on the number of parameters in the corresponding functionals without reducing the quality of the agreement with experimental data.
During the past decades, the CDFT has been successful in describing a large variety of nuclear phenomena, including the equation of state in nuclear matter, ground-state properties of finite nuclei, collective rotational and vibrational excitations, as well as fission landscapes and low-lying spectra of transitional nuclei involving quantum phase transitions.
The readers are referred to Refs.~\cite{Ring1996,Vretenar2005,Meng2006,Paar2007,Niksic2011} for recent reviews.

Part of these successes could only be achieved because the underlying functionals depend only on local densities, i.e., the relativistic Hartree (RH) scheme.
However, a few common problems have been found in the isovector channels of the widely used RH parametrizations.
Since a few years, the RHF calculations are possible for spherical nuclei all over the nuclear chart \cite{Long2006}.
In contrast to the simple Hartree calculations, they reproduce successfully the proton-neutron effective mass splitting \cite{Long2006} (at the Hartree level it can be reproduced only if a isovector-scalar meson is included \cite{Roca-Maza2011}), shell structure evolutions \cite{Long2007,Long2008,Tarpanov2008,Long2009,Moreno-Torres2010}, and spin-isospin resonances \cite{Liang2008,Liang2009,Liang2012}.
However, since the RHF theory introduces non-local potentials $V_{\rm HF}(\mathbf{r}, \mathbf{r}')$, its theoretical framework is much more involved and the simplicity of the Kohn-Sham scheme is lost.
This restricts the applicability of deformed RHF theory to a few very light nuclei \cite{Ebran2011} and it prevents the treatment of effects beyond mean-field.
Therefore, it is highly desirable to find a covariant density functional based on only local potentials, yet keeping the merits of the exchange terms.

Recently, we started from the important observation that one of the successful and widely used RHF parametrizations PKO2 \cite{Long2008} includes only three relatively heavy mesons, $\sigma$, $\omega$, and $\rho$, but no pions.
Since the masses of these mesons are relatively heavy, the zero-range reduction for the corresponding finite-range Yukawa propagators becomes a reasonable approximation.
As a step further, with the help of Fierz transformation \cite{Okun1982,Sulaksono2003}, the Fock terms can be expressed as local Hartree terms.
In such way, a Hartree-Fock equivalent covariant density functional can be derived while the constraints introduced by
the Fock terms of the RHF scheme can be taken into account.
This RHF equivalent parametrization thus obtained is called PKO2-H \cite{Liang2012b}.
For details, the formalism for the effective Lagrangian density as well as the zero-range reduction and Fierz transformation can be found in Ref.~\cite{Liang2012b}.

In the following, we will take the Gamow-Teller resonance (GTR) in $^{90}$Zr as an example to illustrate the ideas of the localized RHF equivalent RPA and the main conclusions.

\textbf{Results and discussion} In the left panel of Fig.~\ref{Fig1}, the strength distribution of GTR in $^{90}$Zr calculated by the particle-hole (\textit{ph}) residual interactions of PKO2-H is shown as solid line.
It is found that, based on the same unperturbed single-particle spectra calculated by the RHF theory (dotted line), the present results by the localized RHF equivalent RPA calculations are quite similar to those by the original RHF+RPA \cite{Liang2008} calculations (dashed line). The slight difference is due to the zero-range approximation.
The experimental peak energy~\cite{Wakasa1997} is well reproduced.

\begin{figure}
  \includegraphics[width=.50\textwidth]{Figa.eps}
  \includegraphics[width=.50\textwidth]{Figb.eps}
  \caption{(Color online) Left panel: Strength distribution of GTR in $^{90}$Zr calculated with the \textit{ph} residual interactions of PKO2-H \cite{Liang2012b} with Lorentzian smearing parameter 1 MeV. The unperturbed excitations and the original RHF+RPA \cite{Liang2008} results with PKO2 \cite{Long2008} are shown as the dotted and dashed lines for comparison. The experimental peak energy \cite{Wakasa1997} is denoted by the arrow.
  Right panel: Diagonal matrix elements of \textit{ph} residual interactions for the $(\nu1g_{9/2})^{-1}(\pi1g_{7/2})^1$ configuration. The total matrix elements are decomposed into the contributions from different channels.
    \label{Fig1}}
\end{figure}

In order to understand the physical mechanisms, we present in the right panel of Fig.~\ref{Fig1} the diagonal matrix elements of the most important GT \textit{ph} configuration $(\nu 1g_{9/2})^{-1}(\pi 1g_{7/2})^1$.
The total matrix elements (full) are decomposed into the contributions from different channels.

In the original RHF+RPA calculations with PKO2, the \textit{ph} residual interactions are provided by the $\sigma$, $\omega$, and $\rho$ mesons.
Both the direct and exchange terms of the \textit{ph} residual interactions must be taken into account.
It has been demonstrated that the total \textit{ph} strengths are essentially determined by the delicate balance between the $\sigma$ and $\omega$ mesons \cite{Liang2008,Liang2012}.
Since these two mesons are isoscalar, they can contribute via the Fock terms only, and correspondingly the effects of $\sigma$ and $\omega$ mesons are repulsive and attractive, respectively.

By taking the zero-range approximation of the finite-range meson-exchange Yukawa propagators and using the Fierz transformation, the \textit{ph} interactions provided by the $\sigma$, $\omega$, and $\rho$ mesons can be mapped into ten possible channels, i.e., the scalar ($S$), vector ($V$), tensor ($T$), pseudoscalar ($PS$), and pseudovector ($PV$) channels with the isoscalar and isovector ($t$) nature.
For the charge-exchange modes, only the isovector channels contribute, since just the direct terms of the \textit{ph} residual interactions are taken into account when the Fierz transformation has been preformed.

Specifically, for the charge-exchange spin-flip modes like GTR, the dominant residual interactions are represented by the matrix elements of the spin-isospin operator $[\boldsymbol\sigma\vec\tau]\cdot[\boldsymbol\sigma\vec\tau]$ between the large components of Dirac spinors.
They are composed of $[\sigma^{ij}\vec\tau]\cdot[\sigma_{ij}\vec\tau]$ and $[\gamma_5\gamma^{i}\vec\tau]\cdot[\gamma_5\gamma_{i}\vec\tau]$ terms with the coefficients
\begin{equation}
    2\alpha^{\rm H}_{tT} = \frac{1}{8}\frac{g^2_\sigma}{m^2_\sigma}
    \qquad\mbox{and}\qquad
    -\alpha^{\rm H}_{tPV} = \frac{1}{8}\frac{g^2_\sigma}{m^2_\sigma} - \frac{1}{4}\frac{g^2_\omega}{m^2_\omega} + \frac{1}{4}\frac{g^2_\rho}{m^2_\rho},
\end{equation}
respectively, where $\sigma^{ij},\gamma^{i},\gamma_5$ are Dirac matrices, $g_x$ and $m_x$ ($x=\sigma,\omega,\rho$) are the meson-nucleon coupling strengths and meson masses in the original RHF parametrization PKO2 \cite{Long2008}, while $\alpha^{\rm H}$ are the Hartree equivalents.
Correspondingly, the net contribution is then approximately proportional to
\begin{equation}\label{Eq2}
    2\alpha^{\rm H}_{tT}-\alpha^{\rm H}_{tPV}
    \approx \frac{1}{4}\left( \frac{g^2_\sigma}{m^2_\sigma} - \frac{g^2_\omega}{m^2_\omega}\right),
\end{equation}
as $g^2_\rho\ll g^2_\omega, m^2_\rho\approx m^2_\omega$.

As shown in the right panel of Fig.~\ref{Fig1}, it is true that, within the present localized RHF equivalent RPA approach, the total \textit{ph} strengths are mainly composed of the positive $tT$ and negative $tPV$ contributions, while the contributions from the $tS$, $tV$, and $tPS$ channels are negligible.
In other words, the total \textit{ph} strengths are determined by nothing but the delicate balance between the $\sigma$ and $\omega$ mesons (see Eq.~(\ref{Eq2})), which is well established in nuclear covariant density functionals.
Therefore, the constraints introduced by the Fock terms of the RHF scheme into the \textit{ph} residual interactions are straight forward and robust.

\textbf{Summary and perspectives} In summary, a new local RHF equivalent covariant density functional was proposed recently in Ref.~\cite{Liang2012b}, where the constraints introduced by the Fock terms of the RHF scheme have been taken into account.
In this way, the advantages of existing RH functionals can be maintained, while the problems in the isovector channel can be solved.
This opens new and interesting perspectives for the development of nuclear local covariant density functionals with proper isoscalar and isovector properties.
To follow this direction, the investigation on the relativistic local density approximation for the Coulomb Fock terms is in progress \cite{Gu}.

\begin{theacknowledgments}
This work is partially supported by the Major State 973 Program 2007CB815000,
National Natural Science Foundation of China under Grant Nos. 10975008, 11105005, 11105006, and 11175002,
the Research Fund for the Doctoral Program of Higher Education under Grant No. 20110001110087,
the Overseas Distinguished Professor Project from Ministry of Education No. MS2010BJDX001,
the DFG Cluster of Excellence "Origin and Structure of the Universe" (www.universe-cluster.de),
and the JSPS Postdoctoral Fellowship for Foreign Researchers.
\end{theacknowledgments}

\bibliographystyle{aipproc}