text:
# Horizontal Muons and a Search for AGN Neutrinos in Soudan 2

## 1 Introduction

The Earth provides an effective shield against the muon component of cosmic ray showers for overburdens exceeding 14,000 meters-water-equivalent (mwe). Consequently, near the horizontal direction ($`|\mathrm{cos}(\theta _z)|<0.14`$) at the site of the Soudan 2 detector (vertical depth 2090 mwe), the muon flux is comprised almost entirely of $`\mu ^\pm `$ initiated by charged current interactions of atmospheric muon neutrinos in the rock surrounding the detector. In addition, astrophysical neutrino sources may contribute to the (horizontal) muon flux. In particular, a high energy ($`>`$ 5 TeV) contribution may arise as the result of neutrinos which originate within astrophysical sources such as Active Galactic Nuclei (AGNs).

Active galaxies are those with extreme variations in brightness corresponding to abnormal emission of large amounts of energy at optical and/or radio wavelengths. Colliding and exploding galaxies also fall into this category. The nuclei of these galaxies may produce both neutral and charged particles. Most of the particles either decay or are absorbed before they can escape the galaxy, leaving only photons and neutrinos to be detected on Earth. Several models predict substantial high-energy neutrino production from AGNs, however there is a theoretical limit which conflicts with many of these models. AGNs may collectively give rise to a diffuse neutrino flux, or they may be detectable as individual point sources. Gamma Ray Burst sources have also been proposed as another possible diffuse source of high energy neutrinos.

Charged current interactions of multi-TeV muon neutrinos ($`\nu _\mu +\overline{\nu }_\mu `$) in the rock surrounding the Soudan 2 detector can be expected to give rise to TeV muons at the detector. Muons with these large energies will very frequently exhibit catastrophic energy loss by bremsstrahlung or pair production as they traverse the tracking calorimeter. This is highly improbable for muons ($`<E_\mu >20`$GeV at the detector) initiated by atmospheric neutrino reactions in the surrounding rock. Thus it is possible to measure the rate of high energy AGN neutrinos.

In this paper we report a measurement of the horizontal muon flux of $`4.01\pm 0.50\pm 0.30\times 10^{13}`$cm-2sr-1s-1 and a negative search for high energy muons from AGNs.

## 2 The Soudan 2 Detector

The detector is located approximately 700 m below the Earth’s surface in an historic iron mine in Soudan, Minnesota, USA. A photograph is shown is Figure 1.

The central detector is an iron sampling calorimeter using drift tubes filled with 85% Ar - 15% CO2 gas as the active medium. The detector is made of 224 modules with dimensions of $`2.7\times 1\times 1.1`$m3. Each module has a mass of $`4.3`$ t. The modules are stacked 2 high for a total height of 5.4 m, and 8 wide for a width of 8 m. The average density is $`1.6`$g/cm3.

The drift tubes are oriented along the North-South direction and are operated in proportional mode. Each end of the tubes is read out by an anode wire and a cathode pad which establish the x (East-West) and y coordinates of the hits. Off line, hits are matched in each half tube, using pulse height to match the cathode and anode pulses. The z (North-South) position of the hit is determined by correlating the drift times.

The ceiling, floor, and walls of the detector hall are covered with an active shield consisting of a two layer aluminum proportional tube array. The shield array resolves charged particle crossings (“hits”) to within four microseconds and within rectangular surfaces of 10 cm $`\times `$ 5 m. For most Soudan 2 analyses, the shield is used to veto events in the main detector which were initiated by particles interacting in the surrounding rock. Use of coincident in-time shield hits in the present analysis is described in Section 4.

The trigger for recording events requires a multiplicity of seven anode pulses or eight cathode pulses in blocks of 16 contiguous readout channels. For long muon tracks this trigger is fully efficient.

A near horizontal muon data event which is (presumably) neutrino-induced is shown in Figure 2. This example is the longest horizontal muon track that was found, with a length of 1464 cm and 149 tube crossings used in the fit. The zenith (azimuth) angle of this track is $`84^{}`$ ($`23^{}`$). Two shield hits can be seen in the top and front views.

The construction and performance of the tracking calorimeter are extensively described in Reference . Of particular interest to this analysis, the single track angular resolution has been measured using the moon shadow to be (0.3 $`\times `$ $`0.3^{}`$) and the electromagnetic energy resolution determined by counting tube crossings has been measured in the ISIS test beam to be

$$\frac{\mathrm{\Delta }E}{E}=\frac{7.0}{\sqrt{E}}+13.5\%$$
(1)
where E is given in GeV and the second term reflects the saturation of the number of hits.

## 3 Horizontal Muons at Soudan 2

Figure 3 shows the vertical muon intensity underground as a function of depth in standard rock, as compiled by Crouch and subsequently updated by the Particle Data Group. Two distinct components are apparent in the fitted curve, which consists of a double exponential plus a constant term. The atmospheric muon rate is observed to fall sharply with slant depth, so that the flatter spectrum of neutrino-induced muons only becomes visible for slant depths greater than 14 kmwe. For angles $`\theta _z`$ with respect to the vertical (zenith angle), the slant depth increases approximately as $`\mathrm{sec}(\theta _z)`$, but depends in detail on both the surface terrain and the rock density and composition.

Based on the vertical muon spectrum and the calculation of the average slant depth at Soudan 2 as shown in Figure 3, we can select a sample of neutrino induced muons by making a cut on the zenith angle near 82. Since the event timing resolution is about 100 ns, it is not possible to ascertain whether a muon is upward- or downward-going. There is thus a two-fold ambiguity in the direction for each measured track and $`\theta _z`$ is defined to be $`90^{}`$.

Calculation of the slant depth requires knowledge of both the rock density and the surface terrain. Differences between standard rock ($`\rho 2.62`$g/cm3) and Soudan rock (3.00 g/cm$`{}_{}{}^{3}>\rho >2.74`$g/cm3) must be taken into account. The topology above and surrounding the Soudan site is shown in Figure 4. The geology of the region is such that the iron deposits are almost vertical, so that the most vertical zenith angles correspond to the highest average density. The majority of the rock is greenstone composed of silicon dioxide (SiO2) and aluminum oxide (Al2O2). Based on Soudan rock effective parameters, a depth of 14 kmwe in standard rock corresponds to 13.2 kmwe in Soudan rock.

Due to the uneven surface terrain, the slant depth was calculated for each track individually. Observed fluxes of cosmic ray induced muons with $`\theta _z<60^{}`$ have been used to determine the zenith angle dependence of the rock density. By considering the azimuthal variation in the slant depth, we infer that systematic uncertainties in the atmospheric muon flux due to incomplete knowledge of the rock density are on the order of 10%.

Upon inclusion of the terrain information, a more refined estimation of the slant depth as a function of zenith angle is obtained, as indicated in Figure 5. The solid curve in that Figure was made assuming a spherical Earth. The dots which appear each one degree in zenith (two degrees in the inset) and are averaged over ten degrees in azimuth reflect the effect of topology variations in the overburden.

## 4 Reconstruction of Neutrino-Induced Muon Tracks

There were approximately $`10^8`$ triggered events in the Soudan 2 detector in the data set considered for this analysis. About half of these were due to cosmic ray muons, and half due to radioactivity in different parts of the detector which satisfied the trigger due to multiplexing. These events were processed using a pattern recognition filter code which identifies muon tracks and discriminates against unphysical background from noise. Event processing took place in three stages. Runs were processed by the mine computer soon after the runs were complete, and sorted into a number of output streams, one of which was a file consisting of events which included a track with a zenith angle greater than 60 degrees. Those files were then subjected to a filter program described below which identified horizontal muon candidates with a zenith angle greater than 78 degrees. There were 4458 events from the output of that filter which were scanned.

The filter program used five cuts which were applied to define a “horizontal muon” sample as listed in Table 1. Tracks were required to have end points consistent with a particle entering and leaving the detector. Acceptable track ends were required to 1) lie outside the fiducial volume which is 50 cm from the edge of the main detector, or 2) lie on a module boundary inside of the main detector and 3) point to a coincident shield hit or 4) point to a portion of the shield where a shield module did not exist. The azimuth cut rejected noisy events reconstructed parallel to the axes of the detector. The zenith cut defined the sample of horizontal muons. The track length cut minimized short tracks such as pions that originated from nuclear interactions within the rock nearby and also low energy cosmic ray muons which had undergone large multiple Coulomb scattering.

The anode-width cut rejected events that were primarily along one anode wire. These were actually along the north-south direction, but were often mis-reconstructed. The efficiency of these cuts for Monte Carlo muons generated isotropically with zenith angles exceeding 82 degrees is 0.56 $`\pm `$ 0.01 as given in Table 2. Note that both the trigger and reconstruction efficiency are lower near zenith angles of 90 degrees.

The pattern recognition software found single horizontal muon tracks, but it also kept a variety of background events. Most background events contained short horizontal tracks accompanied by other activity in the detector. Physicist scanning of all tracks greater than 78 was done to verify the track fits and to eliminate obvious backgrounds. Common backgrounds were due to certain patterns of electronic noise, vertical muons undergoing large radiative stochastic losses, multiple muons, and neutrino interactions within the detector. Tracks with visible multiple scattering (greater than 2 degrees) or with in-time shield hits that were not due to a horizontal muon were also rejected. All events were subjected to two independent scans and residual discrepancies were rescanned. The uncorrelated scan inefficiency as determined from the double scan was less than 1%. The Monte Carlo which was used for the efficiency calculation only generated high energy (100 GeV) muons. This fails to take into account neutrino induced muons with lower energies in the detector which might fail the multiple scattering cut applied during scanning. An independent calculation that took into account the neutrino induced muon energy distribution at the detector gave the result that 5% $`\pm `$ 5% of the events should undergo multiple scattering by more than 2 degrees in the detector, the criterion used during scanning. We therefore used a scan efficiency of 95% for the atmospheric neutrino flux calculation. Note that this inefficiency does not apply to the AGN $`\nu `$ search in the next section.

The track-length cut of 1.75 m (2.2 hadronic interaction lengths) minimized the number of short tracks which arose from processes which were not neutrino-related. The latter backgrounds included pion tracks originating from deep inelastic muon scattering within nearby rock, or else high energy cosmic ray muons from vertical directions which experienced large-angle multiple scattering. This length corresponds to a muon energy threshold of $`400`$MeV, but the restriction on multiple scattering raised this to an effective threshold of $`1.8`$GeV.

The number of accepted events after scanning, with $`\theta _z>78^{}`$, was 1237 events. The slant depth cut of $`14`$kmwe reduced this sample to 65 events. The track length, azimuth and zenith angle distributions of all 1237 events are shown in Figure 6. Figure 7 shows the ($`\theta _z`$, $`\varphi `$) distribution of these events. The depletions at regular intervals on azimuthal projections of these plots are‘ due to the azimuth cuts. The wavy-line contour on Figure 7 shows the 14 kmwe slant depth cut, which is designed to separate the sample into atmospheric muon and neutrino-induced muon candidates. The event with the largest zenith angle is shown in Figure 8.

To estimate the amount of background from atmospheric muons we integrate under the intensity-vs.-slant-depth curve of Figure 3, but without the constant (neutrino-induced muon) term, for slant depths greater than $`14`$kmwe. This procedure yields a background estimate of $`0.1_{0.1}^{+0.2}`$ events in our horizontal muon sample, where the error reflects the uncertainty in the rock density.

Muons which multiple scatter in the Earth above the Soudan 2 detector are a potential source of additional background. Since the overburden above the Soudan 2 detector is so great, it is not appropriate to apply the usual equation for small angle multiple scattering. Instead, a Monte Carlo calculation was performed which propagated muons through the Earth including multiple scattering and calculated the energy loss in steps. With a sample starting with 60 million muons that had $`\theta _z>60`$degrees at the surface of the Earth, the same number that were recorded during our exposure, no additional muons were generated with $`\theta _z>82`$degrees at the detector. The finite statistics of the Monte Carlo introduced an additional uncertainty which corresponds to $`<0.5`$ event background. We conclude that the expected background is less than one event.

The neutrino-induced muon flux, $`\mathrm{\Phi }_\mu `$, can be expressed as

