\section{Introduction}

Recently  a range of ultra-cold gases with appreciable magnetic \cite{Griesmaier2005a,*Bismut2010a,*Pasquiou2011a,Mingwu2011a,Aikawa2012a,Lu2012a}  and electric \cite{Aikawa2010a,Ni2008a} dipoles have been realized in experiments.  The key feature of these gases is that the constituent particles interact via a dipole-dipole interaction (DDI) that is long-ranged and anisotropic. When these systems are polarized by an external field  the DDI causes their shape to distort. The first experimental evidence for such a magnetostrictive effect\footnote{We adopt terminology appropriate to magnetic dipole gases. We note that for polar molecule gases the effects we predict are electrostrictive.} was in the expansion of   dipolar Bose-Einstein condensates (BECs) \cite{Stuhler2005a}  (also see \cite{Yi2003a,Giovanazzi2003a,Giovanazzi2006a,Mingwu2011a}).

Here we examine magnetostrictive effects in normal Bose and Fermi gases with polarized dipoles. In our work we distinguish between two cases: position space magnetostriction  and momentum space magnetostriction, whereby the position and momentum space distributions, respectively, distort relative to the direction of the polarized dipoles. We show that position space magnetostriction is qualitatively similar for Bose and Fermi gases: the gas tends to elongate along the dipole polarization direction to reduce the direct interaction energy. Momentum space magnetostriction arises solely from exchange interactions and distinguishes the quantum statistics of the particles: the Fermi gas momentum distribution elongates along the polarization direction, while the Bose gas instead contracts along this direction.

Previous studies have considered the distortion of ultra-cold gases with DDIs. Work on Bose gases has focused on position space magnetostriction for zero temperature BECs  \cite{Goral2000a,Goral2002a,Eberlein2005a}. For normal Fermi gases Hartree calculations (excluding exchange) were used to study position space magnetostriction in \cite{Goral2001a}. The first theories including exchange interactions for the trapped Fermi gas were based on simple variational approximations    \cite{Miyakawa2008a,Lima2010a,Lima2010b,Endo2010a} (also see \cite{Endo2011a}), and revealed distortion of the momentum distribution.
In this paper we use full Hartree-Fock (HF) calculations to provide a unified quantitative treatment of trapped Bose and Fermi gases. While HF calculations have been applied to trapped Fermi gases \cite{Zhang2009a,Zhang2010a,Baillie2010b}, this is the first application of this theory the trapped dipolar Bose gas (c.f.~Hartree analysis reported in \cite{Bisset2011}, also see \cite{Ticknor2012a}).    Our main concern is in characterizing  magnetostrictive effects,  their relationship to  direct and exchange interactions, and to pair correlations. We support our numerical results by analytical approximations.

\section{Theory}
We consider a gas of  polarized dipolar particles that interact by a long-range dipole-dipole interaction of the form
\begin{align}
U_{\mathrm{dd}}(\mathbf{x})=\frac{C_{\mathrm{dd}}}{4\pi}\frac{1-3\cos^2\theta}{|\mathbf{x}|^3},\label{e:udd}
\end{align}
where the constant $C_{\mathrm{dd}}$ is given by $\mu_0\mu^2_m$ for magnetic dipoles of strength $\mu_m$ and $d^2/\epsilon_0$ for electric dipoles of strength $d$, and $\theta$ is the angle between $\mathbf{x}$ and the polarization axis, which we take to be the $z$ direction. The atoms are confined within a cylindrically symmetric harmonic trap
\begin{equation}
U_{\mathrm{tr}}(\mathbf{x})=\frac{m}{2}\left[\omega_\rho^2(x^2+y^2)+\omega_z^2z^2\right],
\end{equation}
with aspect ratio $\lambda=\omega_z/\omega_\rho$.

Following \cite{Zhang2009a,Zhang2010a,Baillie2010b,Zhang2011a} we introduce the dimensionless interaction parameter $D_t = C_{\mathrm{dd}} N^{1/6}/(4\pi\hbar\omega a_{\mathrm{ho}}\xspace^3)$ where $\omega=(\omega_\rho^2\omega_z)^{1/3}$ and $a_{\mathrm{ho}}\xspace=\sqrt{\hbar/m\omega}$.  When considering thermal effects we characterize the temperature either in terms of the ideal condensation temperature, $T_c^0=\sqrt[3]{N/\zeta(3)}\hbar\omega/k_B$, with $\zeta(s)$ the zeta function, or the Fermi temperature, $T_F^0=\sqrt[3]{6N}\hbar\omega/k_B$. We note that these two degeneracy temperatures relate as $T_c^0\approx0.52T_F^0$.

\subsection{Semiclassical theory}
In the semi-classical approximation \cite{Goral2001a,Zhang2010a,Baillie2010b} the system is described by the Wigner function
\begin{align}
W(\mathbf{x},\mathbf{k})=\frac{1}{\exp([\epsilon(\mathbf{x},\mathbf{k})-\mu]/k_BT)-\eta},\label{e:wigner}
\end{align}
where $\eta=1$ for bosons and $\eta=-1$ for fermions, and $\mu$ is the chemical potential. This approximation furnishes a good description when the temperature is large compared to the trap level spacing ($k_BT\gg \hbar\omega$). For fermions this approximation can be applied at $T=0$ as long as the Fermi energy is sufficiently large (i.e.~$\mu\gg\hbar\omega$) \cite{Zhang2009a}.
For the case of bosons at $T<T_c$, where $T_c$ is the condensation temperature \cite{Glaum2007}, a condensate emerges and the semi-classical theory requires augmentation by a beyond-semi-classical description of the condensate. Here we restrict our attention to normal Bose and Fermi gases.

The position space and momentum space densities are given by
\begin{align}
n(\mathbf{x})=\int  \frac{d\mathbf{k}}{(2\pi)^{3}} W(\mathbf{x},\mathbf{k}),\label{nx}\\
\tilde{n}(\mathbf{k})= \int \frac{d\mathbf{x}}{(2\pi)^3} W(\mathbf{x},\mathbf{k}),
\end{align}
respectively and the total number is given by $N=\int d\mathbf{x} \,n(\mathbf{x})$.

\subsection{Characterization of magnetostriction}
Our primary concern is distortions in the trapped clouds arising from the DDI. It is convenient to characterize this in terms of the rms-widths
\begin{align}
\sigma_{\nu}&\equiv\left[\frac{1}{N}\int \frac{d\mathbf{x} d\mathbf{k}}{(2\pi)^3}\nu^2W(\mathbf{x},\mathbf{k})\right]^{1/2},
\end{align}
with $\nu=\{x,y,z,k_x,k_y,k_z\}$.

In terms of these quantities we define the momentum and position distortion as
\begin{align}
\alpha&\equiv\frac{\sigma_{k_x}}{\sigma_{k_z}},\quad
\beta\equiv\frac{1}{\lambda}\frac{\sigma_{x}}{\sigma_{z}},
\end{align}
respectively, i.e.~as the ratio of the widths in the directions perpendicular ($x$ or $k_x$) and parallel ($z$ or $k_z$) to the dipolar polarization direction.
 The trap anisotropy appears in the definition of $\beta$ so that in the absence of interactions, where the position density has the  anisotropy imposed by the trap, we have $\beta=1$.  In contrast the non-interacting momentum distribution is spherically symmetric. We refer to systems as exhibiting position space magnetostriction when $\beta\ne1$ and momentum space magnetostriction when $\alpha\ne1$.

