text:
# The geometry of atmospheric neutrino production

## 1 Introduction

Atmospheric neutrinos are produced in the hadronic showers generated by cosmic rays in the Earth’s atmosphere. Absorption in the Earth is negligible in the entire relevant energy range, and therefore a detector located near the surface of the Earth will receive a neutrino flux from all directions. In the presence of neutrino oscillations the observed rates of $`\nu `$ interactions will be modified, and the angular distributions distorted. These effects have been in fact measured by Super–Kamiokande  and other atmospheric neutrino detectors  and give clear evidence for the existence of flavor transitions .

The evidence for oscillations is robust, and does not depend on a detailed calculation of the expected neutrino fluxes with and without neutrino oscillations. The strongest evidence comes in fact from the detection of an up–down asymmetry in the angular distribution of muon events and in a small ratio for the rates of $`\mu `$–like and $`e`$–like events. It is possible to predict with very simple considerations that in the absence of oscillations the $`\nu `$ fluxes are approximately up–down symmetric, and the fluxes of electron and muon neutrinos are strictly related because they are produced in the chain decay of the same parent mesons. In order to extract the oscillation parameters from the data it is however important to have detailed predictions for the expected intensity and angular distributions of the no–oscillation fluxes.

In this work I will discuss the angular distributions of the atmospheric neutrino fluxes. These distributions are determined by the presence of neutrino oscillations (if they exist) and a combination of three most important effects:

1. Geomagnetic effects on the primary cosmic rays. Low rigidity particles cannot reach the vicinity of the Earth, and the effect depends on the point of the Earth where cosmic rays arrive (being stronger near the geomagnetic equator), and on their direction (local zenith and azimuth angles).
2. The zenith angle dependence of the neutrino yields. Inclined showers produce more neutrinos that vertical ones because the decay of charged mesons and muons is then more probable.
3. The spherical geometry of the neutrino source volume. This results in an enhancement of the neutrino fluxes from horizontal directions and a suppression for the (up–going and down–going) vertical directions. The effect is exactly up–down symmetric, and is important for sub–GeV neutrinos.

Additional smaller effects are due to the bending of charged particles during shower development, the presence of mountains, and the existence of different air density profiles over different geographycal regions.

The aim of this paper is to discuss in particular the effect of the spherical geometry of the neutrino source region. This effect has been recognized only recently , is not included in the calculations currenty in use  in the analysis of experimental data and is still the object of some controversy.

The work is organized as follows: the next section contains a discussion what is the one–dimensional (1D) approximation used in the first calculations of the atmospheric neutrino fluxes and why the approximation was considered necessary; section 3 and 4 present for completeness very brief qualitative discussion on the geomagnetic effects and of the angular dependence of the neutrino yields; section 5 discusses the geometrical effects related to the shape of the neutrino source volume; section 6 presents the results of a detailed calculations; section 7 gives a summary.

## 2 The one–dimensional approximation

Any cosmic ray shower, produced by a primary particle interacting in an arbitrary point of the atmosphere, with an arbitrary direction can result in one (or more) neutrinos with trajectories that intersect the detector in consideration. Of course only a very small fraction of the cosmic ray showers produce neutrinos that pass inside (or in the vicinity) of a detector, and therefore a full montecarlo calculation of the neutrino fluxes appears as extremely inefficient.

It is because of this difficulty that the first generation calculations of the atmospheric neutrino fluxes have been performed in what is called the one–dimensional approximation. In this approximation one considers only a very small subset of all possible cosmic rays primaries, those that have a trajectory that when continued as a straight line beyond the interaction point intersects the detector. The neutrinos produced in the showers of these primaries are considered as collinear to the primary parent, and also have trajectories that intersect the detector, therefore all $`\nu `$’s generated in a montecarlo calculation can be collected and analysed to estimate the atmospheric neutrino fluxes.

More explicitely, in a 1–D calculation the flux of neutrinos with flavor $`\alpha `$, energy $`E_\nu `$ and direction $`\mathrm{\Omega }`$ observed by a detector located at the position $`\stackrel{}{x}_d`$ is calculated as:

$$\varphi _{\nu _\alpha }(E_\nu ,\mathrm{\Omega }_\nu ,\stackrel{}{x}_d)=\underset{A}{}𝑑E_0\varphi _A[E_0,\mathrm{\Omega }_0\left(\mathrm{\Omega }_\nu \right),\stackrel{}{x}_0(\mathrm{\Omega }_\nu ,\stackrel{}{x}_d)]\frac{dn_{A\nu _\alpha }}{dE_\nu }(E_\nu ;E_0,\mathrm{cos}\theta _0)$$
(1)
the sum is over all primary cosmic rays nuclear species, the quantity $`dn_{A\nu _\alpha }/dE_\nu `$ is the differential yield of neutrinos of flavor $`\alpha `$ from a primary of type $`A`$, energy $`E_0`$ and zenith angle $`\theta _0`$, $`\varphi _A(E_0,\mathrm{\Omega }_0,\stackrel{}{x}_0)`$ is the flux of primary cosmic rays of type $`A`$ , energy $`E_0`$ and direction $`\mathrm{\Omega }_0`$ that reach the Earth at point $`\stackrel{}{x}_0`$. The dependences on the direction and position for the primary flux are the consequence of geomagnetic effects. Since neutrinos travel along straight lines, the interaction position of the primary flux is determined to a good approximation by the neutrino direction $`\mathrm{\Omega }_\nu `$. All down–going neutrinos ($`\mathrm{cos}\theta _\nu >0`$) are produced in showers in the general vicinity of the detector: $`\stackrel{}{x}_0\stackrel{}{x}_d`$, while for up–going trajectories ($`\mathrm{cos}\theta _\nu <0`$) the primary interaction point is in the vicinity of the point where the neutrino enters the Earth. The crucial approximation made in the 1–D approximation, is to consider also the direction of the primary $`\mathrm{\Omega }_0`$ as determined by $`\mathrm{\Omega }_\nu `$. For down–going neutrinos this simply means $`\mathrm{\Omega }_0\mathrm{\Omega }_\nu `$, and for up–going neutrinos, expressing $`\mathrm{\Omega }_0`$ as a zenith and azimuth angle with respect to the local vertical, this means: $`\mathrm{\Omega }_0(\mathrm{cos}\theta _0,\phi _0)=(\mathrm{cos}\theta _\nu ,\phi _\nu )`$.

Equation (1) simplifies enormously the calculation of the atmospheric neutrinos fluxes because now only a very small (indeed infinitesimal) fraction of the cosmic ray showers has to be studied. About half of these showers have trajectories that reach the Earth near the detector point (arriving from all directions in a $`2\pi `$ solid angle), the other half is distributed over the entire surface of the Earth, but a unique direction has to be considered for each point (see fig. 1).

In the equation (1) the problem of the calculation of the neutrino is “factorized” in three parts:

1. Determination of the primary cosmic ray fluxes. These fluxes have of course to be measured experimentally. The primary fluxes have also a time dependence due to the solar modulation. The effect decreases with momentum, and becomes negligible for rigidities $`p/q\stackrel{>}{_{}}100`$ GV.
2. Calculation of the geomagnetic effects. The primary cosmic ray fluxes reaching point $`\stackrel{}{x}_0`$ from the direction $`\mathrm{\Omega }_0`$ can be calculated from the isotropic flux unperturbed by geomagnetic effects (that is measured in the Earth’s magnetic polar regions) with a knowledge of the structure of the geomagnetic field, testing if the trajectories are allowed or forbidden.
3. Calculation of the neutrino yields. The number of neutrinos (and their energy spectrum) produced by a primary of energy $`E_0`$ and zenith angle $`\theta _0`$ can be calculated studying the average development of hadronic showers in the atmosphere.

A calculation using a three–dimensional method can still be considered as divided into the three parts discussed above, and indeed the first two steps are identical. The calculation is however much more complex because it is not anymore sufficient to consider for the neutrino yield simply the inclusive energy distribution of the neutrinos produced in a shower, but also the angular distributions of the neutrinos (and their strong correlation with the energy) are important and have to be calculated. Moreover the simple integral in (1) performed over the primary energy $`E_0`$ keeping the direction fixed has now to be replaced by a more complex convolution that involves not only the energy but also the direction of the primary particle, since the neutrinos with a given trajectory can be produced in showers with a non collinear axis.

## 3 Geomagnetic effects

An assumption common to all atmospheric neutrino calculation, is that the fluxes of cosmic rays at one astronomical unit of distance from the sun, when they are not perturbed by the near presence of the Earth are isotropic, and can be simply described by their energy dependence $`\varphi _A\left(E_0\right)`$. The fluxes reaching the surface of the Earth are however affected by the geomagnetic field and therefore are not isotropic, and have different intensities in different locations. For example it is intuitively clear that low rigidity particles can reach the Earth only traveling parallel to the magnetic field lines arriving near the magnetic poles, while high rigidity particles can reach all points on the Earth’s surface from all possible directions.

Given a map of the magnetic field around the Earth it is a straightforward exercise to compute numerically if a given primary particle three–momentum and position correspond to an allowed or a forbidden trajectory. It is sufficient to integrate the equation of motions of the charged particle in the geomagnetic field, and see if the past trajectory of the particle intersects the Earth’s surface, remains confined to a finite distance from the Earth (forbidden trajectory); or originates from large distances (allowed trajectory).