$$\mathrm{\Phi }_\mu =\frac{N_\mu }{tA\mathrm{\Omega }\epsilon }$$
(2)
where $`N_\mu `$ is the number of observed neutrino induced muons, t is the time exposure of the detector in seconds, A is the effective area of the detector in cm2, $`\mathrm{\Omega }`$ is the solid angle subtended by the detector in steradians, and $`\epsilon `$ is the detection efficiency for the muons in the detector. Both A and $`\mathrm{\Omega }`$ are calculated assuming a uniform acceptance to the right of the contour in Figure 7.

The analyzed data spanned the period from April 14, 1992 to April 24, 2002. The exposure time was calculated from the start and end times of every processed data run. The calculated exposure for this analysis, including corrections for the detector duty cycle and electronics dead time during data taking, is $`2.00\times 10^8`$s. To estimate the effective area for this analysis, Monte Carlo muons were uniformly generated and the effective area for each track calculated. For $`\theta _z>82^{}`$ the average effective area was $`86.7\pm 0.3`$m2. The solid angle was calculated from the acceptance region in Figure 7 to be $`1.77`$sr. The trigger and reconstruction efficiency $`\epsilon `$ = 0.53 (0.556 from Table 2 $`\times `$ 0.95 from scanning.)

Systematic errors arise with the uncertainty in the background of 0.7 events (1%), uncertainties in the scan efficiency (5.6%) and uncertainty in the energy distribution as it affects the efficiency of the $`2^{}`$ multiple scattering cut (5%). Adding these errors in quadrature, we assign an overall systematic uncertainty to the neutrino flux calculation of 7.6%.

The resulting neutrino-induced muon flux for a muon energy threshold of 1.8 GeV is

$$\mathrm{\Phi }_{\nu _\mu }=4.01\pm 0.50(stat)\pm 0.30(sys)\times 10^{13}cm^2sr^1s^1$$
(3)
Upward-going atmospheric neutrino-induced muons have been previously measured by Baksan, IMB, MACRO, Kamiokande, and Super-Kamiokande. However, these experiments have not reported neutrino-induced muon fluxes above the horizon from those regions of solid angle where the overburden might be sufficient to separate atmospheric muons, perhaps due to uncertainties in the surface topology and rock density. The SNO experiment is deep enough that it should be able to see the neutrino-induced muon flux above the horizon. The flux near the horizon has also been measured by LVD and Frejus. The present measurement has aspects in common with that of Frejus, albeit with increased statistics. LVD, which has a much less uniform overburden, did not separate their neutrino-induced signal from background. The measured fluxes closest to $`\theta _Z=90^{}`$ are tabulated in Table 3 along with our estimate of the muon energy threshold used for each analysis.

Super-Kamiokande, and Kamiokande, MACRO, and Soudan 2 as well have analyzed atmospheric neutrino events and concluded that neutrino oscillations modify the zenith angle distribution of the $`\nu _\mu `$ flavor component of the atmospheric neutrino flux. For a neutrino energy of $`40`$ GeV and assuming $`\mathrm{\Delta }m^23.5\times 10^3`$eV2 and maximal ($`\theta _{23}`$) mixing, the probability of oscillation for a $`\nu _\mu `$ from $`82^{}`$ ($`<L>130`$km) is $`2.1\times 10^4`$ while for a $`\nu _\mu `$ from $`98^{}`$ ($`<L>`$ 2000 km) it is $`0.04`$. Our modest statistics and $`90^{}\theta _z`$ ambiguity make an oscillation analysis of these events impractical.

## 5 AGN Neutrino Search

Several models for neutrino production in Active Galactic Nuclei predict large and potentially measurable fluxes of very high energy neutrinos. These neutrinos would then produce high energy muons with energies from a few TeV up to 100s of TeV. In this energy range the dominant energy loss process for the muons in iron is electron pair production, followed by bremsstrahlung. Both processes produce large electromagnetic showers in the Soudan 2 detector which are easily detected and measured.

To understand the response of the tracking calorimeter to these TeV muons, a simulation was performed. High energy muons were propagated through the Soudan 2 detector and the amount of energy they deposited was calculated in a GEANT-based Monte Carlo. Three different muon energies were studied: 5, 20 and 100 TeV. Figure 9 shows the fraction of muons which lost at least a specific amount of energy versus that energy loss for three different muon energies. The results obtained from the simulation agree with an analytic calculation carried out independently. While 60% of the 5 TeV muons lose 5 GeV or more, 91% of the 20 TeV muons and 99% of the 100 TeV muons lose at least 5 GeV in the detector.

The muons identified for the horizontal muon flux measurement were subjected to a predetermined cut of 5 GeV on the amount of energy loss they experienced in the detector. None of the 65 events had visible radiated energy loss greater than this cut.

The largest zenith angle event with substantial visible energy loss is shown in Figure 10. This event has a horizontal muon with a zenith angle of $`80.6^{}`$ and an azimuthal angle of $`133^{}`$ which entered the top of the detector and exits the north wall. Parallel to the muon at the top is an entering electromagnetic shower, presumably from a bremsstrahlung or pair production process along the muon path in the rock above the detector. From Figure 7, the zenith angle ($`80.6^{}`$) tags the event as a high energy cosmic ray muon and not a horizontal muon candidate. Visible energy loss from a cosmic ray muon is not unexpected, because for atmospheric muons at the detector near zenith angle $`80^{}`$, $`<E_\mu >380`$GeV while for atmospheric neutrino induced muons, $`<E_\mu >20`$GeV. The energy of the shower is reconstructed to be 2.2 GeV. There are no such events among the 65 horizontal muon candidates.

Models for neutrinos from AGNs predict a $`\nu `$ flux as a function of energy. For each assumed functional form of the predicted neutrino energy spectrum, an observation or limit on high energy muons can be used to normalize or limit the neutrino flux. It is less meaningful to plot a limit on a neutrino energy distribution for an arbitrary functional form. Therefore we choose to show our limit on the neutrino-induced muon energy distribution. For each muon energy, our non-observation of events with energy loss corresponds to an upper limit on the muon flux.

To provide a comparison of our search with other searches reported in the literature, we show in Table 4 the exposures ($`A\times t\times \mathrm{\Omega }\times \epsilon `$). Our exposure is larger than that of the Frejus experiment, but lower than MACRO or AMANDA. Since the latter two experiments focus on upward-going muons, and the Frejus and Soudan 2 limits use horizontal muons, the limits are complementary for the exotic case involving the highest energy neutrinos for which the Earth’s neutrino opacity is large ($`E_\nu >50`$ TeV).

Based upon the observation of zero events, together with our efficiencies, we calculate 90% CL upper limits for the integral muon flux above three energies, which are shown in Table 5 and Figure 11. Our limit is close to the highest predicted flux; however that flux is not ruled out. The limits are lower for higher energy cut-offs because of the higher probability of detecting a large energy loss.

Other experiments have also used their horizontal muon events to set limits on the flux of the diffuse AGN neutrinos. The Frejus experiment used two methods to set a limit on the flux of high energy neutrinos. One method extrapolated a fit to the atmospheric neutrino energy spectrum obtained previously. Since that fit assumed that one spectral index fit all of their data, using it to deduce a limit for muons arising from a different energy distribution may bias that limit toward small values. The limit was quoted as $`d\varphi /dE_\nu (2.6`$TeV$`)<7.010^{13}`$GeV-1cm-2s-1sr-1. With their second method, Frejus used the absence of muons with large energy loss to constrain the normalization of a number of specific models of the high energy neutrino flux.

The MACRO experiment conducted an AGN search and found 2 candidate events with background of 1.1 expected. They set an upper limit on the muon flux from diffuse neutrinos of $`(1.7\pm 0.2)\times 10^{14}`$cm-2s-1sr-1. Their efficiency was highest for upward-going muons and lowest for horizontal muons. The AMANDA neutrino telescope, located in the ice near the South Pole, has set a limit based on an assumed diffuse $`E_\nu ^2`$ spectrum of $`dN/dE_\nu 10^6E_\nu ^2`$cm-2s-1sr-1GeV-1

As another test for an extraterrestrial component in the horizontal muon events, we have examined them for evidence of point sources. Figure 12 shows the directions of the 65 horizontal muons in galactic coordinates using an Aitoff projection of the galaxy. For each muon, the two possible directions are plotted. We observe that no clusters involving two or more muons within the detector’s angular resolution $`(0.3^{}\times 0.3^{})`$ appear within the plot.

## 6 Summary

In this study we have isolated a sample of 65 neutrino-induced muon candidates, with an estimated background of less than one event. For a slant depth cut 14 kmwe and a muon energy threshold of 1.8 GeV, this corresponds to a muon flux of $`4.01\pm 0.50^{stat}\pm 0.30^{sys}\times 10^{13}`$cm-2sr-1s-1 in the horizontal direction $`0.14<\mathrm{cos}\theta _z<0.14`$. None of the 65 events have a large energy loss in the detector. We set an integral limit on neutrino-induced muons from AGNs and/or other sources from between 1.4 and 2.2 $`\times 10^{14}`$cm-2sr-1s-1 depending on muon energy.

## 7 Acknowledgments

We thank Rick Egeland, Leah M. Goodman and Adrian Ryerson for their help in this analysis. This work was supported by the U.S. Department of Energy, the U.K. Particle Physics and Astronomy Research Council, and the State and University of Minnesota. We also wish to thank the Minnesota Department of Natural Resources for allowing us to use the facilities of the Soudan Underground Mine State Park.


format: text

text:
# Classes of exact wavefunctions for general time-dependent Dirac Hamiltonians in 1+1 dimensions

## Abstract

In this work we construct two classes of exact solutions for the most general time-dependent Dirac Hamiltonian in 1+1 dimensions. Some problems regarding to some formal solutions in the literature are discussed. Finally the existence of a generalized Lewis-Riesenfeld invariant connected with such solutions is discussed.

PACS numbers: 03.65.Ge, 03.65.Pm

The interest in solving problems involving time-dependent systems has attracted the attention of physicists since a long time. This happens due to its applicability for the understanding of many problems in quantum optics, quantum chemistry and others areas of physics -. In particular we can cite the case of the electromagnetic field intensities in a Fabry-Pérot cavity . In fact this kind of problem still represents a line of investigation which attracts the interest of physicists . However, as has been observed recently by Landim and Guedes , the most part of these works deal with nonrelativistic systems. So they tried to bridge this gap, by studying the problem of a fermion under the presence of a time-dependent Lorentz vector linear potential in two-dimensional space-time. For this they used the so called Lewis-Riesenfeld invariant operator, in order to guess the wavefunction which solves the problem. Unfortunately, they just presented a formal solution for the problem under question. In fact, if one consider that in the most part of the works treating the non-relativistic time-dependent systems, the solution is not formal only for a limited number of time-dependence for the potential parameters, it is a little bit strange that in  there is no limitation over those parameters. In view of these comments, as a first example, we are going to complete the program initiated by Landim and Guedes, by solving the problem up to the end, showing that there exists really some restrictions for a complete solution of this problem. After that, we will show that there exists at least two classes of solutions for the time-dependent Dirac equation in 1+1 dimensions. This is done considering the most general combination of Lorentz structures for the potential matrix. Then we show that the problem studied in  is only a particular case of one of those classes. Finally we construct the Lewis-Riesenfeld invariant operators which have our solutions as eigenfunctions.

Let us begin by presenting the Dirac equation in 1+1 dimensions, in the presence of a time-dependent potential for a fermion of rest mass $`m`$,

$$H\psi (q,t)=i\frac{\psi (q,t)}{t},H=\alpha p+\beta m+V_V(q,t)+\beta V_S(q,t)+\alpha \beta V_P(q,t),$$
(1)
where we used $`c=\mathrm{}=1`$ and $`p`$ is the momentum operator. $`\alpha `$ and $`\beta `$ are Hermitian square matrices satisfying the relations $`\alpha ^2=\beta ^2=1,\{\alpha ,\beta \}=0`$. Here we choose $`\alpha =\sigma _1`$and $`\beta =\sigma _3`$, where $`\sigma _1`$ and $`\sigma _3`$ are the $`2\times 2`$ standard Pauli matrices. In this way $`\psi `$ is a two-component spinor. Furthermore the subscripts for the terms of potential denote their properties under a Lorentz transformation: $`V`$ for the time component of the two-vector potential, $`S`$ and $`P`$ for the scalar and pseudoscalar terms, respectively. The absence of the space component of the two-vector potential is due to the possibility of its elimination through a gauge-like transformation.

