text:
# Ultra-high energy cosmic rays: are they isotropic?

## I Introduction

As ultra high energy cosmic rays (UHECR) continue to puzzle physicists and astronomers alike, the basic question of the isotropy of the flux beyond the unobserved GZK cut-off  retains a high priority. A recent analysis by the AGASA group  shows a trend towards isotropy at the highest energies, although some clusters of events (3 doublets and a triplet) are identified with a very low chance probability. This seeming contradiction is still more disturbing as both, bottom-up and top-down production mechanisms generally produce anisotropy to some degree in a natural way. In the present letter, we address the problem of isotropy of supra-GZK cosmic rays ($`E>4\times 10^{19}`$ eV) by using propagation simulations and one- and two-dimensional tests over simulated and existing world data (mainly from AGASA). The significance of the clusters of events detected by AGASA is also discussed in different scenarios.

## II Numerical method and results

Probably the conceptually simplest production models of UHECR are the ones involving bottom-up mechanisms. All of them, require that the sources of the particles group in more or less the same way as luminous matter does. Furthermore, on large scales, luminous matter trace roughly the distribution of cold dark matter (DM), although bias factors have to be taken into account. DM is involved in most of the top-down production mechanisms. It is therefore important to check the expected signature from such a source distribution (which we will call the LLMD - local luminous matter distribution - scenario).

We use a numerical simulation approach to track UHECR propagation through the intergalactic medium and evaluate their arrival distribution. The actual distribution of galaxies is used for the UHECR sources nearer than $`100`$ Mpc . Additionally, the same procedure as in  is used in the description of the intergalactic magnetic field (IGMF): a cell-like spatial structure, with cell size given by the correlation length, $`L_cB_{IGMF}^2(r)`$. The intensity of the IGMF, in turn, scales with luminous matter density, $`\rho _{gal}`$ as $`B_{IGMF}\rho _{gal}^{0.3}(r)`$  and the observed IGMF value at the Virgo cluster ($`10^7`$ G, ) is used as the normalization condition. Note, however, that the IGMF could be ordered and coherent on large scales , in which case the propagation of UHECR should be strongly model dependent . Test particles (protons) are injected at the sources with a spectrum $`dN/dEE^2`$ ($`E<10^{24}`$ eV) and propagated through the intergalactic magnetic field up to the detector on Earth. Energy losses due to redshift, pair production and photo-pion production due to interactions with the cosmic microwave background radiation (CMBR) are also included.

The Aitoff projection of the resultant (2D) arrival probability density is shown in Figure 1 in galactic coordinates. The mask covers the plane of the galaxy, where the actual distribution of galaxies is not well known due to obscuration by dust. The curved, thick line is the celestial equator. Superimposed on the figure are the available events with $`E>4\times 10^{19}`$ eV observed by AGASA (47 events ), Haverah Park (27 ), Yakutsk (24 ) and Volcano Ranch(6 ). The arrival probability contours trace quite well the local large scale structures. The supergalactic plane (SGP), in particular, can be easily distinguished running from North to South approximately along the $`l=135`$ galactic meridian. It is apparent from the figure that, despite some conspicuous clusters in the vicinity of the SGP, the actual observed distribution of UHECR is much more isotropic than what one would expect if their sources aggregate like the luminous matter. Unfortunately, given the non-uniform exposure in declination of the various experiments and the low number statistics involved, it is not trivial to quantify this statement.

The view of the AGASA group  is that supra-GZK events arrive isotropically at Earth. Nevertheless, to complicate things further, three doublets and a triplet within a separation angle of $`2.5^o`$ are also observed.

Independent tests are applied in order to confront AGASA data, at the same level of statistical significance, with two opposite yet plausible scenarios: a completely isotropic UHECR flux and a flux originated from sources that spatially map the large scale distribution of matter inside the GZK-sphere.

Two pools of particles, with 20 million protons each, were constructed: one strictly isotropic, the other obtained from the simulation results depicted in Figure 1. Independent samples are extracted from these reservoirs using the response in declination of the exposure of the AGASA experiment as a selection criteria . The size of each individual sample is equal to the number of events (47) actually observed by AGASA above $`4\times 10^{19}`$ eV.

The most elemental analysis that can be made regarding isotropy is one-dimensional, in right ascension (RA), where other complicating factors like non-uniform exposure in declination and low number statistics are more easily dealt with. Figures 2a and 2b show different forms of visualizing the distribution of events in RA. The shaded bands (figure 2a) in the background correspond to the 68% and 95% confidence levels of the expected distribution of events in RA for a sample of size 47 protons originated in the LLMD scenario. Despite the small size of the samples some features are clearly seen. The largest peak is the signal from the Virgo-Coma line of sight towards the North galactic pole. The opposite half of the SGP (towards the second quadrant in galactic latitude) is responsible for the smaller peak around 30o. The deep depressions surrounding the Virgo peak correspond to the Orion (left) and Local (right) voids, the most prominent structures in our immediate neighborhood, combined with the spurious effect of obscuration of the galaxy distribution due to the galactic plane.

The thick continuous lines in the same figure, correspond to the 68% and 95% confidence levels of the distribution in RA when the incident UHECR flux is isotropic. We can see that, even with so few events, both limits should be distinguishable.

The heavy squares represent the AGASA data (same bin size as for the models above) and are consistent with an isotropic distribution. No signature is seen from the Virgo peak and, furthermore, the most populated bins fall in a region corresponding to the Local Void.

A more quantitative treatment to characterize the anisotropy in RA is the first harmonic analysis . Thus, given a data sample, the amplitude $`r_{1h}=\sqrt{a_{1h}^2+b_{1h}^2}`$ and phase $`\mathrm{\Psi }_{1h}=\text{tan}^1\left(b_{1h}/a_{1h}\right)`$ are calculated, where $`a_{1h}=\frac{2}{N}_{i=1}^N\mathrm{cos}\alpha _i`$, $`b_{1h}=\frac{2}{N}_{i=1}^N\mathrm{sin}\alpha _i`$ and $`\alpha _i`$ is the right ascension of an individual event.

$`r_{1h}`$ and $`\mathrm{\Psi }_{1h}`$ are calculated for $`10^3`$ samples drawn form the isotropic and anisotropic (LLMD) distributions and the results are shown in figure 2b with small dots and crosses respectively. Both cases are very well discriminated in the $`r_{1h}`$-$`\mathrm{\Psi }_{1h}`$ plane. The error box for the first harmonic of AGASA data (calculated by ) is also displayed (hatched region), and is completely consistent with an isotropic UHECR flux. Moreover, the AGASA result by itself, seems completely inconsistent with the LLMD scenario. However, when the phase and amplitudes obtained from other major experiments are considered (large, thick horizontal bars in Figure 2b for Haverah Park -HP- Volcano Ranch -VR- Yakutsk -YK; see ) the picture looks suggestively different, since all the phase observations are clustered inside the same quadrant in RA, covering the right wing of the Virgo peak. That is, despite the fact that every isolated measurement is consistent with isotropy, the observed phases seem to show a systematic enhancement in the direction of the interface between the SGP and the large adjacent Local void. It must also be noted that Haverah Park and Volcano Ranch data behave more like a transition between the isotropic and LLMD scenarios. Three out of four first harmonic phases (HP, YK and VR) include the North galactic pole within one S.D. level, while the forth (AGASA) include it within two S.D.. The exclusion of the observed UHECR events inside the obscuration band, $`b<10^o`$, changes the phase of the AGASA result by only $`6^o`$ (from $`258^o`$ to $`252^o`$) and, therefore, previous conclusions are unchanged by this effect.

Clearly, a two-dimensional analysis of the data would be highly desirably in order to answer questions as simple as whether the data is isotropic or unimodal. One way of doing this, given the small number of events involved and the non-uniformity of the distribution of events in declination due to experimental limitations, is to analyze the normalized eigenvalues $`\tau _1`$, $`\tau _2`$ and $`\tau _3`$ of the orientation matrix $`𝐓`$ of the data. Defining $`𝐓_{𝐢,𝐣}=𝚺_{𝐤=\mathrm{𝟏}}^𝐍𝐯_𝐢^𝐤𝐯_𝐣^𝐤`$, where $`𝐯^𝐤`$ are the $`N`$ unit vectors representing the data over the celestial sphere and assuming $`0\tau _1\tau _2\tau _31`$, the shape, $`\gamma =log_{10}(\tau _3/\tau _2)/log_{10}(\tau _2/\tau _1)`$, and the strength parameter, $`\zeta =log_{10}(\tau _3/\tau _1)`$, can be built . The shape criterion $`\gamma `$ is useful in discriminating girdle-type distributions from clustered distributions. The larger the value of $`\gamma `$ more clustered is the distribution. Uniform, nearly isotropic, distributions have $`\zeta 0`$. Because of the nature of the experimental setup, the observed distribution of UHECR is girdle in nature, regardless of the isotropicity of the UHECR flux. Therefore, in figure 3 we compare the results for $`10^3`$ isotropic (rhombes) and LLMD (circles) samples respectively with the AGASA sample in the $`\gamma `$-$`\zeta `$ plane. It can be seen that the isotropic and LLMD scenarios should be very well separated with the available data, albeit its smallness. AGASA data (thick cross), on the other hand, does not fit either of these scenarios, being an intermediate case.

Figure 4 shows the number of doublets with separation smaller than 2.5o obtained from $`10^4`$ samples of $`47`$ events each, drawn from isotropic and LLMD populations. AGASA observed 3 pairs, which is a large number ( 8% chance probability) for an isotropic UHECR flux, and a rather small number (but still inside the 68% C.L. at 13% probability) if compared with the average value of 5.5 pairs obtained for the anisotropic flux. The situation is analogous for triplets. One triplet was actually observed by AGASA, while $`1.3\pm 3.1`$ is expected for the LLMD model and $`0.02\pm 0.2`$ for the isotropic model. It should also be noted that the AGASA triplet C2 and pair C1 (actually a triplet if Haverah Park data is included), as well as the lower energy cluster BC2 fall on the SGP, on top of a maximum of the arrival probability , strengthening the case for an extragalactic origin inside this structure.

## III Discussion and implications

Different tests have been applied to the analysis of the isotropy of UHECR with $`E>4\times 10^{19}`$ eV. Test samples are drawn from both, an isotropic flux of particles and an anisotropic flux originated in sources with the local luminous matter distribution. Samples used to compare with AGASA data have a size of 47 events and are selected according to the same declination sensitivity as AGASA’s. Our results can be summarized as follows:

1) The comparison between arrival directions of the UHECR above $`4\times 10^{19}`$ eV and the expected arrival probability density, calculated under the assumption of UHECR sources that cluster as the luminous matter in the nearby universe, shows a remarkable degree of isotropy, despite a notorious tendency for clusters to appear on top of large scale structure signatures (Fig.1).

2) AGASA’s RA distribution is consistent with an isotropic distribution (Fig.2).

3) $`10^4`$ simulated experiments equivalent to AGASA show that a set as small as 47 UHECR is enough to separate extremely well isotropic and LLMD scenarios in the amplitude-phase plane; AGASA’s error box is completely consistent with isotropy and inconsistent with LLMD (Fig.2)

4) Nevertheless, the phases of HP, VR, YK and AGASA fall in the same quadrant in phase (Fig.2), which covers the interface between the SGP in the general Virgo direction and the adjacent Local void. The first three experiments include the North galactic pole inside 1 SD, and AGASA at 2 SD.

