format: text

text:
# CMB 𝐵-polarization to map the Large-scale Structures of the Universe

## I Introduction

In the new era of precision cosmology we are entering in, the forthcoming experiments will provides us with accurate data on Cosmic Microwave Background anisotropies. This should lead to accurate determinations of the cosmological parameters, provided the large-scale structures of the Universe indeed formed from gravitational instabilities of initial adiabatic scalar perturbations. It has been soon realized however that even with the most precise experiments, the cosmological parameter space is degenerate when the primary CMB anisotropies alone are considered. Complementary data, that may be subject to more uncontrollable systematics are thus required, such as supernovae surveys (but see ) or constraints derived from the large-scale structure properties. Among the latter, weak lensing surveys are probably the safer, but still have not yet proved to be accurate enough with the present day observations.

Secondary CMB anisotropies (i.e. induced by a subsequent interaction of the photons with the mass or matter fluctuations) offer opportunities for raising this degeneracy. Lens effects are particularly attractive since they are expected to be one of the most important.They also are entirely driven by the properties of the dark matter fluctuations, the physics of which involve only gravitational dynamics, and are therefore totally controlled by the cosmological parameters and not by details on galaxy or star formation rates. More importantly an unambiguous detection of the lens effects on CMB maps would be a precious confirmation of the gravitational instability picture. Methods to detect the lens effects on CMB maps have been proposed recently. High order correlation functions, peak ellipticities or large scale lens induced correlators have been proposed for detecting such effects. All of them are however very sensitive to cosmic variance since lens effect is only a sub-dominant alteration of the CMB temperature patterns. The situation is different when one considers the polarization properties. The reason is that in standard cosmological models temperature fluctuations at small scale are dominated by scalar perturbations. Therefore the pseudo-scalar part, the so called $`B`$ component, of the polarization is negligible compared to its scalar part (the $`E`$ component) and can only be significant when CMB lens couplings are present. This mechanism has been recognized in earlier papers. The aim of this paper is to study systematically the properties of the lens induced B field and uncover its properties.

In section II, we perturbatively compute the lens effect on the CMB polarization $`E`$ and $`B`$ field. This first order equation is illustrated by numerical experiments. Possibility of direct reconstruction of the projected mass distribution is also examined. As it has already been noted a significant fraction of the potential wells that deflect the CMB photons can actually be mapped in local weak lensing surveys. This feature has been considered so far in relation to the CMB temperature fluctuations. We extend in Section III these studies to the CMB polarization exploiting the specificities of the field found in previous section. In particular we propose two quantities that can be built from weak lensing and Cosmic Microwave Background polarization surveys, the average value of which does not vanish in presence of CMB lens effects. Compared to direct analysis of the CMB polarization, such tools have the joint advantage of being less sensitive to systematics –systematic errors coming from CMB mapping on one side and weak lensing measurement on the other side have no reason to correlate!– and so emerge even in presence of noisy data, and of being an efficient probe of the cosmological constant. Indeed the expected amplitude of correlation is directly sensitive to the relative length of the optical bench, from the galaxy source plane to the CMB plane, which is mainly sensitive to the cosmological constant. Filtering effects and cosmic variance estimation of such quantities are considered in this section as well.

## II Lens effects on CMB polarization

### II.1 First order effect

Photons emerging from the last scattering surface are deflected by the large scale structures of the Universe that are present on the line-of-sights. Therefore photons observed from apparent direction $`\stackrel{}{\alpha }`$ must have left the last scattering surface from a slightly different direction, $`\stackrel{}{\alpha }+\stackrel{}{\xi }(\stackrel{}{\alpha })`$, where $`\stackrel{}{\xi }`$ is the lens induced apparent displacement at that distance. The displacement field is related to the angular gradient of the projected gravitational potential. In the following, the lens effect will be described by the deformation effects it induces, encoded in the amplification matrix,

$`𝒜^1`$ $`=`$ $`\left(\begin{array}{cc}1\kappa \gamma _1& \gamma _2\\ \gamma _2& 1\kappa +\gamma _1\end{array}\right)`$ (3)
$`=`$ $`\delta _i^j+\xi _{,i}^j`$ (4)
so that