In the particular case considered by Landim and Guedes , one has

$$H=\alpha p+\beta m+f\left(t\right)q.$$
(2)
As suggested by them, we will use as an Ansatz for the solution the following spinor

$$\psi =\left(\begin{array}{c}M_1\left(t\right)\hfill \\ M_2\left(t\right)\hfill \end{array}\right)e^{i\eta \left(t\right)q},$$
(3)
so that the corresponding coupled equations are

$`i\dot{M}_1`$ $`=`$ $`\left[\left(\dot{\eta }+f\right)q+m\right]M_1+\eta M_2,`$
$`i\dot{M}_2`$ $`=`$ $`\left[\left(\dot{\eta }+f\right)qm\right]M_2+\eta M_1,`$
where the dot denotes differentiation with respect to $`t`$. Imposing that $`\eta \left(t\right)=^tf\left(\lambda \right)𝑑\lambda `$ we eliminate the spatial dependence of the above equations so obtaining

$`i\dot{M}_1`$ $`=`$ $`mM_1+\eta M_2,`$
$`i\dot{M}_2`$ $`=`$ $`mM_2+\eta M_1.`$
Now we simplify even more the above equations if we perform the following identifications

$$M_1G_1e^{imt},M_2G_2e^{imt},$$
(6)
getting

$$\dot{G}_1=\eta _1G_2,\dot{G}_2=\eta _2G_1,$$
(7)
where $`\eta _1\eta e^{2imt}`$ and $`\eta _2\eta e^{\mathrm{\hspace{0.17em}2}imt}`$. By deriving the each one of the above equations it is easy to conclude that one gets the following second-order equations

$$\ddot{G}_i\left(\frac{\dot{\eta }_i}{\eta _i}\right)\dot{G}_i+\eta ^2G_i=0,i=1,\mathrm{\hspace{0.17em}2}.$$
(8)
At this point it is convenient to rescale the $`G_i`$ functions as $`G_i=\sqrt{\eta _i}g_i`$, one gets finally that

$$\ddot{g}_i+\left[\frac{1}{2}\left(\frac{\ddot{\eta }_i}{\eta _i}\right)\left(\frac{\dot{\eta }_i}{\eta _i}\right)^2+\eta ^2\right]g_i=0.$$
(9)
Note that in general these last equations are not somewhat straightforward to solve as asserted in . However we note that for a time exponentially decaying force, these equations takes the form a Schroedinger equation for the Morse potential in the time variable, without the usual boundary conditions of such kind of equation.

¿From now on we will present an extension of the class of time-dependent relativistic systems with exact solutions. For this purpose we begin with the complete Hamiltonian (1), and make the more general Ansatz

$$\psi =\left(\begin{array}{c}M_1\left(t\right)e^{iF_1(q,t)}\hfill \\ M_2\left(t\right)e^{iF_2(q,t)}\hfill \end{array}\right).$$
(10)
Thus we obtain

$`i\dot{M}_1`$ $`=`$ $`\left[\dot{F}_1+V_V+V_S+m\right]M_1+\left(F_2^{}V_P\right)e^{i\left(F_2F_1\right)}M_2,`$
$`i\dot{M}_2`$ $`=`$ $`\left[\dot{F}_2+V_VV_Sm\right]M_2+\left(F_1^{}+V_P\right)e^{i\left(F_2F_1\right)}M_1,`$
where the prime denotes differentiation with respect to $`q`$. As by construction $`M_i`$ do not depend on $`q`$, one sees that it is mandatory to get rid of such a dependence in the above equations. As a consequence two classes of solution emerge, corresponding to $`F_1=F_2`$ or $`F_1=F_2`$.

Let us analyze the class where $`F_1=F_2=F`$. In this case the exponential factor appearing in the equations disappears and, consequently the only way of getting the off-diagonal terms independent of $`q`$, is by imposing that

$$V_P=V_P\left(t\right),F=\theta _1\left(t\right)q.$$
(12)
On the other hand, the effect of this condition over the diagonal terms results that

$$V_S=V_S\left(t\right),V_V=\dot{\theta }_1\left(t\right)q+\theta _2\left(t\right).$$
(13)
Note that this class allows the treatment of systems that are at most linear in the spatial coordinate, and includes the system proposed by Landim and Guedes as a particular case. The general equations to be solved in this class is now given by

$`i\dot{M}_1`$ $`=`$ $`\left(\chi _++m\right)M_1+\eta _{}M_2,`$
$`i\dot{M}_2`$ $`=`$ $`\left(\chi _{}m\right)M_2+\eta _+M_1,`$
with $`\chi _\pm \theta _2\pm V_S`$ and $`\eta _\pm \theta _1\pm V_P`$. The above equations can now be decoupled giving

$`i\ddot{G}_1\left[\chi _++\chi _{}+i\left({\displaystyle \frac{\dot{\eta }_1}{\eta _1}}\right)\right]\dot{G}_1\left[\dot{\chi }_+\left({\displaystyle \frac{\dot{\eta }_1}{\eta _1}}\right)\chi _+i\eta _+\eta _{}+i\chi _+\chi _{}\right]G_1`$ $`=`$ $`0,`$
$`i\ddot{G}_2\left[\chi _++\chi _{}+i\left({\displaystyle \frac{\dot{\eta }_2}{\eta _2}}\right)\right]\dot{G}_2\left[\dot{\chi }_{}\left({\displaystyle \frac{\dot{\eta }_2}{\eta _2}}\right)\chi _{}i\eta _+\eta _{}+i\chi _+\chi _{}\right]G_2`$ $`=`$ $`0,`$
where we made the same definition as above for the $`G_i`$ functions and now $`\eta _1\eta _{}e^{2imt}`$ and $`\eta _2\eta _+e^{\mathrm{\hspace{0.17em}2}imt}`$. It is easy to verify that the above equations recall that appearing in (8) in their particular case. Once again it is in general not solvable, showing that we have got a formal solution, so that if one want to solve the problem until the end, one must to look for particular cases where these equations have explicit solutions.

¿From now on we analyze the class where $`F_1=F_2=F`$. In this new class, the exponential factor appearing in the equations holds and, consequently the only way of getting the off-diagonal terms independent of $`q`$ in the equations (LABEL:3), is by imposing that

$$F^{}=V_P(q,t).$$
(16)
Now, the condition of spatial independency of the diagonal terms implies that

$$V_V=V_V\left(t\right),V_S=V_S(q,t)=\gamma \left(t\right)\dot{F}(q,t),$$
(17)
where $`\gamma \left(t\right)`$ is an arbitrary function of the time. In this case the equation decouple giving simply

$$i\dot{M}_1=\left(V_V+\gamma +m\right)M_1,i\dot{M}_2=\left(V_V\gamma m\right)M_{2.}$$
(18)
Now it is easy to obtain the solution of these equations. These are given by

$`M_1\left(t\right)`$ $`=`$ $`M_1\left(0\right)e^{i\left\{mt+^t\left[V_V\left(\lambda \right)+\gamma \left(\lambda \right)\right]𝑑\lambda \right\}},`$
$`M_2\left(t\right)`$ $`=`$ $`M_2\left(0\right)e^{i\left\{mt^t\left[V_V\left(\lambda \right)\gamma \left(\lambda \right)\right]𝑑\lambda \right\}}.`$
At this point some comments are in order. From above we easily conclude that for this second class of solutions, one really obtains arbitrary non-formal solutions, provided that the remaining integrals can be done. Here it is important to remark that this new class of systems allows arbitrary dependence on the spatial and time variables, in contrast to the previous class of systems.

In what follows, we try to construct an invariant operator of the Lewis and Riesenfeld type , which have the above classes of solutions as eigenfunctions. It could be used, as done in , in order to suggest the form of the spinor presented as an Ansatz in this work. The invariant obviously must satisfy the equation

$$\frac{dI}{dt}=\frac{I}{t}i[I,H]=0,$$
(20)
and can be written as

$$I=𝒜\left(t\right)p+(q,t),$$
(21)
where $`𝒜`$ and $``$ are $`2\times 2`$ matrices, whose elements, in order to guarantee the validity of the above equation (20), must obey the following set of coupled differential equations

$`𝒜_{11}`$ $`=`$ $`𝒜_{22}=const.,𝒜_{12}=𝒜_{21}=const.,`$
$`\mathrm{\Lambda }_{}^{}`$ $`=`$ $`𝒜_{11}V_P^{},\mathrm{\Sigma }_{}=𝒜_{12}\left(m+V_S\right),`$
$`\dot{\mathrm{\Sigma }}_+`$ $`=`$ $`\mathrm{\Lambda }_++\mathrm{\Lambda }_{}𝒜_{11}V_V^{},`$
$`i\dot{\mathrm{\Sigma }}_{}`$ $`=`$ $`2\mathrm{\Lambda }_+i𝒜_{12}V_P^{}i𝒜_{11}V_S^{},`$
$`i\dot{\mathrm{\Lambda }}_+`$ $`=`$ $`2\mathrm{\Sigma }_{}V_P+i\mathrm{\Sigma }_{}^{}+i𝒜_{12}V_V^{}+2\mathrm{\Lambda }_{}\left(m+V_S\right),`$
$`i\dot{\mathrm{\Lambda }}_{}`$ $`=`$ $`i\mathrm{\Sigma }_+^{}+i𝒜_{11}V_P^{}+i𝒜_{12}V_S^{}2\mathrm{\Lambda }_+\left(m+V_S\right),`$
where we defined that $`\mathrm{\Lambda }_\pm \frac{_{12}\pm _{21}}{2}`$ and $`\mathrm{\Sigma }_\pm \frac{_{11}\pm _{22}}{2}`$. On the other hand, we are looking for invariant operators which have the Dirac spinors as eigenfunctions, which implies into further conditions over the elements of the matrices $`𝒜`$ and $``$. This lead us to the following conditions if we consider the eigenvalue equation for the invariant operator acting over the spinor (10) introduced in this work

$`𝒜_{11}(F_1^{}+F_2^{})+_{11}+_{22}`$ $`=`$ $`\chi _1\left(t\right),`$
$`\left(𝒜_{12}F_1^{}+_{21}\right)\left(𝒜_{12}F_2^{}+_{12}\right)`$ $`=`$ $`\chi _2\left(t\right),`$
with the functions $`\chi _{1,2}`$ being arbitrary functions of time. It is not difficult to show that for the first class one can obtain a solution, which includes that proposed in  as a particular case, given by

$`𝒜_{12}`$ $`=`$ $`𝒜_{21}=_{12}=_{21}=0,𝒜_{11}=𝒜_{22}𝒜,_{11}=_{22},`$
$`V_S`$ $`=`$ $`V_S\left(t\right),V_P=V_P\left(t\right),V_V={\displaystyle \frac{\dot{}}{𝒜}}q+\chi \left(t\right).`$
For the second class of solutions however, as can be easily realised from the above set of coupled nonlinear equations (LABEL:nonlinear) and the constraints coming from (LABEL:constraint), it is quite difficult to extract simple solutions. For this reason we were not able to verify if there exist a relativistic invariant which is linear in the momentum, responsible for the generalized solutions here presented. In fact, as it is typical of nonlinear equations, different solutions of it could lead to independent sets of wavefunctions. Finally it is interesting to make some comments, the first one being the observation that the extension of the solutions discussed in this work, by introducing a time-dependent mass is quite simple, as can be seen due to its appearance combined with time-dependent potentials $`V_S`$ and $`V_V`$. On the other hand, the case of the so called Dirac oscillator is not directly solvable from the second class of solutions, this happens because when $`V_P`$ is linear in the coordinate, it implies into the need of a $`V_S`$ quadratic in $`q`$, in order to guarantee the exact solution of the form guessed in this work.

Acknowledgments: The authors are grateful to FAPESP and CNPq for partial financial support.


