text:
# Relations of al Functions over Subvarieties in a Hyperelliptic Jacobian

## 1. Introduction

For a hyperelliptic curve $`C_g`$ given by an affine curve $`y^2=_{i=1}^{2g+1}(xb_i)`$, where $`b_i`$’s are complex numbers, we have a Jacobian $`𝒥_g`$ as a complex torus $`^g/\mathrm{\Lambda }`$ by the Abel map $`\omega `$ \[Mu\]. Due to the Abelian theorem, we have a natural morphism from the symmetrical product $`\mathrm{Sym}^g(C_g)`$ to the Jacobian $`𝒥_g\omega [\mathrm{Sym}^g(C_g)]/\mathrm{\Lambda }`$. As zeros of an appropriate shifted Riemann theta function over $`𝒥_g`$, the theta divisor is defined as

$$\mathrm{\Theta }:=\omega [\mathrm{Sym}^{g1}(C_g)]/\mathrm{\Lambda }$$

which is a subvariety of $`𝒥_g`$. Similarly, it is natural to introduce a subvariety

$$\mathrm{\Theta }_k:=\omega [\mathrm{Sym}^k(C_g)]/\mathrm{\Lambda }$$

and a sequence,

$$\mathrm{\Theta }_0\mathrm{\Theta }_1\mathrm{\Theta }_2\mathrm{}\mathrm{\Theta }_{g1}\mathrm{\Theta }_g𝒥_g$$

Vanhaecke studied the structure of these subvarieties as stratifications of the Jacobian $`𝒥_g`$ using the strategies developed in the studies of the infinite dimensional integrable system \[V1\]. He showed that these stratifications of the Jacobian are connected with stratifications of the Sato Grassmannian. Further Vanhaecke investigated Lie-Poisson structures in the Jacobian in \[V2\]. He showed that invariant manifolds associated with Poisson brackets can be identified with these strata; it implies that the strata are characterized by the Lie-Poisson structures. He also showed that the Poisson brackets are connected with a finite-dimensional integrable system, Henon-Heiles system. Following the study, Abenda and Fedorov \[AF\] investigated these strata and their relations to Henon-Heiles system and Neumann systems.

On the other hand, functions over the embedded hyperelliptic curve $`\mathrm{\Theta }_1`$ in a hyperelliptic Jacobian $`𝒥_g`$ were also studied from viewpoint of number theory in \[C, G, Ô\]. In \[Ô\], Ônishi also investigated the sequence of the subvarieties, and explicitly studied behaviors of functions over them in order to obtain higher genus analog of the Frobenius-Stickelberger relations for genus one case. Though Vanhaecke, Abenda and Fedorov found some relations of functions over these subvarieties explicitly using the infinite universal grassmannians and so-called Mumford’s $`UVW`$ expressions \[Mu\], Ônishi gave more explicit relations on some functions over the subvarieties using the theory of hyperelliptic functions in the nineteenth century fashion \[Ba1, Ba2, Ba3\].

In this article, we will also investigate some relations of functions over the subvarieties based upon the studies of the hyperelliptic function theory developed in the nineteenth century \[Ba2, Ba3, W\]. Especially this article deals with the “sine-Gordon equation” over there.

Modern expressions of the sine-Gordon equation in terms of Riemann theta functions were given in \[\[Mu\] 3.241\],

(1.1)
$$\frac{}{t_\mathrm{P}}\frac{}{t_\mathrm{Q}}\mathrm{log}([2\mathrm{P}2\mathrm{Q}])=A([2\mathrm{P}2\mathrm{Q}][2\mathrm{Q}2\mathrm{P}]),$$

where $`\mathrm{P}`$ and $`\mathrm{Q}`$ are ramified points of $`C_g`$, $`A`$ is a constant number, $`[D]`$ is a meromorphic function over $`\mathrm{Sym}^g(C_g)`$ with a divisor $`D`$ for each $`C_g`$ and $`t_\mathrm{P}^{}`$ is a coordinate in the Jacobi variety such that it is identified with a local parameter at a ramified point $`\mathrm{P}^{}`$ up to constant.

In the previous work \[Ma\], we also studied (1.1) using the fashion of the nineteenth century. In \[W\] Weierstrass defined $`\mathrm{al}`$ function by $`\mathrm{al}_r:=\gamma _r\sqrt{F_g(b_r)}`$ and $`F_g(z):=(x_1z)\mathrm{}(x_gz)`$ over $`𝒥_g`$ with a constant factor $`\gamma _r`$. Let $`\gamma _r=1`$ in this article. As Weierstrass implicitly seemed to deal with it, (1.1) is naturally described by al-functions as \[Ma\],

(1.2)
$$\frac{^2}{v_1^{(g)}v_2^{(g)}}\mathrm{log}\frac{\mathrm{al}_r}{\mathrm{al}_s}=\frac{1}{(b_rb_s)}\left(f^{}(b_s)\left(\frac{\mathrm{al}_r}{\mathrm{al}_s}\right)^2+f^{}(b_r)\left(\frac{\mathrm{al}_s}{\mathrm{al}_r}\right)^2\right).$$

Here $`f^{}(x):=df(x)/dx`$ and $`v^{(g)}`$’s are defined in (2.4). ((1.2) was obtained in the previous article \[Ma\] by more direct computations and will be shown as Corollary 3.3 in this article). We call (1.2) Weierstrass relation in this article.

In this article, we will introduce an “al” function over the subvariety in the Jacobian, $`\mathrm{al}_r^{(m)}:=\sqrt{F_m(b_r)}`$ and $`F_m(z):=(x_1z)\mathrm{}(x_mz)`$ for a point $`((x_1,y_1),`$ $`\mathrm{},`$ $`(x_m,y_m))`$ in the symmetric product of the $`m`$ curves $`\mathrm{Sym}^mC_g`$ $`(m=1,\mathrm{},g1)`$. In \[Mu\], Mumford dealt with $`F_m`$ function (he denoted it by $`U`$) for $`1m<g`$ and studied the properties. Further Abenda and Fedorov also studied some properties of the $`\mathrm{al}_r^{(m)}`$ and $`F_m`$ functions in \[AF\] though they did not mention about Weierstrass’s paper nor the relation (1.2). We will consider a variant of the Weierstrass relation (1.2) to $`\mathrm{al}_r^{(m)}`$ over subvariety in non-degenerated and degenerated hyperelliptic Jacobian.

As in our main theorem 3.1, even on the subvarieties, we have a similar relation to (1.1),

(1.3)
$$\begin{array}{cc}\hfill \frac{}{v_r^{(m)}}\frac{}{v_s^{(m)}}\mathrm{log}\frac{\mathrm{al}_r^{(m)}}{\mathrm{al}_s^{(m)}}& =\frac{1}{(b_rb_s)}\left(\frac{f^{}(b_r)}{(Q_m^{(2)}(b_r))^2}\left(\frac{\mathrm{al}_s^{(m)}}{\mathrm{al}_r^{(m)}}\right)^2+\frac{f^{}(b_s)}{(Q_m^{(2)}(b_s))^2}\left(\frac{\mathrm{al}_r^{(m)}}{\mathrm{al}_s^{(m)}}\right)^2\right)\hfill \\ & +\text{reminder terms}.\hfill \end{array}$$

Here $`Q_m^{(2)}`$ is defined in (2.2). We regard (1.3) or (3.1) as a subvariety version of the Weierstrass relation (1.2). In fact, (1.3) contains the same form as (1.1) up to the factors $`(Q_m^{(2)}(b_t))^2`$ ($`t=r,s`$) and the reminder terms. Thus (1.3) or (3.1) should be regarded as an extension of the sine-Gordon equation (1.2) over the Jacobian to that over the subvariety of the Jacobian.

Further a certain degenerate curve, the remainders in (1.3) vanishes. Then we have a relations over subvarieties in the Jacobian, which is formally the same as the Weierstrass relations (1.2) up to the factors $`(Q_m^{(2)}(b_t))^2`$ ($`t=r,s`$), which means that we can find solutions of sine-Gordon equation over subvarieties in hyperelliptic Jacobian. We expect that our results shed a light on the new field of a relation between “integrability” and a subvariety in the Jacobian, which was brought off by \[V1, V2, AF\].

The author is grateful to the referee for directing his attensions to the references \[AF\] and \[V2\].

## 2. Differentials of a Hyperelliptic Curve

In this section, we will give our conventions of hyperelliptic functions of a hyperelliptic curve $`C_g`$ of genus $`g`$ $`(g>0)`$ given by an affine equation,

(2.1)
$$\begin{array}{cc}\hfill y^2& =f(x)=(xb_1)(xb_2)\mathrm{}(xb_{2g})(xb_{2g+1})\hfill \\ & =Q(x)P(x),\hfill \end{array}$$

where $`b_j`$’s are complex numbers. Here we use the expressions $`Q(x):=Q_m^{(1)}(x)Q_m^{(2)}(x)`$,

(2.2)
$$\begin{array}{cc}\hfill Q_m^{(1)}(x)& :=(xa_1)(xa_2)\mathrm{}(xa_m),\hfill \\ \hfill Q_m^{(2)}(x)& :=(xa_{m+1})(xa_{m+2})\mathrm{}(xa_g),\hfill \\ \hfill P(x)& :=(xc_1)(xc_2)\mathrm{}(xc_g)(xc),\hfill \end{array}$$

where $`a_kb_k`$, $`c_kb_{g+k}`$, $`(k=1,\mathrm{},g)`$ $`cb_{2g+1}`$.

###### Definition 2.1.

\[Ba1, Ba2\]

1. For a point $`(x_i,y_i)C_g`$, we define the following quantities.
2. The unnormalized differentials of the first kind are defined by,

(2.3)
$$dv_k^{(g,i)}:=\frac{Q(x_i)dx_i}{2(x_ia_k)Q^{}(a_k)y_i},(k=1,\mathrm{},g)$$
3. The Abel map for $`g`$-th symmetric product of the curve $`C_g`$ is defined by,

$$v^{(g)}(v_1^{(g)},\mathrm{},v_g^{(g)}):\mathrm{Sym}^g(C_g)^g,$$

(2.4)
$$\left(v_k^{(g)}((x_1,y_1),\mathrm{},(x_g,y_g)):=\underset{i=1}{\overset{g}{}}_{\mathrm{}}^{(x_i,y_i)}𝑑v_k^{(g,i)}\right).$$
4. For $`v^{(g)}^g`$, we define the subspace,

(2.5)
$$\mathrm{\Xi }_m:=v^{(g)}(\mathrm{Sym}^m(C_g)\times (a_{m+1},0)\times \mathrm{}\times (a_g,0))/𝚲.$$

Here $``$ is a complex field and $`𝚲`$ is a $`g`$-dimensional lattice generated by the related periods or the hyperelliptic integrals of the first kind.

The Jacobi variety $`𝒥_g`$ are defined as complex torus as $`𝒥_g:=\mathrm{\Xi }_g`$. As $`\mathrm{\Xi }_m`$ ($`m<g`$) is embedded in $`𝒥_g`$ whose complex dimension as subvariety is $`m`$, the differential forms $`(dv_k^{(g)})_{k=1,\mathrm{},g}`$ are not linearly independent. We select linearly independent bases such as $`v_k^{(m)}:=v_k^{(g)}((x_1,y_1),\mathrm{},(x_m,y_m),(a_{m+1},0),\mathrm{},(a_g,0))`$, $`(k=1,\mathrm{},m)`$ at $`\mathrm{\Xi }_m`$.

$$\mathrm{\Xi }_0\mathrm{\Xi }_1\mathrm{\Xi }_2\mathrm{}\mathrm{\Xi }_{g1}\mathrm{\Xi }_gJ_g$$

For $`(x_1,\mathrm{},x_m)\mathrm{Sym}^m(C_g)`$, we introduce

(2.6)
$$F_m(x):=(xx_1)\mathrm{}(xx_m),$$

and in terms of $`F_m(x)`$, a hyperelliptic $`al`$-function over $`(v^{(m)})\mathrm{\Xi }_m`$, \[Ba2 p.340, W\],

