text:
# Flux of atmospheric muons: Comparison between AIRES simulations and CAPRICE98 data.

## I Introduction

Detailed measurements and studies of the flux of muons in the atmosphere represent a subject of particularly great interest. This is mainly due to the fact that a measurement of the muon flux is an indirect measure of the neutrino flux and can therefore be used to improve the calculation of the atmospheric neutrino flux which in turn is used to compare with the observed neutrino rates in experiments targeted to detect neutrino oscillations. In addition, comparing measurements of fluxes of muons and other particles at different altitudes with simulated data is a powerful tool to check and/or calibrate air shower simulation programs. Such programs are not only essential to predict the atmospheric neutrino flux, but also play an important role in the analysis of data taken at highest energy air shower experiments like Auger AUGERW , AGASA AGASA  or HiRes HiRes . They are also used in Čerenkov light detector experiments like AMANDA AMANDA  or ANTARES ANTARES  to estimate backgrounds from atmospheric neutrinos.

Several calculations of the atmospheric muon flux at different altitudes have been performed in the past naumov ; fluka ; ralph . At ground level more data is available allowing also for simultaneous absolute flux and $`\mu ^+/\mu ^{}`$ flux ratio analysis willi ; wentzicrc ; wentzprd . The comparisons with experimental data show discrepancies in particular at low rigidities. The most significant differences between simulated and real data can be attributed to the following effects:

The normalization of the primary cosmic ray flux. Different measurements disagree in about 10 % or more on the proton and helium spectra between 10 and 50 GeV naumov .

The solar modulation, the geomagnetic effect and the model of the atmospheric profile.

The uncertainties in the particle production models, particularly the hadronic interaction generators muonpaper .

The recent introduction of Ring Imaging Čerenkov counters (RICH) in balloon-borne spectrometers has made it possible to measure in the large proton background also positive muons. The CAPRICE98 experiment with a gaseous RICH yo , measured the flux of positive (negative) muons with momenta in the range 0.3-20 GeV/$`c`$ (0.3-40 GeV/$`c`$), at different altitudes ranging from ground level (885 g/cm2 atmospheric overburden) up to float altitude (5.5 g/cm2).

In this work we make a comparison between the direct measurement of the flux of muons in the atmosphere by the CAPRICE98 experiment yo  and the corresponding simulated values obtained from the air shower simulation program AIRES AIRESManual . This program is a 3D Monte Carlo simulator where the majority of the processes that may undergo the shower particles, are taken into account. AIRES also has the advantage of including effects of the curvature of the Earth and of the geomagnetic field. The AIRES program has the possibility to swap between different hadronic models and it includes links to two well-known external hadronic interaction packages, namely SIBYLL SIBYLL  and QGSJET QGSJET .

This paper is organized as follows: In section II we give some details of the CAPRICE experiment. In section III we report on the method used to simulate the muon flux, using as input experimental values for the cosmic ray fluxes at the top of the atmosphere (TOA). The results of the comparison between experiment and simulations are reported in section IV and finally we summarize our conclusions in section V.

## II The experiment

The analysis presented in this work use the muon fluxes measured by the CAPRICE98 (Cosmic Anti-Particle Ring Imaging Čerenkov Experiment) experiment, at all the atmospheric depths given in table 1 yo . A detailed description of the procedures for muon identification, a careful study of the efficiency of each detectors, the various sources of background, and the rejection criteria and surviving contamination, can also be found in yo .

The balloon was launched from Ft. Summer, New Mexico, USA ($`34.28^{}`$N, $`104.14^{}`$W) on the 28th May 1998 amb99 . The following measurements were performed:

* Muon measurement at ground. Before launch the spectrometer took data at ground level for a period of about 14 hours. The ground level at Ft. Summer is located at an altitude of 1270 m that corresponds to an atmospheric depth of 885 g/cm2.
* Muon measurement during balloon ascending period. During the ascent period of a few hours measurements were done on the flux of particles as a function of momentum and atmospheric depth.
* Muon measurement during balloon floating period. During the 20 hours at float above 35 km, corresponding to approximately 5 g/cm2 of residual atmosphere, the balloon spectrometer recorded data on muons, protons and helium nuclei. This last portion of the flight is important because it provides the data used to estimate the primary flux of proton and helium nuclei at the top of the atmosphere used in the simulations.

At every altitude the experimental data available includes the $`\mu ^+`$ ($`\mu ^{}`$) flux for rigidities ranging from 0.3 GV up to 20 (40) GV. It is important to stress that this is the first experiment able to measure positive muons up to 20 GV. In the previous CAPRICE94 balloon flight boe99b ; muonpaper  the upper limit for $`\mu ^+`$ was 2 GV.

The CAPRICE98 spectrometer accepts particles arriving with an inclination with respect to the vertical axis of less than 20. This characteristic of the instrument needs to be taken into account when performing a simulation. It is important to point out that the axis of the spectrometer remained vertical during the flight.

## III The simulations

### III.1 Flux at the top of the atmosphere

The main input for the simulation of the flux of atmospheric muons at a given altitude is the absolute flux of cosmic rays at the top of the atmosphere. In the calculation we included the fluxes of the following 11 cosmic nuclei: H (protons and deuterium), He (He3 and He4), C, N, O, Ne, Mg, Si, and Fe.

The most important contribution to the total absolute flux at the top of the atmosphere comes from hydrogen and helium nuclei with only small contributions from other nuclei. It is also important to mention that photons and electrons do not contribute significantly to the flux of muons and therefore have not been included in our input.

We have mainly used hydrogen and helium fluxes obtained by the CAPRICE98 experiment emiliano , thus ensuring that the bulk of the input flux is affected by similar systematic errors as all the secondary particles, in particular muons, since all of them are measured with the same apparatus. This implies that when comparing experimental data at a given atmospheric depth with the corresponding simulations, a direct evaluation of the properties of the propagating algorithms is being performed, minimizing the uncertainty due to any inaccuracies in the input flux.

The CAPRICE98 experiment measured the absolute flux of cosmic protons (helium nuclei) with kinetic energies ranging from 3 to 350 GeV (0.9 to 170 GeV/nucleon) emiliano . These absolute fluxes can be adequately parameterized as polynomial functions of the logarithm of the primary energy. The deuterium to proton and He3 to He4 ratios were taken from reference isotopo . For heavier nuclei we have used data from reference gordo .

The absolute TOA fluxes of hydrogen (proton + deuterium), helium (He3 \+ He4), carbon, oxygen, and iron are shown as functions of the energy of the primary particles in figure 1. The lines represent the polynomial functions used as input for the simulations. The CAPRICE98 proton and helium data are also shown. The parameterized representations agree excellently with the experimental data.

It is important to remark that for low energies (below about 3 GeV) the solar modulation was taken into account when fitting the plotted functions, considering the level of solar activity registered at the date of the CAPRICE98 flight.

Notice also that the input fluxes used in the simulations correspond to total absolute fluxes at TOA without any correction due to geomagnetic effects. Such effects are taken into account by the simulation and analysis chain used to evaluate the different simulated observables: a corrective weight is applied to the secondaries generated by each primary entering the atmosphere to effectively take into account the geomagnetic cutoff. Then, the particles are propagated within the atmosphere, also taking into account their deflections due to the geomagnetic field, which is in this case taken as constant.

### III.2 Air shower simulations

The AIRES simulation engine AIRES ; AIRESManual  provides full space-time particle propagation in a realistic environment, taking into account the characteristics of the atmospheric density profile (using the standard US atmosphere atmosfera ), the Earth’s curvature, and the geomagnetic field (calculated for the location and date of the CAPRICE98 flight with an uncertainty of a few percent GF ).

The following particles are taken into account in the AIRES simulations: photons, electrons, positrons, muons, pions, kaons, eta mesons, lambda baryons, nucleons, antinucleons, and nuclei up to $`Z=36`$. Nucleus-nucleus, hadron-nucleus and photon-nucleus inelastic collisions with significant cross-sections are taken into account in the simulation. The hadronic processes are simulated using different models, accordingly to the energy: high energy collisions are processed invoking an external package (SIBYLL 2.1 SIBYLL  or QGSJET01 QGSJET ), while low energy ones are processed using an extension of Hillas splitting algorithm (EHSA) AIRES ; AIRESCCICRC ; AUGERJPG . The threshold energies separating the low and high energy regimes used in our simulations are 200 GeV and 80 GeV for the SIBYLL and QGSJET cases, respectively. The EHSA low energy hadronic model used in AIRES is a very fast procedure, effectively emulating the major characteristics of low energy hadronic collisions. The model is adjusted to retrieve similar results as the high energy hadronic model for energies near the transition thresholds previously mentioned, and the low energy cross sections are calculated from parameterizations of experimental data. A complete discussion on the low energy hadronic models is clearly beyond the scope of this paper. A separate report on this subject will be published elsewhere sjsinpreparation .

AIRES has been successfully used to study several characteristics of high energy showers, including comparisons between hadronic models AUGERJPG ; doqui99 , influence of the LPM effect LPM , muon bremsstrahlung mubrem , and geomagnetic deflections GF  on the shower development. AIRES has also been used to obtain an energy calibration of the AGASA experiment AGASAc , and to study the expected efficiency of the Auger Observatory for detecting quasi-horizontal showers generated by $`\tau `$-neutrinos Pierretau .

In the present work, AIRES has been used to simulate showers with primary energies from 7.5 $`\times `$ 108 up to 1015 eV. In order to accurately simulate the absolute fluxes already described in section III.1, and also to optimize the statistics, we conveniently divided the primary energy range into many subintervals with boundaries chosen so as to have at each of them approximately constant compositions and slope $`\gamma `$ ($`\gamma =d\mathrm{\Phi }/dE_{pr}`$ where $`\mathrm{\Phi }`$ is the flux and $`E_{pr}`$ is the primary energy). The independent sets of simulated showers were generated for each one of the intervals, considering also an isotropic arrival direction distribution with zenith angles ranging from 0 to 89.

The shower simulations performed for this study add up to more than 300 millions of showers, generating particle data files with a total size of about 30 GB, and requiring about 20 days of processing time (using a 1 GHz processor). The generation of such large set of simulated showers proved to be a straightforward computing exercise, where the AIRES system could be easily configured for this particular task, even if it was originally designed to simulate showers with significantly larger energies.

When configuring the simulation program, several aspects have been taken into account to properly set up the input parameters in our case of flux simulation. In particular, we would like to make the following remarks: (i) The low energy hadronic interactions increase their importance as the primary energy decreases. For this reason, we have done a careful setting of the parameters of the EHSA, taking into account experimental results, and comparisons with other models. (ii) The statistical sampling algorithm of AIRES (the so-called thinning) AIRESManual  was completely disabled. This means that all the secondaries generated during the shower development are fully propagated. (iii) All electromagnetic particles (photons, electrons, and positrons) with energies below 100 MeV were discarded. This significantly reduces the processing time required for a given simulation, without altering the propagation of hadrons and muons.

### III.3 Simulating the flux of secondary particles.