\subsection{Hartree-Fock theory\label{s:formalism}}
The HF dispersion relation is \cite{Zhang2010a,Lin2010a,Baillie2010b}
\begin{align}
\epsilon(\mathbf{x},\mathbf{k})=\frac{\hbar^2k^2}{2m}+U_{\mathrm{tr}}(\mathbf{x})+\Phi_D(\mathbf{x})+\eta\Phi_E(\mathbf{x},\mathbf{k}),\label{e:epsilon}
\end{align}
with
\begin{align}
\Phi_D(\mathbf{x})&=\int \frac{d\mathbf{x}' d\mathbf{k}'}{(2\pi)^3}U_{\mathrm{dd}}(\mathbf{x}-\mathbf{x}')W(\mathbf{x}',\mathbf{k}'),\label{e:PhiD}\\
\Phi_E(\mathbf{x},\mathbf{k})&=\int \frac{ d\mathbf{k}'}{(2\pi)^3}\tilde{U}_{\mathrm{dd}}(\mathbf{k}-\mathbf{k}')W(\mathbf{x},\mathbf{k}'),\label{e:PhiE}
\end{align}
the direct and exchange interaction terms, respectively. In Eq.~\eqref{e:PhiE}, $\tilde{U}_{\mathrm{dd}}(\mathbf{k})=C_{\mathrm{dd}} (\cos^2\theta_\mathbf{k}-1/3)$  is the Fourier transform of the dipole-dipole interaction, where $\theta_\mathbf{k}$ is the angle between $\mathbf{k}$ and $\hat{k}_z$.  The HF theory requires solving Eqs.~\eqref{e:wigner}, \eqref{e:epsilon}-\eqref{e:PhiE} in an iterative procedure until the equations are self-consistent. We are interested in solutions for a specified total number, $N$, which requires an additional iterative procedure to adjust $\mu$.

Our system is cylindrically symmetric, since the symmetry axis of the harmonic trap  coincides with the polarization direction. We can exploit this symmetry to simplify the Wigner function to a	function of four	variables $\{\rho,z,k_\rho,k_z\}$, i.e., cylindrical coordinates in $\mathbf{x}$ and $\mathbf{k}$ space.
Our numerical algorithm makes use of Fourier transform techniques and quadrature based on discrete cosine and Bessel functions, and is detailed in \cite{Baillie2010b} (also see \cite{Zhang2010a}). The numerical implementation of HF theory requires extensive computational resources because of the necessity to construct the full Wigner function on grids with suitable resolution to accurately evaluate the DDI.

\subsection{Hartree theory}
We note that the direct interaction potential (\ref{e:PhiD}) can be evaluated in terms of the position space density, i.e.~$\Phi_D(\mathbf{x})=\int  d\mathbf{x}'\,U_{\mathrm{dd}}(\mathbf{x}-\mathbf{x}')n(\mathbf{x}')$, thus if the exchange term is neglected the $\mathbf{k}$ dependence of the dispersion (\ref{e:epsilon}) is  quadratic and Eq.~(\ref{nx}) can be evaluated \cite{Baillie2010b,Bisset2011}
\begin{equation}
n(\mathbf{x})=\frac{1}{\lambda_{\mathrm{dB}}^3}\zeta_{3/2}^\eta\left(e^{[\mu-V_{\mathrm{eff}}(\mathbf{x})]/k_BT}\right),\label{den}
\end{equation}
where  $\zeta_{\alpha}^{\eta}(z)=\sum_{j=1}^{\infty}\eta^{j-1}z^j/j^\alpha$ is the polylogarithm function,  $\lambda_{\mathrm{dB}}=h/\sqrt{2\pi mk_BT}$, and
\begin{equation}
V_{\mathrm{eff}}(\mathbf{x})\equiv U_{\mathrm{tr}}(\mathbf{x})+\Phi_D(\mathbf{x}),\label{Veff}
\end{equation}
is the effective potential. By neglecting exchange we have arrived at a Hartree theory for the system, which involves solving Eqs.~(\ref{den}) and (\ref{Veff}) self-consistently. This theory is computationally much simpler because it does not require the evaluation of the full Wigner function. The Hartree theory provides a baseline for investigating the exchange effects included in the full HF calculations.

\subsection{Analytic approximations}
To develop an analytic approximations for the magnetostrictive behavior we introduce the following ansatz for the HF dispersion relation:
\begin{align}
    \epsilon_A(\mathbf{x},\mathbf{k}) = \frac{\hbar^2}{2m}\left( \kappa_\rho k_\rho^2  + \kappa_z k_z^2\right) +\frac{m\omega_\rho^2}{2}\left( \chi_\rho \rho^2+\chi_z \lambda^2 z^2\right),
\end{align}
where $\kappa_\rho,\kappa_z,\chi_\rho$ and $\chi_z$ are variational parameters to be determined.
\subsubsection{Zero temperature fermions}
We first consider an analytic treatment of $T=0$ fermions. For this system all single particle states up to the Fermi energy are occupied.
We use a Thomas-Fermi ansatz for the Wigner function (after \cite{Miyakawa2008a})
\begin{align}
    W_0(\mathbf{x},\mathbf{k}) &= \Theta[\mu-\epsilon_A(\mathbf{x},\mathbf{k})], \label{e:W0}
\end{align}
where $\Theta[\cdot]$ is the Heaviside step function. In the appendix, we minimize the energy to determine the variational parameters, from which we obtain $\alpha$ and $\beta$. We find, to first order, that
\begin{align}
    \alpha &\approx  1 - c_\alpha D_t, \label{e:alphastep}\\
    \beta &\approx 1 + c_\beta  \lambda^{-2}J'(\lambda^{-2}-1)  D_t, \label{e:betastep}
\end{align}
where $c_\alpha =  \frac{2^{38/3}}{3^{17/6} 175 \pi^2} \approx 0.167$, $c_\beta =  \frac{2^{35/3}}{3^{11/6}35 \pi^2} \approx 1.26$, and the function $J(u)$ is defined in Eq.~(\ref{e:Ju}).
\subsubsection{Finite temperature}
To account for thermal effects in Bose and Fermi gases we use a  Boltzmann ansatz for the Wigner function  (after \cite{Sogo2009a,Endo2011a})
\begin{align}
    W_T(\mathbf{x},\mathbf{k}) &= e^{\left[ \mu -\epsilon_A(\mathbf{x},\mathbf{k})\right]/k_BT}.\label{e:WT}
\end{align}
This approximation fails to account for the quantum statistics (e.g.~phase space densities greater than unity are not prohibited for Fermions), and thus is only applicable at temperatures well above the relevant degeneracy temperatures.\footnote{We are unable to develop simple closed-form analytic results using degenerate forms of the Wigner function [c.f. Eq.~\eqref{e:wigner}]. Nevertheless, such a degenerate form could be treated numerically and would lead to better agreement with the full HF results.}
 In the appendix, we minimize the free energy\footnote{We note that Bose and Fermi statistics are included in the treatment of interactions by adding the appropriate exchange interaction term.} derived from (\ref{e:WT}) to determine the variational parameters and hence the distortion parameters. To first order we find
 \begin{align}
   \alpha &\approx 1 + \frac{\eta}{20\sqrt{\pi}} \left(\frac{\hbar\omega N^{1/3}}{k_BT}\right)^{5/2} D_t   \label{e:alphagaussian},\\
   \beta &\approx 1 +\frac{3}{8\sqrt{\pi}}\lambda^{-2}  J'(\lambda^{-2}-1)  \left(\frac{\hbar\omega N^{1/3}}{k_BT}\right)^{5/2} D_t.  \label{e:betagaussian}
\end{align}
\section{Results}

\subsection{Position space magnetostriction}
 %
In Fig.~\ref{f:beta}(a) we compare the HF position space distortion against the analytic expression (\ref{e:betagaussian}), and Hartree results.
For all results, $\beta$ decreases with increasing dipole interaction strength, i.e.~the gas distorts so as to extend along the direction that the dipoles are polarized.
This distortion increases the number of dipoles in the energetically favorable head-to-tail configuration.
The numerical results terminate at finite DDI due to the system nearing mechanical instability (e.g.~see \cite{Bisset2011}). We observe that at any given value of $D_t$ the position distortion is greatest in near-spherical geometries. The analytic results are independent of whether particles are bosons or fermions, and predict a linear change in $\beta$ with the interaction strength.  The Hartree results are in  good agreement with the full HF solutions for the distortion. This demonstrates that the position space distortion arises from the direct interaction.
\begin{figure}[!tH]
\begin{center}
\includegraphics[width=3.4in]{fig1a.pdf}
\includegraphics[width=3.4in]{fig1b.pdf}
\caption{\label{f:beta}(Color online) Position space distortion for fermions using HF (solid curves) and Hartree (diamonds), and for bosons using HF (dashed curves) and Hartree (circles), also the predictions from   \eqref{e:betastep} (crosses) and \eqref{e:betagaussian} (thin curves) with $\lambda$ = 0.1 [blue (dark gray)], 1 [green (light gray)], and 10 [red (gray)]. (a) Bosons and fermions at $T=1.5T_c^0\approx0.78T_F^0$. (b) Fermions with $D_t=1$. (inset) Fermions at $T=0$.}
\end{center}
\end{figure}

The $T^{-5/2}$ scaling in the finite temperature analytic prediction \eqref{e:betagaussian} suggests that large magnetostriction effects will be obtained at low temperature. For bosons the critical temperature for condensation $T_c$ forms a lower bound for the applicability of our theory.
For fermions our theory could be applied as $T\to0$, however the non-degenerate form of the distribution used to derive the analytic expression \eqref{e:betagaussian} will be unreliable below the Fermi temperature. In Fig.~\ref{f:beta}(b) we compare the finite temperature analytic results to the HF calculations for a range of temperatures below the Fermi temperature. These results indicate that the analytic approximation \eqref{e:betagaussian} significantly overestimates the distortion for $T\lesssim 0.5T_F^0$, and diverges as $T\to0$. We emphasize that this breakdown is not entirely due to the first order expansion used to derive the analytic results, but also because the distribution  \eqref{e:WT} fails to work in the degenerate regime. The $T=0$ analytic expression \eqref{e:betastep} is also shown in Fig.~\ref{f:beta}(b) and provides a qualitative estimate of the HF calculations at $T\to0$. However, for the interaction strength used in Fig.~\ref{f:beta}(b)  ($D_t=1$) corrections to our linearized result are apparent. In the inset to Fig.~\ref{f:beta}(b)  we compare the HF results to the analytic approximation at $T=0$ showing that good agreement is obtained for smaller values of $D_t$.

\begin{figure}[!tH]
\begin{center}
\includegraphics[width=3.4in]{fig2a.pdf}
\includegraphics[width=3.4in]{fig2b.pdf}
\caption{\label{f:alpha}(Color online) Momentum distortion for fermions (solid curves) and bosons (dashed curves) and the predictions from \eqref{e:alphastep} (black cross) and \eqref{e:alphagaussian} (thin black curves) with  $\lambda$ = 0.1 [blue (dark gray)], 1 [green (light gray)], and 10 [red (gray)]. (a) Bosons and fermions at $T=1.5T_c^0\approx0.78T_F^0$, (b) Fermions with $D_t=1$, (inset) Fermions at $T=0$.}
\end{center}
\end{figure}

\subsection{Momentum space magnetostriction}
In Fig.~\ref{f:alpha}(a) we compare the full HF results for the momentum space distortion against the analytic expression (\ref{e:alphagaussian}). For the fermionic system $\alpha$ decreases with increasing dipole interaction strength, i.e.~the gas distorts in momentum space to increase its extent in the  $k_z$ direction. In contrast, for the bosonic system, $\alpha$ increases with increasing dipole interaction strength, thus the system contracts along the $k_z$ direction relative to its radial extent. The Hartree theory (not shown) predicts no momentum distortion (i.e.~$\alpha=1$), thus momentum distortion arises solely from the exchange interactions, which are neglected in the Hartree approximation. As in the previous section, the numerical results terminate at finite DDI due to the system nearing mechanical instability.

In Fig.~\ref{f:alpha}(b)  we consider the Fermi system in the low temperature regime ($T\lesssim T_F^0$). Similar to the position space distortion [see Fig.~\ref{f:beta}(b)] we observe that the analytic approximation \eqref{e:alphagaussian} significantly overestimates the distortion for $T\lesssim 0.5T_F^0$, and diverges as $T\to0$. The $T=0$ analytic expression \eqref{e:alphastep} provides a qualitative estimate of the HF calculations at $T\to0$. In the inset we compare the HF results to the analytic approximation at $T=0$ showing that good agreement is obtained for smaller values of $D_t$.

\subsection{Exchange interactions and second order correlations}\label{SecDirectExchange}
To understand the role of exchange effects it is useful to consider the second order correlation function
\begin{align}
G^{(2)}(\mathbf{x},\mathbf{x}')=\left\langle \hat{\psi}^\dagger(\mathbf{x}) \hat{\psi}^\dagger(\mathbf{x}') \hat{\psi}(\mathbf{x}')\hat{\psi}(\mathbf{x})   \right\rangle,
\end{align}
which describes the probability of finding particles at positions $\mathbf{x}$ and $\mathbf{x}'$, with $\hat{\psi}(\mathbf{x})$ being the field operator. It is convenient to transform this to center of mass $\mathbf{R}=\tfrac{1}{2}(\mathbf{x}+\mathbf{x}')$ and relative $\mathbf{r}=\mathbf{x}-\mathbf{x}'$ coordinates, i.e.~$G^{(2)}(\mathbf{R},\mathbf{r})$, and introduce the  pair correlation function
\begin{equation}
C(\mathbf{R},\mathbf{r})= {G^{(2)}(\mathbf{R},\mathbf{r})} -n(\mathbf{R}+\tfrac{1}{2}\mathbf{r})n(\mathbf{R}-\tfrac{1}{2}\mathbf{r}).
 \end{equation}
Within HF theory, higher order correlation functions decompose into products of first order correlation functions (e.g.~see \cite{Naraschewski1999a,Baillie2010b}) giving
\begin{equation}
C(\mathbf{R},\mathbf{r}) = {\eta} \left|\int \frac{d\mathbf{k}}{(2\pi)^3}W(\mathbf{R},\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{r}}\right|^2.\label{cHF}
\end{equation}
The $\eta$ prefactor of this result shows the important effect of quantum statistics on correlations: For fermions $C(\mathbf{R},\mathbf{r})\le0$ indicates their tendency to anti-bunch, while for bosons $C(\mathbf{R},\mathbf{r})\ge0$ indicates their tendency to  bunch \cite{Jeltes2007a}. Around any given point $\mathbf{R}$ the pair correlations decay as $|\mathbf{r}|$ increases, with this decay length defining the  correlation length.  Equation (\ref{cHF}) shows that the $\mathbf{r}$ dependence of the correlation function is related to the Fourier transform of the momentum distribution. Thus anisotropy in the momentum dependence of $W(\mathbf{R},\mathbf{k})$ will be accompanied by an inverse anisotropy in the  $\mathbf{r}$ dependence of the pair correlation function.

We present a comparison of $C(\mathbf{R},\mathbf{r})$ for a Bose and Fermi gas in Fig.~\ref{f:g2} that demonstrates the features discussed above. The inset to that figure shows the density distribution, and reveals that the correlation length scale ($\sim\lambda_{\mathrm{dB}}$) is much smaller than the spatial extent of the systems.

In terms of the density and pair correlation functions we can evaluate the interaction energy densities
\begin{align}
\mathcal{E}_D(\mathbf{R})&=\frac{1}{2}\int d\mathbf{r}\,U_{\mathrm{dd}}(\mathbf{r})\,n(\mathbf{R}+\tfrac{1}{2}\mathbf{r})n(\mathbf{R}-\tfrac{1}{2}\mathbf{r}),\label{EDden}\\
\mathcal{E}_E(\mathbf{R})&=\frac{1}{2}\int d\mathbf{r}\,U_{\mathrm{dd}}(\mathbf{r})\,C(\mathbf{R},\mathbf{r}),\label{EEden}
\end{align}
for the direct and exchange terms, respectively. Note that the total direct and exchange energies, obtained by integrating the respective energy densities, are given by
\begin{align}
E_D&=\int d\mathbf{R}\,\mathcal{E}_D(\mathbf{R})=\frac{1}{2}\!\int \frac{d\mathbf{x} d\mathbf{k}}{(2\pi)^3} \Phi_D(\mathbf{x})W(\mathbf{x},\mathbf{k}),\\
E_E&=\int d\mathbf{R}\,\mathcal{E}_E(\mathbf{R})=\frac{\eta}{2}\!\int \frac{d\mathbf{x} d\mathbf{k}}{(2\pi)^3} \Phi_E(\mathbf{x},\mathbf{k})W(\mathbf{x},\mathbf{k}).\label{EE}
\end{align}
Equation (\ref{EDden}) shows that the direct energy arises from the density distribution of the system, while the exchange term (\ref{EEden}) arises from correlations, i.e.~with the integrand in Eq.~(\ref{EEden}) only being non-zero over the small range of $\mathbf{r}$ values for which correlations extend   [e.g.~see Fig.~\ref{f:g2}].

The anisotropy of the correlations ensures that the exchange energy is minimized. To see this we note  that the angular part of $U_{\mathrm{dd}}(\mathbf{r})$ is proportional to the $Y_2^0(\theta,\phi)$ spherical harmonic, so that the integral in Eq.~(\ref{EEden}) is sensitive to a quadrupolar moment (i.e.~distortion) of $C(\mathbf{R},\mathbf{r})$ in the $\mathbf{r}$ variable.  For bosons, where the correlation is positive, the correlation peak is prolate:   short range correlations increase the number of pairs in the energetically favorable head-to-tail configuration. In fermions, where the correlation is negative, the correlation dip is oblate:    short range correlations reduce the number of pairs in the energetically unfavorable side-by-side configuration.

Finally, while our focus here has been on second order correlations, because of their relationship to energy, the exchange effects are also apparent in the first order coherence function $G^{(1)}(\mathbf{x},\mathbf{x}')\equiv\langle \hat{\psi}^\dagger(\mathbf{x}')\hat{\psi}(\mathbf{x})\rangle$ which, in center of mass and relative coordinates, relates to the
distribution as
\begin{equation}
G^{(1)}(\mathbf{R},\mathbf{r}) = \int \frac{d\mathbf{k}}{(2\pi)^3}W(\mathbf{R},\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{r}}.\label{G1}
\end{equation}
  This function also exhibits the anisotropy from exchange interactions with a differing sense of distortion for Bose and Fermi gases.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=3.4in]{fig3.pdf}
\caption{\label{f:g2} (Color online) Second order correlation function for $\lambda=10$, $D_t=4$, $N=10^5$, $T=1.1T_c^0\approx 0.57T_F^0$ for Fermi (solid curve) and Bose (dashed curve) for $\mathbf{r}=\rho\hat{\rho}$ (blue, dark grey) $\mathbf{r}=z\hat{z}$ (red, light grey). Inset shows the density profiles for comparison.  }
\end{center}
\end{figure}

\section{Discussion and Conclusions}

In this paper we have examined  magnetostrictive effects in normal Bose and Fermi gases with dipole-dipole interactions.
We have performed  Hartree-Fock calculations for these systems and compared the results against Hartree calculations and analytic approximations that we have developed. Magnetostrictive effects are only apparent for moderately large dipole strengths ($D_t$ of the order of unity), as being realized in current experiments with polar molecules. E.g.~$N=10^5$ KRb molecules (for which  both Bose \cite{Aikawa2010a} and Fermi \cite{Ni2008a} molecules have been realized) with $d=0.5\:\mathrm{D}$ in a  $\omega=2\pi\times100$ s$^{-1}$ trap corresponds to $D_t=4.7$ (higher interaction strengths are possible with more polar molecules, e.g.~NaK which has $d\approx2.72\:\mathrm{D}$ \cite{Wu2012a}). We note that Hartree-Fock theory should provide an accurate description of the system in the weak-coupling regime, e.g.~in the regime of dipolar atom  (magnetic dipole) experiments. For polar molecule experiments, because of the larger dipole strengths, corrections to HF theory may become important, however it is expected that the HF results should be  qualitatively accurate (e.g.~see discussion in Ref.~\cite{Kestner2010a}, also see \cite{Chan2010a,*Liu2011a,*Wang2008a}).  Here we have focused on the case of pure DDIs to allow a more direct comparison of Bose and Fermi systems.  However, we make a few general comments about the effect of contact interactions on the Bose system based on calculations we have performed. First, repulsive contact interactions allow the normal Bose gas to be stable out to higher dipole strengths (e.g.~see \cite{Bisset2011}), where the DDI induced distortions are larger. Second, for a given case (i.e. number of particles, DDI strength and temperature) the addition of repulsive contact interactions tends to lower the density and this in turn causes the DDI induced distortions to slightly decrease.

It should be possible to experimentally verify the effects we have studied in this paper. Spatial distributions have been precisely measured using absorption imaging in cold atom experiments   \cite{Hung2011a,Hung2011b}, and the extension of these techniques to the dipolar gas should allow the position space magnetostrictive behavior to be characterized. The direct measurement of the momentum space effect is possibly more challenging. Releasing the system from the trap requires a model of the expansion dynamics (e.g.~see   \cite{Sogo2009a,Lima2010b,Zhang2011a}) and does not simply reveal the momentum behavior.   Bragg spectroscopy could be used to probe the momentum distribution \textit{in situ} \cite{Stenger1999a,Stamper-Kurn1999a}.  As discussed in this paper  the momentum space magnetostrictive effects in dipolar gases are revealed by the development of anisotropic correlation functions. These could be revealed by a measurement of the  second order (density-density) correlation function, or by measurements of the first order coherence function [c.f.~Eq.~(\ref{G1})]. Such measurements have been performed with high accuracy on alkali Bose gases   \cite{Hung2011a,Hung2011b,Donner2007a}.

\section*{Acknowledgments}  This work was supported by the Marsden Fund of New Zealand (contract UOO0924).


\section{Introduction}
Resonances are ubiquitous in the Physical Sciences and Engineering. In a general sense, a resonance expresses that a system is efficiently coupled to an external excitation within the interval of a certain parameter. When properly controlled, this efficient coupling can be used to measure different phenomena with high sensitivity. For instance, the control of nuclear magnetic resonances in nuclei has enabled a great number of applications in biology and medicine \cite{Gadian1982}. Unfortunately, the lack of control of these resonant phenomena can derive in catastrophic results, one of the most famous being the collapse of the Takoma bridge in the 1940's \cite{Billah1991}. Recently, resonances in plasmonic systems have been gaining a lot of interest \cite{Garcia-Vidal2010} and some techniques to control them have been proposed \cite{Lopez2002, Mojarad2008, preKim2012}. Most resonances can be approximately modelled with the Lorentzian distribution. However, in 1961 a new fundamental kind of resonance was discovered by Fano \cite{Fano1961}. Although it was firstly regarded as a quantum mechanical effect, the Fano resonance is very present in optics and plasmonics. In particular, one of the most important analytical solutions of Maxwell equations in inhomogeneous media, \textit{i.e.} Mie Theory, gives rise to Fano resonances, too \cite{Fano2010}. The Mie Theory \cite{Hergert2012} or the more modern Generalized Lorenz-Mie Theory (GLMT) \cite{GLMT_book} describes the interaction between an incident electromagnetic field propagating in a lossless, homogeneous, isotropic medium and a homogeneous isotropic sphere. The number of applications of this theory is immense. The reader is referred to Gouesbet \textit{et al.} for an extensive review \cite{Gouesbet2011}.

The most appropriate electromagnetic modes for describing the physical phenomena associated with the GLMT are the multipolar modes: $\lbrace \mathbf{A}_{jm_z}^{(m)},  \mathbf{A}_{jm_z}^{(e)} \rbrace$, where $(m)$ and $(e)$ stand for magnetic and electric multipoles. These modes are eigenvectors of the total angular momentum (AM) operator $J^2$ and one of its projections such as $J_z$ with respective values $j$ and $m_z$ \cite{Rose1955, Nora2012}. This set of modes exploit the spherical symmetries of the kind of problems studied by the GLMT. As a consequence, the electromagnetic eigenmodes of the scattered field from a spherical object are precisely the multipolar modes. Experimentally, it is very challenging to prepare and detect a single multipole with an arbitrary value of both $j$ and $m_z$. This is due to the spherical and non-paraxial character of these modes. In this paper we propose a technique to control multipolar resonances of dielectric spheres by engineering the multipolar content of experimentally realizable light beams. In our approach, we use focused cylindrically symmetric beams. These light beams possess a well defined value of $J_z$. In 1992 Allen and co-workers showed that in the paraxial approximation one could find a set of modes which had a well defined value of the AM \cite{Allen1992}. A particular set of paraxial modes with this property are the Laguerre-Gaussian modes (LG$_{q,l}$, with $q$ the radial index and $l$ the azimuthal one). This important finding opened up a whole new field, where the ability to independently control the AM and the polarization components of paraxial beams has allowed innumerable applications in quantum optics, microscopy, biology, optical trapping and astrophysics, just to mention a few \cite{Sonja2008}. The total AM of a vectorial field is composed of an orbital part ($\mathbf{L}$) and a spin part ($\mathbf{S}$), \textit{i.e.} $\mathbf{J}=\mathbf{L}+\mathbf{S}$. Nevertheless, it is important to note that these two components cannot ordinarily be separated. When both of these operators are applied to a Maxwell field, the result from that operation is not a Maxwell field \cite{Lifshitz1982,preIvan2012PRA}.

This paper presents two main findings. First, we show how to experimentally control the multipolar content of cylindrically symmetric modes. We advance that by controlling a component of the AM and the focusing of the incoming beam we can engineer fields which contain just few predetermined multipolar modes. Secondly, we demonstrate that by properly engineering these modes, we can also control the scattered field as well as address single resonances of an arbitrary, homogeneous and isotropic sphere in a lossless, homogeneous and isotropic medium. The second idea falls under the GLMT and is highly relevant to applications where control of the scattered field from a spherical object is desired. This scattering problem is solved by decomposing all the fields involved in the problem (incident $\mathbf{E}^\mathrm{i}$, scattered $\mathbf{E}^\mathrm{sca}$, and interior $\mathbf{E}^\mathrm{l}$) into multipoles and then use the boundary conditions. The three fields of the problem have the same formal expression, with the only difference that the multipolar modes of the scattered field have Hankel functions for the radial amplitudes, while both interior and incident fields have Bessel functions. However, the multipolar decomposition of the fields $\mathbf{E}^\mathrm{sca}$ and $\mathbf{E}^\mathrm{l}$ have all their complex amplitudes modulated by a Mie coefficient. These Mie coefficients are obtained by applying the Maxwell boundary conditions and depend on the AM as well as the parity. Thus, the only technicality of the problem is the projection of the incident field into the basis of multipoles.

The cylindrical symmetry of the incident beam simplifies the problem since the other two fields are required to possess a well defined value of $J_z$, too. In this way, we can solve the problem using a mixture of analytical and numerical techniques. More importantly, by exploiting the symmetries of the system we can get a further insight on the problem only using considerations of conservation of AM. Previous works considered similar problems where spherical metallic particles were excited with Laguerre-Gaussian beams \cite{vandeNes07}. There, the use of purely numerical techniques allowed the authors to distinguish certain properties of the scattered field. Our approach is different and more general. We extend the previous results to excitation with any cylindrically symmetric field and show a clear way of controlling the scattering from spherical particles. The results presented in this paper allow for a better understanding of recent experiments of scattering of silica spheres with Laguerre-Gaussian beams \cite{Gabi2012OL}, as our formulation also includes the paraxial case. In fact, they could also help to expand the computations of the optical forces when the incident beam is not a plane wave \cite{Stilgoe2008} and when the particle is not much smaller than the wavelength \cite{Nieto-Vesperinas2010, Novitsky2011}. Last but not least, our technique to control the multipolar content of cylindrically symmetric beams along with its immediate application to excite single multipolar modes could open up new possibilities in the forefront research of magnetic resonances \cite{Garcia-Etxarri2011, Kuznetsov2012, Miroshnichenko2012, Filonov2012}.

\section{Control of the multipolar content of a focused beam}

The experimental set-up which we propose is as follows. First, we prepare a monochromatic paraxial beam with a well defined $J_z$ controlling its polarization and spatial properties separately. This can be typically achieved with either spatial light modulators or holograms and a set of waveplates and polarisers. After this step, we will focus the resulting beam with a lens which we model as aplanatic \cite{Novotny2006}. The resulting focused beam interacts with the sphere and we analyse the scattering field. From now on, we will express our results in terms of the vector potential in the Coulomb gauge. All electric fields can be recovered using the typical formulas, \textit{i.e.} $\mathbf{E}=ik\mathbf{A}$ \cite{Rose1955}. Hence, the expression of the paraxial beam previously mentioned is $\mathbf{A}^{\mathrm{i}}(\rho,\phi,z=z_1)=F_l(\rho)\exp(il\phi)\exp(ikz_1)\mathbf{e}_p$, where $\{\rho,\phi,z\}$ are the spatial cylindrical coordinates, $z_1$ is the transverse plane right before the lens, $k$ is the modulus of the wavevector, $F_l$ is an arbitrary function representing the complex radial amplitude and $\mathbf{e}_p=(\hat{x}+i\,p\hat{y})/\sqrt{2}$ is a unit polarization vector representing right ($p=+1$) or left ($p=-1$) circular polarization. Here, we are disregarding terms in the amplitude of the order $1/(kd)^2$, where $d$ is the typical size of our light beam, which in our case is approximately the diameter of the focusing lens input aperture. Any paraxial beam can be decomposed as a superposition of this set of modes \cite{Gabi2001}. Note that a paraxial beam with circular polarization has, with a very good approximation, a well defined helicity $\Lambda = \mathbf{J} \cdot \mathbf{P} / \vert \mathbf{P} \vert $, where $\mathbf{P}$ and $\mathbf{J}$ are the linear and angular momentum operators respectively \cite{preIvan2012}. As mentioned before, then the paraxial beam impinges on the back of the microscope objective and it is focused. Microscope objective manufacturers make the objectives so that they behave approximately as aplanatic lenses. In addition to this, they coat the lenses in a way such that the transmission coefficients for the $\hat{\mathbf{s}}$ and $\hat{\mathbf{p}}$ polarized waves \cite{Morse1953} are equal, \textit{i.e.} $t^s(\theta)=t^p(\theta)$. Therefore, specifying the aplanatic model described in \cite{Novotny2006} for this case, and taking into account that this model preserves the helicity and the AM of the incident beam \cite{preIvan2012PRA}, we can write the field at the focus of the lens as
\begin{eqnarray}
\label{Afoc}
& &\mathbf{A}^{\mathrm{foc}}(\mathbf{r}) =  C
\int_{0}^{2\pi} d\phi_k \int_0^{\theta_M} \sin{\theta_k} d\theta_k \sqrt{n \cos \theta_k} \times \nonumber\\
 & & F_l(f \sin\theta_k ) \exp(i (l+p) \phi_k) \mathbf{M}(\theta_k,\phi_k)\mathbf{e}_p \exp(i \mathbf{k} \mathbf{r})
\end{eqnarray}
where $C=f \exp(-ikf)/2\pi$; $f$ is the focal of the lens; $\theta_M=\arcsin(\mathrm{NA})$, with NA the numerical aperture of the lens; $\mathbf{k}=k( \sin\theta_k \cos\phi_k, \sin\theta_k \sin\phi_k, \cos\theta_k)$; $n$ is the index of refraction of the lens; $\mathbf{r}=(\rho \cos(\phi),\rho \sin(\phi),z)$ is the position vector; and $\mathbf{M}(\theta_k,\phi_k)$ is the rotation matrix applied to the polarization vector so that the transversality of the field is fulfilled. Note, that field in Eq. (\ref{Afoc}) fulfils the exact Maxwell equations.

We will now describe how these beams will allow us to control the scattered field off a sphere. As previously mentioned, this problem is solved once the incident field ($\mathbf{A}^{\mathrm{foc}}$ in our case) is decomposed into multipolar modes. This decomposition can be made using the method put forward in \cite{Gabi2008}:

\small{
\renewcommand{\arraystretch}{2.5}
\begin{equation}
\begin{array}{ccl}
&&\mathbf{A}^{\mathrm{foc}}= \displaystyle\sum_{j=| l+p |}^{\infty} i^j (2j+1)^{1/2} C_{jlp}\left[ \mathbf{A}_{j(l+p)}^{(m)}+ip\mathbf{A}_{j(l+p)}^{(e)} \right]\\
&&C_{jlp}= C \displaystyle\int_0^{\theta_M}\sin\theta_k d\theta_k d_{(l+p)p}^j(\theta_k)F_l(f \sin\theta_k) \sqrt{n \cos \theta_k}
\end{array}
\label{decomp}
\end{equation}
\renewcommand{\arraystretch}{1}
\normalsize
As previously stated, the lens does not break the cylindrical symmetry of the beam, meaning that after passing through the lens, the beam still carries the same AM component $J_z=m_z$. Thus, one can see in the decomposition of Eq. (\ref{decomp}) that $m_z$ is the same for all multipoles and corresponds to $m_z=l+p$, which bounds the value of $j\ge|l+p|$. The coefficients of the decomposition $C_{jlp}$ only depend on a one dimensional integral which averages the radial complex amplitude of the input beam ($F_l$) over the reduced rotation matrix $d_{mp}^j(\theta)$, which can be found in \cite{Rose1957}. Finally, note that the beam still has a well defined helicity with value $p$. Indeed, if we apply the helicity operator, which can be written as $\Lambda = (\nabla \times)/k$ in differential form \cite{Tung1985}, to $\mathbf{A}^{\mathrm{foc}}$ in Eq. (\ref{decomp}) we obtain $\Lambda \mathbf{A}^\mathrm{foc} = p \mathbf{A}^\mathrm{foc}$. This is a consequence of the following two relations \cite{Rose1955}:
\begin{equation}
\Lambda \mathbf{A}_{jm_z}^{(m)}= i \mathbf{A}_{jm_z}^{(e)}, \quad \Lambda \mathbf{A}_{jm_z}^{(e)}= -i \mathbf{A}_{jm_z}^{(m)}
\label{TM-TE}
\end{equation}

\begin{figure}[htbp]
\centering\includegraphics[width=12cm]{figure1_v2.eps}
\caption{Multipolar decomposition ($\vert C_{jlp} \vert^2$) for different cases. The insets represent the intensity plots of the modes used for each simulation. The red coloured bars indicate NA=0.25, and the blue ones NA=0.9. The multipolar decomposition of (a) LG$_{0,0}$ and (b) LG$_{0,3}$ is presented. Note that $C_{jlp} $ can be described with very few multipoles when we use a high NA microscope objective to focus the beam.}
\end{figure}
Figure 1 exemplifies one of the main results of this paper. It shows the multipolar decomposition of different focused beams. Two different input beams were used: a Gaussian beam and a Laguerre-Gaussian with optical charge $l=3$, both right circularly polarized. The width of both beams was chosen so that the entrance pupil was filled. We consider two aplanatic lenses with NA's equal to 0.25 and 0.9. Our results show that the higher the NA, the narrower the distribution of multipolar modes is. The other parameter to control the distribution of multipoles is the $J_z$ of the incident beam. As mentioned earlier, the distribution of multipolar amplitudes is zero for values of $j<|l+p|$. This allows us to control the multipolar content of an optical field.

\section{Excitation of single multipolar modes}

The decomposition in Eq. (\ref{decomp}) and the use of the results of Mie scattering are the needed elements to calculate the scattered and interior fields,
\small{
\begin{equation}
\begin{array}{l}
\mathbf{A}^\mathrm{sca}=-\displaystyle\sum_{j=\vert l+p \vert}^{\infty} i^j (2j+1)^{1/2}  C_{jlp} \left[ b_j\mathbf{A}_{j(l+p)}^{(m)}+ipa_j\mathbf{A}_{j(l+p)}^{(e)} \right]  \\
\mathbf{A}^\mathrm{l}=\displaystyle\sum_{j=\vert l+p \vert}^{\infty} i^j (2j+1)^{1/2} C_{jlp} \left[ c_j\mathbf{A}_{j(l+p)}^{(m)}+ipd_j\mathbf{A}_{j(l+p)}^{(e)} \right]
\end{array}
\end{equation}
\normalsize
where $\left\lbrace a_j,b_j,c_j,d_j\right\rbrace$ are the Mie coefficients defined in \cite[Chap. 4]{Bohren1983}. These coefficients only depend on the size parameter of the problem ($x=2 \pi R/\lambda$), and both the relative permeability ($\mu_r$) and permittivity ($\epsilon_r$) of the sphere with respect to the surrounding medium. $R$ is the radius of the particle. Note that the multipolar fields in the expansion of $\mathbf{A}^\mathrm{sca}$ and $\mathbf{A}^\mathrm{l}$ contain different radial functions \cite{Bohren1983}. This completes the solution of the problem, as we can now easily compute the value of the fields in all regions. In particular, we can calculate the scattering efficiency, which is a dimensionless measure of the scattered power by the sphere in the far field \cite[Chap. 3]{Bohren1983}. Using the incident beam given by  Eq. (\ref{decomp}), the following expression is obtained for the scattering efficiency:
\begin{equation}
Q^\mathrm{sca}=\sum_{j={ \vert l+p \vert}}^{\infty} \frac{\left( 2j+1 \right)}{x^2} \vert C_{jlp} \vert ^2 \left( \vert a_j \vert ^2 + \vert b_j \vert ^2 \right)
\label{Csca}
\end{equation}
\begin{figure}[htbp]
\centering\includegraphics[width=12cm]{figure2_v2.eps}
\caption{Suppression of background in the scattering efficiency ($Q^{sca}$). Input beams are a) plane wave and b) LG$_{0,18}$. In a) no aplanatic lens is used, whereas in b) an aplanatic lens with a NA=0.9 is used. The rest of parameters are kept constant for both plots, \textit{i.e.} $R=1.3 \ \mu$m, $n_r=\sqrt{\epsilon_r \mu_r} = 1.5$, $p=1$. The insets indicate the typical profile of the beam used to excite the sphere. The scattering efficiency is represented with a blue continuous line. The Mie coefficients $b_{20}$ and $a_{22}$ are plotted with a red and a green line respectively. In a) we have indicated with a dashed line the position of two particular resonances for these modes. Note that the ordinate axis in a) and b) are different.}
\end{figure}
In comparison, the expression of the scattering efficiency for an incident plane wave \cite[Chap. 10]{Jackson1998}, has all $\vert C_{jlp} \vert ^2$ equal to 1 and the sum starts at $j=1$. This means that in conventional Mie Theory, all the Mie coefficients (and consequently all the multipolar modes) have the same relative weight. In our case, the use of an engineered input beam allow us to control the relative weight of the different multipolar modes (given by $a_j$ and $b_j$) with the use of the multipolar content of the incident field (given by $C_{jlp}$). Then, a tightly focused incident beam will allow us to select a few relevant Mie scattering coefficients. On the other hand, the use of beams with a well defined value of $J_z$ allows us to exclude the first $\vert m \vert$ multipolar modes from the scattering cross section.