text:
# 1 Introduction

## 1 Introduction

Characteristics of nonperturbative QCD make intractable many calculations at low and intermediate energies. Effective quark models therefore retain their usefulness in numerous applications. For example, ref.  used the instanton-extended version 111In ref. , incorporating instanton-induced interactions into the original MIT bag model , was inspired by an analysis  made in the constituent quark model, which found that an effective instanton interaction led to as satisfactory a description of mass splittings of baryons as the conventional approach using one-gluon exchange. of the MIT bag model  in one of many studies of strangeness in nucleons . However, the approach of ref.  contained the so-called linearized approximation, amounting to freezing the baryon radii in their original MIT values. In the present paper we remove this approximation and calculate the effects thereof on the baryon mass splittings, and also on the nucleon strangeness results of ref. . We also explore whether this enables the resolution or alleviation of a long-standing inconsistency between the MIT bag model and nuclear physics: the standard nuclear physics descriptions employ independent nucleons, while the nucleon MIT bag radius is too large for that .

The first approach to consider the instanton-induced interaction within a bag model was due to Kochelev . It is nevertheless important to note that he considered his own bag model , which is somewhat different from the MIT one. As explained in detail in ref. , that line of research  is therefore rather different from the one in ref.  and in the present paper, where we stay as close as possible to the original MIT bag model . The only modification with respect to the MIT model is the inclusion of the instanton-induced interaction . This inclusion is necessary also inside the MIT bag if one allows for the nonvanishing (even if small) probability of penetration of the instanton liquid from the surrounding nonperturbative QCD vacuum, into the “perturbative” MIT bag interior. The instanton density $`n`$ used in this effective instanton-induced interaction inside the bag, is of course reduced with respect to the density in the nonperturbative QCD vacuum: the smaller the probability of this penetration, the larger the reduction. The reduced value of $`n`$ appropriate for the MIT bag interior comes out as a result of our model calculation and fitting.

Besides defining how to incorporate instantons in the MIT bag model and finding the baryon mass shifts caused by the effective instanton-induced quark-quark interaction, ref.  explored the modification of the too large value (required by the “instantonless” model fits, e.g., ) of the strong coupling constant $`\alpha _c`$ used in the supposedly perturbative MIT bag interior. The conclusion of ref.  was that the change in $`\alpha _c`$ was in the desired direction, i.e., it was reduced, but only by about $`6\%`$, which was insufficient to achieve improvement in the consistency of the perturbative description of the bag interior. Actually, it turned out that instanton effects should generally be small in the MIT bag model, since the instanton (plus anti-instanton) density $`n`$ appropriate for the MIT bag interior, was found  to be very much depleted with respect to the instanton density estimated (e.g., by refs. ) for the true, nonperturbative QCD vacuum. The corresponding instanton-induced mass shifts were of the order of only a few MeV. However, the analysis of ref.  (and, consequently, of ref. ) contained an important simplifying assumption: it kept the baryon radii “frozen” in the values obtained by DeGrand et al. . In this way, a full refitting of the bag model parameters (now also including the instanton density $`n`$ inside the bag) was avoided. In fact, this highly nonlinear problem was thereby reduced to solving the four linear equations for the four adjustable model parameters that enter the energy functional linearly: $`\alpha _c`$, the volume-energy density $`B`$, the zero-point energy parameter $`Z_0`$, and a parameter new to the MIT bag model, namely, the instanton density $`n`$. (The quark masses were also not allowed to vary. Ref.  adopted the quark masses of DeGrand et al.  in order to be able to use the results of ref.  and to make a comparison with their results.) The linear equations determining the appropriate value of $`n`$ and the new values of $`\alpha _c`$, $`B`$, and $`Z_0`$ were specified by demanding that the model masses of the proton, neutron, $`\mathrm{\Delta }`$, and $`\mathrm{\Omega }^{}`$, be equal to the empirical masses after the inclusion of the effective instanton-induced interaction. We will call the approach of ref.  the linearized approximation.

In this paper we go beyond this approximation, performing a refitting of baryon masses which allows their radii to vary. It turns out that it leads to larger instanton densities allowed inside the MIT bags, and correspondingly to a stronger share of instantons in the energy balance of the baryon bags, accompanied by decreased $`\alpha _c`$, as well as by acceptable, and for nucleons even highly favorable , changes in the baryon radii. The most important for the present paper is amending the results on the nucleonic scalar strangeness obtained in ref. . On the one hand, larger instanton densities now lead to increased contributions of the instanton–induced interaction to the total scalar strangeness of the nucleon. On the other hand, the basic MIT bag strangeness (8) found by Donoghue and Nappi , if not far from its naive limit, may well still represent the main contribution. If it does, the total nucleon strangeness decreases with diminishing bag radii, which are in turn associated with growing instanton densities.

## 2 Refitting of the baryon bag parameters

Except for removing the linearized approximation, i.e., replacing it with the refitting where the bag radii are not frozen any longer, the incorporation of the instanton effects in the MIT bag follows ref.  closely. The same holds also for other model details, such as the fixed model inputs, the nonstrange and strange quark mass parameters ($`m_u=m_d=0`$ and $`m_s=279`$ MeV, respectively) and quark-antiquark ($`q\overline{q}`$) condensate $`0|\overline{q}q|0=(240\mathrm{MeV})^3`$. Thus, to keep the present paper as concise as possible, we refer to ref.  for all model details and parameters, and to ref.  for the corresponding strangeness calculation. (For detailed technicalities of the latter, ref.  may also be found helpful.)

Here we just recall that the effective instanton-induced interaction $`_I`$, causing the instanton-induced mass shift $`E_I^{}`$ of the baryon $``$, is the sum of the one-, two-, and three-body terms, denoted by $`_1^I`$, $`_2^I`$, and $`_3^I`$, respectively:

$$E_I^{}=|:_I:|=|:_1^I_2^I_3^I:|,$$
(1)
and is defined in detail in ref. . The explicit expressions for the one- and two-body contributions ($`\mathrm{\Delta }M_{}^{(1)}`$ and $`\mathrm{\Delta }M_{}^{(2)}`$, respectively) are also given in ref. .

Before proceeding, let us make two comments regarding our choice of the instanton-induced interaction $`_I`$. It was derived by Nowak et al.  in the framework of random instanton liquid model (RILM). They arrived at the interaction corresponding to the well-known one of Shifman, Vainshtein, and Zakharov (SVZ) , apart from the effects of smearing over the size of an instanton. In the limit of no smearing, it reproduces our chosen  local $`_I`$, which is essentially the same as the SVZ interaction . Since the SVZ interaction is induced by a single (anti-)instanton, our modeling misses multi-instanton effects. Their importance, however, was stressed in, e.g., refs. , putting in doubt the validity of the single-instanton approximation. The caveat is that these effects can be important when baryon bags have diameters larger than average separation of (anti-)instantons, and this will turn out to be the presently relevant situation (since we will find instanton densities inside bags up to one third of the QCD vacuum value of 1 fm-4). Nevertheless, as discussed especially in ref. , we should recall that this interaction was introduced and used  with the aim of capturing the intermediate-range ($`\frac{1}{3}`$ fm) QCD effects, and the interaction we adopted is suitable for that, since the average instanton size is $`\rho \frac{1}{3}`$ fm . Hopefully, it may capture the effects at ranges even a little beyond $`\frac{1}{3}`$ fm, since Nowak et al.  took into account the delocalization of zero modes222Thus, Nowak et al.  took into account the insights of, e.g., refs. , concerning the importance of summing up a large number of interactions with different instantons. The review by Schafer and Shuryak still points out as useful their results and RILM approach in general, observing that interactions among instantons (and hence their correlations) are important but not dominant .. In keeping with the basic idea of the MIT model, one assumes that really long range (i.e., confinement) effects are modeled well by the confining bag boundary.

The second comment is devoted to clarifying our inclusion of the one-body term $`_1^I`$ into the bag model calculations of the instanton-induced contribution (1) to baryon masses. The term $`_1^I`$ has in fact the form of a mass term, and can be thought of as the self-energy, or the effective mass that a quark acquires from the effective interaction caused by the instanton liquid through which quarks move. Now imagine that we are working in some kind of constituent quark model where one from the start uses effective constituent quark masses to parameterize “dressing” by nonperturbative QCD. The self-mass part of the instanton effects would in that case be already included in the constituent mass parameters. Using $`_1^I`$ in the baryon mass calculation would therefore be double-counting, so in that case it must be dropped from Eq. (1). On the other hand, if we employ some approach where one uses the current, Lagrangian quark masses, like in the MIT bag model used presently and in refs. , $`_1^I`$ should be included in the calculation on equal footing with $`_2^I`$ and $`_3^I`$. This procedure was criticized by Dorokhov  on the grounds that in the bag model, the role of the quark constituent mass is played by the single-quark kinetic energy eigenvalue resulting from the boundary condition confining the quarks inside the bag. According to this view, the quark nonperturbative dressing due to the $`_1^I`$ part of the instanton-induced interaction would already be taken into account by the linear bag boundary condition. However, we do not accept this view because this boundary condition serves to incorporate confinement, prohibiting quark separations larger than the bag diameter scale of the order of some 2 fm, while instantons are not responsible for confinement  (contrary to what was thought in early days of instanton physics). Admittedly, this argument is so far only qualitative in the sense that in the model context it is not possible to delineate precisely beyond which scale confinement effects overwhelm instanton effects. Nevertheless, the argument becomes stronger and more precise if one remembers the discussion in the previous passage: there, it was noted that the adopted instanton-induced interaction approximates well the nonperturbative QCD effects at intermediate ranges around $`\frac{1}{3}`$ fm, but not much further than that, and certainly not up to confinement scales of the order of the bag diameter.

As remarked above, the interaction $`_I`$ is actually the same as the well-known SVZ interaction , including the (only seemingly different ) three-body term $`_3^I`$. This term was in fact discussed in ref.  because, at that point, it was not clear whether the contribution of $`_3^I`$ vanished for $`\mathrm{\Lambda }`$, as it did for other baryons. Therefore, ref.  avoided the need to compute the contribution from the complicated-looking $`_3^I`$ by showing that it could contribute only to the mass of the $`\mathrm{\Lambda }`$ and by omitting the $`\mathrm{\Lambda }`$ from the analysis. However, it turns out that the mass shift due to the three-body interaction, if nonzero, must be small for the $`\mathrm{\Lambda }`$ . (In explicit evaluation one can see that all terms in the $`_3^I`$-contribution would cancel in the SU(3)-symmetric limit. This contribution slightly differs from zero only because the strange-quark wave functions differ somewhat from the nonstrange ones.) Neglecting therefore this contribution to $`M_\mathrm{\Lambda }`$, the total instanton-induced mass shift (1) consists of one- and two-body contributions only :

$$E_I^{}=|:_1^I_2^I:|\mathrm{\Delta }M_{}^{(1)}+\mathrm{\Delta }M_{}^{(2)},$$
(2)
for all baryons $``$, including the $`\mathrm{\Lambda }`$. There is hence no need to drop the $`\mathrm{\Lambda }`$ from the analysis, so in this respect, this calculation is slightly more complete than in ref. . However, when we did drop the $`\mathrm{\Lambda }`$ from the present fit in order to check the effects thereof, the results were affected very little.

Therefore, the only significant difference in modeling with respect to ref.  is that in the present paper we want to perform a full refitting of the model parameters, including the variation of the bag radii. Maybe some reader might then object that for each baryon $``$, its bag radius would become a new free parameter, and the number of fitting parameters would become larger than the number of experimental baryon masses $`M_{exp}^{}`$ to be be fitted. Fortunately, this is not so, because each radius $`R_{}`$ of a bag in equilibrium must satisfy the pressure-balance condition. That is, the equilibrium bag radius $`R_{}`$ of the baryon $``$ is fixed by minimizing the bag model mass $`M_{bag}^{}`$,

$$\frac{dM_{bag}^{}}{dR_{}}=0(=N,\mathrm{\Lambda },\mathrm{\Sigma },\mathrm{\Xi },\mathrm{\Delta },\mathrm{\Sigma }^{},\mathrm{\Xi }^{},\mathrm{\Omega }),$$
(3)
and is not a free, adjustable parameter like $`Z_0,\alpha _c,B`$ and $`n`$.

