text:
# Random-phase approximation study of collective excitations in the Bose-Fermi mixed condensate of alkali-metal gases

## I introduction

Collective excitation is one of the most prominent phenomena in quatum many-body systems such as liquid helium, electron gas, nucleus, etc. In the recently developed Bose-Einstein condensates of trapped atomic gases, collective oscillations were the first of the dynamical phenomena discovered. Collective oscillations are characterized by various quantum numbers related to, e.g., the shape of oscillation, or internal degrees of freedom of the constituents such as spin, isospin, etc. The oscillation frequency, the damping width, etc. depend on the interparticle interaction of the constituents, and thus provide a clue to unravel the dynamical correlation of the many-body system.

Study of the properties of trapped neutral atoms has been extended to fermi systems where the occurrence of the Fermi degeneracy was observed, and also to the mixture of bose and fermi particles. The latter system with condensed bosons and degenerate fermions was recently realized experimentally . This system is one of the typical example in which particles obeying different statistics are intermingled. Theoretical studies of the bose-fermi mixed system of cold atomic gases have been done for static properties, for the phase diagram and phase separation, for the stability of the system and for collective excitations.

In Ref.  sum rule approach has been applied for collective excitations in the bose-fermi mixed system. Average excitation energies for the states with multipoles $`L=0,1,2`$ were calculated for both in-phase and out-of-phase modes of the bose and the fermi particles, and the dependence on the bose-fermi interaction strength was studied. The sum rule approach is a powerful technique for collective states in quantum many-body systems, and has been successfully applied to the excitation of bose condensed systems, see Refs. . It does not give, however, direct information on the eigenstates of the system, but rather an average behavior of the strength distribution for the adopted multipole excitation operator. For a more detailed investigation on the dynamical properties of the system, one would need a study of individual eigenstates.

In the present paper we extend the previous study of collective excitations in the bose-fermi mixed system. We calculate full particle-hole type excitations by diagonalization of an RPA type matrix. The motive of this calculation is threefold: first of all, the calculation allows us to study the excitation spectrum and the strength distribution in contrast to the sum rule method which focuses on the strength weighted average of the prescribed multipole operators. We study, for instance, the degree of collectivity of the excitations depending on the strength of boson-fermion interaction, and estimate the damping of the collective excitation albeit within the space of particle-hole excitations. Information on the wave function allows us to calculate observables such as the dynamical structure factor. The latter for bose-condensed systems is now becoming available experimentally by two-photon spectroscopy and is being studied theoretically. Secondly, we compare the results with those obtained from the sum rule approach. This provides a check on the approximation adopted in the calculation such as the model space truncation. The comparison also allows one to examine the structure of the low- and high-lying collective modes which was speculated in  through the mixing angle $`\theta `$ of multipole operators. Third, we can predict the time-dependent behavior of the system for a given external perturbation. This process is actually the one that has been employed in the previous experimental study of collective excitations in BEC. An RPA study of the bose-fermi mixed system has recently been done in , where the response for an external multipole field is formulated in the form of an integral equation. In the present calculation we approach the problem by diagonalizing the RPA matrix, by incorporating the discrete nature of the excitation in an isolated trapped system, and investigate, in particular, the properties of the strength distribution for a combination of bose and fermi operators and for various values of the bose-fermi interaction strengths.

The content of the paper is as follows: in the next section we derive the RPA equation for the bose-fermi mixed system using the equation of motion for particle-hole type excitation operators. The single particle (hole) states are obtained in the mean-field calculation, i.e., by solving the coupled Gross-Pitaevskii and Thomas-Fermi equations. In Section III we first briefly describe the parameters and the numerical procedure employed in the present calculation. We then turn to the detailed studies of the results obtained, including ground state density, single particle states, and strength distribution for each multipole. Comparison with the sum rule calculation is given. Transition densities and dynamical structure factors for some of the collective excitations are presented. We finally consider the time development of the system after an external multipole impulse on the system. Last section is devoted to summary and conclusion.

## II formulation

We consider a dilute spin-polarized bose-fermi mixed system trapped in a spherically symmetric harmonic oscillator potential at $`T=0`$. The system is described by the Hamiltonian

$$\widehat{H}=\widehat{H}_0+\widehat{V}_b+\widehat{V}_{bf}$$
(1)
with

$`\widehat{H}_0`$ $`=`$ $`{\displaystyle d^3r\widehat{\psi }^{}(\stackrel{}{r})\left[\frac{\mathrm{}^2}{2m}\stackrel{}{}^2+\frac{1}{2}m\omega _0^2r^2\right]\widehat{\psi }(\stackrel{}{r})}+{\displaystyle d^3r\widehat{\varphi }^{}(\stackrel{}{r})\left[\frac{\mathrm{}^2}{2m}\stackrel{}{}^2+\frac{1}{2}m\omega _0^2r^2\right]\widehat{\varphi }(\stackrel{}{r})},`$ (2)
$`\widehat{V}_b`$ $`=`$ $`{\displaystyle \frac{1}{2}}g{\displaystyle d^3rd^3r^{}\widehat{\varphi }^{}(\stackrel{}{r})\widehat{\varphi }^{}(\stackrel{}{r}^{})\delta ^{(3)}(\stackrel{}{r}\stackrel{}{r}^{})\widehat{\varphi }(\stackrel{}{r}^{})\widehat{\varphi }(\stackrel{}{r})},`$ (3)
$`\widehat{V}_{bf}`$ $`=`$ $`h{\displaystyle d^3rd^3r^{}\widehat{\varphi }^{}(\stackrel{}{r})\widehat{\psi }^{}(\stackrel{}{r}^{})\delta ^{(3)}(\stackrel{}{r}\stackrel{}{r}^{})\widehat{\psi }(\stackrel{}{r}^{})\widehat{\varphi }(\stackrel{}{r})\text{,}}`$ (4)
where $`\widehat{\varphi }`$ and $`\widehat{\psi }`$ are the boson and the fermion field operators. We take, for simplicity, the mass $`m`$ and the trapping frequency $`\omega _0`$ to be the same for the boson and the fermion. The boson-boson and the boson-fermion interaction strengths for the pseudopotentials, $`g`$ and $`h`$, are related to the $`s`$-wave scattering lengths $`a_{bb}`$ and $`a_{bf}`$ through $`g=4\pi \mathrm{}^2a_{bb}/m`$, $`h=4\pi \mathrm{}^2a_{bf}/m`$, while the fermion-fermion interaction is omitted as we consider a polarized dilute system at low temperature.

According to the Landau picture of quantum liquid, low lying excited states of the system may be described by quasiparticle excitations which have the structure of particle-hole (p-h) type excitations from the ground state. Small amplitude collective oscillations of the system, corresponding to the zero sound, are given by the coherent superposition of these p-h excitations and are well described by RPA for many fermion systems. The corresponding excitations in the bose condensed system have been formulated in the Bogoliubov-de Gennes type equation. The latter is essentially the same as the one obtained in RPA, and we consider below the two types of excitations on the same footing.

In order to formulate RPA for the bose-fermi mixed system, we first determine the ground state in the mean field approximation. The obtained mean field provides us with single particle energies and wave functions. Let us expand the field operators in terms of a complete set of single particle wave functions as

$$\widehat{\varphi }(\stackrel{}{r})=\underset{k}{}\varphi _k(\stackrel{}{r})b_k,\widehat{\psi }(\stackrel{}{r})=\underset{\alpha }{}\psi _\alpha (\stackrel{}{r})a_\alpha ,$$
(5)
where $`b_k`$ and $`a_\alpha `$ are the boson and the fermion annihilation operators in the single particle states specified by the wave functions $`\varphi _k`$ and $`\psi _\alpha `$. They satisfy the standard commutation or anticommutation relations. The single particle states with quantum numbers $`k`$ or $`\alpha `$ are determined by minimizing the energy expectation value. The trial wave function of the system is given by the product of $`N_b`$ bosons in a single state $`\varphi _{k=0}(\stackrel{}{r})`$ and the Slater determinant of $`N_f`$ different fermionic states $`\psi _\alpha (\stackrel{}{r})`$. From the stationary condition of energy for the variation of these wave functions with a given number of bosons $`N_b`$ and fermions $`N_f`$, we obtain the set of coupled equations, i.e., the Gross-Pitaevskii equation for the boson wave function

$$\left[\frac{\mathrm{}^2}{2m}^2+\frac{1}{2}m\omega _0^2r^2+gN_b|\varphi _0(\stackrel{}{r})|^2+h\rho _f(\stackrel{}{r})\right]\varphi _0(\stackrel{}{r})=\mu _b\varphi _0(\stackrel{}{r})$$
(6)
and the mean-field equation for the occupied fermion states

$$\left[\frac{\mathrm{}^2}{2m}^2+\frac{1}{2}m\omega _0^2r^2+hN_b|\varphi _0(\stackrel{}{r})|^2\right]\psi _\alpha (\stackrel{}{r})=ϵ_\alpha ^f\psi _\alpha (\stackrel{}{r}),$$
(7)
where

$$\rho _f(\stackrel{}{r})=\underset{\alpha }{\overset{\text{occupied}}{}}|\psi _\alpha (\stackrel{}{r})|^2$$
(8)
is the ground state density of fermions, $`\mu _b`$ is the boson chemical potential and $`ϵ_\alpha ^f`$ are the single particle energies for the fermion. Set of equations (6) and (7) together with the fermion number condition

$$d^3r\rho _f(\stackrel{}{r})=N_f$$
(9)
constitute a closed set of equations. The single particle wave functions are normalized to unity. In Eq. (6) we approximate the factor $`N_b1`$ by $`N_b`$. Once the single particle wave functions for the occupied states are obtained from Eqs. (6) and (7), the wave functions of the unoccupied states are calculated from a similar set of equations, i.e.,

$$\left[\frac{\mathrm{}^2}{2m}^2+\frac{1}{2}m\omega _0^2r^2+gN_b|\varphi _0(\stackrel{}{r})|^2+h\rho _f(\stackrel{}{r})\right]\varphi _k(\stackrel{}{r})=ϵ_k^b\varphi _k(\stackrel{}{r})$$
(10)
for bosons with $`k0`$, and the same equation (7) for fermions but for unoccupied states. Orthogonality of these wave functions to the occupied ones are automatically satisfied .

For a large number of particles with a smooth mean field potential, the Thomas-Fermi calculation provides a good approximation. Assuming this to be the case, we determine the fermionic ground state density $`\rho _f`$ by using

$$\frac{\mathrm{}^2}{2m}(6\pi ^2\rho _f(\stackrel{}{r}))^{2/3}+\frac{1}{2}m\omega _0^2r^2+hN_b|\varphi _0(\stackrel{}{r})|^2=ϵ_F$$
(11)
together with the equation (6) instead of fully solving the coupled set of equations (6) and (7). The Fermi energy $`ϵ_F`$ is determined by integrating the density $`\rho _f(\stackrel{}{r})`$ so as to give the fermion number. Numerical consistency of this procedure will be shown in the next section. With $`\rho _f`$ and $`\varphi _0`$ so determined, the wave functions $`\varphi _k(k0)`$ and $`\psi _\alpha `$ are obtained as described above.

We assume below that the mean field is spherically symmetric, and the single-particle quantum numbers $`\alpha `$ and $`k`$ involve standard radial and angular momentum quantum numbers $`(n,\mathrm{},m)`$. (Note that the spin quantum number is frozen.) We consider the case in which the fermionic ground state (Slater determinant) is $`m`$-closed and thus is consistent with the spherically symmetric mean field. We denote by $`p`$ ($`h`$) the fermion single particle states unoccupied (occupied) in the Slater determinant, i.e., those above (below) the Fermi energy $`ϵ_F`$. For the bosons at $`T=0`$ the occupied state is the lowest single-particle state $`k=0`$.

