text:
# Many-body symbolic dynamics of a classical oscillator chain

## 1 Introduction

In a historical numerical experiment Fermi Pasta and Ulam  studied the following non-linear model of a 1D crystal described by the Hamiltonian

$$H=\underset{i=1}{\overset{N}{}}\left(\frac{1}{2}m\dot{x}_i^2+V(x_ix_{i1}l)\right),V(x)=\frac{1}{2}m\omega ^2x^2+\frac{1}{4}kx^4.$$
(1)
Contrary to initial expectations of Fermi et al the so-called FPU model (1) behaved in a strong contrast to the ergodic theorem of statistical mechanics, even in quite strongly non-linear regime when one would expect fast relaxation to statistical canonical equilibrium and equipartition of energy from an arbitrary initial state. Instead, FPU model triggered the discovery of non-linear normal modes, the so-called solitons, and indirectly, the whole field of (computational) non-linear dynamics. Later, FPU chain has been studied in the spirit of KAM theory  and various estimates have been given for the critical strength of the dimensionless non-linearity parameter $`kl^2/m\omega ^2`$ required for the motion to become globally stochastic .

Recently, FPU-type models have been studied in the context of energy transport and Fourier heat law in 1D chains of particles . Quite unexpectedly, FPU-type models, namely Hamiltonians of the type (1) with non-linear inter-particle interaction $`V(x)`$ (and no on-site potential), turned out to be anomalous heat conductors, due to slow power law decay of transport (velocity-velocity or current-current) time correlation functions. These results, namely that correlations decay universally as $`C(t)t^{3/5}`$, and correspondingly that Kubo transport coefficient diverges as $`\kappa (L)L^{2/5}`$ as a function of the chain length $`L`$, have later been successfully explained in terms of hydrodynamic arguments and mode-mode coupling theory .

However, a perturbative-like approach such as mode-mode coupling theory should break down when the inter-particle potential has more than one stable position. Indeed, Giardina et al  have performed a series of numerical experiments on the particle chain with an oscillating inter-particle potential $`V(x)=1\mathrm{cos}(x)`$, as well as with the potential $`V(x)=x^2/2+x^4/4`$, and found normal heat conduction with a clean exponential decay of correlations. The key mechanism, which we believe produces non-KAM-like (almost) hyperbolic motion of such a high dimensional Hamiltonian system, are the hyperbolic saddles over which pairs of particles flip from one well to another. This motivated us to study in this paper some fundamental dynamical properties of the simplest version of such a hyperbolic chain with quartic double-well inter-particle potential $`V(x)`$. We present here some intriguing numerical results which suggest existence of a simple (but only approximate) Markov partition with very simple many-body symbolic dynamics of this particle chain. In addition, we find that the transition matrix is formally generated by a Hamiltonian of a certain quantum spin-1/2 chain.

Section 2 introduces the model. In section 3 we define (approximate Markov) partition of phase space and obtain numerical and analytical estimates for the volumes of various cells. Expression for the average time between subsequent jumps among the cells is also obtained. In section 4 we introduce Markovian description of our model with the transition matrix and the flux matrix and numerically check the accuracy of the Markovian property in two independent ways. Section 5 represents the core of the present paper. We numerically calculate the transition matrix and make a formal correspondence between our model and a ferromagnetic Heisenberg quantum spin-$`1/2`$ chain. Using this correspondence we derive various estimates and bounds for the decay rates of the time correlation functions of a class of piece-wise constant functions.

## 2 Model: Inverted Fermi-Pasta-Ulam chain

With the choice for units of mass, time and length of $`m`$, $`1/\omega `$ and $`l`$, respectively, and the new canonical coordinates $`q_i=x_iil`$ and $`p_i=m\dot{x_i}`$, Hamiltonian (1) can be brought to a dimensionless form

$$H=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2}+V(q_iq_{i1}),V(x)=\alpha \frac{x^2}{2}+\frac{x^4}{4}+\frac{\alpha ^2}{4},$$
(2)
where $`\alpha =m\omega ^2/kl^2`$ is the only dimensionless parameter left. In order for the potential minima to be at level 0 (for $`\alpha >0`$) we have added a constant term $`\alpha ^2/4`$ to the potential $`V(x)`$. We can get all different orbits of the system described by the Hamilton function (2) just by varying one parameter $`\alpha `$. At first it would seem that the energy $`E`$ is also a free parameter. But by fixing all the three basic units we have also provided a natural unit for the energy density $`\epsilon =E/N`$. If we scale the dimensionless quantities by a factor $`\xi `$, namely $`(q_i,p_i,\alpha )(\xi q_i,\xi ^2p_i,\xi ^2\alpha )`$, the Hamiltonian is simply multiplied by a factor $`\xi ^4`$. Invariant parameter in this transformation is therefore $`\lambda =\epsilon /\alpha ^2`$. All different nonequivalent Hamiltonians can be obtained just by varying one parameter, namely $`\lambda `$, or equivalently $`\alpha `$ at fixed $`\epsilon `$. Energy density can therefore be fixed without any loss of generality. From now on $`\epsilon =E/N=1`$ is assumed, as well as periodic boundary conditions $`(q_{N+1},p_{N+1})(q_1,p_1)`$. Phase space point will be denoted by $`\mathrm{x}\mathbf{=}\mathbf{(}\mathrm{q}\mathbf{,}\mathrm{p}\mathbf{)}`$. For negative values of the parameter $`\alpha `$ the system has a form of the standard $`\beta `$-FPU model (just in a different parametrization), but in this paper we focus on the positive values of the parameter $`\alpha `$ in which case we call the system Inverted FPU model (IFPU). IFPU model is translationally invariant, therefore besides the energy $`\epsilon `$ there is also a second (trivial) constant of motion, namely the total momentum $`P=_{i=1}^Np_i`$. In all numerical calculations the total momentum has been fixed to $`P=0`$.

For the numerical integration of equations of motion a fourth order symplectic algorithm of Ref.  has been used. This integrator has been checked to be the optimal choice (at least for the model studied here, and at relative accuracy $`10^510^8`$) by making a careful comparison with a number of other symplectic and Runge-Kutta methods.

## 3 Phase space partition and statistical dynamics

If the barrier $`\mathrm{\Delta }U=\alpha ^2/4`$ between the potential minima is sufficiently high the differences $`(q_{i+1}(t)q_i(t))`$ will spend most of the time either around the left or the right minimum, with quite infrequent jumps between the wells. As the motion around the minimum is approximately harmonic the more interesting part will be the transitions between the left and the right potential well, for each pair of neighboring particles, going over the hyperbolic saddle at $`q_{i+1}q_i=0`$. Understanding the motion in phase space of a high-dimensional system is generally very difficult. Some symbolic description of an orbit $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$ is therefore highly desirable. In IFPU model with high potential barrier the choice of a partition of phase space is almost obvious:

For each instant of time $`t`$ we will be interested only in a single binary bit of information $`a_i\{0,1\}`$ for each pair of neighboring particles $`(i,i+1)`$, namely if the neighbors are in the right or left well, $`(q_{i+1}q_i)>0`$ or $`<0`$, we write $`a_i=1`$ or $`a_i=0`$, respectively. Binary positions for the whole chain can be compactly encoded in a binary integer (signature)

$$S=(a_N\mathrm{}a_2a_1)_2=\underset{i=1}{\overset{N}{}}a_i2^{i1}.$$
(3)
Thus we have defined a partition of $`2N`$-dimensional phase space cut by $`N`$ hyper-planes $`q_{i+1}=q_i`$ into $`2^N`$ cells labeled by signatures $`S\{0,1,\mathrm{}2^N1\}`$. By a slight abuse of language, we will use a term signature $`S`$ also to refer to a phase space cell denoted by $`S`$. If the transitions between various signatures are rare, which is obviously the case if the potential barrier $`\mathrm{\Delta }U`$ is high, it seems to be meaningful to concentrate on the (statistical) dynamics of transitions between the signatures. For the periodic boundary conditions the coordinate differences must fulfill the constraint $`_{i=1}^N(q_{i+1}q_i)0`$ which in the limit of infinitely high potential barrier and even $`N`$ translates into the condition $`_{i=1}^Na_i=N/2`$. For odd $`N`$ the number of pairs in the left and the right well should differ by 1. This presents only unnecessary technical complications and from now on we will assume $`N`$ to be even. If the barrier is finite, signatures with different number of pairs in the left and right well are possible to visit. We will call signature $`S`$ to be of order $`i`$ if $`_{i=1}^Na_i=N/2\pm i`$. The number of different signatures $`M_i`$ of order $`i`$ is

$$M_i=\left(\genfrac{}{}{0pt}{}{N}{N/2i}\right)=\frac{N!}{(N/2i)!(N/2+i)!}.$$
(4)
For sufficiently high barrier the system will spend most of its time in signatures of order $`0`$. We will specifically concentrate on signatures of order 0 and treat signatures of higher (mostly 1st) order only as ’tunnels’ for transitions between different signatures of order 0. Figure 2 shows an example of a transition between two order 0 signatures via an intermediate short lived signature of order 1. First we have to know for which values of parameter $`\alpha `$ and size $`N`$ will the description by the signatures of order 0 only be adequate, i.e. will the relative time spent in higher order signatures be negligible.

