text:
# Two-neutrino double beta decay within Fully-Renormalized QRPA: Effect of the restoration of the Ikeda sum rule

## I Introduction

Observation of the neutrinoless double beta decay ($`0\nu \beta \beta `$-decay), violating the total lepton number by two units, would give unambiguous evidence for new physics beyond the Standard Model . For instance, at least one of the neutrinos would have to be a Majorana particle with non-zero mass . The current experimental upper limits on the $`0\nu \beta \beta `$-decay half-life impose stringent constraints, e.g., on the parameters of Grand Unification and super-symmetric extensions of the Standard Model.

Rates of the $`2\nu \beta \beta `$-decay, which is a second order process allowed within the Standard Model, can be calculated within the same nuclear structure models. Thus, the results of the nuclear structure calculations can be directly compared with the corresponding experimental data available for a number of nuclei . Such a comparison provides a very useful test of the models.

The Quasiparticle Random Phase Approximation (QRPA)  has been successfully exploited in nuclear physics to describe properties of the excited states of open-shell nuclei and to calculate intensities of various nuclear reactions, including the double beta decay (see reviews ).

It was shown that the experimental data on the $`2\nu \beta \beta `$-decay rates can be reproduced in QRPA calculations with a sufficiently large strength of the particle-particle interaction . But the proximity of the value to the point of the QRPA collapse questions the reliability of the results. It is known that QRPA collapse occurs due to the use of the quasi-boson approximation (QBA) which violates the Pauli exclusion principle (PEP) and generates too many ground state correlations.

Renormalized QRPA (RQRPA) was formulated in Refs.  to restore PEP in an approximate way. The main goal of the method is to use a self-consistent iteration of the QRPA equation with taking into account quasiparticle occupation numbers in the QRPA ground state. That leads to a modification of the commutation relations for bifermionic operators as compared to the ordinary quasiboson approximation (QBA). At the same time so-called scattering terms (describing transitions of the quasiparticles) are neglected in the Hamiltonian and in the phonon operators. The RQRPA does not collapse for physical values of the particle-particle interaction strength and has been extensively used to calculate the intensities of the double beta decay . It has been also shown that the RQRPA provides better agreement with the exact solution of the many-body problem within schematic models, even beyond the critical point of the standard QRPA (see, e.g.  and references therein).

The self-consistent RQRPA (SCQRPA) is a more complex version of RQRPA to describe the strongly correlated Fermi systems. Within this method one goes a step further beyond the RQRPA. In the SCQRPA at the same time the quasiparticle mean field is changed by minimizing the energy and fixing the number of particles in the correlated ground state of RQRPA instead of the uncorrelated one of BCS as is done in the other versions of the RQRPA. In this way SCQRPA partially overcomes the inconsistency between RQRPA and the BCS approach and is closer to a fully variational theory.

Nevertheless, the main drawback of the modern versions of RQRPA and SCQRPA is the violation of the model-independent Ikeda sum rule (ISR) . A modification of the phonon operator by including scattering terms is needed in order to restore the ISR within RQRPA. The fully-Renormalized QRPA (FR-QRPA) was formulated in Ref.  for even-even nuclei in such a way that it complies with restrictions imposed by the commutativity of the phonon creation operator with the total particle number operator. It was shown analytically that the Ikeda sum rule is fulfilled within the FR-QRPA . Also FR-QRPA is free from the spurious low-energy solutions which would be generated by the scattering terms considered as additional degrees of freedom as suggested in .

The aim of the paper is twofold. First, we would like to describe the FR-QRPA equations in more details (as compared with the original paper ) for a simple case of a Hamiltonian with the separable residual interaction in both particle-hole and particle-particle channels. Second, the first numerical application of FR-QRPA is given to calculate $`2\nu \beta \beta `$-decay intensities and relevant quantities. So far, the full convergence of the FR-QRPA solution has been obtained only for a rather small model space. Nevertheless a comparison of the results obtained within FR-QRPA and SCQRPA can be provided.

## II Basic relationships of the Fully-Renormalized QRPA

Within RPA an excited nuclear state, with angular momentum $`J`$ and projection $`M`$, is created by applying the phonon operator $`Q_{JM}^{}`$ to the vacuum state $`|0_{RPA}^+`$ of the initial, even-even, nucleus:

$$|JM=Q_{JM}^{}|0_{RPA}^+\text{with}Q_{JM}|0_{RPA}^+=0.$$
(1)
As was shown in Ref. , the most appropriate way is to write down the phonon structure in terms of the particle creation and annihilation operators. That allows to fulfill the important principle of the commutativity of $`Q_{JM}^{}`$ with the total particle number operator $`\widehat{A}=\widehat{N}+\widehat{Z}`$. The phonon operator has the following structure:

$$Q_{JM}^{}=\underset{pn}{}\left[x_{(pn,J)}C^{}(pn,JM)y_{(pn,J)}\stackrel{~}{C}(pn,JM)\right],$$
(2)
with $`C^{}(pn,JM)=\left[c_p^{}\stackrel{~}{c}_n\right]_{JM}`$ and $`\stackrel{~}{C}(pn,JM)=()^{JM}C(pn,JM)`$, where $`c_{\tau m_\tau }^+`$ ($`c_{\tau m_\tau }`$) denotes the particle creation (annihilation) operator for protons and neutrons ($`\tau =p,n`$).

Going into the quasiparticle representation, the quasiparticle creation and annihilation operators $`a_{\tau m_\tau }^+`$ and $`a_{\tau m_\tau },(\tau =p,n)`$ can be defined by the Bogolyubov transformation

$$\left(\begin{array}{c}a_{\tau m_\tau }^+\\ \stackrel{~}{a}_{\tau m_\tau }\end{array}\right)=\left(\begin{array}{cc}u_\tau & v_\tau \\ v_\tau & u_\tau \end{array}\right)\left(\begin{array}{c}c_{\tau m_\tau }^+\\ \stackrel{~}{c}_{\tau m_\tau }\end{array}\right),$$
(3)
that leads to the following expression for the phonon operator $`Q_{JM}^{}`$:

$`Q_{JM}^{}={\displaystyle \underset{pn}{}}\left[X_{(pn,J)}\overline{A}^{}(pn,JM)Y_{(pn,J)}\stackrel{~}{\overline{A}}(pn,JM)\right],`$ (4)
(5)
$`\overline{A}^{}=A^{}+\left(u_nv_nB^{}u_pv_p\stackrel{~}{B}\right)/(v_n^2v_p^2)`$ (6)
$`A^{}(pn,JM)=\left[a_p^{}a_n^{}\right]_{JM};B^{}(pn,JM)=\left[a_p^{}\stackrel{~}{a}_n\right]_{JM}`$ (7)
where $`X=u_pv_nxv_pu_ny,Y=u_pv_nyv_pu_nx`$. The bifermionic operators $`\overline{A}^{},\overline{A}`$ now being the basic building blocks of the FR-QRPA automatically contain the quasiparticle scattering terms which, however, are not associated with any additional degree of freedom. That means that there are no spurious low-lying solutions in the present theoretical scheme which would be generated by the scattering terms considered as independent constituents of the phonon operator (as proposed in ).

From this point we can follow the usual way to formulate the RQRPA , substituting $`A`$ by $`\overline{A}`$ everywhere. The forward- and backward-going free variational amplitudes X and Y satisfy the equation:

$$\left(\begin{array}{cc}𝒜& \\ & 𝒜\end{array}\right)\left(\begin{array}{c}X^m\\ Y^m\end{array}\right)=_m\left(\begin{array}{cc}𝒰& 0\\ 0& 𝒰\end{array}\right)\left(\begin{array}{c}X^m\\ Y^m\end{array}\right),$$
(8)
where $`m`$ marks different roots of the QRPA equations for a given $`J^\pi `$,

$`𝒜`$ $`=`$ $`0_{RPA}^+|[\overline{A},[H,\overline{A}^{}]]|0_{RPA}^+,`$ (9)
$``$ $`=`$ $`0_{RPA}^+|[\overline{A},[H,\overline{A}]]|0_{RPA}^+,`$ (10)
and the renormalization matrix $`𝒰_{pn}`$ is

$`𝒰_{pn}`$ $`=`$ $`0_{RPA}^+|[\overline{A}(pn,JM),\overline{A}^{}(p^{}n^{},JM)]|0_{RPA}^+=\delta _{pp^{}}\delta _{nn^{}}𝒟_{pn}.`$ (11)
We use a rather simple, but realistic, Hamiltonian $`H`$ consisting of the quasiparticle mean field $`H_0`$ and the residual separable particle–hole (ph) and particle–particle (pp) interactions:

$`H=H_0+H_{int}^{ph}+H_{int}^{pp},`$ (12)
$`H_0={\displaystyle \underset{\tau =p,n}{}}E_\tau a_\tau ^{}a_\tau ,`$ (13)
$`H_{int}^{ph}=\chi {\displaystyle \underset{M}{}}(1)^M(\beta _{1M}^{}\beta _{1M}^++\beta _{1M}^+\beta _{1M}^{}),`$ (14)
$`H_{int}^{pp}=\kappa {\displaystyle \underset{M}{}}(1)^M(P_{1M}^{}P_{1M}^++P_{1M}^+P_{1M}^{}),`$ (15)
with $`\beta _{1M}^{}=\widehat{J}^1\underset{pn}{}p\sigma n\left[c_p^{}\stackrel{~}{c}_n\right]_{1M}`$, $`P_{1M}^{}=\widehat{J}^1\underset{pn}{}p\sigma n\left[c_p^{}c_n^{}\right]_{1M}`$ and $`J=1`$.

Taking into account the exact (fermionic) expressions for the commutators in (10),(11), one gets the following expressions for the FR-QRPA matrices $`𝒜`$ and $``$:

$`𝒜`$ $`=`$ $`\left[(E_p+E_n)𝒟_{pn}2(E_pE_n)(u_p^2v_p^2+u_n^2v_n^2)_{pn}\right]\delta _{pp^{}}\delta _{nn^{}}`$ (18)
$`+2\chi (u_pv_nu_p^{}v_n^{}+v_pu_nv_p^{}u_n^{})𝒟_{pn}𝒟_{p^{}n^{}}`$
$`2\kappa (u_pu_nu_p^{}u_n^{}\overline{\overline{𝒟}}_{pn}\overline{\overline{𝒟}}_{p^{}n^{}}+v_pv_nv_p^{}v_n^{}\overline{𝒟}_{pn}\overline{𝒟}_{p^{}n^{}}),`$
$``$ $`=`$ $`2(E_pE_n)u_pv_pu_nv_n_{pn}\delta _{pp^{}}\delta _{nn^{}}`$ (22)
$`+2\chi (u_pv_nv_p^{}u_n^{}+v_pu_nv_p^{}u_n^{})𝒟_{pn}𝒟_{p^{}n^{}}`$
$`+2\kappa (u_pu_nv_p^{}v_n^{}\overline{\overline{𝒟}}_{pn}\overline{𝒟}_{p^{}n^{}}+v_pv_nu_p^{}u_n^{}\overline{𝒟}_{pn}\overline{\overline{𝒟}}_{p^{}n^{}}).`$
The renormalization matrices $`𝒟,\overline{𝒟},\overline{\overline{𝒟}}`$ entering (11),(18),(22) can be represented in terms of the relative quasiparticle occupation numbers $`𝒩_p`$ for the level $`p`$ in the RQRPA vacuum:

$`𝒟_{pn}`$ $`=`$ $`1𝒩_n𝒩_p+\left(1v_p^2v_n^2\right)_{pn},`$ (23)
$`\overline{𝒟}_{pn}`$ $`=`$ $`1𝒩_n𝒩_p\left(u_p^2+u_n^2\right)_{pn},`$ (24)
$`\overline{\overline{𝒟}}_{pn}`$ $`=`$ $`1𝒩_n𝒩_p+\left(v_p^2+v_n^2\right)_{pn},`$ (25)
with $`_{pn}=\frac{𝒩_p𝒩_n}{v_n^2v_p^2}`$. In turn, the quasiparticle occupation numbers

$$𝒩_p=\widehat{j}_\tau ^20_{RPA}^+|\underset{m_\tau }{}a_{\tau m_\tau }^{}a_{\tau m_\tau }|0_{RPA}^+;\tau =p,n.$$
(26)
can be expressed in terms of the backgoing amplitudes $`Y`$ of the RQRPA solution (8) . In the calculation we shall use the aproximate expression for $`𝒩_p`$ and $`𝒩_n`$:

$`𝒩_p`$ $``$ $`\widehat{j}_p^2{\displaystyle \underset{n}{}}\left({\displaystyle \underset{J,m}{}}(2J+1)|Y_{pn,J}^m|^2\right)𝒟_{pn}`$ (27)
$`𝒩_n`$ $``$ $`\widehat{j}_n^2{\displaystyle \underset{p}{}}\left({\displaystyle \underset{J,m}{}}(2J+1)|Y_{pn,J}^m|^2\right)𝒟_{pn}`$ (28)
where $`\widehat{j}\sqrt{2j+1}`$. In the present paper we consider only $`J^\pi =1^+`$ contribution to the sums in (28). Along with the modified SCQRPA and FR-QRPA equations for the chemical potential:

$`0_{RPA}^+|\widehat{N}|0_{RPA}^+`$ $`=`$ $`{\displaystyle \underset{n}{}}\widehat{j}_n^2\left(v_n^2+(u_n^2v_n^2)𝒩_n\right)=N,`$ (29)
$`0_{RPA}^+|\widehat{Z}|0_{RPA}^+`$ $`=`$ $`{\displaystyle \underset{p}{}}\widehat{j}_p^2\left(v_p^2+(u_p^2v_p^2)𝒩_p\right)=Z,`$ (30)
a rather complicated set of equations (8)-(30) has to be solved.

It is noteworthy that the renormalization matrices (23) become the same, $`\overline{\overline{𝒟}}_{pn}=\overline{𝒟}_{pn}=𝒟_{pn}`$, in the limit $`_{pn}=0`$ and coinciding with the renormalization matrix of the usual RQRPA (see,e.g., ). Thus, one can argue that the standard versions of RQRPA neglect effectively the differences between the quasiparticle occupation numbers whereas SCQRPA and FR-QRPA take the differences into account.

From now on we follow the usual way of solving RQRPA equations . It is useful to introduce the notation:

$$\overline{X}=𝒰^{1/2}X,\overline{Y}=𝒰^{1/2}Y,$$
(31)
$$\overline{𝒜}=𝒰^{1/2}𝒜𝒰^{1/2},\overline{}=𝒰^{1/2}𝒰^{1/2}.$$
(32)
Then the amplitudes $`\overline{X}`$ and $`\overline{Y}`$ satisfy the equation of usual QRPA:

$$\left(\begin{array}{cc}\overline{𝒜}& \overline{}\\ \overline{}& \overline{𝒜}\end{array}\right)\left(\begin{array}{c}\overline{X}^m\\ \overline{Y}^m\end{array}\right)=_m\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\overline{X}^m\\ \overline{Y}^m\end{array}\right).$$
(33)
Solving the FR-QRPA equations, one gets the fully renormalized amplitudes $`\overline{X}`$, $`\overline{Y}`$ with the usual normalization and closure relations:

$`{\displaystyle \underset{pn}{}}\overline{X}_{(pn,J)}^m\overline{X}_{(pn,J)}^k\overline{Y}_{(pn,J)}^m\overline{Y}_{(pn,J)}^k`$ $`=`$ $`\delta _{km},`$ (34)
$`{\displaystyle \underset{m}{}}\overline{X}_{(pn,J)}^m\overline{X}_{(p_1^{}n_1,J)}^m\overline{Y}_{(pn,J)}^m\overline{Y}_{(p_1^{}n_1,J)}^m`$ $`=`$ $`\delta _{pp_1^{}}\delta _{nn_1},`$ (35)
$`{\displaystyle \underset{m}{}}\overline{X}_{(pn,J)}^m\overline{Y}_{(p_1^{}n_1,J)}^m\overline{Y}_{(pn,J)}^m\overline{X}_{(p_1^{}n_1,J)}^m`$ $`=`$ $`0.`$ (36)
It was shown analytically that the Ikeda sum rule is fulfilled within the FR-QRPA , in contrast to the earlier versions of the RQRPA . The Ikeda sum rule states that the difference between the total Gamow-Teller strengths $`S^{()}`$ and $`S^{(+)}`$ in the $`\beta ^{}`$ and $`\beta ^+`$ channels, respectively, is $`3(NZ)`$ :

$`ISR=S^{()}S^{(+)}=3(NZ),`$ (37)
$`S^{()}={\displaystyle \underset{Mm}{}}\left|1^+M,m|\beta _{1M}^{}|0_{RPA}^+\right|^2,S^{(+)}={\displaystyle \underset{Mm^{}}{}}\left|1^+M,m^{}|\beta _{JM}^+|0_{RPA}^+\right|^2.`$ (38)
With the use of the closure conditions (36), the expressions for $`𝒟_{pn}`$ (23) and the chemical potentials (30), one can show  that

$$ISR=\underset{pn}{}\left|pq_Jn\right|^2(v_n^2v_p^2)𝒟_{pn}=3(NZ).$$
(40)
The inverse half-life of the $`2\nu \beta \beta `$-decay can be expressed as a product of an accurately known phase-space factor $`G^{2\nu }`$ and the second order Gamow-Teller transition matrix element $`M_{GT}^{2\nu }`$:

$$[T_{1/2}^{2\nu }(0_{g.s.}^+0_{g.s.}^+)]^1=G^{2\nu }(g_A)^4|M_{GT}^{2\nu }|^2.$$
(41)
The contribution from the two successive Fermi transitions is safely neglected as they arise from isospin mixing effect . The double Gamow-Teller matrix element $`M_{GT}^{2\nu }`$ for ground state to ground state $`2\nu \beta \beta `$-decay transition acquires the form

$$M_{GT}^{2\nu }=\underset{m_im_f}{}\frac{0_f^+\beta ^{}1_{m_f}^+1_{m_f}^+|1_{m_i}^+1_{m_i}^+\beta ^{}0_i^+}{(\omega ^{m_f}+\omega ^{m_i})/2}.$$
(42)
The sum extends over all $`1^+`$ states of the intermediate nucleus. The index $`i(f)`$ indicates that the quasiparticles and the excited states of the nucleus are defined with respect to the initial (final) nuclear ground state $`|0_i^+`$ ($`|0_f^+`$). The overlap is necessary since these intermediate states are not orthogonal to each other. The two sets of intermediate nuclear states generated from the initial and final ground states are not identical within the considered approximation scheme. Therefore the overlap factor of these states is introduced in the theory as follows:

$$1_{m_f}^+|1_{m_i}^+=\underset{pn}{}[X_{pn}(1^+m_i)X_{pn}(1^+m_f)Y_{pn}(1^+m_i)Y_{pn}(1^+m_f)].$$
(43)
## III Calculation results

In this section we present the $`2\nu \beta \beta `$-decay results obtained within the FR-QRPA for 76Ge, 82Se, 100Mo, 116Cd, 128Te and 130Xe, in comparison with the QRPA and SCQRPA ones. Rather small model bases listed in the Table I are used in order to get full convergence in the FR-QRPA method. The levels are in a vicinity of the Fermi levels and spin-orbit partners are always taken into account. FR-QRPA method is rather sensitive to the differences between occupation probabilities for protons and neutrons entering the denominator in the expression of the bifermionic operators $`\overline{A}^{},\overline{A}`$ (4) and in the expresion for $`_{pn}`$ factor of renormalization matrices (23).

For levels far from the Fermi one, the values for occupation probabilities for protons and neutrons become almost equal. Because of the denominator which appears in the expression of bifermionic operators (4), that causes numerical problems, in particular the method doesn’t converge for large enough values of particle-particle strength. Therefore, the bases are fixed in order to get convergence for a larger interval of particle-particle strength, in particular up to the point of the collapse of the $`2\nu \beta \beta `$-decay matrix elements.

The single particle energies are obtained by using a Coulomb-corrected Woods-Saxon potential with Bertsch parametrization. The proton and neutron pairing gaps are determined phenomenologically to reproduce the odd-even mass differences through a symmetric five-term formula . Then the equations for the chemical potentials (30) are solved for proton and neutron subsystems. The pairing gaps entering the BCS equations are given in the Table I.

The calculation of the QRPA energies and wave functions requires the knowledge of the particle-hole $`\chi `$ and particle-particle $`\kappa `$ strengths of the residual interaction. The value of particle-hole strength $`\chi `$ parameter for each nucleus is fixed in order to reproduce the experimental position of the Gamow-Teller giant resonance in odd-odd intermediate nucleus as obtained from the (p,n) reactions , , . Those values are also given in the Table I. The particle-particle strength $`\kappa `$ is considered as a free parameter.

The numerical results are shown for two groups of nuclei, the nuclei with $`A100`$ and $`A>100`$ respectively. The calculations are done within QRPA, SCQRPA and FR-QRPA in order to show the better stability of the latter method.

The main drawback of the QRPA is the overestimation of the ground state correlations leading to the collapse of the QRPA ground state, near a certain critical interaction strength. Around this point the backward-going RPA amplitudes $`Y_{pn}`$ of the first $`1^+`$ states become overrated, and too many correlations in the ground state are generated with increasing strength of the particle-particle interaction. This phenomenon, as a result of the quasiboson approximation used, leads to QRPA collapse and implies an ambiguous determination of the $`\beta `$ and $`2\nu \beta \beta `$-decay matrix elements.

