text:
# Quasiparticle Decoherence in 𝑑-wave Superconducting Qubits

## Abstract

It is usually argued that the presence of gapless quasiparticle excitations at the nodes of the $`d`$-wave superconducting gap should strongly decohere the quantum states of a $`d`$-wave qubit, making quantum effects practically unobservable. Using a self-consistent linear response non-equilibrium quasiclassical formalism, we show that this is not necessarily true. We find quasiparticle conductance of a $`d`$-wave grain boundary junction to be strongly phase dependent. Midgap states as well as nodal quasiparticles contribute to the conductance and therefore decoherence. Quantum behavior is estimated to be detectable in a qubit containing a $`d`$-wave junction with appropriate parameters.

Among numerous qubit implementations, superconducting ones enjoy long decoherence times because of their gapped electronic excitation spectrum. This fact has recently been confirmed by several striking experiments expt . The key constituents in all of those are Josephson tunnel junctions. One advantage of the tunnel junctions is the exponential dependence of their quasiparticle resistance $`R`$ on temperature $`T`$ (i.e. $`Re^{\mathrm{\Delta }/T}`$ where $`\mathrm{\Delta }`$ is the superconducting gap. Herein $`k_B=\mathrm{}=1`$). The electronic decoherence is therefore exponentially suppressed at low $`T`$. Similar behavior also exists in superconducting point contacts, with the energy of Andreev levels $`ϵ_0(\varphi )=\mathrm{\Delta }\mathrm{cos}\varphi /2`$, replacing $`\mathrm{\Delta }`$ in the exponent LY ; new . Deviation from the exponential dependence is expected at low temperatures new .

Despite their naturally degenerate ground states ilichev01 , desirable for quantum computation, $`d`$-wave qubits dwqubit ; SPQ  are controversial because their quasiparticle spectrum is gapless at the nodes of the order parameter. Moreover, experimentally, the normal resistance extracted from I-V characteristics of $`d`$-wave grain boundary junctions is found to be very small HM ; RC  (with $`RC1`$ps RC ); much smaller than required to observe quantum effects. However, the resistance is measured in the running state of the junctions. The Doppler shift due to large superconducting current in such a state will populate the nodes, enhancing the conductance $`G`$ ($`1/R`$). Moreover, time variation of the phase difference across the junction would effectively phase average $`G`$. As we shall see, midgap states (MGS) can significantly contribute to such an averaged $`G`$, except at very low $`T`$.

The first step is therefore to calculate $`G`$ for a $`d`$-wave grain boundary junction. Most existing methods can study ac properties of a Josephson junction biased with a constant voltage (see e.g. hurd ; zaikin ). This however implies a constant variation of phase difference $`\varphi `$, which is not the case in qubits. In Ref. new, , a self-consistent non-equilibrium quasiclassical technique was developed to calculate linear response of a Josephson junction to an ac voltage with arbitrary frequency. The method was successfully applied to the case of a superconducting point contact. In high $`T_c`$ superconductors, the relatively large $`T_c/E_F`$ ($`E_F`$ is the Fermi energy), makes the quasiclassical approximation only marginally applicable. It nevertheless has proven successful in calculating equilibrium properties of $`d`$-wave superconductors dwave ; amin ; MS . Here, we employ the theory of Ref. new,  to calculate $`G`$, and therefore decoherence, in a $`d`$-wave grain boundary junction.

Let us now briefly describe the technique. More details are provided in Ref. new, . We calculate quasiclassical retarded, advanced, and Keldysh Green’s functions new ; eschrig

$`\widehat{g}^{R,A}=\left(\begin{array}{cc}g^{R,A}& f^{R,A}\\ f^{R,A}& g^{R,A}\end{array}\right),\widehat{g}^K=\left(\begin{array}{cc}g^K& f^K\\ f^K& g^K\end{array}\right),`$ (5)
which are functions of the Fermi velocity $`𝐯_F`$, quasiparticle energy $`ϵ`$, position $`𝐫`$, and time $`t`$. Here $`f^R(𝐯_F,ϵ;𝐫,t)f^R(𝐯_F,ϵ;𝐫,t)^{}`$, etc. In equilibrium, the retarded and advanced Green’s functions can be written in terms of Riccati amplitudes MS  $`a_0^\alpha `$ and $`b_0^\alpha `$ in a way very similar to the conventional method for the Matsubara Green’s functions MS ; amin :

$$g_0^\alpha =s^\alpha \frac{1a_0^\alpha b_0^\alpha }{1+a_0^\alpha b_0^\alpha },f_0^\alpha =s^\alpha \frac{2a_0^\alpha }{1+a_0^\alpha b_0^\alpha },$$
(6)
where $`\alpha =R,A`$ for retarded and advanced functions respectively, and $`s^\alpha =+`$ ($``$) for $`\alpha =R`$ ($`A`$). The subscript “0” denotes equilibrium quantities. The amplitudes satisfy Riccati-type equations

$`𝐯_Fa_0^\alpha `$ $`=`$ $`2iϵ^\alpha a_0^\alpha (a_0^\alpha )^2\mathrm{\Delta }_0^{}+\mathrm{\Delta }_0,`$
$`𝐯_Fb_0^\alpha `$ $`=`$ $`2iϵ^\alpha b_0^\alpha (b_0^\alpha )^2\mathrm{\Delta }_0+\mathrm{\Delta }_0^{},`$ (7)
where $`ϵ^\alpha =ϵ+is^\alpha \eta `$, with $`ϵ`$ and $`\eta `$ being the real and imaginary parts of the quasiparticle energy respectively. In $`d`$-wave superconductors, $`\eta `$ results from both inelastic and impurity scattering processes eta . The boundary conditions are the bulk solutions of (7): $`a_0^\alpha =\mathrm{\Delta }_0/(iϵ^\alpha +s^\alpha \mathrm{\Omega }^\alpha )`$ and $`b_0^\alpha =\mathrm{\Delta }_0^{}/(iϵ^\alpha +s^\alpha \mathrm{\Omega }^\alpha )`$, where $`\mathrm{\Omega }^\alpha =\sqrt{|\mathrm{\Delta }_0|^2(ϵ^\alpha )^2}`$. To calculate $`a_0^R`$ and $`b_0^A`$ ($`b_0^R`$ and $`a_0^A`$), we integrate (7) in the direction of $`𝐯_F`$ ($`𝐯_F`$), starting from the boundary conditions at $`x=\mathrm{}`$ ($`+\mathrm{}`$).

Non-equilibrium behavior emerges when an ac voltage $`V(t)=V_0\mathrm{cos}\omega t`$ is applied across the junction. Choosing a gauge in which the vector potential $`𝐀=0`$, we consider a scalar potential $`\mathrm{\Phi }=\pm (V_0/2)\mathrm{cos}\omega t`$, with the $`+()`$ sign on the left (right) side of the junction. To ensure applicability of the linear response formalism, it is necessary that $`eV_0\mathrm{\Delta },\omega `$ note .

We define the linear response Green’s functions $`\delta \widehat{g}^\alpha =\widehat{g}^\alpha \widehat{g}_0^\alpha `$, related to the linear response amplitudes $`\delta a^\alpha =a^\alpha a_0^\alpha `$ and $`\delta b^\alpha =b^\alpha b_0^\alpha `$ through ($`\alpha =R,A`$)

$`\delta g^\alpha `$ $`=`$ $`2s^\alpha {\displaystyle \frac{\delta a^\alpha b_0^\alpha +\delta b^\alpha a_{0+}^\alpha }{(1+a_{0+}^\alpha b_{0+}^\alpha )(1+a_0^\alpha b_0^\alpha )}},`$ (8)
$`\delta f^\alpha `$ $`=`$ $`2s^\alpha {\displaystyle \frac{\delta a^\alpha \delta b^\alpha a_{0+}^\alpha a_0^\alpha }{(1+a_{0+}^\alpha b_{0+}^\alpha )(1+a_0^\alpha b_0^\alpha )}},`$ (9)
where $`a_{0\pm }^\alpha a_0^\alpha (ϵ\pm \omega /2)`$, etc. We also introduce the anomalous Green’s function $`\delta \widehat{g}^X`$ (with the same matrix form as $`\widehat{g}^R`$) by $`\delta \widehat{g}^K=\delta \widehat{g}^X(_+_{})+\delta \widehat{g}^R_{}\delta \widehat{g}^A_+`$, where $`_\pm \mathrm{tanh}[(ϵ\pm \omega /2)/2T]`$. Correspondingly, we introduce anomalous functions $`\delta a^X`$ and $`\delta b^X`$ which are related to the Green’s functions through

$`\delta g^X`$ $`=`$ $`2{\displaystyle \frac{\delta a^X\delta b^Xa_{0+}^Rb_0^A}{(1+a_{0+}^Rb_{0+}^R)(1+a_0^Ab_0^A)}},`$ (10)
$`\delta f^X`$ $`=`$ $`2{\displaystyle \frac{\delta a^Xa_0^A+\delta b^Xa_{0+}^R}{(1+a_{0+}^Rb_{0+}^R)(1+a_0^Ab_0^A)}}.`$ (11)
The differential equations describing $`\delta a^\alpha `$ are $`𝐯_F\delta a^\alpha =A^\alpha \delta a^\alpha +B^\alpha `$, with bulk boundary conditions $`\delta a^\alpha =B^\alpha /A^\alpha `$, where

$`A^\alpha =\{\begin{array}{cc}2iϵ(a_{0+}^\alpha +a_0^\alpha )\mathrm{\Delta }_0^{},& \alpha =R,A\\ i\omega a_{0+}^R\mathrm{\Delta }_0^{}+b_0^A\mathrm{\Delta }_0,& \alpha =X\end{array}`$ (14)
$`B^\alpha =\{\begin{array}{cc}\delta \mathrm{\Delta }a_{0+}^\alpha a_0^\alpha \delta \mathrm{\Delta }^{}ie\mathrm{\Phi }(a_{0+}^\alpha a_0^\alpha ),& \alpha =R,A\\ a_{0+}^R\delta \mathrm{\Delta }^{}+b_0^A\delta \mathrm{\Delta }ie\mathrm{\Phi }(1+a_{0+}^Rb_0^A).& \alpha =X\end{array}`$ (17)
The corresponding equations for $`\delta b^\alpha `$ can be obtained by applying a $``$-operation to both sides of the above equations. Integrations for $`\delta a^R`$, $`\delta b^A`$, and $`\delta a^X`$ ($`\delta b^R`$, $`\delta a^A`$, and $`\delta b^X`$) are taken in the direction of $`𝐯_F`$ ($`𝐯_F`$), along the quasiclassical trajectories.

The equilibrium order parameter, $`\mathrm{\Delta }_0`$, is calculated self consistently using numerical iteration, while the linear response $`\delta \mathrm{\Delta }\mathrm{\Delta }\mathrm{\Delta }_0`$ is given by new