The MIT bag energy functional $`M_{bag}^{}`$ of the baryon $``$ now depends also on the instanton density $`n`$, because

$$M_{bag}^{}[R_{},Z_0,\alpha _c,B,n]=E_Q^{}+E_0^{}+E_M^{}+E_E^{}+E_V^{}+E_I^{}$$
(4)
now contains the instanton contribution $`E_I^{}`$ (1), in addition to the kinetic energy of the confined quarks $`E_Q^{}`$, the zero-point energy $`E_0^{}`$, the color magnetic energy $`E_M^{}`$, the color electric energy $`E_E^{}`$, and the volume energy $`E_V^{}`$. The expressions for these five latter contributions are given in ref. . (In Eq. (4), the dependence of $`M_{bag}^{}`$ on the quark mass parameters $`m_u,m_d,m_s`$ and the condensate $`0|\overline{q}q|0`$ is not indicated, as they are not adjustable parameters but fixed model inputs.)

In the circumstances explained above, the most practical and numerically tractable way to perform the model fit to the empirical baryon masses, is to pose it as the problem of minimization of the positive definite functional $`F`$

$$F[\{R_{}\},Z_0,\alpha _c,B,n]F_M+F_R,$$
(5)
$$F_M\underset{}{}(M_{exp}^{}M_{bag}^{})^2,$$
(6)
$$F_R\underset{}{}\frac{1}{^2}\left(\frac{dM_{bag}^{}}{dR_{}}\right)^2,$$
(7)
where the both sums run over the baryons $``$ in the ground-state octet ($`N,\mathrm{\Lambda },\mathrm{\Sigma },\mathrm{\Xi }`$) and decuplet ($`\mathrm{\Delta },\mathrm{\Sigma }^{},\mathrm{\Xi }^{},\mathrm{\Omega }`$). Thus, note that in the present fitting procedure all baryon masses enter on an equal footing, whereas ref. , similarly to ref. , chooses some masses somewhat arbitrarily to fix the parameters and predict the other masses.

In the functional $`F`$, the first sum, $`F_M`$, represents the deviation of the bag model masses $`M_{bag}^{}`$ from the experimental baryon masses $`M_{exp}^{}`$. The second sum, $`F_R`$, is a measure of the deviation from the situation of the perfectly satisfied pressure-balance condition. The role of the constant $``$ is just to ensure that the both terms have the same dimension, and we choose the typical baryonic mass scale of 1 GeV to fix its value: $`=1`$ GeV. (Of course, there is some arbitrariness in the choice of the functional $`F`$; for example, we could replace $``$ by $`M_{exp}^{}`$ in each term of the sum $`F_R`$. However, we have checked that varying the scale $``$ does not influence our results significantly, giving us confidence that this arbitrariness is not a problem in practice.)

The functionals defined by Eqs. (5)-(7), namely $`F`$, $`F_M`$ and $`F_R`$, all depend on the model parameters $`Z_0,\alpha _c,B,n`$ and on the set of the bag radii $`\{R_{}\}`$ of the octet and decuplet baryons. The strict approach to the model fitting through the functional minimization would be to pick the initial values of the free parameters, $`Z_0^{(0)},\alpha _c^{(0)},B^{(0)},n^{(0)}`$, and find the equilibrium radii $`R_{}`$ by minimizing the functional $`F_R`$ to, ideally, $`F_R=0`$, where the conditions (3) would be strictly satisfied. Then, the functional $`F_M`$ should be calculated. This two-step process should be repeated with varied values of the free parameters over and over again by some minimization routine (for example based on simplex minimization) till $`F_M`$ is as close to zero as possible. However, this two-step process, where $`F_R`$ would be minimized before each call to $`F_M`$, is computationally rather intractable in practice. Fortunately, it turns out that for the degree of accuracy that is sensible to demand from the MIT bag model (set by $`F_M`$ of the original fit ), it is sufficient to perform the refitting by varying simultaneously $`Z_0,\alpha _c,B,n`$ and the set $`\{R_{}\}`$ to minimize the joint functional $`F`$. Nevertheless, one should accept only those minimizations where $`F_M`$ is the overwhelming share, and $`F_R`$ only a small part of $`F=F_M+F_R`$; otherwise the fit to the experimental masses would be done away from the equilibrium bag radii. Ideally, the aim would be $`F=0`$, but since it is not possible to model all experimental masses exactly, one should look for such parameter values for which $`F`$ is sufficiently small. In the present model it is sensible to demand $`F<3\times 10^3`$ GeV2, since the original MIT bag fit  gives $`F_M=3.2\times 10^3`$ GeV2.

The minimization of $`F`$ by the simplex method , which had already been proved as robust and reliable in earlier applications , has turned out to be very suitable also in the present case.

## 3 Results and discussion

It is necessary to give some thought as to which outputs of the minimization procedure can be accepted as solutions to our problem. The present situation is different than in ref. , where the frozen radius approximation reduced the problem to solving a set of linear equations, so that the solution was unique once we chose which baryon masses would be used to fix the parameters. In the present case, the functional minimization finds many local minima of the functional (5). In which of them the minimization will end up, depends on which part of the parameter space one starts from. Moreover, many of these minima can be acceptable in the sense of sufficiently small value of the minimized functional $`F`$. We thus face the problem of non-uniqueness of the solutions. Fortunately, the smallness of the minimized functional $`F`$ is not the only criterion; clearly, a fit resulting in a good mass spectrum would anyway be unacceptable if it also resulted in physically unacceptable values of the bag radii or fitting parameters. This must always be kept in mind, as the problem was mathematically posed in such a way that it is possible to get an excellent fit to the masses, but with the bag radii and parameters devoid of any physical justification.

### 3.1 Practically instantonless bag

The first thing to check is the limit of the vanishing instanton density $`n`$ inside the bag. This is basically the case of the pure MIT bag model except that our model fitting is done by minimizing the functional (5). This is different from the original MIT bag fitting procedure  where the parameters $`B,Z_0`$ and $`\alpha _c`$ were fixed by singling out three hadrons and constraining their model masses to be the experimental ones. From the model standpoint, it is very satisfying that for $`n=0`$ inside the bag, our different fitting procedure leads to the description of baryons very similar to the original MIT bag fitting procedure . In addition to that, we note that when we depart from the limit of vanishing instanton density inside the bag and finally allow $`n0`$, the minimization of the functional (5) leads, among various outcomes, also to several solutions where the values of $`n`$ are nonvanishing but extremely small, $`n\text{}10^6`$ GeV4. This is practically negligible in comparison with the density $`n_01.6\times 10^3`$ GeV4 = 1 fm-4 estimated reliably (e.g., see refs. ) for the nonperturbative QCD vacuum outside the bag. The resulting baryon masses and radii (as well as values of the variable model parameters) are very close to each other in all these cases of very small $`n`$, and also very similar to the pure MIT case ($`n=0`$), as one would expect, although all those cases are formally different solutions. It is thus clear that they all describe very similarly (“practically uniquely”) the situations when $`n0`$. This shows that non-uniqueness of the solutions is not a problem in practice, and the same happens in the more interesting cases with significant values of $`n`$, discussed in the next subsection.

### 3.2 Appreciable instanton density inside the bag

Let us now discuss the first major interest of this paper: the cases when instanton densities $`n`$ inside the bag are significantly different from zero. Indeed, in most cases we obtained interesting solutions where densities $`n`$ inside the bag are an order of magnitude larger than in the linearized approximation , where333The dimensionless density $`\stackrel{~}{n}`$ used in  and $`n`$ are related by $`n\stackrel{~}{n}\rho ^4`$, where $`\rho `$ is the average instanton radius. Throughout, we have adopted the standard value $`\rho =1/600`$ MeV$`{}_{}{}^{1}1/3`$ fm (e.g., see refs. ). $`n=0.26610^4`$ GeV4. However, for all acceptable minimizations of the functional (5), we find they are still appreciably lower (at least by the factor of 3 or more) than the nonperturbative vacuum density $`n_01`$ fm$`{}_{}{}^{4}=1.610^3`$ GeV4. Therefore, we do not get a description of baryons which would be drastically different from the original MIT bag one , but we do obtain the desirable decrease of $`\alpha _c`$ which is noticeably stronger than the corresponding decrease obtained earlier in the linearized approximation . (In the case depicted in Table 1, $`\alpha _c`$ is by 30% smaller than in ref. .) In the solutions with decreased $`\alpha _c`$, we also observe the decrease of the baryon radii $`\{R_{}\}`$. As mentioned above, this is very desirable from the standpoint of nuclear physics, as explained by, e.g., Brown et al. . Namely, standard nuclear physics descriptions favor the picture of nuclei as made of independent nucleons interacting by effective boson exchange, but the empirical sizes of nuclei indicate that the “standard” MIT nucleon bags with $`R_N1`$ fm are already somewhat too large  for that. For this reason, we give in Table 1 the case with the smallest nucleon radius for which we managed to achieve an acceptable fit. Other physically acceptable solutions have somewhat smaller $`n`$ and somewhat larger radii. Table 2 gives a kind of condensed overview of several representative fits; e.g., the last line in Table 2 summarizes Table 1, the case with the highest $`n`$ which leads to a fit acceptable by all criteria. The general features of the acceptable fits are the following:

a) The values of the functional $`F`$ are around 1.3 to $`1.2\times 10^3`$ GeV2 (out of which only less than a percent is $`F_R`$). This gives the rough limit on the accuracy of reproduction of the mass spectrum within the present model. The average deviation from an experimental baryon mass is 11 MeV. In fact, the predictions for the masses of $`N`$ and $`\mathrm{\Sigma }`$ are the worst. They are too high for $`N`$ and too low for $`\mathrm{\Sigma }`$ by some twenty MeV. The other masses are within 10 MeV from the experimental masses. (In the MIT fit , the $`N`$ mass belongs to those constrained to experimental values to fix the model parameters, but then the $`\mathrm{\Sigma }`$ mass is too low by 45 MeV.) Overall, our approach to fitting of baryon masses gives noticeably smaller sum of squared deviations from the empirical baryon masses, $`F_M`$, than the original MIT bag fit  and the linearized approximation (where $`F_M=3.5\times 10^3`$ GeV2) .

b) Going beyond the linearized approximation and thereby allowing the bag radii to vary, leads to some significant changes with respect to the results in linearized approximation. Notably, Table 1 shows the instanton contributions to baryon energies are an order of magnitude larger than in the linearized approximation . Such instanton contributions are present not only in Table 1, but in large majority of fits, since the instanton densities in most of the presently obtained solutions are an order of magnitude larger than $`n`$ obtained in the linearized approximation . Nevertheless, since the instanton contributions to bag masses are still much smaller than other contributions (except $`E_E^{}`$), the general picture of baryons is not drastically altered with respect to the original MIT bag phenomenology .

c) In most cases this relatively large $`n`$ inside the bag leads to the decrease of $`\alpha _c`$, although there are also some fits with relatively large $`n`$ where $`\alpha _c`$ grows back close to its MIT value . Then, however, such a larger $`\alpha _c`$ is also accompanied by an excessive increase of bag radii. In particular, this yields a nucleon radius even larger than in the MIT case , so that such solutions must be discarded as unacceptable from the point of view of nuclear physics as explained above. The interdependence of the model parameters and the baryon bag radii which minimize $`F`$ is such that $`\alpha _c`$ decreases while $`Z_0`$ and $`B`$ increase with the decrease of the bag radii. This is illustrated in Table 2, which, for four different fits, displays the average octet ($`O`$) and decuplet ($`D`$) radii, $`\overline{R}^O`$ and $`\overline{R}^D`$, for four different fits. The notion of the average multiplet radii $`\overline{R}^O`$ and $`\overline{R}^D`$ is useful since the octet baryon radii are similar to each other, and the decuplet baryon radii are similar to each other. The decuplet radii are also some 10 % larger than the radii of the octet baryons.

### 3.3 Instanton-induced strangeness inside the bag

