format: text

text:
# Parametric Dwarf Spheroidal Tidal Interaction

## 1 Introduction

The Local Group dwarf spheroidal (dS) galaxies may be sensitive probes of the chemical and dynamical formation history of the Milky Way (*cf* Harbeck et al. 2001). They have also been used to constrain dark matter models and early structure formation scenarios in the Universe (*cf* Łokas 2002).

Controversy over whether the dS are equilibrium dynamical systems (*cf* Mateo 1998) or tidally driven non-equilibrium objects (Kuhn and Miller 1989, Kroupa 1997) persists. The primary observational evidence for tidally dominated dynamics comes not only from the dS high velocity dispersions and highly elliptical shapes but also from measurements of a spatially variable velocity dispersion (Kleyna et al. 2002) and significant stellar populations beyond the dS nominal tidal radii (Kuhn et al. 1996, Martínez-Delgado et al. 2001). While there is undisputed evidence that the dS in Sagitarius is tidally disrupted (Ibata et al. 1994) some practitioneers have devised more complex multiparametric models which could formally account for the kinematic data in some of the high M/L dwarfs using dark matter (Kleyna et al. 2002, Wilkinson et al. 2002, Lokas 2002).

A theme of the tidal excitation argument presented here is that the gravitational interaction scale for transfering energy to the dS stellar system from the dS-MW orbit can be much larger than the dS classical tidal radius. In those systems where the stellar crossing time can be comparable to the “period” of the external tide, the simple static arguments for equilibrium break down. Early calculations (Kuhn and Miller 1989 – henceforth KM) demonstrated this possibility using particle-mesh gravitational N-body simulations in idealized circular dS orbits.

This paper develops an analytic formulation of the dS-MW tidal interaction problem. Once again, by treating the response of the dS in terms of galaxy oscillations, but now incorporating the tidal interaction through the Mathieu equation, it is possible to understand the evolution of these tidally interacting systems in non-circular orbits far from the simple resonance condition. As a parametric oscillator we find that a dS stellar system can be “inflated” even when its characteristic pulsational spectrum is not tuned to its MW orbital dynamics.

## 2 Parametric tidal interactions

An oscillator driven by a periodic forcing function can sometimes be described by the Mathieu equation

$$\ddot{x}+\omega _{0}^{}{}_{}{}^{2}\left(1\epsilon \mathrm{cos}2\omega t\right)x=0$$
(1)
Here $`x(t)`$ describes the oscillator amplitude function, $`\omega _0`$ is the oscillator resonant frequency, $`\omega `$ is the driving frequency and $`\epsilon `$ is a small constant. Solutions to equation (1) are most easily developed by a perturbative multi-timescale analysis (Bender and Orzag 1978). One finds that, in general, growing (unstable) oscillatory solutions exist for $`\omega =\omega _0/n`$, $`n`$.

An important property of equation (1) is that even for a weakly damped (high ‘Q’) oscillator a significant, non-zero, frequency range for the driving force leads to growing unstable solutions. For small intrinsic damping this result is independent of the frequency width of the natural resonance. For example, near the first resonance $`\omega `$ between $`\omega _0\left(1\epsilon /4\right)`$ and $`\omega _0\left(1+\epsilon /4\right)`$ leads to instability. The second resonance near $`\omega =\omega _0/2`$ is unstable for a narrower frequency range

$$\frac{1}{24}\epsilon ^2\omega >\omega \omega _0/2>\frac{5}{24}\epsilon ^2\omega $$
(2)
We show here that the internal dynamics of stellar systems in the MW’s gravitational environment can be described as a two-dimensional coupled parametric oscillator. These equations can be solved analytically using multi-scale perturbation methods to explain the dynamical and morphological properties of the Local Group dS.

### 2.1 Tide geometry

Each of the dS orbit the MW at a distance which is long compared to the galactic disk scale length so that it is reasonable to treat the MW as a purely spherical central force problem. We define the MW central force

$$\stackrel{}{F}_{\mathrm{MW}}=\stackrel{}{F}(\stackrel{}{r})=F(r)\widehat{r}$$
(3)
Letting $`\stackrel{}{r}`$ describe the displacement between the MW and dS centers and $`\stackrel{}{\delta r}`$ be the displacement to a dS star from the dS center we write the MW tidal force on the star at lowest order as

$`\stackrel{}{\mathrm{\Delta }F}_{\mathrm{MW}}`$ $`=`$ $`\stackrel{}{F}\left(\stackrel{}{r}+\stackrel{}{\delta r}\right)\stackrel{}{F}\left(\stackrel{}{r}\right)`$ (4)
$`=`$ $`\left[F\left(\stackrel{}{r}+\stackrel{}{\delta r}\right)F\left(\stackrel{}{r}\right)\right]\widehat{r}`$
$`+\mathrm{sin}\left(\mathrm{\Delta }\phi \right)F\left(\stackrel{}{r}\right)\widehat{\phi }`$
Here $`\mathrm{\Delta }\phi `$ is the central angle between a dS star, the MW, and the dS center and $`\widehat{r}`$ and $`\widehat{\phi }`$ are unit vectors in the radial and central angle directions (Figure 1). Since the spatial extent of the dS is small compared to its distance from the galactic center we expand equation (4) using

$`F\left(\stackrel{}{r}+\stackrel{}{\delta r}\right)F\left(\stackrel{}{r}\right)`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}(r)\stackrel{}{\delta r}\widehat{r}`$ (5a)
$`\mathrm{sin}\left(\mathrm{\Delta }\phi \right)`$ $`=`$ $`{\displaystyle \frac{1}{r}}\stackrel{}{\delta r}\widehat{\phi }`$ (5b)
We define a cartesian ($`xy`$) coordinate system centered on the dS which contains the MW and the star so that we can express $`\stackrel{}{\delta r}`$, $`\widehat{r}`$ and $`\widehat{\phi }`$ in term of $`\widehat{x}`$ and $`\widehat{y}`$.

$`\stackrel{}{\delta r}`$ $`=`$ $`x\widehat{x}+y\widehat{y}`$ (6a)
$`\widehat{r}`$ $`=`$ $`\mathrm{cos}\phi \widehat{x}+\mathrm{sin}\phi \widehat{y}`$ (6b)
$`\widehat{\phi }`$ $`=`$ $`\mathrm{sin}\phi \widehat{x}+\mathrm{cos}\phi \widehat{y}`$ (6c)
The MW orbital time dependence is contained in $`r`$ and $`\phi `$. To this order of approximation these equations also describe the three-dimensional stellar problem where $`(x,y)`$ are the projected coordinates in the dS-MW orbital plane. The tidal force perpendicular to the orbital plane yields harmonic motion which is not coupled to the MW orbital motion and is ignorable. Thus the tidal force becomes,

$`\mathrm{\Delta }F_x`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}(r)\left(x\mathrm{cos}^2\phi +{\displaystyle \frac{y}{2}}\mathrm{sin}2\phi \right)+{\displaystyle \frac{F(r)}{r}}\left(x\mathrm{sin}^2\phi {\displaystyle \frac{y}{2}}\mathrm{sin}2\phi \right)`$ (7a)
$`\mathrm{\Delta }F_y`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}(r)\left(y\mathrm{sin}^2\phi +{\displaystyle \frac{x}{2}}\mathrm{sin}2\phi \right)+{\displaystyle \frac{F(r)}{r}}\left(y\mathrm{cos}^2\phi {\displaystyle \frac{x}{2}}\mathrm{sin}2\phi \right)`$ (7b)
The form of the MW potential determines $`\frac{\mathrm{d}F}{\mathrm{d}r}`$ and $`F`$. For a Keplerian potential, it is $`\frac{\mathrm{d}F}{\mathrm{d}r}=2\frac{F(r)}{r}`$ and for the logarithmic case it is $`\frac{\mathrm{d}F}{\mathrm{d}r}=\frac{F(r)}{r}`$. Hence we can write a simplified form of the tide for these two potential forms:

$`\stackrel{}{\mathrm{\Delta }F}_{\mathrm{Kep}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}[(x\mathrm{cos}2\phi {\displaystyle \frac{x}{3}}+y\mathrm{sin}2\phi )\widehat{x}`$
$`+(y\mathrm{cos}2\phi {\displaystyle \frac{y}{3}}+x\mathrm{sin}2\phi )\widehat{y}]`$
$`\stackrel{}{\mathrm{\Delta }F}_{\mathrm{log}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}[(x\mathrm{cos}2\phi +y\mathrm{sin}2\phi )\widehat{x}`$
$`+(y\mathrm{cos}2\phi +x\mathrm{sin}2\phi )\widehat{y}]`$
The numerical factors are of course different for different potential assumptions but the form of the resonant solution remains unchanged. It seems likely that the MW potential is actually logarithmic in the domain of most dS so that we assume the form

$$\varphi =v_{c}^{}{}_{}{}^{2}\mathrm{ln}\left(\frac{r}{r_0}\right)$$
(9)
### 2.2 Solution of the coupled Mathieu equations at first order

#### 2.2.1 Equation of motion for a star near the core

In the dS reference frame a star experiences forces due to the MW tide and the dS potential. The dS potential generates a radial acceleration, $`a_r=GM(r)/r^2`$, where $`M(r)`$ is the mass enclosed within a spherical volume of radius $`r`$ centered on the dS. Near the center of the dS we approximate $`M(r)\rho _04\pi r^3/3`$ so that we can write $`a_r\omega _{0}^{}{}_{}{}^{2}r`$ with $`\omega _{0}^{}{}_{}{}^{2}=4\pi G\rho _0/3`$. In this case the total dS-frame acceleration of a star can be written as

$$\frac{\mathrm{d}^2\stackrel{}{r}}{\mathrm{d}t^2}=\omega _{0}^{}{}_{}{}^{2}\stackrel{}{r}+\frac{\stackrel{}{\mathrm{\Delta }F}}{m}$$
(10)
With additional definitions,

$`k`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}r}}(r)`$ (11a)
$`\epsilon `$ $`=`$ $`{\displaystyle \frac{k}{\omega _{0}^{}{}_{}{}^{2}}}={\displaystyle \frac{\omega _{c}^{}{}_{}{}^{2}}{\omega _{0}^{}{}_{}{}^{2}}}`$ (11b)
where the circular frequency of the dS around the Milky Way is $`\omega _c=v_c/r`$. Combining these terms we finally obtain a recognizable set of coupled equations that describe the behavior of a star in the center of a dS as projected in the $`(xy)`$ plane.

$`\ddot{x}+\omega _{0}^{}{}_{}{}^{2}\left(1\epsilon \mathrm{cos}2\phi \right)x`$ $`=`$ $`\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{sin}2\phi y`$ (12a)
$`\ddot{y}+\omega _{0}^{}{}_{}{}^{2}\left(1+\epsilon \mathrm{cos}2\phi \right)y`$ $`=`$ $`\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{sin}2\phi x`$ (12b)
Notice that the coupling between $`x`$ and $`y`$ motion results from the anisotropic form of the galactic tide where the expansive and compressive tidal contributions are both important in determining the dS dynamical response.

#### 2.2.2 Galaxy oscillations

Individual stars also experience forces due to the dS gravitational perturbations but equations (12) neglect temporal variations (due to $`\rho (t)`$) in the dS potential. Thus there is an additional acceleration term which should be included here at late times in the dS-tide evolution when the temporal variability of the dS density is important.