(2.7)
$$\mathrm{al}_r^{(m)}(v^{(m)})=\sqrt{F_m(b_r)}.$$

Further we introduce $`m\times m`$-matrices,

$$_m:=\left(\begin{array}{cccc}\frac{1}{x_1a_1}& \frac{1}{x_2a_1}& \mathrm{}& \frac{1}{x_ma_1}\\ \frac{1}{x_1a_2}& \frac{1}{x_2a_2}& \mathrm{}& \frac{1}{x_ma_2}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{1}{x_1a_m}& \frac{1}{x_2a_m}& \mathrm{}& \frac{1}{x_ma_m}\end{array}\right),$$

$$𝒬_m=\left(\begin{array}{cccc}\sqrt{\frac{Q(x_1)}{P(x_1)}}& & & \\ & \sqrt{\frac{Q(x_2)}{P(x_2)}}& & \\ & & \mathrm{}& \\ & & & \sqrt{\frac{Q(x_m)}{P(x_m)}}\end{array}\right),𝒜_m=\left(\begin{array}{cccc}Q^{}(a_1)& & & \\ & Q^{}(a_2)& & \\ & & \mathrm{}& \\ & & & Q^{}(a_m)\end{array}\right).$$

###### Lemma 2.2.

1. $$det_m=\frac{(1)^{m(m1)/2}P(x_1,\mathrm{},x_m)P(a_1,\mathrm{},a_m)}{_{k,l}(x_ka_l)},$$

where

$$P(z_1,\mathrm{},z_m):=\underset{i<j}{}(z_iz_j).$$
2. $$_m^1=\left[\left(\frac{F_m(a_j)Q_m^{(1)}(x_i)}{F_m^{}(x_i)Q_m^{(1)}(a_j)(a_jx_i)}\right)_{i,j}\right],$$

where $`F_m^{}(x):=dF_m(x)/dx`$ and $`Q_m^{(1)}(x)=dQ_m^{(1)}(x)/dx`$.
3. (2.8)
$$(𝒬)^1𝒜=\left[\left(\frac{2y_iF_m(a_j)}{F_m^{}(x_i)Q_m^{(2)}(x_i)(a_jx_i)}\right)_{i,j}\right].$$

###### Proof.

(1) is a well-known result \[T\]. Since the zero and singularity in the left hand side give the right hand side as

$$C\mathrm{P}(x_1,\mathrm{},x_m)\mathrm{P}(a_1,\mathrm{},a_m)/\underset{k,l}{}(x_ka_l),$$

for a certain constant $`C`$. In order to determine $`C`$, we multiply $`_{k,l}(x_ka_l)`$ both sides and let $`x_1=a_1`$, $`x_2=a_2`$, $`\mathrm{}`$, and $`x_m=a_m`$. Then $`C`$ is determined as above. (2) is obtained by the Laplace formula using the minor determinant for the inverse matrix. On (3) we note that $`Q_m^{(1)}Q_m^{(2)}=Q(x)`$ in (2.2) and thus $`Q_m^{(1)}(x)\sqrt{P(x)/Q(x)}=y/Q_m^{(2)}`$. Then we obtain (3). ∎

###### Corollary 2.3.

Let $`_{v_i}^{(r)}:=/v_i^{(r)}`$, and $`_{x_i}:=/x_i`$.

(2.9)
$$\left(\begin{array}{c}_{v_1}\\ _{v_2}\\ \mathrm{}\\ _{v_m}\end{array}\right)=2(𝒬_m)^1𝒜_m\left(\begin{array}{c}_{x_1}\\ _{x_2}\\ \mathrm{}\\ _{x_m}\end{array}\right).$$

## 3. Weierstrass relation on $`\mathrm{\Xi }_m`$

The hyperelliptic solution of the sine-Gordon equation over $`𝒥_g`$ related to ramified points $`(a_1,0)`$ and $`(a_2,0)`$ is obtained as (1.1) by Mumford \[Mu\], whose expression in an old fashion is the Weierstrass relation (1.2). Let us consider an extension of the Weierstrass relation (1.2) over $`\mathrm{\Xi }_m`$ as our main theorem. We will give the theorem as follows.

###### Theorem 3.1.

$`\mathrm{al}_r^{(m)}`$ and $`\mathrm{al}_s^{(m)}`$ $`(r,s\{1,2,\mathrm{},m\})`$ over $`\mathrm{\Xi }_m`$ in (2.5) obey the relation,

(3.1)
$$\begin{array}{cc}\hfill \frac{}{v_r^{(m)}}\frac{}{v_s^{(m)}}& \mathrm{log}\frac{\mathrm{al}_r^{(m)}(v^{(m)})}{\mathrm{al}_s^{(m)}(v^{(m)})}\hfill \\ & =\frac{1}{(a_ra_s)}\left(\frac{f^{}(a_r)}{(Q_m^{(2)}(a_r))^2}\left(\frac{\mathrm{al}_s^{(m)}(v^{(m)})}{\mathrm{al}_r^{(m)}(v^{(m)})}\right)^2+\frac{f^{}(a_s)}{(Q_m^{(2)}(a_s))^2}\left(\frac{\mathrm{al}_r^{(m)}(v^{(m)})}{\mathrm{al}_s^{(m)}(v^{(m)})}\right)^2\right)\hfill \\ & +\frac{f^{}(a_{m+1})(\mathrm{al}_r^{(m)}(v^{(m)}))^2(\mathrm{al}_s^{(m)}(v^{(m)}))^2(a_ra_s)}{(a_{m+1}a_r)(a_{m+1}a_s)(\mathrm{al}_{m+1}^{(m)}(v^{(m)}))^4(Q_m^{(2)}(a_{m+1}))^2}\hfill \\ & +\mathrm{}\mathrm{}\hfill \\ & +\frac{f^{}(a_g)(\mathrm{al}_r^{(m)}(v^{(m)}))^2(\mathrm{al}_s^{(m)}(v^{(m)}))^2(a_ra_s)}{(a_ga_r)(a_ga_s)(\mathrm{al}_g^{(m)}(v^{(m)}))^4(Q_m^{(2)}(a_g)^{})^2}.\hfill \end{array}$$

###### Proof.

From (2.7), we will consider the following formula instead of (3.1) without loss of generality,

(3.2)
$$\begin{array}{cc}\hfill \frac{}{v_1^{(m)}}\frac{}{v_2^{(m)}}\mathrm{log}\frac{F_m(a_1)}{F_m(a_2)}& =2\frac{F_m(a_1)F_m(a_2)}{(a_1a_2)}(\frac{f^{}(a_1)}{F_m(a_1)^2(Q_m^{(2)}(a_1))^2}+\frac{f^{}(a_2)}{F_m(a_2)^2(Q_m^{(2)}(a_1))^2}\hfill \\ & +\frac{f^{}(a_{m+1})(a_1a_2)^2}{(a_{m+1}a_1)(a_{m+1}a_2)F_m(a_{m+1})^2(Q_m^{(2)}(a_{m+1}))^2}\hfill \\ & +\mathrm{}\hfill \\ & +\frac{f^{}(a_g)(a_1a_2)^2}{(a_ga_1)(a_ga_2)F_m(a_g)^2(Q_m^{(2)}(a_g))^2}).\hfill \end{array}$$

Before we start the proof, we will comment on our strategy, which is essentially the same as \[Ba3\]. First we translate the words of the Jacobian into those of the curves; we rewrite the differentials $`v_{(r)}^{(m)}`$’s in terms of the differentials over curves as in (3.3). We count the residue of an integration and use a combinatorial trick as in Lemma 3.2. Then we will obtain (3.2).

From (2.8) and (2.9), the derivative $`v`$’s over $`\mathrm{\Xi }_m`$ in (2.5) are expressed by the affine coordinate $`x_i`$’s,

(3.3)
$$\frac{}{v_i^{(m)}}=F_m(a_i)Q_m^{(2)}(a_i)\underset{j=1}{\overset{m}{}}\frac{2y_j}{F_m^{}(x_j)Q_m^{(2)}(x_j)(x_ja_i)}\frac{}{x_j}.$$

The right hand side of (3.2) becomes,

$$\begin{array}{cc}& \frac{^2}{v_1v_2}\mathrm{log}\frac{F_m(a_1)}{F_m(a_2)}=F_m(a_1)Q_m^{(2)}(a_1)\hfill \\ & \underset{j=1i=1}{\overset{m}{}}\frac{2y_j}{(x_ia_1)F_m^{}(x_j)Q_m^{(2)}(x_j)}\frac{}{x_j}\frac{2y_iF_m(a_2)Q_m^{(2)}(a_1)}{F_m^{}(x_i)Q_m^{(2)}(x_i)(x_ia_2)}\frac{(a_1a_2)}{(x_ia_1)(x_ia_2)}.\hfill \end{array}$$

The right hand side is

$$F_m(a_1)F_m(a_2)(_{i=1}^m\frac{1}{F_m^{}(x_i)}[\frac{}{x}\left(\frac{f(x)(a_2a_1)}{(xa_1)^2(xa_2)^2(Q_m^{(2)}(x))^2F_m^{}(x)}\right)]_{x=x_i}$$

$$\underset{k,l,kl}{}\frac{2y_ky_l(a_2a_1)}{F_m^{}(x_k)F_m^{}(x_l)(x_la_1)(x_la_2)Q_m^{(2)}(x_l)(x_ka_1)(x_ka_2)Q_m^{(2)}(x_k)(x_lx_k)}).$$

Then the proof of Theorem 3.1 is completely done due to next lemma. ∎

###### Lemma 3.2.

1. $$\begin{array}{cc}& \underset{i=1}{\overset{m}{}}\frac{1}{F_m^{}(x_i)}\left[\frac{}{x}\left(\frac{f(x)}{(xa_1)^2(xa_2)^2(Q_m^{(2)}(x)^2F_m^{}(x)}\right)\right]_{x=x_i}\hfill \\ & =\frac{2}{(a_1a_2)^2}(\frac{f^{}(a_1)}{F_m(a_1)^2(Q_m^{(2)}(a_1))^2}+\frac{f^{}(a_2)}{F_m(a_2)^2(Q_m^{(2)}(a_1))^2}\hfill \\ & +\frac{f^{}(a_{m+1})(a_1a_2)^2}{(a_{m+1}a_1)(a_{m+1}a_2)F_m(a_{m+1})^2(Q_m^{(2)}(a_{m+1}))^2}\hfill \\ & +\mathrm{}\hfill \\ & +\frac{f^{}(a_g)(a_1a_2)^2}{(a_ga_1)(a_ga_2)F_m(a_g)^2(Q_m^{(2)}(a_g))^2}).\hfill \end{array}$$
2. $$\underset{k,l,kl}{}\frac{2y_ky_l(a_2a_1)}{F_m^{}(x_k)F_m^{}(x_l)(x_la_1)(x_la_2)Q_m^{(2)}(x_l)(x_ka_1)(x_ka_2)Q_m^{(2)}(x_k)(x_lx_k)})=0.$$

###### Proof.

: (1) will be proved by the following residual computations: Let $`C_g^o`$ be the boundary of a polygon representation $`C_g^o`$ of $`C_g`$,

(3.4)
$$_{C_g^o}\frac{f(x)}{(xa_1)^2(xa_2)^2F_m(x)^2(Q_m^{(2)}(x))^2}𝑑x=0.$$

The divisor of the integrand of (3.4) is

$$\underset{i=1}{\overset{2g+1}{}}(b_i,0)4\underset{i=1,2,m+1,m+2,\mathrm{},g}{}(a_i,0)2\underset{i=1}{\overset{m}{}}(x_i,y_i)2\underset{i=1}{\overset{m}{}}(x_i,y_i)+3\mathrm{}$$

We check these poles: First we consider the contribution around $`\mathrm{}`$ point.

$$\begin{array}{cc}\hfill \mathrm{res}_{(x_k,\pm y_k)}& \frac{f(x)}{(xa_1)^2(xa_2)^2F_m(x)^2(Q_m^{(2)}(x))^2}dx\hfill \\ & =\frac{1}{F_m^{}(x_k)}\left[\frac{}{x}\left(\frac{f(x)}{(xa_1)^2(xa_2)^2(Q_m^{(2)}(x))^2F_m^{}(x)}\right)\right]_{x=x_k}.\hfill \end{array}$$

