text:
# Correlation functions for diffusion-limited annihilation, 𝐴+𝐴→0

## Abstract

The full hierarchy of multiple-point correlation functions for diffusion-limited annihilation, $`A+A0`$, is obtained analytically and explicitly, following the method of intervals. In the long time asymptotic limit, the correlation functions of annihilation are identical to those of coalescence, $`A+AA`$, despite differences between the two models in other statistical measures, such as the interparticle distribution function.

The kinetics of nonequilibrium processes, in particular diffusion-limited reactions, have attracted much recent interest . Because of the lack of a comprehensive approach for the study of such systems, models that yield to an exact analysis are of prime importance. In this respect, none have been studied more than diffusion-limited annihilation, $`A+A0`$ , and coalescence, $`A+AA`$ . Known exact results include the time dependence of the particle concentration and the two-point correlation function (for finding two particles at two different points, simultaneously). It has also been shown that the full hierarchy of $`n`$-point correlation functions for the two processes is identical , but explicit expressions for $`n>3`$ are unavailable.

Here we attack the problem of correlation functions for annihilation, using the method of parity intervals (or even/odd intervals) . We recover the identity relation of the $`n`$-point correlation functions for annihilation and coalescence, and we derive explicit expressions, valid in the long time asymptotic limit, for all $`n`$.

Consider the annihilation model, defined on the line $`\mathrm{}<x<\mathrm{}`$. Particles $`A`$ are represented by points which perform unbiased diffusion with a diffusion constant $`D`$. When two particles meet they annihilate instantly. Since the reaction step is infinitely fast, the system models the diffusion-limited annihilation process $`A+A0`$.

An exact treatment of the problem is possible through the method of parity intervals . The key parameter is $`G(x,y;t)`$—the probability that the interval $`[x,y]`$ contains an even number of particles at time $`t`$ . Particles near the edges of an interval may diffuse into or out of the interval, affecting the probability $`G`$. (On the other hand, reactions inside the interval cannot affect its parity.) With this observation in mind, one derives a rate equation for the probability $`G(x,y;t)`$ :

$$\frac{}{t}G(x,y;t)=D(\frac{^2}{x^2}+\frac{^2}{y^2})G(x,y;t).$$
(1)
The annihilation reaction imposes the boundary condition

$$\underset{xy\mathrm{or}yx}{lim}G(x,y;t)=1,$$
(3)
and $`G`$ must also obey the conditions required from a probability density function. If the initial distribution of particles is random, then $`G(x,y;0)=\frac{1}{2}+\frac{1}{2}\mathrm{exp}[2c_0(yx)]`$, where $`c_0`$ is their initial density. In this case we have the additional boundary condition

$$\underset{x\mathrm{}\mathrm{or}y\mathrm{}}{lim}G(x,y;t)=\frac{1}{2}.$$
(4)
From $`G(x,y;t)`$ one can derive the particle concentration:

$$\rho (x;t)=\frac{}{y}G(x,y;t)|_{y=x}.$$
(5)
Let $`\rho _n(x_1,x_2,\mathrm{},x_n;t)`$ be the $`n`$-point density correlation function, for finding particles at each of the locations $`x_1,x_2,\mathrm{},x_n`$ at time $`t`$. The particle concentration, $`\rho (x,t)\rho _1(x,t)`$, represents merely the first term in the hierarchy $`\{\rho _n\}`$, $`n=1,2,\mathrm{}`$.

The correlation functions may be obtained from a generalization of the method of parity intervals, in the following way. Let $`H_n(x_1,y_1,\overline{x_2,y_2},\mathrm{},x_n,y_n;t)`$ be the joint probability that the interval $`[x_1,y_1]`$ contains an even number of particles, $`[x_2,y_2]`$ contains an odd number, etc., (odd intervals are denoted by an overbar), at time $`t`$. The intervals are non-overlapping, and ordered: $`x_1<y_1<\mathrm{}<x_n<y_n`$. Let $`F_n(x_1,y_1,x_2,y_2,\mathrm{},x_n,y_n;t)`$ denote the probability that the total number of particles contained in $`_{i=1}^n[x_i,y_i]`$, is even. Thus,

$`F_1(x_1,y_1;t)`$ $`=`$ $`H_1(x_1,y_1;t)=G(x_1,y_1;t),`$ (7)
$`F_2(x_1,y_1,x_2,y_2;t)`$ $`=`$ $`H_2(x_1,y_1,x_2,y_2;t)+H_2(\overline{x_1,y_1},\overline{x_2,y_2};t),`$ (8)
and, in general, $`F_n`$ is expressible as a sum of $`2^{n1}`$ $`H_n`$ functions, corresponding to the different combinations of interval parities that contribute to a total number of particles that is even. Then, in view of Eq. (5), the $`n`$-point correlation function is given by

$$\rho _n^{\mathrm{anni}}(x_1,\mathrm{},x_n;t)=\frac{(1)^n}{2^{n1}}\frac{^n}{y_1\mathrm{}y_n}F_n(x_1,y_1,\mathrm{},x_n,y_n;t)|_{y_1=x_1,\mathrm{},y_n=x_n}.$$
(9)
The $`H_n`$ satisfy a $`2n`$-dimensional diffusion equation, analogous to Eq. (1), and for similar reasons. However, the boundary conditions of this equation are complicated by the following fact. For $`y_ix_{i+1}`$, a particle moving from the $`i`$-th interval to the $`(i+1)`$-th interval, or vice versa, would flip the parity of the two adjacent intervals. On the other hand, $`F_n`$ satisfies the same diffusion equation as $`H_n`$,

$$\frac{}{t}F_n(x_1,y_1,\mathrm{},x_n,y_n;t)=D(\frac{^2}{x_1^2}+\frac{^2}{y_1^2}+\mathrm{}+\frac{^2}{x_n^2}+\frac{^2}{y_n^2})F_n,$$
(10)
but the boundary conditions are simpler: $`F_n`$ contains also the case where the parity of the intervals $`i`$ and $`(i+1)`$ is flipped, so it is not affected by a particle hopping between the two intervals. If interval $`i`$ is shrunk to zero, Eq. (3) yields the boundary condition

$$\underset{x_iy_i\mathrm{or}y_ix_i}{lim}F_n(x_1,y_1,\mathrm{},x_n,y_n;t)=F_{n1}(x_1,y_1,\mathrm{},\overline{)}x_i,\overline{)}y_i,\mathrm{},x_n,y_n;t),$$
(12)
where we use the notation that crossed out arguments (e.g. $`\overline{)}x_i`$) have been removed. If the endpoints of two adjacent intervals are brought together, the intervals merge, resulting in the boundary condition

$$\underset{y_ix_{i+1}\mathrm{or}x_{i+1}y_i}{lim}F_n(x_1,y_1,\mathrm{},x_n,y_n;t)=F_{n1}(x_1,y_1,\mathrm{},\overline{)}y_i,\overline{)}x_{i+1},\mathrm{},x_n,y_n;t).$$
(13)
Finally, for a random initial distribution of particles, we have

$$\underset{x_1\mathrm{}\mathrm{or}y_n\mathrm{}}{lim}F_n(x_1,y_1,\mathrm{},x_n,y_n;t)=\frac{1}{2}.$$
(14)
The $`F_n`$ are tied together in an hierarchical fashion through the boundary conditions (10a) and (10b): one must know $`F_{n1}`$ in order to compute $`F_n`$. At the root of the hierarchy, $`F_1=G`$ is obtained from Eqs. (1) and (1).

The problem posed by Eqs. (1), (1), (10), (10) is similar to that of diffusion-limited coalescence, $`A+AA`$ . In that case one defines $`E_n(x_1,y_1,\mathrm{},x_n,y_n;t)`$ as the joint probability of finding the intervals $`[x_i,y_i]`$, $`i=1,2,\mathrm{},n`$, empty at time $`t`$. $`E_1(x,y;t)E(x,y;t)`$ satisfies the same equation as Eq. (1):

$$\frac{}{t}E(x,y;t)=D(\frac{^2}{x^2}+\frac{^2}{y^2})E(x,y;t),$$
(15)
with the boundary conditions

$`\underset{xy\mathrm{or}yx}{lim}E(x,y;t)`$ $`=`$ $`1,`$ (17)
$`\underset{x\mathrm{}\mathrm{or}y\mathrm{}}{lim}E(x,y;t)`$ $`=`$ $`0.`$ (18)
Note the difference between the boundary conditions (4) and (18). Likewise, $`E_n`$ satisfies the same equation as Eq. (10):

$$\frac{}{t}E_n(x_1,y_1,\mathrm{},x_n,y_n;t)=D(\frac{^2}{x_1^2}+\frac{^2}{y_1^2}+\mathrm{}+\frac{^2}{x_n^2}+\frac{^2}{y_n^2})E_n,$$
(19)
with boundary conditions analogous to (10a) and (10b),

$`\underset{x_iy_i\mathrm{or}y_ix_i}{lim}E_n(x_1,y_1,\mathrm{},x_n,y_n;t)`$ $`=E_{n1}(x_1,y_1,\mathrm{},\overline{)}x_i,\overline{)}y_i,\mathrm{},x_n,y_n;t),`$ (21)
$`\underset{y_ix_{i+1}\mathrm{or}x_{i+1}y_i}{lim}E_n(x_1,y_1,\mathrm{},x_n,y_n;t)`$ $`=E_{n1}(x_1,y_1,\mathrm{},\overline{)}y_i,\overline{)}x_{i+1},\mathrm{},x_n,y_n;t).`$ (22)
but

$$\underset{x_1\mathrm{}\mathrm{or}y_n\mathrm{}}{lim}E_n(x_1,y_1,\mathrm{},x_n,y_n;t)=0,$$
(23)
instead of (10c). The $`n`$-point correlation function for coalescence is

$$\rho _n^{\mathrm{coal}}(x_1,\mathrm{},x_n;t)=(1)^n\frac{^n}{y_1\mathrm{}y_n}E_n(x_1,y_1,\mathrm{},x_n,y_n;t)|_{y_1=x_1,\mathrm{},y_n=x_n}.$$
(24)
Eqs. (1), (1), (10), (10), and (15) – (19) imply that the solutions for $`F_n`$ and $`E_n`$ are simply related:

$$F_n(x_1,y_1,\mathrm{},x_n,y_n;t)=\frac{1}{2}+\frac{1}{2}E_n(x_1,y_1,\mathrm{},x_n,y_n;t),$$
(25)
provided that the same relation holds also for the initial conditions. Suppose that the initial distribution of particles is random, with initial concentration $`c_0^{\mathrm{anni}}`$ for annihilation, and $`c_0^{\mathrm{coal}}`$ for coalescence. Then $`E_n(x_1,\mathrm{},y_n;0)=\mathrm{exp}[c_0^{\mathrm{coal}}(y_1x_1+\mathrm{}+y_nx_n)]`$, while $`F_n(x_1,\mathrm{},y_n;0)=\frac{1}{2}+\frac{1}{2}\mathrm{exp}[2c_0^{\mathrm{anni}}(y_1x_1+\mathrm{}+y_nx_n)]`$. Thus, the relation (25) is satisfied if $`c_0^{\mathrm{anni}}=\frac{1}{2}c_0^{\mathrm{coal}}`$. Moreover, Eqs. (9) and (24) imply that in this case

$$\rho _n^{\mathrm{anni}}(x_1,\mathrm{},x_n;t)=(\frac{1}{2})^n\rho _n^{\mathrm{coal}}(x_1,\mathrm{},x_n;t),$$
(26)
for all $`n`$. In other words, the n-point correlation functions for annihilation and coalescence are identical, as already found by others .

We now produce explicit expressions for the $`n`$-point correlation functions in the long time asymptotic limit. Recall first the solution for $`E_n`$. For $`n=2`$ the solution is ,

$`E_2(x_1,y_1,x_2,y_2;t)=`$ (27)
$`E(x_1,y_1;t)E(x_2,y_2;t)E(x_1,x_2;t)E(y_1,y_2;t)+E(x_1,y_2;t)E(y_1,x_2;t),`$ (28)
where $`E(x,y;t)`$ is the solution of Eqs. (15), (15). Generally, for $`n2`$ ,

