text:
\section{Introduction}

\noindent Phase space methods for wave propagation, such as these presented here, are in the junction between time-frequency analysis \cite{Co1,Co2,Co3,Fl1,Fl2,Hla,Jaf,Jan,Mal}, microlocal analysis \cite{Ben,Ge1,Ge2,Ji1,Spa,Zha}, PDEs \cite{Bal,Bao,Eng,Fan,Fil} and physics \cite{Gro,Hal,Moy,Tat,Wig} (at least). Distributions from Cohen's class (i.e. Wigner Transforms convolved with various smoothing kernels) have been long well known and used in signal processing contexts, as semi-empirical ways to `clean up' the time-frequency structure of various chirps, e.g. \cite{Fl1}. Microlocal measures have evolved in the last twenty years, in close connection with PDEs, in fact with a very similar job in mind, e.g. \cite{Ge1,Ge2}. The need to extend powerful spectral-density-based ideas from stationary to nonstationary stochastic problems gave birth to various nonstationary spectral densities \cite{ Bal,Fan,Fl2,Hal,Ryz}. These three lines of work all converge to the Wigner Transform (WT), and explain the recent explosion in the literature for WT methods for wave propagation (see e.g.\cite{Ath,Bal,Bao,Ben,Eng,Fan,Fil,Ge2,Hal,Ji1,Ji2,Ji3,Lio,Mar,Mau,Ryz,Spa,Zha}). One soon discovers however that, despite the analytic and asymptotic power of these methods, there has been much less work in computing with WTs (for numerical analysis with WTs see \cite{Arn,Bal,Bao}). This is not an accident, and many people complain that it is not at all practical to compute with WTs, see e.g. \cite{Ben,Jaf}. Computations with Wigner measures are feasible, but this is a different thing, and we discuss it below. In this work we propose a way to go around the WTs `bad' features, and exploit its `good' features in the numerical solution of semiclassical wave propagation problems. To do that, we work with a particular distribution from Cohen's class, which we will call the Smoothed Wigner Transform (SWT).

This paper is a continuation of \cite{Ath}. The two papers are complementary in the following way: in \cite{Ath} the results corresponding to individual numerical experiments were presented in some detail. The objective there was to present direct evidense that `the SWT approach works'. However, as is discussed in section \ref{sec4_2}, the efficiency of the method can only be compared to the more conventional approaches when working with a population of numerical experiments for different values of the semiclassical parameter $\varepsilon$. In this work we focus on this synoptic question, and present clear evidence that the SWT approach is indeed more efficient -- in a well defined sense -- than more conventional approaches. \\

Historically, the WT
\[
W^\varepsilon \left[ {f,g} \right]\left( {x,k} \right)=\int\limits_{y\in
\mathbb{R}^n} {e^{-2\pi iky}f\left( {x+\frac{\varepsilon y}{2}} \right)\bar
{g}\left( {x-\frac{\varepsilon y}{2}} \right)dy}
\]
was introduced for the study of the semiclassical limit of quantum mechanics \cite{Wig}, i.e. the high-frequency Schr\"{o}dinger equation. Consider the system of linear equations

\begin{equation}
\label{eq1}
\displaystyle
\begin{array}{l}
\varepsilon \frac{\partial }{\partial t}u^\varepsilon \left( {x,t}
\right)+L^\varepsilon \left( {x,\varepsilon \partial _x }
\right)u^\varepsilon \left( {x,t} \right)=0,
\\
\left( {u_i^\varepsilon \left( {x,0} \right)} \right)_{i=1,...,d} =\left(
{u_{i;0}^\varepsilon \left( x \right)} \right)_{i=1,...,d}
\end{array}
\end{equation}

\noindent for the complex-valued $d$-component wavefield $u^\varepsilon \left(x,t\right)=\left(u_i ^\varepsilon \left(x,t\right) \right)_{i=1,...,d} $, where $x\in \mathbb{R}^n$  and the $d\times d$ matrix-valued Weyl symbol $L^\varepsilon \left( {x,k} \right)$ of the spatial pseudodifferential operator $L^\varepsilon \left( {x,\varepsilon \partial _x } \right)$ is anti-Hermitian, $L^\varepsilon \left( {x,k} \right)=- \overline { L^\varepsilon \left( {x,k}  \right)}^T $. This is a natural way to quantify the requirement that equation (\ref{eq1}) is a wave equation; indeed many interesting linear wave equations have anti-Hermitian Weyl symbols. The use of the Weyl pseudodifferential calculus with WTs is natural, see e.g. \cite{Gro,Poo}. The small parameter $\varepsilon $ scales the typical wavelength of the wavefunctions; the length scale at which the coefficients vary (or equivalently one over the slopes of the coefficients), as well as the typical propagation distance we're interested in are significantly larger than $\varepsilon $. Many semiclassical problems in quantum mechanics and long distance propagation problems in (classical) continuum mechanics can be cast in the form (\ref{eq1}) \cite{Ge2,Ryz,Whi}. For the remainder of this paper when we refer to a high-frequency wave propagation problem we will mean a problem that can be cast in the form (\ref{eq1}) with anti-Hermitian Weyl symbol.

There can be (at least) two problems associated with the IVP (\ref{eq1}): The \textit{high-frequency-limit problem}: recover the limits of certain observables of the wavefield, such as the $L^2$-norm density,
\[
N^0\left( {x,t} \right):=\mathop {\lim }\limits_{\varepsilon \to 0} \left| {\sum\limits_{i=1}^d {u_i^\varepsilon \left( {x,t} \right)} } \right|^2
\]
 and the energy density
\[
E^0\left( {x,t} \right):=\mathop {\lim }\limits_{\varepsilon \to 0} \,\,\left\langle {L^\varepsilon \left( {x,\varepsilon \partial _x } \right)u^\varepsilon \left( {x,t} \right),u^\varepsilon \left( {x,t} \right)} \right\rangle _{\mathbb{R}^d}.
\]
 On the other hand, the \textit{high-frequency problem} consists in solving for the observables at a given nonzero $\varepsilon \ll 1$; naturally, simplifications making use of the fact that $\varepsilon $ is small can be made, but the idea is not to lose all the $\varepsilon $-dependent
information. The difference is nontrivial; for example it is often the case that limit-observables have singularities (become infinite) while for any $\varepsilon >0$ the same observables stay finite for all spacetime points.

WT methods for high-frequency wave propagation problems consist in studying
the Wigner matrix $W^\varepsilon \left({x,k,t} \right) =  \left[ {W^\varepsilon \left[ {u_i^\varepsilon ,u_j^\varepsilon } \right]\left( {x,k,t} \right)} \right]_{i,j=1,...,d} $ instead of the wavefields $\left( {u_i^\varepsilon \left( {x,t} \right)} \right)_{i=1,...,d} $ themselves. It can be said that this is a good idea for two reasons: first of all, due to `time-frequency density' character of the WT wave-like functions are in principle `nicely' represented in phase space. Moreover, the operator for the evolution in phase space is
(asymptotically) simpler, and becomes a Hamiltonian flow in many cases. In
fact this is a general theme in WT-based methods; there are `two
transforms': one of the `data' $ {u^\varepsilon \left( {x,t} \right)}  \mapsto  {W^\varepsilon \left( {x,k,t} \right)} $, and one of the `dynamics' under which the data evolve.

A formidable tool for high-frequency-limit problems consists in Wigner
measures (WM), \[W^0\left( {x,k} \right):=\mathop {w- \lim }\limits_{\varepsilon
\to 0} \left[ {W^\varepsilon \left[ {u_i^\varepsilon ,u_j^\varepsilon }
\right]\left( {x,k} \right)} \right]_{i,j=1,...,d} \] \cite{Ge1,Lio}. Pretty much all bilinear observables of the wavefields can be expressed (in the limit)
in terms of the WM. Recently computational methods for the high-frequency-limit wave propagation based on the WM have been developed \cite{Ji1,Ji2,Ji3}.

Although often discussed in the context of the WT, the question of smoothing phase-space dynamic models and its implications has only rarely been actually studied in a deterministic context - see however \cite{Co1,Mur}.

The remainder of this paper is structured as follows: in section 2 we review the WM approach for the limit problem; SWTs and the phase-space reformulations of problem (\ref{eq1}) in terms of them are reviewed in section 3. Numerical results on the SWT methods for the high-frequency problem are presented in section 4; conclusions and further work are discussed in section 5.

\section{The WM approach for the high-frequency-limit problem}

