text:
# Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases

## Abstract

We show that one-dimensional binary mixtures of bosons or of a boson and a spin-polarized fermion are Luttinger liquids with the following instabilities: i) For different particle densities, strong attraction between the mixture components leads to collapse, while strong repulsion leads to demixing, ii) For a low-density mixture of two gases of impenetrable bosons (or a spin-polarized fermion and an impenetrable boson) of equal densities, the system develops a gap and exhibits enhanced pairing fluctuations when there is attraction between the components. In the boson-fermion mixture, the pairing fluctuations occur at finite momentum. Our conclusions apply to mixtures both on the continuum and on optical lattices away from integer or fractional commensurability.

The exploration of properties of matter at increasingly lower temperatures and densities has yielded many surprises. The achievement of Bose-Einstein condensation (BEC) in dilute gases of ultracold atoms Anglin02  is one beautiful example. Ever since the observation of BEC, one line of current experimental (and theoretical) studies is on cold gases in very elongated atomic traps. In effectively one-dimensional (1D) systems strong phase fluctuations can occur. These systems are interesting for both large and small $`\gamma `$ values, where $`\gamma =Mg/\mathrm{}^2\rho _0`$ is the dimensionless coupling ratio, $`g`$ is the coupling characterizing the interaction between the atoms in 1D, $`M`$ the atom mass, and $`\rho _0`$ the particle density. For $`\gamma 1`$, a single-component 1D Bose fluid is a quasicondensate with a fluctuating phase but strongly suppressed density fluctuations Petrov00Andersen02 . Large $`\gamma `$ leads to a crossover to the Tonks-Giradeau (TG) regime Petrov00Andersen02 : the bosons become impenetrable, and resemble noninteracting fermions Girardeau60 . The low-energy properties in both regimes can be well described by the Luttinger liquid conceptHaldane ; Cazalilla02 .

These aforementioned regimes make up the state space of the single-component system of 1D repulsive bosons. It may therefore seem surprising that adding a second component to form a mixture of two types of bosons ($`B+B`$) or of a boson plus a spin-polarized fermion ($`B+F`$) can lead to a richer phase diagram. In this Letter, we show that such binary mixtures are Luttinger liquids which can become unstable to gap-opening in one branch of the excitation spectrum, or to demixing, or to collapse.

These instabilities occur at low temperatures by tuning some parameters of the system, signalling a quantum phase transition (QPT) Sachdev99 . Some of these instabilities have direct analogues in higher dimensional systems, while others may be more specific to 1D systems (and are analogues of 1D spin-$`1/2`$ fermionic systems LutherEmery74 ; Gogolin98 ). Dilute cold gases are ideal for observing the instabilities that we predict because, in contrast to conventional solid-state materials, they are cleaner and there is a large degree of control and tunability of parameters such as particle densities and the sign and strength of the interaction. Furthermore, diluteness makes it a priori possible to relate microscopic parameters such as scattering lengths and atomic masses to the parameters of the low-energy effective theory.

Contrary to their 3D counterparts MathewsHall98 ; TruscottSchreck01 ; mix ; Timmermans98 , 1D binary mixtures have only recently attracted interests EckernParedes03 ; Das03 . To our knowledge, there has been no analytic treatment of these systems within the Luttinger liquid framework. This is undertaken here using the harmonic fluid approach, which takes into account phase and density fluctuations in 1D fluids Haldane ; Cazalilla02 . We summarize our results:

Case 1.- We consider a binary mixture of two distinct atomic species (either $`B+B`$ or $`B+F`$) or of two internal states of the same bosonic species. When a sufficiently repulsive interaction exists between the components, the mixture will be unstable against demixing, whereas for sufficient attraction, it will collapse.

Case 2.- We consider a binary mixture in a low-density limit, when the two components are sufficiently close in density and if their sound velocities are similar, the system undergoes a QPT when there is an attraction, no matter how weak, between the two components. The resulting state exhibits a strong pairing tendency. In the case of a $`B+F`$ mixture, the pairing fluctuations occur at the Fermi momentum, unlike the $`B+B`$ case (or $`F+F`$ LutherEmery74 ) where they occur at zero momentum.

Model and calculations: The experimental realization of the quasi one-dimensional systems requires a tight transverse confinement Olshanii98 ; Petrov00Andersen02 , which nowadays can be achieved in various setups traps ; Hansel01 ; Sauer01 . Our discussion below will focus on the properties of homogeneous 1D systems. This should be appropriate for the central region of large harmonic traps  or square-well traps Hansel01 . Alternatively, a toroidal trap Sauer01  can provide the conditions for realizing a homogeneous 1D system. When necessary, we shall also discuss the effect of longitudinal confinement and finite size on the results to be presented below.

Our starting point will be the following Hamiltonian for a mixture of two 1D dilute gases ($`\alpha ,\beta =1,2`$):