At the point $`(a_1,0)`$, noting that the local parameter $`t`$ is given by $`t=\sqrt{(xa_1)}`$ there, we have

$$\mathrm{res}_{(a_1,0)}\frac{f(x)}{(xa_1)^2(xa_2)^2F_m(x)^2(Q_m^{(2)}(x))^2}dx=\frac{2f^{}(a_1)}{(a_1a_2)^2F_m(a_1)^2(Q_m^{(2)}(a_1))^2}.$$

The residue at $`(a_2,0)`$ is similarly obtained. For the points $`(a_k,0)`$ $`(gk>m)`$, we have

$$\mathrm{res}_{(a_k,0)}\frac{f(x)}{(xa_1)^2(xa_2)^2F_m(x)^2(Q_m^{(2)}(x))^2}dx=\frac{2f^{}(a_k)}{(a_ka_1)^2(a_ka_2)^2F_m(a_2)^2(Q_m^{(2)}(a_k))^2}.$$

By arranging them, we obtain (1). (2) is obvious. ∎

As a corollary, we have Weierstrass relation (1.2) which was proved in \[Ma\]:

###### Corollary 3.3.

For $`m=g`$ case, we have the Weierstrass relation for a general curve $`C_g`$,

(3.5)
$$\frac{}{v_r^{(g)}}\frac{}{v_s^{(g)}}\mathrm{log}\frac{\mathrm{al}_r^{(g)}}{\mathrm{al}_s^{(g)}}=\frac{1}{(a_ra_s)}\left(f^{}(a_r)\left(\frac{\mathrm{al}_s^{(m)}}{\mathrm{al}_r^{(m)}}\right)^2+f^{}(a_s)\left(\frac{\mathrm{al}_r^{(m)}}{\mathrm{al}_s^{(m)}}\right)^2\right).$$

Now we will give our finial proposition as corollary.

###### Corollary 3.4.

For a curve satisfying the relations $`c_j=a_j`$ for $`j=m+1,\mathrm{},g`$, $`\mathrm{al}_r^{(m)}`$ and $`\mathrm{al}_s^{(m)}`$ $`(r,s\{1,2,\mathrm{},m\})`$ over $`\mathrm{\Xi }_m`$ in (2.5) obey the relation,

(3.6)
$$\frac{}{v_r^{(m)}}\frac{}{v_s^{(m)}}\mathrm{log}\frac{\mathrm{al}_r^{(m)}}{\mathrm{al}_s^{(m)}}=\frac{1}{(a_ra_s)}\left(\frac{f^{}(a_r)}{(Q_m^{(2)}(a_r))^2}\left(\frac{\mathrm{al}_s^{(m)}}{\mathrm{al}_r^{(m)}}\right)^2+\frac{f^{}(a_s)}{(Q_m^{(2)}(a_s))^2}\left(\frac{\mathrm{al}_r^{(m)}}{\mathrm{al}_s^{(m)}}\right)^2\right).$$

###### Proof.

Since the condition $`c_j=a_j`$ for $`j=m+1,\mathrm{},g`$ means $`f^{}(a_j)=0`$ for $`j=m+1,\mathrm{},g`$, Theorem 3.1 reduces to this one. ∎

Under the same assumption of Corollary 3.4, letting $`A={\displaystyle \frac{2\sqrt{f^{}(a_r)f^{}(a_s)}}{(a_ra_s)Q_m(a_r)Q_m(a_s)}},`$ and