\noindent For an excellent (defining) overview of the WM approach one should refer to \cite{Ge2} and for a comparison with WKB-based asymptotic techniques to \cite{Spa}. WMs are a refinement of microlocal concepts related to microlocal defect measures, H-measures etc \cite{Ge1,Ge2,Lio}. They provide a description for `asymptotic signals'; properly speaking, an asymptotic signal is a family of functions $\left\{ {u^\varepsilon \left( x \right)} \right\}_{\varepsilon >0} $ that in some sense share the same time-frequency structure. A prototype asymptotic signal is given by the WKB ansatz
\[
\left\{ {u^\varepsilon \left( x \right)=A\left( x \right)e^{\frac{2\pi i}{\varepsilon }S\left( x \right)}} \right\}_{\varepsilon >0} ,
\]
and the corresponding WM is
\[
W^0 \left( {x,k} \right) =  \left|{A\left( x \right)} \right|^2\delta \left( {k-\nabla _x S\left( x \right)} \right) ,
\]
 unfolding amplitude and `local wavenumber' (often referred to as `instantaneous frequency' in time-frequency analysis) information in phase space. Naturally, a unique `local wavenumber' doesn't always exist, as e.g. for multi-WKB signals $u^\varepsilon \left( x \right)=\sum\limits_{n=1}^N {A_n \left( x \right)e^{\frac{2\pi i}{\varepsilon }S_n \left( x \right)}} $. In fact an aspect of the WKB ansatz is that the number of `branches' $N$ at each point $\left( {x,t} \right)$ is a `hidden' parameter. This becomes very important when seeking a solution of (\ref{eq1}) in the form $u^\varepsilon \left( {x,t} \right)=\sum\limits_{n=1}^N {A_n \left( {x,t} \right)e^{\frac{2\pi
i}{\varepsilon }S_n \left( {x,t} \right)}} $, because at different points
of spacetime $\left( {x,t} \right)$ there are different unknown numbers of
(nonzero) branches $N=N\left( {x,t} \right)$. The fact that WMs treat this (i.e. the `number of different local wavenumbers') in a more subtle way is an important advantage in the study of high-frequency limits \cite{Spa}.

For multicomponent wavefields the WM is a positive-definite-Hermitian-matrix-valued measure on phase space \cite{Ge2}. To
link matrix-valued WMs with physical quantities in general one must work in the context of a wave problem like (\ref{eq1}). In that case
\[
N^0\left(
x \right)=\int\limits_{k\in \mathbb{R}^n} {tr\left( {W^0 \left( {x,k}
\right)} \right)dk}
\]
is the total (i.e. taking into account all components
of the wavefield) limit norm density and
\[
E^0\left( x
\right)=\int\limits_{k\in \mathbb{R}^n} {tr\left( {L^0\left( {x,k}
\right)W^0 \left( {x,k} \right)} \right)dk}
\]
the total limit energy density [17,35]. Observe that the WT allows the construction of phase space
densities, `elaborating' on the more traditional physical space
densities, i.e.
\[
N^0\left( {x,k} \right)=tr\left( {W^0 \left( {x,k}
\right)} \right), \,\,\,
 E^0\left( {x,k} \right)=tr\left( {L^0\left( {x,k}
\right)W^0 \left( {x,k} \right)} \right).
\]

The derivation of exact equations for the evolution in time of the WT $W^\varepsilon \left( {x,k,t} \right)$ corresponding to a problem of the form (\ref{eq1}) is a well studied topic \cite{Ge2,Lio,Mar,Mau,Moy,Ryz,Wig,Zha}. From them, asymptotic equations governing the evolution of WMs $W^0 \left( {x,k,t} \right)$ can be obtained. In simple cases they consist in decoupled Liouville equations in phase space for scalar phase-space measures \cite{Ge2,Lio}. Therefore, WT-based numerical methods for the high-frequency-limit problem deal with the numerical solution of the Liouville equation in phase space when the initial condition is an appropriate measure. Moment and level set methods have been formulated and used in this connection \cite{Eng,Ji1,Ji2,Ji3}.

A natural question arises here: can't these methods be used for the small
$\varepsilon >0$ problem? So far the answer seems to be negative. The
$\varepsilon \to 0$ weak limit `kills' a lot of fine scale features. If they
are not killed, but included in the numerical problem, we get an overly
heavy problem, that typically is slower to solve than a full numerical
solution of (\ref{eq1}), thus beating its purpose. The need to put the extra data under control clearly emerges; this is possible in terms of the SWT.

\section{The SWT approach for the high-frequency problem}

\subsection{The WT vs. the SWT}

\noindent The Wigner Transform is defined as a sesquilinear mapping,
\[
 W^\varepsilon :f\left( x \right),g\left( x \right)\mapsto W^\varepsilon \left[ {f,g} \right]\left( {x,k} \right).
\]
When $f=g$ it is called the Wigner distribution of $f$, and denoted as $W^\varepsilon \left[ f \right]\left( {x,k} \right)$. The Wigner distribution of a wavefunction $f\left( x \right)$ is a good way to realize a joint breakdown of the wavefunction's $L^2$-norm over space $x$ and wavenumber $k$ - with the caution that it takes on negative values as well.  The books \cite{Fl1,Fol} are two very important sources on the WT and its properties, \cite{Fl1} emphasizing the signal processing point of view, and \cite{Fol} the relation between the WT and pseudodifferential operators.

Cohen's class of distributions \cite{Co2,Co3,Fl1,Jaf,Jan} is defined as the class of all sesquilinear transforms of the form
\[
C\left[ {f,g} \right]\left( {x,k} \right)= \int\limits_{y,u,z} {f\left( {u+\frac{y}{2}} \right)\bar {g}\left( {u-\frac{y}{2}} \right)e^{-2\pi i\left[ {ky+zx-zu} \right]}\phi _C \left( {z,y} \right)dydudz}
\]
for any distribution $\phi _C \left( {z,y} \right)$ (restricting $f,g$ to test functions). An equivalent definition is the class of transforms that results from convolving the WT with a distributional kernel
\[
K_C \left( {x,k} \right)=\mathcal{F}_{\left( {z,y} \right)\to \left( {x,k}
\right)}^{-1} \left[ {\phi _C \left( {z,y} \right)} \right]
\]
\[
C\left[ {f,g} \right]\left( {x,k} \right)=\int\limits_{{x}',{k}'\in
\mathbb{R}^n} {K_C \left( {x-{x}',k-{k}'} \right)\,\,W\left[ {f,g}
\right]\left( {{x}',{k}'} \right)d{x}'d{k}'}.
\]

Among them there are many attractive alternatives to the WT. When the kernel
is chosen to be the Wigner distribution of some function $h$, $K_C \left( {x,k}
\right)=W\left[ h \right]\left( {x,k} \right)$, and if $f=g$, we get the
spectrogram of $f$ with window $h$. Smoothing WTs with an appropriate kernel
tames the interference terms, but doesn't necessarily kill them completely;
there is a balance between smoothing enough and not smoothing too much. In
this connection we define the scaled SWT as the sesquilinear transform

\begin{equation}
\label{eq2}
\begin{array}{c}
\widetilde{W}^{\sigma _x ,\sigma _k ;\varepsilon }\left[ {f,g} \right]\left(
{x,k} \right) = \\
=\frac{2}{\varepsilon \sigma _x \sigma _k
}\int\limits_{{x}',{k}'\in \mathbb{R}^n} {e^{-\frac{2\pi \left| {x-{x}'}
\right|^2}{\varepsilon \sigma _x^2 }-\frac{2\pi \left| {k-{k}'} \right|^2}{\varepsilon \sigma _k^2 }}W^\varepsilon \left[ {f,g}
\right]\left( {{x}',{k}'} \right)d{x}'d{k}'}.
\end{array}
\end{equation} \\

This is a (scaled) WT convolved with a tensor-product Gaussian with space-domain variance proportional to $\varepsilon \sigma _x^2 $ and wavenumber-domain variance proportional to $\varepsilon \sigma _k^2 $. The scaling is selected with the problem (\ref{eq1}) and the semiclassical regime in mind; see Figures 1,2 for intuition on the scale of smoothing. The SWT provides  coarse-scale\footnote{Different authors favor different terms, e.g. slow-scale, coarse-scale, macroscopic, homogenized etc. Each of these terms stands for a description which contains information for the underlying object, understood to be sufficient in terms of certain qualitative criteria, but which is not complete -- usually it is essentially less than complete, thus leading to compression as well. } description of the time-frequency structure of the underlying signal / wavefield; its resolution, i.e. how much information is effectively discarded, can be tuned through the smoothing parameters.

\begin{figure}
\centering
\includegraphics[width=110mm,height=50mm]{WtTanh_s.eps}
\caption{Wigner Transform of the `tanh chirp', $f^{\varepsilon} (x) = A(x) e^{\frac{2\pi i}{\varepsilon} S(x)}$, $A(x)=e^{-25(x-0.5)^2}$, $S(x)=\frac{-1}{5} log \left( e^{10(x-0.5)}+e^{-10(x-0.5)} \right) $
{\footnotesize (Remark: Image quality has to do with -- the very understanable -- arXiv size limits, and does not represent the quality of the output; see e.g. \cite{Ath}}).
}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=110mm,height=50mm]{SWtTanh_s.eps}
\caption{Smoothed Wigner Transform of the `tanh chirp', correspondig to the scaling $\sigma_x,\sigma_k=O(1)$
{\footnotesize (Remark: Image quality has to do with -- the very understanable -- arXiv size limits, and does not represent the quality of the output; see e.g. \cite{Ath}}).
}
\end{figure}

\subsection{Equations for SWTs}

In \cite{Ath} (an early version of which can be found in the arXiv) we presented the derivation of exact equations for the SWT of a
wavefunction. That result essentially stated that, given the multicomponent
problem (\ref{eq1}), its smoothed Wigner matrix satisfies
\begin{equation}
\label{eq3}
\varepsilon \frac{\partial }{\partial t}\widetilde{W}+A+B^T=0,
\end{equation}
where

\[
A=L^\varepsilon \left( {x-\frac{\varepsilon }{4\pi i}\partial _k
-\frac{\varepsilon \sigma _x^2 }{4\pi }\partial _x ,k+\frac{\varepsilon
}{4\pi i}\partial _x -\frac{\varepsilon \sigma _k^2 }{4\pi }\partial _k }
\right)\widetilde{W},
\]
\[
B=\overline{L^\varepsilon} \left( {x+\frac{\varepsilon }{4\pi i}\partial _k
-\frac{\varepsilon \sigma _x^2 }{4\pi }\partial _x ,k-\frac{\varepsilon
}{4\pi i}\partial _x -\frac{\varepsilon \sigma _k^2 }{4\pi }\partial _k }
\right)\widetilde{W}^T,
\]
\noindent with initial data
\[
 \widetilde{W}\left( {x,k,0} \right)=\left[
{\widetilde{W}^{\sigma _x ,\sigma _k ;\varepsilon }\left[
{u_{i;0}^\varepsilon ,u_{j;0}^\varepsilon } \right]\left( {x,k} \right)}
\right]_{i,j=1,...,d}.
\]

We have used the abbreviation $\widetilde{W}=\left[ {\widetilde{W}^{\sigma _x ,\sigma _k ;\varepsilon }\left[ {u_i^\varepsilon ,u_j^\varepsilon } \right]\left( {x,k,t} \right)} \right]_{i,j=1,...,d} $, and, as we discussed earlier, the Weyl PDO calculus. For concreteness we
will now focus on the (scalar) linear Schr\"{o}dinger equation on $\mathbb{R} $, i.e.

\begin{equation}
\label{eq4}
\begin{array}{c}
\varepsilon \frac{\partial }{\partial t}u^\varepsilon \left( {x,t}
\right)-i\frac{\varepsilon ^2}{2}\Delta u^\varepsilon \left( {x,t}
\right)+iV\left( x \right)u^\varepsilon \left( {x,t} \right)=0,\\
u^\varepsilon \left( {x,0} \right)=u_0^\varepsilon \left( x \right)
\end{array}
\end{equation}

\noindent In that case (\ref{eq3}) becomes
\begin{equation}
\label{eq5}
\begin{array}{c}
 \frac{\partial }{\partial t}\widetilde{W}\left( {x,k,t} \right)+\left(
{2\pi k\frac{\partial }{\partial x}+\frac{\varepsilon \sigma _k^2
}{2}\frac{\partial ^2}{\partial x\partial k}} \right)\widetilde{W}\left(
{x,k,t} \right)+ \\
+\frac{2}{\varepsilon} Re\left(
{iV\left( {x-\frac{\varepsilon }{4\pi i}\partial _k -\frac{\varepsilon
\sigma _x^2 }{4\pi }\partial _x } \right)\widetilde{W}\left( {x,k,t}
\right)} \right)=0, \\ \\
\widetilde{W}\left( {x,k,t} \right)=\widetilde{W}^{\sigma _x ,\sigma _k
;\varepsilon }\left[ {u_0^\varepsilon } \right]\left( {x,k,t} \right). \\
 \end{array}
\end{equation}

This is the exact equation, not an asymptotic approximation. Observe that the smoothing parameters $\sigma _x ,\sigma _k $ are now part of the data of the phase-space equation. `Little' smoothing, i.e. a regularized WT, corresponds to $\sigma _x ,\sigma _k \leqslant 1$ in the scaling (\ref{eq2}). In that case the asymptotic phase space equation is

\begin{equation}
\label{eq6}
\begin{array}{l}
 \frac{\partial }{\partial t}\widetilde{W}\left( {x,k,t} \right)+\left(
{2\pi k\frac{\partial }{\partial x}-\frac{{V}'\left( x \right)}{2\pi
}\frac{\partial }{\partial k}} \right)\widetilde{W}\left( {x,k,t} \right)+
\\
+\varepsilon \left( {{V}''\left( x
\right)\frac{\sigma _x^2 }{8\pi ^2}+\frac{\sigma _k^2 }{2}}
\right)\frac{\partial ^2}{\partial x\partial k}\widetilde{W}\left( {x,k,t}
\right)=O\left( {\varepsilon ^2} \right). \\
 \end{array}
\end{equation}

\noindent In particular, the leading order part is a Liouville equation, and is the same as the corresponding leading-order equation for the WT.

The numerical solution of these equations for the SWT are expected to provide a way for a very efficient coarse-scale solver for wave propagation problems. Indeed, in the more standard approach of solving the original wave equation -- in this case equation (\ref{eq4}) -- we need to keep track of the wavefunction `exactly'. That is, even if a coarse-scale amplitude is all we are looking for, in order to be able to predict how the coarse-scale amplitude will evolve, we need as well to keep track of the fine-scale features of the wavefunction -- up to numerical errors of course. The SWT approach on the other hand, enables us to formulate (exactly or asymptoticly) a numerical problem for a coarse-scale description of the wavefield. Thus, we have to keep track of significantly less information, and ultimately we are able to carry out the computations significantly faster.

\section{Numerical results}
\subsection{The numerical scheme}

\noindent The subject of this paper is the formulation of a numerical method suitable for the simulation of semiclassical problems in terms of the SWT, and the investigation of its efficiency as compared to more conventional methods for the same problem. The test cases we work with is the IVP (\ref{eq4}), up to time $t_{max}$, independent of $\varepsilon$, with initial data scaled with $\varepsilon$, for different initial conditions and potentials. The model we discretize is the Liouville equation for the SWT with $\sigma _x ,\sigma _k \leqslant 1$, i.e. the leading order part of (\ref{eq6}). The conclusions of this investigation, based on the numerical results of sections 4.2 and 4.3 are discussed in section 5. An outline of the algorithm follows.

The WT is computed on a Cartesian grid in phase space with the FFT, with complexity $O\left( {N^2\log N} \right)$ \cite{Mal}. The complexity for the computation of the SWT is $O\left( {L^2N^2\log N} \right)$ where $L$ is the number of sampling points needed for the smoothing kernel. When $\sigma _x ,\sigma _k =O\left( 1 \right)$, $L$ is of $O\left( 1 \right)$. Adaptive computation of the SWT that doesn't spend much time on regions of phase space with no energy is also possible \cite{Ozd}, and might be essential for two- and higher dimensional problems.

The Liouville equation is solved numerically with the use of particles, i.e. the numerical implementation of the method of characteristics (related, but not identical, to \cite{ Arn,Fer,Rav}). An initial population of particles is created, so that the SWT can be interpolated up to an error tolerance from its values on them. The trajectory of each particle is computed according to Hamilton's ODEs with a Runge-Kutta solver; the value of the density on each particle remains unchanged in time. The solution at each moment in time is constructed by interpolating the density from its values on the particles.

\subsection{SWT computations vs. solutions of the original problem (\ref{eq4})}
\label{sec4_2}

The results of the SWT method are compared with exact as well as numerical solutions of the corresponding Schr\"{o}dinger equation (\ref{eq4}). Being able to compute faster with the SWT than by numerical methods for the original problem (\ref{eq4}) (or (\ref{eq1}) in general), is essential to our point of view. So far, the method we present here appears to be significantly faster than WT methods and numerical solutions of (\ref{eq4}); however first we have to make precise what exactly we mean by that.

It is well known that as the semiclassical parameter $\varepsilon$ becomes smaller, the numerical solution of (\ref{eq4}) eventually becomes intractable. Saying that a numerical method is {\em slower} or {\em faster} than another in the $\varepsilon \ll 1$ regime we mean the computation time for the `same' problem grows slower or faster in $\varepsilon$ for the one method than for the other. Hence we're looking for the {\em semiclassical complexity} $T=T\left( \varepsilon \right)$, that is the asymptotic dependence of the total computation time on the semiclassical parameter. In the same way, we're intersted in the rate of growth of the degrees of freedom needed for each numerical method $D=D(\varepsilon)$. That way we can compare numerical methods that are implemented in very different ways, and still get a well defined, meaningful result.

As was shown in \cite{Bao}, meaningful numerical solutions by means of finite differences methods for (\ref{eq4}) can be obtained with two different meshing strategies: a finer meshing strategy is needed for the correct $L^2$ approximation of the wavefunction, while a coarser meshing strategy can lead to correct slow-scale observables (but wrong fast-scale information due to numerical dispersion). As we discussed earlier, the SWT method can only recover slow-scale information. Independent numerical experiments, carried out by Kostas Politis (NTUA), lead to findings consistent with those of \cite{Bao}: the semiclassical complexity of standard finite differences methods appears to be $O \left( \varepsilon ^{-3} \right)$ or more, depending on the meshing strategy.

A wavelet method developed by K. Politis shows a rate of growth of computation time around $O(N^{2.5})$, where $N$ is the number of points (i.e. translated versions of the finest-scale wavelets) used; $N$ is proportional to $\varepsilon^{-1}$. Thus wavelet methods seem very promising for semiclassical computations as well. Wavelet based solutions also appear to be better behaved than finite differences in terms of numerical dispersion, and can be accelerated with thresholding.

In the numerical experiments carried out so far, the SWT approach successfully captures the slow-scale observables of the wavefunction with semiclassical complexity no more than $O \left( -\varepsilon ^{-2} log \left( \varepsilon \right) \right) $; see Figures 3, 4. It is in that sense that our preliminary investigation shows that the SWT method is {\em faster} than full solutions of the problem (\ref{eq4}). (In all examples considered the SWT was significantly faster than the wavelet method, although the rates of growth of computation time in $\varepsilon^{-1}$ are relatively close for the two methods). A more thorough study of this question of semiclassical complexity, including more theoretical aspects as well as different numerical methods (e.g. time-splitting spectral methods, which seem to be particularly well adapted to semiclassical computations) is in order, and in fact in progress, in collaboration with other groups as well.

\begin{figure}
\centering
\includegraphics[width=110mm,height=70mm]{CumTime.eps}
\caption{Total computation time against $\varepsilon^{-1}$. All problems are special cases of (4); Problem 1 corresponds to
$u_0^{\varepsilon} (x) = A(x) e^{\frac{2\pi i}{\varepsilon} S(x)} $,
$ A(x)=\frac{1}{4} \left(  tanh\left( 6.87(x+2.42) \right) +1  \right)$
$ \left(  tanh\left( 6.87(2.42-x) \right) +1  \right)$,
$S(x)=-\frac{x^4}{4}-x^2+2x$; $V(x)=0$;
Problem 2 corresponds to $u_0^{\varepsilon} (x) = A(x) e^{\frac{2\pi i}{\varepsilon} S(x)} $, $A(x)$ as above, $S(x)=-\frac{x^4}{4}+2x$; $V(x)=x$;
Problem 3 corresponds to $u_0^{\varepsilon} (x) = e^{-\frac{1+7i}{10\varepsilon}x^2}
+e^{-\frac{0.2+3i}{10\varepsilon}x^2}
+e^{-\frac{0.9-8i}{10\varepsilon}x^2}$; $V(x)=x$;
Problem 4 corresponds to $u_0^{\varepsilon} (x) = e^{-(1+\frac{3i}{\varepsilon})x^2-2x-4}+e^{-(1+\frac{2i}{\varepsilon})x^2-x-1}
+e^{-(1+\frac{i}{\varepsilon})x^2-\frac{2}{3}x-\frac{4}{9}}$; $V(x)=0$. The SWT solutions are found to be valid slow-scale representations of full numerical solutions of Problems 1,2,3 (courtesy of K. Politis), and of the exact solutions of Problem 4. Quantifying in a satisfactory way how good a `slow-scale representation' is has proved somewhat subtle, and is in progress; see \cite{Ath} for qualitative comparisons. The time for the computation of the SWT of the initial condition, although asymptotically important, is very small in this examples.  Certain normalizations were necessary to fit all the data in the same axes; in any case the slopes are not affected.
}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=110mm,height=70mm]{CumDOFs.eps}
\caption{Number of particles used against $\varepsilon^{-1}$.}
\end{figure}

\subsection{SWT vs. WT computations}

WT methods are slower than either full solutions of (\ref{eq4}) or SWT solutions. As compared with the SWT method we propose, the only difference would be another, more complicated initial condition. Apart from being represented with more points, the initial WT is qualitatively different than the SWT. It features many short waves with important amplitude, being a high-frequency 2-d wavefunction itself.

Our numerical experiments indicate that WT methods are in any case one order of magnitude in $\varepsilon^{-1} $ slower, and sometimes even more. One particular phenomenon we observed in those numerical experiments is that, especially in problems with strong interference terms, in order to recover a meaningful approximation of the evolution of the WT many more points are required than in order to represent it initially. That is, if we use just enough particles to approximate the initial WT well, then often our numerical solution evolves to something very noisy; we need to use many more particles than that. This point is in contrast with the respective behavior for the SWT, and in fact we have no precise understanding at this point of how the number of particles for a given WT should be selected. (When this behavior kicks in, the complexity of the WT solution seems to be more than $O \left( \varepsilon ^{-3} \right)$).

\section{Conclusions -- Further work}
\noindent Semiclassical computations are well known to be a very demanding numerical problem. However, it seems that the more precise study of how different methods, in general solving for different things, behave as the semiclassical parameter goes to zero has started only recently. In this work we carry out a preliminary numerical investigation of a recently formulated SWT-based method for semiclassical simulations \cite{Ath}, and find that it is in general faster than full solutions of the original (physical space) PDE, based on numerical experiments, consistently with analysis of the more traditional methods \cite{Bao}. Also, that -- as expected -- SWT based computations are significantly faster and better behaved than WT computations.

As we mentioned earlier, equation (\ref{eq5}) is in fact a whole family of phase-space models, which can then be asymptotically treated to create a hierarchy of approximations for each phase-space model. In this paper we have only worked in a specific asymptotic regime, the semiclassical regime. One reason for this is that it allows direct comparisons with existing WT-based studies; also the history and evolution of phase-space methods draws heavily on semiclassical problems. However it should be clear that the principal reason for the SWT being `faster' is the fact that it is, and remains as time passes, a much simpler function than the WT. In other words {\em homogenization} here comes from the smoothing, and it is an essentially different thing than {\em asymptotics}.

These insights lead to the following conclusion: numerical methods corresponding to the exact equations should be formulated and compared to full solutions for different problems (e.g. `any' linear problem) and various smoothing choices; we can then look for practical, possibly problem specific, `slow-scale solvers'. Thus, further work at hand includes the treatment of different models included in (\ref{eq5}). Numerical methods for such models are needed; parallel implementation might be crucial for higher-dimensional problems.

\section{Acknowledgements}

The author would like to thank Profs. G. Papanicolaou, I. Daubechies, G. Makrakis, and L. Ryzhik for helpful discussions. Also, K. Politis for helpful discussions and for providing numerical results used for comparison. This work has been partially supported by NSF grant DMS-0530865.

\footnotesize


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

Gravitationally bound objects follow Keplerian trajectories in the classical two-body problem. In an ambient gaseous medium, however, orbits evolve under the friction that the gas exerts on the embedded bodies. Gaseous friction changes the orbital parameters through apsidal precession, orbital migration, and eccentricity change.
We will discuss two mechanisms of gas--object interaction \citep[e.g.][]{Thun2016}: (i) hydrodynamic drag, and (ii) gaseous dynamical friction. Both the hydrodynamic drag force and the dynamical friction force oppose the direction of motion of the orbiting body, transferring energy and angular momentum to the gas. As a consequence, the semi-major axis shrinks and the orbital eccentricity may change as well. The resulting orbital evolution, therefore, depends on the properties of the gaseous and the parameters of the system orbit.

The response of gaseous surroundings to the passage of a gravitating object is described by the Bondi-Hoyle-Lyttleton accretion theory \citep{1939PCPS...35..405H,1944MNRAS.104..273B,1952MNRAS.112..195B,2004NewAR..48..843E}. Further, the overdense wake left by a passing object implies a momentum exchange \citep{1971MNRAS.154..141H,1985MNRAS.217..367S,1994ApJ...427..351R,1994A&AS..106..505R,1995A&AS..113..133R,1996A&A...311..817R,1999A&A...346..861R,2016A&A...589A..10T} much like that in collisionless dynamical friction \citep{Chandrasekhar1943}.
\cite{Ostriker1999} derived gaseous dynamical friction force formulae using time-dependent linear perturbation theory.

More broadly, the theory of gaseous dynamical friction has applications throughout astrophysics, where gas and gravitating bodies frequently intermingle. Some particularly rich applications have included the dynamics of young stars and clusters embedded in molecular clouds \citep[e.g.][]{2006ApJ...638..369K,2008ApJ...679L..33K,2010MNRAS.402.1758S,2011ApJ...735...25N,2014ApJ...794..167S,2019ApJ...876..142K,Antoni2019,2019MNRAS.490.5196C}, the interactions and capture of stars by discs in active galactic nuclei \citep{1991MNRAS.250..505S,1993ApJ...409..592A,Narayan2000,2005ApJ...619...30M,2011ApJ...726...28B,2017MNRAS.464..946S,2016MNRAS.460..240K,2018MNRAS.476.4224P,MacLeod2020,2020ApJ...903..133S,2021arXiv211003741M}, and common envelope phases  in which a star engulfs its companion star or planets within its gaseous envelope \citep{1976IAUS...73...75P} and dynamical friction drives orbital inspiral \citep[e.g.][]{Iben1993,Ivanova2013,2007ApJ...661.1192V,2015ApJ...803...41M,Staff2016,MacLeod2017,2018ApJ...853L...1M,2020ApJ...897..130D,2020ApJ...899...77E}.

Following \cite{Ostriker1999}, subsequent authors examined various aspects of orbital evolution both numerically and analytically in different astrophysical contexts \citep[see][for recent discussion of this literature]{MacLeod2017,Antoni2019}. Of particular relevance to our subsequent discussion is work that considered the effect of dynamical friction on orbits. For example, \cite{Salcedo2001} used numerical simulations to study the orbital evolution of objects in different models of gaseous spheres. \cite{Kim2007} and \cite{2008ApJ...679L..33K} investigated circular orbits in a uniform gaseous medium using semi-analytic and hydrodynamic methods, respectively. \cite{Salcedo2019} examined the evolution of eccentric orbits, comparing analytic and numerical results. \cite{Bonetti2020} considered the implications of prograde and retrograde motion in rotating media. Among the shared questions are the characteristic length scale over which gravitational wakes extend (defining the Coulomb logarithm term) and how accelerated, rather than linear, motion affects the forces.

In this paper, we focus particularly on the eccentricity evolution of two-body systems. We adopt a model in which a less massive secondary object orbits in the gaseous envelope of a more massive, extended primary object. We use both (i) a simple numerical integrator to solve the equations of motion, reconstruct the orbital path, and measure the eccentricity evolution; and (ii) a semi-analytic approach to directly estimate the eccentricity evolution from the orbital parameters and drag forces. In section \ref{sec:2}, we introduce the equations of motion and our numerical technique. Then we outline our semi-analytic formalism in section \ref{sec:3}. We demonstrate the application of our formalism in a toy model first, in section \ref{sec:isothermal}, where the primary object is an isothermal sphere with a power-law density profile. Then we use both numerical and semi-analytic techniques to calculate the eccentricity evolution of a Jupiter-like planet within the stellar envelope of a red giant star in section \ref{sec:4}. We examine a series of different polytropic models and show how the eccentricity evolution depends on the polytropic density profile of the envelope. We summarize our conclusions in section \ref{sec:conclusion}.

\section{Method}
\label{sec:2}

\subsection{Equations of Motion}
We start with a simple model of a two-body system where a more massive object (primary) is represented by an extended, static, spherically symmetric gaseous halo; and a less massive companion (secondary) is orbiting within the gas. In this picture, the acceleration of the secondary originates from both the gravitational force of the enclosed mass ($m_\star$) and the external force ($\vec{F}_\mathrm{ext}$) that the gas exerts on it
\begin{equation}
\label{EoM}
    \vec{a} = -G\dfrac{m_\star}{r^2} \dfrac{\vec{r}}{r} + \dfrac{\vec{F}_\mathrm{ext}}{m} \, .
\end{equation}
Here, $G$ is the gravitational constant, $\vec{r}$ is the position and $m$ is the mass of the secondary. We only examine the orbital evolution of the secondary and neglect its back-reaction on the gas. The external forces are considered to be either the hydrodynamic drag force or the gas dynamical friction. The hydrodynamic drag force \citep{Villaver2009} is
\begin{equation}
\label{hyd}
    \vec{F}_\mathrm{hyd} = - \dfrac{1}{2}C_d R^2 \pi \rho v^2 \dfrac{\vec{v}}{v}
\end{equation}
where $R$ is the radius, $\vec{v}$ is the velocity of the secondary and we use $C_d = 1.0$. $\rho$ is the local density in the gas. The gas dynamical friction force \cite{Ostriker1999} is
\begin{equation}
\label{df}
    \vec{F}_\mathrm{df} = - 4\pi(G m)^2  \dfrac{\rho I}{v^2} \dfrac{\vec{v}}{v}
\end{equation}
where $I$ is the Mach-number dependent parameter i.e.~$I_\mathrm{sub} = 0.5 \ln \left[ (1+v/c_s)/(1-v/c_s) \right] - v/c_s$ in the subsonic regime and $I_\mathrm{sup} = 0.5\ln \left( 1 -(v/c_s)^{-2}\right) - \ln \left( r / R\right)$ in the supersonic regime. $c_s$ is the local sound speed in the gas, $r$ is the characteristic size of the system (here we adopt the separation of the secondary object).

One implicit assumption in these expressions is that force depends on the local conditions (those in the vicinity of the orbiting object). In principle, this might not be satisfied if there are changes in quantities like the density over the characteristic length scale of the Coulomb logarithm \citep[see][for a more detailed description of the applicability and limits of this local formalism as compared to hydrodynamic simulations]{Salcedo2019}.

\subsection{Numerical Integration}
We solve Eq.~\eqref{EoM} with a simple leapfrog algorithm applying either the hydrodynamic drag or the dynamical friction as the source of the external force to reconstruct the orbital path. The $\vec{r}_{i}$ position of the object is stepped by the following algorithm:
\begin{align}
    \vec{v}_{i+1/2} &= \vec{v}_i + \vec{a}_i \dfrac{\Delta t}{2} \\
    \vec{r}_{i+1} &= \vec{r}_i + \vec{v}_{i+1/2} \Delta t \\
    \vec{a}_{i+1} &= -\dfrac{G m_\star}{\vec{r}_{i+1}^2} + \dfrac{\vec{F}_\mathrm{ext}}{m}\\
    \vec{v}_{i+1} &= \vec{v}_{i+1/2} + \vec{a}_{i+1} \dfrac{\Delta t}{2}
\end{align}
where $\vec{a}$ is the acceleration and $\Delta t$ is a fixed time-step of the integration. $m_\star = 4\pi \int_0^r \rho(r) r^2 dr$ is the enclosed mass of the envelope at the $r$ position. We integrate the trajectory until the secondary object either (i) approaches a certain fraction of the initial orbital separation e.g.~$10\%$ of the initial semi-major axis, or (ii) the enters the subsonic regime of the dynamical friction. In the dynamical friction simulations, we focus on the eccentricity evolution in the supersonic regime only. After integrating the entire orbital trajectory (up to the above separation criteria), we measure the eccentricity per orbit as $e_i = (r_{a,i} - r_{p,i})/(r_{a,i} + r_{a,i} )$ where $r_{a,i} = \max(r)_i$ and $r_{p,i} = \min(r)_i$ are the local $i$th maxima and minima of the oscillating separation. Similarly, the local semi-major axis is measured as $a_i = (r_{a,i} + r_{p,i})/2$.

\section{Semi-Analytic Model}
\label{sec:3}

The orbital eccentricity changes due to the external forces that the gas exerts on the secondary as it travels through the medium. The gas--object interaction both leads to the dissipation of the specific orbital energy ($\varepsilon$) and the change in the specific angular momentum $(h)$. We estimate the eccentricity following changes in these values to be
\begin{equation}\label{eq:ecc}
    e \approx e_0 + \Delta e  =  \sqrt{1 + \frac{2(\varepsilon_0 + \Delta \varepsilon) (h_0 + \Delta h)^2}{(GM)^2}} ,
\end{equation}
where $\varepsilon_0 = -\dfrac{GM}{2a_0}$, $h_0 = \sqrt{GMa_0(1-e_0^2)}$ are the initial values of $\varepsilon$ and $h$ with $a$ being the semi-major axis and $e$ being the eccentricity. $M = m_\star + m$ is the total mass.
 To simplify Eq.~\eqref{eq:ecc}, we introduce $u_0 = e_0^2 -1$ and $u = e^2 - 1$ variables and express the change as
\begin{equation}\label{eq:du}
    \dfrac{\Delta u}{u_0} \approx \dfrac{\Delta \varepsilon}{\varepsilon_0} + \dfrac{2 \Delta h}{h_0}
\end{equation}
where we neglect quadratic terms in the energy and angular momentum changes as
\begin{equation}
    \frac{(\varepsilon_0 + \Delta \varepsilon) (h_0 + \Delta h)^2}{(GM)^2} \approx \frac{\varepsilon_0 h_0^2}{(GM)^2}\left[ 1 +  \dfrac{\Delta \varepsilon}{\varepsilon_0} + \dfrac{2 \Delta h}{h_0}
\right] \,.
\end{equation}

The source of the energy change is the dissipative external force
\begin{equation}\label{epsilon}
    \Delta\varepsilon = \frac{\vec{v}\cdot \vec{F}_\mathrm{ext}}{m} \Delta t,
\end{equation}
where $\Delta t$ is the time window in which we estimate the transfer of orbital energy to the gas. Similarly, the loss of specific angular momentum comes from the torque that the gas exerts on the companion for $\Delta t$ time
\begin{equation}\label{ha}
    \Delta h = \frac{|\vec{r}\times \vec{F}_\mathrm{ext}|}{m} \Delta t.
\end{equation}

For the approximation of Eq.~\eqref{eq:ecc} to hold, we need to choose a time interval, $\Delta t$, such that $|\Delta \varepsilon/\varepsilon_0|\ll 1$ and $\Delta h/h_0 \ll 1$. In that case the resulting change in eccentricity will be small and the expansion of Eq~\eqref{eq:ecc} will be justified. In practice, we use the orbit as a representative time unit. For systems with slowly-varying orbital elements due to dynamical friction (equivalently where $|F_{\rm ext}/F_{\rm grav}|\ll 1$) changes over an orbital period will be small enough to apply  Eq~\eqref{eq:ecc}.

Additionally, for an eccentric orbit, the gaseous conditions and relative motion of the object through the gas change as a function of orbital phase. Therefore, applying the local conditions at a single position cannot represent the whole orbit \citep[e.g.][]{Salcedo2019}. Instead, we consider the local conditions at the two extremes of periapse and apoapse. Thus,
\begin{equation}
\Delta\varepsilon \approx \Delta\varepsilon_{\rm p} + \Delta\varepsilon_{\rm a}
\end{equation}
and
\begin{equation}
    \Delta h \approx \Delta h_{\rm p} + \Delta h_{\rm a}
\end{equation} where we introduce the lower index notation for periapsis ($p$) and apoapsis ($a$). To do so, we must approximate the time spent at the orbital locations of periapse, $\Delta t_p$, and apoapse $\Delta t_a$.

Given these approximations, we calculate $e$ of Eq.~\eqref{eq:ecc} using Eq.~\eqref{epsilon} and Eq.~\eqref{ha} for hydrodynamic drag and gas dynamical friction forces, respectively in the sections that follow.

\subsection{Hydrodynamic drag}

In case of a pure hydrodynamic drag, the energy dissipation and the change of angular momentum are
\begin{equation}
    \Delta \varepsilon =  - \dfrac{C_d}{2 m}\pi\rho R^2 v^3 \Delta t \; ,
\end{equation}
and
\begin{equation}
    \Delta h = -\dfrac{C_d}{2 m}\pi\rho R^2 |\vec{r}\times\vec{v}^2| \Delta t \; .
\end{equation}
We evaluate these terms only at the peri- and apoapsis where the position is \begin{equation}
    r_\mathrm{p,a} = a(1\mp e)
\end{equation}
and the velocity is
\begin{equation}
    v_\mathrm{p,a} = \sqrt{\dfrac{GM(1\pm e)}{a(1 \mp e)}}.
\end{equation}
 At these locations in the orbit, the position and velocity are perpendicular. We substitute these expressions into Eq.~\eqref{eq:du}:
\begin{equation}
    \dfrac{\Delta u_{\rm p,a}}{u_0} = \alpha \rho a_0 \left( \dfrac{1\pm e_0}{1\mp e_0} - 1 \right)
\end{equation}
where we introduced $\alpha = \pi R^2 \dfrac{C_d}{m} $ and applied the approximation that $\Delta t_{\rm p,a} \approx a/v_{\rm p,a}$. The eccentricity at the peri- and apoapsis is
\begin{equation}
    e_{\rm p,a} = \sqrt{\alpha \rho a_0 [ 1-e_0^2 - (1 \pm e_0)^2  ] -e_0^2 }.
\end{equation}
 This yields a net eccentricity following a complete orbit of $e \approx e_{\rm p} + e_{\rm a} $.

\subsection{Dynamical friction}

Similarly, we outline the effect of gaseous dynamical friction on the eccentricity evolution. The specific orbital energy and the specific angular momentum changes are
\begin{equation} \label{eq:Deps}
    \Delta \varepsilon =  -4\pi G^2 m\rho I \frac{1}{v}  \Delta t
\end{equation}
\begin{equation} \label{eq:Dh}
    \Delta h = -4\pi G^2 m \rho I \frac{|\vec{r}\times\vec{v}|}{v^3} \Delta t
\end{equation}
We evaluate these terms at the apses again, substitute them to Eq.~\eqref{eq:du}:
\begin{equation}
    \dfrac{\Delta u_{\rm p,a}}{u_0} = \beta \rho a_0^3 \xi \left[ \dfrac{1\mp e_0}{1 \pm e_0} -  \left(\dfrac{1\mp e_0}{1 \pm e_0}\right)^2 \right]
\end{equation}
where we introduced $\beta = 8 \pi m I / M^2$ and $\Delta t_{\rm p,a} = \xi a/v_{\rm p,a}$. We note that $\xi \approx 1$ but can be fit by comparing the prediction with the numerical results. Then
\begin{equation}
    e_{\rm p,a} = \sqrt{\beta \rho a_0^3 \xi \left[ \dfrac{(1\mp e_0)^3}{1 \pm e_0} - (1\mp e_0)^2\right] - e_0^2}
\end{equation}
is the  eccentricity at the peri- and apoapsis. Just as in the case of hydrodynamic drag, the net eccentricity after an orbit can be computed as $e \approx e_{\rm p} + e_{\rm a} $.

\section{Results}
\subsection{Isothermal sphere}
\label{sec:isothermal}

We use the numerical simulations and the above semi-analytic formalism to measure and predict the rate of eccentricity change of a secondary object in an isothermal gas sphere of mass $m_\star$ as a model of an extended primary object. The density profile of the gas is chosen to be a power-law function of the radius as $\rho(r) = \rho_0 r^{-\gamma}$ where
\begin{equation}
    \rho_0 = \frac{(3-\gamma) m_\star}{4 \pi \left(R_\star^{3-\gamma} - R_c^{3-\gamma}\right)} \; ,
\end{equation}
$R_c$ and $R_\star$ are the inner and outer edges of the density profile. We use $R_\star = 1 \, \mathrm{au}$ as the unit of distance and the total enclosed mass at $R_\star$ as the unit of mass $m_\star(R_\star) = 1 M_\odot$.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig1.pdf} \vspace{-15pt}
\caption{\label{fig1_isotermal_175} \textit{Top:} Planar view of the orbital path of an orbiter in a isothermal gas sphere with a power-law density profile with $\rho(r) = \rho_0 r^{-1.75}$. The companion spirals inward due to gas dynamical friction (DF with blue) on the left panel, and due to hydrodynamic drag force (HD with green) on the right panel, respectively, until the orbit reaches $0.1 a_0$ in both cases. \textit{Bottom:} Orbital eccentricity as a function of semi-major axis calculated per orbit. The mass ratio is $0.001$, the initial eccentricity is $0.05$, and the initial semi-major axis is $0.94 R_\star$ in both DF and HD simulations. Under the influence of dynamical friction, the orbit becomes more eccentric, while under the influence of the hydrodynamic drag it circularizes. }
\end{figure}
In Figure \ref{fig1_isotermal_175}, we show the orbital paths of a  simulation with a $\gamma = 1.75$ isothermal density profile on the top panels. Here, we chose a secondary with mass $m=10^{-3} M_\odot$ with initial semi-major axis $a = 0.94 \, \mathrm{au}$ and initial eccentricity $e=0.05$. The radius of the secondary object is chosen to be $R = 2\cdot10^{-3} \, \mathrm{au}$ and sound speed is a constant as $c_s = \sqrt{1.4 k_B T / m_\mathrm{H2}}$ with $T = 10^4 K$, $k_B$ being the Boltzmann-constant, and $m_\mathrm{H2}$ being the mass of molecular hydrogen. These initial parameters set the secondary to orbit within the density profile ($r < R_\star$). The top left and right panels show the simulated inspiral of the secondary object under pure dynamical friction (with blue) or hydrodynamic drag (with green), respectively, in the x-y plane of the motion. Solid circles indicate $R_\star$ (the outer cut-off of the gas) and $0.1 a_0$ (i.e.~$10\%$ of the initial semi-major axis) radii. We ran the simulations between $r = 0.09 - 0.94 \, \mathrm{au}$ radius range. The bottom panel shows the eccentricity change during the inspiral. With $\gamma = 1.75$ density profile, the dynamical friction causes the orbit to become more eccentric while hydrodynamic drag circularizes the orbit.

Using the same mass ratio but $a_0 = 0.5 \, \mathrm{au}$, we predict the change of the eccentricity at only the first orbit for both apses as a function of $\gamma$ and $e_0$. We show the ratio between $\Delta e_{\rm p} = e_{\rm p} - e_0$ and $\Delta e_{\rm a} = e_{\rm a} - e_0$ in Figure \ref{gamma_isothermal} where the colors indicate the initial eccentricities of the systems from $e_0 = 0$ up to $1$.

Because forces at periapse tend to circularize the orbit while those at apoapse tend to make it more eccentric, the ratio $\Delta e_{\rm p} / \Delta e_{\rm a}$ is an important indicator of the orbital evolution. We find that $\Delta e_{\rm p} / \Delta e_{\rm a} < 1$ for $\gamma < 3$ which means that the orbits tend to become more eccentric in shallow gas density profiles. By contrast, orbits circularize in cuspy density profiles for $\gamma > 3$. More eccentric the initial orbits are, the more eccentric they get ($\gamma<3$) or the more they circularize ($\gamma>3$).

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig2.pdf} \vspace{-15pt}
\caption{\label{gamma_isothermal} The ratio between the eccentricity change at peri- and apoapses as a function of the slope ($\gamma$) of the isothermal density profile. The ratio of change of the orbital eccentricity is shown for 100 different realizations of a system. Colors indicate the initial eccentricities of the 100 systems. For all eccentricities, $\gamma=3$ is a critical value that divides greater eccentricity change at periapse and apoapse. The signs of $\Delta e_p$ and $\Delta e_a$ imply that when $\gamma<3$ orbits become more eccentric, while when $\gamma>3$ they the become less.}
\end{figure}

We have compared the above semi-analytic results with the numerical simulations and found that $\xi$ scaling parameter can be fit as $\xi \approx {\frac{\pi}{2}} (1-e_0^2)$ in case of dynamical friction. Note that $\xi$ is a fitting parameter that defines the time window $\Delta t = \xi a/ v$ in which we approximate the apses transit in Eq.~\eqref{epsilon} and Eq.~\eqref{ha}. $\xi$ is measured by comparing the discrepancy between the eccentricity change in the numerical simulations versus the eccentricity change in the semi-analytic model with $\xi = 1$. In Figure \ref{fitting}, we show how we measured the deviations at 9 points of initial eccentricities for 3 models ($\gamma = 0,1.5,3.5$).

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig3.pdf} \vspace{-15pt}
\caption{\label{fitting} The fitting parameter $\xi$ on which we approximate the apses transit ($\Delta t = \xi a / v$) as a function of initial eccentricity for $3$ models of power-law density profiles ($\gamma = 0.0,1.5,3.5$). }
\end{figure}