For a qualitative understanding it can be useful to consider the historically important approximation of describing the Earth geomagnetic effect as an exactly dipolar field. In a dipolar field the problem can be solved analytically. All positively particles with rigidity $`R>R_S^+`$ are allowed and the trajectories of all particles with $`R<R_S^+`$ are forbidden because they remain confined to a finite distance from the dipole center 111In this solution it is assumed that the field fills the entire space. A fraction of the allowed trajectories have really a segment “inside” the Earth (with $`r<R_{}`$) and therefore should be considered forbidden also in a dipole field.. The quantity $`R_S^\pm (\stackrel{}{x},\widehat{n})`$ is the Störmer rigidity cutoff :

$$R_S^+(r,\lambda _M,\theta ,\phi )=\left(\frac{M}{2r^2}\right)\left\{\frac{\mathrm{cos}^4\lambda _M}{\left[1+\left(1\mathrm{cos}^3\lambda _M\mathrm{sin}\theta \mathrm{sin}\phi \right)^{1/2}\right]^2}\right\}$$
(2)
where we have made use of the cylindrical symmetry of the problem, $`M`$ is the magnetic dipole moment of field, $`r`$ is the distance from the dipole center, $`\lambda _M`$ the magnetic latitude, $`\theta `$ the zenith angle and $`\phi `$ an azimuth angle measured counterclockwise from magnetic north. For negatively charged particles the cutoff is obtained with the reflection $`\phi \phi +\pi `$, that is exchanging east and west.

The qualitative features that are important for our discussion are the following:

1. For a fixed direction the cutoff rigidity grows monotonically from a vanishing value at the magnetic pole to a maximum value at the magnetic equator.
2. The cutoffs for particles traveling toward magnetic west are higher than those for particles traveling toward magnetic east. Note that geomagnetic effects are the only mechanism discussed in this work that can generate a non–flatness in the azimuthal angle distributions of atmospheric neutrinos.
3. The highest cutoff corresponds to westward going, horizontal particles reaching the surface of the Earth at the magnetic equator ($`\phi =90^{}`$, $`\theta =90^{}`$, $`\lambda _M=0^{}`$). The maximum rigidity cutoff is approximately 60 GV. All geomagnetic effects vanish for particles above the rigidity value.

## 4 The neutrino yields

The “neutrino yield” is the average number of neutrinos of a certain flavor produced by a primary cosmic ray particle. It depends on the primary particle type, energy, and zenith angle. There are also weaker dependences on the azimuth angle of the shower and its geographycal locations due to the effects of the geomagnetic field on the shower development.

In the 1–D approximation the angles of the secondary particles with respect to the primary directions are neglected (or integrated over), and therefore the neutrino yields can be described simply as a set of functions that give the energy distribution of the neutrinos (of different flavors) produced by primaries of given type, energy, and zenith angle.

The one–dimensional neutrino yields have been obtained with the integration (numerical  or even analytical in simplified tratments ) of a set of differential equations that describe the average development of hadronic showers. The most accurate calculations of the yields are however performed with montecarlo techniques. A large number of showers (for primary particles of different type, energy and zenith angle) is generated, and the number, flavor and energy of all produced neutrinos is recorded to compute numerically the neutrino yields as a function of primary energy, direction and particle type.

The authors of the montecarlo calculations of the yields  have also usually made the approximation of considering all secondary particles in the shower as collinear to the primary particle. This can be obtained “rotating” all final state particles so that their momentum is parallel to the projectile (or parent) particle, and neglecting multiple scattering and bending in the geomagnetic field. Strictly speaking this rotation of course implies a small (negligible) deviation from exact conservation of the longitudinal momentum. In the procedure energy is exactly conserved.

The neutrino yields have some important dependence on the zenith angle of the primary particle. The yields grow monotonically when the zenith angle of the primary particle changes from the vertical to the horizontal direction. The reasons for this growth are simple to understand, and originate from two effects. The first effect is simply that an inclined shower develops in air for a longer distance before hitting the ground, therefore the muons produced in high zenith angle showers have more time to decay in air and produce more neutrinos. Essentially all muons that hit the ground rapidly lose all their energy in ionization and radiation processes and decay at rest producing very soft neutrinos (with $`E_\nu <m_\mu /2`$) that can be neglected because are below the thresholds of the existing detectors.

A second effect, contributing also to an enhancement of neutrino production for horizontal showers is due to the fact that the air density is not constant but decreases (approximately exponentially in the stratosphere) with increasing altitude. Therefore the inclination of a trajectory (assuming the same starting point) determines the column density $`X`$ (g cm-2) that corresponds to a fixed length $`L`$ (cm). The column density is largest for vertically down–going particles (zenith angle $`\theta =0`$) and decreases monotonically with increasing zenith angle. The relation between $`L`$ and $`X`$ determines the relative probability of decay or interactions for charged pions (or kaons), and the energy loss of a muon before decay. A smaller $`X`$ (for a fixed $`L`$) implies that the interaction of weakly decaying mesons is suppressed and their decay enhanced; it also means that the muons will decay with a higher energy. For inclined showers the first effect results in a higher number of neutrinos, and the second one in a slightly harder spectrum of neutrinos

The ratio between the vertical and horizontal yields is a function of the neutrino energy. For low $`E_\nu `$ the ratio is close to unity, and there is very little enhancement for the horizontal directions. This can be qualitatively easily understood. Low energy neutrinos are produced in the decay of low energy pions and low energy muons. Because of relativistic effects the decay length of unstable particles is proportional to the particle momentum, for low energy particles the decay is so rapid, that the decay probability is unity, and the energy loss before decay is negligible independently from the direction of the particle. In these circumstances there is no enhancement for the horizontal directions, and the yield is isotropic. With increasing energy the decay length increases, and the effects outlined above start to be significant. As an illustration,the muon decay length is

$$L_\mu =\tau _\mu \frac{p_\mu }{m_\mu }=6.23p_\mu \left(\mathrm{GeV}\right)\mathrm{Km}$$
(3)
The muons are produced at an average altitude of $`17.5`$ Km ($`30`$ Km) for vertical (horizontal) particles with only a weak energy dependence. For $`p_\mu <1`$ GeV, most muons decay independently from their direction. For $`p_\mu =3`$ GeV (10 GeV) the decay of the vertical muons is suppressed by a factor $`2`$ ($`4`$) while the decay probability of horizontal particles remains approximately unity.

These effects are carefully included in the existing calculations, that in fact give neutrino zenith angle distributions that are approximately flat for low $`E_\nu `$ (with distortions that reflect only geomagnetic effects and the angular distribution of the primary cosmic ray flux), and develop a stronger and stronger enhancement on the horizontal with increasing energy.

## 5 Spherical geometry of the source

Most atmospheric neutrinos are produced in a spherical shell of air at an altitude between 10 and 40 km. The approximately spherical geometry of the source region has some very important effects on the angular distribution of the neutrinos and results in an enhancement of the neutrino flux from the horizontal directions. This enhancement has been overlooked in all calculations of the atmospheric neutrino fluxes before the work of Battistoni et al. .

To illustrate the nature of the new effect, let us consider a situation where the geomagnetic effects are absent (and therefore the flux of primary cosmic rays is exactly isotropic), and the zenith angle dependence of the neutrino yield is negligible, that is the number of neutrinos produced by a primary particle does not depend on inclination of the shower in the atmosphere (as it is the case for low neutrino energy). In this situation a one–dimensional calculation predicts an isotropic neutrino flux. This result is incorrect. A realistic (three–dimensional) calculation results in an angular distribution of the neutrino flux that is exactly symmetric for up–down reflections, and for rotations around the vertical axis but that is not isotropic and exhibits an enhancement for the horizontal directions. The enhancement for the horizontal directions depends on the neutrino energy and is more and more marked with decreasing energy.

The origin of this result can appear surprising, and for “pedagogical reasons” it can be instructive to consider two simple problems that have clear analogies with the emission of neutrinos in thre atmosphere. In problem A an observer is located inside a thin spherical shell of stars of radius $`R_s`$; in problem B the observer is inside a spherical cavity in a black body with temperature $`T`$. In both cases we want to compute the angular distribution of the photon flux measured by the observer.

### 5.1 Problem A: a shell of stars (isotropic emission)

Let us consider an observer inside a thin shell composed of isotropically emitting stars (see fig. 2). Let us assume that $`L_0`$ (erg/s) is the average luminosity of the stars, $`n`$ is the surface number density (cm-2) of stars in the shell, and $`R_A`$ (with $`R_A<R_s`$) the distance of the observer from the center of the star shell. The energy flux per unit solid angle $`dF\left(\mathrm{\Omega }\right)/d\mathrm{\Omega }`$ (erg/(cm2 s sr)) measured by the observer can be calculated as:

$$\left(\frac{dF\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}\right)_{\mathrm{iso}}=n\frac{L_0}{4\pi }\frac{1}{\mathrm{cos}\theta _e}=n\frac{L_0}{4\pi }\left[1\left(\frac{R_A}{R_s}\right)^2\left(1\mathrm{cos}^2\theta \right)\right]^{\frac{1}{2}}$$
(4)
where $`\theta `$ (see fig. 2) is a polar (“zenith”) angle around the axis that passes through the center of the spherical shell and the observer.

The solution respects the cylindrical symmetry of our problem (does not depend on the azimuthal angle $`\phi `$) and is up–down symmetric (it remains identical after a reflection $`\theta \pi \theta `$), but neglecting the special case of the observer at the center of the shell, it is not spherically symmetric, and is peaked for horizontal directions ($`\mathrm{cos}\theta =0`$). When the observer is close to the emitting shell (that is when $`R_A`$ approaches $`R_s`$) the peaking is very strong, and the observer sees the shell as a narrow bright disk. As a numerical example: for $`R_AR_{}=6371`$ Km and $`R_s=R_{}+15`$ km, formula (4) predicts a horizontal flux 14.6 times more intense than the vertical flux (see fig. 3).