$$\varphi _m^{(r,s)}(u):=\frac{1}{\sqrt{1}}\mathrm{log}\sqrt{\frac{f^{}(a_r)}{f^{}(a_s}}\frac{Q_m(a_r)}{Q_m(a_s)}\frac{F_m(a_r)}{F_m(a_s)},$$

defined over $`\mathrm{\Xi }_m`$, $`\varphi _m^{(r,s)}`$ obeys the sin-Gordon equation,

(3.7)
$$\frac{}{v_r^{(m)}}\frac{}{v_s^{(m)}}\varphi _m^{(r,s)}=A\mathrm{cos}(\varphi _m^{(r,s)}).$$


format: text

text:
# A first-principles approach to electrical transport in atomic-scale nanostructures

## I Introduction

Molecular- and atomic-scale electronic devices are attracting an ever increasing interest due to the impact they are expected to make in the world of Nanotechnology. The number of experimental and theoretical works in this particular area of research, generically known as molecular electronics, is growing exponentially. The design of devices at the molecular and even atomic scale poses new challenges which require new theoretical and experimental techniques to be developed. Scanning Tunneling Microscopy (STM) is probably the pioneer of the experimental techniques in this research area. It can be used not only in the tunneling regime to image adsorbates, but also in the contact regime to build few-atom nanoscopic contacts. STM can also be used to investigate the electrical properties of nanotubes and DNA molecules with one or both of their ends attached to a suitable electrode. In addition to STM, mechanically controllable break junctions have also revealed themselves as powerful tools to study electrical transport in metallic nanobridges or individual molecules.

The basics to calculate the zero-bias conductance $`G`$ of a nanoscale contact had been established by Landauer’s in his pioneering work long before these systems were commonplace. In Landauer’s formalism $`G`$ is simply given by the quantum mechanical transmission of the electrons around the Fermi energy. The value of this transmission is essentially determined by the region where the number of channels available for conduction is the smallest. In molecular or atomic-scale nanocontacts the region of relevance is the molecule and/or the few atoms forming the nanoscopic bridge between electrodes. The transmission is thus strongly dependent on the particular molecule, the detailed atomic arrangement of the electrodes in the contact region, and the chemical nature of them. Knowing the atomic arrangement of the electrodes or the way the molecule binds to the electrodes is, however, a major problem in itself. Furthermore, even if these important details were known, implementing Landauer’s formalism still requires to know the electronic structure and this is a formidable task as well.

Calculations based on tight-binding or semi-empirical models have been, and still are, very popular since these models capture the atomic-scale character in some detail and are easy to implement. However, they do not allow for structural relaxations to be performed. Most importantly, these simple models, in general, do not yield correct values for the local electronic charges. In other words, the chemical potential is not uniform across the entire system in equilibrium. While imposing local charge neutrality is a straightfoward improvement on these models for metallic nanoconstrictions, there does not exist any simple modification in the case of more complex systems like metal-molecule-metal heterostructures. A way around this problem is to perform self-consistent first-principles calculations which, at least at a mean-field level, guarantee the uniformity of the chemical potential. However, most numerical implementations commonly used to carry out ab-initio electronic calculations are either restricted to finite systems, such as the Gaussian code, or require the infinite system to be periodic such as the SIESTA code. None of these methods are suitable to address the systems here studied which are both infinite and non-periodic. Finally, a perhaps more serious difficulty is the intrinsic non-equilibrium character of electrical transport.

In recent years several proposals have appeared to tackle this problem. Most of all are based upon density functional (DF) theory. In addition to the well-known virtues of the DF theory, it presents the additional advantange that Landauer’s theoretical framework does not need to be modified since DF is still a single-particle description of the many-body problem. In the pioneering works of Lang and co-workers, Tsukada and co-workers, and Guo and co-workers the electrodes were described within the jellium approximation. Jellium models are still used today; for they are convenient in one way: They provide featureless contacts which represent generic situations. As mentioned below this is one major feature of our model too. However the jellium model presents serious drawbacks. How can the differences observed experimentally between, for instance, Al and Au electrodes, be taken into account by means of a jellium?. Furthermore, this approach is not satisfactory when one is trying to describe, for instance, STM experiments where the detailed atomic structure of the tip determines, to a large extent, whether or not the STM can resolve the topography or molecular structure of the adsorbate. Recent approaches, which essentially differ only on the numerical implementation, intend to incorporate the atomic structure of the electrodes in the DF calculation. It is pertinent noting here that in most of these studies a periodic structure beyond a given point within the leads is assumed. Efficient transfer matrix techniques make this reasonable assumption tractable, but it forces one to consider a very specific type of lead (typically a finite-section wire or an infinite surface). Using a jellium model for the electrodes is harmless but it lacks the minimum atomic detail which is crucial to describe, e.g., the contact with molecules. However, employing well-defined specific electrodes is not desirable either since their own electronic structure can interfere with the interpretation of the results. For instance, the appearance of gaps in the conductance close to the Fermi energy for perfectly conducting systems such as Au chains can only be considered an artifact due to the unphysical electrode model. Furthermore, actual nanocontacts are not expected to have high symmetry. As explained below these difficulties are circumvented in our method.

Recently, we have presented an alternative to the ab-initio methods above mentioned. Close in spirit to those presented in Refs. , the main differences and advantages with respect to them are the use of the standard Gaussian98 code to carry out the DF calculation of the relevant transport region and the description of the electrodes bulk by means of appropriate Bethe lattices. The Gaussian98 code provides a versatile method to perfom first-principles calculations of clusters, incorporating the major advancements in the field in terms of functionals, basis sets, pseudopotentials, etc.. On the other hand, the Bethe lattices are two-fold convenient: (i) They reproduce the essential features of the bulk density of states and (ii) their directional self-energies can be easily calculated (see Appendix). In Ref.  these ideas were applied to investigate electrical transport of a C60 molecule in between Al electrodes. Here we describe in detail an improvement to our previous approach which can be summarized in that we now incorporate self-consistently the semi-infinite electrodes into the conductance calculation within the Gaussian98 code. This requires working with Green’s functions from the very start. The method has some resemblance with the cluster Bethe lattice method developed to investigate the physical properties (electronic structure, phonons, etc.) of disordered systems. We illustrate the possibilities of our method by investigating electronic transport in Al and Au nanocontacts. These systems have been the subject of extensive studies in the past by means of tight-binding methods and, more recently, by means of ab-initio methods. We show that, even in these systems where charge transfer is apparently unessential, tight-binding methods may fail to provide a correct quantitative picture. It turns out that the conductance in the case of Al nanocontacts is strongly dependent on the detailed atomic structure. In the case of Au nanocontacts, on the contrary, the results are more universal as confirmed by experiments. Next, we choose a system with somewhat appealing features: Carbon-atom chains. As shown in Ref.  these chains, when contacted by Al electrodes, exhibit a conductance that oscillates with the number of atoms in the chain. Here we address this problem by taking proper account of the binding to the electrodes and investigate how the results depend on the type of electrode (Al or Au).

The rest of the paper is organized as follows. In Section II we discuss the main characteristics of our method. Section III is devoted to an extensive discussion of the results. Finally, we end the paper by summarising the main features of the method and the most remarkable results (Section IV).

## II The Gaussian Embedded Cluster Method (GECM)

In previous work we have presented a method to study transport in atomic-scale and molecular devices which is based on standard quantum chemistry calculations with the Gaussian98 code. This scheme, which has been recently adopted by other groups (see, e.g., Refs. ), is taken here a step further. A DF calculation of the region that includes the molecule or set of atoms forming the contact between electrodes and a significant part of the electrodes is performed (see Fig. 1). As far as transport is concerned, the hamiltonian (or Fock matrix $`\widehat{F}`$) of this central cluster or supermolecule contains the relevant information since it embraces the region with the smallest number of channels for conduction. However, according to the usual theoretical transport schemes, its associated Green’s functions are unsuitable for the evaluation of the current (note that they simply have poles). The retarded(advanced) Green’s functions associated with $`\widehat{F}`$ needs to be extended to include the rest of the semi-infinite electrodes:

$$\widehat{G}^{r(a)}(ϵ)=(ϵ\widehat{I}\widehat{S}\widehat{F}\pm i\delta )^1[ϵ\widehat{I}\widehat{S}\widehat{F}\widehat{\mathrm{\Sigma }}^{r(a)}(ϵ)]^1.$$
(1)
In this expression

$$\widehat{\mathrm{\Sigma }}^{r(a)}(ϵ)=\widehat{\mathrm{\Sigma }}_R^{r(a)}(ϵ)+\widehat{\mathrm{\Sigma }}_L^{r(a)}(ϵ),$$
(2)
where $`\widehat{\mathrm{\Sigma }}_R`$($`\widehat{\mathrm{\Sigma }}_L`$) denotes a self-energy matrix that accounts for the part of the right(left) semi-infinite electrode which has not been included in the DF calculation. $`\widehat{S}`$ is the overlap matrix and $`\widehat{I}`$ is the unity matrix. The added self-energy matrices can only be explicitly calculated in ideal situations, which, in principle, limits the desired applicability of the above procedure. For instance, in Refs.  the authors consider finite-section wires as electrodes. As a result of this choice gaps appear in the conductance of otherwise perfectly conducting central clusters.

In order to overcome this type of problem, we choose to describe the bulk electrode with a Bethe lattice tight-binding model with the appropriate coordination numbers and parameters (see Appendix). The advantage of choosing a Bethe lattice resides in that it reproduces fairly well the bulk density of states of any metallic electrode, avoiding this way the appearance of spurious results. In addition to this the self-energy matrices that appear in Eq. 1 can be calculated iteratively in a simple way (see Appendix for more details). For each atom of the outer planes of the cluster, we choose to add a branch of the Cayley tree in the direction of any missing bulk atom (including those missing in the same plane). In Fig. 1 the directions in which branches are added are indicated by smaller atoms which represent the first atom of the branch in that direction. Assuming that the most important structural details of the electrode are included in the central cluster, the Bethe lattices should have no other relevance than that of introducing a featureless reservoir.

In our present approach the self-consistent process does not stop once the finite central cluster has been solved. Instead, we reformulate the Gaussian98 code to proceed with the self-consistency of the now infinite system. More specifically, once self-consistency for the finite cluster has been almost attained, we calculate the Green’s function as explained above. Then, the density matrix is obtained from the Green’s function according to

$$\widehat{n}=\frac{1}{\pi }_{\mathrm{}}^{ϵ_F}\mathrm{Im}\left[\widehat{G}^r(ϵ)\right]dϵ,$$
(3)
where $`ϵ_F`$ is the Fermi level fixed by the condition of overall charge neutrality in the cluster. The integral in Eq. (3) is calculated along a contour in the complex plane as explained in Refs.  with an efficient automatic numerical integration scheme of Pérez-Jordá et al.. The density matrix is now used to recalculate the matrix elements of the Fock’s operator and the process is repeated until self-consistency is achieved. We note that in this method the standard eigenvalue problem, inherent to the Gaussian98 code, is replaced by the calculation of Green’s functions. In the end the Green’s functions describe an infinite system in a more consistent way than in the method discussed previously by us since it effectively removes finite-size effects in the self-consistency. It is interesting to note that the applicability of this approach, which we hereafter name the Gaussian Embedded Cluster Method (GECM), goes well beyond the present study. In fact it could be a powerful tool whenever an infinite media has to be described (for instance adsorption of molecules on solid or liquid surfaces).

The conductance can now be simply calculated through the expression

$$G=\frac{2e^2}{h}\mathrm{Tr}[\widehat{T}],$$
(4)
where Tr denotes the trace over all the orbitals of the cluster and $`\widehat{T}`$ is the transmission matrix which, in turn, is given by

$$\widehat{T}=\widehat{\mathrm{\Gamma }}_L\widehat{G}^r\widehat{\mathrm{\Gamma }}_R\widehat{G}^a,$$
(5)
where the matrices $`\widehat{\mathrm{\Gamma }}_R`$ and $`\widehat{\mathrm{\Gamma }}_L`$ are given by $`i(\widehat{\mathrm{\Sigma }}_R^r\widehat{\mathrm{\Sigma }}_R^a)`$ and $`i(\widehat{\mathrm{\Sigma }}_L^r\widehat{\mathrm{\Sigma }}_L^a)`$, respectively. In order to single out the contribution of individual channels to the current one can diagonalize the transmission matrix. It turns out (see below) that only few channels give a non-negligible contribution to the current. The symmetry of each channel was identified by looking at its weight on the atomic orbitals of the central atom of the constriction or the C chain.

Finally, it is worth pointing out that if the current in the finite-bias regime has to be obtained, one should simply integrate in energy the above expression with appropriate Fermi distributions functions. Note, however, that a new definition of the density matrix generalized to non-equilibrium needs to be used in the calculation of the Green’s functions. Apart from this, the Landauer-type expression (4) remains valid as long as one does not give up the single-particle description. In this work we are concerned with basic and still open aspects of transport in the systems studied and we will focus on the linear regime.

## III Results

For all the DF calculations we have used the Becke’s three-parameter hybrid functional using the Lee, Yang and Parr correlation functional (B3LYP) together with the semilocal shape consistent pseudopotential (SCPP) and basis sets of Christiansen et al. . We have selected this combination of exchange-correlation functional and pseudopotential for two reasons: First, B3LYP is one of the most accurate and certainly the most popular among the gradient-corrected exchange-correlation functionals; second, Christiansen et al. SCPP provide accurate results for a wide variety of atoms while retaining the simplicity of a minimal basis set . Needless to say that there is no need to restrict the calculations to a minimal basis set since the Bethe lattice can be constructed for any basis set. It is only a matter of computational convenience that we have done it so. Nevertheless, in some cases, we have checked that better basis sets and the use of other exchange-correlation functionals do not modify the main conclusions of our work.

In all cases we have investigated the influence of the number of electrode atoms included in the DF calculation on the conductance. In general, the results do not vary qualitatively with the cluster size, allowing us to extract some general conclusions. However, at a quantitative level, this is not always the case, particularly for Al nanocontacts. Finally we note that, although most of the calculations were carried out taking the bulk inter-atomic distances for the electrode atoms, in some cases we investigated the effects of ab-initio relaxations.

### A Metallic nanocontacts

A complete theoretical study of electrical transport in metallic atomic contacts requires a realistic modelization of the formation process of these nanocontacts. Some structural studies using molecular dynamics for Al and ab-initio relaxations for Al and Na have been reported. This is, however, a problem beyond the scope of this work. Here we consider archetypical atomic structure models that are likely to appear in the last stages of the formation process of atomic contacts before the break-up: Single-atom contacts and atomic wires. More specifically, our first structure consists of two opposite pyramids grown in the (001) direction and “glued” by a single atom \[see Fig. 2\]. Single-atom contacts have been studied in the past with modified tight-binding models. We find to our surprise that, even in this simple case, our ab-initio results are qualitatively different from those obtained with these models, particularly for Al. Our second structure is similar to the previous one, but with a chain of three atoms instead of a single atom \[see Fig. 3\]. Finally we have studied the same chain in between two (111) surfaces with the chain placed on top of a surface atom \[see Fig. 4\]. A similar geometry has been recently studied with ab-initio techniques for Au and Al. Our results agree with what has been reported for Al using a jellium model for the electrodes, but not entirely with what has been obtained for Au.

#### 1 Aluminum

Figure 5 (top panel) shows the conductance vs. energy for a single-atom Al contact. We have considered a $`3s3p`$ basis set and bulk interatomic distances. Curves (a), (b), and (c) correspond to clusters (a), (b), and (c) in Fig. 2, respectively. In all the clusters the contact between electrodes occurs through a single atom, but the number of (001) planes explicitly included in the DF calculation for each pyramid increases from (a) to (c) (remember that Bethe lattices are always attached to the outer planes as in Fig.1).

It is impossible to know the actual atomic structure of the metallic contact in detail unless relaxation calculations are performed, but we do not expect the detailed geometry away from the neck to be important. In fact, as Fig. 5 shows, the conductance does not change significantly from (b) to (c), apart from minor changes in the fine structure. This is a clear indication that, to a good extent, the conductance is determined by the atomic structure in the narrowest region of the neck. However, from our results we see that, at least, the nine central atoms must be explicitly considered in this example. This is in contrast to the conclusions drawn in Ref.  by Cuevas et al. using a modified tight-binding model. Furthermore, the value of $`G`$ around the Fermi level is $`3`$ which is remarkably different from the value they obtained. This discrepancy is due to a combination of facts. First, the hopping parameters that reproduce bulk properties in tight-binding models are not adequate for atoms with low coordination numbers, being these typically smaller than the ab-initio ones. Second, there are non-negligible contributions from next-to-near-neighbor hoppings in Al. This is clearly seen in the bottom panel of Fig. 5 where the contributions to the total current of the main individual channels for the cluster (c) are depicted. As shown in the figure, a channel associated with 2nd-nearest-neighbors hoppings (labelled $`\sigma ^{}`$) can give a contribution of almost one conductance quantum $`G_0`$ ($`G_0=2e^2/h`$) at the Fermi level. The main contribution at $`ϵ_F`$ comes from two degenerate $`p_x,p_y`$-like channels ($`\pi _x,\pi _y`$) that account for almost two conductance quanta $`2G_0`$. (The $`z`$-axis has been chosen along the main symmetry axis of the cluster). In addition there are two channels that have $`sp_z`$ character ($`\sigma `$) with non-negligible contributions at the Fermi energy that add approximately 0.3$`G_0`$ and 0.1$`G_0`$ to the total conductance, respectively. None of this seems consistent with the tight-binding results for a similar geometry.

The conductance of the three-atom chain shows a different behavior from the one in the previous example and presents the same features for the two electrode models considered (see top panels in Figs. 6 and 7). For the first electrode model we have chosen a separation of 4.8Å between apex atoms of the pyramids. The conductance changes appreciably from (a) to (b) in Fig. 6, but it does not so from (b) to (c). (This dependence is similar to that in the single-atom contact). In the second case we have chosen a separation of 9.1Å between planes and we have performed an ab-initio relaxation of the chain atom positions (only the surface layer is included in the DF calculation where the number of atoms increases from 7 to 35). Here the conductance does not depend too much on the number of surface atoms \[see Figs. 4(a), (b), and (c)\]. In all the cases there are oscillations as a function of the energy which, as the bottom panels in Figs. 6 and 7 show, appear mostly in the $`\pi `$ channels. This is reminiscent of the behavior of the transmission in a Fabry-Perot interferometer due to scattering at the interfaces. These results are similar to the ones reported in Ref.  where a jellium model was used to describe the Al electrodes.

In the three cases considered there is no trace of possible conductance quantization. A general trend that can be observed in our results for the chain is that the onset of the transmission through the $`\pi `$ channels occurs close to the Fermi energy. This makes $`G`$ strongly dependent on small variations in the positions of the atoms in the chain and on the atomic structure of the electrodes close to the chain. This might explain why the experimental conductance histograms for Al are much more irregular than those for Au (see below) and other metals. However, there are many open questions regarding the details of the conductance steps for Al that illustrate the necessity of performing both relaxation and conductance calculations at the same time.

#### 2 Gold

We have repeated the conductance calculations for the same structures considered above, but now consisting of Au atoms (we have used a here $`5d6s6p`$ basis set). In principle, only the electron of the $`6s`$ orbital is expected to contribute to the conductance at the Fermi energy which should make the analysis of conductance simpler. As Fig. 8 shows, the conductance around the Fermi energy for the single-atom contact varies little from (a) to (c). We note again that, contrary to the tight-binding predictions, the conductance at the Fermi energy surpasses $`G_0`$ in all the curves. As concerns the contribution of the individual channels we note that the major contribution (almost a conductance quantum) has $`sp_zd_0`$ character ($`\sigma `$). Two degenerate channels ($`\pi _x,\pi _y`$) of $`p_xd_1`$ and $`p_yd_1`$ character (mainly $`d`$) give around 0.25$`G_0`$ quanta each. This should be expected since the $`d`$ orbitals contribute significantly to the density of states in bulk atoms and the number of near neighbors (eight) of the central atom in this cluster is almost the bulk coordination number of an fcc structure (twelve). Of the two channels that give a significant contribution at rather high energies (above 2.0 eV) one has $`sp_zd_0`$ symmetry ($`\sigma `$) and the other corresponds to 2nd-nearest-neighbors.

It is interesting to compare the single-atom contact results for Al and Au. As remarked above, whereas in the case of Al there was a very important contribution to the current at the Fermi level coming from 2nd-nearest-neighbors hoppings, in Au this was only appreciable at high energies. This cannot be understood in terms of the respective atomic radii which are very similar (1.43 Åand 1.45 Åfor Al and Au, respectively), but rather it is due to the character of the wavefunctions at the Fermi level in each case. Namely, while in Al the density of states at $`E_F`$ mainly comes from $`p`$ orbitals that are rather extended, in Au the wavefunction at that energy mainly has an $`s`$ character and is more localised.

On the other hand, the conductance around the Fermi energy for the three-atom Au chain shows clearly an upper limit of $`G_0`$ (see top panels in Figs. 9 and 10), and it does not change qualitatively with the cluster size. Nevertheless, the exact value at the Fermi energy is elusive, changing by as much as 20% from cluster to cluster. We have not been able to verify whether the conductance curves for larger clusters converge to a given one, but all the curves present a characteristic behavior: Above the Fermi energy the conductance is fairly constant while below oscillates and vanishes right above the $`d`$ channel contribution. The channel decomposition analysis is quite simple: A single $`sp_zd_0`$-like channel ($`\sigma `$) contributes around the Fermi energy. Nevertheless, it is still difficult to explain from our results the robustness of the quantization observed in the experiments which does not deviates from $`G_0`$ by more than a small percentage over a large range of stretching force. Large-scale structural studies along with conductance calculations are also desirable here in order to make a precise quantitative comparison with experiments. Recent ab-initio works for Au nanocontacts have partially addressed this problem. However, as already pointed out, the electrode model considered there seems to introduce serious difficulties in the interpretation of their conductance results.

### B Carbon chains

The conductance of C-atom chains attached to Al electrodes has been calculated from first principles in previous works by Lang and Avouris. In their calculations semi-infinite jellium models were used to describe the metal electrodes and a pseudopotential for the C cores. The self-consistent density functional procedure they used (see Refs.  for details) predicted an oscillatory behavior in the conductance at the Fermi level with maxima (minima) in those chains that had an odd (even) number of electrons. This oscillatory behavior differs from that corresponding to a closed shell electronic structure and an $`sp_z`$ hybridization of linear C chains. In this situation each C atom added to the chain provides two $`sp`$ orbitals that contribute to the $`\sigma `$ molecular orbitals (MOs) and two $`p`$ orbitals that contribute to the corresponding $`\pi `$ MOs. With a closed-shell electronic structure and an even number of C atoms there is a partial filled $`\pi `$ shell, and the chain could be considered to be intrinsically a conductor. On the other hand, if the chain has an odd number of C atoms, the $`\pi `$ and $`\sigma `$ shells are completely filled and we would have a semiconductor. Thus, one would expect that chains with an even number of C atoms would provide higher conductances than chains with an odd number of C atoms. However, either if we consider an open shell electronic structure for the C chain, or if the edge C atoms are bonded to a metal surface, the MOs can be filled in a different way, which may lead to an inversion of the aforementioned trend: that is, chains with an odd number of C atoms having a partially filled $`\pi `$ shell while chains with an even number of C atoms having this shell completely filled. The first scenario, i.e., an open-shell electronic structure, is actually the ground-state for isolated linear C chains, while the second one explains the findings of Lang and Avouris for the conductance of C chains attached to Al electrodes. In this second situation the details of the bonding region between the C and metal atoms is of primary importance. This involves knowing the geometry of the contact, the C-metal surface distance, and the electronic structure of the metal under interest. This, however, cannot be easily deduced by using a jellium model to describe the leads.

For this reason, we have applied the method described in the previous sections to these kind of systems. More precisely, we have calculated the conductance for C chains starting from 3 C atoms attached to Al and Au (001) fcc electrodes. The distance between C atoms was fixed to that used in Refs. , i.e., an equal spacing of 2.5 a.u. between C atoms. The description of the electrode surface was reduced to four metal atoms describing a hollow site (see Fig. 11) in the center of which the C chain was attached. The importance of the C-Al surface distance was analyzed by making two sets of calculations for each C chain. In the first set this distance was kept fixed at a reference value of 2Å, while in the second one the distance between the C chain and the metal electrode was that providing the minimum energy of the corresponding cluster. The results are commented next.

#### 1 Al electrodes

The discussion about the conductor or semiconductor character of C chains given in the preceding section can be easily extended to the conductance of these chains once they have been contacted to semi-infinite electrodes. Note that the oscillatory character of the conductance at the Fermi level with the number of C atoms $`n`$ in the chain found in these systems can be understood as a consequence of the oscillatory nature of the conductance with respect to the energy (reflected in Fig. 12). The oscillations in $`G(E)`$ come from two degenerate $`\pi `$ channels, thus explaining the maximum value of 2$`G_0`$ (the only exceptions being the C3 and C4 chains where additional $`\sigma `$ channels appear around the Fermi energy). The number of peaks is given by the number of $`\pi `$ orbitals in the chain which, in turn, depends on the number of C atoms present in the chain. The positioning of the Fermi level near the center of these peaks (partially filled $`\pi `$ shell) or in between peaks (completely filled $`\pi `$ shell) would determine the oscillatory trend with $`n`$. On the other hand, both the exact position of the conductance peaks with respect to the Fermi level and their average size are governed by the distance between the edge of the C chain and the Al electrode.

After inspection of Fig. 12, where we plot the $`G(E)`$ for the two sets of calculations mentioned at the beginning of this section, it is evident that at a distance of 2Åbetween the edge C atom and the Al metal surface we end up with the same situation found by Lang and Avouris , namely, maxima (minima) located at $`n`$ odd (even) (see Fig. 13). Nevertheless, at the equilibrium distance, where the C-Al surface distance is smaller (1Å), and the transferred charge from the metal to the C chain is larger, the situation has changed and now the maxima are located on chains with an even number of C atoms (the only exception being again the C3 chain) as Fig. 13 shows. However, one would expect that the influence of the bonding region, a local effect, would become less important as the size of the C chain increases. This is what is actually found in the conductance of larger number of atoms, such as those presented in Fig. 14, where we plot $`G(E)`$ for $`n`$=14,15 and for $`n`$=30,31 and a distance C-Al surface of 1Å. In both cases the chains with $`n`$ odd give larger values for $`G(ϵ_F)`$ than those with $`n`$ even. The same behavior is found for the same chains attached to Al electrodes at a distance of 2Å, but has not been included in the figures for simplicity. It is worth noting that, strictly speaking at zero temperature, the peaks in $`G(E)`$ and, thus, the oscillations in $`G(ϵ_F)`$ are expected to survive in the thermodynamic limit $`n\mathrm{}`$. This effect will disappear when the temperature becomes of the same magnitude as the peaks width, where $`G(ϵ_F)`$ will be the average value for the odd- and even-$`n`$ chains.

#### 2 Au electrodes

Another key point in the conductance of nanocontacts is the nature of the metal used in the electrodes which has been thoroughly discussed in previous sections. This is specially the case when the chemical bond between the molecule and the metal contact changes markedly as reflected in Fig. 15 as compared to Fig. 12. There we show $`G(E)`$ for the same systems previously analyzed for Al but with Au instead. When we move from Al to Au the bonding between the C chain and the metal surface weakens. This reflects in a sharpening of the conductance peaks and in the amount of charge transferred from the metal to the C chain (compare Figs. 16 and 13). The equilibrium distance between the edge C atom and the Au surface is also typically larger in this case: 1.5Å for all the chains. As for Al, there also appear two $`\pi `$ channels, except for the C3 and C4 chains. However, the narrowing of the conductance peaks makes critical the alignment of the Fermi level. With respect to the oscillations of $`G(ϵ_F`$), we observe the same tendency irrespective of the C-Au distance, that is: maxima located at $`n`$ odd. This is also a consequence of the weaker bond between Au and C and the lesser amount of charge transferred from the former to the latter.