We now define the creation and annihilation multipole operators of the excited state $`|\nu LM`$ with angular momentum quantum numbers $`L,M`$ and an additional quantum number $`\nu `$:

$`Q_{\nu LM}^{}`$ $`=`$ $`F_{\nu LM}^{}+B_{\nu LM}^{},`$ (12)
$`F_{\nu LM}^{}`$ $`=`$ $`{\displaystyle \underset{ph}{}}{\displaystyle \underset{m_pm_h}{}}[X_{phL}^\nu l_pm_pl_hm_h|LM(1)^{l_hm_h}a_{pl_pm_p}^{}a_{hl_hm_h}`$ (14)
$`Y_{phL}^\nu l_hm_hl_pm_p|LM(1)^{l_pm_p}a_{hl_hm_h}^{}a_{pl_pm_p}],`$
$`B_{\nu LM}^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_b}}}{\displaystyle \underset{k0}{}}\left[U_{kL}^\nu b_{kLM}^{}b_0V_{kL}^\nu (1)^{L+M}b_0^{}b_{kLM}\right],`$ (15)
where $`X,Y,U,V`$ are the amplitudes to be determined in RPA, and $`l_pm_pl_hm_h|LM`$ is the Clebsh-Gordan coefficient.

In RPA the operator $`Q`$ is determined so as to satisfy the equation of motion

$$[H,Q_{\nu LM}^{}]\mathrm{}\mathrm{\Omega }_\nu Q_{\nu LM}^{},$$
(16)
where the terms in the l. h. s. which involve different combination of operators $`a^{},a,b^{},b`$ from that in $`Q^{}`$ are neglected by the assumption of RPA. The state

$$|\nu LM=Q_{\nu LM}^{}|0\text{with}Q_{\nu LM}|0=0$$
(17)
is then an approximate eigenstate of the Hamiltonian together with the correlated ground state given by $`|0`$. The amplitudes satisfy the orthonormality condition

$$\underset{ph}{}(X_{phL}^\nu X_{phL}^\nu ^{}Y_{phL}^\nu Y_{phL}^\nu ^{})+\underset{k0}{}(U_{kL}^\nu U_{kL}^\nu ^{}V_{kL}^\nu V_{kL}^\nu ^{})=\delta _{\nu \nu ^{}}.$$
(18)
Substituting Eq. (12) into Eq. (16) we obtain the eigenvalue equation in matrix form:

$`\left(\begin{array}{cccc}\hfill A_{XX}& \hfill A_{XU}& \hfill B_{XY}& \hfill B_{XV}\\ \hfill A_{UX}& \hfill A_{UU}& \hfill B_{UY}& \hfill B_{UV}\\ \hfill B_{XY}& \hfill B_{XV}& \hfill A_{XX}& \hfill A_{XU}\\ \hfill B_{UY}& \hfill B_{UV}& \hfill A_{UX}& \hfill A_{UU}\end{array}\right)\left(\begin{array}{c}X\\ U\\ Y\\ V\end{array}\right)=\mathrm{}\mathrm{\Omega }_\nu \left(\begin{array}{c}X\\ U\\ Y\\ V\end{array}\right),`$ (31)
where submatrices $`A`$ and $`B`$ are given by

$`(A_{XX})_{ph,p^{}h^{}}=[ϵ_{pl_p}^fϵ_{hl_h}^f]\delta _{p,p^{}}\delta _{h,h^{}}\delta _{l_p,l_p^{}}\delta _{l_h,l_h^{}},`$ (32)
$`(A_{XU})_{ph,k}={\displaystyle \frac{h\sqrt{N_b}}{4\pi }}{\displaystyle \frac{\widehat{l_p}\widehat{l_h}}{\widehat{L}}}(1)^{l_h}l_p0l_h0|L0{\displaystyle 𝑑rr^2R_{pl_p}^fR_{hl_h}^fR_{kL}^bR_{00}^b},`$ (33)
$`(A_{UX})_{k,ph}=(A_{XU})_{ph,k},`$ (34)
$`(A_{UU})_{k,k^{}}=[ϵ_{kL}^b\mu _b]\delta _{kk^{}}+{\displaystyle \frac{gN_b}{4\pi }}{\displaystyle 𝑑rr^2R_{kL}^bR_{k^{}L}^bR_{00}^bR_{00}^b},`$ (35)
$`(B_{XY})_{ph,p^{}h^{}}=0,`$ (36)
$`(B_{XV})_{ph,k}={\displaystyle \frac{h\sqrt{N_b}}{4\pi }}{\displaystyle \frac{\widehat{l_p}\widehat{l_h}}{\widehat{L}}}(1)^{l_p}l_p0l_h0|L0{\displaystyle 𝑑rr^2R_{pl_p}^fR_{hl_h}^fR_{kL}^bR_{00}^b},`$ (37)
$`(B_{UY})_{k,ph}=(B_{XV})_{ph,k},`$ (38)
$`(B_{UV})_{k,k^{}}={\displaystyle \frac{gN_b}{4\pi }}(1)^L{\displaystyle 𝑑rr^2R_{kL}^bR_{k^{}L}^bR_{00}^bR_{00}^b}`$ (39)
with $`\widehat{l}=\sqrt{2l+1}`$. In Eq. (36) $`R_{nl}^b(r)`$ and $`R_{nl}^f(r)`$ are the radial parts of the single particle wave functions $`\varphi `$ and $`\psi `$, which are defined by $`\varphi _{nlm}(\stackrel{}{r})=R_{nl}^b(r)Y_{lm}(\theta ,\phi )`$ and $`\psi _{nlm}(\stackrel{}{r})=R_{nl}^f(r)Y_{lm}(\theta ,\phi )`$ where $`Y_{lm}(\theta ,\phi )`$ are the spherical harmonics. We omitted $`O(N_b^1)`$ terms in this calculation. Actually we obtain the same set of equations if we replace $`b_0^{},b_0`$ in Eq. (15) with $`\sqrt{N_b}`$, which is equivalent to the Bogoliubov-de Gennes equation for the bosonic system.

The response of the system to an external field is represented by the transition matrix elements of the relevant operators which connect the ground state and an excited state. The external probe for a collective oscillation is assumed to be long-ranged and is given by the standard multipole operators, i.e.,