The resulting oscillations in the mean stellar orbits are not included above but can be described in a manner which is analogous to polytropic stellar normal mode oscillations. This description of galaxy dynamics in terms of linear oscillations was the basis of KMs prediction that Local Group dS could be tidally inflated even when the stars were contained within the static tidal radius. Miller and Smith (1994) and VanderVoort (1999) have now demonstrated that weakly damped or even growing finite amplitude normal mode galaxy oscillations are a physically interesting feature of realistic galaxy dynamics simulations.

The lowest order oscillation mode should have no radial nodes – as is the case for radial oscillations in stellar polytropes (*cf* Cox 1980). VanderVoort (1999) showed that the total Lagrangian radial displacement summed over all stars in an N-body system exhibits a very similar oscillatory behavior. This is not entirely intuitive given that the collisionless orbital motion of stars in a galaxy occurs over the length scale of the galaxy, while a mass element in a stellar polytrope oscillates only over the small scale of the oscillation displacement amplitude. Nevertheless VanderVoort proves that the fundamental N-body radial eigenmode for the mean Lagrangian radial particle displacement, $`\stackrel{}{\xi }`$, satisfies

$$\frac{\mathrm{d}^2\stackrel{}{\xi }}{\mathrm{d}t^2}=\omega _{0}^{}{}_{}{}^{2}\stackrel{}{\xi }$$
(13)
with $`\omega _{0}^{}{}_{}{}^{2}=W/I=GM(r)\rho (r)rdr/\rho (r)r^4dr`$ and $`\stackrel{}{\xi }\stackrel{}{r}`$. Here $`W`$ and $`I`$ are the potential energy and moment of inertia, while $`\rho (r)`$ and $`M(r)`$ are the galaxy density and enclosed mass distribution and $`G`$ is the gravitational constant. Realistic density distributions lead to frequencies that are of order $`\pi \sqrt{G\rho _0}`$ where $`\rho _0`$ is the central density.

#### 2.2.3 Circular orbit

If we launch the dS on a circular orbit around the MW, then $`\phi =\omega t`$ and equations (12) describe a set of two-dimensional Mathieu equations. We now show that these equations exhibit resonant instability around $`\omega =\omega _0`$. To obtain a valid solution near the resonance, we use a multi-timescale analysis in terms of $`t`$ and $`\tau =\epsilon t`$ and develop the relevant variables in terms of the perturbation parameter $`\epsilon `$.

$`\omega _0`$ $`=`$ $`\omega +\epsilon \omega _1`$ (14a)
$`x(t,\tau )`$ $`=`$ $`x_0(t,\tau )+\epsilon x_1(t)`$ (14b)
$`y(t,\tau )`$ $`=`$ $`y_0(t,\tau )+\epsilon y_1(t)`$ (14c)
The equations obtained from the lowest order ($`\epsilon ^0`$) for $`x_0`$ and $`y_0`$ are simple harmonic oscillators with frequency $`\omega _0`$ as we would expect in a harmonic system. From the lowest order $`x_0`$ and $`y_0`$ solutions and first order equation we can construct a refined solution. Using complex notation we write

$`x_0(t,\tau )`$ $`=`$ $`A(\tau )\text{e}^{\text{i}\omega t}+\text{c.c.}`$ (15a)
$`y_0(t,\tau )`$ $`=`$ $`B(\tau )\text{e}^{\text{i}\omega t}+\text{c.c.}`$ (15b)
Gathering the terms of order $`\epsilon ^1`$, we obtain differential equations for $`x_1`$ and $`y_1`$

$`{\displaystyle \frac{^2x_1}{t^2}}+\omega ^2x_1`$ $`=`$ $`2\omega \omega _1x_0+\omega ^2\mathrm{cos}\left(2\omega t\right)x_0+\omega ^2\mathrm{sin}\left(2\omega t\right)y_02{\displaystyle \frac{^2x_0}{t\tau }}`$ (16a)
$`{\displaystyle \frac{^2y_1}{t^2}}+\omega ^2y_1`$ $`=`$ $`2\omega \omega _1y_0\omega ^2\mathrm{cos}\left(2\omega t\right)y_0+\omega ^2\mathrm{sin}\left(2\omega t\right)x_02{\displaystyle \frac{^2y_0}{t\tau }}`$ (16b)
These are combined with the lowest order solutions for $`x_0`$ and $`y_0`$. We construct a solution for $`A(\tau )`$ and $`B(\tau )`$ using the well-known ansatz (Fredholm alternative) of forcing all harmonic terms which lead to unstable (secular) solutions to vanish. Obtaining a solution with unbounded coefficients implies that there is no stable oscillatory solution for $`x`$ and $`y`$. Thus we force the coefficients of $`\text{e}^{\text{i}\omega t}`$ and $`\text{e}^{\text{i}\omega t}`$ on the RHS of equations (16a) and (16b) to vanish to obtain two differential equations (and their conjugates) for the functions $`A(\tau )`$ and $`B(\tau )`$:

$`2\omega _1A(\tau )+{\displaystyle \frac{\omega }{2}}A^{}(\tau )+{\displaystyle \frac{\omega }{2\text{i}}}B^{}(\tau )2\text{i}{\displaystyle \frac{\mathrm{d}A(\tau )}{\mathrm{d}\tau }}`$ $`=`$ $`0`$ (17a)
$`2\omega _1B(\tau ){\displaystyle \frac{\omega }{2}}B^{}(\tau )+{\displaystyle \frac{\omega }{2\text{i}}}A^{}(\tau )2\text{i}{\displaystyle \frac{\mathrm{d}B(\tau )}{\mathrm{d}\tau }}`$ $`=`$ $`0`$ (17b)
We solve these four coupled differential equations in terms of the real and imaginary parts of $`A`$ and $`B`$, i.e. $`(\mathrm{}\left(A\right),\mathrm{}\left(A\right),\mathrm{}\left(B\right),\mathrm{}\left(B\right))`$. This yields the matrix equation

$$\frac{\mathrm{d}X}{\mathrm{d}\tau }=MX$$
(18)
where,

$$X=\left(\begin{array}{c}\mathrm{}\left(A\right)\\ \mathrm{}\left(A\right)\\ \mathrm{}\left(B\right)\\ \mathrm{}\left(B\right)\end{array}\right)$$
(19)
$$M=\frac{1}{4}\left(\begin{array}{cccc}0& \left(4\omega _1+\omega \right)& \omega & 0\\ \left(4\omega _1\omega \right)& 0& 0& \omega \\ \omega & 0& 0& \left(4\omega _1\omega \right)\\ 0& \omega & \left(4\omega _1+\omega \right)& 0\end{array}\right)$$
(20)
This equation is solved most readily by finding the eigenvectors of $`M`$. Taking

$$H=\sqrt{\frac{\omega +2\omega _1}{\omega 2\omega _1}}$$
(21)
the eigenvectors of M of respective eigenvalues $`\frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}`$, $`\frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}`$, $`\text{i}\omega _1`$ and $`\text{i}\omega _1`$ are

$$\left(\begin{array}{c}H\\ 1\\ 1\\ H\end{array}\right),\left(\begin{array}{c}H\\ 1\\ 1\\ H\end{array}\right),\left(\begin{array}{c}1\\ \text{i}\\ \text{i}\\ 1\end{array}\right),\text{and}\left(\begin{array}{c}1\\ \text{i}\\ \text{i}\\ 1\end{array}\right)$$
(22)
The solution to equation (18) is obtained in terms of four complex constants of integration $`C=(\alpha ,\beta ,\gamma ,\delta )`$ and the eigenvectors as:

$`\alpha \mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\tau \right]`$ $`=`$ $`H\mathrm{}\left(A\right)+\mathrm{}\left(A\right)+\mathrm{}\left(B\right)+H\mathrm{}\left(B\right)`$ (23a)
$`\beta \mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\tau \right]`$ $`=`$ $`H\mathrm{}\left(A\right)+\mathrm{}\left(A\right)+\mathrm{}\left(B\right)H\mathrm{}\left(B\right)`$ (23b)
$`\gamma \text{e}^{\text{i}\omega _1\tau }`$ $`=`$ $`\mathrm{}\left(A\right)\text{i}\mathrm{}\left(A\right)+\text{i}\mathrm{}\left(B\right)+\mathrm{}\left(B\right)`$ (23c)
$`\delta \text{e}^{\text{i}\omega _1\tau }`$ $`=`$ $`\mathrm{}\left(A\right)+\text{i}\mathrm{}\left(A\right)\text{i}\mathrm{}\left(B\right)+\mathrm{}\left(B\right)`$ (23d)
At late times the growing exponential term dominates, and these equations simplify to yield expressions for $`A`$ and $`B`$,

$`A(\tau )`$ $`=`$ $`{\displaystyle \frac{\alpha }{4H}}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\epsilon t\right]\left(1+\text{i}H\right)`$ (24a)
$`B(\tau )`$ $`=`$ $`{\displaystyle \frac{\alpha }{4H}}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\epsilon t\right]\left(H+\text{i}\right)`$ (24b)
Evidently $`A`$ and $`B`$ define the mean position of a star in the dS.

$`x(t)`$ $`=`$ $`\alpha ^{}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\epsilon t\right]\mathrm{cos}\left(\omega t+\theta \right)`$ (25a)
$`y(t)`$ $`=`$ $`\alpha ^{}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\omega ^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\epsilon t\right]\mathrm{sin}\left(\omega t+\theta \right)`$ (25b)
Equations (25a) and (25b) show that instability occurs if $`\omega >\left|2\omega _1\right|`$. Thus if the circular frequency of the orbit satisfies

$$\left|\omega _0\omega _c\right|<\omega _c\epsilon /2$$
(26)
the dS will be tidally excited. The constants $`C`$ determine the initial position and velocity of a given star in the dS. Notice that the growing eigenmode solution involves stars oriented along a line that makes an angle $`\theta `$ with respect to the $`x`$ axis where $`\mathrm{tan}\theta =H.`$ Thus, we find in general that stars expand away from the center of the dS depending on the difference between the dS orbital frequency and $`\omega _0`$, and on the perpendicular distance a star makes to a line in the $`(xy)`$ plane of the dS at angle $`\theta `$ with respect to the line-of-sight toward the MW.

Extending this calculation to second order in $`\epsilon `$ reveals that, unlike the simple one-dimensional case, the next resonance near $`\omega _0/2`$ is stable with no secular amplitude variation.

#### 2.2.4 Elliptical orbit

Equations (25a) and (25b) assume a circular orbit. Their applicability is also restricted, depending upon the central concentration and oscillation frequency $`\omega _0`$ of the dS. For example, the perturbation parameter, $`\epsilon =\omega _{c}^{}{}_{}{}^{2}/\omega _{0}^{}{}_{}{}^{2}`$, can approach unity for a circular orbit in a logarithmic potential. Under these conditions the static tide force alone is strong enough to disrupt the dS. Since the dS are generally not in circular orbits it is important to extend the model to elliptical orbits. In a logarithmic MW potential we can recover the form of the Mathieu equation with the epicycle approximation. We define an “eccentricity” $`e`$ from $`r(t)`$ where $`\kappa =\sqrt{2}\omega `$ is the epicycle frequency for a logarithmic potential. It follows that (Binney and Tremaine 1987),