An inspection of ref.  easily shows that going beyond the linearized approximation, and the above effects thereof, do not change the results of ref.  on the vector, axial-vector and pseudo-scalar strangeness of the nucleon bag: the instanton-induced contributions to them are still vanishing.

In contrast to that, the instanton-induced scalar strangeness is enhanced an order of magnitude over what it was in the linearized approximation , following the increase of the instanton density $`n`$. This is seen in Table 3, which shows the dependence on the instanton density, or on the bag radius associated with this density, of various scalar strangeness components of the nucleon. The instanton-induced contributions due to $`_1^I`$ and $`_2^I`$, respectively denoted by $`N|\overline{s}s|N_{_1^I}`$ and $`N|\overline{s}s|N_{_2^I}`$, comprise the overwhelming share of the total instanton-induced contribution $`N|\overline{s}s|N__I`$. We do not display $`N|\overline{s}s|N_{_3^I}`$, the contribution due to $`_3^I`$, as it contributes only to the third decimal place.

Although our present interest are the instanton-induced contributions, we should also comment on the basic strangeness of the nucleon MIT bag, $`N|\overline{s}s|N_{\mathrm{basic}}`$ (found by Donoghue and Nappi ),

$$N|\overline{s}s|N_{\mathrm{basic}}(\eta 1)0|\overline{q}q|0\frac{4\pi }{3}R_N^3.$$
(8)
It is the product of the nucleon bag volume $`V_N=(4\pi /3)R_N^3`$ and $`0|\overline{q}q|0`$, the expectation value of the $`\overline{q}q`$ scalar condensate in the true, nonperturbative QCD vacuum, but also of the factor $`\eta 1`$ which has unfortunately remained quantitatively undetermined. Its determination is beyond the scope of the present paper. Let us just quote  that the $`\eta `$ ($`0<\eta <1`$) is in general some decreasing function of the bag radius, since $`R_N\mathrm{}`$ corresponds to $`\eta 0`$. The case $`\eta =0`$ is called the naive bag model limit and obviously maximizes the basic bag strangeness (8). This limit was, for definiteness, the only case of the basic bag strangeness $`N|\overline{s}s|N_{\mathrm{basic}}`$ considered in ref. . Although in the present paper even the $`\eta =0`$ limit of $`N|\overline{s}s|N_{\mathrm{basic}}`$ is not so much larger than (now increased) $`N|\overline{s}s|N__I`$ as was the case in ref. , it is still larger by an order of magnitude for all radii displayed in Table 3. The most widely accepted value of the condensate, adopted also in ref.  and the present paper, $`0|\overline{q}q|0=(240\mathrm{MeV})^3`$, leads to $`N|\overline{s}s|N_{\mathrm{basic}}`$ exceeding considerably the empirical value of the total scalar strangeness (determined by, e.g., ref. , from the $`\sigma `$-term estimated from the $`\pi N`$ scattering data and the masses of $`\mathrm{\Xi },\mathrm{\Sigma }`$ and $`\mathrm{\Lambda }`$),

$$N|\overline{s}s|N2.8.$$
(9)
Of course, lower values of $`0|\overline{q}q|0`$ trivially decrease $`N|\overline{s}s|N_{\mathrm{basic}}`$. For example, the choice $`0|\overline{q}q|0=(200\mathrm{MeV})^3`$, as in ref. , amounts to reducing values of $`N|\overline{s}s|N_{\mathrm{basic}}/(1\eta )`$ in Table 3 by the factor $`(200/240)^3=1/1.728`$, but this still gives rather large values. Instantons inside the bag help with that. Admittedly, Eq. (8) shows clearly that $`N|\overline{s}s|N_{\mathrm{basic}}`$ is not directly dependent on instantons and their density $`n`$ inside the bag, but there is an indirect connection: first, the volume factor in Eq. (8) decreases with radii as $`R_N^3`$, and diminishing radii are associated with increasing $`n`$. Second, $`1\eta `$ also falls with $`R_N`$. Thus, even if $`N|\overline{s}s|N_{\mathrm{basic}}`$ in the original, instantonless MIT bag model would be too close to its (too large) naive value, this potential problem would now be alleviated (more strongly than $`R_N^3`$) by lower values of $`R_N`$, occurring at higher $`n`$.

As already stressed, the most interesting effect of presently increased values of the instanton density $`n`$ is the considerable enhancement of the instanton-induced scalar strangeness, and we would like to point out that this enhancement is not due to a favorable choice of fixed input parameters $`m_u,m_d,m_s`$ and $`0|\overline{q}q|0`$. In fact, our adoption of the fixed input parameters of refs.  and  was motivated not only by easiness of comparison with these papers. This choice is also suitable for stressing that the present enhancement (of the instanton-induced strangeness) is not an effect of the choice of the model parameters. This is because the values of the quark mass parameters and of the vacuum quark-antiquark ($`\overline{q}q`$) scalar condensate used in ref.  and in obtaining all presently displayed results ($`m_u=m_d=0`$, $`m_s=279`$ MeV, $`0|\overline{q}q|0=(240\mathrm{MeV})^3`$), actually lead to smaller instanton-induced scalar strangeness than those used in ref.  ($`m_u=m_d=8`$ MeV, $`m_s=200`$ MeV, $`0|\overline{q}q|0=(200\mathrm{MeV})^3`$). That the latter choice  of these inputs gives (at a given instanton density $`n`$) even more enhanced instanton-induced scalar strangeness than that in Table 3, is most easily understood if one notes the role of the characteristic pre-factors, denoted by $`_f`$’s in refs. , appearing in the instanton-induced interaction $`_I`$.

The factor $`_f`$ pertaining to a flavor $`f`$ $`(f=u,d,s)`$, is composed of the corresponding quark mass parameter $`m_f`$, the average instanton size $`\rho \frac{1}{3}\mathrm{fm}`$ , and the $`\overline{q}q`$ condensate $`0|\overline{q}q|0`$, in the following way:

$$_f\frac{1}{m_f\rho \frac{2\pi ^2}{3}\rho ^30|\overline{q}q|0},(f=u,d,s).$$
(10)
Obviously, smaller values of $`m_u,m_d,m_s`$ and $`0|\overline{q}q|0`$ will increase $`_f`$’s, and vice versa.

Let us consider the concrete sets of input parameters, those of refs.  and . Changing $`m_u=m_d`$ from 0 to 8 MeV actually does not have significant influence on the instanton-induced strangeness, since $`m_u`$ and $`m_d`$ are anyway small at the hadronic mass scale (where 8 MeV can be approximated by 0). Nevertheless, the decrease of $`m_s`$ from 279 MeV to 200 MeV is quite important for further increasing the instanton-induced strangeness significantly over the values in Table 3. In fact, the effect thereof is comparable to the effect of the decrease of $`|0|\overline{q}q|0|`$ from $`(240\mathrm{MeV})^3`$ to $`(200\mathrm{MeV})^3`$.

Besides the effect on $`_f`$’s, the change of quark masses changes the quark wave functions, and, more importantly, the quark energy denominators appearing in the calculation of the nucleon strangeness (as can be seen in ref. ). This way, the decrease of $`m_s`$ from 279 MeV to 200 MeV increases still further the instanton-induced strangeness. (Again, the increase of $`m_u=m_d`$ from 0 to 8 MeV is too small to have a significant influence.) The effect of the quark energy denominators and wave functions is not so clearly disentangled as the effect of the $`_f`$–factors, so the explicit calculation is needed to show that the effect is of comparable magnitude. However, the important thing for the present discussion is that this effect changes the nucleon strangeness in the same direction as the $`_f`$–factors.

## 4 Conclusion

To summarize, we first remark that we did not perform a “first-principle”-type calculation of the probability of penetration of the instanton liquid from the surrounding nonperturbative QCD vacuum into the bag. Rather, we performed model fits to baryon masses and these fits showed which values of the instanton densities can be accommodated inside the MIT bag (in a physically acceptable way) and what the effects thereof would be. In the present paper, we went beyond the linearized approximation of ref. , and the bag radii were allowed to vary in the course of parameter fitting, which was performed so that the radii had to satisfy the pressure-balance condition. In this approach, the importance of the instanton-induced interaction allowed to act inside the quark bag is increased in every way, and not only for the baryon mass shifts, the size of which follows the increase of the instanton density inside the bag. We namely found that the MIT bag interior can accommodate instanton densities an order of magnitude larger than found in the linearized approximation . They grow faster than the inverse of the bag volume with decrease of the bag radii returned by the fitting procedure. The growth of the instanton-induced scalar strangeness of the nucleon is even slightly faster than that when the nucleon radius falls. The instanton-induced part of the nucleon strangeness is now an order of magnitude larger than the instanton-induced strangeness found in the linearized approximation . The quantity (8), considered as the basic MIT bag contribution to nucleon strangeness , remains undetermined also in the present work, but we could show that it must be smaller in the instanton-enriched MIT bag than in the original MIT bag. This is good, because this quantity alone has the potential to overshoot strongly the empirical value (9) of $`N|\overline{s}s|N`$. Also, it turns out that allowing for the possibility of instanton densities significantly different from zero inside the MIT bags now enables the favorable, up to 30 % decrease of $`\alpha _c`$. More importantly, it also enables the decrease of the nucleon MIT bag radius by more than 10 %, improving somewhat the consistency of the MIT bag model with nuclear physics.

After so summarizing the concrete results of this paper, we close by making a more speculative comment on how our enabling substantial instanton-liquid densities inside the MIT bag seems to improve the consistency of the model. Let us first note that, contrary to some earlier statements , Schäfer  has recently shown that the instanton liquid model, including the one we use, is not necessarily in conflict with the expansion in large $`N_c`$, the number of QCD colors. Then we recall an observation of Bardeen and Zakharov  concerning the $`N_c`$ scaling of the bag constant $`B`$. Inside hadrons, which can be modeled by the bag model, the quark color fields most probably suppresses instantons and other nonperturbative fluctuations. However, a smooth large-$`N_c`$ limit seems to indicate  that the suppression of these fluctuations, including instantons, is not very strong (in keeping with relatively low values of $`B`$ coming from phenomenological fits). Our modification of the bag, containing considerable instanton-liquid densities inside, is obviously more consistent with the Bardeen and Zakharov’s result  than is the original MIT model, where this suppression is complete.

## Acknowledgments

The support of the Croatian Ministry of Science and Technology contract 0119261 is gratefully acknowledged.

| Baryon $``$ | $`M_{exp}^{}`$ | $`M_{bag}^{}`$ | $`R_{}`$ | $`E_0^{}`$ | $`E_V^{}`$ | $`E_Q^{}`$ | $`E_M^{}`$ | $`E_E^{}`$ | $`E_I^{}`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`N`$ | $`0.938`$ | $`0.959`$ | $`4.365`$ | $`0.526`$ | $`0.176`$ | $`1.404`$ | $`0.128`$ | $`0.000`$ | $`0.033`$ |
| $`\mathrm{\Lambda }`$ | $`1.116`$ | $`1.120`$ | $`4.383`$ | $`0.524`$ | $`0.178`$ | $`1.557`$ | $`0.127`$ | $`0.003`$ | $`0.033`$ |
| $`\mathrm{\Sigma }^+`$ | $`1.189`$ | $`1.172`$ | $`4.529`$ | $`0.507`$ | $`0.197`$ | $`1.513`$ | $`0.094`$ | $`0.003`$ | $`0.061`$ |
| $`\mathrm{\Xi }^0`$ | $`1.315`$ | $`1.306`$ | $`4.475`$ | $`0.513`$ | $`0.190`$ | $`1.688`$ | $`0.109`$ | $`0.003`$ | $`0.047`$ |
| $`\mathrm{\Delta }`$ | $`1.232`$ | $`1.248`$ | $`5.130`$ | $`0.448`$ | $`0.286`$ | $`1.195`$ | $`0.109`$ | $`0.000`$ | $`0.107`$ |
| $`\mathrm{\Sigma }^{}`$ | $`1.385`$ | $`1.388`$ | $`5.073`$ | $`0.453`$ | $`0.277`$ | $`1.370`$ | $`0.096`$ | $`0.003`$ | $`0.095`$ |
| $`\mathrm{\Xi }^{}`$ | $`1.533`$ | $`1.526`$ | $`5.027`$ | $`0.457`$ | $`0.269`$ | $`1.543`$ | $`0.084`$ | $`0.003`$ | $`0.083`$ |
| $`\mathrm{\Omega }^{}`$ | $`1.672`$ | $`1.661`$ | $`4.978`$ | $`0.461`$ | $`0.261`$ | $`1.716`$ | $`0.074`$ | $`0.000`$ | $`0.071`$ |
| $`B^{1/4}=0.150`$ GeV $`Z_0=2.296`$$`\alpha _c=0.394`$$`n=0.51210^3`$ GeV4 | | | | | | | | | |