\subsection{Common Envelope Inspiral}
\label{sec:4}

An ubiquitous feature of simulations of hydrodynamic simulations of common envelope inspiral is the development of moderate eccentricity -- even when these models are initialized in circular orbits. This is seen clearly, for example in Figure 2 of \citet{2012ApJ...746...74R}, Figure 11 of \citet{2012ApJ...744...52P} and Figures 1 of \citet{2016ApJ...816L...9O} and \citet{2019MNRAS.490.3727C}. The analysis of \citet{2019MNRAS.490.3727C} progresses much further and analyzes the drag forces applied to the infalling bodies, and discusses the similarity (in the early inspiral) and departure (after the envelope has been disturbed) of these forces from those predicted by gaseous dynamical friction and an undisturbed profile. Yet, results appear quite varied across the models presented by different groups and initial conditions. For example, the eccentricity seen in the models of \citep{2020A&A...644A..60S} is quite small compared to that of \citet{2016ApJ...816L...9O} when the same code but different initial conditions are employed. To our knowledge, no clear explanation for this eccentricity, or why it differs across simulation models, has been outlined. We explore this question further in this section by examining the eccentricity evolution in a subset of hypothetical model common envelope inspiral episodes.

We simulate the eccentricity evolution of a Jupiter-like planet engulfed by a Sun-like giant star at its late, red giant evolutionary phase. The total mass of the star is $m_\star = 1 M_\odot$ and its radius is $R_\star = 1 \, \mathrm{au}$, the planet has a mass of $m = 10^{-3} M_\odot$ and a radius of $R = 2\cdot10^{-3} \, \mathrm{au}$ (or approximately 4.28 Jupiter radii). The star comprises two structural parts: the inner core which makes up $1/3$ of its total mass and a polytropic envelope which extends up to $R_\star$ and has a polytropic index $n = 1.5$. For the sake of simplicity, we treat the stellar core as a point mass in the center, and construct the envelope solution according to the Lane-Emden equation for $n=1.5$ polytropic index:
\begin{equation}\label{eq:LE}
    \dfrac{1}{r^2} \dfrac{d}{dr} \left( r^2 \dfrac{d\Theta}{d r} \right) = -\Theta^n.