$$F_L(\stackrel{}{r})=\{\begin{array}{cc}\frac{1}{\sqrt{4\pi }}r^2\hfill & (L=0)\hfill \\ r^LY_{L0}(\theta )\hfill & (L0)\hfill \end{array}.$$
(40)
As the system is composed of two kinds of particles, we define operators

$$F_L^b=d^3r\widehat{\varphi }^{}(\stackrel{}{r})F_L(\stackrel{}{r})\widehat{\varphi }(\stackrel{}{r}),F_L^f=d^3r\widehat{\psi }^{}(\stackrel{}{r})F_L(\stackrel{}{r})\widehat{\psi }(\stackrel{}{r})$$
(41)
and their combination

$$F_L^\pm =F_L^f\pm F_L^b,$$
(42)
where $`F_L^\tau `$ with $`\tau =f,b,+,`$ are resepectively called ‘fermionic’, ‘bosonic’, ‘in-phase’ and ‘out-of-phase’ type operators. The transition amplitudes for these operators are calculated in RPA as:

$`0|F_L^\pm |\nu L0=0|[F_L^\pm ,Q_{\nu L0}^{}]|0=0|[F_L^f,F_{\nu L0}^{}]|0\pm 0|[F_L^b,B_{\nu L0}^{}]|0,`$ (43)
$`0|[F_L^f,F_{\nu L0}^{}]|0={\displaystyle \underset{ph}{}}{\displaystyle \frac{1}{\sqrt{4\pi }}}{\displaystyle \frac{\widehat{l_p}\widehat{l_h}}{\widehat{L}}}l_p0l_h0|L0[(1)^{l_h}X_{phL}^\nu +(1)^{l_p}Y_{phL}^\nu ]{\displaystyle 𝑑rr^2R_{hl_h}^fr^\lambda R_{pl_p}^f},`$ (44)
$`0|[F_L^b,B_{\nu L0}^{}]|0={\displaystyle \underset{k}{}}{\displaystyle \frac{\sqrt{N_b}}{\sqrt{4\pi }}}[U_{phL}^\nu +(1)^LV_{phL}^\nu ]{\displaystyle 𝑑rr^2R_{kL}^br^\lambda R_{00}^b},`$ (45)
where $`\lambda =2`$ for $`L=0`$ and $`\lambda =L`$ for $`L0`$. The strength distribution for the multipole operator $`F_L^\tau `$ is then given by

$$\underset{\nu }{}\delta (\mathrm{\Omega }\mathrm{\Omega }_\nu )|0|F_L^\tau |\nu L0|^2.(\tau =f,b,+,)$$
(46)
They measure the collectivity of the excited states with respect to the multipole operator which characterizes the shape (or volume for $`L=0`$) oscillation of the system.

More detailed information on the structure of the collective excitation may well by represented by the transition densities

$`0|\widehat{\rho }^\tau (\stackrel{}{r})|\nu LM=\delta \rho _{\nu L}^\tau (r)Y_{LM}(\theta ,\phi )(\tau =f,b),`$ (47)
$`\delta \rho _{\nu L}^f(r)={\displaystyle \underset{ph}{}}{\displaystyle \frac{1}{\sqrt{4\pi }}}{\displaystyle \frac{\widehat{l_p}\widehat{l_h}}{\widehat{L}}}l_p0l_h0|L0[(1)^{l_h}X_{phL}^\nu +(1)^{l_p}Y_{phL}^\nu ]R_{hl_h}^fR_{pl_p}^f,`$ (48)
$`\delta \rho _{\nu L}^b(r)={\displaystyle \underset{k}{}}{\displaystyle \frac{\sqrt{N_b}}{\sqrt{4\pi }}}[U_{phL}^\nu +(1)^LV_{phL}^\nu ]R_{kL}^bR_{00}^b,`$ (49)
where $`\widehat{\rho }^f(\stackrel{}{r})`$,$`\widehat{\rho }^b(\stackrel{}{r})`$ are the fermion and boson density operators:

$$\widehat{\rho }^f(\stackrel{}{r})=\widehat{\psi }^{}(\stackrel{}{r})\widehat{\psi }(\stackrel{}{r}),\widehat{\rho }^b(\stackrel{}{r})=\widehat{\varphi }^{}(\stackrel{}{r})\widehat{\varphi }(\stackrel{}{r}).$$
(50)
Recent development of the two-photon Bragg spectroscopy allows us to obtain dynamical structure factors (response function) for the excitations in trapped atomic gases. They are related to the transition densities as

$`S_L^\tau (q,\mathrm{\Omega })={\displaystyle \underset{\nu }{}}|_{\nu L}^\tau (q)|^2\delta (\mathrm{\Omega }\mathrm{\Omega }_\nu ),`$ (51)
$`_{\nu L}^\tau (q)={\displaystyle 𝑑rr^2j_L(qr/\mathrm{})\delta \rho _{\nu L}^\tau (r)},`$ (52)
where $`q`$ denotes momentum transferred to the system and $`j_L(qr/\mathrm{})`$ are spherical Bessel functions of order $`L`$. The long wavelength limit of the dynamical structure factor is proportional to the strength distribution for the multipole operators defined above.

## III Calculation

### A Numerical procedure

We consider a boson-fermion mixture of potassium isotopes, i.e., 41K (boson) and 40K (fermion). We take the same value $`m=0.649\times 10^{25}`$kg for the boson and the fermion masses. The boson-boson interaction strength $`g`$ is obtained from the scattering length $`a_{bb}=15.13`$nm for 41K-41K, while the boson-fermion interaction strength $`h`$ is varied. (A negative value for the scattering length $`a_{bf}`$ for 40K-41K has been suggested in , although not well established yet.) Both particles are assumed to be trapped in the spherical oscillator potential with the same trap frequency $`\omega _0=100`$Hz. The oscillator length $`\xi =\sqrt{\mathrm{}/m\omega _0}`$ is 4.03$`\mu `$m for the adopted value of $`\omega _0`$.

Number of bosons is fixed at $`N_b=`$1000, while that of fermions is calculated in the Thomas-Fermi approximation for a fixed Fermi energy $`ϵ_F=(6\times 1140)^{1/3}\mathrm{}\omega _0`$. The latter is determined so that the last oscillator shell with a number of quanta $`N=2n+\mathrm{}`$ is given by $`N_{\mathrm{Fermi}}=17`$ at $`h/g=0`$. This gives the number of fermions $`N_f=1140`$ at $`h/g=0`$. $`N_f`$ is dependent on the boson-fermion interaction strength, and it is assumed that all the subshells given by the quantum numbers $`(n,\mathrm{})`$ are $`m`$-closed, which gives, for instance, $`N_f=1050`$ at $`h/g=8`$ and $`N_f=1227`$ at $`h/g=6`$ for the given value of $`ϵ_F`$.

The single-particle wave functions were obtained from the coupled Gross-Pitaevskii and Thomas-Fermi equations by expanding the wave functions in the harmonic oscillator basis for the given oscillator constant $`\xi `$. We included associated Laguerre functions $`L_n^\alpha `$ up to $`n=30`$ for fermions and $`n=15`$ for bosons.

The excitation energies and the wave functions for each multipole $`L`$ are obtained by diagonalization of the RPA matrix. To construct the particle-hole basis we included single particle states for bosons up to $`n=7`$ for $`\mathrm{}=0`$ and $`n=6`$ for $`\mathrm{}=1,2,3`$, so that seven p-h configurations are included for each $`L=0,1,2,3`$. For fermions, four (particle) states above the Fermi level $`ϵ_F`$ and up to four (hole) states below $`ϵ_F`$ have been included for each $`\mathrm{}`$. The number of fermi p-h configurations is 240 for monopole, 464 for dipole, 672 for quadrupole and 1124 for octupole. The dimension of the RPA matrix is twice the sum of bose and fermi configurations for each multipole. The number of bose configurations is the same for the four multipoles considered in the present calculation.

### B Static properties

In the ground state the fermion density distribution is much broader than that for the boson due to Fermi pressure , although at large negative value of $`h/g`$ the bose-fermi attraction tends to produce a larger overlap of the two kinds of particles as shown in Ref. . As the bose-fermi interaction becomes repulsive the fermions are squeezed out from the central part. This in turn causes a smaller overlap of bosons and fermions, and thus a relatively small net effect of the interaction on the binding energy. For the present choice of parameters the fermions are distributed outside of the boson ‘core’ around $`h/g7`$, forming a ‘shell’ like structure . Note that this behavior would be changed if one adopts different values for $`N_f/N_b`$ and $`g`$, which may be represented by a single parameter $`\alpha `$ introduced in Ref. .

The condition for the validity of the Thomas-Fermi approximation used to obtain the above fermion density distribution may be expressed in terms of the local de Brogile wave length $`\lambda (r)=\mathrm{}/p(r)`$, where $`p(r)=\sqrt{2m(ϵ_FV_{\mathrm{eff}})}`$ with fermion mean field potential $`V_{\mathrm{eff}}(r)=\frac{1}{2}m\omega _0^2r^2+h\rho _b(r)`$. With this quantity, the condition becomes

$$f(r)=\left|\frac{d\lambda (r)}{dr}\right|1.$$
(53)
At $`h/g=8`$ where the fermionic potential may have a most pronounced structure, the value of $`f(r)`$ is $`10^2`$ except around turning points. At the turning points, however, the fermion density produces only a negligible influence on the mean field potential. The validity of the present Thomas-Fermi calculation may also be checked by comparing the fermion density distribution $`\rho _f(r)`$ obtained from Eq. (11) with the one calculation from the single particle wave functions according to Eq. (7). The comparison of the two distribution in the range $`h/g=68`$ shows that they agree within the order of $`10^2`$ over the entire radial range, suggesting the consistency of the present calculation.

In Fig. 1 we show the fermion single-particle energies measured from Fermi energy, $`ϵ_\alpha ϵ_F`$, against orbital angular momentum $`\mathrm{}`$ for three values of the interaction parameter, $`h/g=5.0,1.0,3.0`$. The Fermi energy $`ϵ_F`$ is denoted by the horizontal line at zero energy. Note that $`ϵ_\alpha `$ are given by those of the harmonic oscillator at $`h/g=0`$, $`ϵ_{n,\mathrm{}}=(2n+\mathrm{}+3/2)\mathrm{}\omega _0`$, as we have no direct interaction among fermions. One should note also that the yrast state, i.e., the lowest state for each angular momentum, has no radial node, and the yrare one, the second lowest state, has one radial node, etc. The nodal structure of the particle and hole states around the Fermi surface is responsible for the multipole strength distribution of low energy particle-hole excitations. The figure shows that the states with low orbital angular momentum are much influenced by the bose-fermi interaction, while those with high angular momentum are almost insensitive to the values of $`h/g`$. This is because the additional fermion potential due to the bose-fermi interaction vanishes outside the boson density distribution. For instance, since the yrast single particle wave functions with angular momentum $`\mathrm{}`$ are peaked around $`\sqrt{\mathrm{}}\xi `$, those fermion states with $`\mathrm{}>(R_B/\xi )^2`$, $`R_B`$ being a typical edge radius of the boson distribution, would not be much influenced by the bose-fermi interaction. (In the present case $`R_B3\xi `$, and the above relation gives $`\mathrm{}>9`$.)

### C Energy weighted moments and comparison with the sum rule calculation

RPA calculation provides approximate eigenvalues and wave functions for all the individual eigenstates and is useful to study details of the dynamics of the system. If one is interested in the gross behavior of the strength distribution or an approximate frequency of the oscillation, one may rather consider the $`p`$-th energy weighted moments of the strength distribution,

$$m_p(L,\tau )=\underset{\nu }{}(\mathrm{}\mathrm{\Omega }_\nu )^p|0|F_L^\tau |\nu |^2,$$
(54)
where $`\mathrm{}\mathrm{\Omega }_\nu `$ is the excitation energy of the state $`|\nu `$. These moments can be expressed as ground-state expectation values of multiple commutators of $`F_L^\tau `$ with the Hamiltonian. It is known that this relation, sum rule, is conserved in RPA for some of the moments, i.e., the sum in the r.h.s. of Eq. (54) obtained from RPA eigenstates coincides with the expectation value of the commutator in the HF ground state. Thus a comparison of the two would provide a criterion on the consistency of the RPA calculation. In particular, the first moment (energy weighted sum, EWS) $`m_1`$ is calculated from the double commutator of $`F_L^\tau `$ and the Hamiltonian as

$`m_1(L,\tau )={\displaystyle \frac{1}{2}}[F_L^\tau ,[\widehat{H},F_L^\tau ]]_0,`$ (55)
where $`_0`$ denotes a ground state expectation value. We have checked that the sum rule is satisfied within $`1`$% for all the multipole modes.

One can estimate the average frequency of the collective oscillation based on sum rules. There are several ways to define the average frequency depending on which part of the strength distribution is emphasized. In Ref.  the ratio of the third and the first moments

$$\overline{\omega }=\sqrt{\frac{m_3}{m_1}}$$
(56)
was studied based on sum rules. This definition is advantageous as it can be used to test the validity of the RPA calculation as mentioned earlier. In later discussions we consider also the average frequency calculated by

$$\omega _{av}=\frac{m_1}{m_0}$$
(57)
which has more weight in the low-frequency strength compared with $`\overline{\omega }`$. The two frequencies should coincide when a single collective state exhausts the strength. We discuss later also the width of the strength distribution defined by

$$\sigma =\left[\frac{m_2}{m_0}\left(\frac{m_1}{m_0}\right)^2\right]^{1/2}.$$
(58)
In Fig. 2 we show the average frequency Eq. (56) for $`L=0,1,2`$ and $`\tau =+,`$ given by the sum rule (55) and the one directly calculated in RPA. We find a good agreement of the two calculations in a wide range of the bose-fermi interaction parameter $`h`$ for a fixed $`g`$. There is a discrepancy at very large values of $`|h|`$, especially in the strongly attractive case. In the latter case, an induced instability of the ground state due to the bose-fermi attraction may be close, which suggests that a larger configuration space may be required to satisfy the sum rule. In fact, by increasing the number of particle-hole states, we obtained a better agreement.

### D Distribution of multipole strengths

Now we show the distributions of the multipole strengths. Figures 3,7,10,11 respectively show the strength distributions for $`L=`$0,1,2 and 3, for either the fermionic/bosonic or the in-phase/out-of-phase operators.

#### 1 monopole

In Fig. 3 we show the monopole strength distribution, $`|0|F_0^\tau |\nu |^2`$ against $`\mathrm{\Omega }_\nu `$, for $`\tau =f,b`$ and for three values of $`h/g`$. Here one expects a volume oscillation of boson and/or fermion densities. For $`h=0`$ and for a large number of particles, the bosonic monopole oscillation will be located around $`\sqrt{5}\omega _0`$ as given by the collisionless hydrodynamics, while the fermionic one is concentrated at $`2\omega _0`$. The latter property is due to the almost degenerate values of the relevant particle-hole energies contributing to the monopole oscillation. The present calculation show that this situation persists even for large values of $`h/g`$, although the bosonic frequency is slightly shifted. The fermionic strength is distributed over a few states around $`2\omega _0`$, showing that the induced fermi-fermi interaction via bosons is not strong enough to make a coherent superposition of the fermion particle-hole states. The total monopole strength for fermion is by an order of magnitude larger than that for boson, even though the number of fermions involved in the excitation is smaller because of the Pauli principle, e.g. $`m_1=2590\mathrm{}\omega _0\xi ^4`$ for fermion and $`m_1=419\mathrm{}\omega _0\xi ^4`$ for boson at $`h/g=1`$. This is because of the broader density distribution for fermions.

The strength distribution may suggest that the fermions and bosons are moving independently even for rather large values of $`h/g`$. This is not necessarily the case, however, as one may see in Fig. 4, where we plot the dynamical structure factor for the three cases, i.e., the states at $`\mathrm{\Omega }/\omega _0=2.03`$ for $`h/g=5.0`$, at $`2.01`$ for $`h/g=1.0`$ and at $`1.98`$ for $`h/g=`$-3.0, which carry the largest strengths for the in-phase monopole operator. Although the fermionic strength is far larger than the bosonic one for this state, the mixture of the bosonic component is not small and is peaked at larger values of the momentum transfer $`q`$ as shown by the dashed lines. The latter shows that the bosons oscillate in the inner region although not quite recognizable as far as one studies only the $`F_0`$ distribution.

In Refs.  we suggested that an instablity towards collapse may occur at large negative values of $`h/g`$ and may be signaled by the lowering of $`\overline{\omega }`$. In the present calculation with smaller number of particles this is not apparent in Fig. 2. In Fig. 5 we show the strength distribution at $`h/g=6.65`$. This value is close to the critical value of instability around $`6.7`$ estimated from the Mølmer’s condition  and around $`6.8`$ from the condition by Roth and Feldmeier . We find a lowering of a single state which carries 13.9% of the EWS. The average frequency for this case is $`\overline{\omega }=1.92\omega _0`$ appreciably lower than the value $`2\omega _0`$ for $`h/g=55`$. At more negative values of $`h/g`$ we could not find a stable ground state.

Let us study the character of the low-lying excited states by considering the response to the probe

$$F(\theta )=F^+\mathrm{cos}\theta +F^{}\mathrm{sin}\theta $$
(59)
parametrized by $`\theta `$. In Ref.  the angle $`\theta `$ was determined by minimizing the average frequency $`\overline{\omega }`$ for $`F(\theta )`$ so as to find the character of the probe which favors the low-lying states. Once the value of the parameter, $`\theta _{\mathrm{min}}`$ is determined for a given $`h/g`$, one may consider an operator ‘perpendicular’ to $`F(\theta _{\mathrm{min}}`$), i.e., $`F_{}=F^+\mathrm{sin}\theta _{\mathrm{min}}F^{}\mathrm{cos}\theta _{\mathrm{min}}`$, which may favor the high-lying state. Alternatively, one may maximize the EWS within the model space to determine the value $`\theta _{\mathrm{max}}`$ of the operator to characterize the high-lying states. Figure 6 shows the strength distribution for these three types of operators at $`h/g=5.0`$. The value $`\theta _{\mathrm{min}}=0.12\pi `$ suggests that, at this highly repulsive value of the interaction, a slightly in-phase type oscillation is favored in order to avoid an overlap of bosons and fermions, see the density distribution in Ref . We note that the distribution for the ‘perpendicular’ operator $`F_{}`$ is concentrated and is similar to the one for $`F(\theta _{\mathrm{max}})`$. The latter, however, is an almost bosonic type operator and suggest that the strength distribution alone is not sufficient to characterize the structure of high-lying states. A similar study has been made for a very attractive case $`h/g=6.65`$. In this case we obtain $`\theta _{\mathrm{min}}=0.01\pi `$, i.e., almost in-phase, and the strength distribution is similar to Fig. 5. Thus the in-phase character of the low-lying mode which favors the overlap of the two kinds of particles is clearly seen in the determination of $`F(\theta _{\mathrm{min}})`$.