$`H`$ $`=`$ $`{\displaystyle 𝑑x\underset{\alpha }{}\left[\frac{\mathrm{}^2}{2M_\alpha }_x\mathrm{\Psi }_\alpha ^{}(x)_x\mathrm{\Psi }_\alpha (x)\mu _\alpha \rho _\alpha (x)\right]}`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x\underset{\alpha ,\beta }{}g_{\alpha \beta }\rho _\alpha (x)\rho _\beta (x)}.`$
For bosons, the field operators, $`\mathrm{\Psi }_\alpha (x)`$, obey $`[\mathrm{\Psi }_\alpha (x),\mathrm{\Psi }_\alpha ^{}(x^{})]=\delta (xx^{})`$, commuting otherwise. For a spin-polarized fermion, the field operator obeys anti-commuting relations: $`\{\mathrm{\Psi }_2(x),\mathrm{\Psi }_2^{}(x^{})\}=\delta (xx^{})`$. The density operators $`\rho _\alpha (x)=\mathrm{\Psi }_\alpha ^{}(x)\mathrm{\Psi }_\alpha (x)`$. The inter-atomic interaction potential is $`v_{\alpha \beta }(x)=g_{\alpha \beta }\delta (x)`$, with $`g_{\alpha \beta }=2\mathrm{}\omega _{}a_{\alpha \beta }`$ Olshanii98 , where $`a_{\alpha \beta }`$ are the scattering lengths parametrizing the 3D interaction between species $`\alpha `$ and $`\beta `$, and $`\omega _{}`$ is the transverse confinement frequency.

To study the low-energy properties of the above model, we use the harmonic fluid approach Haldane ; Cazalilla02 , which employs the phase-density representation of the field operators. Phonon fields are introduced for each component: $`\varphi _\alpha (x)`$ for long wavelength phase fluctuations and $`\theta _\alpha (x)`$ for long wavelength density fluctuations. It is assumed that $`_x\theta _\alpha (x)`$ and $`_x\varphi _\alpha (x)`$ are small compared to the equilibrium densities $`\rho _{0\alpha }`$. This approach treats bosons and fermions in 1D on equal footing Haldane ; Cazalilla02 , and does not assume the existence of a BEC, as usual mean-field treatments do. Therefore, it can describe both systems of impenetrable bosons (e.g. the Tonks gas Girardeau60  and others Cazalilla03 ), and quasicondensates Petrov00Andersen02 . Introducing Haldane ; Cazalilla02  $`\mathrm{\Psi }_\alpha (x)\left[\rho _{0\alpha }+_x\theta _\alpha (x)/\pi \right]^{1/2}_me^{im\pi \rho _{0\alpha }x+im\theta _\alpha (x)}e^{i\varphi _\alpha (x)}`$, where the sum over $`m`$ involves only even (odd) integers if the operator is bosonic (fermionic), and $`\rho _\alpha =\left[\rho _{0\alpha }+_x\theta _\alpha (x)/\pi \right]_{m=\mathrm{}}^+\mathrm{}e^{2im\pi \rho _{\mathrm{o}\alpha }x+2im\theta _\alpha (x)}`$, into Eq. (1) one obtains the low-energy Hamiltonian:

$`H_{\mathrm{eff}}={\displaystyle \frac{\mathrm{}}{2\pi }}{\displaystyle 𝑑x\underset{\alpha }{}\left[\frac{v_\alpha }{K_\alpha }\left(_x\theta _\alpha \right)^2+v_\alpha K_\alpha \left(_x\varphi _\alpha \right)^2\right]}`$
$`+{\displaystyle \frac{\mathrm{}}{2\pi }}{\displaystyle 𝑑x\left[2\stackrel{~}{g}_f_x\theta _1_x\theta _2+\stackrel{~}{g}_b\mathrm{cos}2(\theta _1\theta _2+\pi \delta x)\right]}.`$ (2)
In Eq.(Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases), we have retained only terms which can have the most dominant effects at low energies (i.e. marginal or relevant in the renormalization group sense). Its validity is restricted to energies smaller than the chemical potential. The low-energy physics is fully characterized by the phenomenological parameters $`K_{1,2}`$, the sound velocities $`v_{1,2}`$, the couplings $`\stackrel{~}{g}_{\mathrm{f},\mathrm{b}}`$, and $`\delta =\rho _{01}\rho _{02}`$. For small enough $`g_{12}`$, it is possible to relate these phenomenological parameters to the microscopic ones because exchange and correlation effects will be of $`\mathrm{O}(g_{12}^2)`$ in $`\stackrel{~}{g}_{\mathrm{f},\mathrm{b}}`$ (Otherwise, they must be extracted from numerics or experiments Cazalilla02 ). First notice that $`v_{N\alpha }=1/(\mathrm{}\pi \rho _\mathrm{o}^2\kappa _{S\alpha })`$, where $`\kappa _{S\alpha }=\rho _\alpha ^2(\rho _\alpha /\mu _\alpha )`$ is the compressibility Haldane ; Cazalilla02 . Furthermore, in homogeneous systems, Galilean invariance fixes the product $`v_{J\alpha }=v_\alpha K_\alpha =v_{F\alpha }=\mathrm{}\pi \rho _{0\alpha }/M_\alpha `$ Haldane . Thus from the exact solution for a single-component 1D Bose fluid Lieb63 , $`v_{N\alpha }=v_\alpha K_\alpha ^1=v_{F\alpha }\left[18\gamma _\alpha ^1+O(\gamma _\alpha ^2)\right]`$ for $`\gamma _\alpha =M_\alpha g_{\alpha \alpha }/\mathrm{}^2\rho _{0\alpha }1`$ and $`v_{N\alpha }v_{F\alpha }\gamma _\alpha /\pi ^2`$, for $`\gamma _\alpha 1`$. Therefore, $`1K_\alpha <+\mathrm{}`$ for bosons, $`K_\alpha =1`$ in the Tonks limit. For spin-polarized fermions, $`K_21`$, because the s-wave scattering length vanishes thanks to the Pauli principle, leaving the much weaker p-wave channel, which can be neglected. Finally, in the weak coupling limit, $`\stackrel{~}{g}_fg_{12}/(\mathrm{}\pi )`$ and $`\stackrel{~}{g}_b4\pi g_{12}\rho _{01}\rho _{02}/\mathrm{}`$. The sign and strength of $`g_{12}`$ can be controlled using a Feshbach or a confinement induced Olshanii98  resonance. Interestingly, the same Hamitonian, Eq (Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases), also describes a binary mixture in a 1D optical lattice unpub  provided that none of the components is commensurate with the lattice periodicity. This allows for further possibilities of tuning the parameters $`v_\alpha `$, $`K_\alpha `$, and $`g_{12}`$ (e.g. the strength of $`g_{12}`$ can be modified by shifting relative to each other the lattices where each component hops).

The above Hamiltonian, Eq. (Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases), describes a rich variety of one-dimensional binary mixtures:

Case 1.- When the equilibrium densities of the two species are different, i.e. for sufficiently large $`\delta =\rho _{01}\rho _{02}`$, the cosine term in Eq. (Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases) can be neglected (see below). Thus one is left with two Luttinger liquids coupled by $`\stackrel{~}{g}_f_x\theta _1_x\theta _2`$. The normal modes of the system can be found from the equations of motion for $`\theta _1(x,t)`$ and $`\theta _2(x,t)`$. The phase velocity of the normal modes is:

$$v_\pm ^2=\frac{1}{2}(v_1^2+v_2^2)\pm \frac{1}{2}\sqrt{\left(v_1^2v_2^2\right)^2+4(v_{J1})(v_{J2})\stackrel{~}{g}_f^2}.$$
(3)
Hence, the mixture will become unstable provided that

$$\sqrt{v_{N1}v_{N2}}<|\stackrel{~}{g}_f|.$$
(4)
For a repulsive interaction ($`\stackrel{~}{g}_f>0`$), the two components of the mixture repel each other sufficiently strongly to demix. However, for an attractive interaction ($`\stackrel{~}{g}_f<0`$), the system is expected to collapse.

We now examine the linear instability condition, Eq. (4), in several limiting cases. For a $`B+B`$ mixture in the quasicondensate regime ($`\gamma _{1,2}1`$), $`v_{N\alpha }v_{F\alpha }\gamma _\alpha /\pi ^2`$ and $`\stackrel{~}{g}_fg_{12}/(\mathrm{}\pi )`$. Then the mixture is unstable when $`\sqrt{g_{11}g_{22}}<|g_{12}|`$. In higher dimensional systems this condition is obtained from a mean-field theory Timmermans98 , which implicitly assumes the existence of a condensate. It is therefore interesting that it also applies to 1D systems in the quasicondensate regime, where the phase fluctuates. Another interesting limit is a $`B+F`$ mixture where the boson is a quasi condensate. For the fermion, $`v_{N2}v_{F2}=\mathrm{}\pi \rho _{02}/M_2`$. Hence, the system will become unstable below a critical fermion density $`\rho _{02}^{\mathrm{crit}}=M_2g_{12}^2/(\pi ^2\mathrm{}^2g_{11})`$. This agrees with the mean-field result of Das Das03 . Lastly, consider a $`B+F`$ or a $`B+B`$ mixture with impenetrable boson(s) ($`\gamma _\alpha 1`$). Then Eq. (4) implies that the system becomes unstable when $`\left(\rho _{01}\rho _{02}\right)<M_1M_2g_{12}^2/(\pi ^4\mathrm{}^4)`$.

We now briefly discuss the properties of a stable mixture. As found above, at long wavelengths there will be two branches of modes with phase velocities $`v_\pm `$ given by Eq. (3). This result is not qualitatively affected by the presence of a longitudinal confinement, though the actual energies of the modes may. More interestingly, a stable binary mixture will be a Luttinger liquid with all correlations characterized (in the thermodynamic limit and at zero temperature) by the usual power laws Haldane . In a trapped system, the power laws will be accurate near the center of the trap Cazalilla02 . A full discussion of these issues will be given elsewhere unpub .

Case 2: When the components of the mixture are close enough in density, one can no longer neglect the cosine term in Eq. (Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases). To see this, consider the Tonks limit where the bosons behave as free fermions. Then this term describes a backscattering process between the fermions: A fermion of each species is excited from one Fermi point to the other, exchanging a momentum close to $`2\pi (\rho _0\pm \delta )`$. If $`\delta \rho _0=(\rho _{01}+\rho _{02})/2`$ then this process can occur at low energies and needs to be taken into account. To simplify our analysis, we take $`\delta =0`$ ($`\delta 0`$ will be addressed at the end). By a perturbative analysis of the free energy, we see that the correction due to this term becomes larger than the zeroth order term as $`T0`$ for $`K_1+K_2<2`$. To go further, we use the renormalization group (RG) method and take $`K_1=K_2=1`$ and $`v_1=v_2`$ (We have checked that small deviations from these assumptions do not alter the physical picture qualitatively unpub ). The rationale for this is that the lower bound in a boson system is $`K_\alpha =1`$, which is reached for the Tonks gas limit (i.e. by reducing the density, increasing the scattering lengths $`a_{\alpha \alpha }`$ or both Petrov00Andersen02 ). In a $`B+F`$ mixture, however, the fermion already has $`K_21`$, and we only need to tune the boson to the Tonks limit. Furthermore, if the velocites $`v_1`$ and $`v_2`$ can be made made similar 111On the continuum, this is the case for a $`B+B`$ mixture of two internal states of 87Rb MathewsHall98  or a $`B+F`$ mixture of two isotopes of Li TruscottSchreck01  because their masses are equal or similar. On an optical lattice, one can adjust the lattice parameters to have $`v_1v_2`$ even if $`M_1M_2`$., then in the Tonks limit at equal densities, $`v_1=v_2=v_F`$. Thus the Hamiltonian (Instabilities in Binary Mixtures of One-Dimensional Quantum Degenerate Gases) becomes $`H_{\mathrm{eff}}=H_++H_{}`$, where:

$`H_+`$ $`=`$ $`{\displaystyle \frac{\mathrm{}v_+}{2\pi }}{\displaystyle 𝑑x\left[K_+^1\left(_x\theta _+\right)^2+K_+\left(_x\varphi _+\right)^2\right]},`$ (5)
$`H_{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}v_F}{2\pi }}{\displaystyle 𝑑x\left[\left(_x\theta _{}\right)^2+\left(_x\varphi _{}\right)^2\right]}`$ (6)
$`+`$ $`{\displaystyle \frac{\mathrm{}}{2\pi }}{\displaystyle 𝑑x\left[\stackrel{~}{g}_b\mathrm{cos}\sqrt{8}\theta _{}(x)\stackrel{~}{g}_f\left(_x\theta _{}\right)^2\right]}.`$
The new fields $`\theta _\pm (x)=\left(\theta _1(x)\pm \theta _2(x)\right)/\sqrt{2}`$ and $`\varphi _\pm (x)=\left(\varphi _1(x)\pm \varphi _2(x)\right)/\sqrt{2}`$ describe the in-phase (“$`+`$”) and out-of-phase (“$``$”) density and phase fluctuations of the components of the mixture. The fact that in this limit they are decoupled is an analog of spin-charge separation, so ubiquitous in 1D Fermi systems Gogolin98 ; Recati03 .

$`H_+`$ has the form of a Luttinger liquid Hamiltonian characterized by the parameters $`K_+=\left[1+\stackrel{~}{g}_f/v_F\right]^{1/2}\left[1+g_{12}/(\mathrm{}\pi v_F)\right]^{1/2}`$ and $`v_+=v_F\left[1+\stackrel{~}{g}_f/v_F\right]^{1/2}v_F\left[1+g_{12}/(\mathrm{}\pi v_F)\right]^{1/2}`$. Hence, in-phase modes are gapless (in the thermodynamic limit).

$`H_{}`$ is a sine-Gordon model Gogolin98 . The relative importance of the marginal operators, the cosine term and $`(_x\theta _{})^2`$, can be assessed using the RG at weak coupling, where $`g_f=\stackrel{~}{g}_f/v_F1`$ and $`g_b=\rho _0^2\stackrel{~}{g}_b/4\pi ^2v_F1`$. To second order in these couplings, the scaling flow is of Berezinskii-Kosterlitz-Thouless (BKT) type (see e.g. Refs. Gogolin98 ; Sachdev99 ). It is interesting to note that if the sound velocities are tuned to be equal (i.e. $`v_1=v_2`$) in the limit $`K_1=K_2=1`$, then the RG flow proceeds along the separatrix $`g=g_f=g_bg_{12}`$ and the spectrum of the system exhibits an enhanced SU$`(2)`$ symmetry. This is not present in the microscopic Hamiltonian (1), and seems quite striking, especially for a $`B+F`$ mixture. Thus, for small $`g_{12}>0`$, the effective coupling $`g`$ is renormalized to zero as the temperature is decreased and the cosine term leads only to subleading corrections to the low-temperature properties: the out-of-phase modes will behave as a Luttinger liquid with $`K_{}=1`$ and $`v_{}=v_F`$. But when small $`g_{12}<0`$, the effective coupling $`g_b+\mathrm{}`$ with decreasing temperature. The exact solution of the sine-Gordon model Gogolin98  then tells us that a gap opens for the excitations of the out-of-phase mode. At weak $`gg_{12}`$ the gap $`\mathrm{\Delta }|g|^{1/2}e^{1/|g|}`$ Gogolin98 ; Sachdev99 , i.e. exponentially small, and it becomes larger as $`g_{12}`$ increases. Hence, at low enough temperatures only in-phase modes can be excited. In a finite-size system, however, the RG flow to strong coupling will be cut-off at a length scale comparable to the size of the system. Therefore, for weak attraction one should observe that the out-of-phase modes become more “stiff” until the full gap develops.

If $`v_1v_2`$ and/or $`K_1,K_21`$ (with the differences small) the effective SU(2) symmetry is lost and the transition is of BKT type, the gap again being exponentially small but with a different functional form (see e.g. Gogolin98 ; Sachdev99 ). One notable consequence of the development of the gap is that density fluctuations of both components will become correlated because only in-phase modes are allowed at low enough temperatures: specifically, the $`2\pi \rho _0`$ part of density correlation $`\rho _\alpha (x)\rho _\beta (0)|_{2\pi \rho _0}\mathrm{cos}(2\pi \rho _0)x^{K_+}`$ for arbitrary $`\alpha =\beta =1,2`$ and $`\alpha \beta `$. These correlations can be observed by exciting the modes of the system and imaging each species independently.

The aforementioned correlations are also reflected in the “pairing” correlations described by the operator $`\mathrm{\Delta }(x)=\mathrm{\Psi }_1(x)\mathrm{\Psi }_2(x)`$. For a $`B+B`$ mixture, $`\mathrm{\Delta }^{}(x)\mathrm{\Delta }(0)x^{1/K_+}`$, which decays more slowly than the density correlation, since for $`g_{12}<0`$, $`K_+\left[1+g_{12}/(\mathrm{}\pi v_F)\right]^{1/2}>1`$. For a $`B+F`$ mixture, the same correlation function has an interesting oscillation at $`\pi \rho _0`$ (rather than at $`2\pi \rho _0`$): $`\mathrm{\Delta }^{}(x)\mathrm{\Delta }(0)\mathrm{sin}\left(\pi \rho _0x\right)x^{1/K_+K_+/4}`$. Thus, in a $`B+F`$ mixture, the pairing correlations will decay more slowly than the density correlations for $`K_+>2/\sqrt{3}1.16`$.

We now discuss the physical picture underlying the above calculations. On intuitive grounds, an attraction between the two components of the mixture ($`g_{12}<0`$) will lead to the formation of bound pairs, and might eventually lead to collapse. However, in the regime studied above, pairing takes place between hard-core bosons or a hard-core boson and a fermion, which means that as long as the repulsion between bosons of the same species is stronger, only bound pairs containing one particle of each component are possible. This picture is confirmed by considering a strong-coupling limit of the Hamiltonian (1), where first we take $`g_{11},g_{22}+\mathrm{}`$ and then $`g_{12}\mathrm{}`$ such that $`|g_{12}|g_{11},g_{22}`$ 222For a $`B+B`$ mixture, the limit of $`g_{11},g_{22}+\mathrm{}`$ but finite $`g_{12}`$ maps to an integrable model of spinful fermions solved by C. N. Yang. For $`g_{12}<0`$ the fermions form bound pairs. See e.g. M. Takahashi, Thermodynamics of One-Dimensional Solvable Models, Cambridge University Press (Cambridge, 1999). . Moreover, in a $`B+B`$ mixture pairs are also bosons, but their condensation is forbidden in 1D dimension by rigorous theorems Petrov00Andersen02 . In a $`B+F`$ mixture, the pairs are fermions, and this slightly suppresses pairing correlations relative to the $`B+B`$ case. We expect a smooth crossover from the weak coupling regime to the strong-coupling regime where tightly bound pairs form. As the strength of attraction between the two components increases, the size of the gap also grows thus preventing collapse (Fig. 1).

However, the gap can be destroyed when a sufficiently large difference in chemical potential leads to a density difference $`|\delta |>0`$. This transition is of commensurate-incommensurate type Gogolin98 . Small inhomogeneities in the density due to the trap, will not have a strong effect in the pairing provided that the density difference between the two components does not become too large. Otherwise, one can expect phase segregation into pairing and nonpairing regions, with the former occurring mainly near the center of the trap.

In conclusion, we have shown that binary mixtures of 1D quantum gases exhibit some interesting instabilities. In particular, we have found a pairing instability which may be special to one dimension and should have unusual experimental signatures unpub .

We acknowledge inspiring conversations with A. Nersesyan and M. Fabrizio, and useful correspondence with B. Svistunov. A. F. H. thanks the hospitality of DIPC and a grant from EPSRC (U.K.). M. A. C. thanks Gipuzkoako Foru Aldundia for finacial support.


format: text

text:
# Signature of stripe pinning in optical conductivity

## Abstract

The response of charge stripes to an external electric field applied perpendicular to the stripe direction is studied within a diagrammatic approach for both weak and strong pinning by random impurities. The sound-like mode of the stripes described as elastic strings moves to finite frequency due to impurity pinning. By calculating the optical conductivity we determine this characteristic energy scale for both a single stripe and an array of interacting stripes. The results explain the anomalous far-infrared peak observed recently in optical-conductivity measurements on cuprates.

Recent optical conductivity measurements in high-$`T_c`$ cuprates have revealed the existence of a strong peak in the far-infrared regime, which has been attributed to the presence of charge stripes cu1 ; paolo . The typical energy of this peak is much lower than the mid-infrared band detected in previous experiments uchida  and reproduced by numerical simulations of the Hubbard model lorenzana . A similar feature reported for slightly overdoped Bi2Sr2CuO6 paolo2  was explained in terms of a charge-ordering instability of the Fermi liquid caprara . However, these arguments cannot be extended to the underdoped regime, where the Fermi-liquid description breaks down and the charge degrees of freedom are known to order tranquada .

In this paper we propose that this anomalous peak may originate from the pinning of the zero-energy phonon-like mode associated with transverse fluctuations of stripes, which shifts to finite energies in the presence of impurities. We calculate the optical conductivity of striped systems using a diagrammatic approach, which allows us to determine the peak frequency for both weak and strong impurity pinning. In the single-stripe case, we find a peak at a frequency $`\nu `$ given by the ratio between the stripe velocity $`v`$ and a characteristic impurity length scale $`\lambda `$, which is determined by the strength of the pinning. For a stripe array, a second energy scale $`\omega _U`$, proportional to the strength $`U`$ of the stripe-stripe interaction, becomes relevant. In the weak-pinning regime, which is appropriate for describing cuprates, a second peak appears around $`\omega _U`$. When $`\omega _U\nu `$, a resonance arises and the main peak splits. Our calculations yield a good estimate for the peak frequency measured in La cuprates cu1 ; paolo ; lobo  and clarify some peculiar features observed in experiments, such as the peak-sharpening cu1  and splitting paolo  which occur when stripe-stripe interactions are important.

We begin by considering a single line of holes (stripe) embedded in a two-dimensional $`\mathrm{}\times L`$ antiferromagnetic background with lattice constant $`a`$. The transversal dimension $`\mathrm{}`$ correspond to the interstripe distance when the stripe array is explicitly considered (see below). The stripe is oriented along the $`y`$-axis and acts as a domain-wall for the underlying magnetization. The holes are assumed to move in the transverse $`x`$ direction and their dynamics is governed by the $`tJ`$ model. In the long-wavelength limit a wave-like action describes the transverse displacement $`u(y)`$ of the stripe with respect to the equilibrium position prb98 . The stripe velocity $`v=\sqrt{Jt}a`$ is related to the kinetic energy of the holes ($`t`$) and the staggered magnetization ($`J`$) of the surrounding regions. In terms of the dimensionless displacement field $`\psi =u/a`$, the action reads

$$𝒮=\frac{1}{2at}𝑑\tau 𝑑y\left[_0[\psi ]+V[\psi ]\right],$$
(1)
where $`_0[\psi ]=(\psi /\tau )^2+v^2(\psi /y)^2`$ and $`V[\psi ]`$ is the pinning potential. In analogy with the standard problem of pinning of elastic manifolds giamarchi ; cdw1 ; cdw2 , the interaction of the (striped) density of carriers $`\rho (𝐫)`$ with the short-range disorder potential $`U(𝐫)=U_0_i\delta (𝐫𝐫_i)`$ can be written as $`𝑑𝐫U(𝐫)\rho (𝐫)\rho _0U_0_i\mathrm{cos}[2\pi (ux_i)/\mathrm{}]`$, where $`𝐫_i=(x_i,y_i)`$ denotes the position of the defects, and $`\rho _0`$ is the average stripe density giamarchi . We consider that only defects localized near the stripe effectively act to pin it, so that the cosine term can be expanded around the minimum to yeld, within the quadratic approximation, $`V[\psi ]=V_0at_i(\psi (y)\beta _i)^2\delta (yy_i)`$. The impurities pin the local displacement $`\psi (y_i)`$ at a value $`\beta _i`$, with positive strength $`V_0`$, and we assume $`\mathrm{}=k_B=1`$.

In the absence of dynamical fluctuations, the stripe accommodates the impurity potential by adopting a configuration $`\psi _0(y)`$, which is a solution of the static equation of motion. The full, time-dependent solution $`\psi (y,\tau )`$ is described as a fluctuation around the equilibrium configuration $`\psi _0(y)`$, i.e., $`\psi (y,\tau )=\psi _0(y)+\varphi (y,\tau )`$. Introducing the Fourier transform $`\varphi (q,\omega _m)`$ with respect to the momentum $`q`$ and the Matsubara frequency $`\omega _m=2\pi mT`$, the reciprocal-space action for the $`\varphi `$ field becomes

$`𝒮={\displaystyle \frac{1}{2}}{\displaystyle \underset{q,\omega _m}{}}𝒟_0^1(q,\omega _m)|\varphi (q,\omega _m)|^2`$
$`+{\displaystyle \frac{V_0}{2}}{\displaystyle \underset{q,q^{},\omega _m}{}}\varphi ^{}(qq^{},\omega _m)\varphi (q,\omega _m)S(q),`$ (2)
where $`S(q)=(1/L)_ie^{iqy_i}`$ is the form factor, and $`𝒟_0(q,\omega _m)=at(\omega _m^2+v^2q^2)^1`$ is the bare Green function.

An electric field applied in the $`x`$-direction, perpendicular to the stripe, generates a current $`J=𝑑ye(\psi /t)`$. The linear charge density of the stripe is denoted by $`e/a`$. From the Kubo formula for the optical conductivity mahan , $`\sigma (\omega )=e^2i\omega 𝒟(q=0,i\omega _m\omega i\delta )`$ and by analytic continuation in the lower half-plane of the dressed Green function $`𝒟(q,\omega _m)`$ one finds

$$\mathrm{}\sigma (\omega )=e^2\omega \mathrm{}𝒟(q=0,i\omega _m\omega i\delta ).$$
(3)
In the absence of impurities, $`𝒟`$ coincides with the bare Green function $`𝒟_0`$. Its replacement in Eq. (3) yields the expected response of a massless, sound-like mode, $`\mathrm{}\sigma (\omega )=e^2at\delta (\omega )`$. The dressed Green function $`𝒟`$ may be evaluated from the Dyson equation, $`𝒟^1(q,\omega _m)=𝒟_0^1(q,\omega _m)\mathrm{\Sigma }=(\omega _m^2+v^2q^2\mathrm{\Gamma }t)/(at)`$, where we have performed a rescaling $`\mathrm{\Sigma }=\mathrm{\Gamma }/a`$ such that $`\mathrm{\Gamma }`$ has dimensions of energy. The self-energy $`\mathrm{\Gamma }`$ is evaluated by averaging the product of the several $`S(q)`$ factors, which appear in each term of the perturbative series for $`\mathrm{\Gamma }`$, over the random impurity positions mahan . Crossing diagrams, arising from scattering processes by different impurities, correspond to increasing powers of the impurity density $`n_i`$. If $`\mathrm{\Gamma }`$ is evaluated up to order $`V_0^2`$, crossing terms do not appear (Born approximation) and one obtains

$$\mathrm{\Gamma }(\omega _m)=\frac{V_0}{2}n_ia+\left(\frac{V_0}{2}\right)^2n_ia\frac{1}{L}\underset{q}{}D_0(q,\omega _m),$$
(4)
where $`𝒟_0`$ was used for calculating the second-order correction. On defining the frequency $`\omega _0=n_iv=n_ia\sqrt{Jt}`$, the dimensionless quantities $`\stackrel{~}{\omega }=\omega /\omega _0`$, $`G=\mathrm{\Gamma }t/\omega _0^2`$, and the parameter $`\alpha =V_0/(2n_iaJ)`$, Eq. (4) takes the form

$$G=\alpha \frac{i\alpha ^2}{2\stackrel{~}{\omega }}.$$
(5)
According to Eq. (3) the optical conductivity is given by

$$\mathrm{}\sigma (\omega )=\sigma _0\frac{\stackrel{~}{\omega }G^{\prime \prime }}{(\stackrel{~}{\omega }^2+G^{})^2+G^{\prime \prime 2}},$$
(6)
where $`\sigma _0=e^2at/\omega _0`$ and $`G=G^{}+iG^{\prime \prime }`$. Note that as soon as a real part is generated in $`G`$, the delta function in the optical conductivity moves to a finite frequency and its amplitude is determined by $`G^{\prime \prime }`$. Using (5), we find

$$\mathrm{}\sigma (\omega )=\sigma _0\frac{2(\omega /\nu )^2}{4(\omega /\nu )^2[(\omega /\nu )^21]^2+\alpha },$$
(7)
with $`\nu =\sqrt{\alpha }\omega _0`$. We note that $`\mathrm{}\sigma (\omega )0`$ both as $`\omega 0`$ and as $`\omega \mathrm{}`$, and has a peak at $`\omega \nu `$, as shown by the solid line in Fig. 1a. The parameter $`\alpha `$, which controls the perturbative expansion for $`G`$, depends both on the ratio $`V_0/J`$, i.e. the competition between impurity potential and elastic energy, and on the number of impurities per site, $`n_ia`$. Thus the weak-pinning limit, where only first powers of $`\alpha `$ need be considered (as in Eq. (5)), is realized for a large number ($`n_ia1`$) of weak ($`V_0/J1`$) impurity centers. In this case, Eq. (7) shows that the peak frequency scales according to $`\nu \sqrt{V_0tn_ia}`$, moving toward zero energy as the strength of the pinning potential decreases. In the evaluation of Eq. (4), we used the bare Green function $`𝒟_0`$. The divergence of the self-energy (5) at small $`\omega `$ indicates the break down of this crude approximation at low frequencies. Inclusion of the next-order term in the perturbative expansion shows that higher-order diagrams become important for $`\omega <\sqrt{\alpha }\nu /2`$, i.e. at frequencies much smaller than $`\nu `$. It has been argued cdw1  that a better result for $`G`$ may be achieved by using the dressed Green function in Eq. (5), which leads to the self-consistent equation $`G=\alpha +(\alpha ^2/2)(\stackrel{~}{\omega }^2G)^{1/2}`$. However, the resulting expression for $`\sigma (\omega )`$ shows a spurious gap below $`\omega 0.7\nu `$, which is thought to originate from the fact that the self-consistent Born approximation selects arbitrarily a subset of diagrams of higher order in $`\alpha `$ without including crossing terms, whose contribution cannot be excluded if $`n_ia`$ is large. Because we are interested primarily in the evaluation of the peak frequency $`\nu `$, and not on details of the low-frequency behavior of $`\sigma (\omega )`$ fogler , we restrict our analysis to the non-self-consistent Born approximation, which provides a reliable estimate both of the maximum of $`\sigma (\omega )`$ and of its scaling with $`\alpha `$.

A different approach must be considered in the strong-coupling limit $`\alpha 1`$, which corresponds to a small number ($`n_ia1`$) of strong ($`V_0/J1`$) pinning centers. In this case one might expect that the entire perturbative series in $`\alpha `$ should be resummed, but there exists a criterion for selecting the relevant higher-order diagrams. For $`n_ia1`$, it is possible to select the subset of diagrams of first order in $`n_i`$ after averaging over the impurity distribution ($`T`$-matrix approximation) mahan . By explicit summation of the resulting series, using the full Green function $`𝒟`$ while evaluating each diagram, one obtains a self-consistent equation for $`G`$. In the limit $`\alpha \mathrm{}`$, one finds the analytical solution $`\mathrm{}\sigma (\omega )=2\sigma _0\sqrt{\stackrel{~}{\omega }^21}/\stackrel{~}{\omega }^3`$ for $`\stackrel{~}{\omega }>1`$ and $`\mathrm{}\sigma (\omega )=0`$ for $`\stackrel{~}{\omega }<1`$ cdw2 . Thus, $`\mathrm{}\sigma (\omega )`$ displays a gap at $`\omega _0`$ and a maximum at a frequency $`\omega `$ slightly above this value. When the solution is evaluated at arbitrary $`\alpha 1`$, one finds essentially the same features. In this case the requirement of full self-consistency is crucial, and the gap feature arises from the strong pinning of the sound mode. The difference between the weak- and strong-pinning results can be understood by considering that in both cases the peak frequency corresponds to the sound mode of a free string, which is now trapped on a characteristic length scale $`\lambda `$. When few strong impurities are present, $`\lambda `$ is given by the average distance between impurities $`1/n_i`$ and the peak appears at the frequency $`\omega _0`$. However, when the number of impurities increases and $`V_0`$ decreases, the elastic energy cost of accommodating the string on a length scale $`1/n_i`$ is too high, and $`\psi `$ is pinned on the Larking length $`\lambda _{loc}`$ giamarchi ; larkin . Dimensional estimates based on Eq. (1) yield $`\lambda _{loc}v/\sqrt{V_0n_iat}`$, which indeed corresponds to the characteristic frequency $`\nu `$ in the weak-pinning regime.

Until now we have concentrated on the transverse pinning of a single stripe. However, a more realistic description of cuprates requires the study of an array of stripes. By considering a harmonic coupling strength $`U`$ between neighboring stripes, which are labeled by an index $`n=1,\mathrm{}N`$, we obtain the action hass99

$$S=\frac{1}{2at}\underset{n}{}𝑑\tau 𝑑y\left\{[\psi _n]+Ut(\psi _n\psi _{n+1})^2+V[\psi _n]\right\}.$$

Here $`\psi _n=\psi _n(y,\tau )`$ and the fluctuating field $`\varphi _n(y,\tau )\varphi (x=n\mathrm{},y,\tau )`$ is a function of the discrete $`x`$ values. By repeating the same steps as before, we observe that Eq. (2) is still valid, on taking $`q𝐪=(q_x,q_y)`$, $`S(𝐪)=(LN)^1_{ni}e^{iq_xn\mathrm{}+iq_yy_i^n}`$, where $`y_i^n`$ is the impurity coordinate at the $`n`$-th stripe, and $`𝒟_0`$ becomes

$$𝒟_0(𝐪,\omega _m)=\frac{at}{\omega _m^2+v^2q_y^2+4Ut\mathrm{sin}^2(q_x\mathrm{}/2)}.$$
(8)
The system of interacting stripes displays an anisotropic sound-like mode at long wavelengths, with an effective elastic coefficient $`U/a`$ in the $`x`$ direction, and an optical mode at $`q_x\mathrm{}=\pi `$ with energy $`\omega _U=2\sqrt{Ut}`$.

The general procedure used above to address the effect of disorder is extended readily to the case of interacting stripes. We assume that the average density of impurities does not depend on the stripe index $`n`$, i.e. $`n_i^n=n_i`$. In the Born limit, the main difference with respect to the single-stripe case arises from the additional integration in $`q_x`$ of the Green function $`𝒟_0`$ in Eq. (4). By introducing the rescaled variables and $`\overline{\omega }=\omega /\omega _U`$, one obtains

$$G=a+\frac{\alpha ^2}{\pi }\frac{\omega _0}{\omega _U}_0^{\pi /2}𝑑x\frac{1}{\sqrt{\mathrm{sin}^2x\overline{\omega }^2}}.$$
(9)
When $`\overline{\omega }>1`$ the second term in Eq. (9) contributes only to $`G^{\prime \prime }`$. In the limit of vanishing stripe-stripe interaction $`U0`$, or $`\overline{\omega }\mathrm{}`$, it reproduces the single-stripe result Eq. (5). For finite $`U`$, a second energy scale $`\omega _U`$ is relevant for determining the optical conductivity, which is defined by Eq. (3) with $`q`$ replaced by $`𝐪`$. The new feature is the divergence of $`G^{\prime \prime }`$ at $`\omega _U`$, which leads to a second peak around this frequency, where $`\mathrm{}\sigma (\omega )`$ vanishes. The vanishing of the conductivity at $`\omega _U`$ in the weak-pinning regime should indeed be expected: at this energy, neighboring stripes fluctuate in antiphase, and because the total current arises from the sum of the currents of each stripe, the net current in such a case is zero. At $`\omega _U\nu `$ or $`\omega _U\nu `$ the peak at $`\omega _U`$ is almost completely suppressed by the overall decrease of $`\mathrm{}\sigma (\omega )`$ at small and high frequencies, see dashed lines in Figs. 1a and 1e. However, when $`\omega _U\nu `$ the peak at $`\omega \nu `$ splits (Fig. 1c), signaling a resonance between the characteristic mode of the stripe array at the zone boundary and the characteristic frequency of sound-like oscillation of each stripe at the Larkin length. Dimensional estimates indicate that this resonance arises when the average stripe separation coincides with the Larkin length. Moreover, for the stripe array we observe a sharpening of the peak for $`\omega _U>\nu `$ due to the flattening of $`G^{\prime \prime }`$ around $`\nu `$ (dashed line in Fig. 1e). The softening of the peak frequency for $`\omega _U>\nu `$ scales according to $`\alpha ^2\nu /\omega _U`$, so it never moves far away from $`\nu `$.

In the strong-pinning regime the qualitative behavior of $`\mathrm{}\sigma (\omega )`$ for the stripe array is similar to that of the single stripe, and no resonance is observed. Because the interaction with the impurities is the predominant effect, the behavior of $`G`$ is almost entirely dominated by the pinning of the sound mode at $`\omega _0`$. More details about the strong-pinning regime will be presented elsewhere unpu .

We next compare our results with the available experimental data. The velocity of the unpinned sound mode may be estimated by using standard values $`tJ0.1`$ eV for cuprates. The other relevant quantity is $`n_i`$. In these materials, the chemical-substituted dopants, whose two-dimensional density is $`x/a^2`$ ($`x`$ being the doping), simultaneously provide charge carriers in the plane and act as pinning centers. As we explained at the beginning, we consider that on average the impurities within a range $`\mathrm{}`$ contribute to pin each stripe, i.e. $`n_i=x\mathrm{}/a^2`$. Neutron diffraction measurements have shown that in the underdoped regime ($`0.05<x<0.12`$) of La2-xSrxCuO4 (LSCO) $`\mathrm{}=a/2x`$ tranquada , which yields $`n_ia=1/2`$ and $`\omega _00.05`$ eV. Because $`n_ia1`$ corresponds to the weak-pinning regime, we expect that the peak frequency is given by $`\nu =\sqrt{\alpha }\omega _0`$ and thus is reduced with respect to the estimate of $`\omega _0`$, in agreement with experimental data, see Fig. 1b, 1d, 1f cu1 ; paolo . As shown in Fig. 1b-1f, the peak frequency softens as $`x`$ increases. Even though we also expect a general decreasing of $`\nu `$ due to a screening effect on $`V_0`$, a more quantitative comparison with experimental data would require the inclusion of lattice effects, which is beyond the scope of this paper.

Another interesting feature can be explained by our calculations: with increasing doping concentration, stripe-stripe interactions become more important and one should consider the stripe array. Our results predict that the peak splits when the average stripe separation $`\mathrm{}`$ is of the order of the Larkin length, see Fig. 1c. This splitting has been observed: for $`x=0.05`$ doped LSCO, a single peak appears (Fig. 1b) whereas for $`x=0.07`$ the peak splits (Fig. 1d) paolo . Although the value of $`\alpha `$ as a function of doping $`x`$ is in general not known, we may nevertheless estimate the peak frequency where the splitting arises as $`\nu =v/\mathrm{}=0.014`$ eV $`=113`$ cm-1, in excellent agreement with the measured value (Fig. 1d) paolo . Our calculations also allow us to understand the results obtained by Dumm et al. cu1  in Nd-doped LSCO (Fig. 1f): the broad peak observed in LSCO cu1  moves to slightly lower frequency and becomes sharper in the Nd-doped compound, where the stripe array is ordered tranquada . We have shown that when inter-stripe interactions increase, the peak becomes slightly softer, narrower and more intense (Fig. 1e), as observed experimentally. Note that in Ref. cu1  the feature seen in the Nd-free sample in Fig. 1f was interpreted as an anomalous large Drude peak. However, we claim that a stripe signature is still present in these data, even though a partial overlap with the single-particle response is possible. Indeed, resistivity data from Ref. cu1  give $`\sigma (0)810^3`$ $`\mathrm{\Omega }^1`$cm-1, which is much larger than the value obtained by extrapolating the supposed Drude peak. Moreover, a peak at nearly the same frequency is also seen at $`x=0.12`$ in Ref. paolo .

In addition, our results provide the correct order of magnitude for the peak height. The one-dimensional optical-conductivity unit $`\sigma _0\sigma _0^{1d}`$ converts into a three-dimensional one as $`\sigma _0^{3d}=\sigma _0^{1d}/\mathrm{}d=(2e^2/\mathrm{})(t/\omega _0)(x/d)`$, where $`d`$ is the interlayer spacing. At $`x=0.07`$, using $`d6`$ Å, we estimate the maximum of $`\sigma (\omega )`$ to be of order $`10^3\mathrm{\Omega }^1`$cm-1, as experimentally observed paolo .

Note that longitudinal transport, which is possible for half-filled stripes as in the cuprates, was neglected in our approach. However, the experimental observation of the same far-infrared peak in nickelates ni , which exhibit filled stripes, supports our hypothesis that the anomalous peak in the far-infrared regime originates from transverse fluctuations. In fact, our results could be used also to describe the nickelates, but in this case the effective $`t`$ and $`J`$ parameters within a one-band model are not known.

In conclusion, we calculated the contribution from stripe pinning to the optical conductivity in cuprates within a diagrammatic approach accounting for the presence of randomly distributed impurities. We found a peak in the far-infrared range at a frequency which has the same order of magnitude as observed in the experiments. In addition, our results explain the splitting of the peak, as well as the narrowing and increased intensity, which arise when stripe-stripe interactions are relevant.

We acknowledge fruitful discussions with D. Baeriswyl, P. Calvani, S. Caprara, T. Giamarchi, J. Lorenzana, B. Normand and J. Tranquada. We thank the authors of Refs. cu1 ; paolo  for providing us with the experimental data. This work was supported by the Swiss NSF under grant No. 620-62868.00.


text:
# Three-dimensional vortex configurations in a rotating Bose Einstein condensate

## I Introduction

Several experimental groups have produced vortices in Bose Einstein condensates (BEC) mal ; AK ; MCWD ; MCWD2 ; RK ; RBD . One type of experiments consists in imposing a laser beam on the magnetic trap holding the atoms to create a harmonic anisotropic rotating potential. For a prolate trap, it has been observed MCWD2 ; MCWD ; RBD  that when a single vortex exists, the vortex line is not straight along the axis of rotation, but bending. Theoretical works AR ; AJ  establish a simpler expression of the Gross Pitaevskii energy that only depends on the vortex lines. In AJ , it is proved that bending occurs for prolate condensates, but not for oblate ones.

Minimization algorithms GP1 ; CM  have been used to compute local minima of the Gross Pitaevskii energy and provide an evidence of the bending in the same setting as the experiment. Bending (or $`U`$) vortices are described in detail, and multiple vortex configurations are addressed in these studies.

Recently, the ENS group RBD  has further studied configurations with a single vortex line. They have observed planar bent vortices ($`U`$) but also different configurations ($`S`$). They study the length of the line, its deviation from the center and its angular momentum.

In this paper, motivated by the recent experiments at the ENS RBD , we numerically look for local minimizers of the Gross Pitaevskii energy and we want to understand the various vortex configurations observed in the experimental setting: $`U`$ vortices but also $`S`$ vortices. We compute solutions with up to 4 vortices and describe their three-dimensional structure. Different solution branches are followed and the evolution of the corresponding energy and angular momentum are shown. The framework of this study is the case of a prolate condensate where bending is an important phenomenon.

We consider a pure BEC of $`N`$ atoms confined in a harmonic trapping potential rotating along the $`z`$ axis at angular velocity $`\mathrm{\Omega }`$. The equilibrium of the system corresponds to local minima of the Gross-Pitaevskii energy in the rotating frame

$`(\varphi )`$ $`=`$ $`{\displaystyle _𝒟}{\displaystyle \frac{\mathrm{}^2}{2m}}|\varphi |^2+\mathrm{}𝛀(i\varphi ,\varphi \times 𝐱)`$
$`+{\displaystyle \frac{m}{2}}\omega _x^2(x^2+\alpha ^2y^2+\beta ^2z^2)|\varphi |^2+Ng_{3D}|\varphi |^4,`$
where $`g_{3D}=4\pi \mathrm{}^2a/m`$ and the wave function $`\varphi `$ is normalized to unity $`_𝒟|\varphi |^2=1`$. Here, for any complex quantities $`u`$ and $`v`$ and their complex conjugates $`\overline{u}`$ and $`\overline{v}`$, $`(u,v)=(u\overline{v}+\overline{u}v)/2`$.

For numerical applications, it is more convenient to rescale the variables as follows: $`𝐫=𝐱/R`$, $`u(𝐫)=R^{3/2}\varphi (𝐱)`$, where $`R=d/\sqrt{\epsilon }`$ and

$$d=\left(\frac{\mathrm{}}{m\omega _x}\right)^{1/2},\epsilon =\left(\frac{d}{8\pi Na}\right)^{2/5},\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }/(\epsilon \omega _x).$$

In this scaling the Thomas-Fermi limit of $`u`$ is

$$\rho _{\text{TF}}(𝐫)=\rho _0(x^2+\alpha ^2y^2+\beta ^2z^2).$$
(1)
Then, we use the dimensionless energy introduced in AR

$$E(u)=H(u)\stackrel{~}{\mathrm{\Omega }}L_z(u),$$
(2)
with

$$H(u)=_𝒟\frac{1}{2}|u|^2\frac{1}{2\epsilon ^2}\rho _{\text{TF}}|u|^2+\frac{1}{4\epsilon ^2}|u|^4,$$
(3)
$$L_z(u)=i_𝒟\overline{u}\left(y\frac{u}{x}x\frac{u}{y}\right),$$
(4)
defined in the domain $`𝒟=\left\{\rho _{\text{TF}}(𝐫)0\right\}.`$

### I.1 Numerical method

In the present study we compute critical points of $`E(u)`$ by solving the norm-preserving imaginary time propagation of the corresponding equation:

$$\frac{u}{t}\frac{1}{2}\mathrm{\Delta }u+i(\stackrel{~}{𝛀}\times 𝐫).u=\frac{1}{2\epsilon ^2}u(\rho _{\text{TF}}|u|^2)+\mu _\epsilon u,$$
(5)
with $`u=0`$ on $`𝒟`$ and $`\mu _\epsilon `$ the Lagrange multiplier for the norm constraint $`_𝒟|u|^2=1`$. A hybrid 3 steps Runge-Kutta-Crank-Nicolson scheme Orlandi  is used to advance the equation in time:

$$\frac{u_{l+1}u_l}{\mathrm{\Delta }t}=a_l_l+b_l_{l1}+c_l\mathrm{\Delta }\left(\frac{u_{l+1}+u_l}{2}\right),$$
(6)
where $``$ contains the remaining non-linear terms. The corresponding constants for every step ($`l=1,2,3`$) are :

$$\begin{array}{ccc}a_1=8/15,\hfill & a_2=5/12,\hfill & a_3=3/4,\hfill \\ b_1=0,\hfill & b_2=17/60,\hfill & b_3=5/12,\hfill \\ c_1=8/15,\hfill & c_2=2/15,\hfill & c_3=1/3.\hfill \end{array}$$

The resulting semi-implicit scheme is second order time accurate and allows reasonably large time steps, making it appropriate for long time integration. The large sparse matrix linear systems resulting from the implicit terms are solved by an alternating direction implicit (ADI) factorization technique.

For the spatial discretization we use finite differences on a Cartesian uniform mesh with periodic boundary conditions in all directions. To accurately resolve sharp gradients of the variable in presence of vortices, low numerical dissipation and very accurate schemes are required for the spatial derivatives. A sixth-order compact finite difference scheme Lele  with spectral-like resolution was chosen to this end.

### I.2 Physical and numerical parameters

The values of constants in (5) are set to $`\epsilon =0.02,\alpha =1.06,\beta =0.067,`$ corresponding to the experiments of the ENS group CM ; MCWD  ($`m=1.44510^{26}kg`$, $`a=5.810^{11}m`$, $`N=1.410^5`$ and $`\omega _x=1094s^1`$). The angular frequency $`\mathrm{\Omega }`$ will be varied from 0 to the maximum value of $`0.9\omega _x`$, for which no deformation of the condensate has to be taken into account.

Equation (5) is propagated in imaginary time until the evolution of the energy (2) has a gradient in time smaller than $`10^6`$. For the considered range of $`\mathrm{\Omega }`$, the numerical domain is fixed to an elongated box $`(x,y,z)[0.6,0.6]\times [0.6,0.6]\times [8.5,8.5]`$. A refined grid using $`72\times 72\times 510`$ nodes is employed, which is sufficient to achieve grid-independence for all considered numerical experiments.

Different initial conditions are used in to trigger single or multiple vortex configurations and follow the corresponding branches as $`\mathrm{\Omega }`$ is varied. The simplest initial condition assumes a steady-state solution $`u(x,y,z)=\sqrt{\rho _{\text{TF}}(x,y,z)}`$ and is useful to study vortex-free configurations and their degeneracy into multiple vortex configurations when increasing the value of $`\mathrm{\Omega }`$. Initial conditions with vortices are obtained by superimposing to the steady-state a simple ansatz for the vortex. For example, an initial condition with a centered straight vortex of radius $`\epsilon `$ is obtained by imposing

$`u(x,y,z)=\sqrt{\rho _{\text{TF}}}u_\epsilon ,`$ (7)
$`u_\epsilon =\sqrt{0.5\left\{1+tanh\left[{\displaystyle \frac{4}{\epsilon }}\left(r\epsilon \right)\right]\right\}}exp(i\phi ),`$
where $`(r,\phi )`$ are the polar coordinates in the $`(x,y)`$ plane. The 3D shape of the vortex can be easily modified by shifting the center $`r_0`$ of the vortex in successive $`(x,y)`$ planes; for instance, to obtain a planar S shape vortex, the following function can be used:

$$r_0(z)=\{\begin{array}{cc}1+tanh\left[\alpha _v\left(1+\frac{z}{\beta _v}\right)\right]/tanh(\alpha _v),\hfill & z<0\hfill \\ 1+tanh\left[\alpha _v\left(1+\frac{z}{\beta _v}\right)\right]/tanh(\alpha _v),\hfill & z0\hfill \end{array}$$

The constants $`\alpha _v,\beta _v`$ control, respectively, the curvature and the height of the vortex.

We first focus on single vortex configurations and describe later multi vortex configurations.

## II Single vortex lines

We have observed three different types of single vortex configurations as shown in figure 1: planar $`U`$ vortices, planar $`S`$ vortices and non-planar $`S`$ vortices. The $`U`$ vortices are the bent vortices computed in GP1 ; CM  and theoretically studied in AR ; AJ . They are global minimizers of the energy. The $`S`$ configurations were observed experimentally very recently RBD  and are only local minimizers of the energy.

### II.1 $`U`$ vortex

The $`U`$ vortex is a planar vortex formed of 2 parts: the central part is a line which stays on the $`z`$ axis and the outer part reaches the boundary of the condensate perpendicularly. When $`\mathrm{\Omega }`$ increases, the central straight part gets longer (figure 2) and the angular momentum ($`L_z`$) increases to 1 (figure 3).

The $`U`$ vortex is obtained by starting the simulation with an initial condition containing a straight vortex away from the $`z`$ axis. In fact, the $`U`$ vortex lies either in the $`xz`$ or $`yz`$ plane. Starting with an initial condition which is not in one of these plane yields a final state in the $`yz`$ plane, which is the plane closest to the $`z`$ axis.

The shape of the the $`U`$ vortex and its preferred location in the $`yz`$ plane can be analyzed using the approximate energy derived in AR ; AJ : setting the vortex free solution to 0 energy, then the energy of a vortex line $`\gamma `$ can be approximated by

$$_\gamma =_\gamma \rho _{\text{TF}}𝑑lC\mathrm{\Omega }_\gamma \rho _{\text{TF}}^2𝑑z,$$
(8)
where $`C`$ is a constant which depends on the experimental parameters and $`\rho _{\text{TF}}`$ is given by (1). If $`\gamma `$ is not in the $`xz`$ or $`yz`$ plane, then one can construct small perturbations of $`\gamma `$ that preserve $`\rho _{\text{TF}}`$ and lower the energy. This implies that $`\gamma `$ cannot be a critical point of the energy because the gradient is not zero. Of course, if the ellipticity of the cross section is small, the gradient is small, which may allow to observe these configurations.

In order to understand the existence of the straight central part of the $`U`$ vortex, one can also refer to the analysis of AJ : from equation (8) we can infer that a vortex line with a lower energy than the vortex free solution is obtained when the quantity $`\rho _{\text{TF}}C\mathrm{\Omega }\rho _{\text{TF}}^2`$ is negative, i.e. $`C\mathrm{\Omega }\rho _{\text{TF}}>1`$. Let $`\overline{\mathrm{\Omega }}`$ be such that $`C\overline{\mathrm{\Omega }}\rho _0=1`$; it corresponds to the 2d critical velocity for the existence of a vortex in the plane $`z=0`$. For $`\mathrm{\Omega }`$ close to $`\overline{\mathrm{\Omega }}`$, the inner region where $`C\mathrm{\Omega }\rho _{\text{TF}}>1`$ is concentrated near the center of the condensate. In this region, the vortex line has to be straight (see AJ ). This straight part is getting longer as $`\mathrm{\Omega }`$ increases since the region where $`C\mathrm{\Omega }\rho _{\text{TF}}>1`$ is getting bigger. This region corresponds to $`\mathrm{\Omega }>\mathrm{\Omega }_{2d}(z)`$, where $`\mathrm{\Omega }_{2d}(z)`$ is the critical velocity for the existence of a vortex in the 2 dimensional section where $`z`$ is constant. In the outer region, the vortex reaches the boundary using the shortest path.

Figure 3 shows the energy and angular momentum variation with $`\mathrm{\Omega }`$ for the single vortex configurations. The $`U`$ vortices exist only for $`\mathrm{\Omega }`$ bigger than a critical value $`\mathrm{\Omega }_c=0.42\omega _x`$. It is interesting to note that at $`\mathrm{\Omega }_c`$, the energy of the $`U`$ vortex is bigger than the energy of the vortex free solution (we have set to zero the energy of the vortex free solution). A zoom in this region shows that $`\mathrm{\Omega }_c`$ is very close to the angular velocity $`\mathrm{\Omega }_1`$ for which the energy of the vortex free solution is equal to the energy of the $`U`$ vortex.

Figure 3 also shows that the angular momentum $`L_z`$ of the $`U`$ vortex for $`\mathrm{\Omega }=\mathrm{\Omega }_c`$ does not go to 0. This suggests that in fact there could be another $`U`$ solution for $`\mathrm{\Omega }>\mathrm{\Omega }_c`$. Using an ansatz, another type of $`U`$ solution is obtained in CM  which is a saddle point of the energy: it is away from the axis and has lower angular momentum. In AJ , it is proved rigorously that for $`\mathrm{\Omega }`$ small, there is no $`U`$ as a critical point of the energy.

For an initial condition with a straight vortex centered on the $`z`$ axis, if $`\mathrm{\Omega }<0.8\omega _x`$, the straight vortex is unstable and the final configuration is a $`U`$, but if $`\mathrm{\Omega }>0.8\omega _x`$, the straight vortex is stable. This is in agreement with the result of AJ  where the local stability of the straight vortex for $`\mathrm{\Omega }`$ larger is proved.

For small $`\mathrm{\Omega }`$, the $`U`$ vortex disappears and a vortex-free configuration is obtained, while for $`\mathrm{\Omega }`$ larger the $`U`$ vortex degenerates into a three-vortex configuration (described later).

### II.2 $`S`$ vortex

Motivated by the experiments of RBD , we compute new critical points of the energy, which are $`S`$ configurations (see figure 1). Several numerical experiments were performed, starting from different initial conditions containing an ansatz for the $`S`$ vortex (see section I.2).

The planar $`S`$ can be regarded as a $`U`$, with the half-part in the plane $`z<0`$ rotated with respect to the $`z`$ axis by 180 degrees (see figure 4). The non planar $`S`$ are such that the projections of the branches on the $`xy`$ plane are orthogonal, i.e. the rotation of the branches is of 90 degrees. We could check that non planar $`S`$ configurations with an angle between the branches different from 90 degrees do not exist.

As already mentioned for the $`U`$ vortex, stable planar $`S`$ configurations lie either in the $`xz`$ or $`yz`$ plane. As for the $`U`$, this can be explained using the limiting energy obtained in AJ  and considering separately the upper or lower part of the $`S`$. As soon as the cross section is not a disc, if the upper or lower branch of the $`S`$ configuration does not lie in the $`xz`$ or $`yz`$ plane, then the gradient of vortex line energy (8) can never be zero when $`\gamma `$ is varied.

The $`S`$ vortices exist for all values of $`\mathrm{\Omega }`$ while the $`U`$ only exist for $`\mathrm{\Omega }>\mathrm{\Omega }_c`$. When $`\mathrm{\Omega }`$ decreases, the extension of the $`S`$ along the $`z`$ axis goes downwards as shown in figure 5, the angular momentum decreases to 0 (figure 3) and the vortex tends to the horizontal axis. Note that a vortex along the horizontal axis has $`L_z=0`$, but a positive energy. On the other side, when $`\mathrm{\Omega }`$ increases, the $`S`$ gets straighter and it tends to the vertical axis.

The global minimum of the energy is never an $`S`$. But the difference in energy (and angular momentum) between $`U`$ and $`S`$ vortices is very small, as illustrated in figure 3 because an $`S`$ vortex is almost like a $`U`$ with a half-part rotated by 180 degrees.

### II.3 Minimizer with fixed $`L`$

As pointed out in RBD , the minimization problem which is related to the experiments, is rather to minimize $`H`$ (see (3) while fixing $`L_z`$, rather than minimizing $`E=H\mathrm{\Omega }L_z`$. This has been studied in the 2 dimensionnal setting inBS . One can notice that if a given configuration with $`H=h`$ and $`L_z=l`$ minimizes $`E=H\mathrm{\Omega }L_z`$ for some $`\mathrm{\Omega }`$, then $`h`$ minimizes $`H`$ under the constraint that $`L_z=l`$: indeed if $`H^{}=H(u)`$ with $`L_z(u)=l`$, then $`H^{}\mathrm{\Omega }lh\mathrm{\Omega }l`$, since $`(h,l)`$ minimizes $`E`$, and this implies that $`H^{}h`$. Moreover $`\mathrm{\Omega }`$ is the slope to the curve $`H(L_z)`$ at the point $`(h,l)`$ and the property of minimizing $`E`$ that is for all $`h^{}`$, $`l^{}`$,

$$h^{}\mathrm{\Omega }l^{}h\mathrm{\Omega }l$$

implies that the curve $`H(L_z)`$ lies above its tangent at this point.

We have plotted $`H`$ as a function of $`L_z`$. We can check that the curve is convex, and above its tangent, which is consistent with the fact that we have computed minimizers of the energy.

We know that the $`U`$ solution exists for $`\mathrm{\Omega }\mathrm{\Omega }_c`$ and has $`L_z>0.4`$. For $`L_z<0.4`$, we expect that the process of minimizing $`H`$ with fixed $`L_z`$ would produce $`U`$ vortices and the curve $`H(L_z)`$ should be concave in this region. In AJ , we have proved that for $`L_z`$ close to 0, $`HCL_z^{2/3}`$, which is a first indication to the concavity of the curve.

## III Multiple vortices

Multiple vortex configurations are obtained based upon different numerical strategies. The first one is to start the computation from a vortex-free steady state and to abruptly increase $`\mathrm{\Omega }`$ to a very hight value; multiple vortices are thus obtained. The second strategy is to generate an initial condition with vortices as described in section I.2 (the advantage being the control of the shape and initial arrangement of the vortices).

Both techniques are used to follow solution branches with two, three or four vortices in the condensate. Figures 7 and 8 display energy and angular momentum vs $`\mathrm{\Omega }`$ for all studied configurations.

### III.1 3 vortices

When $`\mathrm{\Omega }`$ is increased, the single $`U`$ vortex solution switches to a 3 vortex configuration ($`\mathrm{\Omega }=0.9\omega _x`$). As shown in figure 9a, the configuration is invariant under rotation in a central plane near $`z=0`$ but not near the edges. For large $`\mathrm{\Omega }`$, three-dimensional views show (figure 9 a,b) that there are 2 vortices of similar size and a longer one which is bending near the boundary. For $`\mathrm{\Omega }=0.8\omega _x`$, all vortices display contorted shapes (figure 9c), very similar to those reported in GP1 . Let us point out that the angular momentum of all these 3 vortex configurations is lower than 3 (see figure 8).

When we put as initial condition a configuration with 3 identical $`U`$ vortices at 120o, in the final state, one of them gets a little longer (figure 10a) and the symmetry is lost. This configuration has almost the same energy and angular momentum as the configuration displayed in figure 9b. In exchange, the initial condition with three straight vortices on the $`x`$-axis has its symmetry preserved (figure 10 b), but with a higher energy that the previous one.

When further decreasing $`\mathrm{\Omega }`$, the 3 vortex branch switches to a 2-vortex displaying irregular shapes (figure 11).

### III.2 2 vortices

The two-vortex branch presented in this section was obtained by starting from a vortex-free solution and suddenly increasing $`\mathrm{\Omega }`$ to a value of $`0.8\omega _x`$. The configuration is planar and symmetric, like twice a single $`U`$ vortex, but away from the axis (there is a repulsion between the lines).

When $`\mathrm{\Omega }`$ increases, the lines are almost straight and get closer to each other. This is in agreement with the fact that when $`\mathrm{\Omega }`$ gets large, the straight vortex is a local minimizer of the energy. Hence the bending is no longer the important phenomenon.

We recall that decreasing $`\mathrm{\Omega }`$ from a configuration with 3 vortices, we obtained 2 vortices which are not symmetric, one being longer than the other (figure 11). This configuration has slightly bigger energy than the 2 symmetric vortices.

### III.3 4 vortices

Starting from an initial condition without vortices and increasing $`\mathrm{\Omega }`$ to $`0.86\omega _x`$, we have obtained stable configurations with 4 curved vortices (figure 13 a). When decreasing $`\mathrm{\Omega }`$, this configuration rapidly degenerates into a three-vortex state. For lower $`\mathrm{\Omega }`$ we could obtain stable configurations with four symmetric vortices (figure 13 b), but with bigger energy. The location of the vortices in the plane $`z=0`$ is the same.

We have to point out that with the initial condition of 4 identical vortices, the symmetry is preserved as displayed in figure 13b, which was not the case for 3 vortices.

## IV Concluding remarks

We have studied different vortex configurations in a prolate Bose Enstein condensate by solving the Gross Pitaevskii equation. We have computed $`U`$ and $`S`$ vortices, motivated by the recent experiments of RBD . Our computations involve a parameter $`\epsilon `$, which is small when the number of atoms $`N`$ is large. Decreasing $`\epsilon `$, that is increasing the number of atoms forces the vortex lines to be almost straight in their central part, while for $`\epsilon `$ larger, the central straight part is not so obvious as in some figures of GP1 .

We have found that the $`S`$ vortices are only local minimizers of the energy and exist for all values of the angular velocity $`\mathrm{\Omega }`$, while $`U`$ vortices are global minimizers existing for $`\mathrm{\Omega }\mathrm{\Omega }_c`$. A planar $`S`$ vortex can be regarded as a $`U`$ vortex with a half-part rotated by 180o. Moreover, $`U`$ or planar $`S`$ vortices lie only in the $`xz`$ or $`yz`$ plane while non planar $`S`$ vortices exists only for an angle of 90o between the two branches.

We have followed the branches of solutions when varying $`\mathrm{\Omega }`$ and found configurations with two, three and four vortices.


text:
# Ultraviolet structure in the lensed QSOs 0957+561 1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA Inc under NASA contract NAS5-26555 2footnote 22footnote 2Guest User, Canadian Astronomy Data Centre, which is operated by the Dominion Astrophysical Observatory for the National Research Council of Canada’s Herzberg Institute of Astrophysics

## 1 Introduction and observations

Q0957+561 is a bright gravitationally lensed QSO pair at redshift 1.415. The two images are 6.2” apart, and deep imaging reveals a chief lensing galaxy between them, which lies in a cluster at redshift 0.36. The lensing model and QSO host galaxy morphology have been discussed recently by Keeton et al (2000), and their paper contains a good bibliography of related earlier work. The QSO UV spectra reveal that they have damped Ly$`\alpha `$ absorption and other associated absorbers that arise in the QSO host galaxy (see Michalitsianos et al 1997, and Dolan et al 2000). Repeated UV nuclear spectroscopy has revealed changes and differences in the UV absorbers that are discussed in the above references.

In order to look for extended Ly$`\alpha `$ in the vicinity of the QSO, and to study the structure of the nuclear regions, direct images and slitless spectra were taken with the NUV detector of the Space Telescope Imaging Spectrograph (STIS). Both QSO images were in the field of view. Table 1 summarizes the observations, and also others taken with the same STIS grating, which were obtained from the HST archive. This paper discusses the non-nuclear spectra and UV imaging, which have not been covered in the other work referred to above.

## 2 UV images

The UV images are shown in Figure 1, along with images of stars taken with the same detector. For comparison, several star images, with exposure times and signal levels similar to the QSO images, were taken from the HST archive. While the star images all show deviations from circular symmetry, due to HST tracking jitter, they are all smaller and rounder than the QSO images, as can be seen in Figure 1 and also from the luminosity profiles in Figure 2. Figure 3 shows contours of the two QSO images.

Figure 1 shows both QSOs, rotated to the conventional orientation, the most non-circular star image, and the average image from several stars. The images have been smoothed with a gaussian of $`\sigma `$=2 ‘pixels’ (0.025”) to reduce the granularity of the photon-counting MAMA detector, but this does not significantly affect the spatial resolution (FWHM about 5 pixels) which is determined by the HST optics. The QSOs both have very similar ‘kidney-shaped’ inner profiles, and fainter radial structures, and match each other best with a small relative rotation, depending on whether the cores or radial spokes are fit. No such structure is seen in any star image, and the jitter in the QSO observation does not indicate a difference from the star observations. Thus, we conclude that there is resolved structure in the inner 0.3 arcsec diameter of the QSO images, and that the lensing distortion and rotation within this spatial scale is small. The radial structures are seen out to distances of about 0.4 arcsec, and do show larger differences in extent and relative brightness.

The luminosity profiles are derived from the IRAF ‘ellipse’ task, and plot the semi-major axis of the best fit ellipse to each signal level. The star profile is the mean profile from all the star images. The two spatial axis plots show that these inner profiles are not exponential or power-law. The QSO profiles are resolved at essentially all radii: the kink in the star profiles is caused by the HST PSF structure that shows up most clearly in the unresolved stellar images. We can also see that the B image is larger than the A image in the inner parts (by some 20% in radius, normalised to the peaks). The B image is brighter, but only by 1%, as the peak signal is lower. We discuss the image morphology in the last section.

There is no sign of the lensing galaxy or any other object in the 25 arcsec field of the UV image, as expected for the redshifted SED of a cluster galaxy at redshift 0.36. The only WFPC2 images of the QSOs have 0.1” pixels and show no resolved structure, but do reveal the lensing galaxy.

## 3 Slitless spectra

The dispersed image contains spatially resolved spectra of both QSO images. As these have no wavelength calibration, it was necessary to process the raw image into one with spatial and wavelength axes, separately for each QSO. The wavelength scale was established by using the nuclear emission lines and assuming the QSO emission-line redshift of 1.415. These are described together with the standard long-slit spectra (see Table 1) in the next section. Figure 4 shows the nuclear spectra, and off-nuclear spectra, also described in the next section. As the nuclear UV spectra have been discussed in detail by the papers noted above, this paper does not investigate them further.

The dispersed image was inspected carefully for signs of Ly$`\alpha `$ emission elsewhere in the field. This might arise from emission-line gas at the QSO redshift (i.e. associated with the host galaxy), either on the scale of the host galaxy, or close to the QSO nucleus. The lensing of the QSO and environments means that this may be found anywhere between the two QSO nuclei and at some distance beyond. The background is remarkably uniform, and no significant extra flux is seen in the likely regions, compared with the rest of the image. Close to the nuclei, Ly$`\alpha `$ emission is not noticeably extended beyond the rest of the UV continuum. Thus, the data do not reveal any detached Ly$`\alpha `$ emission line gas around the QSO.

## 4 QSO off-nuclear spectra

The cross-sections of the UV spectra are shown in Figure 5. The resolved width of the nuclear regions cause the slitless spectra to be wider than the 0.2” slit spectra. Further, the extended structure seen in the UV images in Figure 1, will give rise to more off-nuclear signal in the slitless spectra. Both these are seen. The closest comparison is between the two spectra of QSO A with dispersion angles 110o and 114o. There are additional uncertainties in the exact positioning of the nuclear region with the 0.2” slits.

It was found that the regions outside 0.1” in the slit spectra and outside 0.21” in the slitless spectra are free of nuclear contamination, by noting that they do not show the damped Ly$`\alpha `$ absorption. The regions inside that in the slitless spectra (dotted lines in Figure 4), are slightly contaminated, and Figure 3 shows the effect of subtracting a scaled nuclear spectrum just sufficient to remove the damped Ly$`\alpha `$ absorption.

The off-nuclear spectra for components A and B in the slitless data are very similar, extracted as shown in Figure 5. Thus, Figure 3 shows the combined A+B spectra. There is clearly off-nuclear Ly$`\alpha `$, which we discuss in more detail below, but we also note the possible presence of N V (1238-42Å), and an associated UV continuum, which has the same spatial extent as the line emission. Both suggest the presence of a young stellar population. Extraction of off-nuclear spectra from the G430L and G750L observations (see Table 1), shows a noisy broad emission at redshifted C IV (1550Å) which also suggests a distributed young population. At longer wavelengths, the off-nuclear spectra have a flater (i.e. redder) SED, which presumably arises from an older stellar population as well. At the QSO redshift, we can see rest wavelengths out to about 4000Å, but the signal is too weak and noisy to reveal line features.

Figure 6 shows the Ly$`\alpha `$ emission separately for all off-nuclear spectra, and both QSO components. The combination of the spatial extent of the off-nuclear emission-line regions, the slit width (or lack of slit), and the uncertain placement of the nucleus within the slit, make velocity measures of the line rather hard to interpret. However, the width and offset of the lines from the QSO redshift (shown dotted in) indicates that some Ly$`\alpha `$ emission gas has velocities at least of order 1000 km s-1 with respect to the nucleus. The W spectra at all orient angles indicate redshifted emission that corresponds to radial velocities of several thousand km s-1. As noted above, the behaviour in components A and B is similar and suggests that the lensing does not distort the spatially resolved near-nuclear regions very much.

## 5 Discussion

The far-UV rest-frame imaging and spectra of Q0957+561 show that the QSOs are resolved over scales of about 1 arcsec, and that the two lensed components have very similar morphology out to radii of 0.3 arcsec, although component B is larger by about 20%. There are also fainter radial features that extend over diameters of 1.5 arcsec, that do show differences between A and B. There is no sign of the lensing galaxy, or of more extended L$`\alpha `$ in the UV data.

The radio structure is of interest, and known to comparable spatial resolution. Conner, Lehar, and Burke (1992) published VLBI maps that resolve jets from both A and B over some 100mas. Garrett et al (1994) show these in more detail, and there is larger structure mapped with the VLA by Greenfield, Roberts, and Burke (1985). The 100mas jets are at position angles 23o and 17o for A and B, respectively, although they are both near 14o in the inner 10mas, and the A jet curves to PA nearly 90o over 5 arcsec. The radio flux magnification (B/A) is close to 0.75 (see e.g. Press and Rybicki 1998), while the optical values are near to unity. This disagreement is considered to be due to microlensing, or possibly the nearness of B to a lensing caustic.

The UV images reported here do have an extension with position angle about 25o, over the innermost 100mas, which suggests that they are associated with the radio jets. They do also have extensions with similar size, but less sharpness, at angles -60o and -170o, as seen in Figure 3. The best fit between A and B over the inner 200mas radius is with a 4.5o rotation clockwise of B, which is the opposite direction of the radio jet angle differences. As noted, the UV flux magnification is 1.01 and the size is 1.2 (both B/A). Thus, it appears that some of the UV structure is aligned with the radio jets, but the difference in magnification is similar to that at optical wavelengths.

The off-nuclear spectra indicate that there are clouds of emission line gas within the QSO host galaxy with relative velocities of the order of 1000 km s-1, and possibly significantly higher. This suggests association with the damped absorber seen in the nuclear spectra, which indicates outflows of about 2500 km s-1. The off-nuclear Ly$`\alpha `$ flux is about 10-12 erg s-1, (10% of the nuclear line BLR flux), and the equivalent widths similar to the nuclear at about 70Å. The off-nuclear spectra also indicate the presence of a young hot stellar population, as well as a redder population that begins to dominate above rest wavelengths 4000Å.

Further high resolution imaging and spectra should provide more information on the spatial structure, and possibly associate it with changes seen in the nuclear absorption spectra. The inner morphology also may provide further constraints on lensing models for the QSO.

References

Conner S.R., Lehar J., Burke B.F., 1992, ApJ, 387, L61

Dolan J.F., Michalitsionos A.G., Nguyen Q.T., Hill R.J., 2000, ApJ, 539, 111

Garrett M.A., Calder R.J., Porcas R.W., King L.J., Walsh D, Wilkinson P.N., 1994, MNRAS, 270, 457

Greenfield P.E., Roberts D.H., Burke B.F., 1985, ApJ, 293, 370

Keeton, C.R., et al, 2000, ApJ, 542, 74

Michalitsianos A.G. et al, 1997, ApJ, 474, 598

Press W.H., Rybicki G.B., 1998, ApJ, 507, 108

Captions to figures

1. UV images (rest frame 900-1300Å) of the QSOs (upper), compared with stellar images taken at the same observed wavelengths and detector (lower). Images are all 8” on a side, have similar signal levels, and the same display lookup table. The QSOs are oriented N up and E to the right. The left QSO image is component A and the right one is B. The left star is the ‘worst’ one found in the archive, and the right one is the composite of four stars. Note the similar core structure in the QSOs, larger than the stars, and also the radial faint rays emerging from the QSO images.

2. Radial profiles of the QSO images and the mean of the star images, plotted on R1/4 and R scales. The QSO images are extended, with the B image being larger, but with a very similar profile. Formal error bars from the ellipse task are plotted.

3. Contour plots of the two UV images. North is up and East to the left. The extension to the NE is aligned with the VLBI radio jets. The two images align best with a rotation of 4.5o.

4. Spatially resolved UV spectra of Q0957+561. The wavelength scale is the rest frame, and spectra step across the nucleus as shown in the cross sections in Figure 5. As the inner QSO structure (and the extracted spectra) are so similar, the off-nuclear spectra are combined from both A and B components. The central panel shows the two nuclear spectra (scaled down by some 20 times compared with the others). Note the small differences in the damped absorption lines. The two extractions nearest the nucleus are probably contaminated by nuclear light, and are also shown after nuclear subtraction sufficient to remove the Ly$`\alpha `$ absorption.

5. Cross-sections of all spectra in Table 1. The vertical dotted lines show the extraction windows for the spectra in Figure 4. The lower, slit spectra are narrower due to the slit width. The off-nuclear spectra in Figure 6 are extracted using the light outside the inner dotted lines in the upper panel. The different dispersion orientations and uncertainties in the exact slit placement cause the different amounts of off-nuclear light.

6. Panel of Ly$`\alpha `$ emission separately for the two QSO images. These spectra are normalised to a continuum fitted to the entire UV spectrum in all cases. The upper four (dashed) spectra are from the slitless data. Note the Ly$`\alpha `$ profile changes, and the possible appearance of N V away from the nucleus.