In order to show these effects, in Fig. 2(a) we plot the conventional Mie scattering efficiency when the incident field is a plane wave, as well as the values of a pair of high order Mie coefficients ($b_{20}$ and $a_{22}$). In this case, the cross section is dominated by the low order multipoles. The higher order resonances appear as little ripples on top of the large oscillation. This is the so-called ripple structure. However, if an engineered cylindrically symmetric beam is used, the lowest $\vert m \vert$ modes can be removed and thus a particular Mie coefficient of interest can be excited. This is shown in Fig. 2(b), where the scattering efficiency associated to a LG$_{0,18}$ beam is depicted, as well as the same two Mie coefficients plotted in Fig. 2(a).
\begin{figure}[htbp]
\centering\includegraphics[width=12cm]{figure3_v2.eps}
\caption{Effect of a) the input beam and b) the NA of the lens on the excitation of a single resonance in the scattering efficiency ($Q^{sca}$). In a), three different beams are used: a plane wave and two LG's beams with $l=4,8$ focused with a lens with NA=0.25. In b) we use a LG$_{0,18}$ (same as Fig. 2(b)) but we focus it with different NA=0.1, 0.3, 0.5 and 0.9. The beams always fill the entrance pupil of the lens. Note that in a) the ordinate axis is linear and in b) logarithmic. Also, in a) the efficiency produced by each of the LG beams is multiplied by a different factor so that all the curves are comparable in the same plot. In b) the increase of the floor of the plotted lines is not due to an increase of the background, but rather to an amplification of the tails of the resonances.}
\end{figure}
\begin{figure}[htbp]
\centering\includegraphics[width=12cm]{figure4_v2.eps}
\caption{Projection into the helicity basis of the intensity of the interior and total field from a sphere with parameters $\left\lbrace R=1.3 \ \mu \mathrm{m},n_r=1.5 \right\rbrace $. The helicity of the incident beam is always $p=1$. The intensity is plotted in logarithmic scale, where 0 corresponds to the maximum of intensity in a).  Images a), b), e) and f) have been simulated for a $\lambda_{\mathrm{resonance}}=503$ nm, whereas a $2$ nm shift has been introduced in the off-resonance wavelength for the others, \textit{i.e.} $\lambda_{\mathrm{off-resonance}} = 505$ nm. The simulations for the LG$_{0,0}$ beam have been executed with a NA=0.25 lens  in order to be able to excite the resonance of order 20. Note that the scattering is dominated by the off-resonance modes. Calculations with higher numerical apertures give similar results. On the contrary, a lens with a NA=0.9 has been used for the excitation with a LG$_{0,18}$. In both cases, the entrance pupil of the lens was filled. Note that, in contrast with b), the mode of order 20 is greatly enhanced when we tune the wavelength to the resonance.}
\end{figure}