$`r(t)`$ $`=`$ $`r_c\left(1e\mathrm{cos}\left(\kappa t\right)\right)`$ (27a)
$`\phi (t)`$ $`=`$ $`\omega _ct+e\sqrt{2}\mathrm{sin}\left(\kappa t\right)`$ (27b)
Here $`k=\left(v_c/r(t)\right)^2\omega _{c}^{}{}_{}{}^{2}\left(1+2e\mathrm{cos}\left(\kappa t\right)\right)`$ where $`v_c=r_c\omega _c`$ is the circular velocity of the potential. We note that with this convention for $`r(t)`$, the dS is at perigalacticon at $`t=0`$. Replacing $`k(t)=\epsilon (t)\omega _{0}^{}{}_{}{}^{2}`$ and $`\phi (t)`$ in the equations (12) gives, at first order in $`e`$,

$`\ddot{x}+\omega _{0}^{}{}_{}{}^{2}\left(1\epsilon _c\mathrm{cos}\left(2\omega _ct\right){\displaystyle \underset{\pm }{}}b_\pm e\epsilon _c\mathrm{cos}\left(2a_\pm \omega _ct\right)\right)x=`$
$`\omega _{0}^{}{}_{}{}^{2}\left(\epsilon _c\mathrm{sin}\left(2\omega _ct\right)+{\displaystyle \underset{\pm }{}}b_\pm e\epsilon _c\mathrm{sin}\left(2a_\pm \omega _ct\right)\right)y`$ (28a)
$`\ddot{y}+\omega _{0}^{}{}_{}{}^{2}\left(1+\epsilon _c\mathrm{cos}\left(2\omega _ct\right)+{\displaystyle \underset{\pm }{}}b_\pm e\epsilon _c\mathrm{cos}\left(2a_\pm \omega _ct\right)\right)y=`$
$`\omega _{0}^{}{}_{}{}^{2}\left(\epsilon _c\mathrm{sin}\left(2\omega _ct\right)+{\displaystyle \underset{\pm }{}}b_\pm e\epsilon _c\mathrm{sin}\left(2a_\pm \omega _ct\right)\right)x`$ (28b)
where $`a_\pm =1\pm \sqrt{2}/2`$, $`b_\pm =1\pm \sqrt{2}`$ and $`\epsilon _c=\omega _{c}^{}{}_{}{}^{2}/\omega _{0}^{}{}_{}{}^{2}`$. Equations (28) have a form identical to equations (12) where $`\omega `$ is replaced by $`a_\pm \omega _c`$ and $`\epsilon `$ becomes $`b_\pm e\epsilon _c`$. Evidently two sets of resonant frequencies occur when $`\omega _c(1+\sqrt{2}/2)`$ or $`\omega _c(1\sqrt{2}/2)`$ is near $`\omega _0`$. These frequencies are sufficiently different that we treat them as independent resonant solutions. It follows that the form of the growing mode late-time solution for $`x(t)`$ and $`y(t)`$ in an elliptical orbit is

$`x(t)`$ $`=`$ $`\alpha ^{}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\left(a_\pm \omega _c\right)^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\left|b_\pm \right|e\epsilon _ct\right]\mathrm{cos}\left(a_\pm \omega _ct+\theta _\pm \right)`$ (29a)
$`y(t)`$ $`=`$ $`\alpha ^{}\mathrm{exp}\left[{\displaystyle \frac{\sqrt{\left(a_\pm \omega _c\right)^24\omega _{1}^{}{}_{}{}^{2}}}{2}}\left|b_\pm \right|e\epsilon _ct\right]\mathrm{sin}\left(a_\pm \omega _ct+\theta _\pm \right)`$ (29b)
where $`\omega _1`$ has been defined by $`\omega _0=a_\pm \omega _c+b_\pm e\epsilon _c\omega _1`$. Here $`\mathrm{tan}\theta _+=H_+`$ and $`\mathrm{tan}\theta _{}=1/H_{}`$ where $`H_\pm `$ has the form of equation (21) with $`\omega `$ replaced by $`a_\pm \omega _c`$. Comparing to the circular orbit solution we see that the expansion parameter $`\epsilon `$ becomes $`\frac{\omega _{c}^{}{}_{}{}^{2}}{\omega _{0}^{}{}_{}{}^{2}}eb_\pm `$ so that instability occurs when $`\left|\omega _0a_\pm \omega _c\right|<\omega _0\frac{\omega _{c}^{}{}_{}{}^{2}}{\omega _{0}^{}{}_{}{}^{2}}e\left|b_\pm \right|/2`$. This expression is comparable to the circular orbit solution (26) but resonance is “easier” to reach even as the strength of the tidal force compared to the binding force (which is given by $`\epsilon _c=k_c/\omega _{0}^{}{}_{}{}^{2}=\left(\omega _c/\omega _0\right)^21/a_{\pm }^{}{}_{}{}^{2}`$) can be weaker by a factor of 3 or more compared to the circular orbit.

Expanding $`\phi (t)`$ and $`\epsilon (t)`$ to order $`3`$ in $`e`$ yields additional resonant behavior. Analogous to the above expansion we find two new sets of excitation frequencies. The additional terms in the differential equation for $`x(t)`$ and $`y(t)`$ above have the form $`d_{2\pm }e^2\epsilon _c\mathrm{cos}\left(2b_\pm \omega _ct\right)`$ and $`d_{3\pm }e^3\epsilon _c\mathrm{cos}\left(2c_\pm \omega _ct\right)`$ – with corresponding $`\mathrm{sin}()`$ terms as in (28). Here $`c_\pm =1\pm \frac{3}{2}\sqrt{2}`$ and $`d_{n\pm }`$, $`n\{2,3\}`$, are numerical factors of order unity. Thus, the higher order terms lead to excitations at frequencies where $`\omega _0`$ is $`b_\pm \omega _c`$ (order 2) and at $`c_\pm \omega _c`$ (order 3). The growth rates of these modes and instability range for $`\omega _0`$ are smaller by corresponding powers of $`e`$ so we expect to excite higher order modes only as the orbital eccentricity increases.

#### 2.2.5 Projection on the line of sight

An interesting, and perhaps observable signature from the simulations, is the elliptical shape caused by the strong geometric asymmetry in the growing mode solution for the distance of a star away from the center of the dS. Here we consider how the induced elliptical shape varies with respect to the line-of-sight from the MW center-of-force toward the dS center. Consider, for the sake of generality, that the excitation is encountered at a frequency $`a\omega _c`$ where $`a`$ can be chosen from the set $`\{1,a_\pm ,b_\pm ,c_\pm \}`$, with the corresponding exponential term $`E_a(t)`$. Thanks to equations (6b) and (6c), we obtain (where $`t=0`$ corresponds to perigalacticon)

$`\left(x\widehat{x}+y\widehat{y}\right)\widehat{r}`$ $`=`$ $`E_a(t)\mathrm{cos}\left(\left(a1\right)\omega _ct+\theta \right)`$ (30a)
$`\left(x\widehat{x}+y\widehat{y}\right)\widehat{\phi }`$ $`=`$ $`E_a(t)\mathrm{sin}\left(\left(a1\right)\omega _ct+\theta \right)`$ (30b)
In the circular case ($`a=1`$) the dS is expanded along an axis rotated away from the line-of-sight to the MW by an angle $`\theta `$ with $`\mathrm{tan}\theta =H`$. For elliptical orbits the bar appears to rotate. At order $`n`$ the bar formed by the tidal interaction will be turning at a frequency $`n\kappa /2`$. Thus as the dS orbits between successive perigalacticon the bar will turn $`\frac{n}{2}`$ times.

## 3 Numerical Calculations and Experiments

These analytic calculations suggest that a broad range of dS-MW orbit dynamic conditions can lead to parametric tidal dS excitation. Evidently this excitation can occur even when the orbital and dS galaxy resonant frequencies are widely “mistuned.” Although our analytic calculations simplified the force equation terms from the exact spatial and temporal dependence of the dS self-gravity, we can show by numerical methods how many of the properties are retained in more realistic self-gravitating stellar systems.

In the discussion below we assume a system of units where the gravitational constant $`G=1`$. We fix one numeric length unit to a physical scale of $`1\mathrm{kpc}`$ and we’ve taken the dS to have a total initial mass of $`10^6`$ solar masses. With the gravitational constant set to unity, this implies a numerical time unit of $`415\mathrm{Myr}`$. The logarithmic external Milky Way potential has been chosen to have the form $`\varphi (r)=v_{c}^{}{}_{}{}^{2}\mathrm{log}(r/r_0)`$ where we assume a circular velocity $`v_c=220\mathrm{km}.\mathrm{s}^1`$ (or 100 numeric units).

We first consider the problem of a more realistic non-harmonic, but still static, dS potential/density model. This is addressed by numerical integration of the modified two-dimensional Mathieu differential equation. Finally we consider direct N-body calculations to explore the collective and time-dependent effects of the dS potential as it responds to an external tide.

### 3.1 Plummer Potential Solutions

The Plummer model, $`\varphi (r)=M/\sqrt{r^2+b^2}`$ (*cf* Binney and Tremaine 1987), is a rough but useful approximation for a self-gravitating dS stellar distribution. We use this model both to initialize our numerical simulations and in the differential equation integrations below. For example, in this case the local cartesian $`x`$ acceleration of a star can be written as $`\ddot{x}=x\omega _{0}^{}{}_{}{}^{2}\frac{b^3}{(b^2+r^2)^{3/2}}`$ where $`\omega _{0}^{}{}_{}{}^{2}=4\pi \rho _0/3`$ and $`\rho _0`$ is the central density. Notice that the parameter $`b`$ conveniently describes potentials which range from completely unbound ($`b=0`$), to homogenous, purely harmonic ($`b\mathrm{}`$) models. Substituting for the harmonic dS force term in eqs. (12) we obtain

$`\ddot{x}+x\omega _{0}^{}{}_{}{}^{2}{\displaystyle \frac{b^3}{(b^2+r^2)^{3/2}}}x\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{cos}2\phi `$ $`=`$ $`y\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{sin}2\phi `$ (31a)
$`\ddot{y}+y\omega _{0}^{}{}_{}{}^{2}{\displaystyle \frac{b^3}{(b^2+r^2)^{3/2}}}+y\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{cos}2\phi `$ $`=`$ $`x\epsilon \omega _{0}^{}{}_{}{}^{2}\mathrm{sin}2\phi `$ (31b)
We were unable to construct analytic solutions to equations (31), except for the cases $`b=0`$, or $`b\mathrm{}`$. Notice that this more general equation, preserves the form of the two-dimensional parametric oscillator but with a spatially variable “resonance” frequency $`\omega _{p}^{}{}_{}{}^{2}(r)=\omega _{0}^{}{}_{}{}^{2}b^3/(b^2+r^2)^{3/2}`$ and with $`\epsilon \epsilon \omega _{0}^{}{}_{}{}^{2}/\omega _{p}^{}{}_{}{}^{2}`$. For example, stars on circular orbits with initially constant $`r`$ should be tidally excited at the first Mathieu equation resonance condition when the dS orbit satisfies $`\omega =\omega _p(r)`$. We can expect that, to the extent that individual stellar orbits aren’t circular and there is a broad range in $`\omega _p(r)`$, the dS tidal response should be a broad function of its orbital frequency, $`\omega `$. For finite $`b`$ we also note that $`\omega _p(r)<\omega _0`$ so that resonant behavior must occur at lower frequencies than the dS harmonic oscillation frequency.