\end{equation}
The unit system of this expression is defined by $4\pi \int_0^1 \Theta^n(r) r^2 dr = 1$. Accordingly $G=M_\mathrm{unit}=R_\mathrm{unit}=1$, which can be rescaled to any physical dimension.
We note that the inner boundary condition of a point-mass core does not effect the orbital evolution of the planet through its inspiral in the stellar envelope because the planet's acceleration only depends on the gravity of the enclosed mass and the local density of the envelope in our simplified model.

In case of dynamical friction however, the orbital evolution also strongly depends on the system's Mach-number. We only examine the first, rapid evolutionary phase of the orbit where the planet orbits in the supersonic regime i.e.~ $v > c_s$. Here, we evaluate the corresponding $I$ parameter of Eq.\eqref{df} locally as $I_\mathrm{sup} = 0.5\ln ( 1 -[v/c_s]^{-2}) - \ln \left( r / R\right)$. The second term is also know as the Coulomb-logarithm in which $r$ is the separation between the planet and the stellar core. In the first term, we approximate the local sound speed as $c_s = c_0 \sqrt{(1+1/n) \Theta(r)}$ where $c_0 = 29.7 km/s$ is a unit conversion parameter because $\Theta(r)$ is the unit-normed solution of equation \eqref{eq:LE}.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig4.pdf} \vspace{-15pt}
\caption{\label{app1} \textit{Top:} Planar view of the orbital path of a companion in a polytropic envelope of a giant star with $n = 1.5$ polytropic index and $R_\star$ radius. The companion spirals inward due to gas dynamical friction (DF with blue) on the left panel, and due to hydrodynamic drag force (HD with green) on the right panel, respectively, until the orbit reaches the subsonic regime of dynamical friction i.e. where $v(r) = c_s(r)$ (for which $r \approx 0.5 R_\star$). \textit{Bottom:} Orbital eccentricity as a function of semi-major axis calculated per orbit. The mass ratio is $0.001$, the initial eccentricity is $0.05$, and the initial semi-major axis is $0.94 R_\star$ in both DF and HD simulations.}
\end{figure}
In Figure \ref{app1}, we show the orbital paths of the planet on the top panels (similarly to Figure \ref{fig1_isotermal_175}). The left panel shows the case in which the planet's orbit evolves under the effect of dynamical friction only. Similarly, the right panel shows the same hypothetical case where the orbit of the planet shrinks under the influence of hydrodynamic drag only. The simulations were started at $a_0 = 0.94 \, \mathrm{au}$ with $e_0 = 0.05$ initial eccentricity. In both cases, the planet migrates inward and its orbit start to circularize. In case of dynamical friction, we observe a transitional phase in which the circularization stops between $a \approx 0.5 - 0.7 \, \mathrm{au}$ where the planet's orbit becomes more eccentric. These changes reflect the changing steepness of the density profile within the polytropic structure -- near the stellar limb, the scale height is small and the density gradient is very steep, leading to orbital circularization. Deeper in the stellar interior, the scale height can be much larger and the orbit becomes more eccentric under the influence of the shallow density profile. Just before entering the subsonic regime (at which point we stop the simulation), the dynamical friction force drops as the Mach-number approaches unity  $F_\mathrm{df} \propto \ln(1-\mathcal{M}^{-2})$. At this limit, the orbit starts to circularize again that continues in the subsonic regime.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig5.pdf} \vspace{-15pt}
\caption{\label{encmass}\textit{Left:} The enclosed mass of the stellar envelope as a function or radius for different polytropic model. \textit{Right:} The $d \ln \rho / d \ln r$ derivative of the density as a function of radius. A horizontal line compares to the critical value of $-d \ln \rho / d \ln r=3$ identified in the context of power-law density profiles. The simulated range of polytropic indexes varies from $n = 0.1$ up to $n = 3.5$ see the colorbar. Note that simulations adopted a point mass-like stellar core with mass $1/3$ of the total mass.}
\end{figure}
After investigating one specific configuration, we extend our analysis to understand how the observed transition point in the eccentricity evolution depends on the density profile of the stellar envelope. We generate $35$ independent density profiles with polytropic indexes ranging from $n = 0.1$ up to $n = 3.5$ evenly. We show the enclosed masses of the different polytropic envelopes with colored curves as a function of radius on the left panel of Figure \ref{encmass}. (In all figures from Figure \ref{encmass}, colors represent the same polytropic indexes of the corresponding stellar envelopes.) We also show the slope i.e.~$-d\ln \rho /d \ln r$ of the radial density profiles of each polytropic model as a function of radius on the left panel of Figure \ref{encmass}. In all models, we use the same point mass-like stellar core with mass $1/3 M_\odot$.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig6.pdf} \vspace{-15pt}
\caption{\label{rvcm} \textit{Left:} The orbital separation (top) and velocity (bottom) of the planet as a function of enclosed mass for a series different of polytropic models indicated by the colormap. \textit{Right:} The local sound speed (top) and Mach number (bottom) as a function of enclosed mass for the same series of polytropic models. Simulations are stopped either when the planet reach $10\%$ of the initial semi-major axis or enter the subsonic regime in the envelope.}\vspace{-5pt}
\end{figure}
We run $35$ simulations that are shown in Figure \ref{rvcm} using the same initial eccentricity ($e_0 = 0.05$) for the planet but different stellar envelope profiles and different initial semi-major axes. We set the initial semi-major axes in a way that the initial periapsis is at the radius corresponding to $0.95 M_\odot$ enclosed mass. Panels on the right show the separation (top) and the velocity (bottom) of the planet as a function of enclosed mass. Simulations are stopped either when the planet reaches $10\%$ of the initial semi-major axis (e.g.~$ n \lesssim 0.5$ or $3.2 \lesssim n$) or enter the subsonic regime (e.g.~$0.5 \lesssim n \lesssim 3.2$) in the envelope.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig7.pdf} \vspace{-15pt}
\caption{\label{eccentricity}  \textit{Top:} Eccentricity as a function of enclosed masses during the inspiral of the planet in the stellar envelope.
Color scale labels the polytropic index of the $35$ different models. Dots show the measured eccentricities per orbit from the numerical simulations while solid curves show the corresponding semi-analytic prediction. \textit{Bottom:} The differential values of the data from numerical simulation from the top panel i.e.~the eccentricity change between two consecutive orbits in the given polytrope. The horizontal black line separates regimes in which the orbit circularize (below zero) or becomes more eccentric (above zero).}
\end{figure}
Figure \ref{eccentricity} shows our results regarding the eccentricity evolution of a planet's orbit during its inspiral in various models of a stellar envelope. The top panel shows the eccentricity as a function of enclosed mass in the numerical simulation (dots) and in the semi-analytic estimate (solid curves). Here, we calculate $e_{i+1}$ by taking the local $(\varepsilon_0,h_0)$ values of the numerical simulation in Eq.~\eqref{eq:ecc}. Similarly we calculate $\Delta \varepsilon$ of Eq.~\eqref{eq:Deps} and $\Delta h$ of Eq.~\eqref{eq:Dh} under the local conditions in each model. The bottom panel shows the differential values from the top panel i.e.~$\Delta e_i = e_{i+1} - e_i$ of the numerical simulation. The horizontal line indicates $\Delta e = 0$ below which orbits tend to circularizes and above which orbits become more eccentric. Our results show that after initial circularization the orbital eccentricity start to increase in various models before either approaching $0.1 a_0$ or entering the subsonic regime. This is most prominent for polytropes with $n \lesssim 2$ or $n \gtrsim 3$ which simulations never enters the subsonic regime but appears in all models.