One interesting possibility which raises from the ability to freely excite multipolar resonances is that one can address very sharp resonances. The Mie scattering coefficients oscillate from 0 to 1 in a very exotic fashion. This behaviour is given by spherical Bessel functions. Every resonance of a Mie coefficient is given by a zero of a transcendental equation involving Bessel functions of order $j$. These equations have an infinite number of zeros, giving raise to an infinite number of resonances for a given order $j$. In general, it is very difficult to compare the shape of different resonances. Nevertheless, the first peaks of the Mie coefficients do follow a general trend. The higher the order of the Mie coefficient, the sharper the resonance. This implies that the high multipolar resonances of the sphere have smaller linewidths and consequently very large Q factors.  This also means that this sphere will be much more sensitive to the incident radiation. This fact is used in order to excite whispering gallery modes in large spheres, reaching very high Q factors for high values of $j$ \cite{Schiller1991}. In theory, Q could reach values as high as $10^{400}$. Yet, due to radiative losses of the material, the maximum values of Q experimentally achievable are of the order of $Q=10^9$ \cite{Oraevsky2002}.  Although this feature should also be present in the scattering of light from spheres, it is typically hidden due to the background produced by the lower order modes.

As shown in Fig. 3(a), our technique allows to excite a single multipolar resonance with free-space radiation. We fix our attention to the little ripple created by the Mie coefficient $b_{20}$ in the plane wave scattering (also present in Fig. 2(a)). The figure shows how to reduce the background created by the excitation of the lower modes. This has a power drawback though. Increasing the AM of the incident mode reduces the absolute value of the scattering cross section if the NA of the lens is preserved. However, if we increase the NA, the retrieved scattered power increases by many orders of magnitude and the shape of the resonance is maintained. This is shown in Fig. 3(b). Thus, high NA microscope objectives in conjunction with high AM modes can be used to experimentally address single multipolar resonances.