Let us consider a given observing level located at an altitude $`h`$ smaller than the injection altitude $`h_i`$. After simulating the shower in the conditions of section III.1, the secondary particles arriving to the observing level, are processed to estimate the corresponding fluxes, according to the following conditions:

When comparing with CAPRICE98 measurements yo , only particles reaching the observing level with an inclination of less than 20 have been considered. In this paper all particles satisfying that selection criterion are referred as quasi-vertical particles. Unless otherwise specified all the results discussed in the following sections apply to quasi-vertical particles. Notice also that full fluxes refer to fluxes of particles coming from all directions.

When necessary, the selected secondary particles are binned according to their momenta, using the momentum intervals of the experimental data.

The selected particles are weighted to take into account the absolute flux normalization associated with the primary energy subinterval (see section III.1) that corresponds to the shower being analyzed.

The procedure to evaluate the fluxes is repeated at each of the observing altitudes that are considered, using at each altitude an independent set of simulated showers.

## IV Results

### IV.1 General analysis of the simulated flux

The fluxes of secondary particles have been simulated for all altitudes listed in table 1. We have also performed simulations at sea level (1035 g/cm2). In figure 2 the total fluxes of muons, pions, protons, neutrons and helium nuclei are plotted as function of atmospheric depth and altitude above sea level. Electrons, positrons and photons were not propagated in detail in the simulation (see section III) and are therefore not considered in our study. Some of the main characteristics of these particle fluxes show up in this figure. At the highest level (5.5 g/cm2) protons are the most abundant particles because only a small fraction of them have interacted. On the other hand, at altitudes near sea level, the total flux is dominated by muons, which account for more than 96 % of the considered particles. Notice also that, as expected, the muon flux increases with decreasing altitude, until reaching a maximum around $`X=150`$ g/cm2. The muon flux then decreases with decreasing altitude until sea level. For atmospheric depths in the range 90 – 300 g/cm2 the flux of muons is quite constant. The pion flux behaves similarly, reaching its maximum approximately at the same altitude as that for the muons.

We have studied how the different primaries that make up the primary flux at the TOA, already described in section III.1, contribute to the total muon flux at different altitudes. In figure 3 the relative contributions of proton, helium nuclei, and other heavy nuclei (C, N, O, Ne, Mg, Si, Fe) to the muon flux are plotted versus atmospheric depth and altitude above sea level. It is evident that protons give the largest contribution at all altitudes, and their contribution increases when the altitude decreases. The heavy nuclei contribute with a small but not completely negligible fraction that ranges from about 8 % at sea level up to more than 10 % at $`X=5.5`$ g/cm2.

The contributions to the flux of muons with momenta larger than 10 GeV/$`c`$ is also displayed (open symbols) in figure 3. The different fractions are similar to the corresponding ones for the total muon flux, except at high altitudes where the proton fraction is substantially larger in comparison with the previous case.

It is also interesting to analyze how primaries with different energies (for nuclei total kinetic energy) contribute to the muon flux. To this end, we have divided the primary energies $`E_{pr}`$ into four ranges, namely, (1) $`E_{pr}`$ $`<`$ 10 GeV; (2) 10 GeV $`<`$ $`E_{pr}`$ $`<`$ 100 GeV; (3) 100 GeV $`<E_{pr}<1`$ TeV; and (4) $`E_{pr}`$ $`>`$ 1 TeV. The corresponding contributions to the muon flux using the data of our simulations are displayed in figure 4, where the relative contributions of each of the four primary energy categories to the full muon flux (including all arrival directions) are plotted versus atmospheric depth and altitude. As expected, in the energy range 10-100 GeV primaries, that account for the most significant fraction of particles capable of entering into the Earth’s atmosphere, are the ones that contribute most to the muon flux at all altitudes. Lower energy primaries contribute significantly (30 to 40 %) only at high altitudes, but their contribution decreases rapidly for $`X`$ longer than 100 g/cm2. Near sea level, the contributions of energy range 3 and 4 primaries increases significantly with a contribution from 100 GeV -1 TeV primaries of more than 30 % (20 % when selecting only quasi-vertical muons). We have studied the contribution of very low energy primaries ($`E_{pr}`$ $`<`$ 3 GeV) (the results have not been included in figure 4 for simplicity). Their contribution to the flux is always small ($`<1.5`$%) in the entire range of altitudes considered. This implies that any errors in the estimation of the input flux for very low energies are unlikely to have a significant impact on the simulated muon fluxes (notice that 3 GeV is the lowest primary energy measured by CAPRICE98).

To complete our study of how the different primary particles contribute to the flux of muons at different altitudes, we have analyzed the number of muons generated from primaries with given zenith angles. The results of this analysis are displayed in figures 5 and 6.

In figure 5 the relative contributions to the full flux of muons corresponding to primaries with zenith angle, $`\mathrm{\Theta }`$, less than 30 (circles), greater than 30 and less than 60 (squares), and greater than 60 (triangles), are plotted versus atmospheric depth and altitude. The solid (open) symbols correspond to all muons (muons with momentum larger than 10 GeV/$`c`$).

The curves in this figure present a complicated dependency on altitude of the corresponding contributions. At sea level, showers with $`30^{}<\mathrm{\Theta }<60^{}`$ dominate the contribution of the full flux, and the different contributions are independent of the muon momentum. At high altitudes the contribution of very inclined showers is dominating. This is due to the fact that such showers pass through a thicker layer of air and are more developed and producing more muons before reaching observation level, in comparison with vertical ones (see the discussion on angular distributions of muons later in this section). This effect is more important in the case of high energy muons, as indicated by the open symbol curves.

This picture changes dramatically when considering only quasi-vertical muons, as displayed in figure 6. In this case the showers with $`\mathrm{\Theta }<30^{}`$ are the ones that most contribute at all altitudes, and their relative contribution is larger for high energy particles, reaching virtually 100 % when selecting muons with momenta greater than 10 GeV/$`c`$. Notice that in figure 6 the open circles correspond to muons with momenta larger than 1 GeV/$`c`$.

Another very important observable that we have considered in our analysis, is the distribution of arrival directions of muons at different altitudes. As we have already commented in section II, and as explained in detail in reference yo , the normalization of the muon fluxes measured by CAPRICE98 is calculated under the assumption of an isotropic distribution of muons within the acceptance cone (zenith angle less than 20). In order to check the validity of this assumption, we have recorded angular distributions of muons at all the simulated levels and, simultaneously, calculated the fluxes using several aperture cones.

To start with, let us consider the angular distributions of all atmospheric muons. In figure 7 the normalized angular distributions of muons are plotted versus the cosine of the arrival zenith angle in the entire $`[0,1]`$ range for several altitudes namely (a) 5.5 g/cm2, (b) 77 g/cm2,(c) 308 g/cm2, and (d) 885 g/cm2.

We have found that in the entire range of altitudes considered the angular distributions can be accurately fitted by the following function:

$$f(\mathrm{cos}\theta )=U(\mathrm{cos}\theta )^{(\alpha +\beta \mathrm{cos}\theta )}$$
(1)
where $`\alpha `$ and $`\beta `$ are constants, and $`U`$ is a normalization factor determined by the condition $`_0^1f(x)𝑑x=1`$. This function can be conveniently fitted, with $`\alpha `$ and $`\beta `$ as free parameters, using simulated or experimental data.

The constants $`\alpha `$ and $`\beta `$ vary slowly as function of the atmospheric depth of the observing level as illustrated in figure 8. Notice that for $`x600`$ g/cm2 $`\beta 0`$, and therefore in this case the distribution reduces to a power of $`\mathrm{cos}\theta `$.

An important characteristic of the distributions plotted in figure 7 is that in the range $`\theta <20^{}`$ (corresponding to $`\mathrm{cos}\theta >0.94`$) the variations as functions of angle are small, in particular for high altitudes. This implies that the isotropy hypothesis is reasonably justified and at the same time, the fluxes estimated using aperture cones smaller than, say, 20, should not differ significantly.

This property shows up clearly from figures 9, 10 and 11 where the simulated $`\mu ^+`$ and $`\mu ^{}`$ fluxes calculated using different acceptance cones, are plotted versus muon momentum at atmospheric depths 22.6, 308, and 885 g/cm2, respectively. Notice that there are no important differences between the 10 and 20 cases.

From the distributions in figures 10 and 11 it is clear that the high momentum end of the spectra do not significantly depend on the aperture cone used. However, this is not the case for the distributions in figure 9 where it can be seen that the steepness of the distributions for large momenta decreases with aperture angle. This can be explained analyzing the angular distributions of high energy atmospheric muons. In figure 12, the angular distributions for (a) all simulated muons, (b) $`p_\mu >1`$ GeV/$`c`$, (c) $`p_\mu >10`$ GeV/$`c`$, (d) $`p_\mu >100`$ GeV/$`c`$, are shown for four different altitudes. We conclude from these figures that (i) the distributions for the lower altitudes, 5.5 and 77 g/cm2 show a large variation in angular distribution for different energies, (ii) at high altitudes $`\theta =0`$ does no more correspond to the maximum of the distribution (d) which is located rather close to $`\theta =90^{}`$. Therefore, when enlarging the aperture angle an increasing number of muons is accepted, consequently producing a larger flux at high momenta.

It is worthwhile mentioning that this feature of the angular distribution of high energy muons at large altitudes implies that any experiment having a narrow acceptance cone centered at the vertical only accumulates a very small fraction of the muon flux.

### IV.2 Comparison with experimental data

We have compared the fluxes obtained by Monte Carlo simulations with the experimental data available at every altitude. Figures 13, 14 and 15 summarize the results obtained for the ascent, float and ground measurement phases, respectively. In these figures the fluxes of negative and positive muons are plotted as functions of the muon momentum for every altitude considered. In figure 13, that corresponds to the ascent phase of the CAPRICE98 flight, the fluxes corresponding to different altitudes have been multiplied by powers of 100, as indicated in the respective figure caption. The full lines correspond to the AIRES simulations; while the solid squares represent CAPRICE98 data. When available, CAPRICE94 data boe99b ; muonpaper  have been displayed as well (open triangles). Notice that the CAPRICE94 flight corresponds to a different geographical location, and a different time. Therefore the corresponding results, especially for low muon momenta, are not directly comparable to the simulations, that were performed taking into account the CAPRICE98 environment.

The shaded bands drawn together with each distribution illustrate how the flux changes inside the atmospheric depth intervals of table 1. The varying width of these shaded areas can be understood taking into account the behavior of the total muon flux for varying atmospheric depth, represented in figure 2. Near the top of the atmosphere, the muon flux grows with $`X`$, and a visible difference between fluxes evaluated at the minimum and maximum depth of each measurement interval exists as made evident for all momenta in the highest levels plotted in figure 13. In these cases the lower (upper) curve of the shaded region corresponds to the respective minimum (maximum).