$`\kappa `$ $`=`$ $`{\displaystyle \frac{1}{2}}(\xi _{,x}^x+\xi _{,y}^y)`$
$`\gamma _1`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\xi _{,x}^x\xi _{,y}^y)`$
$`\gamma _2`$ $`=`$ $`\xi _{,x}^y=\xi _{,y}^x.`$ (5)
The lens effect affects the local polarization just by moving the apparent direction of the line of sight. Thus, if we use the Stokes parameters $`Q`$ and $`U`$ to describe the local polarization vector,

$$\stackrel{}{P}=\left(\begin{array}{c}Q\\ U\end{array}\right)$$

we can relate the observed polarization $`\widehat{\stackrel{}{P}}`$ to the primordial one by the relation

$$\widehat{Q}(\stackrel{}{\alpha })=Q(\stackrel{}{\alpha }+\stackrel{}{\xi }),\widehat{U}(\stackrel{}{\alpha })=U(\stackrel{}{\alpha }+\stackrel{}{\xi }).$$
(6)
From now on we will denote $`\widehat{x}`$ an observed quantity and $`x`$ the primordial one. $`\stackrel{}{\alpha }^{}=\stackrel{}{\alpha }+\stackrel{}{\xi }`$ is the sky coordinate system for the observer, therefore the amplification matrix $`𝒜`$ is also the Jacobian of the transformation between the source plane and the image plane. We will restrain here our computation to the weak lensing effect so observed quantity will not take into account any other secondary effect. It is very important at this point to note that the lensing effect does not produce any polarization nor rotate the Stokes parameter. In this regime its effect reduces to a simple deformation of the polarization patterns, similar to the temperature maps. This is the mechanism by which the geometrical properties of the polarization field are changed.

To see that we have to consider the *electric* ($`E`$) and *magnetic* ($`B`$) components instead of the Stokes parameters. At small angular scales (we assume that a small fraction of the sky can be described by a plane), these two quantities are defined as,

$`E`$ $``$ $`\mathrm{\Delta }^1\left[\left(_x^2_y^2\right)Q+2_x_yU\right]`$ (7)
$`B`$ $``$ $`\mathrm{\Delta }^1\left[\left(_x^2_y^2\right)U2_x_yQ\right].`$
This fields reflect non-local geometrical properties of the polarization field. The electric component accounts for the scalar part of the polarization and the magnetic one, the pseudo-scalar part: by parity change $`E`$ is conserved, whereas $`B`$ sign is changed. As it has been pointed out in previous papers, lens effects partly redistribute polarization power in these two fields.

We explicit this latter effect in the weak lensing regime where distortions, $`\kappa `$ and $`\gamma _i`$ components are small. This is indeed expected to be the case when lens effects by the large-scale structures are considered, for which the typical value of the convergence field $`\kappa `$ is expected to be $`1\%`$ at 1 degree scale. The leading order effect is obtained by simply pluging (6) in (7) and by expanding the result at leading order in $`\xi `$ , $`\kappa `$, and $`\gamma `$ . Noting that (these calculations are very similar to those done in ),

$`_i\widehat{X}`$ $`=`$ $`\widehat{_kX}\left(\delta _i^k+\xi _{,i}^k\right)`$ (8)
$`_i_j\widehat{X}`$ $`=`$ $`\widehat{_k_lX}(\delta _i^k+\xi _{,i}^k)(\delta _j^l+\xi _{,j}^l)`$
$`+\widehat{_kX}\xi _{,ij}^k`$
we can write a perturbation description of the lensing effect on electric and magnetic components of the polarization. At leading order one obtains:

$`\mathrm{\Delta }\widehat{E}`$ $`=`$ $`\mathrm{\Delta }E+\xi ^i_i\mathrm{\Delta }E2\kappa \mathrm{\Delta }E`$
$`2\delta _{ij}\left(\gamma ^i\mathrm{\Delta }P^j+\gamma _{,k}^iP^{j,k}\right)+O(\gamma ^2)`$
$`\mathrm{\Delta }\widehat{B}`$ $`=`$ $`\mathrm{\Delta }B+\xi ^i_i\mathrm{\Delta }B2\kappa \mathrm{\Delta }B`$
$`2ϵ_{ij}\left(\gamma ^i\mathrm{\Delta }P^j+\gamma _{,k}^iP^{j,k}\right)+O(\gamma ^2),`$
Where we used the fact that $`\widehat{\mathrm{\Delta }X}=\mathrm{\Delta }X+\xi ^i_i\mathrm{\Delta }X`$ at the leading order. The formulas for $`E`$ and $`B`$ are alike. The only difference stands in the $`\delta _{ij}`$ and $`ϵ_{ij}`$ (the latter is the totally antisymmetric tensor, $`ϵ_{11}=ϵ_{22}=0,`$ $`ϵ_{12}=ϵ_{21}=1`$) that reflects the geometrical properties of the two fields. The first three terms of each of these equations represent the naive effect: the lens induced deformation of the $`E`$ or $`B`$ fields. This effect is complemented by an enhancement effect (respectively $`\kappa \mathrm{\Delta }E`$ and $`\kappa \mathrm{\Delta }B`$) and by shear-polarization mixing terms.The latter effects consist in two parts. One which we will call the $`\mathrm{\Delta }`$-term that couples the shear with second derivative of the polarization field. The other one, hereafter the $``$-term, mixes gradient of the shear and polarization. Although terms like $``$ have been neglected in similar computations we cannot do that here a priori. We will indeed show later that these two terms have similar amplitudes.

One consequence of standard inflationary models on CMB anisotropies is the unbalanced distribution of power between the electric ($`E`$) and magnetic ($`B`$) component of its polarization. Adiabatic scalar fluctuations do not induce $`B`$-type polarization and they dominate at small scales over the tensor perturbations (namely the gravity waves). So, even though gravity waves induce $`E`$ and $`B`$ type polarization in a similar amount, *primary* CMB sky is expected to be completely dominated by $`E`$ type polarization at small scales. Then for this class of models the actual magnetic component of the polarization field is generated by the corrective part of eq. (II.1),

$$\mathrm{\Delta }\widehat{B}=2ϵ_{ij}\left(\gamma ^i\mathrm{\Delta }\widehat{P}^j+\gamma _{,k}^i\widehat{P}^{j,k}\right)$$
(10)
This result extends the direct lens effects described in Benabed & Bernardeau who focused their analysis on the lens effect due to the discontinuity of the polarization field in case of cosmic strings. Previous studies of the weak lensing effect on CMB showed that with lensing, the $`B`$ component becomes important at small scales. We obtain here the same result but with a different method; eq. (10) means that the polarization signal $`P`$ is redistributed by the lensing effect in a way that breaks the geometrical properties of the primordial field. Note here that it is mathematically possible to build a shear field that preserves these geometrical properties and that does not create any $`B`$ signal at small scales. We will discuss this problem in Sec. II.3. It also means that $`B`$ directly reflects the properties of the shear map. We will take advantage of this feature to probe the correlation properties of $`B`$ with the projected mass distribution in next sections.

### II.2 Lens-induced $`B`$ maps

We show examples of lens induced $`B`$ maps. These maps have been calculated using “CMBSlow” code developed by A. Riazuelo (see ) to compute primordial polarization maps (we use realizations of standard CDM model to illustrate lens effects). Then various shear maps are applied. We present both true distortions, (obtained by Delaunay triangulation used to shear the $`Q`$ and $`U`$ fields), and the first order calculations given by eq. (10).

Fig. 1 presents the shear effect induced by an isothermal sphere with finite core radius (and the lens edges have been suppressed by an exponential cutoff to minimize numerical noise). The agreement between true distortion (central panel ) and first order formula (right panel) is good. However, a close examination of the maps reveals that some structures in the true map are slightly wider than their counterparts in the first order map. This error is more severe in the center, where the distortion is bigger, which is to be expected since the limits of the validity region of first order calculations are reached.

Fig. 2 shows the $`B`$ field induced by a *realistic* distortion. We use second order Lagrangian dynamics to create a $`2.5\times 2.5`$ degree map that mimics a realistic projected mass density up to $`z=1000`$ and used its gravitational distortion to compute a typical weak lensing-induced $`B`$ map. Again we compare the *exact* effect (i.e. left panel where Delaunay triangulation is used) and the first order formula (middle panel). Right panel shows the difference between the two maps. It reveals the locations where the two significantly disagree. In fact most disagreements are due to slight mismatch of the $`B`$ patch positions, which lead to dipole like effects in this map.

We also show here a comparison of the two parts of the first order formula eq. (10) in order to see which of the $`\mathrm{\Delta }`$ or $``$ terms dominates. It would be more comfortable if one of the two terms was dominant, however, Fig. 3 shows that it is not the case. Even if the $`\mathrm{\Delta }`$-term dominates at low ($`<1000)`$ $`\mathrm{}`$, it is only twice bigger than $``$-one at this scale. The inverse is true for higher ($`30005000)`$ $`\mathrm{}`$s. This can be seen by looking at Fig. 4 where we show the relative amplitudes of the $`\mathrm{\Delta }`$ and $``$ contributions. The $`\mathrm{\Delta }`$ part gives birth to large patches (around $`10^{}`$) while $``$ panel shows a lot more of small features.

### II.3 Direct reconstruction – Kernel problem

The fact that the observable $`B`$ is at leading order proportional to the weak lensing signal invites us to try a direct reconstruction, similar to the lensing mass reconstruction. In fact, we can write

$$\mathrm{\Delta }\widehat{B}=2ϵ_{ij}\left(\gamma ^i\mathrm{\Delta }\widehat{P}^j+\gamma _{,k}^i\widehat{P}^{j,k}\right)\mathrm{F}[\gamma ]$$
(11)
and our reconstruction problem becomes an inversion problem for the operator $`\mathrm{F}`$. Unfortunately, one can prove that this problem has no unique solution. It is due to the fact that $`\mathrm{F}`$ admits a huge kernel, in the sense that, given a polarization map, there is a wide class of shear fields that will conserve a null $`B`$ polarization. The demonstration of this property is sketched in the following.

Since the unlensed polarization is only electric in our approximation, we can describe it by the Laplacian of a scalar field ;

$$E\mathrm{\Delta }\phi \text{ so }\{\begin{array}{c}Q=\left(_x^2_y^2\right)\phi \hfill \\ U=2_x_y\phi \hfill \end{array}.$$
(12)
The same holds for the shear and convergence fields

$$\kappa \frac{\mathrm{\Delta }\psi }{2},\gamma _1=\frac{1}{2}\left(_x^2_y^2\right)\psi \text{ },\gamma _2=_x_y\psi .$$
(13)
Thus we need to know, for a given $`\phi `$ field, whether there is any $`\psi `$ that fulfills the equation

$$\gamma _2\mathrm{\Delta }Q\gamma _1\mathrm{\Delta }U+_i\gamma _2^iQ_i\gamma _1^iU=0.$$
(14)
$`\phi `$ and $`\psi `$ can be written as polynomial decompositions

$`\phi (x,y)`$ $`=`$ $`{\displaystyle \underset{n,l}{}}a_{nl}x^ny^l`$
$`\psi (x,y)`$ $`=`$ $`{\displaystyle \underset{m,k}{}}b_{mk}x^my^k.`$ (15)
Using (15) in (14) we are left with a new polynomial whose coefficients $`c_{ij}`$ are sums of $`a_{nl}\times b_{mk}`$ and have to be all put to zero. With the coefficient equations in hand, it is easy to prove that assuming all the $`b_{mk}`$ coefficient up to $`m+k=N`$ are known and writing the equations $`i+j=(N+1)3,c_{ij}=0`$, we can compute out of all the $`a_{nl}`$ all but three $`b_{mk}`$ with $`m+k=N+1`$. This is somewhat similar to mass reconstruction problems from galaxy surveys where one cannot avoid the mass sheet degeneracy. The situation is however worse in our case since not only constant convergence but also translations and a whole class of $`a_{nk}`$ realization dependent complex deformations are indiscernible. Thus, with the only knowledge of the $`B`$ component of the polarization one cannot, with the first order eq. (10), recover the projected mass distribution.

## III Cross-correlating CMB maps and weak lensing surveys

### III.1 Motivations

Even with the most precise experiments it is clear that clean detections of $`B`$ component will be difficult to obtain. The magnetic polarization amplitude induced with such a mechanism is expected to be one order of magnitude below the electric one. Besides even if we know that there is a window in angular scale where the other secondary effects will not interfere too much with the detection of the lens-induced $`B`$, few is known about removing the foregrounds to obtain clean maps reconstruction algorithms would require.

These considerations lead us to look for complementary data sets to compare $`B`$ with. Although the source plane for weak lensing surveys is much closer than for the lensed CMB fluctuations, we expect to have a significant overlapping region in the two redshift lens distributions, so that weak lensing surveys can map a fair fraction of the line-of-sight CMB lenses. Consequently, weak lensing surveys can potentially provide us with shear maps correlated with $`B`$, but which have different geometrical degeneracy, noise sources and systematics than the polarization field.

The correlation strength between the lensing effects at two different redshifts can be evaluated. We define $`r`$ as the cross-correlation coefficient between two lens planes:

$$r(z_{\text{gal}})=\frac{\kappa \kappa _{\text{gal}}}{\sqrt{\kappa ^2\kappa _{\text{gal}}^2}}.$$
(16)
In a broad range of realistic cases (see tab. 1), $`r40\%`$. To take advantage of this large overlapping we will consider quantity that cross correlates the CMB $`B`$ field and galaxy surveys. Moreover, cross-correlation observations are expected to be insensitive to noises in weak lensing surveys and in CMB polarization maps. This idea has already been explored for temperature maps. We extend this study here taking advantage of the specific geometrical dependences uncovered in the previous section.

### III.2 Definition of $`b_\mathrm{\Delta }`$ and $`b_{}`$.

The magnetic component of the polarization in eq. (10) appears to be built from a pure CMB part, which comes from the primordial polarization, and a gravitational lensing part. It is natural to define $`b`$, in such a way that mimics the $`\mathrm{\Delta }\widehat{B}`$ fonction dependance, by replacing the CMB shear field by the galaxy one.

$`b`$ $`=`$ $`ϵ_{ij}\left(\gamma _{\text{gal}}^i\mathrm{\Delta }\widehat{P}^j+\gamma _{\text{gal},k}^i\widehat{P}^{j,k}\right)`$
$`=`$ $`ϵ_{ij}\left(\gamma _{\text{gal}}^i\mathrm{\Delta }P^j+\gamma _{\text{gal},k}^iP^{j,k}\right)+O(\kappa ^2).`$
In the following, we will label local lensing quantities, such as what one can obtain from lensing reconstruction on galaxy surveys, with a gal index. This new quantity can be viewed as a guess for the CMB polarization $`B`$ component if lensing was turned on only in a redshift range matching the depth of galaxy surveys. The correlation coefficient of this guess with the true $`\mathrm{\Delta }B`$ field, that is $`\mathrm{\Delta }\widehat{B}b`$, is expected to be quadratic both in $`P`$ and in $`\gamma `$ and to be proportional to the cross-coefficient $`r`$.

For convenience, and in order to keep the objects we manipulate as simple as possible, we will not exactly implement this scheme, as it will lead to uneven angular derivative degrees in the two terms of resulting equations. We can, instead, decompose the effect in the $`\mathrm{\Delta }`$ and $``$-part. These two are not correlated, since their components do not share the same degrees of angular derivation<sup>1</sup><sup>1</sup>1generically, a random field and its derivative at the same point are not correlated. . Hence, we can play the proposed game, considering the two terms of eq. (10) as if they were two different fields, creating two guess-quantities that should correlate independently with the observed $`B`$ field. Following this idea we build $`b_\mathrm{\Delta }`$ as,

$`b_\mathrm{\Delta }`$ $``$ $`ϵ_{ij}\gamma _{\text{gal}}^i\mathrm{\Delta }\widehat{P}^j`$
$`=`$ $`ϵ_{ij}\gamma _{\text{gal}}^i\mathrm{\Delta }P^j+O(\kappa ^2)`$
which corresponds to the $`\mathrm{\Delta }`$-term in eq. (10). The amplitude of the cross-correlation between $`\mathrm{\Delta }B`$ and $`b_\mathrm{\Delta }`$ can easily be estimated. At leading order, we have

$$\mathrm{\Delta }\widehat{B}b_\mathrm{\Delta }=2ϵ_{ij}ϵ_{kl}\gamma ^k\gamma _{\text{gal}}^i\mathrm{\Delta }P^l\mathrm{\Delta }P^j.$$
(19)
The corresponding $``$ correlation is

$$\mathrm{\Delta }\widehat{B}b_{}=2ϵ_{ij}ϵ_{kl}_m\gamma ^k_n\gamma _{\text{gal}}^i_mP^l_nP^j$$
(20)
where we have defined

$$b_{}ϵ_{ij}_k\gamma _{\text{gal}}^i_k\widehat{P}^j.$$
(21)
Fig. 4 shows numerical simulations presenting maps of first order $`\mathrm{\Delta }\widehat{B}`$, its $`\mathrm{\Delta }`$ and $``$ contributions and the corresponding guess maps one can build with a low $`z`$ shear map. The similarities between the top maps and the bottom maps are not striking. Yet, under close examination one can recognize individual patterns shared between the maps. This is confirmed by the computation of the correlation coefficient between the maps, that shows significant overlapping, between 50% and 15%, depending correlation and filtering strategy. The calculations hereafter will evaluate the theoretical correlation structure between maps given in figs. 4-b and 4-g & h.

For galaxy surveys, the amplification matrix is,

$`𝒜_{\text{gal}}^1(\stackrel{}{\alpha })\text{Id}={\displaystyle _0^{z_{\text{gal}}}}d\chi w_{\text{gal}}(\chi )`$ (22)
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^{\frac{3}{2}}}}\delta (\stackrel{}{k})\mathrm{e}^{\mathrm{i}\left(k_r\chi +\stackrel{}{k}_{}D(\chi )\stackrel{}{\alpha }\right)}`$
$`\times \left(\begin{array}{cc}1+\mathrm{cos}(2\varphi _k_{})& \mathrm{sin}(2\varphi _k_{})\\ \mathrm{sin}(2\varphi _k_{})& 1\mathrm{cos}(2\varphi _k_{})\end{array}\right)`$ (25)
where $`\delta (k)`$ is the Fourier transform of the density contrast at redshift $`z(\chi )`$, $`w`$ is the lens efficiency function, $`\mathrm{D}`$ is the angular distance, and $`\varphi _k_{}`$is the position angle of the transverse wave-vector $`k_{}`$in the $`k_{}=(k_x,k_y)`$ plane. Assuming a Dirac source distribution the efficiency function is given by

$$w_{\text{gal}}(z)=\frac{3}{2}\mathrm{\Omega }_\mathrm{o}\frac{\mathrm{D}_z\mathrm{D}_{zz_{\text{gal}}}}{a\mathrm{D}_{z_{\text{gal}}}}.$$
(26)
Note that the Fourier components $`\delta (k)`$ include the density time evolution. They are thus proportional to the growth factor in the linear theory. The time evolution of these components is much more complicated in the nonlinear regime (see ).

Then, $`b_{\mathrm{}}`$ is

$$b_{\mathrm{}}(\stackrel{}{\alpha })=^{\chi _{\text{gal}}}𝒟(\chi ,\stackrel{}{l},\stackrel{}{k})\stackrel{~}{E}(l)\delta (k)𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k}_{})$$
(27)
with the integration element defined as,

$$𝒟(\chi ,\stackrel{}{l},\stackrel{}{k})=\mathrm{d}\chi w_{\text{gal}}(\chi )\frac{\mathrm{d}^3k}{(2\pi )^{3/2}}\frac{\mathrm{d}^2l}{2\pi }\mathrm{e}^{\mathrm{i}\left[k_r\chi +\left(\stackrel{}{k}_{}\mathrm{D}(\chi )+\stackrel{}{l}\right)\stackrel{}{\alpha }\right]},$$

(it actually depends on the position of the source plane through the efficiency function $`w(z)`$) and where $`\mathrm{}`$ stands for either $`\mathrm{\Delta }`$ or $``$. The geometrical kernel $`𝒢^{\mathrm{Ker}}`$ is given by (using eq. (12))

$`𝒢_\mathrm{\Delta }^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k})`$ $``$ $`l^2\mathrm{sin}2\left(\varphi _k\varphi _l\right)`$ (28)
$`𝒢_{}^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k})`$ $``$ $`lk\mathrm{cos}\left(\varphi _k\varphi _l\right)\mathrm{sin}2\left(\varphi _k\varphi _l\right).`$ (29)
This function contains all the geometrical structures of the $`\mathrm{\Delta }`$ and $``$ terms. We can write the same kind of equation for $`\mathrm{\Delta }\widehat{B}`$. Then, the cross-correlation is