In Fig.1 and Fig.2 the dependence of the energy of the first excited Gamow-Teller state in daughter nuclei is plotted versus the $`\kappa `$ parameter. Hereafter, the dashed line corresponds to the QRPA case, the dotted line represents the SCQRPA case and the solid line describes FR-QRPA case. For all studied nuclei the collapse of the first excited state is shifted to higher values of $`\kappa `$ for each method and the stability increases in the FR-QRPA case. In Fig.3 and Fig.4 the $`2\nu \beta \beta `$-decay matrix elements as a function of the particle-particle strength $`\kappa `$ are shown. The calculations are done for all nuclei within the three metods. The horizontal dashed line indicates the experimental values taken from .

For all nuclei there is a similar behaviour in the sense that QRPA and SCQRPA collapse a bit earlier than the FR-QRPA does. Although the chosen bases are rather small, the new effects we intend to emphasise as the differences between QRPA extensions are evident. The FR-QRPA method offers considerably less sensitive dependence of $`M_{GT}^{2\nu }`$ on $`\kappa `$ and shifts the collapse to larger values of particle- particle strength.

According to the definition (37) $`S^{}`$ ($`S^+`$) is the total summed Gamow-Teller $`\beta ^{}`$($`\beta ^+`$) transition strength from the ground state of an even-even nucleus. In Fig.5 and Fig.6 we plot the relative $`\beta ^{}`$ strength, $`S^{}/3(NZ)`$ for the mother nucleus (left side) and relative $`\beta ^+`$ strength, $`S^+/3(NZ)`$ for the daugther (right side), for $`A100`$ and $`A>100`$ respectively as a function of particle-particle interaction parameter in order to show the magnitude and the nature of violation of ISR.

Now we would like to discuss the conservation of the Ikeda sum rule $`ISR=S_{}S_+=3(NZ)`$ in the FR-QRPA framework and to compare with the previous calculations for QRPA ans SCQRPA. We didn’t include the calculations for RQRPA because SCQRPA goes beyond and brings more improvements than RQRPA, especially for Ikeda sum rule.

Finally we combine the data of previous plots and show the ratio of SCQRPA and FR-QRPA sum $`ISR/3(NZ)`$ as a function of $`\kappa `$. In the QRPA the Ikeda sum rule is exactely conserved as long as all spin-orbit partners of the single-particle orbitals are included. In the other extended versions of QRPA the sum rule is violated with a degree of deviation lying between $`17\%`$ (RQRPA) and $`3\%`$ (SCQRPA) . In our study, following the analytical calculation of , we have shown numerically that Ikeda sum rule is exactly fullfiled within FR-QRPA formalism.

## IV Conclusions

In summary, the first calculation of the $`2\nu \beta \beta `$-decay matrix elements within the recently proposed Fully-Renormalized QRPA, which fulfills Ikeda sum rule exactly, are presented. The considered nuclear model includes the separable Gamow-Teller residual interaction. The subject of interest is the effect of the restoration of the Ikeda sum rule on the $`2\nu \beta \beta `$-decay observable for $`A=76,82,100,116,128,130`$ systems. Within the present work we arrived to the folowing important conclusions:

i) The SCQRPA violates the Ikeda sum rule. This phenomenon has been indicated in the previous studies , but the degree of violation we obtained is less than in the other calculations because we did not include all multipolarities.

ii) In the limit when the difference between proton and neutron quasiparticle occupation numbers is neglected the FR-QRPA coincides with SCQRPA.

iii) From a comparison of FR-QRPA with SCQRPA results we conclude that the effect of the restoration of the Ikeda sum rule is important in the range of large value of particle-particle strength beyond the point of collapse of the standard QRPA.

It is worth to mention that the FR-QRPA approach is sensitive to the precise evaluation of the proton and neutron quasiparticle occupation numbers. Due to the limitation of the approximate expression given in (28) (motivated by a similarity to the SQRPA approach) the convergence of the FR-QRPA is achieved only for relatively small model space. However, even for such a model space the differences among the standard QRPA, SCQRPA, FR-QRPA approach are evident for $`\kappa `$ close to the point the standard QRPA breaks up. There is a hope that for a proper ansatz of the RPA ground state the FR-QRPA approach can work also for a large model space. This is the subject of our further study.

###### Acknowledgements.

This work was supported in part by the Landesforschungsschwerpunktsprogramm Baden-Wuerttemberg ”Low Energy Neutrino Physics”. L.P. and V.R. would like to thank the Graduiertenkolleg ”Hadronen im Vakuum, in Kernen und Sternen” GRK683 and IKYDA02 project for support. The work of F. Š. was supported in part by the Deutsche Forschungsgemeinschaft (436 SLK 17/298) and by the VEGA Grant agency of the Slovac Republic under contract No. 1/0249/03.


format: text

text:
# On the complexity of curve fitting algorithms

## Abstract

We study a popular algorithm for fitting polynomial curves to scattered data based on the least squares with gradient weights. We show that sometimes this algorithm admits a substantial reduction of complexity, and, furthermore, find precise conditions under which this is possible. It turns out that this is, indeed, possible when one fits circles but not ellipses or hyperbolas.

In many applications one needs to fit a curve described by a polynomial equation

$$P(x,y;\mathrm{\Theta })=0$$

(here $`\mathrm{\Theta }`$ denotes the vector of unknown parameters) to experimental data $`(x_i,y_i)`$, $`i=1,\mathrm{},n`$. In this equation $`P`$ is a polynomial in $`x`$ and $`y`$, and its coefficients are either unknown parameters or functions of unknown parameters. For example, a number of recent publications  are devoted to the problem of fitting quadrics $`Ax^2+Bxy+Cy^2+Dx+Ey+F=0`$, in which case $`\mathrm{\Theta }=(A,B,C,D,E,F)`$ is the parameter vector. The problem of fitting circles, given by equation $`(xa)^2+(yb)^2R^2=0`$ with three parameters $`a,b,R`$, also arises in practice .

It is standard to assume that the data $`(x_i,y_i)`$ are noisy measurements of some true (but unknown) points $`(\overline{x}_i,\overline{y}_i)`$ on the curve, see  for details. The noise vectors $`e_i=(x_i\overline{x}_i,y_i\overline{y}_i)`$ are then assumed to be independent gaussian vectors with zero mean and a scalar covariance matrix, $`\sigma ^2I`$. In this case the maximum likelihood estimate of $`\mathrm{\Theta }`$ is given by the orthogonal least squares fit (OLSF), which is based on the minimization of the function

$$(\mathrm{\Theta })=\underset{i=1}{\overset{n}{}}d_i^2$$
(1)
where $`d_i`$ denotes the distance from the point $`(x_i,y_i)`$ to the curve $`P(x,y;\mathrm{\Theta })=0`$.

Under these assumptions the OLSF is statistically optimal – it provides estimates of $`\mathrm{\Theta }`$ whose covariance matrix attains its Rao-Cramer lower bound . The OLSF is widely used in practice, especially when one fits simple curves such as lines or circles. However, for more general curves the OLSF becomes intractable, because the precise distance $`d_i`$ is hard to compute. In those cases one resorts to various alternatives, and the most popular one is the algebraic fit (AF) based on the minimization of

$$_\mathrm{a}(\mathrm{\Theta })=\underset{i=1}{\overset{n}{}}w_i[P(x_i,y_i;\mathrm{\Theta })]^2$$
(2)
where $`w_i=w(x_i,y_i;\mathrm{\Theta })`$ are suitably defined weights. The choice of the weight function $`w(x,y;\mathrm{\Theta })`$ is important. The AF is known  to provide a statistically optimal estimate of $`\mathrm{\Theta }`$ (in the sense that the covariance matrix will attain its Rao-Cramer lower bound) if and only if the weight function satisfies

$$w(x,y;\mathrm{\Theta })=a(\mathrm{\Theta })/P(x,y;\mathrm{\Theta })^2$$
(3)
for all points $`x,y`$ on the curve, i.e. such that $`P(x,y;\mathrm{\Theta })=0`$. Here $`P=(P/x,P/y)`$ is the gradient vector of the polynomial $`P`$, and $`a(\mathrm{\Theta })>0`$ may be an arbitrary function of $`\mathrm{\Theta }`$ (in practice, one simply sets $`a(\mathrm{\Theta })=1`$). Any other choice of $`w`$ will result in the loss of accuracy, see . We call $`w(x,y;\mathrm{\Theta })`$ a gradient weight function if it satisfies (3) for all $`x,y`$ on the curve $`P(x,y;\mathrm{\Theta })=0`$. The AF (2) with a gradient weight function $`w(x,y;\mathrm{\Theta })`$ is commonly referred to as the gradient weighted algebraic fit (GRAF). It was introduced in the mid-seventies  and recently became standard for polynomial curve fitting, see, for example, .

Even though the GRAF is much cheaper than the OLSF, it is still a nonlinear problem requiring iterative methods. For example, in a popular reweight procedure  one uses the $`k`$-th approximation $`\mathrm{\Theta }^{(k)}`$ to compute the weights $`w_i=w(x_i,y_i;\mathrm{\Theta }^{(k)})`$ and then finds $`\mathrm{\Theta }^{(k+1)}`$ by minimizing (2) regarding the just computed $`w_i`$’s as constants. Note that if the parameters $`\mathrm{\Theta }`$ are the coefficients of $`P`$, then (2), with fixed weights, becomes a quadratic function in $`\mathrm{\Theta }`$, and its minimum can be easily found. Another algorithm is based on solving the equation $`_\mathrm{\Theta }_\mathrm{a}(\mathrm{\Theta })=0`$, i.e.

$$P_i^2_\mathrm{\Theta }w_i+2w_iP_i_\mathrm{\Theta }P_i=0$$
(4)
for which various iterative schemes could be used. In the case of fitting quadrics, for example, the most advanced algorithms are the renormalization method , the heteroscedastic error-in-variables method  and the fundamental numerical scheme . In all these algorithms, one needs to evaluate $`𝒪(n)`$ terms at each iteration. Therefore, the complexity of those algorithms is $`𝒪(kn)`$, where $`k`$ is the number of iterations. Moreover, each algorithm requires access to individual coordinates $`x_i,y_i`$ of the data points at each iteration. These difficulties can be sometimes avoided in a remarkable way, as we show next.

Suppose we need to fit circles given by equation

$$P(x,y)=(xa)^2+(yb)^2R^2=0.$$

Then we have

$$P(x,y;\mathrm{\Theta })^2=4(xa)^2+4(yb)^2=4P(x,y)+4R^2$$
(5)
hence $`P(x,y;\mathrm{\Theta })^2=4R^2`$ for all the points $`(x,y)`$ lying on the circle $`P(x,y)=0`$, and we can set $`w(x,y;\mathrm{\Theta })=1/R^2`$. Therefore

$`_\mathrm{a}(a,b,R)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}R^2\left[x_i^2+y_i^22ax_i2by_i+a^2+b^2R^2\right]^2`$ (6)
$`=`$ $`R^2[z_1+az_2+bz_3+a^2z_4+b^2z_5+abz_6+cz_7+acz_8+bcz_9+c^2n]`$
where we denoted $`c=a^2+b^2R^2`$ for brevity, and