On the other hand, in the cases of measurement intervals located close to the ground level, the muon flux decreases with $`X`$, and therefore the positions of the distributions corresponding to the minimum and maximum depth of the respective interval interchange with respect to the highest ones discussed in the previous paragraph. In these cases there are no important variations at the high energy end of such distributions.

For $`X`$ ranging from about 70 g/cm2 to about 200 g/cm2 the muon flux remains approximately constant. In this case the curves corresponding to minimum, maximum and average depth for a given level overlap. As a consequence, any errors due to uncertainties in the altitude of the balloon and/or the atmospheric model used in the simulations should be very small at these intermediate depths.

The relative differences between simulated and experimental data, are shown in figure 16 as function of the atmospheric depth. For each atmospheric depth, the relative difference is a weighted average over momenta, where the weight is the inverse of the relative error on each measured point. The error bars have been calculated taking into account both the experimental and Monte Carlo errors. To estimate these last ones we have taken into account only the uncertainty derived from the uncertainty in the flux at TOA, whose average over the energies with measurements is around 10%. The Monte Carlo statistics is very large so statistical fluctuations of mean values can be neglected. We have not attempted in this work to do a detailed analysis of other uncertainties that can affect the simulated data.

The relative differences are generally compatible with zero within one standard deviation of the experimental flux. At float altitude the simulated flux for $`\mu ^+`$ is 38 $`\pm `$ 18 % above the measured one. At ground level the simulated flux is approximately 20 % lower than the experimental one.

The global averages of the relative differences, represented in figure 16 as dotted lines, are essentially zero (there is a slight positive tendency of about 5 % in the $`\mu ^{}`$ case, that is compatible with zero well within $`1\sigma `$).

The ground level data is worth of a special analysis because in this case the flux was measured with high statistics, allowing for narrower momentum bins. The comparison between simulated and experimental data indicate that the simulations predict a flux that is smaller than the experimental one, especially for high muon momenta, even if such difference is compatible with zero within 2$`\sigma `$. There are several possible reasons for these differences, as already mentioned in the introduction. We discuss here some of our results.

The differences between simulated and experimental data at ground, amounting to about 20%, suggest that there are systematic errors on that level since the statistical errors are very small. The relative importance of primaries with energies above 100 GeV (see figure 4) is larger for ground data than for data at higher altitudes, suggesting that primary energy dependent systematic errors could have a non negligible effect on the muon data at ground. This will be the subject of a detailed study to appear later.

The study of the $`\mu ^+/\mu ^{}`$ flux ratio can also provide important information to test the models used in the simulations. We consider only the data at ground level where contamination corrections are negligible, and the experimental data is abundant enough to give small errors. In figure 17 the flux ratio at ground is plotted versus muon momentum. The squares (triangles) correspond to the CAPRICE98 (CAPRICE97) experimental data, while the shaded area represents the simulation results including an approximate estimation of modeling errors coming from different uncertainties (geomagnetic effects, input flux, hadronic models, etc.). It is worthwhile mentioning that the statistics of the CAPRICE97 data set is about four times larger than that of the CAPRICE98 case, resulting, as is evident from the plot, in smaller error bars. When comparing the results of the simulations with the experimental data, we find in general an acceptable agreement, with relative differences always located within 2$`\sigma `$. However, in the approximate range 3-40 GeV/c the simulations give ratios that are about 10% lower than the experimental results. Such differences are of the same order of magnitude than the ones reported in a similar study performed using various hadronic models wentzicrc , giving a qualitative indication of how much these models can change the simulated $`\mu ^+/\mu ^{}`$ flux ratio. A complete study of the impact of different hadronic models, as well as other factors that could influence the final simulated ratio is beyond the scope of this work, and will be addressed in a future publication.

It is also interesting to study the change with altitude of the flux of muons with a given momentum. In figure 18 the absolute fluxes of muons corresponding to each of the experimental momentum bins are plotted as functions of the atmospheric depth.

Those curves indicate that there exists a very good agreement between our simulations with AIRES and the experimental data in all the considered cases, with no evident bias of any kind. We also note that at each momentum bin both experimental and simulated curves have their maxima approximately at the same altitude.

## V Conclusions

An in depth analysis of the absolute flux of muons at different atmospheric depths has been performed. This includes a comparative study of fluxes measured by the CAPRICE98 and CAPRICE94 experiment with the corresponding simulations with AIRES. Additionally, the most important correlations between the simulated observables and different input parameters are discussed.

Many of the results presented in the previous sections indicate clearly that the strategy followed to perform the simulations is capable of producing acceptable results. We could also verify that the following conditions cannot be simplified without altering to some extent the simulated fluxes:

* Input fluxes that include H, He, and other heavy nuclei, with energies up to a maximum value of at least 10 TeV.
* Full 3d propagation of particles.
* Wide range of primary zenith angles, especially for simulation of full fluxes.

The fluxes simulated with AIRES present a very good agreement with the corresponding ones measured experimentally, at all altitudes and muon momenta considered. With the exception of a few isolated points, the relative differences between simulated and experimental data are always smaller than 20 %. In most cases the relative differences are compatible with zero within error bars. Global averages of relative differences are technically zero, and there are no evident bias that could indicate that the simulation algorithms are not performing adequately.

Detailed comparisons between measured and simulated atmospheric muon data are surely useful to tune the simulation program but requires experimental measurements substantially more accurate than the ones available at present. In this direction, the CAPRICE re-flight project circella  constitutes a very important effort capable of providing improved measurements that will allow the production of more detailed comparative analysis of real and simulated data.

## VI Acknowledgments

One of the authors (S. J. S.) is indebted to Prof. Per Carlson and colleagues at the Department of Physics, KTH, Sweden, for their kind hospitality. This work was partially supported by CONICET and Fundación Antorchas of Argentina. P. H. thanks Roberto Liotta for useful discussion and the Swedish Foundation for International Cooperation in Research and Higher Education for support. E. M. was supported by the Foundation BLANCEFLOR Boncompagni–Ludovisi, née Bildt.


format: text

text:
# Fits to data on polarised structure functions and spin asymmetries with power law corrections

## Abstract

We have compared polarized parton densities determined in the NLO QCD fits to polarized structure functions and spin asymmetries. We consider models of such distributions based on MRST 99 and MRST 2001 fits to non-polarized data. Simple power law corrections corresponding to higher twists are taken into account and their importance is analyzed. The role of positivity conditions for parton densities and their influence on the values of $`\chi ^2`$ is discussed.

PACS numbers: 12.38.-t, 13.88.+e, 14.20.Dh

In  we got parton densities using NLO QCD fit from the data on polarized structure function $`g_1`$ and spin asymmetries. We have used in our fits functional form for parton densities corresponding to unpolarized distributions obtained in the Martin, Roberts, Stirling and Thorne (MRST 98) fit . The results for integrated parton densities were very similar in both cases: for ones gotten from the polarized structure functions and spin asymmetries. In that comparison we have not taken into account power law corrections connected with higher twists contributions. It has been pointed out  that these corrections can play a role in the determination of parton densities. It is an aim of the present paper to consider the influence of the power law corrections in the different methods of extracting polarized parton distributions from fits to the experimental data.

As before we will use all available experimental data on polarized structure functions and corresponding spin asymmetries . Functional form of our polarized parton densities will be taken from corresponding unpolarized parton densities from the fits made by Martin, Roberts, Stirling and Thorne (called MRST 99  and MRST 2001 ). We will use the same method of extracting such densities as in our previous papers , where the experimental data for spin asymmetries were used to obtain polarized parton densities. Many experimental groups gave also experimental results for polarized structure functions $`g_1`$ measured for proton, neutron and deuteron. The asymptotic behaviour of our polarized parton distributions is determined (up to the condition that the corresponding spin densities are integrable) by the fit to unpolarized data. Experiments on unpolarized targets provide information on the spin averaged quark and gluon densities $`q(x,Q^2)`$ and $`G(x,Q^2)`$ inside the nucleon.

Our quark distributions (at $`Q^2=1\mathrm{GeV}^2`$) are parametrized as follows:

$$\mathrm{\Delta }q_i(x)=A_ix^{\lambda _i}(1x)^{\eta _i}(1+ϵ_i\sqrt{x}+\mu _ix).$$
(1)
The similar form one has for gluon density:

$$\mathrm{\Delta }G(x)=B_Gx^{\lambda _G}(1x)^{\eta _G}(1+ϵ_G\sqrt{x}+\mu _Gx).$$
(2)
The values of constants $`\lambda _i`$ and $`\eta _i`$ are given in .

We have for total quark distributions:

$`\mathrm{\Delta }u`$ $`=`$ $`\mathrm{\Delta }u_v+2\mathrm{\Delta }\overline{u},`$
$`\mathrm{\Delta }d`$ $`=`$ $`\mathrm{\Delta }d_v+2\mathrm{\Delta }\overline{d},`$ (3)
$`\mathrm{\Delta }s`$ $`=`$ $`2\mathrm{\Delta }\overline{s},`$
and for axial charges:

$`\mathrm{\Delta }\mathrm{\Sigma }`$ $`=`$ $`\mathrm{\Delta }u+\mathrm{\Delta }d+\mathrm{\Delta }s,`$
$`a_8`$ $`=`$ $`\mathrm{\Delta }u+\mathrm{\Delta }d2\mathrm{\Delta }s,`$ (4)
$`a_3`$ $``$ $`g_A=\mathrm{\Delta }u\mathrm{\Delta }d.`$
In order to determine the unknown parameters in the expressions for polarized quark and gluon distributions we calculate the spin asymmetries (starting from initial value $`Q^2`$ = 1 $`\text{GeV}^2`$) for measured values of $`Q^2`$ and make a fit to the experimental data on spin asymmetries for proton, neutron and deuteron targets. The spin asymmetry $`A_1(x,Q^2)`$ can be expressed via the polarized structure function $`g_1(x,Q^2)`$ as:

$$A_1(x,Q^2)\frac{(1+\gamma ^2)g_1(x,Q^2)}{F_1(x,Q^2)}=\frac{g_1(x,Q^2)}{F_2(x,Q^2)}[2x(1+R(x,Q^2))],$$
(5)
(In the case of SLAC data  $`g_1/F_1`$ values were given experimentally) where $`R=[F_2(1+\gamma ^2)2xF_1]/2xF_1`$, whereas $`\gamma =2Mx/Q`$ ($`M`$ stands for proton mass). We will take the value of $`R`$ from the  (of course formula for $`R`$ determined experimentally already includes corrections from higher twists). In calculating $`g_1(x,Q^2)`$ and $`F_2(x,Q^2)`$ in the next to leading order we use procedure described in  (power law corrections are not included). Having calculated the asymmetries according to eq.(5) for the value of $`Q^2`$ obtained in experiments we can make a fit to asymmetries on proton, neutron and deuteron targets.