In our model polytropes, we measure the first transition points where the initial circularization turns into eccentricity growth as a function of semi-major axis, Figure \ref{app3}. We observe that our numerical integrations and semi-analytic predictions are in very close agreement, and trace a relationship that defines the fractional critical radius as a function of the envelope's polytropic index. We note, however, that polytropes with a different core mass fraction would lead to a different result. We note that our result can be fit with the linear relation $n(a/R_\star) = -4.53 a/R_\star + 4.54$.  Finally, we compare this relationship to a criterion based on the local density slope, $d\ln \rho / d\ln r = -3$. We find that this approximation, while not as accurate as the full semi-analytic theory (because it neglects the dependence of the dynamical friction drag force on the local Mach number), provides useful context for the resulting eccentricity evolution.

Given these findings, we argue that the ubiquitous presence of eccentricity in simulated common envelope inspirals can be traced to gaseous dynamical friction on the extended envelope. In the early inspiral, the steep gradient of the outer envelope ($-d\ln \rho / d\ln r  \gtrsim3$) damps any orbital eccentricity,  while in the later inspiral, eccentricity can be enhanced by interaction with the comparatively homogeneous envelope interior ($-d\ln \rho / d\ln r \lesssim 3$).

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig8.pdf} \vspace{-15pt}
\caption{\label{app3}The transition points from circularization to eccentricity growth in different polytropes as function of semi-major axis as a fraction of the stellar radius. Both numerical results (with red $\times$-s) and semi-analytic predictions (with blue dots) are shown together with a linear fit (black solid line) to the numerical results. The orange curve shows the transition points corresponding to $d \ln \rho / d \ln r = 3$ in Figure \ref{encmass}.}
\end{figure}

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