#### 2 dipole

Dipole strength distributions for the in-phase and the out-of-phase operators are shown in Fig. 7. The strong peak in the in-phase strength distribution corresponds to the center-of-mass oscillation of the whole system. For many-particle systems confined in a common oscillator potential the center-of-mass motion is decoupled from other (intrinsic) degrees of freedom of the system, resulting in the oscillation with the same frequency $`\omega _0`$. This relation can be represented by the commutation relation of the dipole operator and the Hamiltonian, and should hold also within RPA . In the present numerical calculation, however, the state at $`\mathrm{\Omega }=\omega _0`$ does not exhaust the whole strength for the in-phase oscillation; the largest deviation occurs at $`h/g=1`$ having 80% of the EWS. This is because the single-particle model space in our numerical calculation is not sufficient to completely decouple the center of mass motion, especially at strong coupling cases where the deviation of the single particle potential from the oscillator becomes large. The $`m_1`$ sum rule is, nevertheless, almost satisfied as mentioned earlier, and is determined by the total number of bose and fermi particles, where the latter number is dependent on the value of $`h/g`$. Smaller value of the sum rule percentage at $`h/g=1`$ is mainly due to the appearance of the almost degenerate state ($`\mathrm{\Delta }\omega /\omega _010^5`$) at this energy.

All other dipole oscillations should be orthogonal to the center-of-mass motion and thus are out-of-phase type in character. One may note that at repulsive values of $`h`$ the dominant part of the out-of-phase strength is located below $`\omega _0`$, although the strength carried by each individual state is not very large. As discussed in Ref. , this may reflect the ground state density distribution which favors the out-of-phase oscillation by making the overlap of boson and fermion distributions smaller. Figure 8 shows transition densities for the out-of-phase type oscillations at $`h/g=5.0,1.0,3.0`$ which carry the largest strengths; $`\mathrm{\Omega }/\omega _0=0.700`$ for $`h/g=5.0`$, $`0.956`$ for $`h/g=1.0`$ and $`1.14`$ for $`h/g=`$-3.0. Bosons and fermions move in the opposite directions so as to make the overlap smaller around the surface of the boson distribution. One can observe that the bosonic transition density is large and robust, while the fermionic one is feeble and is spread over the whole system. This situation may be analogous to the soft dipole mode speculated in neutron-rich nuclei wherein the protons oscillate almost free in the sea of neutrons. Dynamical structure factors corresponding to these states are shown in Fig. 9. Here the bosonic and the fermionic responses are shown together with the one for a hypothetical out-of-phase type probe. By changing the momentum transfer one would find a structure corresponding to the oscillation of the two kinds of particles.

It was suggested in  that for a very strong repulsion between bosons and fermions with sufficient number of particles, the system may become unstable towards a phase separation of the two kinds of particles. Out-of-phase dipole strength distribution in Fig. 7 indeed show a softening of the strength distribution at large value of $`h/g`$. We cannot conclude from the present calculation, however, if this tendency is related to the mentioned instability.

#### 3 quadrupole

One finds from the quadrupole strength distribution in Fig. 10 that the fermionic strength is split into low- and high-energy parts, while the bosonic one is concentrated. The higher strengths for the fermions are due to the 2$`\mathrm{}\omega _0`$ excitation of the particle from the occupied single-particle states and are similar in character to the bosonic excitation. In contrast, the lower fermion strengths come from the matrix elements of the quadrupole operator $`F_2`$ which re-orient the single particle states within the same oscillator shell of $`N=2n+\mathrm{}`$ around the Fermi surface. This transition normally involves a change of nodes by one, see Fig. 1, and the corresponding strengths are smaller than the higher ones. One may note that the quadrupole bosonic strength is located slightly higher than $`\sqrt{2}\omega _0`$ as expected for the collective oscillation in the large $`N`$ limit, which may be due to a rather small number of particles in the present calculation.

We note that the strength distribution tends to be fragmented at large (repulsive) value of the interaction parameter $`h/g`$ as seen in Fig. 10. This is because the average potential deviates appreciably from the harmonic oscillator potential, giving rise to a dispersion in the fermion particle-hole energies, see Fig. 1. Strong fermion-boson interaction would then be responsible to scatter the bosonic strengths, too.

#### 4 octupole

Octupole strength distribution is shown in Fig. 11. Here again the fermionic strength is split into $`\mathrm{}\omega _0`$ and $`3\mathrm{}\omega _0`$ regions due to the character of the multipole operator $`F_3`$. We note also that the frequency of the bosonic oscillation is much larger than the hydrodynamic value $`\sqrt{3}\omega _0`$.

We find a striking change in the fermionic strength distributions depending on the sign of the bose-fermi interaction. While for a strongly repulsive case ($`h/g=5.0`$) the region between the low- and high-lying states is almost filled up by small strengths, the one for a strongly attractive case ($`h/g=3.0`$) shows a large gap in the strength distribution. This may be traced back to the structure of the single-particle states in Fig. 1. The large strengths around $`3\omega _0`$ are due to high orbital angular momentum states with $`\mathrm{}_p\mathrm{}_h=3`$ which stay robust against $`h/g`$ as discussed earlier. Smaller strengths due to low $`\mathrm{}_p`$ and $`\mathrm{}_h`$ states are, on the other hand, sensitive to the interaction, and the correponding particle-hole energies are below (for $`h/g>0`$) or above (for $`h/g<0`$) the unperturbed value $`3\omega _0`$. This may be understood by noting that the particles passing through the center (small $`\mathrm{}`$ orbit) feel the shallow (for $`h/g>0`$) or deep (for $`h/g<0`$) potentials and are thus more easy or hard to excite. For the strengths around $`\omega _0`$ different combinations of orbits are involved, and the strengths are broadened by the deviation from the oscillator potential.

### E Time evolution

To study the collective oscillation of the cold atomic gases, the current experiments on BEC exert a time-dependent field on the system and then observe the time development of the shape and the size of the condensate. This procedure generally excites a number of normal modes, and one may in principle be able to resolve each mode by performing Fourier transform as far as the time duration is long enough. To simulate the situation we consider a time evolution of the system after one applies an external weak pulse of the step-function type, i.e., $`\mathrm{\Delta }VF_L^\tau `$ for $`\mathrm{\Delta }tt0`$. The calculation is performed within the lowest order perturbation theory. ($`\mathrm{\Delta }t`$ is chosen as $`\mathrm{\Delta }t=10^2\omega _0^1`$.)

#### 1 monopole

First we consider a perturbation of monopole type

$$\mathrm{\Delta }V=m\omega _0\mathrm{\Delta }\omega \sqrt{4\pi }F_0^\tau $$
(60)
with $`\mathrm{\Delta }\omega \omega _0`$ which is applied to the ground state of the system for a short period $`\mathrm{\Delta }t`$. The effective frequency of the external oscillator potential for fermions and/or bosons is then changed into $`\omega _0\pm \mathrm{\Delta }\omega `$ depending on the choice of $`\tau `$. Time dependence of the rms radius of bosons and fermions is plotted in Fig. 12 at $`h/g=3`$ for the fermionic external field, $`\tau =f`$. Frequency of the short period oscillation is related to the average frequency of the strength distribution, e.g., $`\omega _{av}=1.99\omega _0`$. That of the long period one, on the other hand, reflects the width of the distribution. For instance, the slow decrease of the envelope for $`\sqrt{r_f^2}`$ would imply a very narrow width of the strength distribution. The bottom of Fig. 12 shows the Fourier transform of $`r_f^2(t)r_f^2_0`$ which recovers the sharp peak structure of the strength distribution in Fig 3. Time evolution of the bosonic radius shows a rather regular modulation, and its beat frequency is estimated about $`0.125\omega _0`$, which is consistent with the half width $`\sigma /2=0.127\omega _0`$.

#### 2 dipole

For the dipole case we take the perturbing potential

$$\mathrm{\Delta }V=m\omega _0^2\mathrm{\Delta }z\sqrt{\frac{4\pi }{3}}F_1^\tau ,$$
(61)
where $`\mathrm{\Delta }z(\xi )`$ measures the shift of the centers of the oscillator potentials felt by fermions and/or bosons. In Fig. 13 we show the time dependence of the center-of-masses of bosons and fermions for the out-of-phase external dipole impulse, $`\tau =`$, at $`h/g=5`$. The frequency of the oscillation is consistent with the average frequency $`0.90\omega _0`$ for the dipole strength distribution. The amplitude of the oscillation shows an irregular behavior, suggesting that the strength distribution does not have a simple structure. The lowest part of the figure shows the Fourier transform of the relative distance of the fermion and boson center-of-mass positions, $`N_fz_fN_bz_b`$, which is sufficient to recover the behavior of the original strength distribution given in Fig. 7.

## IV Summary and Conclusion

In the present paper we performed an RPA calculation of the polarized bose-fermi mixed system of alkali-metal gases at zero temperature. We solved a coupled Gross-Pitaevskii-Thomas-Fermi equations to obtain the ground state for several values of the boson-fermion interaction strength with a fixed value of the boson-boson interaction. Single particle states of bose and fermi particles are calculated based on the mean field produced by the obtained ground state density. The density distribution constructed from the occupied single-particle orbits agrees well with the original one, suggesting the self-consistency of the calculation. We calculated and diagonalized the RPA matrix to obtain excitation energies and wave functions for multipoles $`L=0,1,2,3`$. We then calculated strength distributions for bosonic/fermionic and in-phase/out-of-phase multipole operators, and also transition densities and dynamical structure factors for some of the collective states.

We first calculated energy weighted moments of the strengths to study average behavior of the strength distribution. Comparison of the first moments with the sum rule predictions for $`L=0,1,2`$ shows that the RPA configuration space is sufficient to within 1%. A glance at the monopole distribution suggests that fermions and bosons are moving rather independently for the adopted parameters , which turns out not to be the case if we study the dynamical structure factors. We find, in fact, that the long wavelength oscillation of fermions in the surface is coupled to the internal short-wavelength oscillation of bosons. For a strong attractive boson-fermion interaction, the calculation suggests a softening of in-phase monopole oscillation. We studied also the structure of the low-lying out-of-phase type dipole modes. Lowering of the energy is related to the ground state density distribution of bosons and fermions which favors the out-of-phase oscillation. Transition densities for these modes suggest that the bosons oscillate on their way through the cloud of fermions, which is analogous to the soft dipole mode discussed in neutron-rich nuclei. We also calculated the strength distributions for quadrupole and octupole modes, and studied, in particular, the origin of the fragmentation of fermionic strengths. Finally we considered time dependent behavior of the trapped bose-fermi system after an impulse of a multipole external field. By Fourier transforming the time-dependent oscillating behavior of the system, we could recover the gross structure of the strength distribution.