## IV Concluding remarks

In summary, we have developed a methodology to self-consistently calculate ab-initio transport properties in atomic-scale systems based upon the Gaussian98 code. This opens the doors to the use of the most standard quantum chemistry tool to the study of transport through molecules, an interdisciplinary subject of increasing interest. We have chosen to study two systems that illustrate the capabilities of our approach: metallic constrictions of simple and noble metals and C chains with a reactive electrode (Al) and an “inert” (Au) electrode. In the first case, we have shown that assuming local charge neutrality, as commonly done in semi-empirical methods, may lead to qualitatively incorrect results. On the other hand, the system metal-C chain-metal illustrated how the chemistry of the contact may determine electrical transport and, therefore, the incorrectness of ascribing one or another behavior to a given molecule leaving apart the specific electrode/molecule chemistry. One should speak, instead, of electrical transport characteristics of the whole system, electrode-molecule-electrode.

###### Acknowledgements.

Part of this work was supported by the Spanish CICYT under Grants Nos. 1FD97-1358 and PB96-0085, and by the Generalitat Valenciana under Grants No. GV00-151-01 and GV00-095-2. Discussions with L. Pastor-Abia, J. M. Pérez-Jordá, J. C. Sancho, N. Agraït, P. Serena, A. Hasmy, D. Sánchez-Portal, C. Untiedt, F. Guinea, and J. J. Sáenz are acknowledged. We are also grateful to D. A. Papaconstantopoulos for useful correspondence concerning the tight-binding parameters used in band structure calculations.

## V appendix

In this appendix we discuss how selfenergies for Bethe Lattices (BL) with no symmetry can be calculated. Symmetry can be broken due to either the spatial atomic arrangement or the orbitals on the atoms that occupy each lattice site, or both. When no symmetry exists the selfenergy in an arbitrary direction cannot be obtained by rotating that for a given direction as done in Ref. . Instead, the following procedure has to be followed. The method is valid for any basis set or lattice. Let $`\tau _𝐢`$ be the $`N`$ nearest-neighbor directions of the lattice we are interested in and $`\widehat{V}_{\tau _𝐢}`$ the interatomic interaction matrix in these directions. The selfenergies associated to each direction have to be obtained from the following set of $`2N`$ coupled selfconsistent equations,

$$\widehat{\mathrm{\Sigma }}_{\tau _𝐢}=\widehat{V}_{\tau _𝐢}\left[E\widehat{\mathrm{𝟏}}\widehat{E}_0(\widehat{\mathrm{\Sigma }}_{\overline{T}}\widehat{\mathrm{\Sigma }}_{\overline{\tau _𝐢}})\right]^1\widehat{V}_{\tau _𝐢}^{}$$
(7)
$$\widehat{\mathrm{\Sigma }}_{\overline{\tau _𝐢}}=\widehat{V}_{\overline{\tau _𝐢}}\left[E\widehat{\mathrm{𝟏}}\widehat{E}_0(\widehat{\mathrm{\Sigma }}_T\widehat{\mathrm{\Sigma }}_{\tau _𝐢})\right]^1\widehat{V}_{\overline{\tau _𝐢}}^{},$$
(8)
where $`i=1,\mathrm{},N`$ and $`\overline{\tau _𝐢}=\tau _𝐢`$. $`E\widehat{\mathrm{𝟏}}`$ is the energy times the identity matrix, $`\widehat{E}_0`$ is a diagonal matrix containing the orbital levels, $`\widehat{V}_{\tau _𝐢}`$ is the interatomic interaction in the $`\tau _𝐢`$ direction, and $`\widehat{\mathrm{\Sigma }}_T`$ and $`\widehat{\mathrm{\Sigma }}_{\overline{T}}`$ are the sums of the selfenergy matrices entering through all the Cayley tree branches attached to an atom and their inverses, respectively, i.e.,

$$\widehat{\mathrm{\Sigma }}_T=\underset{i=1}{\overset{N}{}}\widehat{\mathrm{\Sigma }}_{\tau _𝐢}$$
(10)
$$\widehat{\mathrm{\Sigma }}_{\overline{T}}=\underset{i=1}{\overset{N}{}}\widehat{\mathrm{\Sigma }}_{\overline{\tau _𝐢}}.$$
(11)
This set of $`2N`$ matricial equations has to be solved iteratively. It is straightforward to check that, in cases of full symmetry, it reduces to the single equation discussed in .