In this paper, we examined the eccentricity evolution of an object orbiting in an extended gaseous medium due to hydrodynamic drag or gas dynamical friction. We built a simple numerical integrator to calculate the orbital path of the companion under the influence of the gravity of the primary and the frictional forces exerted by the gas. We measured the eccentricity evolution during the inspiral of the companion and compared the results with the prediction of a semi-analytic approach.
Some key findings of our study are:
\begin{enumerate}
    \item Drag forces applied to the system at periapse tend to make orbits more circular, while those applied at apoapse tend to make orbits more eccentric.
    \item In all centrally-concentrated mass distributions, hydrodynamic drag causes orbital eccentricity to decrease because the drag increases with increasing velocity (e.g. at periapse in an eccentric orbit), see equation \eqref{hyd} and Figure \ref{fig1_isotermal_175}.
    \item Because the gaseous dynamical friction drag force decreases with increasing velocity in the supersonic regime, equation \eqref{df}, whether orbits become more or less eccentric under the influence of gaseous dynamical friction depends on the density profile.
    \item We find that the critical value for a radial power-law density distribution $\rho \propto r^{-\gamma}$ is $\gamma=3$, where lower values of $\gamma<3$ drive orbital eccentricity increase while higher values of $\gamma>3$ drive orbital eccentricity decrease (Figure \ref{gamma_isothermal}).
    \item Both the sign and rate of change of orbital eccentricity can be be accurately predicted by our semi-analytic theory of section \ref{sec:3} when coupled with the numerical coefficient $\xi \approx \pi/2 (1-e^2)$, as reported in Figure \ref{fitting}.