1.5cm

Table 1: The fit for the input quark masses $`m_u=m_d=0`$, $`m_s=279`$ MeV and the quark-antiquark vacuum condensate $`0|\overline{q}q|0=(240\mathrm{MeV})^3`$. We display the separate energies $`E_X^{}`$ ($`X=0,V,Q,M,E,I`$) contributing to $`M_{bag}^{}`$, the mass of the baryon bag, to be compared with the corresponding experimental baryon mass $`M_{exp}^{}`$ in the first column. The output values of the bag model parameters $`B,Z_0,\alpha _c`$ and $`n`$ are given in the lowest part of Table 1. All the masses and energies are given in GeV, and the bag radii $`R_{}`$ in inverse GeV, while $`Z_0`$ and $`\alpha _c`$ are dimensionless.

| $`n\times 10^3`$ | $`\alpha _c`$ | $`B\times 10^4`$ | $`Z_0`$ | $`\overline{R}^O`$ | $`\overline{R}^D`$ | $`F\times 10^3`$ |
| --- | --- | --- | --- | --- | --- | --- |
| \[GeV4\] |  | \[GeV4\] |  | \[GeV-1\] | \[GeV-1\] | \[GeV2\] |
| $`0.290`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}0.485}`$ | $`4.031`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}1.865}`$ | $`5.0`$ | $`5.6`$ | $`1.14`$ |
| $`0.310`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}0.474}`$ | $`4.188`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}1.930}`$ | $`4.9`$ | $`5.4`$ | $`1.18`$ |
| $`0.398`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}0.437}`$ | $`4.612`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}2.114}`$ | $`4.7`$ | $`5.2`$ | $`1.26`$ |
| $`0.512`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}0.394}`$ | $`5.058`$ | $`\mathrm{\hspace{0.17em}\hspace{0.17em}2.296}`$ | $`4.4`$ | $`5.1`$ | $`1.29`$ |

1.5cm

Table 2: Brief overview of some typical fits. (The fit given in Table 1 is one example of them.) The values of functional $`F`$ in the last column shows good quality of the fits. The interdependence of the adjustable bag parameters ($`n`$, $`\alpha _c`$, $`B`$, $`Z_0`$) and the bag radii is summarily depicted utilizing the average octet and decuplet radii, $`\overline{R}^O`$ and $`\overline{R}^D`$, respectively.

| $`n\times 10^3`$ \[GeV4\] | $`R_N`$ \[GeV-1\] | $`N|\overline{s}s|N_{_1^I}`$ | $`N|\overline{s}s|N_{_2^I}`$ | $`N|\overline{s}s|N__I`$ | $`\frac{N\left|\overline{s}s\right|N_{\text{basic}}}{\left(1\eta \right)}`$ |
| --- | --- | --- | --- | --- | --- |
| 0.290 | 4.994 | 0.22 | 0.09 | 0.31 | 7.21 |
| 0.310 | 4.909 | 0.24 | 0.10 | 0.34 | 6.85 |
| 0.398 | 4.658 | 0.29 | 0.15 | 0.44 | 5.85 |
| 0.512 | 4.365 | 0.36 | 0.22 | 0.58 | 4.82 |

Table 3: Dependence of the scalar strangeness of the nucleon on the instanton density $`n`$, or on the bag radius $`R_N`$ associated with this density. For comparison, we recall that in the linearized approximation  the instanton-induced strangeness from $`_1^I`$ was $`N|\overline{s}s|N_{_1^I}=0.035`$, whereas the contribution from $`_2^I`$ was $`N|\overline{s}s|N_{_2^I}=0.023`$ (at $`R_N=5.00`$ GeV-1$`1`$ fm). In the last column, the choice $`\eta =0`$ maximizes the basic bag strangeness contribution of ref. ; this is nevertheless only the so-called naive bag model limit, and in fact $`\eta `$ remains undetermined. The fixed model inputs ($`m_u,m_d,m_s`$ and $`0|\overline{q}q|0`$) are the same as in Table 1 and Table 2, and are discussed in detail in the main text.


text:
# A model differential equation for turbulence

## I Model Equation

In this letter, we present a model equation for the energy spectrum of isotropic homogeneous turbulence,

$$E_t=\frac{1}{8}(k^{11/2}E^{1/2}(E/k^2)_k)_k+f\nu k^2E,$$
(1)
where $`t`$ is time, $`k`$ is the absolute value of the wavenumber, $`\nu `$ is the kinematic viscosity coefficient, $`f`$ is an external forcing and the energy spectrum $`E(k,t)`$ is normalised so that the kinetic energy density is $`E𝑑k`$. The third term on the RHS of this equation is obvious and describes the linear process of viscous dissipation. The first term contains a nontrivial model of the nonlinear processes in turbulence which rests on three basic properties111The parts 1 and 2 were inspired by paper by Kulsrud and Anderson Kulsrud and Anderson (1992). However, they used a 1st order in $`k`$ equation which does not allow thermodynamic equilibria and which fails to describe spreading of the spectrum toward lower $`k`$.:

1. The total energy density $`E𝑑k`$ is conserved by the inviscid dynamics. The characteristic time of the spectral energy redistribution for local interaction of turbulent scales is of order of the vortex turnover time, $`1/\sqrt{k^3E}`$.

2. The steady state in forced turbulence corresponds to a constant energy cascade though the inertial range of scales which is described by the Kolmogorov spectrum,

$$E=CP^{2/3}k^{5/3},$$
(2)
where $`P`$ is the energy flux (constant in $`t`$ and $`k`$) and $`C`$ is the Kolmogorov constant. Experimental measurements Sreenivasan (1995) give $`C=1.6\pm 0.17`$. As we will see below, equation (1) has an exact solution of form (2) with $`C=(24/11)^{2/3}1.68`$ ($`C`$ can be changed by tuning the constant factor in the first term).

3. When the wave-number range is truncated to a finite value and both forcing and dissipation are absent, turbulence reaches a thermodynamic equilibrium state characterized by equipartition of energy over the wave-number space Kraichnan and Montgomery (1980). In terms of the one-dimensional energy spectrum this means $`Ek^2`$ which is obviously a steady state solution of the equation (1) for $`f=\nu =0`$. Mathematical simplicity and respect to the above basic properties of Navier-Stokes turbulence make model (1) useful for practical and numerical applications. We will now analyze solutions of (1) in greater detail to find other properties that are predicted by this model.

## II Stationary Solutions

Let us consider steady-state spectra in the inertial range. For $`f=\nu =0`$, we have the following general time-independent solution,

$$E=Ck^2(Pk^{11/2}+Q)^{2/3},$$
(3)
where $`C=(24/11)^{2/3}1.68`$ and $`P`$ and $`Q`$ are arbitrary constants. For $`Q=0`$, this gives the pure Kolmogorov cascade solution (2), whereas for $`P=0`$ this is a pure thermodynamic equilibrium. For the general solution, both the constant flux of energy $`P=\frac{1}{8}k^{11/2}E^{1/2}(E/k^2)_k0`$ and a thermodynamic part ($`Q0`$) are present; they appear as an interesting nonlinear combination and not just as a linear superposition because equation (1) is nonlinear. Thus, one can refer to solution (3) with finite $`P`$ and $`Q`$ as a warm cascade to distinguish it from the pure Kolmogorov solution which could be viewed as a cold cascade.

Let us suppose that turbulence is produced near some scale $`k=k_0`$ so that $`f(k)>0`$ only in a finite range in the vicinity of $`k_0`$. Suppose that $`f(k)=0`$ to the left of this range (at large scales) and in a large inertial range to the right which ends at a very high $`kk_d`$ where viscosity $`\nu `$ or some other dissipation mechanisms ($`f(k)<0`$) become important. Then, up-scale of the forcing there will be a pure thermodynamic solution with $`P=0`$ and $`Q0`$ because there is no dissipation or forcing assumed to be present near $`k=0`$ to absorb or generate a finite energy flux. In the inertial range down-scale of the forcing there will be a constant flux cascade solution. This solution typically takes the form of a pure Kolmogorov (cold) cascade and extends down to the dissipation range where the energy flux is absorbed. Typically, the solution only penetrates a finite distance into the dissipation range and adapts itself until it provides sufficient dissipation to absorb the supplied flux. In the presence of dissipation, the model does not develop structure at arbitrarily high $`k`$ as it would in the inviscid case. The qualitative features of the steady state are independent of the detailed form chosen for the dissipation. Figure 1 shows the steady state solutions obtained numerically for several different choices of the dissipation.

However, if the dissipation is not sufficiently strong, the solution can penetrate far enough into the dissipation range to reach the maximal wave-number which necessarily exists in any numerical solution. If one imposes a zero flux condition at the right end of the computational interval, the energy flux is reflected from the maximal wave-number leading to greater values of $`E`$ in the dissipative range. Such a cascade stagnation acts to enhance the dissipation rate and thereby to adjust it to the energy flux to be absorbed. This phenomenon is common in numerical simulations of turbulence and is usually called the bottleneck phenomenon Falkovich (1994). Figure 2 shows a numerically obtained steady state for the dissipation function,

$`\nu (k)`$ $`=`$ $`\nu _0(kk_D)^2k>k_D`$
$`0k<k_D`$
with $`\nu _0=1.0\times 10^5`$ and $`k_D=500`$. The bottleneck phenomenon is clearly seen as an energy “pile up” over the cold cascade solution near the dissipative scale. In our model, the bottleneck phenomenon is described by the warm cascade solutions; in particular the theoretical curve in figure 2 is computed by taking $`P14.5`$, $`Q1.5\times 10^9`$ in equation (3). The relative importance of the “thermal” effects with respect to the cascade grows as one moves from larger to smaller scales; in extreme cases the spectrum can be nearly pure Kolmogorov near the forcing range and almost purely thermodynamic near the dissipative scale.

Because of the nonlinearity, solutions for given forms of forcing and dissipation are usually hard to find analytically and one has to use numerics. However, a lot of insight about the qualitative behavior of the system can be gained from considering stationary solutions (3) in an inertial range $`k_1<k<k_2`$ and fixing the spectrum at its boundaries, $`E(k_1)=E_1,E(k_2)=E_2`$. This kind of the boundary conditions roughly models the forcing and the dissipation effects outside of the inertial range. This gives

$`P=[(E_2/Ck_2^2)^{3/2}(E_1/Ck_1^2)^{3/2}]/(k_2^{11/2}k_1^{11/2})`$ (4)
$`Q=[k_2^{5/2}(E_2/C)^{3/2}k_1^{5/2}(E_1/C)^{3/2}]/(k_2^{11/2}k_1^{11/2})`$ (5)
Thus, the sign of $`P`$ is opposite to the sign of $`(E_2/E_1k_2^2/k_1^2)`$ and can be either positive or negative depending of the spectrum steepness with thermodynamic $`k^2`$ solution been a borderline case for which $`P=0`$. Constant $`Q`$ can also be either positive or negative with Kolmogorov $`5/3`$ been the borderline slope. It is convenient to think of the $`Q<0`$ solutions as negative temperature states.

## III Non-stationary Solutions

So far we concentrated on the stationary solutions but how do these solutions form? We consider from now on the inviscid case. The Kolmogorov $`5/3`$ spectrum is of a finite capacity type in that it contains only a finite amount of energy at the high $`k`$ end. Let us take an initial condition which is compactly supported and force the system by imposing a constant flux boundary condition across the left end of the computational interval. Owing to the finite capacity of the Kolmogorov solution, an infinitely remote dissipative scale must be reached in a finite time. The solution has a nonlinear front at $`k=k_{}(t)`$ and this front accelerates explosively, reaching $`k=\mathrm{}`$ at a finite singular time which we shall denote by $`t_{}`$. We can equally consider the decaying case where an initial distribution of energy, compactly supported at large scales, is allowed to spread under the action of the nonlinearity without any external forcing. Provided there is sufficient energy in the initial condition, the right front still reaches $`k=\mathrm{}`$ within a finite time but in addition a second front propagates to the left spreading the spectrum to large scales. This second front does not exhibit singular behaviour in finite time. Hence the decaying case is broadly similar to the forced case as far as large $`k`$ behaviour is concerned since the initial concentration of energy at large scales acts as an effective forcing for the right front.

It is well known (see for example Lacey et al. (1982)) that the solutions of nonlinear diffusion equations with compactly supported initial data often have the property of remaining compactly supported during the time evolution. This turns out to be the case here. The left and right nonlinear fronts actually correspond to the left and right extrema of the support of the solution. Since these points must be determined as part of the solution, we are, in principle, required to solve a moving boundary problem with two free boundaries. In order that the problem be well-posed, we require an additional moving boundary condition on each interfaceElliott and Ockendon (1982). For the problem under consideration, the appropriate condition for the right front is

$$\frac{dk_{}^{(R)}}{dt}=\underset{kk_{}^{(R)}}{lim}k^{\frac{11}{2}}E^{\frac{1}{2}}\frac{}{k}\left(k^2E\right).$$
(6)
A similar condition holds for the left front.

This formula can be obtained as follows. Consider the time interval $`[0,t]`$ with $`t<t^{}`$. For $`0<\tau <t`$ the right interface is given by the curve $`k_{}^{(R)}(\tau )`$ as shown in figure 3. Taking, $`ϵ<<1`$ and assuming that the solution goes continuously to zero at the right interface we define the curve $`k_ϵ(\tau )`$ by the condition $`E(k_ϵ(\tau ),\tau )=ϵ`$ for $`0<\tau <t`$. We now integrate around the contour, $`\mathrm{\Gamma }`$, in $`(k,t)`$ plane as shown in the figure,

$$_\mathrm{\Gamma }\left(\frac{E}{T}+\frac{P}{k}\right)𝑑k𝑑\tau =0.$$
(7)
We obtain

$$_\mathrm{\Gamma }\left(Edk+Pd\tau \right)=0.$$
(8)
Proceed as follows

$$_\mathrm{\Gamma }E\left(dk+\stackrel{~}{P}d\tau \right)=0,$$
(9)
where

$$\stackrel{~}{P}=k^{\frac{11}{2}}E^{\frac{1}{2}}\frac{}{k}\left(k^2E\right)$$
(10)
$$_{k_ϵ(t)}^{k_{}^{(R)}(t)}E𝑑k+_{k_{}^{(R)}(\tau )}E\left(dk+\stackrel{~}{P}d\tau \right)+_{k_{}^{(R)}(0)}^{k_ϵ(0)}E𝑑k+_{k_ϵ(\tau )}E\left(dk+\stackrel{~}{P}d\tau \right)=0.$$
(11)
Now $`E=0`$ on $`k_{}^{(R)}(\tau )`$ and $`E=ϵ`$ on $`k_ϵ(\tau )`$. Using these facts we obtain

$$k_ϵ(t)k_ϵ(0)+_{k_ϵ}\stackrel{~}{P}𝑑\tau =_{k_ϵ(0)}^{k_{}^{(R)}(0)}\frac{E}{ϵ}𝑑k_{k_ϵ(t)}^{k_{}^{(R)}(t)}\frac{E}{ϵ}𝑑k.$$
(12)
We now take $`ϵ0`$ so that $`k_ϵ(\tau )k_{}^{(R)}(\tau )`$. The integrands on the RHS are bounded by 1 and therefore give no contribution. We are left with

$$k_{}^{(R)}(t)=k_{}^{(R)}(0)_0^t\stackrel{~}{P}(k_{}^{(R)}(\tau ),\tau )𝑑\tau ,$$
(13)
which in conjunction with (10) yields (6) upon differentiation.

Let us now look for a self similar solution of equation (1) taking the following form

$$E=(t_{}t)^aF(\eta );\eta =k/k_{},k_{}=c(t_{}t)^b,$$
(14)
where $`a,b`$ and $`c`$ are constants. Clearly, $`b`$ must be negative since we require that $`k_{}\mathrm{}`$ as $`tt_{}`$. Substituting (14) into (1) with $`f=\nu =0`$ we find that the $`t`$ dependence drops out of the equation if $`a=23b.`$ We then have the following equation for $`F`$,

$$(3b+2)F+b\eta F^{}=\frac{C^{3/2}}{8}(\eta ^{11/2}F^{1/2}(F/k^2)^{})^{},$$
(15)
where prime means differentiation with respect to $`\eta `$. Equation (15) defines a one-parameter family of self-similar solutions. Note that substitution of the form (14) into the moving boundary condition, (6) yields the same similarity relations. Thus the moving boundary condition is consistent with the similarity ansatz but does not provide any additional constraints. In particular it does not tell us which member of the 1 parameter family of solutions is selected by the PDE.

The solution near the front tip can be found by expanding $`F`$ in series with respect to small $`(1\eta )`$; in the leading order we have

$$F=\frac{16b^2}{C^3}(1\eta )^2$$
(16)
which gives for the spectrum

$$E=\frac{16b^2}{k_{}^{}{}_{}{}^{3}(t_{}t)^2}\left(1\frac{k}{k_{}}\right)^2.$$
(17)
Notice that the quadratic decay of the solution as one approaches the tip is exactly the decay required in order to ensure that the speed of the tip as given by (6) is finite and non-zero.

We are interested in solutions which behave like a power law far behind the front. That is $`Ek^x`$ as $`k0`$. The relations (14) then imply that $`x=a/b`$. The pure Kolmogorov spectrum, $`x=5/3`$, therefore requires $`b=3/2`$, corresponding to what one might consider to be normal scaling in the wake of the front.

We performed numerical simulations of the forced solutions of equation (1) with compact initial data and constant input flux on the left. We use the numerical solution to check the development of a self-similar front with tip of form (16) and to determine which value of $`b`$ is selected. The results are shown in figure figure 4. The scaling parameter, $`b`$ and the singular time, $`t_{}`$, are most conveniently extracted from the relation

$$k_{}\left(\frac{dk_{}}{dt}\right)^1=\frac{1}{b}(t_{}t),$$
(18)
which allows one to calculate $`b`$ and $`t_{}`$ from a linear fit of the data near $`t=t_{}`$ as shown in figure 5. We find that $`t_{}=0.0799`$ and $`b=1.748`$ which corresponds a significantly steeper than Kolmogorov slope, $`x=1.856`$. The singular time, $`t_{}`$ depends on the choice of initial conditions but the anomalous scaling exponent does not. In particular we verified that the same value of $`b`$ is obtained for the decaying case. Such anomalous scaling behavior whereby the exponent of the solution in the wake of the nonlinear front is larger than the Kolmogorov value has been observed before. Examples include MHD wave turbulence Galtier et al. (2000) and optical turbulenceConnaughton et al. (2003). We expect that this anomaly is a general property of finite capacity systems rather than a property of our model and therefore may also be present in the Navier-Stokes equations.

## IV Origin of the Transient Spectrum

For the model (1) the origin of the anomaly can be traced to the question of existence of a solution of the similarity equation (15) which has the correct behaviour both for the wake and the tip. Written in terms of $`x`$ rather than $`b`$ and re-scaled to get rid of constants, the similarity equation is

$$\frac{2}{x3}\left(\eta \frac{dF}{d\eta }+xF\right)=\frac{d}{d\eta }\left(\eta ^{\frac{11}{2}}\sqrt{F}\frac{d}{d\eta }\left(\eta ^2F\right)\right).$$
(19)
We require that this equation have a solution which behaves as $`\eta ^x`$ as $`\eta 0`$ and also satisfies the boundary condition, (6), at the front tip, $`\eta 1`$. Such a solution is not typical and actually exists only for one value of $`x`$. In particular, such a solution does not exist for $`x=5/3`$, the Kolmogorov value. The structure of the problem can be studied by introducing a new independent variable, $`s=\mathrm{log}\eta `$ and a pair of dependent variables, $`f(s)`$, $`g(s)`$ defined by

$$F=25\eta ^3f^2,\frac{dF}{d\eta }=\frac{25}{3}\eta ^4fg.$$
(20)
The purpose of this transformation is to “autonomise” the equation222Simpler autonomising transformations can be certainly be found. The form of the transformation has been chosen in order to make more convenient the positions of the equilibrium points of the resulting dynamical system.. That is, we remove the explicit dependence on $`\eta `$ from the equation. When this is done, equation (19) is equivalent to the following autonomous first order system:

$`{\displaystyle \frac{df}{ds}}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(f+g\right)`$ (21)
$`f{\displaystyle \frac{dg}{ds}}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left(5f^2+6fg9g^2+{\displaystyle \frac{10}{x3}}\left(3f+xg\right)\right),`$
and methods of phase plane analysis can be applied. The associated dynamical system has three equilibria P1 $`=(0,0)`$, P2 $`=(0,10/3(x3))`$ and P3 $`=(1,1)`$. Note that P1 and P2 are singular points of the original equations, (21).

Let us now attempt to identify the trajectory in the $`(f,g)`$ plane describing the physical solution. We note that the required trajectory must remain in the quadrant $`f>0`$, $`g<0`$. The wake is at $`\eta 0`$, or $`s\mathrm{}`$. In the wake, $`F\eta ^x`$ and $`F^{}x\eta ^{x1}`$. It follows that $`f\eta ^{3x}`$ and $`gx\eta ^{3x}`$. The value of $`x`$ is less than 3 so both $`f`$ and $`g`$ go to zero as $`\eta 0`$. The wake is therefore at $`(0,0`$ and should be reached as $`s\mathrm{}`$.

Now consider the tip. We know that the front tip is at $`\eta =1`$, or $`s=0`$.Both $`F`$ and $`F^{}`$ go to zero as $`\eta 1`$ but it follows from (6) that the product $`F^{1/2}F^{}`$ remains finite. The front tip therefore lies at some finite point on the negative $`g`$ axis. In fact the tip lies at P2. To show this, we show that the moving boundary condition is satisfied there. Translating (6) into the self-similar variables we obtain

$$b\tau ^{b1}=\tau ^{\frac{5}{2}b+\frac{1}{2}a}\underset{\eta 1}{lim}\eta ^{\frac{7}{2}}\frac{1}{\sqrt{F}}\frac{dF}{d\eta }.$$
(22)
Balancing the $`\tau `$ dependence gives us back the self-similarity condition, $`a=23b`$. We now carefully do the various changes of variables and rescalings to obtain

$$b=\underset{s0}{lim}\frac{3}{5}h.$$
(23)
Noting the similarity relations (14), we see that this condition is satisfied identically at the point $`P_2`$.

The critical value of $`x`$ occurs when the unstable manifold of $`P_1`$ intersects the stable manifold of $`P_2`$ thus forming a connection between the two singular points of the original equations. This is then the only trajectory which can satisfy the required conditions to describe both the tip and the wake. A numerical approximation to the unstable manifold of the origin as $`x`$ is varied through the critical value is shown in figure 6. The corresponding front profiles, converted back into the self-similar variables are shown in figure 7.

In practice, any model should include dissipation so that this self-similar solution above will be valid only until the front tip meets the dissipation scale. After this, the transient slope gets replaced in the inertial range by the stationary cascade solution, with or without bottleneck depending on the dissipation, as discussed above.

In summary, we presented a simple model in which basic properties of the Navier-Stokes turbulence are built-in: dimensionality and scaling, Kolmogorov and thermodynamic spectra. This model allows to obtain new predictions about the “warm” cascade states which are a nonlinear mixture of the cascade and thermodynamic solutions and which describe the bottleneck phenomenon. The model also allowed us to study the self-similar dynamics of the finite-time formation of the steady cascade states.

We thank Ildar Gabitov and Alan Newell for helpful discussions.