$`\mathrm{\Delta }\widehat{B}b_{\mathrm{}}(\stackrel{}{\alpha })`$ $`=`$ $`2{\displaystyle ^{\chi _{\text{gal}}}}𝒟(\chi _{\text{gal}},\stackrel{}{l}_{\text{gal}},\stackrel{}{k}_{\text{gal}})`$
$`\times {\displaystyle ^{\chi _{\text{cmb}}}}𝒟(\chi _{\text{cmb}},\stackrel{}{l}_{\text{cmb}},\stackrel{}{k}_{\text{cmb}})𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{gal}},\stackrel{}{k}_{\text{gal}})`$
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{cmb}},\stackrel{}{k}_{\text{cmb}})\delta (\stackrel{}{k}_{\text{gal}})\delta (\stackrel{}{k}_{\text{cmb}})\stackrel{~}{E}(\stackrel{}{l}_{\text{gal}})\stackrel{~}{E}(\stackrel{}{l}_{\text{cmb}}).`$
The completion of this calculation requires the use of the small angle approximation,

$`\delta (\stackrel{}{k}_{\text{gal}})\delta (\stackrel{}{k}_{\text{cmb}})`$ $`=`$ $`P(k)\delta ^3(k_{\text{gal}}+k_{\text{cmb}})`$
$`P(k_{})\delta ^2(k_{\text{gal}_{}}+k_{\text{cmb}_{}})\delta (k_{\text{gal}_r}+k_{\text{cmb}_r})`$
which implies

$$\stackrel{}{k}_{\text{gal}}=\stackrel{}{k}_{\text{cmb}}=\stackrel{}{k}$$
(32)
and after the radial components have been integrated out,

$$\chi _{\text{gal}}=\chi _{\text{cmb}}=\chi .$$
(33)
We also define the $`C_E(l)`$ as the angular power spectrum of the $`E`$ field,

$$\stackrel{~}{E}(\stackrel{}{l}_{\text{gal}})\stackrel{~}{E}(\stackrel{}{l}_{\text{cmb}})=C_E(l)\delta ^2(l_{\text{gal}}l_{\text{cmb}})$$
(34)
Eventually one gets,

$`\mathrm{\Delta }\widehat{B}b_{\mathrm{}}(\stackrel{}{\alpha })`$ $`=`$ $`2{\displaystyle ^{z_{\text{gal}}}}d\chi w_{\text{gal}}w_{\text{cmb}}{\displaystyle \frac{\mathrm{d}^2k\mathrm{d}^2l}{(2\pi )^4}}`$
$`\times C_E(l)P(k)𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k})^2`$
Then, integrating on the geometrical dependencies in $`𝒢_{\mathrm{}}^{\mathrm{Ker}}`$, we have

$`\mathrm{\Delta }\widehat{B}b_\mathrm{\Delta }(\stackrel{}{\alpha })`$ $`=`$ $`2{\displaystyle ^{z_{\text{gal}}}}d\chi w_{\text{gal}}w_{\text{cmb}}`$ (36)
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}k\mathrm{d}l}{2(2\pi )^2}}kl^5C_E(l)P(k)`$
$`=`$ $`\mathrm{\Delta }E^2\kappa \kappa _{\text{gal}},`$
and

$`\mathrm{\Delta }\widehat{B}b_{}(\stackrel{}{\alpha })`$ $`=`$ $`{\displaystyle ^{z_{\text{gal}}}}d\chi w_{\text{gal}}w_{\text{cmb}}`$ (37)
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}k\mathrm{d}l}{2(2\pi )^2}}k^3l^3C_E(l)P(k)`$
$`=`$ $`{\displaystyle \frac{1}{2}}(\stackrel{}{}E)^2\stackrel{}{}\kappa \stackrel{}{}\kappa _{\text{gal}},`$
implying that, ignoring filtering effects, we are able to measure directly the correlation between lensing effect at $`z_{\text{cmb}}`$ and any $`z_{\text{gal}}`$ a weak lensing survey can access. Since $`\mathrm{\Delta }\widehat{E}=\mathrm{\Delta }E(1+O(\kappa ))`$ we get, for the $`\mathrm{\Delta }`$ type quantity,

$`\mathrm{\Delta }\widehat{E}^2`$ $`=`$ $`\mathrm{\Delta }E^2\left(1+O(\kappa )\right)^2`$
$`=`$ $`\mathrm{\Delta }E^2\left(1+O\left(\kappa ^2\right)\right).`$
The same holds for $`.`$ We are then able to construct two quantities insensitive to the normalization of CMB and $`\sigma _8`$

$`𝒳_\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }\widehat{B}b_\mathrm{\Delta }(\stackrel{}{\alpha })}{\mathrm{\Delta }\widehat{E}^2\kappa _{\text{gal}}^2}}={\displaystyle \frac{\kappa \kappa _{\text{gal}}}{\kappa _{\text{gal}}^2}}`$
$``$ $`r\sqrt{{\displaystyle \frac{\kappa ^2}{\kappa _{\text{gal}}^2}}}.`$
and

$`𝒳_{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }\widehat{B}b_{}(\stackrel{}{\alpha })}{(\stackrel{}{}\widehat{E})^2(\stackrel{}{}\kappa _{\text{gal}})^2}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\stackrel{}{}\kappa \stackrel{}{}\kappa _{\text{gal}}}{\stackrel{}{}\kappa _{\text{gal}}^2}}`$
$``$ $`{\displaystyle \frac{1}{2}}r_{}\sqrt{{\displaystyle \frac{\kappa ^2}{\kappa _{\text{gal}}^2}}}.`$
We implicitly defined $`r_{}`$ like $`r`$ but with $`\kappa `$ instead of $`\kappa `$

$$r_{}(z_{\text{gal}})=\frac{\stackrel{}{}\kappa \stackrel{}{}\kappa _{\text{gal}}}{\sqrt{(\kappa )^2(\kappa _{\text{gal}})^2}}.$$
(41)
We will see in Sect. III.4 that they behave very much alike. This result is to be compared with the formula for $`\mathrm{cos}(\theta _g)`$ established in  where the obtained quantity was going like $`r\sqrt{\kappa ^2}`$. These calculations however have neglected the filtering effects that may significantly affect our conclusions. These effects are investigated in next section.

### III.3 Filtering effects

In above section we conduct our calculations assuming no filtering. Obviously we have to take it into account! We will show here that the results obtained before hold, in certain limits, when one adds filtering effects.

In the following, we consider, for simplicity, top-hat filters only. It is expected that other window functions will show very similar behaviors and this simplification does not restrain the generality of our results. Let us call $`W(x)`$ the top-hat filter function in Fourier space

$$W(x)2\frac{\mathrm{J}_1(x)}{x}.$$
(42)
$`\mathrm{J}_1`$ is the first $`\mathrm{J}`$-Bessel function. We will also define $`W_i(x)`$ a general function

$$W_i(x)2\frac{\mathrm{J}_i(x)}{x}$$
(43)
where $`\mathrm{J}_i`$ is the $`i^{\mathrm{th}}`$ $`\mathrm{J}`$-Bessel function, so that $`W=W_1`$. Then, if $`X(\stackrel{}{\alpha })`$ is the value of any quantity $`X`$ at position $`\stackrel{}{\alpha }`$ on the sky, its top-hat filtered value can be computed as

$$X_{(\theta )}(\stackrel{}{\alpha })=\frac{\mathrm{d}^2k}{2\pi }\stackrel{~}{X}_kW(k\theta )\mathrm{e}^{\mathrm{i}\stackrel{}{k}\stackrel{}{\alpha }},$$
(44)
where $`\stackrel{~}{X}`$ is $`X`$ Fourier transform. In the following we will note $`X_{(\theta )}`$ the filtered quantity at scale $`\theta `$.

The tricky thing for $`\mathrm{\Delta }\widehat{B}b_{\mathrm{}}`$ is that the CMB part and the low-redshift weak lensing part are *a priori* filtered at different scale. For $`\mathrm{\Delta }\widehat{B}`$, which is a measured value, its pure CMB part and its weak lensing part are filtered at the same scale $`\theta `$. Hence, $`\widehat{B}`$ reads,

$`\mathrm{\Delta }\widehat{B}(\stackrel{}{\alpha })_{(\theta )}`$ $`=`$ $`2{\displaystyle ^{\chi _{\text{cmb}}}}𝒟(\chi ,\stackrel{}{l},\stackrel{}{k})\stackrel{~}{E}(l)\delta (k)`$ (45)
$`\times \left[𝒢_\mathrm{\Delta }^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k}_{})+𝒢_\mathrm{\Delta }^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k}_{})\right]W\left(|\stackrel{}{k}_{}\mathrm{D}+\stackrel{}{l}|\theta \right)`$
*A contrario* $`b_{\mathrm{}}`$ is a composite value. The CMB part is still filtered at $`\theta `$ whereas the weak lensing part (which comes from a weak lensing survey of galaxies) is filtered independently at another scale which we denote $`\theta _{\text{gal}}`$. It implies that,

$`b_{\mathrm{}}(\stackrel{}{\alpha })_{(\theta )}`$ $`=`$ $`2{\displaystyle ^{\chi _{\text{gal}}}}𝒟(\chi ,\stackrel{}{l},\stackrel{}{k})\stackrel{~}{E}(l)\delta (k)`$ (46)
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k}_{})W(k\mathrm{D}\theta _{\text{gal}})W(l\theta ).`$
Taking filtering into account, the cross-correlation coefficient becomes,