### 3.1 Estimating the fractional volumes of higher order signatures

We will now make a rough theoretical estimate for the time spent in signatures of different orders. For an ergodic system this time is proportional to the measure (volume) of the phase space covered with the signatures in question. Let us designate by $`t_i`$ the time spent in signatures of order $`i`$. We have

$$t_i(E)_{\mathrm{order}i}\frac{dS}{|H|}=\frac{d\mathrm{\Gamma }_i(E)}{dE},$$
(5)
where the region of integration is the part of the energy surface $`H=E`$ intersecting with the full set of signatures of order $`i`$ and $`\mathrm{\Gamma }_i(E)`$ is the total volume of the phase space in the signatures of order $`i`$ and with energy less than $`E`$. First, let us make an estimate for $`\mathrm{\Gamma }_0(E)`$. Potential around both minima is to the lowest order harmonic, therefore we can roughly say $`\mathrm{\Gamma }_0(E)\mathrm{\Gamma }(E)^N`$, where $`\mathrm{\Gamma }(E)`$ is the volume for one harmonic oscillator, that is $`\mathrm{\Gamma }(E)=2\pi E/\sqrt{\alpha }`$. For the signatures of order $`i`$ the argument is very similar. In this case there are $`2i`$ particles more in one potential well than in the other. Because the trivial condition $`_{i=1}^N(q_{i+1}q_i)0`$ must still be satisfied, the equilibrium positions of pairs will not be at $`\pm \sqrt{\alpha }`$ any more but will be shifted for $`\mathrm{\Delta }q`$ in order to keep the center of mass of all pairs at 0. This gives for the shift $`\mathrm{\Delta }q=2i\sqrt{\alpha }/N`$ and contributes $`\epsilon _i=\alpha \mathrm{\Delta }q^2=4i^2\alpha ^2/N^2`$ to the energy density, provided the cubic part of the potential is negligible, i.e. $`2i/N1`$. Harmonic approximation for the volume $`\mathrm{\Gamma }_i(E)`$ can therefore still be used, but now at energy $`EN\epsilon _i`$. Since we are interested only in the dependences for large $`N`$ we can omit (’cancel’) energy derivatives coming from (5), so that we have

$$\frac{t_i}{t_0}\frac{\mathrm{\Gamma }(\epsilon \epsilon _i)^N}{\mathrm{\Gamma }(\epsilon )^N}=\left(1\frac{\epsilon _i}{\epsilon }\right)^N\mathrm{exp}\left(4i^2\frac{\alpha ^2}{N\epsilon }\right).$$
(6)
We have fixed $`\epsilon 1`$ so that the relative fraction of signatures of order 1, $`t_1/t_0`$, falls exponentially in $`\alpha ^2/N`$. We must stress that the exponential dependence is expected only when $`4\alpha ^2/N^21`$ and that the overall prefactor in front of an exponential function could still depend on $`N`$. This is nicely confirmed by the numerical data in figure 3.

Grouping the signatures by their order is meaningful since members of each group have the same phase space volume and the same minimal potential energy. The minimal potential energy of order 0 signatures is $`\epsilon _0=0`$, of order 1 is $`\epsilon _1=4\alpha ^2/N^2`$ and of order $`i`$ is $`\epsilon _ii^2\epsilon _1`$. The fraction of time spent in higher order signatures is just a power of relative time spent in order 1 signatures, $`(t_i/t_0)(t_1/t_0)^{i^2}`$ and this can be reduced by increasing $`\alpha `$.

Of course, for the transitions to be possible at all, the barrier height must be smaller than the total energy $`E`$. This yields the condition

$$\alpha ^2/4<N$$
(7)
which is just the opposite of the condition to keep the $`t_1/t_0`$ small. Nevertheless, we can still reduce the fraction of higher order states to around $`\mathrm{exp}(16)10^6`$, which is small. But the condition (7) can be problematic for small $`N`$. Indeed, for the smallest nontrivial case of $`N=4`$ we have numerically found that there are no transitions between signatures for $`\alpha 2.8`$. For smaller $`\alpha `$ we do find transitions, but there the system exhibits non-ergodic behaviour. Condition (7), although fulfilled, seems to be too weak and is still keeping the phase space cut into non-connected parts. The smallest system worth studying regarding transitions is therefore $`N=6`$.

In all numerical work $`\alpha `$ and $`N`$ have been chosen so as to keep the fractional volume of higher order signatures small. We have therefore studied the dynamics on a set of $`M_0`$ signatures of order 0 exhibiting transitions thru short-lived order 1 signatures as shown in figure 2. The most important physical scale that characterizes such a transport is the average time $`\tau `$ between subsequent transitions.

### 3.2 Ergodicity and average transition time $`\tau `$ between signatures of order 0

Though the main body of this paper is concerned with the transport and decay of time-correlations, i.e. with the system’s mixing property, we should first mention, that we have also checked the weaker dynamical property of ergodicity directly. We have employed a method of Robnik et al  of comparing the rate of visiting of different phase space cells (in out case, order 0 signatures) with the rate for a fully random dynamics. The results turned to be fully consistent with (uniformly) ergodic behaviour of our IFPU model in the high-barrier regime discussed above.

Now we define the time scale $`\tau `$ to be an average time from the point when the orbit enters a certain order 0 signature $`S`$ to the point where the same orbit enters the next order 0 signature $`S^{}`$ ($`S^{}S`$) concluding the transition from $`S`$ to $`S^{}`$. The average is taken over many different orbits with microcanonically distributed initial conditions, or over one very long orbit if we assume ergodicity. We obtain an approximate functional form for the dependence of $`\tau `$ on $`N`$ and $`\alpha `$ by similar arguments as for the relative volume of higher order signatures. The probability for a jump will be estimated by the ratio between the phase space volume of the set of states just before a jump $`\mathrm{\Gamma }_\mathrm{t}`$ and the volume of the set of equilibrium states $`\mathrm{\Gamma }_{\mathrm{eq}}`$. The approximate volume $`\mathrm{\Gamma }_{\mathrm{eq}}(E)`$ of the phase space of equilibrium states has been derived in the previous subsection and is just $`\mathrm{\Gamma }_{\mathrm{eq}}(E)\mathrm{\Gamma }(E)^N=(2\pi /\sqrt{\alpha })^NE^N`$. In the state just before the transitions there should be more than $`\mathrm{\Delta }U`$ of energy in one pair of particles. This energy will be denoted by $`a\alpha ^2`$, where $`a`$ is for now some unspecified constant. Certainly $`a`$ must be greater than $`1/4`$, but because our oscillators are coupled we allow for the possibility that the energy of a pair must actually be bigger than the barrier height. Phase volume $`\mathrm{\Gamma }_\mathrm{t}(E)`$ for a state with one pair having the energy $`a\alpha ^2`$ and the other $`N1`$ pairs residing in the minima is

$`\mathrm{\Gamma }_\mathrm{t}(E)={\displaystyle _0^{Ea\alpha ^2}}𝑑\epsilon \epsilon ^{N1}\left({\displaystyle \frac{2\pi }{\sqrt{\alpha }}}\right)^{N1}(N1)(\mathrm{\Gamma }(E)\mathrm{\Gamma }(\epsilon ))`$
$`\left({\displaystyle \frac{2\pi }{\sqrt{\alpha }}}\right)^Na\alpha ^2(Ea\alpha ^2)^{N1}.`$ (8)
Probability for a jump is proportional to the quotient of derivatives $`d\mathrm{\Gamma }_\mathrm{t}(E)/d\mathrm{\Gamma }_{\mathrm{eq}}(E)`$. Reciprocal average time between transitions $`1/\tau (N,\alpha )`$ is proportional to the above quotient and to the frequency of oscillation around potential minimum. Using also $`E=N`$, we finally get

$$\tau (N,\alpha )\frac{A}{\sqrt{\alpha }}\frac{N}{a\alpha ^2}\frac{1}{(1a\alpha ^2/N)^N},$$
(9)
where $`A`$ is some overall numerical factor. We have determined numerical parameters $`a`$ and $`A`$ by fitting the numerical data for different lengths $`N`$ in the range $`N=6,\mathrm{},30`$. Both parameters do depend on $`N`$ for small chains, but are, of course, independent of $`\alpha `$. The value of $`A`$ is between $`0.02`$ and $`0.05`$, whether the $`a`$ has dependence $`a=1/4+𝒪(1/N)`$. In the limit of large chains $`a`$ has the value $`1/4`$ predicted by our simple physical picture. This picture can also be confirmed by looking at the local energy density just before the jump. Indeed, local energy of a pair is exactly $`a\alpha ^2`$ for all chain sizes. So, the parameter $`a`$ is not just some fitting parameter in equation (9) but it has a clear physical interpretation. In figure 4 we depict the estimate of $`\tau `$ (9) together with the numerical data. Worth noticing is also an interesting scaling law, although only approximate. Global dependence of $`\tau (N,\alpha )`$ is such that we can write approximately