It is elementary to deduce equation (4). Only stars in the surface element $`dA`$ subtended by the solid angle $`d\mathrm{\Omega }`$ contribute to the energy flux from that solid angle interval with a total of $`ndA`$ stars (see fig. 2); each star contributes an energy flux $`L_0/\left(4\pi \mathrm{}^2\left(\theta \right)\right)`$, where $`\mathrm{}\left(\theta \right)`$ is the distance between the observer and the sources in the direction $`\theta `$. Therefore:

$$\left(dF\right)_{\mathrm{iso}}=ndA\frac{L_0}{4\pi \mathrm{}^2}=n\frac{d\mathrm{\Omega }\mathrm{}^2}{\mathrm{cos}\theta _e}\frac{L_0}{4\pi \mathrm{}^2}$$
(5)
In the second equality we have written explicitely the expression for the area of the element of shell seen in the solid angle $`d\mathrm{\Omega }`$: $`dA=d\mathrm{\Omega }\mathrm{}^2/\mathrm{cos}\theta _e`$. The area of the shell element obviously scales as $`\mathrm{}^2`$, and this factor cancels exactly the factor $`\mathrm{}^2`$ that takes into accout the reduction with distance of the apparent luminosity of the stars, however the expression for $`dA`$ contains also a geometrical term $`\left(\mathrm{cos}\theta _e\right)^1`$ that takes into account the orientation of the surface element of the shell with respect to the line of sight. The “emission” angle $`\theta _e`$ is the angle between the normal to the source surface (toward the center of the shell) and the line of sight from the source to the observer. From elementary geometry one obtains:

$$\mathrm{cos}\theta _e=\sqrt{1\left(\frac{R_A}{R_s}\right)^2\left(1\mathrm{cos}^2\theta \right)}$$
(6)
It is this factor that is responsible for the strong enhancement of the flux from the horizontal direction.

It is interesting to note that the up–going and down–going vertical fluxes ($`\mathrm{cos}\theta =\pm 1`$) are independent from the position of the observer inside the shell, while the horizontal flux depends on the observer position, and grows when the observer approaches the emitting spherical shell. It follows that the total (angle integrated flux) is a function of the observer position:

$$F_{\mathrm{tot}}^{\mathrm{iso}}=_{\left[4\pi \right]}𝑑\mathrm{\Omega }\left(\frac{dF\left(\mathrm{\Omega }\right)}{d\mathrm{\Omega }}\right)_{\mathrm{iso}}=nL_0\frac{1}{r}\mathrm{sinh}^1\left[\frac{r}{\sqrt{1r^2}}\right]$$
(7)
(where $`r=R_A/R_s`$) and grows monotonically from a value $`nL_0`$ for an observer at the center of the shell to a divergent value when $`r1`$. The divergence is connected to the fact that an observer exactly on the shell must be (in the approximation we have used of a shell of vanishing thickness) infinitesimally close to a star.

### 5.2 Problem B: a cavity in a blackbody (emission $`\mathrm{cos}\theta _e`$)

The solution of problem B, is of course well kown, we actually live in a (very large) cavity whose walls have a temperature of 2.7 Kelvin. An observer inside a spherical “cavity” in a black body (and in fact in a cavity of arbitrary shape and dimension) observes an isotropic black body spectrum independently from its position. It can be instructive to deduce this well known result from the same considerations used in the previous discussion. The surface of a black body emits energy per unit surface and unit solid angle at a rate:

$$\frac{dL_{bb}}{d\mathrm{\Omega }_e}=\sigma T^4\frac{\mathrm{cos}\theta _e}{\pi }$$
(8)
where $`\sigma `$ is the Stefan–Boltzmann constant. The emission is not isotropic but decreases linearly with the cosine of the emission angle with respect to the normal to the surface element. The energy flux from the direction $`\mathrm{\Omega }`$ for the observer of problem B can be calculated as:

$$\left(dF\right)_{bb}=dA\frac{dL_{bb}}{d\mathrm{\Omega }}\frac{1}{\mathrm{}^2}=\left(\frac{d\mathrm{\Omega }\mathrm{}^2}{\mathrm{cos}\theta _e}\right)\left(\sigma T^4\frac{\mathrm{cos}\theta _e}{\pi }\right)\left(\frac{1}{\mathrm{}^2}\right)=\frac{\sigma T^4}{\pi }d\mathrm{\Omega }$$
(9)
Note that in this case there is a cancellation not only for the factor $`\mathrm{}^2`$, but also for the orientation factor $`\mathrm{cos}\theta _e`$. If the normal to emitting surface has a large angle with the line of sight, a larger surface ($`\left(\mathrm{cos}\theta _e\right)^1`$) can contribute to the emission, but the large angle emission is suppressed by a factor $`\mathrm{cos}\theta _e`$. Combining the two effects one has an exact cancellation. The demonstration has been actually general, and is valid for a cavity of arbitrary shape or dimension and for an arbitrary position of the observer.

### 5.3 Atmospheric neutrinos

We can now come back to the problem of atmospheric neutrinos. The source volume of atmospheric neutrinos is also with very good approximation a spherical shell with the observers placed inside, very close to the inner radius of the source volume. However an element of atmosphere emits neutrinos with an angular distribution that is neither isotropic, nor linear in zenith angle, therefore we can expect a neutrino flux at sea level with an angular distribution with a form that is intermediate between the solution for problem A (isotropic emission in the source and strong horizontal peaking of the observed flux) and problem B (emission $`\mathrm{cos}\theta _e`$ in the source and isotropic observed flux).

Let us consider an isotropic flux of cosmic rays that impinges on the Earth atmosphere. The quantity $`C_0`$ that is the number of cosmic rays absorbed per unit time by an element of unit area of the atmosphere (units cm-2 s-1) can be calculated as:

$$C_0=_{\left[2\pi \right]}𝑑\mathrm{\Omega }\varphi _0\mathrm{cos}\theta _0=\pi \varphi _0$$
(10)
(the integration is over one hemisphere). Note that the absorption rate of cosmic rays is proportional to $`\mathrm{cos}\theta _0`$, with $`\theta _0`$ the angle between the normal to the surface element and the cosmic ray direction. If the cosmic rays produce an average number $`n_\nu `$ of neutrinos per showers, then each element of the atmosphere is also a source of neutrinos, with an emission rate $`S_\nu =C_0n_\nu `$ (again in units cm-2 s-1).

The angular distribution of the neutrino emission from a “patch” of atmosphere is not easy to calculate and depends on the average angle between neutrino and primary particle. Knowing this distribution it is possible to calculate the observable flux at sea level, following the same steps outlined above.

The results for two limiting cases are easy to obtain. If the neutrinos in a shower are produced collinearly with the primary trajectory ($`\mathrm{\Omega }_\nu =\mathrm{\Omega }_0`$) then the angular distribution of the neutrino emission from a surface element of the atmosphere is $`\mathrm{cos}\theta _\nu `$, because it simply reflects the angular distribution of the absorbed primary flux (see equation (10)). From equation (9) (neglecting the angular dependence of the neutrino yield) it follows that the observed neutrino flux is isotropic.

The opposite limiting case is obtained when then neutrinos are emitted quasi isotropically in the primary cosmic ray showers. This is approximately true only for very low neutrino energy, ($`E_\nu \stackrel{<}{_{}}10`$ MeV). In this case all memory of the direction of the primaries is erased, and the neutrino emission from an element of the atmosphere is isotropic. From equation (4) it follows that the angular distribution of the observed neutrino fluxes is very sharply peaked for the horizontal directions. Considering the altitude distribution of the neutrino production points (that is related to the density profile of the atmosphere and the value of the hadronic cross sections) it can be estimated that for the very low energy neutrinos that are emitted quasi–isotropically, the flux on the horizontal plane is more than an order of magnitude more intense than in the vertical directions (see fig. 3).

In general neutrinos are emitted in a cone with an axis that corresponds to the primary particle direction and an opening angle that shrinks with increasing neutrino energy. We can therefore expect that the emission of the neutrinos from the atmosphere is in general intermediate, between the two extremes corresponding to isotropic emission, and an emission $`\mathrm{cos}\theta _e`$. The two limiting cases are approached for very low neutrino energy, when the neutrino emission is very poorly correlated to the primary direction, and for high energy neutrinos when because of the Lorentz boost, the neutrinos are emitted quasi–collinearly with the primary particle.

These expectations can be verified with more detailed calculations. As an illustration, let us consider the simplified problem where (i) the atmosphere is formed by a single thin layer at an height $`h`$ (the radius of the emitting layer is therefore $`R_s=R_{}+h`$); (ii) the primary cosmic ray flux is exactly isotropic; (iii) the neutrino production is independent of the zenith angle of the primary particle; and (iv) the angle $`\theta _{0\nu }`$ between each neutrino and its primary particle has a fixed value $`\theta _{0\nu }=\alpha `$ (and the emission has cylindrical symmetry around the axis defined by the primary particle trajectory). The problem is to compute the the neutrino flux received from an observer located at sea level. This simplified problem contains essentially all the interesting geometry of atmospheric neutrino production . Defining $`r=R_{}/R_s`$, and $`S_\nu =\pi \varphi _0n_\nu `$ the number of neutrinos emitted per unit time and unit surface by the emitting layer, the solution can be written as:

$$\varphi _\nu (\mathrm{\Omega }_\nu ,r,\alpha )=\frac{S_\nu }{2\pi }\left[\frac{F_\alpha \left(y\right)}{y}\right]_{y=\sqrt{1r^2\left(1\mathrm{cos}^2\theta _\nu \right)}}$$
(11)
with (for $`\alpha <90^{}`$):