format: text

text:
# On the amelioration of quadratic divergences

## Abstract

Once massless quadratically divergent tadpole diagrams are discarded, because they contain no intrinsic scale, it is possible to convert other divergences into logarithmic form, using partial fraction identities; this includes the case of quadratic divergences, as has been applied to the linear sigma model. However the procedure must be carried out with due care, paying great attention to correct numerator factors.

In QED one knows that the formally quadratically divergent photon vacuum polarization graph is reduced to a logarithmically divergent singularity by invoking gauge invariance. However, in a linear $`\sigma `$ model field theory, involving spinless pion and $`\sigma `$ mesons, no principle like gauge invariance can be invoked; yet the quadratic divergence can be tamed and converted into a logarithmic one. This fact was first noted for a quark-level SU(2) model, by using the $`2\mathrm{}`$-dimensional regularization lemma (DRL) :

$$I\overline{d}^4p\left[\frac{m^2}{(p^2m^2)^2}\frac{1}{p^2m^2}\right]=\underset{\mathrm{}2}{lim}\frac{i(m^2)^\mathrm{}1}{(4\pi )^{\mathrm{}}}[\mathrm{\Gamma }(2\mathrm{})+\mathrm{\Gamma }(1\mathrm{})]=\frac{im^2}{16\pi ^2},$$
(1)
since  the combination $`\mathrm{\Gamma }(2\mathrm{})+\mathrm{\Gamma }(1\mathrm{})=\mathrm{\Gamma }(3\mathrm{})/(1\mathrm{})1`$, as $`\mathrm{}2`$. In this note we shall show that lemma (1) holds for all regularization schemes which set the massless quadratically divergent tadpole integral $`d^4p/p^2`$ equal to zero, for the reason that it contains no intrinsic mass scale (dimensional analysis). This therefore includes dimensional regularization , analytic and zeta function regularization, Lodder log  and Pauli-Villars regularizations .

To prove the result (1), start with the partial fraction identity

$$\frac{m^2}{(p^2m^2)^2}\frac{1}{p^2m^2}=\frac{1}{p^2}\left[\frac{m^4}{(p^2m^2)^2}1\right],$$

and integrate the rhs over $`\overline{d}^4p`$. Making a Wick rotation to $`q^2=p^2,\overline{d}^4p=iq^2dq^2/16\pi ^2`$ and dropping the tadpole term, we remain with

$$I=i\frac{m^4}{16\pi ^2}_0^{\mathrm{}}\frac{dq^2}{(q^2+m^2)^2}=\frac{im^2}{16\pi ^2}$$
(2)
coinciding with the rhs of (1). An alternative way of reaching this conclusion is to apply the Karlson  trick (advocated also by Schwinger):

$$\frac{d}{dm^2}\frac{d^4p}{p^2m^2}=\frac{d^4p}{(p^2m^2)^2},$$
(3)
since this serves to eliminate a possible quadratic tadpole infinity. Thus one may readily verify, without having to make any statement about quadratic divergences tadpole infinities, that

$$dI/dm^2=2m^2\overline{d}^4p/(p^2m^2)^3=i/16\pi ^2.$$

The only question/problem is what happens when we integrate this up in $`m^2`$. The regularization (1) presumes that no constant $`\mathrm{\Lambda }^2`$ intrudes—amounting to setting the tadpole contribution to zero.

At this point we want to inject a note of caution; this concerns numerators of divergent Feynman integrals. Consider the vacuum polarization integral

$`\mathrm{\Pi }^{\mu \nu }(k)`$ $`=`$ $`ie^2\mathrm{Tr}{\displaystyle \frac{[\gamma ^\mu ((p+k)\gamma +m)\gamma ^\nu (p\gamma +m)]}{((p+k)^2m^2)(p^2m^2)}\overline{d}^2\mathrm{}p}`$ (4)
$`=`$ $`c\eta ^{\mu \nu }+(k^2\eta ^{\mu \nu }k^\mu k^\nu )f(k^2).`$ (5)
Gauge invariance requires $`c=0`$ (and dimensional regularization will guarantee it). Thus whatever the chosen regularization method, we must insist upon $`\mathrm{\Pi }^{\mu \nu }(0)`$ or $`c=0`$. In dimensional regularization one strictly finds

$$c_2\mathrm{}2^{\mathrm{}}\frac{(1/\mathrm{}1)p^2+m^2}{(p^2m^2)^2}\overline{d}^2\mathrm{}p2^{\mathrm{}}[(1/\mathrm{}1)\mathrm{\Gamma }(1\mathrm{})\mathrm{\Gamma }(2\mathrm{})/\mathrm{}](m^2)^\mathrm{}10.$$
(6)
The delicate cancellation only works because the correct dimensional factors have been included in the numerator of (4). However if one were blindly to set $`\mathrm{}=2`$ at the start of eq. (5), one would find

$$c_4\frac{4m^22p^2}{(p^2m^2)^2}\overline{d}^4p=2\frac{\overline{d}^4p}{p^2}\left[1+\frac{m^4}{(p^2m^2)^2}\right],$$

so dropping the tadpole term directly in 4-D would leave us with a finite but nonzero contribution $`2m^4\overline{d}^4p/[p^2(p^2m^2)^2]m^2`$, which is undesirable! It is therefore vitally important to incorporate the correct factors in the numerators of such integrals, before jumping to conclusions. If in any doubt one should apply a regularization scheme which respects the symmetries that one holds dear—such as dimensional regularization.

To conclude, we can with appropriate care interpret the regularized quadratically divergent integral as

$$16i\pi ^2\overline{d}^4p/(p^2m^2)=m^2[\mathrm{ln}(m^2/M^2)+C],$$

where $`M`$ is some mass scale and $`C`$ is some indeterminate constant. This form is of course consistent with the once $`m^2`$-differentiated integral:

$$16i\pi ^2\left(m^2\frac{d}{dm^2}1\right)\frac{\overline{d}^4p}{p^2m^2}=16i\pi ^2I=m^2,$$
(7)
and it obeys the DRL, eq (1); but the naive ultraviolet cutoff integral $`d^4p/(p^2m^2)=\mathrm{\Lambda }^2`$ violates all decent regularization schemes discussed above.

The fact that dimensional regularization regards all one-loop divergences as simple poles in $`\mathrm{}`$ means that, with appropriate care, one should be able to convert them all into ‘logarithmic’ form by suitable discarding all integrals of the type $`d^4p(p^2)^n`$, at least for $`n>2`$, because they contain no intrinsic scale. Taming such divergences may have deeper repercussions; in particular, one of the reasons for invoking the elegant idea of supersymmetry is that it provides for a nice cancellation between boson and fermion loop contributions (including the numbers of each), with consequences for fine-tuning. However if these loop contributions are effectively logarithmic, the argument for supersymmetric cancellations is weakened considerably and one requires instead that $`_i(1)^{F_i}g_i^2m_i^2[\mathrm{ln}m_i^2+C]`$ be finite, rather than $`_i(1)^{F_i}g_i^2`$, with all $`g_i`$ normally taken to be equal and where the summation is over all spin and charge states.

It should be emphasized that we have not actually eliminated the quadratic divergence, we have simply transmogrified it. \[Thus $`\lambda \varphi ^4`$ theory still has a dilatation current anomaly, or anomalous dimension for the $`\varphi `$-field, which occurs to order $`\lambda ^2`$ when $`\varphi `$ is massless; but a massive $`\varphi `$ will contribute to the self-mass renormalization to order $`\lambda `$ when the tadpole loop contains a mass.\] Since the above results apply to one-loop diagrams one might well ask if the tadpole subtraction technique can be extended to higher loops. To show that this can be done, we shall exemplify the process by considering Yukawa theory,

$$=Z_\varphi [(\varphi )^2(\mu ^2\delta \mu ^2)\varphi ^2]/2+Z_\psi \overline{\psi }[i\gamma .m+\delta m]\psi +Z_gg\overline{\psi }\gamma _5\psi \varphi +Z_\lambda \lambda \varphi ^4,$$

focusing on the fermion loop contributions up to order $`g^4`$ to the meson self-energy, which are the main source of quadratic infinities. (We have taken a pseudoscalar model rather than a scalar model because it renders the calculations more transparent in the soft meson momentum limit). Specifically we shall examine the $`g^4`$ overlapping divergence as a good test case of the procedure.

Before doing so, let us take note of a couple of one-loop results as they are needed at the next perturbative order. The meson self-energy to order $`g^2`$ is given by

$$\mathrm{\Pi }^{(2)}(k)=ig^2\mathrm{Tr}\overline{d}^4p\gamma _5\frac{1}{\gamma .(p+k)m}\gamma _5\frac{1}{\gamma pm}=4ig^2\frac{[p(p+k)m^2]\overline{d}^4p}{[(p+k)^2m^2][p^2m^2]}.$$
(8)
For simplicity, renormalise this at $`k=0`$, allowing the self-mass counterterm to be massaged from quadratic to logarithmic by tadpole subtraction:

$$\delta \mu ^2=\mathrm{\Pi }^{(2)}(0)=4ig^2\frac{\overline{d}^4p}{p^2m^2}4ig^2\frac{m^2\overline{d}^4p}{p^2(p^2m^2)}.$$

The wave-function renormalization $`Z_\varphi `$, being logarithmic, requires no treatment; nor is it needed for the next part of the argument. The vertex correction is also logarthmic but is required to the next order when considering the overlapping divergent graph, so let us just note that at zero momentum transfer,

$$\mathrm{\Gamma }^{(2)}(p)=ig^2\gamma _5\frac{\overline{d}^4p}{(qp)^2(q^2m^2)},$$
(9)
from which it follows that to order $`g^2`$ one can take

$$Z_g=1ig^2\frac{\overline{d}^4q}{q^2(q^2m^2)}$$

as the coupling constant counterterm. Observe that we have neglected the mass of the meson, as that has little bearing on the ultraviolet behaviour: the difference $`1/(k^2\mu ^2)1/k^2=\mu ^2/k^2(k^2\mu ^2)`$ has improved high-energy convergence; non-vanishing $`\mu `$ just muddies the argument to come.

Now concentrate on the overlapping $`g^4`$ order fermion loop contribution to the meson self-energy. Including the vertex counterterm, it is given by

$`\mathrm{\Pi }_{\mathrm{overlap}}^{(4)}(0)`$ $`=`$ $`g^4\mathrm{Tr}{\displaystyle \overline{d}^4p\left[\gamma _5\frac{1}{\gamma pm}\gamma _5\frac{1}{\gamma pm}(\frac{\overline{d}^4q}{(qp)^2}\gamma _5\frac{1}{\gamma qm}\gamma _5\frac{1}{\gamma qm}+i(Z_g1))\right]}`$ (10)
$`=`$ $`4g^4{\displaystyle \frac{\overline{d}^4p}{p^2m^2}\frac{\overline{d}^4q}{q^2m^2}\left(\frac{1}{q^2}\frac{1}{(qp)^2}\right)}4g^4{\displaystyle \frac{f(p^2)\overline{d}^4p}{p^2m^2}}.`$ (11)
The vertex renormalization has guaranteed that the $`q`$-integral is finite, but the remaining $`p`$-integral is potentially quadratically infinite and it is this part that we wish to convert to logarithmic by tadpole subtraction. The trick is to notice that $`f(p^2)`$ can be converted into the dispersive form