$$z_1=(x_i^2+y_i^2)^2,z_2=4x_i(x_i^2+y_i^2),\mathrm{}$$

are some expressions involving $`x_i`$ and $`y_i`$ only.

The minimization of (6) is still a nonlinear problem requiring iterative methods , but it has obvious advantages over the reweight procedure described above and other generic methods for solving the equation (4). First of all, the values of $`z_1,\mathrm{},z_9`$ only need to be computed once, and then the cost of minimization of (6) will not depend on $`n`$ anymore. Thus, the complexity of this algorithm is $`𝒪(n)+𝒪(k)`$, where $`𝒪(n)`$ is the cost of evaluation of $`z_1,\mathrm{},z_9`$ and $`𝒪(k)`$ is the cost of some $`k`$ iterations spent on the subsequent minimization of $`_\mathrm{a}(a,b,R)`$. Moreover, once the values of $`z_1,\mathrm{},z_9`$ are computed and stored, the coordinates $`x_i,y_i`$ can be destroyed. Practically, $`z_1,\mathrm{},z_9`$ can be computed “on-line”, when the data are collected. The minimization procedure per se can be implemented “off-line”, without storage of (or access to) the data points. The quantities $`z_1,\mathrm{},z_9`$ here play the role of sufficient statistics.

Inspired by the above example, we might say that the problem of fitting a polynomial curve $`P(x,y;\mathrm{\Theta })=0`$ admits a reduction of complexity if there are $`\mathrm{}`$ functions $`z_j(x_1,y_1,\mathrm{},x_n,y_n)`$, $`1j\mathrm{}`$, with $`\mathrm{}`$ being independent of $`n`$ and $`\mathrm{\Theta }`$, and a gradient weight function $`w(x,y;\mathrm{\Theta })`$ such that

$$_\mathrm{a}=f(z_1,\mathrm{},z_{\mathrm{}};\mathrm{\Theta })$$
(7)
i.e. $`_\mathrm{a}`$ is a function of $`z_1,\mathrm{},z_{\mathrm{}}`$ and $`\mathrm{\Theta }`$ only.

This definition does not suggest how to find the functions $`z_1,\mathrm{},z_{\mathrm{}}`$ in practical terms, though. Since $`_\mathrm{a}`$ is given by (2) with $`P(x_i,y_i;\mathrm{\Theta })`$ being a polynomial in $`x_i,y_i`$, then the most natural (if not the only) way to construct the functions $`z_1,\mathrm{},z_{\mathrm{}}`$ is to express the gradient weight function (3) in the form

$$w(x,y;\mathrm{\Theta })=\underset{k=1}{\overset{K}{}}C_k(\mathrm{\Theta })D_k(x,y)$$
(8)
where $`C_k`$ are functions of the parameter vector $`\mathrm{\Theta }`$ alone, and $`D_k`$ are functions of $`x`$ and $`y`$ only (here the number of terms, $`K`$, must be independent of $`\mathrm{\Theta }`$). Indeed, suppose that the representation (8) is found. Since $`P^2`$ is a polynomial in $`x,y`$, we can expand it as

$$P^2(x,y)=\underset{p,q}{}c_{p,q}x^py^q$$

where $`c_{p,q}=c_{p,q}(\mathrm{\Theta })`$ denote its coefficients. Now the function $`_\mathrm{a}`$ can be evaluated as

$`_\mathrm{a}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{K}{}}}{\displaystyle \underset{p,q}{}}C_k(\mathrm{\Theta })c_{p,q}(\mathrm{\Theta }){\displaystyle \underset{i=1}{\overset{n}{}}}x_i^py_i^qD_k(x_i,y_i)`$
$`=`$ $`{\displaystyle \underset{k=1}{\overset{K}{}}}{\displaystyle \underset{p,q}{}}C_k(\mathrm{\Theta })c_{p,q}(\mathrm{\Theta })z_{k,p,q}`$
where

$$z_{k,p,q}=\underset{i=1}{\overset{n}{}}x_i^py_i^qD_k(x_i,y_i)$$

The values of $`z_{k,p,q}`$ depend on the data $`x_i,y_i`$ only, hence we obtain the desired representation (7). Therefore, (8) implies (7). We believe that the converse is also true, i.e. the conditions (7) and (8) are actually equivalent, but we do not attempt to prove that.

Motivated by the above considerations, we adopt the following definition: the problem of fitting a polynomial curve $`P(x,y;\mathrm{\Theta })=0`$ admits a reduction of complexity if the gradient weight function (3) can be expressed in the form (8).

As we have seen, the problem of fitting circles admits a reduction of complexity (and so does the simpler problem of fitting lines). Now if the problem of fitting ellipses and/or hyperbolas admitted a reduction of complexity as defined above, we would be able to dramatically improve the known GRAF algorithms . Unfortunately, this is impossible – there are deep mathematical reasons which prevent a reduction of complexity in the case of ellipses, hyperbolas, and parabolas.

In this paper we find general conditions on the polynomial $`P(x,y;\mathrm{\Theta })`$ under which the problem of fitting the curve $`P(x,y;\mathrm{\Theta })=0`$ allows a reduction of complexity. It turns out that lines and circles satisfy these conditions, but ellipses, hyperbolas, and parabolas do not. Our results thus demonstrate (in a rigorous mathematical way) that fitting noncircular conics is an intrinsically more complicated problem than fitting circles or lines.

For convenience, let us denote

$$Q(x,y;\mathrm{\Theta }):=P(x,y;\mathrm{\Theta })^2=(P/x)^2+(P/y)^2$$

Clearly, $`Q(x,y;\mathrm{\Theta })`$ is itself a polynomial in $`x`$ and $`y`$. Our subsequent arguments will involve some facts from complex analysis. We will treat $`x`$ and $`y`$ as complex, rather than real, variables.

Theorem. The problem of fitting curves $`P(x,y;\mathrm{\Theta })=0`$ admits a reduction of complexity (as defined above) under the condition that the system of polynomial equations