$$F_\alpha (y)=\{\begin{array}{cc}0\hfill & \text{for }y<\mathrm{sin}\alpha \text{,}\hfill \\ & \\ \mathrm{cos}\alpha y+\frac{2}{\pi }\sqrt{1\mathrm{cos}^2\alpha y^2}\frac{2}{\pi }\mathrm{cos}\alpha y\mathrm{tan}^1\left[\frac{\mathrm{cos}\alpha y\sqrt{1\mathrm{cos}^2\alpha y^2}}{\mathrm{cos}^2\alpha +y^21}\right]\hfill & \text{for }y<|\mathrm{sin}\alpha |\text{ ,}\hfill \\ & \\ 2\mathrm{cos}\alpha y\hfill & \text{for }y>\mathrm{sin}\alpha \text{.}\hfill \end{array}$$
(12)
For $`\alpha >90^{}`$ one has to substitute: $`F_\alpha \left(y\right)=F_{\pi \alpha }\left(y\right)`$. This solution can be easily checked numerically with a montecarlo method writing a simple program with few lines of computer code. Some examples of the solution are shown in fig. 4 and fig. 5. For $`\theta _{0\nu }0`$ the observed neutrino flux is isotropic with a value $`\varphi _\nu =\varphi _0n_\nu `$ where $`\varphi _0`$ is the isotropic primary flux and $`n_\nu `$ the average number of neutrino produced by a primary particle. When the angle $`\theta _{0\nu }`$ grows the observed flux is suppressed for the vertical direction and is enhanced on the horizontal plane. The suppression and enhancement become more important when the angle $`\theta _{0\nu }`$ increases. For example (see fig. 4) if the neutrinos are produced at an height $`h=20`$ km and are emitted with an angle $`\theta _{0\nu }=10^{}`$, $`30^{}`$, or $`60^{}`$ with respect to the primary particle direction, the observed neutrino flux at sea level has a horizontal/vertical ratio of 1.45, 3.25, and 8.5.

It is important to observe that the enhancement for the horizontal directions for a fixed angle $`\theta _{0\nu }`$ depends on the ratio $`r=R_{}/R_s`$ (that is on the altitude of the neutrino production), and is stronger when the altitude decreases. On the contrary the suppression for the vertical direction does depend on the angle $`\theta _{0\nu }`$ but is independent of the ratio $`R_{}/R_s`$ (see equation (12) and fig. 5). An important consequence of this is the fact that the average (angle integrated) neutrino flux observed at sea level depends on the altitude of neutrino production, and grows when the average altitude of production decreases. This can easily understood qualitatively observing that for finite $`\theta _{0\nu }`$ a fraction of the produced neutrinos will “miss” the Earth because are emitted on trajectories that do not intersect her surface. The fraction neutrinos “lost in space” will decrease if the production points are closer to the surface.

It is less intuitive to note that for large $`\theta _{0\nu }`$ the flux of neutrinos in a full 3–D calculations can be higher than what is obtained in a collinear approximation. This can be checked integrating over angle equation (11) and it is possible because the enhancement on the horizontal plane can be stronger than the suppression on the vertical. This does not represent a violation of “unitarity” since it can happen only for a small range of radii just below the source region in the atmosphere, while in most of the volume inside the Earth for larger $`\theta _{0\nu }`$ the flux is suppressed.

From fig. 3 and fig. 5 one can see that the height of neutrino production is very important in determining the enhancement for the horizontal directions. To compute the geometrical effects for atmospheric neutrinos it is not possible to approximate the emission volume as a thin surface. The results for a thick shell can of course be obtained integrating over a continuos distributions of thin shells with different radii, The assumption of emission from an altitude of $`15`$–20 km can also give a reasonable first order approximation.

In summary this “geometrical effect” can be very important for low neutrino energy. In a realistic calculation it must be of course be combined with the other effects: geomagnetic effects and the zenith angle dependence of the neutrino yields to determine the angular distribution of atmospheric neutrinos.

It has been stated that the enhancement on the horizontal is the result of the inclusion in a 3–D calculation of primary particles with “grazing” trajectories that if continued would not intersect the Earth’s surface and are not included in a 1–D calculation. These grazing trajectories give a neglible contribution to the horizontal neutrino fluxes, and are not the source of the enhancement. The large flux from the horizontal directions can be understood as the effect of the production of “horizontal” neutrinos from the more numerous “vertical” showers after emission with an appropriate angle with respect to the shower axis. Note than an isotropic flux implies that most of the interacting tracks are in fact vertical (the absorption rate is $`\mathrm{cos}\theta _0`$), and therefore the opposite effects of “vertical” showers producing “horizontal” neutrinos and vice versa are of different size and do not cancel each other. Note also that the enhancement for the horizontal directions is linked with a suppression of the vertical fluxes.

A second 3–D calculation of the atmospheric neutrino fluxes (Tserkovnyak et al. ) has been recently made public. In this calculation the enhancement of the neutrino fluxes for the horizontal direction is not present, and the predicted angular distributions of the neutrinos are very similar to those calculated in a 1–D approach. The reasons for the discrepancy with the Battistoni et al work  are not easily understandable from a simple reading of the papers. Also an early attempt at a 3–D calculation of Lee and Koh  did not find a flux enhancement for horizontal directions. As we are trying to illustrate in this work the enhancement of the neutrino fluxes in the horizontal directions is the consequence of simple geometry and must be present in a correctly performed 3–D calculation. The absence of this effect in  and  is therefore evidence of the existence of errors it the calculation method used for the prediction.

## 6 A complete three-dimensional calculation

To illustrate the points discussed in the previous section with a concrete example we have performed a detailed calculation of the neutrino fluxes with a three–dimensional method.

The method of the calculation is conceptually very simple. The space around the Earth has been divided into two regions: an outer region ($`r>R_{}+\mathrm{\Delta }R`$ where $`R_{}=6371.2`$ Km is the Earth’s average radius and $`\mathrm{\Delta }R=80`$ km) where the matter density is assumed to be vanishingly small and only a magnetic field is present, and an inner region ($`R_{}<r<R_{}+\mathrm{\Delta }R`$) where particles interact not only with the magnetic field but also with air. The air is modeled as having a spherically symmetric density profile $`\rho \left(r\right)`$. The value of $`\mathrm{\Delta }R`$ was determined calculating that less that one percent of the cosmic rays that graze the Earth with a distance of closest approach equal to $`\mathrm{\Delta }R`$ interact in the residual high altitude atmosphere. The result of the calculation are independent of $`\mathrm{\Delta }R`$ (if it is sufficiently large), but the calculation become inefficient for $`\mathrm{\Delta }R`$ too large,

The fluxes of cosmic rays entering the spherical surface at $`R_s=R_{}+\mathrm{\Delta }R`$ were determined using a model for the isotropic flux in the absence of geomagnetic effects (we used the results of primary flux of the Bartol calculation ), and taking into account the effects of the geomagnetic field (we used the IGRF field model ) in the propagation of the particles. This was obtained generating a uniform and isotropic flux that enters the spherical surface $`R_s`$, and rejecting the particles that correspond to forbidden trajectories.

The primary particles are then propagated inside the inner shell, taking into account the bending in the magnetic field and the interactions with the air nuclei. When an inelastic interaction occurs, a set of final state particles is generated with a realistic distributions, of multiplicity, flavor, energy and transverse momentum. In this calculation we have used a simple model of the hadronic interactions due to Hillas . A fraction of approximately 3% of the primary cosmic rays that enter the inner shell exit without interacting.

The subsequent development of the shower is a standard one. The trajectories of charged particles are propagated along curved trajectories because of the presence of the magnetic field. When neutrinos are produced their trajectories (simple straight lines) are studied. The trajectories intersect $`n`$ times the surface of the Earth with $`n=0`$ or $`n=2`$. If $`n=0`$ the neutrino is discarded. If $`n=2`$ two neutrinos (corresponding to down–going and up–going trajectories) are “collected” and their position and direction are recorded.

Some of the results of our calculation are collected in six figures (from fig. 6 to fig. 11). For illustration (and debugging purposes) we have performed the calculation three times.

* A first calculation (represented by the thick) histograms was performed using a fully 3–D method. This includes geomagnetic effects on the primary cosmic ray fluxes, and in the shower development the inclusion of realistic $`p_{}`$ distributions in hadronic interactions, and the correct treatment of the kinematics in particle decays.
* A second calculation (represented by thin histograms) was performed to reproduce the 1–D algorithms. The geomagnetic effects are calculated for the primary cosmic ray particles (exactly as in the previous case) but not for the shower development (particles travel along straight lines for $`r<R_{}+80`$ Km.), and all final state particles are collinear with the projectile (or parent) particle. This is achieved modeling the interactions and particle decays exactly as in the previous case, (including therefore transverse momentum), and performing as a last step a rotation of all the 3–momenta of the final state particles so that they become parallel to the projectile (for interactions) or parent (for decays) particle.
* A third calculation (represented by dashed histograms) was performed neglecting the geomagnetic effects on the primary flux (therefore considering exactly isotropic primary fluxes) and using the 1–D algorithms outlined in the previous point.

Since the atmospheric neutrino fluxes depend on the geographycal position of the detector, we have divided the Earth’s surfaces into five equal area regions, according to the geomagnetic latitude: two polar regions, two intermediate latitude regions and an equatorial region, and have collected separately all neutrinos generated in the montecarlo calculation that land in the different regions (each neutrino is recorded in two separate regions, once as down–going and once as up–going). The results in the two (north and south) polar regions and the two intermediate latitude regions are essentially undistinguishable, and therefore only the results for the north regions are presented here.

The six figures (from 6 to 11) with the results of our calculation are organized as follows: three figures are for $`\mu `$–like events and three for $`e`$–like events; in each set the three figures correspond to the three different magnetic latitude regions: polar ($`\mathrm{sin}\lambda _M=[0.6,1]`$ where $`\lambda _M`$ is the magnetic latitude), intermediate ($`\mathrm{sin}\lambda _M=[0.2,0.6]`$) and equatorial ($`\mathrm{sin}\lambda _M=[0.2,+0.2]`$). In each figure four panels show the zenith angle distributions for the event rates in four different neutrino energy intervals. To compute the rates the $`\nu `$ fluxes have been convoluted with the neutrino cross sections model described in . Note that the zenith angle is the neutrino one, and no experimental smearing or inefficiency has been included.

The motivation to perform the three different calculations are that a comparison between the three allows to put in evidence the different effects that determine the shape of the zenith angle distributions of the neutrino fluxes. In fact:

1. In the calculation (iii) (1–D approximation with no geomagnetic effects) the only source of a non–flatness in the zenith angle dependence of the event rates are the neutrino yields.
2. In calculation (ii) (1–D approximation) the zenith angle distributions reflect both the geomagnetic effects on the primary cosmic rays, and the zenith angle dependence of the neutrino yields.
3. Finally for the full three–dimensional calculation (i) the zenith angle distributions are determined by a combination of all three effects: the geomagnetic effects, the neutrino yields and the spherical geometry of the source.

Note that in the calculation (iii) where geomagnetic fields are neglected and the primary cosmic ray fluxes are considered as isotropic the neutrino event rates are independent from the detector position, therefore the results (described by the dashed histograms) are equal (within the statistical errors of the montecarlo calculations) for the same event type ($`\mu `$ or $`e`$) and same energy range in each of the three regions that correspond to different magnetic latitude. For the other two calculations the neutrino event rates do depend on the detector position and therefore the results are different in the three magnetic latitude regions.

Let us consider first the zenith angle distributions of the $`\nu `$ fluxes in calculation (iii) (1–D with no geomagnetic effects, represented by the dashed histograms in fig. 611). For $`E_\nu \stackrel{<}{_{}}1`$ GeV, the $`\nu `$ fluxes are essentially isotropic, while with increasing energy the fluxes start to develop an enhancement for horizontal directions. This enhancement is due to the larger neutrino yields in more inclined showers, and becomes more marked when the energy increases as discussed in section 4. Note that the enhancement caused by the neutrino yields is more marked for $`e`$–like events, reflecting the fact that the effect of the shower inclination is especially important for muon decay, that is the source of essentially all $`e`$–like events and of only a fraction (less or equal to one half) of the $`\mu `$ events. Note also that this effect is exactly up–down symmetric, and all dashed histograms remain equal (within statistical errors) for a reflection ($`\mathrm{cos}\theta _z\mathrm{cos}\theta _z`$).

The zenith angle distributions of calculation (ii) (1-D, represented by thin histograms) very clearly display the effect of the geomagnetic cutoffs. This is reflected in two effects: the calculated event rates now depend on the detector position and are highest in the polar region, and lowest in the equatorial region (the $`y`$ axis is an absolute scale in fig. 611), and the up–down symmetry is broken. For example in the calculations for the polar region (fig. 6 and fig. 9) the up–going ($`\mathrm{cos}\theta _z<0`$) rates are smaller that the down–going ones ($`\mathrm{cos}\theta _z>0`$) reflecting the fact that the geomagnetic cutoffs above the detector (in the magnetic polar region) are lower than average, while up–going trajectories are produced everywhere on the Earth. On the contrary in the calculations for the equatorial region (fig. 8 and fig. 10) the up–going rates are larger than the down–going ones that are suppressed because of the high geomagnetic cutoffs in the equatorial region. The geomagnetic effects become negliglible when the neutrino energy becomes larger than a few GeV, reflecting the fact that higher energy neutrinos are produced in the showers of higher energy primary particles that are not affected by the geomagnetic effects.

Finally the 3–D calculation (thick histograms) clearly exhibits the contribution from all three sources of zenith angle dependence. A new and remarkable feature is present for $`E_\nu `$ below few GeV’s, and is the sharp enhancement for the horizontal directions, whose origin we have discussed in section 5.

It should be noticed that the “geometric” enhancement for the horizontal plane in a realistic calculation that includes the geomagnetic effects is stronger than what can be estimated assuming that the primary spectrum is isotropic. This can be easily understood qualitatively noting that the value of the geomagnetic rigidity cutoff (see eq. 2) for a fixed position and azimuth grows monotonically with zenith angle from a minimum value for the vertical direction to a maximum value on the horizontal plane. The enhancement on the horizontal plane, as discussed in the previous section, can be understood as the consequence of the fact that there are more “vertical” showers producing “horizontal” $`\nu `$’s than vice versa. This is true for an isotropic flux, when the source of the vertical/horizontal asymmetry is due to the fact that the zenith angle distribution of the showers is $`\mathrm{cos}\theta _0`$ (with $`\theta _0`$); the asymmetry between vertical and horizontal showers becomes more important when the geomagnetic effects are included, since in this case the horizontal primary flux is also suppressed by an higher rigidity cutoff. For a realistic calculation the size of the “geometric” horizontal enhancement also depends on the position of the detector, reflecting the effect discussed above, and the fact that the energy spectrum of the $`\nu `$ flux observed in different location changes being softer (harder) at high (low) magnetic latitude, and the angle $`\theta _{0\nu }`$ is correlated with the neutrino energy.

The “spherical geometry” enhancement can be easily distinguished from the “neutrino yield” enhancement since it is most important at low neutrino energy reflecting the higher average neutrino–primary angle, and vanishes as the angle $`\theta _{0\nu }`$ shrinks with increasing $`E_\nu `$. Conversely the “neutrino yield” enhancement becomes more important as the neutrino energy increases.

## 7 Conclusions

In this work we have discussed the different sources that together with neutrino oscillations determine the zenith angle distributions of the atmospheric neutrino fluxes. The source volume of the atmospheric neutrinos has the shape of a spherical shell, and a simple consequence of this geometry is the existence of an enhancement of the flux intensity for the horizontal directions. This enhancement is stronger when the neutrinos are emitted with a large average angle with respect to the primary particle direction, and therefore it is most important at low energy, and disappears for $`E_\nu \stackrel{>}{_{}}3`$–10 GeV when the neutrino emission is approximtely collinear with the shower axis.

The detection of this enhancement is not easy, since the angular resolution for sub–GeV neutrinos is poor. This is true for a detector such as Super–Kamiokande that can measure only the charged lepton produced in a quasi–elastic scattering (such as $`\nu _\mu n\mu ^{}p`$), since the angle between the neutrino and the final state charged lepton is large, but it is also true  for higher resolution detectors (such as Soudan–2 or Icarus) that are capable of detecting the recoil protons. Even in this case the angular resolution will be limited by the fact that the spectator nucleons (that are undetected) can absorb a non negligible 3–momentum, and therefore the vector sum of the recoil nucleon and charged lepton momenta does not correspond exactly to the initial neutrino momentum.

This difficulty in the detection of the enhancement does not imply that the effect is negligible when the atmospheric neutrino data is analysed to extract information about the oscillation parameters. Several effects can be significant. The most important one is a significant modification of the pathlength distributions of the neutrino events. These distributions are controlled by the $`\nu `$ zenith angle distributions that can very significantly distorted by the geometrical effects discussed here. Note also that the height distribution of the horizontal neutrinos for a fixed value of the zenith angle, is also modified with respect to a 1–D calculation. For example in a 3–D calculation most of the sub–GeV horizontal neutrinos are actually produced in less inclined (and deeper developing showers) after emission at large angles with respect to the shower axis, and are therefore produced at lower altitude with respect to the 1–D case. This effect is second order with respect to the distortion of the zenith angle distribution but must also be taken into account. As an example, even if the detector resolution is so poor that the shape of the neutrino angular distribution is not measurable, in the presence of $`\nu _\mu \nu _\tau `$ oscillations the pathlength distribution determines the average (zenith angle integrated) suppression of the $`\mu `$–like event rate for a given set of oscillation parameters.

The geometrical effects can also be a non negligible correction for the calculation of the $`\mu /e`$ ratio. The enhancement for the horizontal directions depend in fact on the average angle $`\theta _{0\nu }`$ between neutrino and primary particle and on the average altitude of the neutrino production. Both these quantities are slightly different for neutrinos produced directly in meson decay or in the chain decay of a muon, and are therefore also slightly different for electron and muon (anti)–neutrinos.

Also the absolute observable rates for sub–GeV neutrino events depend on these effects. The absolute rate is not important in the determination of the oscillation parameters, but the consistency of data and prediction also for the absolute rate is certainly desirable.

The improved calculation of the atmospheric $`\nu `$ fluxes is now an active field of research, strongly stimulated by the great importance of the discovery of the existence of flavor transitions for these neutrinos. New improved calculations will be important in the detailed interpretation of the existing and future data on atmospheric neutrinos, and in the unbiased determinations of the $`\nu `$ oscillation parameters (or other parameters describing the new physics of the flavor transitions). New measurements of the primary cosmic rays have been recently performed , and have to be included in a new fit for the primary cosmic rays, and at the same a significant effort has been dedicated to the developoment of improved modeling of hadronic interactions in cosmic ray showers. Because of these developments in the field we postpone to a future work a quantitative discussion of a new estimate of the oscillation parameters. The newly recognizied geometric effects that we have discussed in this work will be important in these future calculations.

Acknowledgments Special thanks to Giuseppe Battistoni for several very useful discussions and to Takaaki Kajita for encouragement. I’m grateful to Tom Gaisser, Teresa Montaruli and Todor Stanev for reading the manuscript and giving suggestions for improvements. I also acknowledge useful discussions with Ed Kearns, Todd Haines and M. Honda.


format: text

text:
# Steady-state selection in driven diffusive systems with open boundaries

## Abstract

We investigate the stationary states of one-dimensional driven diffusive systems, coupled to boundary reservoirs with fixed particle densities. We argue that the generic phase diagram is governed by an extremal principle for the macroscopic current irrespective of the local dynamics. In particular, we predict a minimal current phase for systems with local minimum in the current–density relation. This phase is explained by a dynamical phenomenon, the branching and coalescence of shocks, Monte-Carlo simulations confirm the theoretical scenario.

\]