As usual in our previous papers the value of $`a_3`$ is not constrained in such fit. We will also do not fix $`a_8`$ but we put a constraint on its value. Simply we will add it as an extra experimental point (from hyperon decays one has $`a_8=0.58\pm 0.03`$, but we enhance an error to 3$`\sigma `$, i.e. to 0.1). It is also possible to use experimentally determined $`F_2(x,Q^2)`$, for example by SMC group from CERN, where the power law corrections are included in the fit . It is also possible to use directly experimental results for polarized structure functions $`g_1`$ given by many experimental groups (the problem is that $`g_1`$ could be determined in a different way in different experiments). We will try to use these three methods of obtaining polarized quark and gluon densities taking into account very simple ($`h/Q^2`$) power law (with no $`x`$ dependence) correction to the NLO QCD expression for the polarized structure function $`g_1(x,Q^2)`$ ($`g_1(x,Q^2)g_1(x,Q^2)+h/Q^2`$). Simple $`x`$ dependence of coefficient $`h(x)`$, as well as terrms of order $`1/Q^4`$, give negligible corrections to $`\chi ^2`$.

Table 1

| $`fit`$ | $`\mathrm{\Delta }u`$ | $`\mathrm{\Delta }d`$ | $`\mathrm{\Delta }s`$ | $`\mathrm{\Delta }\mathrm{\Sigma }`$ | $`g_A`$ | $`\mathrm{\Delta }G`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`A99`$ | $`0.74`$ | $`0.48`$ | $`0.16`$ | $`0.10`$ | $`1.22`$ | $`0.47`$ |
| $`B99`$ | $`0.73`$ | $`0.49`$ | $`0.17`$ | $`0.08`$ | $`1.22`$ | $`0.49`$ |
| $`C99`$ | $`0.78`$ | $`0.48`$ | $`0.14`$ | $`0.16`$ | $`1.26`$ | $`0.43`$ |
|  | $`h_{1p}`$ | $`h_{1n}`$ | $`\chi ^2`$ | $`\chi _0^2`$ | $`\chi _{(np)}^2`$ | $`\chi _{0(np)}^2`$ |
| $`A99`$ | $`0.04`$ | $`0.02`$ | $`165.04`$ | $`177.58`$ | $`153.83`$ | $`156.45`$ |
| $`B99`$ | $`0.03`$ | $`0.03`$ | $`159.35`$ | $`165.94`$ | $`152.01`$ | $`154.02`$ |
| $`C99`$ | $`0.0008`$ | $`0.008`$ | $`150.77`$ | $`150.86`$ | $`149.34`$ | $`149.58`$ |

In the Table 1 we present our results for integrated parton densities with corresponding $`\chi ^2`$ for the fit (subscript $`0`$ corresponds to the fit with no higher twists corrections, whereas $`(np)`$ stands for the fit with no positivity conditions taken into account). Power law corrections are taken in the simple form (suggested in ) $`h_{1p}/Q^2`$ for $`g_1^p`$ and $`h_{1n}/Q^2`$ for $`g_1^n`$ (and corresponding combination for $`g_1^d`$). Fit A99 corresponds to the fit to experimental data on polarized structure functions, B99 to $`g_1`$ but with $`F_2`$ taken from formula from  and C99 to the fit to spin asymmetries where power law corrections are negligible. In the Table 1 we also give $`\chi ^2`$ for corresponding fits without power law corrections. By comparing corresponding columns we see how the introduction of these corrections reduces $`\chi ^2`$. It is seen from Table 1 that integrated parton densities do not differ much for the three fits (It is also true in the case of fits without power law corrections) however there is essential increase in first two models relatively to the third. It seems that the increase in $`\chi ^2`$ of fits A99 and B99 in comparison with fit C99 is connected with positivity conditions assumed for the parton distributions.

To show the influence of positivity conditions we present also the values for $`\chi ^2`$ for the fits without positivity conditions taken into account with power corrections and without such corrections. The values of integrated quark densities do not change significantly (the strongest change is in gluon contributions which become positive) in the fits without positivity conditions comparing with those where positivity conditions for parton densities were assumed. From Table 1 we see that in the models without positivity conditions the effect of taking into account power low corrections is rather small, i.e. 2.6 in the first ($`h_{1p}=0.02\pm 0.01`$, $`h_{1n}=0.01\pm 0.04`$) and 2 in the second model ($`h_{1p}=0.02\pm 0.02`$, $`h_{1n}=0.03\pm 0.04`$) for 2 degrees of freedom. That means that in those fits the parameters $`h_{1p}`$ and $`h_{1n}`$ are at the border of being relevant. It is not very surprising that the fits with assumed positivity conditions (that means with additional restrictions on the form of the fitted parton densities) give $`\chi ^2`$ higher then fits where the positivity conditions are not assumed. On the other hand we probably should not expect in NLO QCD very strong violation of positivity conditions. The values for the integrated parton densities given in Table 1 for the third model can be compared with the values obtained using MRST 98 model:

$`\mathrm{\Delta }u`$ $`=`$ $`0.77,`$
$`\mathrm{\Delta }d`$ $`=`$ $`0.59,`$
$`\mathrm{\Delta }s`$ $`=`$ $`0.20,`$
$`\mathrm{\Delta }\mathrm{\Sigma }`$ $`=`$ $`0.02,`$ (6)
$`g_A`$ $`=`$ $`1.36`$
$`\mathrm{\Delta }G`$ $`=`$ $`0.01`$
$`\chi ^2`$ $`=`$ $`150.47.`$
Our fit to spin asymmetries where 431 points were considered (without averaging over $`Q^2`$) was discussed in  and similar values for integrated parton densities were obtained.

In Table 2 we present in a similar way as in Table 1 results obtained from the fit corresponding to the fit MRST 2001.

Table 2

| $`fit`$ | $`\mathrm{\Delta }u`$ | $`\mathrm{\Delta }d`$ | $`\mathrm{\Delta }s`$ | $`\mathrm{\Delta }\mathrm{\Sigma }`$ | $`g_A`$ | $`\mathrm{\Delta }G`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`A01`$ | $`0.79`$ | $`0.32`$ | $`0.04`$ | $`0.47`$ | $`1.14`$ | $`56.6`$ |
| $`B01`$ | $`0.81`$ | $`0.35`$ | $`0.06`$ | $`0.40`$ | $`1.16`$ | $`42.6`$ |
| $`C01`$ | $`0.86`$ | $`0.38`$ | $`0.05`$ | $`0.44`$ | $`1.24`$ | $`29.9`$ |
|  | $`h_{1p}`$ | $`h_{1n}`$ | $`\chi ^2`$ | $`\chi _0^2`$ | $`\chi _{(np)}^2`$ | $`\chi _{0(np)}^2`$ |
| $`A01`$ | $`0.06`$ | $`0.00`$ | $`199.73`$ | $`217.44`$ | $`155.10`$ | $`158.51`$ |
| $`B01`$ | $`0.05`$ | $`0.02`$ | $`188.88`$ | $`197.73`$ | $`152.77`$ | $`154.41`$ |
| $`C01`$ | $`0.0005`$ | $`0.002`$ | $`158.40`$ | $`158.41`$ | $`148.14`$ | $`148.23`$ |

The values of $`\chi ^2`$ corresponding to the model with positivity condition for parton densities in this case are much higher then for the solution corresponding to MRST 99 presented in Table 1. For the spin asymmetries one has $`\chi ^2`$=158.41, which have to be compared with 150.75. These solutions have more singular behaviour for valence $`u`$ quark and very different gluon behaviour for small $`x`$. Therefore it is more difficult to fit (with positivity conditions for quark densities) the experimental data and $`\chi ^2`$ for fits A01 and B01 and C01 is higher. On the other hand the solutions without assuming positivity conditions are not very different from the solutions corresponding to MRST 99 ($`\chi ^2`$ corresponding to spin asymmetries 148.23 is smaller then 149.58). The reduction in $`\chi ^2`$ coming from power law corrections for solutions with assumed positivity conditions for quark densities is significant (for B01 $`h_{1p}=0.05\pm 0.02`$, $`h_{1n}=0.02\pm 0.04`$) but similarly to the situation presented in Table 1 for the solutions without assuming positivity conditions the changes in $`\chi ^2`$ caused by these corrections are rather small (for A01$`(np)`$ $`h_{1p}=0.03\pm 0.02`$, $`h_{1n}=0.01\pm 0.04`$ and for B01$`(np)`$ $`h_{1p}=0.02\pm 0.02`$, $`h_{1n}=0.01\pm 0.04`$). It is difficult to draw a conclusion that taking into account our simple power law corrections of the form $`h/Q^2`$ is very important for the fits. The integrated quark densities for different fits do not differ much (like in Table 1 there are some changes for integrated glouon distributions). In general taking into account positivity conditions for parton densities, in spite of relatively big changes in $`\chi ^2`$, does not change significantly the integrated parton densities.

There are changes when we compare solutions corresponding to MRST 98, MRST 99 and MRST 2001. From the eq., Table 1 and Table 2 we see that $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ increase, $`\mathrm{\Delta }s`$ decrease and as consequence $`\mathrm{\Delta }\mathrm{\Sigma }`$ increase. The values of $`\mathrm{\Delta }\mathrm{\Sigma }`$ are not very different from $`g_8=0.58`$. We have the situation when the solutions give different integrated parton densities and the corresponding $`\chi ^2`$ values are very close (especially when we take solutions corresponding the case when we do not assume positivity conditions for parton densities). There is a dependence on the form of assumed parton densities but that does not influence strongly $`\chi ^2`$ values. It is tempting to choose (in spite of higher $`\chi ^2`$ value when positivity condition for quark densities are assumed) solution corresponding to MRST 2001 (recent fit to unpolarized data) where integrated parton densities give relatively high value of $`\mathrm{\Delta }\mathrm{\Sigma }`$ and small value of $`\mathrm{\Delta }s`$ (or at least consider it seriously).

We have compared polarized parton distributions, corresponding to MRST 99 and MRST 2001 determination, using three different methods of obtaining polarized parton densities. Simplest power law correction corresponding to higher twists have been taken into account. As expected these power law corrections are negligible in the fits to spin asymmetries. At the first sight they seem to be important in the first two models. For comparison we have also considered the solutions where positivity conditions for parton densities have not been assumed. We observe strong increase in $`\chi ^2`$ in models A99, A01 and B99, B01 (fits to $`g_1`$) connected with positivity conditions. Reduction in $`\chi ^2`$ by taking into account power law corrections is significant in the models where positivity conditions were assumed and rather marginal in models without positivity conditions for parton densities. Hence, it seems that the determination of power law corrections is not very reliable. There are some differences in integrated parton densities corresponding to the models MRST98, MRST 99 and MRST 2001. Within the definite model three different methods of fitting give very similar results for integrated quark densities (but not very much for a gluon one) also in the case when positivity conditions for parton densities are not assumed. The latest model (C01) even if not by $`\chi ^2`$ (with slightly higher $`\chi ^2`$ value with positivity conditions for parton densities assumed) values is preferred by interpretation reasons.


text:
# First-principles thermodynamics of transition metals: W, NiAl, PdTi

## I Introduction