$$E_n(x_1,y_1,\mathrm{},x_n,y_n;t)=\underset{p=1}{\overset{(2n1)!!}{}}\sigma _pE(z_{1,p},z_{2,p};t)E(z_{3,p},z_{4,p};t)\mathrm{}E(z_{2n1,p},z_{2n,p};t),$$
(29)
where $`z_{1,p},z_{2,p},\mathrm{},z_{2n,p}`$ is an ordered permutation, $`p`$, of the variables $`x_1,y_1,\mathrm{},x_n,y_n`$, such that

$$z_{1,p}<z_{2,p},z_{3,p}<z_{4,p},\mathrm{},z_{2n1,p}<z_{2n,p},\mathrm{and}z_{1,p}<z_{3,p}<z_{5,p}\mathrm{}<z_{2n1,p}.$$
(30)
There are exactly $`(2n1)!!=13\mathrm{}(2n1)`$ such permutations. $`\sigma _p`$ is $`+1`$ for even permutations (permutations that require an even number of exchanges between pairs of variables), or $`1`$ for odd permutations. Alternatively, the $`E_n`$ may be obtained through the recursion relation:

$`E_n(x_1,y_1,\mathrm{},x_n,y_n;t)=`$ $`+{\displaystyle \underset{j=1}{\overset{n}{}}}E(x_1,y_j;t)E_{n1}(\overline{)}x_1,y_1,\mathrm{},x_j,\overline{)}y_j,\mathrm{},x_n,y_n;t)`$ (32)
$`{\displaystyle \underset{j=2}{\overset{n}{}}}E(x_1,x_j;t)E_{n1}(\overline{)}x_1,y_1,\mathrm{},\overline{)}x_j,y_j,\mathrm{},x_n,y_n;t),`$
then $`\rho _n`$ may be computed through the relation (24), or (9) and (25).

Consider the long-time asymptotic limit, where

$$E(x,y;t)=\mathrm{erfc}\left(\frac{yx}{\sqrt{8Dt}}\right).$$
(33)
Then, the long-time asymptotic $`n`$-point correlation function is:

$$\rho _n(x_1,y_1,\mathrm{},x_n,y_n;t)\underset{t\mathrm{}}{}(\rho )^n\underset{p=1}{\overset{(2n1)!!}{}}\sigma _pC(z_{1,p},z_{2,p};t)C(z_{3,p},z_{4,p};t)\mathrm{}C(z_{2n1,p},z_{2n,p};t),$$
(35)
where