We can explore this solution space numerically. To distinguish stable and unstable solutions as a function of orbital frequency $`\omega `$ we numerically integrate eqs. (31) using a fourth-order Runge-Kutta algorithm. For example, according to the analytic solution, with $`b\mathrm{}`$, $`\omega _0=4`$, and $`\epsilon =0.1`$ we should obtain solutions for $`r(t)=\sqrt{x^2+y^2}`$ which are growing only for dS orbits with circular frequencies $`\omega \omega _0`$. To verify this we integrate an ensemble of solutions with $`\omega `$ between 0 and 8 over $`t=0`$ to $`50`$ and compute the mean square radius over time for each solution. Here $`x(0)=1`$ and the stellar orbits were chosen with a zero initial velocity condition to yield radial orbits. Figure 2 plots this mean solution radius for each dS orbital frequency.

Figure 2 illustrates several features of the analytic parametric solution. First, instability occurs when the orbital frequency is close to the dS resonant frequency ($`\omega _p`$) and, unlike in the one-dimensional parametric oscillator, the next resonance where the orbital frequency is one-half of the oscillation frequency is stable. This agrees with our analytic result. It is perhaps surprising that a circular orbit with half the orbital frequency of the dS resonance does not lead to instability – since the tide force has 180-degree symmetry. This is a consequence of the compressive and expansive parts of the tensor tide interaction.

The instability frequency width is determined by the tidal strength and form of the MW potential (through $`\epsilon `$ in eq. (26)). With $`\omega _0=\omega _c=4`$ and $`\epsilon =0.1`$ we expect a resonance width of 0.4 – which is confirmed in figure 2. The growth rate of $`r(t)`$, as determined by eq. (25), is also reproduced in the numerical calculation. For $`\omega _1=0`$ we expect a stellar orbit to expand like $`\mathrm{exp}\omega _0\epsilon t/2.=\mathrm{exp}0.2t`$ which by $`t=50`$ yields the amplitude plotted here.

The other extreme, where $`b=0`$, corresponds to a dS which is completely unbound, i.e. where the dS self-gravity is unimportant. For example, at late times a dS might be tidally distended so that its self-gravity becomes negligible. A naive interpretation of the MW effect here might be to conclude that unbound dS member stars should be tidally accelerated to large distances from the dS center-of-mass. Figure 3 shows the logarithm of the rms stellar orbit radius for the unbound ($`b=0`$) case with the same initial conditions, range of dS-MW orbital frequencies, and integration parameters as in fig. 2. We see that for static and slowly rotating tides the stellar system does indeed rapidly blow up – saturating our numerical dynamic range before $`t=50`$ for orbits near $`\omega =0`$. Analytic solutions using computer aided symbolic manipulation imply that the leading growing mode term in the solution varies as $`\mathrm{exp}\left[\sqrt{\epsilon \omega _{0}^{}{}_{}{}^{2}\omega ^2}t\right]`$. Interestingly, as this shows, for dS orbits with short enough periods satisfying

$$\omega >\sqrt{\epsilon }\omega _0=\sqrt{k}$$
(32)
localized near the dS center-of-mass, despite being unbound to the dS.

Intermediate potential models ($`0<b<\mathrm{}`$) yield unstable orbits depending on the initial conditions of the stellar orbit. The dS orbital frequency which leads to instability of a dS star depends on the initial radius of the star with respect to the dS center-of-mass in units of $`b`$, and the dS central density which determines $`\omega _0`$. The previous two cases have shown that dS stars which begin their orbits at normalized distances of 0 and $`\mathrm{}`$ develop secular instability for dS-MW orbits that satisfy, respectively, $`\omega \omega _0`$ and $`\omega 0`$. Intermediate cases lead to intermediate resonant frequencies.

The top panel of figure 4 plots the orbital expansion of a star starting from $`r/b=1`$ using the same integration parameters ($`\omega _0=4`$) as the previous Runge-Kutta solutions. Secular growth now occurs at lower frequencies near $`\omega 2.8`$. In a Plummer potential-density model the half-mass radius occurs at $`1.4b`$. The lower panel of figure 4 demonstrates for this radius how the dS resonant orbital frequency continues to decrease for stars at larger normalized dS distances. Evidently stars near the half-mass radius are tidally excited when $`\omega =1.1`$. The initial position of a star within the dS potential helps to determine its fate in the time-dependent MW tide.

These solutions describe a system where the tide amplitude $`\epsilon `$ and the orbital frequency $`\omega `$ are independent. A logarithmic MW potential satisfies $`\epsilon \omega _{0}^{}{}_{}{}^{2}=\omega ^2`$. Thus for increasing orbital frequency, $`\epsilon `$ in eqs. (31) is also increasing. Qualitatively we find that this increases the orbital frequency domain over which resonance expands the dS. Systems which are resonant with $`\omega \omega _p`$ for fixed $`\epsilon `$ in the above discussion are now tidally disrupted for $`\omega \omega _p`$.

#### 3.1.1 Elliptical Orbits

Elliptical orbits are also described by eqs. (31) by taking $`\epsilon \epsilon (t)`$ and using the actual form of $`\phi (t)`$ and $`\epsilon (t)`$ for an elliptical orbit. Here we assume a logarithmic form for the MW potential so that $`\epsilon `$ and $`\omega `$ are not independent.

$`\phi (t)`$ $`=`$ $`1+e\sqrt{2}\mathrm{sin}\kappa t`$ (33a)
$`\epsilon (t)`$ $`=`$ $`{\displaystyle \frac{\omega ^2}{\omega _{0}^{}{}_{}{}^{2}}}{\displaystyle \frac{1}{(1e\mathrm{cos}\kappa t)^2}}`$ (33b)
where $`\kappa =\sqrt{2}\omega `$ is the epicyclic frequency.

We compute the mean distance of a star from the dS center as we did in the previous section. We assume a constant perigalacticon distance with $`e=1\omega /2`$, take $`\omega _0=6.7`$, $`x(0)/b=1.4`$, and let $`\omega `$ vary from 1.0 to 2.0 so that we can compare with the results of section 3.2.3. Figure 5 shows that a broad range of circular frequencies lead to parametric excitation and, as we will see, this general behavior also appears in more realistic self-gravitating systems. We also note that a circular orbit ($`\omega =2`$) can be less dispersive than an elliptical orbit ($`\omega <2`$) even if the tide is always stronger in the circular case.

### 3.2 Direct N-body Calculations

In order to account for the self-gravity and collective dynamics of many stars within the dS we have used a direct N-body calculation. The numerical simulations were computed using the modified treecode v1.4 (Barnes 1990). Typically between 1000 and 10,000 particles were introduced into simulations which included several different “external” forces designed to study variously dS internal oscillations and the tidal coupling of the stellar system to an imposed large-scale logarithmic external MW potential.

The effects of numerical viscosity and small spatial scale potential fluctuations can be important in these calculations and we used a range of integration time step and potential “softening” parameter to understand this dependence. Most of the simulations were done with a time step of 0.008 units. We explored the effects of the interparticle potential softening parameter over values ranging from 0.025 to 0.3. Significant damping of oscillations appears to be minimal in these models with softening parameter $`a>0.1`$.

With a crossing time of about 0.2 time units the relaxation timescale for $`N=1000`$ simulations was at least 2 time units and much larger for softer potentials. Many of our results only depend on measuring differential changes between simulations, e.g. to determine the frequency dependence in the particle ejection rate between models with different orbital frequencies. Some results are computed from relatively long model simulation runs.

Several different initial particle configurations were generated although most of the results we present here used the treecode Plummer model realization (*cf* Aarseth et al. 1974 as implemented by Barnes 1990). The rms particle radius was typically $`0.1\mathrm{kpc}`$ and the Plummer parameter was typically $`b=0.16`$ for a softening parameter of $`a=0.1`$. To achieve this an initially isotropic velocity dispersion of about $`1\mathrm{km}.\mathrm{s}^1`$ was used and allowed to relax before imposing any external time-dependent tidal forces.

We have not tried to reproduce any of the dS in detail since their MW orbital uncertainties and the non-uniqueness of the model solutions make detailed comparisons difficult to interpret. Nevertheless we believe these N-body calculations are most representative of a dS like Draco or Ursa Minor.

#### 3.2.1 Oscillations

The parametric galaxy model described by eqs. (12) neglects the effect of the oscillating potential of the dS. Thus, we can expect normal mode galaxy oscillations to dominate even the growing parametric tidal modes if the system’s normal mode frequencies are ever excited by the MW orbit. To explore this numerically we need to understand the normal mode oscillation spectrum of our dS system.

Based on previous experience with particle-mesh calculations (KM) involving $`10^5`$ particles or more we expected free galaxy oscillations to be readily generated by small deviations from equilibrium in the initial conditions of the particle configuration. For example, Miller and Smith (1999) had difficulty damping and suppressing such oscillations. Even though the total energy in our simulations of isolated dS was conserved to within 0.03 %, we initially had difficulty detecting oscillations. To excite an oscillation spectrum we applied an impulsive radial stretching acceleration at $`t=0`$ and looked for the resulting ringing. From a temporal fourier analysis of projections of the particle motion onto the fundamental eigenmode, $`\stackrel{}{\xi }(r)=\stackrel{}{r}`$, we measure the oscillation spectrum.

In detail we compute a timeseries by summing over particles (labeled with $`i`$)

$$c_j=\frac{1}{N}\underset{i}{}\stackrel{}{r_i}(t_j)\stackrel{}{v_i}(t_j)$$
(34)
The logarithm of the temporal fourier transforms of $`c_j`$ from two different simulation runs are indicated in Fig. 6. For an homogeneous sphere, the power spectrum should peak at the oscillation frequency $`\omega _0`$. Increasing the softening parameter both changes the peak frequency in the power spectrum and the peak amplitude response to the excitation. The virial equation (13) yields a reasonable estimate of the observed oscillation peak frequency. For example with softening parameter of $`a=0.1`$ and using the observed N-body density distribution we compute a frequency of $`\omega _0=3.54`$, which is consistent with the peak in Fig. 6 (lower panel).

The shorter timestep minimizes the integration errors and the intrinsic numerical damping noise (Figure 7). Similarly by increasing the softening parameter we can decrease the effects of integration error from the short range particle interactions. Figure 6 confirms that as the softening parameter increases there is an enhancement in the fundamental mode amplitude.

It is difficult to compute the true physical damping of galaxy oscillations but comparing our results with Miller and Smith (1999) suggests that the direct N-body calculations have a larger numerical damping than particle-cell approximations. We expect the dS instability growth rates from these N-body simulations to be underestimated compared to the analytic and Runge-Kutta solutions.

#### 3.2.2 Parametric resonance

In order to demonstrate Eqs. (25) for a dS in a circular orbit, we decouple the strength of the tidal force from its time dependence. We isolate the dependence of the tide amplitude and frequency by fixing the amplitude while independently setting the orbital frequency ($`\omega `$). This is analogous to the approach KM used in their particle-mesh calculations. It allows the frequency dependence of the resonance to be distinguished from the effect of a growing static tide strength as the dS galactocentric distance is decreased in order to increase $`\omega `$. Evidence of this resonant behavior is illustrated in Figure 8 which shows a snapshot of the particle positions for three different models corresponding to $`\omega =\omega _p/2`$, $`\omega _p`$, and $`2\omega _p`$. The direction of the instantaneous tide is indicated by the line on each plot.