A recurrent problem in the investigation of many-body systems far from equilibrium is posed by the coupling of a driven particle system with locally conserved particle number to external reservoirs with which the system can exchange particles at its boundaries. In the presence of a driving force a particle current will be maintained and hence the system will always remain in an nonequilibrium stationary state characterized by some bulk density and the corresponding particle current. While for periodic boundaries the density is a fixed quantity, the experimentally more relevant scenario of open boundaries naturally leads to the question of steady-state selection, i.e. the question which stationary bulk density the system will assume as a function of the boundary densities . In the topologically simplest case of quasi one-dimensional systems this is of importance for the understanding of many-body systems in which the dynamic degrees of freedom reduce to effectively one dimension as e.g. in traffic flow , kinetics of protein synthesis , or diffusion in narrow channels .

Within a phenomenological approach this problem was first addressed in general terms by Krug  who postulated a maximal-current principle for the specific case where the density $`\rho ^+`$ at the right boundary to which particles are driven is kept at zero. The exact solution of the totally asymmetric simple exclusion process (TASEP)  for arbitrary left and right boundary densities $`\rho ^{},\rho ^+`$ confirms and extends the results by Krug. The complete phase diagram comprises boundary-induced phase transitions of first order between a low-density phase and a high density phase, and a second-order transition from these phases to a maximal current phase . Analysis of the density profile  provides insight into the dynamical mechanisms that lead to these phase transitions and shows that the phase diagram is generic for systems with a single maximum in the current-density relation . Experimental evidence for the first-order transition is found in the process of biopolymerization for which the TASEP with open boundaries was originally invented as a simple model , and more directly also in recent measurements of highway traffic close to on-ramps . Renormalization group studies  indicate universality of the second order phase transition.