$$\tau (N,\alpha )\frac{f(4N/\alpha ^2)}{\sqrt{\alpha }}\mathrm{exp}(\alpha ^2/4),$$
(10)
where $`f(x)`$ is some parameter-free function. The shape of this function can be seen in figure 5.

We have also numerically computed statistical distribution of transition times $`\tau _n`$, i.e. the time intervals between entering different subsequent signatures of order 0 such that $`\tau =\tau _n`$, for a fixed $`N`$ and $`\alpha `$. This distribution is shown in figure 6 and is clearly seen to be exponential.

There is a discrepancy only for extremely small times. This is a consequence of the duration of intermediate unstable order 1 state that occurs between order 0 states (see figure 2). The fraction of order 1 states was in this case $`t_1/t_0=0.0013`$ which is the same as the time scale of the discrepancy in the figure. Our distribution should be a convolution of distributions of $`t_0`$ and of $`t_1`$. But the precise explanation of this discrepancy is the following: we obtain too many of “short” jumps because of the way we measure the time when the system has arrived into a new order 0 signature. We mark the system as arriving into a new signature exactly at the top of the barrier, when the difference $`(q_{i+1}q_i)`$ changes the sign and not when it relaxes to the bottom of the well. This allows for some ‘fake transitions’ in the regions of phase space where two or more 0 order signatures almost touch with a saddle-like region in between (for instance near coordinate origin $`q_i0`$). There it is possible that the system, which although it gets over the barrier, subsequently relaxes into some nearby well because it does not have the ‘right’ momentum vector.

## 4 Markov process and its transition matrix

We decided to describe dynamics with a sequence of binary signatures instead of a full orbit $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$ with the hope of some simplification. We are hence looking for some simple statistical description of the transitions between signatures. Stochastic system is fully described by giving conditional transition probabilities $`P(S,S^{},S^{\prime \prime },\mathrm{};t,t^{},\mathrm{})`$. The simplest kind of a stochastic process is called Markov process. For a Markov process the probabilities $`P(S,S^{};t)`$ contain entire information about the system. We will shortly write a matrix $`𝐏(t)`$ instead of $`P(S,S^{};t)`$. For now let us define $`P(S,S^{};t)`$ on the full space of $`2^N`$ signatures although we will later be interested only in transition probabilities among $`M_0`$ signatures of order 0. The transition matrix $`𝐏`$ must satisfy the condition

$$𝐏(t+t^{})=𝐏(t)𝐏(t^{}).$$
(11)
The matrix $`𝐏(t)`$ contains therefore some redundant information. We can write instead

$$𝐏(t)=\mathrm{exp}(𝐅t),$$
(12)
where for a Markov process the matrix $`𝐅`$ is time independent. The matrix element $`F(S,S^{})`$ is a probability flux from the signature $`S`$ to the signature $`S^{}`$, and can be defined by the limit

$$𝐅=\underset{t0}{lim}\frac{𝐏(t)𝐈}{t}.$$
(13)
The conservation of probability imposes the condition $`_S^{}F(S,S^{})=0`$ for the matrix elements.

By a statistical description, using $`𝐏(t)`$ instead of $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$, we lose information about the details within a single cell of a partition. In other words, we study the dynamics on the space of functions that are constant over one phase space cell. This can be formalized by introducing characteristic functions $`B_S(\mathrm{x}\mathbf{)}`$ over signatures $`S`$