\end{enumerate}
We apply this theory of eccentricity evolution under the influence of gaseous dynamical friction to the orbital evolution of engulfed objects in common envelope phases. Here we consider an example of a Jupiter-like planet interacting with the envelope of a Sun-like star at its late, red giant evolutionary phase. The hydrostatic mass distribution of the stellar envelope has a steep density gradient (lower temperature and smaller scale height) near the surface and a shallower density gradient (higher temperature and larger scale height) in the deep interior. We show that as orbiting objects pass through these mass distributions, the experience circularization in the outer envelope, and eccentricity excitation in the inner envelope. For polytropes of varying index, $n$, we demonstrate the eccentricity evolution and the inflection between eccentricity decrease in the outer envelope and increase in the inner envelope (Figure \ref{eccentricity}). The inflection between eccentricity decrease and growth can be modeled accurately with our semi-analytic model, or approximately by finding the radius within the stellar model where $d \ln \rho / d\ln r = -3$ (Figure \ref{app3}).

Conceptually, our results provide a framework for understanding the evolution of eccentricity in objects being dragged inward in gaseous distributions. In particular, we demonstrate that the development of orbital eccentricity in global hydrodynamic simulations of common envelope phases is indeed realistic, rather than being an artifact of numerics or initial conditions. Our semi-analytic model adopts the Mach-number dependent coefficients of \citet{Ostriker1999}, but could equally be extended to coefficients of dynamical friction that depend on the local density gradient or other properties \citep[e.g.][]{MacLeod2017,2020ApJ...897..130D}. Our results are further suggestive that the emergence of objects from common envelope phases with moderate, non-zero eccentricities may be a natural consequence of the physics of gaseous dynamical friction. There are many other applications of our results, including the formation of Thorne-Zytkow objects \citep[][]{1975ApJ...199L..19T} and the migration of stars in accretion flows around black holes.

\vspace{-15pt}
\section*{Acknowledgements}

The project was sponsored by a Fulbright Student Grant through the Hungarian - American Commission for Educational Exchange. This work was supported by the National Science Foundation under Grant No. 1909203 and by the Black Hole Initiative - which is supported by grants from JTF and GBMF.

\section*{Data Availability}

The python scripts used to create and analyze our simulations and models are available from the corresponding author upon request and publicly available from http://galnuc.elte.hu/Eccentricity\_Evolution\_in\_Gaseous\_Dynamical\newline \_Friction.zip

\vspace{-15pt}
\bibliographystyle{mnras}

redpajama_set_name: RedPajamaArXiv