5) The comparison of AGASA data with $`10^4`$ simulated data sets from isotropic and LLMD fluxes on the $`\gamma `$-$`\zeta `$ (shape-strength) plane show that the former is an intermediate case, more clustered than isotropic samples but less than LLMD (Fig.3).

6) The number of pairs observed by AGASA is too large for an isotropic flux, but it is within the 68% CL for LLMD flux (Fig.4).

While a first order interpretation of the AGASA data certainly points to an isotropic flux of UHECR, consideration of the first harmonic analysis of other data sets and of two-dimensional tests over the AGASA data itself, as well as expected numbers of doublets for isotropic and anisotropic samples, point to a more complicated, intermediate picture with a certain degree of mixture of both limiting cases.

We can envisage at least three scenarios in which such a result could be obtained:

1) The sources involve bottom-up mechanisms associated with luminous matter but some of the events are scattered in the intergalactic medium such that we observe the composition of a diffuse and a direct component .

2) The sources involve bottom-up mechanisms associated with luminous matter but there is a large local magnetic structure, like a magnetized galactic wind, which isotropize the UHECR flux upon traversing the galactic halo . As the energy of the particles increases, and as long as they all have the same mass, the degree of isotropization should decrease making the galactic pole visible.

3) The sources involve top-down mechanisms associated with dark matter whose distribution roughly associates with the LLMD. In this case, the observed flux is the composition of an extragalactic component, whose signature is not very different from that of the LLMD, and a component originated in the halo of our own galaxy.  showed that, under general conditions, the halo component would dominate the extragalactic flux by at least two orders of magnitude. This is only true, however, in the unrealistic case of dark matter uniformly distributed in intergalactic space. Nevertheless, dark matter aggregates strongly and tends to be overabundant, by factors of $`10^2`$, in the center of galaxy clusters when compared to its abundance in the halos of isolated galaxies. It can therefore be shown that, in a sample of 47 events, and assuming Virgo as the only source of extragalactic events, 3-7 events should originate in Virgo and arrive inside a solid angle of approximately the size of the cluster. This could give rise to a slight anisotropy that correlates with the SGP when combined with the almost isotropic flux originated in a large galactic halo. Note, however, that the solid angle does not need to point exactly in the direction of Virgo, depending on the large scale structure of the intervening magnetic field.

Obviously, more high quality, high energy data from HiRes and Auger are badly needed.

Acknowledgments. The author benefited from discussions with Alan. A. Watson, Peter Biermann, Todor Stanev, Eun-Joo Ahn and Torsten Ensslin. This work is partially supported by the Brazilian agencies FAPESP and CNPq.


format: text

text:
# Money and Goldstone modes

## I Introduction

In classical equilibrium theory in economics , agents submit their demand-versus-price functions to a “central agent” who then determines the relative prices of goods and their allocation to individual agents. The absolute prices are not fixed, so the process does not determine the value of money, which merely enters as a fictitious quantity that facilitates the calculation of equilibrium. Thus, traditional equilibrium theory does not offer a fundamental explanation of money, perhaps the most essential quantity in a modern economy.

Money is essentially a dynamical phenomenon, since its role is intimately related to the temporal sequence of events. Suppose an agent has apples and wants oranges. He might have to sell his apples to another agent before he buys oranges from a third agent: hence money is needed for the transaction, supplying liquidity. It stores value between transactions. If everybody could get together at a common market-place, there would be no need for money.

A “search-theoretic” approach to monetary economics has recently been proposed, , in which agents may be either money traders, producers, or commodity traders. The agents randomly interact with each other, and they decide whether or not to trade based on “rational expectations” about the value of a transaction. After a transaction the agent changes into one of the two other types of agents. This theory has a steady state where money circulates. As in other equilibrium theories, this theory does not describe a dynamics leading to the steady state, of sufficient detail, say for one to simulate it if one would so desire. The value of money comes about by: (1) pre-defining the ratio of agents which are money traders; and (2) requiring each money-trader to spend all his money at each transaction. Actually, Trejos and Wright  found their monetary equilibrium to be unstable with respect to a small perturbation of the expected value of money, causing the system to slide to a stable fixed point where money has no value!

Within our picture  there is no restoring force for the value of money—which could in principle be anything, with no consequences for the actual physical transactions that takes place. The value of money represents a “continuous symmetry”. If, at some point, the value of money was globally redefined by a certain factor, this would have no consequences whatsoever. Thus, in order to arrive at a specific value of money, the continuous symmetry must be broken.

The principle of continuous symmetries, and, particularly, the spontaneous breaking of those symmetries, is hugely important in physics. Indeed, it is the basic principle behind our understanding of the properties of fundamental particles. It is also responsible for magnetic order and superconductivity. As an example, consider, for instance, a lattice of interacting atoms forming a crystal. The crystal’s physical properties, including its energy, are not affected by a uniform translation $`X`$ of all atoms. The position of the crystal represents a continuous symmetry. It could in principle be anywhere. Nevertheless, starting with an arbitrary position of the individual atoms, the dynamics of the atoms interacting with their neighbors fix the actual position of the entire crystal! We say that the continuous symmetry is “dynamically broken”. Note that if we shift the crystal there is no restoring force. Also, if the system is subjected to random noise, the lack of restoring force, expressed by the underlying continuous symmetry, leads to large positional fluctuations. The atoms exhibit “bounded rationality”: they simply adjust to the positions of their neighbors, wherever they might me. In order to calculate their dynamics one does not need any information of the ultimate equilibrium (if one exists) of the entire system.

In economics, the agents represent the atoms of the economy. The value that an agent assigns to money corresponds to the position of an atom. The neighbors represent the subsystem of agents with which an individual agent interact, in one way or another. The local actions of the agents determine the value of money for each agent at each instant of time. The average value of money, corresponding to the “center of mass” for the atomic chain, is fixed by the collective dynamics of all the agents.

Figure 1 illustrates the concept of continuous symmetries in the two cases. The physical properties of the chain of atoms are invariant under a uniform shift of all atoms. The economy is invariant under a uniform shift of all prices. The actual positions, including the “center of mass” in both cases are fixed by the dynamics.

It is hard to envision a picture of an economy which does not derive its properties starting from a reasonably realistic and “simulatable” model of the behavior of the individual agents. In the final analysis, even equilibrium theory must have a dynamical underpinning if, in any sense, it aims at modeling a dynamical world.

It is not so important for our considerations how the agents actually behave, i.e., to what extent their behavior might be considered rational, as long as they somehow adjust their behavior in response to their local environment. To be specific, we study a network of vendors and buyers, each of whom have a simple local optimization strategy. Whenever a transaction is considered, the agent must decide the value of the goods and services in question, or, equivalently, the value of money relatively to that of the goods and services he intends to buy or sell. He will associate that value to his money that he believes will maximize his utility. Thus, the value of money is a “strategic variable,” that the agent in principle is free to choose as he pleases. However, if he makes a poor choice he will loose utility. In any case, the dynamics of his action, together with those of other agents, will break the continuous symmetry.

For simplicity, we assume that agents are rather myopic: they have short memories, and they take into account only the properties of their “neighbors,” i.e., the agents with which they interact directly. They have no idea about what happens elsewhere in the economy.

Despite the bounded rationality of these agents, the economy self-organizes into an equilibrium state. Since we define the dynamics explicitly, we are, however, able to treat the nature of this relaxation to the equilibrium state, as well as the response of the system to perturbations, and to noise-induced fluctuations around the equilibrium. These phenomena are intimately related to the dynamics of the system, and cannot be discussed within any theory concerned only with the equilibrium situation. The value of money is fixed by a “bootstrap” process: agents benefit from assigning a specific value to money, despite this value’s global indeterminacy.

## II The Model

In our model, we consider $`N`$ agents, $`n=1,2,\mathrm{},N`$, placed on a one-dimensional lattice with periodic boundary conditions. This geometry is chosen in order to have a simple and specific way of defining who is interacting with whom. The geometry is not important for our general conclusions concerning the principles behind the fixation of prices.

We assume that agents cannot consume their own output, so in order to consume they have to trade, and in order to trade they need to produce. Each agent produces a quantity $`q_n`$, of one good, which is sold at a unit price $`p_n`$, to his left neighbor $`n1`$. He next buys and consumes one good from his neighbor to the right, who then buys the good of his right neighbor, etc., until all agents have made two transactions. This process is repeated indefinitely, say, once per day.

Our model is a simple extension of Jevons’  example of a three agent, three commodity economy with the failure of the double coincidence of wants, i.e., when only one member of a trading pair wants a good owned by the other. A way out of the paradox of no trade where there is gain to be obtained by all is to utilize a money desired by and held by all. Originally this was gold, but here we show that the system dynamics can attach value to “worthless” paper money.

For simplicity, all agents are given utility functions of the same form

$$u_n=c(q_n)+d(q_{n+1})+I_n(p_nq_np_{n+1}q_{n+1}).$$
(1)
The first term $`c`$, represents the agent’s cost, or displeasure, associated with producing $`q_n`$ units of the good he produces. The displeasure is an increasing function of $`q`$, and $`c`$ is convex, say because the agent gets tired. The second term $`d`$, is his utility of the good he can obtain from his neighbor. Its marginal utility is decreasing with $`q`$, so $`d`$ is concave.

An explicit example is chosen for illustration and analysis,

$$c(q_n)=aq_n^\alpha \text{ , }d(q_{n+1})=bq_{n+1}^\beta ,$$
(2)
where $`a`$ is $`\frac{1}{2}`$ and of dimension \[utility per (unit of $`q_n`$)α\], $`b`$ is 2 and of dimension \[utility per (unit of $`q_{n+1}`$)β\], and $`\alpha `$ and $`\beta `$ are chosen as 2 and $`\frac{1}{2}`$. The specific values of $`a`$, $`b`$, $`\alpha `$, and $`\beta `$ are not important for the general results, as long as $`c`$ remains convex and $`d`$ concave.

The last term represents the change in utility associated with the gain or loss of money after the two trades. Notice that the dimension of $`I_n`$ is \[utility per unit of currency\], i.e., the physical interpretation is the value of money.

Each agent has knowledge only (indirectly through the prices that were charged) about the utility functions of his two neighbors, as they appeared the day before. The agents are monopolistic, i.e., agent $`n`$ sets the price of his good, and agent $`n1`$ then decide how much $`q_n`$, he will buy at that price. This amount is then produced and sold—there is no excess production. The goal of each agent is to maximize his utility, by adjusting $`p_n`$ and $`q_{n+1}`$, while maintaining a constant amount of money. Money has value only as liquidity.

The agents aim to achieve a situation where the expenditures are balanced by the income:

$$p_nq_np_{n+1}q_{n+1}=0.$$
(3)
This condition is important for our considerations. The agent envisions an ongoing process, which will be repeated many times. In real life, he will eventually deal with very many agents. If he spends more money than he earns, he will eventually run out of money, or reach a ceiling for debt, causing loss of utility thereafter. If he uses less than he earns, he will be hoarding money, which gives no pleasure, but with a loss of utility. We assume that there is no preference for consumption now rather than later. Perhaps he will over-spend for a while, but eventually he will try to balance his budget.

Note that when the value of money is fixed at a global value, $`I_n`$ = $`I`$, the agents optimize their utility by charging a price

$$p=2^{\frac{1}{3}}I^1$$
(4)
and selling an amount