$`\mathrm{\Delta }\widehat{B}_{(\theta )}b_{\mathrm{}(\theta ,\theta _{\text{gal}})}`$ $`=`$ $`2{\displaystyle ^{z_{\text{gal}}}}d\chi w_{\text{gal}}w_{\text{cmb}}`$
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}^2k\mathrm{d}^2l}{(2\pi )^4}}C_E(l)P(k)𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l},\stackrel{}{k})`$
$`\times W(k\mathrm{D}\theta _{\text{gal}})W(l\theta )W(|\stackrel{}{k}\mathrm{D}+\stackrel{}{l}|\theta ).`$
It can be shown (from the summation theorems of the Bessel functions) that,

$`W_1(|\stackrel{}{k}\mathrm{D}+\stackrel{}{l}|\theta )=`$ (48)
$`{\displaystyle \underset{i=1}{}}iW_i(k\mathrm{D}\theta )W_i(l\theta )(1)^i{\displaystyle \frac{\mathrm{sin}i(\varphi _k\varphi _l)}{\mathrm{sin}(\varphi _k\varphi _l)}}`$
It is then possible to break the $`W\left(|\stackrel{}{k}\mathrm{D}+\stackrel{}{l}|\theta \right)`$ into a sum of $`W_i(k\mathrm{D}\theta )W_i(l\theta )`$ with coefficients that depend on the geometrical properties of our problem. Integrating over the geometrical dependencies of $`𝒢_{\mathrm{}}^{\mathrm{Ker}}`$, leaves us with only a few non vanishing terms in our sum,

$$d\varphi \mathrm{sin}^2(2\varphi )\frac{\mathrm{sin}(i\varphi )}{\mathrm{sin}\varphi }=\{\begin{array}{cc}\pi \hfill & i\text{ }=\text{ }1\text{ or }i=3\hfill \\ 0\hfill & \text{elsewhere}\hfill \end{array},$$
(49)
for the $`\mathrm{\Delta }`$-term and

$$d\varphi \mathrm{cos}\varphi \mathrm{sin}^2(2\varphi )\mathrm{sin}(i\varphi )\mathrm{sin}\varphi =\{\begin{array}{cc}\pi /2& i=1\hfill \\ 3\pi /4& i=3\hfill \\ \pi /4& i=5\hfill \\ 0& \text{elsewhere}\hfill \end{array},$$
(50)
for the $``$-term. Each term can be computed exactly, and it turns out that the terms built from $`W_i,i>1`$ are always negligible compared to the ones coming from $`W_1`$. It implies that we can safely ignore the $`W_3`$ and $`W_5`$ in both $`\mathrm{\Delta }`$ and $``$ expressions, therefore it is reasonable to assume that $`W\left(|\stackrel{}{k}\mathrm{D}+\stackrel{}{l}|\theta \right)=W(k\mathrm{D}\theta )W(l\theta ).`$ It is expected that other windows, in particular the Gaussian window function, share similar properties. Then, taking into accounts the filtering effects, the equations for the cross-correlations reduce to

$$\mathrm{\Delta }\widehat{B}_{(\theta )}b_{\mathrm{\Delta }(\theta ,\theta _{\text{gal}})}=\mathrm{\Delta }E_{(\theta )}^2\kappa _{(\theta )}\kappa _{\text{gal}(\theta _{\text{gal}})}$$
(51)
and

$$\mathrm{\Delta }\widehat{B}_{(\theta )}b_{(\theta ,\theta _{\text{gal}})}=\frac{1}{2}E_{(\theta )}^2\kappa _{(\theta )}\kappa _{\text{gal}(\theta _{\text{gal}})},$$
(52)
so that our correlation coefficients can be written,

$$𝒳_{\mathrm{\Delta }(\theta ,\theta _{\text{gal}})}=r_{(\theta ,\theta _{\text{gal}})}\sqrt{\frac{\kappa _{(\theta )}^2}{\kappa _{\text{gal}(\theta _{\text{gal}})}^2}}$$
(53)
and

$$𝒳_{(\theta ,\theta _{\text{gal}})}=\frac{1}{2}r_{(\theta ,\theta _{\text{gal}})}\sqrt{\frac{\kappa _{(\theta )}^2}{\kappa _{\text{gal}(\theta _{\text{gal}})}^2}}.$$
(54)
The results obtained in eqs. (III.2-III.2) are thus still formally valid. Actually, eqs. (53-54) simply tell that filtering effects can simply be assumed to act independently on the lensing effects and on the primary Cosmic Microwave Background maps. We are left with two quantities that only reflects the line-of-sight overlapping effects of lensing distortions.

### III.4 Sensitivity to the cosmic parameters

We quickly explore here the behavior of $`𝒳_{\mathrm{}}`$ in different sets of cosmological parameters. These quantities only depend on weak lensing quantities. Ignoring the $`\mathrm{\Omega }_0`$ dependence in the angular distances and growing factor, one would expect $`\kappa ^2`$ to scale like $`\mathrm{\Omega }_0^2`$. Yet, because of the growth factor, the convergence field exhibits a weaker sensitivity to $`\mathrm{\Omega }_0`$. Assuming $`\mathrm{\Lambda }=0`$ and a power law spectrum, we know from  that $`\kappa _{\text{gal}}^2\mathrm{\Omega }_0^{1.66}`$ for $`z_{\text{gal}}=1`$. The same calculation leads to $`\kappa _{\text{cmb}}\kappa _{\text{gal}}\mathrm{\Omega }_0^{1.68}`$, $`(\kappa _{\text{gal}})^2\mathrm{\Omega }_0^{1.91}`$ and $`\stackrel{}{}\kappa _{\text{cmb}}\stackrel{}{}\kappa _{\text{gal}}\mathrm{\Omega }_0^{1.915}`$. Then, in this limit, the quantities $`𝒳_{\mathrm{}}`$ have a very low dependence on $`\mathrm{\Omega }_0`$ :

$$𝒳_\mathrm{\Delta }\mathrm{\Omega }_0^{0.02}\text{ and }𝒳_{}\mathrm{\Omega }_0^{0.005}.$$

Eventually, the $`𝒳_{\mathrm{}}`$ quantities should exhibit a seizable sensitivity to $`\mathrm{\Lambda }`$; changing $`\mathrm{\Lambda }`$ increases or reduces the size of the optic bench and accordingly the overlapping between $`\kappa _{\text{cmb}}`$ and $`\kappa _{\text{gal}}`$.

Figs. 5 and 6 present contour plots of the amplitude of $`𝒳_\mathrm{\Delta }`$ and $`𝒳_{}`$ in the $`(\mathrm{\Omega }_0,\mathrm{\Lambda })`$ plane for CDM models. They show the predicted low $`\mathrm{\Omega }_0`$ sensitivity and the expected $`\mathrm{\Lambda }`$ dependency. Both figures are very alike. This is due to the fact that the dominant features are contained in the efficiency function dependences on the angular distances.

### III.5 Cosmic variance

In previous sections we looked at the sensitivity of observable quantities which mixed galaxy weak lensing surveys and CMB polarization detection. It is very unlikely that both surveys will be able to cover, with a good resolution and low foreground contamination, a large fraction of the sky. It seems however reasonable to expect to have at our disposal patches of at least a few hundreds square degrees. The issue we address in this section is to estimate the cosmic variance of such a detection in joint surveys in about 100 square degrees.

The computation of cosmic variance is a classical problem in cosmological observation . A natural estimate for an ensemble average $`X`$ is its geometrical average. If the survey has size $`\mathrm{\Sigma }`$ then,

$$\overline{X}=\frac{1}{\mathrm{\Sigma }}_\mathrm{\Sigma }\mathrm{d}^2\alpha X(\stackrel{}{\alpha })$$
(55)
For a compact survey with circular shape of radius $`\mathrm{\Xi }`$ we formally have,

$$\overline{X}=\frac{\mathrm{d}^2k}{2\pi }\stackrel{~}{X}(\stackrel{}{k})W(k\mathrm{\Xi }).$$
(56)
For sake of simplicity this is what we use in the following but we will see that the shape of the survey has no significant consequences.

Taking $`\overline{X}`$ as an estimate of $`X`$ (the ensemble average of $`X`$), leads to an error of the order $`\sqrt{\overline{X}^2\overline{X}^2}`$which usually scales like $`1/\sqrt{\mathrm{\Sigma }}`$ if the survey is large enough.

When we are measuring $`𝒳_{\mathrm{}}`$ on a small patch of the sky, we are apart from the statistical value by the same kind of error. We can neglect the errors on $`\mathrm{\Delta }\widehat{E}^2`$, $`(\widehat{E})^2`$, $`(\kappa _{\text{gal}})^2`$and $`\kappa _{\text{gal}}^2`$; those may not be the dominant source of discrepancy and can even be measured on wider and independent samples. The biggest source of error is the measure of $`\mathrm{\Delta }\widehat{B}b_{\mathrm{}}`$. It is given by,

$$C_{\mathrm{}}=\sqrt{\left(\overline{\mathrm{\Delta }\widehat{B}b_{\mathrm{}}}\overline{\mathrm{\Delta }\widehat{B}}\overline{b_{\mathrm{}}}\right)^2\overline{\mathrm{\Delta }\widehat{B}b_{\mathrm{}}}\overline{\mathrm{\Delta }\widehat{B}}\overline{b_{\mathrm{}}}^2}.$$
(57)
The computation of (57) is made easier if we write explicitly the geometrical average as a summation over $`N`$ measurement points ($`N`$ can be as large as we want),

$$\overline{X}=\frac{1}{N}\underset{i=1}{\overset{N}{}}X(\theta _i),$$
(58)
we then developed (57), and replace the ensemble average of the summation sign by the geometrical average over the survey size. We are left with a sum of correlators containing $`8`$ fields taken at 2, 3 and 4 different points. The calculations can be carried out analytically if we assume that all our fields follow Gaussian statistics, which is reasonable at the scale we are working on. In that case indeed, we can take advantage of the Wick theorem to contract each of the $`8`$ fields correlators in products of $`2`$ points correlation functions. By definition, (57) contains only connected correlators, moreover the ensemble averages $`\mathrm{\Delta }\widehat{B}`$ and $`b_{\mathrm{}}`$ vanish, therefore only a small fraction of correlators among all the possible combination of the $`8`$ fields survive. We can use a simple diagrammatic representations to describe their geometrical shape. All the non vanishing terms in $`C_{\mathrm{}}`$ are given in Fig. 7. Each line between two vertex represents a $`2`$ points correlation function such as $`X(\stackrel{}{\alpha }_1)X(\stackrel{}{\alpha }_2)`$, and the different symbols at the vertex correspond to different $`X`$ fields (the cross stands for $`\mathrm{\Delta }P`$, the dot for $`\gamma _{\text{cmb}}`$, and the open dot stands for $`\gamma _{\text{gal}}`$). The $`𝒜`$-terms represent terms where the two top (and the two bottom) $`\mathrm{\Delta }B`$ and $`b_{\mathrm{}}`$ are taken at the same point, but top and bottom fields are not at the same place. The $``$-terms are three points diagrams: the top $`\mathrm{\Delta }B`$ and $`b_{\mathrm{}}`$ are at the same point whereas the right and left bottom vertexes are at two different locations. The $`𝒞`$ terms are four-points diagrams, where each vertex is at a different point. To illustrate our notations, let us write $`_{2c}^{\mathrm{}}`$ as an example,

$$\begin{array}{ccc}_{2c}^{\mathrm{}}& =& \gamma _{\text{cmb}}(\stackrel{}{\alpha }_1)\gamma _{\text{gal}}(\stackrel{}{\alpha }_2)\gamma _{\text{gal}}(\stackrel{}{\alpha }_3)\gamma _{\text{cmb}}(\stackrel{}{\alpha }_1)\times \\ & & \mathrm{\Delta }P(\stackrel{}{\alpha }_1)\mathrm{\Delta }P(\stackrel{}{\alpha }_1)\mathrm{\Delta }P(\stackrel{}{\alpha }_2)\mathrm{\Delta }P(\stackrel{}{\alpha }_3)\end{array}$$

We only focus on the calculation of the $`𝒜`$ terms because we can use the approximation that

$$𝒜𝒞.$$
(59)
Indeed, in perturbative theory, if the survey is large enough, the $`n`$-points correlation functions naturally dominates over the $`n+1`$-points correlation function. This is true as long as the local variance is much bigger than the autocorrelation at survey scale and we assume the surveys are still large enough to be in this case.

The general expression for any $`𝒜`$ diagram is

$`𝒜_i^{\mathrm{}}`$ $`=`$ $`4{\displaystyle ^{\text{cmb}}}𝒟(\chi _{\text{cmb}1},\stackrel{}{l}_{\text{cmb}1},\stackrel{}{k}_{\text{cmb}1})𝒟(\chi _{\text{cmb}2},\stackrel{}{l}_{\text{cmb}2},\stackrel{}{k}_{\text{cmb}2})`$
$`\times {\displaystyle ^{\text{gal}}}𝒟(\chi _{\text{gal}1},\stackrel{}{l}_{\text{gal}1},\stackrel{}{k}_{\text{gal}1})𝒟(\chi _{\text{gal}2},\stackrel{}{l}_{\text{gal}2},\stackrel{}{k}_{\text{gal}2})`$
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{cmb}1},\stackrel{}{k}_{\text{cmb}1_{}})𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{cmb}2},\stackrel{}{k}_{\text{cmb}2_{}})`$
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{gal}1},\stackrel{}{k}_{\text{gal}1_{}})𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{gal}2},\stackrel{}{k}_{\text{gal}2_{}})_i\stackrel{}{k}_i\stackrel{}{l}_j`$
$`\times W(|\stackrel{}{k}_{\text{cmb}1_{}}\mathrm{D}+\stackrel{}{l}_{\text{cmb}1}|\theta )W(|\stackrel{}{k}_{\text{cmb}2_{}}\mathrm{D}+\stackrel{}{l}_{\text{cmb}2}|\theta )`$
$`\times W(k_{\text{gal}1_{}}\mathrm{D}_1\theta _{\text{gal}})W(l_{\text{gal}1}\theta )W(k_{\text{gal}2_{}}\mathrm{D}_2\theta _{\text{gal}})W(l_{\text{gal}2}\theta )`$
$`\times W\left(|\stackrel{}{k}_{\text{gal}1_{}}\mathrm{D}_1+\stackrel{}{l}_{\text{gal}1}+\stackrel{}{k}_{\text{cmb}1_{}}\mathrm{D}_1+\stackrel{}{l}_{\text{cmb}1}|\mathrm{\Xi }\right)`$
$`\times W\left(|\stackrel{}{k}_{\text{gal}2_{}}\mathrm{D}_2+\stackrel{}{l}_{\text{gal}2}+\stackrel{}{k}_{\text{cmb}2_{}}\mathrm{D}_2+\stackrel{}{l}_{\text{cmb}2}|\mathrm{\Xi }\right)`$
where $`_i`$ gives the 2-point correlations associated with the lines of the diagram. For example :

$`_1`$ $`=`$ $`\delta (\stackrel{}{k}_{\text{gal}1})\delta (\stackrel{}{k}_{\text{cmb}1})\delta (\stackrel{}{k}_{\text{gal}2})\delta (\stackrel{}{k}_{\text{cmb}2})`$ (61)
$`\times \stackrel{~}{E}(l_{\text{gal}1})\stackrel{~}{E}(l_{\text{gal}2})\stackrel{~}{E}(l_{\text{cmb}1})\stackrel{~}{E}(l_{\text{cmb}2})`$
We explicit in the following the computation of $`𝒜_1^{\mathrm{}}`$. The other terms follow the same treatment or can be neglected. The lines in the $`𝒜_1^{\mathrm{}}`$ diagram give us the relations

$`\stackrel{}{k}_{\text{cmb}1}`$ $`=`$ $`\stackrel{}{k}_{\text{gal}1}=\stackrel{}{k}_1`$
$`\stackrel{}{k}_{\text{cmb}2}`$ $`=`$ $`\stackrel{}{k}_{\text{gal}2}=\stackrel{}{k}_2`$ (62)
$`\stackrel{}{l}_{\text{cmb}1}`$ $`=`$ $`\stackrel{}{l}_{\text{cmb}2}=\stackrel{}{l}_{\text{cmb}}`$
$`\stackrel{}{l}_{\text{gal}1}`$ $`=`$ $`\stackrel{}{l}_{\text{gal}2}=\stackrel{}{l}_{\text{gal}}`$
Then, using these relations and the small angular approximation, we have :

$`𝒜_1^{\mathrm{}}`$ $`=`$ $`4{\displaystyle ^{\text{gal}}}d\chi _1d\chi _2w_{\text{cmb}1}w_{\text{gal}1}w_{\text{cmb}2}w_{\text{gal}2}`$
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}^2k_1\mathrm{d}^2k_2}{(2\pi )^4}}{\displaystyle \frac{\mathrm{d}^2l_{\text{gal}}\mathrm{d}^2l_{\text{cmb}}}{(2\pi )^4}}`$
$`\times C_E(l_{\text{gal}})C_E(l_{\text{cmb}})P(k_1)P(k_2)`$
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{cmb}},\stackrel{}{k}_1)𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{cmb}},\stackrel{}{k}_2)`$
$`\times 𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{gal}},\stackrel{}{k}_1)𝒢_{\mathrm{}}^{\mathrm{Ker}}(\stackrel{}{l}_{\text{gal}},\stackrel{}{k}_2)`$
$`\times W(|\stackrel{}{k}_1\mathrm{D}+\stackrel{}{l}_{\text{cmb}}|\theta )W(|\stackrel{}{k}_2\mathrm{D}+\stackrel{}{l}_{\text{cmb}}|\theta )`$
$`\times W(k_1\mathrm{D}\theta _{\text{gal}})W(k_2\mathrm{D}\theta _{\text{gal}})W^2(l_{\text{gal}}\theta )`$
$`\times W^2(|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|\mathrm{\Xi }).`$
We apply the decomposition of $`W_1\left(|\stackrel{}{k}\mathrm{D}(\chi )+\stackrel{}{l}|\theta \right)`$ we used in eq. (48). The geometry of our problem is the same and the result (49) still holds for the terms in $`W_1\left(|\stackrel{}{k}_1\mathrm{D}(\chi _1)+\stackrel{}{l}_{\text{cmb}}|\theta \right)`$ and $`W_1\left(|\stackrel{}{k}_2\mathrm{D}(\chi _2)+\stackrel{}{l}_{\text{cmb}}|\theta \right)`$. This however is not true for $`W_1^2\left(|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|\mathrm{\Xi }\right)`$ for which the application of the re-summation theorem does not bring any simplification. Then, neglecting all the $`W_3`$ parts and after integration on the $`\varphi _{k_i}`$, for the $`\mathrm{\Delta }`$-term, we have,

$`𝒜_1^\mathrm{\Delta }`$ $`=`$ $`{\displaystyle ^{\text{gal}}}d\chi _1d\chi _2w_{\text{cmb}1}w_{\text{gal}1}w_{\text{cmb}2}w_{\text{gal}2}`$
$`\times {\displaystyle }{\displaystyle \frac{\mathrm{d}k_1\mathrm{d}k_2}{(2\pi )^2}}{\displaystyle \frac{\mathrm{d}^2l_{\text{gal}}\mathrm{d}^2l_{\text{cmb}}}{(2\pi )^4}}l_{\text{gal}}^4l_{\text{cmb}}^4k_1k_2`$
$`\times C_E(l_{\text{gal}})C_E(l_{\text{cmb}})P(k_1)P(k_2)`$
$`\times W^2(|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|\mathrm{\Xi })\mathrm{cos}^22\left(\varphi _{l_{\text{cmb}}}\varphi _{l_{\text{gal}}}\right)`$
$`\times W(k_1\mathrm{D}\theta _{\text{gal}})W(k_2\mathrm{D}\theta _{\text{gal}})W^2(l_{\text{gal}}\theta )`$
$`\times W(k_1\mathrm{D}\theta )W(k_2\mathrm{D}\theta )W^2(l_{\text{cmb}}\theta ).`$
Note that for the evaluation of the $``$ part, using the same kind of method, we obtain the same equation as eq. (III.5) where $`l_{\text{gal}}^4l_{\text{cmb}}^4`$ is replaced by $`l_{\text{gal}}^2l_{\text{cmb}}^2k_1^2k_2^2/2`$.

We can get rid of the remaining $`W^2\left(|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|\mathrm{\Xi }\right)`$ with another approximation. The power spectrum $`C_E(l)`$ favors large values of $`l`$ whereas $`W^2(|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|\mathrm{\Xi })`$ will be non-zero for $`|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|1/\mathrm{\Xi }`$. Then for typical survey size of about one hundred square-degrees, $`|\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}|l_i`$ and we can assume $`\stackrel{}{l}_{\text{gal}}\stackrel{}{l}_{\text{cmb}}`$ and $`\stackrel{}{l}_{\text{gal}}+\stackrel{}{l}_{\text{cmb}}=\stackrel{}{ϵ}`$. In this limit, $`\mathrm{cos}^22\left(\varphi _{l_{\text{cmb}}}\varphi _{l_{\text{gal}}}\right)=1`$ and $`𝒜_1^{\mathrm{}}`$ can be written

$`𝒜_1^\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{l\mathrm{d}l}{(2\pi )^2}l^8C_E^2(l)W^4(l\theta )\frac{\mathrm{d}^2ϵ}{2\pi }l^8W_1^2(ϵ\mathrm{\Xi })}`$
$`\times \left[{\displaystyle ^{\text{gal}}}d\chi w_{\text{cmb}}w_{\text{gal}}{\displaystyle \frac{k\mathrm{d}k}{2\pi }P(k)W(k\mathrm{D}\theta )W(k\mathrm{D}\theta _{\text{gal}})}\right]^2`$
which is essentially the cosmic variance of $`\mathrm{\Delta }E^2`$, for the $`\mathrm{\Delta }`$ part and of $`(E)^2`$ for the $``$ one (where $`l^8`$ in eq. (III.5) is replaced by $`l^4k_1^2k_2^2/2`$). Finally we have,

$`{\displaystyle \frac{𝒜_1^\mathrm{\Delta }}{B_{(\theta )}b_{\mathrm{\Delta }(\theta ,\theta _{\text{gal}})}^2}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{\Sigma }}}{\displaystyle \frac{dll^9C_E^2(l)W_1^4(l\theta )}{\left(dll^5C_E(l)W_1^2(l\theta )\right)^2}}`$
$``$ $`\text{Cosmic variance of }\mathrm{\Delta }E^2`$
where $`\mathrm{\Sigma }=\pi \mathrm{\Xi }^2`$ in case of a disc shape survey. We show in Fig. 8 numerical results for a $`100\mathrm{deg}^2`$ survey although the numerical calculations were done with a Gaussian window function instead of a top-hat.

Numerically, for $`\theta =10^{}`$, we get

$$\frac{𝒜_1^{\mathrm{}}}{B_{(\theta )}b_{\mathrm{}(\theta ,\theta _{\text{gal}})}^2}\frac{(3.7\%)^2}{\mathrm{\Sigma }/100\text{ deg}^2}.$$
(67)
We expect that for the same reasons, the $`𝒜_2^{\mathrm{}}`$ terms will be dominated by the weak lensing variance. Yet a correct evaluation here is harder to reach. We have made this estimation within the framework of a power law $`P(k)`$. With this simplification in hand, we can write for $`𝒜_{2n}^{\mathrm{}}`$ (we focus only the $`\mathrm{\Delta }`$ part, but the same discussion holds for the $``$ observable.)

$`{\displaystyle \frac{𝒜_{2n}^\mathrm{\Delta }}{B_{(\theta )}b_{\mathrm{\Delta }(\theta ,\theta _{\text{gal}})}^2}}`$ $`=`$
$`{\displaystyle \frac{1}{r^2}}{\displaystyle \mathrm{d}^2k_1\mathrm{d}^2k_2P(k_1)P(k_2)\mathrm{cos}^2\left(\varphi _{k_1}\varphi _{k_2}\right)}`$
$`\times {\displaystyle \frac{W_1^2(k_1\theta )W_1^2(k_2\theta _{\text{gal}})W_1^2(|\stackrel{}{k}_1+\stackrel{}{k}_2|\mathrm{\Xi })}{\left[\mathrm{d}^2kP(k)W_1(k\theta )W_1(k\theta _{\text{gal}})\right]^2}}.`$
The last integral behaves essentially like the cosmic variance of $`\kappa ^2`$. More exactly, it goes like $`1/\sqrt{2}`$ this variance. It should even be smaller, because of the extra $`\mathrm{cos}^2`$ factor. We evaluated this cosmic variance using the ray-tracing simulations described in . These simulations provide us with realistic convergence maps (for the cosmological models we are interested in) with a resolution of 0.1’, and a survey size of 9 square degrees. The sources have been put at a redshift unity, and the ray-lights are propagated through a simulated Universe whose the density field has been evolved from an initial CDM power spectrum. The measured cosmic variance of $`\kappa _{(\theta )}\kappa _{(\theta _{\text{gal}})}`$ is about $`3\%`$ (see Table 2) when filtered at scales $`\theta _{\text{gal}}=5^{}`$ and $`\theta =10^{}`$ for a $`\mathrm{\Omega }_0=0.3`$ cosmology.

An estimation of $`𝒜_{2n}^\mathrm{\Delta }`$ is then given by,

$$\frac{𝒜_{2n}^\mathrm{\Delta }}{B_{(\theta )}b_{\mathrm{\Delta }(\theta ,\theta _{\text{gal}})}^2}\left(\frac{2.12\%}{r}\right)^2\frac{1}{\mathrm{\Sigma }/100\text{ deg}^2}.$$
(70)
Since $`r_{}`$ is very comparable to $`r`$, we very roughly estimate $`𝒜_{2n}^{}`$

$$\frac{𝒜_{2n}^{}}{B_{(\theta )}b_{(\theta ,\theta _{\text{gal}})}^2}\left(\frac{2.12\%}{r}\right)^2\frac{1}{\mathrm{\Sigma }/100\text{ deg}^2}.$$
(71)
The same considerations gives

$`{\displaystyle \frac{𝒜_{2n}^{\mathrm{}}}{B_{(\theta )}b_{\mathrm{}(\theta ,\theta _{\text{gal}})}^2}}`$ $`=`$ $`{\displaystyle \frac{(2.12\%)^2}{\mathrm{\Sigma }/100\text{ deg}^2}}.`$ (72)
There is no $`r`$ dependency here; the diagram cross-correlates $`\kappa _{\text{cmb}}`$ and $`\kappa _{\text{gal}}`$.

We can approximate the remaining $`𝒜`$-terms. They should be smaller than the former. We have

$`𝒜_{3n}^{\mathrm{}}`$ $``$ $`{\displaystyle \frac{1}{r_{\mathrm{}}^2}}{\displaystyle \frac{(2.12\%\times 3.7\%)^2}{\mathrm{\Sigma }/100\text{ deg}^2}}B_{(\theta )}b_{\mathrm{}(\theta ,\theta _{\text{gal}})}^2`$
$``$ $`𝒜_{2n}^{\mathrm{}}`$
and

$`𝒜_{3c}^{\mathrm{}}`$ $``$ $`{\displaystyle \frac{(2.12\%\times 3.7\%)^2}{\mathrm{\Sigma }/100\text{ deg}^2}}B_{(\theta )}b_{\mathrm{}(\theta ,\theta _{\text{gal}})}^2`$
$``$ $`𝒜_{2c}^{\mathrm{}}.`$
Then, only the $`𝒜_1^{\mathrm{}}`$ and $`𝒜_2^{\mathrm{}}`$ terms (boxed on Fig. 7) contribute substantially to the cosmic variance of $`𝒳_{\mathrm{}}`$. Since $`𝒜_1^{\mathrm{}}`$ and $`𝒜_2^{\mathrm{}}`$ are respectively the cosmic variance of $`\mathrm{\Delta }E^2`$ (resp. $`(\stackrel{}{}E)^2)`$ and of $`\kappa ^2`$ (resp. $`(\stackrel{}{}\kappa )^2`$), we can write the variance of $`𝒳_{\mathrm{}}`$ as

$`\mathrm{C}osVar(𝒳_\mathrm{\Delta })`$ $`=`$
$`\mathrm{C}osVar\left(\mathrm{\Delta }E^2\right)+\left({\displaystyle \frac{1+r^2}{2r^2}}\right)\mathrm{C}osVar\left(\kappa ^2\right).`$
and

$`\mathrm{C}osVar(𝒳_{})`$ $`=`$
$`\mathrm{C}osVar\left((\stackrel{}{}E)^2\right)+\left({\displaystyle \frac{1+r_{}^2}{2r_{}^2}}\right)\mathrm{C}osVar\left((\stackrel{}{}\kappa )^2\right).`$
Table 3 presents numerical results for various filtering scenarii and models.

The two quantities, $`b_\mathrm{\Delta }`$ and $`b_{}`$, lead to similar cosmic variance that are rather small. Obviously it would be even better to use $`b=b_\mathrm{\Delta }+b_{}`$. For such a quantity the resulting cosmic variance for the cross-correlation coefficient should even be smaller, by a factor $`\sqrt{2}`$, although a detailed analysis is made complicated because of the complex correlation patterns it contains.

## IV Conclusion

We have computed a first order mapping that describes, in real space, the weak lensing effects on the CMB polarization. In particular we derived the explicit mathematical relation between the primary CMB polarization and the shear field at leading order in lens effect. It demonstrates that a $`B`$-component of the polarization field can be induced by lens couplings. We have shown however that the $`B`$-map alone cannot lead to a non-ambiguous reconstruction of the projected mass map.

Nonetheless, the $`B`$-component can potentially exhibit a significant correlation signal with local weak lensing surveys. This opens a new window for detecting lens effects on CMB maps. In particular, and contrary to previous studies involving the temperature maps alone, we found that such a correlation can be measured with a rather high signal to noise ratio even in surveys of rather modest size and resolution. Anticipating data sets that should be available in the near future, ($`100\text{deg}^2`$ survey, with $`5^{}`$ resolution for galaxy survey and $`10^{}`$ Gaussian beam size for CMB polarization detection), we have obtained a cosmic variance around $`8\%`$. Needless is to say that this estimation does not take into account systematics and possible foreground contaminations. It shows anyway that Cosmic Microwave Background polarization contains a precious window for studying the large scale mass distribution and consequently putting new constraints on the cosmological parameters.

In this paper we have investigated specific quantities that would accessible to observations. They both would permit to put constraint on the cosmological constant. The simulated maps we presented here are only of illustrative interest. We plan to complement this study with extensive numerical experiments to validate our results (in particular on the cosmic variance), and explore the effect of realistic ingredients we did not include in our simple analytical framework, a shear non-gaussianity, lens-lens coupling and so forth.

###### Acknowledgements.

We thank B. Jain, U. Seljak and S. White for the use of their ray-tracing simulations. KB and FB thank CITA for hospitality and LvW is thankful to SPhT Saclay for hospitality. We are all grateful to the TERAPIX data center located at IAP for providing us computing facilities.


text:
# Instabilities of vortices in binary mixture of trapped Bose-Einstein condensates: Role of collective excitations with positive and negative energies

## I Introduction

Recent observations of the quantized vortices in a dilute ultracold gas of $`{}_{}{}^{87}Rb`$ atoms  are spectacular evidences of the superfluid properties of atomic gases below transition temperature for formation of Bose-Einstein condensate (BEC). First experimental results were obtained, using method suggested by Williams and Holland , in a binary mixture of the hyperfine states of $`{}_{}{}^{87}Rb`$  and more recent experiments  have demonstrated vortices and vortex lattices in an optically stirred single component condensate.

Gross-Pitaevskii (GP) equation  is a widely used approximation to describe dynamics of the quantized vortices in the superfluid component of Bose gases at temperatures below critical. The GP equation falls into a general class of Hamiltonian nonlinear systems. The theory of stability of stationary solutions (equilibria) of such systems is well established at the moment, see e.g. . However, a certain gap between formal mathematical knowledge and applications to physical examples still exists and BEC field is not an exception. One of the controversial examples in the BEC context is complex or imaginary eigenfrequencies in the spectrum of elementary excitations. Existence of these frequencies was first pointed out by Bogoliubov in his original work , for more recent references see, e.g. . Unjustified negligence by the nonselfadjoitness of the Bogoliubov equations often lead to the association between frequencies and energies in a way standard for quantum mechanics based on selfadjoint operators, thus admitting possibility of complex energies in a conservative system.

However, rigorously proved theorem allowing clear physical interpretation of complex frequencies is known from the general theory of Hamiltonian systems . It states that excitations with complex frequencies can appear only as a result of the resonance between two excitations with positive and negative energies . The roots of this theorem go back to 19th century  and it has been used for interpretation of instabilities with complex eigenvalues in plasma physics  and fluid dynamics . Motivated by recent experiments on observation of vortices in binary mixture of condensates  we will demonstrate that these vortices can have complex frequencies in their spectrum, thereby giving good practical ground for selfconsistent theoretical interpretation of this phenomenon in the BEC context.

Properties of vortices in the single component magnetically trapped ultracold gases have been subject to the intensive theoretical investigations, see e.g. \[9,17-32\], which significantly extended classical works  dealing with spatially unbounded case. The properties of unit vortices in the two component condensates have also been studied, in parallel and independently from this work, by Garcia-Ripoll and Pérez-Garcia . Richness of the dynamics of the two-condensate system and different approaches to the problem have led only to few overlaps which are outlined where appropriate.

In the next section we introduce coupled GP equations and briefly describe their general properties. Then, in section III, we derive Bogoliubov equations for excitations, clearly specifying differences between frequency and energy spectra of the excitations. We also explain scenario of appearance of the excitations with complex frequencies and show that they have zero energies. In section IV we verify validity of general results presented in Sec. III using perturbation theory in the limit of weak interaction and direct numerical study of Bogoliubov equations. In Sections IV, V we describe how long term dynamics of unit and higher order vortices can be interpreted using linear Bogoliubov theory.

## II Gross-Pitaevskii equations

Studies of superfluid mixtures using coupled GP equations have long history and attracted significant recent activities, see e.g.  and references therein. Following these works we assume that wave functions $`\psi _{1,2}`$ of two-species condensate inside an axial harmonic trap obey equations

$`i\mathrm{}_t\psi _1={\displaystyle \frac{\mathrm{}^2}{2m}}\stackrel{}{}^2\psi _1+{\displaystyle \frac{1}{2}}m\mathrm{\Omega }^2(r^2+\sigma ^2z^2)\psi _1`$ (1)
$`+(u_{11}|\psi _1|^2+u_{12}|\psi _2|^2)\psi _1,`$ (2)
$`i\mathrm{}_t\psi _2={\displaystyle \frac{\mathrm{}^2}{2m}}\stackrel{}{}^2\psi _2+{\displaystyle \frac{1}{2}}m\mathrm{\Omega }^2(r^2+\sigma ^2z^2)\psi _2`$ (3)
$`+(u_{22}|\psi _2|^2+u_{21}|\psi _1|^2)\psi _2,`$ (4)
where for simplicity we have neglected by possible differences of atomic masses $`m_{1,2}=m`$ and trap frequencies, $`\mathrm{\Omega }_{1,2}=\mathrm{\Omega }`$, $`\sigma _{1,2}=\sigma `$, $`\stackrel{}{}=\stackrel{}{i}_x_x+\stackrel{}{i}_y_y+\stackrel{}{i}_z_z`$, $`r^2=x^2+y^2`$. Coefficients $`u_{ij}=4\pi \mathrm{}^2a_{ij}/m`$ characterise intra- and inter-species interaction with corresponding two-body scattering lengths $`a_{11}a_{22}`$ and $`a_{12}=a_{21}`$.

At this point we introduce dimensionless time and space variables $`\stackrel{~}{t}=\mathrm{\Omega }t`$ and $`(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{z})=(x,y,z)/a_{ho}`$ and normalisation for the wave functions $`\stackrel{~}{\psi }_{1,2}=a_{ho}^{3/2}\psi _{1,2}/\sqrt{N_1}`$, where $`a_{ho}=\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$ is the harmonic oscillator strength and $`N_{1,2}=𝑑V|\psi _{1,2}|^2`$ are the numbers of particles. We will consider the quasi-2D model to simplify our numerical study. This approximation was previously used in several works, see e.g. , and it is applicable not only for pancake traps, $`\sigma 1`$, but also captures main qualitative features of spherical traps. To make further reduction of Eqs. (2) we redefine the wave functions once more:

$$\stackrel{~}{\psi }_{1,2}=\left[\frac{\sigma }{2\pi }\right]^{1/4}\mathrm{\Psi }_{1,2}(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{t})e^{\sigma \stackrel{~}{z}^2/4}e^{i\mu _{1,2}\stackrel{~}{t}i\sigma \stackrel{~}{t}/2},$$
(5)
where $`\mu _{1,2}`$ are the chemical potentials. Dropping tilde we find that equations for $`\mathrm{\Psi }_{1,2}`$ and $`\mathrm{\Psi }_{1,2}^{}`$ can be put into Hamiltonian form

$$i_t\stackrel{}{\psi }+\widehat{\eta }\frac{\delta H}{\delta \stackrel{}{\psi }^{}}=0,$$
(6)
$$\stackrel{}{\psi }=\left[\begin{array}{c}\mathrm{\Psi }_1\\ \mathrm{\Psi }_1^{}\\ \mathrm{\Psi }_2\\ \mathrm{\Psi }_2^{}\end{array}\right],\widehat{\eta }=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],\frac{\delta H}{\delta \stackrel{}{\psi }^{}}=\left[\begin{array}{c}\delta H/\delta \mathrm{\Psi }_1^{}\\ \delta H/\delta \mathrm{\Psi }_1\\ \delta H/\delta \mathrm{\Psi }_2^{}\\ \delta H/\delta \mathrm{\Psi }_2\end{array}\right],$$