$$\rho =\rho _1(x;t)=\{\begin{array}{cc}1/\sqrt{8\pi Dt}\hfill & \text{annihilation,}\hfill \\ 1/\sqrt{2\pi Dt}\hfill & \text{coalescence,}\hfill \end{array}$$
(36)
$$C(z_1,z_2;t)=\{\begin{array}{cc}1\hfill & (z_1,z_2)=(x_k,y_k),\hfill \\ \mathrm{erfc}(\xi _{lk})\hfill & (z_1,z_2)=(x_k,x_l),\hfill \\ e^{\xi _{lk}^2}\hfill & (z_1,z_2)=(x_k,y_l),\hfill \\ e^{\xi _{lk}^2}\hfill & (z_1,z_2)=(y_k,x_l),\hfill \\ \sqrt{\pi }\xi _{lk}e^{\xi _{lk}^2}\hfill & (z_1,z_2)=(y_k,y_l),\hfill \end{array}$$
(37)
and we used the notation $`\xi _{lk}=(x_lx_k)/\sqrt{8Dt}`$. For example, for $`n=2,3`$, we get the long-time asymptotic expressions:

$$\frac{\rho _2(x_1,x_2;t)}{\rho ^2}=1e^{2\xi _{21}^2}+\sqrt{\pi }\xi _{21}e^{\xi _{21}^2}\mathrm{erfc}(\xi _{21}),$$
(39)
$`{\displaystyle \frac{\rho _3(x_1,x_2,x_3;t)}{\rho ^3}}`$ $`=`$ $`1e^{2\xi _{21}^2}e^{2\xi _{32}^2}e^{2\xi _{31}^2}+2e^{\xi _{21}^2\xi _{32}^2\xi _{31}^2}`$ (43)
$`+\sqrt{\pi }\xi _{21}(e^{\xi _{21}^2}e^{\xi _{32}^2\xi _{31}^2})\mathrm{erfc}(\xi _{21})`$
$`+\sqrt{\pi }\xi _{32}(e^{\xi _{32}^2}e^{\xi _{21}^2\xi _{31}^2})\mathrm{erfc}(\xi _{32})`$
$`+\sqrt{\pi }\xi _{31}(e^{\xi _{31}^2}e^{\xi _{21}^2\xi _{32}^2})\mathrm{erfc}(\xi _{31}).`$
In summary, we have confirmed the fact that the infinite hierarchies of $`n`$-point correlation functions for coalescence and annihilation are identical, using the method of parity intervals. The simplicity of our approach allowed us to obtain explicit expressions for the long-time asymptotic limit, given in Eqs. (33). We note that our results are not restricted to long times. Indeed, for the case of a random distribution of particles, such that the initial concentration for annihilation is half that of coalescence, the identity holds at all times. In this case, explicit expressions for the correlation functions (valid at all times) can be obtained by using the full solution of Eqs. (15), (15),

$$E(x,y;t)=\mathrm{erfc}(\frac{yx}{\sqrt{8Dt}})\frac{1}{2}e^{2Dc_0^2t}\{e^{c_0(yx)}[1\mathrm{erf}(\frac{yx+4Dc_0t}{\sqrt{8Dt}})]e^{c_0(yx)}[1+\mathrm{erf}(\frac{yx4Dc_0t}{\sqrt{8Dt}})]\},$$
(44)
instead of the asymptotic expression of Eq. (33).

Remarkably, the particle distributions in coalescence and annihilation differ, despite the correspondence of the $`n`$-point correlation functions. The probability density function, $`p(x)`$, for the distance $`x`$ between two neighboring particles illustrates this difference. For large $`x`$, $`p(x)e^{x^2}`$ for coalescence, while $`p(x)e^x`$ for annihilation. Evidently, the complete hierarchy of $`n`$-point correlation functions is not sufficient to determine an infinite-particle system uniquely, and $`p(x)`$ cannot be computed from a knowledge of the $`\rho _n`$. A study of $`p(x)`$ and $`\{\rho _n\}`$ in finite systems might illuminate this curious phenomenon.

###### Acknowledgements.

Funding of this project by the NSF (PHY-9820569) is gratefully acknowledged.


format: text

text:
# Do Hadronic Charge Exchange Reactions Measure Electroweak 𝐿=1 Strength?

## I Introduction

Charge exchange reactions ($`A_Z,A_{Z\pm 1}`$) induced by hadronic projectiles are a powerful tool for probing spin-isospin degrees of freedom in nuclei . The spin-isospin parts of the operators that mediate charge exchange reactions such as ($`p,n`$) are the same as those involved in the corresponding processes induced by electromagnetic and weak interactions. As a result, the matrix elements that describe hadronic charge exchange reactions are closely related to those that describe the rates of $`\beta `$ decay or the cross sections of reactions induced by neutrinos. It would be fortunate if this relationship were quantitatively accurate, since it is often difficult to study the leptonic processes directly. For example, the range of excitation energy kinematically accessible in a $`\beta `$ decay transition does not encompass the majority of the allowed (Gamow-Teller) strength and the experimental study of neutrino induced reactions is difficult.

A promising direction of future activity is to determine leptonic strengths for otherwise inaccessible nuclides by studying charge exchange reactions using radioactive (secondary) beams in inverse kinematics . This would provide nuclear properties important for problems of nuclear physics, particle physics, astrophysics and cosmology . One could clarify the relationship between the spatial properties of nuclear halo systems and the nature of soft multipole modes . One could also determine the strength of neutrino-nucleus interactions needed to describe the chemical evolution of the Universe, especially the abundances of the light elements  and the products of $`r`$-process nucleosynthesis , and to calibrate terrestrial detectors of supernova neutrinos .

However, it is not obvious a priori that the correlation of charge exchange and leptonic matrix elements is sufficiently close for this purpose. In contrast to leptonic processes, hadronic reactions involve operators that have an additional radial dependence and are subject to distortion by complex nuclear potentials; the medium renormalization of effective operators and the contributions of multi-step processes introduce additional uncertainties. Therefore it was an important advance to establish that there is an approximate proportionality between the cross section of charge exchange reactions at very forward angles leading to the Gamow-Teller (GT) excitations (with transferred $`T=1,L=0,J=S=1`$) and the transition strength $`B`$(GT) determined by intrinsic nuclear matrix elements . The proportionality has been confirmed for strong transitions in a variety of nucleon and nucleus induced charge exchange reactions; a more detailed analysis is needed for weak GT processes . The model-independent character of the relation between the $`L=0`$ cross section for reactions induced by 12C projectiles and the GT strength was clarified by a theoretical analysis  based on a sensitivity function which identified the important part of the target transition density in momentum space.

In contrast to GT transitions, the first forbidden matrix elements of weak processes explicitly include orbital degrees of freedom. The corresponding nuclear response in the $`L=1`$ channel is associated with the states forming the spin-dipole and giant-dipole charge-exchange resonances (SDR and GDR). This excitation was discovered  and studied on different targets  mostly using ($`p,n`$) reactions. Recently the energy splitting of the $`L=1`$ charge exchange resonances (GDR and SDR) was determined . There is very little information about the quantitative relationship of charge exchange cross sections and leptonic strength for these excitations. Here we take a first step in providing this information by studying the relationship between first forbidden strength and ($`p,n`$) cross sections, both calculated from the same wave functions. We show that, at the same level of accuracy as for the GT case, for strong transitions one can expect an approximate proportionality between the observed cross sections of charge exchange reactions populating spin-dipole states and the corresponding nuclear transition probabilities. We follow the general approach that was successfully applied to GT excitations by Osterfeld et al. , extended to describe $`L=1`$ transitions and the effects of the real part of the optical potential.

## II Theoretical framework

In order to understand the relationship between the charge exchange cross section and the nuclear response strength, we need to examine the effects specific to the excitation of the SDR. For this purpose we consider the influence of distortion in a simple eikonal approximation (EA). In many cases, even at rather low energy, the EA gives a good qualitative description of the reaction cross section in the SDR region. We obtain the relevant wave functions and transition densities from the shell model.

The strength $`B_J`$ of the SDR in the long wavelength limit and the transition form-factors $`F_J(q)`$ are calculated in terms of the elements of the single-particle density matrix of the given transition $`if`$ from the ground state, $`\rho _{fi}(\nu ,\nu ^{})=f|a_\nu ^{}^{}a_\nu |i`$,

$$B_J(if)=\left|\underset{\nu \nu ^{}}{}\rho _{fi}(\nu ,\nu ^{})(\nu ^{}rO_J\nu )\right|^2,$$
(1)
$$F_J(q)=\underset{\nu \nu ^{}}{}\rho _{fi}(\nu ,\nu ^{})(\nu ^{}j_1(qr)O_J\nu ),$$
(2)
where $`(\nu ^{}j_1(qr)O_J\nu )`$ is the product of the radial matrix element of the spherical Bessel function and the reduced matrix element of the charge exchange spin-dipole operator $`O_{JM}=\{\sigma Y_1(\widehat{𝐫})\}_{JM}\tau ^\pm `$. In the limit of low momentum transfer, $`qR1`$, the squared transition form-factor (2) is related directly to the strength of the SDR,

$$F_J^2(q)\frac{q^2}{9}B_J(if).$$
(3)
In the case of the GT resonance a direct proportionality between the experimentally measured cross sections of charge exchange reactions at very forward angles and nuclear matrix elements was confirmed by a number of studies  at a level of accuracy of 10-15%. From the viewpoint of the underlying physics, this important result is based on a single-step mechanism for the process, a simple bare operator which does not include orbital degrees of freedom, and the dominance of the central spin-isospin interaction $`V_{\sigma \tau }`$ over a broad range of energies. It is a priori unclear whether these features pertain to the SDR case ($`\mathrm{\Delta }L=1`$). At the maximum of the differential cross section for the SDR, $`qR1`$. The $`q`$-dependence of the hadronic operator and the effects of tensor forces may show up at the larger $`q`$. As a result, there may not be a simple relationship between the cross section and the nuclear response strength.

In the SDR case, the direct part of the $`(p,n)`$ reaction amplitude (the exchange part will be discussed later) is

$$T_{fi}^{dir}=d^3𝐫\chi _f^{()}(𝐤_f,𝐫)\frac{d^3𝐪^{}}{(2\pi )^3}F_J(q^{})V_{JM}(𝐪^{})\mathrm{exp}(ı𝐪^{}𝐫)\chi _i^{(+)}(𝐤_i,𝐫),$$
(4)
where the effective operators for the channels with angular momenta $`J^\pi =0^{},1^{}`$ and $`2^{}`$ contain contributions from central and tensor forces,

$$V_{JM}(𝐪^{})=V_{JM}^c(𝐪^{})+V_{JM}^t(𝐪^{}),$$
(5)
$$V_{JM}^c(𝐪^{})=4\pi ı\sqrt{\frac{2}{2J+1}}t_{\sigma \tau }^c(q^{})\{\sigma Y_1^{}(\widehat{𝐪}^{})\}_{JM},$$
(6)
$$V_{JM}^t(𝐪^{})=4\pi ı\sqrt{\frac{2}{2J+1}}t_\tau ^t(q^{})[\{\sigma Y_1^{}(\widehat{𝐪}^{})\}_{JM}3(𝝈\widehat{𝐪}^{})\{\widehat{𝐪}^{}Y_1^{}(\widehat{𝐪}^{})\}_{JM}].$$
(7)
In Eqs. (6) and (7), $`t_{\sigma \tau }^c(q^{})`$ and $`t_\tau ^t(q^{})`$ are, respectively, the central and tensor components of the nucleon-nucleon $`t`$-matrix, see Franey and Love , that are responsible for spin-isospin transfer. The excitation of different $`J`$-components of the SDR proceeds via different combinations of the amplitudes of the nucleon-nucleon effective interaction. For the $`0^{}`$ part, the tensor interaction can be combined with the central one using the relations

$$\{\widehat{𝐪}^{}Y_1^{}(\widehat{𝐪}^{})\}_{00}=\frac{1}{\sqrt{4\pi }},\{\sigma Y_1^{}(\widehat{𝐪}^{})\}_{00}=\frac{(𝝈\widehat{𝐪}^{})}{\sqrt{4\pi }}.$$
(8)
These identities produce the combination $`t^l(q)=t_{\sigma \tau }^c(q)2t_\tau ^t(q)`$ which is just the spin-longitudinal component of the nucleon-nucleon $`t`$-matrix. For the $`1^{}`$ part, the operator $`\{\widehat{𝐪}^{}Y_1^{}(\widehat{𝐪}^{})\}_{1M}=0`$, and the amplitude in Eq. (4) becomes proportional to $`t^{tr}(q)=t_{\sigma \tau }^c(q)+t_\tau ^t(q)`$, which is the spin-transverse component of the nucleon-nucleon $`t`$-matrix. For the $`2^{}`$ part, both components contribute, and the amplitude in Eq. (4) has a more complicated form.

The functions $`\chi ^{()}`$ and $`\chi ^{(+)}`$ in the amplitude (4) are optical-model wave functions describing the motion of the initial proton and final neutron in the optical potential of the target nucleus. To disentangle the nuclear transition form-factor from the observed cross section, one needs to unravel the intrinsic nuclear dynamics masked by the distorted waves $`\chi ^{(\pm )}`$.

## III Distortion factor

The effective operators (6) and (7) in the reaction amplitude (4) are evaluated at the value $`𝐪^{}`$ of the local momentum transfer that corresponds to the charge exchange event. However, because of the distortion by the optical potential, $`𝐪^{}`$ does not coincide with the asymptotic momentum transfer $`𝐪=𝐤_i𝐤_f`$. In the absence of distortion (the plane wave approximation) we would have

$$[\chi _f^{()}(𝐤_f,𝐫)\chi _i^{(+)}(𝐤_i,𝐫)]_{PW}=\mathrm{exp}(ı𝐪𝐫)$$
(9)
so that the integration over $`𝐫`$ could be performed explicitly resulting in $`\delta (𝐪𝐪^{})`$. For distorted waves this is no longer true.

In the EA, the product of two optical-model wave functions in Eq. (9) can be estimated by

$$\chi _f^{()}(𝐤_f,𝐫)\chi _i^{(+)}(𝐤_i,𝐫)=\mathrm{exp}(ı𝐪𝐫)D(𝐫_{}),$$
(10)
where the distortion factor $`D(𝐫_{})`$ is defined by

$$D(𝐫_{})=\mathrm{exp}\left[\frac{ı}{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\frac{dz}{v}U_{opt}(z,𝐫_{})\right].$$
(11)
In the spirit of the eikonal approximation, the longitudinal momentum is still preserved, while the distortion is effective in the plane perpendicular to the trajectory. In Eq. (11), the optical potential is different in the initial and final channels. In first order we account for this difference by assuming a fast single-step process which leads to

$$\frac{U_{opt}}{v}=\frac{1}{2}\left(\frac{U_{opt}^i}{v_i}+\frac{U_{opt}^f}{v_f}\right).$$
(12)
For a square well potential of depth $`U_0`$ with a sharp boundary at $`r=R`$ the distortion factor can be calculated analytically:

$$D(r)=\mathrm{exp}\left(ı\frac{2U_0}{\mathrm{}v}\sqrt{R^2r^2}\right),$$
(13)
and $`D(r)=1`$ for $`r>R`$. This approximation might be insufficient in the region of minimum of the cross section where the details of the potential shape are essential. However, near the maximum, which is our region of interest, the cross section is insensitive to the diffuseness of the optical-model potential. We have checked this point by varying the diffuseness in a DWIA calculation.

The exchange part of the reaction amplitude was estimated in the “fixed $`𝐐`$” approximation . This approximation has been found  to improve with increasing energy and multipolarity of the excitation for even-state forces and mainly $`S=0`$ natural parity transitions. We need to deal with $`S=1`$ transitions and include tensor forces. Recently, the fixed $`𝐐`$ approximation for exchange was discussed in detail in a review article ; it is more accurate for central interactions than for the tensor interaction so that the latter is of most concern. Fortunately, the approximation for tensor forces is best when $`Q`$ is not far from $`k_i`$, which is the case for the relatively high energies considered here ($`k_i`$ is greater than the Fermi momentum in the nucleus) and for small angles. It is also relatively less important for the transitions we consider where there are large central and direct tensor contributions . Baker et al.  note that for the intermediate energy ($`p,p^{}`$) isovector excitation of $`1^{}`$ and $`2^{}`$ states in 40Ca, the cross sections obtained with the fixed $`𝐐`$ approximation were usually within 30% of the results of calculations done with exchange treated exactly. As shown in Fig. 1, our eikonal calculations with the fixed $`𝐐`$ approximation are within 20% of DWBA calculations with exact exchange, consistent with the ($`p,p^{}`$) result. This level of agreement seems adequate for our purposes: the overall effect of the approximation is simply to change slightly the scale of the cross sections for a given $`J`$ from that shown in Fig. 2. If the effect were different for different transitions, it would show up in the overall scatter of the points about the average line which is small for strong transitions.

In the fixed $`𝐐`$ approximation  the exchange momentum $`𝐐`$ coincides in the laboratory frame with the initial momentum $`𝐤_i`$. For the exchange amplitude we then obtain

$$T_{fi}^{ex}=\sqrt{2}\left[(\stackrel{~}{t}_{\sigma \tau }^c(k_i)+\stackrel{~}{t}_\tau ^t(k_i))𝝈3\stackrel{~}{t}_\tau ^t(k_i)(𝝈\widehat{𝐤}_i)\widehat{𝐤}_i\right]f|𝐎^{ex}(𝐪)|i,$$
(14)
where the central $`\stackrel{~}{t}_{\sigma \tau }^c`$, and tensor, $`\stackrel{~}{t}_\tau ^t`$, interactions are defined by Franey and Love . In Eq. (14) the effective exchange operator

$$𝐎^{ex}(𝐪)=\underset{j}{}\mathrm{exp}(ı𝐪𝐫_j)D(𝐫_j)𝝈_j\tau _j^{}$$
(15)
includes the distortion factors $`D(r)`$ specific for each nucleon inside the nuclear matrix element.

## IV Sensitivity function

In Eq. (4) the integration over $`𝐫`$ is not well defined at large distances. It is convenient to single out a no-distortion contribution proportional to $`\delta (𝐪𝐪^{})`$ by using the decomposition $`D(r)1+[D(r)1]`$. The first term describes the plane wave contribution and the second the effects of distortion. Since $`\left|D(r)\right|1`$, the distortion term reduces the plane wave contribution. A convenient form of Eq. (4) can be obtained for transitions to $`0^{}`$ and $`1^{}`$ states by writing it as

$$T_{fi}^{dir}=\frac{\{\sigma \mathrm{T}^{(J)}\}_{JM}}{\sqrt{2J+1}},$$
(16)
where

$$\mathrm{T}_m^{(J)}=\mathrm{T}_m^J(PW)+_0^{\mathrm{}}𝑑q^{}S_m^J(q,q^{})F_J(q^{})$$
(17)
is the amplitude describing the excitation of the SDR with the longitudinal, $`m=0`$, or transverse, $`m=\pm 1`$, relative proton-neutron spatial oscillations. In Eq. (17), $`\mathrm{T}_m^{(J)}(PW)`$ is the plane wave contribution, and we have introduced the sensitivity function

$$S_m^J(q,q^{})=\frac{2\sqrt{2}}{\pi }q^2d^3𝐫\mathrm{exp}(ı𝐪𝐫)(D(𝐫_{})1)Y_{1m}(\widehat{𝐫})j_1(q^{}r)$$

$$\times \{\begin{array}{c}t^l(q^{})\mathrm{for}J=0,\\ t^{tr}(q^{})\mathrm{for}J=1\end{array}.$$
(18)
The sensitivity function in Eq. (18) characterizes the range of $`q^{}`$ which contribute importantly to the charge exchange cross section for a given asymptotic momentum transfer $`q`$.

## V Example: 12C$`(p,n)^{12}`$N reaction

As an example of application of the method we performed numerical calculations for the 12C$`(p,n)^{12}`$N reaction to compare with experimental data  for the excitation of spin-dipole states at a proton energy of 135 MeV.

For 12C, with the optical-model potential of Ref. , $`D(r)`$ varies smoothly inside the nucleus. Near the surface it changes rapidly from its value at the center $`D(0)0.5`$, to the value of 1. It is then a good approximation to write the exchange matrix element in Eq. (14) as

$$f|𝐎^{ex}(𝐪)|iD(r_0)f|\underset{j}{}\mathrm{exp}(ı𝐪𝐫_j)𝝈_j\tau _j^{}|i.$$
(19)
The result is not very sensitive to a particular choice of the reference point $`r_0`$; we used $`r_0=0`$. This approximation underestimates the exchange part. Near the maximum of the cross section its contribution is not significant. It becomes important at large angles where the difference of distortion along different trajectories is noticeable; however, the cross section at large angles is small. We mention parenthetically that, opposite to the results of  for even-state forces, the exchange contribution reduces the cross section. This is related to our inclusion of odd-state forces whose net effect is small after exchange is accounted for.

In our calculations, the wave functions and transition densities for the spin-dipole states were obtained using a harmonic oscillator basis including the orbitals of $`p,sd`$ and $`pf`$ shells that form the $`3\mathrm{}\omega `$ model space necessary for the description of the $`L=1`$ excitations. The calculations were performed with the WBN residual interaction and a harmonic oscillator parameter of $`1.64(A/A1)^{1/2}`$ fm . The cross sections for the 12C($`p,n`$)12N reaction leading to the $`1^{}`$ state at $`E_x=1.8`$MeV and the $`2^{}`$ state at $`E_x=4.3`$ MeV were calculated as the sum of the direct and exchange amplitudes, Eqs. (7) and (14). The results are shown in Fig. 1 together with data from Ref. . For comparison, a calculation with the DW81 code is also presented. The calculations give similar cross section shapes near the maximum and significantly overestimate the magnitude of the cross section. The results are very similar for other excited states. They are also similar to the distorted wave results obtained for the same transition in Ref.  using a $`1\mathrm{}\omega `$ model space and the MK interaction.

The systematics of the cross sections at their maximum divided by the calculated $`\beta `$ decay strengths are shown in Fig. 2 for $`0^{},\mathrm{\hspace{0.17em}1}^{}`$ and $`2^{}`$ states. As seen from Fig. 2, there is an approximate proportionality between the cross section at the maximum and the spin-dipole strength, accurate to within 10-15%, for states with strength $`B_J>0.1\mathrm{fm}^2`$. This is the same level of proportionality as for GT ($`L=0`$) excitations at very forward angles. One may ask whether the validity of this conclusion is affected by the poor agreement in the magnitude of the cross section for the 1.8 MeV $`1^{}`$ state. We would argue that this is not the case: since the wave functions are sufficiently complex, they provide a reasonable sample of possible behavior with respect to the operators involved. Furthermore, one might expect proportionality to fail for such weak transitions.

## VI Discussion

It is not clear a priori that the high degree of proportionality shown in Fig. 2 should occur. The cross section involves an integral of the transition form factor over a range of $`q^{}`$ while the value of $`B_J`$ is determined by evaluating the form factor at $`q^{}0`$. To examine what leads to the observed proportionality, we return to Eq. (17). Two factors determine the ($`p,n`$) cross section: the transition form-factor $`F_J(q)`$ and the sensitivity function $`S_m^J(q,q^{})`$. In Fig. 3 we show the transition form-factors for different $`1^{}`$ states normalized to the same maximum value in order to compare their shapes. The shapes are very similar near the maximum but differ at higher momentum transfer $`q`$. If the region of high $`q^{}`$ does not contribute significantly in the integration over $`q^{}`$ in Eq. (17), the integrals for different form-factors will be proportional.

Samples of the imaginary parts of the sensitivity functions are shown in Figs. 4 and 5; the real parts have very similar shapes and are typically a factor of two smaller in magnitude. A general remark should be made about the $`q^{}`$-dependence at small $`q^{}`$. Since $`D(𝐫_{})`$ does not depend on the longitudinal coordinate $`z`$, the $`z`$-component of the local momentum transfer $`𝐪^{}`$ must coincide with the $`z`$-component of the asymptotic momentum transfer $`𝐪`$. When the absolute value of $`q^{}`$ is smaller than $`q_z`$, this condition cannot be fulfilled at any angles of $`𝐪^{}`$, and the sensitivity function must be equal to zero.

As noted above, the sensitivity function for $`J=0`$ is proportional to the spin-longitudinal component of the effective interaction, and that for $`J=1`$ to the spin-transverse one. We, therefore, expect a different $`q^{}`$-dependence reflecting the different behavior of $`t^l(q^{})`$ and $`t^{tr}(q^{})`$. For $`0^{}`$ states the projection $`m=0`$ dominates, corresponding to the spin-longitudinal behavior of the reaction amplitude for $`J=0`$; the sensitivity function for $`m=1`$ is smaller by an order of magnitude. At the small scattering angle corresponding to the peak of the cross section, $`\theta =4.3^{}`$, the momentum transfer $`𝐪`$ is almost parallel to the initial proton momentum $`𝐤_i`$, thus enhancing the $`m=0`$ component. For $`J=1`$ the picture is different, as is seen in Fig. 5. Projections $`m=0`$ and $`m=\pm 1`$ give comparable contributions.

Given the nature of the sensitivity functions it is clear why the cross sections and $`B_J`$ are closely proportional. For both $`J=0`$ and $`J=1`$ the main contribution comes from the peak region where the transition form factors have the same shape, leading to the observed proportionality.

In summary, our results imply that there will be an approximate proportionality of the observed cross section at the maximum of the charge exchange reaction exciting spin-dipole modes and the leptonic strength. This supports the possibility of using such reactions for extracting leptonic strengths of astrophysical interest. Having established here the basic apparatus to examine this issue, it will next be important to consider transition densities for heavier nuclides, so as to determine whether their shapes are similar enough that cross sections and $`B_J`$ strengths will be proportional. It will also be important to examine the nature of the sensitivity functions for heavier nuclei, to ascertain whether they remain concentrated in a relatively small range of $`q^{}`$ where the transition densities are similar.

ACKNOWLEDGMENTS

The authors are thankful to W. G. Love for a helpful discussion. The shell model information necessary for the calculations was kindly provided by B. A. Brown. V.F.D. is thankful to T. Nakamura for assistance and discussions and to the NSCL for the hospitality and support. This work was supported by the National Science Foundation.


text:
# Anomalous Rashba spin splitting in two-dimensional hole systems

## Abstract

It has long been assumed that the inversion asymmetry-induced Rashba spin splitting in two-dimensional (2D) systems at zero magnetic field is proportional to the electric field that characterizes the inversion asymmetry of the confining potential. Here we demonstrate, both theoretically and experimentally, that 2D heavy hole systems in accumulation layer-like single heterostructures show the opposite behavior, i.e., a decreasing, but nonzero electric field results in an increasing Rashba coefficient.

Spin degeneracy of electron and hole states in a solid stems from the inversion symmetry in space and time. If the spatial inversion symmetry is broken, then there is a splitting of the single particle states even at magnetic field $`B=0`$ (Ref. ). In quasi two-dimensional (2D) semiconductor structures, the bulk inversion asymmetry (BIA) of the underlying crystal structure (e.g., a zinc blende structure), and the structure inversion asymmetry (SIA) of the confining potential contribute to the $`B=0`$ spin splitting.  While BIA is fixed, the so-called Rashba spin splitting  due to SIA can be tuned by means of external gates that change the electric field $`𝐄`$ in the sample.  The $`B=0`$ spin splitting is of significant current interest both because of its fundamental importance as well as its possible device applications.

For many years it has been assumed that the Rashba spin splitting in 2D systems is proportional to the electric field that characterizes the inversion asymmetry of the confining potential.  In single heterostructures, where SIA is the dominant source of spin splitting, the electric field is determined by the density-dependent self-consistent potential.  One would thus expect that the spin splitting decreases with density, although for 2D electron systems this effect may partly be compensated by many-particle effects that tend to increase the Rashba spin splitting for low densities.  Here we demonstrate, both theoretically and experimentally, that 2D heavy hole systems in accumulation layer-like single heterostructures show the opposite behavior, namely, a decreasing, but nonzero electric field results in an increasing Rashba coefficient. Contrary to electrons, however, exchange-correlation effects in the low-density regime decrease the spin splitting. We show that this surprising result is a consequence of heavy hole–light hole (HH-LH) coupling in 2D hole systems. We obtain good qualitative agreement between calculated and measured spin splittings in 2D hole systems in GaAs heterostructures where the density and the spin splitting are varied by means of an external gate. Our results are applicable to many systems, as most III-V semiconductors have essentially the same band structure that is underlying our investigation.

To lowest order of the wave vector $`𝐤`$ and electric field $`𝐄`$, the SIA spin splitting of electron states in the $`\mathrm{\Gamma }_6^c`$ conduction band is given by the Rashba term

$$H_{6c}^{\mathrm{SO}}=\alpha 𝐤\times 𝐄𝝈.$$
(1)
Here $`𝝈=(\sigma _x,\sigma _y,\sigma _z)`$ denotes the Pauli spin matrices and $`\alpha `$ is a material-specific prefactor.  We assume $`𝐄=(0,0,E_z)`$. Treating the off-diagonal $`𝐤𝐩`$ coupling between electron and hole states by third order Löwdin perturbation theory,  we obtain for the Rashba coefficient $`\alpha _\lambda `$ of the lowest electron subband $`\lambda =1`$

$$\alpha _1=eP^2a\left(\frac{1}{\mathrm{\Delta }_{11}^{cl}}\frac{1}{\mathrm{\Delta }_{12}^{cl}}\frac{1}{\mathrm{\Delta }_{11}^{cs}}\frac{1}{\mathrm{\Delta }_{12}^{cs}}\right),$$
(2)
where $`P`$ is Kane’s momentum matrix element  and $`\mathrm{\Delta }_{\lambda \lambda ^{}}^{\nu \nu ^{}}_\lambda ^\nu _\lambda ^{}^\nu ^{}`$ with $`_\lambda ^c`$, $`_\lambda ^h`$, $`_\lambda ^l`$, and $`_\lambda ^s`$ the energy of the $`\lambda `$th electron, HH, LH, and split-off subband, respectively.  The numerical prefactor $`a`$ depends on the geometry of the confining quantum well (QW). In an infinitely deep rectangular QW we have $`a=256/(81\pi ^2)`$. According to Eq. (1) we obtain a spin splitting $`\pm \alpha E_zk_{}`$ of the subband dispersion $`(𝐤_{})`$ that is proportional to the electric field $`E_z`$ and is linear in the in-plane wave vector $`𝐤_{}=(k_x,k_y,0)`$. A detailed analysis reveals that spin splitting of electron states depends on the electric field $`E_v`$ in the valence band that differs from the electric field $`E_c`$ in the conduction band by the contributions of the interfaces.  However, the important point here is that in a single heterostructure both $`E_v`$ and $`E_c`$ are determined by the self-consistent Hartree potential.

For hole systems in the $`\mathrm{\Gamma }_8^v`$ valence band (point group $`T_d`$), the dominant contribution to Rashba spin splitting is given by the term

$$H_{8v}^{\mathrm{SO}}=\beta 𝐤\times 𝐄𝐉,$$
(3)
where $`\beta `$ is a system-dependent prefactor and $`𝐉=(J_x,J_y,J_z)`$ are the angular momentum matrices for $`j=3/2`$. We neglect here the small corrections in $`H_{8v}^{\mathrm{SO}}`$ due to the $`𝐤𝐩`$ coupling to remote bands such as the higher $`\mathrm{\Gamma }_8^c`$ and $`\mathrm{\Gamma }_7^c`$ conduction bands.  Quantum confinement reduces the symmetry from the cubic point group $`T_d`$ to $`D_{2d}`$. Therefore, the four-fold degeneracy of the $`\mathrm{\Gamma }_8^v`$ band ($`T_d`$) is lifted and we obtain two-fold degenerate subspaces transforming according to $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_6`$ of the point group $`D_{2d}`$ (Refs. ). The irreducible representation $`\mathrm{\Gamma }_7`$ corresponds to the LH states ($`z`$ component of angular momentum $`j_z=\pm 1/2`$) whereas $`\mathrm{\Gamma }_6`$ corresponds to the HH states ($`j_z=\pm 3/2`$). To lowest order in $`k`$ the effective Rashba Hamiltonian for the LH states is the same as in Eq. (1) for electron states.  For HH states, however, spin splitting is mediated by a coupling to the LH states so that, to lowest order, spin splitting of HH states is of third order in $`𝐤`$ (Ref. ). Neglecting anisotropic corrections in the Hamiltonian we have

$$H_h^{\mathrm{SO}}=\beta ^hE_z(\sigma _+k_{}^3+\sigma _{}k_+^3),$$
(4)
with $`\sigma _\pm =1/2(\sigma _x\pm i\sigma _y)`$ and $`k_\pm =k_x\pm ik_y`$. Treating the off-diagonal HH-LH coupling by third order Löwdin perturbation theory  we obtain for the Rashba coefficient $`\beta _\lambda ^h`$ of the lowest HH subband $`\lambda =1`$

$$\beta _1^h=ia\gamma _3(\gamma _2+\gamma _3)\left[\frac{1}{\mathrm{\Delta }_{11}^{hl}}\left(\frac{1}{\mathrm{\Delta }_{12}^{hl}}\frac{1}{\mathrm{\Delta }_{12}^{hh}}\right)+\frac{1}{\mathrm{\Delta }_{12}^{hl}\mathrm{\Delta }_{12}^{hh}}\right]$$
(5)
where $`\gamma _2`$ and $`\gamma _3`$ are the Luttinger parameters.  In an infinitely deep rectangular QW we have $`a=64/(9\pi ^2)`$. We see from Eq. (5) that the Rashba spin splitting of HH states depends not only on the electric field $`E_z`$ but also on the separation between the HH and LH subbands. A decreasing separation gives rise to an increasing Rashba coefficient $`\beta _1^h`$. The factor $`\gamma _3(\gamma _2+\gamma _3)`$ in Eq. (5) refers to a quantum structure grown in the crystallographic direction . The expressions for other growth directions are similar, but the other terms in Eq. (5) remain unchanged. We remark that for typical hole densities only the lowest HH subband is occupied.

The electric field $`E_z`$ that enters into the Rashba Hamiltonian depends on the charges in the system. We will show now that accumulation layer-like single heterostructures behave rather differently with respect to changes of the 2D charge density as compared to other quasi 2D semiconductor structures.

In a rectangular QW, a small density $`N`$ and a small asymmetry imply that the properties of the system are controlled by the effective potential steps at the interfaces, i.e., changes in $`N`$ or $`E_z`$ have a minor effect in this regime. In an inversion layer-like heterostructure, we always have a band bending of the order of the fundamental gap so that, for small densities, the Hartree potential and $`E_z`$ are determined by the space charges due to the given concentration of ionized majority impurities in the system. For accumulation layer-like systems, on the other hand, it was shown by Stern,  that the space charge layer is controlled by the much smaller concentration of minority impurities in the system. Thus, even for a small 2D density, the dominant contribution to the Hartree potential stems from the charges in the 2D system itself. Therefore, over a wide range of densities $`N`$, the electric field $`E_z`$ is proportional to $`N`$. In single heterostructures, the subband separations are approximately proportional to $`E_z`$. Using the triangular well approximation  we have, for the subband energies $`_\lambda ^\nu `$ measured from the corresponding bulk band edge, $`_\lambda ^\nu E_z^{2/3}`$ which implies $`_\lambda ^\nu N^{2/3}`$ and $`\beta _\lambda ^hN^{4/3}`$. Therefore, we can expect from Eqs. (4) and (5) that accumulation layer-like 2D HH systems show a Rashba spin splitting that increases when $`N`$ and $`E_z`$ are reduced. On the other hand, the coefficient (2) is essentially independent of $`N`$ and $`E_z`$ because the energy gaps $`\mathrm{\Delta }_{\lambda \lambda ^{}}^{\nu \nu ^{}}`$ are always of the order of the fundamental gap.

In order to validate these qualitative arguments we present next the results of realistic, fully self-consistent subband calculations.  We use an $`8\times 8`$ multiband Hamiltonian  that includes the lowest conduction band $`\mathrm{\Gamma }_6^c`$, the topmost valence band $`\mathrm{\Gamma }_8^v`$ and the split-off valence band $`\mathrm{\Gamma }_7^v`$. The simpler $`4\times 4`$ Luttinger Hamiltonian,  taking into account only the band $`\mathrm{\Gamma }_8^v`$, gives essentially the same results. We have checked that higher conduction bands have a minor influence. Many-particle effects are taken into account based on a density-functional approach.  From these calculations we obtain the difference $`\mathrm{\Delta }N=N_+N_{}`$ between the spin subband densities $`N_\pm `$ as a function of the total density $`N=N_++N_{}`$.

In Fig. 1(a) we show $`\mathrm{\Delta }N/N`$ calculated as a function of $`N`$ for a 2D HH system in the accumulation layer of a GaAs-Al0.5Ga0.5As single heterostructure on a (001) GaAs substrate.  From $`N=5\times 10^{11}`$ cm-2 to $`1\times 10^{10}`$ cm-2 the parameter $`r_s`$, the Coulomb energy to Fermi energy ratio, increases from 4.3 to 17. Therefore, one can expect that many-particle effects are quite important in this regime of densities $`N`$. Indeed, we find that taking into account exchange-correlation (solid lines) reduces $`\mathrm{\Delta }N/N`$ as opposed to a calculation without exchange-correlation (dashed lines). This behavior, which is opposite to 2D electron systems,  can be traced back to the fact that exchange-correlation increases the subband spacings  so that the Rashba coefficient $`\beta _1^h`$ is reduced, in agreement with Eq. (5).

It is convenient to characterize our numerical results in terms of an effective Rashba coefficient $`\beta _1^hE_z`$ (Ref. ). Assuming that the spin-split HH subband dispersion is approximately of the form $`_\pm ^h(k_{})=\mu _hk_{}^2\pm \beta _1^hE_zk_{}^3`$, where $`\mu _h`$ (times $`2/\mathrm{}^2`$) is the reciprocal effective mass, we have

$$\beta _1^hE_z=\sqrt{\frac{2}{\pi }}\mu _h\frac{N\left(\sqrt{N+\mathrm{\Delta }N}\sqrt{N\mathrm{\Delta }N}\right)+\mathrm{\Delta }N\left(\sqrt{N+\mathrm{\Delta }N}+\sqrt{N\mathrm{\Delta }N}\right)}{6N^2+2\mathrm{\Delta }N^2}.$$
(6)
Figure 1(b) shows that $`\beta _1^hE_z/\mu _h`$ increases when $`N`$ is reduced.

For comparison, we have calculated $`\mathrm{\Delta }N/N`$ for a 2D electron system in the accumulation layer of a Ga0.47In0.53As-Al0.47In0.53As single heterostructure  \[Fig. 2(a)\]. Here spin splitting is given by Eq. (1). Therefore, the spin-split subband dispersion is approximately of the form $`_\pm ^c(k_{})=\mu _\mathrm{c}k_{}^2\pm \alpha _1E_zk_{}`$, and we obtain similarly to Eq. (6) (Ref. )

$$\alpha _1E_z=\sqrt{2\pi }\mu _c\left(\sqrt{N+\mathrm{\Delta }N}\sqrt{N\mathrm{\Delta }N}\right).$$
(7)
In Fig. 2(b) it can be seen that, in contrast to Fig. 1(b), the spin splitting coefficient $`\alpha _1E_z/\mu _c`$ decreases rapidly with decreasing $`N`$. We remark that unlike the HH system in Fig. 1, exchange correlation has only a weak influence on the electron system in Fig. 2 (Ref. ).

To further analyze our numerical results, we define an effective electric field $`E_z=_zV_H(z)`$, where $`V_H(z)`$ is the Hartree potential without the effective potential due to the position-dependent band edges. In an accumulation layer the contribution of the space charge layer to $`V_H(z)`$ is very small.  It follows then, by partial integration of the Poisson equation, that $`E_z=e/(2\epsilon \epsilon _0)N`$, where $`\epsilon `$ is the dielectric constant and the expectation value refers to the 2D charge density that gives rise to $`V_H(z)`$. Using these values for $`E_z`$ and $`\mu _c=89`$ eVÅ2 we obtain $`\alpha _134.3`$ eÅ2 independent of $`N`$, consistent with Eq. (2). This implies that in Fig. 2(b) the drastic change of $`\alpha _1E_z/\mu _c`$ merely reflects the change of the electric field $`E_z`$. On the other hand, the weak variation of $`\beta _1^hE_z/\mu _h`$ in Fig. 1(b) indicates that the “bare” Rashba coefficient $`\beta _1^h`$ increases by a factor of 250 when $`N`$ is lowered from $`5\times 10^{11}`$ to $`1\times 10^{10}`$ cm-2 (Ref. ). This is in good, qualitative agreement with the analytical model discussed above that predicts an increase of $`\beta _1^h`$ by a factor of $`50^{4/3}184`$. Note that for low densities the third order perturbation approach, that underlies Eqs. (2) and (5), breaks down because the subbands are merging together so that higher order corrections become important. These higher order terms are fully taken into account in our numerical calculations.  If the density $`N`$ is reduced below $`10^{10}`$ cm-2 the Hartree potential and spin splitting are ultimately controlled by the fixed concentration of minority impurities.

It is interesting to compare the spin splittings in accumulation layers with those in QW’s where $`E_z`$ is tuned externally, e.g., by means of gates.  The dotted lines in Figs. 1 and 2 show the calculated results for a $`200`$ Å wide rectangular QW where the external electric field $`E_z^{\mathrm{ext}}`$ was chosen according to $`E_z^{\mathrm{ext}}(N)=e/(2\epsilon \epsilon _0)N`$. In an electron system (Fig. 2) we obtain spin splittings very close to the results for the accumulation layer. In particular, we have $`\alpha _130.6`$ eÅ2 independent of $`N`$. Similarly, for a 2D HH system in a QW (Fig. 1) and $`N1\times 10^{11}`$ cm-2 we obtain $`\beta _1^h7.54\times 10^6`$ eÅ4. \[For larger densities higher order corrections in $`_\pm ^h(k_{})`$ become important. \] Since in QW’s the subband spacings are essentially determined by the QW width (i.e., are independent of $`N`$) this is consistent with Eq. (5). These calculations also indicate that for 2D HH systems in a QW, spin splitting becomes negligible in the regime of low densities,  which is due to the fact that spin splitting of $`_\pm ^h(k_{})`$ is proportional to $`k_{}^3`$. However, for 2D HH systems in single heterostructures, spin splitting can be very important in the low-density regime. We note that inversion layers give results similar to QW’s, but the specific numbers depend on the details of the doping profile.

In order to reinforce our conclusions, we present next a comparison between measured and calculated spin splittings in a GaAs-Al0.3Ga0.7As single heterostructure grown on a nominally undoped (311)A GaAs substrate with a weak $`p`$-type background doping.  A back gate was used to tune the density $`N`$ from $`1.8\times 10^{10}`$ to $`4.2\times 10^{10}`$ cm-2. To measure the spin subband densities $`N_\pm `$, the Shubnikov-de Haas (SdH) oscillations at low magnetic fields $`B`$ were examined  at a temperature $`T50`$ mK (see inset of Fig. 3). The frequencies $`f_{\mathrm{SdH}}`$ of these oscillations are a measure of the zero-$`B`$ spin splitting.  In Fig. 3(a) we present the measured and calculated spin subband densities exhibiting remarkably close agreement. Fig. 3(b) shows $`\beta _1^hE_z/\mu _h`$ determined by means of Eq. (6). On average, $`\beta _1^hE_z/\mu _h`$ increases as the density is reduced. Taking into account the orders-of-magnitude change that we have for $`\beta _1^hE_z/\mu _h`$ in QW’s and for $`\alpha _1E_z/\mu _c`$ in electron systems, the agreement between experiment and theory is quite satisfactory.  We wish to emphasize that it is indeed the anomalous enhancement of the Rashba coefficient in 2D HH systems in accumulation layer-like single heterostructures that allows us to experimentally resolve the spin splitting in this density regime; data on QW samples with comparable densities reveal no measurable spin splitting.

Support from NSF, DOE, and Humboldt Foundation are gratefully acknowledged.


text:
# The Globular Cluster Systems of Ellipticals and Spirals

## 1. Spatial Distribution

For the Milky Way (Minniti 1995) and other spiral galaxies (Forbes, Brodie & Larsen 2001) a consenus is emerging that the bulk of the metal-rich globular clusters (GCs) in spirals are associated with the bulge rather than the disk component, and the metal-poor GCs with the halo. A similar situation appears to be present in ellipticals with the red GCs associated with the elliptical bulge/spheroid and the blue GCs with the (hot gas) halo (Forbes, Brodie & Larsen 2001).

Globular cluster systems typically reveal a surface density distribution that is near constant close to the galaxy centre and which falls off rapidly like a power-law at larger radii. By fitting a King core profile to the GC density profile, a system core radius can be defined. Early-type galaxies reveal a trend for more massive galaxies to have larger core radii (Forbes et al. 1996). The Milky Way GC system, with a core of $``$ 1 kpc, seems to obey this trend.

## 2. Numbers and Specific Frequency

For small ellipticals, the number of GCs scales with the galaxy luminosity, while the most massive ellipticals appear to have a slightly steeper scaling relation (Djorgovski & Santiago 1992). One possible explanation for this comes from the work of McLaughlin (1999). This model requires a universal efficiency factor which is defined to be the ratio of mass in GCs to that of the initial gas mass available to form those GCs. For most ellipticals the initial gas supply has been converted into stars and is hence well represented by the current stellar mass (or luminosity). As a consequence these galaxies have a near constant SN $``$ 3. For high luminosity galaxies, McLaughlin argues that the current stellar mass grossly underestimates the initial gas supply and so they appear to be rich in GCs when they should be considered as poor in field stars. However by including the current gas content (as indicated by the presence of extensive X-ray gas), along with the current stellar mass, the GC systems of such galaxies can be understood as part of a continuous trend. McLaughlin’s relation can successfully explain the wide variation in GC number with luminosity.

Fig. 1 shows the McLaughlin relation with a universal efficiency factor of 0.2%. The data for early-type galaxies scatters uniformly about the relation. I have also included the location of three well-studied spiral galaxies (using the bulge luminosity), i.e. M104 (NGC 4594), M31 and the Milky Way (data from Forbes et al. 2000 and Larsen et al. 2001). The GC systems of the three spirals are consistent with the early-type galaxy number–luminosity trend.

## 3. Metallicities and Abundances

Historically, one difference between the GCs in spirals and ellipticals was the mean metallicity of the two peaks. The Milky Way and M31 have GC subpopulations with metallicities of \[Fe/H\] $``$ –1.5 and –0.5 (Barmby et al. 2000) while ellipticals were thought to have mean GC \[Fe/H\] $``$ –1.0 and 0.0 (Harris 1991). Recently, two developments have caused us to reasses the mean GC metallicity in ellipticals – both of which make the GCs in ellipticals more metal–poor. The first effect is the use of more accurate transformations from optical colours to \[Fe/H\]. For example, the new transformation of Kissler–Patig et al. (1998) converts a typical V–I = 1.05 to \[Fe/H\] = –1.07, where the old Galactic–based transformation would give \[Fe/H\] $``$ –0.5. The second effect is that the more accurate Galactic extinction values of Schlegel et al. (1998) tend to be larger on average by up to AV $``$ 0.1, than the traditionally-used Burstein & Heiles (1984) values. Thus extinction-corrected GC colours are now bluer than previously, and more metal–poor when transformed. If these two effects are taken into account, the two GC subpopulations in ellipticals have mean metallicities of \[Fe/H\] $``$ –1.5 and –0.5 which is similar to that of the M31 and Milky Way GC systems. Thus there appears to be very little difference between the mean metallicity of the two subpopulations in spirals and giant ellipticals.

The GC metallicity distributions for a wide range of Hubble types and luminosities is shown in Fig. 2 (including the ‘Local Group Elliptical’ – see Forbes et al. 2000). They are remarkably similar. All galaxies appear to possess a population of metal-poor (\[Fe/H\] $``$ –1.5) and presumably old GCs, while all bulge systems also have some metal-rich GCs (\[Fe/H\] $``$ –0.5). Hilker (this conference) confirms the bimodality in NGC 3311 with the same peaks.

There has been some debate over the years about the exact nature of the GC mean metallicity vs mass relation. To investigate this in more detail, we have created the largest sample to date of (37) early-type galaxies from the literature with bimodal GC colour distributions. This is therefore larger than the Forbes & Forte (2001), Brodie (this conference) or Kundu (this conference) samples. We also include M104, M31 and the Milky Way. We examine the correlation of colour (metallicity) with galaxy velocity dispersion (mass) in Fig. 3.

The Spearman rank correlation test finds that both the red and blue GCs are correlated with galaxy velocity dispersion with probabilities of 99.9% and 99.2% respectively. A least squares fit gives a positive slope for the red GCs at the 4$`\sigma `$ level and 2.5$`\sigma `$ for the blue GCs. The mean colors of the GC subpopulations in the three spirals are consistent with the best fit relation.  We also show the galaxy color – velocity dispersion relation which has the same slope as for the red GCs.

The existence of a GC color (metallicity) – galaxy mass relation, indicates a common chemical enrichment history for the red GCs and the host galaxy (Forbes & Forte 2000). Additionally, this relation has a similar slope to the galaxy bulge metallicity–mass relation. Thus the red GCs and the galaxy (bulge) stars appear to have formed in the same star formation event with fairly similar ages and metallicities for both early and late type galaxies. The red GCs and bulges of spirals follow these same trends and hence a similar evolutionary history.

The situation for the blue GCs is not quite so clear. The slope and hence statistical significance of a relation between colour and mass is less than for the red GCs. If the slope is indeed zero, i.e. there exists a constant mean colour of about (V–I)o = 0.95, then it indicates that the blue GCs formed before the potential well of the galaxy did. This suggests a ‘pregalactic’ origin (Forbes, Brodie & Grillmair 1997). We note that the model of GC formation at cosmological reionization proposed by Cen (2001) predicts no strong correlation between blue GC mean colour (metallicity) and galaxy mass. Alternatively, a relation with galaxy mass would indicate that they formed in situ in the current potential well.

In terms of abundances, it appears that from the few GC spectra available GCs in ellipticals have normal, ie Galactic, abundance ratios. This implies \[Mg/Fe\] $``$ +0.3. These supersolar abundances suggest either a very short star formation timescale or an IMF that is skewed to high mass stars. It will be interesting to extend this work to GCs in low luminosity galaxies for which the field stars have lower \[Mg/Fe\] ratios (Terlevich & Forbes 2001). From the limited data available it seems that GCs in spirals and ellipticals have the same abundance ratios.

## 4. Sizes

An advantage of HST over ground-based telescopes for the study of GCs is the ability of HST to measure GC sizes. The largest study to date is that of Larsen et al. (2001). In Fig. 4 we show the effective radius for the blue and red GC subpopulations in a wide range of galaxies, i.e. from spiral to elliptical to cD galaxy. The figure shows that the general trend of the blue GCs being $``$20% larger than the red ones exists over this large range in Hubble types. What ever process is causing this difference (eg destruction on different orbits or initial conditions) it is doing so irrespective of the host galaxy.

## 5. Ages and Kinematics

Determining the ages, and to a lesser extent kinematics, of GCs is perhaps the key to discriminating between the different formation models. Although controversial it appears that the metal-rich GCs in our Galaxy may be $``$ 2 Gyrs younger than the metal-poor ones. Observations to date of GCs in ellipticals suggest that both subpopulations are old and possibly coeval – a difference of a few Gyrs can not be ruled out.

Kinematics from large numbers of GCs are now available for one elliptical (M49) and two cD galaxies (M87, NGC 1399). Kinematics also hold important clues about formation mechanisms, although we note that both the merger and multi-phase collapse scenarios would probably expect the blue GCs to show higher V/$`\sigma `$ than the red ones. See Bridges (this conference) for further details.

## 6. Concluding Remarks

In terms of the spatial distribution, numbers, chemical properties, and sizes the GCs in spirals and ellipticals are remarkably similar. It is perhaps too early to say about the ages and kinematics. However we should take seriously the possibility that GC systems formed in a similar way in all galaxies. We suggest the following formation timeline. About 15 Gyrs ago, GCs with \[Fe/H\] $``$ –1.5 formed in all galaxies. They are associated with galaxy halos. A few Gyrs later, metal-rich \[Fe/H\] $``$ –0.5 GCs form along with galaxy bulges in a clumpy gaseous collapse. This is the main epoch of galaxy field star formation (see Forbes, Brodie & Grillmair 1997 for more details). In addition, recent galaxy mergers form elliptical galaxies with moderate SN values (ie $``$ 3).

## 7. Acknowledgments

I thank my SAGES colleagues M. Beasley, J. Brodie, C. Grillmair, J. Huchra, M. Kissler-Patig, S. Larsen. Part of this research was funded by the Ian Potter Foundation and NSF grant AST 9900732.

## References

Djorgovski, S., Santiago, B., 1992, ApJ, 391, L85

Burstein, D., Heiles, C., 1984, ApJS, 54, 33

Cen, R., 2001, astro-ph/0101197

Forbes, D. A., et al. 1996, AJ, 467, 126

Forbes, D. A., Brodie, J. P., Grillmair, C. J., 1997, AJ, 113, 1652

Forbes, D. A., et al. 2000, A&A, 358, 471

Forbes, D. A., Forte, J. C., 2001, MNRAS, 322, 257

Forbes, D. A., Brodie, J. P., Larsen, S. S., 2001, ApJ, submitted

Harris, W. E., 1991, ARA&A, 29, 543

Kissler-Patig, M., et al. 1998, AJ, 115, 105

Larsen, S. S., Forbes, D. A., Brodie, J. P., 2001, MNRAS, in press

Larsen, S. S., et al. 2001, AJ, in press

McLaughlin, D., 1999, AJ, 117, 2398

Minniti, D., 1995, AJ, 109, 1663

Schlegel, D. J., Finkbeiner, D. P., Davis, M., 1998, ApJ, 500, 525

Terlevich, A. I., Forbes, D. A., 2001, MNRAS, submitted


text:
# Complex poles and oscillatory screened potential in QED plasmas

## Abstract

We discuss the mechanism for the oscillatory behavior of the static interparticle potential in a degenerate electron plasma. This behavior, observed in metallic alloys, is commonly referred to as ’Friedel oscillations’, and its origin associated to the Kohn singularity. We show that, although this interpretation is adequate for large distances, the oscillations at short distances originate from a complex pole of the in-medium photon propagator in the complex q-plane, which exists aside the (purely imaginary) Debye pole. Such short-range oscillations can be physically discriminated if they remain at temperatures for which Friedel oscillations have already disappeared. This is suggested by finite temperature calculations in non-electromagnetic plasma models showing a similar pole structure.

It is well-known that wave propagation and particle interaction in vacuum are substantially modified by medium effects inside a plasma at a given density and temperature . In the case of a QED plasma, ordinary transverse photons (with respect to the wavelength vector) acquire an ’effective mass’, and are governed by a dispersion relation which substantially differs from the vacuum wave equation. In addition to these photons, there appear also longitudinal modes with different propagation properties . These modes arise from the poles of the photon propagator in the medium, and give raise to a number of phenomena in plasma and solid state physics. It is then important to analyze the structure of the photon propagator and its physical consequences. In this letter, we will concentrate on this study for the static photon propagator at zero temperature in the complex q-plane (where q is the momentum of the photon), and examine some of the consequences of our findings when considering the interparticle potential in the plasma.

We calculate the quasi-photon propagator in the RPA approximation to the dielectric function, and include a local field correction, as given by Ichimaru  (we have also performed calculations with other parameterizations of this correction). The expression of the static potential between two point-like impurities of charge $`e`$ is obtained from this propagator, and reads :

$$V(r)=\frac{e^2}{4\pi ^2r}\text{Im}_{\mathrm{}}^{\mathrm{}}𝑑q\frac{qe^{iqr}}{D(q)}$$
(1)
where $`D(q)=q^2ϵ(\omega =0,q)`$, $`ϵ(\omega ,q)`$ is the dielectric function and $`r`$ the interparticle distance.

The denominator $`D(q)`$ contains the function $`f(q)=\mathrm{log}|\frac{q+2p_f}{q2p_f}|`$, with $`p_f`$ the Fermi momentum of the electrons. This function can be written, for real values of $`q`$, as follows :

$$f(q)=\{\begin{array}{c}f_1(q)2\mathrm{arg}\mathrm{tanh}(\frac{q}{2p_f}),|\frac{q}{2p_f}|<1\\ f_2(q)2\mathrm{arg}\mathrm{tanh}(\frac{2p_f}{q}),|\frac{2p_f}{q}|<1\end{array}$$
(2)
We define the functions $`D_1(q)`$ and $`D_2(q)`$, which are obtained by replacing $`f(q)f_1(q)`$ and $`f(q)f_2(q)`$, respectively, in $`D(q)`$. Then, Eq. (1) becomes

$$V(r)=\frac{e^2}{4\pi ^2r}\mathrm{𝐼𝑚}\left[_{\mathrm{}}^{2p_f}𝑑q\frac{qe^{iqr}}{D_2(q)}+_{2p_f}^{\mathrm{}}𝑑q\frac{qe^{iqr}}{D_2(q)}+_{2p_f}^{2p_f}𝑑q\frac{qe^{iqr}}{D_1(q)}\right]$$
(3)
We next consider the analytical continuation of the function $`D_1(q)`$ to the complex q-plane, and apply the residue theorem in order to transform the last integration in Eq. (3) into the integration domain of the first two integrals in the same equation. After this operation, and with the help of the symmetry properties of the photon propagator due to the Onsager relations, Eq. (1) can be cast under the form :

$$V(r)=V_P(r)+V_C(r)$$
(4)
Here, $`V_P(r)`$ is the contribution of the residues of $`1/D_1(q)`$ in the region $`Im(q)>0`$. We have also defined

$$V_C(r)=\frac{e^2}{2\pi ^2r}\mathrm{𝐼𝑚}\left[_{2p_f}^{\mathrm{}}𝑑qqe^{iqr}\left(\frac{1}{D[f_2(q)]}\frac{1}{D[f_2(q)+i\pi ]}\right)\right]$$
(5)
The analytical structure of $`1/D_1(q)`$ on the upper half-plane is shown in Fig. 1, where we have also drawn the integration contour used to perform the integration. As we discuss later, in addition to the well-known ’Debye pole’, lying on the imaginary axis, *there exist complex poles*, which are responsible for strong oscillations of the interparticle potential at short distances.

In fact, while the oscillatory behavior of the screened potential created by a ionic impurity in metallic alloys is known since the works of Friedel  (see also ), its algebraic decay with distance has been associated to the existence of the Kohn singularity  in the photon self-energy, induced by the sharp edge of the degenerate electron distribution: when considering energy and momentum conservation in the collisions between the electrons and soft quasi-photons, we can easily see that only a fraction of the electrons can interact with quasi-photons and, therefore, screening of a Fermi gas at T=0 K is less effective than screening at finite temperature. At the threshold of the interaction, i. e. at the threshold of electron-hole creation, the momentum of the quasi-photon is equal to the diameter of the Fermi sphere. It must be noticed that the photon self-energy presents a singularity at this point. The presence of a singularity permits to obtain, with the help of Lighthill’s method, the asymptotic form ($`r\mathrm{}`$) of the potential, as an expansion in terms of the form $`\mathrm{cos}(2p_fr)`$ and $`\mathrm{sin}(2p_fr)`$, damped as negative powers of $`r`$, and enhanced by powers of $`\mathrm{log}(4p_fr)`$ . One has to remember, for further discussion, that Lighthill’s method consists essentially in replacing the Fourier (anti) transform of a non-analytical function by the (anti) transform of the (formal) expansion of that function around its singular points.

In Fig. 2 we have plotted the result for the first, dominant term ($`\mathrm{cos}(2p_f)/r^2`$), of Lighthill’s expansion (dotted line) for the potential around an impurity placed into an electron gas with ***This is a typical value corresponding to metals. $`r_s=3`$, where $`r_s`$ is the characteristic interparticle distance in units of the Bohr’s radius. In the same figure we have also plotted for comparison (solid line) the result of a numerical calculation of $`V(r)`$, as obtained by direct integration from Eq. (1). As we see, the Lighthill expansion converges to the true result for large values of $`r`$. An even better approximation to the exact result is given by $`V_C(r)`$, as defined above, and represented by the dashed line. The reason for this good agreement is that $`V_P(r)`$ exponentially goes to zero when $`p_fr1`$, and thus only $`V_C(r)`$ is important for long distances. Of course, the advantage of the Lighthill’s method for large distances is to provide a systematic way to construct corrections to the simple formula we used here, and to give an analytical expression for these terms, although they become rapidly cumbersome as the order of the approximation increases.

However, as we go to shorter distances, the situation becomes quite different. The reason is that, in this range, the dominant contribution to the form of the potential is given by the poles of the photon propagator in the complex plane, which induce exponentially-damped oscillations. This is clearly seen in Fig. 3, where we have represented the different contributions to the potential for the same value $`r_s=3`$, but for a shorter distance range. As before, the exact result (obtained from direct numerical integration) has also been plotted. In order to calculate $`V_P(r)`$, we have to locate numerically the zeros of $`D_1(q)`$ in the upper half-plane, while $`V_C(r)`$ is calculated by numerical integration. As we already mentioned, we find a zero lying on the imaginary-$`q`$ axis, which is commonly referred to as the ’Debye pole’. In addition to this, there exists a genuinely complex pole, with both a real and an imaginary part, which gives raise to exponentially-damped oscillations. As we can see, the contribution from the poles $`V_P(r)`$, represented by the dashed-dotted line, gives a much closer result to the exact line than the Lighthill term. We have also plotted, for the sake of completeness, the contribution from $`V_C(r)`$, which accounts for the difference between the poles term and the exact potential We have verified that the addition of $`V_C(r)`$ to $`V_P(r)`$ reproduces the exact potential, within machine-size precision numbers. . Thus, for this range of distances, the potential is dominated by the pole contribution. As a matter of fact, for the present value of $`r_s`$ the contribution of the Debye and complex poles to $`V_P(r)`$ are of the same order.

As the density is changed, the relative importance of these two contributions to $`V_P(r)`$ changes. This can be understood from Fig. 4, where we have shown the evolution of the real and imaginary parts of the complex pole, and the imaginary part of the Debye pole. For large densities ($`r_s2`$), $`V_P(r)`$ reduces to the well-known Debye screened potential, since the contribution arising from the complex pole is strongly damped. For larger values of $`r_s`$, however, the situation is reversed, since $`V_P(r)`$ is dominated by the complex pole and, therefore, this term alone constitutes a good approximation to the exact potential.

Let us summarize our results. We have studied the static interparticle potential in a QED electron plasma at zero temperature. This can be obtained from the knowledge of the dielectric function (or equivalently, of the photon polarization), which was taken from the RPA approximation, with local field corrections included. For the latter, we adopted the formulae given by Ichimaru, although we have performed also calculations with the formulae of Ref. , with very similar results. We found that the potential shows an oscillatory behavior for the whole distance range. This behavior has been traditionally referred to as Friedel oscillations, and its origin attributed to the Kohn singularity. We have shown that one can get a much better insight into this behavior through the study of the analytic structure of the photon propagator in the complex-$`q`$ plane. We find that, in addition to the well-known Debye pole on the imaginary axis, there exists a complex pole, which accounts for the oscillations of the potential at short ($`r1\AA `$) distances.

A similar pole structure, and the associated short-range oscillations, has been found in other (non electromagnetic) plasma models describing the one-pion exchange  and the exchange of other mesons in nuclear matter . In these cases, one can easily show that the separation of the screened potential into a short-range component and a Friedel-like component is not merely a mathematical artifice, but in fact can have physically observable consequences. Indeed, the phenomenon we have revisited (Friedel oscillations) and the new one obtained here come from two different effects. Friedel oscillations arise from the non-analytical behavior of the Fourier transform of the potential created by an impurity, while our results come from the zeros of the dispersion relation, in the complex plane, of the quasi-photon momentum. Using Maxwell equations, it can be easily obtained the charge distribution from the potential of the system. It comes out that the form of the charge and potential distribution are qualitatively identical. The presence of an impurity, then, leads to the apparition of a spatial static structure of charges in the medium, which is characteristic of a highly-interacting system. As long as the temperature increases, the kinetic energy of particles increases, and this structure breaks down. At the same time, the non-analytic behavior of the photon self-energy disappears at finite temperature, and therefore Friedel oscillations are exponentially damped with temperature.

This effect has been explicitly shown in  for the one-pion exchange. We also found that the complex pole is more stable, and remains present (as well as the associated short-range oscillations), for higher temperatures. This fact should permit to experimentally distinguish between both phenomena.

Even if the immediate extrapolation of these results to the present case might not be justified, the similarity of the analytic structure of the dressed boson propagators in both systems at $`T=0`$ suggest that a similar behavior can be expected for the electronic plasma. If such a thermal behavior is genuine, it could probably be experimentally tested in heated liquid metals, although a calculation of the characteristic temperatures at which Friedel and short-range oscillations disappear will be necessary.

Consequently, in order to get a better understanding of this phenomena and their possible consequences in plasma physics, it would be interesting to perform a study similar to the present one at non-zero temperature. Unfortunately, to do this we need a local-field correction formula valid at finite temperature. To our knowledge, such formula has not been given yet in the literature.

Acknowledgments

This work has been partially supported by the Spanish Grants DGES PB97-1432 and AEN99-0692.


text:
# Penetration of dynamic localized states in DC-driven Josephson junction ladders by discrete jumps

## I Introduction

The subject of large-amplitude anharmonic dynamics in lattices has received widespread attention over the past decade. In particular, intense theoretical focus has centered on so-called intrinsic localized modes (ILMs), also known as discrete breathers, with the result that many of their properties are now well understood. These excitations result from the interplay between nonlinearity and discreteness, and they can be highly localized in perfect lattices, with or without external driving. They can occur in a variety of different lattices: recent experiments have reported vibrational ILMs in a quasi-1D charge-density wave system, spin-wave ILMs in a quasi-1D antiferromagnetic system, and discrete breathers in Josephson junction (JJ) ladders.

The latter systems are noteworthy, in that arrays of coupled JJs have served for many years as reliable laboratory systems for studying diverse nonlinear phenomena. The nonlinear dynamics are particularly rich. A single “small” JJ subject to an applied constant DC bias current can be mapped onto the problem of a damped pendulum driven by a constant torque, with the dynamical degree of freedom being the Josephson phase difference. There are thus two qualitatively different states, namely a static (superconducting) state and a dynamic (whirling or resistive) state, with the latter producing a readily measured voltage $`V\dot{\phi }`$ across the whirling junction. When several junctions are assembled to form a regular array, such as the ladder shown in Fig. 1, they become inductively coupled. In the coupled system, junctions in the superconducting state can also exhibit steady-state librations, when JJs in the whirling state are present. In view of the mapping onto the pendulum problem, JJ ladders share features with lattices of nonlinearly coupled electric dipole rotors, driven by an external monochromatic AC electric field, but with the important simplification that they can be driven with purely DC bias currents.

Figure 1 sketches an anisotropic JJ ladder consisting of small JJs of two types, “horizontal” and “vertical,” which are, respectively, perpendicular and parallel to the applied bias current (arrows). The anisotropy arises from the different areas of the horizontal and vertical junctions and is characterized by the parameter $`\eta =A_h/A_v=I_c^h/I_c^v`$, where $`I_c^h`$ and $`I_c^v`$ are the the critical currents for each type of junction.

References  report experimental observations of various discrete breathers in ladders driven by a homogeneously applied DC bias current, represented by the dashed arrows in Fig. 1. For these states, the localized voltage patterns have a simple structure involving only two nonzero steady-state voltages (rotational frequencies). The breathers were found to be stable in the limit of small coupling ($`\eta 1`$) and for bias currents $`I_{\mathrm{ext}}I_c^v`$. For the case of large 2-D JJ arrays subject to a homogeneously applied DC bias, more complicated inhomogeneous states, with meandering voltage patterns, have also been reported.

Here, we study the dynamics of a JJ ladder with an external DC bias current applied at only one edge (solid arrows, Fig. 1). For increasing bias ($`I_{\mathrm{ext}}I_c^v`$) and over a wide range of anisotropies, we find by direct numerical simulations that the dynamic state expands into the ladder one cell at a time, by a sequence of abrupt jumps. This behavior is in marked contrast to the well-known cases of long JJs and JJ parallel arrays ($`\eta =\mathrm{}`$), where the entire system abruptly switches to the resistive state at a particular value of the DC bias. It is also different than the breather case, since all of the junctions within the boundary of this localized dynamic state whirl, and at each expansion the number of different frequencies (voltages) grows. The sequence of $`I`$-$`V`$ characteristics and threshold currents can be modeled analytically, yielding very good agreement with the numerical results.

## II Numerical Simulations

We consider a ladder with a large but finite number of cells $`N`$. The ladder’s state is specified by the time-dependent Josephson phases $`\{\phi _i\}`$, $`\{\psi _i\}`$, and $`\{\stackrel{~}{\psi }_i\}`$ for the vertical, upper horizontal and lower horizontal junctions, respectively, where $`i`$ denotes the cell. We have found in our simulations that the symmetry condition $`\stackrel{~}{\psi }_i=\psi _i`$ holds for the phenomena to be discussed here. The ladder dynamics are then determined by the coupled nonlinear equations of motion obtained in Refs.  and :

$`\widehat{L}(\phi _i)`$ $`=`$ $`\gamma _i+{\displaystyle \frac{1}{\beta _L}}[\phi _{i1}2\phi _i+\phi _{i+1}+2(\psi _i\psi _{i1})],`$ (2)
$`i=2,\mathrm{},N1,`$
$`\widehat{L}(\psi _i)`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\phi _i\phi _{i+1}2\psi _i),i=1,\mathrm{},N,`$ (3)
where the operator $`\widehat{L}(\phi )\ddot{\phi }+\alpha \dot{\phi }+\mathrm{sin}(\phi )`$. The equations for the vertical junctions at $`i=1`$ and $`N`$ are

$`\widehat{L}(\phi _1)`$ $`=`$ $`\gamma _1+{\displaystyle \frac{1}{\beta _L}}(\phi _2\phi _1+2\psi _1),`$ (4)
$`\widehat{L}(\phi _N)`$ $`=`$ $`\gamma _N+{\displaystyle \frac{1}{\beta _L}}(\phi _{N1}\phi _N2\psi _{N1}).`$ (5)
Equations (2) and (4) describe each junction within the resistively and capacitively shunted junction (RCSJ) model, and the unit of time is the inverse of the plasma frequency $`\omega _J\sqrt{2eI_c/C\mathrm{}}`$. Since each junction’s critical current and capacitance scale with the area, $`\omega _J`$ is independent of the anisotropy parameter $`\eta `$, as is the effective damping constant $`\alpha 1/(\omega _JRC)`$. The normalized bias current $`\gamma _i`$ is defined as $`I_{i,\mathrm{ext}}/I_c^v`$. The inductive coupling between junctions is determined by the parameter $`\beta _L2\pi LI_c^v/\mathrm{\Phi }_0`$, where $`L`$ is the self-inductance of a single cell and $`\mathrm{\Phi }_0=hc/2e`$ is the elementary flux quantum. Coupling beyond that described by $`\beta _L`$ is not included.

We performed direct numerical integration of the equations of motion for ladders with $`N=20`$ cells, using a fifth-order Gear predictor-corrector algorithm, for a range of anisotropies: $`\eta =`$ 0.5, 1.0, 2.0, 3.0, and 5.0. The arrays were underdamped, with $`\alpha =0.1`$, and we used a moderate value of the coupling parameter $`\beta _L=0.5`$. The external DC bias was applied at one edge, i.e. $`\gamma _1=\gamma `$ and all other $`\gamma _i=0`$. To simulate the $`I`$-$`V`$ curves, we started with all phases at zero and gradually increased the external bias $`\gamma `$ from zero to 50, in increments of 0.005. When junctions were present in the whirling state, the MD time scale was set by the time-average period of the fastest rotating phase. For a given value of gamma, we waited for at least 100 of these reference periods before averaging, in order to avoid transients, following which we computed the time-averages $`\dot{\phi }_i`$ and $`\dot{\psi }_i`$ over at least 100 additional reference periods. These averages are proportional to the average voltages across the junctions. The current was then incremented, with the initial phase configuration being that from the preceding MD time step. In all runs, the time step was 1/200 of the reference period.

Our simulated $`I`$-$`V_i`$ curves for an anisotropy of $`\eta =2`$ are shown in Fig. 2. The most striking finding is the occurrence of extremely sharp voltage jumps. At each of the corresponding threshold currents $`\gamma _n^{thr}`$, a new cell is added to the ladder’s dynamic state. With the applied current below the first threshold $`\gamma _1^{thr}`$, all junctions are in a static (zero voltage) state. When $`\gamma `$ exceeds $`\gamma _1^{thr}`$, the first vertical junction and its adjacent top and bottom horizontal junctions abruptly switch into the dynamic state, with all other junctions remaining in the zero voltage state—the rotating phases are confined to the first cell. As the bias is increased further, all three average voltages for this 1-cell dynamic state increase linearly until the next threshold current $`\gamma _2^{thr}`$ is reached, at which point the dynamic state suddenly expands into the second cell, accompanied by sharp changes of the voltages in the first cell. This process continues, yielding successive transitions from $`n`$-cell dynamic localized states to $`(n+1)`$-cell dynamic states. The distribution of threshold currents and voltage ratios depends on the ladder’s anisotropy. Over the range $`0<\gamma <50`$, the ladder reached a 3-cell state for $`\eta =0.5`$ and 1, a 4-cell state for $`\eta =2`$, and a 5-cell state for $`\eta =3`$ and 5. In the following, these states will be termed $`n`$-cell edge states.

An $`n`$-cell edge state is in striking contrast to an n-cell breather. The breather occurs away from the ladder’s edge and is homogeneously driven by same DC bias current ($`I_{\mathrm{ext}}I_c`$) applied at every cell, whereas the edge states are driven by a DC bias ($`I_{\mathrm{ext}}I_c`$) applied at just one edge. The edge states have a richer internal structure—all of the junctions within an edge state are in a nonzero voltage state (see Fig. 3), whereas in a breather state, all of the horizontal junctions are in the zero voltage state except for those on the breather’s boundary. Moreover, all of a breather’s vertical junction phases rotate at the same average frequency, whereas the $`n`$-cell edge state exhibits a peculiar distribution of average frequencies. This frequency (voltage) distribution depends on both $`n`$ and the ladder’s anisotropy. For example, our simulations for the $`\eta =2`$, 3-cell edge state of Fig. 3(d) yield the ratios given in second column of Table I. Comparison with the third column shows that they are in very good agreement with analytic predictions derived below.

The superconducting state forming ahead of the $`n`$-cell edge state is also unusual. Fig. 3 gives snapshot images of the Josephson phase distribution for several values of the applied DC current bias, for the anisotropy $`\eta =2`$. In panel (a), the current is just below the first threshold, and one sees a single Josephson vortex in the superconducting part of the ladder. The remaining panels (b)–(e) show the phases just after a new cell is added to the dynamic state. At the threshold currents $`\gamma _n^{thr}`$, the superconducting state becomes unstable, and the vortex jumps into the next cell as the edge state expands. Our simulations show that in general the superconducting state is sensitive to the anisotropy. Thus for rather small values of $`\eta 1`$, there are no vortices trapped in the superconducting part of the ladder over the range of currents studied. For these cases, the Josephson phases of the vertical junctions in the superconducting portion of the ladder simply decrease with distance from the boundary of the resistive portion, corresponding to the Meissner state of the superconductor. With increasing anisotropy, single Josephson vortices appear in the superconducting portion, as in Fig. 3. For large anisotropies ($`\eta 5`$) more complex vortex trains are observed, and we also find that the penetration of the edge state changes the nature of the superconducting vortex state, rather than simply pushing it ahead as for $`\eta =2`$. A detailed discussion of the superconducting state will be given elsewhere (Ref. ).

## III Theoretical Analysis

The unusual voltage distributions in the $`n`$-cell edge states can be explained analytically by making use of Kirchhoff’s laws, applied to the time-average currents (normalized to $`I_c^v`$) and corresponding dimensionless voltages in each cell. The key assumption is the coexistence of the resistive and superconducting states in different portions of the ladder. For a cell $`i`$ within an $`n`$-cell edge state, current conservation gives $`I_i^v+I_i^h=I_{i1}^h`$, while the voltage condition is $`I_i^v\frac{2}{\eta }I_i^hI_{i+1}^v=0`$. Combining these yields an equation for just the horizontal currents:

$$I_{i+1}^h+I_{i1}^h2\left(1+\frac{1}{\eta }\right)I_i^h=0.$$
(6)
This equation and the two from which it was derived apply to all cells $`i`$ in $`1in`$, provided we define $`I_0^h\gamma ,I_{n+1}^v0`$, and $`I_{n+1}^hI_n^h`$, in order to take proper account of the $`n`$-cell edge state’s boundaries.

Equation (6) is readily solved by substituting $`I_i^h=\lambda ^i`$, which yields two roots, namely $`\lambda 1+\frac{1}{\eta }+\sqrt{(1+\frac{1}{\eta })^21}`$ and $`1/\lambda `$. Hence the general solution of Eq. (6) is $`I_i^h=c_1\lambda ^i+c_2\lambda ^i`$, where the constants $`c_1`$ and $`c_2`$ are obtained from the above definitions of $`I_0^h`$ and $`I_{n+1}^h`$. With the horizontal currents thus determined, the vertical currents may be computed from $`I_i^v=I_{i1}^hI_i^h`$. The currents are then converted into the average voltages via $`V_i^v=I_i^v/\alpha `$ and $`V_i^h=I_i^h/(\alpha \eta )`$. The resulting voltage distribution within an $`n`$-cell edge state is ($`1in`$):

$$V_i^v=\frac{\gamma (1\lambda )(\lambda ^{i1}\lambda ^{2n+1i})}{\alpha (\lambda ^{2n+1}+1)},$$
(7)
$$V_i^h=\frac{\gamma (\lambda ^i+\lambda ^{2n+1i})}{\alpha \eta (\lambda ^{2n+1}+1)}.$$
(8)
Equations (7) and (8) give the predicted voltage ratios in Table I, which for $`\eta =2`$ are seen to be in very good agreement with our MD simulations. Indeed, we find that for all of the values of $`\eta `$ studied, the predicted $`I`$-$`V`$ curves are in very good agreement with the MD curves, such as those of Fig. 2. Only the values of the current thresholds for the jumps are left undetermined by these equations.

We can also predict the distribution $`\{\gamma _n^{thr}\}`$ of threshold currents for each $`\eta `$ by assuming that the superconducting state associated with the $`(n1)`$-cell edge state becomes unstable and converts to the $`n`$-cell edge state when the current $`I_n^h`$ reaches a depinning current $`I_{dp}`$, which we take to be independent of $`n`$. This yields an expression for the threshold currents

$$\gamma _n^{thr}=I_{dp}\frac{\mathrm{cosh}[(n\frac{1}{2})\mathrm{ln}\lambda ]}{\mathrm{cosh}(\frac{1}{2}\mathrm{ln}\lambda )}.$$
(9)
The ratios $`\gamma _n^{thr}/I_{dp}`$ predicted by Eq. (9) for $`n=`$ 2, 3, 4 and 5 are plotted versus $`\eta `$ in Fig. 4. To compare with the MD results, we fit $`I_{dp}`$ to the first observed MD threshold current for each $`\eta `$, namely $`\gamma _1^{thr}=`$ 1.295, 1.510, 2.040, 2.525, and 4.010, for $`\eta =`$ 0.5, 1.0, 2.0, 3.0, and 5.0, respectively. The resulting MD ratios are shown by the circles in Fig. 4 and agree well with the predictions.

## IV Hysteresis

Despite their rich structure of frequency ratios, the $`n`$-cell edge states are found to be highly stable in our simulations, for the case of increasing current. However, we also find notable hysteresis effects when the simulations are started with an $`n`$-cell edge state for large $`n`$ and the applied DC current is gradually decreased to zero. Figure 5 is representative of the hysteretic behavior encountered. The threshold currents and stability properties for the sequence of down-conversions $`\{nn1\}`$ are very different than for the increasing-current case. In particular, we observe resonant steps, switching processes, and nonlinear regions in the $`I`$-$`V`$ curves. We believe that all of these features arise from the resonant interaction between the n-cell edge states and other excitations, both localized and delocalized, as will be discussed elsewhere.

## V Summary

In summary, our numerical simulations have revealed unusual localized dynamic states in anisotropic JJ ladders subject to a DC bias current at one edge. Increasing the bias causes these states to expand by adding single cells in a sequence of sudden jumps, giving rise to a diverse set of voltage distributions and sharp changes in the $`I`$-$`V`$ curves. This behavior occurs for a wide range of parameters and should be observable through the $`I`$-$`V`$ characteristics or by direct visualization using low temperature scanning laser microscopy techniques.

###### Acknowledgements.

We thank A. V. Ustinov and S. Flach for useful discussions. J. B. Page gratefully acknowledges the Max Planck Institute for the Physics of Complex Systems, Dresden, for their support and hospitality.