First-principles calculations have established an excellent record for describing the ground state properties of materials at zero temperature over almost twenty years. More recently, the full description of the electronic structure allows for accurate calculation of phonon frequenciesbaroni  and the thermodynamic properties derivable from themfreeE . These methods have been applied to insulating systems and elemental metals. In this paper, we extend the methodology to investigate transition metals and compounds that have a low symmetry or phonon anomalies. Previous authors have looked at elemental metals Al, Li, Naquong , Agxie , Cunarasimhan , and Wdebernardi  and found that the quasiharmonic approximation holds most of the way to melting. We extend this work to alloys, which introduces some new features. As our prototype materials, we study W, NiAl and PdTi, which have many similarities but each of which has aspects of peculiar interest.

NiAl is a high melting point alloy which finds applications in aerospace and turbine blade design. Its electronic structure is characterized by Kohn anomalies in the phonon spectrumhuang:nial , arising from a Fermi surface nesting effect. A huge amount of work has been done on the system using empirical potentialsackland:mrs ; gao ; mishin:2002 ; mishin:1997 ; yan , which cannot describe the Kohn anomaly, and it is interesting to examine whether this neglect has serious structural consequences. Crystallographically, NiAl has the cubic $`B2`$ (CsCl) crystal structure which is stable at all temperatures. This high symmetry simplifies the electronic structure calculation enormously.

PdTi provides a sharp contrast to NiAl. It is a high temperature shape memory alloy, isostructural with NiAl at high temperatures in the $`B2`$ phase, but undergoing a martensitic phase transition to the tetragonal $`B19`$ phase at 810 K and predicted to undergo a further transition to the monoclinic $`B19^{}`$ at very low temperatureshuang:2002b . The dynamically stabilized high temperature phase and low symmetry low temperature phase each provide a challenging test for theoretical prediction.

Tungsten is another high melting point transition metal. In addition to studying its own thermal expansion, we use this to investigate the effects of anisotropic strain on the phonon contribution to the free energy. This enables us to make a simple test of some approximations required in the calculation of PdTi.

The difference between the materials being studied is highlighted by the phonon spectra for the $`B2`$ phase: NiAl and PdTi both adopt the $`B2`$ structure at high temperature while tungsten is a refractory metal and has a bcc structure at all temperature. The bcc structure can be regarded as a special case of $`B2`$ where both atoms are the same. NiAl has phonon anomalies due to Fermi surface nesting, while $`B2`$ structure PdTi is dynamically stabilized: it has negative frequency phonons (i.e. mechanically unstable) at low temperatures. Such a situation leads to temperature dependent modes, which cannot be treated within the quasiharmonic approximation and would need to be dealt with separately either in realzhong  or reciprocaldrummond  space.

The quasiharmonic approximationquasi  to the free energy assumes all phonons can be treated as simple harmonic oscillators with a frequency dependent on the volume of the material. The static lattice energy ($`T=0`$) is evaluated at a range of volumes, and for each of these the phonon dispersion relation is calculated. The quasiharmonic contribution to the free energy at finite temperature is then evaluated from Bose-Einstein statistics of the phonons at a given volume. From this the pressure can be evaluated for all conditions of $`T`$ and $`V`$, giving an equation of state. Alternately in the zero pressure case, the relationship between $`T`$ and $`V`$ (i.e. the thermal expansion) can be determined. Since larger volumes typically lead to lower phonon frequencies, the entropy is typically lowered by expansion. Thus larger volume (entropy) is favoured at higher temperatures, and thermal expansion occurs.

There are, of course, other contributions to thermal expansion in ordered alloys (antisite defect creation, vacancy creation, coupling of phonons to lattice parameters, anharmonic effects due to finite atomic displacements) and it is interesting to test how much of the observed expansion can be attributed to quasiharmonic effects.

## II Theory

Phonon dispersions can be evaluated using density-functional perturbation theory (DFPT) baroni:1987  and finite displacement methodsprac . Here we use the former. Previous work on isotropic equations of state for cubic materials covers much of the theory, here we investigate the practical requirements for lower symmetry materials, in particular what approximations can be made to made calculations tractable.

Virtually all implementations of density-functional theory (DFT) based ab initio codes take as input the unit cell vectors a, b and c and atomic positions. We define a matrix B comprised of the three unit cell vectors.

$$B=(\text{a, b, c})$$

The volume of the unit cell is

$$V=|B|=(\text{a}\times \text{b})\text{c}.$$

For a cubic material, the quasiharmonic approximation uses the phonon frequencies calculated at different volume to obtain the free energy as a function of two variables,

$$F(V,T)=E(V)k_BT\mathrm{ln}Z(V)$$
(1)
where $`E(V)`$ is the cold curve energy at the volume $`V`$ (i.e. the energy of a system with the atoms placed on their lattice sites), $`k_B`$ is the Boltzmann’s constant and $`Z`$ is the vibrational partition function. We assume that for all temperatures of interest the contribution of electronic excitations is negligible - generally true for equation of statesdswift , but not for transport properties such as heat conduction.

From $`F(V,T)`$ one can obtain the equation of state of the single phase of the system:

$$P=\left(\frac{\delta F}{\delta V}\right)_T$$
(2)
and the volume thermal expansion:

$$\alpha _V=\frac{1}{V}\left(\frac{\delta V}{\delta T}\right)_P$$
(3)
These can be evaluated for a given structure regardless of whether an alternate phase of lower Gibbs free energy exists.

For non-cubic materials, the free energy is a function of all independent cell parameters, and the linear thermal expansion is a tensor quantity. Calculating the free energy across such a multidimensional space is impractical, so we examine a further approximation: that the effect of increasing temperature is equivalent to a negative pressurefoot  At each pressure, the internal coordinates must be re-relaxed to their equilibrium position within the unit cell, since the phonon spectrum is only well defined by an expansion about such an equilibrium.

Thus, we proceed as follows:

1/ Evaluate the equilibrium structure at $`P=0`$, by minimizing $`U(𝐚,𝐛,𝐜,\{u_i\})U_0`$ where $`u_i`$ are the internal atomic coordinates.

2/ Evaluate the force constants and hence the phonon dispersion relation for the relaxed structure, and the partition function $`Z(V)`$ for these oscillators. $`V=|B|=(\text{a}\times \text{b})\text{c}`$.

3/ Evaluate Gibbs free energy at a range of temperatures $`F(V,T)=U_0kT\mathrm{ln}Z+PV`$.

4/ Evaluate the equilibrium structure at P, by minimizing enthalpy $`U(𝐚,𝐛,𝐜,\{u_i\})+PVH_0(V)`$ where $`u_i`$ are the internal atomic coordinates.

5/ Evaluate the phonon dispersion relation for this relaxed structure, and the partition function $`Z`$ for these oscillators.

6/ Evaluate Gibbs free energy at fixed V for a range of temperatures $`F(V,T)=U_0kT\mathrm{ln}Z+PV`$.

7/ Repeat 4, 5, 6 for a range of positive and negative pressures.

Once a coarse grid of $`V`$-points has been generated (for thermal expansion, as few as 3 different pressures suffices) a finer grid can be generated without further expensive ab initio calculation by interpolating force constantsdswift . A dense grid of $`T`$-points can easily be generated at a given $`V`$.

In terms of speed and accuracy there is little to choose between DFPT or finite displacements: the former involves calculation of second derivatives and the later requires large supercells. We note that for finite displacements, once the eigenvectors are calculated at $`P=0`$, a single calculation of restoring forces with all atoms displaced suffices at other volumesfreeE .

The allowed occupation states of each phonon mode of frequency $`\omega `$ are $`ϵ_n(\omega )=(n+1/2)\mathrm{}\omega `$ where $`n`$ is the the number of phonons populating that mode. The canonical partition function is therefore

$$\xi (\omega )\underset{n=0}{\overset{\mathrm{}}{}}e^{\beta ϵ_n(\omega )}=1/(2\mathrm{sinh}(\mathrm{}\omega /2k_BT)).$$
(4)
The total partition function is then:

$$Z=\mathrm{\Pi }_\omega \xi (\omega )$$
(5)
For calculating thermal expansion under constant (zero) pressure boundary condition, the essential quantity is the free energy. The mean vibrational free energy:

$$k_BT\mathrm{ln}Z(V)=k_BT_0^{\mathrm{}}g(\omega ,V)\mathrm{ln}(\xi (\omega ,T))𝑑\omega ,$$
(6)
introducing the phonon density of states $`g(\omega ,V)`$ to make the dependence of density of states on volume explicit.

Phonon frequencies are evaluated at a dense set of $`q`$-points in the first Brillouin zone using the force constants matrix obtained by fourier transform of the “exact” DFPT frequencies. The density of phonon states is calculated by integration using the tetrahedron methodlehmann . This provides the $`Z(V)`$ required in step 2/ above.

The structure obtained at a particular temperature can be found by minimizing the free energy with respect to the lattice vectors and internal coordinates of the atoms

$$\frac{dF}{dB_{ij}}=0;\frac{dF}{d𝐮_𝐤}=0$$
(7)
These derivatives can be split into cold curve ($`E(B,𝐮_𝐤)`$) and phonon ($`F_{phon}`$) contributions.

To make calculations for complex structures tractable, we make the assumption that the phonon part of the free energy is a function of volume only, and is independent of $`𝐮_𝐤`$. We examine this approximation in the case of anisotropic deformation of tungsten and find it to be good. Intuitively, this can be understood if the force constants vary linearly with separation: we might expect that volume conserving shears and changes in $`𝐮_𝐤`$ will stiffen some force constants and weaken others, raising some frequencies and lowering others to give an overall cancellation in the small strain limit. By contrast, volume increases weaken all force constants and give a systematic decrease in mean frequency. Hence the conditions above become:

$`{\displaystyle \frac{E}{B_{ij}}}|_{𝐮_𝐤}=0`$ (8)
$`{\displaystyle \frac{F_{phon}}{V}}|_{𝐮_𝐤}=0`$ (9)
$`{\displaystyle \frac{E}{𝐮_𝐤}}|_V=0`$ (10)
## III Results

### III.0.1 Anisotropic strain - Tungsten

We use tungsten as a benchmark to investigate the effects of anisotropic strain on the phonon contribution to the free energy. Tungsten has a bcc crystal structure. The calculations were done with the generalized gradient approximation (GGA) and the equilibrium lattice parameter was computed to be 2.905 Å, comparable to the experimental result of 2.887 Å. Phonon calculations were done at three different volumes: $`V_0`$, 1.01$`{}_{}{}^{3}V_{0}^{}`$ and 1.02$`{}_{}{}^{3}V_{0}^{}`$. At each volume, phonon dispersion was done with the bcc structure and tetragonal strained structures with four different $`c/a`$ ratios, i.e. 0.96, 0.98, 1.02 and 1.04.