$`H`$ $`={\displaystyle }dxdy(|\stackrel{}{}\mathrm{\Psi }_1|^2+|\stackrel{}{}\mathrm{\Psi }_2|^2`$ (9)
$`+(\widehat{V}\mu _1)|\mathrm{\Psi }_1|^2+(\widehat{V}\mu _2)|\mathrm{\Psi }_2|^2`$
$`+{\displaystyle \frac{g}{2}}[\beta _{11}|\mathrm{\Psi }_1|^4+\beta _{22}|\mathrm{\Psi }_2|^4+2\beta _{12}|\mathrm{\Psi }_1|^2|\mathrm{\Psi }_2|^2]).`$
where $`H`$ is the energy functional (or Hamiltonian), $`\stackrel{}{}=\stackrel{}{i}_x_x+\stackrel{}{i}_y_y`$, $`\widehat{V}=r^2/4`$ is the harmonic potential, $`g`$ is the interaction parameter $`g=8\sqrt{\pi \sigma }N_1a_{11}/a_{ho}`$, $`\beta _{12}=a_{12}/a_{11}`$, and $`\beta _{22}=a_{22}/a_{11}`$. $`\beta _{11}=1`$ and it is left in the equations for the sake of the symmetry.

Invariancies of $`H`$ with respect to the infinitesimal rotations and two parameter phase transformation

$$(\mathrm{\Psi }_1,\mathrm{\Psi }_2)(\mathrm{\Psi }_1e^{i\varphi _1},\mathrm{\Psi }_2e^{i\varphi _2}),$$
(10)
result in conservation of the total angular momentum and of the total number of particles in each component.

Radially symmetric stationary states of the condensate (equilibria) can be presented in the form

$$\mathrm{\Psi }_j=A_j(r)e^{iL_j\theta },j=1,2$$
(11)
where $`\theta `$ is the polar angle and $`A_j`$ are real functions. Using method suggested in  only states with vortex in one of the components can be created and therefore we will consider below only cases with $`L_2>0`$ and $`L_1=0`$. Functions $`A_j(r)`$ were found numerically using Newton method. Chemical potentials $`\mu _{1,2}`$ were found from the normalization conditions

$$2\pi r𝑑rA_1^2=1,2\pi r𝑑rA_2^2=\frac{N_2}{N_1}N.$$
(12)
## III Frequency and energy spectra of collective excitations

### A Bogoliubov equations and frequency spectrum

To study spectrum of BEC at equilibrium we linearise Eqs. (6) using substitutions

$$\mathrm{\Psi }_j=(A_j(r)+f_j(r,\theta ,t))e^{iL_j\theta },$$
(13)
where $`f_j`$ are small and complex. Assuming that excitation are periodic in $`\theta `$ with period $`2\pi `$ we expand $`f_j`$ into Fourier series:

$$f_j=\underset{l}{}\left(U_{lj}(r,t)e^{il\theta }+V_{lj}^{}(r,t)e^{il\theta }\right),$$
(14)
where $`l=0,\pm 1,\pm 2,\mathrm{}`$ Then

$$i_t\stackrel{}{W}_l+\widehat{\eta }\widehat{}_l\stackrel{}{W}_l=0$$
(15)
is the set of linear partial differential equations resulting from the substitution of (13), (14) into Eq. (6). Here $`\stackrel{}{W}_l=(U_{l1},V_{l1},U_{l2},V_{l2})^T`$,

$$\widehat{}_l=\left[\begin{array}{cccc}\widehat{}_{l,1}& g\beta _{11}A_1^2& g\beta _{12}A_1A_2& g\beta _{12}A_1A_2\\ g\beta _{11}A_1^2& \widehat{}_{l,1}& g\beta _{12}A_1A_2& g\beta _{12}A_1A_2\\ g\beta _{12}A_1A_2& g\beta _{12}A_1A_2& \widehat{}_{l,2}& g\beta _{22}A_2^2\\ g\beta _{12}A_1A_2& g\beta _{12}A_1A_2& g\beta _{22}A_2^2& \widehat{}_{l,2}\end{array}\right],$$

is a self-adjoint operator and

$`\widehat{}_{l,j}={\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}r{\displaystyle \frac{}{r}}+{\displaystyle \frac{1}{r^2}}(L_j+l)^2+\widehat{V}\mu _j`$ (16)
$`+2g\beta _{jj}A_j^2+g\beta _{12}A_j^{}^2;j=1,2;j^{}=2,1.`$ (17)
Frequency spectra of $`\widehat{\eta }\widehat{}_l`$ are discrete, providing $`\widehat{V}0`$, therefore phonons, strictly speaking, are absent in the trap geometry. In accord with standard terminology , all spatially bounded elementary excitations can be called collective excitations (or collective modes) of an equilibrium under consideration. Linearised equations for excitations in a Bose gas, similar to Eq. (15), were first derived by Bogoliubov and expansion (14) was first applied in the context of the vortex excitations by Pitaevskii . To find frequency spectrum and collective modes we need to solve set of the eigenvalue problems

$$\widehat{\eta }\widehat{}_l\stackrel{}{w}_{ln}=\omega _{ln}\stackrel{}{w}_{ln}.$$
(18)
$`\widehat{\eta }\widehat{}_l`$ are non-self-adjoint operators and therefore complex frequencies are not forbidden. If $`\stackrel{}{w}_{ln}`$ is an eigenvector of $`\widehat{\eta }\widehat{}_l`$ with eigenvalue $`\omega _{ln}`$ it is selfevident that $`\stackrel{}{w}_{ln}^{}`$ is also an eigenvector with eigenvalue $`\omega _{ln}^{}`$ and it can be shown that $`\widehat{\eta }\widehat{}_l`$ has eigenvectors $`\stackrel{}{w}_{ln}=\widehat{\tau }\stackrel{}{w}_{ln}`$ and $`\stackrel{}{w}_{ln}^{}`$ with eigenvalues $`\omega _{ln}`$ and $`\omega _{ln}^{}`$, respectively. Here

$$\widehat{\tau }=\left[\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right].$$

Thus spectrum of $`\widehat{\eta }\widehat{}_l`$ can be obtained by reflection of the spectrum of $`\widehat{\eta }\widehat{}_l`$ with respect to the line $`Re\omega =0`$ in the plane $`(Re\omega ,Im\omega )`$. In other words, it means that purely real or purely imaginary frequencies of the elementary excitation exist in pairs and complex frequencies exist in quartets and that

$$Tr\{\widehat{\eta }(\widehat{}_l+\widehat{}_l)\}=\underset{n}{}(\omega _{ln}+\omega _{ln})=0.$$
(19)
Any equilibrium state of the condensate is spectrally stable  if its spectrum is real. If there is at least one frequency with negative imaginary part then corresponding collective mode will grow in time destabilising the equilibrium, which is called spectrally unstable.

### B Biorthogonality

Eigenmodes of $`\widehat{\eta }\widehat{}_l`$ are biorthogonal to the modes of the adjoint eigenvalue problem , i.e.