The tight-binding parameters, which include only nearest-neighbors interactions, used in these calculations are given in Table I. The Table reports data not only for the metals taken as electrodes in the present work (Al and Au) but also for two additional metals (Ti and W) commonly used in experiments and/or calculations. All were obtained through fittings to the electronic bulk band structures calculated by including second or even third nearest-neighbors interactions. In the case of hexagonal-close-packed Ti we took as nearest-neighbors six out of plane and six in-plane neighbors, as actual interatomic distances differ in less than a 2%. The densities of states on bulk atoms are shown in Fig. (17) and (18). Although some of the features of the actual DOS are not reproduced (as it commonly occurs in the Bethe lattice approximation) the overall results are satisfactory. The major discrepancy is noted for Ti since in this metal the Fermi level in the crystalline case lies in a valley of the DOS. In the calculations reported in this work the $`sp`$ basis was used for Al. We have described the electronic structure of Au by means of the $`spd`$ basis (see Fig. 17 and Table I). In the case of Au and Ti, the reduced $`sd`$ basis gives an excessively narrow conduction band. Finally, we note that in the case of Ti there are no major differences between the DOS for the fcc and hcp lattices.


text:
# Laser-induced collective excitations in a two-component Fermi gas

## I Introduction

Fermionic atoms have already been cooled well below the Fermi temperature . Once the quantum regime has been achieved, the predicted BCS superfluid transition for two different hyperfine states with attractive interactions between them would be the next goal to achieve. There are several methods proposed for probing the superconducting gap . Some of them are based on the use of laser light. For instance, it has been proposed to measure the superfluid coherence and the existence of the superfluid gap by the absorption of almost on-resonance laser light coupled between the two hyperfine states  or between one of them and another atomic hyperfine state . Also the use of off-resonant (or very low intensity) laser light  has been proposed. There both hyperfine components are coupled to other excited states not involved in the pairing, and the excited state populations remain negligible.

In this paper we review and analyze the linear density response of the system to the different laser excitations in the collisionless regime where the lifetime of the quasiparticles is much longer than the period of the applied field. The motivation for studying the response is two-fold. First, in the above-mentioned probing schemes for the gap a below-gap density response (e.g. Anderson- Bogoliubov phonon) would be undesireable because the schemes rely on the absence of below-gap excitations. Therefore it is of interest to clarify in what kind of laser probing schemes density reponses can be avoided. Second, the density response itself may be the object of study and the laser excitation schemes that can produce a response are of interest.

Off-resonant laser excitation caused by two intersecting waves which form a density grating in space and time (Bragg spectroscopy) has been used to induce collective excitations in atomic Bose-Einstein condenstates . In case of near-resonant laser excitations, we transform the non-diagonal atom-light interaction Hamiltonian into a form of a density perturbation. We show that the initial Hamiltonian of the interacting two-component Fermi gas is preserved by this transformation in case of a contact interaction and therefore a standard linear density response calculation for a BCS-system can be applied. In general, an intensity grating in space and time is needed for providing momentum and energy transfer for density perturbations. We show that when the laser light is coupled between the two paired hyperfine components, phonons cannot be excited even when the light forms an intensity grating.

In section II, the Hamiltonian for a two-component Fermi gas is given and the density response is defined. The laser probing schemes considered in this article are summarized in section III. In section IV we transform the perturbation caused by the probing lasers into a density perturbation by appropiate rotation of the coupled states. The results are presented in detail in section V and summarized in section VI. The linear response calculation is presented in the Appendix.

## II The system

The Hamiltonian for the two-component Fermi gas, including only 2-body interactions between the atoms reads