We calculated phonon density of states (DOS) at 15 different structures corresponding to tetragonal strains ($`e_{xx}`$; $`e_{yy}=e_{zz}`$; $`e_{xy}=e_{yz}=e_{zx}=0`$). The phonon dispersion of bcc tungsten at equilibrium lattice parameter along with experimental results is shown in Fig. 1, while the phonon dispersions under tetragonal strains and a larger volume are shown in Fig. 2(a) and Fig. 2(b), respectively. The calculated phonon frequencies are in good agreement with experiment (Fig. 1) and with a recent calculationdebernardi . Fig. 2(a) and (b) show that compared to volumetric strain, volume-conserving shear strains have negligible contribution to phonon free energy.

The free energy was evaluated at each strain state. In Fig. 3 we show the thermal expansion calculated using the three isotropic strains. The linear thermal expansion adopts the form from Ref. amer, . The agreement between theoretical and experimental results is very good.

Fig. 4 shows the phonon free energy with different strain. As expected, the phonon contribution to the free energy is only weakly affected by volume-conserving anisotropic strain. Thus it is a good approximation to neglect this term: a given volume-conserving strain contributes about 5% as much as equivalent volume strain.

### III.0.2 Anomalous phonons - NiAl

The phonon spectrum of $`B2`$ NiAl is characterized by Kohn anomalies. As a consequence, to describe the phonon dispersion properly requires a thorough DFPT $`q`$-point sampling or (equivalently) large supercells in direct space methods. In Fig. 5, we show the convergence of the result to experimental observation with $`q`$-point sampling density.

Despite the anomalies, we find that the thermodynamic properties are almost independent of $`q`$-point sampling aside from a rigid shift of the free energy. This does not affect the thermal expansion, but stabilizes the structure slightly against alternate structures because of the increased entropy of the low frequency phonons.

The entire calculated phonon dispersion relations at two different volumes are shown in Fig. 6(a). All phonon frequencies increase under compression. Hence the mode Grüneisen parameters, $`\gamma _j(\text{q})`$, shown in Fig. 6(b), and defined as

$$\gamma _j(\text{q})=\frac{\omega _j(\text{q})}{V}\frac{V}{\omega _j(\text{q})}.$$

are positive throughout the whole BZ. There is no anomalous thermal expansion. Integration of free energies derived from phonon spectra gives the thermal expansion shown in Fig. 7. The excellent agreement with experimental data suggests that both DFPT phonons and the quasiharmonic approximation are valid in this case.

### III.0.3 Complex structure - PdTi

For PdTi we calculate phonon spectra for two phases, the low temperature $`B19^{}`$ and the ambient temperature $`B19`$ phase. In both cases, the equilibrium cold curve structure is determined for a particular pressure by relaxing internal parameters and lattice constants simultaneouslyhuang:2002b . It is then assumed that the effect of thermal expansion on lattice parameter ratios and internal parameters is equivalent to that of (negative) pressure, such that the phonon spectrum need only be evaluated once at each volume, and the free energy can be interpolated as a function of volume only, with the free parameter relaxation done on the cold-curve structure. The results for anisotropic deformation in tungsten (Fig. 4) suggest that even if the effects of pressure and thermal expansion on free parameters are not equivalent, the consequent change in thermal free energy due to non-volumetric strains will be small.

The $`B19^{}`$ phase has stable phonons throughout the entire BZ (Fig. 8(a)) but the $`B19`$ phase (Fig. 8(b)) has some imaginary frequency modes: we assume these contribute to the free energy like free particlesdswift . Assuming the phonons in the $`B19^{}`$ phase to remain harmonic up to and above the phase transition, this enables us to calculate the relative free energies of the two phases, and determine for each temperature not only the lattice parameters but also which is the stable phase. This gives rise to the thermal expansion shown in Fig. 9 with a discontinuity marking the phase transition.

The calculated phase transition temperature for the $`B19`$-$`B19^{}`$ phase, on the basis of equal quasiharmonic free energies, is 140K. This can be compared with another estimate based on treating the $`B19`$ phase as a barrier between equivalent $`B19^{}`$ variants, which barrier (about 0.0007eV/atom) can continually be crossed when the temperature reaches $`\mathrm{\Delta }E/k_B=9K`$ in the higher temperature phase. The discrepancy between these arises because one concentrates on dynamics of a particular soft mode, while the other considers the static average over all modes: experimentally, the lower of the two is expected.

## IV Discussion

We have calculated the thermal expansion properties of W, NiAl and PdTi using the quasiharmonic approximation and DFPT finding that agreement with experiment is excellent. This shows that the harmonic phonon free energy is by far the dominant effect in thermal expansion, and that other effects such as thermal defects, electron free energy, coupling of phonons to lattice parameters, and anharmonic effects due to finite atomic displacements can safely be neglected.

These three materials are chosen to represent represent an increasingly challenging test to this methodology. The accurate thermal expansion in tungsten was expected on the basis of previous workquong ; xie , and allowed for the demonstration that anisotropic strain has only a second order effect on the free energy. Similarly the results on NiAl show that the method can be applied to alloys just as effectively as with elements and showed that the presence of phonon anomalies also has little effect on thermal expansion.

Without these two results, the PdTi calculation would have been intractable, but armed with the knowledge that the anisotropic temperature effect can be treated as a pressure effect on the cold curve, allowing the vibrational partition function to be treated as a function of volume only, calculation and minimizations of the full $`F(B_{ij},𝐮_𝐤,T)`$ becomes tractable. In addition, implementing the previously-described treatment of soft phononsdswift ; drummond  allow us to predict the thermal expansion of both $`B19^{}`$ and the dynamically stabilized $`B19`$ PdTi crystal, in addition to estimating the phase transition temperature for $`B19`$-$`B19^{}`$. To our knowledge, these quantities have yet to be measured experimentally and such measurement will provide a sensitive test to our methods.

In sum, this paper represents a significant theoretical advance in the type of materials whose thermal expansion can be calculated from ab initio simulation.

This work was carried out under grant AFOSR/MURI F49620-98-1-0433. The calculations were performed on the SGI Origin 2k/3k at ARL MSRC. GJA would also like to thank the Fulbright Foundation for support.


text:
# A Possible Mechanism of Biological Memories in terms of Quantum Fluids

## Abstract

A mechanism of memories, especially biological memories, is studied in terms of quantum fluids. Two-dimensional flows in central potentials $`V_a(\rho )=a^2g_a\rho ^{2(a1)}`$ ($`a0`$ and $`\rho =\sqrt{x^2+y^2}`$) have zero-energy eigenstates that degenerate infinitely for all $`a`$. It is shown that stable standing waves constructed from the zero-energy flows are confined in various types of polygons which can be the minimum units of memory systems. Vortex patterns awoken in the units by stimuli correspond to the memories of the stimuli. This memory system is not a system for preserving memories as usual but that for awaking memories. The system has interesting properties; (i) the absolute economy as for the energy consumption, (ii) the infinite variety for a huge number of memories, (iii) the perfect recovery of the system from any disturbances by stimuli, and (iv) the large flexibility in the construction of the system. A process for thinking is also proposed in terms of this memory system.

1. Introduction

Researches for the mechanism of brains are extensively going on from many different experimental and theoretical standpoints, and many interesting results have been reported. At present, however, these approaches are still not enough to understand the mechanism of brains from the viewpoint of a fundamental theory of physics. The first step of the understanding seems to be the elucidation of the mechanism of memories in brains as a physical process. The insides of our bodies are mainly occupied by many kinds of fluids, and the total unification of many organs are also preserved mainly by such fluids. It is not a peculiar idea that the mechanism in brains is governed by some kind of flows of such fluids. In this sense it is quite natural that the dynamics of flows (hydrodynamics) is one of the fundamental dynamics in our bodies. On the other hand the dynamics in brains does not seem to be a completely deterministic dynamics like classical dynamics, but some kind of probability seems to play an important role in the mechanism of thinking. Considering the coexistence of the probability and the fundamental principle of superposition in wave dynamics, we may invoke some kind of quantum process related to hydrodynamics to describe the dynamics in brains. Such a hydrodynamical approach of quantum mechanics was vigorously examined in the early stage of the development of quantum mechanics \[1-8\], and then some fundamental properties of the quantum flows like vortices were investigated \[9-13\]. Recently interesting solutions that have exactly zero-energy eigenvalue and degenerate infinitely have been found in two-dimensional quantum flows \[14-16\]. From the viewpoint of the energy consumption in the construction and the preservation of the mechanism of brains such zero-energy flows can be very useful tools. Furthermore the infinite degeneracy can provide the infinite variety of vortex patterns \[15-17\] corresponding to a huge number of memories. Note that the vortices have recently been observed in many phenomena such as condensed matters \[18-20\], non-neutral plasma \[21-24\] and Bose-Einstein gases \[25-29\]. Here we show a possibility of the mechanism of memories in terms of the zero-energy flows. It is shown that the memory system proposed here is not a system for preserving memories as usual but that for awaking memories, and it has quite suitable properties to represent biological memory systems. 2. A short review on zero-energy flows in two-dimensional Schr$`\ddot{\mathrm{o}}`$dinger equations

In this section we briefly study two-dimensional (2D) quantum flows and vortices. It has been shown that 2D Schr$`\ddot{\mathrm{o}}`$dinger equations with central potentials $`V_a(\rho )=a^2g_a\rho ^{2(a1)}`$ ($`a0`$ and $`\rho =\sqrt{x^2+y^2}`$), of which eigenvalue problems with the energy eigenvalue $``$ are written as

$$[\frac{\mathrm{}^2}{2m}\mathbf{}+V_a(\rho )]\psi (x,y)=\psi (x,y)$$
(1)
where $`\mathbf{}=^2/x^2+^2/y^2,`$ have zero-energy ($`=0`$) eigenstates \[14-16\]. Note here that in this equation the mass $`m`$ and the coordinate vector $`(x,y)`$ can represent not only those of the single particle but those of the centre of mass for a many particle system as well . It has also been shown that, as far as the zero-energy eigenstates ($`\psi _0`$) are concerned, the Schr$`\ddot{\mathrm{o}}`$dinger equations for all $`a`$ can be reduced to the following equation in terms of the conformal transformations $`\zeta _a=z^a`$ with $`z=x+iy`$;

$$[\frac{\mathrm{}^2}{2m}\mathbf{}_ag_a]\psi _0(u_a,v_a)=0,$$
(2)
where $`\mathbf{}_a=^2/u_a^2+^2/v_a^2,`$ using the variables defined by the relation $`\zeta _a=u_a+iv_a`$ . That is to say, the zero-energy eigenstates for all the different numbers of $`a`$ are described by the same plane-wave solutions in the $`\zeta _a`$ space. Furthermore it is easily shown that the zero-energy states degenerate infinitely. Let us consider the case for $`a>0`$ and $`g_a>0`$. Putting the function $`f_n^\pm (u_a;v_a)e^{\pm ik_au_a}`$ with $`k_a=\sqrt{2mg_a}/\mathrm{}`$ into (2), where $`f_n^\pm (u_a;v_a)`$ are polynomials of degree $`n`$ ($`n=0,1,2,\mathrm{}`$), we obtain the equation for the polynomials