$$q=2^{\frac{2}{3}}$$
(5)
at that price. This is the monopolistic equilibrium. In general, of course, with a different set of strategies there is no reason that the dynamics should lead to the monopolistic equilibrium, or to any stationary state for that sake. Note that the resulting quantities, Eq. (5), are independent of the value of money, which thus represents a continuous symmetry. There is nothing in the equations that fixes the value of money and the prices. Mathematically, the continuous symmetry expresses the fact that the equations for the quantities are “homogeneous of order one.” The number of equations is one less than the number of unknowns, leaving the value of money undetermined. This is the continuous symmetry which eventually is broken by the dynamics.

Agent $`n`$ tries to achieve his goal, of maximizing his utility and keeping constant his amount of money, by estimating the amount of goods $`q_n`$ his neighbor will order at a given price, and the price $`p_{n+1}`$ his other neighbor will charge at the subsequent transaction. Knowing that his neighbors are rational beings like himself, he is able to deduce the functional relationship between the price $`p_n`$, he demands and the amount of goods $`q_n`$, that will be ordered in response to this. Furthermore, he is able to make a guess at the size of $`p_{n+1}`$, based on the previous transaction with his right neighbor. This enables him to decide what the perceived value of money should be, and hence how much he should buy and what his price should be. This process is then continued indefinitely, at times $`\tau =1,2,3,\mathrm{}`$

This defines the game. The strategy we investigate contains the assumption made by each agent that the other agents do not change their valuation of money $`I`$, between their two daily transactions, and hence they maximize their utility accordingly.

The process is initiated by choosing some initial values for the $`I`$s. They could, e.g., be related to some former gold standard. In fixing his price at his first transaction of day $`\tau `$, agent $`n`$ exploits the knowledge he has of his neighbors’ utility functions, i.e., he knows that the agent to the left will maximize his function with respect to $`q_{n,\tau }`$

$$\frac{u_{n1,\tau }}{q_{n,\tau }}=0.$$
(6)
His knowledge of the functional relationship, between the amount of goods $`q_{n1,\tau }`$, ordered by agent $`n1`$ at time $`\tau `$ and the price $`p_{n,\tau }`$, set by agent $`n`$, allows agent $`n`$ to gauge the effect of his price policy. Lacking knowledge about the value of $`I_{n1,\tau }`$, agent $`n`$ instead uses the value from the last transaction $`I_{n1,\tau 1}`$. Substituting the expression for $`q_{n,\tau }`$, found from solving Eq. (6), into Eq. (1) and maximizing with respect to $`p_n`$ and $`q_{n+1}`$ we find (for details see )

$$I_{n,\tau +1}=\left(I_{n1,\tau }^4I_{n,\tau }^2I_{n+1,\tau }\right)^{\frac{1}{7}}$$
(7)
which is a weighted geometric average of the value the two neighbors and agent $`n`$ prescribed to their money the previous day. Using this value of $`I_n`$, agent $`n`$ can fix his price $`p_n`$ and decide which quantity $`q_{n+1}`$, he should optimally buy. This simple equation completely specifies the dynamics of our model. The entire strategy can be reduced to an update scheme involving only the value of money—everything else follows from this. Thus, the value of money can be considered the basic strategic variable!

Even though there is no utility in the possession of money, as explicitly expressed by Eq. (3), the strategies and dynamics of the model nevertheless leads to a real utility, given by the third term in Eq. (1), being assigned to money. The dynamics in this model is driven by the need of the agents to make estimates about future transactions. In a sense, this models the real world where agents are forced to make plans about the future, based on knowledge about the past—and, in practise, only a very limited part of the past.

In the steady state, where the homogeneity of the utility functions give $`I_n=I_{n+1}`$, we retrieve the monopolistic equilibrium equations, Eqs. (4) and (5). Of course, this is somehow accidental. Other more or less rational strategies will lead to different values of money.

### A Stability analysis

Taking the logarithm and introducing $`h_{n,\tau }=\mathrm{ln}(I_{n,\tau })`$ yields the linear equation:

$$h_{n,\tau +1}=\frac{4}{7}h_{n1,\tau }+\frac{2}{7}h_{n,\tau }+\frac{1}{7}h_{n+1,\tau }.$$
(8)
This describes a Markov process. Clearly, any $`h=constant`$ solves this equation, in accordance with the fact that $`h`$ is a continuous symmetry. Thus, in order to get solutions, one must specify initial conditions. A simulation of 1000 agents with random initial values for the variable $`h`$ (sampled from a uniform distribution on the interval ) is shown below for illustration.

Now assume that $`h_{n,\tau }`$ is a slowly varying function of $`(n,\tau )`$ and that we may think of it as the value of a differentiable function $`h(x,t)`$ in $`(x,t)=(n\delta x,\tau \delta t)`$. Then, expanding to first order in $`\delta t`$ and second order in $`\delta x`$, we find the diffusion equation

$$\frac{h(x,t)}{t}=D\frac{^2h(x,t)}{x^2}v\frac{h(x,t)}{x},$$
(9)
where $`D=\frac{5}{14}\frac{(\delta x)^2}{\delta t}`$ and $`v=\frac{3}{7}\frac{\delta x}{\delta t}`$. The generator $`T`$, of infinitesimal time translations is defined by

$$\frac{h(x,t)}{t}=Th(x,t).$$
(10)
Taking the lattice Fourier transformation, the eigenvalues of $`T`$ are found to be $`\lambda _k=k^2Dikv`$, where the periodic boundary condition yields $`k=\frac{2\pi }{N}l`$; $`l=0,1,\mathrm{},N1`$. The damping time for each mode $`k`$, is given by $`t_k=(k^2D)^1`$. The only mode that is not dampened has $`k=0`$, and is the soft “Goldstone mode”  associated with the broken continuous symmetry with respect to a uniform shift of the logarithm of prices in the equilibrium:

All prices can be changed by a common factor, but the amount of goods traded will remain the same; as already noted by Marx  when discussing gold as the unit in which prices are expressed, i.e., the measure of value. The rest of the modes are all damped, and hence the system eventually relaxes to the steady state.

Thus, despite the myopic behavior of the agents, the system manages to organize itself into the equilibrium. But in contrast to equilibrium theory, we obtain the temporal relaxation rates towards the equilibrium, as well as specific absolute values for the individual prices. The value of money is fixed by the history of the dynamical process, i.e., by the initial condition combined with the actual strategies of the bounded rational agents.

Note that the steady states are all Nash equilibria: the agents can not change their behavior without loosing utility. In the physics analogy, it requires energy for the atoms to shift away from the equililibrium, as long as their neighbors remain fixed.

Thus, there is a continuous set of perfectly stable Nash equilibria. In contrast, the “search theoretical” approach  yields a single Nash equilibrium, with a well defined value of money. The anchoring of the value is imposed by requiring the pre-determined number of money traders, somewhat artificial and unrealistic, we believe, to spend all their money at each transaction. If the condition for an agent changes he will change his value of money. If he, say invents an easier way of production, he will lower his price, sell more, and buy more goods, the effect being deflation propagating through the system, as described by the solution to Eq. (9) for a delta-function initial condition . Likewise, if production becomes more difficult for a single agent, or if he increases his lust for consumption, this will cause an inflationary wave in the the whole system. For all agents, except the agent with changed production capability and his neighbors, the changes in utility are transient effects. In the steady state all these agents will produce and consume the same amount of goods as before the change, but at a different price level.

How much money is needed to run the economy? Since the agents merely seek to preserve a constant amount of money, it does not really matter how much money they have, as long as they have sufficient liquidity to bridge the interval between buying and selling. One could even imagine an economy with no real money at all, as long as the seller would provide sufficient credit. If, at some point, it were decided by the government that dollar bills and coins can not be used for payment, this would probably just amount to a nuisance. We would be forced to use credit and debit cards and the like, but everything would work just as before. One could easily model this situation. The money would merely be numbers on accounts.

## III Conclusion:

Here we considered a simple toy model with monopolistic agents. In general, economy deals with complicated heterogeneous networks of agents, with complicated links to one another, representing the particular “games” they play with one another. We submit that the general picture remains the same. At each trade, the agents evaluate the value of money, by analyzing their particular local situation, and act accordingly. The prices charged by the agents will be constrained by those of the interacting agents. The geometry and topology of realistic networks is much more complicated, of course, but the continuous symmetry remains. The simple equation Eq. (8) would be replaced by a much more complicated avereraging over many interacting neighbors. In any case, there is a continuous symmetry related to a global change of prices, which is broken by the agents’ actions.

It would be interesting to study the formation and stability of markets where very many distributed players are interested in the same goods, but not generally interacting directly with one another. Modifications of this network model may also provide a toy laboratory for the study of the effects of the introduction of the key financial features of credit and bankruptcy as well as the control problems posed by the governmental role in varying the money supply.

Here the underlying principle for assigning value to money was the dynamical breaking of the continuous symmetry, leading to large fluctuations in the system when subjected to noise. For stock markets, it is often assumed that there is a fundamental value to be discovered by the market. In reality, there could easily be so much uncertainty about this value, that there would be no real restoring force for price fluctuations within a wide range. Effectively, we would again have a continuous symmetry, and consequent large price fluctuations above and below what could be expected from a theory based on rational expectations.


text:
# REFERENCES

d=3 random field behavior near percolation

F. C. Montenegro

Departamento de Fisica, Universidade Federal de Pernambuco, 50670-901 Recife PE, Brasil

D. P. Belanger and Z. Slanič

Department of Physics, University of California, Santa Cruz, California 95064

J. A. Fernandez-Baca

Solid State Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-6393

## Abstract

The highly diluted antiferromagnet $`Mn_{0.35}Zn_{0.65}F_2`$ has been investigated by neutron scattering for $`H>0`$. A low-temperature ($`T<11`$ K), low-field ($`H<1T`$) pseudophase transition boundary separates a partially antiferromagnetically ordered phase from the paramagnetic one. For $`1<H<7`$ T at low temperatures, a region of antiferromagnetic order is field induced but is not enclosed within a transition boundary.

Diluted uniaxial antiferromagnets (AF) in applied fields are ideal random-field Ising model (RFIM) systems. $`Mn_{0.35}Zn_{0.65}F_2`$ is a three dimensional ($`d=3`$) system with a magnetic concentration $`x`$ close to the percolation threshold, $`x_p=0.25`$. Its behavior may be contrasted with $`Fe_xZn_{1x}F_2`$, which differs magnetically from $`Mn_xZn_{1x}F_2`$ only in the nature and strength of the anisotropy, which is dipolar in the latter and is dominated by a much larger single-ion anisotropy in the former. Samples of $`Fe_xZn_{1x}F_2`$ have been investigated with $`x=0.25`$, $`0.27`$ and $`0.31`$. The first two do not order even with $`H=0`$ and show spin-glass-like behavior at all $`H`$. The $`x=0.31`$ sample shows AF long-range order (LRO) at low $`H`$ and spin-glass-like behavior at higher fields. In $`Mn_xZn_{1x}F_2`$, ac susceptibility measurements indicate a spin-glass-like clustering at low temperatures for samples with $`0.2<x<0.35`$. We present evidence suggesting that frozen spin-glass-like clusters also affect the stability of the AF long-range order in $`Mn_{0.35}Zn_{0.65}F_2`$.

The neutron scattering experiments were performed at the Oak Ridge National Laboratory High Flux Isotope Reactor using a two-axis configuration with a monochromated $`14.7`$ meV beam. Further experimental details have been previously reported. For simplicity, the transverse line shapes used in fits of the data for $`H>0`$ are of the mean-field RFIM form

