text:
# 1 Introduction

## 1 Introduction

Let us briefly review the present state of the quark mass problem, considering that, modern physics has enough experimental data to be certain in existence of quarks. Quarks acquire their masses in a process which is related to the Higgs sector and goes back to the original Higgs’s work of 1964 . Quarks are interacting with the Higgs boson (bosons) proportionally to masses of quarks. The extensive review of the Higgs boson search and theories presented in . The Higgs boson is believed to have the mass $`115\text{ GeV}`$  or perhaps even up to $`200\text{ GeV}`$ . In view of indeterminacy of a theoretical position, both the number of bosons and their masses remain unknown and therefore the concrete mechanism of quark mass formation is also unclear. As twenty years ago  and now  practically the same opinion is stated:

“… we do not at all understand the lepton and quark mass spectrum …”

“The fact, that the mass spectrum of quarks and leptons envelops 5 orders of magnitude, is a mystery from any point of view… Even with universal recognition of standard model, these problems seem very far from resolution” . Theory urges us to employ the concept of a quark mass with a great deal of discretion, because only bound quark systems (i. e. particles) are observable. Depending on the concrete character of quark interaction, we can expect either positive or negative binding energy. The negative binding energy is analogous to mass defect in nucleus when the mass sum of the component quarks is more than the mass of the relevant hadron. The positive binding energy signifies that in addition to quark masses there is an additional energy in a hadron, for example, the energy of color field, which provides the quark confinement. The notion of quark mass and its value depends upon the theoretical model used in a given experimental or theoretical work. With respect to quark mass modern physics developed a number of notions: there exist bare quark mass, physical quark mass, constituent quark mass and effective quark mass. Effective quark mass is inherently related to the impulse used to calculate it  and we won’t consider it here. An interacting field theory has a well-defined notion of the rest mass of a single particle state, which is called physical mass. There can also be discrete bound states of the particle-antiparticle pairs.

The mass parameter in the QCD Lagrangian is considered to be a bare quark mass . The various estimations within the QCD theory using the $`\overline{MS}`$-schemes and lattice simulations show that the bare masses of the light $`u`$\- and $`d`$-quarks are less than $`10\text{ MeV}`$ . The estimation of the bare mass ratios for the $`u`$-, $`d`$\- and $`s`$-quarks is $`1:2:40`$ respectively. Therefore, it is accepted, that the fundamental Lagrangian of strong interaction does not contain any hint on a symmetry with respect to quark masses .

Approximate constituent quark masses can be obtained from experimental data of the hadron masses . Constituent quark is deemed to be formed from the bare quark and its “coat”, i.e. surrounding gluon field energy, the sea quarks and all other unsolved hadron structure problems. First of all, the “coat” notion makes the concept of a constituent quark mass fuzzy, because the quark’s incapability of solitary existence makes the mass of it s “coat” dependent on a particular environment, i.e. from the quark structure of hadron and the quark-gluon dynamics.

On the other hand, the mass discreteness of the light hadrons with equal quantum numbers is several times greater than pion masses. It is difficult to understand why a light constituent quarks of pions have to dress in additional heavy “coat”. At the same time significant discreteness of the hadron mass spectrum makes it possible to hope that the quark-gluon field energy is also subject to quantization.

On the whole, after the forty years the birth of the quark model, the solution of the quark mass problem seems to be far still. It is also apparent that both the hadron internal energetic structure and quark-quark interaction dynamics is impossible to fully understand without the knowledge of quark mass quantization and the precise values of their spectrum. These reasons and the author’s conviction in the existence of the analytic relation for quark masses gave rise to the present work. Hereinafter under a quark mass we shall understand a physical rest mass of the single particle state of some interacting quantum field.

## 2 Quark properties and experimental data

As a first stage of the study we shall ascertain the fact that there is both theoretical reason and experimental data in favor of certain relation between quark masses. However, let’s begin with a few words about flavor as it is the only known quantum number, which is obviously connected with quark mass. All other quantum numbers is strictly quantized and are either constant or periodic by quark generations. Let’s ask ourselves, is quark mass quantized? It evidently is, for all available energetic information about matter indicates that quantum effects and energy discrete levels became dominant as an object dimensions decrease. The mass spectrum of the hadrons with different flavors support the idea of quark mass quantization. Therefore another question looms: is quark mass strictly quantized? According to field theory, the answer is yes. The flavor quantum number is essentially a reflection of quark’s internal energy—its physical mass. Let’s see how the flavor changes in the basic weak interactions of quarks.

$`Q_n+W^\pm `$ $``$ $`Q_{n+1}`$ (1)
$`Q_n+e^\pm `$ $``$ $`Q_{n+1}+\nu `$ (2)
These are transition reactions between quarks with neighboring flavors. It is irrelevant whether real or virtual bosons the quarks are interacting with. The point is that both the quark’s mass and charge are changed. Therefore, two successive transformations will result in double mass transformation with the charge remaining constant—that is, flavor will be the only quantum number that will change.

We shall now examine theoretical computations of bare quark masses and certain mesons . The values of bare quark masses are known with one digit precision at best. Let’s select for the analysis the charged mesons consisting of two neighboring quarks $`Q_n`$ and $`Q_{n+1}`$ with minimal mass in their ground state. In the table 1 these masses are given together with their relations 111The table 1 data are taken from the Particle Data Group listing .

Table 1. Masses of bare quarks and certain mesons.

| quark | quark mass | quark mass | meson | meson mass | meson mass |
| --- | --- | --- | --- | --- | --- |
|  | MeV | ratio |  | MeV | ratio |
| $`d`$ | 3–9 | - | - | - | - |
| $`u`$ | 1–5 | $`0.5`$ ($`u/d`$) | $`\pi ^\pm `$($`du`$) | 139.57 | - |
| $`s`$ | 75–170 | $`40`$ ($`s/u`$) | $`K^\pm `$($`us`$) | 493.64 | 3.54 ($`K^\pm `$/$`\pi ^\pm `$) |
| $`c`$ | 1150–1350 | $`10`$ ($`c/s`$) | $`D_s^\pm `$($`sc`$) | 1968.5 | 3.99 ($`D_s^\pm `$/$`K^\pm `$) |
| $`b`$ | 4000–4400 | $`3.4`$ ($`b/c`$) | $`B_c^\pm `$($`cb`$) | 6400 | 3.25 ($`B_c^\pm `$/$`D_s^\pm `$) |
| $`t`$ | 174300 | $`42`$ ($`t/b`$) | - | - | - |

The bare quark mass ratios do not show any regularity but the meson mass ratios definitely present certain degree of uniformity. The mean value of their mass ratios is 3.6. Assuming that the neighboring quark mass ratio is actually a constant of 3.6 we can obtain the following quark masses (starting from the mean table value for the $`b`$-quark): 4200($`b`$), 1167($`c`$), 324($`s`$), 90($`u`$) and 25($`d`$). The sum of the neighboring quark masses will be about 17% less than the meson masses for all considered mesons except $`D_s^\pm `$ ($`25\%`$). This sounds promising, for after a proportional increase of the quark masses we’ll have a match for three of the four examined mesons. Thus, taking the $`b`$-quark mass value equal to $`5000\text{ MeV}`$ we’ll have for the mesons $`\pi ^\pm `$, $`K^\pm `$ and $`B_c^\pm `$ the mass values of 137, 493 and $`6390\text{ MeV}`$ respectively.

Other examples that meson mass ratio is roughly constant are given in the table 2.

Table 2. Meson masses and mass ratios

| meson | meson mass, MeV | meson mass ratio |
| --- | --- | --- |
| $`\omega `$ | 782 | - |
| $`J`$ | 3097 | 3.96 $`(J/\omega )`$ |
| $`{\rm Y}(1S)`$ | 9460 | 3.05 $`({\rm Y}/J)`$ |
| $`\phi `$ | 1019.4 | - |
| $`\Psi (2S)`$ | 3686 | 3.61 $`(\Psi /\phi )`$ |
| $`{\rm Y}(2S)`$ | 10023 | 2.72 $`({\rm Y}(2S)/\Psi )`$ |

In contrast to table 1 the mesons given here have a hidden flavor (that is, a quark-antiquark pair). Thus, we can claim that this experimental data set also supports the idea of approximate constancy of the neighboring quarks mass ratio ($`3.6`$).

So then, there is certain multiplicative pattern in mass transformation between quark flavors. This probably means that quarks are in a similar state in these mesons with various masses. This is quite hard to explain considering that quark and meson masses differ by two orders of magnitude. The light quarks of a pion should be relativistic whereas the heavy quarks of a heavy meson are not very relativistic even when taken at their bare masses. There would be no such contradiction if the quark masses themselves have the scale similarity with the same factor that charged mesons do. Then the quarks should not be very relativistic in all these mesons.

## 3 Construction of the quark mass model

Now we will attempt to build a quark transformation model and to calculate the quark mass transformation factor. We should investigate the quark interaction process that is essentially to guess the new key moments of the strong interaction. A couple of simple question to start things off: first of all—what is the character of quark-antiquark interaction, and second of all—what is the result of that interaction? It is rather easy to answer the second question.

This interaction produces either a meson based on the quark-antiquark pair (e.g. a vector meson) or the complete annihilation of the pair with the birth of photons or smaller mass quarks and other particles. The formation of vector meson has only deferred the quark-antiquark annihilation. There are seemingly two extremes—either the pair is present (meson) or it disappears (annihilation). But this “yes or no” situation can be avoided by introducing the concept of a discrete annihilation spectrum or discrete coupling of a quark-antiquark pair. It is not much of a novelty actually. A good old discrete bound states of an interacting field theory, which are lower than the doubled physical mass of the quark. Bound quark-antiquark states are well known (e.g., their paired states are presented in table 2). We are approaching the main question now—what is the binding energy value of quark-antiquark bound state and what quantum laws operate by the levels of binding energy? It is still an open question now as we poorly know even the quark masses. It is similar a movement on the closed circle and our chance only to guess. Then new variant of the answer is below stated.

Proposition #1

The quark-antiquark pair of the same flavor $`n`$ can, instead of the complete annihilation, annihilate partly to the bound state of a harmonic oscillator with the mass equal to $`m_n\frac{4}{\pi }`$.

The partial annihilation means that the quark pair with total mass $`2m_n`$ are involved in such an interaction with the following equation for their bound state with total energy $`m_{sum}`$:

$$m_{sum}=2m_n_0^\pi \frac{sin(x)}{\pi }𝑑x=\frac{4}{\pi }m_n$$
(3)
This process (let us name it the harmonic quark oscillation) can be also interpreted as a harmonic oscillation of the quark-antiquark pair over space-time from complete separation of quarks to their complete annihilation. So for example, the time-dependent quadratic mass term can be noted as the equation:

$$m^2(\tau )=(2m_n)^2cos^2(\omega \tau )$$
(4)
The time-averaged mass will also be equal to $`m_n\frac{4}{\pi }`$. In eq. (4) we are not considering the transition of the mass into other energy forms and these oscillations should be forbidden by the law of energy conservation. Nonetheless in the virtual meaning with the help of uncertainty relation we have a right to consider them to a certain boundary frequency $`\omega `$.

The equation (3) is also corresponds to a stable state of the quark pair with a certain mass defect value ($`36\%`$). This state could be similar to quarkonium solution of the wave equations with quarks mutually revolving on some stable orbit, except that it has the very high binding energy.

At this stage, it’s important to examine the stability of the harmonic oscillation. The process defined by (4) would not be stable for exactly the same reason which has allowed not to forbid it. Namely, that we are already at the energy fluctuation level comparable to the quark mass in the context of uncertainty relation, which means that the other fluctuations would immediately destroy the harmonic oscillator. But what will follow the destruction? There are at least two ways to overcome this instability. The first one is that an unstable annihilating oscillator is fast damped and then quark pair goes to mentioned stationary state, which corresponds to equation (3). The second way is to accept the proposition two.

Proposition #2

The labile symmetric state of the harmonic annihilating oscillator based on a quark-antiquark pair of the same flavor can be broken in one or another way and would subsequently pass into the asymmetric state based on quarks of different kinds, i.e. consisting of a quark-antiquark pair with neighboring flavors.

To clarify this proposition let us write one of the possible reaction
of $`u`$-antiquark birth.

$$(s\overline{s})_{ho}+e^{}=s+\overline{u}+\nu $$
(5)
where $`(s\overline{s})_{ho}`$ is a harmonic oscillator based on a strange quark-antiquark pair and $`\nu `$ is a neutrino. Reactions similar to (5) can be written down for any quark harmonic oscillator, both with actual and virtual electrons. Doubling the number of participants we can write the reactions with quarks only:

$$2(s\overline{s})_{ho}=s+\overline{s}+u+\overline{u}$$
(6)
This reaction demonstrates the birth of the quark-antiquark $`u`$-pair from two harmonic oscillators based on strange quark-antiquark pairs. The main point here is the possibility of such reactions with neither extra participants nor energy addition. Note that, it is weak interaction with participation of other particles that is responsible for the symmetry breaking in one harmonic oscillator with the subsequent transition into the asymmetric state.

Now we can write down a recurrent relation between quarks of neighboring flavors, considering that a harmonic state with the energy value of $`\frac{4}{\pi }m_n`$ generates two quarks with masses $`m_n`$ and $`m_{n1}`$:

$$m_{n1}=\frac{4}{\pi }m_nm_n=m_n(\frac{4}{\pi }1)=0.27323954m_n$$

$$\frac{m_n}{m_{n1}}=\frac{\pi }{(4\pi )}=3.6597924\stackrel{\mathrm{def}}{=}MQ$$
(7)
Thus, having introduced two propositions we managed to obtain an equation, which links together the masses of neighboring quarks. Hereafter using the word quarks, we shall mean the harmonic quarks, i.e. the quarks, which subject to the equation 7.

Let’s summarize the implications of the equation (7):

* Quarks form a rigid chain and their masses are all bound together
* The quarks have a logarithmically equidistant mass spectrum
* The multiplication factor $`MQ`$ between neighboring quark masses is equal to $`3.66`$ (which is practically equal to the mean meson mass ratio (3.6) mentioned above)
* A quark flavor uniquely determines the transformation number

To stress this last point, the (7) can be written as

$$m_n=m_0\left[\frac{\pi }{4\pi }\right]^n$$
(8)
where $`m_0`$ is a hypothetical initial quark mass, and the quark order number $`n`$ is essentially its flavor. Thus we have:

| n | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| quark | $`d`$ | $`u`$ | $`s`$ | $`c`$ | $`b`$ | $`t`$ | $`b^{}`$ |

It’s worthwhile to note that the total mass of any neighboring quark pair is precisely equal to the mass of the harmonic oscillator based on the quark with greater mass. Moreover, with respect to a given quark its two neighbors can be considered as an upward excitation (the quark with greater mass) and a downward excitation (the quark with lesser mass), with the equal electric charges of both excitations.

Finally, we have the chance to calculate actually the genuine quark masses and afterwards to verify the validity of the original propositions using ample experimental data.

## 4 Harmonic quark masses

### 4.1 Quark masses estimation

We can estimate the harmonic quark masses assuming that the charged meson mass is a total mass of the valence quarks composing it. Table 3 presents these estimations for each of the four stable mesons from the table 1.

Table 3. Harmonic quark masses, MeV (zero approximation)

| meson | $`d`$ | $`u`$ | $`s`$ | $`c`$ | $`b`$ | $`t`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`\pi ^\pm `$($`ud`$) | 29.9 | 109.6 | 401 | 1468 | 5373 | 19665 |
| $`K^\pm `$($`su`$) | 28.9 | 105.9 | 387.7 | 1419 | 5193 | 19005 |
| $`D_s^\pm `$($`cs`$) | 31.5 | 115.4 | 422.4 | 1546 | 5658 | 20708 |
| $`B_c^\pm `$($`bc`$) | 28.0 | 102.5 | 375.3 | 1373 | 5026 | 18396 |

The dispersion of these estimations does not exceed 12%. It is worth mentioning that $`B_c^\pm `$ meson mass value is measured with 8% inaccuracy and the estimation for meson $`D_s^\pm `$ supposedly produce an overestimation222The reasons for this overestimation are beyond the scope of the present paper and will be explained in further publications.. The pion- and kaon-estimated masses differ by 3.5%. This difference has a natural explanation if the sizes $`\pi ^\pm `$ and $`K^\pm `$ are approximately equal. The only difference of their structure is in $`d`$ and $`s`$-quarks and the localization of the $`d`$-quark, which is different by an order of mass magnitude from the $`s`$-quark will demand an additional kinetic energy. That’s why the pion-estimated quark masses will be greater.

### 4.2 Bounds of harmonic quark masses

Keeping in mind that quark masses in the model are rigidly bound, the heavier mesons we use for the quark mass determination, the more accurate it will be. Therefore we would start with the heavy charged two-quark meson with the mass measured exact enough, for example, $`b^\pm `$ with the mass $`5279.0\pm 0.5\text{ MeV}`$. The kaon-estimated (table 3) total mass of $`b`$ and $`u`$-quarks is $`5299\text{ MeV}`$. Provided that binding energy in the meson is not greater than the $`u`$-quark mass, in other words not greater than 50% of its mass per single quark (a binding energy in harmonic oscillator is only 36% from mass of the everyone quark) we have that $`b^\pm `$ mass is an upper bound for the harmonic $`b`$-quark mass. The lower bound can be determined as follows. Let’s assume that the value of binding energy in a meson is positive but less than a pion producing threshold (about $`140\text{ MeV}`$). Hence the lower bound of a harmonic $`b`$-quark mass is equal to $`5279140109=5030\text{ MeV}`$. Now we can calculate bounds for all other quarks using equation (8). The calculated data are given in the table 4.

Table 4. Bounds and mean values of the harmonic quark masses, MeV

|  | $`d`$ | $`u`$ | $`s`$ | $`c`$ | $`b`$ | $`t`$ |
| --- | --- | --- | --- | --- | --- | --- |
| lower bound | 28.0 | 102.6 | 375 | 1374 | 5030 | 18409 |
| upper bound | 29.4 | 107.7 | 394 | 1442 | 5279 | 19320 |
| mean value | 28.7 | 105.2 | 384 | 1408 | 5154 | 18864 |

We may see that the bound values for quark masses situate inside same data of tables 3 and the mean values are almost equal to kaon-calculated values. The fact that $`\pi ^\pm `$ mass exceeds the sum of the upper mass bounds of the $`u`$ and $`d`$ by $`2.5\text{ MeV}`$ means that pion binding energy is positive. It means that quarks are completely separated in $`\pi ^\pm `$ and there is an additional energy of several MeVs. What kind of energy it is? The $`u`$ and $`d`$-quarks of a pion charged alike and will repel each other. The repulsion force will be the more effective the less distance is between the quarks. On the other hand, the color field would act on the longer distances. Thus the additional kinetic energy would be maximal at the color and Coulomb forces equality and would transform to Coulomb energy on the small distances and to color-field energy on large distances. The value of this additional energy is naturally corresponds to the electromagnetic split in hadronic duplets and triplets. The complete separation of the quarks in $`\pi ^\pm `$ and the structural similarity of the stable charged mesons discovered earlier makes it possible to suppose that the quarks in other charged mesons (in particular $`K^\pm `$ and $`b^\pm `$) are also completely separated.

### 4.3 Electromagnetic split is taken into account

To calculate the mean value of electromagnetic split the author used well-known mass differences of the charged and neutral particles up to $`1350\text{ MeV}`$, i.e. of pions, kaons (with $`J=0,1`$) and basic baryon octet. The seven differences altogether. The mean split value is $`4.13\pm 1.47\text{ MeV}`$. Decreasing by this the mass values of $`\pi ^\pm `$, $`K^\pm `$ and $`b^\pm `$ and using equation 7 we obtain the quark masses to a first approximation (table 5). The inaccuracy value given depends only on electromagnetic split measurement inaccuracy and the meson mass. (For the $`b^\pm `$ its mass measurement inaccuracy was also considered).

Table 5. Harmonic quark masses (MeV)

determined with electromagnetic split consideration

(first approximation)

| meson | $`d`$ | $`u`$ | $`s`$ | $`c`$ | $`b`$ | $`t`$ | error,% |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $`\pi ^\pm `$ | 29.066 | 106.374 | 389.31 | 1424.8 | 5214.4 | 19083.7 | $`\pm `$1.08 |
| $`K^\pm `$ | 28.706 | 105.058 | 384.49 | 1407.2 | 5149.9 | 18847.5 | $`\pm `$0.30 |
| $`b^\pm `$ | 28.815 | 105.456 | 385.95 | 1412.5 | 5169.4 | 18919.0 | $`\pm `$0.030 |

There are two positive moments to emphasize here. First of all, pion- and kaon-calculated quark masses differ only by 1.2%, which is three times less than the same difference of masses to a zero approximation and, moreover, is in accordance with the inaccuracy calculated. And second of all, the most accurate values of the harmonic quark masses, calculated with the $`b^\pm `$ meson ($`0.03\%`$ inaccuracy), fall inside the range of pion- and kaon-calculated values. Hereafter the phrase “harmonic quark masses” would relate to these most accurate values.

## 5 First results

We will now some demonstrate the ability of harmonic model to solve real-world issues. Even the most superficial examination of the harmonic quark masses produce two impressive results:

* The total mass of $`u`$ and $`d`$ quarks is extremely close to the mass of neutral pion.
* $`u`$-quark mass is extraordinarily close to the muon mass.