$$\stackrel{}{w}_{ln},\stackrel{}{a}_{ln^{}}=0,$$
(20)
where $`nn^{}`$ and $`\stackrel{}{a}_{ln}`$ obey

$$(\widehat{\eta }\widehat{}_l)^{}\stackrel{}{a}_{ln}=\widehat{}_l\widehat{\eta }\stackrel{}{a}_{ln}=\omega _{ln}^{}\stackrel{}{a}_{ln},$$
(21)
and $`\stackrel{}{b},\stackrel{}{c}=2\pi _k_0^{\mathrm{}}r𝑑rb_k^{}c_k`$ for any $`\stackrel{}{b}`$, $`\stackrel{}{c}`$. Factor $`2\pi `$ is introduced to mimic integration over $`\theta `$. The key feature of our model, originating in its Hamiltonian structure, is that transformation linking $`\stackrel{}{w}_{ln}`$ and its adjoint $`\stackrel{}{a}_{ln}`$ can be found in explicit and simple form. If $`\stackrel{}{w}_{ln}`$ is an eigenmode of $`\widehat{\eta }\widehat{}_l`$ with frequency $`\omega _{ln}`$ then it can be checked that $`\widehat{\eta }\stackrel{}{w}_{ln}^{}`$ and $`\widehat{\eta }\stackrel{}{w}_{ln}`$ are eigenmodes of $`\widehat{}_l\widehat{\eta }`$ with eigenvalues, respectively, $`\omega _{ln}^{}`$ and $`\omega _{ln}`$. The mode adjoint to $`\stackrel{}{w}_{ln}`$ is $`\widehat{\eta }\stackrel{}{w}_{ln}^{}`$, therefore if $`Im\omega _{ln}0`$ then biorthogonality condition (20) implies

$$\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}=0.$$
(22)
If $`\omega _{ln}`$ is real, then $`\stackrel{}{w}_{ln}`$ is also real and $`\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}0`$. Normalisation constant can always be chosen in such a way that

$$|\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}|=2.$$
(23)
Convenience of making the left hand side of Eq. (23) equal to $`2`$ will become clear below, when Eqs. (19), (20) for energies of elementary excitations are derived. Eq. (23) makes it explicit that inner product $`\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}`$ can be either positive or negative. This point often remains silent if one derives conditions similar to (23) as part of the diagonalisation procedure of the second-quantised Hamiltonian disregarding eigenmodes with negative and zero values of $`\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}`$ . The fact that inner product of a real eigenmode with its adjoint can be negative, is different from standard quantum mechanics based on the selfadjoint operators. These difference can have series of consequences and one of them is that energy levels are not necessarily linked to the eigenfrequencies according to the standard rule $`\mathrm{}\omega `$, see below. Note, that origin of the nonselfadjoitness in our case is the nonlinearity of GP equations. If particles in the condensate do not interact, $`g=0`$, then $`\widehat{\eta }\widehat{}_l`$ become diagonal and self-adjoint.

### C Energy spectrum

As a prelude to calculation of energies of the elementary excitations it is instructive to introduce notion of the nonlinear stability, i.e. stability under the full nonlinear dynamics. According to Dirichlet’s theorem , nonlinear stability is ensured if an equilibrium state under consideration is either minimum or maximum of the functional $`H`$. Note here, that excitations change number of particles in the equilibrium state, therefore it was convenient to introduce energy functional $`H`$, which is actually the so-called modified energy , i.e. it is the energy functional for Eq. (1) modified by addition of the number of particles integrals, see terms proportional to $`\mu _j`$ in Eq. (9). Let us assume that $`\stackrel{}{w}_l(r)e^{il\theta }`$ is a small initial perturbation of an equilibrium state of the condensate, then

$$H=H_0+\frac{1}{2}\stackrel{}{w}_l,\widehat{}_l\stackrel{}{w}_l+\mathrm{},$$
(24)
where $`H`$ is the energy of the perturbed equilibrium and $`H_0`$ is calculated at the exact equilibrium. The equilibrium is nonlinearly stable if eigenvalue problems

$$\widehat{}_l\stackrel{}{\beta }_{lm}=\alpha _{lm}\stackrel{}{\beta }_{lm}$$
(25)
have all their eigenvalues either negative or positive, except zero eigenvalues generated by continuous symmetries. Spectral instability implies nonlinear one, and nonlinear stability implies spectral one, but not vise versa . In the nonrotating traps all the higher order states of the condensate including vortices are not the local extrema of the energy . Let us stress, however, that their nonlinear instability can not be guaranteed by this fact alone and requires separate consideration.

If $`\stackrel{}{w}_l`$ in Eq. (24) is an eigenmode $`\stackrel{}{w}_{ln}`$ of $`\widehat{\eta }\widehat{}_l`$ then

$$ϵ_{ln}=\frac{1}{2}\stackrel{}{w}_{ln},\widehat{}_l\stackrel{}{w}_{ln}=\frac{1}{2}\omega _{ln}\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}$$
(26)
measures energy of this mode. Assuming that $`\omega _{ln}`$ is real and using biorthogonality conditions (23) one gets

$$ϵ_{ln}=\omega _{ln}sign\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}.$$
(27)
If $`Im\omega _{ln}0`$ then Eq. (22) implies

$$ϵ_{ln}=0.$$
(28)
Energy in physical units is given by $`ϵ_{ln}\mathrm{}\mathrm{\Omega }`$.

It is readily demonstrated that

$$\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln}=\widehat{\tau }\stackrel{}{w}_{ln},\widehat{\eta }\widehat{\tau }\stackrel{}{w}_{ln}=\stackrel{}{w}_{ln},\widehat{\eta }\stackrel{}{w}_{ln},$$
(29)
therefore modes with frequencies $`\omega _{ln}`$ and $`\omega _{ln}=\omega _{ln}`$ have their energies equal in sign and value.

Eigenfunctions $`\stackrel{}{\beta }_{ml}`$ are orthogonal $`\stackrel{}{\beta }_{lm},\stackrel{}{\beta }_{lm^{}}=\delta _{mm^{}}`$ and form complete basis. Therefore $`\stackrel{}{w}_{ln}`$ in Eq. (26) can be expanded in terms of $`\stackrel{}{\beta }_{lm}`$, which gives $`ϵ_{ln}=_m\alpha _{lm}\stackrel{}{w}_{ln},\stackrel{}{\beta }_{lm}^2`$. Thus an equilibrium can have collective modes carrying energy with opposite signs only if it is not a local extremum of the energy functional.

### D Resonances and instabilities, crossings and avoided crossings

Instabilities with imaginary eigenvalues have been found to dominate dynamics of homogeneous BEC with attractive interaction  and mixture of two condensates with repulsive interaction . Higher order vortices and dark solitons in single component trapped condensates can be unstable with respect to modes with complex frequencies . In our problem instabilities with complex eigenvalues are the most important ones also. Therefore it is desirable to have criterion or theorem allowing interpretation and prediction of these instabilities. Such theorem is actually known from the general theory of the Hamiltonian dynamical systems  and, in our context, it can be reformulated as: (i) Sign of the energy of a collective excitation preserves as parameters vary as long as there is no frequency resonance with another excitation. (ii) Condensate at equilibrium can lose spectral stability as parameters vary only in two ways: either via frequency resonance of two elementary excitations with positive and negative energies or by resonance at zero frequency. Proof of these results  is based on the fact, that transition from spectral stability to instability can not violate energy conservation law: $`_tH=0`$.

From this theorem and preceding considerations one can conclude that complex eigenvalues in the spectrum of the vortices in Bose-Einstein condensate can appear only due to mutual annihilation of the collective excitations with positive and negative energies. The latter can coexist only for nonlinearly unstable equilibria. The theorem does not forbid for frequencies of two excitations with either the same or opposite signs of the energies simply cross each other without change of the spectral stability. In our example it typically happens when corresponding eigenmodes remain orthogonal at the exact resonance. If, however, eigenmodes of the colliding excitations start to compete for the same direction in the functional phase space and become degenerate at the resonance then crossing is not a generic scenario. If the signs of the energies of the excitations are opposite, then quartet of complex frequencies appears upon passing the resonance. Alternatively, if the signs are the same, the exact resonance can not be achieved and it becomes replaced by the so-called avoided crossing of the energy levels .

Figs. 1(a),(b) show numerically calculated frequency and energy spectra of the collective excitations with $`l=\pm 1`$ for spectrally stable (left panel) and spectrally unstable (right panel) unit vortices $`(L_1=0,L_1=1)`$. The negative energy excitation is clearly seen in the spectrally stable situation. At the transition threshold to spectral instability this excitation and another one, having energy with the same absolute value but opposite sign, annihilate each other. This transition is accompanied by appearance of the zero energy excitations. Examples of frequency and energy evolution under the parameter variation resulting in instabilities, crossings and avoided crossing are given in the section IV. Figs. 1(c) show spectra of $`\widehat{}_{\pm 1}`$ illustrating that there are no qualitative changes in these spectra after appearance of complex frequencies.

Possibility of observation of the collective modes with complex frequencies has caused some concerns and discussions . However, initial perturbations of the equilibrium state having nonzero projections on the adjoint mode $`\widehat{\eta }\stackrel{}{w}_{ln}^{}`$ with $`Im\omega _{ln}^{}>0`$ will lead to the ultimate growth of the corresponding mode $`\stackrel{}{w}_{ln}`$ because it has $`Im\omega _{ln}<0`$. The consequences of this growth is in any way diminished by the fact that the energy of this mode is zero. Wealth of references and examples of instabilities with complex eigenvalues existing in other physical contexts can be found in .

### E Goldstone and dipole modes

Infinitesimal variations of $`\varphi _j`$, see Eq. (10), generate two zero-energy eigenmodes (Goldstone modes) $`(A_1,A_1,0,0)^T`$ and $`(0,0,A_2,A_2)^T`$ belonging to the null-eigenspace of $`\widehat{\eta }\widehat{}_0`$. Harmonic trapping modifies spectrum in such a way that $`\widehat{\eta }\widehat{}_{\pm 1}`$ acquire couple of parameter independent eigenvalues $`\omega =\pm 1`$ with eigenfunctions

$$\stackrel{}{w}_{\pm 1d}=\left[\begin{array}{c}\\ \omega \left(\frac{dA_1}{dr}\frac{1}{r}L_1A_1\right)+\frac{1}{2}rA_1\\ \\ \omega \left(\frac{dA_1}{dr}\pm \frac{1}{r}L_1A_1\right)\frac{1}{2}rA_1\\ \\ \omega \left(\frac{dA_2}{dr}\frac{1}{r}L_2A_2\right)+\frac{1}{2}rA_2\\ \\ \omega \left(\frac{dA_2}{dr}\pm \frac{1}{r}L_2A_2\right)\frac{1}{2}rA_2\end{array}\right],$$
(30)
where $`\omega `$ can take values $`\pm 1`$ for both eigenmodes. $`\stackrel{}{w}_{\pm 1d}`$ are often called dipole modes and $`ϵ_{\pm 1d}=1`$. Eqs. (30) generalise expressions previously derived for single-species condensates . Existence of the dipole modes can also be associated with the Kohn’s theorem .

## IV Collective excitations of unit vortices: $`L_1=0,L_1=1`$

We begin our analysis considering weakly interacting condensate: $`g1`$, $`N1`$. In this limit potential energy due to harmonic trapping $`\widehat{V}`$ strongly dominates over the interaction energy, which allows to make explicit calculations of the frequency and energy spectra of the elementary excitations. Calculations of the excitation spectra of the vortices and dark solitons in the singly component weakly interacting condensates have been recently done by several groups of authors . However, these calculations lack analysis of the energy sign in a sense explained in the preceding section.

We substitute asymptotic expansions $`\mu _{1,2}=\mu _{1,2}^{(0)}+g\mu _{1,2}^{(1)}+O(g^2)`$, $`A_{1,2}=A_{1,2}^{(0)}+gA_{1,2}^{(1)}+O(g^2)`$ into the stationary ($`_t=0`$) version of Eqs. (6) and derive recurrent system of linear equations. In the zero approximation we have two uncoupled harmonic oscillator problems with eigenmodes

$$A_1^{(0)}=\frac{1}{\sqrt{2\pi }}e^{r^2/4},A_2^{(0)}=\sqrt{\frac{N}{\pi }}\frac{r}{2}e^{r^2/4}.$$
(31)
Using solvability condition of the first order problem we find asymptotic expressions for the chemical potentials: $`\mu _1=1+g(2\beta _{11}+n\beta _{12})/(8\pi )+O(g^2)`$, $`\mu _2=2+g(n\beta _{22}+\beta _{12})/(8\pi )+O(g^2)`$.

Then we expand operators $`\widehat{}_l`$, eigenmodes $`\stackrel{}{w}_l`$ and frequencies $`\omega _l`$ into the series $`\widehat{}_l=\widehat{}_l^{(0)}+g\widehat{}_l^{(1)}+O(g^2)`$, $`\stackrel{}{w}_l=\stackrel{}{w}_l^{(0)}+g\stackrel{}{w}_l^{(1)}+O(g^2)`$, $`\omega _l=\omega _l^{(0)}+g\omega _l^{(1)}+O(g^2)`$. After substitution into Eq. (18) in the first order we find standard equation

