text:
# Frequency analysis of 𝛿 Scuti and RR Lyrae stars in the OGLE-1 database

## 1. Introduction

The Optical Gravitational Lensing Experiment (Udalski et al. 1992) is devoted to search for dark matter in our galaxy using microlensing phenomena. As a by-product of this program an extensive photometry of variable stars in the Galactic Bulge has been accumulated (Udalski et al. 1994, 1995, 1996, 1997). We have performed a systematic frequency analysis of RR Lyrae and $`\delta `$ Scuti stars of this sample. The results for monoperiodic variables are presented elsewhere (Poretti 2001). Here, we discuss the main properties of the identified multiperiodic pulsators.

## 2. Search for multiperiodicity

The OGLE-1 database contains 215 RR Lyrae and 53 $`\delta `$ Scuti variables. The data span $`900`$ days, with typically 130-150 I-band measurements per star. As a first step of analysis, the lightcurve is fitted with the Fourier sum of the form

$$\mathrm{m}_\mathrm{I}(t)=\mathrm{A}_0+\underset{i=1}{\overset{N}{}}\mathrm{A}_i\mathrm{cos}(\frac{2\pi }{\mathrm{P}_1}t+\varphi _i)$$

with primary pulsation period $`\mathrm{P}_1`$ being adjusted as well. Next, the search for secondary periodicities is performed. Two different metods are applied:

a) we compute Fourier power spectrum of the residuals of the fit

b) we supplement the fitting formula with additional cosine term with a trial period $`\mathrm{P}_2`$. We then fit the data for different values of $`\mathrm{P}_2`$, searching for the value that reduces the dispersion of the fit in a significant way.

Both methods yield the same results. As a third step, the Fourier fit with two identified frequencies and their detectable linear combinations is performed. The search for additional periodicities is then repeated. The process is stopped when no significant term appears.

## 3. Canonical double-mode pulsators

We have identified 6 ”normal” double-mode variables, five among $`\delta `$ Scuti stars and one (BW7 V30) among RR Lyrae stars. They are listed in Table 1. The period ratios are very typical and correspond to the ratios of the first two radial overtones. The only exception is BWC V82 with $`\mathrm{P}_2/\mathrm{P}_1=0.55`$, where a higher order radial mode (third overtone) has to be invoked.

## 4. Blazhko effect: nonradial oscillations in RR Lyrae stars

In 37 RR Lyrae variables a different type of multiperiodicity is found – additional peaks very close to the primary pulsation frequency are present. Their amplitudes are usually below 0.06 mag. The secondary frequencies are well-resolved within our dataset and are not due to a secular period variability. Their beating with the primary (radial) pulsation results in an apparent long-term amplitude and phase modulation, a phenomenon referred to as Blazhko effect.

This new multiperiodic behaviour comes in two different flavours: we see either a single secondary peak, forming a doublet with the primary frequency (1 RRc and 27 RRab stars) or a pair of secondary peaks, which together with the primary frequency form an equidistant triplet centered on the primary peak (1 RRc and 8 RRab stars).

While the frequency triplet can result from periodic amplitude and/or phase modulation of a purely radial pulsation, such a process cannot produce a doublet. The observed period ratios ($`\mathrm{P}_2/\mathrm{P}_1=0.951.02`$) are not compatible with excitation of two radial modes. The unavoidable conclusion is that a secondary component of the doublet must correspond to a nonradial mode of oscillation.

In Table 2 we present the statistics of variables with closely spaced frequencies for the Galactic Bulge sample, compared with the numbers derived for the LMC (Alcock et al. 2000; Welch et. al. 2002). In both populations close doublets/triplets occur more frequently in fundamental-mode pulsators than in overtone pulsators. Interestingly enough, the fraction of multiperiodic RRab stars is two times higer in the Galactic Bulge than in the LMC. It is tempting to speculate that the difference in the incidence rate is related to difference in metallicity of the two populations. This hypothesis can be tested when photometry of RR Lyrae stars in the SMC is analysed.

With the sample of 35 multiperiodic RRab stars, we can discuss the group properties of this type of variables. In Fig. 1 we plot amplitudes of the Galactic Bulge RRab stars as a function of their periods. The presence of secondary frequencies has no apparent effect on the pulsation amplitude. We note that close frequency doublets are detected with roughly the same probability at all periods represented in this sample. The occurrence of triplets, on the other hand, seems to be limited to $`\mathrm{P}<0.6`$ day.

Fig. 2 shows the frequency separation $`\mathrm{\Delta }\mathrm{f}=\mathrm{f}_2\mathrm{f}_1`$ ($`\mathrm{f}=1/\mathrm{P}`$) for multiperiodic RRab variables identified in this work. For 80% of all doublets $`\mathrm{\Delta }\mathrm{f}`$ is positive, corresponding to secondary frequency being higher than the primary one. Identical distribution of $`\mathrm{\Delta }\mathrm{f}`$ has also been found for RRab stars in the LMC (Welch et al. 2002). The negative values of $`\mathrm{\Delta }\mathrm{f}`$ occur only in a narrow range of periods between 0.49 day and 0.53 day. Curiously, this particular period range seems to be avoided by stars with triplets.

The frequency separation in the triplets tends to be slightly smaller than in the doubles. For both triplets and doublets the separation $`\mathrm{\Delta }\mathrm{f}`$ is significantly smaller than in the overtone RR Lyrae variables of the LMC (Alcock et al. 2000).

Acknowledgements. This work has been supported by Polish KBN grants 2 P03D 002 20 and 5 P03D 012 20.

## References

Alcock, C., Allsman, R., Alves, D. R., et al. 2000, ApJ, 542, 257

Poretti, E. 2001, A&A, 371, 986

Udalski, A., Szymański, M., Kałużny, J. & Kubiak, M. 1993, Acta Astron., 42, 253

Udalski, A., Kubiak, M., Szymański, M., et al. 1994, Acta Astron., 44, 317

Udalski, A., Olech, A., Szymański, M., et al. 1995, Acta Astron., 45, 433

Udalski, A., Olech, A., Szymański, M., et al. 1996, Acta Astron., 46, 51

Udalski, A., Olech, A., Szymański, M., et al. 1997, Acta Astron., 47, 1

Welch, D., Kovács, G., Cook, K., et al. 2002, these proceedings


format: text

text:
# Probing colour separate singlet states in 𝑒⁺𝑒⁻ annihilation at 𝑍⁰ pole

## 1 Introduction

In high energy strong interaction processes, e.g., $`e^+e^{}hadrons`$, the confinement mechanism demands that at the end of the PQCD cascade the partons should be connected, phenomenologically, as the strings in Lund model, or clusters in Webber model,  from which the hadronization starts. This ‘topology’ of the parton system on the hard and soft interface is beyond PQCD approach.

In the ‘standard models’ of parton cascade and hadronization, , the $`N_C\mathrm{}`$ approximation($`N_C`$ is the number of colour) is applied. This approximation leads to the result that at the end of the PQCD cascade, there is a direct correspondence between the parton states and the string/cluster chain states. Such a simple but phenomenological approximation conceals the complexity of the colour connections on the hard and soft interface. It had not yet met any real challenge mainly because the ‘standard models’,  which employ the large $`N_C`$ approximation can reproduce most of the experimental data especially the global properties until now.

On the other hand, the straightforward ‘group calculus’ (i.e., direct product and reduction of the representations of the $`SU_C(3)`$ group) shows that the colour structure of a multi-parton state is copious and complex as long as $`N_C=3`$(see, e.g., , ). This implies that the colour connections of the multi-parton system on the hard and soft interface can be of many kinds. The traditional one, which is employed in the ‘standard models’ and will be called colour singlet chain states hereafter, is only one kind among them. The PQCD calculations show that, projected onto the colour space of PQCD final state, such kind of colour connection can never have the probability 1 as long as $`N_C=3`$ and that when parton number becomes larger the probability is reduced. Another kind of colour connections is colour separate singlet states, with the distinguish feature that several gluons form colour singlet clusters and can hadronize independently. Since only 3 kinds of colour(8 kinds of ‘double colours’), when parton number is large111as calculated by JETSET, the average numbers of gluons are: $`\sqrt{s}=91GeV,<N>6;\sqrt{s}=200GeV,<N>9;\sqrt{s}=1TeV,<N>17`$, it is inevitably that more and more gluons have identical colours, and the colour separate states could have large probability.

We should notice that, from our previous works, the above two kinds of colour configurations, colour singlet chains with various orders for the gluons in the chain, colour separate singlets with different gluons combined into singlet clusters, are not orthogonal to each other as long as $`N_C=3`$ and within PQCD framework. They belong to different complete sets of bases for the PQCD final state colour space. Hence PQCD sees no differences among them. On the other hand, different colour connections mean different string/cluster states. This implies that different colour connections can lead to differences in handron states. Furthermore, the cross section of any colour connection we calculated in PQCD by projecting it onto the PQCD colour space could have no realistic meanings. Their relative probabilities should depend on the unknown NPQCD mechanism. To analysize these issues, one needs:

1) To confirm if the colour singlet chain states and the colour separate singlet states really do not equal to each other by seeing if they lead to the differences in hadron states for the same parton state and hadronization scheme.

2)If they do not equal, to see how does nature chooses among them (maybe) by the(unknown) NPQCD laws.

Because of the lack of the knowledge of NPQCD, we can not answer these questions from the QCD Lagrangian, but employing phenomenological models. The straightforward way is to modify the event generators(e.g., JETSET as discussed following): Only the colour connections changed from colour singlet chains to colour separate singlets while the other parts (parton cascade and hadronization scheme, etc) kept. Then one can use the modified generators to produce the hadron states to see if there are some deviations from the corresponding traditional results. Further step is to compare with the data especially those the unmodified generators can not explain, to extract the relative weights of all the connections.