$$f(p^2)p^2_{m^2}^{\mathrm{}}\frac{ds}{s(sp^2)}=\mathrm{ln}(1\frac{p^2}{m^2}),$$

whereupon the overlap integral can be re-expressed as

$$\mathrm{\Pi }_{\mathrm{overlap}}^{(4)}(0)\overline{d}^4p𝑑s\left[\frac{1}{sm^2}\left(\frac{1}{p^2m^2}+\frac{1}{sp^2}\right)\frac{1}{s(p^2m^2)}\right].$$

Having manoeuvred the quadratic infinity into suitable form, we can apply the amelioration procedure, $`\overline{d}^4p/(p^2M^2)\overline{d}^4pM^2/p^2(p^2M^2)`$, and end up with

$$\mathrm{\Pi }_{\mathrm{overlap}}^{(4)}(0)\frac{\overline{d}^4p}{p^2}\frac{ds}{sm^2}\left[\frac{m^4}{s(p^2m^2)}+\frac{s}{sp^2}\right].$$

Thus we have achieved our goal of conversion of quadratic to logarithmic in both of the ‘internal momentum’ integrals. This does not imply that the logarithmic infinities remain at first order like $`\mathrm{log}(M^2/m^2)`$; indeed we encounter higher powers of logs in the overlapping example above but the important point is that the quadratic infinity has been ameliorated. Of course, this resulting divergence needs as ever to be subtracted by a new fourth order counterterm, but this is just standard renormalization fare and there is no problem with that. For higher loop integrals we are optimistic that the same procedure should work, provided one replaces the logarithmic dispersive integral $`f`$ by the more general Feynman parametric form ($`\alpha `$ represents a set of Feynman parameters),

$$f(p^2)p^2𝑑\alpha /[p^2M^2(\alpha )],$$

and follows the steps above.


text:
# Understanding deterministic diffusion by correlated random walks

## I Introduction

Deterministic diffusion is a prominent topic in the theory of chaotic dynamical systems and in nonequilibrium statistical mechanics . To achieve a proper understanding of the mechanism of deterministic diffusion with respect to microscopic chaos in the equations of motion, a number of simple model systems was proposed and analyzed. An interesting finding was that, for a simple one-dimensional chaotic map on the line, the diffusion coefficient is a fractal function of a control parameter . Diffusion coefficients exhibiting similar irregularities are also known for more complicated sytems such as sawtooth maps , standard maps , and Harper maps , the latter two models being closely related to physical systems like the kicked rotor, or a particle moving in a periodic potential under the influence of electric and magnetic fields. More recently, related irregularities in the diffusion coefficient were discovered in Hamiltonian billiards such as the periodic Lorentz gas , and for a particle moving on a corrugated floor under the influence of an external field . That this phenomenon is not specific to diffusion coefficients is suggested by results on other transport coefficients such as chemical reaction rates in multibaker maps , electrical conductivities in the driven periodic Lorentz gas , and magnetoresistances in antidot lattices  which, again, exhibit strongly irregular behavior under parameter variation.

In this paper we wish to contribute to the further analysis and understanding of irregular transport coefficients in simple model systems. We use certain forms of the Green-Kubo formula as starting points for expansions of the diffusion coefficient in terms of correlated random walks. In Section 2 we exemplify this approach for the piecewise linear map on the line studied in Refs. . In Section 3 we apply a suitably adapted version of it to diffusion in the periodic Lorentz gas. In both cases our approximations converge quickly to the precise diffusion coefficient as obtained from other numerical methods. We argue that for both models this approximation scheme provides a physical explanation for the irregular structure of these diffusion coefficients in terms of dynamical correlations, or memory effects. In Section 4 we summarize our results and relate them to previous work on irregular and fractal transport coefficients.

## II Deterministic diffusion in one-dimensional maps on the line

Probably the most simple models exhibiting deterministic diffusion are one-dimensional maps defined by the equation of motion

$$x_{n+1}=M_a(x_n)$$
(1)
where $`a\text{}`$ is a control parameter and $`x_n`$ is the position of a point particle at discrete time $`n`$. $`M_a(x)`$ is continued periodically beyond the interval $`[0,1)`$ onto the real line by a lift of degree one, $`M_a(x+1)=M_a(x)+1`$. We assume that $`M_a(x)`$ is anti-symmetric with respect to $`x=0`$, $`M_a(x)=M_a(x)`$. Let

$$m_a(x):=M_a(x)\text{mod}\mathrm{\hspace{0.33em}1}$$
(2)
be the reduced map related to $`M_a(x)`$. This map governs the dynamics on the unit interval according to $`x_n=m_a^n(x)`$. Let $`\rho _n(x)`$ be the probability density on the unit interval of an ensemble of moving particles starting at initial conditions $`xx_0`$. This density evolves according to the Frobenius-Perron continuity equation

$$\rho _{n+1}(x)=_0^1𝑑y\rho _n(y)\delta (xm_a(y)).$$
(3)
Here we are interested in the deterministic diffusion coefficient defined by the Einstein formula

$$D(a)=\underset{n\mathrm{}}{lim}<(x_nx_0)^2>/(2n)$$
(4)
with $`x_n`$ governed by $`M_a(x)`$, where the brackets denote an average over the initial values $`x_0`$ with respect to the invariant probability density on the unit interval, $`\mathrm{}:=_0^1𝑑x\rho _a^{}(x)\mathrm{}`$. Eq.(4) can be transformed onto the Green-Kubo formula for maps ,

$$D(a)=\frac{1}{2}j^2(x_0)+\underset{k=1}{\overset{\mathrm{}}{}}j(x_0)j(x_k),$$
(5)
where the jump velocity

$$j(x_n):=[x_{n+1}][x_n],$$
(6)
$`[x]`$ being the largest integer less than $`x`$, takes only integer values and denotes how many unit intervals a particle has traversed after one iteration starting at $`x_n`$. The map we study as an example is defined by

$$M_a(x)=\left\{\begin{array}{cc}\hfill ax,& 0<x\frac{1}{2}\hfill \\ \hfill ax+1a,& \frac{1}{2}<x1\hfill \end{array}\right\},$$
(7)
where the uniform slope $`a>2`$ serves as a control parameter. The Lyapunov exponent of this map is given by $`\lambda (a)=\mathrm{ln}a`$ implying that the dynamics is chaotic. A sketch of this map is shown in Fig. 1 (a). Fig. 1 furthermore depicts the parameter-dependent diffusion coefficient of this map in the regime of $`2a4`$ as calculated in Refs.  by means of a numerical implementation of analytical methods. In these references numerical evidence was provided that the irregularities on a fine scale are reminiscent of an underlying fractal structure.

To analyze this fractal diffusion coefficient, Eq.(5) forms a suitable starting point because it distinguishes between two crucial contributions of the dynamics to the diffusion process, which are (1) the motion of the particle on the unit interval $`x_n\text{mod}\mathrm{\hspace{0.33em}1}`$ generating the invariant density $`\rho _a^{}(x)`$, and (2) the integer jumps from one unit interval to another one as related to $`j(x_n)`$. In fact, in Refs.  it was shown that both parts are independent sources of fractality for the diffusion coefficient. However, in these references Eq. (5) was only discussed in the limit of infinite time. In this paper we suggest to look at the contributions of the single terms in this expansion, and to analyze how the Green-Kubo formula approaches the exact diffusion coefficient step by step. We do this by systematically building up hierarchies of approximate diffusion coefficients. These approximations should be defined such that the different dynamical contributions to the diffusion process are properly filtered out. Another issue is how to evaluate the single terms of the Green-Kubo expansion on an analytical basis. In this section we show how to do this for the one-dimensional model introduced above, in the next section we apply the same idea to billiards such as the periodic Lorentz gas.

We start by looking at the first term in Eq.(5). For $`a<4`$ the absolute value of the jump velocity $`j(x_n)`$ is either zero or one. Assuming that $`\rho _a^{}(x)1`$ for $`a2`$ and cutting off all higher order-terms in Eq. (5), the first term leads to the well-known random walk approximation of the diffusion coefficient

$$D_0(a)=\frac{1}{2}_0^1𝑑xj_a^2(x)$$
(8)
which in case of the map Eq. (7) reads

$$D_0(a)=(a2)/(2a).$$
(9)
This solution is asymptotically correct in the limit of $`a2`$. More generally speaking, the reduction of the Green-Kubo formula to the first term only is an exact solution for arbitrary parameter values only if all higher-order contributions from the velocity autocorrelation function $`C(n):=<j(x_0)j(x_n)>`$ are strictly zero. This is only true for systems of Bernoulli type . Conversely, the series expansion in form of Eq. (5) systematically gives access to higher-order corrections by including higher-order correlations, or memory effects. This leads us to the definition of two hierarchies of correlated random walk diffusion coefficients:

(1) Again, we make the approximation that $`\rho _a^{}(x)1`$. We then define

$$D_n^j(a):=\frac{1}{2}_0^1𝑑xj^2(x)+\underset{k=1}{\overset{n}{}}_0^1𝑑xj(x)j(x_k),n>0,$$
(10)
with $`D_0(a)`$ given by Eq. (8). Obviously, this series cannot converge to the exact $`D(a)`$.

(2) By using the exact invariant density in the averages of Eq.(5) we define

$$D_n^\rho (a):=\frac{1}{2}_0^1𝑑x\rho _a^{}(x)j_a^2(x)+\underset{k=1}{\overset{n}{}}_0^1𝑑x\rho _a^{}(x)j_a(x)j_a(x_k),n>0,$$
(11)
here with $`D_0^\rho (a)=\frac{1}{2}_0^1𝑑x\rho _a^{}(x)j_a^2(x)`$, which of course must converge exactly. The approximations $`D_n^j(a)`$ and $`D_n^\rho (a)`$ may be understood as time-dependent diffusion coefficients according to the Green-Kubo formula Eq. (5). According to their definitions, $`D_n^j(a)`$ enables us to look at contributions coming from $`j(x_n)`$ only, whereas $`D_n^\rho (a)`$ assesses the importance of contributions resulting from $`\rho _a^{}(x)`$. The rates of convergence of both approximations give an estimate of how important higher-order correlations are in the different parameter regions of the diffusion coefficient $`D(a)`$.

Similarly to Eq. (9), the next order $`D_1(a)`$ can easily be calculated analytically and reads

$$D_1(a)=\left\{\begin{array}{cc}\hfill (a2)/(2a),& 2<a1+\sqrt{3}\text{and}2+\sqrt{2}<a4\hfill \\ \hfill 3/23/a2/a^2,& 1+\sqrt{3}<a3\hfill \\ \hfill 1/2+3/a2/a^2,& 3<a2+\sqrt{2}\hfill \end{array}\right\}.$$
(12)
Further corrections up to order $`n=3`$ were obtained from computer simulations, that is, an ensemble of point particles was iterated numerically according to Eq. (1). All results are contained in Fig. 1 showing that this hierarchy of correlated random walks generates a self-affine structure which resembles, to some extent, the one of the well-known Koch curve. Fig. 1 (a) illustrates that this structure forms an important ingredience of the exact diffusion coefficient $`D(a)`$ thus explaining basic features of its fractality. Indeed, a suitable generalization of this approach in the limit of time to infinity leads to the formulation of $`D(a)`$ in terms of fractal generalized Takagi functions .

Fig. 1 (b) depicts the results for the series of $`D_n^\rho (a)`$ up to order $`n=2`$ as obtained from computer simulations. This figure illustrates that there exists a second source for an irregular structure related to the integration over the invariant density, as was explained above. In the Green-Kubo formula Eq.(5) both contributions are intimately coupled with each other via the integration over phase space.

We now perform a more detailed analysis to reveal the precise origin of the hierarchy of peaks in Fig. 1. For this purpose we redefine Eq. (10) as