$$(\widehat{\eta }\widehat{}_l^{(0)}\omega _l^{(0)})\stackrel{}{w}_l^{(1)}=(\omega _l^{(1)}\widehat{\eta }\widehat{}_l^{(1)})\stackrel{}{w}_l^{(0)}.$$
(32)
$`\widehat{\eta }\widehat{}_l^{(0)}`$ is diagonal and self-adjoint and all its eigenmodes and eigenvalues can be found explicitly. Then using solvability condition for Eq. (32) one can find corrections for all frequencies and coefficients in the linear superposition of the zero approximation eigenmodes. Computer algebra makes technical realization of this plan by a straightforward exercise. We will present and analyse explicit analytical results only for the operator $`\widehat{\eta }\widehat{}_1`$ considering vicinity of the spectral point $`\omega =1`$, because it contains information about origin of spectral instabilities of unit vortices. Equivalent analysis of $`\widehat{\eta }\widehat{}_1`$ near $`\omega =1`$ gives the same results.

$`\widehat{\eta }\widehat{}_1^{(0)}`$ has three unit frequencies, $`\omega _1^{(0)}=1`$. Corresponding eigenmodes are

$`\stackrel{}{b}_1=(0,1,0,0)^T{\displaystyle \frac{r}{2\sqrt{\pi }}}e^{r^2/4},`$ (33)
$`\stackrel{}{b}_2=(0,0,1,0)^T{\displaystyle \frac{1}{\sqrt{2\pi }}}e^{r^2/4},`$ (34)
$`\stackrel{}{b}_3=(0,0,0,1)^T{\displaystyle \frac{r^2}{4\sqrt{\pi }}}e^{r^2/4}.`$ (35)
Solvability conditions for Eq. (32) lead to the characteristic determinant of the three by three matrix. We find that one of the three frequencies is associated with dipole mode $`\stackrel{}{w}_{1d}`$ and that the other two are

$`\omega _1^\pm =1+{\displaystyle \frac{g}{32\pi }}\left(3\beta _{12}N\beta _{22}\pm \sqrt{R}\right)+O(g^2),`$ (36)
$`R(3\beta _{12}+N\beta _{22})^28\beta _{12}^2(N+1).`$ (37)
Corresponding unnormalised eigenmodes are

$`\stackrel{}{w}_1^\pm =`$ $`\sqrt{2N}\left(N\beta _{22}\beta _{12}\pm \sqrt{R}\right)\stackrel{}{b}_1`$ (40)
$`+\sqrt{2}\left(3N\beta _{22}+[14N]\beta _{12}\sqrt{R}\right)\stackrel{}{b}_2`$
$`+4N\left(\beta _{12}\beta _{22}\right)\stackrel{}{b}_3+O(g).`$
Pair of frequencies (36) becomes complex if $`R<0`$, signaling spectral instability of the unit vortex. If $`R>0`$, then

$`\stackrel{}{w}_1^\pm ,\widehat{\eta }\stackrel{}{w}_1^\pm =`$ (41)
$`4\left[R(N1)\pm \sqrt{R}\left\{(15N)\beta _{12}+N(3+N)\beta _{22}\right\}\right]`$ (42)
and one can check that

$$ϵ_1^\pm =\pm \omega _1^\pm ,$$
(43)
i.e. modes with frequencies $`\omega _1^+`$ and $`\omega _1^{}`$ have, respectively, positive and negative energies, see Fig. 2(c). One can also verify biorthogonality $`\stackrel{}{w}_1^\pm ,\widehat{\eta }\stackrel{}{w}_1^\pm =0=ϵ_1^\pm `$ for $`R<0`$, see Eq. (22), and $`\stackrel{}{w}_1^\pm ,\widehat{\eta }\stackrel{}{w}_1^{}=0`$ for $`R>0`$, see Eq. (23).

Assuming $`\beta _{ij}>0`$ and rewriting instability condition $`R<0`$ in the form

$$\frac{\beta _{22}}{\beta _{12}}<\frac{1}{N}(2\sqrt{2(N+1)}3),$$
(44)
one can see that unit vortices are more stable if intraspecies interaction of atoms in the vortex containing part of the condensate is somewhat larger than interspecies interaction. The choices of scattering lengths corresponding to the experiment  are: $`\beta _{22}=1.03/0.97`$, $`\beta _{12}=1/0.97`$ ($`\beta _{22}>\beta _{11}`$) and $`\beta _{22}=0.97/1.03`$, $`\beta _{12}=1/1.03`$ ($`\beta _{22}<\beta _{11}`$). The former case corresponds to the vortex in the spin state $`\{F=1,m_f=1\}`$ of $`{}_{}{}^{87}Rb`$ and the latter to the vortex in the state $`\{2,2\}`$. These two states will be called – state $`1`$ and state $`2`$. It is clear that for $`N=1`$ Eq. (44) predicts instability for the vortex in the state $`2`$ and stability for the vortex in the state $`1`$, which supports results of the experimental observations .

Fig. 2 shows frequency resonance accompanied by the simultaneous mutual annihilation of excitations with positive and negative energies happening at some critical value of $`\beta _{12}`$, see Eq. (44). In fact it models transition from the situation with vortex in the state 1 to the case with vortex in the state 2. Performing numerical studies for wide range of parameters, outside the validity region of analytical considerations, we have not been able to find regions of spectral instability of the unit vortex in state $`1`$. Contrary, existence of instabilities of the vortex in state $`2`$ due to exactly the same scenario, which is predicted in the weak interaction limit, can be readily demonstrated, see Fig. 3.

It is important to stress, that if two condensates are decoupled, $`\beta _{12}=0`$, then the negative energy mode $`\stackrel{}{w}_1^{}`$ belongs to the vortex containing component and the positive energy mode $`\stackrel{}{w}_1^+`$ belongs to the vortex free component. Because of this separation instability is not possible for any values of $`g`$, which agrees with spectral stability of unit vortices in the singly component case reported in . Thus we can conclude that the vortex instability in our example has essentially two-component nature and its analog may also exist in the case when second component is a non-condensate one.

Criterion similar to (44), but without corresponding energy analysis, has also been independently obtained in  using two mode approach, i.e. condensate wave functions $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ have been presented as linear superposition of modes (31) with time dependent coefficients and GP equations have been reduced to the set of ordinary differential equations for these coefficients . However, this method fails taking into account an eigenmode proportional to $`r^2`$, see $`\stackrel{}{b}_3`$ in (34), which makes an important contribution to the expressions for frequencies. Therefore Eqs. (36) and (44) are different from the corresponding results presented in .

To illustrate crossings and avoided crossings in the spectrum of unit vortices we show in Fig. 4 positive part of the frequency spectra of $`\widehat{\eta }\widehat{}_{\pm 1}`$ and $`\widehat{\eta }\widehat{}_{\pm 2}`$. One can see numerous points, which on the first glance can be interpreted as crossings in the frequency spectrum. However, under the close investigation those of them which are marked by the open circles, turns to be the avoided crossing of the excitations with equal energy signs, see Fig. 5.

## V Drift of unit vortices

The unstable modes of the vortex in the state $`2`$ are of the dipole type, i.e. $`|l|=1`$, therefore their growth leads to the initial displacement of the vortex from the trap center. Displaced vortex carries less of the total angular momentum, compare to the momentum of the vortex positioned at the trap center. Lack of the angular momentum is compensated by two factors. First, vortex acquires a nonzero tangential velocity and therefore its trajectory is actually a spiral, which is similar to the dynamics of an optical vortex displaced from the center of a gaussian beam . Second, angular momentum, and vortex itself, become gradually transferred into the second condensate component, which was first demonstrated in .

If drift instability is absent, then dissipative effects still can result in the vortex drift, which was predicted for the single component condensates by several authors . Presence of the second condensate opens a channel for the energy and angular momentum transfer from the vortex into the vortex free component. It can be seen that vector $`\stackrel{}{b}_1`$, being excited by the instability, provides a channel for this transfer. Thus drift instability can also be interpreted as due to dissipation of the energy and momentum by the vortex free condensate component.

In the limits $`N1`$ and $`N1`$ our model can be approximately considered as a single component condensate with ($`N1`$) or without ($`N1`$) vortex. Unit vortex and ground state of the single component condensate are known to be spectrally stable. Therefore drift instability disappears in both limits, see Fig. 3. Increase of $`g`$ for fixed $`N`$ also results in the suppression of the instability, see Fig. 3, which means that not only relative, but also absolute increase of the number of atoms in the vortex free component stabilizes the condensate.

## VI Drift and splitting of higher order vortices

Considering higher order vortices one can expect to find instability scenario resulting in their splitting into unit vortices. This scenarion is expected to be due to growth of the collective modes with $`|l|>1`$. However, as we will see below the drift instability linked to the dipole like modes, $`|l|=1`$, also can be presented. It leads to displacement of the whole vortex from the trap center without splitting, at least at the onset of the instability.

Both drift and splitting instabilities of the higher-order vortices appear as a result of the frequency resonances of the elementary excitations with negative and positive energies, similar to the case of unit vortices. $`N1`$ corresponds to the vortex free condensate and therefore both instabilities disappear in this limit. In the limit $`N1`$ only drift instability is suppressed and one can recover periodic in $`N`$ bands of the instabilities with complex frequencies and $`|l|>1`$ similar to the results reported for higher order vortices in a singly component condensate . Note, that splitting itself was not explicitly demonstrated in . It is also interesting to note that higher-order vortices in the free, $`\widehat{V}=0`$, single component condensate are spectrally stable . Thus splitting can be considered as induced by the trapping. The vortex free condensate component plays crucial role in the drift instability, but will the latter one be presented without trapping or not remains an open problem.

### Double vortices: $`L_1=0,L_2=2`$

Considering double vortex we have found that it can be unstable with respect to the $`|l|=1,2`$ excitations. The vortex in the state $`1`$ has been found surprisingly stable. One has to take relatively small values of $`g`$ and large $`N`$ to find splitting instability. Vortex in the state $`2`$ is more unstable in a sense that splitting exists already for $`N1`$, see Fig. 6. As it is evident from Fig. 6 either drift or splitting instability can dominate vortex dynamics. If drift instability is dominant, then vortex first gets displaced from the trap center and only then splits into the unit ones, see Fig. 7. The latter happens due to the curved background which breaks cylindrical symmetry with respect to the vortex axis. After the splitting vortices remain close to each other and move towards the condensate periphery. Results of the numerical simulation of GP equations (1) presented in this Section were obtained starting from equilibrium states perturbed by random noise with amplitude $`0.05A_{1,2}`$.

The dynamics is quite different when splitting instability is dominant, see Fig. 8. In this case unit vortices appear straight at the onset of the instability development and spiral out of the condensate center being always positioned symmetrically with respect to it. After a certain period of time vortices move back to the trap center and condensate state close to the initial one is restored, see Fig. 8, then the cycle is repeated with gradually worsening degree of periodicity.

During development of the instability angular momentum and vorticity become partially transferred into the second condensate. Analysis of the transverse profiles of the phases corresponding to the density profiles shown in Fig. 8, has revealed that black spots appearing in the second condensate are indeed unit vortices, not the density holes without topological structure.

### Triple vortices: $`L_1=0,L_2=3`$

Reach variety of beautiful vortex lattices can be found considering instabilities of vortices of the order 3 and higher. This reachness can be understood in terms of spatial profiles of the unstable collective modes. E.g. triple vortex has been found unstable with respect to the perturbations with $`|l|=1,2,3`$. All components of the $`|l|=2`$ excitations are equal to zero at the trap center. Therefore one can expect that growth of this mode will develop into a spatial structure preserving vortex at the trap center. Contrary some of the components of the $`|l|=3`$ modes have humps at the center. Therefore their growth should repel all unit vortices out of the center, which leads to the breaking of the triple vortex into a triangular structure of the unit vortices moving away from the trap center. Both instabilities can be found for the same parameters values and can have close growth rates. Therefore the winning mode is selected through the process of complex competition, see Figs. 9,10.

## VII Summary

We have described general approach to stability of equilibrium in BEC using Bogoliubov theory and GP equation. Biorthogonality conditions (22), (23) and correspondence between frequency and energy spectra of the elementary excitations (27), (28) have been derived selfconsistently from first principals revealing several novel and conceptually important aspects originating in nonseldjoitness of the Bogoliubov operator.

It has been demonstrated that frequency resonances of the excitations with positive and negative energies can lead to their mutual annihilation and appearance of the collective modes with complex frequencies and zero energies. Conditions for the avoided crossing of energy levels have also been discussed. General theory has been verified both numerically and analytically in the weak interaction limit considering an example of vortices in a binary mixture of condensates.

Growth of excitations with complex frequencies leads to the two main scenarios of the instability development. First one is the spiraling of unit and double vortices out of the condensate center to its periphery. Second scenario is the splitting of the double and higher order vortices into unit ones. Absolute and/or relative increase of the number of particles in the vortex free condensate component have been found to have stabilizing effect.

###### Acknowledgements.

Author is particularly grateful to referees for their critical and very helpful comments. He also acknowledges discussions with S.M. Barnett and W.J. Firth. Numerical part of the work was significantly speeded up due to access to the computer equipment obtained via U.K. EPSRC grant GR/M31880 and assistance of G. Harkness and R. Martin.