$$\delta \mathrm{\Delta }(\omega )=2e(\mathrm{\Delta }_0/\omega )\mathrm{\Phi }(\omega ).$$
(18)
$`\delta \mathrm{\Delta }`$ satisfies the self-consistency relation up to an error of O$`(j_c/j_{\mathrm{c},\mathrm{bulk}})`$, where $`j_c`$ ($`j_{\mathrm{c},\mathrm{bulk}}`$) is the Josephson (bulk) critical current density. In $`d`$-wave grain boundary junctions, $`j_c/j_{\mathrm{c},\mathrm{bulk}}`$ is smallest for a $`0^{}45^{}`$ junction. Application of (18), instead of common self-consistent iterative methods eschrig , is a clear advantage of the present technique.

We consider an asymmetric $`0^{}45^{}`$ misoriented $`d`$-wave grain boundary junction, with perfect transparency and no roughness. To clearly see the effect of the nodes and bound states, it is useful to define an angle resolved conductance $`G_\theta =\mathrm{Re}[\delta j_\theta /V_0]S`$, where $`S`$ is the the area of the junction and

$$\delta j_\theta =\frac{ev_FN_F}{4}_{ϵ_c}^{ϵ_c}𝑑ϵ\mathrm{Tr}\left[\widehat{\tau }_3\delta \widehat{g}^K\right]_{𝐯_F=v_F\widehat{𝐧}}$$
(19)
is the contribution to the ac current density from a trajectory in $`\theta `$ direction (cf. Fig. 1b). Here $`\widehat{\tau }_3`$ is the third Pauli matrix, $`\widehat{𝐧}`$ is a unit vector in the direction of the trajectory ($`\widehat{𝐧}\widehat{𝐱}=\mathrm{cos}\theta `$), $`N_F`$ is the density of states at the Fermi surface, and $`ϵ_c`$ is an energy cutoff, which in our numerical calculations is taken to be 10$`T_c`$. (The results however are insensitive to the exact value of $`ϵ_c`$.) The total conductance is given by $`G=(d\theta /2\pi )G_\theta \mathrm{cos}\theta `$.

Figure 1a shows $`G_\theta `$ for two values of phase difference $`\varphi `$ across the junction. At $`\varphi =0`$, $`G`$ is maximum at $`\theta =22.5^{}`$ and $`67.5^{}`$ (paths $`\alpha `$ and $`\beta `$ in Fig. 1b). These are the directions along which the quasiparticles see different signs of the $`d`$-wave order parameter on either side of the boundary. Zero energy Andreev bound states formed in these directions are responsible for the large $`G`$ new . Notice the asymmetry of $`G_\theta `$ with respect to $`\theta =0`$ at $`\varphi =0`$. This asymmetry leads to an ac current parallel to the boundary hurd . The current vanishes at $`\varphi =\pi /2`$ and is not dissipative (it is perpendicular to the electric field), therefore not central to our discussion.

A nonzero phase difference splits the bound states, moving them away from zero energy. Maximum splitting occurs at $`\varphi =\pi /2`$. Nodal directions dominate the quasiparticle current at this phase difference (cf. Fig. 1a). Notice orders of magnitude difference in the conductance maxima between the two phase differences.

Figure 2a displays $`\varphi `$-dependence of the total $`G`$ at two temperatures. Sharp peaks at $`\varphi =0`$ and $`\pi `$, 2–3 orders of magnitude larger than the minimum (at $`\varphi =\pi /2`$), result from the MGS. The conductance spikes are much narrower at $`T=0.01T_c`$ than at $`0.1T_c`$ and become extremely sharp at even lower $`T`$. As a result, the influence of the MGS is small at low $`T`$ (see below). Figure 2b shows the $`T`$-dependence of $`G`$ at $`\varphi =0`$ and $`\pi /2`$. While the two curves coincide at high $`T`$, they diverge as $`T`$ is lowered. The conductance at $`\varphi =\pi /2`$ is suppressed at low $`T`$. This is because the gapless quasiparticles on one side of the junction see a gapped order parameter on the other side bruder . As $`T0`$, however, it starts to saturate to a value almost proportional to $`\eta `$ (see Fig. 3b). Similar behavior was also found for the case of superconducting point contacts new .

To reduce decoherence, the qubit should be operated at low $`T`$. Figure 3 shows the $`\omega `$ and $`\eta `$ dependence of $`G`$ at $`T=10^4T_c`$. While strong $`\omega `$-dependence of $`G`$ exists at $`\varphi =0`$, it is almost frequency independent when $`\varphi =\pi /2`$ (Fig. 3a). The $`\eta `$-dependence of $`G`$ at $`\varphi =\pi /2`$, on the other hand, is close to linear (Fig. 3b). The positive slope can be attributed to the broadening of the MGS with $`\eta `$, which enhances their overlap near zero energy, increasing their contribution to $`G`$.

We now try to calculate the decoherence time in a $`d`$-wave qubit, assuming a double well potential $`U(\varphi )`$ with minima at $`\pm \varphi _0\pm \pi /2`$. (A practical example of such a system is given in Ref. SPQ, .) We first write down a Hamiltonian which, in classical regime, reproduces the above calculated $`G`$. At low frequencies, $`G`$ shows slow $`\omega `$-dependence for almost all $`\varphi `$ (except for $`\varphi =0,\pi `$, which do not contribute to decoherence—see below). One can therefore assume our system to be coupled to an Ohmic heat bath leggett . To get a $`\varphi `$-dependent $`G`$, however, it is necessary to consider nonlinear coupling to the bath caldeira ; smirnov  \[see Eq. (24) below\]. We therefore write the Hamiltonian as

$$H=\frac{Q^2}{2C}+U(\varphi )F(\varphi )X,$$
(20)
where $`Q=2eP_\varphi `$ ($`P_\varphi `$ momentum conjugate to $`\varphi `$) is the charge stored on the junction, and $`X`$ is a heat bath operator. To find the relation between the nonlinear coupling function $`F`$ and $`G`$, we first study the classical behavior of the above Hamiltonian. The equation of motion is

$`(C/4e^2)\ddot{\varphi }=_\varphi U(\varphi )+_\varphi F(\varphi )X,`$ (21)
Assuming small back-action of the system on the heat bath, one can use linear response theory smirnov  to write

$`X(t)=X^{(0)}(t)+{\displaystyle 𝑑t^{}D(tt^{})F[\varphi (t^{})]},`$ (22)
where $`X^{(0)}`$ is the unperturbed bath operator, $`D(tt^{})=i[X^{(0)}(t),X^{(0)}(t^{})]_{}\theta (tt^{})`$ is the bath’s retarded Green’s function, and $`\mathrm{}`$ denotes equilibrium statistical averaging. For an Ohmic heat bath $`D(t)=\alpha \delta ^{}(t)`$ \[in Fourier space $`D(\omega )=i\alpha \omega `$\], with $`\alpha `$ being the dissipation coefficient and $`\delta ^{}(t)`$ defined by $`𝑑t^{}\delta ^{}(tt^{})z(t^{})=\dot{z}(t)`$.

Eq. (21) therefore leads to a stochastic equation that resembles the classical Langevin equation

$`(C/4e^2)\ddot{\varphi }+\alpha \left(_\varphi F\right)^2\dot{\varphi }+_\varphi U=\xi _\varphi ,`$ (23)
where $`\xi _\varphi =(_\varphi F)X^{(0)}(t)`$ is the fluctuation force xi . The second term in (23) gives dissipation and therefore corresponds to the $`\varphi `$-dependent $`G`$ caldeira : $`G(\varphi )=4e^2\alpha \left(_\varphi F\right)^2.`$ As a result

$`F(\varphi )={\displaystyle \frac{1}{2e\sqrt{\alpha }}}{\displaystyle _0^\varphi }𝑑\varphi ^{}\sqrt{G(\varphi ^{})}.`$ (24)
Notice the nonlocal dependence of $`F`$ on $`G`$. The spike in $`G`$ at $`\varphi =0`$ (cf. Fig. 2a) becomes extremely sharp as $`T0`$, with little contribution to $`F(\varphi )`$. This, although easily justified numerically, can be understood by making an analogy to the case of a superconducting point contact new : $`G(\varphi )(\mathrm{\Delta }/T)e^{\mathrm{\Delta }|\varphi |/T}`$ near the spike ($`\varphi `$ measured from the center of the spike), leading to a contribution to $`F(\varphi )`$ proportional to $`\sqrt{T/\mathrm{\Delta }}`$, which is negligible as $`T0`$. Numerical integration shows an even smaller effect.

In the quantum regime, one can truncate the Hilbert space to the left and right degenerate states, $`|L,|R`$, of the double-well potential, with $`L|F(\varphi )|R0`$ and $`R|F(\varphi )|RL|F(\varphi )|LF(\varphi _0).`$ Equation (20) then gives the effective two-state Hamiltonian of the system (at the degeneracy), coupled to the heat bath, as

$$H_{\mathrm{eff}}=\frac{\delta E}{2}\sigma _xF(\varphi _0)X\sigma _z,$$
(25)
where $`\delta E`$ is the energy splitting between the two lowest energy states of the system. The dephasing rate is proportional to the spectrum of heat bath fluctuations, $`S(\omega )=\mathrm{Im}[D(\omega )]\mathrm{coth}(\omega /2T),`$ taken at the resonance frequency of the two-level system \[see, e.g., Eq. (3.11) in Ref. leggett, \]. Then, the dephasing time $`\tau _\phi `$ due to coupling to the Ohmic heat bath is obtained from

$`\tau _\phi ^1=\alpha F(\varphi _0)^2\delta E\mathrm{coth}{\displaystyle \frac{\delta E}{2T}}.`$ (26)
Numerical calculation at $`T=10^4T_c`$ gives $`_0^{\pi /2}\sqrt{G}𝑑\varphi 0.38\sqrt{G_0}`$, for $`\eta =0.01T_c`$. Substituting into (24) and (26) and restoring $`\mathrm{}`$, we find $`\tau _\phi ^10.023R_QG_0\delta E/\mathrm{}`$, where $`R_Q=h/(2e)^26.45`$ k$`\mathrm{\Omega }`$ is the quantum resistance.

For quantitative estimation of $`\tau _\phi `$, we need to know the value of $`G_0`$. We extract $`G_0`$ from the Josephson critical current density whose value is available from experiment HM : $`j_c10^210^4`$ A/cm2. Our calculations at $`T=0.05T_c`$ (close to 4.2K where most experiments are performed) show a critical current $`I_c0.08G_0T_c/e`$, almost independent of $`\eta `$. For a submicron junction of area $`S0.01\mu `$m2, we obtain (taking $`T_c100`$ K) $`G_010^510^3\mathrm{}`$. Assuming $`\delta E/h1`$ GHz, we find $`\tau _\varphi 1100`$ ns (qubit quality factor $`Q1100`$). Smaller $`\eta `$ will result in a larger $`\tau _\varphi `$. It is therefore desirable to use materials with low disorder and junctions with small roughness. Depending on the parameters, it is possible to observe signatures of quantum behavior (e.g. coherent tunneling), although the decoherence time may not be long enough for quantum computation.