$$S(q)=\frac{A}{q^2+\kappa ^2}+\frac{B}{(q^2+\kappa ^2)^2},$$
(1)
where $`\kappa `$ is the inverse correlation length for fluctuations. We do not fit data in the Bragg scattering region $`|q|<0.008`$ reciprocal lattice units (rlu).

Figure 1 shows AF (100) transverse scans obtained by heating with $`H>0`$ after cooling in zero field (ZFC) and the fits to Eq. 1 for $`H=0.35`$ T. A striking feature of these scans is the relatively small change in the line shape with temperature when compared with the $`H=0`$ scans reported previously, particularly for $`T>T_c(H)`$, where $`T_c(H)`$ is a transition-like temperature. The fits to the data are essentially identical whether $`B`$ is fixed to zero or allowed to vary, indicating that the squared-Lorentzian term in Eq. 1 is unimportant. The results for $`\kappa `$, which reflect the unusual line shape behavior, are shown in Fig. 2. The values for $`H=0.35`$ T drop well below the values for $`H=0`$ which are concentration-gradient limited, indicating that Eq. 1 does not describe the line shapes very well for $`H>0`$. The clear minimum in $`\kappa `$ vs. $`T`$ and Bragg scattering intensities vs. $`T`$ (not shown) do indicate that the system is trying to undergo a transition, but it is questionable that one successfully occurs. For $`H=0.75`$ T, when $`B`$ is fixed to zero, the values of $`\kappa `$ are artificially much smaller than the instrumental resolution. The results when $`B`$ is allowed to vary are those shown in Fig. 2. The minimum in $`\kappa `$ vs. $`T`$ is extremely shallow although a transition-like region is indicated by a peak in the scattering at $`|q|=0.008`$ rlu vs. $`T`$ as well as the rapid decline in the Bragg scattering intensity. Either a transition does not actually take place or very few spins are involved. From temperatures at which the maxima in the off-Bragg scattering at $`|q|=0.008`$ rlu and the maxima in the temperature derivative of the Bragg scattering with respect to temperature occur, the AF pseudophase boundary is determined. The boundary is shown by the points with horizontal error bars in the inset of Fig. 3. The points with vertical error bars indicate the magnetic field at which an intensity change is observed in the Bragg scattering and at which the peak of the off-Bragg scattering occurs as the field is changed while $`T`$ is kept constant. The transition region again appears quite broad. An example of this at $`T=5K`$ is shown in Fig. 3 which also indicates substantial hysteresis when comparing ZFC data and data obtained upon cooling in the field (FC) for both $`q=0`$ and $`|q|=0.008`$ rlu and is consistent with earlier magnetization measurements.

For $`2<H<7`$ T a region of AF LRO is indicated by the $`q=0`$ intensities in Fig. 4. The order is most intense at $`H5`$ T. There is no peak in the off-Bragg ($`|q|=0.008`$ rlu) scattering for $`T>5`$ K, as exemplified for $`H=5`$ T in the inset of Fig. 4. This indicates that there is no phase transition boundary associated with this AF order and that the AF order is most likely field induced. Within a cluster the sublattice with an excess of spins will preferentially order along the field direction. The other sublattice will then order in the opposite direction as a result of the exchange interaction. At fields above $`5`$ T, the AF order weakens. This contrasts the behavior observed in $`Mn_{0.5}Zn_{0.5}F_2`$, where a spin-flop phase with a clear transition boundary occurs above the AF one.

Previous magnetization measurements were made of the pseudotransition boundary. The inset in Fig. 3 shows the peak positions of $`d(MT)/dT`$ for $`H<0.5T`$ as the starred data points. The widths of these peaks are very large. For example, the half width at half maximum at $`H=0.35T`$ is more than $`1`$ K. These widths are much larger compared to measurements at larger $`x`$. Note that the positions do not coincide with the boundary determined with neutrons and even the curvature is different. For the neutron boundary, $`TT_c(H)H^{2/\varphi }`$ with $`\varphi =1.4`$ as in RFIM systems further from $`x_p`$. In contrast, the magnetization data are fit with $`\varphi =3.4`$, the typical spin-glass exponent. Above $`H=0.5`$ T, it was not possible to reliably determine a peak position. These results reinforce the scenario in which the boundaries do not represent a true transition to AF LRO, but rather an AF transition that is strongly interfered with by spin-glass-like frozen clusters. The transition at $`T=11`$ K for $`H=0`$, in contrast, seems much more normal, though spin-glass-like behavior is evident below $`T=7`$ K. Slow relaxation seems evident in both the neutron scattering and magnetization for $`H>0`$, probably a result of clustering induced by the alignment of domains by the field. At low $`H`$, hysteresis is observed well above $`T_c(H)`$ and clustering is indicated by a strong deviation from Curie-Weiss behavior below $`T=22`$ K. The small anisotropy in the $`Mn_xZn_{1x}F_2`$ system allows the field to align clusters more easily than in $`Fe_xZn_{1x}F_2`$ and this effect certainly must account for the differences between these two systems.

This work has been supported by DOE Grant No. DE-FG03-87ER45324 and by ORNL, which is managed by Lockheed Martin Energy Research Corp. for the U.S. DOE under contract number DE-AC05-96OR22464. One of us (F.C.M.) also acknowledges the support of CAPES, CNPq, FACEPE and FINEP (Brazilian agencies).


text:
# Exact bounds on the ground-state energy of the infinite-𝑈 Hubbard model

\[

## Abstract

We give upper and lower bounds for the ground-state energy of the infinite-$`U`$ Hubbard model. In two dimensions, using these bounds we are able to rule out the possibility of phase separation between the undoped-insulating state and an hole-rich state.

\]

The problem of phase separation (PS) in strongly correlated systems is one of the most debated subject, especially after the seminal paper by Emery and co-workers , who also pointed out its relevance in the framework of high-temperature superconductivity. However, after several years of intensive numerical and analytical investigation, there is no general consensus even for the qualitative features of the possible instability in the Hubbard and $`tJ`$ models .

In this work we focus on the Hubbard model on a square lattice of $`L`$ sites in the infinite-$`U`$ limit, and prove that there is no PS in the low doping limit: our result being one of the few rigorous statements in this field. The Hamiltonian is defined by

$$=t\underset{i,j,\sigma }{}𝒫_Gc_{i,\sigma }^{}c_{j,\sigma }𝒫_G,$$
(1)
where the $`i,j`$ denotes nearest neighbor sites, and the Gutzwiller projector $`𝒫_G`$ enforces the constraint of no double occupancy on each site.

An upper bound on the ground-state (GS) energy can be easily obtained by use of a variational wavefunction. A very simple choice is given by the fully polarized ferromagnetic (FM) state of energy $`E_{FM}=2t_𝐤_{i=1,d}\mathrm{cos}k_i`$, where $`d`$ is the dimensionality of the lattice and the sum over $`𝐤`$ is restricted to the lowest hole energy orbitals.

A lower bound to the GS energy is instead more involved and can be derived as follows. The GS of the Hamiltonian $``$ with $`M`$ holes can be generally written as

$$|\psi =\underset{𝐑_1\mathrm{}𝐑_M}{}\mathrm{\Psi }(𝐑_1\mathrm{}𝐑_M)|𝐑_1\mathrm{}𝐑_M,$$
(2)
where $`|𝐑_1\mathrm{}𝐑_M`$ are normalized states representing suitable superpositions of spin states with the constraint of having the holes placed in the sites labeled by $`(𝐑_1\mathrm{}𝐑_M)`$. The sum is over all distinct hole configurations. The matrix elements $`K_{ij}`$ of the hopping operator $`\widehat{K}`$ on the reduced Hilbert space defined by the hole configurations are non vanishing only if a single hole hops at nearest neighbor distance. Therefore, $`K_{ij}`$ has no diagonal elements. For instance, the element corresponding to the hopping $`𝐑_1𝐑_{}^{}{}_{1}{}^{}`$ is $`𝐑_{}^{}{}_{1}{}^{}\mathrm{}𝐑_M|\widehat{K}|𝐑_1\mathrm{}𝐑_M`$ and its modulus is always smaller or equal than $`t`$. By construction, the lowest eigenvalue $`E_0`$ of this matrix is the GS energy of the infinite-$`U`$ Hubbard model.

Now we prove that a lower bound of the GS energy is obtained by considering the matrix with all the non vanishing entries equal to $`t`$. In fact, for the hermitian matrix with only non diagonal elements $`K_{ij}`$, GS energy $`E_0<0`$ and with components $`\psi _i=\mathrm{\Psi }(𝐑_1\mathrm{}𝐑_M)`$ we have:

$`E_0`$ $`=`$ $`|{\displaystyle \underset{i,j}{}}\psi _i^{}K_{ij}\psi _j|{\displaystyle \underset{i,j}{}}|K_{ij}||\psi _i||\psi _j|`$ (3)
$``$ $`{\displaystyle \underset{i,j}{}}B_{ij}|\psi _i||\psi _j|E_B`$ (4)
Here, the matrix $`B`$ describes a gas of $`M`$ hard-core bosons (HCB) with nearest neighbor hopping $`B_{ij}=t`$ for the non zero entries, which satisfies $`|B_{ij}||K_{ij}|`$. The last inequality in (4) follows from the Perron-Frobenius theorem  and $`E_B`$ is the lowest eigenvalue of $`B`$. This proves that a lower bound on the energy of the infinite-$`U`$ Hubbard model is given by the GS energy of an HCB gas with density equal to the doping of the infinite-$`U`$ Hubbard model: $`E_BE_0E_{FM}`$. These bounds clearly hold in any dimension $`d`$.

It is worth noting that the lower bound is not at all trivial, indeed i) the Hilbert space for a system of $`M`$ HCB is much smaller than the Hilbert space of a system of spin-1/2 with $`M`$ holes, ii) in general it is not true that, for a given Hamiltonian $``$, the energy of the bosonic GS is lower than the energy of the fermionic GS. This is actually the case for the spin-1/2 infinite-$`U`$ Hubbard model but it does not hold in general (it is true if all the off-diagonal matrix elements of $``$ for the bosonic model are negative). For example let us consider the two-dimensional $`tJ`$ model

$`_{tJ}=t{\displaystyle \underset{i,j,\sigma }{}}𝒫_Gc_{i,\sigma }^{}c_{j,\sigma }𝒫_G`$ (5)
$`+J{\displaystyle \underset{i,j}{}}\left(S_iS_j{\displaystyle \frac{1}{4}}n_in_j\right),`$ (6)
where $`S_i=\frac{1}{2}_{\sigma ,\sigma ^{}}c_{i,\sigma }^{}\tau _{\sigma ,\sigma ^{}}c_{i,\sigma ^{}}`$, being $`\tau _{\sigma ,\sigma ^{}}`$ the Pauli matrices, and $`n_i=_\sigma c_{i,\sigma }^{}c_{i,\sigma }`$. If we consider $`L=16`$ and $`2`$ holes, the fermionic GS has an energy lower than the corresponding bosonic GS  for $`J0.2t`$. Indeed for this Hamiltonian, even for bosons, the off-diagonal terms are not all negative definite. In Table I we report the GS energies of the $`tJ`$ model for different value of $`J`$ both for fermions and bosons. Remarkably Lanczos diagonalizations on finite clusters show that, even at finite $`J`$, the hole-hole correlations in the $`tJ`$ model are quite similar to the corresponding HCB results .