There have been some works along this route(e.g.,  and refs. therein). Recently we constructed a colour separate singlet modelby modifying the standard JETSET7.4 code. Hereafter, we use CS(Colour Separate) for the modified JETSET, and CC(Colour singlet Chain) for the standard. However, the results of the CS seems an ambiguity: Though the colour connection is different from that of CC, the observables describing the hadron states averaged among unbiased events have no significant differences from those of th CC provided that a T-measure to weight the colour separate states (Hence ‘short strings’ are favoured(see, e.g.,  and refs. therein). Considering the uncertainties in the model, somehow we can say that for unbiased events, these two kinds of colour connections can be thought of as leading to the same hadron states.

The clues to understand this ambiguity are: The variables describing ‘global properties’ are mainly decided by PQCD, not sensitive to the NPQCD properties, especially when averaged for unbiased events, as indicated by the hyperthesis local parton-hadron duality. Hence jet shape variables, inclusive momentum distributions ’should’ be similar for CS and CC, when averaged for UNBIASED events. This indicate that we must look into special events. The baryon anti-baryon rapidity correlations which are expected to be more sensitive to the NPQCD properties however, are also not able to expose the differences of the two colour connections. So other more sensitive observables need to be looked for. Another reason lies on the colour separate model we constructed . Only the singlet ‘glueball’ cluster and its neighbouring 2 strings of the CS are different from those of CC, so the difference is somehow ‘localized’ in a region of the phase space. This tells that the sensitive observables are more possibly defined in a ’window’ of the phase space. The above arguments will be much more clear after considering the following sections.

In this paper we selected two-jet like events, which will be described following. First, taking the advantage of generators, we select two-jet like events in the parton level, i.e., by requiring the energy of all the gluons(including medium state gluons) smaller than a certain value. This part of events shows that CS and CC lead to hadron states with significant differences when measure a certain observable(section 2). Second, in section 3, we select 2-jet like events according to the transverse momenta of the hadrons, which can be done in experiment. This part of events shows that the colour separate singlet state can be ’probed’: For these events, the CS also lead to significant differences from the CC, and comparing with data we can extract the relative rate between the two kinds of colour connections . All the numerical analysize are performed with the default JETSET7.4 (for CC) and the modified generator(for CS). The numerical results are for the $`Z^0`$ pole and the selected two-jet like events are of the rate $``$ one per mill.

## 2 Two-jet like events at parton level and a sensitive variable

The parton level selected events and its properties are described in . Here we define

$$R=\underset{i}{}|p_{Ti}|,|y_i|<y_0$$
(1)
. in which the $`p_{Ti}`$ and $`y_i`$ are the transverse momentum and rapidity of the hadrons concerned. This observable can ’accumulate’ the various differences between CS from CC. Its distribution for CS and CC is shown in Fig. 1

The large differences confirm us that, for the events we select, the CS and CC colour connections really lead to large differences in hadron states. Hence we can say that they do not equal to each other from the point of view of hadronization/NPQCD. The differences can not be smeared by tuning various parameters in the generator.

## 3 Selecting two-jet like events at hadron level

The discussions of the above section can be considered as a ideal experiment, which confirm us that the different colour connections really leads to different hadron states. Now let’s see if we can select such kind of events from experiments. By comparing the data and the predictions of different colour connections, we can extract the relative weight of them.

It is easy to understand that no direct correspondence to the parton level selected two-jet like events from the point of view of experiment. The crucial reason is that the set of events with the energies of all the gluons in a event smaller than a certain value have no definite meaning. It is affected by, e.g., the cut-off of the parton cascade. At the same time, It is difficult to select these events with high purity by the jet shape variables. However, we can see from Fig.2 that the $`p_T`$ distribution of final state hadrons for these events is significant different from that of the unbiased events. Especially in the region $`p_T<0.5GeV/c`$, these events can take a large part. Hence from now on we define our two-jet like events as those events whose final state hadrons all have the transverse momenta smaller than $`0.5GeV/c`$. This can be measured in experiments and have definite physical meaning.

The most important thing is that, the newly-defined events have very the same properties as the parton level selected events.

Now we find that for the new-defined two-jet like events, there is large differences for the R distribution. We can select these events, and compare the data and theoretical prediction for variouis observables, the rate of different colour connections can be extracted.

## 4 Discussions

The conclusion of last section is very clear. One thing to emphasize is that, The two-jet like events we defined in last section have never been selected and compare with generator results until now. So when taking the comparison to extract the rate of different colour connections, one should worry about the validity of the generators. So what we have proposed here is also a very interesting thing to do to test the generator.

The two-jet like events and the observable R looks similar as those in  at the first sight, But there are crucial differences. In that paper, for investigating hadronization, the less gluons the better, however we would like events with many soft gluons so that CS will have large probability. In that paper, they take a sum of the vector $`p_T`$ to demonstrate different correlations among different models, while we take the sum of the absolute value of the $`p_T`$ to accumulate the differences.


text:
# Modern Shell-model Calculations

## 1 Introduction

The shell model is the basic framework for nuclear structure calculations in terms of nucleons. To make it a truly microscopic approach demands, however, that the model-space effective interaction $`V_{\mathrm{eff}}`$ be derived from the nucleon-nucleon ($`NN`$) interaction. This has long been indeed a topical problem in nuclear structure theory, the first steps in this direction dating back to the mid 1960s. The two basic ingredients which come into play in what have become known as “realistic” shell-model calculations are the free $`NN`$ potential and the many-body methods for deriving the effective interaction in the chosen model space. ¿From the early 1980s on there has been much progress in this field. On the one hand, high-quality $`NN`$ potentials have been constructed which reproduce quite accurately the $`NN`$ scattering data and are suitable for application in nuclear structure. On the other hand, an accurate calculation of the Brueckner $`G`$ matrix has become feasible while the so-called folded-diagram series for the effective interaction $`V_{\mathrm{eff}}`$ can be summed up to all orders using iterative methods. These improvements have brought about a revival of interest in realistic shell-model calculations.

In this context, the main question is how accurate a description of nuclear structure properties can be provided by effective interactions derived from the free $`NN`$ potential. In recent years, we have concentrated our efforts on this problem studying a number of nuclei around doubly magic 100Sn, 132Sn, and 208Pb.2-9 Our attention has been focused on nuclei with few valence particles or holes since they provide the best testing ground for the basic ingredients of shell-model calculations, especially as regards the matrix elements of the two-body effective interaction. The results of our calculations have turned out to be in remarkably good agreement with experiment for all the nuclei considered (about 20), providing evidence that realistic effective interactions are able to describe with quantitative accuracy the spectroscopic properties of complex nuclei. A brief survey of the present status of realistic shell-model calculations based on this comprehensive study has been given in Ref. 10.

While we have also performed calculations starting from different $`NN`$ potentials, all the results summarized in Ref. 10 have been obtained by making use of effective interactions derived from the Bonn-A free $`NN`$ potential by means of a $`G`$-matrix formalism including renormalizations from both core polarizations and folded diagrams.

The use of the $`G`$ matrix has long been proved to be a valuable tool to overcome the difficulty posed by the strong repulsive core contained in all modern $`NN`$ potentials. As mentioned above, an accurate calculation of it can be carried out which is based on the Tsai-Kuo method. It should be recalled, however, that the $`G`$ matrix is model-space dependent, which makes its actual calculation rather involved. Very recently, a completely new approach has been proposed which consists in deriving from the $`NN`$ potential a low-momentum effective potential $`V_{lowk}`$, defined within a given cut-off momentum $`\mathrm{\Lambda }`$. This is a smooth potential which may be used directly in nuclear structure calculations.

The main aim of this paper is to report on realistic shell-model calculations performed by using as initial input the CD-Bonn potential and deriving the effective interaction by means of both the $`G`$-matrix and the $`V_{lowk}`$ formalisms. Our calculations concern the 132Sn neighbors with two, three, and four valence protons, namely 134Te, 135I, and 136Xe.

Our presentation is organized as follows. In Sect. 2 we give an outline of our calculations. In Sect. 3 we present the results obtained by using the $`V_{lowk}`$ as well as the $`G`$-matrix and compare them with the experimental data. Some concluding remarks are given in Sect. 4.

## 2 Outline of calculations

As already mentioned in the Introduction, we present in this paper the results of calculations performed by using a new approach to the derivation of $`V_{\mathrm{eff}}`$ as well as those obtained by the usual $`G`$-matrix formalism.

A description of our procedure for deriving $`V_{\mathrm{eff}}`$ within the framework of a $`G`$-matrix folded-diagram formalism can be found in Refs. 15 and 16. We only outline here the essentials of the method to point out the main differences between this approach and that based on the use of $`V_{lowk}`$.

The $`G`$ matrix is defined by the integral equation

$$G(\omega )=V+VQ_2\frac{1}{\omega Q_2TQ_2}Q_2G(\omega ),$$
(1)
where $`V`$ represents the $`NN`$ potential, $`T`$ denotes the two-nucleon kinetic energy, and $`\omega `$ is an energy variable (the so-called starting energy). The two-body Pauli exclusion operator $`Q_2`$ prevents double counting, namely the intermediate states allowed for $`G`$ must be outside of the chosen model space. Thus the Pauli operator $`Q_2`$ is dependent on the model space, and so is the $`G`$ matrix. The operator $`Q_2`$ is specified, as discussed in Ref. 17, by a set of three numbers ($`n_1,n_2,n_3`$), each representing a shell-model orbital. We employ a matrix inversion method to calculate the above $`G`$ matrix. This method gives the solution of Eq.(1) as

$$G=G_F+\mathrm{\Delta }G.$$
(2)
The “free” $`G_F`$ matrix is

$$G_F=V+V\frac{1}{e}G_F,$$
(3)
and the Pauli correction term $`\mathrm{\Delta }G`$ is given by

$$\mathrm{\Delta }G=G_F\frac{1}{e}P_2\frac{1}{P_2(\frac{1}{e}+\frac{1}{e}G_F\frac{1}{e})P_2}P_2\frac{1}{e}G_F,$$
(4)
with $`e(\omega T)`$ and $`P_2=1Q_2`$. The calculation of $`G_F`$ is straightforward, as it does not contain the Pauli projection operator. Then $`\mathrm{\Delta }G`$ can be calculated by performing some matrix operations in the model space $`P_2`$.

Using the above $`G`$ matrix we can calculate $`V_{\mathrm{eff}}`$. This can be written schematically in operator form as

$$V_{eff}=\widehat{Q}\widehat{Q^{}}\widehat{Q}+\widehat{Q^{}}\widehat{Q}\widehat{Q}\widehat{Q^{}}\widehat{Q}\widehat{Q}\widehat{Q}+\mathrm{},$$
(5)
where $`\widehat{Q}`$ (referred to as $`\widehat{Q}`$-box) is a vertex function composed of irreducible linked diagrams, and the integral sign represents a generalized folding operation. $`\widehat{Q^{}}`$ is obtained from $`\widehat{Q}`$ by removing terms first order in $`V`$. In our calculations we take the $`\widehat{Q}`$-box to be composed of $`G`$-matrix diagrams through second order. After the $`\widehat{Q}`$-box is calculated, $`V_{\mathrm{eff}}`$ is obtained by summing up the folded-diagram series of Eq.(5) to all orders by means of the Lee-Suzuki iteration method. This last step can be performed in an essentially exact way for a given $`\widehat{Q}`$-box.

A detailed description of the derivation of $`V_{lowk}`$ as well as a discussion of its main features can be found in Refs. 13 and 20. The basic idea of this approach is to reduce a realistic $`NN`$ potential $`V_{NN}`$ to an effective low-momentum potential $`V_{lowk}`$ by integrating out the high-momentum modes. The $`V_{lowk}`$ constructed in this way is confined within a cut-off momentum $`\mathrm{\Lambda }`$, and reproduces exactly the deuteron binding energy and the low-energy phase shifts (up to $`E_{Lab}=2\mathrm{}^2\mathrm{\Lambda }^2/M`$) of $`V_{NN}`$. It has the desirable feature of being a smooth potential and as such it is suitable to be used directly in shell-model calculations. The $`V_{\mathrm{eff}}`$ is derived by following the same procedure as that outlined above, except that the $`G`$ matrix is replaced by $`V_{lowk}`$. More explicitly, we first calculate the $`\widehat{Q}`$-box including diagrams up to second order in $`V_{lowk}`$ and then obtain $`V_{\mathrm{eff}}`$ by summing up the folded-diagram series using the Lee-Suzuki iteration method.

In our study of 134Te, 135I, and 136Xe we assume that 132Sn is a closed core and let the valence protons occupy the five single-particle (s.p.) orbits $`0g_{7/2}`$, $`1d_{5/2}`$, $`1d_{3/2}`$, $`2s_{1/2}`$, and $`h_{11/2}`$. As regards the energy spacings between the five s.p. levels, we take three of them from the experimental spectrum of 133Sb. As for the $`s_{1/2}`$ state, we adopt the value $`ϵ_{1/2}`$=2.8 MeV, which reproduces the experimental energy of the $`\frac{1}{2}^+`$ level at 2.15 MeV in 137Cs. For the shell-model oscillator parameter $`\mathrm{}\omega `$ we use the value 7.88 MeV, as obtained from the expression $`\mathrm{}\omega =45A^{1/3}25A^{2/3}`$ for $`A=132`$.

In the $`G`$-matrix calculations we have fixed ($`n_1,n_2,n_3`$)=(16,28,55) for the neutron orbits, and ($`n_1,n_2,n_3`$)=(11,21,55) for the proton orbits. As regards the calculations using $`V_{lowk}`$ as the input interaction, an important question is how to choose the momentum cut-off $`\mathrm{\Lambda }`$. According to the considerations made in Refs. 13 and 20, we have used in our calculation the value $`\mathrm{\Lambda }=2.05`$ fm-1. We have verified that the results are practically insensitive to changes in $`\mathrm{\Lambda }`$ within reasonable limits, namely 1.9-2.2 fm-1.

In both the $`G`$-matrix and the $`V_{lowk}`$ approaches the effective interaction has been derived from the CD-Bonn free $`NN`$ potential which fits the world $`NN`$ data below 350 MeV with a $`\chi ^2`$/datum of 1.02. Once the effective interaction has been derived, the shell-model calculations are carried out employing the OXBASH code.

## 3 Results

We present here the results of our realistic shell-model calculations for the $`N=82`$ isotones 134Te, 135I, and 136Xe. As already pointed out in the preceding sections, these results have all been obtained starting from the CD-Bonn $`NN`$ potential. In Figs. 1-3 we report the experimental23-25 and calculated spectra of these three nuclei. The latter have been obtained by using as input interaction the usual Brueckner $`G`$-matrix and the low-momentum $`NN`$ potential $`V_{lowk}`$.

The spectra of 134Te and 135I obtained by using the Bonn-A $`NN`$ potential and the $`G`$-matrix formalism have been presented in a previous paper, where a detailed comparison between theory and experiment is made. While we refer

the reader for details to the above paper, we emphasize here the very good agreement with experiment of both the calculated spectra. A measure of the overall agreement between the calculated and experimental energy levels is provided by the $`rms`$ deviation $`\sigma `$. For the spectra obtained using $`V_{lowk}`$ the values of $`\sigma `$ are 110, 73, and 111 keV for 134Te, 135I, and 136Xe, respectively, while for the $`G`$-matrix calculations they are 163, 101, and 154 keV. This comparison shows that the $`V_{lowk}`$ results are rather similar to, and actually slightly better than, the $`G`$-matrix ones.

## 4 Concluding remarks

In this paper, we have presented some results of realistic shell-model calculations which have been performed within the framework of a new approach to the derivation of the effective interaction. This is based on the use of a low-momentum $`NN`$ potential $`V_{lowk}`$ obtained by integrating out the high-momentum part of the original $`V_{NN}`$ potential containing a strong repulsive core. This $`V_{lowk}`$ is a smooth potential and can be used directly in nuclear structure calculations. To assess the practical value of this new approach, we have also performed traditional $`G`$-matrix calculations and compared the $`V_{lowk}`$ and $`G`$-matrix results between them and with experimental data. It turns out that not only the former are in very good agreement with experiment, they are even somewhat better than the latter.

On the above grounds we may conclude that the $`V_{lowk}`$ approach provides an advantageous alternative to the usual $`G`$-matrix formalism in realistic shell-model calculations. We are currently using this approach for further nuclear structure studies. It is worth mentioning, however, that it is being profitably used also in other fields, as for instance the Fermi liquid theory for nuclear matter.

## Acknowledgments

This work was supported in part by the Italian Ministero dell’Università e della Ricerca Scientifica e Tecnologica (MURST) and by the U.S. DOE Grant No. DE-FG02-88ER40388.


text:
# Concentrations and Assembly Histories of Dark Matter Halos

## 1 METHOD

We investigate the connection between halo density profiles and their mass assembly histories, using a structural merger tree constructed from a high-resolution N-body simulation of a flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_m=0.3`$, $`h=0.7`$ and $`\sigma _8`$ = 1.0, whose evolution has been simulated with the ART code . The trajectories of $`256^3`$ cold dark matter particles are followed within a cubic, periodic box of comoving size 60$`h^1`$Mpc from redshift $`z=40`$ to the present. We use distinct halo catalogs at 36 output times spaced between $`z=7`$ and $`0`$. NFW density profiles , $`\rho _{\mathrm{NFW}}(r)=\rho _\mathrm{s}/[(r/R_\mathrm{s})\left(1+r/R_\mathrm{s}\right)^2]`$, are measured for each halo with more than 200 particles, corresponding to halos more massive than $`2.2\times 10^{11}h^1\mathrm{M}_{}`$. For each halo at $`z=0`$, we identify its full merger history, and determine concentration parameters $`c_{\mathrm{vir}}R_{\mathrm{vir}}/R_\mathrm{s}`$ for all progenitors, thus providing a structural merger tree for each of $`3000`$ halos.

## 2 RESULTS

Figure 1 (left) shows the history of mass growth for the major progenitors of several different halos, spanning a range of masses and concentration parameters. Massive halos tend show substantial mass accumulation up to late times, while the growth curves for less massive halos tend to flatten out earlier. By examining a range of full mass assembly histories for our sample of halos, we find that both average mass accretion histories and mass accretion histories for individual halos are well characterized by a simple function: $`M(a)=M_\mathrm{o}e^{\alpha z},a=(1+z)^1`$. Fits to this equation are shown in Figure 1 (left) for representative individual halos. The single free parameter $`\alpha `$ can be related to a characteristic epoch for formation, $`a_\mathrm{c}`$, defined as the expansion scale factor $`a`$ when the logarithmic slope of the accretion rate, $`\mathrm{d}\mathrm{log}M/\mathrm{d}\mathrm{log}a`$, falls below some specified value, $`S`$; the functional form implies $`a_\mathrm{c}=\alpha /S`$. The same formation epoch can be defined equivalently for any “observing” epoch $`z_\mathrm{o}`$ of that halo, by replacing $`a`$ by $`a/a_\mathrm{o}`$, in which case $`a_\mathrm{c}=a_\mathrm{o}\alpha /S`$. Thus at any such observing redshift, with scalefactor $`a_\mathrm{o}=1/(1+z_\mathrm{o})`$ and mass $`M_\mathrm{o}=M(z_\mathrm{o})`$, the mass growth is fit by

$$M(a)=M_\mathrm{o}\mathrm{exp}\left[a_\mathrm{c}S\left(\frac{a_\mathrm{o}}{a}1\right)\right].$$
(1)
This implies that, for any halo whose mass accretion trajectory follows this functional form, the characteristic formation time is the same regardless of the redshift $`z_\mathrm{o}`$ at which the halo is observed.

We find that the concentration of a halo, defined as $`c_{\mathrm{vir}}R_{\mathrm{vir}}/R_\mathrm{s}`$, is tightly correlated with the characteristic formation epoch as defined above, and that this relation holds at all redshifts when, properly scaled by $`a_\mathrm{o}`$:

$$c_{\mathrm{vir}}=c_1a_\mathrm{o}/a_\mathrm{c},$$
(2)
where $`c_1`$ is the typical concentration of halos whose formation time is at the time of measurement, $`a_c=a_\mathrm{o}`$. Figure 1 (right) shows that this formula provides a good description of the observed correlation between concentration and formation time for halos at all masses and redshifts. As mentioned, the typical formation time is a function of mass, but there is significant scatter in $`a_\mathrm{c}`$ for a fixed mass. The relation defined by Eq. 2 is able to account for the complete mass and redshift dependence of $`c_{\mathrm{vir}}`$, and for the scatter in $`c_{\mathrm{vir}}`$ measured for fixed mass halos. This correlation is has important consequences for galaxy formation modeling. For further details on this work, see .


text:
# Semiclassical analysis of a two-electron quantum dot in a magnetic field: dimensional phenomena

## Abstract

While the dynamics for three-dimensional axially symmetric two-electron quantum dots with parabolic confinement potentials is in general non-separable we have found an exact separability with three quantum numbers for specific values of the magnetic field. Furthermore, it is shown that the magnetic properties such as the magnetic moment and the susceptibility are sensitive to the presence and strength of a vertical confinement. Using a semiclassical approach the calculation of the eigenvalues reduces to simple quadratures providing a transparent and almost analytical quantization of the quantum dot energy levels which differ from the exact energies only by a few percent.

PACS numbers: 73.21.La, 03.65.Sq, 75.75.+a, 05.45.Mt

Current nanofabrication technology allows one to control the size and shape of quantum dots . Due to the confinement of the electrons in all three spatial directions the energy spectrum is quantized creating excellent experimental and theoretical opportunities to study controlled single-particle and collective dynamics at the atomic scale. For example, depending on the experimental setup, the spectrum of a quantum dot displays shell structure  or follows predictions of random matrix theory (for a review see ). Furthermore, it becomes possible to trace the transition from a quantum mechanical to an almost classical regime.

Few-electron quantum dots have attracted special attention , since they may provide a natural realization of a quantum bit . The simplest quantum dot with the essential features of more complex systems contains two electrons. Experimental data, including transport measurements  and spin oscillations in the ground state under a perpendicular magnetic field , have been explained quantum mechanically as a result of the interplay between the two-dimensional lateral confinement potential, electron correlations and the magnetic field . While in these experiments effects of the third spatial dimension are somewhat hidden, they naturally come into play with a tilted magnetic field . However, even with a perpendicular magnetic field the vertical confinement has at least two important consequences which will be worked out in the following: first, it changes the magnetic moment and susceptibility with respect to the 2D results, second, the generically non-separable 3D dynamics becomes separable for certain values of the magnetic field. We will use a semiclassical description which offers a simple and accurate approach to explore the effects of dimensionality in quantum dots. In contrast to a circular (2D) two-electron quantum dot whose classical dynamics is always separable and therefore regular, the corresponding 3D-system with axial symmetry is in general a non-integrable problem with typical features of mixed dynamics (regular/chaotic).

The Hamiltonian for the 3D two-electron quantum dot reads

$$H=\underset{j=1}{\overset{2}{}}\left\{\frac{1}{2m^{}}(𝐩_j\frac{e}{c}𝐀_j)^2+\frac{m^{}}{2}[\omega _0^2(x_j^2+y_j^2)+\omega _z^2z_j^2]\right\}+V_C+H_{\mathrm{spin}},$$
(1)
where $`V_C=\alpha /|𝐫_1𝐫_2|`$ is the Coulomb energy ($`\alpha =e^2/(4\pi \epsilon \epsilon _0)`$) and $`H_{\mathrm{spin}}=g^{}(𝐬_1+𝐬_2)𝐁`$ describes the Zeeman energy. Here $`m^{}`$ and $`g^{}`$ are the effective electron mass and $`g`$-factor, respectively, and $`\epsilon `$ is the dielectric constant. The confining potential is approximated with a 3D axially-symmetric harmonic oscillator and $`\mathrm{}\omega _z\mathrm{}\omega _0`$ are the energy scales of confinement in the $`z`$-direction and in the $`xy`$-plane, respectively. For the typical voltage $`1`$ V applied to the gate, the confining potential is some eV deep which is large compared to the few meV of the confining frequency . Hence, the electron wave function is localized close to the minimum of the well which always can be approximated by a parabolic potential. In real samples the electron-electron interaction is usually screened. However, the pure Coulomb interaction should suffice to understand the main features of the system. For the perpendicular magnetic field $`(𝐁z)`$ we choose a gauge described by the vector $`𝐀=[𝐁\times 𝐫]/2=\frac{1}{2}B(y,x,0)`$. Introducing the relative and center-of-mass coordinates $`𝐫=𝐫_1𝐫_2`$, $`𝐑=\frac{1}{2}(𝐫_1+𝐫_2)`$, the Hamiltonian, Eq.(1), can be separated into the center-of-mass (CM) $`H_{\mathrm{CM}}`$ and relative-motion (RM) $`H_{\mathrm{rel}}`$ terms: $`H=H_{\mathrm{CM}}+H_{\mathrm{rel}}+H_{\mathrm{spin}}`$. The solution to the CM-Hamiltonian is well known  and the effect of the Zeeman energy has been discussed in . In the following we will concentrate on the dynamics of $`H_{\mathrm{rel}}`$.

For our analysis it is convenient to use cylindrical scaled coordinates, $`\stackrel{~}{\rho }=\rho /l_0`$, $`\stackrel{~}{p}_\rho =p_\rho l_0/\mathrm{}`$, $`\stackrel{~}{z}=z/l_0`$, $`\stackrel{~}{p}_z=p_zl_0/\mathrm{}`$, where $`l_0=(\mathrm{}/\mu \omega _0)^{1/2}`$ is the characteristic length of the confinement potential with the reduced mass $`\mu =m^{}/2`$. The strength parameter $`\alpha `$ of the Coulomb repulsion goes over to $`\lambda =2\alpha /(\mathrm{}\omega _0l_0)`$. Using the effective mass $`m^{}=0.067m_e`$, the dielectric constant $`\epsilon =12`$, which are typical for GaAs, and the confining frequency $`\mathrm{}\omega _0`$=3 meV, we obtain $`\lambda 3`$. Hereafter, for the sake of simplicity, we drop the tilde, i.e. for the scaled variables we use the same symbols as before scaling.

In these variables the Hamiltonian for the relative motion takes a particular simple form (in units of $`\mathrm{}\omega _0`$)

$$ϵ\frac{H_{\mathrm{rel}}}{\mathrm{}\omega _0}=\frac{1}{2}\left[p_\rho ^2+\frac{m^2}{\rho ^2}+p_z^2+\left(\frac{\omega _\rho }{\omega _0}\right)^2\rho ^2+\left(\frac{\omega _z}{\omega _0}\right)^2z^2+\frac{\lambda }{\sqrt{\rho ^2+z^2}}\right]\frac{\omega _L}{\omega _0}m,$$
(2)
where $`m=l_z/\mathrm{}`$, $`\omega _L=eB/2m^{}c`$ is the Larmor frequency and

$$\omega _\rho =(\omega _L^2+\omega _0^2)^{1/2}$$
(3)
is the effective confinement frequency in the $`\rho `$-coordinate which depends through $`\omega _L`$ on the magnetic field.

Due to the cylindrical symmetry, the $`z`$-component $`l_zp_\varphi `$ of the angular momentum is conserved and the motion in $`\varphi `$ is separated from the motion in the $`(\rho ,z)`$-plane. Since the Coulomb term couples the two coordinates, the problem is in general non-integrable which is reflected in the Poincaré sections shown in Fig.1 for increasing magnetic field. The chosen ratio $`\omega _z/\omega _0=3`$ is of the same order of magnitude as in the experiment . For $`\omega _L=0`$ and small values of $`m`$ the motion is mainly chaotic (see Fig.1a). With the magnetic field the frequency of oscillations along the $`\rho `$-coordinate can be controlled which leads to qualitatively different dynamical situations (Fig.1b-d). For equal effective confinement frequencies $`\omega _\rho ^2=\omega _z^2`$, the Hamiltonian Eq.(2) becomes separable in spherical coordinates and the dynamics is integrable (Fig.1c). For two other limiting cases, the dynamics is nearly integrable, namely in the limit $`m\mathrm{}`$ and for $`\omega _z\mathrm{}`$. The latter case represents a two-dimensional quantum dot, classically, we have $`p_z,z0`$ in this limit.

The semiclassical quantization of the circular 2D quantum dot is particularly simple since it reduces to a 1D WKB quantization of the $`\rho `$-motion due to the separability of the problem. For given $`m`$ and $`p_z=z=0`$ the momentum $`p_\rho `$ determined from Eq.(2) enters the action integral

$$I_\rho =\frac{\mathrm{}}{2\pi }p_\rho 𝑑\rho =\frac{\mathrm{}}{\pi }_{\rho _{\mathrm{min}}}^{\rho _{\mathrm{max}}}|p_\rho |𝑑\rho ,$$
(4)
with the turning points $`\rho _{\mathrm{min}}`$, $`\rho _{\mathrm{max}}`$ as the positive roots of equation $`p_\rho (\rho )=0`$. The WKB quantization conditions

$$I_\rho (ϵ)=\mathrm{}(n_\rho +\frac{1}{2}),n_\rho =0,1,\mathrm{},m=0,\pm 1,\mathrm{}$$
(5)
determine the energy levels. For non-interacting electrons ($`\lambda =0`$) the analytical calculation of the action integral leads to the (quantum mechanically exact) eigen-energies

$$ϵ=\sqrt{1+\left(\frac{\omega _L}{\omega _0}\right)^2}(2n_\rho +|m|+1)\frac{\omega _L}{\omega _0}m,$$
(6)
which are the well known Fock-Darwin energies . For $`\lambda 0`$, we calculate the action integral Eq.(4) numerically with a few iterations to determine the quantum eigenvalues.

The energy spectra for non-interacting and interacting electrons are shown in Fig.2. In the interacting case the semiclassical result, although not exact (the error is less than $`1\%`$), reproduces very well the quantum mechanical results .

Turning now to the 3D quantum dot we have seen that the dynamics is separable for $`\omega _z^2=\omega _\rho ^2\omega _{L}^{}{}_{}{}^{2}+\omega _0^2`$ and the Hamiltonian Eq. (2) in scaled spherical coordinates takes the form

$$ϵ=\frac{1}{2}\left\{p_r^2+\left(\frac{\omega _z}{\omega _0}\right)^2r^2+\frac{\lambda }{r}+\frac{(𝐥/\mathrm{})^2}{r^2}\right\}\frac{\omega _L^{}}{\omega _0}m.$$
(7)
In this case the square of the total angular momentum $`𝐥^2`$ is an additional integral of motion. Therefore, the classical dynamics reduces again to a one-dimensional, radial problem. Using Eq.(7) and calculating the action integral for the radial motion analogous to that in Eq.(4) (i.e. with $`r`$ instead of $`\rho `$), we obtain the energy levels from the standard WKB quantization conditions

$`I_r(ϵ)=\mathrm{}(n_r+\frac{1}{2}),|𝐥|=\mathrm{}(l+\frac{1}{2}),`$ (8)
$`n_r,l=0,1,\mathrm{},m=0,\pm 1,\mathrm{},\pm l.`$ (9)
Note that it is only the magnetic field which generates the spherical symmetry of the problem and therefore its separability leading to three good quantum numbers $`n_r`$, $`l`$ and $`m`$.

In the general case of an axially symmetric 3D quantum dot we have non-integrable motion and a semiclassical quantization is neither straight forward nor does it give results which allow for a simple understanding of the dynamics. Therefore, we make use of the fact that in real samples the confining potential in the $`z`$-direction is much stronger than in the $`xy`$-plane which allows us to analyze the 3D non-integrable system with the ’removal of resonances’ method (RRM) . To lowest order the RRM consists of averaging the Hamilton function over the fastest angle of the unperturbed motion $`(\lambda =0)`$ after rewriting coordinates and momenta in terms of action-angle variables $`(J_\rho ,J_z,\theta _\rho ,\theta _z)`$:

$`\rho ^2={\displaystyle \frac{\omega _0}{\omega _\rho }}\left(2j_\rho +|m|2\sqrt{j_\rho (j_\rho +|m|)}\mathrm{cos}2\theta _\rho \right),`$ (10)
$`z^2={\displaystyle \frac{2j_z\omega _0}{\omega _z}}\mathrm{sin}^2\theta _z,`$ (11)
and $`p_\rho =\dot{\rho }`$, $`p_z=\dot{z}`$. Here, $`j_z=J_z/\mathrm{}`$ and $`j_\rho =J_\rho /\mathrm{}`$. If $`\omega _z>\omega _\rho `$ one averages over the angle $`\theta _z=\omega _zt`$. As a result, the motion effectively decouples into an unperturbed motion in the $`z`$-coordinate governed by the potential $`(\omega _z/\omega _0)^2z^2/2`$ and into the relative motion in the $`\rho `$-coordinate governed by the effective potential

$$V_{\mathrm{eff}}(\rho ,j_z)=\frac{1}{2}\left(\frac{\omega _\rho }{\omega _0}\right)^2\rho ^2+\frac{m^2}{2\rho ^2}+\frac{\lambda }{\pi \rho }K\left(2\frac{\omega _0}{\omega _z}\frac{j_z}{\rho ^2}\right),$$
(12)
where $`K`$ is the first elliptic integral. Hence, the effective Hamiltonian reads

$$ϵ=\frac{p_\rho ^2}{2}+V_{\mathrm{eff}}\frac{\omega _L}{\omega _0}m+\frac{\omega _z}{\omega _0}j_z.$$
(13)
Applying a similar procedure as in the 2D case, we calculate the action integral numerically. The momentum $`p_\rho `$ is determined from Eq.(13) and the turning points $`\rho _{\mathrm{min}}`$, $`\rho _{\mathrm{max}}`$ are as usual the (positive) roots of the equation $`p_\rho (\rho )=0`$. Finally, the WKB-quantization conditions

$`I_\rho (ϵ)=\mathrm{}(n_\rho +\frac{1}{2}),j_z=n_z+\frac{1}{2},`$ (14)
$`n_\rho ,n_z=0,1,2,\mathrm{},m=0,\pm 1,\pm 2,\mathrm{},`$ (15)
determine the energy levels.

Comparing the exact results for eigen-energies for the spherical case $`\omega _z/\omega _\rho =1`$ we found good agreement even for large values of the magnetic field (Fig.3a) although RRM is expected to work best for $`\omega _\rho /\omega _z<1`$. Without magnetic field we have $`\omega _\rho /\omega _z=1/3`$ which means that the motion in $`z`$ and $`\rho `$ approximately decouples justifying the widely used 2D approximation. This is also reflected in the small difference between 2D and 3D results (compare Fig.2b with Fig.3a at $`\omega _L=0`$). Turning on the magnetic field increases the coupling of the dynamics in $`\rho `$ and $`z`$ which allows the two electrons eventually to access the full 3D space. As a consequence, the electrons can avoid each other more effectively and the Coulomb interaction has a smaller effect on the 3D spectrum than on the 2D spectrum which is most clearly visible for the $`m=0`$ energies, see Fig.3a. We can understand this effect quantitatively by averaging the elliptic integral in Eq.(12) over the unperturbed ($`\lambda =0`$) motion in $`\rho `$. It gives rise to an effective charge in the Coulomb interaction $`V_C\lambda _{\mathrm{eff}}/2\rho `$, where

$$\lambda _{\mathrm{eff}}=\frac{2\lambda }{\pi ^2}_0^\pi K(\frac{\omega _\rho /\omega _z}{1+|m|\sqrt{1+2|m|}\mathrm{cos}2\theta _\rho })𝑑\theta _\rho $$
(16)
for $`n_\rho =n_z=0(j_\rho =j_z=1/2)`$. The 3D energy quantized with this effective charge for the repulsion is close to the full interaction (dotted line in Fig.3a).

The effective charge $`\lambda _{\mathrm{eff}}/\lambda `$ as a function of $`\omega _\rho /\omega _z`$ for different $`m`$ is shown in Fig.3b. The maximum repulsion at $`\omega _\rho /\omega _z=0`$ corresponds with $`\omega _z\mathrm{}`$ to the 2D case. The 3D case without magnetic field starts for our parameters $`\omega _\rho /\omega _z=1/3`$ at some value $`\lambda _{\mathrm{eff}}/\lambda <1`$ which decreases further for increasing $`\omega _\rho /\omega _z`$, i.e., increasing magnetic field. This explains quantitatively through the effective charge the difference of the effect of a magnetic field on a quantum spectrum in 2D and 3D cases. However, this difference becomes weaker for larger $`m`$ as it is seen in Fig.3b.

Although the ground state as a function of the magnetic field is formed piecewise by levels of increasing $`m`$ and alternating singlet-triplet character (see Fig.3a) the magnetic properties of the ground state nevertheless reveal the dimensional difference between 2D and 3D. At temperature $`T=0`$ the dot is in the ground state and the magnetic moment and the magnetic susceptibility are defined by $`\mu _{\mathrm{mag}}=E_{\mathrm{gr}}/B`$ and $`\chi =\mu _{\mathrm{mag}}/B`$, respectively. Both quantities exhibit discontinuities as a function of the magnetic field due to the symmetry changes (with respect to $`m`$ and spin). We find that these discontinuities shift when going from the 2D quantum dot to the 3D case as shown in Fig.4.

By relaxing the restriction of two dimensions for a quantum dot and working in the physical three-dimensional space we have investigated physical examples of non-integrable systems close to integrability. For this situation the RRM method is naturally justified, since the confining frequencies in quantum dots obey the condition $`\mathrm{}\omega _z\mathrm{}\omega _0`$. The WKB-approach provides a simple and transparent way to calculate the spectrum of the 3D two-electron quantum dot even for small values of the quantum numbers. We have found that at specific values of the magnetic field $`\omega _L^{}=\sqrt{\omega _z^2\omega _0^2}`$ an axially-symmetric quantum dot exhibits spherical symmetry and its dynamics becomes completely separable with three integrals of motion and three corresponding quantum numbers. We have shown that the confinement in the $`z`$-direction, neglected in the 2D description of quantum dots, does have an influence on the spectrum and consequently also on magnetic properties of the dot. In fact, by changing the confining frequency in the $`z`$-direction only slightly one can increase or decrease the magnetic moment and the susceptibility, i.e. one can control the magnetic properties of the two-electron quantum dot.

Figure captions:

Fig.1. Poincaré surfaces of sections $`z=0`$, $`p_z>0`$ of the relative motion for the axially-symmetric 3D two-electron quantum dot ($`\omega _z/\omega _0=3`$, $`\lambda =3`$, $`m=0`$, $`ϵ=5`$) in the magnetic field for: (a) $`\omega _L=0`$, (b) $`\omega _L/\omega _0=2.5`$, (c) $`\omega _L/\omega _0=\sqrt{8}`$ and (d) $`\omega _L/\omega _0=3.3`$. The section (c) indicates that for the corresponding value of the magnetic field the system is integrable.

Fig.2. The energy spectrum of a circular 2D quantum dot (in units $`\mathrm{}\omega _0`$) as a function of the ratio $`\omega _L/\omega _0`$ for $`n_\rho =0`$ and $`m=0,\mathrm{},9`$ in the cases: (a) $`\lambda =0`$ and (b) $`\lambda =3`$.

Fig.3. (a) The comparison between energy levels (in units $`\mathrm{}\omega _0`$) of the axially-symmetric 3D quantum dot with $`\omega _z/\omega _0=3`$ and $`\lambda =3`$ for $`n_\rho =n_z=0`$ and $`m=0,\mathrm{},9`$ obtained using the RRM (full lines) and exact results for the spherical case (circles). The inset shows a good agreement between the RRM and the exact results. The dashed and dotted lines display the energy level with $`m=0`$ for the 2D and 3D cases with $`\lambda _{\mathrm{eff}}`$ at $`\omega _0/\omega _z=0`$ and $`1/3`$, respectively. (b) The dependence of the effective strength of the Coulomb interaction $`\lambda _{\mathrm{eff}}/\lambda `$ on the ratio $`\omega _\rho /\omega _z`$.

Fig.4. Magnetic moments $`\mu _{\mathrm{mag}}`$ (a) in the units of effective Bohr magneton $`\mu _B^{}=(m_e/m^{})\mu _B`$ and the magnetic susceptibility $`\chi `$ (b) for the 2D (dashed lines) and 3D (full lines) cases as a function of the magnetic field strength (in $`\omega _L/\omega _0`$-units). We use the same parameters as in Figs.2,3.


text:
# Unzipping double-stranded DNA with a force: scaling theory

## Abstract

A double stranded DNA molecule under the stress of a pulling force acting on the strand terminals exhibits a partially denatured structure or can be completely unzipped depending the magnitude of the pulling force. A scaling argument for relationships amongst basic length scales is presented that takes into account the heterogeneity of the sequence. The result agrees with our numerical simulation data, which provides a critical test of the power laws in the unzipping transition region.

The characteristics of the structure of double stranded DNA (dsDNA) have been recently probed through mechanical manipulation at a microscopic level by applying an external force. From a theoretical modelling perspective, we have just begun to understand these experimental results, which are usually provided in the form of manipulating force verses extension. Most of the current theoretical work is based a homopolymer model, in which the binding forces between nucleotides in various types are treated equally . One of the theoretical efforts is to reproduce the force-extension phase diagram observed experimentally, together with other interesting conformational properties that characterize the unzipping transition. It is certainly a challenging task to understand the system-dependent and the universal features in these complex systems, as the pairing along the two DNA strands is in general heterogeneous rather than homogeneous.

The present article serves two purposes. Firstly, a physically direct scaling argument is employed to analyze a simplified heterogeneous dsDNA model, which exhibits an unzipping transition under the external pulling force acting upon the terminal end. The usefulness of some well-tested scaling arguments in polymer theory is demonstrated, both in terms of a Gaussian model and of a freely jointed bond model that correctly addresses the finite extensibility of each strand. Secondly, based on a continuous description of the problem , numerical simulations are presented to test some of the obtained power laws, in particular those near the unzipping transition. A special attention is paid to display the sequence-dependent features verses averaged, sequence-independent features.

To start with, we consider two bounded DNA strands as shown schematically in Fig. 1, where some segments interact with each other predominantly with repulsions by displaying excursion loops, and some segments predominantly with attractions in a mostly bound form. To simplify the consideration, we assume here that the interaction energy per base pair when they are confined in the bounding distance can either be $`+ϵ`$ (repulsive) or $`ϵ`$ (attractive). Furthermore, we assume that the averaged free energy per base pair is $`|f_0|`$ when the dsDNA is below the melting temperature with no additional pulling force.

When an external pulling force is applied to the two terminal nucleotides (or two terminal groups) of the dsDNA, we assume that on average $`M`$ base units can be pulled out. Hence, the reduced free energy of the entire chain reads

$$\beta =(NM)\beta |f_0|+\frac{D}{2}\frac{Z^2}{Ma^2}\beta FZ\sqrt{M}\beta ϵ.$$
(1)
where the first term describes the fact that there are $`NM`$ base pairs that still remain in the bound state on average. The three additional terms describe the effects of the pulling process that results in $`M`$ separated base pairs and an average distance $`Z`$ between the two terminals. One of the essential assumptions here is that the polymer segment in the stretched portion obeys a Gaussian distribution, an assumption that will be questioned below. In order to make a more precise comparison with a previous result, a prefactor $`D/2`$ is explicitly introduced in the second term, to make it exact even at the level of numerical coefficients. The third term is simply a potential energy reduction of the terminal pairs when an external pulling field is applied.

The last term in Eq. 1 describes the binding energy reduction of the pulled portion. When the dsDNA is being pulled, two possibilities might occur near the end of the separated portion. In the event that this portion is connected to an originally bound segment, a stronger force is required to unzip this locked portion. In the event that this portion is connected to a segment that is an excursion loop, the dsDNA would spontaneously unwind the excursion loop. Hence under a fixed pulling force and distance $`Z`$, $`M`$ may vary and the entire separated portion of the $`M`$ base pairs would adjust itself so that it has a net sum of a lower energy. Since the binding energy is fluctuating between $`+\sqrt{M}ϵ`$ and $`\sqrt{M}ϵ`$ randomly the average energy that the pulled portion can acquire is $`\sqrt{M}ϵ`$, where the square root is a result of the randomness in the system.

Now the free energy has to be minimized with respect to $`Z`$ and $`M`$ as adjustable parameters. Minimizing with respect to $`Z`$ leads to

$$Z=a^2\beta FM/D$$
(2)
which has a linear form in $`F`$ and reflects the fact that the pulled portion is near a highly stretched limit. Inserting Eq. 2 into Eq. 1, we have

$$=N|f_0|+tM\sqrt{M}ϵ$$
(3)
where we have defined a reduced force parameter

$$t=|f_0|\beta F^2a^2/(2D).$$
(4)
The above equation is our final expression for the free energy of the system. We see that when $`t0`$, there is no finite minimum in the free energy as $`M`$ tends to increase unboundedly; as a consequence the dsDNA is entirely unzipped. When $`t>0`$ however, there is an equilibrium state, where, by minimizing the above expression with respect to $`M`$, the average number of unzipped base pairs is

$$M=(ϵ/t)^2/4$$
(5)
and the free energy is

$$=N|f_0|ϵ^2/t/4.$$
(6)
The unzipping transition takes place when $`t=0`$, which is equivalent to saying that there exists a critical unzipping force,

$$F_\mathrm{c}=\sqrt{2D|f_0|/(\beta a^2)}.$$
(7)
The main conclusion drawn here, namely, Eqs. 5, 6, and 7, is consistent with a much more involved derivation. The numerical coefficient of the critical force is exact, in view of the fact that the original involved terms in the free energy are exact to start with. The free energy for the entire system is no longer self-averaging, as the last term in Eq. 6 is not proportional to $`N`$. Indeed, from the physical perspective, the equilibrium state of the pulled portion is established regardless of the total length of the dsDNA, as long as the pulled portion is still shorter. The independence of $`N`$ in the last term is merely a statement that for different $`N`$, the equilibrium state of the pulled portion remains essentially the same.

This simple treatment for the unzipping transition of the dsDNA produces several relationships between the basic length scales and the transition properties. However, there is one important feature of the stretched portion that has been ignored both above and in Ref. , that might become a pitfall of the entire derivation. The usage of a Gaussian entropic term for the stretched portion in Eq. 1 can be traced back to the spring-bead model commonly used in polymer modelling. As a spring can be stretched without a bound under a strong pulling force, the current model physically allows the pulled segment to stretch beyond the natural length $`Ma`$, as demonstrated in Eq. (2). Hence, it becomes desirable to incorporate the effect of finite extensibility in the above discussion.

It has been previously demonstrated by Grosberg and Khokhlov  that a better representation for a stretched polymer is the freely jointed bond model that contains a rigid bond that models a polymer monomer unit. Instead of the second and thrid terms in Eq. 1, the contribution to the free energy due to the entropic penalty in the strong stretched limit ($`D=3`$) can be shown to be

$$_{\mathrm{stretch}}=\frac{M}{\beta }\mathrm{ln}[\frac{\mathrm{sinh}(\beta Fa)}{\beta Fa}]$$
(8)
and the separation distance can be directly calculated,

$$Z=Ma[\mathrm{coth}(\beta Fa)1/(\beta Fa)].$$
$`(9)`$
Replacing the the second and the third terms in Eq. 1, we have

$$=(NM)f_0+\frac{M}{\beta }\mathrm{ln}[\frac{\mathrm{sinh}(\beta Fa)}{\beta Fa}]\sqrt{M}ϵ$$
$`(10)`$
The linear dependence of the stretching free energy on $`M`$ essentially renders the same scaling theory as discussed above. Indeed, the free energy expression in Eq. 6 remains the same and the only change is the definition of $`t`$, which now reads

$$t=|f_0|+\frac{1}{\beta }\mathrm{ln}[\frac{\mathrm{sinh}(\beta Fa)}{\beta Fa}].$$
$`(11)`$
We note that this expression recovers the same form as in Eq. 4 in the small $`\beta Fa`$ limit.

Hence, when the finite extensibility is taken into account in the dsDNA model, the critical force can be implicitly solved through finding the root for $`t=0`$ in Eq. 11. The power laws in Eqs. 5 and 6 remain the same.

The rest of this article is devoted to the discussion of numerical simulations that can be used to address system specific features in comparison with these universal features. Returning to the Gaussian description of the conformation for both strands of DNA, we write a more rigorous expression for the probability functional in one dimension, after Lubensky and Nelson (LN):

$$P_0\mathrm{exp}\{\frac{1}{2a^2}_0^N[(\frac{dz}{ds})^2]𝑑s\beta _0^NW[z(s),s]𝑑s\}$$
$`(12)`$
where $`z(s)`$ describes the distance between the base pair labeled $`s`$ in the dsDNA, and $`W`$ is a strong bounding potential, which depends on the (disordered) sequence of nucleotides on both strands. In a Gaussian chain model, the Kuhn length $`a`$ in this reduced probability function is $`\sqrt{2}`$ times the Kuhn length of each single strand of DNA . The unzipping probability is then simply

$$P_{\mathrm{unzip}}[z(s),\{s\}]=P_0[z(s),\{s\}]\mathrm{exp}\{\beta Fz(N)\}$$
$`(13)`$
with the assumption that the pulling force is directed along the $`z`$ direction and acting on the $`N`$th monomer. The $`\{s\}`$ dependence in these functions stresses the fact that the interaction between the base pairs depends on a given sequence.

A mathematically similar problem to the above probability distribution function arises in the study of random copolymer localization in an external field, specifically, at a sharp interface.The basic idea of our numerical simulations is to simulate the disordered heterogeneousness in the system by generating an ensemble of dsDNA models that have explicit bounding interaction dependence. The calculation of the probability function of the LN model is then translated into solving a Schrödinger type equation. For each generated sequence, the Schrödinger equation is solved and thus yields an exact result with no arbitrary approximations.

We have used basically the same numerical procedure described in Ref.  by computing the probability function, $`\psi (z,s)`$, of finding the $`s`$th base pair having a separation distance $`z`$. Now, the calculation of $`\psi (z,s)`$ in a potential field can be simplified to solving a time-dependent Schrödinger equation,

$$\frac{\psi (z,s)}{s}=[\frac{a^2}{2}\frac{d^2}{dz^2}+\beta W(z,s)]\psi (z,s)$$
$`(14).`$
where the “wave function” $`\psi `$ is subject to the initial condition $`\psi (z,0)=1`$. The probability function for the terminal pair to have a separation distance $`z`$, $`q(z,\{s\})`$, is given by $`q(z,\{s\})=\psi (z,N)`$. At this stage, $`q(z,\{s\})`$ has an explicit sequence dependence.

For further numerical development we must choose a function form for the interaction potential and we write

$$W[z,\{s\}]=\zeta (s)ϵw(z)$$
$`(15)`$
with $`w(z)`$ being modelled by a Gaussian potential well $`w(z)=\mathrm{exp}((z/a)^2)/\sqrt{\pi }`$. The coefficient $`\zeta (s)`$ represents the binding energy of different pairs with $`ϵ`$ assumed here positive. In the remainder of this article we describe numerical results of having chosen $`\beta ϵ=1`$.

There are in total 10 possible pairs of nucleotides in DNA, corresponding to 10 different values of $`\zeta (s)`$, where guanine-thymine and adenine-cytosine pairings are known to be attractive. Whether or not the disordered sequence of bases along DNA is truly random is a still debatable question. However, for the present purpose, it suffices to assumed that $`\zeta (s)`$ takes a statistically independent random value for each $`s`$. In particular, we assume $`\zeta (s)`$ is a random number uniformly ranging from $`1`$ to $`1`$ as a theoretical abstraction. The numerical data for a biased randomness will be discussed elsewhere.

Figure 2 shows the relationship between the force and the separated distance $`\overline{z}`$ for 10 different, randomly generated disorder sequences. In producing the simulation data, a finite system with a limited size of $`z<1600`$ was considered; this finite-size effect shows up in all curves as they cannot exceed a separation distance of 1600. A length of $`N=32000`$ units has been used here, which guarantees a minimum finite-size effect due to the modelled DNA length. The step-like functions in this plot are a signature of the sequence disorderliness. This figure is a clear indication that for a specific dsDNA model (and indeed for real dsDNA) a variety of force-distance curves can be obtained with the unzipping transition occur near an averaged critical force. The critical force for each sequence, however, varies. For reference, the anticipated critical force, averaged over the sequence ensemble, is around $`\beta F_ca=0.3724`$ from an accurate estimate of the bound free energy $`f_0`$ of this system (see Eq. 7).

Though each sequence might show strong individual characteristics, the universal behavior around which the individual properties fluctuate, can be obtained by taking an average over the sequence ensemble. In order to do so, we have produced $`10^3`$ sample DNA models, each being examined through the solution of the Schrödinger equation. The average separation, $`<\overline{z}>`$, as a function of the pulling force $`F`$, is presented in Fig. 2 as the thick dark curve, where the average $`<\mathrm{}>`$ represent sampling in the disorder ensemble.

The same curve is also presented in Fig. 3 by circles as a function of $`F_\mathrm{c}F`$, in a double logarithmic scale. We see that within the region of $`[0.03,0.02]`$, $`<\overline{z}>`$ indeed follows a power behavior,

$$<\overline{z}>(F_\mathrm{c}F)^{q_1},$$
$`(16)`$
with a fitted slope of $`q_1=0.194`$, directly verifying the scaling exponent of 2 in Eq. 5. Exception of the scaling can be seen in smaller $`F_\mathrm{c}F`$ region where the saturation of $`<\overline{z}>`$ to a finite value occurs. As a direct consequence of the compound  finite size effect shown for various sequences in Fig. 2, the averaged $`<\overline{z}>`$ is more drastically affected by the finite size of the system.

LN have quoted two other power laws. Translating the averages of $`M`$ into the averages of $`z`$ permits us to examine these power laws:

$$\delta z_1^2<\overline{z}^2><\overline{z}>^2(F_\mathrm{c}F)^{q_2}$$
$`(17)`$
and

$$\delta z_2^2<\overline{z^2}><\overline{z}^2>(F_\mathrm{c}F)^{q_3}$$
$`(18)`$
The square and diamond symbols in Fig. 3 represent our simulation data for $`\delta z_1^2`$ and $`\delta z_2^2`$. By fitting straight lines to these curves in the region $`[0.03,0.02]`$, we obtained $`q_2=3.99`$ and $`q_3=2.94`$, which agree well with analytical determination of these exponents.

To summarize, we have carried out a scaling analysis of the unzipping transition phenomenon of a heterogeneous dsDNA model. Our main conclusions, on the free energy and stretching separation, agree with the results drawn by previous authors. These results compare favorably well with the numerical simulation data, which, for the first time, demonstrate the relationship between the sequence-dependent and universal features.

Financial support of this work is provided by the Natural Science and Engineering Research Council of Canada and the grant associated with the Premier’s Research Excellence Award of Ontario, Canada.


text:
# Soft Photon Spectrum in Orthopositronium and Vector Quarkonium Decays

## 1 Contradiction between Low’s theorem
and Ore-Powell differential decay rate

We will concentrate on the orthopositronium decay into three photons. The discussion also applies to vector quarkonium ($`J/\psi `$, $`\mathrm{{\rm Y}}`$,…).

From the requirement of gauge invariance, and because of the quantum numbers of the initial and final particles involved (neutral self-conjugate bosons), Low’s theorem predicts that the decay amplitude must vanish linearly when the energy of one of the photons is going to zero

$$\left(o\text{-}Ps\gamma \gamma \gamma \right)\stackrel{\omega 0}{}𝒪\left(\omega \right)$$
(1)
with $`\omega `$ the energy of one of the photon. The squared modulus of the amplitude therefore behaves as $`𝒪\left(\omega ^2\right)`$ for small photon energy. Aside from that, the three-photon phase-space alone (i.e. with a constant decay amplitude) gives a differential rate as

$$\frac{d\mathrm{\Gamma }\left(o\text{-}Ps\gamma \gamma \gamma \right)}{dx}|_{PhaseSpace}x$$
(2)
with $`x=2\omega /M`$ the reduced photon energy, $`M`$ the orthopositronium mass. Combining the amplitude with the phase-space, one finds that the low-energy end of the photon spectrum must behave as

$$\frac{d\mathrm{\Gamma }\left(o\text{-}Ps\gamma \gamma \gamma \right)}{dx}\stackrel{x0}{}x^3$$

No matter the model used to compute the orthopositronium decay rate, the differential rate must exhibit this $`x^3`$ behavior for small photon energies. As we have already pointed out in , the lowest order decay amplitude, as found by Ore and Powell, leads to a differential decay rate in contradiction with the analytical requirement of Low’s theorem. Their model is based on the formula :

$$\mathrm{\Gamma }\left(o\text{-}Ps\gamma \gamma \gamma \right)=\frac{1}{3}\left|\psi \left(0\right)\right|^2\left(4v_{rel}\sigma \left(e^{}e^+\gamma \gamma \gamma \right)\right)_{v_{rel}0}$$
(3)
with $`v_{rel}`$ the relative velocity of the $`e^+e^{}`$ in their center-of-mass frame, and $`\psi \left(0\right)`$ the positronium Schrödinger wavefunction at zero separation ($`m`$ is the electron mass)

$$\left|\psi \left(0\right)\right|^2=\frac{\alpha ^3m^3}{8\pi }$$
(4)
This formula states that in first approximation, the positronium decay rate can be computed from the static limit of the scattering cross section $`e^+e^{}\gamma \gamma \gamma `$. Equivalently, it is found from the squared modulus amplitude for an $`e^+e^{}`$ pair at rest into $`\gamma \gamma \gamma `$. Summed over photon polarizations, this is easily shown to be

$$\underset{polarizations}{}\left|\left(\left(e^+e^{}\right)_{v_{rel}=0}\gamma \gamma \gamma \right)\right|^2=\frac{\left(1x_1\right)^2}{x_2^2x_3^2}+\frac{\left(1x_2\right)^2}{x_1^2x_3^2}+\frac{\left(1x_3\right)^2}{x_1^2x_2^2}$$

which behaves as a constant when one of the $`x_i`$ is vanishing (as can be seen by using energy-momentum conservation $`x_1+x_2+x_3=2`$), while it should vanish as $`x_i^2`$. In turn, the well-known differential rate inherits an incorrect analytical behavior

$$\frac{d\mathrm{\Gamma }\left(o\text{-}Ps\gamma \gamma \gamma \right)}{dx_1}=\frac{2\alpha ^6m}{9\pi }\mathrm{\Omega }\left(x_1\right)$$

where the spectrum function is

$`\mathrm{\Omega }\left(x_1\right)`$ $`={\displaystyle _{1x_1}^1}𝑑x_2\left|\left(\left(e^+e^{}\right)_{v_{rel}=0}\gamma \gamma \gamma \right)\right|_{x_3=2x_1x_2}^2`$ (5)
$`={\displaystyle \frac{2\left(2x_1\right)}{x_1}}+{\displaystyle \frac{2\left(1x_1\right)x_1}{\left(2x_1\right)^2}}+4\left[{\displaystyle \frac{\left(1x_1\right)}{x_1^2}}{\displaystyle \frac{\left(1x_1\right)^2}{\left(2x_1\right)^3}}\right]\mathrm{ln}\left(1x_1\right)`$
$`={\displaystyle \frac{5}{3}}x_1+𝒪\left(x_1^2\right)\text{ near }x_1=0`$
In the Ore-Powell model, the photon energy spectrum vanishes only linearly near zero, instead of the required $`\mathrm{\Omega }\left(x_1\right)=𝒪\left(x_1^3\right)`$.

Remarks

For completeness, recall that it is this differential rate that gives the total width

$$\mathrm{\Gamma }\left(o\text{-}Ps\gamma \gamma \gamma \right)=\frac{2\left(\pi ^29\right)}{9\pi }\alpha ^6m\text{ since }_0^1𝑑x_1\mathrm{\Omega }\left(x_1\right)=\pi ^29$$

Up to color factors, wavefunctions, coupling constants, the present analysis can be repeated for the quarkonia decay modes

$`{\displaystyle \frac{d\mathrm{\Gamma }\left(V\gamma \gamma \gamma \right)}{dx_1}}`$ $`={\displaystyle \frac{64}{3}}e_Q^6\alpha ^3{\displaystyle \frac{\left|\varphi _0\right|^2}{M^2}}\mathrm{\Omega }\left(x_1\right)`$
$`{\displaystyle \frac{d\mathrm{\Gamma }\left(Vggg\right)}{dx_1}}`$ $`={\displaystyle \frac{160}{81}}\alpha _S^3{\displaystyle \frac{\left|\varphi _0\right|^2}{M^2}}\mathrm{\Omega }\left(x_1\right)`$
$`{\displaystyle \frac{d\mathrm{\Gamma }\left(Vgg\gamma \right)}{dx_1}}`$ $`={\displaystyle \frac{128}{9}}e_Q^2\alpha \alpha _S^2{\displaystyle \frac{\left|\varphi _0\right|^2}{M^2}}\mathrm{\Omega }\left(x_1\right)`$
where $`V`$ is the $`1^{}`$ vector bound state made of the $`Q\overline{Q}`$ pair, $`M`$ is the mass of $`V`$, $`\varphi _0`$ the (unknown) quarkonium configuration space wavefunction at zero separation and $`e_Q`$ the heavy quark electric charge in units of the electron one. All these decay spectra do violate the basic requirement of analyticity.

## 2 Contradiction between Low’s theorem
and current positronium decay models

The model used by Ore-Powell, based on the factorized formula (3), may seem a bit naive, especially in view of the enormous amount of work done by various groups (see for example , , , , and references quoted there). Nevertheless, it is very illustrative of the current approaches in its treatment of intermediate states. Indeed, models derived from Bethe-Salpeter analyses, or from QED non-relativistic effective theory (NRQED, ), always connect the process of annihilation of bound charged particles to that of scattering of real, asymptotic charged particles. The difficulty with such approaches is thereby apparent: asymptotic and bound charged particles have drastically different radiation properties: the former exhibit bremsstrahlung-type radiations, while the later do not radiate zero energy photons (for very low-energy photons, a positronium state is just a neutral, self-conjugate boson, hence it does not radiate in that limit).

Bremsstrahlung radiations typically lead to a Laurent expansion for the amplitude as

$$^\mu (\omega ,\mathrm{})=𝒪\left(1/\omega \right)+𝒪\left(1\right)+𝒪\left(\omega \right)$$

What Low’s theorem state is that both the terms of $`𝒪\left(1/\omega \right)`$ and $`𝒪\left(1\right)`$ must disappear . While the cancellation of $`𝒪\left(1/\omega \right)`$ terms is automatic from selection rules, that of $`𝒪\left(1\right)`$ terms is much more delicate, requiring a non-perturbative treatment of the binding energy. Typically, the $`𝒪\left(\omega \right)`$ term is of the form

$$^\mu (\omega ,\mathrm{})\stackrel{\omega 0}{}\omega \left(\frac{M^2}{M^24m^2}+\text{regular terms as }M2m\right)$$
(6)
Obviously, a perturbative expansion in the binding energy $`M2m`$ is mathematically inconsistent with the soft-photon expansion: if the limit $`M2m`$ is taken before $`\omega 0`$, spurious $`𝒪\left(1\right)`$ terms arise. Because the basis of current computations is a perturbative expansion in the binding energy $`M2m`$, computed as relativistic and radiative corrections to the Ore-Powell result, one can expect that spurious radiations affect the corrections presented in the literature.

As a comment, notice that the present considerations apply to amplitudes. For orthopositronium, any term less singular than $`𝒪\left(1/\sqrt{\omega }\right)`$ would lead to an IR finite decay rate, because of (2). In other words, unphysical terms can lead to IR finite contributions to the decay rate. Computations presented in the literature sometimes explicitly exhibit such a bad behavior .

To illustrate how the Low’s theorem is implemented in the presence of charged intermediate states, let us give an example: $`K_S^0e^+e^{}\gamma `$. The process is modeled by a charged pion loop. The resulting amplitude can be found in many places (see for example , , ), and it behaves exactly as predicted by Low’s theorem (i.e. as $`\omega ^3`$ near $`\omega =0`$), even if in this case $`M_K>2m_\pi `$, i.e. the intermediate charged pion pair can be on-shell.

Now, let us imagine that one is willing to compute the decay rate for $`K_S^0e^+e^{}\gamma `$ by assuming that intermediate on-shell $`\pi ^+\pi ^{}`$ dominates. The decay process is then factorized as $`K_S^0\pi ^+\pi ^{}`$ times $`\pi ^+\pi ^{}e^+e^{}\gamma `$. This is exactly the approximation done to get the Ore-Powell result: $`o`$-$`Pse^+e^{}`$ times $`e^+e^{}\gamma \gamma \gamma `$. However, the soft-photon spectrum of the factorized approximation is completely wrong, being in contradiction with Low’s theorem.

The approximation done was too stringent. To get the correct answer, one must also consider processes like $`K_S^0\pi ^+\pi ^{}\gamma `$ times $`\pi ^+\pi ^{}\left(\gamma \right)e^+e^{}\left(\gamma \right)`$ (where the photon is disconnected). These bremsstrahlung processes interfere destructively with the previous ones, giving a finite vanishing complete amplitude in the soft-photon limit. In the framework of dispersion relations , this is a simple application of the Cutkosky rule to get the absorptive part of an amplitude. Similarly, we state that the reason why the Ore-Powell result fails to exhibit a correct soft-photon spectrum is because some contributions to the amplitude are missed (like $`o`$-$`Pse^+e^{}\gamma `$ times $`e^+e^{}\left(\gamma \right)\gamma \gamma \left(\gamma \right)`$, see , ).

Higher order corrections to the Ore-Powell result   are similarly incomplete. Trying to take binding energy effects into account with models based on any kind of factorization $`o`$-$`Pse^+e^{}`$ times $`e^+e^{}\gamma \gamma \gamma `$ is hopeless in view of Cutkosky rules. The fact that the on-shell intermediate state method ”works” (i.e. is not IR divergent at the decay rate level) for orthopositronium is not a valid argument, since the amplitudes produced by that method fails to fulfill a basic requirement of QED, namely Low’s theorem.

## 3 Conclusions and Perspectives

From a theoretical perspective, as shown in , the introduction of additional contributions (i.e. process like $`o`$-$`Pse^+e^{}\gamma `$ times $`e^+e^{}\left(\gamma \right)\gamma \gamma \left(\gamma \right)`$) to the positronium decay amplitude is unavoidable if one is willing to fulfill the basic requirement of Low’s theorem. As we have explained, the $`\omega ^3`$ low-energy spectrum is a consequence of the properties of the positronium or quarkonium ”as seen from far away”, i.e. as a neutral point-like self-conjugate bosonic particle. Such a particle does not radiate zero-energy photons, and the resulting photon spectrum must exhibit a $`\omega ^3`$ shape near $`\omega =0`$. This is equally true for $`J/\psi `$, $`\mathrm{{\rm Y}}`$…, and our solution is naturally extended to quarkonium theory.

On the practical side, it appears quite obvious that the violation of Low’s theorem is very small in positronium decay. In other words, the missing contributions are subleading. As discussed in , we can expect them to introduce corrections of the order of the binding energy $`E_B=M2m`$, i.e. $`\alpha ^2`$, or beyond. Nevertheless, even if very small, those corrections are relevant to the current theoretical considerations. Indeed, it is at the $`\alpha ^2`$ level that some discrepancies have been found among experiments, the so-called orthopositronium lifetime puzzle , , . What we claim is that no definite answer to that puzzle could be given at present. Indeed, the additional contributions could turn out to be less suppressed than usually thought, and the current theoretical result for the $`\alpha ^2`$ correction is not fully reliable.

This state of affair is to be contrasted to the quarkonium case. There, the missing contributions could become sizeable since the binding energy is non-negligible (see (6)). In other words, the relevance of the Ore-Powell result for quarkonia is doubtful for any precision calculation.

For instance, the photon spectrum in inclusive quarkonium decay into hadrons + photon  will clearly be affected by the missing contributions. Indeed, the modification of the spectrum needed at low energy to fulfill Low’s theorem will affects the spectrum also at high energy (since when the photon has its maximum energy, one of the gluon’s energy can go to zero, see the integration ranges in (5)).

Also, the missing contributions should be crucial to solve the $`\rho \pi `$ puzzle , . The so-called $`14\%`$ rule is obtained from the ratio of the leptonic mode of the $`J/\psi `$ and $`\psi \left(2S\right)`$ (which is essentially the ratio of wavefunctions at zero separation). There is no missed contribution for the leptonic modes. On the other hand, at least $`12`$ additional contributions need to be considered for three-gluon modes, and those will depend on the binding energy, which is quite different for $`J/\psi `$ and $`\psi \left(2S\right)`$. Those additional contributions, which arise already at the lowest order, are essential to enforce Low’s theorem through their destructive interferences with the standard factorized Ore-Powell ones.

Finally, the previous remark also shows that the extraction of $`\alpha _S`$  from quarkonia branching fractions will be affected, at least partially, by the additional contributions.

In conclusion, we have shown that the implications of gauge invariance and analyticity, in the form of constraints on the low-energy end of photon spectrum are not met by current bound state decay models. Restoring a correct behavior in that low-energy region could lead to potentially interesting advances in both QED and QCD bound state description.

Acknowledgments: We are very pleased to acknowledge useful discussions with Gabriel Lopez Castro, Jean-Marc Gérard, Stéphanie Trine and Jacques Weyers. C. S. acknowledges financial support from FNRS and IISN (Belgium).