$$[\mathbf{}_a\pm 2ik_a\frac{}{u_a}]f_n^\pm (u_a;v_a)=0.$$
(3)
Note that from the above equations we can easily see the relation $`(f_n^{}(u_a;v_a))^{}=f_n^+(u_a;v_a)`$ for all $`a`$ and $`n`$. General forms of the polynomials have been obtained by using the solutions in the $`a=2`$ case (2D parabolic potential barrier (2D PPB)) . Since all the solutions have the factors $`e^{\pm ik_au_a}`$ or $`e^{\pm ik_av_a}`$, we see that the zero-energy states describe stationary flows . Taking account of the direction of incoming flows in the $`\zeta _a`$ plane that is expressed by the angle $`\alpha `$ to the $`u_a`$ axis, the general eigenfunctions with zero-energy are written as arbitrary linear combinations of eigenfunctions included in two infinite series of $`\left\{\psi _{0n}^{\pm (u)}(u_a(\alpha );v_a(\alpha ))\right\}`$ for $`n=0,1,2,\mathrm{}`$, where

$$\psi _{0n}^{\pm (u)}(u_a(\alpha );v_a(\alpha ))=f_n^\pm (u_a(\alpha );v_a(\alpha ))e^{\pm ik_au_a(\alpha )},$$
(4)
$`u_a(\alpha )=u_a\mathrm{cos}\alpha +v_a\mathrm{sin}\alpha `$, and $`0<\alpha 2\pi `$. (For the details, see the sections II and III of Ref. 15.) We briefly comment on states having non-zero energy eigenvalues. Energy eigenvalues are in general complex values, and the states with complex energies can be understood to be resonance states which are not stable and decay in some time scale . Exact solutions not only for zero-energy but for non-zero energies as well have been solved in 2D PPB . Note here that we consider the motions of the $`z`$ direction perpendicular to the $`xy`$ plane to be free motions represented by $`e^{\pm ik_zz}`$. We should, therefore, notice that the total energies $`E_T`$ of the flows are given by $`E_T=E_z`$, where $`E_z`$ are the energies of the free motions in the $`z`$ direction. We shall, however, not take account of the motions in the $`z`$ direction in the following discussions on the motions in the $`\zeta _a`$ plane.

Note that the zero-energy eigenfunctions can not be normalized as same as those in scattering processes . Actually it has been shown that all the states for zero and non-zero energies are eigenstates in the conjugate spaces of the Gel’fand triplets \[13-16\]. This fact means that the probability density ($`\rho (t,x,y)=|\psi (t,x,y)|^2`$) and the probability current ($`𝒋(t,x,y)=\mathrm{Re}[\psi (t,x,y)^{}(i\mathrm{})\psi (t,x,y)]`$) lose the meanings. We, however, see that the velocity defined by $`𝒗=𝒋(t,x,y)/\rho (t,x,y),`$ can have the well-defined meaning, because the ambiguity arising from the normalization is cancelled in the definition. The velocity was widely discussed in early stage of quantum mechanics \[1-8\]. It is also well-known that vortices appear at the zero points of the density, that is, the nodal points of the wave functions \[9-13\]. Let us here remember the fact that the zero-energy states have the infinite degeneracy, and therefore we can construct wave functions having the nodal points at almost arbitrary positions in terms of linear combinations of the zero-energy states. Actually various vortex patterns have been studied \[15-17\]. It should be emphasized that we can study the vortex problems for all $`V_a(\rho )`$ as a problem with the constant potential $`V_a=g_a`$ in the $`\zeta _a`$ plane, because the conformal transformations do not change fundamental properties of vortices such as the numbers of vortices and the strengths of vortices $`\mathrm{\Gamma }`$. In quantum mechanical phenomena it is shown that the circulation $`\mathrm{\Gamma }=_C𝒗𝑑𝐬`$ is quantized such that $`\mathrm{\Gamma }=2\pi l\mathrm{}/m`$, where $`l`$ is an integer \[9-14\]. Vortex patterns in the $`(x,y)`$ plane can be obtained by the inverse transformations of the conformal mappings . 3. Standing waves and vortex patterns

Let us start with the study of zero-energy standing waves that will play an important role in the following discussions. Considering the relation $`(f_n^{}(u_a;v_a))^{}=f_n^+(u_a;v_a)`$, we can easily make two infinite series of standing waves in terms of the zero-energy series $`\left\{\psi _{0n}^{\pm (u)}\right\}`$ as follows;

$`\psi _{0n}^{+\mathrm{S}(u)}`$ $`={\displaystyle \frac{1}{2}}(\psi _{0n}^{+(u)}+\psi _{0n}^{(u)})=f_n^{\mathrm{Re}(u)}\mathrm{cos}k_au_a(\alpha )f_n^{\mathrm{Im}(u)}\mathrm{sin}k_au_a(\alpha ),`$
$`\psi _{0n}^{\mathrm{S}(u)}`$ $`={\displaystyle \frac{1}{2}}(\psi _{0n}^{+(u)}\psi _{0n}^{(u)})=i[f_n^{\mathrm{Re}(u)}\mathrm{sin}k_au_a(\alpha )+f_n^{\mathrm{Im}(u)}\mathrm{cos}k_au_a(\alpha )],`$ (5)
where $`f_n^{\mathrm{Re}(u)}`$ and $`f_n^{\mathrm{Im}(u)}`$, respectively, stand for the real and imaginary parts of $`f_n^+`$. In order to simplify our argument we pick up only four standing waves that are constructed in terms of the linear combinations of the flows with the lowest order, that is, $`e^{\pm i(k_uu_a\pm k_vv_a)}`$ corresponding to the choices of $`\alpha =\theta _a,\pi \pm \theta _a`$ and $`2\pi \theta _a`$, where $`k_u=k_a\mathrm{cos}\theta _a`$ and $`k_v=k_a\mathrm{sin}\theta _a`$ with $`0<\theta _a2\pi `$, such that

$$\psi _{cc}=\mathrm{cos}k_uu_a\mathrm{cos}k_vv_a,\psi _{ss}=\mathrm{sin}k_uu_a\mathrm{sin}k_vv_a,\psi _{sc}=i\mathrm{sin}k_uu_a\mathrm{cos}k_vv_a,\psi _{cs}=i\mathrm{cos}k_uu_a\mathrm{sin}k_vv_a.$$
(6)
These wave functions can be confined in quadrangles of the $`\zeta _a`$ plane, which are surrounded by rigid fences described by infinite potentials, because we can take zero as boundary values on the sides of the quadrangles. For example $`\psi _{cc}=\mathrm{cos}k_uu_a\mathrm{cos}k_vu_a`$ becomes zero on the lines where either $`k_uu_a=(2m+1)\pi /2`$ or $`k_vv_a=(2n+1)\pi /2`$ ($`m`$ and $`n`$ are integers) is satisfied. In the cases where $`a`$ are positive integers the quadrangles in the $`\zeta _a`$ plane are represented by $`4a`$-gons in the original $`xy`$ plane. Now we get the standing wave states confined in the special regions which are described by the quadrangles in the $`\zeta _a`$ plane. It should be noticed that the quadrangles have a wide variety of the sizes which are determined by the numbers of wave lengths contained in the quadrangles. That is to say, the number $`N_u`$ for the $`u`$ direction and that $`N_v`$ for the $`v`$ direction are free parameters, where $`N_u`$ and $`N_v`$ are positive integers or positive half-integers.

Now let us discuss vortex patterns made from the zero-energy flows. It is known that vortex patterns are one of topological properties of flows, and then they can be stable in many physical phenomena. Following the above consideration, we can prepare many different standing waves depending on the sizes of quadrangles. The standing waves, of course, have no vortex, because $`𝒋=0`$. Let us study an example of a quadrangle where a standing wave $`\psi _{cc}`$ with $`N_u=1/2`$ and $`N_v=3/2`$ is confined. We would like to investigate what happens when a flow is put into the quadrangle. Let us put a different standing wave, for example, the standing wave with the wave number vector $`\stackrel{}{q}=(q_u,q_v)`$, which is described by the wave function $`\psi _{sc}=i\mathrm{sin}q_uu_a\mathrm{cos}q_vv_a`$, into the quadrangle prepared. Note that in the process where the new flow is added we can use the free motions of the $`z`$ direction. We see that the ranges of $`u`$ and $`v`$ in the quadrangle are given by $`|u|<u_B\pi /2k_u`$ and $`|v|<v_B3\pi /2k_v`$. (Hereafter we omit the suffix $`a`$ from $`u_a`$, $`v_a`$, $`k_a`$ and etc.) The total wave function is written as

$$\mathrm{\Psi }=\mathrm{cos}k_uu\mathrm{cos}k_vv+iC\mathrm{sin}q_uu\mathrm{cos}q_vv,$$
(7)
where $`C`$ is taken to be a real number for the simplicity of the present discussion. In order that the inserted wave $`\psi _{sc}`$ is stable in the quadrangle, the wave function of the inserted wave must vanish on the boundaries of the quadrangle, that is, on the lines fulfilling $`|u|=u_B`$ or $`|v|=v_B`$. This constraint gives us the relations

$$q_u=2(L_u+1)k_u\mathrm{and}q_v=(2L_v+1)k_v/3,$$
(8)
where $`L_u`$ and $`L_v`$ are zero or positive integers. Let us investigate vortices in the quadrangle. From the fact that $`\psi _{cc}`$ has two nodal lines fulfilling $`v=\pm v_B/3`$ in the $`\zeta `$ plane and on the other hand $`\psi _{sc}`$ has nodal lines fulfilling $`u=n\pi /q_u`$ or $`v=(2m+1)\pi /2q_v`$ ($`n`$ and $`m`$ are integers), we see that vortices appear at the nodal points of the total wave function $`(u_V,v_V)`$, where

$$u_V=n\pi /q_u=\pm nu_B/(L_u+1),\mathrm{and}v_V=\pm v_B/3.$$
(9)
Considering the constraint $`|u|<u_B`$, we obtain the following results;

in the case of $`L_u=0`$ two vortices appear at $`(0,\pm v_B/3)`$ for $`n=0`$,

in the case of $`L_u=1`$ six vortices appear at $`(0,\pm v_B/3)`$ for $`n=0`$ and $`(\pm u_B/2,\pm v_B/3)`$ for $`n=1`$, and generally

$`2(L_u+2)`$ vortices appear on the lines of $`v=\pm v_B/3`$ for $`L_u>0`$. In general cases the following results are found out: (I)In this vortex-formation process we find out the selection rule for the inserted wave numbers $`\stackrel{}{q}`$ such that only the waves with $`q_u=(L_u+1)k_u/N_u`$ and $`q_v=(2L_v+1)k_v/2N_v`$ can format some stable vortex-patterns in the above example. In general the ratios of $`q_u/k_u`$ and $`q_v/k_v`$ must be rational numbers. (II) From the observation of the number of the vortices we can read off the number $`L_u`$ (or $`L_v`$). Furthermore we can add one more important fact: (III) In the processes putting standing waves into the quadrangle the number $`C`$ that represents the magnitude of the inserted wave play no role in the determination of the vortex patterns. 4. Mechanism of memories