Figure 9 shows the effect of changing $`\omega `$ while keeping the tidal amplitude (determined by $`\epsilon `$) fixed. Here the total number of particles ejected from the dS is plotted versus time and frequency $`\omega `$. Recall that $`\omega _p`$ is approximately 2 numerical units and the tide amplitude corresponds to $`\epsilon =0.5`$. Successive curves diplaced upward in this plot show how the ensemble of models have ejected more particles at later times.

Several points are illustrated by fig. 9. The most important is that, as expected from the Mathieu equation solutions, the resonant behavior (measured by particle ejection from the dS) extends to frequencies relatively far from the effective resonant frequency $`\omega _p`$. According to our solution to Eq. (25) and for our choice of $`\epsilon `$ we expect an instability frequency range $`\mathrm{\Delta }\omega /\omega _p=\epsilon /2=1/4`$ which is consistent with the numerical results. It is also notable that there is no second order resonance at $`\omega _p/2`$, consistent with the analytic predictions of the coupled Mathieu equation model.

We also find in the N-body solutions that for dS circular frequencies $`\omega 1.5`$ the ejected stars migrate toward $`r\mathrm{}`$ as the dS evolves. For $`\omega 1.6`$ the ejected stars (more distant than 1 kpc from the core) appear to remain in a bounded region of space, at least over the duration of these simulations. This result is anticipated by eq. 32. In these simulations $`k=2`$ so the cutoff frequency should be about $`\omega =1.4`$, in fair agreement with the N-body calculation.

#### 3.2.3 Elliptical orbits

More realistic orbital calculations also demonstrate how the dynamical tide interaction affects the internal dynamics of a dwarf spheroidal. We adapted the numerical calculations described above to describe a self-gravitating ensemble of masses orbiting within a logarithmic external potential. We use a 1024 point dS galaxy characterized by oscillation frequency, $`\omega _07`$ numeric units. Setting this dS onto a circular orbit at $`50\mathrm{kpc}`$ from the galactic center yields a stable system with 75% of the initial stars and with a particle loss rate of less than 50 stars over 250 time units.

To generate a family of elliptical orbits we launched the dS at $`50\mathrm{kpc}`$ along the $`x`$ axis with variable speed $`v_0v_c`$ along the $`y`$ axis in a constant galactic potential $`\varphi (r)=v_{c}^{}{}_{}{}^{2}\mathrm{log}(r)`$. Thus in our units when $`\omega _c=2`$ the orbit is circular and in general the eccentricity, $`e`$, varies like $`e=1\omega _c/2`$. Thus, the perigalacticon of all runs is exactly $`50\mathrm{kpc}`$ and the mean tidal force on the dS is largest for the circular case ($`e=0`$, $`\omega _c=2`$) and decreases with $`\omega _c`$ or as the eccentricity $`e`$ increases.

Figure 10 shows how the number of particles ejected from the dS varies with orbital ellipticity and time. Here the orbit time is measured from the injection of the equilibrated dS into the logarithmic potential. The broad frequency response of the dS in an elliptical orbit is confirmed from the Runge-Kutta solutions (Figure 5). It is interesting that elliptical orbits, with weaker tides than a circular orbit can lose more stars. Eccentric variable tides can eject stars where a stronger constant tide cannot.

### 3.3 Velocity Dispersion

The velocity dispersion of our parametrically excited dS galaxies varies by a large factor over their orbits and can be a strong function of the extent of the dS over which the dispersion is calculated. For example a general feature of these models is that the dispersion increases outward from the core of the dS. Nevertheless it is possible to inflate the dispersion by an order of magnitude or more depending on the dS resonant and orbital conditions. For example, in Fig. 11 we plot the instantaneous line-of-sight velocity dispersion for the model above with orbital eccentricity of 0.5. While this is a relatively high eccentricity orbit, it is not particularly close to resonance (with a slow growth rate) and we plot the dispersion after the growing mode is established. It is clear that at late times the dispersion increases by more than a factor of 10 yielding virial mass estimates exaggerated by a factor of 100. Figure 12 shows this model when the dispersion was 10 times larger than its initial configuration. This system is not particularly close to dissolution and it has lost only $`300`$ stars out of 734 when it was launched in its orbit. It also seems to retain its central concentration and surface brightness profile.

### 3.4 Interpreting galaxy morphology from N-body calculations

Our analytic parametric oscillator calculations and numerical integrations of more realistic Plummer model systems showed that resonant excitation of dS stars should produce elliptical systems with their long axis rotated by an angle $`\theta `$ away from the direction toward the force center in the plane of the orbit. For a circular or elliptical orbit this angle is calculated from $`\mathrm{tan}\theta =H`$ (section 2.2.3). At the first parametric resonance this angle is 45 degrees. Further from the resonance condition this angle varies from an orientation along the center-of-force direction to being perpendicular to it. The existence of a nonzero dS bar angle is evidence of parametric resonance.

We found in section 2.2.5 that this angle also describes the dS elongation in an elliptical orbit at perigalacticon. During the orbit the bar of a dS should rotate with respect to the center-of-force direction with increasing frequency, depending on the order of the resonance in orbital ellipticity, $`e`$.

We have computed the direction of the dS bar from the moment of inertia tensor of the dS stars. Figure 13 shows a series of snapshots from the simulation with $`e=0.22`$. In this figure the center-of-force is always to the left and $`\theta (t)`$ is the angle between the two solid lines which intersect at the center of the dS. Figure 14 shows how the angle varies with time for this calculations. Each spike in the graph corresponds to one rotation of the bar. We find that the rotation rate increases with increasing eccentricity (as it should) as higher order modes in powers of $`e`$ are excited.

The bar rotation indicated in Fig. 14 is not strictly sinusoidal with the epicycle period. We also note that at perigalacticon our dS systems tend to have the bar directed toward the force center ($`\theta =0`$). Recall from our discussion of the Plummer model that not all stars in the dS are characterized by the same resonance frequency $`\omega `$, so that stars near the center of the dS have larger frequencies than stars at larger distances. Thus for given orbital circular and epicyclic frequencies the response of many of the dS stars is non-resonant. Thus near perigalacticon many of the stars (especially near the core of the dS) simply respond to the tide like any non-resonant fluid system. Our moment of inertia calculation is thus controlled by contributions from both resonant and non-resonant stellar components. This can lead to a more complex behavior of the bar pointing direction than our model predicts in detail.

The outer parts of the dS (see Fig. 13) are curved as if $`H`$ depends on distance from the dS. This is anticipated from the Plummer calculations since the effective $`\omega _p`$ decreases outward. The curvature in the bar at large distances from the dS raises the possibility that snapshot observations of real tidally excited stellar systems may provide useful constraints on the form of the dS potential.

## 4 Comparing the parametric oscillator model with observations

Our earlier comparisons of dS dynamics and morphology with a simpler resonant tidal excitation model depended on the commensurability of the dS oscillation frequency (and damping width) and its MW orbital driving frequency to generate elliptical systems with large velocity dispersions. The parametric oscillator model presented here more accurately describes the dS-MW interaction and is a more precise explanation for dS properties. In particular their elliptical shape and sometimes large velocity dispersions is a natural consequence of this model.

Parametric dS galaxy oscillations effectively increase the dS gravitational interaction cross-section with the MW. Properly accounting for the time-dependent tidal interaction of a dS with the MW, even when the dynamical frequency of the dS ($`\omega _0`$) is not commensurate with its orbital or epicyclic frequency around the MW, is critical to describing the dS dynamics. We have shown how the frequency domain near the dS fundamental frequency which describes resonance is not characterized by the galaxy oscillation damping time but by the tide amplitude ($`\epsilon `$). For a logarithmic MW potential in a circular orbit the fractional resonant-growth frequency domain is simply $`\omega _{c}^{}{}_{}{}^{2}/\omega _{0}^{}{}_{}{}^{2}`$ (the square of the ratio of the circular orbit frequency and the internal galaxy oscillation frequency). Equation 26 implies that any dS with oscillation frequency less than about 1.4$`\omega _c`$ will be tidally disrupted.

As we noted above, elliptical orbits can be even more likely to disrupt dS stellar systems because of their rich epicyclic harmonic structure. Eccentric orbits lead to disruption of “stiffer” dS with internal frequencies as large as 1.7, 2.4, and 3.1 (or larger) times the circular orbital frequency – corresponding to 1st, 2nd, and 3rd order terms in orbital eccentricity. Note that parametric resonance increases the mean dS central distance of its member stars and their velocities (and velocity dispersions). The exponential growth time depends on how close the dS is to a resonance condition but as eq. (25) shows the growth times of, for example, the velocity dispersion can be comparable to the orbit period. In general we expect larger growth rates as the ratio of circular to internal dS frequency increases toward unity.

Despite the undetermined orbital characteristics of the dS, it is interesting to compare our best estimate of $`\epsilon `$ for each of the MW dwarf spheroidals with their kinematic properties. It is our contention that the dS kinematics are not dominated by dark matter halos, but by their MW tidal interaction. Thus we estimate their mass from their luminosity and not their velocity dispersion, since the dS are not in virial equilibrium. Core (half-light) radii ($`r_c`$) and galactocentric distances ($`d`$) are used from the compilation by Mateo (1998). A pulsation frequency $`\omega _0`$ is computed from the VanderVoort approximation $`\omega _0=\pi \sqrt{G\rho _0}`$. We compute the central density $`\rho _0`$ from the observed total luminosity (assumed equal to the total dS mass in solar units) and the expression $`\rho _0=\frac{3M(0.64)^3}{4\pi r_c^3}.`$ This is exact when the density distribution equals a Plummer model and is a reasonable empirical approximation for the MW dS. We also take $`\omega _c=v_c/d`$ (with $`v_c=220\mathrm{km}.\mathrm{s}^1`$) as the circular frequency around the MW. Adequate for our purposes, we note that the parametric growth rate increases with $`\epsilon `$ and the ratio $`A=\omega _c/\omega _0`$. If the dynamics of the dS are dominated by tidal interactions we expect observed dS virial “M/L” ratios to increase with $`A`$. Figure 15 plots $`A`$ against M/L. The general trend of increasing “M/L” with $`A`$ is in good agreement with the parametric model, especially given the uncertainty in the dS orbital parameters.

Our “snapshot” knowledge of the Local Group Dwarf Spheroidals means that we cannot measure their orbital ellipticity (although see Kuhn 1993 for a discussion of this point) nevertheless, knowledge of the dS radial velocities and galactocentric distances is sufficient to suggest that at least the nearest dS must be tidally inflated by parametric resonance.

Another interesting and perhaps observationally verifiable consequence of parametric excitation is that the long axis of the dS bar should, in general, be inclined with respect to the separation vector between the MW center and the dS. Orbits which are nearly circular should yield a dS with the leading edge closer to the MW than the trailing edge and inclined at an angle of 45 degrees. Elliptical orbits produce a rotating bar whose phase with respect to the radius vector direction depends on how close the dS is to resonance and the proximity of the dS in its orbit to perigalacticon.

## 5 Summary and Conclusions

We have demonstrated an analytic model of time dependent dS-MW tidal interactions which extends the Mathieu Equation for parametric oscillations to two dimensions. Several important general conclusions follow from this model:

1. Secular instability can cause dS star orbits and velocities to grow, depending on the ellipticity and circular frequency of the MW orbit and the resonant frequency of the dS.
2. Exponential growth can occur over a broad range of frequencies, not just near the characteristic dS harmonic frequency. Growth times can be comparable to the MW orbital period.
3. Direct integration of more realistic, non-harmonic dS potentials, yield the same qualitative behavior predicted by the Mathieu equation.
4. Self-gravitating N-body calculations also confirm the existence of broad resonance conditions that lead to stellar ejection from the dS system, elliptical bar formation, and inflated (non-virial) velocity dispersions.
5. Model expectations for the parametric growth rate in MW dwarf galaxies have been contrasted with the observed M/L measurements for the 8 Milky Way dS. The expected trend of increasing M/L with $`\omega _c/\omega _0`$ is confirmed.

Considerable attention has been given to finding dark matter density distributions which will both stabilize the luminous components of the dS so they survive MW tidal encounters, but not so much that they don’t allow stars to be ejected to populate “tidal tails”. These tails are now clearly observed by several methods (Kuhn et al. 1996, Smith et al. 1997, Martinez-Delgado et al 2001) and are an observational hurdle that any dS model must pass. While it appears that there may be dark matter models which cannot be ruled out (*cf* Mayer et al. 2002) here we have developed a consistent model of MW-dS interactions which does not require any invisible dS mass component and which appears to account for the morphology and dynamics of the dS. The prediction of extratidal dS stars in the simple KM resonance model were later observed. We can hope that the considerably refined predictions of the parametric oscillator model (for example the inclined bar) may also be empirically established.

Acknowledgements

We’re grateful to Josh Barnes for help using his treecode v1.4.

JJF is grateful to the IfA for support in many forms during his internship from the ENS – and to his grand-father who left us while this work was being done.


text:
# 1 Introduction

## 1 Introduction

The Landau theory of second-order phase transitions? and its realization for superconductors, the Ginzburg-Landau (GL) gauge theory,? can be classified as belonging to the most illuminating theoretical achievements in XXth-century physics. The basic concepts? advanced in these theories often find applications in interdisciplinary research. The microscopic Bardeen-Cooper-Schrieffer (BCS) theory? makes it possible to calculate the parameters of the GL theory. Thus the phenomenology of superconductivity can be reliably derived once the parameters of the microscopic Hamiltonian are specified. Such a scheme ensures that there is no missing link between the microscopic theory and the material properties of the superconductors. On the other hand Bernoulli effect is one of the oldest effects known in the physics which in as a result of the energy conservation can be realized strictly speaking only for superfluids. For superconductor the Bernoulli theorem gives a relation between current density and the current induced contact potential difference at constant in thermal equilibrium electrochemical potential.

The purpose of the present paper is to systematize the results for the GL theory of disordered and anisotropic gap superconductors and to point out what new information for the fundamental parameters of superconductors can be extracted using the Bernoulli effect. We are describing in short how the experiments can be done and how the effective mass of Cooper pairs, the volume density of charge carriers and London penetration depth can be extracted from investigation of the Bernoulli effect. The proposed experiments require only standard low frequency electronic measurements. In parallel we analyze a simple methodical derivation of the London penetration depth for clean superconductors. We will start our analysis with the GL equation, then we will analyze disorder renormalization and experimental accessibility of the effective Cooper pair mass, we will analyse the Bernoulli effect and temperature dependence of the penetration depth. Finally we will conclude which experimental development we consider as most perspective.

## 2 Ginzburg-Landau equation

Consider now the gradient terms of the Ginzburg-Landau (GL) theory. A general expression for the tensor of the squared coherence lengths $`(\widehat{\xi }^2)_{\alpha \beta }`$ has been derived by Pokrovsky and Pokrovsky?:

$`(\widehat{\xi ^2})_{\alpha \beta }=`$ $`{\displaystyle \frac{\mathrm{}^2}{\left(4\pi k__\mathrm{B}T_\mathrm{c}\right)^2}}{\displaystyle \frac{1}{\chi ^2}}(\zeta _{3,0}v_\alpha v_\beta \chi ^2`$ (1)
$`+2x_c\zeta _{3,1}v_\alpha v_\beta \chi \chi +x_c^2\zeta _{3,2}v_\alpha v_\beta \chi ^2),`$
where

$$\zeta _{k,l}\zeta _{k,l}(x_c+1/2),\zeta _{k,l}(z)=\underset{n=0}{\overset{\mathrm{}}{}}(n+z)^k(n+1/2)^l,x_c=\frac{\mathrm{}}{2\pi k_\mathrm{B}\tau (T_c)},$$
(2)
$`\tau (T_c)`$ is the electron scattering time, $`v_\alpha `$ are the components of the velocity, and $`\mathrm{}`$ means averaging on the Fermi surface. In terms of $`\widehat{\xi ^2}`$ the GL equation for the space-dependent order parameter $`\mathrm{\Xi }(𝐫)`$ can be written as

$$\left(𝐃\widehat{\xi ^2}𝐃+ϵ+\frac{b}{a_0}\left|\mathrm{\Xi }\right|^2\right)\mathrm{\Xi }=0,$$
(3)
where

$`𝐃`$ $`={\displaystyle \frac{}{𝐫}}{\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}𝐀(𝐫),`$ $`\mathrm{\Phi }_0`$ $`=2\pi \mathrm{}/e^{},`$ (4)
$`|e^{}|`$ $`=2|e|,`$ $`ϵ`$ $`={\displaystyle \frac{TT_\mathrm{c}}{T_\mathrm{c}}}.`$ (5)
Equation (3) is the extremum (minimum) condition for the GL free energy

$$F[\mathrm{\Xi },𝐀]=𝑑𝐫\left\{\frac{\mathrm{}^2}{2}[𝐃\mathrm{\Xi }(𝐫)]^{}\stackrel{~}{M}^1𝐃\mathrm{\Xi }(𝐫)+a(T)|\mathrm{\Xi }(𝐫)|^2+\frac{1}{2}b|\mathrm{\Xi }(𝐫)|^4\right\},$$
(6)
where

$$(\stackrel{~}{M}^1)_{\alpha \beta }=\frac{2a_0}{\mathrm{}^2}(\widehat{\xi ^2})_{\alpha \beta },\text{ and }a(T)=ϵa_0.$$
(7)
Using the eigenvalues $`\xi _\alpha `$ and $`\stackrel{~}{M}_\alpha `$ of the tensors $`\widehat{\xi ^2}`$ and $`\stackrel{~}{M}`$ we can introduce the temperature-dependent GL coherence lengths and penetration depths, respectively ($`0<ϵ1`$):

$$\xi _{\alpha ,\mathrm{GL}}(T)\xi _\alpha /\sqrt{ϵ},\lambda _{\alpha ,\mathrm{GL}}(T)\lambda _\alpha /\sqrt{ϵ},$$
(8)
which satisfy the following GL relations

$$\frac{1}{\lambda _\alpha ^2}=\frac{(a_0/b)e^2}{c^2\epsilon _0}\frac{1}{\stackrel{~}{M}_\alpha }=\frac{2a_0^2e^2}{\mathrm{}^2\epsilon _0c^2b}\xi _\alpha ^2=8\pi ^2\mu _0\frac{a_0^2\xi _\alpha ^2}{b\mathrm{\Phi }_0^2}=\frac{8\pi ^2\mu _0}{\mathrm{\Phi }_0^2}T_c\mathrm{\Delta }C\xi _\alpha ^2,$$
(9)
where in Gaussian units $`\mu _0=4\pi `$ and $`\mathrm{\Delta }C=a_0^2/T_cb`$ is the jump of the specific heat at $`T_c`$ per unit volume. The explicit formulae for the GL coefficients $`a_0`$ and $`b`$ are given in Ref.?

## 3 Isotropic alloys. Disorder renormalization of the Cooper pair mass

Let us illustrate now the operation of the above general expressions on the important for the applications case of dirty isotropic alloys. In this case we have

$`v_\alpha v_\beta `$ $`={\displaystyle \frac{1}{3}}v_\mathrm{F}^2\delta _{\alpha \beta },`$ $`\chi _𝐩`$ $`=\mathrm{const},`$
$`(\widehat{\xi ^2})_{\alpha \beta }`$ $`=\xi ^2\delta _{\alpha \beta }`$ $`\stackrel{~}{M}_{\alpha \beta }`$ $`=\stackrel{~}{M}\delta _{\alpha \beta }.`$ (10)
Using the identity

$$\zeta _{3,0}+2x_c\zeta _{3,1}+x_c^2\zeta _{3,2}=\zeta _{1,2}$$
(11)
one obtains the classical result due to Gor’kov,see the textbook Ref. 8

$$\xi ^2=\frac{1}{3}\left(\frac{\mathrm{}v_\mathrm{F}}{4\pi k__\mathrm{B}T_\mathrm{c}}\right)^2\zeta _{1,2}=\xi _{\mathrm{clean}}^2\frac{\stackrel{~}{M}_{\mathrm{clean}}}{\stackrel{~}{M}(x)}.$$
(12)
Here we wish to recall some notation used in different works on the physics of superconuctivity (see, e.g., Ref. 9):

$`\xi _{\mathrm{clean}}`$ $`={\displaystyle \frac{\xi _0}{4}}{\displaystyle \frac{2\pi }{\gamma }}\sqrt{{\displaystyle \frac{7\zeta (3)}{12}}}={\displaystyle \frac{\xi _0}{4}}{\displaystyle \frac{2\mathrm{\Delta }(0)/k__\mathrm{B}T_\mathrm{c}}{\sqrt{\mathrm{\Delta }C/C_n(T_\mathrm{c})}}}=0.738\xi _0\xi _{0\mathrm{P}},`$
$`\xi _0`$ $`={\displaystyle \frac{\mathrm{}v_\mathrm{F}}{\pi \mathrm{\Delta }(0)}}={\displaystyle \frac{\gamma }{\pi ^2}}{\displaystyle \frac{\mathrm{}v_\mathrm{F}}{k__\mathrm{B}T_\mathrm{c}}}=0.1805{\displaystyle \frac{\mathrm{}v_\mathrm{F}}{k__\mathrm{B}T_\mathrm{c}}},`$
$`2\mathrm{\Delta }(0)`$ $`={\displaystyle \frac{2\pi }{\gamma }}k__\mathrm{B}T_\mathrm{c}=3.53k__\mathrm{B}T_\mathrm{c},\gamma =1.781\mathrm{},`$
$`{\displaystyle \frac{\xi ^2}{\xi _{\mathrm{clean}}^2}}`$ $`={\displaystyle \frac{\lambda _{\mathrm{clean}}^2}{\lambda ^2}}={\displaystyle \frac{\stackrel{~}{M}_{\mathrm{clean}}}{\stackrel{~}{M}(x)}}={\displaystyle \frac{\zeta _{1,2}\left(x_c+\frac{1}{2}\right)}{7\zeta (3)}}{\displaystyle \frac{1}{1+\frac{14\zeta (3)}{\pi ^2}x_c}}={\displaystyle \frac{1}{1+\frac{\xi _{0\mathrm{P}}}{l}}},`$ (13)
$`\xi _{0\mathrm{P}}`$ $`={\displaystyle \frac{7\zeta (3)}{2\pi \gamma }}\xi _0=0.752\xi _0={\displaystyle \frac{7\zeta (3)}{2\pi ^3}}{\displaystyle \frac{\mathrm{}v_\mathrm{F}}{k__\mathrm{B}T_\mathrm{c}}}`$
$`=0.1357{\displaystyle \frac{\mathrm{}v_\mathrm{F}}{k__\mathrm{B}T_\mathrm{c}}}={\displaystyle \frac{2\sqrt{21\zeta (3)}}{\pi ^2}}\xi _{\mathrm{clean}}=1.018\xi _{\mathrm{clean}},`$
and finally $`l=v_\mathrm{F}\frac{\tau }{2}.`$