Realistic junctions, in general, suffer from finite reflectivity, roughness, and faceting note3 . The effect of roughness, to some extent, is similar to that of $`\eta `$; it broadens the MGS, increasing their contribution to $`G`$ as they split golubov . Imperfect transparency also affects the MGS in a nontrivial way, influencing $`G`$. Presence of a subdominant order parameter at the junction, although yet unjustified experimentally NV , can enhance the quasiparticle resistance by opening a gap at the nodes. A small bulk size can also have a similar effect by quantizing momentum along the nodal directions. Most of the above mentioned effects can be studied within the framework of the present formalism and are the subject of further investigations.

Other sources, such as fluctuations of the external fields, coupling of bulk nodal quasiparticles to the electromagnetic field produced by the qubit (due to e.g. spontaneous currents), coupling to nuclear spins or paramagnetic impurities, background charge fluctuations, etc. will also contribute to the decoherence.

It is our pleasure to thank Ya. Fominov, A.A. Golubov, W.N. Hardy, Y. Imry, A. Maassen van den Brink, N. Schopohl, P.E.C. Stamp, and A.M. Zagoskin for helpful discussions, A.Ya. Tzalenchuk for experimental information, and A.J. Leggett for pointing out the necessity for nonlinear coupling to the heat bath.


format: text

text:
# A New Technique to improve RFI suppression in Radio Interferometers

## 1. Introduction

Radio interferometers provide high-resolution, high-dynamic range radio images of celestial objects. It is well known that radio interferometers are less susceptible to radio frequency interference (RFI) compared to single dish. However, the dynamic range of the images in many cases is limited by interference indicating the need to improve RFI suppression.

In this paper, we describe a new technique to improve RFI suppression in interferometers 111To our knowledge a technique similar to that discussed in this paper was used earlier in Mauritius Radio Telescope data reduction (Sachdev, Udayashankar 1995). . The basic idea is to replace the averaging done at the output of the correlator with a multirate filtering process. A brief discussion of RFI suppression in existing interferometers is given in Sec. 2. Sec. 3 describes the proposed technique, gives the design of the multirate filter and presents simulation results. The limitations of the proposed technique are discussed in Sec. 4 and the effects of the multirate filtering process in the visibility and spatial domain are discussed in Sec. 5.

## 2. RFI suppression in existing Radio Interferometers

### 2.1. Correlated and Uncorrelated RFI

In the case of interferometric observations interference can be broadly classified as correlated and uncorrelated RFI. Consider a two element interferometer near a RFI source (Fig. 1). The voltages $`V_1`$ and $`V_2`$ at the outputs of the two antennas each consists of three components — the system noise ($`V_{1sys}`$, $`V_{2sys}`$), the source noise222The source noise is that component of the sky noise which is not ‘resolved’. ($`V_{1source}`$, $`V_{2source}`$) and the RFI component ($`V_{1rfi}`$, $`V_{2rfi}`$). For the moment we will not consider fringe stopping, which is the case when the source is at the north pole. Also we consider that the RFI is not moving with respect to the interferometer. Interferometers usually use complex correlators, which are constituted of ‘real’ and ‘imaginary’ parts. For the discussion in this paper, we consider a ‘real’ correlator. But the analysis is equally applicable to the ‘imaginary’ part of the correlator as well. The mean value of the correlator output, $`V_c`$, is given by

$$V_c=<V_1V_2>=<V_{1source}V_{2source}>+<V_{1rfi}V_{2rfi}>,$$
(1)
where $`<>`$ denotes time average. If the separation between antennas (baseline length) is $`>`$ a few 100 kms, which is the case for the VLBI (very long baseline interferometry), and if the RFI source is near one of the antennas, then the second antenna will not pick up the RFI (ie say $`V_{1rfi}=`$ 0). This means that $`V_c`$ does not have any contribution from the RFI source or that the RFI is uncorrelated. If this condition is not true then the RFI is correlated (for this discussion we consider that the bandwidth of the system is very small so that there is no bandwidth decorrelation). In most connected interferometers RFI is correlated in many baselines. For the VLBI, most terrestrial RFI is uncorrelated, but satellite interference can be correlated in most ‘short’ baselines. The technique discussed in this paper to suppress interference deals with correlated RFI. Note that the RMS (root mean square) value of the correlator output is proportional to the geometric mean of the total powers at the inputs of the correlator if the source noise power is much smaller than the system noise power. The correlator input powers have contributions from the RFI even if it is uncorrelated. Thus uncorrelated RFI also limits the dynamic range of the image.

### 2.2. Fringe stopping and correlated RFI

For simplicity consider that the baseline of the two element interferometer is in the East-West direction. The output of the two antennas are connected to a correlator with cables of equal length (Fig. 2). The response of the interferometer $`P(\theta )`$ is given by (Thompson, Moran, Swenson, 2001)

$$P(\theta )=P_{ant}(\theta )\text{cos}\left(\frac{2\pi D\text{sin}(\theta )}{\lambda }\right),$$
(2)
where $`P_{ant}(\theta )`$ is the cross response of the two antennas, $`D`$ is the baseline length and $`\lambda `$ is the operating wavelength. $`\theta `$ is measured with respect to the zenith, which we equate to the hour angle ($`HA`$) in radians for simplicity. A source moving through the interferometer response produces “fringes” at the output of the correlator. The fringe frequency $`f_{fringe}`$ in the present case is given by,

$$f_{fringe}=\frac{D}{\lambda }\text{cos}(HA)\frac{\text{d}(HA)}{\text{d}t}\text{cos}(\delta ),$$
(3)
where $`\delta `$ is the declination of the source and $`\frac{\text{d}(HA)}{\text{d}t}`$ is in units of radians per sec. The general equation for the fringe frequency is (Thompson et al. 2001)

$$f_{fringe}=\frac{\text{d}(HA)}{\text{d}t}u\text{cos}(\delta ),$$
(4)
where $`u`$ is the projected baseline in the East-West direction in units of $`\lambda `$. The above equation shows that the fringe frequency due to a source depends on the projected baseline and earth’s rotation rate ($`=\frac{\text{d}(HA)}{\text{d}t}`$). Now consider a terrestrial, stationary, correlated RFI. Since the position of the RFI does not change with respect to the interferometer, there will be no fringe due to the interference.

For observations with the interferometer, the interferometer is ‘phased’ in the direction of the interested radio source through out the observations. This process of phasing the interferometer (which includes both delay compensation and fractional phase compensation) is called ‘fringe stopping’. Fringe stopping essentially makes the correlator output constant due to a (point) source at the phase center. Sources away from the phase center produce fringes with frequencies equal to the difference between $`f_{fringe}`$ of the source at the phase center and that of the source away from the phase center. Thus, in the above case, larger the difference between the $`HA`$ of a source and that at the phase center the larger will be its ‘differential’ fringe frequency. This also means that the correlated RFI now produces a fringe with a frequency equal to $`f_{fringe}`$ of the source at the phase center because of fringe stopping.

### 2.3. Fringe averaging and RFI suppression

Usually in interferometric observations, one is interested in mapping sources within the primary beam of the antenna. Typically the differential fringe frequencies due to these sources is much smaller than $`f_{fringe}`$ of the source at the phase center. The maximum differential fringe rate is given by

$$\mathrm{\Delta }f_{fringe}\frac{\text{d}(HA)}{\text{d}t}\frac{\mathrm{\Omega }_{antenna}}{\mathrm{\Omega }_{synth}}\frac{\text{d}(HA)}{\text{d}t}\frac{D}{D_{antenna}},$$
(5)
where $`\mathrm{\Omega }_{antenna}`$ and $`D_{antenna}`$ are the beam and diameter of each element of an interferometer, $`\mathrm{\Omega }_{synth}\lambda /D`$ is the synthesized beam, $`\lambda `$ is the operating wavelength and $`D`$ is the baseline length. Consider, for example, $`D`$ = 10 km and $`D_{antenna}`$ = 25 m. The differential fringe frequency is $``$ 0.005 Hz while $`f_{fringe}`$ of the source at the phase center is $``$ 3.5 Hz, for $`\lambda =`$ 0.21 m. The low differential fringe frequency allows one to average the output of the correlator for several seconds. Since the RFI produces a fringe (due to fringe stopping process) this averaging reduces the amplitude of the RFI at the correlator output thus making the interferometer less susceptible to RFI (see Perley 2002 for a quantitative estimate of RFI attenuation in Very Large Array due to fringe averaging). An additional suppression of broadband RFI occurs due to bandwidth decorrelation, which is a bonus over the former effect. We will not discuss the bandwidth decorrelation effect any further here (see Thompson et al. 2001).

## 3. A new technique to improve RFI suppression in Interferometers

The averaging at the correlator output can be thought of as a moving average filtering process. The response of such a filter is shown is Fig. 4. The stop-band attenuation (i.e. the attenuation of frequencies beyond the cutoff frequency) of the filter determines the amount of RFI suppression achieved in this process if the fringe due to RFI is beyond the cutoff frequency of the filter. As seen from Fig. 4, stop-band attenuation in the worst case is 13.5 dB. We propose that RFI suppression in interferometers can be improved by using a filter with high stop-band attenuation at the output of the correlator. In the following sections we show how such filters can be designed and present results of the simulation done to tests its performance.

### 3.1. Multirate filtering

Present day correlators are designed using digital techniques. The typical bandwidth of the interferometer is a few tens of MHz or more. Sampling the signals at Nyquist rate implies that the data rate at the output of the multiplier will be of the order of the bandwidth. As mentioned above the output of the correlator can be averaged for several secs. This means that the required cutoff frequency of the filter that replaces the simple averaging should be a few tenths of a Hz. Thus the normalized (normalized with half the Nyquist frequency of the baseband) cutoff frequency of the filter is $`<10^6`$. A direct realization of a digital filter is not practical in such situations. However, such filters are realizable using multirate filtering technique. This technique essentially implements the required filter after reducing the sampling frequency (see Fig. 3; Bellanger 2000).

The basic filter used for the multirate filtering design is a digital elliptic filter. Using Matlab, a digital elliptic filter of order 10 with normalized cutoff frequency of 0.1 and stopband attenuation of 100 dB is designed. The passband ripple is chosen as 0.1 dB. The elliptic filter is realized using IIR (infinite impulse response) filter sections and its response is shown in Fig. 3. The sampling rate reduction is done in two ways: (1) using polyphase filtering technique; (2) by low-pass filtering the data using a moving average filter and then resampling the data at half the number of points averaged. The second type of filter is economical and simpler to implement in hardware. For testing the performance of the filter the sampling rate is reduced by 100. Thus the effective cutoff frequency of the filter is $`10^3`$ and we believe that this could be scaled further down by appropriately resampling the data. The responses of the multirate filters are shown in Fig. 4.