In this Letter we develop a dynamical approach to generic driven one-component systems with several maxima in the current-density relation, and we show that a novel phase of rather unexpected nature appears: For a certain range of boundary densities the steady state carries the minimal current between two maxima even though both boundary densities support a higher current. We shall refer to this phase as minimal current phase. More generally, we claim that the current always obeys the extremal principle

$`j`$ $`=`$ $`\underset{\rho [\rho ^+,\rho ^{}]}{\mathrm{max}}j(\rho )\text{ for }\rho ^{}>\rho ^+`$ (1)
$`j`$ $`=`$ $`\underset{\rho [\rho ^{},\rho ^+]}{\mathrm{min}}j(\rho )\text{ for }\rho ^{}<\rho ^+.`$ (2)
To understand the origin of this extremal principle we first note that in the absence of detailed balance stationary behavior cannot be understood in terms of a free energy, but has to be derived from the system dynamics. For definiteness consider a driven lattice gas with hard-core repulsion. At $`\rho =1`$ no hopping can take place and hence the current vanishes. Two local maxima can arise as the result of sufficiently strong repulsion between nearest neighbor particles as opposed to the pure on-site repulsion of the usual TASEP which leads to a single maximum (Fig. 1).

At first sight one might not expect such a little change in the interaction radius of the particles to affect the phase diagram. However, the extremal principle (1), (2) predicts that the full phase diagram (Fig. 2) generically consists of seven distinct phases, including two maximal current phases with bulk densities corresponding to the respective maxima of the current  and the minimal current phase in a regime defined by

$$j(\rho ^+),j(\rho ^{})>j(\rho _{min});\rho ^{}<\rho _{min}<\rho ^+.$$
(3)
Here the system organizes itself into a state with bulk density $`\rho _{bulk}`$ corresponding to the local minimum of the current. As in the maximal current phases no fine-tuning of the boundary densities is required.

In order to understand the basic mechanisms which determine the steady state selection we first follow Ref.  and consider the collective velocity

$$v_c=j^{}(\rho )$$
(4)
which is the velocity of the center of mass of a local perturbation in a homogeneous, stationary background of density $`\rho `$ (Fig. 3(a)). A second quantity of interest is the shock velocity

$$v_s=\frac{j_2j_1}{\rho _2\rho _1}$$
(5)
of a ‘domain wall’ between two stationary regions of densities $`\rho _{1,2}`$ (Fig. 3(b)). Notice that both velocities may change sign, depending on the densities.

These velocities are sufficient to understand the phase diagram of systems with a single maximum in the current . Further to these notions we need to introduce here the concepts of coalescence and branching of shocks. A single large shock (with a large density difference $`\rho _2\rho _1`$) may be understood as being composed of subsequent smaller shocks with narrow plateaux at each level of density (Fig. 4). In the usual asymmetric exclusion process these shocks travel with different relative speeds such that they coalesce and form a ‘bound state’ equivalent to a single shock . In the present situation, however, the minimum in the current-density relation leads to more complicated dynamics. A closer investigation of Eq. (5) shows that depending on $`\rho _{1,2}`$ a single shock may branch into two distinct shocks, moving away from each other, discussed below.

With these observations the dynamical origin of the phase transition lines can be understood by considering the time evolution of judiciously chosen initial states. Because of ergodicity, the steady state does not depend on the initial conditions and a specific choice involves no loss of generality. We turn our attention to a line $`\rho ^+=c`$ with $`\rho _{min}<c<\rho _2^{}`$ in the phase diagram which crosses the minimal current phase. Along this line it is convenient to consider an initial configuration with a shock with densities $`\rho ^{}`$ and $`\rho ^+`$ on the left and on the right respectively, which is composed of many narrow subsequent shocks at various levels of intermediate densities (Fig. 4).

(i) Let us start with the point with equal boundary densities in which case the system evolves into a steady state with the same bulk density $`\rho _{bulk}=\rho ^{}=\rho ^+`$.

(ii) Lower $`\rho ^{}`$ slightly below $`\rho ^+`$ with just a single shock separating both regions. According to (5) the shock travels with speed $`v_s=(j^+j^{})/(\rho ^+\rho ^{})>0`$ to the right, making the bulk density equal to $`\rho ^{}`$. At the same time, small disturbances will according to (4) also drift to the right, as $`v_c=j^{}(\rho ^{})>0`$ in this region, thus stabilizing the single shock.

(iii) Now, lower $`\rho ^{}`$ below $`\rho _{min}`$. While the shock velocity $`v_s`$ is still positive, so that one expects the shock to move to the right, the collective velocity $`v_c=j^{}(\rho ^{})<0`$ indicates that disturbances will spread to the left. This discrepancy marks the failure of a single shock scenario. In order to resolve it, consider for simplicity $`\rho ^{}=\rho _1^{}`$ and return to the picture with many subsequent shocks at each density level between $`\rho ^{}`$ and $`\rho ^+`$ (Fig.4). Eq. (5) shows that small shocks below $`\rho _{min}`$ will move to the left, while those above $`\rho _{min}`$ will move to the right***Actually, the small shock between the levels $`\rho _k^{},\rho _k^+`$ will be stable only if $`v_c(\rho _k^+)<v_s<v_c(\rho _k^{})`$, see Eq.(4,5). Some of small shocks satisfy the above condition and some do not. This brings about additional structure to the resulting shock, which will be discussed elsewhere . The leftmost of the left-moving shocks will merge in a single one, and so will the rightmost of the right-moving shocks. The result are two single shocks and moving different directions and thus expanding the region with the density $`\rho _{bulk}=\rho _{min}`$. The system is in the minimal current phase. This picture is well supported by the Monte-Carlo simulations shown in Fig. 5, demonstrating the branching of a single shock into two distinct shocks moving in opposite directions. The minimal current phase will persist for any left boundary density in the range $`\rho ^{}[\stackrel{~}{\rho }_1,\rho _{min}]`$. Notice that the change of bulk density is continuous across the point $`\rho ^{}=\rho _{min}`$ to the minimal current phase, so the transition is of the second order.

(iv) As we lower $`\rho ^{}`$ below $`\stackrel{~}{\rho }_1`$, the shock velocity $`v_s=(j_{min}j(\rho ^{}))/(\rho _{min}\rho ^{})>0`$ becomes positive. The shock is moving to the right, leading to a low density phase with bulk density $`\rho _{bulk}=\rho ^{}`$ which drops discontinuously from $`\rho _{bulk}=\rho _{min}`$ at $`\rho ^{}=\stackrel{~}{\rho }_1+0`$ to $`\rho _{bulk}=\rho ^{}`$ at $`\rho ^{}=\stackrel{~}{\rho }_10`$. The system undergoes a first order phase transition. On the transition line the shock performs an unbiased random walk, separating coexisting regions of densities $`\rho _{min}`$ and $`\rho ^{}`$ respectively.

(v) Let us start again from $`\rho ^{}=\rho ^+`$ and now increase $`\rho ^{}`$. Until one reaches $`\rho ^{}=\rho _2^{}`$, the collective velocity $`v_c=j^{}(\rho ^{})>0`$ is positive, leading to $`\rho _{bulk}=\rho ^{}`$.