In Fig. 4 we show the interior and the total electric field (scattered plus incident) surrounding the sphere for two different incident beams. Two sets of four plots corresponding to two different incident beams (LG$_{0,0}$ and LG$_{0,18}$) at the same combination of $\left\lbrace R=1.3 \  \mu m,n_r=1.5 \right\rbrace $ are presented. Both incident fields have a well defined helicity $\Lambda^+$ (or $p=1$) in the focus of the lens. Thus, the resonant behaviour of the scattered field will be specially clear when the opposite component of the helicity $\Lambda^-$ is examined, since the contribution due to the incident field is 0. Indeed, if we compare the resonant excitation shown in Fig. 4(f) with the off-resonant behaviour of Fig. 4(h), it can be seen that the intensity differ in two orders of magnitude. Note that as it was shown in Fig. 2, the resonance depicted here is driven by the $b_{20}$ Mie coefficient, hence the single excited multipolar mode is $\mathbf{A}_{20,19}^{(m)}$. Another consequence of the single excitation of a multipolar mode in Fig. 4(e-h) is the fact that in Fig. 4(e-f) the two helicity components of the scattered field have the same power. This is due to the fact that the we are exciting a single magnetic multipole $\mathbf{A}_{20,19}^{(m)}$ and that the multipoles considered in Eq. (\ref{TM-TE}) are a superposition of multipoles with a well defined helicity \cite{preIvan2012PRA}. Finally, it can be observed that unlike Fig. 4(e-h), where the field can be described with a single mode, many multipoles are needed to describe the interferences seen in Fig. 4(a-d). In fact, the pattern-like observed in Fig. 4(a-h) has drawn a lot of attention since 2004, where it was first characterized as photonic nanojet \cite{Chen2004}. Since then, nanojets have been extensively studied and applied in many different fields \cite{Heifetz2009}.

\section{Conclusion}

We have demonstrated how to address multipolar resonances and how to control the scattered field of a spherical object by engineering the incident light. This technique allows to excite single multipolar resonances and also, to choose resonances with very high Q factors even for spheres with sizes of the order of the wavelength of the incident radiation. Besides the fundamental control of material particles that this technique grants, we foresee applications in the fields of cytometry and dark field microscopy, where the control of the AM of the illuminating beams has not been  explored until now.

\section*{Acknowledgements}
This work was funded by the Australian Research Council Discovery Project DP110103697. G.M.-T. is the recipient of an Australian Research Council Future Fellowship (project number FT110100924). We are thankful to I. Fernandez-Corbaton for some very useful discussions.

\end{document}


\section{Introduction}

The idea of extra dimension was proposed to unify different forces. The representative one is the Kaluza-Klein (KK) theory, which unifies the 4-dimensional (4D) gravity and electromagnetism by a 5-dimensional Einstein's theory with a spatial extra dimension. And the string/M-theory also requires that there are 6/7 extra dimensions. But the extra dimensions in these theories are so extremely small that can not be detected by the present experiments.

When the Arkani-Hamed-Dimopoulos-Dvali (ADD) \cite{Antoniadis1998} and the Randall-Sundrum (RS) \cite{Randall1999a,Randall1999b} brane-world models were put forward, the extra dimension theory was revived. One of the reasons is because the scale of extra dimensions can be millimetre or even infinite, which increases the possibility to be detected in experiment. Another one is because the theory opens a new way to solve the long-standing hierarchy problem and the cosmology problem \cite{ArkaniHamed1998rs,Cosmological2000,coscon2009,Cosmoligicalbrane2010}.