Let us consider now the relevance of these bounds on the GS energy in the two dimensional case of the infinite-$`U`$ Hubbard model. The GS energy of an HCB gas can be computed numerically by quantum Monte Carlo on finite systems . Moreover the T-matrix formalism in the low density limit allows to obtain the following analytic expression :

$$\frac{E_B(\delta )}{L}4t\delta 2\pi t\frac{\delta ^2}{\mathrm{log}\delta },$$
(7)
with $`\delta =M/L`$.

If PS occurs between a hole-free region (with $`\delta =0`$) and a hole-rich phase (with $`\delta =\delta _c`$) then for finite size systems the function

$$ϵ(\delta )=\frac{E_0(\delta )E_0(0)}{Mt}$$
(8)
displays a minimum at $`\delta =\delta _c`$  while, in the thermodynamic limit, $`ϵ(\delta )`$ would be constant in the range $`0\delta \delta _c`$.

In Fig. (1) we show the previously discussed bounds on $`ϵ(\delta )`$. Since $`E_B(0)=E_{FM}(0)=E_0(0)=0`$, $`ϵ(\delta )`$ for the infinite-$`U`$ Hubbard model lies in between $`ϵ_B(\delta )`$ and $`ϵ_{FM}(\delta )`$. For any $`finite`$ density the curvature of both curves is positive ruling out the possibility to have PS starting from $`\delta =0`$. Notice that a previous lower bound to the GS energy given by Trugman  just corresponds to $`ϵ(\delta )=4`$, and then it would not allow to draw any conclusion on the occurrence of PS in this model.

In conclusion we have obtained a rigorous result that rules out PS in the infinite-$`U`$ Hubbard model starting from the zero doping limit. Of course our findings do not exclude that PS could take place between two finite densities or at finite $`U>0`$.

This work has been partially supported by INFM and MURST (COFIN99).


text:
# Calculation of Ground State Energy for Confined Fermion Fields

## 1 Introduction

Nowadays a number of methods is developed for evaluation of the Casimir energy for systems with quantized fields under nontrivial boundary conditions . The most known of them are the Green function method, the zeta function approach, the application of contour integration, the multiple scattering method, and the direct mode summation with thermal regularization  (for recent discussion see ). An important field of application for these techniques is an investigation of ground state energy for the quark bag models with account of one-loop corrections due to the filled Dirac’s negative energy sea . In this paper we propose an approach for explicit calculation of ground state energy for the models with confined fermion field coupled to the boson (scalar) field in a certain spatial region, so called hybrid bag models . This method allows to extract the singular terms from the divergent sums of eigenvalues $`\omega _n`$ and obtain the finite Casimir energy of fermionic sea as the result of absorbing singularities into the “dressed” model parameters. The realization of a such strategy requires, as a rule, an including of contact terms into the “classical” expression for the energy of system. The form of these terms is determined by their dependence on the geometric parameters . By means of this method, we renormalize the ground state energy of the hybrid chiral bag model in $`(1+1)D`$ and study it as a function of the bag’s size. It’s shown that this function has unique minimum what means that there exists a stable configuration for the considered model.

## 2 Regularization

The ground state energy of fermion field is defined as the vacuum expectation value of the Hamiltonian:

$$0|H|0E_0=\frac{1}{2}\text{ }\left(\underset{\omega _n>0}{}\omega _n\underset{\omega _k<0}{}\omega _k\right).$$
(1)
Provided that the spectrum $`\omega _n`$ possesses the symmetry $`\omega _n\omega _n`$, the vacuum energy is

$$E_0=\underset{\omega _n>0}{}\omega _n.$$
(2)
This sum diverges. On of the ways to regularize it is to use the exponential cutoff  what yields

$$E_0^{exp}(\tau )=\underset{\tau 0}{lim}\left(\underset{\omega _n>0}{}\omega _n\text{e}^{\tau \omega _n}\right)=\underset{\tau 0}{lim}\frac{d}{d\tau }\underset{\omega _n>0}{}\text{e}^{\tau \omega _n},$$
(3)
where $`\tau =\stackrel{~}{\tau }\mu ^1`$ and $`\mu `$ has the dimension of mass. This regularization can be applied for calculation of the Casimir energy in several simple situations , but usually only for systems with an explicitly known spectrum $`\omega _n`$. In this paper we’ll try to show that one can succeed in calculations even if this is not the case. Suppose that an asymptotical expansion of $`\omega _n`$ for $`n1`$  can be found for a model with an unknown spectrum. This expansion reads

$$\omega _n=\underset{i=1}{\overset{\mathrm{}}{}}\mathrm{\Omega }_in^i=\mathrm{\Omega }_1n+\mathrm{\Omega }_0+\frac{\mathrm{\Omega }_1}{n}+O\left(\frac{1}{n^2}\right).$$
(4)
These are the three leading terms in (4) which determine the divergencies in (2). To begin, consider the case when the expansion (4) contains only two terms (what really takes place, e. g. for the free massless fermion field in (1+1)-bag):

$$\omega _n^{(1)}=\mathrm{\Omega }_1n+\mathrm{\Omega }_0.$$
(5)
Then the regularized energy is

$$E_1^{exp}(\tau )=\underset{\tau 0}{lim}\frac{d}{d\tau }\underset{n=1}{}\text{e}^{\tau (\mathrm{\Omega }_1n+\mathrm{\Omega }_0)}\omega _0=\underset{\tau 0}{lim}\frac{d}{d\tau }\left(\text{e}^{\tau \mathrm{\Omega }_0}\sigma _1(\tau )\right)\omega _0,$$
(6)
$$\sigma _1(\tau )=\underset{n=1}{}\text{e}^{\tau \mathrm{\Omega }_1n}=\frac{1}{\text{e}^{\tau \mathrm{\Omega }_1}1}.$$
(7)
The term $`\omega _0`$ in (6) is written separately because the expansion (4) is not valid for it. By virtue of

$$\frac{1}{\text{e}^x1}=\underset{k=0}{}\frac{B_k}{k!}x^{k1},$$
(8)
where $`B_k`$ are the Bernoulli numbers, one gets

$$E_1^{exp}(\tau )=\frac{1}{\tau ^2\mathrm{\Omega }_1}+\frac{\mathrm{\Omega }_1}{12}+\frac{\mathrm{\Omega }_0}{2}+\frac{\mathrm{\Omega }_0^2}{2\mathrm{\Omega }_1}\omega _0.$$
(9)
This expression contains the quadratic divergence $`E^{quad}(\tau )=\frac{1}{\tau ^2\mathrm{\Omega }_1}=\frac{\mu ^2}{\stackrel{~}{\tau }^2\mathrm{\Omega }_1}`$ depending on the arbitrary mass $`\mu `$ and the geometric parameters of the bag from $`\mathrm{\Omega }_1`$.

Consider now the case when

$$\omega _n^{(2)}=\mathrm{\Omega }_1n+\mathrm{\Omega }_0+\frac{\mathrm{\Omega }_1}{n}.$$
(10)
In this situation the regularized energy is

$$E_2^{exp}(\tau )=\underset{\tau 0}{lim}\frac{d}{d\tau }\sigma _2(\tau )\omega _0,$$
(11)
$$\sigma _2(\tau )=\underset{n=1}{}\text{e}^{\tau \left(\mathrm{\Omega }_1n+\mathrm{\Omega }_0+\frac{\mathrm{\Omega }_1}{n}\right)}=\text{e}^{\tau \mathrm{\Omega }_0}\underset{n=1}{}\text{e}^{\tau \mathrm{\Omega }_1n}\text{e}^{\tau \mathrm{\Omega }_1/n}.$$
(12)
One easily see that only the two leading terms in the expansion of $`\text{e}^{\tau \mathrm{\Omega }_1/n}`$ in powers of $`\tau `$ yield a non-vanishing contribution to $`E_2^{exp}(\tau )`$ as $`\tau 0`$. Then we get

$$\sigma _2(\tau )=\text{e}^{\tau \mathrm{\Omega }_0}\underset{n=1}{}\text{e}^{\tau \mathrm{\Omega }_1n}\left(1\frac{\tau \mathrm{\Omega }_1}{n}\right)+O(\tau ^2)=$$

$$=\sigma _1(\tau )+\tau \mathrm{\Omega }_1\text{e}^{\tau \mathrm{\Omega }_0}\underset{n=1}{}\frac{1}{n}\text{e}^{\tau \mathrm{\Omega }_1n}+O(\tau ^2),$$
(13)
where $`\sigma _1(\tau )`$ is already obtained in (7). Using the known relation

$$\underset{n=1}{}\frac{1}{n}\text{e}^{\alpha n}=\mathrm{ln}(1\text{e}^\alpha )$$

we find

$$\sigma _2(\tau )=\sigma _1(\tau )+\tau \mathrm{\Omega }_1\text{e}^{\tau \mathrm{\Omega }_0}\left(\mathrm{ln}\tau \mathrm{\Omega }_1\frac{\tau \mathrm{\Omega }_1}{2}\right)+O(\tau ^2).$$
(14)
So the regularized energy reads

$$E_2^{exp}(\tau )=\frac{1}{\tau ^2\mathrm{\Omega }_1}+\mathrm{\Omega }_1\mathrm{ln}\tau \mu +\mathrm{\Omega }_1\mathrm{ln}\frac{\mathrm{\Omega }_1}{\mu }+\frac{\mathrm{\Omega }_1}{12}+\frac{\mathrm{\Omega }_0}{2}+\frac{\mathrm{\Omega }_0^2}{2\mathrm{\Omega }_1}+\mathrm{\Omega }_1\omega _0.$$
(15)
The contribution of the terms of order $`O(1/n^2)`$ ($`E^{fin}`$) is finite and can be found in any particular case. The divergent parts

$$E^{quad}(\tau )=\frac{1}{\tau ^2\mathrm{\Omega }_1}\text{and}E^{log}(\tau )=\mathrm{\Omega }_1\mathrm{ln}\tau \mu $$

is to be compensated using a certain renormalization scheme, which will be discussed later.

Let us compare the result obtained above with another method of regularization — the $`\zeta `$-function approach. In this case the energy is regularized as:

$$E_1^{zeta}(s)=\underset{n}{}\omega _n\underset{s1}{lim}\mu ^{1+s}\underset{n}{}\omega _n^s,$$
(16)
where the arbitrary mass $`\mu `$ is included in order to restore the correct dimension. We assume that this mass equals $`\mu `$ in (16). For $`\omega _n^{(1)}=\mathrm{\Omega }_1n+\mathrm{\Omega }_0`$ one has

$$E_1^{zeta}(s)=\underset{s1}{lim}\mu ^{s+1}\underset{n=1}{}(\mathrm{\Omega }_1n+\mathrm{\Omega }_0)^s\omega _0=\underset{s1}{lim}z_1(s)\omega _0,$$
(17)
where

$$z_1(s)=\mu ^{s+1}\underset{1}{}(\mathrm{\Omega }_1n+\mathrm{\Omega }_0)^s=$$

$$=\underset{1}{}\frac{1}{(\mathrm{\Omega }_1n)^s}\left(1\frac{s\mathrm{\Omega }_0}{\mathrm{\Omega }_1n}+\frac{s(s+1)}{2!}\left(\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }_1n}\right)^2+O(1/n^3)\right).$$
(18)
It is clear that the contribution of terms $`O(1/n^3)`$ vanishes as $`s1`$, so the regularized sum $`E_1^{zeta}(s)`$ is determined by the leading three terms in the expansion (18). Then one finds:

$$z_1(s)=\mu ^{s+1}\left(\mathrm{\Omega }_1^s\zeta (s)\frac{s\mathrm{\Omega }_0}{\mathrm{\Omega }_1^{s+1}}\zeta (s+1)+\frac{\mathrm{\Omega }_0^2}{2\mathrm{\Omega }_1^{s+2}}s(s+1)\zeta (s+2)\right).$$
(19)
The values of $`\zeta (z)`$ analytically continued to the total real axis are known :

$$\zeta (0)=\frac{1}{2},\zeta (1)=\frac{B_2}{2}=\frac{1}{12}.$$
(20)
For $`\zeta (z)`$ in the vicinity of $`1`$ we have

$$\underset{z1}{lim}\zeta (z)=\frac{1}{z1}+C,$$
(21)
where $`C=0.5772156649\mathrm{}`$ is the Euler constant. Thus taking the limit $`s1`$ one obtains

$$E_1^{zeta}(s)=\frac{\mathrm{\Omega }_1}{12}+\frac{\mathrm{\Omega }_0}{2}+\frac{\mathrm{\Omega }_0^2}{2\mathrm{\Omega }_1}\omega _0.$$
(22)
This expression reproduces explicitly the finite part of (9). The absence of the divergent term is due to the analytical continuation for $`\zeta (z)`$. Taking into account the next term of the expansion in powers of $`1/n`$ (10) we find

$$E_2^{zeta}(s)=$$

$$=\underset{s1}{lim}\mu ^{s+1}\underset{n=1}{}\left(\mathrm{\Omega }_1n+\mathrm{\Omega }_0+\frac{\mathrm{\Omega }_1}{n}\right)^s\omega _0=\underset{s1}{lim}z_2(s)\omega _0,$$
(23)
where

$$z_2(s)=z_1(s)\frac{s\mathrm{\Omega }_1\mu ^{s+1}}{\mathrm{\Omega }_1^{s+1}}\zeta (s+2)$$
(24)
(we write down only the terms which yield a non-vanishing contribution). The second term in (24) can be found using (21) and the expansion $`x^\epsilon =1+\epsilon \mathrm{ln}x+O(\epsilon ^2)`$, $`\epsilon =s+1`$, $`\epsilon 0`$. Then

$$\underset{s1}{lim}s\mathrm{\Omega }_1\left(\frac{\mu }{\mathrm{\Omega }_1}\right)^{s+1}\zeta (s+2)=$$

$$=\mathrm{\Omega }_1\underset{\epsilon 0}{lim}(1+\epsilon )\left(1\epsilon \mathrm{ln}\frac{\mathrm{\Omega }_1}{\mu }\right)\left(\frac{1}{\epsilon }+C\right)=\mathrm{\Omega }_1\left(\frac{1}{\epsilon }1+C\mathrm{ln}\frac{\mathrm{\Omega }_1}{\mu }\right).$$
(25)
Therefore the regularized energy is of the form

$$E_2^{zeta}(s)=\frac{\mathrm{\Omega }_1}{12}+\frac{\mathrm{\Omega }_0}{2}+\frac{\mathrm{\Omega }_0^2}{2\mathrm{\Omega }_1}\omega _0+\mathrm{\Omega }_1\left(\mathrm{ln}\frac{\mathrm{\Omega }_1}{\mu }+1\right)\mathrm{\Omega }_1\left(\frac{1}{\epsilon }+C\right),$$
(26)
what coincides with (15) (except of the quadratic singularity $`E_2^{quad}`$) provided that

$$\frac{1}{\tau }=\mu \gamma \text{e}^{1/\epsilon },\text{ln}\gamma =C.$$
(27)
If one needs to take into account the contribution of terms of order $`O\left(\frac{1}{n^2}\right)`$, then their sum $`(E^{fin})`$ should be added to (15) and (26). Therefore it is shown that the exponential and $`\zeta `$-function regularizations provide the equivalent results for regularized Casimir energy in $`(1+1)D`$.

The singular part is:

$$E^{div}(\tau )=\frac{1}{\tau ^2\mathrm{\Omega }_1}+\mathrm{\Omega }_1\mathrm{ln}\tau \mu .$$
(28)
It can be removed, for example, by means of the redefinition of the “bare” constants from the initial Lagrangian. It is very important to note that the way of renormalization is prescribed by the dependence of the extracted singularities on the geometric parameters (in our case the only geometric parameter is the bag’s size $`R`$). We collect all singularities with the similar dependence from $`R`$, and then find the “contact term” to be redefined .

Let us describe this procedure in the simplest case, i. e., MIT bag model with massive fermions in $`(1+1)D`$ . The Lagrangian reads

$$_{MIT}=\theta (|x|<R)(i\overline{\psi }\widehat{}\psi m_F\overline{\psi }\psi B)+\theta (|x|>R)(i\overline{\psi }\widehat{}\psi M_F\overline{\psi }\psi ),$$
(29)
where $`B`$ is the so-called bag constant which characterizes an excess of the energy density inside a hadron compared to the energy of the nonperturbative vacuum. Taking the limit $`M_F\mathrm{}`$ in the exterior region $`|x|>R`$ one gets the “bag”, which is just the segment of real axis $`[R,R]`$. The boundary conditions

$$(\pm i\gamma ^1+1)\psi (\pm R)=0$$
(30)
lead to the spectrum

$$\omega _n=\sqrt{\left(\frac{\pi }{2R}n+\frac{\pi }{4R}\right)^2+m_F^2}.$$
(31)
Assume that the fermion mass $`m_F`$ is small and expand the energy in powers of it. Then (with the accuracy up to $`m_F^4`$) one has

$$\omega _n=(\mathrm{\Omega }_1n+\mathrm{\Omega }_0)+\frac{m_F^2}{2(\mathrm{\Omega }_1n+\mathrm{\Omega }_0)}+O(m_F^4),$$
(32)
$`\mathrm{\Omega }_1=\frac{\pi }{2R}`$, $`\mathrm{\Omega }_0=\frac{\pi }{4R}`$. For $`n1`$ the expansion (4) can be found, where

$$\mathrm{\Omega }_1=\frac{m_F^2}{2\mathrm{\Omega }_1}=\frac{m_F^2R}{\pi }$$
(33)
and for $`n=0`$ we have

$$\omega _0=\mathrm{\Omega }_0+\frac{m_F^2}{2\mathrm{\Omega }_0}=\mathrm{\Omega }_0+2\mathrm{\Omega }_1.$$
(34)
The divergent part of the energy ( according to (28)) is

$$E^{div}(\tau |R)=\left(\frac{1}{\tau ^2\pi }+\frac{m_F^2}{2\pi }\mathrm{ln}\tau \mu \right)2R=2B^{}(\tau )R.$$
(35)
In accordance with the approach of , the renormalization is performed by means of the redefinition of the bare bag constant $`B`$ in the Lagrangian (29):

$$B=B_0B^{}(\tau ).$$
(36)
The remaining parts of energy (15) are finite and can be found explicitly. The contribution of terms $`O\left(n^2\right)`$ in the expansion (4) is determined by

$$E^{fin}(R)=\frac{m_F^2}{2}\underset{n=1}{}\left(\frac{1}{\mathrm{\Omega }_1n+\mathrm{\Omega }_0}\frac{1}{\mathrm{\Omega }_1n}\right)=$$

$$=\frac{m_F^2}{2\mathrm{\Omega }_1}\underset{n=1}{}\frac{1}{n(2n+1)}=\frac{2m_F^2R}{\pi }(1\mathrm{ln}2).$$
(37)
Therefore, the renormalized energy of the fermionic sea for $`(1+1)D`$ massive MIT bag model as the function of the bag’s size $`R`$ and renormalized bag constant $`B_0`$ reads (up to the terms of order $`m_F^4`$):

$$E_{MIT}(R)=2B_0R\frac{\pi }{48R}+\frac{m_F^2R}{\pi }\left(12\mathrm{ln}2+\mathrm{ln}\frac{\pi }{2R\mu }\right)+O(m_F^4).$$
(38)
Taking the limit $`m_F0`$ (massless MIT bag model), one gets from (38) the well-known result of . Note, that this configuration is unstable and tends to $`R0`$. It is possible to make it stable by adding one valence fermion into the lowest level $`n=0`$. Then the energy reads

$$\stackrel{~}{E}_{MIT}(R)=2B_0R+\frac{11}{48}\frac{\pi }{R}+\frac{m_F^2R}{\pi }\left(32\mathrm{ln}2+\mathrm{ln}\frac{\pi }{2R\mu }\right)+O(m_F^4).$$
(39)
## 3 Two-phase Hybrid Bag Model

Now we are ready to study the hybrid bag model, in which the fermion (“quark”) field interacts with scalar (“meson”) field. Compared to the models considered in , here is no phase of the massless fermions. The Lagrangian reads

$$=i\overline{\psi }\widehat{}\psi +\frac{1}{2}\text{ }_\mu \phi ^\mu \phi \theta (|x|<R)\left(\frac{1}{2}\text{ }m_F[\overline{\psi },\text{e}^{i\gamma _5\phi }\psi ]B\right)$$

$$\theta (|x|>R)(V(\phi )+\frac{1}{2}\text{ }M_F[\overline{\psi },\text{e}^{i\gamma _5\phi }\psi ]).$$
(40)
The commutator in the fermionic forms provides the charge-conjugation symmetry. For the sake of simplicity, we consider the scalar field $`\phi (x,t)`$ in the mean-field approximation (MFA) , i. e., assume it to be a $`c`$number function of the space-time variables; suppose also that it is independent of the temporal coordinate: $`\phi (x,t)=\phi (x)`$. Besides this, we will take into account only the leading terms of the perturbative expansion in the chiral coupling constant . In other words, the fermionic mass $`m_F`$ is assumed to be small inside the segment $`|x|<R`$ while the mass $`M_F`$ in the exterior region $`|x|>R`$ is infinitely large. Thus the fermion field at $`|x|>R`$ vanishes and the equations of motion for $`|x|<R`$ are

$$(i\widehat{}m_F\text{e}^{i\gamma _5\phi })\psi (x,t)=0,$$
(41)
$$\phi ^{^{\prime \prime }}=i\frac{m_F}{2}[\overline{\psi },\gamma _5\text{e}^{i\gamma _5\phi }\psi ]_{sea},$$
(42)
where in r. h. s. of (42) the v. e. v. of the axial current is taken accordingly to MFA. For $`|x|>R`$, the scalar field is determined by the non-linear equation

$$\phi ^{\prime \prime }(x)=V_\phi ^{}(\phi ).$$
(43)
The solution of the equations (41, 42) have been firstly found by Sveshnikov and Silaev and studied in detail in , so we will discuss it here only briefly. Let us assume that the scalar field is an odd soliton function and for $`|x|>R`$ it has the form:

$$\phi (x)=\pi \left(1A\text{e}^{mx}\right),x>0$$
(44)
where $`m`$ is the meson mass, and $`\phi (x)`$ for $`x<0`$ is determined by the oddness. The boundary conditions are

$$(\pm i\gamma ^1+\text{e}^{i\gamma _5\phi (x)})\psi (\pm R)=0,$$

$$\phi (\pm R\pm 0)=\phi (\pm R0),$$