$`H_0={\displaystyle 𝑑𝐫\underset{\alpha =,}{}\left[\psi _\alpha ^{}(𝐫)\left(\frac{^2}{2m}+V_\alpha (𝐫)\mu _\alpha \right)\psi _\alpha (𝐫)\right]}`$ (1)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha ,\beta =,}{}}{\displaystyle 𝑑𝐫𝑑𝐫^{}g_{\alpha \beta }(𝐫𝐫^{})\psi _\alpha ^{}(𝐫)\psi _\beta ^{}(𝐫^{})\psi _\beta (𝐫^{})\psi _\alpha (𝐫)},`$ (2)
where the two hyperfine states trapped are denoted by $``$ and $``$. When the interaction $`g_{}(𝐫𝐫^{})`$ is attractive, the system has been predicted to undergo a BCS transition  associated with the appearance of an order parameter. We assume the system to be below the critical temperature and the Hamiltonian $`H_0`$ treated within the BCS approximation.

The density response of the system to perturbations, caused by e.g. probe lasers or the trapping fields, is studied in the linear response regime . The density perturbation is a linear function of an external potential $`U(𝐫,t)`$ expressed through the density-density response function $`\chi (𝐫,t;𝐫^{},t^{})`$:

$$\delta n(𝐫,t)=\delta n_{}(𝐫,t)+\delta n_{}(𝐫,t)=𝑑𝐫𝑑t\chi (𝐫,t;𝐫^{},t^{})U(𝐫,t),$$
(3)
where the density operators are defined in the usual way $`n_{\alpha =,}(𝐫)=\psi _\alpha ^{}(𝐫)\psi _\alpha (𝐫)`$. Fourier transforming the previous expression one obtains

$$\delta n(𝐤,\omega )=\chi (𝐤,\omega )U(𝐤,\omega ).$$
(4)
The poles of the reponse function $`\chi `$ give the collective modes of the system while its modulus gives the spectral weight of the modes.

It is known that a perturbation of the form $`U(𝐫,t)[\psi _{}^{}(𝐫)\psi _{}(𝐫)+\psi _{}^{}(𝐫)\psi _{}(𝐫)]`$ leads in the long wavelength limit to appearance of the Anderson-Bogoliubov phonon , that is, a collective mode with energy below the superconducting gap.

In this paper we study the reponse function $`\chi (𝐤,\omega )`$ and potential $`U(𝐤,\omega )`$ for a perturbation $`H^{}`$ created by different laser probing schemes. Most of the schemes have been considered in the literature as proposed techniques for observing the superconducting gap.

## III Probing schemes

### A Laser coupling between two states

a) Coupling the two paired states. Let us assume that the states $``$ and $``$ are coupled by light . Using the Rotating-Wave-Approximation (RWA) , the interaction of the laser light with the matter fields can be described by a time-independent Hamiltonian in which the detuning $`\delta `$ plays the role of a chemical potential and a Rabi frequency $`\mathrm{\Omega }(𝐫)`$ characterises the local strength of the interaction. The perturbation Hamiltonian is $`H^{}=H_\mu +H_T`$, where the transfer Hamiltonian ($`H_T`$) and $`H_\mu `$ are given by

$`H_T={\displaystyle 𝑑𝐫\mathrm{\Omega }(𝐫)\psi _{}^{}(𝐫)\psi _{}(𝐫)}+\mathrm{\Omega }^{}(𝐫)\psi _{}^{}(𝐫)\psi _{}(𝐫)`$ (5)
$`H_\mu ={\displaystyle \frac{\delta }{2}}{\displaystyle 𝑑𝐫\psi _{}^{}(𝐫)\psi _{}(𝐫)}{\displaystyle \frac{\delta }{2}}{\displaystyle 𝑑𝐫\psi _{}^{}(𝐫)\psi _{}(𝐫)}`$ (6)
In principle, one can imagine a direct coupling of the states $``$ and $``$ by one field, in which case $`\delta `$ and $`\mathrm{\Omega }`$ correspond directly to that field. In practice, for coupling of the two hyperfine states which have an energy difference far from any laser frequency one may prefer a Raman transition. In that case the Rabi frequencies and $`\delta `$ depend on the parameters from both lasers, especially $`\mathrm{\Omega }\mathrm{\Omega }_1^{}\mathrm{\Omega }_2`$.

b) Coupling to a non-paired state. When one of the states $`\alpha =\{`$ or $`\}`$ is coupled to some excited hyperfine state $`e`$ , the perturbation reads

$`H_T={\displaystyle 𝑑𝐫\mathrm{\Omega }(𝐫)\psi _\alpha ^{}(𝐫)\psi _e(𝐫)}+\mathrm{\Omega }^{}(𝐫)\psi _e^{}(𝐫)\psi _\alpha (𝐫)`$ (7)
$`H_\mu ={\displaystyle \frac{\delta }{2}}{\displaystyle 𝑑𝐫\psi _\alpha ^{}(𝐫)\psi _\alpha (𝐫)}{\displaystyle \frac{\delta }{2}}{\displaystyle 𝑑𝐫\psi _e^{}(𝐫)\psi _e(𝐫)}.`$ (8)
c) Far off-resonant light. In  the coupled light was proposed to be almost on-resonant, that is, atomic population is transfered in the excitation process. In  coupling of the paired states with some excited states was assumed to be done by far off-resonant or very weak intensity light meaning the excited state population remains negligible and the excited state can be eliminated from the problem. Far off-resonance coupling can be treated by starting with an initial Hamiltonian of the same form as (8), then after adiabatically eliminating the excited states one arrives at $`\frac{|\mathrm{\Omega }|^2}{\delta }\psi _\alpha ^{}\psi _\alpha `$, c.f. section IV.

### B Space and time dependence of the perturbation

d) Spatial variation of the intensity. The common way of creating intensity gratings used e.g. for creating collective excitations in Bose-Einstein condensates  is to use two intersecting waves. When two intersecting waves with the same polarization but different wave-vector and frequency are coupled to the same two-level system, the total Rabi frequency is $`\mathrm{\Omega }=\mathrm{\Omega }_1e^{i𝐤_\mathrm{𝟏}𝐫i\omega _1t}+\mathrm{\Omega }_2e^{i𝐤_\mathrm{𝟐}𝐫i\omega _2t}`$. After making the RWA for the frequency $`\omega _1`$, $`\mathrm{\Omega }=\mathrm{\Omega }_1e^{i𝐤_\mathrm{𝟏}𝐫}+\mathrm{\Omega }_2e^{i𝐤_\mathrm{𝟐}𝐫+i\omega _{12}t}`$ where $`\omega _{12}=\omega _1\omega _2`$.

Also the beam profile may vary spatially in the cloud size scale. The Rabi frequency would then be for instance $`\mathrm{\Omega }(𝐫)e^{r^2/(2\sigma ^2)}`$ for a Gaussian beam shape.

As will be shown in the following section, the essential feature considering a possible density response is whether $`|\mathrm{\Omega }|^2`$ is time and space dependent. Clearly, for Bragg spectroscopy as considered above this is true. For the Gaussian beam profile, $`|\mathrm{\Omega }|^2`$ is dependent on position but not on time. In contrast, for a single laser $`\mathrm{\Omega }=|\mathrm{\Omega }|e^{i𝐤𝐫i\omega t}`$ and for a Raman excitation $`\mathrm{\Omega }|\mathrm{\Omega }_1||\mathrm{\Omega }_2|e^{i𝐤_\mathrm{𝟏}𝐫i\omega _1t}e^{i𝐤_\mathrm{𝟐}𝐫+i\omega _2t}`$, thus in both cases $`|\mathrm{\Omega }|^2`$ is a constant.

e) Time dependece of the perturbation. As discussed above, Bragg scattering schemes lead to a time-varying intensity grating, that is, the energy $`\omega _{12}`$ is transferred to the system via a perturbation that is proportional to the intensity.

Another possible source of time-dependence is the turning on of the perturbation. As is usually done in the linear response theory , we consider smooth turning on of the perturbing potential.

Direct intensity modulation of the probing laser(s) would be one more source of time-dependence.

## IV Transforming $`H_T`$ into a density perturbation

In order to use linear response theory we have to convert all the perturbing Hamiltonians of the previous section into a density perturbation.

a) The perturbation term $`H_T`$ involves products of the field operators of the different atomic states. In the mean field approach they correspond to the Fock terms $`<\psi _{}^{}(𝐫)\psi _{}(𝐫)>`$ which are zero. Therefore no linear response theory can be directly applied.

A perturbation of the type (note that the two-component $`\psi `$ here is not the same as that in the standard BCS-theory, used for instance in the Appendix)

$`\left(\begin{array}{cc}\psi _{}^{}(𝐫)& \psi _{}^{}(𝐫)\end{array}\right)\left(\begin{array}{cc}\delta /2& \mathrm{\Omega }(𝐫)\\ \mathrm{\Omega }^{}(𝐫)& \delta /2\end{array}\right)\left(\begin{array}{c}\psi _{}(𝐫)\\ \psi _{}(𝐫)\end{array}\right)`$ (14)
$``$ $`\psi ^{}(𝐫)W\psi (𝐫)`$ (15)
can be digonalized for any $`\delta `$ and $`\mathrm{\Omega }`$ by an appropiate rotation in the space of the two states $``$ and $``$, $`\psi =𝒰`$$`\stackrel{~}{\psi }`$ and $`𝒰`$=$`R_z(\alpha )R_y(\beta )R_z(\gamma )`$ with $`R_j(\theta )=\mathrm{exp}(i\sigma _j\theta )`$, $`j=\{x,y,z\}`$ and $`\sigma _j`$ is the corresponding Pauli spin matrix. The transformation matrix is

$`𝒰=\left(\begin{array}{cc}e^{i(\alpha +\gamma )/2}\mathrm{cos}\beta /2& e^{i(\alpha \gamma )/2}\mathrm{sin}\beta /2\\ e^{i(\alpha \gamma )/2}\mathrm{sin}\beta /2& e^{i(\alpha +\gamma )/2}\mathrm{cos}\beta /2\end{array}\right).`$ (18)
The rotated perturbation matrix is $`𝒰^{}`$W$`𝒰`$. The off-diagonal terms are

$$e^{i\gamma }((\delta /2)\mathrm{sin}(\beta )e^{i\alpha }\mathrm{\Omega }^{}\mathrm{sin}^2(\beta /2)+\mathrm{\Omega }\mathrm{cos}^2(\beta /2)e^{i\alpha })$$
(19)
and the complex conjugate. Off-diagonal terms that are zero can be obtained by a rotation with the Euler angles $`\{\gamma =0,\mathrm{cos}\alpha =\frac{\mathrm{\Omega }+\mathrm{\Omega }^{}}{2|\mathrm{\Omega }|},\mathrm{tan}\beta =\frac{2|\mathrm{\Omega }|}{\delta }\}`$ which yields

$`\left(\begin{array}{cc}\sqrt{(\frac{\delta }{2})^2+\mathrm{\Omega }\mathrm{\Omega }^{}}& 0\\ 0& \sqrt{\mathrm{\Omega }\mathrm{\Omega }^{}+(\frac{\delta }{2})^2}\end{array}\right)`$ (22)
for the perturbation matrix. This procedure is extensively used in describing the interaction of laser light with a two-level system. The form $`|\mathrm{\Omega }|^2/\delta `$ for far off-resonant light comes from (22) for $`|\delta |>>|\mathrm{\Omega }|`$.

One has to check whether the rotated states feel the same two-body interaction that the non-rotated ones, that is, whether the rotated states are still described by a Hamiltonian of the initial form which implies BCS pairing. Assuming, for simplicity, that $`\psi _{}`$ and $`\psi _{}`$ see the same potential $`V_\alpha \mu _\alpha `$ , we obtain $`\psi _{}^{}(𝐫)\psi _{}(𝐫)+\psi _{}^{}(𝐫)\psi _{}(𝐫)\stackrel{~}{\psi }_{}^{}(𝐫)\stackrel{~}{\psi }_{}(𝐫)+\stackrel{~}{\psi }_{}^{}(𝐫)\stackrel{~}{\psi }_{}(𝐫)`$. The kinetic energy term transforms as $`\psi _{}^{}(𝐫)\psi _{}(𝐫)+\psi _{}^{}(𝐫)\psi _{}(𝐫)\stackrel{~}{\psi }_{}^{}(𝐫)\stackrel{~}{\psi }_{}(𝐫)+\stackrel{~}{\psi }_{}^{}(𝐫)\stackrel{~}{\psi }_{}(𝐫)`$ if the momentum k transfered by the laser is small. The extra terms coming in the transformation of the kinetic energy are $`k,k^2`$. The momentum of the laser and the recoil energy are very small compared to the momentum of the atoms participating in the process (typically close to the Fermi-surface, $`k_F`$).

The interaction term in the Hamiltonian transforms $`\psi _{}^{}(𝐫)\psi _{}^{}(𝐫^{})\psi _{}(𝐫^{})\psi _{}(𝐫)\stackrel{~}{\psi }_{}^{}\stackrel{~}{\psi }_{}^{}\stackrel{~}{\psi }_{}\stackrel{~}{\psi }_{}`$ in case of a contact interaction $`g(𝐫𝐫^{})\delta (𝐫𝐫^{})`$. As $`\stackrel{~}{\psi }_\alpha `$ is a linear combination of the two species $``$ and $``$, one might expect interaction terms between the same species for the new matter fields (interactions between both $`\stackrel{~}{\psi }_{}`$ and $`\stackrel{~}{\psi }_{}`$, and between $`\stackrel{~}{\psi }_{}`$ and $`\stackrel{~}{\psi }_{}`$). But the interactions of the type $`\stackrel{~}{\psi }_\alpha \stackrel{~}{\psi }_\alpha `$ are forbidden by the fermionic behaviour (commutation rules are preserved by the rotation). The requirement of a contact interaction is obvious here: $`\stackrel{~}{\psi }_\alpha (𝐫)\stackrel{~}{\psi }_\alpha (𝐫)=0`$, but $`\stackrel{~}{\psi }_\alpha (𝐫)\stackrel{~}{\psi }_\alpha (𝐫^{})`$ can be non-zero. As a summary, since the fermionic field commutator relations are preserved by the rotation, for contact interaction $`\stackrel{~}{H}_0=H_0`$.

We have thus shown that the system can be transformed in such a way that one can study the collective mode spectra by linear response theory for a BCS system for a perturbation of the type $`U(n_{}n_{})`$, where $`U(𝐤,\omega )=\left[\sqrt{\mathrm{\Omega }(𝐫,t)\mathrm{\Omega }(𝐫,t)^{}+(\frac{\delta }{2})^2}\right]`$, and the rest of the Hamiltonian still corresponds to the BCS Hamiltonian. The main difference to the standard linear density response treatment for an BCS system is the minus sign between $`n_{}`$ and $`n_{}`$. We will show that this leads to zero overall response in the homogeneous case.

b) When diagonalizing the interaction Hamiltonian of Eq.(8) in the same way as in a), the initial Hamiltonian $`H_0`$ is not preserved. We thus obtain a perturbation $`U(n_\alpha n_e)`$ with $`U(𝐤,\omega )`$ as given above, but with the initial Hamiltonian modified, $`\stackrel{~}{H}_0H_0`$.

c) In the case of coupling to an excited state by far-detuned light the rotation does not modify the initial Hamiltonian because at the limit $`|\delta |>>|\mathrm{\Omega }|`$ the transformed states are very close to the initial ones, $`\psi _\alpha \stackrel{~}{\psi }_\alpha `$ and $`\psi _e\stackrel{~}{\psi }_e`$, moreover the excited state population is assumed to be negligible. Thus the light-atom interaction Hamiltonian is in the form of a density perturbation, and the initial Hamiltonian is preserved.

In all these cases, $`U(𝐤,\omega )`$ can be either a constant $`U(𝐤,\omega )=U^{}`$ or dependent on $`𝐤`$ and $`\omega `$. The latter happens when the intensity $`|\mathrm{\Omega }|^2`$ is space and time dependent. Note that the k-dependence $`\mathrm{\Omega }e^{i𝐤𝐫}`$ is not sufficient, that is, the laser momentum given in a single-laser or Raman two-level excitation (near-resonant or off-resonant) is not enough to transfer momentum through a density perturbation.

## V Results

### A GRPA for the case of light coupled between the two paired states

The response of the system to the perturbation $`U(k,\omega )(n_{}n_{})`$ in a homogenous system is calculated in the Appendix following the derivation in. The result is that for the homogeneous case, the responses of two components are equal but with opposite signs and no collective modes are excited at any frequencies.

The response of one of the spin components is $`L_{11}=\frac{(a+b)(1+g_0R)}{(1g_0R)(g_0+a+b)}`$, where we use the integrals a, b, c, d, R as defined in . The reponse of the spin component has the factor (a+b) which at T=0 in the low q,$`\omega `$ region is $`a+bN(ϵ_F)c_B^2q^2/(12\mathrm{\Delta }^2)`$ where $`N(ϵ_F)`$ is the density of states at the Fermi surface and the sound velocity $`c_B^2=v_F^2/3`$. This factor is very small and it makes the response for one spin component negligible even without the cancellation effect. This is because the superfluid cannot participate in the relative motion of the two spin components. Above $`T_c`$, b vanishes and the response of one of the spins is not negligible anymore, as the integral $`a`$ is some fraction of $`N(ϵ_F)`$. However, due to the cancellation the overall response is zero.

We have considered the homogeneous case, that is, when the trapping potential treated within the local density approximation is sufficient. The strongly trapped case can be approached by using multipole expansions : the overall response is not zero for spin-dipole excitations in the trapped case above $`T_c`$. The difference to our results is a consequence of the geometry, as in the homogeneous case we get the same spatially uniform response in opposite directions for both spin species. As $`T0`$, one gets vanishing response also for the inhomogeneous case , at least at energies below the gap, because one needs to break pairs in order to have relative motion of the two components.

### B Density response for different probing schemes

a)Coupling the paired states. The perturbation reads $`U(n_{}n_{})`$ for the rotated states, where $`U(𝐤,\omega )=[\sqrt{|\mathrm{\Omega }(\text{r},t)|^2+\delta ^2}]`$. When the light is far off-resonant, one can take the limit $`U(𝐤,\omega )[|\mathrm{\Omega }(𝐫,t)|^2/\delta ]`$. In the case of homogeneous laser intensity $`\mathrm{\Omega }=|\mathrm{\Omega }|e^{i𝐤𝐫}`$, $`U\delta (\omega )\delta (\text{k})`$. There is no density reponse simply because U does not provide the momentum and energy for a collective excitation.

A different behaviour arises when using a Bragg scattering scheme to provide an intensity grating. Assuming $`|\mathrm{\Omega }_1|=|\mathrm{\Omega }_2|`$ ($`\mathrm{\Omega }_1,\mathrm{\Omega }_2`$ as introduced in section III d)) one obtains $`|\mathrm{\Omega }|^2=|\mathrm{\Omega }_1|^2(2+2\mathrm{cos}(\omega _{12}t𝐤_{\mathrm{𝟏𝟐}}𝐫))|\mathrm{\Omega }_1|^22\mathrm{cos}(\omega _{12}t𝐤_{\mathrm{𝟏𝟐}}𝐫)`$ and the density response $`\chi (\omega ,𝐤)`$ has to be analyzed at $`\chi (\omega _{12},𝐤_{\mathrm{𝟏𝟐}})`$. The density response $`\chi `$ gives zero in the homogeneous case for all temperatures, as one component annihilates the response of the other, and for the inhomogeneous case well below $`T_c`$  as considered in the previous subsection.

For a general case where the beam profile of the laser light gives a spatially varying intensity of the form $`|\mathrm{\Omega }(𝐫)|^2e^{\frac{r^2}{\sigma ^2}}`$, the far off-resonant perturbation reads $`U(𝐤,\omega )=\frac{|\mathrm{\Omega }(k)|^2}{\delta }e^{\frac{\sigma ^2k^2}{4}}`$. If $`\sigma `$ is some fraction $`f`$ of the trap size $`R_{TF}`$, one can estimate the momentum given by $`p=\frac{\mathrm{}}{\sigma }`$. A collective mode with an energy of the order of the gap energy (e.g. half the gap) has momentum $`q\sqrt{3}\frac{\mathrm{\Delta }_G}{v_F}`$  and the ratio between them is $`\frac{p}{q}=\frac{\mathrm{}\omega _T}{f\mathrm{\Delta }_G}`$. For typical parameters, this could be a non-negligible number, however, this alone is not enough to produce a density response because of the lack of the time dependence, i.e. $`U(\omega ,𝐤)\delta (\omega )f(𝐤)`$.

b) Coupling to an excited state. When diagonalizing the interaction Hamiltonian in this case the initial Hamiltonian was modified. Therefore the linear response calculation for a BCS system cannot be directly applied. One may guess, however, that the arguments about the perturbation $`U`$ being time- and space-dependent will hold also in this case even when the response function $`\chi `$ would have a different form. Therefore we do not expect a density response except when a Bragg spectroscopy scheme is used.

c) Far off-resonant coupling. When both paired states are coupled to some excited states and the light is far detuned, the density perturbation potential is proportional to $`|\mathrm{\Omega }|^2/\delta `$. No collective mode of the system is excited if there is no intensity variation in space and time. When considering intensity modulations of the type $`|\mathrm{\Omega }|^2|\mathrm{\Omega }_1|^22\mathrm{cos}(\omega t𝐤𝐫)`$, an Anderson-Bogoliubov phonon can be excited .

## VI Conclusions

We have reviewed the density response of some typical laser excitation schemes. We have shown that since the BCS Hamiltonian for a contact interaction is preserved under a rotation, most of the considered laser excitations can be expressed in terms of a perturbation acting on the density. In this form, the perturbation potential is proportional to $`|\mathrm{\Omega }|^2`$ where $`\mathrm{\Omega }`$ is the (effective) Rabi frequency. Therefore, even when the laser light provides momentum and energy ($`\mathrm{\Omega }e^{i𝐤_𝐋𝐫}e^{i\omega _Lt}`$), the transformed potential acting as a density perturbation is not time- and space-dependent. This leads to absence of a density response whenever $`|\mathrm{\Omega }|^2`$ is a constant spatially and temporally. This makes many proposed laser-probing schemes well suited for observing the superconducting gap since they do not induce below-gap collective excitations.

For Bragg scattering, $`|\mathrm{\Omega }|^2f(𝐫,t)`$. In this case Anderson-Bogoliubov phonons can be excited, in general. The exception is the case when the laser(s) couple between the two paired components of the gas. We have shown that, in the homogeneous case, the density response becomes zero because the contributions of the components cancel each other. In a harmonic trap, spin-dipole response is predicted for temperatures near $`T_c`$. Therefore, the presence or absence of low-energy collective excitations under a perturbation of the type $`U(\omega ,k)(n_{}n_{})`$ could be used to observe whether the trapped system can be approximated by a homogeneous system (local density approximation) or whether the trapping effects are dominant.

## VII Appendix

The formalism used is based on the functional differentiation technique as described in  following the derivation in . The response function $`\chi `$ is viewed as a functional derivative of the one-particle matrix Green’s functions $`G`$ with respect to the external field $`U(𝐫,t)`$. To allow for pairing in the superconducting phase, it is convenient to work with a single-particle Green’s function G given by a $`2\times 2`$ matrix defined as

$`G(1,2)T\mathrm{\Psi }(1)\mathrm{\Psi }^{}(2)=\left(\begin{array}{cc}G_{}(1,2)& F(1,2)\\ F^{}(1,2)& G_{}(2,1)\end{array}\right),`$ (25)
where $`\mathrm{\Psi }=(\psi _{}\psi _{}^{})^T`$ and $`1(𝐫_\mathrm{𝟏},\tau _1)`$. Imaginary times (Matsubara formalism) are used so that one can deal with finite temperatures. In the absence of external fields, the equal time ($`\tau _2=\tau _1^+`$) single-particle Green’s function components reduce to $`F(𝐫)\psi _{}(𝐫)\psi _{}(𝐫)`$, the s-wave order parameter, and $`G_,(𝐫)=n_,(𝐫)`$.

¿From the equation of motion of the Green’s function one gets the generalized Dyson equation for G(1,2) in terms of $`G_0`$, the non-interacting single-particle Green’s function, the matrix self-energy $`\mathrm{\Sigma }(1,2)`$ which is evaluated in the pairing approximation (Hartree-Fock-Bogoliubov), and W(1) which is the the external perturbing field matrix.

The density response matrix is obtained in RPA by taking the functional derivative of the Green’s function with respect to the external field U. One can define the three-point correlation function $`L(1,2,5)\sigma _3\frac{\delta G(1,2)}{\delta U(5)}`$, whose limit $`L(1,2)L(1,1^+,2)`$ will give the density-density response function $`\chi (1,2)=L_{11}(1,2)+L_{22}(1,2)`$. Here $`\sigma _3`$ is the third Pauli matrix.

When deriving the Dyson equation for G(1,2) with respect to U in order to get the three-point correlation function L, both the self energy matrix and the external field matrix W contribute. The lowest-order (single-bubble) result  for $`L`$ is given by

$$L^0(1,2,5)=\sigma _3𝑑\overline{3}𝑑\overline{4}G(1,\overline{3})\frac{\delta W(\overline{3},\overline{4})}{\delta U(5)}G(\overline{4},2).$$
(26)
For the Anderson-Bogoliubov phonon (see ) the pertubation matrix $`W(1)=U(1)\sigma _3`$, and $`L_{AB}^0(1,2,5)=\stackrel{~}{G}(1,5)\stackrel{~}{G}(5,2)`$, where we have introduced $`\stackrel{~}{G}\sigma _3G`$. In our case $`L^0(1,2,5)=\stackrel{~}{G}(1,5)G(5,2)`$ because the perturbation matrix $`W(1)=U(1)`$ due to the minus sign in $`U(n_{}n_{})`$. ¿From now on, we denote by the subindex AB the quantities that are calculated for the perturbation matrix $`W(1)=U(1)\sigma _3`$.

It is useful to rewrite the GRPA integral equations in terms of irreducible two-particle Green’s functions $`\overline{L}_{ij}(1,2,5)`$, and for the homogeneous system, to Fourier transform in order to solve coupled equations. The result reads

$$L_{ij}(𝐪,i\omega _n)=\overline{L}_{ij}(𝐪,i\omega _n)+\frac{\overline{L}_{ABij}(𝐪,i\omega _n)g_{}(𝐪)\overline{L}_{ll}(𝐪,i\omega _n)}{1g_{}(𝐪)\overline{L}_{ABll}(𝐪,i\omega _n)},$$
(27)
where equal indices mean summation over the possible values. We denote by the subindex AB the three-point correlation functions for the Anderson-Bogoliubov phonon type of perturbation.

The sum of diagonal terms reduces to

$`\chi (𝐪,i\omega _n)`$ $`=`$ $`L_{11}(𝐪,i\omega _n)+L_{22}(𝐪,i\omega _n)`$ (28)
$`=`$ $`{\displaystyle \frac{\overline{L}_{ii}(𝐪,i\omega _n)}{1g_{}(𝐪)\overline{L}_{ABll}(𝐪,i\omega _n)}}.`$ (29)
This shows that all linear perturbations have the same poles as the Anderson-Bogoliubov phonon (plus possibly some additional ones), but the spectral weight depends on the trace of the irreducible two-particle Green’s functions $`\overline{L}_{ii}(𝐪,i\omega _n)`$ of the specific perturbation.

For a contact interaction ($`g_{}(𝐪)=g_0`$), the equation for the irreducible correlation function  reduces to a set of linear algebraic equations

$$\overline{L}_{ij}(𝐪,i\omega _n)=L_{ij}^0(𝐪,i\omega _n)L_{ABiklj}^0(𝐪,i\omega _n)g_0\overline{L}_{kl}(𝐪,i\omega _n),$$
(30)
where we have defined the four-index tensor $`L_{ABiklj}^0`$, indicating components of the two factor matrices. Defining column vectors $`\overline{L}=\left(\begin{array}{c}\overline{L}_{11}\overline{L}_{12}\overline{L}_{21}\overline{L}_{22}\end{array}\right)^T`$ and a $`4\times 4`$ matrix $`L_{(AB)mn}^0`$ as in , Eq.(30) reduces to

$$\overline{L}=[I+g_0L_{AB}^0]^1L^0.$$
(31)
As discussed in  one can get $`L_{ijkl}^0`$ by Matsubara frequency summations, (Eq.(4.27) in ) and using the symmetry properties then reduce just to 6 independent elements. For the weak-coupling limit and contact interaction, the independent elements are reduced to four: a, b, c, d which are integrals defined in  (Eq.A8). $`L_m^0`$ can be calculated from the $`4\times 4`$ matrix by $`L_m^{0T}=\left(\begin{array}{c}L_{1ll1}^0L_{1ll2}^0L_{2ll1}^0L_{2ll2}^0\end{array}\right)`$.

Now we derive the $`4\times 4`$ matrix for our case when $`W(1)=U`$ by reconsidering the symmetry properties of the matrix elements in order to derive the vector $`L_m^0`$ that we insert in Eq.(31). The sum of the first and last component of $`\overline{L}_m`$ contribute to Eq.(29) and characterize the spectral weight and perhaps some additional poles of the response function. Such a perturbation gives from (26) the lowest order correlation function $`L^0(1,2,5)=\stackrel{~}{G}(1,5)G(5,2)`$ (c.f. $`L_{AB}^0(1,2,5)=\stackrel{~}{G}(1,5)\stackrel{~}{G}(5,2)`$). We calculate the $`4\times 4`$ matrix $`L_{mn}^0`$ using the symmetry properties of the new $`L_{ijkl}^0`$ elements, in the same fashion as in  (see appendix A there) obtaining

$`L_{mn}^0`$ $`=`$ $`\left(\begin{array}{cccc}a& c& c& b\\ c& d& b& c\\ c& b& d& c\\ b& c& c& a\end{array}\right).`$ (36)
It leads to $`L_m^{0T}=\left(\begin{array}{cccc}a+b& 0& 0& ab\end{array}\right)`$. Multiplying $`L_m^0`$ by $`[I+g_0L_{ABmn}^0]^1`$ as in (31) one gets $`\overline{L}_m^T=\left(\begin{array}{cccc}\frac{a+b}{g_0+a+b}& 0& 0& \frac{a+b}{g_0+a+b}\end{array}\right)`$. This means $`\overline{L}_{11}=\overline{L}_{22}`$, and the spectral weight in (29) vanishes and $`\chi =0`$.

Inserting the well-known result  $`\overline{L}_{ABll}=\frac{2R}{1+g_0R}`$ into Eq.(29), one obtains the spin density response $`L_{11}=\frac{(a+b)(1+g_0R)}{(1g_0R)(g_0+a+b)}`$, where $`R=A+\frac{4c^2}{1+g_0B}`$ and $`Bb+d`$, $`Aab`$ are defined in .

In the case of a general perturbation

$$W=U\left(\begin{array}{cc}1& 0\\ 0& c_u\end{array}\right),$$
(37)
where $`c_u`$ is any real number, the single-bubble result is

$$L_{ijkl}^0=\stackrel{~}{G}_{ij}\left(\left(\begin{array}{cc}1& 0\\ 0& c_u\end{array}\right)_{ks}G\right)_{sl}.$$
(38)
By calculating the symmetry properties of the elements of $`L_{mn}`$, one obtains $`L_m^{0T}=\left(\begin{array}{cccc}a+c_ub& cc_uc& cc_uc& c_uab\end{array}\right)`$ which leads to

$$\overline{L}_{11}\overline{L}_{22}=\frac{(c_u+1)(a+b)}{1+g_0(a+b)}$$
(39)
and

$$\overline{L}_{11}+\overline{L}_{22}=\frac{(1c_u)R}{1+g_0R}.$$
(40)
The reponse for one spin would be

$$L_{11}=\frac{r_u}{(1+g_0R)(1+g_0B)}$$
(41)
where $`r_u=g_0AB+g_02c^2\frac{(1c_u)+2g_0(a+b)}{1+g_0(a+b)}+\frac{(a+c_ub)+g_0(a^2+c_ub^2+db(1+c_u))}{1+g_0(a+b)}`$.