{In the Randall-Sundrum model, the Standard Model (SM) fields are confined on the brane, and gravity propagates in the whole space-time. This idea is motivated by the string theory, in which the SM fields are excitation modes of open string whose ends lie on the brane, and graviton is a closed string excitation having no-end points, thus it does not have to be attached to the brane. Moreover, the $q-$form fields also comprise the excitation of a closed string \cite{stringvol.1,stringvol.2}. So it is reasonable to introduce the $q-$form fields in the RS (or RS-like) brane, and see the localization of these fields on different branes, including the thin and thick branes \cite{Rubakov1,Rubakov2,ThickBraneDewolfe,Thickbrane2000,
ThickBrane2001,ThickBrane2002,ThickBrane2003,Blochbrane2004,
ThickBraneBazeia2006,ThickBrane20071,ThickBrane20082,ThickBrane20093,
KeYang2009Weyl,ThickBraneZhongYuan2011JHEP,ThickBraneKeYang2011PRD,
ThickBraneCriticalGravityLiu2012,BraneFR2012}.

On the other hand, although free q-form fields in the bulk are just equivalent to scalar or vector fields by a duality in 4D space-time, they are new types of particles in the higher space-time.

In this work, we consider the branes with $p$ spacial dimensions in the $D=p+2$ dimensional space-time, which are called $p-$branes as in the string/M-theory. The realistic world is the $3-$brane, and the higher dimensional branes with $p\geq5$ ($D\geq7$) may also have realistic applications if the branes have 3 infinite large dimensions (which are those we can feel) and $p-2$ finite size dimensions with topology $S^1 \times S^1 \times ... \times S^1 = T^{p-2}$ and small enough radius. Thus we will investigate the localization of $q-$form fields on these $p-$branes.}

The introduction of extra dimensions implies that there are KK modes accompanying all the particles on the brane \cite{WarpedPassages}. These KK modes carry the messages of the extra dimensions. If the KK modes can be localized on the brane, there is hope to detect them and prove the existence of the extra dimensions, which is one main meaning of  localization. The massless KK modes exactly stand for the particles on the brane, and surely must be localized for a realistic brane-world model. And the next heavier massive modes are the key clue.

In Refs.~\cite{LocalizationPRDLanglois2003,
LocalizationGravityJHEP2005,
Liu2007WeylVo,Liu2008WeylPT,LocalizationMO2008,JHEPadsGuo2009,
LocalizationFermionSplit2010,
LocalizationFuPRD2011,PRDefromedLiu,LocalizationCastro2011,
LocalizationZhao2011JHEP,LocalizationWithoutScalar2011JHEP,
LocalizationWaveBrane2012,LocalizationFermion2012}, the localization of the scalar, vector and fermion fields has been investigated. In this paper, we try to find one general method for localizing any form ($q-$form) field. For $q=0,1$ the fields are just the scalar and vector ones, and for higher $q$ they stand for new particles in the higher space-time \cite{Quevedo2010ui-eprintv1}. In the RS model, the localization of $2-$form and $3-$form fields has been investigated \cite{PRLKR2002,Mukhopadhyaya2004,QformRS,qRSdilation-eprintv2}, and it was shown that it is possible to find some signal about the fields coupled with a dilaton at Tev scale. Then the authors discussed the localization of the $q-$form fields on a domain wall in a RS-like scenario \cite{qdual2010-eprintv3}.

While our aim is to find the general method for localizing the $q-$form fields, and then apply this to some typical branes. In the following, we first study how to obtain the effective potential of KK modes for the $q-$form fields in Sec. \ref{SecLocalize}. Then, with some typical $p-$brane models in Sec. \ref{pbranes}, we investigate the condition for localizing the $q-$form fields and the effect of the dimension of the $p-$branes. Finally, we give a simple conclusion in Sec. \ref{conclusion}.

\section{How to localize $q-$form fields}
\label{SecLocalize}

We consider a free antisymmetric $q-$form field $X_{M_1M_2...M_q}$ on $p-$branes in $D-$dimensional space-time with $D=p+2$, and treat the field as a small perturbation around the background. The action for the field is
\begin{eqnarray}\label{qaction}
S_q=\int d^D x\sqrt{-g}\;Y_{M_1M_2... M_{q+1}}Y^{M_1M_2... M_{q+1}},
\end{eqnarray}
where $Y_{M_1M_2...M_{q+1}}$ is the field strength defined as $Y_{M_1M_2...M_{q+1}}=\partial_{[M_1}X_{M_2...M_{q+1}]}$.
Because the $(D-1)-$form and higher form fields have no role to play in the brane \cite{QformRS}, we only analyze the cases $q=0,1,2,...,p$.

The line-element of the $D-$dimensional space-time is assumed as
\begin{eqnarray}\label{line-element}
 ds^2&=&\text{e}^{2A(z)}\big(\eta_{\mu\nu}dx^\mu dx^\nu
          + dz^2\big)
\end{eqnarray}
with $z$ denoting the extra dimension vertical to the brane. Then the equations of motion can be derived from the action (\ref{qaction}) and the conformal metric (\ref{line-element}):
\begin{eqnarray}
 \partial_{\mu_1} ( \sqrt{-g}\;Y^{\mu_1\mu_2...\mu_{q+1}})
 +\partial_z(\sqrt{-g}\;Y^{z\mu_2...\mu_{q+1}})&=& 0, \label{equ1}\\
 \partial_{\mu_1} ( \sqrt{-g}\;Y^{\mu_1...\mu_qz})&=& 0.
\end{eqnarray}