$$\phi ^{}(\pm R\pm 0)=\phi ^{}(\pm R0).$$
(45)
The fermionic spectrum possesses the symmetry $`\nu \nu `$, where $`\nu ^2=k^2+m_F^2`$. The corresponding unitary transformations of the wave function are $`\chi i\gamma _1\chi `$ (here $`\chi =\text{e}^{i\gamma _5\phi /2}\psi `$), therefore the v. e. v. of the axial current in the r. h. s. of (42) is equal to zero and the solution of (42) for $`|x|<R`$ appears to be the linear function $`\phi (x)=2\lambda x`$ (for detailed discussion of the self-consistent solution for this system see ). The eigenvalues $`\omega _n`$ can be obtained from the equation

$$\left(1\text{e}^{2ikR}\frac{m_F+i(\nu +k)}{m_F+i(\nu k)}\right)\left(1\text{e}^{2ikR}\left(\frac{\nu k}{\nu +k}\right)\frac{m_Fi(\nu +k)}{m_Fi(\nu k)}\right)=$$

$$=\left(1\text{e}^{2ikR}\frac{m_Fi(\nu +k)}{m_Fi(\nu k)}\right)\left(1\text{e}^{2ikR}\left(\frac{\nu k}{\nu +k}\right)\frac{m_F+i(\nu +k)}{m_F+i(\nu k)}\right),$$
(46)
where $`\nu =\omega \lambda `$. Provided that the signatures of $`\omega _n`$ and $`\nu _n`$ are the same for all $`n`$, the fermionic Casimir energy is determined by (2) with the replacement $`\omega _n\nu _n`$. For the sake of being definite we regularize this expression by means of the exponential cutoff (3):

$$E_0^{exp}(\tau )=\underset{\tau 0}{lim}\left(\underset{\nu _n>0}{}\nu _n\text{e}^{\tau \nu _n}\right)=\underset{\tau 0}{lim}\frac{d}{d\tau }\underset{\nu _n>0}{}\text{e}^{\tau \nu _n}.$$

The equation (46) can be written as

$$m_F\mathrm{sin}2Rk+k\mathrm{cos}2Rk=0.$$
(47)
Using the expansion in powers of $`m_F`$ for $`k`$

$$k=\stackrel{~}{k}_0+m_F\stackrel{~}{k}_1+m_F^2\stackrel{~}{k}_2+O(m_F^3)$$
(48)
one obtains

$$\nu =\stackrel{~}{k}_0+m_F\stackrel{~}{k}_1+m_F^2\left(\stackrel{~}{k}_2+\frac{1}{2\stackrel{~}{k}_0}\right)+O(m_F^3).$$
(49)
One can solve (47) in any order, what yields (up to $`O\left(m^3\right)`$),

$$\nu _n=\frac{\pi }{2R}n+\frac{\pi }{4R}+\frac{2}{\pi }m_F(1+m_FR)\frac{1}{(2n+1)}\frac{16Rm_F^2}{\pi ^3(2n+1)^3}+O(m_F^3).$$
(50)
Expanding $`\nu _n`$ in powers of $`1/n`$ one gets

$$\nu _n=\mathrm{\Omega }_1n+\mathrm{\Omega }_0+\frac{\mathrm{\Omega }_1}{n}+O\left(\frac{1}{n^2}\right),$$
(51)
where

$$\mathrm{\Omega }_1=\frac{\pi }{2R},\mathrm{\Omega }_0=\frac{\pi }{4R},\mathrm{\Omega }_1=\frac{m_F}{\pi }(1+m_FR),$$
(52)
what differs from the massive MIT bag model by the term $`m_F/\pi `$ in $`\mathrm{\Omega }_1`$ which is independent of the bag’s size $`R`$. For $`n=0`$ we have

$$\nu _0=\mathrm{\Omega }_0+2\mathrm{\Omega }_1\frac{16Rm_F^2}{\pi ^3}.$$
(53)
Now one can use the formulae for the regularized energy (15). The finite fermionic energy should be of the form

$$E_F=E_{Cas}+2BR+\mathrm{\Lambda },$$
(54)
where $`B`$ and $`\mathrm{\Lambda }`$ are the bag and cosmological constants, and $`E_{Cas}`$ is the divergent fermionic Casimir energy. Hence the renormalization requires the redefinition of $`B`$ (36) and $`\mathrm{\Lambda }`$:

$$\mathrm{\Lambda }=\mathrm{\Lambda }_0\frac{m_F}{\pi }\mathrm{ln}\tau \mu .$$
(55)
The contribution $`E^{fin}`$ of the non-singular terms $`O\left(\frac{1}{n^2}\right)`$ reads

$$E^{fin}=\mathrm{\Omega }_1\underset{n=1}{}\frac{1}{n(2n+1)}+\frac{16Rm_F^2}{\pi ^3}\underset{n=0}{}\frac{1}{(2n+1)^3}=$$

$$=\mathrm{\Omega }_12(1\mathrm{ln}2)+\frac{16Rm_F^2}{\pi ^3}A,$$
(56)
where $`A=1.051799\mathrm{}`$. Therefore the finite energy of the fermionic sea is (provided that $`\mathrm{\Lambda }_0=0`$):

$$E(R)=2B_0R\frac{\pi }{48R}+$$

$$+\frac{m_F}{\pi }(1+m_FR)\left[\mathrm{ln}\frac{\pi }{2\mu R}+12\mathrm{ln}2\right]+\frac{16Rm_F^2}{\pi ^3}A.$$
(57)
A total energy of this system contains also the contribution of the scalar field. In order to find it we use the boundary condition for $`\phi `$ and it’s derivative (45) what yields

$$2\lambda =\frac{\pi m}{mR+1}.$$
(58)
By virtue of the virial theorem in the external region $`|x|>R`$ one obtain $`\frac{1}{2}\text{ }\phi ^2(x)=V(\phi )`$, so the scalar field energy reads

$$E_\phi (R)=\frac{1}{2}\text{ }_R^R𝑑x\phi ^2(x)+(_{\mathrm{}}^R+_{\mathrm{}}^R)dx\phi ^2(x)=\frac{\pi ^2m}{mR+1}.$$
(59)
Note that the representation of the scalar field in the form (44) is valid only at the distances much larger than the soliton size, i. e. of order $`m^1`$. This restriction allows us to postulate (in the framework of our model) the following relation between the meson mass $`m`$ and the bag’s radius:

$$mR=K_\pi ,K_\pi 1.$$

Thus the contribution of the scalar field to the total energy of the bag is

$$E_\phi (R)=\frac{\pi ^2K_\pi }{(1+K_\pi )R}$$
(60)
and the total energy consists of the fermionic and the bosonic parts:

$$E_{tot}(R)=E_F(R)+E_\phi (R).$$
(61)
In terms of the dimensionless variables the total energy $`=E/m_F`$ reads $`(x=m_FR)`$:

$$(x)=2B_1x\frac{\pi }{48x}+$$

$$+\frac{1}{\pi }(1+x)\left[\mathrm{ln}\frac{\pi m_1}{2x}+12\mathrm{ln}2\right]+\frac{16x}{\pi ^3}A+\frac{\pi ^2K_\pi }{(1+K_\pi )x}.$$
(62)
It is easy to see that this energy has the unique minimum and the configuration is stable.

## 4 Conclusion

The present approach to renormalization of the infinite ground state energy of the models with confined fermion fields in $`(1+1)D`$ is based on the analytical regularization of the divergent sums  (e. g., the exponential cutoff, or $`\zeta `$-function regularization), the extraction of the singular terms in the form of poles at the point of the physical value of regularization parameter, and subsequent absorbing of these singularities by means of redefinition of bare constants from the initial Lagrangian. In this framework the singularities are isolated unambiguously using the scheme analogous to MS in QFT. The remaining parts are finite and demonstrate the non-trivial dependence on the geometrical parameters of the model. The generalization of the proposed framework to the higher dimensions appears to be a separate problem and will be the subject of the future investigations .

## 5 Acknowledgments

Author would like to express his gratitude to Prof. K. A. Sveshnikov for fruitful discussions and critique and Dr. O. V. Pavlovsky for discussions and many useful remarks. Author thanks the Abdus Salam ICTP in Trieste for warm hospitality during 2 month in 1999 when a part of this work was done. This work is partially supported by RFBR (00-15-96577) and Young scientists Fellowship of Nuclear Physics Institute, Moscow.


text:
# Structure and dynamics of a model glass: influence of long-range forces

## I Introduction

The connection between the structure and the dynamics of a given system is of particular interest in glasses since the structural entities at the origin of some specific vibrational features like the Boson Peak (BP) are not yet well identified. For instance in silica, the archetype of network forming glasses, the existence of the BP has been connected alternatively to the presence of SiO4 tetrahedra , ”domains”  or voids . In molecular dynamics (MD) simulations the most important ingredient is the force field used to mimic the interparticle interactions: the aim is to reproduce as closely as possible the experimental vibrational data and then, by looking at the atomic positions, to give an interpretation in terms of structure. In that sense the so-called ”BKS” potential  has been widely used to describe silica and the results obtained on the structure , the dynamical  and the thermal properties  compare well with experiments. This pairwise classical potential is the combination of a short-range term and the long-range Coulomb interaction (CI) which is obviously the most important ingredient of the force field . One remarkable property of this potential is that CI forces alone are sufficient to reproduce the tetrahedral covalent bonding network of the silica atoms without the addition of extra terms such as angular terms. Our idea in this paper is to vary the amplitude of the CI (this is not strictly equivalent to varying the range) in order to study the variation of the structure and the dynamics of what can be seen as a toy silica glass. With this procedure we can study the properties of a system going from a strong network glass (the silica glass described by the usual BKS potential) to a fragile two-component soft sphere glass when the CI is turned off. The question we would like to address is: how do the structure and the dynamics change along this path ? The aim is to study the connection between structural features like the First Sharp Diffraction Peak (FSDP) or the existence of a tetrahedral network and the BP, this connection being the topic of contradictory interpretations in the literature .

## II Simulations

We performed molecular dynamics simulations for microcanonical systems containing 216 silicon and 432 oxygen atoms confined in a cubic box of edge length $`L=21.48`$ Å, which corresponds to a mass density of $`2.18`$g/cm3 very close to the experimental value of $`2.2`$g/cm3. Periodic boundary conditions were used to limit surface effects. The particles interact via the BKS potential whose standard functional form is, for two particles $`i`$ and $`j`$ (which can be either Si or O):

$$U(r_{ij})=\frac{q_iq_je^2}{r_{ij}}+A_{ij}\mathrm{exp}(B_{ij}r_{ij})\frac{C_{ij}}{r_{ij}^6}$$
(1)
where $`r_{ij}`$ is the interparticle distance, $`e`$ the charge of an electron and the parameters $`A_{ij}`$, $`B_{ij}`$ and $`C_{ij}`$ are fixed as follows: $`A_{SiO}=18003.7572`$ and $`A_{OO}=1388.773`$ eV;$`B_{SiO}=4.87318`$ and $`B_{OO}=2.76`$ Å-1; $`C_{SiO}=133.5381`$ and $`C_{OO}=175.0`$ eVÅ6. The partial charges used to compute the standard CI are given by $`q_{Si}=2.4`$ and $`q_O=1.2`$. This original form contains an unphysical divergence at very short distances which we have corrected by adding a short range repulsive term ($`1/r_{ij}^{40}`$).
In order to vary the strength of the CI (which is calculated using the Ewald summation method) we simply multiply the partial charges by a parameter $`\lambda `$ which is fixed at the beginning of each simulation:

$$q_i\lambda q_i$$
(2)
A simulation starts from a liquid configuration at 7000K which is equilibrated during 50000 time steps (35 ps). This system is then cooled to zero temperature at a quench rate of $`2.3\times 10^{14}`$ K/s which is obtained by removing the corresponding amount of energy from the total energy of the system at each iteration. The aim here is not to obtain fine details of the structure of real silica which explains the use of a rather fast quenching rate compared to the rates that can nowadays be achieved. In any case computer quenches are still far away from those used in experiments. The quenched glass samples are finally relaxed in the microcanonical ensemble during 100000 time steps (70 ps) after which the temperature has reached $``$ 2K. We then average the structural informations over the last 20000 time steps to avoid transient configurations and we diagonalize the dynamical matrix corresponding to the final configuration to obtain the vibrational density of states.
This procedure has been repeated for eleven values of $`\lambda `$ going from $`\lambda =0.1`$ to $`\lambda =1.1`$. The ideal ”$`\lambda =0`$” system would consist of a soft-sphere like arrangement for the oxygen atoms (similar to the systems we have studied earlier ) together with a random arrangement of the silicon atoms in the interstices of the oxygen structure. Because of the weak attractive O-O interaction present in the short-range part of the BKS potential, a small Coulomb interaction is necessary in order to stabilize this system: we found that a value of $`\lambda `$ around 0.065 is sufficient which explains why we started our investigation at $`\lambda =0.1`$. For this lowest $`\lambda `$ value we have checked on the radial pair distribution functions as well as on other structural characteristics that the structure was indeed very close to the ideal arrangement described above. On the other hand, strong values of $`\lambda `$ lead to an artificial local densification of the samples which are non homogeneous (holes start to form in the structure): for this reason the maximum value of $`\lambda `$ is limited to $`1.1`$. Of course the standard BKS results are recovered for $`\lambda =1`$.

## III Results and discussion

In order to examine the “quality” of the tetrahedral network we first study the tetrahedral O-Si-O angle $`\theta `$ which should be ideally equal to $`109.47^{}`$ in a perfect tetrahedron as well as the Si-O-Si angle $`\varphi `$ which measures the relative position and orientation of two neighboring SiO4 tetrahedra. To numerically determine the bond angle O-Si-O (resp. Si-O-Si), we determine for each Si (resp. O) atom the two nearest O (resp. Si) atoms and we calculate the angle between the two corresponding segments Si-O (resp. O-Si), the result being averaged over all the Si (resp. O) atoms. The variation of $`\theta `$ and $`\varphi `$ as a function of $`\lambda `$ is plotted in figure 1. For $`\lambda 0.5`$ it is obvious that the SiO4 tetrahedra are not formed and that the network is not connected: the system remains analogous to a soft sphere glass up to a threshold value of about 0.6. This can also be seen in the inset of figure 1 where the Si coordination number is plotted versus $`\lambda `$. This quantity has been calculated by the integration of the Si-O radial pair distribution functions up to a cut-off distance of 2.3 Å which corresponds to the first minimum for all the values of $`\lambda `$. This coordination number is greater than 5.5 for $`\lambda =0.1`$ and is close to 4 for $`\lambda 0.60.7`$. Nevertheless for $`\lambda 0.5`$ $`\theta `$, the tetrahedral angle, is almost constant and close to the angle in a perfect tetrahedron (indicated by the black arrow) while $`\varphi `$, the angle between the tetrahedra, is continuously increasing with $`\lambda `$. This indicates that the CI has a direct influence on the medium-range order and once the tetrahedra (the building blocks) are created, they form a connected network whose characteristics depend on the strength of the CI. In agreement with a previous study  the mean Si-O-Si angle is close to $`150^{}`$ for $`\lambda =1.0`$ (standard BKS description) slightly larger than the experimental value of $`\varphi `$ ($`144^{}`$  indicated by the hollow arrow). The calculated value closest to $`144^{}`$ corresponds to $`\lambda =0.7`$ but this value is not sufficient to recover other structural features of vitreous silica as we will see below, therefore this is not an argument to modify the standard BKS charges.
A second way to analyze the connectivity of the network is to monitor the distribution of $`n`$-fold rings: an $`n`$-fold ring is defined as the shortest closed path of alternating Si-O bonds. Therefore an $`n`$-fold ring consists of $`2n`$ alternating Si-O bonds. The ring distribution as a function of ring size is plotted in figure 2 for 3 values of $`\lambda `$: 0.2, 0.5 and 1. For $`\lambda =0.2`$ the ring distribution is typical of an open structure with a majority of small rings (the estimated number of large rings should be taken with a pinch of salt since finite size effects due to the limited length of the simulation box bias the results). For $`\lambda =0.5`$ the distribution becomes peaked for rings with a specific size (5 and 6) which extend over distances greater than the inter-atomic distances: molecules are formed in the system. Finally for $`\lambda =1.0`$ we find the usual distribution of rings, peaked at sixfold rings and almost symmetric around $`n=6`$ , typical of the silica network.
A complementary way of analyzing the structure is to compute the static structure factor S(q) as a function of $`\lambda `$. We used the standard way of calculating S(q)  and focused our attention on a common feature in glasses: the First Sharp Diffraction Peak. This pre-peak is related to medium-range structures in the glass and its origin has been the subject of vigorous debate over a number of years. In figure 3 we represent the amplitude of the FSDP as a function of $`\lambda `$ and in the inset we show the shape of S(q) for $`\lambda `$ = 0.2, 0.5 and 1 in the vicinity of the FSDP, together with the experimental curve . With decreasing values of $`\lambda `$ the amplitude of the FSDP decreases and one can say that for $`\lambda 0.4`$ this feature does not exist anymore in S(q). Indeed the residual value of the amplitude for small values of $`\lambda `$ is not significant and is an artifact of the Fourier transform as can be seen in the inset. Together with the results shown in figure 1, this indicates that the FSDP appears when the tetrahedra (already formed) arrange themselves into a connected network.
All this structural information shows that below a certain amplitude of the Coulomb interactions the structure of our toy silica glass has lost all the medium or extended range features and is similar to a soft-sphere glass. From this study it appears that the transition value of $`\lambda `$ between a fragile and a network glass is around 0.5. Of course the only proper way to describe the covalent bonding network of silica atoms would be to use quantum mechanical calculations and therefore no precise physical meaning should be attached to this threshold value $`\lambda =0.5`$. This value is linked to the specific form of the BKS potential which does not contain 3-body terms and in which the Coulomb interactions are sufficient to build the SiO4 tetrahedra. Nevertheless the good quality of this “simple” potential has been demonstrated in a recent study which has shown that SiO2 glass samples generated with the BKS potential are good starting configurations for ab initio quantum mechanical molecular dynamics calculations. One interesting consequence of our results is that one should be cautious when screened Coulomb interactions are used to simulate real silica glasses. Indeed screened interactions are in a sense similar to “weak” interactions like the ones we have simulated here with the small values of $`\lambda `$ and mistreat the sum in reciprocal space which is made in the Ewald summation method. We show here that the screening should not be too strong in order to simulate a realistic silica glass.

In parallel with the study of the structure we have calculated the vibrational spectrum as a function of $`\lambda `$. During the diagonalization of the dynamical matrix we did not calculate the eigenvectors (for computer time reasons) but used the eigenvalues (always positive) to compute the vibrational density. We were of course interested in the Boson Peak which appears to be a characteristic of glasses even though it can also exist in some crystals . This peak, seen in Raman and inelastic neutron scattering spectra, reflects an enhancement in the density of states compared to the Debye distribution and is located at frequencies around $`\nu _B1.5`$ THz in silica. When varying $`\lambda `$, the energy (and the pressure since we consider a microcanonical ensemble) of our system varies appreciably and this is reflected in the variation of the largest frequency of the whole vibrational spectrum, $`\nu _{\mathrm{max}}`$, which increases from $`\nu _{\mathrm{max}}`$ 20 THz for $`\lambda =0.1`$ to $`\nu _{\mathrm{max}}`$ 40 THz for $`\lambda =1.0`$. As expected, also the shape of the spectrum varies with $`\lambda `$: for $`\lambda `$ small we have basically one band since the oxygen and silicon atoms are almost identical whereas for large values of $`\lambda `$ we observe a clear separation between an acoustic and an optical band typical of the formation of well defined Si-O bonds. But the Boson peak is in all cases clearly present. However due to the $`\lambda `$ dependence of the whole spectrum, it is not legitimate to compare the BPs directly and therefore we have first determined the normalized density of states $`g(f)`$ using reduced frequencies $`f=\nu /\nu _{\mathrm{max}}`$. The normalization is such that $`_0^1g(f)𝑑f=1`$. Then we have used a standard way to extract the BP from the vibrational spectrum which is to plot $`g(f)/f^2`$ because in the Debye approximation $`g(f)f^2`$ at low frequency. This curve has then been fitted by a “generalized Lorentzian” :

$$I=I_0f^n\frac{1}{[f^2+f_0^2]^m}$$
(3)
which has been used to extract the position of the BP and the corresponding height in a straightforward manner. In figure 4 we show the $`g(f)/f^2`$ data together with the fitting curves for three typical values of $`\lambda `$ and in the inset we report the position and the intensity of the BP as a function of $`\lambda `$. As expected, due to the normalization of the $`g(f)`$ curves, the variations of the position and the intensity have opposite trends. Apart from a first decrease of the intensity between $`\lambda =0.1`$ and $`\lambda =0.2`$ that we cannot interpret simply (at very small values of $`\lambda `$, there are some artifacts due to the form of the potential as indicated above), one observes a net sharpening of the BP when $`\lambda `$ becomes larger than $`0.5`$. This seems to indicate that a strong BP is linked to the existence of tetrahedral units or to the presence of an FSDP and therefore supports previous interpretations for the origin of the BP in silica . Moreover the global variation of the BP characteristics with $`\lambda `$ confirms the recent experimental observation that the BP increases when going from fragile to strong glasses . Of course these results, based on a comparison between reduced spectra, should be taken with care: in addition to the increase of the bandwidth when increasing $`\lambda `$ there is also a deformation of the whole spectrum which is difficult to take into account simply. Finally, if some significance can be given to the apparent decrease of the BP intensity with increasing $`\lambda `$ for $`\lambda >0.7`$, one could conclude that, once the tetrahedral network is formed, the increase of the medium range order leads to a decrease of the BP intensity. Of course, such a statement needs to be validated by further investigations.

## IV Conclusion

In conclusion, we have performed classical molecular dynamics simulations in a toy silica glass model in which we have changed the intensity of the long-range Coulomb forces. With this procedure it is possible to analyze the structure and the dynamics of a fragile glass (without Coulomb interaction) and a strong glass (with Coulomb interaction) within the same system. We have shown that the formation of a tetrahedral network, the existence of a First Sharp Diffraction Peak in the static structure factor, a characteristic peaked ring size distribution, as well as a relatively sharp Boson peak, are correlated to the strength of the Coulomb forces: all these features are absent in the fragile glasses and appear once the Coulomb forces exceed a potential-dependent threshold. Of course our study is limited to a very specific model and other ways of creating locally ordered units could have been explored. Nevertheless it sheds light on what is really happening in a system which can be considered as either a fragile or a strong glass. By tuning one single parameter (which therefore appears to be the predominant ingredient in the simulation of glasses) one is able to connect structural features like the FSDP to vibrational specificities like the Boson peak and give answers to questions previously raised by experimental studies.

Part of the numerical calculations were done at the computer center CINES (Centre Informatique National de l’Enseignement Superieur) in Montpellier.