### 3.2. Results of the simulation

A simulation was done in Matlab to determine the improvement in RFI rejection of a two element interferometer by implementing the multirate filtering process. Fig. 5 shows the block diagram of the interferometer implemented for the simulation. The fringe generator adds a fringe to the source noise. The antenna blocks are power combiners, combining the source noise, the two independent system noises and the RFI. All noises are Gaussian random variables. The response of the receiver is realized using the receiver filters. Fringe stopping is done before multiplying the two antenna outputs. The relative delay between the two antennas is made Zero for the simulation.

The designed multirate filters and the moving average filter are used as post-detection filters to study their relative performance. A sine wave with normalized frequency of 0.65 is used as RFI. The receiver bandwidth in normalized frequency is selected as 0.2 centered at 0.625. The receiver noise ($`T_{sys}`$) and source noise powers are made equal to 0.2 in arbitrary units. Fig. 6 shows the spectra of the two receiver outputs. The source noise is multiplied with a fringe of normalized frequency $`3\times 10^3`$. If we consider that this fringe freq should translate to say 3 Hz, then the Nyquist frequency will be 2.0 KHz. This choice of the fringe frequency is limited by the practicality of the simulation. The outputs of the correlator for different types of post-detection filters are shown in Fig. 7. The remarkable ($`>`$ 40 dB) improvement in the RFI suppression compared to the case when a moving average filter is used as post-detection filter is evident from this figure. Comparing the spectrum of the correlator output with a moving average post-detection filter and multirate post-detection filter shows that the spectral values below the cutoff freq of the filter are almost the same for both cases. This means that the amplitudes of the desired fringe frequencies are not affected by the multirate filtering process.

The effectiveness of the technique is also determined by estimating the spectral power of the RFI at the output of the correlator with different types of post-correlation filters. To estimate the spectral power of RFI a power spectrum of the correlator output is first made. The difference of the spectral values of the channels where the RFI is expected and the average of the values of the two nearby channels is then taken. If the RFI power is comparable to the noise fluctuations in the spectrum then the difference can also be negative. We assigned a value of zero for the RFI spectral power if the estimate is negative. Fig. 6 shows the spectral power of RFI vs the interference-to-noise ratio (INR) of the RFI at the output of the receiver. The INR is the ratio of the RFI power to the noise power (receiver noise + source noise) at the output of the receiver. The remarkable suppression of the spectral power of the RFI with a multirate filter in post-detection stage is clear from Fig. 6.

## 4. Limitations of the technique

The technique discussed in this paper to improve RFI suppression of interferometers is effective only for correlated RFI (see Sec. 2). Usually terrestrial RFI is correlated on baselines of length less than a few 100 kms. Correlated RFI produces fringes at the output of the interferometer. The fringe frequency depends on the baseline length as well as the source declination for a given hour angle. Effective suppression of RFI occurs if the fringe frequency is much above the cut-off frequency of the multirate post-detection filter (see Sec. 3). However, if the baselines are small ($`<`$ a few kms) then the fringe frequency will be small and it becomes difficult to make filters that can attenuate these frequencies. Thus the technique presented here gives good RFI rejection for intermediate ($``$ a few kms to a few 100 kms) baselines. The fringe frequency also becomes smaller for sources near the north celestial pole (declinations near 90), which means that this technique will not work efficiently for this case as well.

Often RFI will have a strong carrier and associated modulating frequency components. The modulation results in a wider spectrum of the RFI at the output of the correlator. Thus only those components below the cut-off frequency of the multirate filter will be attenuated.

The analysis and simulation presented in this paper consider that the RFI is not time variable. However, in reality RFI is time variable. The response of the multirate filter to transient RFI needs to be studied. This will be presented elsewhere.

## 5. Effects of multirate filtering on the visibility and spatial domain

Till now the effect of the multirate filtering technique is discussed in relation with the time series at the output of the correlator. Any processing of the correlator output has an equivalent effect on the visibility domain and correspondingly on the spatial domain (i.e. positions in sky). The effective position of a stationary RFI is the north celestial pole. Thus suppressing the RFI using the multirate filtering process is equivalent to introducing a spatial null at the north celestial pole. In fact, the spatial null is not confined to the north celestial pole alone but extends along the spatial domain where the differential fringe frequencies are larger than the cutoff frequency of the multirate filter. The attenuation due to the spatial null will be proportional to the attenuation of the RFI due to the multirate filter.

In the visibility domain, the multirate filtering process essentially convolves the true visibility with the filter response along $`uv`$ tracks traced by pairs of antennas in an interferometer. Depending on the filter response, the amplitude and phase of the spatial frequencies of interest are modified. Since the filter is implemented digitally its response is known well and is identical for all the correlators. With the knowledge of this filter response, we believe a correction similar to the primary beam correction can be done to take out the effect of filtering. A detailed discussion on this will be presented elsewhere.

## 6. Acknowledgment

DAR thanks Rick Fisher for all useful discussions during the course of this work. The basic idea described in this paper stemmed from a discussion DAR had with Prof. A. Pramesh Rao, National Center For Radio Astrophysics, Pune, India.

## References

Bellanger, M., 2000, Digital Processing of Signals, John Wiley & Sons, Ltd, Chichester, UK.

Perley, R., 2002, EVLA Memo No. 49, NRAO Library, Socorro.

Sachdev, S., Udayashankar, N., 1995, BASI, 23, 574.

Thompson, A. R., Moran, J. M., Swenson, Jr., G. W., 2001, Interferometry and Synthesis Technique, John Wiley & sons, inc., New York.

## Discussion

B. Carlson: Could notch filters be used instead?

Anish Roshi: No, using notch filter instead of low-pass filter to suppress the fringe due to RFI has two disadvantages. (a) Since the fringe frequency is a function of baseline length the characteristics of the required notch filter for each baseline is different. This difference in filter characteristics can introduce baseline based effects thus may increase closure errors. (b) Man made RFI in general are modulated carriers. The modulation results in amplitude modulation of the fringe due to the interference. Using a low-pass filter, instead of notch filter, suppresses most of the sidebands of the ‘modulated fringe’ thus suppressing the RFI better.

D. Ferris: Replacing the normal rectangular integration function with the impulse response of your filter can be expected to increase the variance of the data, as well as smearing along time tracks in the u-v plane. Do you have estimates of these effects for the filters you presented?

Anish Roshi: Regarding variance, there will not be an increase in variance due to multirate filtering if the noise equivalent bandwidth of the filter is smaller than that of the conventional averaging filter. Regarding the smearing in the u-v plane: in the conventional analysis (see Thompson et. al 2001, p205) of the smearing of visibility in the u-v plane due to averaging, the averaging time is determined such that the u-v track in the u-v plane should not cross the cell size in the u-v plane, which is determined by the maximum field-of-view to be mapped. The condition thus obtained on the integration time is identical to that given by Eq. (5) in the text. The integration time is just the inverse of the maximum differential fringe frequency in Eq (5). Thus, as long as the multirate filter retains the maximum differential fringe frequency given by Eq. (5), there will not be any smearing in u-v plane for those sources within the field-of-view.


text:
# Importance of Perturbed Gravitational Potentials in Differentially Rotating Newtonian Stars

## 1 Introduction

Concerning the oscillations of spherical stars, Cowling (1941) started to use an approximation that neglect the perturbed gravitational potential to rewrite the system of equations into easier form. By using this approximation (later it is called the Cowling approximation), he found that there are three types of oscillations; the high frequency acoustic mode (p-mode), the low frequency gravity mode (g-mode), and the intermediate fundamental mode (f-mode).

In the last two decades, the Cowling approximation has been applied to general relativistic compact stars. In spherical stars, McDermott, Van Horn, and Scholl (1983) proposed the scheme in which all components of the metric perturbations are neglected. On the other hand, Finn (1988) proposed the formulation which neglects most components of the metric perturbations except for the $`(r,t)`$ component. The accuracies of these two approaches were compared by Lindblom and Splinter (1990) for mainly p-modes, and it was found that generally the approximation neglecting all components of the metric perturbations is better than the other one, while there is a possibility that the latter becomes better when the g-modes are considered. Hereafter, “the Cowling approximation” means the one that neglects all components of the metric perturbations.

It is also shown that the Cowling approximation is reliable in slowly rotating general relativistic stars (Kojima 1997). Furthermore, Ipser and Lindblom (1992) proposed the formulation of applying the Cowling approximation to rapidly rotating relativistic stars. The accuracy of the Cowling approximation in rotating stars was studied by Yoshida and Kojima (1997). They compared the two sets of solutions obtained by using the procedure produced by Ipser and Lindblom (1992) and by solving the exact equations for the f- and p-mode oscillations of nearly spherical stars with slow rotation. According to their results, solutions of the Cowling approximation agree well with the exact solutions for slowly rotating and highly relativistic models. Recently, Yoshida and coworkers have used the Cowling approximation to find neutral points of the f-mode oscillations of rapidly rotating polytropes (Yoshida and Eriguchi 1997, 1999, 2001; Yoshida, Rezzolla, Karino, and Eriguchi 2002). Also they checked the accuracy of the approximation in relativistic stars and found the conditions that the approximation becomes accurate (Yoshida and Eriguchi 1997, 1999).

Our knowledges about the Cowling approximation to the oscillations of rotating stars thus far can be summarized as follows:
The Cowling approximation becomes more accurate
(1)when the compressibility of the star is sufficiently high,
(2)when the gravity of the star is sufficiently strong.
It may be natural that the Cowling approximation becomes more accurate if the effect of the gravity does not work effectively to the fluid mass. For the simple oscillations such as f-mode, the amplitude of the mode becomes a simple monotonically increasing function of a stellar radius, and takes its maximum at the surface. Here, when we consider the density distribution of a compressible star, the density of surface region becomes relatively low. In other words, there are only small masses in the region where the mode amplitude gets large, in the centrally condensed star. Also, effectively, the strong gravitation in general relativity (GR) can have the same effect. Consequently, the effects of compressibility and strong gravity make the accuracy of the Cowling approximation better.
and
(3)when the rotation of the star is sufficiently rapid, when the star rotates rigidly.
Here, notice that the effect of differential rotation on the approximation has not studied, therefore we will focus on it.

From the astrophysical points of view, we should take additional effects into account. For example, some stars are rotating differentially. In particular, newly born neutron stars may rotate differentially (see e.g. Liu and Lindblom 2001, Shibata and Uryu 2000, Dimmelmeier, Font and Müller 2002). Young neutron stars have been considered to be astrophysically interesting objects because they can be candidates of sources of gravitational waves due to the rotational instability (Chandrasekhar 1970, Friedman and Schutz 1978ab).However, if stars rotate differentially, the eigenfunctions of oscillations as well as the equilibrium configurations will be deformed considerably (see e.g. Eriguchi and Müller (1985) for Newtonian equilibrium stars; Komatsu, Eriguchi, and Hachisu (1988) for relativistic equilibrium stars; Karino, Yoshida, and Eriguchi (2001) for eigenfunctions of oscillations of Newtonian stars.) If one would like to apply the Cowling approximation to such differentially rotating stars, it is necessary to estimate the applicability of the Cowling approximation to them. In this paper we will study the accuracy of the Cowling approximation in rotating Newtonian stars by comparing the results obtained by the Newtonian Cowling approximation with those by solving the exact equations for the perturbed equations.