The 2% difference between the length $`\xi ,`$ introduced by Ginzburg and Landau,?,? and $`\xi _{0\mathrm{P}}`$ introduced by Pippard,? is experimentally inaccessible. So is the several percent difference between the exact microscopic result (13) and its Padé approximant shown in Fig. 1. The Padé approximant for the mass renormalization can be easily derived? within framework of the Pippard-Landau theory, i.e. using the generalization of the local GL theory and nonlocal Pippard electrodynamics. These are generalized? by inserting the Pippard kernel between the two gradients in the GL expression for the free energy (6):

$$𝑑\mathrm{r}𝑑\mathrm{r}^{}\underset{\alpha ,\beta }{}[D_\alpha \mathrm{\Xi }(𝐫)]^{}\delta (𝐑)D_\beta \mathrm{\Xi }(𝐫^{})$$
(14)
$$𝑑\mathrm{r}𝑑\mathrm{r}^{}[𝐃\mathrm{\Xi }(𝐫)]^{}\frac{3𝐑𝐑}{4\pi \xi _{0\mathrm{P}}R^4}\mathrm{e}^{R/\xi _\mathrm{P}}𝐃\mathrm{\Xi }(𝐫^{}),$$

where $`𝐑=𝐫𝐫^{},`$ $`\xi \xi _\mathrm{P},`$ and

$$\frac{1}{\xi _\mathrm{P}}=\frac{1}{\xi _{0\mathrm{P}}}+\frac{1}{l},\xi _{\mathrm{GL}}(T)\frac{\xi _\mathrm{P}}{\sqrt{|ϵ|}}.$$
(15)
## 4 Experimental accessibility of the effective Cooper pair mass $`(T=0)`$

The tensor introduced in Eq. (7) acquires dimension of mass if we renormalize the GL order parameter so as to have dimension of density. The total charge carrier density $`n_{\mathrm{tot}}`$ is a fundamental theoretical notion. Thus it is of principle interests whether this quantity is experimentally accessible through the equilibrium thermodynamic properties of the vortex-free Meissner-Ochsenfeld phase. An approach to this problem has been outlined in Ref. 11 and is based on the current-induced contact-potential difference at the surface of a superconductor? (we call it the London-Hall effect):

$$\mathrm{\Delta }\phi =_\mathrm{H}\frac{B^2}{2\mu _0}.$$
(16)
For superconductors characterized with local electrodynamics $`\mathrm{\Delta }\phi `$ is a manifestation of the energy conservation law and the Bernoulli theorem? for charged superfluids in thermodynamic equilibrium:

$$\mathrm{\Delta }\phi =\frac{1}{2\epsilon _0c^2}_\mathrm{H}\lambda ^2(T)j^2.$$
(17)
For $`T=0`$ this equation reads as one fluid energy conservation low

$$\frac{1}{2}M_{\mathrm{CP}}v^2+e^{}\mathrm{\Delta }\phi =0,$$
(18)
where $`M_{\mathrm{CP}}`$ is the effective mass of Cooper pairs. For anisotropic superconductors $`\lambda `$ in Eq. (17) corresponds to the direction of the current density $`𝐣`$ at the superconductor surface. All radio frequencies are too small compared to the typical gap parameters of superconductors. Hence $`\mathrm{\Delta }\phi `$ can be measured using a lock-in amplifier in conjunction with a suitable low-noise preamplifier. The superconductor surface under investigation, as shown in Fig. 2, plays the role of one of the plates of the capacitor with capacitance $`c`$ and resistance $`r_c`$.

The condition under which $`\mathrm{\Delta }\phi `$ can be directly measured reads

$$Cf_{\text{lock-in}}\frac{1}{R_C}+\frac{1}{R_{\text{lock-in}}},$$
(19)
where $`f_{\text{lock-in}}`$ and $`R_{\text{lock-in}}`$ are, respectively, the operating frequency and the internal resistance of the lock-in amplifier. Typical voltages are $``$nV but for high-$`T_\mathrm{c}`$ superconductors the Bernoulli signals will be considerably stronger. The key technological problem lies in the the properties of the superconductor-insulator interface—its quality should be comparable to that of the samples used for investigation of electric effects in superconductors.

For pure crystals of elemental metals the total charge density is accessible trough Hall-effect measurements in the normal phase, but for dirty superconductors only the London-Hall effect gives such a possibility. Determination of the Hall constant by the current-induced contact-potential difference provides a tool for determination of the superfluid density at zero temperature $`n(0),`$ $`en_{\mathrm{tot}}=e^{}n(0)=1/_\mathrm{H}.`$ Thus knowing $`n(0)`$ and $`\lambda (0)`$ one can determine the effective mass of the Cooper pairs $`M_{\mathrm{CP}}`$:

$$\frac{1}{\lambda ^2(0)}=\frac{n(0)e^2}{M_{\mathrm{CP}}c^2\epsilon _0},M_{\mathrm{CP}}=\frac{e^{}\lambda ^2(0)}{c^2\epsilon _0_\mathrm{H}}.$$
(20)
All the three important parameters, $`_\mathrm{H}`$, $`M_{\mathrm{CP}},`$ and $`\lambda (T),`$ can be determined by measuring the contact potential difference in thin ($`d_{\mathrm{film}}\lambda (T)`$) and thick ($`d_{\mathrm{film}}\lambda (T)`$) superconducting films. We should mention that strong enhancement of the effective mass due to disorder is valid for arbitrary band anisotropy if the averaged gap does not vanish:

$$(\widehat{\xi }_{\mathrm{dirty}}^2)_{\alpha \beta }\frac{\tau }{T_\mathrm{c}}\frac{\chi ^2}{\chi ^2}v_\alpha v_\beta ,(x1).$$
(21)
As it was commented earlier,? the effective mass of the Cooper pairs was first measured by Fiory et al.? using electrostatic doping of YBa2Cu3O7-δ films. This pioneering experiment confirmed that the effective mass shows weak temperature dependence. Certainly, there are other methods? for absolute determination of $`M_{\mathrm{CP}},`$ e.g., from the surface Hall effect.? For the layered cuprates the problem of determining the Cooper pair mass is related to problem of the vortex charge.? The aim of the model experiment suggested in Ref. 15 was to check whether the Bernoulli potential of the pancake vortices creates an electrical polarization necessary to explain the Hall effect in the vortex creep regime. The supposed experimental setup is sketched in Fig. 3. Up to now the is only one experimental hint concerning the electrostatic polarization around the vortices, the NMR and NQR study by Kumagai, Nozaki and Matsuda.?

Finally we wish to mention that electric polarization of a metal-oxide-superconductor plane capacitor by a sudden light impulse can give a the condensation energy

$$\frac{\mathrm{\Delta }\phi }{_\mathrm{H}}=\frac{B_c^2(T)}{2\mu _0}.$$
(22)
We have to take into account that electric polarization is related to work function, workfunction is related to the chemical potential which is the Gibbs free energy per particle and in such a way we have an electric method for measurement of the thermodynamic critical magnetic field $`B_c(T).`$

## 5 What has to be done. A short working program for the experiment

After the analysis of the theory let us describe what has to be done experimentally for creation of Cooper pair mass spectroscopy based on the Bernoulli effect. Among the other electric field effect in superconductors such as surface Hall effect and electrostatic modulation of the kinetic inductance the Bernoulli effect looks simplest and experimental research is better to start with him. The main advantage is that in this case the electric field is very small and there are no restrictions related to the brake-trough voltage of the insulating layer. It is necessary insulator to be good only for small voltages. These investigations can be a by-product of preliminary research of the superconductor insulator interfaces in further trial to make superconducting field effect transistors. The first step will be qualitatively observation of the Bernoulli effect in already existing superconductor structures. A numerical calculation of the current distribution in the framework of London electrodynamics can convert this observation in quantitative measurement. The simplest possible set-up is a superconducting nanobridge in which an appropriate gate can be performed. In this field effect transistor type structure? current in the strip can be excited by capacitor connections id source and drain area. In such a way generation of the current harmonics by the contacts can be avoided.

Later on a systematic investigation can be started. The main purpose is to try to use already prepared big number of superconductor films without any destroying patterning to be carried out. The gold Bernoulli electrodes should be evaporated on a thin SrTiO3 substrate. The investigated film should be affixed on this plate. A small coil, or a solenoid with diameter of order of 1 mm, for example, will excite the eddy currents. The detector circuit should be switched between the circular and ring electrodes patterned on the insulator plate, cf. Fig. 2. The cross talk between detector and drive coil should be annulated by a variable mutual inductance. The harmonic current should be applied to the drive coil and the Bernoulli signal will be detected at second harmonics using a low noise preamplifier before the loch-in. A resonance technique additionally can use the high-Q factor to suppres the noise and signal at basic frequency.

In parallel a numerical simulation of the experimental setup will allow the result of measurements to be delivered as material constants. After a systematic investigation of the collected films would be possible to extract fundamental information of the doping dependence of the effective mass of Cooper pairs in cuprate CuO2 superconductors and the almost temperature independent London-Hall coefficient.

Finally the experimental setup could be elaborated as a commercial device (superfluid densitometer) intended for fast contactless investigation of the Bernoulli effect in superconductors. Such a device can find application for the monitoring of the quality of the thin superconducting films used in second-generation superconducting cables, for example. The investigation of other electric field effects in the superconductors will be the next step of the investigation of electric field effect in superconductors.

The measurement of the Bernoulli signal in irradiated thin films can become a routine tool for investigation of these disordered superconductors.

Except for the Cooper pair mass spectroscopy and charge carriers densitometry Bernoulli effect, i.e. measurement of the contact potential difference can become unique tool for determination of the condensation energy especially for type-II superconductors. For sudden heating we can use laser light illumination of the set-up for the Bernoulli measurements.

## 6 Temperature dependence of the penetration depth and optical mass

Let us attempt now a simple derivation of the temperature dependence of the penetration depth for clean superconductors. Our objective is to show that the result is in agreement with the formula for the coherence length, Eq. (1).

For a superconductor at $`T=0`$, consider the current response to a small space-homogeneous variation of the vector potential $`𝐀(t).`$ The Fermi surface is shifted in the momentum space as a rigid object; see, e.g., Ref. 6. In this respect there is no principal difference between a superconductor and a normal metal at optical frequencies. The difference is only that the optical approach, $`\omega \tau 1,`$ can be used for the static response of the clean superconductor. The electromagnetic field brings about a small change of the electron momentum $`𝐐=e𝐀(t)`$ and the shift in the momentum space, $`𝐩𝐩+𝐐,`$ creates in turn a small shift of the electron kinetic energies:

$$\epsilon _𝐩\epsilon _{𝐩+𝐐}\epsilon _𝐩+𝐯_𝐩𝐐+\frac{1}{2}𝐐m_𝐩^1𝐐,$$
(23)
$$m_𝐩^1=\frac{𝐯_𝐩}{𝐩}=\frac{^2\epsilon _𝐩}{𝐩^2}.$$
(24)
Thus for the increase of the kinetic energy density we obtain:

$`w_{\text{kin}}`$ $`={\displaystyle \frac{1}{2}}𝐐m^1𝐐n_c={\displaystyle \frac{1}{2}}𝐯_{\mathrm{dr}}m𝐯_{\mathrm{dr}}n_c,`$ (25)
$`={\displaystyle \frac{1}{2}}\epsilon _0c^2{\displaystyle \underset{\alpha ,\beta =1}{\overset{3}{}}}A_\alpha (\widehat{\lambda }^2(0))_{\alpha \beta }A_\beta ,`$
where $`𝐣=en_c𝐯_{\mathrm{dr}},`$ $`𝐯_{\mathrm{dr}}=m^1𝐐,`$ and

$$n_c=2_{\epsilon _𝐩<E_\mathrm{F}}\frac{d𝐩}{(2\pi \mathrm{})^3}$$

$$m^1=\frac{{\displaystyle \frac{d𝐩}{(2\pi \mathrm{})^3}m_𝐩^1\theta (E_\mathrm{F}\epsilon _𝐩)}}{{\displaystyle \frac{d𝐩}{(2\pi \mathrm{})^3}\theta (E_\mathrm{F}\epsilon _𝐩)}}=\frac{{\displaystyle _{\epsilon _𝐩=E_\mathrm{F}}}{\displaystyle \frac{dS_𝐩}{(2\pi \mathrm{})^3v_𝐩}}𝐯_𝐩𝐯_𝐩}{{\displaystyle _{\epsilon _𝐩<E_\mathrm{F}}}{\displaystyle \frac{d𝐩}{(2\pi \mathrm{})^3}}}.$$
(26)
Whence the expressions for the penetration depth and the optical mass tensor read

$`(\widehat{\lambda }^2(0))_{\alpha \beta }`$ $`={\displaystyle \frac{e^2n_c}{\epsilon _0c^2}}(m^1)_{\alpha \beta }={\displaystyle \frac{e^2}{\epsilon _0c^2}}2\nu _\mathrm{F}v_\alpha v_\beta ,`$ (27)
$`(m^1)_{\alpha \beta }`$ $`={\displaystyle \frac{2\nu _\mathrm{F}}{n_c}}v_\alpha v_\beta ,`$ (28)
$`\nu _\mathrm{F}`$ $`={\displaystyle \frac{d𝐩}{(2\pi \mathrm{})^3}\delta (E_\mathrm{F}\epsilon _𝐩)}.`$ (29)
The comparison with the Bernoulli effect considered in Sec. 3.1 tells us that the effective mass of the Cooper pairs in the clean limit is just twice the optical mass:

$$M_{\mathrm{CP}}=2m,n(T=0)=\frac{1}{2}n_c,$$
(30)
and neglecting some subtleties, e.g., the appearance of hole pockets, the total number of electrons corresponds to the total number of charge carriers $`n_c`$ and $`n_{\mathrm{tot}}.`$

Formally, the effect of nonzero temperature reduces? to introducing an additional multiplier $`r_d`$ in the averaging of the velocity-velocity tensor at the Fermi surface, i.e.

$$(\widehat{\lambda }^2(T))_{\alpha \beta }=\frac{e^2}{\epsilon _0c^2}2\nu _\mathrm{F}v_\alpha v_\beta r_d\left(\frac{\mathrm{\Delta }_𝐩}{2\pi k__\mathrm{B}T}\right),$$
(31)
where

$$r_d\left(\frac{\mathrm{\Delta }_𝐩}{2k__\mathrm{B}T}\right)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{2\pi k__\mathrm{B}T\mathrm{\Delta }_𝐩^2}{\left\{\mathrm{\Delta }_𝐩^2+\epsilon _n^2\right\}^{3/2}}\{\begin{array}{cc}1,\hfill & \text{ for }k__\mathrm{B}T|\mathrm{\Delta }_𝐩|\hfill \\ 7\zeta (3)\left(\frac{\mathrm{\Delta }_𝐩}{2\pi k__\mathrm{B}T}\right)^2,\hfill & \text{ for }|\mathrm{\Delta }_𝐩|k__\mathrm{B}T\hfill \end{array},$$

where $`\epsilon _n=(2n+1)\pi k__\mathrm{B}T.`$ The BCS multiplier $`r_d`$ can be easily derived using Matsubara Green function method?, averaging over the Fermi surface also can be easily derived using this technique. For superconductors with strong coupling corrections we suggest an interpolation formula which can be used before a detailed theory to be developed. Having the solution of the Eliashberg equations $`Z_{n,𝐩}`$, $`\stackrel{~}{\epsilon }_{n,𝐩}`$ and $`\stackrel{~}{\mathrm{\Delta }}_{n,𝐩}`$ we can perform in Eq. (6) the substitution $`\epsilon _n\stackrel{~}{\epsilon }_{n,𝐩}=Z_{n,𝐩}\epsilon _n`$ and $`\mathrm{\Delta }_{n,𝐩}\stackrel{~}{\mathrm{\Delta }}_{n,𝐩}`$. The procedure used by Kogan? is related to further substitution in Eq. (6) $`\stackrel{~}{\epsilon }_{n,𝐩}\stackrel{~}{\epsilon }_{n,𝐩}^{}=\stackrel{~}{\epsilon }_{n,𝐩}+1/\tau _𝐩`$. The function

$$r_d(z)=z^2\underset{n=0}{\overset{\mathrm{}}{}}\left[z^2+\pi ^2\left(n+1/2\right)^2\right]^{3/2}\{\begin{array}{cc}7\zeta (3)z^2/\pi ^2,\hfill & z1\hfill \\ 1,\hfill & z1\hfill \end{array},$$
(32)
gives implicitly the BCS expression? for the current response (for a simple derivation using temperature Green’s functions see, for example, Ref. 3)

For multiband superconductors, like MgB2, the general formula for the penetration depth reads

$`(\widehat{\lambda }^2(T))_{\alpha \beta }`$ $`=2\nu _\mathrm{F}{\displaystyle \frac{e^2}{\epsilon _0c^2}}{\displaystyle \underset{b}{}}c_br_d\left({\displaystyle \frac{\mathrm{\Delta }_b(T)}{2\pi k__\mathrm{B}T}}\right)v_\alpha v_\beta _b,`$
$`={\displaystyle \underset{b}{}}c_b(\widehat{\lambda }_b^2(0))_{\alpha \beta }r_d\left({\displaystyle \frac{\mathrm{\Delta }_b(T)}{2\pi k__\mathrm{B}T}}\right),`$ (33)
where $`\lambda _b^2(0)`$ is the contribution of the $`b`$th band at $`T=0.`$ A careful fit of $`\lambda (T)`$ for MgB2 could make it possible to separate the influence of the $`\pi `$\- and $`\sigma `$-bands.

Returning to the expression for the kinetic energy density, Eq. (25), let us introduce the momentum of the pairs $`𝚷=2𝐐,`$ and the small-gap approximation

$$r_d(𝐩)7\zeta (3)(|\mathrm{\Xi }|\chi _𝐩/2\pi k__\mathrm{B}T)^2.$$

Then for a clean superconductor one obtains

$`w_{\mathrm{kin}}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle _{\epsilon _𝐩<E_\mathrm{F}}}{\displaystyle \frac{2d𝐩}{(2\pi \mathrm{})^3}}𝐐m_𝐩^1𝐐r_d(𝐩)`$ (34)
$`={\displaystyle \frac{1}{2}}(𝚷\mathrm{\Xi })^{}\stackrel{~}{M}^1(𝚷\mathrm{\Xi })=a_0[𝐃\mathrm{\Xi }(𝐫)]^{}\widehat{\xi }^2𝐃\mathrm{\Xi }(𝐫)`$
$`={\displaystyle \frac{7\zeta (3)\mathrm{}^2\nu _\mathrm{F}}{\left(4\pi k__\mathrm{B}T_\mathrm{c}\right)^2}}[𝐃\mathrm{\Xi }(𝐫)]^{}𝐯𝐯\chi ^2𝐃\mathrm{\Xi }(𝐫),`$
where we have used the replacement

$$𝚷\mathrm{\Xi }i\mathrm{}𝐃\mathrm{\Xi }(𝐫)=\left(i\mathrm{}\frac{}{𝐫}e^{}𝐀(𝐫)\right)\mathrm{\Xi }(𝐫).$$
(35)
Equation (34) agrees with Eq. (1) for $`x_c=0`$ which is perhaps the simplest validation of the gradient terms in the GL expansion. In this way, for temperatures slightly below $`T_\mathrm{c}`$ we obtain

$$\stackrel{~}{M}_{\mathrm{clean}}^1=\frac{7\zeta (3)\nu _\mathrm{F}}{8\pi ^2(k__\mathrm{B}T_\mathrm{c})^2}𝐯𝐯\chi ^2,$$
(36)
$$\widehat{\lambda }_{\mathrm{clean}}^2(T)\frac{4e^2\nu _\mathrm{F}}{\epsilon _0c^2}\frac{\chi ^2𝐯𝐯\chi ^2}{\chi ^4}\frac{T_\mathrm{c}T}{T_\mathrm{c}}.$$
(37)
For isotropic-gap superconductors, $`\chi =\mathrm{const},`$ and arbitrary band anisotropy this equation together with Eq. (27) gives

$$\widehat{\lambda }_{\mathrm{clean}}^2(0)=\frac{1}{2}\frac{d}{dT}\widehat{\lambda }_{\mathrm{clean}}^2(T)|_{T_\mathrm{c}}.$$
(38)
## 7 Conclusions

We demonstrated that temperature dependence of the penetration depth? in clean superconductors can be presented and an inserting of an additional BCS factor $`r_d`$ in the formula for far infrared skin depth of the metal. We conclude that London-Hall coefficient is temperature- and disorder independent, and the Bernoulli effect is an adequate method for determination of $`_\mathrm{H}.`$ Given the penetration depth and the London-Hall coefficient at $`T=0`$ we can determine the effective mass of the superfluid particles. Measuring the Bernoulli effect for thin and thick films, see Eq. (16) and Eq. (17) we conclude that London penetration depth $`\lambda (T)`$ also can determined by purely electrostatic measurements. This hydrodynamic mass $`M_{\mathrm{CP}}`$ which parameterizes the relation between $`_\mathrm{H}`$ and $`\lambda (0)`$ is in the clean limit the optical mass extrapolated to zero frequency. But for disordered superconductors, irradiated thin films, for example, investigation of the disorder dependence of the Bernoulli signal can give an important information for the disorder renormalization of the effective Cooper pair mass. On Fig. 1 is presented the result for an isotropic gap superconductor but for a different gap anisotropies the curves will be different. In such a way we are coming to the conclusion that the impurity and temperature dependence of the Bernoulli signal is an important property which can be used to check the comparison between the theory and the experiment. For high-$`T_\mathrm{c}`$ superconductors the Bernoulli effect gives the possibility to evaluate the vortex charge.? We emphasize that Bernoulli effect and other electric field effects in superconductors, e.g., the surface Hall effect,? gives a complete set of parameters describing the low-frequency current response of superconductors. We conclude also that bulk condensation energy of a superconductor can also be determined by electrostatic polarization created by a sudden heating due to Eq. (22).

Acknowledgments

This work was partially supported by the Flemish Government Programme GOA.

References