$$B_S(\mathrm{x}\mathbf{)}\mathbf{=}\mathbf{\{}\begin{array}{cc}1\hfill & \mathrm{x}\mathbf{}\mathrm{S}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\mathbf{.}$$
(14)
Characteristic functions $`B_S`$ span the space of all observables which are piece-wise constant over the signatures. Every piece-wise constant observable $`W(\mathrm{x}\mathbf{)}`$ can be expanded over the base functions $`B_S`$ as $`\mathrm{W}(\mathrm{x}\mathbf{)}\mathbf{=}\mathbf{}_\mathrm{S}\mathrm{w}_\mathrm{S}\mathrm{B}_\mathrm{S}\mathbf{(}\mathrm{x}\mathbf{)}`$ with the expansion coefficients $`w_S`$ given by

$$w_S(t)=\frac{B_S\mathrm{W}(t)_\mathrm{E}}{B_S_\mathrm{E}},$$
(15)
where the brackets $`_\mathrm{E}`$ denote a microcanonical phase space average over the energy surface $`H=E`$. The matrix elements of a transition matrix $`𝐏`$ are nothing but the correlation functions between the characteristic functions

$$P(S,S^{};t)=B_S^{}(t)B_S_\mathrm{E}/B_S_\mathrm{E}.$$
(16)
As the transition matrix propagates the probabilities between the signatures we can write the time dependent vector $`\mathrm{w}(t)=(w_S(t),w_S^{}(t),\mathrm{})`$ as

$$\mathrm{w}(t)=𝐏(t)\mathrm{w}(0)=\mathrm{exp}(𝐅t)\mathrm{w}(0).$$
(17)
Similarly, the autocorrelation function of the observable $`W`$ is

$$\mathrm{W}(t)\mathrm{W}(0)_\mathrm{E}=\mathrm{w}^T\mathrm{exp}(𝐅t)\mathrm{w}.$$
(18)
We can see that the behaviour of correlation functions is determined by the spectrum of the (Markov) propagator $`\mathrm{exp}(𝐅t)`$.

From now on we assume that the fraction of higher order signatures is very small $`\eta t_1/t_01`$ and we restrict the matrices $`𝐏(t)`$ and $`𝐅`$ on the $`M_0`$-dimensional subspace of signatures of order 0, which essentially contain all non-vanishing matrix elements. Numerically we calculated the transition matrix $`𝐏(t)`$ from one very long orbit as an average

$$P(S,S^{};t)=\frac{\delta _{S^{},S(t^{}+t)}\delta _{S,S(t^{})}_\mathrm{t}^{}}{\delta _{S,S(t^{})}_\mathrm{t}^{}},$$
(19)
where the brackets $`_\mathrm{t}^{}`$ denote a time average over the orbit and $`\delta _{S,S(t)}`$ has value 1 if the orbit is in a signature $`S`$ at time $`t`$ and 0 otherwise. For small $`t`$ the average can be performed, and by using expression (13), yields for $`F(S,S^{})`$, $`(SS^{})`$

$$F(S,S^{})=\frac{n(S,S^{})}{t_S},$$
(20)
where $`n(S,S^{})`$ is the number of direct transitions between the order 0 signatures $`S`$ and $`S^{}`$ and $`t_S`$ is the total time spent in the signature $`S`$. Here a small comment regarding higher order signatures is in order. We saw in figure 2 that between each order 0 states there is a short intermediate order 1 state of duration $`\tau t_1/t_0`$. So there are no direct transitions between order 0 states and all fluxes $`F(S,S^{})`$ (20) among $`M_0`$ order 0 signatures are strictly zero. As we have decided to study only (indirect) transitions between order 0 states (for the parameter values where the fraction of higher order states is negligible) we must somehow circumvent this difficulty. One solution is to calculate the derivative of the transition matrix $`𝐏(t)`$ in equation (13) not at time $`t=0`$ but at the time $`t`$ with $`\tau t_1/t_0t\tau `$. In this case the number $`n(S,S^{})`$ in (20) is the number of jumps between order 0 states $`S`$ and $`S^{}`$ with one short intermediate order 1 state. Still more simple and convenient solution is to replace the orbit $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$ by a sequence of times $`t_i`$ when the orbit hits the boundary of a cell of order 0 coming from outside. Orbit $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$ is therefore replaced by a sequence of pairs $`\mathrm{},(S_i,t_i),(S_{i+1},t_{i+1}),\mathrm{}`$ which tell you that at time $`t_i`$ orbit $`\mathrm{x}\mathbf{(}\mathrm{t}\mathbf{)}`$ came to the order 0 signature $`S_i`$, then until time $`t_{i+1}`$ was in signature $`S_i`$ or any higher order signature and at time $`t_{i+1}`$ arrived into the next order 0 signature $`S_{i+1}`$ ($`S_{i+1}S_i`$), etc. Then we define $`n(S,S^{})`$ as the number of subsequent pairs $`(S,S^{})`$ in the sequence $`\{S_i\}`$, and $`t_S=_i^{S_i=S}(t_{i+1}t_i)`$, hence the flux matrix is calculated as

$$F(S,S^{})=\frac{_i^{(S,S^{})=(S_i,S_{i+1})}1}{_i^{S_i=S}(t_{i+1}t_i)}.$$
(21)
This definition is used in the numerical calculation of this paper.

### 4.1 Tests for the Markov process

Before describing the process of transitions as Markovian we must check if this is at all permissible. Checking the rigorous conditions for the partition to signatures to be Markovian, see e.g. , is a considerable mathematical problem, at least to our judgement. On the other hand, heuristic arguments and numerical results support the idea, that for sufficiently large transition times $`\tau `$ the process will indeed be Markovian. The system has positive Lyapunov exponents.<sup>3</sup><sup>3</sup>3We have numerically checked that the distribution of Lyapunov exponents is approximately linear, as is typically the case for sufficiently chaotic systems . Vaguely this means that the system quickly “forgets” its history. If the average time between the jumps $`\tau `$ is longer than the Lyapunov time, we can expect transitions to be Markovian. However, the most straightforward test is to check explicitly the composition formula (11).

Condition $`𝐏(2t)=𝐏^2(t)`$ is trivially fulfilled for small and for large times $`t`$, $`t\tau `$ and $`t\tau `$, respectively. We checked it for time $`t=\tau `$ by calculating the quantity

$$\sigma =\frac{𝐏(2\tau )𝐏^2(\tau )_\mathrm{E}}{𝐏(2\tau )_\mathrm{E}},$$
(22)
where $`𝐀_\mathrm{E}=(_{i,j}|A_{i,j}|^2)^{1/2}`$ is the Euclidean norm of a matrix $`𝐀`$ (all results have also been re-checked by using the spectral norm giving almost identical results). As the matrix $`𝐏(t)`$ will be calculated numerically from one very long but finite orbit it will have statistical error $`\sigma _\mathrm{T}`$ due to the finite number of simulated transitions. In addition, $`\sigma `$ will also have a contribution $`\sigma _\mathrm{M}`$ from a systematic error because the process may not be precisely Markovian. We have assumed $`\sigma ^2=\sigma _\mathrm{M}^2+\sigma _\mathrm{T}^2`$ and made an estimate for the statistical error $`\sigma _\mathrm{T}`$ as a square root of a number of average transitions per matrix element $`n`$. It is

$$\sigma _\mathrm{T}=c\frac{1}{\sqrt{n}}=c\frac{N}{2}\sqrt{\frac{M_0}{N_0}},$$
(23)
where $`N_0`$ is the total number of all transitions between the signatures of order 0 and $`c`$ is some unspecified numerical constant which has been determined from the numerical data. Its value turned out to be $`c0.7`$. Results for the dependence of a systematic error $`\sigma _\mathrm{M}`$ are shown in figure 7.

The error $`\sigma _\mathrm{M}(\tau )`$ goes to zero with increasing $`\tau `$, as expected. Even more interesting is the fact that the error also seems to decrease with increasing $`N`$ at constant $`\tau `$. This is interesting as it means that with increasing $`N`$ we do not have to keep the fraction of higher order states small in order for the process to stay Markovian. <sup>4</sup><sup>4</sup>4In such a case one could bring higher order signatures back and describe a complete hierarchy of transitions between signatures of various orders as a Markov process on $`2^N`$ dimensional space.

Whether a given system can be described as a Markov process can be checked also by calculating the transition matrix $`𝐏(t)`$ and comparing it with the exponential formula (12) as a function of time. $`𝐏(t)`$ can be calculated numerically directly from definition, for instance by using the formula (19). This definition would have a meaning also if the transitions were not Markovian. We then compare the correlation functions between the characteristic functions (16) and the matrix elements of the propagator $`S|\mathrm{exp}(𝐅t)|S^{}`$. These calculations should agree only if the flux matrix $`𝐅(t)=\mathrm{log}𝐏(t)/t`$ were time independent and this would signal that the transitions are described by a Markov process. Numerical results for this test are shown in figure 8.

Correlation functions calculated in both ways have been computed for the characteristic function $`B_S|S`$ on the signature $`S=54`$ and for the macroscopic observable with the vector $`w_S=S\text{ mod }2`$, i.e. characteristic function on a union of half of all signatures with the first bit $`a_1=1`$, denoted by $`|s1`$. It can be seen from the figure that the autocorrelation function of the observable $`|s1`$ begins to fall as $`\mathrm{exp}(\lambda _1t)`$, where $`\lambda _1`$ is the biggest nontrivial eigenvalue of $`𝐅`$, very early. This is a consequence of relatively large overlap between the macroscopic state $`|s1`$ and the eigenstate $`|v1`$ corresponding to the eigenvalue $`\lambda _1`$. We will show in subsection 10 that $`|s1|v_1|^21/N`$ while $`|S|v_1|^21/M_0`$, so $`|s1|v_1|^2|S|v_1|^2`$ (for $`N1`$).

## 5 The form and the spectral properties of the flux matrix

Something can be said immediately about the spectrum of $`𝐅`$. Conservation of probability implies all eigenvalues $`\lambda _i`$ to be smaller or equal to zero. Further, because the system is believed to be ergodic, there should be a non-degenerate constant eigenvector with the corresponding eigenvalue $`\lambda _0=0`$. This is all the consequence of the Perron-Frobenius form of the matrix $`𝐏`$. For finite $`N`$ the spectrum of the flux matrix is discrete and the correlation functions will fall asymptotically as $`\mathrm{exp}(\lambda _1t)`$ where $`\lambda _1`$ is the largest nonzero eigenvalue. An interesting question is what happens in the thermodynamic limit $`N\mathrm{}`$. If there is a spectral gap we will have an exponential decay, otherwise some anomalous behaviour can occur.

### 5.1 Numerical calculation of the flux matrix

Non-vanishing matrix elements of $`𝐅`$ are expected only between the signatures that differ just by two bits provided that $`t_1/t_01`$. In other words, two bits $`a_i=1`$ and $`a_j=0`$ only switch their position. Such a transition will be called a (bit) jump of length $`d`$, where $`d=|ji|`$ is the distance between the two bits involved. Example of a such transition has been shown in figure 2. We have numerically calculated complete matrices $`𝐅`$ for sizes $`N=6,8,10`$ and $`12`$ and different $`\tau `$ (or $`\alpha `$), in order to check the above hypothesis and to obtain some knowledge about the sizes and the structure of matrix elements. The only matrix elements which are really non-vanishing are those of the bit jumps, and what is more, the value of the matrix element depends only on the distance $`d`$ of the jump and not on the particular signatures involved. This practically means that in addition to the diagonal elements we have only $`N/2`$ different matrix elements in the flux matrix $`𝐅`$. This number must be compared to the total number of matrix elements $`M_0^2`$ which grows exponentially with $`N`$. We have checked this also for the matrices $`𝐅`$ of sizes up to $`N=30`$ but in this case we have compared only the jumps of length $`d`$ with different number of bits 1 on the sites between the two jump sites. We confirmed that the probability fluxes of jumps depend only on the length $`d`$. It is therefore meaningful to introduce dimensionless coefficients $`c_d`$ which are ratios between the matrix elements of a jump of length $`d`$ and a jump of length $`d=1`$. Therefore $`c_11`$, by definition, and the rest of $`N/21`$ coefficients $`c_d`$ together with the average transition time $`\tau `$ is all that we need in order to specify the flux matrix $`𝐅`$ completely. We can write

$$𝐅=\frac{1}{\tau }\frac{2N2}{N^2(c_{N/2}/2+_{d=1}^{N/21}c_d)}𝐂,$$
(24)
where the matrix $`𝐂`$ is just a matrix involving the coefficients $`c_d`$ only. Numerical values of $`c_d`$ for different $`N`$ and $`\tau `$ are listed in table 1. It can be seen that for not too big $`\tau `$ (e.g. $`\tau =11`$ and $`N=10`$) the coefficients $`c_d`$ decrease monotonically with the distance $`d`$. For bigger $`\tau `$ (e.g. $`\tau =621`$ and $`N=10`$), this dependence ceases to be monotonic but becomes slightly well-shaped with a minimum at around $`dN/4`$.

There also seems to be a general trend that with increasing $`N`$, at a constant fraction $`\eta t_1/t_0`$, the coefficients $`c_d`$ all approach 1. This can be explained by the following heuristic argument. To keep the fraction $`t_1/t_0`$ of higher order states constant with increasing $`N`$ we must increase the barrier height $`\mathrm{\Delta }U`$. As a consequence, the transition time $`\tau `$ increases and correspondingly also the average time $`\tau t_1/t_0`$ spent in the non-equilibrium intermediate state of order 1 in between the states of order 0. At a constant $`t_1/t_0`$ the transition time $`\tau `$ depends approximately exponentially on the chain length $`N`$ (see figure 4). On the other hand, the relaxation time for the non-equilibrium energy distribution at the moment of transition grows perhaps only linearly with $`N`$ (inversely proportional to the sound speed), but certainly slower than exponentially. For sufficiently large $`N`$ the relaxation time will therefore be smaller than the time spent in the intermediate order 1 cell. Initial locally non-equilibrium distribution at the beginning of a jump will therefore relax before the transition to the new order 0 state, which will therefore all be equally probable. Because of this we can assume that the coefficients $`c_d`$ will not depend on $`d`$ in the thermodynamic limit. This hypothesis still needs further verification as the limit of high $`N`$ could not be tested due to the fast growth of the dimension $`M_0`$ with $`N`$. Despite of that, the numerical results are consistent and point in the right direction (see table 1). We will show now, that our flux matrix has a specially appealing form which can be described using a formal connection to a certain Hamiltonian of a quantum spin-$`1/2`$ chain. Based on this we will also deduce the essential spectral properties of the flux matrix.

### 5.2 Correspondence between the Markov process and the quantum spin chain

We can make the announced correspondence by connecting the eigenvalue problem for the flux matrix $`𝐅`$ with the eigenvalue problem for the quantum Hamilton operator $`\widehat{\mathrm{H}}`$ over $`M_0`$-dimensional Hilbert space with the same matrix elements as $`𝐅`$. This correspondence would generally be of no particular use since the Hamiltonian $`\mathrm{H}`$ does not represent any simple quantum system. In our example of IFPU chain this is not the case, since $`𝐅`$ has a particularly nice and simple form, with non-vanishing elements only between the signatures connected by a bit jump where signatures are sequences of binary bits, 0 and 1. This should immediately remind us of Heisenberg chains of quantum spin-$`1/2`$ particles. Let us write the Hamilton function for one-dimensional ferromagnetic Heisenberg spin chain of length $`N`$, with Pauli variables $`\stackrel{}{\sigma }_j=(\sigma _j^\mathrm{x},\sigma _j^\mathrm{y},\sigma _j^\mathrm{z})`$ and coupling constants $`J(d)`$

$$\widehat{\mathrm{H}}=\frac{1}{2}\underset{j,d=1}{\overset{N}{}}J(d)\stackrel{}{\sigma }_j\stackrel{}{\sigma }_{j+d}=\frac{1}{2}\underset{j,d=1}{\overset{N}{}}J(d)\left\{\sigma _j^\mathrm{z}\sigma _{j+d}^\mathrm{z}+\frac{1}{2}\left(\sigma _j^+\sigma _{j+d}^{}+\sigma _j^{}\sigma _{j+d}^+\right)\right\},$$
(25)
where periodic boundary condition is assumed $`\sigma _{N+j}\sigma _j`$ and $`J(d)J(Nd),J(0)0`$. Operators $`\sigma _j^\pm =\sigma _j^\mathrm{x}\pm \mathrm{i}\sigma _j^\mathrm{y}`$ are the standard raising and lowering operators. Quantum state of the spin chain, an eigenstate of $`\sigma _j^\mathrm{z}`$, will be denoted by the signature $`|S`$, with the obvious interpretation. Bit $`a_j=1`$ or $`0`$ denotes spin $`j`$ up or down, respectively. Now we can calculate the matrix elements of the Hamilton operator (25). For the off diagonal elements we get

$$S|\widehat{\mathrm{H}}|S^{}=\{\begin{array}{cc}2J(d)& \text{if there is a jump of length }d\text{ connecting }S\text{ and }S^{}\hfill \\ 0& \text{otherwise}\hfill \end{array}.$$
(26)
Diagonal elements are

$$S|\widehat{\mathrm{H}}|S=J(N/2)\left[N/22s_{N/2}(S)\right]\underset{d=1}{\overset{N/21}{}}J(d)\left[N2s_d(S)\right],$$
(27)
where $`s_d(S)`$ is a number of different signatures that can be reached from the signature $`S`$ with a jump of length $`d`$. For any signature $`S`$ we can write an identity

$$\underset{d=1}{\overset{N/2}{}}s_d(S)=\left(\frac{N}{2}\right)^2.$$
(28)
By denoting $`s=M_0N^2/(2N2)`$, the following equality is also valid

$$\underset{i=1}{\overset{M_0}{}}s_d(S_i)=\{\begin{array}{cc}s& d=1,\mathrm{}N/21\hfill \\ s/2& d=N/2\hfill \end{array}.$$
(29)
The ground state for the ferromagnetic Heisenberg Hamiltonian (25) can be found immediately for any $`J(d)>0`$. We have to keep in mind that the Hilbert space is in our case spanned just by $`M_0`$ signatures of order 0 and not by all possible signature states as is usual for the quantum spin chains. Order $`i`$ of the signature is simply an eigenvalue of the $`\mathrm{S}_\mathrm{z}`$ component of the total spin, namely $`\mathrm{S}_\mathrm{z}=_j\sigma _j^\mathrm{z}/2=i`$. But the operator $`\mathrm{S}_\mathrm{z}`$ commutes with the Hamilton function and the Hilbert space is therefore a direct sum of subspaces labelled by eigenvalues of $`\mathrm{S}_\mathrm{z}`$. We choose $`M_0`$ dimensional Hilbert subspace with $`\mathrm{S}_\mathrm{z}(=i)=0`$. Ground state denoted by $`|0`$ is

$$|0=\frac{1}{\sqrt{M_0}}\underset{i=1}{\overset{M_0}{}}|S_i,$$
(30)
with the energy

$$E_00|\widehat{\mathrm{H}}|0=N\left(\underset{d=1}{\overset{N/21}{}}J(d)+\frac{1}{2}J(N/2)\right).$$
(31)
Now that we have the matrix elements of $`\mathrm{H}`$ and the ground state, we can see that if we prescribe $`J(d)`$ to equal $`c_d`$ multiplied by the constant prefactor in front of the matrix $`𝐂`$ in equation for $`𝐅`$ (24), we can formally write our flux matrix in terms of the Heisenberg spin Hamiltonian (25)

$$𝐅=𝐇+E_0𝐈.$$
(32)
Transition probabilities $`\mathrm{P}(\mathrm{S},\mathrm{S}^{};\mathrm{t})`$ can now be written as

$$𝐏(t)=\mathrm{exp}(𝐅t)=A(t)\mathrm{exp}(\beta 𝐇),$$
(33)
where we wrote $`\beta =t`$ and $`A(t)=\mathrm{exp}(E_0t)`$. Expression for $`𝐏`$ has the same form as the density operator for the quantum canonical distribution at temperature $`1/t`$. Time dependence of $`𝐏`$ is therefore the same as the dependence of the canonical distribution on cooling. Understanding time dynamics of the IFPU model is equivalent to the cooling of the ferromagnetic Heisenberg spin chain or its imaginary time dynamics. In the limit $`t\mathrm{}`$ the matrix elements of $`𝐏`$ go towards the constant value $`1/M_0`$ and in the corresponding Heisenberg spin chain any non-equilibrium distribution relaxes to the ground state $`|0`$.

For higher lying states we numerically solved the eigenvalue problem for the flux matrix at various $`\tau `$ and sizes up to $`N=12`$. One such example of a numerical spectrum is shown in figure 9.

We have said that with increasing $`N`$, at constant $`t_1/t_0`$, we expect coefficients $`c_d`$ to approach a constant value $`c_d=1`$ independent of $`d`$. This case corresponds to the simplest Heisenberg spin chain with uniform coupling. Apart from the constant factor the Hamilton function (25) is in this case $`\mathrm{H}=2𝐒^2+3N/2`$ and the matrix $`𝐅`$ reads

$$𝐅=\frac{2}{N^2\tau }\left(2𝐒^2\frac{N(N+2)}{2}𝐈\right),$$
(34)
where $`𝐒^2=1/4(\stackrel{}{\sigma }_1+\mathrm{}+\stackrel{}{\sigma }_N)^2`$. Eigenvalues of $`𝐒^2`$ are $`S(S+1)`$, where $`S`$ is the quantum number of the total spin with values from $`S=0,\mathrm{},N/2`$. The Hilbert space is composed of order 0 states with $`\mathrm{S}_\mathrm{z}=0`$. If we denote by $`\lambda _j`$ the eigenvalues of operator (34) and by $`n(j)`$ the corresponding multiplices we have

$`\lambda _j`$ $`=`$ $`{\displaystyle \frac{4j(N+1j)}{N^2\tau }}`$
$`n(j)`$ $`=`$ $`{\displaystyle \frac{(N+12j)N!}{(N+1j)!(j)!}}j=0,\mathrm{},N/2.`$ (35)
Particularly interesting are the first two eigenvalues $`\lambda _0`$ and $`\lambda _1`$. The biggest, non degenerate eigenvalue $`\lambda _0=0`$ belongs to the ground state $`|0`$. The largest nontrivial eigenvalue is $`\lambda _1=4/N\tau `$ with the multiplicity $`n(1)=N1`$. In the thermodynamic limit and keeping $`\tau `$ constant, it goes to zero as $`1/N`$. The correlation functions therefore decay slower than exponential in the thermodynamic limit. From the figures 9 and 10 it can also be deduced that with increasing $`\tau `$ the spectrum is indeed approaching the spectrum for the $`c_d1`$ (35), as predicted. What is more, the eigenvalue $`\lambda _1`$ seems to be monotonically approaching the limiting case ($`c_d1`$) from above. We can therefore conclude by observation that the correlation functions of observables which can be spanned by $`|S`$ decay slower than $`\mathrm{exp}(4t/N\tau )`$.

Now that we expect the spectrum to be close to the one for $`c_d1`$ we can explain the decay of the correlation functions for one particle state $`|s1`$ and the one signature state $`|S`$ in figure 8. State vector for $`|S`$ has only one component different from 0. Eigenvector $`|v1`$ corresponding to the eigenvalue $`\lambda _1`$ has on the other hand all components approximately of the same order, i.e. $`1/\sqrt{M_0}`$. Square of the scalar product is therefore $`|S|v_1|^21/M_0`$. One particle state $`|s1`$ is proportional to the $`|s1\sigma _1^+\sigma _1^{}|0`$. Because we act on the state $`|0`$ with the eigenvalue $`S=N/2`$ with the operator $`\sigma _1^+\sigma _1^{}`$, we can decompose the state $`|s1`$ as a linear combination of states with the eigenvalue $`S=N/2`$ and $`S=N/21`$. State $`|s1`$ can therefore be written as a sum of $`N`$ states, one of which is also $`|v1`$. If we assume all the expansion coefficients to be of the same order we immediately obtain $`|s1|v_1|^21/N`$.

## 6 Discussion and conclusion

We divide the phase space of a $`N`$-body Hamiltonian, namely the inverted FPU model, into $`2^N`$ cells that are uniquely tagged by a binary number, called a signature (3). Signatures are ordered according to the absolute difference in the number of particle pairs/bits in the left/0 and right/1 potential well. For a sufficiently high potential barrier an approximate Markovian description on space spanned by order 0 signatures is possible. The accuracy of the Markovian property has been checked numerically to be increasing with increasing chain length $`N`$ at a constant fraction of higher order states. In the flux matrix $`𝐅`$ we have non-vanishing matrix elements only between the signatures connected by an exchange of a pair of bits. This enables for the formal correspondence between the IFPU and ferromagnetic Heisenberg quantum spin-$`1/2`$ chains. Understanding time dependence of a Markov transition matrix is equivalent to cooling or imaginary time dynamics of a quantum spin chain. If we increase the chain length $`N`$ at a fixed fraction of higher order states, the transition rates are expected to become independent of a jump length $`d`$. Such trend is confirmed by a numerically calculated matrices $`𝐅`$. In this limit exact eigenvalues of a flux matrix can be obtained explicitly. The largest nontrivial eigenvalue in this case is $`\lambda _1=4/N\tau `$ and goes to $`0`$ in the thermodynamic limit. There is also a numerical evidence that the largest nontrivial eigenvalue in a general case (for a non-constant transition rate) is strictly larger than $`4/N\tau `$. Therefore, the correlation functions of the observables which can be spanned by functions $`|S`$ decay asymptotically slower than $`\mathrm{exp}(4t/N\tau )`$.

Unfortunately these results do not imply directly the behaviour of more general correlation functions, such as current-current correlation which is needed in order to understand the transport properties. We have made some numerical calculations of current-current time correlation functions of IFPU chain for different sizes $`N`$ and different parameters. As expected, the decay of correlations for observables that vary on time scales smaller than $`\tau `$ (e.g. current) is faster than the decay of correlation functions spanned by the piece-wise constant basis $`|S`$. Transition from the algebraic (anomalous conductivity) to the exponential (normal conductivity) decay of current autocorrelation function occurs rather abruptly at $`\alpha 3`$ (for $`N20`$).

Finally, though we suggested several interesting properties of the above model, we must admit that most of our conclusions are based on numerical evidence. Therefore, we believe that it should be a challenging and not impossible future task to try to provide more rigorous justifications (and perhaps proofs) of our results. Establishing rigorous asymptotic Markovian property would enable one to systematically code and enumerate all the many-body (unstable, hyperbolic) periodic orbits and to use them explicitly in a classical or semi-classical trace formulae, for example to calculate the transport coefficients directly. We feel that IFPU chain may become a useful toy model of a chaotic field theory .

## Acknowledgement

Financial support by the Ministry of Science and Technology of Slovenia is gratefully acknowledged.

## References


format: text

text:
# Generalization of a 3-D resonator model for the simulation of spherical enclosures11footnote 1Accepted for publication in Applied Signal Processing, 2000

## 1 Introduction

The feedback delay network (FDN) of order $`N`$, as depicted in fig. 1 for $`N=4`$, is the multivariable generalization of the recursive comb filter, and it has been widely used to simulate the late reverbation of an enclosure . The FDN with a diagonal matrix can be used to simulate a box with perfectly-reflecting walls. Energy absorption at the walls and in air can be accounted for by cascading the delay lines with lowpass filters. In regular rooms, some energy gets transferred from one mode to another due to non-mirror reflection at the walls. This effect is called diffusion  and is encompassed by the non-diagonal elements of the feedback matrix. As observed in the time domain, diffusion produces a gradual increase in density of the impulse response. The delay lengths of a FDN are sometimes chosen with number-theoretic considerations in order to minimize the overlapping of echoes, as it was done with classic reverberation structures . A more physical criterion is based on the following observation. In a rectangular enclosure, the distribution of normal modes can be obtained as the composition of (infinite) harmonic series, each series being associated with the spatial direction of propagation of the plane wave fronts supporting the modes. For instance, the longitudinal size of a rectangular box is associated with a low-pitch mode and with all its multiples. Since any harmonic series of resonances can be reproduced by means of a recursive comb filter, a reference FDN can be constructed as a parallel connection of comb filters or, in other words, with a diagonal feedback matrix. For the rectangular enclosure, the delay lengths can be computed exactly from the geometry of the room . Specifically, given a limited number of delay units $`N`$, it is possible to determine the $`N`$ delay lengths such that the model reproduces all the resonances of the room up to a limit frequency that depends on $`N`$ and the dimensions of the room. Conversely, the number $`N`$ of delay units can be determined in such a way that the limit frequency coincides with the Schroeder frequency , above which the exact reproduction of resonances is not perceptually relevant and it is sufficient to keep the modal density constant.

The fact that different modes are differently excited according to the position of the sound source can be taken into account by a proper choice of the $`𝐛`$ coefficients in fig. 1. Similarly, the pickup point should affect the choice of the $`𝐜`$ coefficients.

A model where normal modes, absorptions, and diffusion are all represented in the FDN structure has been called “Ball within the Box” (BaBo), since it can be visualized as a box of mirrors containing a single diffusing object (the Ball) represented by the matrix . Usually the model is completed with the explicit computation of early reflections that uses a geometric description of the acoustical scene.

Non-rectangular enclosures usually do not have an even distribution of resonances. In some relevant cases, however, the modal distribution can be calculated in closed form from the geometric specification of the enclosure. In particular, this paper deals with the spherical resonator, whose resonances can be found by computing the local extremal points of Bessel functions. The spherical Bessel functions tend to cosine functions for large values of the argument . A prior realization of the spherical resonator, depicted in fig. 2, exploited the fact that the extremal points are asymptotically equidistant, using recursive comb filters with feedback high-pass filters to reproduce the medium- and high-frequency resonances . On the other hand, a set of low-frequency resonances were individually reproduced by tuned second-order resonant filters. Such prior realization was successfully experimented in the AML $`|`$ Architecture and Music Laboratory, a museum installation where the visitor can experience how shapes such as, e.g. a tube, a cube or a sphere imprint a specific signature on the sounds. Informal reports from many listeners convinced us that it is indeed possible to identify basic shapes from the kind of resonance distribution they display.

Whereas the models used within the AML are specific to each shape, we try here, starting from the BaBo model, to design a single model valid for a few shapes. The BaBo model was initially designed for rectangular shapes but we extend it to the simulation of non-rectangular ones, in the hope that we can even feature a “shape control handle”. This paper reports the extension of the BaBo model to spherical enclosures and the comparison of audible results with recordings made through acoustical resonators. We also indicate how the simulation of cylinders of various lengths and radii might be achieved by the same structure with the proper parameter design.

## 2 Rectangular resonator model

The BaBo model provides parametric control over the geometric and physical properties of a rectangular enclosure . The kernel of the model is a feedback delay network where the delay lines have length in seconds given by

$$d_{l,m,n}=\frac{2}{c\sqrt{(l/X)^2+(m/Y)^2+(n/Z)^2}}$$
(1)
where $`c`$ is the speed of sound and $`l,m,n`$ are triplets of small positive integers sharing no common (nontrivial) divisor. Triplets with two zeros correspond to axial modes. Triplets with one zero correspond to tangential modes. Triplets with no zero correspond to oblique modes .

If the feedback matrix is diagonal we have a parallel of comb filters and this corresponds to a perfectly reflecting enclosure. In this case each comb filter represents the triplets $`(l,m,n)`$, $`(2l,2m,2n)`$, $`(3l,3m,3n)`$, etc., thus giving a perfectly harmonic series of resonances. In fact, each harmonic series of resonances is provided by a closed plane wave path propagating back and forth along a precise direction in space. If there is diffusion at the walls the plane wave fronts get scattered along many directions at every reflection. This diffusion phenomenon is reproduced in the BaBo model by non-diagonal matrix coefficients. If the elements of the matrix have all the same magnitude we have perfect, so called Lambertian, diffusion. Since the comb filters feed each other as the wavefronts scatter along several directions, the harmonicity of each comb series is somewhat broken in a diffusive enclosure.

Whereas the BaBo structure is well suited for imitating resonators characterized by overlapping series of harmonic resonances, it seems difficult to adapt it for the simulation of non-rectangular enclosures. Indeed, this is possible with moderate extra effort, as we will explain in section 5.

## 3 Acoustics of the sphere

The modal frequencies $`f_{ns}`$ of a sphere are proportional to the roots $`z_{ns}`$ of

$$j_n^{}(x)=0,n=0,1,\mathrm{}$$
(2)
where $`j_n`$ is the spherical Bessel function222The spherical Bessel function of order 0 is just the popular sinc function $`\mathrm{sin}x/x`$. of order $`n`$ and $`z_{ns}`$ is the $`s^{th}`$ root of the $`j_n^{}`$ function. The theoretical resonance frequencies are

$$f_{ns}=\frac{c}{2\pi a}z_{ns}$$
(3)
where $`c`$ is the propagation speed of sound and $`a`$ is the radius of the sphere .

The spherical Bessel functions $`j_n`$ can be studied by reference to the cylindrical Bessel functions $`J_n`$, thanks to the relation

$$j_n(x)=\sqrt{\frac{\pi }{2x}}J_{n+\frac{1}{2}}(x).$$
(4)
The function $`j_n^{}(x)`$ has the following property:

$$j_n^{}(x)=\frac{1}{2n+1}[nj_{n1}(x)(n+1)j_{n+1}(x)].$$
(5)
By substituting (4) in (5) we get

$$j_n^{}(x)=\frac{1}{2n+1}\sqrt{\frac{\pi }{2x}}[nJ_{n\frac{1}{2}}(x)(n+1)J_{n+\frac{3}{2}}(x)].$$
(6)
Some roots $`z_{ns}`$ of equation (2), i.e., the zeros of (6), are listed in Table 1. A closer look at the set of roots shows that they are not uniformly distributed, unlike the longitudinal resonances in tubes or between parallel boundaries. The roots are wider apart at low frequencies than at high frequencies. This effect is stronger for higher values of $`n`$ but, as it is clearly visible from Table 2, any series of roots tend to be periodic in $`\pi `$ for high values of $`s`$. This can be interpreted as dispersion at low frequencies and will give us a hint on how to implement the spherical resonator.

Table 3 shows the actual frequency values of the resonances of a sphere having radius $`a=0.188\text{m}`$, filled with air at the temperature333For the sake of accuracy, we take into account the fact that the propagation speed $`c`$ is temperature dependent according to the following relation:
$$c=331.8\sqrt{\frac{t+273}{273}}(m/s)$$
(7) where $`t`$ is the temperature $`{}_{}{}^{}\text{C}`$ of the air contained in the shape. $`t=23^{}\text{C}`$. This particular enclosure will be used throughout the paper as a test case.

It should be noticed that the Bessel functions of any order give rise to a resonance at dc, except for order $`n=1`$. In the design procedure of section 6 we neglect this singularity and we treat all the Bessel orders as if they would provide a resonance at dc. This approximation produces inaudible effects, as humans are insensitive to very low frequencies.

Relation (3) yields the theoretical resonant frequencies of the sphere. Since the envelope of the sphere might be vibrating as well as dissipating some acoustic energy, these frequencies should be corrected for the effects that occur at the boundary .

## 4 Measurements

Data are available from 3 experiments: Moldover et al.  display a spectrum measured in an argon-filled thick metal shell, and we have measured resonances in a rigid plastic shell as well as in an inflatable plastic ball.

In the experiment conducted on the metal shell, many roots $`z_{ns}`$ can be identified . The series of roots $`z_{1s}`$ can be measured but the amplitudes are very small. The series $`z_{0s}`$ shows larger amplitudes but the largest amplitudes are found at $`z_{22}`$, $`z_{42}`$ and $`z_{62}`$. Due to this observation, these three important resonances were explicitely modelled by second-order filters in the early implementation of fig. 2 . In general, in the experiment done with the metal shell, the series $`z_{2p,s}`$ are well represented, whereas the series $`z_{2p+1,s}`$ have lower amplitudes. The relative difference in amplitude between the resonances is due to different sensitivity to wall absorption and to the position of transducers during the measurement.

The resonant frequencies of the spherical loudspeaker depicted in fig. 3 were measured. In this loudspeaker 12 transducers are mounted on a spherical ABS plastic enclosure 444We thank Speaker Array Logic for providing the loudspeaker.. Several resonances could be identified with sufficient accuracy (see figures 45 and Tab. 4).

The resonances of an inflatable plastic ball having a diameter of $`0.67m`$ were also measured by posing the plastic ball onto a small loudspeaker. The loudspeaker was playing test signals through the ball and a microphone recorded the sound filtered by the ball (fig. 14) at the temperature of $`23^{}\text{C}`$ degrees. The position of the microphone was chosen in order to balance the amplitude of the various resonances. A spurious resonance showed up at $`200\text{Hz}`$, probably due to the imperfect coupling between the loudspeaker and the ball, so that it was necessary to prefilter the excitation signal with a parametric notch equalizer tuned at $`200\text{Hz}`$.

Fig. 15 shows the frequency response measured with a white noise generator and a spectrum analyzer. We found that the low-frequency resonances are systematically sharper than the theoretical values (see table 5). We assume that these deviations are due to the compliance offered by the plastic boundary, which can not be considered as a rigid wall at low frequency. The prominent resonances could be identified with $`f_{11}`$, $`f_{22}`$, $`f_{32}`$, $`f_{42}`$, $`f_{52}`$, $`f_{62}`$, $`f_{72}`$ and $`f_{92}`$.

## 5 Spherical resonator model

The first attempt to imitate spherical resonators, as shown in figure 2, has provided audible results that allow users to distinguish simulated spheres from other shapes such as cube, quader or tube just from their specific sound colour. This seems to confirm experiments on object recognition performed with blind subjects . From the AML model, we have learnt that the fundamental frequencies ($`f_{11},f_{22},f_{32},\mathrm{},f_{n2}`$) should be modelled with a high accuracy. The periodic spacing of resonances at higher frequencies is an effective parameter but it is difficult to tune. If all the resonances were spaced according to the asymptotic spacing $`\pi `$ of the $`z_{ns}`$ roots, they would imitate a cylinder rather than a sphere. On the other hand, the spacing at low frequencies is much larger than the asymptotic spacing (Tab. 2), so the spacing was tuned empirically at an intermediate value. In the early experiment, the $`z_{ns}`$ were grouped according to the number $`s`$ of the roots rather than to the order $`n`$ of the function because we noticed that the spacing in the $`z_{ns}`$ series was more regular along $`s`$ than along $`n`$.

Building on the experience gathered with this first model, we try here to design a more systematic one. Consider the perfectly-reflecting rectangular box and its representation in the BaBo model. Let us see how the model can be extended to spherical enclosures.

Indeed we would like to consider a perfectly reflecting sphere as a parallel connection of non-harmonic comb filters. The resonances of the $`n^{th}`$ comb filter correspond to the local extremal points of the $`n^{th}`$ order spherical Bessel function.

Given the resonance frequencies, we can sketch the ideal phase response of the $`n^{th}`$ comb filter loop, since the loop phase has to be equal to a multiple of $`2\pi `$ in order to sustain the mode associated with a resonance. Fig. 6 schematically shows with bold crosses the desired phase response at the positions of resonances, as it is generally found for any Bessel function order. A straight line fits the crosses fairly well for all Bessel orders, so that we can construct a realizable target phase response that is the sum of a linear contribution (realizable by a pure delay) and a piecewise linear contribution (that can be approximated by a stable allpass filter). The scheme of one of the inharmonic comb filters is shown in fig. 7, where the allpass filter and delay line are the main components of the feedback loop. In a practical implementation, one should also include a lowpass filter that accounts for frequency-dependent losses in a modal series. Such a filter can be designed starting from measured modal decay times. Alternatively, its parameters can be left open to user adjustment. In any case, we are neglecting losses at this stage, as we are mainly interested in the frequency distribution of resonances rather than in their relative strength.

Indeed, in the allpass filter design stage, we do not need to approximate all the piecewise linear dashed line of fig. 6. Instead, we just have to make sure that the designed phase response passes as close as possible to the points represented as small circles. Fig. 8.a shows with crosses the phase response of the $`0`$ order feedback loop at the resonance points. A monotonic phase curve interpolating those points can be obtained as the sum of a linear ramp and a nonlinear residual, also shown in fig. 8.a with dots and circles, respectively. According to the scheme of fig. 7, the linear component is given by a delay line whose length is equal to the slope of the linear ramp, and the nonlinear residual can be provided by an allpass filter.

A second design trade off appears here: up to which number $`S`$ should the $`f_{ns,0sS}`$ be approximated? Measurements and computations show that it is rather difficult to identify $`f_{ns,s>1}`$ resonances since slight inaccuracies suffice to change assignements. The experience gained in the first model (fig. 2) indicates that $`S`$ can not be too small if the “sense of roundness” has to be maintained in the medium and high frequency range. The order chosen for the allpass filter will determine the number $`S`$ of resonances that will be accurately imitated.

The nonlinear phase curve can be roughly approximated by a couple of linear segments, a low-frequency slope and a high-frequency slope. With this observation, the allpass filter for the $`0`$ order Bessel function can be designed by placing the poles on the unit circle according to these two slopes. Fig. 8.b shows a zero-pole distribution that gives the two-slope phase response interpolating the small circles, which is depicted in solid line in fig. 8.a. Figures 9 to 12 show the phase responses and pole-zero plots for the same sphere as in fig. 8, but with a Bessel function order going from 1 to 4. Table 6 shows the designed parameters for radius equal to $`0.188\text{m}`$ and $`0.32\text{m}`$, and Bessel functions going from 0 to 2. Reasonably shaped allpass phase responses can be obtained just by controlling two parameters, the argument of the first (low-frequency) pole, and the angular distance between any couple of contiguous poles.

## 6 Design Procedure and Examples

We have written a design procedure that finds the combinations of allpass filter and delay length for a set of sphere radii by iterative optimization over the position of the knee and the second, largest slope of the phase curve. In other words, the procedure iteratively changes the angle of the low-frequency pole and the relative distance between all the other equidistant poles. The relative contribution of the delay line to the loop phase response is also subject to optimization. The distance of the poles from the centre is kept fixed as it is mainly responsible for the magnitude of the ripples in the phase response. This pole radius plays a minor role because we are not interested in the phase values between the positions of resonance peaks and, therefore, we are not interested in a particular value of the phase ripples. A criterion for choosing the initial inter-pole distance is that the asymptotic distance between the peaks associated with a Bessel function is equal to a period of the ripple. The error that the iterative design procedure tries to minimize is the sum of the square phase differences at the target resonance points, weighted by some monotonically decreasing function that allows a more accurate positioning of resonances in low frequency.

If the goal is having a spherical resonator model that can be controlled in its radius, the design procedure should be run several times, each for a selected value of radius. In this way, a table can be constructed with a complete set of parameters for the resulting inharmonic comb filters. Notice that, even though a sixth-order allpass filter has six coefficients, the observation made in section 5 allows us to implement it as three second-order sections that can be controlled by two parameters: the angle of the first pole, and the angular distance of the following poles. The results of the design procedure for two distinct radii are reported in table 6. Since the precise position of resonances is of some importance only in the first few thousands Hz (say, $`4\text{kHz}`$, we see that a low order allpass filter is adequate for small spheres, e.g. order $`6`$ or less for radii smaller than $`0.5\text{m}`$. However, for larger spheres, the filter order should be increased in order to provide a decent approximation at least under the first kHz.

Fig. 13.a shows the frequency response of the parallel connection of dispersive comb filters, here designed for radius $`0.188\text{m}`$. The crosses represent the ideal modal positions for Bessel functions of order 0 to 4. In order to have a good match between the first resonance of each comb and its theoretical position, we properly shaped a weighting function to be used in the iterative optimization procedure. Maximum weight is used around the first resonance, while the following resonances become gradually less important. Psychoacoustic investigations should be conducted in order to better understand if the approximations introduced can be perceived and if they affect the perceived object shape. However, informal listening seems to indicate that significant deviations from the theoretical partial positions can be tolerated without loosing the “sense of roundness”. For instance, the sharpness of resonances measured in the plastic ball, and reported in Table 5, does not seem to prevent the listener from identifying the enclosure as a sphere.

Fig. 13.b shows how the comb filters from order 0 to 3 separately contribute to the response of fig. 13.a, which is obtained by pure summation of the comb outputs. Around some resonances different modes coming from different Bessel series interact with each other, and the local result is either a magnification or an attenuation of the peak. As well as with actual enclosures where the shape of the frequency response is dependent on the positions of exciter and pickup, with the FDN we can vary the shape of the response, without moving the resonances, just by changing the input and output coefficients , indicated as $`b_i`$ and $`c_i`$ in fig. 1.

## 7 Trials, Discussion, and Tuning

The ability of a model to simulate actual objects can be checked by comparing sounds played through the real objects with sounds processed by the model. Let us consider two different objects: a cylinder and a ball. Their shapes are far enough from each other to allow for a clear difference between processed sounds.

Figure 14 shows the setup used to record sounds processed by actual objects. The plastic ball is layed onto a small loudspeaker fed by a CD-player. Since the coupling between loudspeaker and sphere provoques a strong resonance at $`200Hz`$, a notch filter is inserted between the sound source and the loudspeaker. The microphone picks up the sound radiated in the recording room. In a similar fashion, sounds played in a tube closed at one end are recorded. Since the transducers and the room are not completely neutral, a third recording is made that serves as a reference: the sound played through the loudspeaker alone.

In order to compare the sounds played through the sphere to those played through the tube, some kind of tuning is necessary. Tuning a tube is easy but tuning a sphere is definitely not, because its resonances are not harmonic from each other. Several tuning methods have been considered. One could tune the lowest resonance of the sphere to that of the tube but that would lead to sounds located in different registers. We have chosen a tuning method that adjusts the brightness of the sounds. Another possible method would be to ajust the volume of the objects. This last method is not suitable here because a long tube with a small diameter will sound dramatically different from a short tube with a large diameter but having the same volume. This tuning is however meaningful for the comparison of a cube and a sphere. The brightness of the sounds played through the tube has been adjusted by modifying the length of the tube until it empirically fits the brightness of the sounds played through the sphere.

Some listeners easily recognise the shape of the object just by listening to their sounds . Most listeners agree that each object has a very peculiar sound colour but they fail at describing the differences and at telling one shape from the other. Many of them are capable of making a decision after a few learning trials, as if the knowledge had been available in their brain for a long time (perhaps since childhood) but that this knowledge had to be reactivated.

Similar tests were performed with sounds processed through the computer models. The FDN model of the sphere produces sounds with a refined reverberation, but the typical character of the sphere is less prominent. When comparing with FDN simulations of cubes, it is more difficult to tell the shape that is currently being simulated even though prior exposure still helps significantly. The fact that a sophisticated model performs worse than a simpler model for specific cases must be interpreted thinking that the initial model of the sphere (fig. 2) was constructed for a specific radius by fine tuning of a few prominent resonances. On the other hand, the FDN model aims at achieving generality and complex behaviour with a compact parametric structure. However, further research in sound perception has to be pursued in order to understand which are the salient characteristics of a sound spectrum that bring us its shape signature. When such results become available, the design procedure of the FDN spherical model will be improved further.

### 7.1 Including deviations from the ideal, rigid sphere

The FDN has been optimized according to the theory of the sphere but, in order to compare it with the plastic ball, deviations from the theoretical tuning had to be implemented. Such deviations (see table 5) can be introduced in our model just by moving the theoretical resonance positions in the procedure for designing the allpass filters. This can be done fairly easily if the deviations are small, otherwise it can be difficult to assign a certain resonance to an inharmonic comb series. Alternatively, one can start with the filters designed for the ideal sphere and adjust the position of the first pole, as we did to obtain the frequency response of fig. 16, which should be compared with fig. 15. In the feedback loop, we used second-order FIR filters (exhibiting a one-sample delay) to simulate the faster attenuation of higher modes. Moreover, a first-order lowpass filter has been cascaded with the whole structure in order to resemble the lowpass characteristic of fig. 15.

This tuning of the model improves the simulation but it is not yet enough to allow for a reliable detection of the simulated shape from the sound itself. Again, it seems that the physical theory of the sphere is not sufficient to drive the design stage, and perceptual issues should be considered. To this end, the issue of identifying a shape from its audible signature should be investigated more deeply in the future.

## 8 Conclusion and Further Research

We have shown that the basic structure of the BaBo model can also be used to simulate spherical geometries, just by using properly designed allpass filters within the filtering blocks $`H_i`$ of fig. 1. The resulting model might be called “Ball within the Ball” (BaBa). We have shown that a simple design procedure and, possibly, some manual tweaking, allow an efficient structure to be realized that can be tuned either like an ideal sphere or like a real one. An open question is nevertheless how accurate this design has to be, since it is still unsettled which are the relevant parameters that affect the “perceived roundness” of the object.

An important aspect of the model is that very few parameters are added to the BaBo model to control the allpass filters of the spherical model. Even deviations from ideal boundary conditions are reasonably achieved by moving only the position of one pole per inharmonic series. So far, we have simulated inharmonic series given by Bessel functions of order ranging from 0 to 6. In many cases, higher order inharmonic series are needed to achieve realism, but the fundamental resonance of those series seems to be most important while the higher resonances get drowned in the dense mixture of resonances from other series and their position is out of the bandwidth of perceived inharmonicity. So, we suggest to implement higher-order series as harmonic comb filters tuned to the fundamental frequency of that series. Some form of delay interpolation  might be needed to position such fundamental frequencies with sufficient accuracy.

The same structure used for spherical resonators might be as well used for cylindrical resonators. In this case there is a harmonic series given by longitudinal modes, superimposed with inharmonic series of resonances whose positions are determined by the extremal points of cylindrical Bessel functions. Therefore, the feedback delay network should have the first delay line accounting for the longitudinal modes, and the remaining lines, cascaded with properly-designed allpass filters, accounting for the inharmonic, transversal modal series.

If the inharmonic series, each corresponding to a Bessel function of a certain order, are recreated by comb filters having a delay line and an allpass filter in the feedback loop, it is conceivable to control the degree of “roundness” of the enclosure by changing the relative contribution to the overall phase response given by the delay and by the allpass filter. Namely, if all the delays are increased we gradually move from a sphere to a cube with rounded faces. This continuous shape control, as well as the extension of the BaBo model to cylindrical shapes will be covered in future research.