The plan of this paper is as follows. In the next section, we will briefly summarize the scheme to solve the linearized perturbed equations for differentially rotating stars with and without the gravitational perturbation. Our numerical results will be presented in Section 3. In Section 4, the accuracy of the Cowling approximation will be discussed. In the last section, we will summarize our present study and make a brief comment on the applicability of the Cowling approximation to GR. Throughout this paper, the ordinary spherical polar coordinate system $`(r,\theta ,\phi )`$ is used.

## 2 Method

The motivation of this paper is to estimate the accuracy of the Cowling approximation. For this purpose, we will solve the linearized equations of perturbed fluid in Newtonian gravity and numerically obtain the eigenvalues and eigenfunctions of oscillation modes. Then we will compare them with the results computed by neglecting the gravity perturbation, $`\delta \varphi `$, from the basic equations (the Cowling approximation). The method of solving the full system of the basic equations for the oscillations of rotating stars is shown in Karino, Yoshida, Yoshida and Eriguchi (2000). Hence, in this section, we will briefly explain the solving method.

The equilibrium states of axisymmetric polytropes can be obtained by the method proposed by Eriguchi and Müller (1985). The polytropic relation is expressed as

$$p=K\rho ^{1+\frac{1}{N}},$$
(1)
where $`p`$, $`\rho `$, $`K`$ and $`N`$ are the pressure, the density, the polytropic constant and the polytropic index, respectively. In this paper, we compute polytropes with $`N=`$ 0.5, 1.0, and 1.5. Such relatively stiff polytrope can be used as the models of neutron stars.