(vi) As soon as $`\rho ^{}`$ crosses the point $`\rho ^{}=\rho _2^{}`$, the sign of collective velocity $`v_c`$ changes and the overfeeding effect  occurs: a perturbation from the left does not spread into the bulk  and therefore further increase of the left boundary density does not increase the bulk density. The system enters the maximal current phase II through a second-order transition.

Using analogous arguments one constructs the complete phase diagram (Fig. 2) and obtains the extremal principle (1, 2). The velocities (4), (5) which determine the phase transition lines follow from the current-density relation. This behavior can be checked with Monte Carlo simulations. A model with two maxima of the current is a TASEP with next nearest neighbor interaction defined by the bulk hopping rates (see A26,A27 in )

$`\mathrm{0\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}0}`$ $``$ $`\mathrm{0\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0}\text{ with rate }1+\delta `$ (6)
$`\mathrm{1\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}0}`$ $``$ $`\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}0}\text{ with rate }1+ϵ`$ (7)
$`\mathrm{0\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}1}`$ $``$ $`\mathrm{0\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}1}\text{ with rate }1ϵ`$ (8)
$`\mathrm{1\hspace{0.33em}1\hspace{0.33em}0\hspace{0.33em}1}`$ $``$ $`\mathrm{1\hspace{0.33em}0\hspace{0.33em}1\hspace{0.33em}1}\text{ with rate }1\delta `$ (9)
with $`|ϵ|<1;|\delta |<1`$. For sufficiently strong repulsive interaction $`1ϵ1`$ the current at half-filling is strongly suppressed, which brings about two maxima structure in the current-density relation (Fig. 1). The limit $`ϵ=1`$ leads to $`j_{min}=0`$. For negative $`ϵ`$ (attractive interaction) the current-density relation always has a single maximum. The other parameter, $`\delta `$, is responsible for the particle-hole symmetry. $`\delta =0`$ corresponds to symmetrical graph $`j(\rho )=j(1\rho )`$. $`\delta 0`$ breaks the particle-hole symmetry in favor of larger particles current ( $`\delta >0`$) or larger vacancies current ( $`\delta <0`$). In particular, $`\delta 0`$ is responsible for different heights of two maxima on Fig. 1.

The injection at the left boundary site 1 and extraction of particles at the right boundary site $`L`$ is chosen to correspond to coupling to boundary reservoirs with densities $`\rho ^\pm `$ respectively. Along the line $`\rho ^+=\rho ^{}`$ the stationary distribution is then exactly given by the equilibrium distribution of an 1-D Ising model with boundary fields and the bulk field such that the density profile is constant with density $`\rho =\rho ^+=\rho ^{}`$ . The current $`j=(1+\delta )\mathrm{\hspace{0.17em}0100}+(1+ϵ)\mathrm{\hspace{0.17em}1100}+(1ϵ)\mathrm{\hspace{0.17em}0101}+(1\delta )\mathrm{\hspace{0.17em}1101}`$ as a function of the density can be calculated exactly using standard transfer matrix techniques. The exact graph is shown in Fig. 1 for specific values of the hopping rates.

We performed Monte Carlo simulations for systems of sizes $`L`$ from 100 to 1000. Densities and currents were averaged over at least 50$`L`$ rounds, and averaged over 7 different histories. We consider the bulk density $`\rho _{bulk}`$ as order parameter. In analogy to equilibrium phase transitions a singularity (jump discontinuity) in the first derivatives of $`\rho _{bulk}`$ signals a first (second) order phase transition. In the finite system we determined the location of the first (second) order transition by the appearance of peak (jump) in the first derivatives of the bulk density $`\rho _{bulk}(\rho ^+,\rho ^{})`$ with respect to $`\rho ^+`$ and $`\rho ^{}`$. As initial state we chose either the empty or the completely filled lattice, whichever gave the faster convervence. Despite finite size effects, precise analysis of which requires further investigations, the overall agreement of the simulated phase diagram with the predicted one is very good already for $`L=150`$ (Fig.2) We have chosen the set of parameters $`ϵ`$ and $`\delta `$, which shows distinct minimum between two maxima. For the other choices of $`ϵ,\delta `$, the phase diagram can be obtained in the same way, e.g. by current-density relation and extremal principles Eqs.(1,2)..

We conclude that the interplay of density fluctuations and shock diffusion, coalescence and branching resp. as described above represents the basic mechanisms which determine the steady state selection (1,2) of driven diffusive systems with a nonlinear current–density relation. A surprising phenomenon is the occurrence of the self-organized minimal current phase. Since little reference is made to the precise nature of the dynamics we argue that the phase diagram is generic and hence knowledge of the macroscopic current-density relation of a given physical system is sufficient to calculate the exact nonequilibrium phase transition lines which determine the density of the bulk stationary state. For systems with more than two maxima in the current the interplay of more than two shocks has to be considered.

A fundamentally different scenario seems likely in systems with long-range bulk correlations, which, however, is untypical in one dimension. Perhaps the most interesting open question which one might address in a similar fashion is the question of steady-state selection in systems with two or more local conservation laws. This could provide insight into the mechanisms responsible for boundary-induced spontaneous symmetry breaking  in 1-D driven diffusive systems with open boundaries.

Acknowledgements This work was partially supported by DAAD and CAPES within the PROBRAL programme. G.M.S. thanks H. Guiol for useful discussions.


text:
# Ferroelectricity and structure of BaTiO3 grown on YBa2Cu3O7-δ thin films

## I Introduction

In recent years, there has been much interest in the preparation of superconducting field effect transistors (SuFETs) with a (Ba,Sr)TiO3 thin film as gate insulator and YBa2Cu3O7-δ (YBCO) as source-drain channel . More recently the development of high frequency superconducting devices increased the interest of the community on YBCO/ferroelectric heterostructures . From these studies it became clear that it is possible to grow high quality perovskite ferroelectric thin films on top of YBCO electrodes. This motivated us to investigate the room temperature ferroelectric and structural properties of YBCO/BaTiO3/Au heterostructures.

The objective of this work is to analyse the domain wall contribution to the dielectric properties of high quality BaTiO3 (BTO) thin films. So far, similar studies have only been done on Pt/PbZrxTi1-xO3 (PZT)/Pt heterostructures .

## II Experiment

The 100 nm thick YBCO layers and the BTO films were sputtered in-situ in pure O2 atmosphere. The top Au electrodes were thermally evaporated through a shadow mask. Before deposition, the YBCO target was pre-sputtered for 1 hour and the BTO target for 30 minutes. After this time the dc-voltage of the plasma did not decrease anymore showing that the targets reached a stable state. For both processes the substrate target distance was about 20 mm. The heterostructures have been prepared on MgO substrates because of the better High Frequency properties of this material in comparison to SrTiO3. The YBCO ground electrode was prepared following the method already described in details in Ref. . The difference in the present work is that the post annealing was performed after the deposition of the BTO layer. The stoichiometric BTO targets were prepared at the University of Mainz. The BTO films were deposited by rf-sputtering in on axis geometry without a magnetron in O2 pressures between 0.5 and 1.5 mbar. This was needed in order to prevent re-sputtering. During deposition, the rf power was regulated at 130 W resulting in a typical bias voltage of 30 V. The best films were deposited at 700 C. At the end of the process, the bi-layers were cooled down to room temperature during 1 hour and annealed in 1 bar oxygen pressure. The electrode contact size varied between 0.5 mm2 and 1 mm2. X-ray diffraction patterns in Bragg-Brentano and four-circle geometry were obtained using a Philips X’Pert-MPD and a STOE STADI4 diffractometer, respectively. The capacitance and the hysteresis loops were measured with a HP4284A LCR-meter and a self made Sawyer-Tower circuit, respectively.

## III Results and Discussion

### A X-ray diffraction

The $`\theta `$/2$`\theta `$-scan shown in Fig. 1 reveals only $`c`$-axis oriented growth of the BTO films. However as shown by detailed TEM studies there may still exist small amounts of $`a`$-axis oriented domains embedded on a $`c`$\- axis matrix even if by x-ray diffraction no splitting caused by $`a`$-and $`c`$-axis oriented domains can be found . The $`c`$-axis length of the BTO thin films in our work sputtered at 0.7 mbar is determined to be 3% larger in comparison with the bulk material and coincides with the value of the work of Abe et al. . On the other hand, by using a deposition pressure of 1.5 mbar we obtain the same $`c`$-axis length like the bulk material. Hence the lattice expansion at a lower pressure may result from atomic defects due to the high kinetic energy of the incident species. This explanation is in accordance with our observation that the lattice expansion occurs only for BTO but not for SrTiO3 thin films. As barium ions are much heavier than strontium ions they can induce more serious damage upon impact on the previously deposited film. We expect therefore the lattice parameters to be also dependent on the deposition technique, as already reported by Srikant et al. .

In order to characterise the in-plane orientation of the YBCO/BTO heterostructures we performed a $`\varphi `$-scan of the YBCO (103) reflex (see Fig. 2) and a two-dimensional scan of the plane $`\mathrm{}=1`$ for BTO (see Fig. 3). For both films there are 8 peaks, indicating a fourfold symmetry and an in-plane alignment of YBCO’s and BTO’s crystallographic axes parallel and rotated 45 relative to the MgO lattice. The rotated alignment can be explained by the fact that 2 lattice constants of MgO equals 3 diagonal lattice parameters of YBCO. Since BTO grows on YBCO epitaxial it reproduces this 45 rotation.