Here let us consider the mechanism of memories in terms of the vortex-formation process. In the present framework we consider that every vortex pattern corresponds to a memory. In this viewpoint a quadrangle confining a standing wave is not a memory corresponding to a special item, but the system in which a vortex pattern is awoken by a flow induced by a stimulus. A quadrangle can awake different vortex patterns corresponding to different wave numbers expressed by the relations like $`q_u=(L_u+1)k_u/N_u`$ in the example of the section 3. We may visualize the following process; a stimulus corresponding to an item is taken by some organs and converted to flows. The flows reach at boundaries of many quadrangles with different sizes ($`4a`$-gons in the real space). Every quadrangle analyses the flow in terms of the standing waves given by (6). Notice that this analysis for the inserted flow is nothing but Fourier analysis in the $`\zeta _a`$ plane, of which boundaries are given by the lines $`u=\pm u_B`$ for $`u`$ and $`v=\pm v_B`$ for $`v`$. If the flow does not contain the waves with the wave numbers fitting to the selection rule (I), the flow cannot make any stable vortex pattern. Only the flows containing the wave numbers fitting to the rule (I) make stable vortex patterns in the quadrangles. Thus one stimulus sent to many quadrangles awakes various vortex patterns in those quadrangles. We may consider that the set of these vortex patterns made in the different quadrangles represents the memory corresponding to the stimulus. As noted in (III), the magnitude of the stimulus does not affect the vortex formation at all.

Let us see the characteristic features of this memory system. From the mechanical standpoint the fundamental properties of a quadrangle are determined by the parameters of the potential inside the quadrangle, that is, $`a`$ and $`g_a`$, where $`a`$ determines the form ($`4a`$-gon) of the $`xy`$ plane and $`g_a`$ does the fundamental wave number $`k_a`$. It must be noticed that further steps are needed for making a quadrangle confining a standing wave. One is the division of the fundamental wave number between two directions $`u_a`$ and $`v_a`$, following the relation $`k_a^2=k_u^2+k_v^2`$, that is, $`\stackrel{}{k}_a=(k_u,k_v)`$ is the wave number vector in the $`\zeta _a`$ plane. In order to determine the actual lengths of the sides of the quadrangle, we need one more step for the decision of two numbers $`N_u`$ and $`N_v`$ which determine the lengths of the sides of the quadrangle. Summarizing the above argument, the following three steps are needed to set up a quadrangle;

(1) the choice of $`a`$ and $`g_a`$ (the determination of the potential inside the quadrangle),

(2) the division of $`k_a`$ between $`k_u`$ and $`k_v`$ (the determination of the wave vector $`\stackrel{}{k}_a`$ confined

in the quadrangle), and

(3) the choice of the numbers of the wave lengths contained in the quadrangle, $`N_u`$ and $`N_v`$

(the determination of the lengths of the sides of the quadrangle). Through these processes a quadrangle confining a definite standing wave is determined. We see that in these processes we have five parameters, that is, $`a`$, $`g_a`$, $`N_u`$, $`N_v`$ and one more parameter $`\theta _a`$ to determine $`\stackrel{}{k}_a`$ such that $`k_u=k_a\mathrm{cos}\theta _a`$ and $`k_v=k_a\mathrm{sin}\theta _a`$. The angle parameter $`0<\theta _a2\pi `$ is a continuous free-parameter, and therefore the two wave numbers $`k_u`$ and $`k_v`$ can be also continuous numbers limited by $`k_a`$. This fact means that from one fundamental wave number $`k_a`$ we can make quadrangles having arbitrary values as for the ration of the lengths of two perpendicular sides, which is given by tan$`\theta _a`$. Taking also account of the free parameters $`N_u`$ and $`N_v`$, we have large freedom to the determination of the sizes and forms of quadrangles. We can, thus, provide quadrangles which can correspond to almost arbitrary inserted wave number vectors $`\stackrel{}{q}`$ even if $`\stackrel{}{q}`$ must be satisfied by the selection rule given in (I) to make stable vortex patterns in the quadrangles.

Let us continue to discuss about the mechanism for making memories. As already noted, the quadrangles are not the mechanism for preserving memories but that for awaking memories. One quadrangle characterized by a set of the four parameters $`k_a,\theta _a`$, $`N_u`$ and $`N_v`$ can fundamentally produce infinite numbers of different vortex patterns corresponding to the differences of the wave numbers of the inserted waves expressed by the relations $`q_u=L_uk_u/N_u`$ and $`q_v=L_vk_v/N_v`$, where $`L_u`$ and $`L_v`$ are positive integers or positive half-integers. It should once more be emphasized that one quadrangle awakes only one pattern for a stimulus but the pattern corresponding to the stimulus is generally different from others. Considering the limit of the discrimination of the vortex patterns, each quadrangle may have a certain maximum values $`(L_{\mathrm{max}})`$ for $`L_u`$ and $`L_v`$. In such a case, for the system where $`N_q`$ number of quadrangles cooperate we can roughly estimate the number ($`N_{VP}`$) of the different vortex patterns that can be distingushed by the system as $`N_{VP}=_{i=1}^{N_q}(L_{u,\mathrm{max}}\times L_{v,\mathrm{max}})_i.`$ Provided that the mean number of $`(L_{u,\mathrm{max}}\times L_{v,\mathrm{max}})_i`$ is 10 and $`N_q`$ is a large number, we obtain a huge number of the distinguishable vortex patterns counted to be $`N_{VP}10^{N_q}`$. Even if the mean number is 2 and $`N_q`$ is not a very large number, e.g. 100, we still have a huge number $`2^{100}10^{30}`$. From this fact we may consider that a total system consisting of a huge number of quadrangles is divided into many small blocks, and each block plays a fairly independent role in the total system.

We would like here briefly to note on the role of the higher order standing waves given in (5). Taking account of the freedom introduced in linear combinations of those standing waves, we will find out that the boundary conditions are expressed by some complicated relations for $`u_a`$ and $`v_a`$. On the other hand, from (5) we see that all the terms of the standing waves have one of the four factors $`\mathrm{cos}k_uu_a`$, $`\mathrm{sin}k_uu_a`$, $`\mathrm{cos}k_vv_a`$ and $`\mathrm{sin}k_vv_a`$. Therefore, as far as the discussions of vortex patterns are concerned, we can perform arguments similar to those done in the lowest cases. In actual situations, however, the fences of quadrangles ($`4a`$-gons in real spaces) will be made of some membranes in living beings. They are, of course, not perfectly rigid. In such cases the wave functions do not generally vanish at the boundaries of the quadrangles. Boundaries are, therefore, not given by strict quadrangles but by some deformed forms, and the values of wave functions on the boundaries are not strictly 0 but some values determined by the continuity conditions of wave functions. The infinite degeneracy of the standing waves can possibly play a role to solve the problem. That is to say, we will be able to provide various forms corresponding to various boundary conditions in terms of the infinite degeneracy. The vortex positions must also be changed. By gathering the terms with the same factor (cos or sin) discussions similar to those done in the lowest order cases will be also performed here. This infinite degeneracy and the potential parameters $`a`$ and $`g_a`$ will bring a large flexibility in the construction of the minimum units of the memory system. 5. Interesting properties of the system composed of the zero-energy flows

We have used only the zero-energy states for making the memory system. We would like to point out four interesting properties of the present system composed of the zero-energy flows ; (i) the absolute economy as for the energy consumption in the construction of the memory system and also the preservation of the system, where no flow in the $`z`$ direction should be taken, and therefore the total energy $`E_T=0`$. (ii) the almost infinite variety of the vortex patterns provided for a huge number of memories, and (iii) the perfect recovery of the memory system from shocks induced by stimuli, because the zero-energy ($`E_T=0`$) states have no time dependence at all, and therefore the standing waves in the minimum units of the system are perfectly recovered even if any disturbances induced by flows decaying or diffusing in some time scale are taken place . These properties are suitable for interpreting the mechanism of memories in brains. That is to say, stimuli taken in sense organs are analyzed in terms of the states including both of zero- and nonzero-energy states in the conjugate spaces of Gel’fand triplets, and after some time scale only the zero-energy flows selected by the rule (I) remain in the minimum units and make stable vortex patterns. The vortex patterns are sent to other places by the free motions of the inserted waves in the direction perpendicular to the $`xy`$ plane, and then every unit perfectly recover the initial state. We may point out one more important property; (iv) the large flexibility due to the infinite degeneracy and the parameters $`a`$ and $`g_a`$ in the potentials is provided in making the minimum units like the quadrangles. 6. A comment on thinking processes

Finally we would like to present a conjecture about thinking processes. At present we do not know how the mechanism of memories presented here does work in thinking processes (algorithm of thinking processes). We shall here point out a possibility that the memory system proposed here plays a role to distinguish accustomed stimuli, which have been repeatedly experienced in the past, from many stimuli included in an item caught by sense organs. We may consider that the same stimulus repeatedly given in our experience provokes organizations in brains to make a memory system consisting of many minimum units like quadrangles (possibly neurons in our brains) for selecting the very stimulus among many stimuli. Some systems like autonomic nervous system are made in the womb of mother. Such a memory system provided for the selection of some special stimuli will be made in a rather small block of the total system.

Now we can consider the following process for thinking: When an item is caught by some sense organs, the stimuli taken by the organs will be gathered in a special place of brains, where the stimuli are converted into signals that fit to the detection by the memory system. From the place the signals are delivered to many blocks, each of which consist of a certain number of minimum units. The signal sent to one of the blocks are analyzed by each unit in the block in terms of the standing waves, and then every unit selects the standing waves, following the selection rule (I). If there exist wave numbers fitting to the rule, the unit make a vortex pattern. Thus the signal sent to the block is converted to a set of vortex patterns. This set represents a recognition of an accustomed experience. Thus the stimuli caught by the organs are represented by a set of such recognitions . The recognitions characterized by various vortex patterns are once more delivered to a certain number of new blocks (call them higher blocks) and gathered together in each block. In this process the transportations will be carried out by using the free motions of the direction perpendicular to the standing waves. It should be noticed that the vortex patterns which are a kind of topological property can be very stable in such transportations. Then they proceed to the next step carried out in the higher blocks, where interactions between vortices take place. Some vortices will be created and/or annihilated , and a new vortex pattern appears in each block. Such new patterns produced in the higher blocks will be sent to more higher blocks, where the higher step is performed. The same processes will be repeated, until the brain gets a set of vortex patterns that is recognized to be the final answer. Such hierarchy may be the base of our thinking processes. As we know, flows of fluids are very much flexible and easily deformed, but not easy to control them perfectly. Experiences and trainings in daily life must play a very important role in the construction of this hierarchy. At present we have so many unkown problems such as the definite correspondence between the vortex patterns and the recognitons, the dynamics in the higher blocks, and so on. We are still at the entrance of the understanding. The view given here, however, is possibly an idea for the understanding of thinking in terms of real physical processes.