In order to investigate the localization of the $q-$form field on a $p-$brane, we first choose a gauge $X_{\mu_1...\mu_{q-1}z}=0$ using the gauge free $\delta X_{M_1M_2...M_q}=\partial_{[M_1}\Lambda_{M_2...M_q]}$ with $\Lambda_{M_2M_3...M_q}$ an antisymmetric tensor \cite{Mukhopadhyaya:2007jn}, and make a decomposition for the $X_{\mu_1\mu_2...\mu_q}$:
\begin{eqnarray}\label{KK1}
X_{\mu_1\mu_2...\mu_q}(x^\mu,z)=\sum_n
\hat{X}^{(n)}_{\mu_1\mu_2...\mu_q}(x^\mu)U^{(n)}(z) \text{e}^{(2q-p)A/2},
\end{eqnarray}
so the field strength becomes
\begin{eqnarray}
Y_{\mu_1\mu_2...\mu_{q+1}}(x^\mu,z)&=&\sum_n \hat{Y}^{(n)}_{\mu_1\mu_2...\mu_{q+1}}(x^\mu)
U^{(n)}(z) \text{e}^{(2q-p)A/2},\label{KK10}\\
Y_{\mu_1\mu_2...\mu_qz}(x^\mu,z)&=&\sum_n\frac{1}{q+1}\hat{X}^{(n)}_{\mu_1\mu_2...\mu_q}(x^\mu)
\left(U^{(n)}{'}+\frac{2q-p}{2}A'U^{(n)}\right)U^{(n)}(z) \text{e}^{(2q-p)A/2},~~~~~\label{KK11}
\end{eqnarray}
where $\hat{Y}_{\mu_1...\mu_{q+1}}(x^\mu)=\partial_{[\mu_1}\hat{X}_{\mu_2...\mu_{q+1}]}(x^\mu)$ is the field strength on the brane.

Substituting the relations (\ref{KK10}) and (\ref{KK11}) into the equations of motion (\ref{equ1}), we can get the following Schr\"{o}dinger-like equation for the KK modes of the $q-$form field:
 \begin{eqnarray}
\big[ -\partial^2_z+ V(z)\big]U^{(n)}(z)
  =m_n^2 U^{(n)}(z),
  \label{SchEq}
\end{eqnarray}
where $m_n$ are the masses of the KK modes, and the effective potential $V(z)$ takes the form
\begin{eqnarray}\label{V}
V(z)=\frac{(p-2q)^2}{4}A'^2(z) +\frac{p-2q}{2}A''(z)
\label{Veff}
\end{eqnarray}
with the prime standing for the derivative with respect to $z$.

Furthermore, provided the orthonormality condition
\begin{equation}\label{NormCondi}
  \int U^{(m)}U^{(n)}dz=\delta_{mn},
\end{equation}
the effective action of the $q-$form field on the brane could be obtained:
\begin{eqnarray}
S_{\text{eff}} = \sum_n \int d^{D-1}x\sqrt{-\hat{g}}
\bigg(\hat{Y}^{(n)}_{\mu_1\mu_2...\mu_{q+1}}\hat{Y}^{(n)\mu_1\mu_2...\mu_{q+1}}
+\frac{1}{q+1}m_n^2\hat{X}^{(n)}_{\mu_1\mu_2...\mu_q}\hat{X}^{{(n)}\mu_1\mu_2...\mu_q}
\bigg).
\end{eqnarray}

For $m_0^2=0$, we get the zero mode for the $q-$form field from (\ref{SchEq}):
\begin{equation}\label{U0}
  U_0\propto \text{e}^{\frac{p-2q}{2}A},
\end{equation}
which results in $X^{(0)}_{\mu_1\mu_2...\mu_q}(x^\mu,z)=
\hat{X}^{(0)}_{\mu_1\mu_2...\mu_q}(x^\mu)$ (see from (\ref{KK1})), and it is irrelevant to the extra dimension.
It means that the KK mode with $m^2=0$ (the zero mode) exactly reproduces the particle on the brane. It surely should be localized on the realistic branes. And we can judge this with the help of the  orthonormality condition (\ref{NormCondi}). So the messages of the extra dimension only can be disclosed by the massive KK modes. We could investigate the feature of them using the Schr\"{o}dinger-like equation (\ref{SchEq}). At last the core of the problem is focused on the effective potential in (\ref{SchEq}), which is decided by three factors, i.e., the parameters $q$, $p$, and the background space-time. And this is the clue of the following work.

\section{Localization of $q-$form fields on different $p-$branes}
\label{pbranes}

In the section we try to investigate the localization of the $q-$form fields on different types of flat $p-$branes, i.e., the thin RS $p-$branes, the thick ones with finite and infinite extra dimension, and the pure geometrical ones. So we first give the solutions of the $p-$branes. The effective potentials on these branes mainly have two types: the volcano-like type and the P\"{o}schl-Teller (PT)-like one. For the former type there exist only continuous massive KK modes, but for the later one some massive bound KK modes appear. The effect of the parameters $p$ and $q$ on the localization will also be discussed.

\subsection{Thin $p-$brane model}

Firstly, we consider thin $p-$brane (RS brane) in $D-$dimensional AdS space-time with a negative cosmological constant $\Lambda$. The action of this system is
\begin{eqnarray}
  S=\frac{1}{2\kappa_D^2}\int d^Dx\sqrt{-g}\big[R-p\;\Lambda\big]
        + \int d^{D-1}x \sqrt{-g^{(b)}}(-\sigma),
\end{eqnarray}
where $g^{(b)}_{\mu\nu}$ is the induced metric on the brane, $\sigma$ is the brane tension, and $\kappa_D^2=8 \pi G_D$ with $G_D$ the D-dimensional Newton constant. Here we let $\kappa_D=1$. Then the Einstein equation reads:
\begin{equation}
  G_{MN}+\frac{\Lambda}{2}p\;g_{MN}
  =-\sigma\delta^\mu_M\delta^\nu_N\eta_{\mu\nu}\text{e}^A\delta(z),
\end{equation}
where $G_{MN}=R_{MN}+\frac{1}{2}g_{MN}R$ is the Einstein tensor. Then the solution of the thin flat $p-$brane is
\begin{eqnarray}
  A(z)&=&-\ln{\big(1+\sqrt{\frac{-\Lambda}{p+1}}|z|\big)},\label{WarpRS}\\
  \sigma&=&2\;p\sqrt{\frac{-\Lambda}{p+1}}.
\end{eqnarray}

So in this brane model, the form of the potential (\ref{V}) is expressed as
\begin{eqnarray}
  V(z)=-\frac{(p-2q)(p-2q+2)\Lambda}{p+1}\frac{1}{4(1+k_1|z|)^2}
       -\frac{k_1(p-2q)}{1+k_1|z|}\delta(z)
\end{eqnarray}
with $k_1=\sqrt{\frac{-\Lambda}{p+1}}$. As $\Lambda$ is negative, the potential is volcano-like for $q<p/2$. And there is a zero mode for the $q-$form field with the wave function (\ref{U0}). We can check whether it can be localized on the brane through (\ref{NormCondi}):
\begin{eqnarray}
  \int U_0^2dz=\int (1+k_1 |z|)^{2q-p}dz=\frac{2}{k_1(2q-p+1)}(1+k_1 |z|)^{2q-p+1},
\end{eqnarray}
from which it is clear that, for $q<\frac{p-1}{2}$ the integral is finite, hence, the zero mode can be localized on the brane. For example, for $p=3$, only the scalar field $(q=0)$ can be localized. While for higher dimension case, the vector field can be as well. { And the $2-$form field can also not be localized on the $3-$brane, such as the $2-$form Kalb-Ramond (KR) field. This result corresponds to that in Refs.~\cite{PRLKR2002,Mukhopadhyaya2004,LocalizationFuPRD2011,PRDefromedLiu}. It is shown that the KR field can be considered as a source of torsion \cite{PRLKR2002,Mukhopadhyaya2004}, however, there is no phenomenon about the torsion of our space-time ($3-$brane). Here, we can give a possible reason for this, i.e., the KR field should be localized on the higher dimensional brane ($p\geq5$). And for the following brane model, we will obtain the same result for the KR field.}

\subsection{Thick flat $p-$branes model}

In this subsection, we will investigate the $q-$form field on thick Minkowski $p-$branes, which are generated by one or two scalar fields.

\subsubsection{Thick flat $p-$brane with infinite extra dimension}

Considering the thick $p-$brane generated by one scalar field $\phi$, the action is given by
\begin{equation}
S=\int d^Dx \sqrt{-g}
 \big[\frac{1}{2\kappa_D^2}R -\frac{1}{2}(\partial\phi)^2-V(\phi)\big].
\end{equation}
 With the metric
\begin{equation}\label{line-element-y}
ds^2=\text{e}^{2A}\eta_{\mu\nu}dx^\mu dx^\nu+dy^2
\end{equation}
the equations of motion can be get:
\begin{eqnarray}
  \phi'^2&=&-p\;A'',\\
  -2V&=&p\;(p+1)\;A'^2+p A'',\\
  dV/d\phi&=&\phi''+(p+1)A'\phi',
\end{eqnarray}
where the prime denotes the derivative with respect to $y$. Then we use the superpotential method to solve these equations. With the introduction of the superpotential function $W(\phi)$ defined as $\phi'=dW/d\phi$, the above second-order differential equations are reduced to the following one-order differential ones:
\begin{eqnarray}
  A'=-\frac{1}{p}W,~~~
  V=-\frac{1}{2}\left[\frac{p+1}{p}W^2-(\frac{d W}{d\phi})^2\right].
\end{eqnarray}
If we suppose that the superpotential is the sine-Gordon type, i.e., $W(\phi)=\frac{3}{2}bc\sin{\sqrt{\frac{2}{3b}}\phi}$ with $b$ and $c$ both constants, the solution of this flat thick brane can be obtained:
\begin{eqnarray}
  \phi(y)&=&\sqrt{6b}\arctan{\tanh{(cy/2)}},\\
  A(y)&=&\frac{3b}{2p}\ln{\text{sech}(cy)}\label{Solu-sineGorden}.
\end{eqnarray}

With the coordinate transformation $dz=\text{e}^{-A}dy$, the warp factor tends to $A(z)\rightarrow-\ln{(\frac{3bc}{2p}z+1)}$ when $z \rightarrow \pm\infty$. Then the behavior of the potential $V(z)$ at $z=0$ and at infinity is
\begin{eqnarray}
  V(z=0)&\rightarrow&-\frac{3b^2c}{4p}(p-2q),\\
  V(z\rightarrow\pm\infty)&\rightarrow&0.
\end{eqnarray}
For $q<p/2$, the potential is also a volcano-like one. So there is a zero mode for the $q-$form field. With the orthonormality condition (\ref{NormCondi})
\begin{eqnarray}
  \int U_0^2dz\rightarrow\int(1+\frac{3bc}{2p}z)^{-(p-2q)}dz,
\end{eqnarray}
we see that, for $q<\frac{p-1}{2}$, the zero mode can be localized on the brane, which is the same with the thin brane case.

\subsubsection{Thick flat $p-$brane with finite extra dimension}
\label{flatfinite}

Now let's turn to another kind of thick flat $p-$brane, which is generalized by two scalar fields $\phi$ and $\pi$. The action of this system is
\begin{eqnarray}
  S=\int d^D x \sqrt{-g}\bigg[\frac{1}{2\kappa_D^2}R-\frac{1}{2}\phi'^2-\frac{1}{2}\pi'^2
  -V(\phi,\pi)\bigg].
\end{eqnarray}
Then, with the conformally flat metric (\ref{line-element}), the equations of motion read
\begin{eqnarray}
  \phi'^2+\pi'^2&=&p(A'^2-A''),\\
  -2\text{e}^{2A}V&=&p^2A'^2+p A'',\\
  dV/d\phi&=&\text{e}^{-2A}\big[\phi''+p A'\phi'\big],\\
  dV/d\pi&=&\text{e}^{-2A}\big[\pi''+p A'\pi'\big].
\end{eqnarray}
Using the superpotential method, we can find that the following one-order differential equations are the solutions of the above ones:
\begin{eqnarray}
  A'=-\frac{1}{p}W, ~~~~\pi=\sqrt{p}\;A,~~~~
  V=-\frac{1}{2}\text{e}^{-\frac{2}{\sqrt{p}}\pi}
  \left[W^2-(\frac{\partial W}{\partial \phi})^2\right].
\end{eqnarray}
Supposing the form of the superpotential as $W(\phi)=va\phi(1-\frac{\phi^2}{3v^2})$ with $v$ and $a$ positive constants, we obtain the following solution:
\begin{eqnarray}
  \phi&=&v \tanh{(az)},\\
  A&=&-\frac{v^2}{3p}\left[\ln{\cosh^2{(az)}}+\frac{1}{2}\tanh^2{(az)}\right],\\
  \pi&=&\sqrt{p}\;A.
\end{eqnarray}
For this solution, the physical extra dimension $y$ is finite, which can be seen from the coordinate transformation $\int dy=\int \text{e}^{A}dz\propto \int \text{e}^{-k_2z}dz$ with $k_2$ a positive constant. And we have discussed the finity of the extra dimension in Ref.~\cite{PRDefromedLiu} in 5D case.

For this brane-world model, the effective potential $V(z)$ reads
\begin{eqnarray}
  V(z)=\frac{a^2v^2(p -2 q)}{36p^2}
   \bigg[-18 p\;\text{sech}^4{(az)} + (p-2q)v^2\tanh^2{(az)} \big[2 + \text{sech}^2(az)\big]^2 \bigg],
   \label{VzFiniteED}
\end{eqnarray}
whose asymptotic behavior is
\begin{eqnarray}
 V(z=0)&=&-\frac{a^2 v^2(p - 2 q)}{2 p},\\ V(z\rightarrow\infty)&\rightarrow&\frac{a^2v^4(p - 2 q)^2}{9 p^2}.
\end{eqnarray}
For $q<p/2$, the potential is a PT-like one, and there exist a zero mode and finite massive bound KK modes. We plot the shapes of the potential for different $q,p$ in Fig.~\ref{FigVfinite}.

\begin{figure*}[htb]
\begin{center}
\includegraphics[width=6cm]{V_FiniteExtra_q.eps}
~~~~\includegraphics[width=6cm]{V_FiniteExtra_p.eps}
\end{center}
 \caption{The shapes of the effective potential $V(z)$ (\ref{VzFiniteED}) for $v=5,a=1$. The left one is for different $q$ with $p=7$, and the right one is for different $p$ with $q=1$.}
 \label{FigVfinite}
\end{figure*}

For the zero mode, we have
\begin{equation}
  \int U_0^2dz=\int \text{e}^{(p-2q)A}dz\propto\int\text{e}^{-k_2(p-2q)z}dz.
\end{equation}
So the zero mode can be localized on the brane if $q<p/2$. For example, for the case of $p=3$, both the scalar field $(q=0)$ and the vector field $(q=1)$ can be localized on the brane.

While for the massive KK modes, we mainly focus on the effect of the parameters $p$ and $q$ on them. For the scalar field ($q=0$ ), the dimension of the brane $p$ will not affect the effective potential; so the mass spectra of the bound scalar KK modes are the same on different $p-$branes. But for other $q-$form fields, the depth of the effective potential increases with $p$, and so does the number of the massive bound KK modes. For example, for the vector field ($q=1$), the mass spectra are
\begin{eqnarray}
  m^2&=&(0, 6.20)\cup (7.71, \infty), ~~~~~~~~~~~~~~~~~~~\text{for}~~~p=3,\\
  m^2&=&(0, 12.93, 21.49)\cup (25, \infty), ~~~~~~~~~~~~\text{for}~~~p=5,\\
  m^2&=&(0, 15.80, 27.30, 34.05)\cup (35.43, \infty), \text{for}~~~p=7,
\end{eqnarray}
from which we see that the number of the massive bound KK modes increases with the dimension of the $p-$brane.

Moreover, on the same $p-$brane, the depth of the potential decreases with $q$, thus the number of the massive bound KK modes also decreases. For example, on the $7-$brane, there are more massive bound KK modes of the scalar field ($q=0$) than that of the higher form fields, which can be seen from the following mass spectra:
\begin{eqnarray}
  m^2&=&(0, 22.96,41.72,56.04,65.51)\cup (69.44, \infty), ~~~~\text{for}~~~q=0,\\
  m^2&=&(0, 15.80, 27.30, 34.05)\cup (35.43, \infty), ~~~~~~~~~~~~\text{for}~~~q=1,\\
  m^2&=&(0, 8.61)\cup (12.76, \infty), ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\text{for}~~~q=2.
\end{eqnarray}

\subsubsection{Pure geometrical thick $p-$brane}

At last we consider the pure geometrical thick $p-$brane. For the $D-$dimensional Weyl-integrable manifold $M_D^W$, which is specified by a pair of $(g_{MN}, \omega)$ with $g_{MN}$ the $D-$dimensinal metric and $\omega$ a Weyl scalar function, the action of the Kaluza-Klein theory is \cite{KeYang2009Weyl,BarbosaCendejas2006,Arias2002,BC2007WeylGravity}
\begin{equation}\label{actionW}
  S_D^W=\int_{M_D^W}d^Dx\sqrt{-g}\;\text{e}^{-\frac{p}{2}\omega}
  \left[\frac{1}{2\kappa_D^2}R-\frac{3}{2}\xi(\nabla\omega)^2-3U(\omega)\right],
\end{equation}
where $\xi$ is the coupling constant and $U(\omega)$ is the self-interaction potential of the Weyl scalar function $\omega$. By performing a conformal transformation $\hat{g}^{MN}=\text{e}^{-\omega}g^{MN}$, the action (\ref{actionW}) is mapped into the Einstein frame:
\begin{eqnarray}\label{actionR}
  S^R_D=\int d^Dx\;\sqrt{-\hat{g}}\left[\frac{1}{2\kappa_D^2}\hat{R}-\frac{3}{2}\xi(\hat{\nabla}\omega)^2
  -3\hat{U}(\omega)\right]
\end{eqnarray}
with $\hat{U}=\text{e}^{\omega}U$. Then we can obtain one solution in a compact manifold along the extra dimension with the range $-\frac{\pi}{2}\leq ky\leq\frac{\pi}{2}$:
\begin{eqnarray}
  \omega&=&-\frac{2c(p+1)}{p}\ln{\cos{(k y})},\\
  A&=&c\;\ln{\cos(k y)},
\end{eqnarray}
where $c=\frac{12p^2}{(p+1)[p-12\xi(p+1)]}$, and $k$ is a positive constant. Here, in order to get a finite warp factor $\text{e}^{2A}$ everywhere, we need $c>0$. Further, in order for the photon to spend infinite time from the brane to the boundary, we require $c\ge 1$. The solution in 5-dimensional space-time can be found in Refs.~\cite{Liu2008WeylPT,BarbosaCendejas2006}.

Through the coordinate transformation $dz=\text{e}^{-A}dy$, the effective potential (\ref{V}) can be rewritten in $y$-coordinate as
\begin{eqnarray}\label{VW}
  V(z(y))&=&\text{e}^{2A}\left[ \frac{(p-2q)^2+2(p-2q)}{4}(\partial_y A)^2
  +\frac{p-2q}{2}\partial_{y}^2 A\right]\nonumber\\
      &=&\frac{1}{4}c(p-2q)k^2\cos^{2(c-1)}(ky)\big[c(2+p-2q)\sin^2{(ky)}-2\big],
\end{eqnarray}
from which we can see that the potential $V(z(y))$ has different types of shapes decided by the parameter $c$,
and it is the same for $V(z)$. We plot the shapes of the $V(z(y))$ and $V(z)$ for different $c$ in Fig.~\ref{figVWeyl}.
Note that, for $c\ge 1$, the range of the $z$ is infinite; while for $0<c < 1$,
the range of the $z$ is finite, i.e., $-z_c <z<z_c$ with $z_c=\sqrt{\pi}~\Gamma(\frac{1-c}{2})/2\Gamma(1-\frac{c}{2})$.

\begin{figure*}[htb]
\begin{center}
\subfigure[$c>0~(-\pi/2<ky<\pi/2)$]{\label{figVWeyla}
\includegraphics[width=6cm]{Vy_Weyl.eps}}
\subfigure[$c\ge 1~ (-\infty<z<\infty)$]{\label{figVWeylb}
\includegraphics[width=6cm]{Vz1_Weyl.eps}}
\subfigure[$\frac{2}{2+p-2q}{\le}c<1 ~ (-z_c <z<z_c)$]{\label{figVWeylc}
\includegraphics[width=6cm]{Vz2_Weyl.eps}}
\subfigure[$0<c<\frac{2}{2+p-2q} ~ (-z_c <z<z_c)$]{\label{figVWeyld}
\includegraphics[width=6cm]{Vz3_Weyl.eps}}
\end{center}
 \caption{The shapes of the effective potential $V(z(y))$ and $V(z)$ for pure geometrical $p-$brane with $k=1$, $p=3$, $q=1$, and $c=1.4$ for the thick lines, $c=1$ for the thin solid lines, $c=0.76$ for the dashed lines, $c=0.67$ for the dotted lines, and $c=0.1$ for the dot dashed lines.}
 \label{figVWeyl}
\end{figure*}

From (\ref{VW}) it is clear that, for $c>1$ and $q<p/2$, $V(y=0)=-{ck^2}(p-2q)/{2}<0$ and $V(ky\rightarrow\pm\pi/2)\rightarrow 0_{+}$, and the effective potential is a volcano-like one. So there exists a zero mode for $q<p/2$. With the orthonormality condition (\ref{NormCondi})
\begin{eqnarray}\label{zeroMW}
\int_{-z_b}^{z_b}U_0^2dz&=&2\int_{0}^{z_b} \text{e}^{(p-2q)A}dz\nonumber\\
&=&\int_0^{\frac{\pi}{2k}} \text{e}^{(p-2q-1)A}dy
=\int_0^{\frac{\pi}{2k}} \cos^{w}{(ky)} dy\nonumber\\
&=&
  \frac{\sqrt{\pi}~\Gamma(\frac{1+w}{2})}{2\Gamma(1+\frac{w}{2})}
  - \frac{\cos^{1+w}{(\pi/2)}}{1+w},
\end{eqnarray}
where $w=c(p-2q-1)$, it is found that when $q<\frac{p-1+\frac{1}{c}}{2}$ the above integral is finite, hence, the zero mode can be localized on the brane. Considering $c>1$, the final result is that the $q-$form fields with $q<p/2$ have a localized zero mode, which is different from that on the thin $p-$brane and on the one with infinite extra dimension.

In fact, for the zero mode of the $q-$form field, the $\int U_0^2dz$ in the physical coordinate system becomes $\int \text{e}^{(p-2q-1)A}dy$; it is always finite for $q=(p-1)/2$ and finite $y$. So the zero mode of the fields with $q=(p-1)/2$ can always be localized on the brane with finite extra dimension. For the pure geometrical thick brane with compact extra dimension considered here, the $q-$form field with $q=(p-1)/2~(<p/2)$ can also be localized on the brane.

While for the case $c=1$, the potential is a PT-like one for $q<p/2$ (see Fig. \ref{figVWeylb}), whose behavior is $V(y=0)=-k^2(p-2q)/2 <0$ and $V(ky\rightarrow\pm\pi/2)\rightarrow k^2(p-2q)^2/4>0$, so there must exist a zero mode and some massive bound KK modes. From (\ref{zeroMW}), it is clear that the zero mode can be localized for $q<p/2$. And the number of the massive bound KK modes of the $q-$form field increases with the dimension of the $p-$brane. Here, the parameter $p$ also effects the number of the massive bound scalar KK modes, which is different with that on the thick flat $p-$brane with finite extra dimension discussed in section \ref{flatfinite}.

For $\frac{2}{2+p-2q}<c<1$ and $q<p/2$, $V(y=0)=-c(p-2q)k^2/2<0$ and $V(ky\rightarrow\pm\pi/2)\rightarrow\infty_+$, so the effective potential is an infinite potential well (see Fig. \ref{figVWeylc}). Thus, there are a bound zero mode and infinite number of massive bound KK modes, all these KK modes can be localized on the brane.

For $c=\frac{2}{2+p-2q}$ and $q<p/2$, we will get a negative PT-like potential (see Fig. \ref{figVWeylc}), and the zero mode can be localized on the brane, which can be seen from (\ref{zeroMW}).

We do not consider the case of $0<c<\frac{2}{2+p-2q}$, for which the potential approaches negative infinity at the boundary (see Fig. \ref{figVWeyld}).

\section{Conclusion}
\label{conclusion}

In this paper, we first investigated one general method of localizing the $q-$form fields on $p-$branes. We found that the KK modes of the $q-$form fields satisfy the Schr\"{o}dinger-like equation. Thus, the characters of the KK modes are decided by the behavior of the effective potential, which depends on the parameters $q$ and $p$ and the background geometry.

Then we applied this to some typical $p-$brane models. It was found that for the thin Minkowski $p-$brane the effective potential is a volcano-like one, which results in that there exist a zero mode and continuous massive KK modes, and only for $q<(p-1)/2$ the zero mode of the $q-$form field can be localized on the brane. And it is the same for the thick Minkowski $p-$brane with infinite extra dimension. While for the thick Minkowski $p-$brane with finite extra dimension, there is a PT-like potential. So a zero mode and some massive bound KK modes can exist. For the fields with $q<p/2$, the zero mode can be localized on the brane. The number of the massive bound KK modes for the same $q-$form field increases with the dimension of the $p-$brane, except for the scalar field ($q=0$). Moreover, on the same $p-$brane, the $q-$form fields with higher $q$ have less number of the massive bound KK modes.

For the pure geometrical thick $p-$brane, the shape of the effective potential is decided by the value of the parameter $c$. For $c>1$ and $q<p/2$, the effective potential is a volcano-like one, and the condition for localizing the zero mode of the $q-$form field is $q<p/2$, which is different from that on the thin $p-$brane and on the one with infinite extra dimension. For $c=1$ and $q<p/2$, the effective potential is a PT-like one, and only the fields with $q<p/2$ have a localized zero mode. The number of the massive bound KK modes for any $q-$form field, including the scalar field, increases with the dimension of the $p-$brane. And for $\frac{2}{2+p-2q}<c<1$ and $q<p/2$, the effective potential is an infinite potential well, and there exist a localized zero mode and infinite number of massive bound KK modes. For $c=\frac{2}{2+p-2q}$ and $q<p/2$, there is a negative PT-like potential, and the zero mode can also be localized on the brane. However, we only need the case of $c\ge 1$ in order to secure a fully predictive initial value Cauchy problem \cite{Mannheim2005}.

\section*{Acknowledgement}

This work was supported by the Program for New Century Excellent Talents
in University, the National Natural Science Foundation of China (No. 11075065),
the Doctoral Program Foundation of Institutions of Higher Education of China
(No. 20090211110028), the Huo Ying-Dong Education Foundation of
Chinese Ministry of Education (No. 121106), and the Fundamental Research Funds for the Central Universities (No. lzujbky-2012-207).