$$\{\begin{array}{c}P(x,y)=0\\ Q(x,y)=0\end{array}$$
(9)
has no solutions, real or complex.

Before we prove our theorem, we shall show how to use it. For the problem of fitting circles, we have already computed $`Q=4P+4R^2`$, see (5), hence the system (9) has indeed no solutions for nondegenerate circles (for which $`R0`$).

When using the theorem, the following invariance property will be helpful. Let $`(x,y)(\stackrel{~}{x},\stackrel{~}{y})`$ be a transformation of the $`xy`$ plane that is a composition of translations, rotations, mirror reflections and similarities (the latter are defined by $`(x,y)(cx,cy)`$ for some $`c0`$). Denote by $`\stackrel{~}{P}(\stackrel{~}{x},\stackrel{~}{y})`$ the polynomial $`P`$ in the new coordinates $`\stackrel{~}{x},\stackrel{~}{y}`$. Then the system (9) has a solution (real or complex) if and only if the corresponding system

$$\{\begin{array}{c}\stackrel{~}{P}(\stackrel{~}{x},\stackrel{~}{y})=0\\ \stackrel{~}{Q}(\stackrel{~}{x},\stackrel{~}{y})=0\end{array}$$

has a solution, real or complex. Here $`\stackrel{~}{Q}=\stackrel{~}{P}^2`$. This simple fact, which can be verified directly by the reader, allows us to simplify the polynomial $`P(x,y)`$ before applying the theorem.

Consider the problem of fitting ellipses and hyperbolas. By using a translation and rotation of the $`xy`$ plane we can always reduce the polynomial $`P`$ to a canonical form $`ax^2+by^2+c=0`$ (with $`ab`$ and $`abc0`$). Then $`Q=4a^2x^2+4b^2y^2`$ and we arrive at a system of equations

$$\{\begin{array}{c}ax^2+by^2+c=0\\ a^2x^2+b^2y^2=0\end{array}$$

It is easy to see that it always has a solution

$$x=\pm \sqrt{\frac{bc}{a(ab)}},y=\pm \sqrt{\frac{ac}{b(ab)}}$$

(note that $`x`$ or $`y`$ may be an imaginary number, which is allowed by our theorem). Therefore, the problem does not admit a reduction of complexity.

If our curve is a parabola, then we can use its canonical equation $`y=cx^2`$ for $`c>0`$, hence $`P=ycx^2`$ and $`Q=4c^2x^2+1`$. Here again we have a common zero of $`P`$ and $`Q`$ at the point $`x=𝐢/2c`$ and $`y=1/4c`$. Thus, no conic sections (except circles) satisfy the conditions of our theorem.

We now prove our theorem. Since $`w(x,y;\mathrm{\Theta })`$ must be a gradient weight function, the requirement (8) is equivalent to

$$\frac{1}{Q(x,y)}=\underset{k=1}{\overset{K}{}}C_k(\mathrm{\Theta })D_k(x,y)\mathrm{whenever}P(x,y)=0$$
(10)
(here we incorporated the factor $`a(\mathrm{\Theta })`$ into the coefficients $`C_k(\mathrm{\Theta })`$, for convenience). We emphasize that the left identity in (10) does not have to hold on the entire $`xy`$ plane, it only has to hold on the curve $`P(x,y)=0`$. If we denote that curve by $``$, then (10) can be restated as

$$\frac{1}{Q(x,y)}=\underset{k=1}{\overset{K}{}}C_k(\mathrm{\Theta })D_k(x,y)\mathrm{whenever}(x,y)$$
(11)
The functions $`D_k(x,y)`$ in (10) cannot be arbitrary, they must be easily computable, i.e. available in the machine arithmetics. That is, they must be combinations of elementary functions – polynomials, exponentials, logarithms, trigonometric functions, etc. In that case $`D_k(x,y)`$ are analytic functions of $`x`$ and $`y`$. Therefore, they have analytic extensions to the two-dimensional complex plane $`\text{C}|\text{ }|\text{ }^2`$. We note that they do not need be entire functions, i.e. analytic everywhere in $`\text{C}|\text{ }|\text{ }^2`$, they may have some singularities. For example, the function $`(1+x^2+y^2)^1`$ is analytic in $`\text{I}\text{R}^2`$ but has singularities in $`\text{C}|\text{ }|\text{ }^2`$, e.g. the point $`x=𝐢`$ and $`y=0`$ is its singularity. Also, those extensions maybe multivalued functions (examples are $`\mathrm{ln}x`$ or $`\sqrt{x}`$).

Now, the following function will also be analytic in $`\text{C}|\text{ }|\text{ }^2`$:

$$G(x,y)=1Q(x,y)\underset{k=1}{\overset{K}{}}C_k(\mathrm{\Theta })D_k(x,y)$$

since it is a combination of analytic functions. By (11), it vanishes on the curve $``$ in the real $`xy`$ plane. Consider the subset $`𝒵\text{C}|\text{ }|\text{ }^2`$ defined by the equation $`P(x,y)=0`$, where $`x`$ and $`y`$ are treated as complex variables. Note that $``$ is a curve on the two-dimensional manifold $`𝒵`$. We will prove that the function $`G(x,y)`$ vanishes on the entire $`𝒵`$.

We can assume that $`P(x,y)`$ is an irreducible polynomial (otherwise we can apply our argument to each irreducible factor of $`P`$). Then $`𝒵`$ is an algebraic variety, hence it admits a complex parametrization (a complex coordinate, $`z`$), and the restriction of the function $`G`$ onto $`𝒵`$ will be an analytic function of $`z`$. It is known in complex analysis that if an analytic function $`G(z)`$, $`z\text{C}|\text{ }|\text{ }`$, vanishes on a one-dimensional curve in C$`|`$ $`|`$ , then it is identically zero on C$`|`$ $`|`$ , hence $`G(z)0`$ for all $`z\text{C}|\text{ }|\text{ }`$. In our case the curve on which $`G`$ vanishes is $``$ (and we assume, of course, that it is a nondegenerate curve for all the relevant values of the parameter $`\mathrm{\Theta }`$). Hence, $`G`$ vanishes on the entire $`𝒵`$, and therefore

$$G(x,y)=0\mathrm{whenever}(x,y)𝒵$$
(12)
On the other hand, if the system of equations (9) has a complex solution $`(x,y)`$, then (12) would be impossible, since any solution of (9) lies on the manifold $`𝒵`$ (because $`P(x,y)=0`$), and at the same time $`Q(x,y)=0`$ implies $`G(x,y)=1`$. Therefore, if the system (9) has a solution (real or complex), then the representation (8) cannot possibly exist.

It remains to show that if the system (9) has no solutions, then the representation (8) is possible, and hence our problem indeed admits a reduction of complexity. Assuming that (9) has no solutions, we will construct the representation (8) in the simplest, polynomial form:

$$w(x,y;\mathrm{\Theta })=\underset{p,q}{}w_{p,q}(\mathrm{\Theta })x^py^q$$
(13)
the degree of this polynomial being independent of the parameter $`\mathrm{\Theta }`$. Consider a polynomial equation

$$P(x,y)U(x,y)+Q(x,y)W(x,y)=1$$
(14)
where $`U(x,y)`$ and $`W(x,y)`$ are unknown polynomials. A classical mathematical theorem, Hilbert’s Nullstellensatz , says that the equation (14) has polynomial solutions $`U(x,y)`$ and $`W(x,y)`$ if and only if $`P(x,y)`$ and $`Q(x,y)`$ have no common zeroes in $`\text{C}|\text{ }|\text{ }^2`$, i.e. whenever the system (9) has no complex solutions, which is exactly what we have assumed. Note that since $`P`$ and $`Q`$ depend on $`\mathrm{\Theta }`$, then so do $`U`$ and $`W`$, but we suppress this dependence in the equation (14).

Now the polynomial $`W(x,y)`$ solving (14) gives us the weight function $`w(x,y;\mathrm{\Theta })=W(x,y)`$, and it is easy to see that

$$W(x,y)=1/Q(x,y)\mathrm{whenever}P(x,y)=0$$

Technically, the theorem is proved, but we make a further practical remark. Suppose we know that the system (9) has no solutions, so that the problem admits a reduction of complexity. In this case we need to find the polynomial $`W(x,y)`$ solving (14) in an explicit form, in order to determine the weight function $`w(x,y;\mathrm{\Theta })`$. To this end we describe a finite and relatively simple algorithm for computing the coefficients $`w_{pq}`$ of the polynomial $`W`$. We substitute the expansions

$$W(x,y)=\underset{p,q}{}w_{p,q}x_i^py_i^q\mathrm{and}U(x,y)=\underset{p,q}{}u_{p,q}x_i^py_i^q$$

into the identity (14) and then equate the terms on the left hand side and those on the right hand side with the same degrees of the variables $`x,y`$. This gives a linear system of equations for the unknown coefficients $`w_{pq}`$ and $`u_{pq}`$. This might be a large system (its size depends on the degrees of $`U`$ and $`W`$), but it is a linear system whose solution can be found by routine matrix methods. If the assumed degrees of $`U`$ and $`V`$ are high enough, then the above system is always solvable by the so called effective Nullstellensatz, see . By solving that system we can obtain explicit formulas for the coefficients $`w_{pq}`$ and $`u_{pq}`$. In fact, we only need the coefficients of $`W`$, not $`U`$. Lastly, we remark that those coefficients will be rational functions of the coefficients of the polynomial $`P(x,y)`$, hence they will be easily computable.

Acknowledgement. N. Chernov is partially supported by NSF grant DMS-0098788 and N. Simányi is partially supported by NSF grant DMS-0098773.


text:
# Picosecond time-resolved two-dimensional ballistic electron transport

## Abstract

Time-resolved transport of ballistic electrons in a two-dimensional electron gas has been measured with a resolution of less than $`5ps`$. This was accomplished by using picosecond electrical pulses to launch electrons from the emitter of a transverse magnetic focusing structure and optoelectronically sampling the collector voltage. Both plasma resonances and the ballistic transport signal are clearly resolved. The transit time appears to be somewhat longer than expected from simple Fermi velocity considerations.

Time-dependent transport measurements in two-dimensional electron gas (2DEG) systems facilitate investigations of the charge distribution and its collective dynamicsButtiker (2000). Techniques implementing fast sampling oscilloscopes have been used to detect edge-magnetoplasmon excitations in the time domain at temperatures as low as $`0.3K`$ and magnetic fields as high as $`13T`$Ashoori et al. (1992); Ernst et al. (1996). Typical time resolution achieved in these experiments was on the order of $`100ps`$. To measure the dynamic response of a mesoscopic device, a topic of significant theoretical interestButtiker (2000); Zozoulenko and Blomquist (2003), one encounters additional experimental difficulties. Ballistic electron transport times in mesoscopic devices fabricated from AlGaAs/GaAs 2DEGs can easily range from a few to hundreds of picoseconds. To fully characterize the device response in this regime, it is advantageous to provide picosecond scale excitation and detection while still working in the environment of low temperatures and high magnetic fields. In this letter, we report the direct measurement of picosecond time-resolved ballistic electron transport in a 2DEG transverse magnetic focusing device. Our approach utilizes an ultrafast waveguide and integrated detector coupled to the device under test.

For high frequency investigations of 2DEG transport, it is convenient to have a way to distinguish the electron signal from other coupling between the device nodes. Previous experiments investigating edge-magnetoplasmons have used gates to selectively deplete the 2DEG and thus determine its contribution to the overall signalErnst et al. (1996). In our experiment we chose to work in a transverse magnetic focusing geometry as shown in Fig.1a. In a conventional experiment (DC), a small current between contacts $`(1)`$ and $`(2)`$, is passed through the emitter point contact (labelled E in Fig.1b) and the voltage between contacts $`(3)`$ and $`(4)`$ is sensed across the collector point contact (labelled C in Fig.1b). By applying a perpendicular magnetic field, ballistic electrons injected into the base of the device can be deflected, by the Lorentz force, into the collector either directly or via skipping trajectories along the device boundary. As the magnetic field is swept, this leads to a periodic change in collector voltage with peaks at fields meeting the condition $`B=nB_{focus}`$ for integer n where $`B_{focus}=2\mathrm{}k_F/eL`$van Houten et al. (1989). Here $`L`$ is the point contact spacing and $`k_F`$ is the Fermi wave vector determined by the 2DEG density. This signal is only observed in one direction of magnetic field as the opposite polarity deflects electrons away from the collector. The periodicity in B along with the one sided spectrum provides a means to distinguish the ballistic electron signal from the background.

The focusing device of Fig.1a was fabricated from an $`Al_{0.3}Ga_{0.7}As/GaAs`$ modulation doped heterostructure. The 2DEG formed at the interface was located $`700\AA `$ below the sample surface. The device was patterned using electron beam lithography and a shallow etch to define the point contacts ($`0.7\mu m`$ wide with $`3.7\mu m`$ center-to-center spacing). Alloyed Ni-Au-Ge contacts were made followed by a $`1000\AA `$ deep mesa etch to remove the 2DEG outside the region of interest. This device isolation was essential as bulk 2DEG magnetoplasmons are known to couple strongly to waveguide structuresShaner and Lyon (2002). The processed device was without detectable carriers before illumination and had an electron density of $`2.3\times 10^{11}cm^2`$ after illumination. The estimated mean free path of the electrons in this device was $`9\mu m`$. The illumination in this case originated from stray laser light (used for driving Auston switches) striking the device mesa. When the laser was turned off, the carrier density remained stable and there were no noticeable changes when it was left on for the duration of the experiment.

A separate waveguide device was fabricated that utilized low temperature grown GaAs (LT-GaAs) Auston switches Auston (1983) for both pulse generation and detection. Due to the short ($`1ps`$) recombination time in LT-GaAs, it is an excellent photoconductor for ultrafast applicationsKorn et al. (2002). The waveguide conductors were made using $`200\AA `$ Ti $`/2000\AA `$ Au in a liftoff process. The launch (L) and detection (D) points on the waveguide are shown in Fig.1c, along with the relevant device dimensions. The waveguide substrate was a multilayer structure comprised of a $`2\mu m`$ LT-GaAs epilayer transferred onto a $`200\mu m`$ thick sapphire plate that was subsequently bonded to silicon. Holes were etched completely through the silicon to hold single mode optical fibers that illuminate the switchesShaner and Lyon (2002). Indium bumps of $`2\mu m`$ height were evaporated on the connection points (enumerated in Fig.1 on both devices). The pieces were then brought together using a flip-chip bonder to make the electrical connections between the waveguide and the focusing device. Epoxy was used to fill the remaining gap of $`23\mu m`$ between the two devices.

Pump and probe optical fibers were attached to the back of the waveguide in order to illuminate the launch and detection switches. The sample was mounted inside a small superconducting magnet and measurements were performed at $`4.2K`$ in liquid helium. When a DC bias was applied across the launch switch, as shown in Fig.1c, an optical pulse sent down the pump fiber generated a lifetime limited photocurrent in the LT-GaAs. This current produced a voltage transient on the waveguide. Using a fiber dispersion compensator and $`250fs`$ wide $`790nm`$ wavelength optical pulses produced by a mode locked Ti:Sapphire laser, the minimum pulse width attained on a similar calibration waveguide (same conductor spacing and substrate) was $`1.8ps`$. On the detection side, an optical pulse received from the probe fiber electrically closed the detector switch and sampled the collector voltage with similar lifetime limited resolution.

To perform measurements, a DC bias of $`15V`$ was applied across the launch switch while one end of the center waveguide conductor (connected to contact (1) at the other end) and contact (2) were held at signal ground. The pump beam was chopped at a frequency of $`11Hz`$ allowing for differential lock-in detection of the voltage between the output of the detector switch and the opposing transport connection (4). An optical delay line was inserted into the pump beam path and the delay between probe and pump was stepped by $`1ps`$ between magnetic field sweeps from $`330mT`$ to $`330mT`$. Time averaged optical powers of $`300\mu W`$ and $`1mW`$ were used for the pump and probe respectively.

Zero delay (t=0) denotes the time when the electrical pulse generated on the waveguide reaches contact (1) of the focusing device. Its location (delay line position) was determined by careful measurements of the optical fiber lengths. At the conclusion of this series of experiments, the fibers were removed from the device and the location of zero delay confirmed in the calibration waveguide. Figure 2 shows magnetic field sweeps at selected delay times. For times up to $`t=30ps`$, no focusing signal can be seen. At $`t=50ps`$, the focusing signal is clearly developed and the focusing oscillations are seen to decay in amplitude for longer delay times.

The time dependence of the collector voltage at zero magnetic field is shown in Fig.3a. An initial pulse is seen in the collector response that coincides with the expected arrival time of the electrical pulse at the emitter to within $`1ps`$. We interpret this signal as coupling between the emitter and collector nodes. Although this coupling interferes with extracting a transport signal, it is useful in that it gives both an indication of the time dependence of the emitter voltage and a reference point in the data for when the electrical pulse reaches the focusing device. Even though pulses with $`1.8ps`$ width are launched on the waveguide, this signal is broadened to $`5ps`$ which serves as an upper limit on the exciting pulse duration. Also, because the waveguide did not have an impedance matched termination at the left handed end (Fig.1c), reflections are seen at later times which could launch additional electrons. However, there is a window of about $`50ps`$ in which there are no major reflections. The structure in this zero field response is reproducible at all times (within $`5\mu V`$ noise level - about the thickness of the trace in Fig.3a) and the signal at times earlier than those displayed was zero.

To further isolate the transport response, the raw data was filtered so that only those parts of the response which vary with magnetic field remain. This removes the coupling signal that is largely independent of magnetic field (a constant offset of each field scan). Filtered scans are assembled in the greyscale plot of Fig.3b with time along the x-axis and magnetic field along the y-axis. Here, the light areas correspond to electrons entering the collector. The ballistic electron signal is quite pronounced in the form of horizontal white streaks at the focusing peak locations in magnetic field. This is further emphasized by the time averaged plot to the right of the main image. A weakly damped magnetoplasmon oscillation is also seen in the time response beginning at about $`t=70ps`$. This resonance appeared in both positive and negative magnetic field and its frequency varied from $`50GHz`$ at zero field, to approximately $`85GHz`$ at $`B=175mT`$. Although not displayed in the data of Fig.3b, these oscillations decayed away after approximately $`250ps`$. A less pronounced set of similar but higher frequency oscillations can be seen from $`t=20ps`$ to $`t=40ps`$.

The time response of the ballistic electrons travelling directly from the emitter to collector, displayed as the upper trace in Fig.3c, was obtained by averaging a $`5mT`$ wide horizontal slice at the first focusing maximum of the filtered response in Fig.3b at $`B=0.043T`$ (indicated by the arrow). A similar slice of the filtered data, taken at $`B=0.043T`$ is shown as the lower trace in Fig.3c. We interpret the absence of any structure in the time response at this negative field to mean that the peaks in the ballistic electron response are due to changes in the emitter output, and not an artifact of the magnetoplasmon resonance modulating the signal.

The first peak of the ballistic electron time response is centered at $`t=52ps`$, with an onset at about $`30ps`$. Given the $`50ps`$ window before the arrival of waveguide reflections, the time-resolved ballistic electron signal for $`t<80ps`$ ($`50ps`$ after the detection of the first electrons) will not be affected by the reflections. The structure in the ballistic electron signal (Fig.3c) around $`120ps`$ and $`170ps`$ may arise from these waveguide reflections. However, the first peak in Fig.3c occurs at too early a time to be affected and we attribute it to the ballistic electrons launched by the electrical pulse at $`t=0`$.

The width of the first ballistic electron peak in Fig.3c is approximately $`35ps`$, which is significantly broader than the initial electrical pulse. Group velocity dispersion was calculated to be only $`2ps`$ for a $`1ps`$ pulse injected into the transport region. A wide point contact produces an electron beam with approximately a $`cos(\alpha )`$ angular distribution where $`\alpha `$ is the angle from perpendicular injectionBoggild et al. (1999). Despite this, at fields meeting the focusing condition a high percentage of the electrons converge on the collectorBoggild et al. (1999). These electrons, launched at different angles, traverse paths of different lengths and result in a distribution of transit times. A classical analysisBoggild et al. (1999) was used to model the resulting time response to a $`1ps`$ emitter current pulse. Assuming all carriers move at the Fermi velocity, and using the geometry of Fig.1b, a response width of approximately $`18ps`$ was predicted. While this simple model does not fully explain the width of the measured response, it nonetheless indicates that angular dispersion of the emitter point contact beam contributes heavily to the overall response width.

The modeled ballistic response was peaked at $`28ps`$ with an onset at $`17ps`$. However, in Fig.3c the ballistic electron signal is instead seen to be peaked at $`52ps`$. In Fig.2 and Fig.3b it is also clear that we see no magnetic focusing until $`t>30ps`$. Precise agreement with the calculated $`28ps`$ base transit time is not expected since there are other possible delays. For example, an additional delay of up to $`10ps`$ could be expected in order to traverse the point contacts (including a reduced velocity in the point contact). At this time we do not have a full explanation of why these time-resolved measurements give a significantly longer delay than that expected from the simple Fermi velocity considerations.

To investigate the effects of electrical pulse amplitude, the device of Fig.1 was also operated in the reverse direction by biasing the detector switch and exciting it with the pump beam to make it function as a pulse launcher. In this configuration, the center waveguide conductor and one of the ground planes (connections 1 and 2 respectively in Fig.1c) form direct electrical connections to the focusing device and were used to measure the time averaged focusing signal. Figure 4 compares the focusing spectrum measured using both continuous wave (CW) and pulsed laser excitation. A current of approximately $`100nA`$ was measured as being drawn from the DC voltage source at $`100\mu W`$ laser power (either CW or pulsed).

In general, the focusing peaks in the pulsed case degraded both in contrast and amplitude at higher laser powers. This is a signature of increased scattering due to electron heatingHornsey et al. (1993) and can be reproduced in similar magnetic focusing devices (using conventional CW current drive) by running them at high currents. It was found that reducing the switch bias at fixed optical power had the same effect as lowering the optical power at fixed switch bias. This confirmed that the heating effects were directly dependent on electrical pulse amplitude (a combination of switch bias and pump laser power) and are not arising from laser heating.

With the hot electron effects in mind, measurements were taken in the time-resolved configuration (Fig.1) at pump laser powers as low as $`20\mu W`$ (due to differences in switch efficiency, $`300\mu W`$ on the launch switch corresponded to approximately $`80\mu W`$ on the detection switch). There was no change in the time response of the focusing signal as a function of launching power over this range. This demonstrated that, compared to other broadening mechanisms, the emitter pulse voltage amplitude did not contribute heavily to the ballistic electron response of Fig.3c.

We have directly time-resolved ballistic electron transport in a transverse magnetic focusing device (fabricated from an $`Al_{0.3}Ga_{0.7}As/GaAs`$ heterostructure 2DEG) with $`<5ps`$ time resolution. Picosecond electrical pulses, delivered to the device by a waveguide, injected electrons through the emitter point contact and their arrival at the collector was time-resolved optoelectronically. A clear ballistic transport signal was observed along with magnetoplasmon oscillations in the device. The measured transport time was $`52ps`$, nearly twice as long as the expected time-of-flight for ballistic electrons travelling at the Fermi velocity through the focusing orbit. Additional delays inherent in the device have not been measured independently and the time difference cannot be completely explained at this time. The measured response width of the ballistic electron signal was found to be approximately $`35ps`$. After investigating various sources of response broadening, including group velocity dispersion and hot electron injection, the angular dispersion of the emitter point contact was found to contribute most heavily to the response width (accounting for approximately $`18ps`$ in calculations). Other sources of response broadening have not yet been identified.

###### Acknowledgements.

E.A. Shaner was supported by a National Defense Science and Engineering graduate fellowship. This work was supported in part by the U.S. Army Research Office under grant $`\mathrm{\#}`$DAAG55-98-1-0270.


text:
# The origin of the difference in the superconducting critical temperatures of the 𝛽_"H" and 𝛽_"L" phases of (BEDT-TTF)2I3

## Abstract

Incommensurate lattice fluctuations are present in the $`\beta _\text{L}`$ phase ($`T_c1.5`$ K) of ET2I3 (where ET is BEDT-TTF - bis(ethylenedithio)tetrathiafulvalene) but are absent in the $`\beta _\text{H}`$ phase ($`T_c7`$ K). We propose that the disorder in the conformational degrees of freedom of the terminal ethylene groups of the ET molecules, which is required to stabilise the lattice fluctuations, increases the quasiparticle scattering rate and that this leads to the observed difference in the superconducting critical temperatures, $`T_c`$, of the two phases. We calculate the dependence of $`T_c`$ on the interlayer residual resistivity. Our theory has no free parameters. Our predictions are shown to be consistent with experiment. We describe experiments to conclusively test our hypothesis.

It has long been known that the application of hydrostatic pressure, $`P`$, to $`\beta `$-ET2I3 has a dramatic effect on the superconducting critical temperature, $`T_c`$. At ambient pressure $`T_c1.5`$ K but when the applied pressure reaches $`P1`$ kbar a discontinuous increase in $`T_c`$ ($`7`$ K) is observed. The low $`T_c`$ state ($`P1`$ kbar) is denoted the $`\beta _\text{L}`$ phase and the high $`T_c`$ state ($`P1`$ kbar) is labelled the $`\beta _\text{H}`$ phase. When the pressure on the $`\beta _\text{H}`$ phase is decreased the material does not return to the $`\beta _\text{L}`$ phase but rather $`T_c`$ is seen to further increase. Below $`T130`$ K the resistivity of the $`\beta _\text{H}`$ phase is found to undergo a discontinuous decrease while no such anomaly is found in the $`\beta _\text{L}`$ phase Ginodman . Incommensurate lattice fluctuations have been observed in the $`\beta _\text{L}`$ phase but they are absent in the $`\beta _\text{H}`$ phase below $`T130`$ K Emge . The incommensurate lattice fluctuations are stabilised by variations in the conformational ordering of the terminal ethylene groups of the ET molecules and thus can only exist in the presence of disorder Ravy . For a recent review of this phenomenology see Ishiguro .

Although the change in resistivity between the $`\beta _\text{L}`$ and $`\beta _\text{H}`$ phases and the incommensurate lattice fluctuations have been studied both theoretically Ravy  and experimentally Ginodman ; Emge ; Ishiguro , no explanation of the change in $`T_c`$ has been forthcoming. In this paper we propose that difference in the $`T_c`$’s of the two phases is due to disorder, we will show that current data is consistent with our explanation and suggest experiments which could clearly determine whether or not this is the correct explanation of the difference in $`T_c`$.

In a recent paper it was shown BenRHMdisorder  that intrinsically non-magnetic disorder decreases the $`T_c`$ of a wide range of the superconducting salts of ET in line with the Abrikosov–Gorkov (AG) formula:

$`\mathrm{ln}\left({\displaystyle \frac{T_{c0}}{T_c}}\right)=\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\mathrm{}}{4\pi k_BT_c}}{\displaystyle \frac{1}{\tau }}\right)\psi \left({\displaystyle \frac{1}{2}}\right),`$ (1)
where $`T_{c0}`$ is the superconducting transition temperature of a pure sample, $`1/\tau `$ is the quasiparticle scattering rate and $`\psi (x)`$ is the digamma function. There are two scenarios compatible with this observation: either there is opposite spin pairing and the impurities induce localised magnetic moments, or else there is a finite angular momentum pairing state (most probably d-wave pairing). In this paper we will not discuss which of these scenarios is realised, but rather make use of this observation of the effect that increasing $`1/\tau `$ has on $`T_c`$.

In a quasi-two dimensional metal, such as $`\beta `$-ET2I3, the residual interlayer resistivity is given by

$`\rho _0={\displaystyle \frac{\pi \mathrm{}^4}{2e^2m^{}ct_{}^2}}{\displaystyle \frac{1}{\tau }}`$ (2)
where $`m^{}`$ is the quasiparticle effective mass, $`c`$ ($`=15.291`$ Å for $`\beta `$-ET2I3 Ishiguro ) is the interlayer lattice constant and $`t_{}`$ is interlayer hopping integral. Substituting (2) into (1) one finds that

$`\mathrm{ln}\left({\displaystyle \frac{T_{c0}}{T_c}}\right)=\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{e^2m^{}ct_{}^2}{2\pi ^2\mathrm{}^3k_BT_c}}\rho _0\right)\psi \left({\displaystyle \frac{1}{2}}\right).`$ (3)
from quantum oscillation Wosnitza et al. Wosnitza\_ratio\_kIBr  found that $`m^{}/m_e=4.2\pm 0.2`$, where $`m_e`$ is the electronic rest mass, $`t_{}/E_F=1/(175\pm 10)`$, where $`E_F`$ is the Fermi energy, and the Fermi wavevector, $`k_F=(3.4\pm 0.2)\times 10^9`$ cm-1. Thus taking $`E_F=\mathrm{}^2k_F^2/2m^{}`$ one finds that $`t_{}=0.59\pm 0.08`$ meV. $`T_{c0}=7.75`$ K is found from fast electron irradiation experiments BenRHMdisorder ; Forro .

In figure 1 we plot $`T_c`$ as a function of $`\rho _0`$. We stress that there are no free parameters in this plot once the data from Wosnitza et al. and Forro et al. BenRHMdisorder ; Forro  is considered. For comparison we also plot the data of Tokumoto et al. Tokumoto  who deliberately induced impurities in $`\beta _\text{H}`$-ET2I3 by fabricating the alloy $`\beta `$-ET2(I3)1-x(IBr2)x. The data is consistent with our theory although clearly more data is required to properly test our prediction.

We propose that the disorder in the conformational degrees of freedom of the terminal ethylene groups of the ET molecules which is required to stabilise incommensurate lattice fluctuations in the $`\beta _\text{L}`$ phase Ravy , but absent in the $`\beta _\text{H}`$ phase, increases the quasiparticle lifetime. (Either by inducing localised magnetic moments or causing variation in the site energy, depending on which of the scenarios proposed in BenRHMdisorder  is realised). Note that conformational disorder of the terminal ethylene groups of $`\kappa `$-ET2Cu\[N(CN)2\]Br can be controlled by varying the cooling rate and this leads to variations in $`T_c`$ which are well described by the AG formula BenRHMdisorder .

Unfortunately, we are not aware of any reports of the resistivity of $`\beta _\text{L}`$-ET2I3 with which to test our hypothesis. For example Ginodman et al. Ginodman  reported the change in resistance observed between the $`\beta _\text{H}`$ and $`\beta _\text{L}`$ phases, but did not report the size of their crystals. Taking the relative dimensions of their sample as a fitting parameter we find good agreement with our prediction (see figure 1; based on $`A/d=0.4`$ cm, where $`A`$ is the cross-sectional area of the sample and $`d`$ is the length of the sample in the c-axis, this is actually a fairly typical relative dimension for samples of this material). However, as we do not know the size of the sample this does not represent a very stringent test of our hypothesis. The resistance data of Ginodman et al. is well described by the Fermi liquid form $`R(T)=R_0+AT^2`$ in both the $`\beta _\text{L}`$ and $`\beta _\text{H}`$ phases with the same value of $`A`$ in both phases. This suggests that major difference between the two phases is the scattering rate due to impurities, in support of our hypothesis.

Systematic measurements of the resistivity of $`\beta `$-ET2I3 in both the $`\beta _\text{L}`$ and $`\beta _\text{H}`$ phases are required to test our hypothesis. If this were done for several samples with varying amounts of disorder due to random imperfections produced during the fabrication process then it would be possible to map out the entire curve of figure 1 for the $`\beta _\text{H}`$ phase. Clearly only the lower part of the curve could be mapped out for the $`\beta _\text{L}`$ phase because of the disorder in the terminal ethylene groups.

Acknowledgements

I would like to thank Ross McKenzie and Jochen Wosnitza for useful conversations. This work was supported by the Australian Research Council.


text:
# Anomalous Evolution of the Dwarf Galaxy HiPASS J1321-311footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., (AURA), under NASA Contract NAS 5-26555.

## 1 Introduction

Dwarf galaxies are the most common type of galaxy in the Universe. In cosmological models, such as those invoking cold dark matter, the first and most frequently formed gravitationally bound objects resemble present-day dwarfs. Many studies have been made of these galaxies, allowing us a glimpse into the conditions and evolution of galaxies in the youthful Universe. While the stellar populations of dwarf galaxies in the Local Group have been extensively investigated (see Mateo 1998 and references within; Grebel 1997), this represents only one environment in which these types of galaxies are known to exist. Still open is the question of how, and if, dwarf galaxies evolve from initial gas clouds to dwarf irregular galaxies to dwarf spheroidal galaxies (e.g., Grebel, Gallagher, & Harbeck 2003). Furthermore, we need to establish whether dwarf galaxies outside the Local Group have similar relationships between star formation histories and other global and environmental properties.

Numerous studies have explored the global properties of dwarf galaxies in nearby galaxy groups and clusters. However, detailed investigations of star formation histories derived from resolved stellar populations are limited to the nearest galaxy groups, which require the capabilities of the Hubble Space Telescope (HST). Thus we now have Wide-Field Planetary Camera 2 (WFPC2) color-magnitude diagrams (CMDs) for significant samples of dwarf galaxies in the M81 (e.g., Karachentsev et al. 2002a) and the Centaurus A (Cen A; Karachentsev et al. 2002b) groups. These and similar data suggest diverse star formation histories among nearby dwarfs, similar to those found in the Local Group. Dwarf spheroidal and transition galaxies show dominant red giant branches (RGBs), while the luminous stellar populations in irregulars consist of a mixture of the RGB stars along with more massive stars and their evolved descendants. Galaxies where the old RGB is weak are rare, but this condition is seen in the dwarf irregular galaxies Leo A (Tolstoy et al. 1998; Schulte-Ladbeck et al. 2002) and Holmberg IX (Karachentsev et al. 2002a; Marakova et al. 2002). This suggests that larger samples of resolved dwarf galaxies will reveal rarer or shorter lived evolutionary paths.

One property that makes the dwarf galaxies in Cen A interesting to study is the fact that so many of them are gas-rich. No dwarfs in the Local Group are young in the sense that they still contain more mass in gas than in stars. However, some of those that have been found in a blind Hi survey of the nearby Cen A group are (Hi Parkes All-Sky Survey \[HiPASS\]; Banks et al. 1999). After a number of new candidate dwarf galaxies were discovered, we undertook follow-up observations of a select number of them using the HST/WFPC2, giving us the opportunity to study the resolved stellar populations and thus the evolutionary histories in gas-rich Cen A group dwarfs. In this letter we discuss the unusual nature of one of these systems, HiPASS J1321-31. In a future paper we discuss its physical properties along with the other Cen A dwarfs observed is this program, HiPASS J1337-39 and HiDEEP J1337-33 (Grossi et al. 2003). We find that the color-magnitude diagram (CMD) of HiPASS J1321-31 shows a significantly different stellar population from the other two dwarf galaxies surveyed, as well as most Local Group dwarfs. In the following we present the CMD of HiPASS J1321-31 and discuss the possible origins of its “unique” stellar populations.

## 2 Observations and Reductions

HiPASS J1321-31 was observed with the WFPC2 on 2001 June 12 with four exposures in F555W at 1200 sec each and four exposures in F814W at 1300 sec each, with the galaxy centered on the WFC3 (Figure 1). The images for each filter and chip were combined to create a deeper image of the galaxy. daophot, allstar, and allframe were used to obtain the profile-fitted photometry of the stars presented in this letter. The point spread functions (PSFs) used in the reductions were kindly communicated to us by P. B. Stetson (see Stetson et al. 1998). Due to the lack of stars on the PC, we chose not to reduce that chip. Limits in $`\chi `$ ($`\chi >2.0`$) and sharpness ($`0.1<\mathrm{sharp}<0.1`$) were placed on the fits of the stars to remove any objects with questionable photometry. CTE corrections were made using the equations from Whitmore, Heyer, & Casertano (1999). After the aperture corrections were calculated from the unsaturated brightest stars in each chip, the photometry was calibrated to the standard $`V`$ and $`I`$ system using the equations in Holtzman et al. (1995).

To check the reliability of the photometry and confirm that the red plume apparent in the CMD was “real,” the data were also reduced using DoPHOT and independently using daophot/allstar with different PSFs. In both cases, the red plume was still present. We also checked the location of the red plume stars against any defects in the WFC3, such as warm pixels, and found the photometry was unaffected by these difficulties.

## 3 Color-Magnitude Diagram

Figure 2a shows the CMD for the WFC3, which was centered on HiPASS J1321-31. The CMD for the field (WFC2 and WFC4) is seen in Figure 2b. Figure 2a shows a clear red plume of luminous stars. These stars are scattered across the optical extent of the galaxy (see Figure 1). We know that HiPASS J1321-31 has a heliocentric velocity of 572 km/s according to the Hi survey of Banks et al. (1999). This is consistent with it being a member of the Cen A Group and proves to be important in interpreting the CMD. In the following we discuss three possible interpretations for the red plume: A RGB, an asymptotic giant branch (AGB), or a red supergiant branch (RSGB).

### 3.1 The Red Plume

One might think that the red plume seen in Figure 2a is the RGB of HiPASS J1321-31. Assuming this to be true, the tip of the red giant branch (TRGB) would be at $`m_I22.6`$ mag. In this case using this value for the TRGB, we estimate the distance to the galaxy to be about 2.0 Mpc, assuming $`M_{I,\mathrm{TRGB}}=4.05`$ (Da Costa & Armandroff 1990) and $`E(\mathrm{B V})=0.062\pm 0.01`$ (Schlegel, Finkbeiner, & Davis 1998). This estimate places HiPASS J1321-31 just outside the Local Group in the direction of the Cen A group of galaxies. The distance to the giant elliptical galaxy Cen A is $`3.9\pm 0.3`$ Mpc (Harris, Harris, & Poole 1999), but HiPASS J1321-31 is located closer in the sky to M83, which has a distance of $`4.5\pm 0.3`$ Mpc (Thim et al. 2003). This interpretation would place HiPASS J1321-31 somewhere between the Milky Way and M83. This is in conflict with the radial velocity found for the Hi cloud associated with the galaxy, which is consistent with it being a member of the more distant Cen A group. This apparent conflict could be resolved if the Hi cloud were not associated with the stellar galaxy, but such a chance coincidence would be unique and is highly unlikely.

However, the biggest challenge to the RGB interpretation is the fact that the red plume is not well populated even well below its tip. For globular clusters and dwarf galaxies, the number of stars increases as one moves fainter along the RGB. Figure 3 shows the luminosity function for HiPASS J1321-31. The luminosity function has been Gaussian-smoothed using Eq. A1 in Sakai, Madore, & Freedman (1996). In the magnitude range from $`24.0<m_I<22.0`$, the luminosity function is relatively flat. We used the data in Saviane, Held & Piotto (1996) to create a luminosity function for the RGB in Tucana (the dashed line in Figure 3). The tip of the Tucana RGB was shifted to match the tip of the red plume in HiPASS J1321-31 and was scaled so that the upper luminosity functions are reasonably well-matched. The relatively constant luminosity function of the HiPASS J1321-31 red plume is in contrast to the increasing luminosity function of the Tucana RGB. The combination of distance difficulties and flat luminosity function argue that the red plume is not a RGB.

A second possibility is that the red plume is an AGB. A strong AGB population is seen in dwarf galaxies like IC 1613 (Cole et al. 1999) and Sextans A (Dohm-Palmer et al. 2002). In these galaxies, the AGB extends around 0.5 mag in $`I`$ above the TRGB and then turns redward. This differs from the behavior of the red plume in HiPASS J1321-31. It extends at least two magnitudes above any possible TRGB and begins its assent from the bluer part of the possible RGB. Thus, it is implausible that the red plume consists of AGB stars.

Ruling out the RGB and AGB leaves us with the likelihood that these stars make up the RSGB. Such stars originate from intermediate mass stars in post-main sequence helium core burning evolutionary phase, the same phases that produce pronounced “blue loops” at low metallicities. A beautiful example of this phenomenon is seen in Sextans A, where stars with ages of $`<1`$ Gyr mark the red and blue edges of loop evolution (Dohm-Palmer et al. 1997). However, for cases where the evolved stars have ages near 1 Gyr the red branch of the evolutionary loop dominates (see Fig. 3 of Dohm-Palmer et al.). Thus in this model HiPASS J1321-31 would be in something resembling a post burst phase where the star formation rate was enhanced $`1`$ Gyr in the past, giving rise to an unusual number of evolved stars with initial masses slightly less than $`3M_{}`$.

We searched the literature for other dwarf galaxies with similar CMDs as HiPASS J1321-31. In the HST snapshot survey of dwarf galaxies (Karachentsev et al. 2002a, b, c; 2003a, b, c) there are several with red plumes similar to HiPASS J1321-31. They are all found to extend from the bluer part of the RGB and trend slightly to the red as one moves to higher luminosity along the branch. We investigated the number of stars along the red plumes in a couple of these galaxies and found that they remain relatively constant, similar to what we see in the luminosity function for HiPASS J1321-31 (see Figure 3). Therefore, as noted by Karachentsev et al., it is likely these red plumes are RSGBs, including the one in HiPASS J1321-31, arising from intermediate mass stars.

We can test the intermediate mass RSG hypothesis by seeing if the presence of a normal RGB is consistent with our data. Following the approach adopted in the Karachentsev et al. papers, we used the edge-detection filter equation in Sakai, Madore, & Freedman (1996) to search for the TRGB. Applying the equation to those stars between $`0.8<VI<2.0`$ results in an edge being detected at $`m_I=24.65\pm 0.11`$ mag, where the error includes the photometric uncertainty ($`\pm 0.06`$), the uncertainty in the TRGB estimate ($`\pm 0.06`$), the uncertainty of the photometric zeropoint ($`\pm 0.05`$), and the uncertainty in the aperture corrections ($`\pm 0.05`$). This agrees with the visual impression of the point where the luminosity function begins to rise in Figure 3, and is a clear signature of the beginning of the RGB. Assuming the TRGB to be at $`M_I=4.05`$ (Da Costa & Armandroff 1990), $`E(\mathrm{B V})=0.062\pm 0.01`$ (Schlegel, Finkbeiner, & Davis 1998), and $`R_V=3.1`$, the distance modulus estimated from the TRGB is $`(mM)_0=28.59\pm 0.13`$. The distance to HiPASS J1321-31 is thus $`5.2\pm 0.3`$ Mpc. This would place the galaxy beyond M83, but it would still be a member of the Cen A group, in good agreement with the heliocentric velocity of the Hi cloud.

## 4 What Does This Mean?

What are the astrophysical implications of such a pronounced RSGB in HiPASS J1321-31 and other dwarfs? The difficulty in understanding the RSGB in HiPASS J1321-31 is the apparent deficiency of corresponding BSGB stars. Dwarf galaxies such as Sextans A, IC 1613, and those observed in the HST snapshot survey for dwarf galaxies, all exhibit a strong population of blue stars reaching nearly as high as the RSGB. These blue stars are made up of an uncertain combination of upper main sequence and evolved intermediate mass stars with ages $`<300`$ Myr. The comparable luminosity of the red and blue plume in Sextans A and IC 1613 then reflects relatively high levels of ongoing star formation. The absence of bright blue stars in HiPASS J1321-31 cannot be due to observational problems, for we see them in the HST data of the other HiPASS dwarf galaxies (Grossi et al. 2003).

One possible scenario is that HiPASS J1321-31 experienced an episode of enhanced star formation, and possibly a mild starburst, less than a Gyr ago when the stars now populating the RSGB were formed. Using the Bertelli et al. (1994) and Girardi et al. (2002) isochrones, we found that evolved stars with masses around 2-3 $`M_{}`$ and low metallicities on the red branch of the core-helium burning phase would be dominant, while the corresponding blue loop stars would be too faint to be detected in our observations. The presence of the faint blue hump implies that star formation continued for a few hundred million years at a lower rate. Figure 2a shows a number of theoretical isochrones (Girardi et al. 2002) for $`Z=0.0004`$. While this gives some idea of the age range of the RSGB, more detailed modeling will be discussed in Grossi et al. (2003). In addition to an absence of luminous blue plume stars, we did not find any Hii regions in narrow-band H$`\alpha `$ images taken with the WIYN telescope. HiPASS J1321-31 has a high amount of gas in it ($`M_{HI}/L_B=5.0`$; Table 1). If it did experience a significant burst of star formation in its past, this could have spurred further star formation throughout the galaxy. Yet, there is no evidence of recent or continuous star formation. An interesting question then is not only why the star formation rose in the past but also why it now appears to have fallen?

We find the tip of the RSGB to be around $`m_I=22.6`$ mag yielding an absolute magnitude of $`M_I=6.1`$ mag assuming $`M_I=4.05`$ for the TRGB. The Karachentsev et al. (2002a, b, c; 2003a, b, c) galaxies with well-populated normal RSGBs have tips from $`M_I=7.0`$ up to $`8.7`$. So the faintest RSGB tip in these other galaxies is still about one magnitude brighter than that for HiPASS J1321-31. This suggests that its RSGB derives from an older stellar population since the younger stars which would populate the brighter end of the RSGB are gone, i.e., evolved away, and the corresponding BSGB stars are missing as well. Therefore, HiPASS J1321-31 shows a star formation history that is uncommon in known dwarf galaxies.

Searching through the Karachentsev et al. (2002a, b, c; 2003a, b, c) papers, we have found only two galaxies with similar CMDs to HiPASS J1321-31 whose observed properties we list in Table 1. None of the Local Group dwarfs have this feature in their CMDs. ESO 444-G084 (Karachentsev et al. 2002b), another galaxy in the Cen A group, shows a possible RSGB which is less populated than the one in HiPASS J1321-31, along with a fainter population of blue stars. Another nearby dwarf galaxy KK 65 (Karachentsev et al. 2003a) has a CMD with what looks to be a dispersed RSGB. A faint population of blue stars can also be seen, but the difference between the tip of the RSGB and the blue plume is not as great as that in HiPASS J1321-31. From Table 1 we can see that each of these galaxies are gas-rich, but HiPASS J1321-31 stands apart due to its low luminosity. Another dwarf galaxy, UGC 8833, found in the Canes Venatici I cloud has a similar blue plume and RGB to HiPASS J1321-31, but it lacks a RSGB (Karachentsev et al. 2003b). The physical properties of these two galaxies are similar, except that HiPASS J1321-31 has a much higher Hi mass-to-light ratio (see Table 1). This comparison helps illustrate that HiPASS J1321-31 must have had a period of increased star formation in its rather recent past in order to create its RSGB.

## 5 Summary & Conclusions

We present the CMD for the gas-rich dwarf galaxy HiPASS J1321-31 as observed by the HST/WFPC2. A bright red plume of stars along with a much fainter ($`2`$ mag) population of blue stars indicate that the galaxy has experienced a peculiar star formation history. Three explanations for this red plume were discussed. We rule out the possibilities of it being a RGB or AGB, and conclude that the red plume is composed of core helium-burning post-main sequence stars with ages $``$ 1 Gyr, possibly as young as 0.4 Gyr. These are likely to be related to an epoch when the galaxy experienced an increase in the star formation rate. This event was not associated with the birth of the galaxy since an older RGB is present — HiPASS J1321-31 is not a young galaxy. Several puzzles remain, such as the source of the possible starburst and how the large gas content and unusual evolutionary history might be linked (e.g., could this dwarf consists of tidal debris? see Marakova et al. 2002).

This research was supported in part by NASA through grant number GO-09115.07 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. MG and RFM acknowledge the support of the UK Particle Physics and Astronomy Research Council. Thanks to P. B. Stetson for sharing his PSFs for the WFPC2 and his data reduction programs.