A fourfold symmetry was obtained by the use of SrTiO3 substrates. The orientation matrix gives for our films the in-plane lattice constants $`a=b=(0.3994\pm 0,0002)`$ nm for BTO. This value is in excellent agreement with the result on bulk materials. For YBCO we measured $`a=0.3847`$ nm, $`b=0.3851`$ nm. These values coincide well with the average between $`a`$ and $`b`$ in bulk materials. This is due to the usual mixing of these two orientations observed in thin films done on standard substrates. The results above show, that the $`c`$\- axis distortion of BTO cannot be explained by the requirement to conserve the unit cell volume, according to the Poisson ratio. This indicates that the high energy of barium at the sputtering process might be indeed responsible for the observed lattice distortion.

The mosaic spreads of the different layers were determined by measuring the rocking curves of BTO (002) and YBCO (005). For BTO we obtain a FWHM (full width at half maximum) of 0.38 and for YBCO 0.32, indicating a good epitaxial growth of the heterostructures. The small FWHM obtained in our BTO films can be explained by a reduced lattice mismatch between YBCO and BTO. For comparison, Wills et al. got for $`a`$-axis oriented BTO on LaAlO3 substrates a FWHM of 0.6 .

### B Ferroelectric hysteresis loops

The non-remanent component of polarisation is caused by non-ferroelectric ionic and electronic polarisability as well as field induced reversible ferroelectric domain wall motion. The remanent component is induced by switching ferroelectric domains. Using a Sawyer-Tower circuit, we obtain a remanent polarisation of $`P_\mathrm{r}=1.25\mu `$C/cm2 and a coercive field of $`E_\mathrm{c}=30`$ kV/cm (Fig. 4). The remanent polarisation is consistent with other results obained by rf-sputtering . But Yoneda et al. even observed for 30 monolayers thick BTO films a ferroelectric hysteresis up to 600 C . Measuring the dielectric constant of our 300 nm thick BTO films, from room temperature to 200 C, did not indicate any sign of a ferroelectric phase transition, as well. This result is rather surprising since bulk BTO is known to have a Curie temperature of about 120 C . This increase in the curie temperature could be due to internal stresses induced by the substrate.

### C Contribution of the domain walls to the permittivity

The major research in ferroelectric thin films emphasises the nucleation, growth and switching of domains because these features are very important for nonvolatile memory applications. In contrast, contributions to the dielectric or the piezoelectric constant at sub-switching fields due to domain wall displacements and vibrations were investigated mainly for single crystals or ceramics. We follow here the procedure of Taylor et al. by extending the models for the interaction and pinning of domain walls with randomly distributed pinning centres in magnetic materials to the dielectric displacement in ferroelectric thin films . Below we will address the dielectric response of a BTO thin film both in terms of the weak-field (Rayleigh law) and the logarithmic frequency dependence. The Rayleigh law is valid for low-field conditions, where no nucleation of domains occur and the average structure of the domain walls drifts within the sample as the field is cycled . In order to achieve the linear Rayleigh behaviour for the dielectric constant we cycled the sample about 107 times prior to the measurement with an amplitude below the coercive voltage. Figure 5 shows the linear dependence of the capacitance versus ac voltage amplitude $`V_0`$. From the capacitance at 100 Hz, and from the geometry of our electrodes we deduce a dielectric constant $`ϵ=300`$. In accordance to the Rayleigh law the weak-field dependence can be described by the formula $`C=C_{\mathrm{init}}+aV_0`$. $`C_{\mathrm{init}}`$ represents the intrinsic lattice and reversible domain wall movement. The second term $`aV_0`$ with the Rayleigh coefficient $`a`$ is responsible for the weak-field hysteresis due to the pinning and depinning of domain walls. According to Néel it is caused by the irreversible displacement of the walls between two meta-stable states when the applied field is large enough to overcome the barrier. In Fig. 6 we show the logarithmic dependence of the capacitance on frequency. As for magnetic systems this logarithmic behaviour is generally attributed to domain wall motion across randomly distributed pinning centres . In order to determine whether the reversible $`C_{\mathrm{init}}`$ or the irreversible Rayleigh coefficient $`a`$ is responsible for this behaviour, we extracted them for all frequencies from the ac field dependence of the capacitance. The result indicating the logarithmic frequency dependence of both parameters is presented in Fig 7. Using a Sawyer-Tower circuit we finally measured the non-saturating hysteresis loop (Rayleigh loop) of our samples at sub-switching fields (see Fig. 8). The full line representing the calculated hysteresis curve is based on the LCR-meter measurement by using the formula $`P=(ϵ_{init}+\alpha E_0)E\pm \alpha /2(E_0^2E^2)`$, where $`+`$ stands for the decreasing and $``$ for the increasing field. $`ϵ_{\mathrm{init}}`$ and $`\alpha `$ are derived from $`C_{\mathrm{init}}`$ and $`a`$ respectively. $`E_0`$ designates the maximum field reached during the cycle. While the slopes of the loops coincide, the area of the calculated hysteresis curve is underestimated in comparison to the measured one. This fact suggests that the observed hysteresis is not only due to the Rayleigh-type irreversible contribution of the dielectric constant, but by the losses having their origin in the series resistance and parallel conductivity of the capacitor. This explanation is confirmed by the fact that the calculated curve is only based on the imaginary and not on the real part of the impedance.

So far there is no universal behaviour for the frequency response of capacitors based on insulating perovskites. While Taylor et al. observed, for Pt/PZT/Pt capacitors, in the range from 10-1 Hz to 104 Hz an exponential law between capacitance and frequency , Zafar et al. used for Pt/(Ba,Sr)TiO3/Pt between 102 Hz and 106 Hz a power law to describe their data .

While the logarithmic frequency dependence is in accordance to the contribution of domain walls, the power law dependence is known as the Curie-von Schweidler law. Obviously further studies are necessary to determine the correct relationship between the dielectric constant and the frequency. Maybe our data and the results of Taylor et al, fit better to the logarithmic frequency dependence because the films exhibit a remanent polarisation and thus a higher domain wall contribution in comparison to non-ferroelectric (Ba,Sr)TiO3 films. The relative irreversible domain wall contribution to the total capacitance can be obtained by the ratio $`aV/C`$. If we apply an sinusoidal electrical field with a frequency of 100 Hz and an amplitude of $`E_0=5`$ MV/cm we obtain a ratio of 3%. This value is much smaller than the 21% calculated by Taylor et al. for 1.3 $`\mu `$m thick PZT thin films. Since the irreversible part of the dielectric constant determines the remanent polarisation it is reasonable for BTO to obtain a lower contribution of irreversible domain wall movement than for PZT because the remanent polarisation of PZT thin films is typically in the range of 40 $`\mu `$C/cm2. This is more than one magnitude larger in comparison with our BTO thin films, revealing a remanent polarisation of 1.25 $`\mu `$C/cm2. Damjanovic et al. determined the Rayleigh parameters describing the piezoelectric response of PZT and BTO ceramics by applying an external ac pressure to the samples . The Rayleigh law for the piezoelectric behaviour can be described by the formula $`d=d_{\mathrm{init}}+\beta X_0`$ with the piezoelectric constant $`d`$, the amplitude $`X_0`$ of the applied pressure and the Rayleigh coefficient $`\beta `$. In contrast to our results for thin films they observed, for the coarse grain (average grain size 27 $`\mu `$m) BTO ceramic, a higher relative irreversible contribution $`\beta X/d`$ than for soft PZT. Thus the potential barrier for a domain wall displacement is not dominated by the material, but by the grain size and the type of domain walls because 180 domain walls influence the dielectric constant but not the piezoelectric coefficient . Demartin concluded for fine grain samples (0.7 $`\mu `$m grain size) of BTO that domain walls are clamped to a considerable degree by the high internal pressure  and hence they contribute less to the piezoelectric coefficient. In contrast to coarse grain ceramics (26 $`\mu `$m grain size), the domain walls created in fine grain ceramics are not sufficient to relieve the stress produced by the phase transition, resulting in a high internal stress and a weak activity of the domain walls. The decrease in dielectric constant with decreasing grain size (from 50 nm to 12 nm) arises also from the absence of 180 domain boundaries in smaller grains . On the other hand, it is well known that polycrystalline BTO with a grain size of approximately 1 $`\mu `$m possesses a much higher permittivity at room temperature than a single crystal BTO due to the higher contribution of domain walls . The much larger film thickness of 1.3 $`\mu `$m chosen by Taylor and Damjanovic  yields an enhanced activity of the domain walls and a higher Rayleigh coefficient in comparison to our 300 nm thick BTO films.

The grain size and the stresses due to the phase transition determine the dielectric constant in ceramics. But ferroelectric thin films are additionally influenced by the existence of interfaces and additional stresses, caused by the lattice mismatch between film and substrate.

## IV Conclusion

The small lattice mismatch and chemical compatibility between YBCO and BTO allow a good crystal structure for heterostructures of both materials. X-ray diffraction patterns display a high crystalline order of both YBCO and BTO whereas their in-plane orientation is influenced by the choice of MgO as substrate material. The sub-switching ac electric field dependence of the permittivity may be described in terms of the Rayleigh law. At a field of 5 MV/m the total dielectric constant contains 3% contribution of irreversible domain wall movement. In addition to $`C_{init}`$ the irreversible Rayleigh parameter representing irreversible domain wall movement shows also a linear logarithmic frequency dependence of permittivity. This behaviour can be explained by pinning of domain walls. In comparison to conventional Pt as base electrode for ferroelectric memories, YBCO offers the advantages of better fatigue behaviour.

The authors gratefully acknowledge support by the Materialwissenschaftlichen Forschungszentrum (MWFZ), Mainz, the Training and Mobility of young researchers (TMR) program of the European Union under grant number ERBFM-BICT972217 and the German Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie (BMBF) under project number 13N6916.