Recollecting that neutral pion is the lightest existing hadron and the truly neutral particle it seems that neutral pion with great share of probability is a stationary harmonic oscillator based on quark-antiquark u-pair. The pion mass ($`134.97\text{ MeV}`$) differs from the oscillator mass ($`134.27\text{ MeV}`$ by 0.5% only. It’s likely that the relativistic correction would account for that.

What is to $`u`$-quark and muon, the difference between them is only $`0.2\text{ MeV}`$, with the $`u`$-quark mass value inaccuracy of $`0.03\text{ MeV}`$ which is seven times less than the difference. The obvious conclusion is that muon is a successful attempt of Nature to explicitly fix the single $`u`$-quark mass state as a lepton suppressing color and fractional charge. As muon charge is greater by $`\frac{1}{3}e`$ than the $`u`$-quark charge, the Coulomb energy of muon should also be greater. If electron mass is practically Coulomb energy then the charge in $`\frac{1}{3}e`$ add to the quark mass $`1`$/3 from the electron mass. Hence, the difference in $`0.2\text{ MeV}`$ has an easy explanation. Thus the mystery of muon mass is naturally solved. It is strong result and good argument for harmonic model. Furthermore, the $`MQ`$ factor works well not only with quarks but with pions too, thereby supporting the idea of similar structural formation of the mesons with different flavors:

* $`m_{\pi ^0}MQ=493.99\text{ MeV}`$, which is only $`0.34\text{ MeV}`$ greater than the mass of the first meson containing strangeness ($`K^\pm `$)
* $`m_{\pi ^\pm }MQ^2=1869.4(\mathrm{MeV})`$, which is precisely equal to the mass value of the first meson containing charm ($`D^\pm `$)

In the last case, the additional energy of the pion quarks is also transformed by the factor $`MQ^2`$, which leads to the effectively increase of the $`D^\pm `$ meson mass. This could be one of the above-mentioned reasons of the harmonic quark masses overestimation when calculated with $`D^\pm `$ and $`D_s^\pm `$ mesons.

Acknowledgement

I thank V. Teplov for the text correction and help with Latex.


format: text

text:
# 1 Introduction

## 1 Introduction

Scalar field cosmological models are of great importance in modern cosmology. The dark energy is attributed to the dynamics of a scalar field, which convincingly realizes the goal of explaining current accelerating expansion of universe generically using only attractor solutions . Furthermore a scalar field can drive an accelerated expansion and thus provides possible models for cosmological inflation in the early universe . In particular, there have been a number of studies of spatially homogeneous scalar field cosmological models with an exponential potential. They are already known to have interesting properties; for example, if one has a universe containing a perfect fluid and such a scalar field, then for a wide range of parameters the scalar field mimics the perfect fluid, adopting its equation of state . These scaling solutions are attractors at late times . The inflation models and other cosmological consequences of multiple scalar fields have also been considered .

The scale-invariant form makes the exponential potential particularly simple to study analytically. There are well-known exact solutions corresponding to power-law solutions for the cosmological scale factor $`at^p`$ in a spatially flat Friedmann-Robertson-Walker (FRW) model . More generally the coupled Einstein-Klein-Gordon equations for a single field can be reduced to a one-dimensional system which makes it particularly suitable for a qualitative analysis . Recently, adopting a system of dimensionless dynamical variables , the cosmological scaling solutions with positive or negative exponentials have been studied . In general there are many scalar fields with exponential potentials in supergravity, superstring and the generalized Einstein theories, thus multiple potentials may be more interesting. In the previous paper , We studied the stability of cosmological scaling solutions in an expanding universe model with multiple scalar fields with positive or negative exponential potentials. A phase-space analysis of the spatially flat FRW models shows that there exist cosmological scaling solutions which are the unique late-time attractors and successful inflationary solutions driven by multiple scalar fields with a wide range of each potential slope parameter $`\lambda `$. It is assumed that there is no direct coupling between potentials. Multiple cross-coupling exponential potentials arise in many occasions, for instance, from compactifications of vacuum Einstein gravity on product spaces . Indeed they are a natural outcome of the compactification of higher dimensional theories down to 3+1 dimensions. With this in mind it is worth investigating such potential in a bit more detail.

In this paper, we first study a system of dimensionless dynamical variables of multiple scalar fields with a positive or negative cross-coupling exponential potential. We obtain the scaling solutions and analyze their stability. There still exist cosmological scaling solutions which are the unique late-time attractors. In this model we then introduce a barotropic fluid to the system. We discuss the physical consequences of these results.

## 2 Cross-coupling Exponential Potential

We consider $`n`$ scalar fields $`\varphi _i`$ with a cross-coupling potential

$$V=V_0\mathrm{exp}\left(\underset{i=1}{\overset{n}{}}\lambda _i\kappa \varphi _i\right),$$
(1)
where $`\kappa ^28\pi G_N`$ is the gravitational coupling and $`\lambda _i`$ are dimensionless constants characterising the slope of the potential. Further we assume all $`\lambda _i0`$ since we can make them positive through $`\varphi _i\varphi _i`$ if some of them are negative. The evolution equation of each scalar field for a spatially flat FRW model with Hubble parameter $`H`$ is

$$\ddot{\varphi _i}+3H\dot{\varphi _i}\lambda _i\kappa V=0,$$
(2)
subject to the Friedmann constraint

$$H^2=\frac{\kappa ^2}{3}\left(\underset{i=1}{\overset{n}{}}\frac{1}{2}\dot{\varphi }_i^2+V\right).$$
(3)
Defining $`(n+1)`$ dimensionless variables

$$x_i=\frac{\kappa \dot{\varphi _i}}{\sqrt{6}H},y=\frac{\kappa \sqrt{|V|}}{\sqrt{3}H},$$
(4)
the evolution equations (2) can be written as an autonomous system:

$`x_i^{}`$ $`=`$ $`3x_i\left(1{\displaystyle \underset{j=1}{\overset{n}{}}}x_j^2\right)\pm \lambda _i\sqrt{{\displaystyle \frac{3}{2}}}y^2,`$ (5)
$`y^{}`$ $`=`$ $`y{\displaystyle \underset{j=1}{\overset{n}{}}}\left(3x_j^2\lambda _j\sqrt{{\displaystyle \frac{3}{2}}}x_j\right),`$ (6)
where a prime denotes a derivative with respect to the logarithm of the scalar factor, $`N\mathrm{ln}a`$, and the constraint equation (3) becomes

$$\underset{i=1}{\overset{n}{}}x_i^2\pm y^2=1.$$
(7)
Throughout this paper we will use upper/lower signs to denote the two distinct cases of $`\pm V_0>0`$. $`x_i^2`$ measures the contribution to the expansion due to the field’s kinetic energy density, while $`\pm y^2`$ represents the contribution of the potential energy. We will restrict our discussion of the existence and stability of critical points to expanding universes with $`H>0`$, i.e., $`y0`$. Critical points correspond to fixed points where $`x_i^{}=0`$ and $`y^{}=0`$, and there are self-similar solutions with

$$\frac{\dot{H}}{H^2}=3\underset{i=1}{\overset{n}{}}x_i^2.$$
(8)
This corresponds to an expanding universe with a scale factor $`a(t)`$ given by $`at^p`$, where

$$p=\frac{1}{3_{i=1}^nx_i^2}.$$
(9)
The system (5) and (6) has at most one $`n`$-dimensional sphere $`S`$ embedded in $`(n+1)`$-dimensional phase-space corresponding to kinetic-dominated solutions, and a fixed point $`A`$, which is a kinetic-potential-scaling solution listed in Table 1.

In order to study the stability of the critical points, using the Friedmann constraint equation (7) we first reduce Eqs.(5) and (6) to $`n`$ independent equations

$$x_i^{}=\left(\lambda _i\sqrt{\frac{3}{2}}3x_i\right)\left(1\underset{j=1}{\overset{n}{}}x_j^2\right).$$
(10)
Substituting linear perturbations $`x_ix_i+\delta x_i`$ about the critical points into Eqs.(10), to first-order in the perturbations, gives equations of motion

$$\delta x_i^{}=2\left(\lambda _i\sqrt{\frac{3}{2}}3x_i\right)\underset{j=1}{\overset{n}{}}(x_j\delta x_j)3\left(1\underset{j=1}{\overset{n}{}}x_j^2\right)\delta x_i,$$
(11)
which yield $`n`$ eigenvalues $`m_i`$. Stability requires the real part of all eigenvalues being negative.

$`S`$: $`_{i=1}^nx_i^2=1`$, $`y=0`$. These kinetic-dominated solutions always exist for any form of the potential, which are equivalent to stiff-fluid dominated evolution with $`at^{1/3}`$ irrespective of the nature of the potential. Then Eqs.(11) become

$$\delta x_i^{}=2\left(\lambda _i\sqrt{\frac{3}{2}}3x_i\right)\underset{j=1}{\overset{n}{}}(x_j\delta x_j),$$

which yield $`n`$ eigenvalues: one of them, say $`m_1`$, does not vanish, $`m_1=\sqrt{6}(_{i=1}^n\lambda _ix_i\sqrt{6})`$; the remains of them vanish. Thus the solutions are marginally stable for $`_{i=1}^n(\lambda _ix_i)>\sqrt{6}`$. For the special case $`\lambda _i=\lambda `$, using the constraint equation (7) we find $`\sqrt{6}/(n\lambda )<_{i=1}^nx_i/n(_{i=1}^nx_i^2/n)^{1/2}=1/\sqrt{n}`$. That is, if each scalar field has an identical-slope potential, there exist stable points only for $`\lambda ^2>6/n`$.

$`A`$: $`x_i=\frac{\lambda _i}{\sqrt{6}}`$, $`y=\sqrt{\pm (1\frac{1}{6}_{i=1}^n\lambda _i^2)}`$. The potential-kinetic-scaling solution exists for sufficiently flat $`_{i=1}^n\lambda _i^2<6`$ positive potentials or steep $`_{i=1}^n\lambda _i^2>6`$ negative potentials. The power-law exponent, $`p=\frac{2}{_{i=1}^n\lambda _i^2}`$, depends on parameter $`\lambda _i`$. From Eqs.(11) we find the eigenvalues

$$m_i=\frac{1}{2}\left(6\underset{j=1}{\overset{n}{}}\lambda _j^2\right).$$

Thus the scaling solution is always stable when this point exists for a positive potential, which corresponds to the power-law inflation in an expanding universe when $`_{i=1}^n\lambda _i^2<2`$. However, this solution is unstable for a negative potential.

The different regions of $`\lambda _i`$ parameter space lead to different qualitative evolution. As an example we consider the cosmologies containing $`n`$ scalar fields with the cross-coupling potential $`\lambda _i=\lambda `$. For the sufficiently flat ($`\lambda ^2<6/n`$) positive potential, these kinetic-dominated solutions are unstable and the kinetic-potential-scaling solution is the stable late-time attractor. Hence generic solutions start in the former and approach the later at late times. For the steep ($`\lambda ^2>6/n`$) positive potential, there exists a stable kinetic-dominated regime, in which each points are the late-time attractors. Hence generic solutions start in kinetic-dominated solution and approach the stable regime. For the flat sufficiently ($`\lambda ^2<6/n`$) negative potential, only these kinetic-dominated solutions exist which are unstable scaling solutions. For the steep ($`\lambda ^2>6/n`$) negative potential, the kinetic-potential-scaling solution is unstable and there exists a stable kinetic-dominated regime. Hence generic solutions start in a kinetic-dominated regime or the kinetic-potential-scaling solution and approach the stable kinetic-dominated regime at late times.

## 3 Plus a Barotropic Fluid

We now consider multiple scalar fields with the cross-coupling potential (1) evolving in a spatially flat FRW universe containing a fluid with barotropic equation of state $`P_\gamma =(\gamma 1)\rho _\gamma `$, where $`\gamma `$ is a constant, $`0<\gamma 2`$, such as radiation ($`\gamma =4/3`$) or dust ($`\gamma =1`$). The evolution equation for the barotropic fluid is

$$\dot{\rho _\gamma }=3H(\rho _\gamma +P_\gamma ),$$
(12)
subject to the Fridemann constraint

$$H^2=\frac{\kappa ^2}{3}\left(\underset{i=1}{\overset{n}{}}\frac{1}{2}\dot{\varphi }_i^2+V+\rho _\gamma \right).$$
(13)
We define another dimensionless variable $`z\frac{\kappa \sqrt{\rho _\gamma }}{\sqrt{3}H}`$. The evolution equations (2) and (12) can then be written as an autonomous system:

$`x_i^{}`$ $`=`$ $`3x_i\left(1{\displaystyle \underset{j=1}{\overset{n}{}}}x_j^2{\displaystyle \frac{\gamma }{2}}z^2\right)\pm \lambda _i\sqrt{{\displaystyle \frac{3}{2}}}y^2,`$ (14)
$`y^{}`$ $`=`$ $`y\left(3{\displaystyle \underset{i=1}{\overset{n}{}}}x_i^2+{\displaystyle \frac{3\gamma }{2}}z^2\sqrt{{\displaystyle \frac{3}{2}}}{\displaystyle \underset{i=1}{\overset{n}{}}}\lambda _ix_i\right),`$ (15)
$`z^{}`$ $`=`$ $`{\displaystyle \frac{3}{2}}z\left(\gamma +2{\displaystyle \underset{i=1}{\overset{n}{}}}x_i^2+\gamma z^2\right),`$ (16)
and the constraint equation becomes

$$\underset{i=1}{\overset{n}{}}x_i^2\pm y^2+z^2=1.$$
(17)
Critical points correspond to fixed points where $`x_i^{}=0`$, $`y^{}=0`$ and $`z^{}=0`$, and there are self-similar solutions with

$$\frac{\dot{H}}{H^2}=3\underset{i=1}{\overset{n}{}}x_i^2\frac{3\gamma }{2}z^2.$$
(18)
This corresponds to an expanding universe with a scale factor $`a(t)`$ given by $`at^p`$, where

$$p=\frac{2}{6_{i=1}^nx_i^2+3\gamma z^2}.$$
(19)
The system (14)-(16) has at most one $`n`$-dimensional sphere $`S`$ embedded in $`(n+2)`$-dimensional phase-space corresponding to kinetic-dominated solutions, a fixed point $`A`$ which is a kinetic-potential-scaling solution, a fixed point $`B`$ which is a fluid-dominated solution, and a fixed point $`C`$ which is a fluid-potential-kinetic-scaling solution listed in Table 2.

$`S`$: $`_{i=1}^nx_i^2=1`$, $`y=0`$, $`z=0`$. These kinetic-dominated solutions always exist for any form of the potential, which are equivalent to stiff-fluid dominated evolution with $`at^{1/3}`$ irrespective of the nature of the potential. The linearization of system (14)-(16) about these fixed points yields

$`\delta x_i^{}`$ $`=`$ $`2\left(\lambda _i\sqrt{{\displaystyle \frac{3}{2}}}3x_i\right){\displaystyle \underset{j=1}{\overset{n}{}}}(x_j\delta x_j),`$
$`\delta z^{}`$ $`=`$ $`{\displaystyle \frac{3}{2}}(2\gamma )\delta z,`$
which indicate that the solutions are marginally stable for $`_{i=1}^n(\lambda _ix_i)>\sqrt{6}`$ and a stiff fluid $`\gamma =2`$.

$`A`$: $`x_i=\frac{\lambda _i}{\sqrt{6}}`$, $`y=\sqrt{\pm (1\frac{1}{6}_{i=1}^n\lambda _i^2)}`$, $`z=0`$. The potential-kinetic-scaling solution exists for sufficiently flat $`_{i=1}^n\lambda _i^2<6`$ positive potentials or steep $`_{i=1}^n\lambda _i^2>6`$ negative potentials. The power-law exponent, $`p=\frac{2}{_{i=1}^n\lambda _i^2}`$, depends on the slope of the potential. The linearization of system (14)-(16) about this critical point yields $`(n+1)`$ eigenvalues

$`m_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(6{\displaystyle \underset{j=1}{\overset{n}{}}}\lambda _j^2\right),`$
$`m_z`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(3\gamma {\displaystyle \underset{j=1}{\overset{n}{}}}\lambda _j^2\right).`$
Thus the scaling solution is stable for a positive potential with $`_{i=1}^n\lambda _i^2<3\gamma `$, which corresponds to the power-law inflation in an expanding universe when $`_{i=1}^n\lambda _i^2<2`$.

$`B`$: $`x_i=0`$, $`y=0`$, $`z=1`$. The fluid-dominated solution exists for any form of the potential, corresponding to a power-law solution with $`p=2/3\gamma `$.

$`\delta x_i^{}`$ $`=`$ $`3\delta x_i+(3\gamma \sqrt{6}\lambda _i)\delta z,`$
$`\delta z^{}`$ $`=`$ $`3\gamma \delta z,`$
which indicate that the solution is never stable.

$`C`$: $`x_i=\sqrt{\frac{3}{2}}\frac{\gamma \lambda _i}{_{j=1}^n\lambda _j^2}`$, $`y=\sqrt{\frac{3}{2}\frac{(2\gamma )\gamma }{_{i=1}^n\lambda _i^2}}`$, $`z=\sqrt{1\frac{3\gamma }{_{i=1}^n\lambda _i^2}}`$. The fluid-potential-kinetic-scaling solution exists for a positive potential with $`_{i=1}^n\lambda _i^2>3\gamma `$. The power-law exponent, $`p=2/3\gamma `$, is identical to that of the fluid-dominated solution, depends only on the barotropic index $`\gamma `$ and is independent of the slope $`\lambda _i`$ of the potential. The linearization of system (14)-(16) about the fixed point yields

$`\delta x_i^{}`$ $`=`$ $`3(2\gamma )x_i{\displaystyle \underset{j=1}{\overset{n}{}}}(x_j\delta x_j)3\left((1{\displaystyle \frac{\gamma }{2}})(1{\displaystyle \underset{j=1}{\overset{n}{}}}x_j^2)+{\displaystyle \frac{\gamma }{2}}y^2\right)\delta x_i`$
$`+(\sqrt{6}\lambda _i3\gamma x_i)y\delta y,`$
$`\delta y^{}`$ $`=`$ $`3(2\gamma )y{\displaystyle \underset{j=1}{\overset{n}{}}}(x_j\delta x_j)\sqrt{{\displaystyle \frac{3}{2}}}y{\displaystyle \underset{j=1}{\overset{n}{}}}(\lambda _j\delta x_j)`$
$`+\left({\displaystyle \frac{3}{2}}(2\gamma ){\displaystyle \underset{j=1}{\overset{n}{}}}x_j^2\sqrt{{\displaystyle \frac{3}{2}}}{\displaystyle \underset{j=1}{\overset{n}{}}}(\lambda _jx_j)+{\displaystyle \frac{3\gamma }{2}}{\displaystyle \frac{9\gamma }{2}}y^2\right)\delta y,`$
which yield (n+1) eigenvalues

$`m_1`$ $`=`$ $`{\displaystyle \frac{3(2\gamma )}{4}}\left(1+\sqrt{1{\displaystyle \frac{8\gamma (_{i=1}^n\lambda _i^23\gamma )}{_{i=1}^n\lambda _i^2(2\gamma )}}}\right),`$
$`m_2`$ $`=`$ $`{\displaystyle \frac{3(2\gamma )}{4}}\left(1\sqrt{1{\displaystyle \frac{8\gamma (_{i=1}^n\lambda _i^23\gamma )}{_{i=1}^n\lambda _i^2(2\gamma )}}}\right),`$
$`m_3`$ $`=`$ $`\mathrm{}=m_{n+1}={\displaystyle \frac{3}{2}}(2\gamma ).`$
Thus the scaling solution is stable for a positive potential with $`_{i=1}^n\lambda _i^2>3\gamma `$.

The different regions in the ($`\gamma ,\lambda _i`$) parameter space lead to different qualitative evolution. For the sufficiently flat ($`_{i=1}^n\lambda _i^2<3\gamma `$) positive potential, $`S`$, $`A`$ and $`B`$ exist. Point $`A`$ is the stable late-time attractor. Hence generic solutions begin in a kinetic-dominated regime or at the fluid-dominated solution and approach the kinetic-potential-scaling solution at late times. For the intermediate ($`3\gamma <_{i=1}^n\lambda _i^2<6`$) positive potential, all critical points exist. Point $`C`$ is the stable late-time attractor. Hence generic solutions start in a kinetic-dominated regime, at the kinetic-potential-scaling solution or at the fluid-dominated solution and approach the stable fluid-kinetic-potential-scaling solution. For the steep ($`_{i=1}^n\lambda _i^2>6`$) positive potential, $`S`$, $`B`$ and $`C`$ exist. Point $`C`$ is the stable late-time attractor. Hence generic solutions start in a kinetic-dominated regime or at the fluid-dominated solution and approach the stable fluid-kinetic-potential-scaling solution. For the sufficiently flat ($`_{i=1}^n\lambda _i^2<3\gamma `$) negative potential, the kinetic-dominated solution $`S`$ and the fluid-dominated solution $`B`$ exist, which are unstable. For the intermediate ($`3\gamma <_{i=1}^n\lambda _i^2<6`$) negative potential, the kinetic-dominated solution $`S`$ and the fluid-dominated solution $`B`$ exist, which are unstable. For the steep ($`_{i=1}^n\lambda _i^2>6`$) negative potential, $`S`$, $`A`$ and $`B`$ exist. Point $`A`$ is the stable late-time attractor. Hence generic solutions start in a kinetic-dominated regime or at the fluid-dominated solution and approach the stable kinetic-potential-scaling solution at late times.

## 4 Conclusions and Discussions

We have presented a phase-space analysis of the evolution for a spatially flat FRW universe containing $`n`$ scalar fields with a positive or negative cross-coupling exponential potential. In particular, for the $`\lambda _i=\lambda `$ case, we find that in the expanding universe model with a sufficiently flat ($`\lambda ^2<6/n`$) positive cross-coupling potential the only power-law kinetic-potential-scaling solution is the late-time attractor. It is more difficult to obtain assisted inflation in such models since the fields with cross-coupling exponential potential tend to conspire to act against one another rather than assist each other. However, steep ($`\lambda ^2>6/n`$) negative cross-coupling potential has kinetic-dominated solutions with $`at^{1/3}`$, some of which are the late-time attractors. It can be known that the kinetic energy of each field tends to be equal via their effect on the expansion at late times.

Then we have extended the phase-space analysis of the evolution to a realistic universe model with a barotropic fluid plus $`n`$ scalar fields with a positive or negative cross-coupling exponential potential. We have shown that for the sufficiently flat ($`_{i=1}^n\lambda _i^2<3\gamma `$) positive cross-coupling potential, the kinetic-potential-scaling solution is the stable late-time attractor. The energy density of the scalar fields dominates at late times. Moreover, for the steep ($`_{i=1}^n\lambda _i^2>6`$) positive cross-coupling potential, the fluid-kinetic-potential-scaling solution is the stable late-time attractor. However, a negative cross-coupling potential has no stable scaling solutions.

## Acknowledgements

This project was in part supported by NNSFC under Grant Nos. 10175070 and 10047004 as well as by NKBRSF G19990754.


text:
# Bertini intra-nuclear cascade implementation in Geant4

## I INTRODUCTION

The intra-nuclear cascade model (INC) was first proposed by Serber in 1947 serber47 . He noticed that, in particle-nuclear collisions the deBroglie wave-length of the incident particle is comparable to or shorter than the average intra-nucleon distance. Hence, the justification for describing the interactions in terms of particle-particle collisions.

The INC has been successfully used in the Monte Carlo simulations at intermediate energy region since Goldberger made first calculations by hand in 1947 goldberger48 . First computer simulations were done by Metropolis et al. in 1958 metropolis58 . Standard methods in INC implementations were formed when Bertini published his results in 1968 bertini68 . An important addition was exciton model introduced by Griffin in 1966 griffin66 .

Our presentation describes the implementation of a Bertini INC model in the Geant4 hadronic physics frameworks geant4collaboration03 . Geant4geant4  is a Monte Carlo particle detector simulation toolkit, having applications also in the medical field and space sciences. Geant4 provides a flexible framework for modular implementation of various kinds of hadronic interactions. Geant4 exploits advanced Software Engineering techniques and Object Oriented technology to achieve the transparency of the physics implementation and to this way provide the possibility of validating the physics results more easily.

At level-2, the hadronic model frameworks are based on concepts of physics processes, cross-section and models for final state generation. While the process is a general concept, models are allowed to have restrictions in process type, material, element and energy range. Several models can be utilized by one process class; for instance, a process class for inelastic scattering can use distinct models for different energies.

Process classes utilize model classes to determine the secondaries produced in the interaction and to calculate the momenta of the particles. Here we preset a collection of such models providing medium-energy intra-nuclear cascade treatment.

## II GEANT4 CASCADE MODEL

In inelastic particle-nucleus collisions a fast phase ($`10^{23}10^{22}s`$) of INC results in a highly exited nucleus, and is followed by fission and pre-equilibrium emission. A slower ($`10^{18}10^{16}s`$) compound nucleus phase follows with evaporation. A Boltzman equation must be solved to treat the physical process of the collision in detail.

The intra-nuclear cascade model developed by Bertini bertini68 ; bertini69 ; bertini71  solves the Boltzmann equation on the average. The model has been implemented in several codes such as HETC alsmiller90 . Our model is based on re-engineering of INUCL code titarenko99a . Models included are a Bertini intra-nuclear cascade model with exitons, a pre-equilibrium model, a simple nuclear explosion model, a fission model, and a evaporation model.

The nuclear model consist of a three-region approximation to the continuously changing density distribution of nuclear matter within nuclei. Relativistic kinematics is applied throughout the cascade. The cascade is stopped when all the particles which can escape the nucleus, have done so. Then conformity with the energy conservation law is checked.

### II.1 Model limits

Particles treated are nucleons, pions, photons, and nuclear isotopes. A bullet particle can be a proton, neutron or pion. The range of targets allowed is arbitrary.

The necessary condition of validity of the INC model is $`\lambda _B/v<<\tau _c<<\mathrm{\Delta }t`$, where $`\lambda _B`$ is the deBroglie wavelength of the nucleons, $`v`$ is the average relative nucleon-nucleon velocity and $`\mathrm{\Delta }t`$ is the time interval between collisions. The physical foundation becomes approximate at energies less than about $`200MeV`$, and there needs to be supplemented with a pre-equilibrium model. Also, at energies higher than 5-10 GeV the INC picture breaks down. Our implementation in Geant4 has been tested with bullet kinetic energy between 100 MeV and 3 GeV.

### II.2 Intra-nuclear cascade model

The basic steps of the INC model are summarized below:

1. The spatial point, where the incident particle enters, is selected uniformly over the projected area of the nucleus.
2. Total, free particle-particle cross-sections and region-dependent nucleon densities are used to select the path length for the projectile particle.
3. The momentum of a struck nucleon, the type of reaction, and the four momentum of the reaction products are determined.
4. the exciton model is updated as the cascade proceeds.
5. If Pauli’s exclusion principle allows and $`E_{particle}>E_{cutoff}`$ = 2 MeV, step (2) is performed to transport the products.

After INC, the residual excitation energy of the resulting nucleus is used as input for a non-equilibrium model. The schematic presentation of Bertini’s INC is shown in Fig. 1

### II.3 The nuclear model

Some of the basic features of the nuclear model are the following:

* The nucleons are assumed to to have a Fermi gas momentum distribution. The Fermi energy is calculated in a local density approximation i.e. it is made radius dependent with Fermi momentum $`p_F(r)=(\frac{3\pi ^2\rho (r)}{2})^{\frac{1}{3}}`$.
* The nucleon’s binding energies (BE) are calculated using a mass formula. A parametrization of the nuclear binding energy uses a combination of the Kummel mass formula, and experimental data. Also, a asymptotic high temperature mass formula is used, for cases where it is not possible to use experimental data.

#### II.3.1 Initialization

The initialization phase fixes the nucleus radius and momentum according to the Fermi gas model.

If the target is Hydrogen (A = 1), a direct particle-particle collision is performed, and no nuclear modeling is used.

If $`1<A<4`$, a nuclei model consisting one layer with radius of 8.0 fm is created.

For $`4<A<11`$, a nuclei model is composed of three concentric spheres $`i=\{1,2,3\}`$ with radii

$$r_i(\alpha _i)=\sqrt{C_1^2(1\frac{1}{A})+6.4}\sqrt{log(\alpha _i)}$$

where $`\alpha _i=\{0.01,0.3,0.7\}`$ and $`C_1=3.3836A^{1/3}`$

If $`A>11`$, nuclei are modeled with three concentric spheres as well. The sphere radii are then defined as:

$$r_i(\alpha _i)=C_2\mathrm{log}(\frac{1+e^{\frac{C_1}{C_2}}}{\alpha _i}1)+C_1$$

where $`C_2=1.7234`$.

The potential energy for nucleon N is

$$V_N=\frac{p_F^2}{2m_N}+BE_N(A,Z)$$

where $`p_F`$ is a Fermi momentum and BE a binding energy.

Impulse distribution in each region follows Fermi distribution with zero temperature.

$$f(p)=cp^2$$
(1)
where

$$_0^{p_F}f(p)𝑑p=n_porn_n.$$
(2)
Here $`n_p`$ and $`n_n`$ are the number of protons or neutrons in a region and $`p_F`$ is momentum corresponding the Fermi energy

$$E_F=\frac{p_F^2}{2m_N}=\frac{\mathrm{}^2}{2m_N}(\frac{3\pi ^2}{v})^{\frac{2}{3}},$$
(3)
which depends on the density $`n/v`$ of particles, and which is different for each particle and each region.

#### II.3.2 Pauli’s exclusion principle

Pauli’s exclusion principle forbids interactions where the products would be in occupied states. Following the assumption of a completely degenerate Fermi gas, the levels are filled from the lowest level. The minimum energy allowed for a collision product corresponds to the lowest unfilled level of system, which is the Fermi energy in the region. So in practice, Pauli exclusion principle is taken into account by accepting only secondary nucleons which have $`E_N>E_F`$.

#### II.3.3 Cross-sections and kinematics

Path lengths of nucleons in the nucleus are sampled according to the local density and free nucleon-nucleon cross-sections. Angles after collisions are sampled from experimental differential cross-sections. Tabulated total reaction cross-sections are calculated by Letaw’s formulation letaw83 ; letaw93 ; pearlstein89 . For nucleon-nucleon cross-sections, parameterizations based on the experimental energy and isospin dependent data. The parameterization described in barashenkov72  is used.

For pion the INC cross-sections are provided to treat elastic collisions, and inelastic channels: $`\pi ^{}`$n $``$ $`\pi ^0`$n, $`\pi ^0`$p $``$ $`\pi ^+`$n and $`\pi ^0`$n $``$ $`\pi ^{}`$p. Multiple particle production is also implemented.

The S-wave pion absorption channels $`\pi ^+`$nn $``$ pn, $`\pi ^+`$pn $``$ pp, $`\pi ^0`$nn $``$ X , $`\pi ^0`$pn $``$ pn, $`\pi ^0`$pp $``$ pp, $`\pi ^{}`$nn $``$ X , $`\pi ^{}`$pn $``$ nn , and $`\pi ^{}`$pp $``$ pn are implemented.

### II.4 Pre-equilibrium model

The Geant4 cascade model implements the exciton model proposed by Griffin griffin66 . In this model nuclear states are characterized by the number of exited particles and holes (the exitons). INC collisions give rise to a sequence of states characterized by increasing exciton number, eventually leading to a equilibrated nucleus. For practical implementation of the exciton model we use level density parameters from ribansky73  and the matrix elements from kalbach78 .

In the exciton model the possible selection rules for particle-hole configurations in the course of the cascade are: $`\mathrm{\Delta }p=0,\pm 1`$ $`\mathrm{\Delta }h=0,\pm 1`$ $`\mathrm{\Delta }n=0,\pm 2`$, where p is the number of particle, h is number of holes and n = p + h is the number of exitons.

The cascade pre-equilibrium model uses target excitation data, and exciton configurations for neutron and proton to produce the non-equilibrium evaporation. The angular distribution is isotropic in the frame of rest of the exciton system.

The parameterizations of the level density used, are tabulated both with their A and Z dependence and including a high temperature behavior. The nuclear binding energy is using a smooth liquid high energy formula.

### II.5 Break-up models

Fermi break-up is allowed only in some extreme cases, i.e. for light nuclei ($`A<12`$ and $`3(AZ)<Z<6`$ ) and if $`E_{excitation}>3E_{binding}`$ . A simple explosion model decays the nucleus into neutrons and protons and decreases exotic evaporation processes.

The fission model is a phenomenological model using potential minimization. Binding energy parametrization is used and some features of the fission statistical model are incorporated as infong69 .

### II.6 Evaporation model

The statistical theory for particle emission from exited nuclei remaining after INC was originally developed by Weisskopf weisskopf37 . This model assumes complete energy equilibration before particle emission, and re-equilibration of excitation energies between successive evaporation emissions. As a result, the angular distribution of emitted particles is isotropic.

The Geant4 evaporation model for cascade implementation adapts the widely used computational method developed by Dostrowski dostrovsky59 ; dostrovsky60 . The emission of particles is computed until the excitation energy fall below a cutoff. If a light nucleus is highly exited, a Fermi break-up model is executed. In addition, fission is performed when the fission channel is open. The main chain of evaporation is followed until $`E_{excitation}`$ falls below Ecutoff = 0.1 MeV. The evaporation model ends with a $`\gamma `$ emission chain, which is followed until Eexcitation $`<`$ E$`{}_{cutoff}{}^{}{}_{}{}^{\gamma }`$ = 10-15 MeV.

## III IMPLEMENTATION

The Bertini cascade model is implemented in the Geant4 hadronic physics framework. The Bertini cascade source code module is located in directory Geant4/source/processes/hadronic/models/cascade/cascade. All the models are used collectively through interface method ApplyYourself defined in a class G4CascadeInterface. A Geant4 track (G4Track) and a nucleus (G4Nucleus) are given as a parameters.

Currently implementation is not optimized for speed. The typical timing behavior for different target nuclons and bullet energies is presented in Fig. 2.

## IV RESULTS

We have tested the physics performance of Bertini cascade models in the first Geant4 5.0 implementation for proton induced reactions for energies from 100 MeV to 3 GeV. Detailed comparisons with experimental data has been made in energy range 160 – 800 MeV.

The double differential cross-sections for neutrons at angles of $`7.5^{}`$, $`30^{}`$, $`60^{}`$, and $`150^{}`$ are presented in Figs. 3-6. Proton incident particle at 256 MeV energy is hitting an iron target and all Bertini sub-models are used. Data were taken from 256iron .

Figures 7 and 8 show double differential cross-section of neutrons resulting from 597 MeV protons impinging on lead. Data were taken from 597lead

Finally, an example of pion production physics performance. Double differential cross-sections for $`\pi ^+`$ are given for angels from $`22.5^{}`$ to $`135^{}`$ in Figs. 9-12.

The agreement with experimental data is found to be from moderate to relatively good. Data fo the 585 MeV pion comparisons were taken from 585lead .

## V CONCLUSION

We have released a Bertini INC model in Geant4. Exitons, pre-equilibrium, nucleus explosion, fission, and evaporation are modeled. Particles treated are $`\gamma `$, $`\pi `$, n, p, and nuclear isotopes, and incident nucleons and pions can be treated. We have tested the code in a energy range 60 MeV - 3 GeV and found a reasonable agreement with experimental data. Tuning of performance and models are planned for future Geant4 releases. We also plan to extend the upper energy limit of the models.


text:
# A Comparison of Measured Crab and Vela Glitch Healing Parameters with Predictions of Neutron Star Models

## 1 Introduction

Soon after the discovery of radio pulsars (Hewish et al., 1968), they were identified as rotating, highly magnetized neutron stars (NSs) (Gold, 1968, 1969; Ostriker & Gunn, 1969). It was observationally established that pulsar periods slowly and monotonically increase with time as a result of magnetic braking and rotational energy loss. Early in 1969, it was noticed that the Vela pulsar suddenly increased its angular velocity (Reichley & Downs, 1969; Radhakrishnan & Manchester, 1969). These distinct and sudden increases in rotational frequency, known as glitches, were later regularly observed in the Vela and Crab pulsars and, more infrequently, in other pulsars (e.g., Shemar & Lyne 1996; Lyne, Shemar, & Smith 2000). The Crab (PSR B0531+21) and Vela (PSR B0833$``$45) pulsars are the two best-studied glitching pulsars since they have each been observed to glitch a number of times and are bright and easy to monitor. These pulsars have the most extensive and accurate glitch data published and currently provide the best test of the physical mechanisms by which pulsar glitches occur.

There are currently two competing models that are well-accepted to explain how pulsars glitch. In simplest terms, both models treat the NS as a two-component body, with a superfluid interior core surrounded by a rigid external crust (e.g., Ruderman 1972). These two components are dynamically weakly coupled through the magnetic field.

In the starquake glitch model (Ruderman, 1976; Baym & Pines, 1971; Alpar et al., 1996), a slight equatorial oblateness in the crust can be formed if the NS is born rapidly spinning. As the pulsar slows down via magnetic braking, the deformation cannot be supported through centrifugal pressure from the core, and the crust can suddenly crack under gravity. The subsequent reduction in oblateness reduces the moment of inertia, resulting in a sudden increase in the rotational frequency, seen as a glitch.

Alternatively, in the vortex unpinning model of glitches (Anderson & Itoh, 1975; Alpar et al., 1984, 1993), angular momentum is stored in vortices of superfluid which migrate outward from the core as the star slows down. These vortices can become pinned to nuclei in the transition region between the core and outer crust, thereby preventing angular momentum transfer to the crust. A differential rotation develops between the core and crust until a catastrophic unpinning of the vortices occurs. Angular momentum is suddenly transferred to the crust, spinning it up. Since the crust is tightly coupled to the external magnetic field, this sudden spin-up in the rotation period is observed as a glitch.

Within the starquake model, the time evolution of a pulsar’s rotational frequency after the occurrence of a glitch (Wang et al., 2000) can be fit by the equation:

$$\mathrm{\Omega }(t)=\mathrm{\Omega }_0(t)+\mathrm{\Delta }\mathrm{\Omega }_{t=0}\left[1Q(1e^{t/\tau })\right]+\mathrm{\Delta }\dot{\mathrm{\Omega }}_pt.$$
(1)
Here $`\mathrm{\Omega }_0(t)`$ represents an extrapolation of the pre-glitch frequency evolution, $`\mathrm{\Delta }\mathrm{\Omega }_{t=0}`$ is the magnitude of the change in the rotational frequency at the time of the glitch ($`t=0`$), $`\tau `$ is a characteristic exponential healing or recovery time after the glitch, and $`Q`$ is the fraction of the initial frequency change that is eventually recovered (known as the healing parameter). The last term accounts for any permanent (not recovered) change $`\mathrm{\Delta }\dot{\mathrm{\Omega }}_p`$ in the frequency derivative after the glitch. The parameter $`Q`$ in Equation 1 can be determined observationally (Shapiro & Teukolsky, 1983) from the jump in $`\mathrm{\Omega }`$ and its first and second time derivatives:

$$Q=\frac{\mathrm{\Delta }\dot{\mathrm{\Omega }}_{t=0}^2}{\mathrm{\Delta }\ddot{\mathrm{\Omega }}_{t=0}\mathrm{\Delta }\mathrm{\Omega }}_{t=0}.$$
(2)
In the starquake model, $`Q`$ is related to the moments of inertia of the components of the star according to (Pines, Shaham, & Ruderman, 1974):

$$Q=\frac{I_{\mathrm{core}}}{I_{\mathrm{total}}},$$
(3)
where $`I_{\mathrm{core}}`$ and $`I_{\mathrm{total}}`$ are respectively the moment of inertia of the superfluid core and the moment of inertia of the entire star. We have constructed a series of models of the NS interior using representative high-density equations of state (EOSs) in order to test whether the calculated moments of inertia in the models satisfy the prediction of the starquake model for Crab and Vela glitches. We have compared the model results with the measured values of $`Q`$ from Crab and Vela glitches to see if they are consistent.

## 2 Equations of State

We have constructed NS models using seven representative EOSs of superdense matter from which parameterized model stars could be produced. For densities $`\rho {}_{}{}^{}{}_{}{}^{<}1.6\times 10^{14}`$ g cm-3 (crustal densities), we have used three EOSs for three separate density regimes (see, e.g., Shapiro & Teukolsky 1983). These three EOSs are described by Feynman, Metropolis, & Teller (1949), Baym, Pethick, & Sutherland (1971), and Baym, Bethe, & Pethick (1971). The seven high-density EOSs, which dominate the macroscopic characteristics of the model stars, are briefly described below. Further details about each EOS can be found in the sources referenced.

### 2.1 The BJW Equation of State

For densities above the nuclear density, a combination of EOSs presented by Canuto (1975) was used for the modeling. For densities up to $`5\times 10^{15}`$ g cm-3, the EOS described by Bethe & Johnson (1974) was employed. This EOS is an improvement upon the work of Reid (1968) which includes a repulsive nucleonic core arising from the exchange of vector mesons in a hyperonic liquid. At higher densities, the EOS of Walecka (1974) was used, although the NS models produced in this density regime were beyond the stability limit. In this EOS, nucleons interact attractively via exchange of scalar mesons and repulsively via exchange of more massive vector mesons. Summaries of these EOSs are also presented by Shapiro & Teukolsky (1983). We call this EOS BJW.

### 2.2 The FPS Equation of State

Lorenz, Ravenhall, & Pethick (1993) refer to an EOS calculated by Friedman & Pandharipande (1981) for high-density neutron and nuclear matter. This EOS employs the microscopic V14 two-body potential, a three nucleon interaction (TNI) potential, and uses hypernetted chain techniques. This EOS was modified by fitting the microscopic interaction of Friedman & Pandharipande to a Skyrme-like energy density function (Skyrme, 1959). The essential feature of the Skyrme model is a two-body interaction that has the spatial character of a two-body delta function plus derivatives. We refer to this EOS as FPS.

### 2.3 GWM Equation of State

This EOS is based on a variant of the theory of nuclear field coupling (Zimanyi & Moszkowski, 1990) in which the scalar field is coupled to the derivative of the nucleon field. A modification proposed by Glendenning, Weber, & Moszkowski (1992), called the hybrid derivative coupling model, replaces the purely derivative coupling of the scalar field to baryons and vector mesons with a Yukawa point and derivative coupling to baryons and both vector fields. The coupling model is consistent with the experimentally inferred binding energy of lambda hyperons in nuclear matter. The resulting EOS, called GWM here, maintains equilibrium between all baryons to convergence and leptons. The existence of hyperons significantly softens this EOS.

### 2.4 HKP Equation of State

This EOS, proposed by Haensel, Kutschera, & Proszynski (1981), supposes that a minimal condition for the properties of high-density matter is that it must produce the observed values for the saturation properties of nuclear matter. The relativistic mean field theory of Serot (1979) is assumed in which nucleon interactions are governed by the exchange of neutral scalar and neutral vector mesons, pions, and rho mesons. This theory is an extension of the relativistic mean field theory of Walecka (1974). The authors assume a conservative range of possible nuclear saturation densities, the mean of which is $`2.8\times 10^{14}`$ g cm-3. The EOS used here, called HKP, is based on this nuclear saturation value.

### 2.5 WFF Equation of State

Wiringa, Fiks, & Fabrocini (1988) describe a microscopic EOS of dense nuclear matter constrained by nucleon-nucleon scattering data. The interaction includes a two-nucleon Urbana $`v_{14}`$ (UV14) potential and a TNI three-body potential term (Lagaris & Pandharipande, 1981a, b). The three-body potential includes a repulsive term, the primary effect of which is a reduction in the intermediate-range attraction of the two-nucleon potential, and an attractive term, which becomes negligible at high densities. We refer to this EOS as WFF.

### 2.6 APR Equations of State

Akmal, Pandharipande, & Ravenhall (1998) describe a set of realistic EOSs based on the Argonne $`v_{18}`$ (A18) two-nucleon interaction (Wiringa, Stoks, & Schiavilla, 1995), calculated using variational chain summation methods. This is supplemented with a relativistic boost correction term (Forest, Pandharipande, & Friar, 1995) and a TNI term based on the Urbana IX model (Pudliner et al., 1995). The two cases of pure neutron matter (PNM) and symmetric nuclear matter (SNM) are separately considered. The latter is composed of equal numbers of neutrons and protons in beta equilibrium. The two EOSs used to construct the models presented here were derived from a parabolic fit to and extrapolation of the tabulated nucleon density and energy values given for PNM and SNM by Akmal, Pandharipande, & Ravenhall (1998). We refer to the two EOSs based on PNM and SNM as APR(p) and APR(s), respectively.

## 3 The Modeling Procedure

To construct parameterized NS models for each EOS, we first chose a central density $`\rho _c\rho (0)`$ for each stellar model. The EOS and the relativistic expressions for stellar structure (Oppenheimer & Volkoff, 1939) were used to calculate the pressure gradient and incremental mass contained in concentric shells as the model iterated outward from the center of the star. The relativistic structure equations are:

$$\frac{dm(r)}{dr}=4\pi r^2\rho (r)$$

$$\frac{dP(r)}{dr}=\left(\rho (r)+\frac{P(r)}{c^2}\right)\frac{G\left[m(r)+4\pi r^3P(r)/c^2\right]}{r^2\left[12Gm(r)/rc^2\right]},$$
(4)
where $`m(r)`$ represents the mass of the star internal to a distance $`r`$ from the center of the star. When the pressure reached zero, $`P(R)=0`$, the edge of the star was reached ($`r=R`$) and the relevant macroscopic parameters could be read off.

To calculate the relativistic moment of inertia $`I(r)`$ interior to a radius $`r`$, expressions given by Arnett & Bowers (1977) were used which account for the Lense-Thirring frame-dragging effect (e.g., Glendenning 2000) (see discussion below):

$$\frac{dI(r)}{dr}=\frac{8\pi }{3}\frac{[\rho (r)+P(r)/c^2]r^4e^{\varphi (r)}}{[12Gm(r)/rc^2]^{1/2}}\frac{\overline{\omega }(r)}{\mathrm{\Omega }}.$$
(5)
Here $`\mathrm{\Omega }`$ is the observed stellar rotational frequency and $`\overline{\omega }(r)`$ is the angular velocity of the star at a distance $`r`$ relative to the angular velocity of the rotating local inertial frame, which is dragged at angular velocity $`\omega (r)`$:

$$\overline{\omega }(r)=\mathrm{\Omega }\omega (r).$$
(6)
The term $`e^{\varphi (r)}`$ in Equation 5 translates $`\overline{\omega }(r)`$ as measured from infinity to $`\overline{\omega }(r)`$ as measured in the local inertial frame at $`r`$. The following auxiliary relations were used to determine values of $`\varphi (r)`$ and $`\overline{\omega }(r)`$ throughout the star ($`r<R`$) in the moment of inertia calculation:

$$\frac{d\varphi (r)}{dr}=\left(\frac{1}{\rho (r)c^2+P(r)}\right)\frac{dP(r)}{dr}$$
(7)
$$j(r)=\left(1\frac{2Gm(r)}{rc^2}\right)e^{\varphi (r)}$$
(8)
$$\frac{d}{dr}\left(r^4j(r)\frac{d\overline{\omega }(r)}{dr}\right)+4r^3\frac{dj(r)}{dr}\overline{\omega }(r)=0.$$
(9)
An arbitrary value of $`\overline{\omega }(0)`$ was chosen and the following boundary conditions imposed in the models:

$$\frac{d\overline{\omega }(0)}{dr}=0$$
(10)
$$\overline{\omega }(R)=\mathrm{\Omega }\frac{R}{3}\frac{d\overline{\omega }(R)}{dr}.$$
(11)
This latter condition (in which a value for $`\mathrm{\Omega }`$ is determined at the end of the model iteration) was used to scale the arbitrarily chosen $`\overline{\omega }(0)`$ in such a way as to yield the proper observed $`\mathrm{\Omega }`$ in a subsequent modeling pass. This subsequent pass produced the correct $`\overline{\omega }(r)`$ for all $`r<R`$. The condition of finite $`\varphi (0)`$ was also imposed.

Milne’s centered algorithm (Harrison et al., 1965) was used for the iteration. Since there is a transition region between the superfluid core and the rigid crust which consists of a mixture of nuclei and free baryons, there is no clear division where the core-crust transition occurs. We chose a transition density near the nuclear density ($`\rho _{\mathrm{transition}}2.4\times 10^{14}`$ g cm-3) in order to separate the two components in our models. From this we obtained $`I_{\mathrm{core}}I(r_{\mathrm{transition}})`$ and $`I_{\mathrm{total}}I(R)`$ for each model NS and could compare the ratio of these to the measured values of $`Q`$ from glitch observations. In previous work, some authors (e.g., Datta & Alpar 1993; Shapiro & Teukolsky 1983) have assumed a slightly lower core-crust transition density than the one used here. This lower transition density reduces the size of the crust, thereby increasing the fractional moment of inertia of the core. Although slight differences in the choice of the transition density do not affect our model results, a significant reduction would lead to a more stringent test of the starquake model (i.e., it would further limit the range of observed $`Q`$ that would be consistent with the predictions of the starquake model). Thus, our choice of $`\rho _{\mathrm{transition}}`$ is a conservative one for this purpose.

For angular rotation speeds comparable to those of the Crab and Vela pulsars ($`\mathrm{\Omega }190`$ rad s-1 and $`70`$ rad s-1, respectively), the macroscopic star parameters are affected by first-order rotation effects, but not significantly so by second-order effects.

The Lense-Thirring frame dragging effect is a first-order effect which scales with rotational speed as $`(\mathrm{\Omega }/\mathrm{\Omega }_c)`$ (Arnett & Bowers 1977), where $`\mathrm{\Omega }_c=(GM/R^3)^{1/2}`$, the critical value above which centrifugal mass shedding would occur for the rotating star. For realistic NS models, $`\mathrm{\Omega }_c13000`$ rad s-1. This first-order effect is therefore $``$ 1-2% for the Crab and Vela pulsars and is taken into consideration in the calculation of the relativistic moment of inertia in the models.

NS mass shifts and deformations in sphericity from rotation (which would complicate the simple spherical stellar modeling presented here) are second order effects, scaling as $`(\mathrm{\Omega }/\mathrm{\Omega }_c)^2`$ (Arnett & Bowers 1977). This term would be significant for millisecond pulsars, but for the Crab pulsar, the second-order effect is $`10^4`$, and it is even lower for Vela. Thus, the second-order effect can be safely ignored when considering the macroscopic parameters of our NS models.

## 4 Results

We have produced a set of parameterized NS models using the seven high-density EOSs described above. In order to decide which models in each series were appropriate to use in the comparison with the observed values of $`Q`$, we used a realistic NS mass range as a reasonable constraint. Models with a NS mass outside of the range were not used in the comparison.

Thorsett & Chakrabarty (1999) have presented a statistical study of NS masses for the known binary radio pulsar population for which useful mass constraints could be derived. They find that, statistically, the distribution of pulsar masses falls in the range $`M=1.35\pm 0.04M_{}`$. This is consistent with the initial NS mass function from simulations of supernovae, which predict that NSs are formed with $`M{}_{}{}^{}{}_{}{}^{>}1.2M_{}`$ (Timmes, Woosley, & Weaver, 1996). For our comparison and analysis, we consider only the mass range $`M=1.4\pm 0.2M_{}`$, which we assume to be a conservative estimate of the realistic NS mass range.

Figure 1 shows NS mass $`Mm(R)`$ as a function of central density $`\rho _c\rho (0)`$ for our models constructed using the seven different EOSs described above. NS branches with $`dM/d\rho _c>0`$ satisfy a necessary stability criterion (e.g., Arnett & Bowers 1977) and are indicated with solid bold lines. For all EOSs used, there are stable model stars in the mass range $`1.4\pm 0.2M_{}`$. Figure 2 shows NS mass $`M`$ as function of radius $`R`$ for the seven EOSs used. Stable NS branches are again indicated with solid bold lines. None of the models in the stable branches violates the required stability condition $`M/R<4c^2/9G`$ for static relativistic stars (Glendenning, 2000).

Figure 3 shows the ratio of the moment of inertia of the core, $`I_{\mathrm{core}}`$, to the total moment of inertia, $`I_{\mathrm{total}}`$, as a function of NS mass $`M`$. For stable models in the mass range of interest ($`1.4\pm 0.2M_{}`$), this ratio is confined to $`I_{\mathrm{core}}/I_{\mathrm{total}}{}_{}{}^{}{}_{}{}^{>}0.55`$. Starquakes cannot be responsible for producing glitches in a pulsar with measured values of $`Q`$ consistently outside this range unless the pulsar mass is significantly smaller than expected.

Glitch measurements (including measurements of $`Q`$) for the Crab and Vela pulsars have been published in various places in the literature (see Table 1 for references). However, no comprehensive current listing of all measured values of $`Q`$ for Crab and Vela glitches exists. Table 1 is a compilation of all measured values of $`Q`$ published to date for glitches from the Crab and Vela pulsars. Some are different measurements of the same glitch. This complete set of measurements can be used to derive a range of observed $`Q`$ for each pulsar. From the 21 measurements of $`Q`$ for Crab glitches, a weighted mean of the values yields $`Q=0.72\pm 0.05`$111For measurements of $`Q`$ without a quoted uncertainty, an uncertainty of 0.1 was simply assumed in the calculation of the weighted mean.. This is only slightly smaller than the unweighted mean of $`Q=0.83`$. A range of $`Q{}_{}{}^{}{}_{}{}^{>}0.7`$ encompasses the observed distribution for the Crab pulsar. A weighted mean of the 11 measurements of $`Q`$ for Vela glitches yields $`0.12\pm 0.07`$, while an unweighted mean gives $`Q=0.18`$. All estimates for Vela agree that $`Q`$ is small, with a likely range $`Q{}_{}{}^{}{}_{}{}^{<}0.2`$.

## 5 Discussion

The glitch behavior of the Crab (Alpar et al., 1994, 1996) and Vela (Alpar et al., 1993; Chau et al., 1993) pulsars has been extensively studied previously. Below we discuss and compare our model results with the observed range of $`Q`$ in the context of this previous work.

### 5.1 Vela

Alpar et al. (1993) and Chau et al. (1993) have found that their model predictions of vortex unpinning in Vela are completely consistent with the observed Vela glitch characteristics. They deduce that the fractional crustal moment of inertia must be greater than 2.6% for Vela (Chau et al., 1993). Datta & Alpar (1993) have used a corresponding lower limit of 3.4% for Vela in an attempt to rule out soft EOSs. Link, Epstein, & Lattimer (1999) employ a lower limit of 1.4% for a similar purpose. These estimates imply an upper limit to the fractional core moment of inertia of $`0.98`$, which tends to disfavor APR(s) models with mass $`M{}_{}{}^{}{}_{}{}^{>}1.5M_{}`$, but cannot independently constrain the Vela mass since this condition is satisfied for a wide range of model masses (see Figure 3). The conclusions made by Alpar et al. (1993) and Chau et al. (1993) that the Vela mass is probably less than $`1.4M_{}`$ (and is more likely closer to $`1.2M_{}`$) cannot be confirmed or ruled out. Lyne (1992) explains one reason why starquakes cannot be the glitch mechanism for Vela: the required oblateness is not sustainable given the large size of Vela glitches ($`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }10^6`$). After $``$ 100 years, the oblateness would reach zero, and therefore glitches could only have been sustained for $``$ 1% of the current Vela age. The time interval between Vela glitches is also inconsistent with the much longer intervals predicted by the starquake model (Alpar et al., 1996) unless a solid core model is invoked for Vela (e.g., Canuto & Chitre 1973; Pines, Shaham, & Ruderman 1974). Starquakes would also produce a thermal energy dissipation during the large Vela glitches which is expected to be observable as a change in X-ray luminosity soon after the glitch occurs. X-ray observations of Vela show no such signal down to limits of less than a few percent in the fractional change in flux (Seward et al., 2000; Helfand, Gotthelf, & Halpern, 2001). The vortex unpinning model does not have this energy dissipation problem for Vela (Alpar, 1995). The low values of $`Q`$ measured for Vela glitches are additional evidence that the predictions of the starquake model do not match observations; otherwise the implied Vela mass would be too low ($`M{}_{}{}^{}{}_{}{}^{<}0.5M_{}`$ for $`Q{}_{}{}^{}{}_{}{}^{<}0.2`$ in our models). Our model results confirm previous conclusions that starquakes cannot be a feasible glitch mechanism for the Vela pulsar.

### 5.2 Crab

Alpar et al. (1994, 1996) have studied the behavior of Crab glitches and have determined that the starquake model is at least partially responsible for them. The magnitudes of Crab glitches, combined with the observed glitch rate, support this notion. Lyne (1992) indicates that the small size of Crab glitches ($`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }10^8`$) allows for small changes in oblateness which would not significantly deplete the oblateness over the current lifetime of the Crab. Link, Franco, & Epstein (1998) and Franco, Link, & Epstein (2000) also show that the permanent post-glitch offsets in the period derivative seen for the Crab pulsar can be accounted for by starquakes: the net torque on the star increases through shearing effects. Our model results are also consistent with the starquake interpretation for Crab glitches. The fractional moment of inertia values in our models are comparable to observed $`Q`$ values for Crab glitches ($`Q{}_{}{}^{}{}_{}{}^{>}0.7`$) if $`M{}_{}{}^{}{}_{}{}^{>}0.15M_{}`$, which seems likely. Alpar et al. (1994) set a lower limit of 0.2% for the fractional crustal moment of inertia. This is not constraining in our models since the corresponding fractional core moment of inertia range $`I_{\mathrm{core}}/I_{\mathrm{total}}{}_{}{}^{}{}_{}{}^{<}0.998`$ is easily satisfied. Our results support previous suggestions that starquakes could be responsible for Crab glitches.

## 6 Conclusions

Using parameterized NS models produced from seven representative EOSs of superdense matter, we find that the fractional moment of inertia of the core component of the model stars is $`I_{\mathrm{core}}/I_{\mathrm{total}}{}_{}{}^{}{}_{}{}^{>}0.55`$ for all stable configurations in the assumed realistic NS mass range $`1.4\pm 0.2M_{}`$ (Figure 3). This ratio, which is predicted by the starquake model to equal the glitch healing parameter $`Q`$, is not consistent with the observed range $`Q{}_{}{}^{}{}_{}{}^{<}0.2`$ for Vela glitches (see Table 1), unless the Vela pulsar mass is unrealistically small ($`M{}_{}{}^{}{}_{}{}^{<}0.5M_{}`$). This confirms results from previous studies of the Vela pulsar which indicate that starquakes are not the cause of Vela glitches (Lyne, 1992; Alpar et al., 1993; Chau et al., 1993; Alpar, 1995). The much larger values of $`Q{}_{}{}^{}{}_{}{}^{>}0.7`$ seen for Crab glitches (Table 1) are consistent with the moment of inertia values of our models for realistic masses, as predicted by the starquake model. These results support previous conclusions from the analysis of Crab glitch behavior (Alpar et al., 1994, 1996; Link, Franco, & Epstein, 1998; Franco, Link, & Epstein, 2000) in which starquakes have been proposed as the Crab glitch mechanism. Repeated and accurate measurements of $`Q`$ for other glitching pulsars in the future (if obtainable) could be used as a simple test of the starquake glitch model and may help resolve whether differences between Crab-like and Vela-like glitches can be understood on evolutionary grounds (e.g., Alpar 1995).

We thank the referee for helpful comments to improve the manuscript and for suggesting the inclusion of additional equations of state in the modeling. MD was partially supported by a grant of the Polish Committee for Scientific Research.


text:
# Fast and robust two–qubit gates for scalable ion trap quantum computing

## Abstract

We propose a new concept for a two–qubit gate operating on a pair of trapped ions based on laser coherent control techniques. The gate is insensitive to the temperature of the ions, works also outside the Lamb–Dicke regime, requires no individual addressing by lasers, and can be orders of magnitude faster than the trap period.

Trapped ions constitute one of the most promising systems to implement scalable quantum computation. PhysicsTodayMay2003  In an ion trap quantum computer qubits are stored in long-lived internal atomic states. A universal set of single and two qubit gates is obtained by manipulating the internal states with lasers, and entangling the ions via the motional states CiracZoller1 . During the last years a remarkable experimental progress in building an ion trap quantum computer has allowed to realize two–qubit gates Wineland95 ; Wineland1 ; Wineland2 ; Blatt  and also to prepare entangled states Wineland3 ; Wineland4 ; Wineland5 . The ultimate challenge is now the development of scalable ion trap quantum computing. Scalability is based on storing a set of ions, and moving ions independently, in particular to bring together pairs of ions to perform a two-qubit gate Scalable-a ; Scalable-b . Basic steps towards this goal have already been demonstrated experimentally papermotionWineland .

An important question to be addressed is to identify the current limitations of the two–qubit gates with trapped ions (given the fact that one–qubit gates are significantly simpler with those systems). The ideal scheme should: (i) be independent of temperature (so that one does not need to cool the ions to their ground state after they are moved to or from their storage area); (ii) require no addressability (to allow the ions to be as close as possible during the gate so as to strengthen their interaction), and (iii) be fast (in order to minimize the effects of decoherence during the gate, and to speed up the computation). This last property has been identified PhysicsTodayMay2003  as a key limitation: in essentially all schemes suggested so far CiracZoller1 ; Proposals-b ; Proposals-c ; Proposals-d ; ProposalsMilburn ; Proposals-a  one has to resolve spectroscopically the motional sidebands of the ions with the exciting laser, which limits the laser intensity and therefore the gate time.

The two–qubit gate between pairs of ions analyzed below solves the problem of speed by using mechanical effects instead of spectral methods to couple the motion and internal states of the ions. In this way the new limits on the time of the quantum gate are those of laser control, which can be orders of magnitudes faster than the present limits dictated by trap design. Thus, the method proposed here is a significant step forward towards fast and efficient scalable quantum computations with trapped ions.

Below we will first study the dynamics of a pair of trapped ions under the influence of short laser pulses with varying directions. We will prove that there exist certain laser pulse sequences which perform a phase gate on the two qubits, while leaving the motional state unchanged. We illustrate this with two protocols for laser pulses: (i) a sequence of four pulses which gives a gate time of $`T=1.08/\nu `$ with $`\nu `$ the trap frequency, and (ii) a protocol which allows us to perform a gate in a time $`TN_P^{2/3}/\nu `$ where $`N_P`$ is the number of laser pulses. Finally, we will complement our study of the gate dynamics with an analysis of possible errors, which includes fluctuations of the intensity or the duration of the pulses, and temperature. The gate will be shown to be extremely robust to these perturbations.

We consider two ions in a one–dimensional harmonic trap, interacting with a laser beam on resonance. The Hamiltonian describing this situationCiracZoller1  can be written as $`H=H_0+H_1`$, where $`H_0=\nu _ca^{}a+\nu _rb^{}b`$ describes the motion in the trap and

$`H_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }(t)}{2}}\sigma _1^+e^{i\eta _c(a^{}+a)+{\scriptscriptstyle \frac{1}{2}}\eta _r(b^{}+b)}+`$ (1)
$`{\displaystyle \frac{\mathrm{\Omega }(t)}{2}}\sigma _2^+e^{i\eta _c(a^{}+a){\scriptscriptstyle \frac{1}{2}}\eta _r(b^{}+b)}+h.c.`$
Here, $`\nu _c=\nu `$ and $`\nu _r=\sqrt{3}\nu _c`$ are the frequencies of the center of mass and stretching mode, respectively; $`a`$ and $`b`$ are the corresponding annihilation operators, and $`\eta _c=\eta /\sqrt{2}`$ and $`\eta _r=\eta \sqrt[4]{4/3}`$ are proportional to the Lamb–Dicke parameter, $`\eta `$. Note that the Rabi frequency $`\mathrm{\Omega }`$ is the same for both ions, since we have not assumed individual addressing. Also notice that replacing $`\eta `$ with $`(\eta )`$ is equivalent to reversing the direction of the laser beam.

In the following we will consider two different kind of processes: (i) Free evolution, where the laser is switched off ($`\mathrm{\Omega }=0`$) for a certain time; (ii) Sequences of pairs of very fast laser pulses, each of them coming from opposite sides. If we denote by $`\delta t`$ the duration of a pulse and by $`\mathrm{\Omega }`$ the corresponding Rabi frequency, we are interested in the limit $`\delta t0`$ with $`\mathrm{\Omega }\delta t=\pi `$. Processes (i) and (ii) will be alternated (See Fig. 1a): at time $`t_1`$ a sequence of $`z_1`$ pulses is applied, followed by free evolution until at time $`t_2`$ another sequence of $`z_2`$ pulses is applied followed by free evolution and so on. The $`z_k`$ are integer numbers, whose sign indicates the direction of the laser pulses.

For a pulse sequence, consisting of kicks interspersed with free harmonic time evolution (Fig. 1), we write the evolution operator as $`𝒰=𝒰_c𝒰_r`$, where $`𝒰_{c,r}=_{k=1}^NU_{c,r}(\mathrm{\Delta }t_k,z_k)`$ has contributions of the center–of–mass and relative motions,

$`U_c(t_k,z_k)`$ $`=`$ $`e^{i2z_k\eta _c(a+a^{})(\sigma _1^z+\sigma _2^z)}e^{i\nu _c\mathrm{\Delta }t_ka^{}a},`$
$`U_r(t_k,z_k)`$ $`=`$ $`e^{iz_k\eta _r(b+b^{})(\sigma _1^z\sigma _2^z)}e^{i\nu _r\mathrm{\Delta }t_kb^{}b}.`$
The integers $`z_k`$ indicate the direction of the initial pulse in the sequence of pairs of very fast laser pulses, each of them coming from opposite sites.

In order to fully characterize $`U`$, we only have to investigate its action on states of the form $`|i_1|j_2|\alpha _c|\beta _r`$, where $`i,j=0,1`$ denote the computational basis, and $`|\alpha `$ and $`|\beta `$ are coherent states. This task can be easily carried out once we know the action of $`𝒰=_{k=1}^NU(\varphi _k,p_k)`$ on an arbitrary coherent state $`|\alpha `$, where

$$U(\varphi _k,p_k)=e^{ip(a+a^{})}e^{i\varphi _ka^{}a}.$$

We obtain $`𝒰|\alpha =e^{i\xi }|\stackrel{~}{\alpha }`$, where

$`\stackrel{~}{\alpha }`$ $`=`$ $`\alpha e^{i\theta _N}i{\displaystyle \underset{k=1}{\overset{N}{}}}p_ke^{i(\theta _k\theta _N)},`$
$`\xi `$ $`=`$ $`{\displaystyle \underset{m=2}{\overset{N}{}}}{\displaystyle \underset{k=1}{\overset{m1}{}}}p_mp_k\mathrm{sin}(\theta _k\theta _m)\mathrm{}\left[\alpha {\displaystyle \underset{k=0}{\overset{N}{}}}p_ke^{i\theta _m}\right],`$
with $`\theta _k=_{m=1}^k\varphi _m`$.

The crucial point is to realize that if $`_{k=1}^Np_ke^{i\theta _k}=0`$ the motional state $`|\alpha `$ after the evolution is the same as if there was only free evolution \[Fig. 1a\], and a global phase $`\xi `$ appears which does not depend on the motional state. Translating this result to the operators $`𝒰_c|\alpha `$ and $`𝒰_r|\beta `$, we obtain the conditions

$`C_c`$ $``$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}z_ke^{i\nu t_k}=0,`$ (3a)
$`C_r`$ $``$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}z_ke^{i\sqrt{3}\nu t_k}=0.`$ (3b)
If these commensurability conditions are satisfied, the motional state of the ion will not depend on the qubits and the evolution operator will be given by

$$𝒰(\mathrm{\Theta })=e^{i\mathrm{\Theta }\sigma _1^z\sigma _2^z}e^{i\nu _cTa^{}a}e^{i\nu _rTb^{}b}.$$
(4)
The value $`T`$ is the total time required by the gate and

$$\mathrm{\Theta }=4\eta ^2\underset{m=2}{\overset{N}{}}\underset{k=1}{\overset{m1}{}}z_kz_m\left[\frac{\mathrm{sin}[\sqrt{3}\nu \mathrm{\Delta }t_{km}]}{\sqrt{3}}\mathrm{sin}(\nu \mathrm{\Delta }t_{km})\right],$$
(5)
where $`\mathrm{\Delta }t_{km}=t_kt_m`$. Therefore, if Eqs. (3) are fulfilled, and $`\mathrm{\Theta }=\pi /4`$ we will have a controlled–phase gate (which is equivalent to a controlled–NOT gate up to local operations) which is completely independent of the initial motional state, i.e. there are no temperature requirements.

It is straightforward to show that for any value of the time $`T`$ it is always possible to find a sequence of laser pulses which implements the gate, and therefore the gate operation can be, in principle, arbitrarily fast. The search for a sequence of pulses may be done numerically, or even semianalytically. In the following we give two simple protocols which are not optimized in order to reduce the number of pulses.

The first protocol (Protocol I) requires the least number of pulses to produce the gate in a fixed time $`T1.08(2\pi /\nu )`$. The recipe is illustrated in Fig. 1, which provides the phase space plots for the evolution of the motional state. The sequence of pulses is defined as $`(z_n/N,t_n)=\{(\gamma ,\tau _1),(1,\tau _2),(1,\tau _2),(\gamma ,\tau _1)\}`$. Here $`0<\gamma =\mathrm{cos}(\theta )<1.0`$ is a real number, which may be introduced by tilting both lasers a small angle $`\theta `$ with respect to the axis of the trap, so that no transverse motion is excited. It is always possible to find a solution to Eq. (3) with $`\tau _10.538(4)(2\pi /\nu )>\tau _2>0`$. The results for the performance of the gate are summarized in Fig. 2. As shown in Fig. 2(a), for realistic values of the Lamb–Dicke parameter Wineland1  we only need to apply the sequence of pulses one or two times to implement a phase gate.

The second protocol (Protocol II) performs the gate in an arbitrarily short time $`T`$. The pulses are now distributed according to $`(z_n/N,t_n)=\{(2,\tau _1),(3,\tau _2),(2,\tau _3),(2,\tau _3),(3,\tau _2),(2,\tau _1)\}`$. The whole process takes a time $`T=2\tau _1`$ and requires $`N_p=|z_n|=14N`$ pairs of pulses. As Fig. 2 shows, the number of pulses increases with decreasing time as $`N_pT^{3/2}`$.

In order to study the main potential limitations of our scheme, we define the error of the gate $`E`$ in terms of the gate fidelity Nielsen  as $`E=1\mathrm{Tr}_{\mathrm{mot}}\left\{Q_{\mathrm{mot}}\rho _{\mathrm{mot}}Q_{\mathrm{mot}}^{}\right\}`$. Here $`\mathrm{Tr}_{\mathrm{mot}}`$ and $`\mathrm{Tr}_{\mathrm{int}}`$ denote traces over motional and internal degrees of freedom, and $`Q_{\mathrm{mot}}=\mathrm{Tr}_{\mathrm{int}}\left\{𝒰(\pi /4)U_{\mathrm{real}}^{}\right\},`$ depends on $`U_{\mathrm{real}}`$, the gate performed in the presence of imperfections.

We now turn to a discussion of the possible sources of errors. A limiting factor for the gate is the anharmonicities of the restoring forces. The more pulses we apply, the larger the relative displacement of the ions, as Fig. 3(b) shows. When the ions become too close to each other, the increasing intensity of the Coulomb force can lead to a breakdown of the harmonic approximation which is implicit in Eq. (1). In order to analyze this effect, we have made a perturbation analysis for $`\nu T1`$ and found that such an anharmonicity causes an error $`E|0.4a_0/d|^2/(2\pi \nu T)`$, where $`a_0`$ is the ground state size of the external potential and $`d`$ is the ion separation in equilibrium. For typical parameters and imposing an error $`E10^4`$ we obtain $`\nu T10^3`$. A similar analysis could be applied to study anharmonicities of the trap itself.

In addition, we have studied the influence of errors in the laser pulses of our scheme. Up to now, our analytical calculations assumed that the intensity of the laser is very large during each pulse, and that therefore one may neglect the influence of the trap during this process. To validate this assumption we have simulated numerically a system of two ions with only one vibrational mode. We have used the exact sequences developed above to produce the phase gate using only eight laser pulses. In Fig. 4(a) we plot the error of the gate as a function of the duration of the laser pulse, $`\tau =\pi /2\mathrm{\Omega }`$. The longer the pulse, the more important the effect of the trap, and the larger the error. But even for relatively long pulses, we obtain a fidelity which is comparable to the results obtained in current setups Wineland1 ; Blatt ; Wineland2 . We have also studied the influence of noise in the intensity of the laser pulses, or, what is equivalent, random errors in its duration. The larger the amplitude of the error the lower the fidelity of the gate, as Fig. 4(b) shows.

As mentioned before, the scheme is insensitive to temperature. If the commensurability condition (3) is not perfectly satisfied due to, for example, errors in timing of laser pulses, or misalignment of the lasers, then the corresponding contribution to the gate error is

$$E=(C_1^4+C_2^4+C_1C_26)/8$$

with

$`C_1`$ $`=`$ $`\mathrm{exp}\left[(1/2+k_bT/\mathrm{}\nu _c)|2\eta _cC_c|^2\right]`$ (6)
$`C_2`$ $`=`$ $`\mathrm{exp}\left[(1/2+k_bT/\mathrm{}\nu _r)|\eta _rC_r|^2\right],`$ (7)
which is a smooth function of temperature $`T`$.

Finally, we would like to make some remarks regarding the experimental implementation of this scheme. First, it is not necessary to kick the atoms using pairs of counter-propagating laser beams. The same effect (i.e. a change of sign in $`\eta `$) may also be achieved in current experiments by reverting the internal state of both ions simultaneously. One then only needs a laser beam aligned with the trap to kick the atoms, and another laser orthogonal to the axis of the trap to produce the NOT-gate.

The second and more important remark is that it is possible to avoid errors in the laser pulses by using more sophisticated kicking methods. One possibility consists in using STIRAP STIRAP-a ; STIRAP-b . Only one of the qubit states would be connected by two on–resonance laser beams to a third atomic state, $`|e^{}`$. In the first part of the kicking process, the Rabi frequencies of the lasers $`\mathrm{\Omega }_a`$ and $`\mathrm{\Omega }_b`$ are adiabatically switched on an off respectively. The momenta of both laser beams should be different, so that as we slowly proceed from $`\mathrm{\Omega }_a/\mathrm{\Omega }_b0`$ to the opposite regime $`\mathrm{\Omega }_b/\mathrm{\Omega }_a0`$, the ions in the state $`|1`$ are completely transferred to the new dark state $`|e^{}`$ and get a kick $`|1e^{i(\stackrel{}{k}_a\stackrel{}{k}_b)\stackrel{}{x}}|e^{}`$. Next we must change the sense of the laser beams ($`\stackrel{}{k}_{a,b}\stackrel{}{k}_{a,b}`$), and perform the adiabatic transfer from $`\mathrm{\Omega }_b/\mathrm{\Omega }_a0`$ to $`\mathrm{\Omega }_a/\mathrm{\Omega }_b0`$. The total transformation will be $`|0|0`$, $`|1e^{2i(\stackrel{}{k}_a\stackrel{}{k}_b)\stackrel{}{x}}|1`$. The advantages of this method are: (i) the system remains all the time in a dark state, avoiding spontaneous emission; (ii) the process is insensitive to fluctuations of the intensity; (iii) the duration of the pulse need not be precisely adjusted, and (iv) the intensity of the laser need not be the same for both ions.

Summing up, in this work we have developed a new concept of two-qubit quantum gate for trapped ions, in which the trap frequency poses no longer a limitation on the speed of the gate. Rather than performing sideband transitions which weakly couple the internal and the motional states of the atoms, we suggest to push the atoms resonantly during very short times and along different directions. The limitations in that case come from: (i) the anharmonicities of the restoring force that the ions experience when pushed far away from each other, and (ii) the ability to control the laser pulses. The first limitation still allows to perform the gates in a time which is three orders of magnitude smaller than the one imposed by the trap frequency. The second one can be overcome by using adiabatic passage techniques which make these laser excitations tolerant against laser imperfections. In any case, the rapid experimental progress in laser control with very short pulses indicates that it may be possible soon to perform quantum gates with a very high speed. In addition, our scheme is independent of the temperature, requires no addressability, and works beyond the Lamb-Dicke regime.

J.I.C. thanks D. Leibfried for discussions and NIST for the hospitality during his stay in Boulder. P.Z. thanks the Innsbruck ion trap group for discussions. Part of this work was supported by the EU IST project RESQ, the project TOPQIP, the DFG (Schwerpunktprogramm Quanteninformationsverarbeitung) and the Kompetenznetzwerk Quanteninformationsverarbeitung der Bayerischen Staatsregierung. Research at the University of Innsbruck is supported by the Austrian Science Foundation, EU Networks and the Institute for Quantum Information.