We will use the following rotation law (j-constant rotation law):

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_0A^2}{(r\mathrm{sin}\theta )^2+A^2},$$
(2)
where $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }_0`$ are the angular velocity and the angular velocity on the rotation axis. The parameter $`A`$ represents a degree of differential rotation. A small value of $`A`$ corresponds to a strongly differential rotation, while in the limit of $`A\mathrm{}`$ the rotation law becomes that of rigid rotation.

The eigenfrequencies and eigenfunctions of the modes can be obtained by solving the perturbed fluid equations. Here, we assume small perturbations and linearize the system of equations. The linearized quantities are expanded as follows:

$$\delta f(r,\theta ,\phi ,t)=\underset{m}{}\mathrm{exp}(i(\sigma tm\phi ))f_m(r,\theta ),$$
(3)
where $`f`$ and $`\delta f`$ are a certain physical quantity and its Eulerian perturbation, respectively. In this expansion, $`\sigma `$ denotes the eigenvalue of the oscillation mode and $`m`$ is the azimuthal wave number. Here, we only consider the $`m=2`$ f- and r-mode for two reasons: (1) from the astrophysical point of views, the $`m=2`$ mode is the most interesting because it has the fastest growth rate in rapidly rotating stars, and (2) it has been said that the accuracy of the Cowling approximation becomes the worst for the $`m=2`$ mode. Actually, the eigenfrequencies of radial and quasi-radial-modes have been studied by nonlinear computations and it has been shown that the frequency obtained by the Cowling approximation is considerably different from the one given by full GR calculations even in the case with rapid rotations (Font, Goodale, Iyer, Miller, Rezzolla, Seidel, Stergioulas, Suen and Tobias 2002).

On the other hand, in the Cowling approximation, the Eulerian perturbation of the gravitational potential is dropped from the full system and solved the equations with the remaining terms.

## 3 Numerical Results

The eigenvalues of the $`m=2`$ f-modes for the $`N=1.0`$ polytrope are plotted in Fig. 1 as a function of $`T/|W|`$, where $`T`$ is the rotational kinetic energy and $`W`$ is the gravitational potential energy of the star. Three curves of the different types correspond to the different values of the parameter $`A`$; $`A=0.5,1.0`$, and $`\mathrm{}`$ (= rigid rotation). For each $`A`$, two curves are plotted: the upper curves correspond to the results of the full computations and the lower to the approximated solutions. In the case of $`A=1.0`$ and 1.5 including the perturbation of gravity, when the rotation of the star get rapid, the so-called neutral point where $`\sigma =0`$ appears. This neutral point is corresponding to the critical point of secular instability due to gravitational radiation reaction (Chandrasekhar 1970, Friedman & Schutz 1978ab). On the other hand, the curves given by the Cowling approximation have smaller, finite values. This means that we can not estimate the critical points of secular instability when we apply the Cowling approximation.

The differences of the eigenfrequencies between the results of the full perturbed systems and the Cowling approximation are shown in Fig. 2. From this figure, it is clear that the accuracy of the Cowling approximation becomes better as the stellar rotation gets faster for rigidly rotating polytropes. For differentially rotating polytropes, however, the accuracy of the approximation becomes worse and worse as the stellar rotation becomes faster.

The comparison between the approximated results and full perturbed results have been done for rigidly rotating polytrope by Jones, Andersson and Stergioulas (2002). They have carried out time evolving computations of the linearized perturbed equations of rapidly and rigidly rotating Newtonian polytrope with the Cowling approximation, and compared their results about f- and r-modes with the previous results. They have checked the accuracy of the Cowling approximation and confirmed certain amount of errors of eigenfrequencies between them. Their results for rapidly rotating stars (Figure 4 in Jones et al. 2002) is well-matched with our result for rigid rotation case (Figure 1), although they have computed upper branch of f-mode and our results are for lower branch. This similarity shows our result is consistent.

There are two reasons why the accuracy of the Cowling approximation becomes worse for differential rotation. The first reason is related to the density distribution of the equilibrium stars. In Figs. 3 and 4, equi-density contours of the rigidly and differentially rotating stars with $`N=1.0`$ and $`r_\mathrm{p}/r_\mathrm{e}=0.6`$ are shown, respectively. Here $`r_\mathrm{p}`$ and $`r_\mathrm{e}`$ are the polar radius and the equatorial radius, respectively. By comparing Fig. 3 and Fig. 4, in the rapidly and differentially rotating star, the dense region near the center extends to the outer region. Remember that the amplitude of the f-mode oscillation is significant in the outer region of the star. Then, in the differentially rotating star, the gravitational potential perturbation comes to play a relatively important role, because the region with high-amplitude is dense. Therefore, the Cowling approximation which neglects the perturbed gravitational potential from the basic equations gives less accurate solutions.

The second factor is the behavior of the eigenfunction due to differential rotation. In Figs. 5 and 6, the density perturbation $`\delta \rho `$ is shown as a function of the relative distance from the rotation axis for rigidly (Fig. 5) and differentially (Fig. 6) rotating $`N=1.0`$ polytropes. In the upper panels of these figures denote the eigenfunctions for the almost spherical configurations whose axis ratio is $`r_\mathrm{p}/r_\mathrm{e}=0.99`$, while the lower panels show those for the rapidly rotating stars with $`r_\mathrm{p}/r_\mathrm{e}=0.60`$. The eigenfunctions in those figures are normalized by the maximum amplitude of the density perturbation in the star. From Fig. 5, the eigenfunctions of rigidly rotating polytropes are the monotonically increasing functions of the distance from the rotation axis and their amplitudes become the maximum at the surface. This behaviors are common with the previous results given in the cases of spherical stars and/or rigidly rotating stars. In Fig. 6, however, the behavior of the eigenfunctions drastically change. When the stellar rotation is slow, the qualitative behaviors of the eigenfunctions that increase monotonically with the stellar radii are the same as those of rigidly rotating polytropes. When the stellar rotation is sufficiently rapid, the maximum values of the eigenfunctions appear in the inner region, not at the surface. This result means that in rapidly and differentially rotating stars the amplitudes of the modes in the high density region are rather large and that the effects of the gravitational perturbation become important.

In order to see the dependence of the accuracy of the Cowling approximation on the equations of state in differentially rotating stars, the eigenvalues of the $`m=2`$ f-modes in differentially rotating ($`A=1.0`$) stars, with several polytropic indices are also computed. The eigenvalues as functions of $`T/|W|`$ are shown in Fig. 7. Three curves of different types correspond to different models with $`N=0.5`$, $`N=1.0`$, and $`N=1.5`$. Similarly as shown in Fig. 1, the upper sequences of the pair curves show the eigenvalues for the full perturbed systems, while the lower sequences correspond to the results of the Cowling approximation.

In Fig. 8, the differences of the eigenfrequencies between the solutions of the full systems and those of the Cowling approximation are shown as in Fig. 2. For rigidly rotating stars, the high compressibility makes the accuracy of the approximation better, because the density distributions in highly compressible star is considerably condensed to the central region, so the density will be significantly low in the outer region where the amplitudes of the oscillations become large. In Fig. 7 and 8, almost the same tendency with the previous one can be observed. For $`N=1.5`$ polytrope, the matter of the star is more centrally condensed compared with the smaller values of $`N`$. Hence, even if the star is rotating differentially with certain degree, the accuracy of the Cowling approximation does not become inaccurate. On the other hand, for $`N=0.5`$ model, since the density of the star does not drop significantly in the outer region, the accuracy of the Cowling approximation becomes considerably inaccurate mainly due to the effect of differential rotation.

## 4 Discussion

### 4.1 f-mode

In the previous section, we have compared the results given by the system of the full perturbed equations with that in the Cowling approximations, and we have found that the accuracy of the approximation becomes worse when we consider rapidly and differentially rotating stars. And the degree of declination of the accuracy depends on the strength of differential rotation. Two reasons for this tendency can be considered. The first reason will be the prolongational effect of the equilibrium density distribution along equatorial plane, by the presence of differential rotation. The second one will be the change that the maximum points of the eigenfunctions appear in the interior regions of the star, not surface regions.

In this section, we will quantitatively estimate these effects of differential rotation. For this purpose, we will compute the gravity perturbations numerically in Newtonian gravity. As mentioned above, since the reliability of the Cowling approximation becomes better when the density of the fluid components is sufficiently low in the region where the amplitudes of the oscillation get relatively larger, in order to estimate the relation between the density in an equilibrium and the amplitudes of the potential perturbations , we define a new quantity as follows:

$$W_\varphi \rho \delta \overline{\varphi }_m\mathrm{exp}(i(\sigma tm\phi ))𝑑V,$$
(4)
where, $`\delta \overline{\varphi }=\delta \varphi /\delta \varphi _{\mathrm{max}}`$ is the perturbed gravitational potential normalized by the maximum amplitude of $`\delta \varphi `$ in the star. 111 Additionally, here, since it is assumed that the equilibrium configurations are axisymmetric and the perturbations have the dependency as $`\mathrm{exp}(i(\sigma tm\phi ))`$ in the azimuthal direction, the integral in the $`\phi `$-direction, $`\mathrm{exp}(i(\sigma tm\phi ))𝑑\phi `$, need not to be carried out in the present context. In this meaning, we express this quantity as $`\overline{W}_\varphi `$.  These quantities, $`W_\varphi `$, for several sequences of $`N=1.0`$ polytropes are shown in Fig. 9, for $`A=1.0`$, and 0.5. The horizontal axis denotes $`T/|W|`$. Here, actually plotted values, $`\overline{W}_\varphi `$, are normalized by the masses of the equilibrium models:

$$M_0=\rho 𝑑V.$$
(5)
Introducing such a quantity and comparing Figs. 9 with Fig. 2, we can see that the behavior of $`W_\varphi `$ becomes consistent with the absolute values of $`\delta \sigma `$ in terms of the function of $`A`$ and $`T/|W|`$. These tendencies of $`W_\varphi `$ strongly suggest that our conjecture described in previous sections are correct; i.e. the accuracy of the Cowling approximation is related with the values of $`W_\varphi `$.

### 4.2 r-mode

In these days, from astrophysical point of view, r-mode oscillations in neutron stars are noticed. A lot of studies which suggest that r-modes can be unstable in every rotating star and the accompanied gravitational waves can be detected by improved interferometer as like LIGO II has been done (Andersson 1998, Friedman and Morsink 1998, for recent works, see the review of Andersson 2002 and its references).

Since the eigenfunctions of the r-mode oscillations are surface-standing, they may not suffer from the effect of extension of the central dense regions. Additionally the gravity perturbation is much smaller than the fluid perturbations for r-modes (about the feature of r-mode oscillations, Saio 1982, Yoshida and Lee 2002, for the eigenvalue of r-modes in differentially rotating stars, see Karino et al. 2001). Because of these reasons, the Cowling approximation may be accurate in the calculations of r-mode oscillations. In order to prove this, we have carried out the same computation with f-modes for the $`m=2`$ r-mode. The results are shown in Figure 10. The plot is in the same manner with Figure 1. From this result, it can be recognized that the approximation will give good results for r-mode oscillations. This conclusion is also comparable with the assertion of Jones et al. (2002).

## 5 Concluding Remarks

In the present study, the accuracy of the Cowling approximation has been investigated for rapidly and differentially rotating stars in Newtonian gravity. We have adopted a polytropic equation of state with polytropic indices $`N=0.5`$, 1.0, and 1.5 as the stellar models and estimated the accuracy of the approximation for the $`m=2`$ f- and r-mode oscillations.

For the f-mode in rigidly rotating stars, the tendency that the accuracy gets better for rapidly rotating stars has been reproduced. On the other hand, when the angular momentum distribution of the star is not uniform, the accuracy of the Cowling approximation becomes worse when the rotation rate gets higher. The accuracy of the approximation also depends on the equation of state. For softer equations of state, the accuracy of the approximation becomes relatively better even if the effect of differential rotation exists, because of the central condensation of the matter. Contrary, for stiffer equations of state, which are usually used to describe the neutron stellar matter, the accuracy of the approximation becomes worse for rapidly and differentially rotating stars. The reasons why the accuracy becomes worse can be considered as follows. First, when the angular momentum distribution is not uniform, the equilibrium configurations are highly deformed and the high density region extends to the outer part of the equatorial region. Second, when the angular momentum distribution is not uniform, the maximum amplitude of the eigenfunction appears well inside the star where the equilibrium density is rather high. It is possible to quantitatively confirm that the oscillations become large in the high density region, by using a new quantity, $`W_\varphi `$, which is the integration of the density weighted by the eigenfunctions of the oscillations.

As discussed in Introduction, it is extremely difficult to solve oscillations of general relativistic rotating compact stars. This is mainly because it is difficult to treat the perturbations of the metric, i.e. the perturbed gravitational potentials. Thus it is desirable to apply the Cowling approximation to general relativistic configurations and obtain reasonably accurate approximated solutions.

Can we say anything about the applicability of the Cowling approximation to general relativistic configurations ? It is widely believed that the accuracy of the Cowling approximation is better for general relativistic configurations, because the strong gravity makes the effective central condensation of the stellar matter higher. This effect of gravitation is similar to the effect of the increasing of the compressibility. As shown in the present paper, the ascent of the compressibility achieves good accuracy of the approximation. Therefore, in GR, it is expected that the accuracy of the Cowling approximation would become better than that for Newtonian gravity.

However, also in the case of GR, the same problem accompanied with differential rotation will happen. When we apply the Cowling approximation in differentially rotating relativistic stars, the accuracy of the approximation will depend on the positive effect of strong gravity and the negative effect of differential rotation. The detailed study must be carried out to discuss about the accuracy of such an approximation in differentially rotating compact stars. As mentioned in the last section, however, the Cowling approximation works very well in the calculations of r-mode oscillations (also see, Yoshida & Lee 2002). Hence, it is good idea to apply the approximation to obtain the eigenfrequencies of r-modes in general relativity. Also it is shown that the accuracy of the Cowling approximation becomes better for higher modes (Jackson 1995), therefore to compute those higher modes by using this approximation is one of the realistic usages of the Cowling approximation.

## Acknowledgments

The present author thanks Yoshiharu Eriguchi for useful discussion, comments, and continuous encouragements. He would like to thank Keisuke Taniguchi and Koji Uryu for their careful reading and useful comments. He thanks Yasufumi Kojima for his comments. This work was in part supported by JSPS Research Fellowship for Young Scientists.


text:
# Fundamental and vortex solitons in a two-dimensional optical lattice

## Abstract

Fundamental and vortex solitons in a two-dimensional optically induced waveguide array are reported. In the strong localization regime, the fundamental soliton is largely confined to one lattice site, while the vortex state comprises of four fundamental modes superimposed in a square configuration with a phase structure that is topologically equivalent to the conventional vortex. However, in the weak localization regime, both the fundamental and vortex solitons spread over many lattice sites. We further show that fundamental and vortex solitons are stable against small perturbations in the strong localization regime. OCIS codes: 190.0190.

Discrete optical spatial solitons are localized modes of nonlinear waveguide arrays that form when discrete diffraction is balanced by nonlinearity. When a low intensity beam is focused into a waveguide array, the propagating field spreads over many sites (due to optical tunneling) exhibiting a typical discrete diffraction pattern with the intensity mainly concentrated in the outer lobes . However, at sufficiently high power, the beam self-traps to form a localized state.

Discrete solitons in waveguide arrays were first predicted to exists as solutions to the discrete nonlinear Schrödinger equation and later observed in AlGaAs waveguide arrays. This experimental observation stimulated much new research such as diffraction management, soliton interaction and beam steering, and discrete solitons in a two-dimensional (2D) photorefractive optical lattice . The experiments in  are particularly interesting since the 2D waveguide array there is formed optically, thus it is very versatile and easily tunable. This allowed a host of 2D localization phenomena to be observed in . However, many questions on discrete solitons in 2D optical lattices still remain open. For instance, stability properties of such solitons have not been carefully studied. In addition, new localization states such as those analogous to vortex solitons in a homogeneous waveguide have not been explored.

In this Letter, fundamental and vortex solitons in a 2D optical lattice are found and their stability analysed. In the strong localization regime, the fundamental soliton is confined largely on one lattice site with a uniform phase. The vortex soliton comprises of four fundamental modes located at the bottoms of the optical potential in a square configuration with a phase structure that is topologically equivalent to the conventional vortex. By winding around the zero intensity position along any simple closed curve, the phase of the vortex state acquires $`2\pi `$ increment. We call this structure a vortex cell. When the localization is weak, both fundamental solitons and vortex-cells spread over many lattice sites. In the strong localization regime, we show that the fundamental soliton as well as the vortex cell are both linearly and nonlinearly stable under weak perturbations.

We start with the $`2D`$ nonlinear Schrödinger equation

$$i\frac{\psi }{z}+\left(\frac{^2}{X^2}+\frac{^2}{Y^2}\right)\psi V\psi +|\psi |^2\psi =0,$$
(1)
where $`V=V_0\left(\mathrm{cos}^2X+\mathrm{cos}^2Y\right)`$ is the optical lattice potential, and $`V_0`$ is its intensity. Such a potential can be obtained by optically interfering two pairs of laser beams . Transverse distances $`X`$ and $`Y`$ in (1) have been normalized by the lattice period $`D`$, and distance $`z`$ normalized by $`4n_0D^2/\pi \lambda `$, where $`n_0`$ is the refractive index, and $`\lambda `$ is the wavelength of the beam . For typical values $`D=9\mu m`$, $`n_0=2.3`$ and $`\lambda =0.5\mu m`$ , the unit distance $`z`$ in (1) corresponds to 0.5mm in physical distance. Without the lattice potential, solitons would suffer collapse under small perturbations . However, as we shall show here, optical lattices can suppress the collapse of fundamental solitons and vortex-cells.

Eq. (1) conserves two quantities: the power $`P=_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}|\psi |^2𝑑X𝑑Y`$, and the energy $`E`$:

$$E=_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\left\{|\psi |^2\frac{1}{2}|\psi |^4+V|\psi |^2\right\}𝑑X𝑑Y.$$
(2)
We look for stationary solutions of the form $`\psi (X,Y,z)=e^{i\mu z}u(X,Y)`$ where $`\mu `$ being the propagation constant of the soliton. Then $`u(X,Y)`$ satisfies

$$\frac{^2u}{X^2}+\frac{^2u}{Y^2}Vu+|u|^2u=\mu u.$$
(3)
Solutions to this equation can be obtained by a Fourier-iteration method as was suggested by Petviashvili .

Fundamental solitons. A fundamental soliton of Eq. (3) has a single main hump sitting at a bottom of the potential, say $`(X,Y)=(\pi /2,\pi /2)`$. Two examples corresponding to propagation constants $`\mu =0`$ and $`0.88`$ with $`V_0=1`$ are displayed in Fig. 1(c,d). We see that for small $`\mu `$ (Fig. 1c), the beam is largely confined on one lattice site, while at higher $`\mu `$ (Fig. 1d), it spreads over many lattice sites. To quantify these solitons, we calculate the dependence of the normalized power $`P`$ on the propagation constant $`\mu `$ for $`V_0=1`$ as displayed in Fig. 1a. When $`\mu \mathrm{}`$, $`P`$ approaches a constant 11.70. This is apparently because in this limit, the fundamental soliton is highly localized, thus it approaches the lattice-free fundamental-soliton state, which has critical power $`P_c11.70.`$ As $`\mu `$ goes to a cut-off value which is approximately 0.95, $`P`$ appears to go to infinity. In this limit, the fundamental state becomes uniformly distributed in space. Thus, this cut-off value should be the boundary of the band gap in the linear-wave spectrum . When $`\mu 0.72`$, $`P`$ is minimal. When $`\mu >0.72`$, $`dP/d\mu >0`$. In this region, the Vakhitov-Kolokolov (VK) theorem suggests that the soliton is linearly unstable due to the presence of a purely-real unstable eigenvalue $`\sigma `$ in the linearized equation . We have confirmed this instability by numerically simulating the linearized version of (1) around the above soliton. The results for unstable eigenvalues are shown in Fig. 1b. On the other hand, when $`\mu <0.72`$ ($`V_0=1`$), the VK instability is absent, and the soliton state is linearly stable.

For the two-dimensional self-focusing case, collapse is an important issue. The above linear stability analysis does not guarantee that the fundamental soliton will not collapse under small perturbations. In the study of collapse, the energy $`E`$ plays an important role. In the absence of the lattice potential or when the potential is harmonic, the soliton collapses if its energy is negative. In case the energy is positive, however, the soliton collapses only if it is strongly perturbed. For the optical lattice, we have calculated the energy $`E`$ of fundamental solitons at various values of $`\mu `$ and plotted the results in Fig. 1(a). The energy is found to be always positive. Thus we can expect that this state is able to withstand small perturbations without collapse. To confirm the above expectations on stability, we numerically study the nonlinear evolution of the fundamental soliton under small perturbations by directly simulating Eq. (1) with initial condition

$$\psi (X,Y,z=0)=u(X,Y)[1+ϵu_p(X,Y)],$$
(4)
where $`ϵ1`$, and $`u_p(X,Y)`$ is the initial perturbation. We first take $`u_p`$ to be white noise. A large number of simulations with small $`ϵ`$ and various realizations of random-noise perturbations have been performed, and we have found that for $`V_0=1`$, if $`\mu <0.72`$, the fundamental soliton is indeed stable against white noise perturbations; when $`\mu >0.72`$, the soliton is unstable. To study the nonlinear evolution process, we now take $`u_p=1`$. For $`V_0=1`$, $`\mu =0`$ and $`ϵ=\pm 0.01`$, nonlinear evolutions are plotted in Fig. 2a. We see that the perturbed soliton only oscillates weakly around the fundamental-soliton state, meaning that the soliton is both linearly and nonlinarly stable. On the other hand, at $`V_0=1`$ and $`\mu =0.88`$ (where the soliton is linearly unstable), the dynamics is different as two scenarios are identified: (i) at higher input power ($`ϵ>0`$), the perturbed state relaxes into a $`z`$-periodic bound state (Fig.2b,c); (ii) at lower input power ($`ϵ<0`$), the perturbed state decays into linear Bloch waves (Fig. 2b,d) (similar scenarios can be found in  for a different system). Collapse is not observed in either scenarios.

Vortex cell. In addition to the fundamental solitons, we have numerically found vortex solitons as well. Two examples with $`V_0`$=1, $`\mu =0`$ and $`0.82`$ are shown in Fig. 3 (b,c,d). At $`\mu =0`$ (strong localization regime), the vortex state comprises of four fundamental solitons superimposed in a square configuration with a phase structure which is topologically equivalent to a conventional vortex \[see Fig. 3b,c\]. By winding around the center along any closed curve, the phase of the vortex acquires $`2\pi `$ increment, thus, we name it a vortex cell. At $`\mu =0.82`$ (weak localization regime), the vortex cell spreads out to more lattice sites and becomes more intricate, as can be seen in Fig. 3d. But its phase structure is almost the same as with $`\mu =0`$. We should point out that these vortex cells are different from conventional vortices without optical lattice on a major aspect: their densities and phases depend on both $`r`$ and $`\theta `$. The normalized power and energy diagrams versus propagation constant for these vortex cells at $`V_0=1`$ are shown in Fig. 3(a). This figure is similar to that for fundamental solitons \[Fig. 1(a)\], except that both $`P`$ and $`E`$ here are about four times larger. Unstable eigenvalues $`\sigma `$ of vortex-cells are determined by simulating the linearized equation (1) around vortex-cells. The results are shown in Fig. 4a for $`V_0=1`$. We see that vortex cells experience an oscillatory instability for $`\mu >7.8`$ and become stable for $`\mu <7.8`$. These vortex cells also suffer the VK instability in the region $`\mu >0.73`$ where $`dP/d\mu >0`$ \[see Fig. 3a\]. However, the oscillatory instability is much stronger as it occurs over a wider region and has a higher growth rate. Fig. 3(a) shows that vortex cells also have positive energy. Thus, if a vortex cell is linearly stable, it should be able to resist collapse under small perturbations . However, if it suffers the linear oscillatory instability discussed above, this instability could result in power exchange from one part of the cell to another so that the intensity at some small spots becomes high, triggering local collapse to occur. We have observed this scenario numerically. An example is shown in Fig. 4(b,c,d) which displays the development of a vortex cell with $`\mu =0`$ and $`V_0=1`$ when it is amplified by 1% initially.

We next discuss the effect of varying the potential strength $`V_0`$ on the formation and stability of fundamental and vortex solitons. For this purpose, we have chosen $`V_0=1.5`$, and repeated most of the above calculations. The results are summarized as follows: (i) fundamental solitons and vortex-cells exist at higher values of $`\mu `$ (up to approximately 1.35); (ii) at the same value of $`\mu `$, both the fundamental and vortex solitons have lower power ($`P`$) and higher values of energy ($`E`$); (iii) the oscillatory instability suffered by vortex cells is reduced (it could even be completely suppressed when $`V_0`$ becomes even larger). Thus, we conclude that higher lattice potentials stabilize both fundamental and vortex solitons.

In conclusion, we have studied new types of fundamental and vortex solitons in a two-dimensional optical lattice potential, and shown that both solitons are stable in the strong localization regime.

The work of J.Y. was supported by NSF and NASA.


text:
# References

Superconducting Phase Diagram of Li Metal in NearlyHydrostatic Pressures up to 67 GPa

Shanti Deemyad and James S. Schilling

Department of Physics, Washington UniversityC. B. 1105, One Brookings Dr., St. Louis, MO 63130

Abstract. The dependence of the superconducting transition temperature $`T_c`$on nearly hydrostatic pressure has been determined to 67 GPa in an ac susceptibility measurement for a Li sample embedded in helium pressure medium. With increasing pressure, superconductivity appears at 5.47 K for 20.3 GPa, $`T_c`$rising rapidly to $``$14 K at 30 GPa. The $`T_c\left(P\right)`$-dependence to 67 GPa differs significantly from that observed in previous studies where no pressure medium was used. Evidence is given that superconductivity in Li competes with symmetry breaking structural phase transitions which occur near 20, 30, and 62 GPa. In the pressure range 20 - 30 GPa, $`T_c`$is found to decrease rapidly in a dc magnetic field, the first evidence that Li is a type I superconductor.

Of the 52 elemental solids known to be superconducting, 23 enter this state only if compressed under sufficiently high pressures . The latest confirmed member of this ever growing family is the alkali metal Li, an element devoid of superconductivity at ambient pressure to temperatures as low as 4 mK . Following an early report of possible pressure-induced superconductivity in Li near 7 K , recent electrical resistivity studies by Shimizu et al  to 48 GPa (480 kbar) followed by ac susceptibility and resistivity studies by Struzhkin et al  to 80 GPa confirmed the onset of a superconducting state for pressures above 20 GPa. The observed dependences of $`T_c`$ on pressure in these three experiments, however, are in poor agreement. This may result from the fact that no pressure medium was used, the ultrahard diamond anvils  or boron-nitride spacers  pressing directly onto the Li sample, subjecting it to shear stress and plastic flow. Shear-stress effects on $`T_c(P)`$ are well known from studies on such diverse superconducting materials as organic metals , high-$`T_c`$ oxides , MgB2 , and Re metal ; $`\beta `$Hg, a shear-stress-induced body-centered tetragonal modification of ordinary rhombohedral $`\alpha `$Hg, exhibits distinctly different superconducting properties . In a substance like Li, where a multitude of potential phases lie very close in energy , the shear stress may be sufficient to induce structural phase transitions.

Shimizu et al  reported that $`T_c`$ decreases only slightly in moderate magnetic fields, a field of 30,000 Oe being required to completely suppress Li’s superconductivity at 34 GPa. This would indicate that Li is a type II superconductor, an unexpected result for an elemental superconductor, particularly for a simple (s,p-electron) metal, like Li.

Neaton and Ashcroft  have obtained the counterintuitive result that under sufficient compression the electronic structure of fcc Li departs radically from free-electron-like behavior. The anomalous increase in the magnitude of the pseudopotential (and the electron-ion interaction) in the fcc phase under pressure leads not only to the possibility of a superconducting state in Li, one where $`T_c`$ might increase significantly with pressure, but also to possible structural transitions to phases with reduced symmetry, such as a Li-ion pairing phase  or the reduced-symmetry $`hR`$1 or $`cI`$16 phases observed by Hanfland et al  for Li near 180 K. In Li, therefore, superconductivity and symmetry breaking phase transitions are expected to compete with each other  and possibly lead to an anomalous variation in $`T_c`$ with pressure. Other superconducting properties, such as whether type I or type II, values of the critical field(s), gap size, etc., may also be anomalous and should be determined.

Ab initio electronic structure calculations by Christensen and Novikov  predict that under compression fcc Li should exhibit superconductivity where $`T_c`$ increases rapidly with pressure, reaching values as high as 50-70 K. An increase in $`T_c`$ under pressure is highly anomalous for a simple (s,p-electron) metal, like Li. For all known ambient-pressure simple-metal superconductors, $`T_c`$ is found to decrease under pressure because lattice stiffening effects normally dominate over electronic effects . The accurate determination for Li of the intrinsic dependence of $`T_c`$ on pressure is, therefore, of considerable importance.

In this paper we present an extensive determination of the superconducting phase diagram, $`T_c(P),`$ to 67 GPa for a Li sample surrounded by helium, the most hydrostatic pressure medium known. $`T_c`$ is derived from the magnetic signature in the ac susceptibility at 1000 Hz and 3 Oe rms. Over the pressure range to 67 GPa, the pressure dependence $`T_c(P)`$ differs significantly from the results of earlier studies where no pressure medium was used . After the onset of superconductivity at 5.47 K for 20.3 GPa, $`T_c`$ rises rapidly to $``$ 14 K at 30 GPa, followed by an abrupt change in the sign of $`dT_c/dP`$, signalling a likely phase transition from fcc to a lower symmetry structure. This transition at 30 GPa appears to illustrate the competition between the superconducting state and symmetry breaking phase transitions suggested above. In addition, $`T_c`$ is found to decrease rapidly with applied dc magnetic field, giving the first evidence that Li is a type I superconductor.

Because Li reacts readily with oxygen, water vapor, and nitrogen, all steps of the sample preparation and pressure cell loading were carried out in an Ar-gas glove box with continuous purification. High pressures were generated using a diamond-anvil cell (DAC) made from nonmagnetic CuBe. Further experimental details of the DAC and ac susceptibility measurement have been reported previously .

In the present experiment the miniature Li samples were cut from foil (Alfa Aesar 99.9%) and placed in a 250 $`\mu `$m dia. hole drilled through the center of a preindented rhenium gasket (see Fig. 1). Tiny ruby spheres allow the determination  of the pressure in situ to $`\pm `$ 0.2 GPa at 20 K. Before sealing the gasket hole shut with the opposing diamond anvils, the hole is flooded with liquid helium which acts as pressure medium. The pressure was changed at temperatures in the range 150 - 180 K; the Li sample was kept at temperatures below 180 K during the entire high-pressure experiment. After the conclusion of the experiment, the Li sample and rhenium gasket were carefully examined. The preindented gasket thickness was found to be reduced from the original 70 $`\mu `$m to 35 $`\mu `$m which is still greater than that of the Li sample (20 - 25 $`\mu `$m). That neither the diamond anvils nor the gasket wall pressed directly onto the Li sample can be seen in Fig. 1 (lower right) and can also be inferred from the irregular surface and unaltered shape of the sample after the experiment.

In the present experiment, no superconducting transition was observed for Li in the measured temperature region 4 - 60 K for applied pressures of 2, 11.5, 19.2, and 19.7 GPa. However, as seen in Fig. 2, at a pressure of 20.3 GPa a clear superconducting transition appears at 5.47 K which increases rapidly with pressure to $``$ 14 K at 30.2 GPa before falling and rising again above 50 GPa. The magnitude of the transition is consistent with 100% shielding. In Fig. 3 we summarize the results of our two extensive experiments determining the superconducting phase diagram of Li metal. In the first experiment (run 1) the pressure was first increased in 17 steps to 67 GPa (solid circles) and then decreased in 13 steps to ambient pressure (open circles), both in monotonic fashion. In run 2 (solid triangles) the pressure was increased in 8 steps monotonically to 36.5 GPa; a blockage in the cooling line prevented further experimentation. The reproducibility of the data for runs 1 and 2 is seen to be excellent; the $`T_c(P)`$-dependence for increasing or decreasing pressure shows only a small pressure hysteresis.

We now discuss in more detail the salient features exhibited by the superconducting phase diagram in Fig. 3. The sudden appearance of superconductivity in Li at 5.47 K for 20.3 GPa is likely due to a transition from the low-temperature Rh6 phase to the fcc structure, as suggested previously . Increasing the pressure above 20 GPa, $`T_c`$ is seen to rise rapidly at the rate $`dT_c/dP+0.9`$ K/GPa, comparable to previous findings , to a value near 14 K at 30 GPa. Above this critical pressure the derivative abruptly changes sign to $`dT_c/dP`$ -0.4 K/GPa, giving clear evidence for a second pressure-induced structural phase transition at 30 GPa. This conclusion is reinforced by the reversible broadening of the superconducting transition in the region 30 - 36 GPa. In view of the diffraction results on Li at 180 K by Hanfland et al , we tentatively identify this second transition as fcc to cI16, although a transition from fcc to another phase, such as the intermediate hR1 phase, is certainly possible. The data near 50 GPa is not sufficiently dense to allow one to speculate whether the minimum in $`T_c(P)`$ near this pressure is characteristic for a single phase or arises from a structural phase transition. The abrupt disappearance of superconductivity above 62 GPa, however, signals a further phase transition, perhaps from cI16 to some unknown phase. Diffraction experiments on Li at temperatures below 50 K under nearly hydrostatic pressure conditions are clearly needed to arrive at an unequivocal structure assignment in the low-temperature region.

Whatever the exact nature of the structural phase transitions indicated by the data in Fig. 3 for pressures above 20 GPa, it is almost certain that they are from the high symmetry fcc structure to phases of lower symmetry, such as cI16 or hR1. The highly anomalous increase in $`T_c`$ with pressure between 20 and 30 GPa in the fcc phase agrees with the trend indicated in the calculation of Christensen and Novikov  and bears witness to the large enhancement in the pseudopotential and the electron-ion interaction with pressure predicted earlier by Neaton and Ashcroft . We speculate that with increasing pressure the electron-ion interaction in Li is enhanced to the extent that at 30 GPa a structural transition from fcc to a lower symmetry structure finally occurs. This transition removes, or at least relieves, the anomalous enhancement of the pseudopotential ocurring in the fcc phase, allowing $`T_c`$ to decrease with pressure due to dominant lattice stiffening effects, as in all canonical simple-metal superconductors . As the pressure is increased further, however, the story repeats itself and $`T_c`$ passes through a minimum at 50 GPa and increases again due to a renewed enhancement of the pseudopotential, only to fall again above 62 GPa following a further symmetry lowering phase transition. To test this scenario, detailed electronic structure calculations for the phases actually present at low temperatures would be very useful.

The anomalous behavior of $`T_c(P)`$ for the “simple” s,p-metal Li arises from the close proximity of the electron cores of neighboring ions under high compression. This situation will not be relieved, and Li allowed to revert to its former free-electron-like behavior, before such astronomically high pressures (far higher than those in the present experiment!) are applied as to break up the atomic shell structure itself, resulting ultimately in a Thomas-Fermi electron gas.

To further characterize the superconducting state of Li, in run 2 we subjected the sample to dc magnetic fields up to 260 Oe. In the inset to Fig. 4,$`T_c`$ at 21.9 GPa pressure is seen to decrease from 5.95 K to 5.21 K when a field of only 200 Oe is applied, yielding the initial critical field slope $`(dH_c/dT)_{T_c}`$ -270 Oe/K. Since the applied fields are insufficient to directly determine the critical field at 0 K, $`H_c(0)`$, from the present data, we use the expression $`H_c(0)=\frac{1}{2}T_c(dH_c/dT)_{T_c},`$ derived from the standard empirical relation $`H_c(T)=H_c(0)[1(T/T_c)^2]`$ , to obtain the estimate from the 21.9 GPa data that $`H_c(0)`$ -$`\frac{1}{2}`$(5.95 K)(-270 Oe/K) $`=`$ 800 Oe. Critical field data $`H_c(T)`$ at eight further pressures were obtained, five of which are shown in Fig. 4. In the pressure range 20 to 24 GPa the values of $`(dH_c/dT)_{T_c}`$, $`H_c(0)`$ and $`T_c`$ for Li are comparable to those found for canonical type I simple-metal superconductors like Hg and Pb where $`(dH_c/dT)_{T_c}`$ -198.2 Oe/K and -237.3 Oe/K, $`H_c(0)`$ 412 Oe and 803 Oe, and $`T_c`$ 4.15 K and 7.195 K, respectively . For $`P<`$ 30 GPa, the evidence that Li is a type I superconductor is thus quite compelling, although we cannot exclude the possibility of very weak type II behavior. All known ambient-pressure simple-metal (s,p electron) superconductors are type I if in a pure and strain-free condition.

For pressures slightly above versus slightly below 30 GPa, it is interesting to note that the initial slope $`(dH_c/dT)_{T_c}`$ increases approximately twofold, likely signalling a comparable increase in $`H_c(0)`$. An abrupt change in $`(dH_c/dT)_{T_c}`$ and/or $`H_c(0)`$ would be consistent with the occurrence of a structural phase transition at 30 GPa. Since for type I superconductivity we have the relation $`H_c(0)\frac{1}{2}[\mathrm{\Delta }(0)]^2N(E_f)`$ , such an increase in $`H_c(0)`$ could arise from an increase in the magnitude of either the superconducting gap $`\mathrm{\Delta }(0)`$ or the density of states $`N(E_f)`$. Alternatively, the phase transition at 30 GPa may create lattice defects which lead to weak type II behavior and an increase in the (upper) critical field.

Shimizu et al  reported from resistivity measurements at 34 GPa, where $`T_c`$ 7 K, that $`H_c(0)`$ 30,000 Oe, a value 6$`\times `$ larger than the maximum value of $`H_c(0)`$ estimated in the present experiment. This high value of $`H_c(0)`$ clearly points to type II superconductivity for their Li sample. Since in their experiment the direct contact of the Li sample with the diamond anvils and stiff gasket generated relatively large shear stresses, it is conceivable that the resulting plastic deformation led to a high density of defects with a sharp reduction in the electronic mean-free-path, thus promoting type II behavior.

Acknowledgments. The authors thank J. Hamlin, A.-K. Klehe, V.V. Struzhkin, V. Tissen, and especially S. Solin for critical reading the original manuscript. Experimental support by J. Hamlin is acknowledged. One of the authors (JSS) is grateful to N.W. Ashcroft and J.B. Neaton for stimulating discussions. Research is supported by NSF grant DMR-0101809.

Figure Captions

Fig. 1. (upper) Reflected-light photograph into Ar-gas glove box of Au-plated rhenium gasket preindented with diamond anvil (0.5 mm dia. culet). In 250 $`\mu `$m dia. hole are seen for run 2 at ambient pressure the Li sample at 6 o’clock plus clusters of ruby spheres. (lower) Transmitted-light photograph of hole containing Li sample at (left) ambient pressure and (right) 30 GPa. All three photographs are to same relative scale.

Fig. 2. Superconducting transition of Li in the real part of the ac susceptibility (1000 Hz, 3 Oe rms) for 8 values of the pressure. All data are to same scale but are shifted vertically for clarity. All data are from run 1 for increasing pressure, except data at 20.3 GPa (run 2) and 43.0 GPa (run 1, decreasing pressure).

Fig. 3. Superconducting phase diagram of Li metal under nearly hydrostatic pressure: run 1 increasing pressure ($``$), run 1 decreasing pressure ($``$), run 2 increasing pressure ($`\mathrm{}`$). The error in pressure is $`\pm `$ 0.2 GPa. $`T_c`$ is determined by superconducting midpoint to $`\pm `$ 50 mK; vertical “error bars” give temperatures of the superconducting onset and completion. No superconducting transition is observed above 4 K for pressures below 20 GPa or at 67 GPa. Dashed lines are guides to eye with slopes $`dT_c/dP`$ +0.9 K/GPa (left) and -0.4 K/GPa (right).

Fig. 4. Data points ($`,\mathrm{}`$) give to 220 Oe the critical magnetic field $`H_c`$ versus temperature at various pressures for superconducting Li. Solid, dashed, and dotted lines are obtained using the empirical expression $`H_c(T)=H_c(0)[1(T/T_c)^2]`$ to allow fits to the present data for Li in comparison to ambient-pressure data for Pb. For $`P=`$ 28.6, 31.8, and 36.5 GPa we find $`H_c(0)`$ 2400, 5000, and 4500 Oe, respectively. Inset shows ac susceptibility data at 21.9 GPa for applied magnetic fields of 0 and 200 Oe.