$$D_n^j(a)=_0^1𝑑xj(x)J_n(x)\frac{1}{2}_0^1𝑑xj^2(x),n>0,$$
(13)
where we have introduced the jump velocity function

$$J_n(x):=\underset{k=0}{\overset{n}{}}j(x_k).$$
(14)
From Eq. (6) it follows $`J_n(x)=[x_{n+1}]`$, that is, this function gives the integer value of the displacement of a particle starting at some initial position $`x`$. In Fig. 2 we depict some representative results for $`J_1(x)`$ under variation of the control parameter $`a`$. Because of the symmetry of the map we restrict our considerations to $`0<x<0.5`$. Eq. (13) tells us that the product of this function with $`j(x)`$ determines the diffusion coefficient $`D(a)`$. The shaded bar in Fig. 2 marks the subinterval in which $`j(x)=1`$, whereas $`j(x)=0`$ otherwise, thus an integration over $`J_1(x)`$ on this subinterval yields the respective part of the diffusion coefficient. One can now relate the four diagrams (a) to (d) in Fig. 2 to the functional form of $`D_1^j(a)`$ in Fig. 1 (a) thus understanding where the large peak in $`D_1^j(a)`$ for $`2.732<a<3.414`$ comes from: For $`a<2.732`$, $`J_1(x)`$ does not change its structure and the interval where particles escape to other unit intervals increases monotonously, therefore $`D(a)`$ increases smoothly. However, starting from $`a=2.732`$ particles can jump for the first time to next nearest neighbors within two time steps, as is visible in $`J_1(x)`$ taking values of $`2`$ for $`x`$ close to $`0.5`$. Consequently, the slope of $`D(a)`$ increases drastically leading to the first large peak around $`a=3`$. Precisely at $`a=3`$, backscattering sets in meaning that particles starting around $`x=0.5`$ jump back to the original unit interval within two time steps, as is reminiscent in $`J_1(x)`$ in form of the region $`J_1(x)=0`$ for $`x`$ close to $`0.5`$. This leads to the negative slope in $`D(a)`$ above $`a=3`$. Finally, particles starting around $`x=0.5`$ jump to nearest neighbor unit intervals by staying there during the second time step instead of jumping back. This yields again a monotonously increasing $`D(a)`$ for $`a>3.414`$. Any higher-order peak for $`D_n^j(a),n>1,`$ follows from analogous arguments. Thus, the source of this type of fractality in the diffusion coefficient is clearly identified in terms of the topological instability of the function $`J_1(x)`$ under variation of the control parameter $`a`$. Indeed, this argument not only quantifies two previous heuristic interpretations of the structure of $`D(a)`$ as outlined in Refs. , it also explains why, on a very fine scale, there are still deviations between these results and the precise location of the extrema in $`D(a)`$, cp. to the “overhang” at $`a=3`$ as an example. The obvious reason is that contributions from the invariant density slightly modifying this structure are not taken into account.

Looking at the quantities $`J_n(x)`$ furthermore helps us to learn about the rates of convergence of the approximations $`D_n^j(a)`$ to $`D(a)`$ at fixed values of $`a`$, as is illustrated in Fig. 3 (a) to (d). Here we have numerically calculated $`J_n(x)`$ at $`a=3.8`$ for $`n=0,1,2,3`$. Again, the shaded bar indicates the region where $`j(x)=1`$ enabling $`J_n(x)`$ to contribute to the value of the diffusion coefficient according to Eq. (13). In fact, $`J_n(x)`$ may also be interpreted as the scattering function of an ensemble of particles starting from the unit interval, since it sensitively measures the final position to which a particle moves within $`n`$ time steps under variation of its initial position $`x`$. One can clearly see that, with larger $`n`$, $`J_n(x)`$ develops more and more discontinuities eventually leading to a highly singular and irregular function of $`x`$. Integration over further and further refinements of $`J_n(x)`$ determines the convergence of the series of $`D_n^j(a)`$ to a fixed value $`D_{\mathrm{}}^j(a)`$, cp. to Fig. 1 (a).

We remark that a suitable integration over the functions $`J_n(x)`$ leads to the definition of fractal so-called generalized Takagi functions, which can be calculated in terms of de Rham-type functional recursion relations . In a way, the integration over jump velocity functions such as the one shown in Fig. 3 (d) is similar to the integration over Cantor set structures leading to Devils’ staircase-type functions . Our results presented so far thus bridge the gap between understanding the coarse functional behavior of $`D(a)`$ on the basis of simple random walk approximations only, and analyzing its full fractal structure in terms of Takagi-like fractal forms, in combination with an integration over a complicated non-uniform invariant density. We now show that essentially the same line of argument can be successfully applied to more physical dynamical systems such as particle billiards.

## III Deterministic diffusion in billiards

The class of two-dimensional billiards we want to discuss here is described as follows: A point particle undergoes elastic collisions with obstacles of the same size and shape whose centers are fixed on a triangular lattice. There is no external field, thus the equations of motion are defined by the Hamiltonian $`H=mv^2/2`$ supplemented by geometric boundary conditions as induced by the scatterers. A standard example is the periodic Lorentz gas for which the scatterers consist of hard disks of radius $`R`$, see Fig. 4 . In the following we choose $`m=1`$, $`v=1`$, $`R=1`$, and as a control parameter we introduce the smallest inter disk distance $`w`$ such that the lattice spacing of the disks is $`2+w`$. $`w`$ is related to the number density $`n`$ of the disks by $`n(w)=2/[\sqrt{3}(2+w)^2]`$. At close packing $`w=0`$ the moving particle is trapped in a single triangular region formed between three disks, see Fig. 4, part (1). For $`0<w<w_{\mathrm{}}=4/\sqrt{3}2=0.3094`$, the particle can move across the entire lattice, but it cannot move collision-free for an infinite time. For $`w>w_{\mathrm{}}`$ the particle can move ballistic-like in form of arbitrarily far jumps between two collisions.

The diffusion coefficient for this particle billiard can be defined by the two-dimensional equivalent of the Einstein formula Eq.(4) reading

$$D(w)=\underset{t\mathrm{}}{lim}<(𝐱(t)𝐱(0))^2>/(4t),$$
(15)
where, again, the average is taken over the equilibrium distribution of particles with position coordinates $`𝐱(t)`$. It can be proven that in the regime of $`0<w<w_{\mathrm{}}`$ the parameter-dependent diffusion coefficient $`D(w)`$ exists . The full parameter-dependence of this function was discussed particularly in Ref.  showing that, on a fine scale, $`D(w)`$ is again an irregular function of the parameter $`w`$ similarly to the diffusion coefficient of the one-dmensional map $`D(a)`$ as discussed above. Whether $`D(w)`$ is indeed fractal, or maybe $`C^1`$ but not $`C^2`$ in contrast to the one-dimensional map discussed above, is currently an open question. The main issue we want to focus on in this section are quantitative approximations for the full parameter dependence of $`D(w)`$, and to check for the importance of memory effects. A first simple analytical approximation for the diffusion coefficient was derived by Machta and Zwanzig in Ref. . This solution was based on the assumption that diffusion can be treated as a Markovian hopping process between the triangular trapping regions indicated in Fig. 4. For random walks on two-dimensional isotropic lattices the diffusion coefficient then reads $`D=\mathrm{}^2/(4\tau )`$, where $`\mathrm{}=(2+w)/\sqrt{3}`$ is the smallest distance between two centers of the traps, and $`\tau ^1`$ is the average rate at which a particle leaves a trap. This rate can be calculated by the fraction of phase space volume being available for leaving the trap divided by the total phase space volume of the trap thus leading to the Machta-Zwanzig random walk approximation of the diffusion coefficient

$$D_{\mathrm{MZ}}(w)=\frac{w(2+w)^2}{\pi [\sqrt{3}(2+w)^22\pi ]}.$$
(16)
Indeed, this approximation is precisely the billiard analogue to the one-dimensional random walk diffusion coefficient for maps Eqs.(8),(9). Similarly, this approximation is asymptotically exact only for $`w0`$, as is shown in comparison to computer simulation results in Fig. 5 .

In Ref. , Eq. (16) was systematically improved by including higher-order correlations. Two basic approaches were presented both starting from the idea of Machta and Zwanzig of looking at diffusion in the Lorentz gas as a hopping process on a hexagonal lattice of “traps” with frequency $`\tau ^1`$. This picture was quantified by introducing a simple symbolic dynamics for a particle moving from trap to trap as indicated in Fig. 4. Let us follow a long trajectory of a particle starting with velocity $`𝐯`$ parallel to the $`y`$-axis, cp. part (1) in Fig. 4. For each visited trap we label the entrance through which the particle entered with $`z`$, the exit to the left of this entrance with $`l`$, and the one to the right with $`r`$, see part (2) in Fig. 4. Thus, a trajectory in the Lorentz gas is mapped onto words composed of the alphabet $`\{z,l,r\}`$. One can now associate transition probabilities to these symbol sequences reading, for a time interval of $`2\tau `$, $`p(z),p(l),p(r)`$. $`p(z)`$ corresponds to the probability of backscattering, whereas $`p(l)=p(r)=(1p(z))/2`$ indicates forward scattering. For a time interval of $`3\tau `$, we have nine symbol sequences each consisting of two symbols leading to the probabilities $`p(zz)`$, $`p(zl)`$, $`p(zr)`$, $`p(lz)`$, $`p(ll)`$, $`p(lr)`$, $`p(rz)`$, $`p(rl)`$, $`p(rr)`$, and so on. Beside this hierarchy of conditional probabilities defined on a symbolic dynamics there is a different type of probability that a particle leaving a trap jumps without any collision directly to the next nearest neighbor trap. As was shown in Ref. , Eq. (16) can be corrected by analytically including all these probabilities. Alternatively, lattice gas computer simulations were performed by using the probabilities as associated to the symbol sequences. The heuristic corrections of Eq. (16) led to a satisfactorily explanation of the overall behavior of $`D(w)`$ on a coarse scale, however, the convergence was not exact. The lattice gas simulations, on the other hand, were converging exactly, however, here a proper analytical expression for the diffusion coefficient approximations in terms of the associated probabilities was not available.

In analogy to the procedure as outlined for the one-dimensional map, that is, starting from a suitable Green-Kubo formula, we will now define a third approximation scheme which we expect to be generally applicable to diffusion in particle billiards. Compared to the two existing approaches mentioned above the advantage of the new method is two-fold, namely (1) that by using the set of symbolic probabilities the respective Green-Kubo formula can be evaluated according to an analytical scheme, and (2) that the resulting approximations converge exactly to the computer simulation results.

In the Appendix we prove that, starting from the Einstein formula Eq. (15), quite in analogy to the one-dimensional case a Green-Kubo formula can be derived which is defined for an ensemble of particles moving on the hexagonal lattice of traps depicted in Fig. 4. The result reads

$$D(w)=\frac{1}{4\tau }𝐣^2(𝐱_0)+\frac{1}{2\tau }\underset{n=1}{\overset{\mathrm{}}{}}𝐣(𝐱_0)𝐣(𝐱_n).$$
(17)
Here $`𝐣(𝐱_n)`$ defines jumps at the $`n`$th time step in terms of the lattice vectors $`\mathbf{}_{\alpha \beta \gamma \mathrm{}}`$ associated to the respective symbol sequence of the full trajectory on the hexagonal lattice, cp. to Fig. 4. Let us start with $`𝐣(𝐱_0)=\mathbf{}/\tau ,\mathbf{}:=(0,\mathrm{})^{}`$. The next jumps are then defined by $`𝐣(𝐱_1)=\mathbf{}_\alpha /\tau ,\alpha \{l,r,z\}`$, and so on. The averages indicated in Eq. (17) by the brackets are calculated by weighting the respective scalar products of lattice vectors with the corresponding conditional probablities $`p(\alpha \beta \gamma \mathrm{})`$. Note that Eq. (17) is the honeycomb lattice analogue to the Green-Kubo formula derived by Gaspard for the Poincaré-Birkhoff map of the periodic Lorentz gas . The Poincaré-Birkhoff version is very efficient for numerical computations, however, according to its construction it fails to reproduce the Machta-Zwanzig approximation Eq. (16). Consequently, it does not appear to be very suitable for diffusion coefficient approximations of low order. More details will be discussed elsewhere . We remark that, in terms of using a symbolic dynamics, there is also some link between Eq. (17) and respective diffusion coefficient formulas obtained from periodic orbit theory .

We now demonstrate how Eq. (17) can be used for systematic improvements of the diffusion coefficient on the lattice of traps by including dynamical correlations: As in case of one-dimensional maps we start by looking at the first term in Eq. (17) and cut off all higher order contributions. Obviously, the first term is again the random walk expression for the diffusion coefficient on the hexagonal lattice of traps which, by including the respective solution for the jump frequency $`\tau ^1`$, boils down to Eq. (16). For calculating higher-order corrections we now define the hierarchy of approximations

$$D_n(w)=\frac{l^2}{4\tau }+\frac{1}{2\tau }\underset{\alpha \beta \gamma \mathrm{}}{}p(\alpha \beta \gamma \mathrm{})\mathbf{}\mathbf{}(\alpha \beta \gamma \mathrm{}),n>0,$$
(18)
with $`D_0(w)`$ given by Eq. (16). To our knowledge yet there is no method available to analytically calculate the conditional probabilities $`p(\alpha \beta \gamma \mathrm{})`$. Our following evaluations are therefore based on the data presented in Ref.  as obtained from computer simulations. In terms of the formal probabilities it is now easy to calculate the solution for the first order approximation at time step $`2\tau `$ to

$$D_1(w)=D_0(w)+D_0(w)(13p(z)).$$
(19)
For a comparison of this formula with the simulation data $`D(w)`$ see Fig. 5. We remark that the corresponding solution in Ref.  as obtained from a heuristical correction of Eq. (16) reads $`D_{1,\mathrm{MZ}}(w)=D_0(w)3(1p(z))/2`$. Indeed, one can show that the new formula Eq.(19) is closer to $`D(w)`$ for large $`w`$, whereas for small $`w`$ the previous approximation is somewhat better. It is straightforward to calculate the two approximations of next highest order to

$$D_2(w)=D_1(w)+D_0(w)\left(2p(zz)+4p(lr)2p(ll)4p(lz)\right)$$
(20)
and to

$`D_3(w)`$ $`=`$ $`D_2(w)+D_0(w)[p(llr)+p(llz)+p(lrl)+p(lrr)+p(lzl)+p(lzz)+p(rll)+p(rlr)+p(rrl)`$ (23)
$`+p(rrz)+p(rzr)+p(rzz)+p(zll)+p(zlz)+p(zrr)+p(zrz)+p(zzl)+p(zzr)`$
$`2(p(lll)+p(lrz)+p(lzr)+p(rlz)+p(rrr)+p(rzl)+p(zlr)+p(zrl)+p(zzz))].`$
All results are shown in Fig. 5 demonstrating that the series of approximations defined by Eq. (18) converges quickly and everywhere to the simulation results. Our new scheme thus eliminates the deficiency of the semi-analytical approximation proposed in Ref.  that was based on heuristically correcting the Machta-Zwanzig approximation Eq. (16). By comparing this new scheme with the lattice gas simulations of the same reference, on the other hand, it turns out that the rate of convergence of the lattice gas approach is still a bit better. In any case, all these three methods unambiguously demonstrate that for achieving a complete understanding of the density-dependent diffusion coefficient in the periodic Lorentz gas it is unavoidable to take higher-order correlations, or the impact of memory effects, properly into account.

We finally remark that a very good low-order approximation for the diffusion coefficient can already be obtained by combining Eq. (19) with the probability of collisionless flights $`p_{\mathrm{cf}}(w)`$ mentioned above, that is, by taking into account the possibility of next nearest neighbor jumps. The correction of $`D_0(w)`$ as given by Eq. (16) according to collisionless flights only was already calculated in Ref.  and read $`D_{0,\mathrm{cf}}(w)=D_0(w)(1+2p_{\mathrm{cf}}(w))`$. Adding now the second term of Eq. (19) to this expression by just following the Green-Kubo scheme yields

$$D_{1,\mathrm{cf}}(w)=D_0(w)(2+2p_{\mathrm{cf}}(w)3p(z)).$$
(24)
This solution is also depicted in Fig. 5 and shows that this approximation indeed significantly improves Eq. (19) for large $`w`$ yielding a function that is qualitatively and quantitatively very close to $`D(w)`$. We know of no better approximation for $`D(w)`$ based on information such as $`p(z)`$ and $`p_{\mathrm{cf}}(w)`$ only. The succesful application of Eq. (24) suggests that collisionless flights form an important mechanism to understand the full diffusive dynamics of this billiard. However, somewhat surprisingly they are not explicitly contained neither in the Green-Kubo expansion Eq. (18) nor in the lattice gas simulations of Ref. . In both cases exact convergence is achieved by following the hierarchy of symbol sequence probabilities only, where collisionless flights are not apparent.

## IV Summary and conclusions

In this paper we suggested a general scheme of how to understand the structure of parameter-dependent deterministic diffusion coefficients in simple model systems. The important point was to find suitable Green-Kubo formulas to start with. We used the fact that the class of models we studied here was defined on periodic lattices, and respectively we discretized the dynamics of the moving particles according to these lattices. As two different examples we analyzed a simple one-dimensional map on the line as well as the periodic Lorentz gas. In both cases we recovered the respective well-known random walk formulas for the diffusion coefficient as zero-order approximations in our Green-Kubo approach. We then calculated higher-order terms according to Green-Kubo thus systematically including higher-order correlations. As much as possible this was performed analytically, alternatively in combination with data obtained from computer simulations. Our results provided clear evidence that a proper understanding of the parameter-dependent diffusion coefficients in both models can only be achieved by taking strong memory effects into account. In this respect our research appears to be somewhat related to the findings of long-time tails in the velocity autocorrelation function of simple model systems such as Lorentz gases, and to the existence of non-analyticities in the density expansion of transport coefficients, which are both consequences of strong correlations in the dynamics of the moving particles; see Ref.  for a nice review. In case of the one-dimensional map our approach enables a detailed understanding of the dynamical mechanism generating the most pronounced irregularities in the fractal diffusion coefficient. In case of the periodic Lorentz gas this scheme straightforwardly generalizes the Machta-Zwanzig random walk formula in terms of systematic higher-order approximations that converge exactly to the simulation results.

For one-dimensional maps the infinite time limit of this approach, though not the intermediate level as quantitatively discussed here, was already worked out in Refs.  leading to the definition of the diffusion coefficient in terms of fractal so-called generalized Takagi functions. For the periodic Lorentz gas an analogous generalization would be desirable as well. We furthermore remark that the approximation scheme as presented in this paper was already successfully applied (1) to the one-dimensional climbing sine map, where there is a complicated transition scenario between normal and anomalous diffusion , and (2) to the so-called flower shape billiard, where the hard disks of the Lorentz gas are replaced by obstacles of flower shape . In both cases the resulting parameter-dependent diffusion coefficients, as far as existent, are much more complicated functions than in the corresponding models discussed above, yet our scheme yields systematic explanations of the structure of these functions in terms of strong dynamical correlations. We thus expect that using Green-Kubo formulas this way provides a general access road to understanding deterministic diffusion in low-dimensional periodic arrays of scatterers, possibly also in view of experimental results on systems such as antidot lattices, ratchets, and Josephson junctions.

Acknowledgments
We thank T.Harayama for interesting discussions, L.Matyas for a careful reading of the manuscript, and Chr.Dellago for originally having generated the data set that was used for the Lorentz gas.

## A Green-Kubo formula on the hexagonal lattice

In this Section we derive Eq. (17) which yields the diffusion coefficient $`D(w)`$ for the periodic Lorentz gas via a generalization of the Machta-Zwanzig picture. That is, we look at diffusion as a higher-order Markov process on a hexagonal lattice of traps with lattice spacing $`\mathbf{}`$ where particles hop with frequency $`\tau ^1`$ from trap to trap, cp. to Fig. 4. The time $`t`$ is suitably rewritten in terms of the escape time $`\tau =(\pi /6w)(\sqrt{3}/2(2+w)^2\pi )`$  as $`t=n\tau ,n\mathrm{}`$. Let $`𝐱_n`$ be the position of the moving particle at time step $`n`$. We then write $`𝐱_n𝐗_n+\stackrel{~}{𝐱}_n`$, where $`𝐗_n`$ denotes the position of the trap on the hexagonal lattice in which the particle is situated at time step $`n`$. This vector can be expressed by a suitable combination of lattice vectors. For example, one may choose a sum of $`\mathbf{}_{\alpha \beta \gamma \mathrm{}}`$ as introduced in Section III, see Fig. 4, a more precise definition is not important here. Correspondingly, $`\stackrel{~}{𝐱}_n`$ is the distance from the nearest trap center to the actual position of the particle in the Wigner-Seitz cell.

The Einstein formula Eq. (15) then reads

$$D(w)=\underset{n\mathrm{}}{lim}(𝐗_n+𝚫\stackrel{~}{𝐱}_n)^2/(4n\tau )$$
(A1)
where $`𝚫\stackrel{~}{𝐱}_n:=\stackrel{\mathbf{~}}{𝒙}_n\stackrel{\mathbf{~}}{𝒙}_0`$. Multiplying out the nominator we get

$$<𝐗_n^2+2𝐗_n𝚫\stackrel{~}{𝐱}_n+𝚫\stackrel{~}{𝐱}_n^2>$$
(A2)
According to its definition the last term is bounded, $`𝚫\stackrel{~}{𝐱}_n^2<const`$. To the second term we apply the Hölder inequality  yielding

$$|<𝐗_n𝚫\stackrel{~}{𝐱}_n>|\sqrt{<|𝐗_n|^2>}\sqrt{<|𝚫\stackrel{~}{𝐱}_n|^2>}<const\sqrt{<|𝐗_n|^2>}.$$
(A3)
Consequently, in the limit of infinite time only the first term in Eq. (A2) contributes to the positive diffusion coefficient $`D(w)`$ of Eq. (A1) leading to

$$D(w)=\underset{n\mathrm{}}{lim}<𝐗_n^2>/(4n\tau ).$$
(A4)
Starting from this Einstein formula on the hexagonal lattice it is now straightforward to derive Eq. (17). $`𝐗_n^2=X_n^2+Y_n^2`$ tells us that essentially there are two one-dimensional parts, thus we are back to the respective derivation of Eq. (5) in Ref. . To make this precise, let us write

$$X_n=\underset{k=0}{\overset{n1}{}}j(x_k),$$
(A5)
which is in analogy to Eq. (14). However, here $`j(x_k):=X_{k+1}X_k`$ is strictly speaking no jump velocity, since it only denotes the distance a particle jumps within a time interval $`\tau `$, whereas $`\tau =1`$ in Eq. (6). Multiplying out the nominator of Eq. (A4) in terms of Eq. (A5) leads to an equation for $`D(w)`$ in form of velocity autocorrelation functions $`C(k,l):=<j(x_k)j(x_l)>=𝑑x𝑑y\rho ^{}(x,y)j(x_k)j(x_l)`$, where $`\rho ^{}(x,y)`$ is the equilibrium density of the periodic Lorentz gas. Translational invariance implies $`C(k,l)=C(0,lk)`$ which is easily shown by substitution combined with conservation of probability according to the Frobenius-Perron equation of the billiard. Summing up all contributions $`C(0,k)`$ obtained from the multiplication, doing the same for the component $`Y_n`$, and putting both results together yields Eq. (17).