text:
# Semiclassical trace formulas for noninteracting identical particles

## I Introduction

In the semiclassical limit of quantum mechanics, the periodic orbits of the corresponding classical system play a special role in determining the spectral properties of the quantum system. This fundamental fact has been a dominant theme in modern semiclassical physics and was pioneered by Gutzwiller , Balian and Bloch , Strutinsky and Magner  and Berry and Tabor . One of the central results which emerged from this work is the representation of the density of states in terms of classical periodic orbits. Such representations are referred to as trace formulas. Semiclassical analysis based on the use of trace formulas is now common in many areas of physics . Besides providing a natural framework for studying the quantum manifestations of classical chaos , such analysis has been used in the study of nuclei , atoms , metal clusters , molecules , chemical systems , spins , Casimir effects  and tunneling . Trace formulas have also become a prominent analytical tool in the study of mesoscopic systems . New directions continue to be explored .

Despite the vast utility of trace formulas, their use in the few-body or many-body context has received little attention. Although trace formulas are applicable to interacting many-body systems, more effort has gone into developing semiclassical descriptions of single-particle dynamics in an appropriate mean field. One impressive exception is the application of the Gutzwiller trace formula to the study of two-electron atoms and related three-body systems . The main difficulty of applying the theory is that periodic orbits must be found for the interacting many-body system. One approach to this problem has been proposed in Ref.  which develops a particle number expansion of the trace formula.

In this paper, we consider the extension of the Gutzwiller trace formula to systems of noninteracting identical particles. The semiclassical analysis of such systems is quite subtle. As we shall discuss below, if the particles are noninteracting, there is a continuous symmetry. (Discrete symmetries in semiclassical trace formulas are discussed in Refs.  and continuous symmetries in Refs. .) This symmetry has a profound consequence; although an $`N`$-particle system in $`d`$ dimensions can be thought of as a single-particle system in $`Nd`$ dimensions, one cannot simply apply the Gutzwiller trace formula since the presence of a continuous symmetry implies the periodic orbits of the full phase space are not isolated, but rather occur in continuous families.

Recently, we presented a semiclassical formalism for the density of states of two noninteracting identical particles based on an asymptotic analysis of convolution integrals that arise in the decomposition of the semiclassical two-particle density of states . In principle, this approach can be generalised to more than two identical particles. In this paper, we pursue a different approach which uses the full phase space rather than the individual phase spaces of each particle. We show that the formalism for continuous symmetries  can be used to find many-body trace formulas. This approach recovers our previous results, but can also be more easily generalised to arbitrary particle numbers. In addition, it is conceptually cleaner than the convolution method since spurious end-point contributions from convolution integrals do not arise and therefore need not be explained away. The most important difference is that the convolution method cannot be used when there are interactions between the particles whereas the analysis of this paper can be extended to include interactions .

It is also important to understand the effect of particle symmetry on the semiclassical structure of many-body trace formulas. For noninteracting identical particles, there are coexisting discrete and continuous symmetries. While Ref.  considers the symmetry-reduced trace formula due to the discrete permutational symmetry, it is assumed the periodic orbits are isolated, which is only true if the particles are strongly interacting (although there is a brief discussion of the noninteracting case). We include the appropriate continuous symmetries to determine the trace formulas for the bosonic and fermionic densities of states.

This paper is organised as follows. In section II, we study the case of two noninteracting identical particles. We first provide the necessary background material in sections II.1-II.3 and then give the semiclassical formulation in the full phase space. Section III considers the extension to $`N`$ identical noninteracting particles. The symmetry decomposition of the $`N`$-particle density of states is examined in section IV. The results of a numerical study of the three-particle cardioid billiard are then presented to illustrate and test the results of the paper. We finish the paper with a conclusion and several appendices.

## II Two Noninteracting Identical Particles

### II.1 Quantum density of states

The quantum Hamiltonian for two identical noninteracting particles, $`a`$ and $`b`$, is

$$\widehat{H}=\widehat{h}(\widehat{z}_a)+\widehat{h}(\widehat{z}_b),$$
(1)
where $`\widehat{z}_{a/b}`$ denote the set of operators $`(\widehat{x}_{a/b},\widehat{p}_{a/b})`$ and $`\widehat{h}`$ is a one-particle Hamiltonian. The full Hamiltonian (1) is invariant under the unitary transformation $`\widehat{U}`$ which exchanges $`a`$ and $`b`$. We define the single-particle energies and eigenstates by

$$\widehat{h}|j=\epsilon _j|j.$$
(2)
Then, the two-particle energies and eigenstates are $`E_{ij}=\epsilon _i+\epsilon _j`$ and $`|ij`$ so that

$$\widehat{H}|ij=E_{ij}|ij.$$
(3)
Accordingly, the one and two-particle densities of states are

$`\rho _1(\epsilon )`$ $`=`$ $`{\displaystyle \underset{j}{}}\delta (\epsilon \epsilon _j),`$
$`\rho _2(E)`$ $`=`$ $`{\displaystyle \underset{i,j}{}}\delta (EE_{ij}),`$ (4)
and they are related by the convolution identity

$$\rho _2(E)=\left(\rho _1\rho _1\right)(E).$$
(5)
A useful result is the relation between the density of states and the trace of the energy Green function or resolvent. We define $`g(E)=\text{Tr}(\widehat{G}(E))`$, where $`\widehat{G}(E)=1/(E\widehat{H})`$ is the one-sided Fourier transform of the quantum propagator. In terms of the resolvent,

$$\rho (E)=\frac{1}{\pi }\text{Im}\left\{g(E+iϵ)\right\},$$
(6)
and this applies for either the one or two-particle density of states as long as we use the appropriate resolvent on the right-hand side. In the limit $`ϵ0^+`$, the exact density of states is recovered . Henceforth, the $`iϵ`$ will be implicit.

### II.2 Symmetry decomposition

The most interesting aspect of the existence of identical particles is the fact that only certain states are occupied, the fully symmetric ones if the particles are bosons or the fully antisymmetric ones if the particles are fermions. It is important to understand how the above discussion decomposes when we consider the separate densities of symmetric and antisymmetric states. Although not absolutely necessary for the present discussion, it will be useful for later to introduce projection operators. As mentioned, the Hamiltonian (1) is invariant under exchange of the particles $`a`$ and $`b`$, an operation we denote by $`\varsigma `$ (leaving the particles unchanged we denote by $`\iota `$). There is a two-element discrete group comprised of these operations and the representation of these group elements in the Hilbert space (i.e. the quantum operators that exchange the particles) are $`\widehat{U}`$ and $`\widehat{I}`$. Both of these operators commute with $`\widehat{H}`$. This is a simple group with two irreducible representations (irreps) which we identify as the bosonic (symmetric) representation and the fermionic (antisymmetric) representation. Given an arbitrary state with components belonging to both irreps, we can project out the portion belonging to each irrep through the use of the projection operators

$$\widehat{P}_\pm =\frac{1}{2}\left(\widehat{I}\pm \widehat{U}\right),$$
(7)
where the $`\pm `$ refer to the bosonic and fermionic irreps, respectively.

In terms of these projection operators, the bosonic and fermionic densities of states are given as

$$\rho _\pm (E)=\text{Tr}\left(\widehat{P}_\pm \delta (E\widehat{H})\right).$$
(8)
The sum of the bosonic and fermionic densities is the complete two-particle density of states, $`\rho _2(E)`$. The difference is given by $`\text{Tr}(\widehat{U}\delta (E\widehat{H}))`$ and expressing the trace in the energy eigenbasis,

$`\rho _+(E)\rho _{}(E)`$ $`=`$ $`{\displaystyle \underset{i,j}{}}ij|\widehat{U}\delta (E\widehat{H})|ij`$ (9)
$`=`$ $`{\displaystyle \underset{i,j}{}}ji|ij\delta (EE_{ij})`$
$`=`$ $`{\displaystyle \underset{j}{}}\delta (E2\epsilon _j),`$
where we have used the fact that $`\widehat{U}`$ exchanges the state labels in the second line and the fact that $`E_{jj}=2\epsilon _j`$ in the third. The final line we recognise as $`\rho _1(E/2)/2`$ and thereby conclude

$$\rho _\pm (E)=\frac{1}{2}\left(\rho _2(E)\pm \frac{1}{2}\rho _1\left(\frac{E}{2}\right)\right).$$
(10)
### II.3 Review of semiclassical formulation

It is common to decompose the semiclassical density of states into smooth and oscillatory components. For the one-particle density

$$\rho _1^{\text{sc}}(\epsilon )=\overline{\rho }_1(\epsilon )+\stackrel{~}{\rho }_1(\epsilon ),$$
(11)
where $`\overline{\rho }`$ and $`\stackrel{~}{\rho }`$ denote the smooth and oscillating components, respectively. There is an extensive literature on this decomposition . We adopt the point of view that one can simply use the leading-order contribution of each component. We do not consider the subtle issues related to the asymptotic nature of this decomposition (see, for example, Refs. ). For analytic potentials in $`n`$ dimensions, the leading-order term for the smooth density of states is

$$\overline{\rho }_1(\epsilon )\frac{1}{(2\pi \mathrm{})^n}\text{d}z\delta (\epsilon h(z)),$$
(12)
where $`z`$ collectively denotes the $`2n`$ classical phase space coordinates and $`h(z)`$ is the classical Hamiltonian (for an exception to this general result, see Refs. ). There are corrections to (12) involving derivatives of the delta function in the integrand. The first correction is $`O(\mathrm{}^2)`$. For a two-dimensional billiard, the analogous expression is

$$\overline{\rho }_1(\epsilon )\frac{\alpha 𝒜}{4\pi }\pm \frac{\sqrt{\alpha }}{8\pi \sqrt{\epsilon }}+𝒦\delta (\epsilon ),$$
(13)
where $`\alpha =2m/\mathrm{}^2`$, $`𝒜`$ and $``$ refer to the area and perimeter, respectively and the $`\pm `$ refers to Neumann and Dirichlet boundary conditions, respectively. The third term

$$𝒦=\frac{1}{12\pi }\text{d}l\kappa (l)+\frac{1}{24\pi }\underset{i}{}\frac{\pi ^2\theta _i^2}{\theta _i}$$
(14)
is the average curvature integrated along the boundary with corrections due to corners with angles $`\theta _i`$. It does not actually contribute to the density of states, but rather to its first integral and this term will be used in section V. There are also corrections involving powers and derivatives of the curvature (see Refs.  for more exhaustive studies). Similar results hold for higher dimensional billiards (see Ref. ).

The oscillating component can be written as

$$\stackrel{~}{\rho }_1(\epsilon )\frac{1}{\pi }\text{Im}\left\{\underset{\gamma }{}A_\gamma (\epsilon )\mathrm{exp}i\left(\frac{S_\gamma (\epsilon )}{\mathrm{}}\sigma _\gamma \frac{\pi }{2}\right)\right\},$$
(15)
where $`\gamma `$ labels the periodic orbits of the system, $`S_\gamma `$ is the classical action integral along the orbit and $`\sigma _\gamma `$ is a topological index  counting the caustics in phase space encountered by the orbit. $`A_\gamma `$ is the amplitude of the periodic orbit and depends on the specific nature of the system, for example, whether the orbit is isolated or not and if so on its stability. For the case of isolated periodic orbits,

$$A_\gamma (\epsilon )=\frac{1}{i\mathrm{}}\frac{T_\gamma ^0(\epsilon )}{\sqrt{\left|det(\stackrel{~}{M}_\gamma I)\right|}},$$
(16)
where $`T_\gamma ^0`$ is the primitive period of the orbit and $`\stackrel{~}{M}_\gamma `$ is the $`2(n1)\times 2(n1)`$ symplectic stability matrix on any Poincaré section to which the orbit is transverse. Its eigenvalues give the stability exponents of the orbit.

As stated above, the density of states for two noninteracting particles is the autoconvolution of the one-particle density of states (5). Formally, the semiclassical two-particle density of states is the autoconvolution of (11), that is,

$$\rho _2^{\text{sc}}(E)=\left(\rho _1^{\text{sc}}\rho _1^{\text{sc}}\right)(E)=\overline{\rho }_2(E)+\stackrel{~}{\rho }_2(E),$$
(17)
where

$`\overline{\rho }_2(E)=\left(\overline{\rho }_1\overline{\rho }_1\right)(E),`$ (18a)
$`\stackrel{~}{\rho }_2(E)=2\left(\overline{\rho }_1\stackrel{~}{\rho }_1\right)(E)+\left(\stackrel{~}{\rho }_1\stackrel{~}{\rho }_1\right)(E).`$ (18b)
The mixed term $`2\left(\overline{\rho }_1\stackrel{~}{\rho }_1\right)(E)`$ also belongs to the oscillating component of $`\rho _2^{\text{sc}}(E)`$. This is because an asymptotic endpoint analysis of the convolution integral results in an oscillatory function as shown in Ref.  where all components have been evaluated and given explicit semiclassical interpretations in terms of one and two-particle dynamics which support this decomposition. We also showed above how the difference between the bosonic and fermionic densities is given by the one-particle density of states. Formally, we may write (analogous to Eq. (10)),

$$\rho _\pm ^{\text{sc}}(E)=\frac{1}{2}\left(\rho _2^{\text{sc}}(E)\pm \frac{1}{2}\rho _1^{\text{sc}}\left(\frac{E}{2}\right)\right).$$
(19)
The formal results (17-19) can also be understood from a semiclassical analysis in the full phase space. This analysis is not only more fundamental, it is necessary if one wants to include interparticle interactions since the particle dynamics then become coupled and we can no longer make use of calculations that involve the individual one-particle phase spaces. In the following sections, we derive trace formulas for the full and symmetry-reduced densities of states from semiclassical calculations in the full two-particle phase space. Since we are mainly concerned with the extension of the Gutzwiller theory, we focus on the fluctuating part of the density of states. However, since the smooth part is important in constructing the complete density of states, we have provided a discussion of the two-particle Thomas-Fermi term (and the associated symmetry decomposition based on the theory of Ref. ) in Appendix A. To calculate the fluctuating part of the density of states, we need to find all periodic orbits in the full phase space at a specified energy $`E`$.

### II.4 Two-particle dynamics in the full phase space

The two identical particles, $`a`$ and $`b`$, evolve independently in their own one-particle configuration space which we denote as having dimension $`d`$ so that the one-particle phase spaces are of dimension $`2d`$. The full two-particle configuration space has dimension $`2d`$ and the corresponding phase space is of dimension $`4d`$. We reserve the symbol $`z`$ to collectively denote these $`4d`$ phase space coordinates and will use $`z=(z_a,z_b)`$ where $`z_{a/b}`$ denote the $`2d`$-dimensional one-particle phase space coordinates of each particle. We recall that dynamics in the full phase space consists of each particle evolving separately in its own phase space. The dynamics of $`z`$ are defined through one-particle dynamics by $`\mathrm{\Phi }_tz=(\varphi _tz_a,\varphi _tz_b)`$, where $`\varphi _t`$ is the flow for one particle. The (noninteracting) two-particle Hamiltonian is $`H(z)=h(z_a)+h(z_b)`$, where $`h(z_{a/b})`$ is a one-particle Hamiltonian.

We seek periodic orbits with phase space coordinates $`z^{}`$ such that $`\mathrm{\Phi }_Tz^{}=z^{}`$ for some period $`T`$. This is possible if the two particles are on (generally distinct) periodic orbits with the same period. In general, two arbitrary periodic orbits will have different periods. However, there is a parameter which we can vary, namely the way in which the total energy is partitioned between the two particles. Generally, we can find an energy $`E_a`$ (and $`E_b=EE_a`$) such that the two periods are the same. We will assume henceforth that there is only one energy $`E_a`$ for which there is a solution. (This assumption can be relaxed at the cost of heavier notation.) There is another way to have a periodic orbit in the full phase space; one particle can evolve dynamically on a periodic orbit with all of the energy while the other is stationary at a fixed point of the potential. This is discussed later.

#### II.4.1 Dynamical periodic orbits

If *both* particles are on periodic orbits, we call the full phase space periodic orbit a dynamical periodic orbit. We first note that such orbits occur in continuous families. To see this, imagine that a full phase space periodic orbit consists of one particle on a periodic orbit $`\gamma _a`$ and the other particle on a distinct periodic orbit $`\gamma _b`$ (see Fig. 1) and that the energy partition is such that both orbits have the same period $`T`$. We have complete freedom in specifying which points on the respective orbits we choose as initial conditions. Given that we define $`t=0`$ to be when particle $`b`$ is at some specified point on $`\gamma _b`$, we can vary the position of particle $`a`$ on $`\gamma _a`$. By changing its initial position along the orbit, we map out a continuous family of congruent periodic orbits.

This can be formalised as follows. We note that in addition to the total Hamiltonian $`H`$, there is a second constant of motion $`J=h(z_a)`$ in involution with $`H`$. It generates time translations of particle $`a`$ while leaving particle $`b`$ fixed. (In fact, $`J`$ can be chosen as any linear combination of $`h(z_a)`$ and $`h(z_b)`$ as long as it is independent of $`H`$.) Flows generated by $`J`$ are denoted by $`\mathrm{\Psi }_\theta `$ and are mapped in the full phase space as follows: $`\mathrm{\Psi }_\theta z=(\varphi _\theta z_a,z_b)`$. The symmetry parameter $`\theta `$ is conjugate to $`J`$ and has the dimension and interpretation of time. However, since it only measures the evolution of particle $`a`$, it is not time in the usual sense and we will follow the notation of Ref.  in denoting the paramater by $`\theta `$. A combination of flows in $`H`$ and $`J`$ is

$$\mathrm{\Phi }_t\mathrm{\Psi }_\theta z=(\varphi _t\varphi _\theta z_a,\varphi _tz_b)=(\varphi _{t+\theta }z_a,\varphi _tz_b).$$
(20)
Since $`\mathrm{\Psi }_\theta `$ and $`\mathrm{\Phi }_t`$ commute and separately conserve both $`H`$ and $`J`$, the surface mapped out by these flows has constant $`H`$ and $`J`$ (i.e. $`H=E`$ and $`J=E_a`$) . Starting at some point on the full phase space periodic orbit, flows in $`H`$ and $`J`$ map out a two-dimensional torus. This means there is a 1-parameter degenerate family of periodic orbits (the other dimension is parameterised by time and is present even in the case of isolated orbits). Therefore, we cannot use the Gutzwiller trace formula for isolated orbits since it will give a spurious infinity. Due to the continuous family, there is one fewer stationary phase integrals to be done in evaluating the trace so that this family of orbits contributes $`O(1/\sqrt{\mathrm{}})`$ more strongly than an isolated orbit and the calculation of its amplitude must be performed carefully.

For the present, we assume that there are no symmetries other than $`J`$ so that all periodic orbits of the one-particle phase space are isolated. The flow directions generated by $`H`$ and $`J`$ are stable as are the two directions transverse to the constant $`H`$ and $`J`$ surfaces. Thus, there are four directions of neutral stability in phase space. The remaining $`(4d4)`$ directions decompose into separate subspaces of dimension $`(2d2)`$ within each of which there are the standard symplectic possibilities for stability.

In general, the leading-order contribution of one $`f`$-parameter family of orbits (generated by Abelian symmetry) to the resolvent is

$$\begin{array}{cc}\hfill \stackrel{~}{g}_\mathrm{\Gamma }(E)& =\frac{1}{i\mathrm{}}\frac{1}{\left(2\pi \mathrm{}\right)^{f/2}}\frac{T_\mathrm{\Gamma }^0V_\mathrm{\Gamma }^0}{\left|\frac{𝚯}{𝐉}\right|_\mathrm{\Gamma }^{1/2}\left|det\left(\stackrel{~}{M}_\mathrm{\Gamma }I\right)\right|^{1/2}}\hfill \\ & \times \mathrm{exp}i\left(\frac{S_\mathrm{\Gamma }(E)}{\mathrm{}}(\mu \delta )_\mathrm{\Gamma }\frac{\pi }{2}f\frac{\pi }{4}\right).\hfill \end{array}$$
(21)
This contribution is $`O(1/\mathrm{}^{f/2})`$ stronger than an isolated periodic orbit. As mentioned above, every constant of motion implies one fewer stationary phase integrals and therefore $`f`$ fewer powers of $`\sqrt{\mathrm{}}`$ in the prefactor. For a similar reason, there is an additional phase factor of $`f\pi /4`$. The total contribution to the resolvent is a sum over all families of periodic orbits, the capital $`\mathrm{\Gamma }`$ indicating that these are indeed families and not isolated orbits as in the more familiar Gutzwiller trace formula. In our case, the sum over $`\mathrm{\Gamma }`$ can be expressed as a double sum over $`\gamma _a`$ and $`\gamma _b`$ indicating the periodic orbits on which the particles are evolving. We now describe the various factors in (21) and explain what they are in the present situation for which $`f=1`$.

The volume term, $`T_\mathrm{\Gamma }^0V_\mathrm{\Gamma }^0`$ is the integral over the flows generated by $`H`$ and $`J`$, $`_\mathrm{\Gamma }\text{d}t\text{d}\theta `$, integrated over the periodic orbit family. The time integral gives the period of the family $`T_\mathrm{\Gamma }=T_{\gamma _a}(E_a)=T_{\gamma _b}(E_b=EE_a)T`$ while the $`\theta `$ integral gives $`V_\mathrm{\Gamma }=T_{\gamma _a}(E_a)`$ since a flow in $`J`$ by that amount returns particle $`a`$ to where it began. (Hence, the intial phase space coordinate is mapped back to itself under the dynamics.) However, there can be discrete symmetries such that a combination of flows in $`H`$ and $`J`$ for less than $`T`$ restores the intial conditions. This situation occurs when one or both particles is on a repetition of some primitive periodic orbit. To see this, suppose that particle $`a`$ is on the $`n_a`$’th repetition of its orbit while particle $`b`$ is on the $`n_b`$’th repetition of its orbit. Then, the torus is partitioned into $`n_an_b`$ equivalent segments and the primitive volume term is $`T_\mathrm{\Gamma }^0V_\mathrm{\Gamma }^0=T_\mathrm{\Gamma }V_\mathrm{\Gamma }/n_an_b`$. However, the full periods are defined through the primitive periods by $`T_{\gamma _a}(E_a)=n_aT_{\gamma _a}^0(E_a)`$ and similarly for particle $`b`$. Thus, $`T_\mathrm{\Gamma }^0V_\mathrm{\Gamma }^0=T_{\gamma _a}^0(E_a)T_{\gamma _b}^0(E_b=EE_a)`$, which is the product of the primitive periods.

$`\stackrel{~}{M}_\mathrm{\Gamma }`$ is a $`(4d4)\times (4d4)`$ matrix linearising motion on a reduced surface of section. Specifically, it is the section at constant $`(H,J,x_a,x_b)`$ where $`x_{a/b}`$ are chosen so that the dynamics are transverse to the surface on which they are both constant. In our case, this section is simply the direct product of the normal Poincaré surfaces of section for each of the two motions (where one would specify the one-particle energy and some fixed coordinate in each case). As a result, $`\stackrel{~}{M}_\mathrm{\Gamma }`$ has a block diagonal structure since there is no coupling between the two particle spaces. We conclude that $`det(\stackrel{~}{M}_\mathrm{\Gamma }I)=det(\stackrel{~}{M}_{\gamma _a}I)det(\stackrel{~}{M}_{\gamma _b}I)`$, where $`\stackrel{~}{M}_{\gamma _a/\gamma _b}`$ are the stability matrices of the separate one-particle periodic orbits and $`I`$ is the appropriately dimensioned unit matrix on both sides of the equality.

The anholonomy term $`\left(\mathrm{\Theta }/J\right)_\mathrm{\Gamma }`$ measures the amount by which orbits that are periodic in the symmetry-reduced dynamics fail to be periodic in the full phase space. Suppose we vary the value of $`J`$ infinitesmally while keeping the total energy fixed; in our case this amounts to a slight change of the energy partition between the two particles. The periodic orbit is launched as before with the same initial conditions except for $`p_{}`$ (the momentum conjugate to $`x_{}`$) which must be changed appropriately to effect the change in $`J`$. After the original period $`T`$, an initial phase space coordinate will not be mapped back to where it began, but rather infinitesmally close to this initial condition. A flow in $`H`$ for some extra amount of time $`\mathrm{\Delta }t`$ and a flow in $`J`$ by an extra amount $`\mathrm{\Delta }\theta `$ (or vice-versa since the flows commute) closes the orbit in the full phase space. The factor $`\mathrm{\Theta }/J`$ is simply the ratio $`\mathrm{\Delta }\theta /\mathrm{\Delta }J`$ (in the limit that $`\mathrm{\Delta }J0`$). ($`\mathrm{\Theta }`$ is capitalised to stress that $`J`$ and $`\theta `$ can also be used as labels of families of surfaces, in which case, this factor can be interpreted as a Jacobian for a change of label from $`J`$ to $`\theta `$.) Recall that the value of $`J=h(z_a)`$ is just the energy of particle $`a`$. If $`J_aJ_a+\mathrm{\Delta }J_a`$, then $`E_bE_b\mathrm{\Delta }J_a`$, since the total energy is fixed. $`\gamma _a`$ now has a perturbed period, $`T+\mathrm{\Delta }T_{\gamma _a}=T+T_{\gamma _a}^{}\mathrm{\Delta }J_a`$ while $`\gamma _b`$ now has a perturbed period $`T+\mathrm{\Delta }T_{\gamma _b}=TT_{\gamma _b}^{}\mathrm{\Delta }J_a`$, where the primes denote differentiation with respect to energy:

$$T_{\gamma _a}^{}=\frac{\text{d}T_{\gamma _a}(\epsilon )}{\text{d}\epsilon }|_{\epsilon =E_a},T_{\gamma _b}^{}=\frac{\text{d}T_{\gamma _b}(E\epsilon )}{\text{d}\epsilon }|_{\epsilon =E_a}.$$
(22)
Let $`z^{}=(z_a^{},z_b^{})`$ and $`z`$ denote the initial and final phase space coordinates, respectively. Then, after the original period $`T`$,

$$\mathrm{\Phi }_Tz^{}=z=(\varphi _{\mathrm{\Delta }T_{\gamma _a}}z_a^{},\varphi _{\mathrm{\Delta }T_{\gamma _b}}z_b^{}).$$
(23)
We need to find $`(\mathrm{\Delta }t,\mathrm{\Delta }\theta )`$ which map $`z`$ back to $`z^{}`$. Using Eq. (20), the condition for a periodic orbit $`\mathrm{\Phi }_{\mathrm{\Delta }t}\mathrm{\Psi }_{\mathrm{\Delta }\theta }z=z^{}`$ implies $`\mathrm{\Delta }t=\mathrm{\Delta }T_{\gamma _b}`$ and $`\mathrm{\Delta }\theta =\mathrm{\Delta }T_{\gamma _a}\mathrm{\Delta }T_{\gamma _b}`$ so that

$$\frac{\mathrm{\Theta }}{J}=T_{\gamma _a}^{}+T_{\gamma _b}^{}.$$
(24)
The action $`S_\mathrm{\Gamma }(E)=S_{\gamma _a}(E_a)+S_{\gamma _b}(E_b=EE_a)`$ is the action of the periodic orbits in the family (all orbits in $`\mathrm{\Gamma }`$ have the same action because of symmetry). Finally, we discuss the phase indices. $`\mu _\mathrm{\Gamma }`$ is determined from the dynamics in the symmetry-reduced surface of section in the same way as for isolated orbits in the usual Gutzwiller trace formula and following the same logic as above, $`\mu _\mathrm{\Gamma }=\sigma _{\gamma _a}+\sigma _{\gamma _b}`$. $`\delta _\mathrm{\Gamma }`$ is defined as the number of positive eigenvalues of $`\left(\mathrm{\Theta }/J\right)_\mathrm{\Gamma }`$ . In this case, the anholonomy term is simply a scalar and therefore $`\delta _\mathrm{\Gamma }=1`$ if the Jacobian is positive and $`\delta _\mathrm{\Gamma }=0`$ if the Jacobian is negative. We conclude that the contribution to the resolvent from one family of dynamical orbits is

$$\stackrel{~}{g}_\mathrm{\Gamma }^\text{d}(E)=\frac{i\mathrm{}}{\sqrt{2\pi \mathrm{}}}\left(\frac{T_{\gamma _a}^0(E_a)}{i\mathrm{}}\frac{\mathrm{exp}i\left(\frac{S_{\gamma _a}(E_a)}{\mathrm{}}\sigma _{\gamma _a}\frac{\pi }{2}\right)}{\left|det\left(\stackrel{~}{M}_{\gamma _a}I\right)\right|^{1/2}}\right)\left(\frac{T_{\gamma _b}^0(E_b)}{i\mathrm{}}\frac{\mathrm{exp}i\left(\frac{S_{\gamma _b}(E_b)}{\mathrm{}}\sigma _{\gamma _b}\frac{\pi }{2}\right)}{\left|det\left(\stackrel{~}{M}_{\gamma _b}I\right)\right|^{1/2}}\right)\frac{\mathrm{exp}i\left(\delta _\mathrm{\Gamma }\frac{\pi }{2}\frac{\pi }{4}\right)}{\sqrt{\left|T_{\gamma _a}^{}+T_{\gamma _b}^{}\right|}}.$$
(25)
As mentioned above, we assumed that there is only one energy partition such that both particles have the same period. This will be the case when the period is a monotonic function of energy, which is a typical situation. If the period is a more complicated function of energy, there may be further solutions and if so then one must have a sum over $`(\gamma _a,\gamma _b)`$ for each possible solution of this condition, but we suppress this possibility for notational simplicity.

We obtained Eq. (25) in Ref.  by doing a stationary phase analysis of the direct autoconvolution of Eq. (15). (The phase index $`\nu `$ in Eq. (18) of Ref.  has a different definition than $`\delta _\mathrm{\Gamma }`$ in Eq. (25), but the overall phase is consistent in the two formulas.) The condition of stationary phase immediately implied that the energy must be partitioned so that the periods of the two orbits are the same. The stationary phase integral then introduces a factor of $`\sqrt{\mathrm{}}`$ as well as the sum of the second derivatives of the actions with respect to energy evaluated at the stationary phase energy. This is precisely the first derivatives of the periods with energy. Thus, we have shown how these two different approaches yield consistent results.

Either particle can execute any number of repetitions of its primitive orbit. If particle $`a`$ executes $`l_{\gamma _a}`$ repetitions and particle $`b`$ executes $`l_{\gamma _b}`$ repetitions, then the energy must be partitioned so that $`l_{\gamma _a}T_{\gamma _a}(E_a)=l_{\gamma _b}T_{\gamma _b}(E_b)`$. As discussed above, the volume term remains unchanged as the square of the primitive periods due to the discrete symmetry in the family of orbits. The action of this orbit is $`l_{\gamma _a}S_{\gamma _a}(E_a)+l_{\gamma _b}S_{\gamma _b}(E_b)`$, and similarly for the phase index. The separate monodromy matrices are raised to the appropriate power. The anholonomy term follows from the discussion above Eq. (24) as $`l_{\gamma _a}T_{\gamma _a}^{}(E_a)+l_{\gamma _b}T_{\gamma _b}^{}(E_b)`$. Assuming that the periods are monotonic in the energy, the index $`\delta _\mathrm{\Gamma }`$ is unchanged. Otherwise, it depends on the sign of $`l_{\gamma _a}T_{\gamma _a}^{}+l_{\gamma _b}T_{\gamma _b}^{}`$. Therefore, the previous expression continues to apply, but the actions, phase indices and periods are multiplied by the appropriate value of $`l_{\gamma _a/\gamma _b}`$ and the separate monodromy matrices are raised to the appropriate power. This dependence can also be understood to be already contained in the definition of the various orbit properties.

We observe that the amplitude of (25) is proportional to the product of the amplitudes for the single-particle dynamics. The trace formula for two noninteracting particles contains an additional prefactor of $`i\mathrm{}/\sqrt{2\pi \mathrm{}}`$, a factor involving the derivatives of the periods with respect to energy (and the associated phase index $`\delta `$) and an additional phase factor of $`\pi /4`$. This result generalises to cases where the amplitudes are not given by (16). We simply replace the single-particle amplitudes in large brackets by the equivalent ones for the system under consideration. This can be understood by noting that the only coupling between the particles is as we have described and any further symmetry can be handled within the single-particle phase spaces. This conclusion can also be understood in the convolution picture by simply using the appropriate single-particle amplitudes when doing the stationary phase analysis .

#### II.4.2 Symmetry decomposition: dynamical pseudo-periodic orbits (DPPOs)

As discussed initially by Gutzwiller  and later in more generality by Robbins , in the presence of a discrete symmetry, the fluctuating density of states can be decomposed among the various irreducible representations (irreps). For the two-particle case, this is simply the symmetric (bosonic) and antisymmetric (fermionic) cases. To evaluate the separate densities of states, one must calculate $`g_\pm (E)=\text{Tr}(\widehat{P}_\pm \widehat{G}(E))`$ using the projection operators in (7). The first term of the projection operator results in the standard sum over dynamical periodic orbits (25). There is a factor of $`1/2`$ which indicates that this contribution is simply divided evenly between the symmetric and antisymmetric spectra. It is the second term of the projection operator which requires careful analysis.

The oscillating part of $`\text{Tr}(\widehat{U}\widehat{G}(E))`$ can be expressed in terms of orbits on which particles begin at a point in phase space, evolve for some time $`T`$, are then exchanged using the classical analogue of $`\widehat{U}`$ with the net result that the particles are returned to their initial conditions. We call these orbits pseudo-periodic to distinguish them from the (standard) dynamical periodic orbits discussed earlier. We first define the symplectic mapping $`u`$ corresponding to classical particle exchange as $`u(z_a,z_b)=(z_b,z_a)`$. It has the property that $`u^2`$ is the identity mapping. The combination of time evolution for time $`t`$ and particle exchange maps a phase space point $`z^{}=(z_a^{},z_b^{})`$ to $`z=u\mathrm{\Phi }_tz^{}=(\varphi _tz_b^{},\varphi _tz_a^{})`$. To find orbits which are periodic under these combined operations, we require phase space coordinates $`z^{}`$ and periods $`T`$ such that $`z^{}=u\mathrm{\Phi }_Tz^{}`$. Applying this combined operation twice, we find that $`z^{}=\mathrm{\Phi }_{2T}z^{}`$. This is just the condition for a periodic orbit of period $`2T`$ in the full phase space without particle exchange. So we conclude that the initial coordinate $`z^{}`$ must be on a periodic orbit of the full phase space. However, this condition is still more restrictive since the above considerations also imply that after time $`T`$, particle $`a`$ must be where particle $`b`$ began and vice-versa. This is only possible if the two particles are traversing the *same* periodic orbit, with the same energy and furthermore are exactly half a period out of phase. We shall call this a type-1 dynamical pseudo-periodic orbit (DPPO). There is also the degenerate case where both particles begin and evolve together. This shall be called a type-0 DPPO and is discussed below.

Therefore, the set of possible pseudo-periodic orbits is much more restricted than the set of standard periodic orbits since we have only contributions when both particles are executing the same dynamics. Furthermore, these orbits are isolated and do not come in a 1-parameter family. The existence of families for the standard periodic orbits is due to the freedom in specifying the relative phases of the two motions. We no longer have this freedom. This immediately implies that contributions from pseudo-periodic orbits will be weaker by $`\sqrt{\mathrm{}}`$ because there is one more stationary phase integral to do than for the standard periodic orbits. (This can also be understood from the fact that particle exchange does not conserve the separate energies and so does not commute with $`J`$.) Therefore, the usual Gutzwiller trace formula applies and we use it to determine the actions, periods and stabilities of these isolated pseudo-periodic orbits.

Consider an arbitrary periodic orbit $`\gamma `$ of the one-particle phase space with period $`T_\gamma `$ and choose some arbitrary initial condition on it which we shall call $`z_a^{}`$. To have a pseudo-periodic orbit in the full phase space, we can begin at $`z^{}=(z_a^{},z_b^{}=\varphi _{T_\gamma /2}z_a^{})`$. If we flow for a time $`T_\gamma /2`$ and then exchange the particles, we map $`z^{}`$ onto itself (see Fig. 2). Therefore, the set of pseudo-periodic orbits in the full two-particle phase space is one-to-one with the set of standard periodic orbits in the one-particle phase space. The periods of the pseudo-periodic orbits in the full phase space are one-half of the periods of the corresponding standard periodic orbits in the one-particle phase space. Nevertheless, when evaluating the trace integral we must integrate over all initial conditions on the orbit and this gives a full factor of $`T_\gamma ^0`$ in the amplitude. The actions and phase indices for the pseudo orbit are the same as for the standard orbit; although we integrate for only half the time, both particles are in motion and between them, they execute one full motion of the periodic orbit. The stability matrix in the full phase space requires careful analysis. Let $`\stackrel{~}{M}_\gamma `$ be the stability matrix of the full periodic orbit $`\gamma `$ of the one-particle phase space and $`\stackrel{~}{M}_\gamma ^{}`$ be the stability matrix of the pseudo-periodic orbit $`\gamma ^{}`$ in the full phase space. It is shown in Appendix B that $`det(\stackrel{~}{M}_\gamma ^{}I)=4det(\stackrel{~}{M}_\gamma I)`$, where on each side of the equation $`I`$ is understood to be an appropriately dimensioned unit matrix. We conclude that the contribution of this orbit to the oscillating part of $`\text{Tr}(\widehat{U}\widehat{G}(E))`$ is

$$\stackrel{~}{g}_\varsigma ^\text{d}(E)=\frac{1}{2i\mathrm{}}\frac{T_\gamma ^0}{\sqrt{\left|det(\stackrel{~}{M}_\gamma I)\right|}}\mathrm{exp}i\left(\frac{S_\gamma }{\mathrm{}}\sigma _\gamma \frac{\pi }{2}\right),$$
(26)
where all classical quantities are evaluated at the single-particle energy $`E/2`$. (Recall the symbol $`\varsigma `$ denotes the group element that exchanges the particles.) Apart from the energy dependence and the factor of two in the denominator, this contribution is the same as the one from the corresponding primitive orbit for the single-particle density of states \[Eqs. (15) and (16)\].

As mentioned above, there is also the situation where both particles start at the same point on the orbit and evolve together. Interchanging them at the end trivially returns them to the same coordinates. This pseudo-orbit has action $`2S`$, but should not be confused with the standard dynamical orbit where the two particles start at independent points on the orbit and therefore occur in a 1-parameter family. The fact that we interchange the particles at the end ensures that the pseudo-orbit described here is isolated and does not occur in a family. The two types of orbits share the same action, but the standard orbit has a larger amplitude due to the different $`\mathrm{}`$ prefactor and will tend to dominate. This situation of coexisting contributions with the same action is analogous to a potential system with a reflection symmetry where there is a boundary orbit which contributes to both the identity term in the density of states and also to the reflection term. The only difference here is that the two types of orbits contribute with different powers of $`\mathrm{}`$.

The analysis of the contribution of the type-0 DPPO is similar to above. We state without proof that its amplitude is simply the same as the double repetition of the orbit, again divided by an overall factor of two as discussed in Appendix B. This pattern continues for higher repetitions, where for odd multiples of the action the particles start $`T_\gamma /2`$ out of phase while for even multiples they start in phase and interfere with stronger (in an $`\mathrm{}`$ sense) contributions from the standard dynamical orbits. Apart from the energy dependence and the factor of two in the denominator, the sum over repetitions is the same as for the single-particle density of states.

Thus, we see that the contribution of the pseudo-periodic orbits to the bosonic and fermionic densities of states is precisely the same as the fluctuating density of states of the one-particle spectrum except that it is to be evaluated at half the total energy (since the total energy is partitioned equally between the two particles) and should also be divided by an overall factor of two. In conclusion,

$$\stackrel{~}{\rho }_\pm (E)=\frac{1}{2}\left(\stackrel{~}{\rho }_2(E)\pm \frac{1}{2}\stackrel{~}{\rho }_1\left(\frac{E}{2}\right)\right),$$
(27)
consistent with (10).

### II.5 One-particle dynamics in the full phase space

We now discuss contributions to the resolvent from periodic orbits in the full phase space where one particle executes dynamics while the other particle remains stationary. In particular, suppose that particle $`a`$ is stationary at some point in phase space while particle $`b`$ evolves dynamically on a periodic orbit. We call this a heterogeneous periodic orbit. The structure of such orbits is qualitatively different for potential systems and billiards.

For analytic potentials, the stationary particle must be at some extremum of the potential with zero momentum. In this case, the full heterogeneous orbit is isolated in phase space since a flow in $`J=h(z_a)`$ does not map an initial condition $`z^{}`$ to any new phase space point $`z`$. Therefore, we can use the Gutzwiller trace formula for isolated orbits. In billiards, the stationary particle has zero momentum, but it can be anywhere in the billiard. So rather than being isolated, the heterogeneous orbits appear as $`d`$-dimensional families. This means that we can use the formalism of Ref.  to calculate the amplitude of these orbits.

The symmetry decomposition for heterogeneous orbits is trivial. Since the two particles are executing completely different dynamics, the combination of time evolution and particle exchange, as above, can never return the particles to their initial conditions. This requires an equivalence of the two motions. Thus, the contribution from heterogeneous orbits is simply divided evenly between the symmetric and antisymmetric representations and belongs to the $`\stackrel{~}{\rho }_2(E)`$ term of Eq. (27).

#### II.5.1 Analytic potentials

Suppose particle $`b`$ traverses a periodic orbit $`\gamma `$ with action $`S_\gamma `$, primitive period $`T_\gamma ^0`$, stability matrix $`\stackrel{~}{M}_\gamma `$ and topological index $`\sigma _\gamma `$. Particle $`a`$ is assumed to be stationary at a potential minimum with energy $`E_a=0`$. At the minimum, the potential is locally harmonic with $`d`$ frequencies $`\omega _j`$. As explained above, the full heterogeneous orbit is isolated and so we can simply use the Gutzwiller trace formula for isolated orbits. The only required information is the monodromy matrix in the phase space of particle $`a`$ since $`det(\stackrel{~}{M}_\mathrm{\Gamma }I)=det(M_aI)det(\stackrel{~}{M}_\gamma I)`$, where $`\stackrel{~}{M}_\mathrm{\Gamma }`$ is the $`(4d2)\times (4d2)`$ stability matrix of the full heterogeneous orbit and $`M_a`$ is the $`2d\times 2d`$ monodromy matrix of particle $`a`$. Since the dynamics of particle $`a`$ are locally harmonic, we can use the result for a $`d`$-dimensional harmonic oscillator (see Appendix C), $`\sqrt{|det(M_aI)|}=_{j=1}^d2\mathrm{sin}(\omega _jT_\gamma (E)/2)`$. The phase index of this motion is simply $`d`$, one for each transverse harmonic degree of freedom. Thus, the contribution of one heterogeneous orbit to the resolvent is

$$\begin{array}{cc}\hfill \stackrel{~}{g}_\mathrm{\Gamma }^\text{h}(E)& =\frac{1}{i\mathrm{}}\frac{T_\gamma ^0(E)}{\sqrt{\left|det\left(\stackrel{~}{M}_\gamma I\right)\right|}_{j=1}^d2\mathrm{sin}\left(\frac{\omega _jT_\gamma (E)}{2}\right)}\hfill \\ & \times \mathrm{exp}i\left(\frac{S_\gamma (E)}{\mathrm{}}\sigma _\gamma \frac{\pi }{2}d\frac{\pi }{2}\right),\hfill \end{array}$$
(28)
where we have retained the symbol $`\mathrm{\Gamma }`$ to denote the full heterogeneous orbit and $`\gamma `$ to stress that this is the contribution from the situation where only one particle is evolving dynamically. There is also an identical contribution from the situation where particle $`b`$ is fixed while particle $`a`$ evolves dynamically. As before, repetitions can be understood to be implicit in the definitions of the action, period, phase index and stability matrix.

One can also consider extrema other than potential minima, such as saddles or potential maxima. We can expand the $`d`$-dimensional potential around an extremum $`x_0`$ as

$$V(xx_0)=\frac{1}{2}\left(\underset{j=1}{\overset{d_+}{}}\omega _j^2\xi _j^2\underset{j=d_++1}{\overset{d}{}}\omega _j^2\xi _j^2\right),$$
(29)
where the $`\xi `$ measure the deviations of $`x`$ from $`x_0`$. In general, there are $`d_+`$ stable directions and $`d_{}=dd_+`$ unstable directions. Then, the expression (28) is still valid, but the energy of the dynamically evolving particle is replaced by $`EV(x_0)`$, the phase factor $`d\pi /2`$ replaced by $`d_+\pi /2`$ and the $`\mathrm{sin}(\omega _jT_\gamma /2)`$ replaced by $`\mathrm{sinh}(\omega _jT_\gamma /2)`$ for the unstable directions. Finally, we note that for smooth potentials, the dynamical orbits give the leading-order contribution to $`\stackrel{~}{\rho }_2(E)`$ while the hetero-orbits give corrections of higher order in $`\mathrm{}`$.

#### II.5.2 Billiard systems

As mentioned above, heterogeneous orbits in a $`d`$-dimensional billiard occur in $`d`$-dimensional families and we may therefore use Eq. (21) with $`f=d`$ to determine the appropriate trace formula. We first consider a two-dimensional billiard ($`d=2`$), although the result is easily generalised. For this case, the orbits do not appear as three-tori, but rather have the topology of $`\times 𝒮^1`$, where $``$ denotes the billiard domain and $`𝒮^1`$ is the one-torus associated with the dynamics of the evolving particle $`b`$ on the periodic orbit $`\gamma `$. There are two constants of the motion: $`J_1=p_{x_a}`$ and $`J_2=p_{y_a}`$ ($`𝐉=(J_1,J_2)`$) and these generate flows $`𝚯=(x_a,y_a)`$. Clearly,

$$det\left(\frac{𝚯}{𝐉}\right)=det\left(\begin{array}{cc}\frac{x_a}{p_{x_a}}& \frac{x_a}{p_{y_a}}\\ & \\ \frac{y_a}{p_{x_a}}& \frac{y_a}{p_{y_a}}\end{array}\right)=\frac{x_a}{p_{x_a}}\frac{y_a}{p_{y_a}},$$
(30)
since the off-diagonal elements vanish due to the fact that the $`x`$ and $`y`$ motions are uncoupled. After particle $`b`$ has traversed the primitive orbit $`n_\gamma `$ times, $`x_a/p_{x_a}=y_a/p_{y_a}=n_\gamma T_\gamma ^0(E)/m`$, where $`T_\gamma ^0(E)`$ is the primitive period of the orbit and $`m`$ is the mass of the particle. (The minus sign indicates that a backwards flow is required to close the orbits in the full phase space.) This immediately implies that the phase index $`\delta 0`$. The stability matrix defined in (21) in this case is simply the stability matrix of the motion of particle $`b`$. The volume for a family of such orbits is the area of the billiard and combining all of the factors, the leading-order contribution of a family of heterogeneous orbits $`\mathrm{\Gamma }`$ to the two-particle density of states is

$$\stackrel{~}{\rho }_\mathrm{\Gamma }^\text{h}(E)=\frac{\alpha 𝒜}{4\pi ^2}\frac{\mathrm{cos}\left(\frac{S_\gamma (E)}{\mathrm{}}\sigma _\gamma \frac{\pi }{2}\frac{\pi }{2}\right)}{n_\gamma \sqrt{\left|det(\stackrel{~}{M}_\gamma I)\right|}}.$$
(31)
We obtained this expression in Ref.  by doing a direct energy convolution integral of the first term of Eq. (13) with Eq. (15). This once again underlines the equivalence of the two methods. This naturally extends to the higher order terms of (13) through a more careful analysis of the surface corrections, but we do not pursue this analysis here. Also, this result generalises to $`d`$ dimensions as

$$\begin{array}{cc}\hfill \stackrel{~}{\rho }_\mathrm{\Gamma }^\text{h}(E)& =\frac{1}{\pi \mathrm{}}\left(\frac{\mathrm{}\alpha }{4\pi }\right)^{d/2}\frac{T_\gamma ^0(E)\mathrm{\Omega }_d}{\sqrt{\left|det(\stackrel{~}{M}_\gamma I)\right|}\left(n_\gamma T_\gamma ^0(E)\right)^{d/2}}\hfill \\ & \times \mathrm{cos}\left(\frac{S_\gamma (E)}{\mathrm{}}\sigma _\gamma \frac{\pi }{2}d\frac{\pi }{4}\right),\hfill \end{array}$$
(32)
where $`\mathrm{\Omega }_d`$ is the $`d`$-dimensional volume of the billiard. We stress that this is $`O(1/\mathrm{}^{d/2})`$ stronger than an isolated orbit, this factor arising from the fact that this class of orbits appear in $`d`$-dimensional families. The contribution from hetero-orbits is also $`O(1/\mathrm{}^{(d1)/2})`$ stronger than the contribution from dynamical orbits. Thus, for billiards, hetero-orbits give the leading-order contribution to $`\stackrel{~}{\rho }_2(E)`$ while dynamical orbits give corrections of higher order in $`\mathrm{}`$.

## III Several Noninteracting Identical Particles

We now consider the extension to $`N`$ identical particles. The smooth term can be written as an $`(N1)`$-fold convolution integral of the single-particle smooth terms and can also be understood as a single integral in the $`N`$-particle phase space. At this point, we say no more about the smooth term and refer the reader to the appendices for further discussions. For the oscillating term, there are again two possibilities. Either all of the particles are evolving dynamically, or a subset of them is stationary at various potential extrema (or anywhere in a billiard). For the first situation, the discussion closely parallels the two-particle case. The only nontrivial quantity to determine is the anholonomy matrix $`\left(𝚯/𝐉\right)_\mathrm{\Gamma }`$. We consider the case $`N=3`$, but this result readily generalises.

### III.1 Dynamical orbits

In an obvious extension of the notation, there are three single-particle phase spaces with coordinates $`z_a`$, $`z_b`$ and $`z_c`$ so that the full three-particle phase space has coordinates $`z=(z_a,z_b,z_c)`$ and the total Hamiltonian $`H(z)=h(z_a)+h(z_b)+h(z_c)`$. Two other constants of motion which are in involution with $`H`$ are $`J_a=h(z_a)`$ and $`J_b=h(z_b)`$ and they generate time translations of particles $`a`$ and $`b`$, respectively, while having no effect on the other particles. Flows generated by $`H`$, $`J_a`$ and $`J_b`$ are denoted by $`\mathrm{\Phi }_t`$, $`\mathrm{\Lambda }_{\theta _a}`$ and $`\mathrm{\Psi }_{\theta _b}`$, respectively. If $`\varphi `$ is a single-particle flow, then flows in the full phase space are mapped as follows:

$`\mathrm{\Phi }_t(z_a,z_b,z_c)`$ $`=`$ $`(\varphi _tz_a,\varphi _tz_b,\varphi _tz_c),`$
$`\mathrm{\Lambda }_{\theta _a}(z_a,z_b,z_c)`$ $`=`$ $`(\varphi _{\theta _a}z_a,z_b,z_c),`$
$`\mathrm{\Psi }_{\theta _b}(z_a,z_b,z_c)`$ $`=`$ $`(z_a,\varphi _{\theta _b}z_b,z_c).`$ (33)
The periodic orbits of the full phase space (at a given total energy $`E`$) can be found from the one-particle periodic orbits by balancing the energy partition among the three particles (i.e. varying $`J_a`$ and $`J_b`$ while holding $`H`$ fixed) so that all the one-particle periodic orbits have the same period. (The result is a three-particle periodic orbit in the full phase space.) Imagine a slight departure from this equilibrium situation so that $`J_aJ_a+\mathrm{\Delta }J_a`$ while holding $`J_b`$ and $`H`$ fixed. Then,

$$\begin{array}{ccc}E_a\hfill & & E_a+\mathrm{\Delta }J_a\hfill \\ E_b\hfill & & E_b\hfill \\ E_c\hfill & & E_c\mathrm{\Delta }J_a\hfill \end{array}\begin{array}{ccc}T_a\hfill & & T_a+\mathrm{\Delta }T_a\hfill \\ T_b\hfill & & T_b\hfill \\ T_c\hfill & & T_c+\mathrm{\Delta }T_c\hfill \end{array}$$
(34)
where $`\mathrm{\Delta }T_a=T_a^{}\mathrm{\Delta }J_a`$ and $`\mathrm{\Delta }T_c=T_c^{}\mathrm{\Delta }J_a`$, the primes denoting differentiation with respect to energy. The initial condition $`z^{}=(z_a^{},z_b^{},z_c^{})`$ with these modified energies (but each particle still on its periodic orbit at that modified energy) is not on a periodic orbit of the full phase space. However, it is on a generalised periodic orbit. That is, the trajectory can be made to close with additional flows in $`(H,J_a,J_b)`$. Imagine a flow in $`H`$ for the original period $`T`$. The orbits of particles $`a`$ and $`c`$ will fail to close by the amount by which their period is longer (or shorter) due to the changed energy: $`\mathrm{\Phi }_Tz^{}=(\varphi _{\mathrm{\Delta }T_a}z_a^{},z_b^{},\varphi _{\mathrm{\Delta }T_c}z_c^{})`$. Additional flows in $`(H,J_a,J_b)`$ close the trajectory. First, a flow in $`H`$ by the amount $`\mathrm{\Delta }T_c`$ returns particle $`c`$ to $`z_c^{}`$: $`\mathrm{\Phi }_{\mathrm{\Delta }T_c}\mathrm{\Phi }_Tz^{}=(\varphi _{\mathrm{\Delta }T_a+\mathrm{\Delta }T_c}z_a^{},\varphi _{\mathrm{\Delta }T_c}z_b^{},z_c^{})`$. The condition for a periodic orbit $`\mathrm{\Lambda }_{\mathrm{\Delta }\theta _a}\mathrm{\Psi }_{\mathrm{\Delta }\theta _b}\mathrm{\Phi }_{\mathrm{\Delta }T_c}\mathrm{\Phi }_Tz^{}=z^{}`$ immediately implies

$`\mathrm{\Delta }\theta _a`$ $`=`$ $`(T_a^{}+T_c^{})\mathrm{\Delta }J_a,`$
$`\mathrm{\Delta }\theta _b`$ $`=`$ $`T_c^{}\mathrm{\Delta }J_a.`$ (35)
We get a similar result from a deviation in $`J_b`$ (holding $`J_a`$ and $`H`$ fixed) and conclude that

$$\left(\frac{𝚯}{𝐉}\right)=\left(\begin{array}{cc}T_a^{}+T_c^{}& T_c^{}\\ T_c^{}& T_b^{}+T_c^{}\end{array}\right).$$
(36)
The determinant is $`T_a^{}T_b^{}+T_b^{}T_c^{}+T_c^{}T_a^{}`$ and is invariant under a permutation of the indices. Note that we could have chosen the two generators $`J_a`$ and $`J_c`$ and followed through the analogous calculation. In that calculation, the anholonomy matrix would be modified by permuting $`b`$ and $`c`$ in Eq. (36). Therefore, the eigenvalues of $`𝚯/𝐉`$ are not invariant. But, since the determinant is invariant, so too is the number of positive eigenvalues which determines the phase index $`\delta `$. Therefore, the final result is invariant. For $`N>3`$ particles, this generalises to

$$det\left(\frac{𝚯}{𝐉}\right)=\left(\underset{p=1}{\overset{N}{}}T_p^{}\right)\left(\underset{p=1}{\overset{N}{}}\frac{1}{T_p^{}}\right),$$
(37)
where $`T_p`$ is the period of the orbit on which particle $`p`$ is residing. This can be shown by induction.

The other factors which go into the trace formula are simple to determine; the discussion is very similar to the two-particle case and so we refrain from going into great detail. For $`N`$ particles, flows in $`H`$ and $`𝐉=(J_1,\mathrm{},J_{N1})`$ map out an $`N`$-dimensional torus. This means there are $`(N1)`$-parameter families of periodic orbits in the full phase space. The total action is the sum of all the single-particle actions and similarly for the total phase index $`\mu `$. The monodromy matrix is defined holding all of the single-particle energies constant in such a way that it is block diagonal among the various single-particle motions. The volume of the periodic orbit family is the product of the primitive periods. (To see this, recall that the volume term $`T_\mathrm{\Gamma }V_\mathrm{\Gamma }=_\mathrm{\Gamma }\text{d}t\text{d}\theta _1\text{d}\theta _2\mathrm{}\text{d}\theta _{N1}`$ and that the primitive volume should only count distinct configurations.) Using Eq. (21) with $`f=(N1)`$, we conclude that the contribution to the resolvent from one family of dynamical periodic orbits is

$$\stackrel{~}{g}_\mathrm{\Gamma }^\text{d}(N,E)=\frac{1}{i\mathrm{}}\frac{1}{(2\pi i\mathrm{})^{(N1)/2}}\left\{\underset{p=1}{\overset{N}{}}\frac{T_p^0(E_p)\mathrm{exp}i\left(\frac{S_p(E_p)}{\mathrm{}}\sigma _p\frac{\pi }{2}\right)}{\sqrt{\left|det(\stackrel{~}{M}_pI)\right|}\sqrt{\left|T_p^{}(E_p)\right|}}\right\}\frac{\mathrm{exp}i\delta _\mathrm{\Gamma }\frac{\pi }{2}}{\sqrt{\left|_{p=1}^N\frac{1}{T_p^{}(E_p)}\right|}}.$$
(38)
In Eq. (38), we have used the label $`p`$ rather than the more cumbersome $`\gamma _p`$ to refer to the periodic orbit on which particle $`p`$ resides. We will continue to do this for the remainder of the section. The phase factor $`\delta _\mathrm{\Gamma }`$ is the number of positive eigenvalues of the $`(N1)\times (N1)`$ matrix $`(𝚯/𝐉)_\mathrm{\Gamma }`$. If all of the particles are on distinct orbits, then there are $`N!`$ congruent but distinct full phase space orbits, corresponding to the choice of which particle to assign to which orbit. If there is more than one particle on the same orbit, then the number of combinatoric possibilities is accordingly modified. We take this combinatoric factor to be implicit in the sum over orbits and do not explicitly account for it here.

### III.2 Hetero-orbits

The other possibility is that some of the particles are not evolving dynamically, but rather are stationary in a billiard or at potential extrema. Imagine that $`M`$ particles are evolving dynamically and $`(NM)`$ are fixed at extrema. Then, these heterogeneous orbits come in $`(M1)`$-fold families. In the special case where the nonevolving particles are stationary at potential minima,

$$\stackrel{~}{g}_\mathrm{\Gamma }^\text{h}(M,N,E)=\stackrel{~}{g}_\mathrm{\Gamma }^\text{d}(M,E_e)\left\{\underset{p=M+1}{\overset{N}{}}\frac{\mathrm{exp}\left(i\frac{\pi d}{2}\right)}{_{j=1}^d2\mathrm{sin}\left(\frac{\omega _{j_p}T}{2}\right)}\right\}.$$
(39)
The evolving particles share the energy $`E_e=E_{p=M+1}^NV(x_p)`$, where $`x_p`$ denote the positions of the stationary particles. We recall that $`d`$ is the dimension of the one-particle dynamics and the $`\omega _{j_p}`$ denote the $`d`$ local harmonic frequencies around the minimum at which particle $`p`$ resides. As in the two-particle case, if a particle is at a saddle or maximum, we replace the phase $`d\pi /2`$ by $`d_+\pi /2`$, where $`d_+`$ denotes the number of stable directions and replace the $`\mathrm{sin}`$ in the amplitude by $`\mathrm{sinh}`$ for the unstable directions. Again, there are distinct but congruent heterogeneous orbits in which different particles are chosen to be on different orbits or extrema, but we refrain from an explicit discussion of the combinatoric possibilities.

Next suppose that $`(NM)`$ particles are stationary in a $`d`$-dimensional billiard. In addition to the $`(M1)`$ independent generators that exist for the potential system, there are $`(NM)d`$ generators $`𝐉_𝐪=(𝐩_1,\mathrm{},𝐩_{NM})`$. The associated flows are $`𝚯_𝐪=(𝐪_1,\mathrm{},𝐪_{NM})`$. (Each $`𝐩`$ and $`𝐪`$ is $`d`$-dimensional.) Since the generators associated with the stationary particles also generate new orbits, the dimensionality of the orbit families is $`f=(M1)+(NM)d`$. The volume term $`T_\mathrm{\Gamma }^0V_\mathrm{\Gamma }^0=_\mathrm{\Gamma }\text{d}t\text{d}\theta _1\mathrm{}\text{d}\theta _{M1}\text{d}𝐪_1\mathrm{}\text{d}𝐪_{NM}=T_1^0(E_1)\mathrm{}T_M^0(E_M)\mathrm{\Omega }_d^{(NM)}`$. The phase index $`\delta _\mathrm{\Gamma }`$ is the number of positive eigenvalues of the $`f\times f`$ matrix $`(𝚯/𝐉)_\mathrm{\Gamma }`$ which has a block-diagonal structure; one block is the anholonomy associated with the evolving particles analogous to Eq. (36) and one block is the anholonomy associated with the stationary particles analogous to Eq. (30). Thus, the contribution to the resolvent from a family of billiard hetero-orbits is

$$\stackrel{~}{g}_\mathrm{\Gamma }^\text{h}(M,N,E)=\stackrel{~}{g}_\mathrm{\Gamma }^\text{d}(M,E)\left\{\underset{p=M+1}{\overset{N}{}}\frac{\mathrm{\Omega }_d\mathrm{exp}\left(i\frac{\pi d}{4}\right)}{(\frac{2\pi \mathrm{}T}{m})^{d/2}}\right\}.$$
(40)
In Eqs. (39-40), $`T`$ is the global period (recall that the energies of all the dynamically evolving particles have been partitioned so that all of the periodic orbits have a common period) and $`\delta _\mathrm{\Gamma }0`$ if $`M=1`$. As in the two-particle case, hetero-orbits are more important in billiards than in smooth potentials. Their leading-order contribution to $`\stackrel{~}{\rho }_N(E)`$ is $`O(1/\mathrm{}^{(NM)(d1)/2})`$ stronger for billiards and $`O(\mathrm{}^{(NM)/2})`$ weaker for potentials than the corresponding contribution from the dynamical orbits.

We now make some final comments. The above expressions apply for any of the particles executing multiple repetitions of its primitive orbit provided the energy is partitioned among the dynamically evolving particles so that all single-particle periodic orbits have a common period. Then, the various orbit properties which appear in the formulas are understood to be those for the repeated orbit. The formulas written above only account for the contribution of a single family of orbits. The oscillatory part of the resolvent is a sum over all families: $`\stackrel{~}{g}(E)=_\mathrm{\Gamma }\stackrel{~}{g}_\mathrm{\Gamma }(E)`$. We mention that Eqs. (38-40) can also be obtained from convolution integrals by doing a stationary phase analysis of the $`N`$-particle dynamical term and taking appropriate combinations of saddle-point and end-point contributions from the various cross-term integrals. However, the approach outlined above is more illuminating since it reveals the underlying structure of the periodic orbit families. The many-particle trace formulas involve only properties of periodic orbits of the one-particle phase space. Thus, after studying a one-particle system, one can immediately work out the details of the many-particle system. This parallels the situation in quantum mechanics where the problem of $`N`$ noninteracting particles in a potential is a simple extension of the one-particle problem.

## IV Symmetry decomposition of the $`N`$-particle density of states

If the system consists of $`N`$ identical particles, it is invariant under $`S_N`$, the permutation group of $`N`$ identical particles. This group has many different irreps for $`N>2`$, but we only consider the one-dimensional bosonic/fermionic irreps which are fully symmetric/antisymmetric under particle exchange. We first introduce the projection operators

$$\widehat{P}_\pm =\frac{1}{N!}\underset{\tau }{}(\pm 1)^{n_\tau }\widehat{U}_\tau ,$$
(41)
where $`\pm `$ refer to the bosonic/fermionic irreps, respectively. The sum is over the group elements $`\tau `$ of $`S_N`$ which denote particular permutations of the particles, $`\widehat{U}_\tau `$ is the representation of the group element in the Hilbert space (i.e. the quantum operator which exchanges the particles), $`n_\tau `$ is the number of 2-particle exchanges required to obtain $`\tau `$ and the factor $`(\pm 1)^{n_\tau }`$ is a group character. For fermions, the sign of the character depends on the number of times two particles must be interchanged. As before, we need to evaluate $`g_\pm (E)=\text{Tr}(\widehat{P}_\pm \widehat{G}(E))`$ and therefore $`\text{Tr}(\widehat{U}_\tau \widehat{G}(E))`$ for each $`\tau `$. This is a class function, only depending on the cyclic structure of $`\tau `$.

Consider a permutation and break it up into cycles . For $`N`$ particles, $`\tau `$ can be decomposed uniquely into mutually commuting cycles; in each of these cycles, a subset of the particles is being permuted. An $`n`$-cycle is a permutation in which only $`n`$ of the particles are being permuted. In particular, a 1-cycle corresponds to an individual particle being left alone, a 2-cycle corresponds to two particles being exchanged with each other and so on. A general permutation $`\tau `$ may consist of cycles of various sizes and also may have several cycles of the same size. In general, for a given $`\tau `$, there are $`\nu _1`$ 1-cycles, $`\nu _2`$ 2-cycles and so on. Then, the cycle structure of a class of permutations can be given as a set of integers $`(\nu _1,\nu _2,\mathrm{},\nu _N)`$. This set $`𝝂`$ labels the conjugacy classes. Two permutations with the same $`𝝂`$ belong to the same class and thus have the same value of $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$. The analysis of the previous section can be understood as being the special case of the identity element. To decompose the full density of states, one needs to determine both the smooth and oscillating contributions to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$. The smooth contribution is discussed in Appendix D. In this section, we examine the oscillating contribution.

### IV.1 Dynamical cycles

Consider first the case for which all particles are evolving dynamically. A group element $`\tau `$ consists of $`m_\tau `$ cycles, a given cycle $`k`$ consisting of interchanging $`n_k`$ particles. As in the two-particle case, particle interchange does not commute with all of the single-particle energies and so we do not expect periodic orbit families of dimension $`(N1)`$. However, for each cycle, there is a generator $`J_k`$ which is the sum of the single-particle Hamiltonians of the particles involved in that cycle and is preserved under the action of the group element $`\tau `$. These generators commute with each other and with the total Hamiltonian $`H`$. However, this is not an independent set since $`_kJ_k=H`$. There are $`(m_\tau 1)`$ independent commuting generators other than the full Hamiltonian and so we expect periodic orbit families of this dimensionality contributing to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$.

We seek structures in the full phase space which are invariant under the combined operations of time evolution (for time $`T`$) generated by $`H`$ and particle exchange as specified by $`\tau `$. Clearly, this is only possible if all particles of a given cycle $`k`$ are on the *same* periodic orbit, $`\gamma _k`$. They must all have the same energy, which we shall call $`E_k`$ and then $`J_k=n_kE_k`$. For example, imagine that particles $`a`$, $`b`$ and $`c`$ constitute a 3-cycle. Starting with particle $`a`$ at some arbitrary point on a periodic orbit $`\gamma `$ of the one-particle phase space, particle $`b`$ an amount $`T_\gamma /3`$ ahead of it and particle $`c`$ an amount $`T_\gamma /3`$ behind. Then, after a time $`T=T_\gamma /3`$, $`ab`$, $`bc`$ and $`ca`$. However, the group element $`\tau =(acb)`$ maps $`ac`$, $`cb`$ and $`ba`$ simply undoes this change and the original configuration is restored. Such a cycle is shown at the left of Fig. 3. We then imagine that for every cycle comprising $`\tau `$, there is a train of particles with identical energies traversing a periodic orbit of the one-particle phase space. Each particle completes $`(1/n_k)`$th of the full motion on the periodic orbit.

The logic then is as follows. We assign each cycle a periodic orbit $`\gamma _k`$. (We will henceforth label the orbit properties using the subscript $`k`$ rather than the more cumbersome $`\gamma _k`$.) We partition the energy (i.e. the values of $`J_k`$) so that the values of the periods $`T_k/n_k`$ are all the same; this quantity we denote by $`T`$. After time $`T`$ and permutation $`\tau `$, the resulting structure is guaranteed to be globally periodic in the full phase space. Such an orbit comes in an $`(m_\tau 1)`$ degenerate family which can be understood as follows. For each cycle, it is enough to specify the initial condition of one particle after which we know the initial conditions of all the other particles. We choose the initial condition of the first particle arbitrarily for the first cycle. The first particle of the other $`(m_\tau 1)`$ cycles can then begin anywhere on their respective orbits (this constituting the dimensionality of the family). We also understand this from the fact that starting at the arbitrary initial condition, flows generated by any of the $`(m_\tau 1)`$ generators $`J_k`$ map out a surface of this dimensionality. Together with a flow in $`H`$, the periodic orbit surface is a torus of dimension $`m_\tau `$.

For the symmetry decomposition (involving the dynamical orbits) of a two-particle system, it was noted that there were contributions from higher multiples. For instance, one could start both particles at the same point on a single-particle orbit, let them evolve for a full period and then interchange them. There is an analogous structure in the $`N`$-particle case. We can allow the particles to execute a fraction $`l_k/n_k`$ of an orbit as depicted in Fig. 4. As before, the additional factor $`l_k`$ can be absorbed into the definitions of the various classical parameters.

The contribution of an $`m_\tau `$-torus of orbits can now be inferred from our previous work. The only detail is in the determination of $`\left(𝚯/𝐉\right)`$. It is as in (37), but with the understanding that the sum/product over orbits should be replaced by a sum/product over cycles. These become equivalent in the identity contribution which was considered there. Also, since the anholonomy term measures deviations away from global periodicity arising from a change in the energy partition (now among the cycles), $`T_p^{}`$ should be replaced by $`T_k^{}/n_k^2`$. A factor of $`1/n_k`$ comes from the fact that the energy of the cycle must be divided evenly among the $`n_k`$ particles belonging to that cycle. A second factor of $`1/n_k`$ comes from the fact that the orbit has time $`T_k/n_k`$ for the anholonomy to evolve. (Note that if this orbit is a multiple repeat, then it is understood that $`T_k^{}=l_kT_k^0^{}`$, where $`T_k^0`$ is the primitive period.) The entire contribution should also be divided by $`_kn_k`$ arising from the monodromy matrix as discussed in Appendix B. This last fact is the generalisation of the factor of $`1/2`$ appearing as a prefactor in the second term of Eq. (27) for the two-particle case. Therefore, the contribution from a family of dynamical cycles to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$ can be written as

$$\stackrel{~}{g}_\tau ^\text{d}(m_\tau ,E)=\frac{1}{i\mathrm{}}\frac{1}{(2\pi i\mathrm{})^{(m_\tau 1)/2}}\left\{\underset{k=1}{\overset{m_\tau }{}}\frac{T_k^0(E_k)\mathrm{exp}i\left(\frac{S_k(E_k)}{\mathrm{}}\sigma _k\frac{\pi }{2}\right)}{\sqrt{\left|det(\stackrel{~}{M}_kI)\right|}\sqrt{\left|T_k^{}(E_k)\right|}}\right\}\frac{\mathrm{exp}i\delta \frac{\pi }{2}}{\sqrt{\left|_{k=1}^{m_\tau }\frac{n_k^2}{T_k^{}(E_k)}\right|}},$$
(42)
where $`\stackrel{~}{M}_k`$ is the stability matrix for a full cycle $`k`$ (cf. Appendix B). Notice that the contribution of the group element for which *all* of the particles belong to the same cycle is proportional to $`\stackrel{~}{\rho }_1(E/N)`$ <sup>1</sup><sup>1</sup>1The division by $`N`$ of the energy argument simply states that the total energy must be *evenly* divided among all of the particles. The set of orbits corresponding to this cycle is clearly the same as the set of orbits of the one-particle dynamics (almost by definition) and the amplitudes and actions are the same as the one-particle case since the $`N`$ particles *collectively* execute one complete motion (or a multiple repetition) of the periodic orbit. This is observed for $`N=3`$ in section V.2..

### IV.2 Hetero-cycles

It is also possible that $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$ has a contribution from cycles where some particles are fixed (either at extrema of the potential or anywhere in a billiard) while others are evolving dynamically. Let $`s`$ denote the number of cycles that are stationary and $`e`$ the number of cycles that are evolving dynamically. Then, $`s+e=m_\tau `$. To have such a contribution to the oscillating component, group elements must consist of two or more cycles since those that consist of only one $`n_k`$-cycle will either contribute to Eq. (42) if they are dynamical cycles or to the smooth part (cf. Appendix D) if they are stationary cycles. In addition, we require at least one cycle to involve particles which are evolving dynamically ($`e1`$) *and* at least one cycle to involve particles which are stationary ($`s1`$). Thus, hetero-cycles are cycles for which $`1e<N`$ and $`1s<N`$.

For potentials, the dimension of a family of orbits is then $`(e1)`$ since only the generators associated with dynamical cycles generate new orbits. The stationary cycles simply contribute their monodromy matrices and phase indices and otherwise play no essential role. Eq. (42) holds for the particles which are evolving dynamically, but $`m_\tau `$ is replaced by $`e`$ and the energy associated with the $`e`$ dynamical cycles, $`E_e`$, is the total energy minus the sum of the potential energies of the stationary particles. For a potential minimum, the contribution of one family of hetero-cycles to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$ is

$$\stackrel{~}{g}_\tau ^\text{h}(e,m_\tau ,E)=\stackrel{~}{g}_\tau ^\text{d}(e,E_e)\left\{\underset{k=e+1}{\overset{m_\tau }{}}\frac{\mathrm{exp}\left(i\frac{\pi d}{2}\right)}{_{j=1}^d2\mathrm{sin}\left(\frac{\omega _{j_k}T_k}{2}\right)}\right\},$$
(43)
where the $`\omega _{j_k}`$ denote the local frequencies around the potential minimum at which the particles of cycle $`k`$ reside. If this cycle of particles is actually at a saddle or a maximum, the final factor is modified as in the discussion below Eq. (29).

For billiards, the previous relation holds for the dynamical cycles, but the product over stationary cycles is modified. As explained below, the dimension of the orbit families is $`[(e1)+sd]`$ since the generators associated with stationary cycles also generate new orbits. If there are $`s`$ stationary 1-cycles, these generators and their flows are $`𝐉=(𝐩_1,\mathrm{},𝐩_s)`$ and $`𝚯=(𝐪_1,\mathrm{},𝐪_s)`$, respectively. (There are $`sd`$ components since each $`𝐩_i`$ and $`𝐪_i`$ is $`d`$-dimensional.) In fact, this is true regardless of the number of particles belonging to the stationary cycle. At first, this may seem incorrect since longer cycles will introduce additional generators because they involve more particles. However, this larger set of generators is not an independent set. To see this, recall that the particles involved in a stationary cycle can be anywhere in the billiard. If the cycle is not a 1-cycle, but rather an $`n_k`$-cycle, the combined operations of time evolution and particle exchange will not restore the initial configuration unless all the particles involved in that cycle possess the same phase space coordinates. More formally, a stationary $`n_k`$-cycle possesses a set of generators, $`𝐉=(𝐩_1,\mathrm{},𝐩_{n_k}`$) and associated flows $`𝚯=(𝐪_1,\mathrm{},𝐪_{n_k}`$), where each $`𝐩_i`$ and $`𝐪_i`$ has $`d`$ components. However, after the specification of a single ($`𝐩_i`$,$`𝐪_i`$) pair, all others are uniquely determined: $`𝐩_1=𝐩_2=\mathrm{}=𝐩_{n_k}`$ *and* $`𝐪_1=𝐪_2=\mathrm{}=𝐪_{n_k}`$. Thus, one independent set of generators is $`𝐉=𝐩_i/n_k`$, where $`𝐩_i`$ is the momentum of the $`i`$th particle of the stationary cycle. The factor of $`1/n_k`$ comes from the fact that the momentum of the cycle must be equally partitioned among the $`n_k`$ particles belonging to that cycle. We note here that it is not necessary for stationary particles of distinct cycles to have the same phase space coordinates. Thus, for a $`d`$-dimensional billiard, the contribution to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$ from a family of hetero-cycles is

$$\stackrel{~}{g}_\tau ^\text{h}(e,m_\tau ,E)=\stackrel{~}{g}_\tau ^\text{d}(e,E)\left\{\underset{k=e+1}{\overset{m_\tau }{}}\frac{\mathrm{\Omega }_d\mathrm{exp}\left(i\frac{\pi d}{4}\right)}{\left(\frac{2\pi \mathrm{}n_kT_k}{m}\right)^{d/2}}\right\}.$$
(44)
Eqs. (43-44) are the most general formulas of the paper. We allow for any amount of particle permutation and any number of particles can be evolving while the rest are stationary. Each cycle can involve an arbitrary repetition of the primitive motion. As before, if $`e=1`$, then $`\delta 0`$. Hetero-cycles in billiards are $`O(1/\mathrm{}^{sd/2})`$ stronger than in smooth potentials. For potentials \[billiards\], hetero-cycles are $`O(\mathrm{}^{s/2})`$ weaker \[$`O(1/\mathrm{}^{(d1)s/2})`$ stronger\] than dynamical cycles. Thus, the most significant structures (in an $`\mathrm{}`$ sense) are dynamical cycles for smooth potentials and hetero-cycles for billiards.

## V Numerical Example: The Three-Particle Cardioid Billiard

To illustrate the use of the trace formulas derived above, we study a system of three noninteracting identical particles in a two-dimensional cardioid billiard. In a billiard, classical orbits possess simple scaling properties. For instance, the action and period of an orbit $`\gamma `$ are

$`S_\gamma (\epsilon )`$ $`=`$ $`\sqrt{2m\epsilon }L_\gamma ,`$
$`T_\gamma (\epsilon )=S_\gamma ^{}(\epsilon )`$ $`=`$ $`{\displaystyle \frac{\sqrt{2m}L_\gamma }{2\sqrt{\epsilon }}}={\displaystyle \frac{\mathrm{}\sqrt{\alpha }}{2\sqrt{\epsilon }}}L_\gamma .`$ (45)
For this reason, it is natural and convenient to analyse the length spectrum of the various trace formulas. Thus, for our analysis, we shall compare Fourier transforms of quantum spectra with their semiclassical approximations in the reciprocal space of orbit lengths, $`L`$. In reciprocal $`L`$-space, we expect signals at the lengths of the full periodic orbits of the three-particle system. In the subsequent analysis, peaks in the various length spectra are identified with particular periodic orbits of the full classical phase space. We first consider the total (three-particle) density of states for the cardioid system and then study its decomposition among the irreps of $`S_3`$.

### V.1 Total density of states

#### V.1.1 Quantum mechanics

The analogue of Eqs. (II.1) and (5) for the quantum three-particle density of states is

$`\rho _3(E)`$ $`=`$ $`{\displaystyle \underset{i,j,k}{}}\delta (E(ϵ_i+ϵ_j+ϵ_k))`$ (46)
$`=`$ $`\rho _1(E)\rho _1(E)\rho _1(E).`$
In fact, this relation applies even if the particles are not identical where the full density is still the convolution of the three distinct single-particle densities. We construct the three-particle spectrum by adding the energies of the one-particle spectrum. (The billiard has a reflection symmetry which implies that all the single-particle states are either even or odd; this symmetry should not be confused with the symmetry due to particle exchange.) In the subsequent analysis, we work exclusively with the odd-parity one-particle spectrum. We include the first 500 single-particle energies which allows us to construct the first 19,317,062 energy levels representing all three-particle energies less than $`2.8148\times 10^3`$. (The spectrum is truncated at $`E_{\text{max}}=2ϵ_1+ϵ_{500}`$ to ensure there are no missing levels.) It is possible to improve the resolution in $`L`$-space by truncating the spectrum at a higher energy. But, this would require a precise spectrum since there is a rapid increase in the number of three-particle levels with energy and errors accumulate.

#### V.1.2 Weyl term

Using the identity contribution from Appendix D, the smooth three-particle density of states is just the 2-fold convolution integral of the smooth single-particle density of states:

$`\overline{\rho }_3(E)=\overline{\rho }_1(E)\overline{\rho }_1(E)\overline{\rho }_1(E).`$ (47)
For a two-dimensional billiard, we use Eq. (13) for $`\overline{\rho }_1(E)`$. After performing the necessary integrations (ignoring terms $`O(1/\mathrm{}^3)`$), the three-particle smooth term is found to be

$$\begin{array}{cc}\hfill \overline{\rho }_3(E)& =\frac{\alpha ^3𝒜^3}{128\pi ^3}E^2\frac{\alpha ^{5/2}𝒜^2}{32\pi ^3}E^{3/2}\hfill \\ & +\frac{3}{2}\alpha ^2\left(\frac{𝒜^2}{128\pi ^2}+\frac{𝒜^2𝒦}{16\pi ^2}\right)E.\hfill \end{array}$$
(48)
For the odd-parity single-particle spectrum of the cardioid, $`𝒜=3\pi /4`$, $`=6`$, and $`𝒦=3/16`$. Some of the contributions of the higher-order terms of $`\overline{\rho }_3(E)`$ can be calculated, but it is formally meaningless to include them since there are corrections of the same relative order in $`\mathrm{}`$ that are not known. The terms that are $`O(\sqrt{\alpha ^3E})`$ and $`O(\alpha E^0)`$ can be computed numerically.

#### V.1.3 Hetero-orbits

For three particles in a two-dimensional billiard, there are two types of heterogeneous orbits. The first type occurs when one particle is on a periodic orbit while the other two particles are stationary. These orbits come in 4-parameter families. The trace formula is obtained by using Eq. (40) with $`M=1,N=3`$. For the situation where particles $`a`$ and $`b`$ are stationary and particle $`c`$ evolves on the orbit $`\gamma `$, the leading-order contribution to $`\stackrel{~}{\rho }_3(E)`$ is

$$\begin{array}{cc}\hfill \stackrel{~}{\rho }_3^{\text{h1}}(E)& =\frac{\alpha ^{3/2}𝒜^2E^{1/2}}{8\pi ^3}\underset{\gamma }{}\frac{(L_\gamma ^0/L_\gamma ^2)}{\sqrt{\left|det(\stackrel{~}{M}_\gamma I)\right|}}\hfill \\ & \times \mathrm{cos}\left(\sqrt{\alpha E}L_\gamma \sigma _\gamma \frac{\pi }{2}\pi \right).\hfill \end{array}$$
(49)
The second type of hetero-orbit arises from the situation where only one particle is stationary while the other two evolve on periodic orbits. For instance, particle $`c`$ is stationary while particle $`a`$ evolves on $`\gamma _a`$ and particle $`b`$ evolves on $`\gamma _b`$. Using formula (40) with $`M=2,N=3`$, we conclude the leading-order contribution to $`\stackrel{~}{\rho }_3(E)`$ from these hetero-orbits is

$$\begin{array}{cc}& \stackrel{~}{\rho }_3^{\text{h2}}(E)=\frac{\alpha ^{5/4}𝒜E^{1/4}}{(2\pi )^{5/2}}\underset{\gamma _a,\gamma _b}{}\left(\underset{p=a,b}{}\frac{L_{\gamma _p}^0}{\sqrt{\left|det(\stackrel{~}{M}_{\gamma _p}I)\right|}}\right)\hfill \\ & \times \left(L_{\gamma _a}^2+L_{\gamma _b}^2\right)^{3/4}\mathrm{cos}\left(\sqrt{\alpha E}L_\mathrm{\Gamma }\sigma _\mathrm{\Gamma }\frac{\pi }{2}\frac{3\pi }{4}\right),\hfill \end{array}$$
(50)
where $`L_\mathrm{\Gamma }=\sqrt{L_{\gamma _a}^2+L_{\gamma _b}^2}`$ and $`\sigma _\mathrm{\Gamma }=(\sigma _{\gamma _a}+\sigma _{\gamma _b})`$.

For the total density of states, both formulas are multiplied by a factor of $`3`$ since there are three identical contributions depending on the choice of which particle is evolving and which is stationary. Higher order contributions can be obtained using the convolution formalism and the results are given in Appendix E.

#### V.1.4 Dynamical orbits

To use formula (38), we must first determine the energies that satisfy the following conditions:

$`T_{\gamma _a}(E_a)=T_{\gamma _b}(E_b)=T_{\gamma _c}(E_c),`$ (51)
$`E_a+E_b+E_c=E.`$
This leads to a simple linear system which can be solved to give

$`E_i=\left({\displaystyle \frac{L_{\gamma _i}^2}{L_{\gamma _a}^2+L_{\gamma _b}^2+L_{\gamma _c}^2}}\right)E`$ (52)
for $`i=a,b,c`$. We can now proceed to compute each of the quantities involved in formula (38). The anholonomy term (cf. Eq. (36))

$`det\left({\displaystyle \frac{𝚯}{𝐉}}\right)`$ $`=`$ $`T_{\gamma _a}^{}T_{\gamma _b}^{}+T_{\gamma _b}^{}T_{\gamma _c}^{}+T_{\gamma _c}^{}T_{\gamma _a}^{}`$ (53)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2\alpha }{16}}{\displaystyle \frac{L_\mathrm{\Gamma }^8}{L_{\gamma _a}^2L_{\gamma _b}^2L_{\gamma _c}^2E^3}}.`$
In addition, $`\text{Tr}\left(\frac{𝚯}{𝐉}\right)<0`$ and this implies the phase factor $`\delta _\mathrm{\Gamma }0`$. Then, the three-particle dynamical term can be written as

$$\begin{array}{cc}\hfill \stackrel{~}{\rho }_3^\text{d}(E)& =\frac{\alpha }{\left(2\pi \right)^2}\underset{\gamma _a,\gamma _b,\gamma _c}{}\frac{1}{L_\mathrm{\Gamma }}\left(\underset{p=a,b,c}{}\frac{L_{\gamma _p}^0}{\sqrt{\left|det(\stackrel{~}{M}_{\gamma _p}I)\right|}}\right)\hfill \\ & \times \mathrm{cos}\left(\sqrt{\alpha E}L_\mathrm{\Gamma }\sigma _\mathrm{\Gamma }\frac{\pi }{2}\frac{\pi }{2}\right),\hfill \end{array}$$
(54)
where $`L_\mathrm{\Gamma }=\sqrt{L_{\gamma _a}^2+L_{\gamma _b}^2+L_{\gamma _c}^2}`$ and $`\sigma _\mathrm{\Gamma }=(\sigma _{\gamma _a}+\sigma _{\gamma _b}+\sigma _{\gamma _c})`$.

#### V.1.5 Numerics

We first mention that for billiards, it is common to express the density of states in terms of the wavenumber $`k`$, where $`\epsilon =k^2/\alpha `$ so that $`\rho (k)=2k\rho (\epsilon )/\alpha `$. This is convenient since $`k`$ is conjugate to the periodic orbit lengths $`L`$. Therefore, our numerical results will be quoted as functions of $`k`$ with the understanding that these functions have been converted to the $`k`$ domain from the energy domain using the Jacobian relation above. This will always be the case when the argument is $`k`$. As well, for all numerical comparisons, $`\alpha `$ (and $`\mathrm{}`$) are set to unity.

We compare the Fourier transform of the oscillatory part of the density of states

$`\stackrel{~}{F}_3^{\text{sc}}(L)`$ $`=`$ $`\{\stackrel{~}{\rho }_3(k)\}`$ (55)
$`=`$ $`\{3\stackrel{~}{\rho }_3^{\text{h1}}(k)+3\stackrel{~}{\rho }_3^{\text{h2}}(k)+\stackrel{~}{\rho }_3^\text{d}(k)\}`$
with its quantum mechanical analogue which we define to be

$`\stackrel{~}{F}_3^{\text{qm}}(L)=\{\rho _3(k)\overline{\rho }_3(k)\}.`$ (56)
In Eq. (56), the first term is the quantum three-particle density of states, $`\rho _3(k)=_I\delta \left(kk_I\right)`$, where the superindex $`I`$ denotes a triplet of integers ($`i,j,k`$). The subtracted term is the smooth density of states as determined from Eq. (48). The oscillatory part has contributions from hetero-orbits (49-50) and dynamical orbits (54). In all formulas, $`\gamma _p`$ are periodic orbits in the fundamental domain (i.e. the half-cardioid) and $`L_{\gamma _p}^0`$ are their primitive lengths. Orbit properties are discussed in Refs.  and some of the shorter geometrical orbits are shown in Fig. 5. The stability matrices of the Gutzwiller amplitudes are computed using the standard procedure for the stability of free-flight billiards (see, for example, Ref. ). We define the Fourier transform

$$\{\rho (k)\}=_{\mathrm{}}^{\mathrm{}}\text{w}(k)\mathrm{exp}(ikL)\rho (k)\text{d}k$$
(57)
as a function of the conjugate variable $`L`$. Here, $`\text{w}(k)`$ is the three-term Blackman-Harris window function

$$\text{w}(k)=\{\begin{array}{ccc}_{j=0}^2a_j\mathrm{cos}\left(2\pi j\frac{kk_0}{k_fk_0}\right)\hfill & & k_0<k<k_f\hfill \\ 0\hfill & & \text{otherwise}\hfill \end{array}$$
(58)
with $`(a_0,a_1,a_2)=(0.42323,0.49755,0.07922)`$. We choose $`k_0`$ and $`k_f`$ so that the window function goes smoothly to zero at the first and last eigenvalues of the three-particle spectrum. Numerical integration of (55) and (56) using this integral operator is displayed in Fig. 6. In the semiclassical transform, a total of 212 periodic orbits (including multiple repetitions) were used.

We observe good agreement between quantum and semiclassical results for $`L<7`$. In fact, it is difficult to distinguish between the two curves. For this reason, we plot the difference between them in Fig. 7. Clearly, the errors are small with respect to individual peak heights. Furthermore, the errors are largely due to hetero-orbit contributions. This can be understood by considering the first three structures in $`L`$-space. The first structure ($`L2.60`$) is due to a type-1 hetero-orbit where two particles are stationary and one evolves on $`\gamma =\frac{1}{2}(^{}`$2a). The second structure ($`L3.67`$) is from a type-2 hetero-orbit where one particle is stationary and two particles evolve independently on the same orbit $`\gamma =\frac{1}{2}(^{}`$2a). The third structure arises from the interference between a type-1 hetero-orbit ($`L4.62`$) where one of the three particles is on $`\gamma =\frac{1}{2}(^{}`$4b) and a dynamical orbit ($`L4.50`$) where all three particles evolve independently on $`\gamma =\frac{1}{2}(^{}`$2a). We see that the first and third structures have similar errors and thus conclude that the error introduced from the dynamical term is much smaller than that from the heterogeneous terms. All other $`L`$-space structures arise from the interference of many orbits and can be accounted for in a similar manner. For $`L>7`$, the discrepancies are more significant and mostly due to the problematic orbits $`\gamma =\text{4a}`$ and $`\gamma =\frac{1}{2}(^{}10\text{b})`$ that are not well isolated in phase space and pass close to the cusp of the cardioid (cf. Table 1). These orbits have inaccurate Gutzwiller amplitudes for reasons explained in Refs. .

### V.2 Symmetry decomposition

Due to the identical nature of the particles, the eigenstates of $`\widehat{H}`$ can be classified according to the irreps of $`S_3`$, the permutation group of three identical particles. Each group element belongs to one of three classes ($`(3,0,0)`$,$`(1,1,0)`$,$`(0,0,1)`$) based on the cycle structure of that element. Thus, there are also three irreps. These are the symmetric (trivial) irrep $`A_+`$, the antisymmetric irrep $`A_{}`$ and the two-dimensional mixed-symmetry irrep $``$. ($`S_N`$ always possesses exactly two one-dimensional irreps regardless of the size of $`N>1`$.) The character table for $`S_3`$ is given below. Numbers in front of class labels indicate the number of elements in that class.

#### V.2.1 Quantum mechanics

The total three-particle density of states can be decomposed into symmetry-reduced densities $`\rho _{}(E)`$, each belonging to an irrep $``$ of $`S_3`$:

$$\rho _3(E)=\rho _+(E)+\rho _{}(E)+\rho _{}(E).$$
(59)
Each partial density may be obtained by projection $`\rho _{}(E)=\text{Tr}(\widehat{P}_{}\delta (E\widehat{H}))`$, where the operator $`\widehat{P}_{}`$ projects onto the irrep $``$ <sup>2</sup><sup>2</sup>2For a discete group $`G`$, the operator $`\widehat{P}_{}`$ that projects onto the irrep $``$ is $`\widehat{P}_{}=(d_{}/|G|)_g\chi _{_{}}(g)\widehat{U}_g^{}`$, where the sum is over all group elements $`gG`$, $`d_{}`$ is the dimension of the irrep, $`|G|`$ is the order of the group, $`\chi _{_{}}(g)`$ is the character of the group element $`g`$ in the irrep $``$ and $`\widehat{U}_g`$ is the operator that transforms $`\mathrm{\Psi }`$ as prescribed by the group element $`gG`$. Permutation operators are unitary ($`S_N`$ is a unitary group). For $`A_\pm `$ irreps, this reduces to the operator of Eq. (41) and for the irrep $``$ of $`S_3`$, using the information provided in Table 2, $`\widehat{P}_{}=\frac{2}{6}(2\widehat{I}\widehat{U}_{\tau _1}\widehat{U}_{\tau _2})`$, where $`\tau _{1/2}(0,0,1)`$ and $`\widehat{I}`$ denotes the identity operator.. Expressing the trace in the energy eigenbasis as in (9), the symmetry-reduced densities are

$$\rho _\pm (E)=\frac{1}{6}\left[\rho _3(E)\pm \frac{3}{2}\rho _1\left(\frac{E}{2}\right)\rho _1(E)+\frac{2}{3}\rho _1\left(\frac{E}{3}\right)\right],$$
(60)
$`\rho _{}(E)={\displaystyle \frac{2}{6}}\left[2\rho _3(E){\displaystyle \frac{2}{3}}\rho _1\left({\displaystyle \frac{E}{3}}\right)\right].`$ (61)
To understand how the cross term arises in Eq. (60), consider the contribution from $`\tau =(a)(bc)(1,1,0)`$:

$`{\displaystyle \underset{i,j,k}{}}ijk|\widehat{U}_\tau \delta (E\widehat{H})|ijk`$ $`=`$ $`{\displaystyle \underset{i,j,k}{}}ikj|ijk\delta (EE_{ijk})`$ (62)
$`=`$ $`{\displaystyle \underset{i,j}{}}\delta (EE_{ijj})`$
$`=`$ $`{\displaystyle \underset{i,j}{}}\delta (E(\epsilon _i+2\epsilon _j))`$
$`=`$ $`{\displaystyle \frac{1}{2}}\rho _1\left({\displaystyle \frac{E}{2}}\right)\rho _1(E).`$
The other contributions can be found in a similar manner. We could compute each partial density separately for comparison with the numerics, but it is more illuminating to isolate the contribution from each symmetry class by inverting the above system of equations: (we ignore the identity class which reproduces the total density of states)

$`\rho _{(1,1,0)}(E)`$ $``$ $`\rho _+(E)\rho _{}(E)`$ (63)
$`=`$ $`{\displaystyle \frac{1}{2}}\rho _1\left({\displaystyle \frac{E}{2}}\right)\rho _1(E),`$
$`\rho _{(0,0,1)}(E)`$ $``$ $`\rho _+(E)+\rho _{}(E){\displaystyle \frac{1}{2}}\rho _{}(E)`$ (64)
$`=`$ $`{\displaystyle \frac{1}{3}}\rho _1\left({\displaystyle \frac{E}{3}}\right),`$
where $`\rho _{(1,1,0)}(E)`$ and $`\rho _{(0,0,1)}(E)`$ denote the densities belonging to the class of two-particle and three-particle exchanges, respectively. From Eq. (64), we note that the contribution of the longest cycle is directly related to the one-particle density of states as discussed at the end of section IV.1. In the following sections, we discuss the semiclassical decomposition of each partial density into smooth and oscillatory components:

$`\rho _{}^{\text{sc}}(E)=\overline{\rho }_{}(E)+\stackrel{~}{\rho }_{}(E).`$ (65)
#### V.2.2 Stationary cycles

If all particles being permuted are fixed, the cycles are stationary and contribute to $`\overline{\rho }_{}(E)`$. Using the results of Appendix D, it can be shown that the smooth densities belonging to each irrep are given by Eqs. (60-61), but with $`\rho `$ replaced by $`\overline{\rho }`$. Thus,

$`\overline{\rho }_{(1,1,0)}(E)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\overline{\rho }_1\left({\displaystyle \frac{E}{2}}\right)\overline{\rho }_1(E)`$ (66)
$`=`$ $`{\displaystyle \frac{1}{2}}[{\displaystyle \frac{\alpha ^2𝒜^2}{16\pi ^2}}E{\displaystyle \frac{\alpha ^{3/2}𝒜(1+\sqrt{2})}{16\pi ^2}}\sqrt{E}`$
$`+`$ $`{\displaystyle \frac{\alpha }{4\pi }}({\displaystyle \frac{\sqrt{2}^2}{16}}+3𝒜𝒦)],`$
$`\overline{\rho }_{(0,0,1)}(E)`$ $`=`$ $`{\displaystyle \frac{1}{3}}\overline{\rho }_1\left({\displaystyle \frac{E}{3}}\right)`$ (67)
$`=`$ $`{\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{\alpha 𝒜}{4\pi }}{\displaystyle \frac{\sqrt{3\alpha }}{8\pi }}{\displaystyle \frac{}{\sqrt{E}}}+3𝒦\delta (E)\right].`$
Note that in Eq. (66), we have ignored terms that are $`O(1/\mathrm{})`$ since some of the contributions at this order cannot be calculated exactly. These terms can be computed numerically, but are insignificant for our analysis.

#### V.2.3 Hetero-cycles

The leading-order cycles are the three 1-cycles of the identity class. If one 1-cycle is stationary and the other two 1-cycles are dynamical ($`s=1,e=2`$), then the result from Eq. (44) is identical to $`\stackrel{~}{\rho }_\mathrm{\Gamma }^\text{h}(2,3,E)`$. If instead two 1-cycles are stationary and one 1-cycle is dynamical ($`s=2,e=1`$), then Eq. (44) reduces to $`\stackrel{~}{\rho }_\mathrm{\Gamma }^\text{h}(1,3,E)`$. There are three contributions of each type.

The first correction is from permutations $`\tau (1,1,0)`$ that consist of one 1-cycle and one 2-cycle. There are two such contributions. The first one is from hetero-cycles for which the 1-cycle is stationary and the 2-cycle is dynamical. Using formula (44) such that $`k=1`$ is the 1-cycle and $`k=2`$ is the 2-cycle ($`n_1=1,n_2=2,J_2=2E_2=H=EE_2=E/2`$), the result has the structure of the leading-order term of $`\frac{1}{2}\stackrel{~}{\rho }_1(E/2)\overline{\rho }_1(E)`$. There is also the situation for which the 2-cycle is stationary and the 1-cycle is dynamical. Using Eq. (44) such that $`k=1`$ is the dynamical cycle (just a standard periodic orbit of the one-particle phase space) and $`k=2`$ is the 2-cycle, the result has the structure of the leading-order term of $`\frac{1}{2}\overline{\rho }_1(E/2)\stackrel{~}{\rho }_1(E)`$. Group elements $`\tau (0,0,1)`$ consist of single 3-cycles. Therefore, there are no contributions from this class. To summarise, we have shown that

$$\begin{array}{c}\stackrel{~}{\rho }_\pm ^\text{h}(E)=\frac{1}{6}[3\stackrel{~}{\rho }_3^{\text{h1}}(E)+3\stackrel{~}{\rho }_3^{\text{h2}}(E)\hfill \\ \hfill \pm 3(\frac{1}{2}\stackrel{~}{\rho }_1\left(\frac{E}{2}\right)\overline{\rho }_1(E)+\frac{1}{2}\overline{\rho }_1\left(\frac{E}{2}\right)\stackrel{~}{\rho }_1(E))],\end{array}$$
(68)
$`\stackrel{~}{\rho }_{}^\text{h}(E)=2\left[\stackrel{~}{\rho }_3^{\text{h1}}(E)+\stackrel{~}{\rho }_3^{\text{h2}}(E)\right].`$ (69)
#### V.2.4 Dynamical cycles

The leading-order contribution to $`\stackrel{~}{\rho }_\pm ^\text{d}(E)`$ comes from the identity element $`\iota =(a)(b)(c)`$ which consists of three 1-cycles ($`m_\tau =3;J_1=h_a,J_2=h_b,J_3=h_c`$; $`_kJ_k=H`$). Thus, there are 2 independent commuting generators other than $`H`$ and so we expect periodic orbit families of dimension 2. Using Eq. (42) and the fact that 1-cycles are equivalent to periodic orbits of the one-particle phase space, we find the leading-order term of $`\stackrel{~}{\rho }_\pm ^\text{d}(E)`$ is $`\stackrel{~}{\rho }_3^\text{d}(E)/6`$.

The next contribution is from permutations $`\tau (1,1,0)`$. There are three elements in this class each consisting of one 1-cycle and one 2-cycle ($`m_\tau =2`$; $`k=1,n_1=1`$; $`k=2,n_2=2`$). Then, for $`\tau =(ab)(c)`$, $`J_1=h_c,J_2=h_a+h_b`$ and similarly for the other elements in this class. Thus, there is only one independent generator (other than $`H`$) and we expect one-dimensional families. Using Eq. (42), we find this contribution has the structure of a two-particle density. The 1-cycles ($`k=1`$) are assigned to $`\gamma _1`$ and the 2-cycles ($`k=2`$) to $`\gamma _2`$, where $`\gamma _{1/2}`$ are any periodic orbits of the one-particle phase space. Then, all cycle properties are those of the corresponding orbit (cf. the 1-cycle and 2-cycle of Fig. 3; note the repetition $`l_2/n_2=2/2=1`$ which denotes the case where the particles of the 2-cycle evolve together is not shown). Multiple repetitions of the 2-cycle are either fractions (if $`l_2`$ is odd) and correspond to type-1 DPPOs *or* integers (if $`l_2`$ is even) and are type-0 DPPOs (cf. the classification used in section II.4.2). The generators $`J_1=n_1E_1=E_1`$ and $`J_2=n_2E_2=2E_2`$ are the energies of the particles involved in the 1-cycle and 2-cycle, respectively (particles of the 2-cycle *each* have energy $`E/2`$). Thus, the final form is structurally equivalent to $`\frac{1}{2}\stackrel{~}{\rho }_1(E/2)\stackrel{~}{\rho }_1(E)`$.

The two group elements $`\tau (0,0,1)`$ each consist of one 3-cycle ($`m_\tau =1,k=1,n_1=3`$) which implies there are no generators independent of $`H`$ and thus the orbits are isolated. As before, cycle properties can be mapped to those of an orbit of the one-particle phase space (cf. the 3-cycle shown Fig. 4; $`l_k/n_k=l_1/n_1=1/3,2/3,3/3`$; higher repetitions $`l_1/n_1=l_1/3`$ would have action $`l_1S`$, phase index $`l_1\sigma `$, stability matrix $`\stackrel{~}{M}^{l_1}`$, where $`S,\sigma ,\stackrel{~}{M}`$ are the properties of the primitive orbit to which the cycle is assigned). The energy $`E_1`$ in Eq. (42) is the energy of each particle involved in the 3-cycle $`(k=1)`$ and since $`H=J_1=3E_1=E`$, it follows that $`E_1=E/3`$. Thus, the result has the structure of a one-particle trace formula, but it is evaluated at $`E/3`$ and has a cycle structure prefactor of 1/3. Including the prefactors from the projection operator, we conclude that

$$\stackrel{~}{\rho }_\pm ^\text{d}(E)=\frac{1}{6}\left[\stackrel{~}{\rho }_3^\text{d}(E)\pm \frac{3}{2}\stackrel{~}{\rho }_1\left(\frac{E}{2}\right)\stackrel{~}{\rho }_1(E)+\frac{2}{3}\stackrel{~}{\rho }_1\left(\frac{E}{3}\right)\right],$$
(70)
$`\stackrel{~}{\rho }_{}^\text{d}(E)={\displaystyle \frac{2}{6}}\left[2\stackrel{~}{\rho }_3^\text{d}(E){\displaystyle \frac{2}{3}}\stackrel{~}{\rho }_1\left({\displaystyle \frac{E}{3}}\right)\right].`$ (71)
We stress that even though the correction terms have structures equivalent to one- and two-particle densities, they are in fact contributions from the dynamical cycles of the full three-particle phase space.

#### V.2.5 Trace formulas for the two symmetry classes

Combining the results of Eqs. (68-71), the fluctuating densities for the two nontrivial symmetry classes are

$`\stackrel{~}{\rho }_{(1,1,0)}(E)`$ $`=`$ $`\stackrel{~}{\rho }_+(E)\stackrel{~}{\rho }_{}(E)={\displaystyle \frac{1}{2}}\stackrel{~}{\rho }_1\left({\displaystyle \frac{E}{2}}\right)\stackrel{~}{\rho }_1(E)`$ (72)
$`+`$ $`\left[{\displaystyle \frac{1}{2}}\overline{\rho }_1\left({\displaystyle \frac{E}{2}}\right)\stackrel{~}{\rho }_1(E)+{\displaystyle \frac{1}{2}}\stackrel{~}{\rho }_1\left({\displaystyle \frac{E}{2}}\right)\overline{\rho }_1(E)\right]`$
$`=`$ $`\stackrel{~}{\rho }_{(1,1,0)}^\text{d}(E)+\left[\stackrel{~}{\rho }_{(1,1,0)}^{\text{h1}}(E)+\stackrel{~}{\rho }_{(1,1,0)}^{\text{h2}}(E)\right]`$
$`=`$ $`\stackrel{~}{\rho }_{(1,1,0)}^\text{d}(E)+\stackrel{~}{\rho }_{(1,1,0)}^\text{h}(E)`$
and

$`\stackrel{~}{\rho }_{(0,0,1)}(E)`$ $`=`$ $`\stackrel{~}{\rho }_+(E)+\stackrel{~}{\rho }_{}(E){\displaystyle \frac{1}{2}}\stackrel{~}{\rho }_{}(E)`$ (73)
$`=`$ $`{\displaystyle \frac{1}{3}}\stackrel{~}{\rho }_1\left({\displaystyle \frac{E}{3}}\right).`$
The leading-order term of $`\stackrel{~}{\rho }_{(1,1,0)}^\text{h}(E)`$ is given by

$$\begin{array}{cc}& \stackrel{~}{\rho }_{(1,1,0)}^\text{h}(E)=\frac{\alpha 𝒜}{4\pi ^2}\underset{\gamma }{}\frac{(L_\gamma ^0/2L_\gamma )}{\sqrt{\left|det\left(\stackrel{~}{M}_\gamma I\right)\right|}}\hfill \\ & \times \mathrm{cos}\left(\sqrt{\alpha E}L_\gamma \sigma _\gamma \frac{\pi }{2}\frac{\pi }{2}\right)\hfill \\ & +\left\{L_\gamma ^0\sqrt{2}L_\gamma ^0,L_\gamma \frac{L_\gamma }{\sqrt{2}}\right\}.\hfill \end{array}$$
(74)
The first term of (74) is the contribution from two particles being stationary at the same point in the billiard (i.e. a stationary 2-cycle) while the third particle evolves on a periodic orbit (i.e. a dynamical 1-cycle). The second term is the contribution from one particle being stationary (i.e. a stationary 1-cycle) while the other two particles evolve on a periodic orbit (i.e. a dynamical 2-cycle). Higher-order contributions from hetero-cycles can be worked out. These are included in the numerics, but we do not write them out explicitly here (see Appendix E). The contribution from the dynamical cycles as determined above can be written as

$$\begin{array}{cc}& \stackrel{~}{\rho }_{(1,1,0)}^\text{d}(E)=\frac{\alpha ^{3/4}}{(2\pi )^{3/2}E^{1/4}}\underset{\gamma _1,\gamma _2}{}\left(\underset{i=1}{\overset{2}{}}\frac{L_{\gamma _i}^0}{\sqrt{\left|det\left(\stackrel{~}{M}_{\gamma _i}I\right)\right|}}\right)\hfill \\ & \times \left[2(2L_{\gamma _1}^2+L_{\gamma _2}^2)\right]^{1/4}\mathrm{cos}\left[\frac{\sqrt{\alpha E}}{\sqrt{2}}L_{12}\sigma _{12}\frac{\pi }{2}\frac{\pi }{4}\right],\hfill \end{array}$$
(75)
where $`L_{12}=\sqrt{2L_{\gamma _1}^2+L_{\gamma _2}^2}`$ and $`\sigma _{12}=(\sigma _{\gamma _1}+\sigma _{\gamma _2})`$. To understand how this result is obtained, recall the structure of the dynamical cycles in this class. Each full cycle consists of one 1-cycle and one 2-cycle. The total energy is partitioned among the three particles such that the periods of the cycles are the same. Suppose the 1-cycle and 2-cycle are associated with the orbits $`\gamma _1`$ and $`\gamma _2`$, respectively. The energies $`E_1,E_2`$ are determined from the periodicity condition

$`T_{\gamma _1}(E_1)={\displaystyle \frac{1}{2}}T_{\gamma _2}\left(E_2\right)T.`$ (76)
Using the usual relations for actions and periods in a billiard (V), one can show that

$$E_1=\left[\frac{2L_{\gamma _1}^2}{2L_{\gamma _1}^2+L_{\gamma _2}^2}\right]E,E_2=\left[\frac{(L_{\gamma _2}^2/2)}{2L_{\gamma _1}^2+L_{\gamma _2}^2}\right]E,$$
(77)
where $`E_1`$ is the energy of the particle of the 1-cycle, $`2E_2`$ is the total energy of the particles involved in the 2-cycle (each of them has energy $`E_2`$ since their energies must be equal) and $`E=E_1+2E_2`$ is the total energy of the three-particle system. The 2-cycles are similar to the DPPOs of a two-particle system (cf. section II.4.2) and we shall use the same classification scheme for all 2-cycles. The trace formula for $`\stackrel{~}{\rho }_{(0,0,1)}(E)`$ is a one-particle trace formula except that $`L_\gamma ^0\sqrt{3}L_\gamma ^0`$ and $`L_\gamma \frac{L_\gamma }{\sqrt{3}}`$.

#### V.2.6 Numerics

We first consider the class $`(1,1,0)`$. We compare numerically the length spectrum of the dynamical cycles

$$\stackrel{~}{F}_{(1,1,0)}^{\text{sc}}(L)=\{\stackrel{~}{\rho }_{(1,1,0)}^\text{d}(k)\}$$
(78)
with its quantum analogue

$$\stackrel{~}{F}_{(1,1,0)}^{\text{qm}}(L)=\{\rho _{(1,1,0)}(k)\overline{\rho }_{(1,1,0)}(k)\stackrel{~}{\rho }_{(1,1,0)}^\text{h}(k)\}.$$
(79)
We construct the first 241,080 levels of $`\rho _{(1,1,0)}(k)`$ using the first 1000 single-particle energies. The smooth term is computed from Eq. (66) using the billiard parameters (as above) for the odd spectrum. Trace formulas are computed using geometrical orbits with length $`L<10`$. The result is shown in Fig. 8.

We now examine some of the $`L`$-space structures of Fig. 8. The first peak ($`L3.18`$) is due to the dynamical cycle where both 1- and 2-cycles are on the primitive orbit $`\gamma ^0=\frac{1}{2}(^{}`$2a). The particle of the 1-cycle completes one full motion on the orbit while the particles of the 2-cycle each traverse half the orbit and are then exchanged. The second peak ($`L4.17`$) occurs because of a dynamical cycle where the 1- and 2-cycles are on primitve orbits $`\gamma _1^0=\frac{1}{2}(^{}`$2a) and $`\gamma _2^0=\frac{1}{2}(^{}`$4b), respectively. For the third peak ($`L4.50`$), the 1-cycle is as in the first case, except the 2-cycle is type-$`0`$. Thus, as the particle of the 1-cycle completes one full motion on $`\gamma ^0=\frac{1}{2}(^{}`$2a), the particles of the 2-cycle each traverse the full orbit and are then exchanged. This is summarised in Table 3 where some of the dynamical cycles in this class are listed. In each case, the energies are divided according to Eq. (77). For $`L<7`$, there are a total of 27 dynamical cycles. All structures in $`L`$-space can be accounted for in a similar manner and can be checked systematically by noting that due to the energy division between the particles, we expect peaks at positions $`L=\sqrt{L_{\gamma _1}^2+(n_{\gamma _2}^\pm L_{\gamma _2}^0)^2/2}`$. The 2-cycles are type-0 and type-1 for even ($`n_{\gamma _2}^+`$) and odd ($`n_{\gamma _2}^{}`$) integer repetition indices, respectively.

The discrepancies between quantum and semiclassical results are due to the problematic orbits mentioned above. The structure at $`L5`$ is poorly reproduced due to the hetero-cycles involving 1-cycles that are stationary and 2-cycles that are dynamical and evolving on the problematic orbits $`\gamma ^0=4\text{a}`$ and $`\gamma ^0=\frac{1}{2}(^{}`$10b). (The length spectrum of the hetero-cycles is shown in Fig. 9.) All other discrepancies are due to purely dynamical cycles and these are summarised in Table 4.

We next consider the class $`(0,0,1)`$. We numerically compare the Fourier transform of the dynamical 3-cycles

$`\stackrel{~}{F}_{(0,0,1)}^{\text{sc}}(L)=\{\stackrel{~}{\rho }_{(0,0,1)}^\text{d}(k)\}`$ (80)
with its quantum analogue

$`\stackrel{~}{F}_{(0,0,1)}^{\text{qm}}(L)=\{\rho _{(0,0,1)}(k)\overline{\rho }_{(0,0,1)}(k)\}.`$ (81)
We use the first 1000 energies of $`\rho _{(0,0,1)}(k)`$. The smooth term is computed from Eq. (67) and trace formulas are computed using geometrical orbits with length $`L<11`$. The result is shown in Fig. 10. We now identify some of the peak structures with one or several of the orbits listed in Fig. 5.

The first peak ($`L1.5`$) can be identified with a type-$`1`$ dynamical 3-cycle consisting of all three particles evolving on the orbit $`\gamma ^0=\frac{1}{2}(^{}`$2a) with the same energy and exactly $`T_\gamma /3`$ out of phase and each completing one-third of the full motion on the orbit and finally being permuted as specified by $`\tau =(acb)`$. The third peak ($`L3`$) is due to a type-$`2`$ dynamical 3-cycle where all three particles evolve on the orbit $`\gamma ^0=\frac{1}{2}(^{}`$2a) as above, except that each particle completes two-thirds of the full motion on the orbit before being exchanged according to $`\tau =(abc)`$. The peak at $`L4.5`$ is from a type-$`0`$ dynamical 3-cycle consisting of all three particles starting and evolving together in phase on $`\gamma ^0=\frac{1}{2}(^{}`$2a), but each completing one full motion on the orbit and then being trivially exchanged as presribed by any group element $`\tau (0,0,1)`$. As an example of a higher multiple cycle, consider the first repetition of the type-$`1`$ cycle mentioned above. It is the same as before except that each particle completes one and one-third of the motion on the orbit before being permuted. This is summarised in Table 5 where some of the dynamical 3-cycles are listed. (The peak at $`L4.25`$ is completely undetected by the trace formula since it arises from a diffractive orbit. Such orbits require a separate analysis since they are not included in the standard Gutzwiller theory ). As before, the discrepancy that occurs at $`L5`$ arises from the two orbits $`\gamma ^0=4\text{a}`$ and $`\gamma ^0=\frac{1}{2}(^{}`$10b).

All structures in $`L`$-space can be accounted for in a similar manner. As a systematic check, recall that each dynamical 3-cycle can be mapped one-to-one with a periodic orbit $`\gamma `$ of the one-particle phase space. If the orbit has length $`L_\gamma =n_\gamma L_\gamma ^0`$, where $`n_\gamma `$ is a repetition index, it is mapped to a 3-cycle where each particle executes a fraction $`n_\gamma /3`$ of the full motion on $`\gamma `$. We can then write $`n_\gamma /3=i+j/3`$, where $`i=\text{int}(n_\gamma /3)`$ (i.e. integer part of $`n_\gamma /3`$) and $`j/3`$ ($`j=0,1,2`$) is the remainder. If $`j0`$, then the orbit with length $`L_\gamma `$ is associated with the $`i`$th repetition of a type-$`j`$ dynamical 3-cycle. If $`j=0`$, then it is the $`(i1)`$th repetition. To determine peak positions, we recall that all particles of the 3-cycle have the same energy $`E/3`$ and thus we expect peaks at lengths $`L=n_\gamma L_\gamma ^0/\sqrt{3}`$. (Recall that a billiard orbit with length $`L_\gamma `$ has action $`S_\gamma (\epsilon )/\mathrm{}=\sqrt{\alpha \epsilon }L_\gamma `$.)

## VI Conclusion

We began with the case of two noninteracting identical particles evolving dynamically on periodic orbits and explained how the time-translational symmetry leads to families of periodic orbits in the full phase space. Using the trace formula for continuous symmetries , we obtained a trace formula for the two-particle resolvent consistent with the dynamical term of the semiclassical two-particle density of states . We also proved identities for the symmetry-reduced densities of states (see also Appendix A.1) which were stated without proof in our previous work . Dynamical pseudo-periodic orbits (DPPOs) were defined and it was shown that their contribution to the semiclassical exchange term has the form of a one-particle trace formula. We also introduced two-particle heterogeneous periodic orbits in the full phase space. We discussed how the structure of these orbits is different in billiards and in analytic potentials and the explicit contribution of such orbits to the two-particle resolvent was determined.

We have also demonstrated that the approach used in this paper yields results that are consistent with those of the convolution method . In the convolution picture, one is faced with the asymptotic analysis of many convolution integrals and the further issue of spurious contributions from them. The full phase space formalism, on the other hand, more easily generalises to arbitrary particle numbers. It is also more illuminating since it reveals the underlying structure of the periodic orbit families. (One useful property of the convolution approach occurs for billiards where significant higher-order corrections from hetero-orbits can be explicitly calculated (cf. Appendix E)). Most importantly, the convolution formalism does not accomodate interactions and it is necessary to use the full phase space if interactions are to be included. The symmetric/antisymmetric resolvent for the case of noninteracting particles can be expressed as a sum of resolvents, one for each element of the permutation group. As shown above, trace formulas for the oscillatory components can be written as products over cycles, where each cycle is assigned to a periodic orbit of the one-particle phase space. The results of the paper were applied to a specific problem: three noninteracting identical particles in a two-dimensional cardioid billiard. We found that our semiclassical analysis correctly reproduced the quantum results and we explained how these results could be understood in terms of classical structures in the full phase space.

We have assumed that the single-particle dynamics is free of any continuous symmetry. If there are additional symmetries, then they must also be properly accounted for in the theory. If we restrict ourselves to the noninteracting case, then essentially the only difference from what we have presented above is that $`𝐉`$ and $`𝚯`$ become higher dimensional. (An example would be two particles in a disk billiard. In this case, there are *four* independent constants of the motion, any two of $`\{E,E_a,E_b\}`$ and any two of $`\{L_z,L_{z_a},L_{z_b}\}`$. Thus, the periodic orbits occur in 3-parameter families.) One future goal is to consider separately the important zeroth-order problem of harmonic oscillator potentials. The harmonic oscillator has a higher degree of symmetry than we are accounting for here ($`SU(d)`$ in $`d`$ dimensions). This project would require using the theory of Ref.  which derives a trace formula for systems with more general symmetries including non-Abelian cases.

As mentioned above, a major advantage of the formalism presented here (as compared to the convolution formalism of Ref. ) is that it can be extended to include interactions. Any interaction between the particles destroys the periodic orbit families described above and replaces them with a discrete set of isolated orbits. For weak interactions, we can use a perturbative method  which is applicable to any situation where continuous symmetries are broken. This program will be explored in a future publication . One could also study (zero-range) point interactions. Such interactions are often considered as corrections to mean-field approximations. Semiclassically, point-interactions can be understood using the formalism of diffractive orbits . One can imagine that such an interaction leaves the periodic orbit families (of the noninteracting system) largely unchanged, but introduces qualitatively new diffractive orbits. We plan to explore this scenario in future work.

###### Acknowledgements.

We thank Rajat Bhaduri, Stephen Creagh, Randy Dumont and David Goodings for useful discussions. This work was financially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).

## Appendix A Two-particle Thomas-Fermi calculation

We first discuss the smooth two-particle density of states and its decomposition into bosonic and fermionic densities. Using the identity $`\delta (Eh(z_a)h(z_b))=\text{d}\epsilon \delta (\epsilon h(z_a))\delta (E\epsilon h(z_b))`$, we can show that the leading-order smooth term for the two-particle density of states is the autoconvolution of the leading-order smooth term of the one-particle density of states (12). We could verify this term-by-term in the expansion of $`\overline{\rho }_2(E)`$, but we can do it more efficiently for all terms as follows.

We use the partition function $`Z(\beta )=\text{Tr}(\mathrm{exp}(\beta \widehat{H}))`$ which is the Laplace transform of the density of states. It is convenient to work with the Wigner transform which is defined for an arbitrary operator $`\widehat{A}`$ as

$$\widehat{A}_\text{W}(z)=\text{d}xq+\frac{x}{2}\left|\widehat{A}\right|q\frac{x}{2}\mathrm{exp}\left(i\frac{px}{\mathrm{}}\right),$$
(82)
in terms of which the trace is

$$\text{Tr}\{\widehat{A}\}=\frac{1}{(2\pi \mathrm{})^n}\text{d}z\widehat{A}_\text{W}(z).$$
(83)
The trace of a product of two (but not more) operators is given by

$$\text{Tr}\{\widehat{A}\widehat{B}\}=\frac{1}{(2\pi \mathrm{})^n}\text{d}z\widehat{A}_\text{W}(z)\widehat{B}_\text{W}(z).$$
(84)
The Wigner transform of the evolution operator, $`\mathrm{exp}(\beta \widehat{H})_\text{W}(z)`$, can be written as an asymptotic expansion in powers of $`\mathrm{}`$, the first few terms of which are typically retained and used as the smooth approximation to the partition function. Taking the inverse Laplace transform then gives the smooth density of states. In particular, the leading-order term of $`\mathrm{exp}(\beta \widehat{H})_\text{W}(z)`$ is $`\mathrm{exp}(\beta \widehat{H}_\text{W}(z))`$, where the Wigner transform of the quantum Hamiltonian $`\widehat{H}_\text{W}(z)`$ is simply the classical Hamiltonian which we have denoted by $`H(z)`$. (There are corrections to this if the Hamiltonian is not of the kinetic plus potential form.) The inverse Laplace transform of this expression yields the leading-order smooth term (12).

For two independent particles, the full quantum Hamiltonian is the sum of one-particle Hamiltonians and since these are functions of independent phase space variables,

$$\mathrm{exp}(\beta \widehat{H})_\text{W}(z)=\mathrm{exp}(\beta \widehat{h})_\text{W}(z_a)\mathrm{exp}(\beta \widehat{h})_\text{W}(z_b).$$
(85)
Thus, the smoothed two-particle partition function is simply the product of smooth one-particle partition functions. By the Laplace convolution theorem, this implies that the smoothed two-particle density of states is the autoconvolution of the smoothed one-particle density of states. This same argument can be made for the exact density of states as an alternate proof of Eq. (5).

### A.1 Symmetrised $`2`$-particle Thomas-Fermi term

The bosonic and fermionic partition functions are $`Z_\pm (\beta )=\text{Tr}(\widehat{P}_\pm \mathrm{exp}(\beta \widehat{H}))`$, where $`\widehat{P}_\pm `$ are the projection operators defined in Eq. (7). The leading-order term is just the two-particle partition funtion. The next term $`\text{Tr}(\widehat{U}\mathrm{exp}(\beta \widehat{H}))`$ requires slightly more analysis and can be evaluated directly using Eq. (84). We begin by finding $`\widehat{U}_\text{W}(z)`$.

It is shown in Ref.  that for a one-particle system with a symmetry axis through the coordinate $`q`$, the Wigner transform of the reflection operator is

$$\left(\widehat{}\right)_\text{W}(z)=\pi \mathrm{}\delta (q)\delta (p),$$
(86)
where $`p`$ is the momentum conjugate to $`q`$. We map our problem onto that one as follows. First, suppose that the one-particle system is one-dimensional and define the Jacobi coordinates: $`q=q_aq_b`$, $`p=(p_ap_b)/2`$, $`Q=(q_a+q_b)/2`$ and $`P=p_a+p_b`$. Then, exchanging $`a`$ and $`b`$ is equivalent to reflecting in $`q`$ so that the variable $`q`$ in the above equation is replaced by $`q_aq_b`$ and $`p`$ is replaced by the conjugate momentum $`(p_ap_b)/2`$. Then, $`\widehat{U}_\text{W}(z)=2\pi \mathrm{}\delta (q_aq_b)\delta (p_ap_b)`$. If the one-particle system is higher dimensional, then $`\widehat{U}`$ is the product of one such inversion in every component. All of them are independent so that the final result is the product of the individual ones (for the same reason that Eq. (85) is multiplicative). The final result is

$$\widehat{U}_\text{W}(z)=(2\pi \mathrm{})^d\delta (z_az_b),$$
(87)
where the delta function represents the product of all $`2d`$ delta functions (two for each component). Equivalent results can be found in Ref. .

Then,

$$\begin{array}{cc}& \text{Tr}\left(\widehat{U}\mathrm{exp}(\beta \widehat{H})\right)\hfill \\ & =\frac{1}{(2\pi \mathrm{})^{2d}}\text{d}z\widehat{U}_\text{W}(z)\mathrm{exp}(\beta \widehat{H})_\text{W}(z)\hfill \\ & =\frac{1}{(2\pi \mathrm{})^d}\text{d}z_a\mathrm{exp}(2\beta \widehat{h})_\text{W}(z_a),\hfill \end{array}$$
(88)
where we used the delta functions from $`\widehat{U}_\text{W}(z)`$ in Eq. (87) to do the integrals over the $`z_b`$ variables and the multiplicative property of the Wigner functions as in Eq. (85). The first few terms of Eq. (88) give the smooth approximation to the one-particle partition function evaluated at $`2\beta `$. Under the inverse Laplace transform, this becomes $`\overline{\rho }_1(E/2)/2`$ and we conclude

$$\overline{\rho }_\pm (E)=\frac{1}{2}\left(\overline{\rho }_2(E)\pm \frac{1}{2}\overline{\rho }_1\left(\frac{E}{2}\right)\right),$$
(89)
consistent with (10).

## Appendix B Stability matrix of pseudo-periodic orbits

We first prove that $`det(\stackrel{~}{M}_\gamma ^{}I)=4det(\stackrel{~}{M}_\gamma I)`$, where $`\stackrel{~}{M}_\gamma ^{}`$ is the stability matrix of a pseudo-periodic orbit $`\gamma ^{}`$ of the full two-particle phase space and $`\stackrel{~}{M}_\gamma `$ is the stability matrix of the corresponding periodic orbit $`\gamma `$ in the one-particle phase space. This admits various generalisations which are used in the main discussion. The (type-1) dynamical pseudo-periodic orbit (DPPO) consists of both particles evolving for half of the single-particle period $`T_\gamma /2`$ followed by the symplectic mapping $`u`$ that exchanges the two particles.

We define coordinates as follows (cf. Fig. 11). For particle $`a`$, we define an initial section $`\mathrm{\Sigma }_a`$ such that the phase space flow is transverse to it and all points on the section are at equal energy. We define a coordinate pointing along the orbit which we call $`\eta _a`$. Without loss of generality, we can take $`\eta _a/t1`$. We also define a coordinate transverse to the constant $`h_a`$ surface (but in the phase space of particle $`a`$) which we call $`\kappa _a`$. If we consider it to take the values of $`h_a`$, then it is canonically conjugate to $`\eta _a`$ and has zero time derivative under the flow since $`h_a`$ is conserved. The remaining $`(2d2)`$ coordinates for particle $`a`$ lie on the section $`\mathrm{\Sigma }_a`$ and will be collectively denoted by $`\xi _a`$. As the flow evolves, changes in the $`\xi _a`$ coordinates are described by the $`(2d2)\times (2d2)`$ symplectic stability matrix (for the one-particle dynamics) $`\stackrel{~}{N}_a`$. Similarly, we define $`\mathrm{\Sigma }_b`$, $`\eta _b`$, $`\kappa _b`$, $`\xi _b`$ and $`\stackrel{~}{N}_b`$ for particle $`b`$. We also need a way to connect coordinates on $`\mathrm{\Sigma }_a`$ to those on $`\mathrm{\Sigma }_b`$; we will take them to be such that they are connected by parallel transport so that, for example, the stable and unstable manifolds are mapped onto each other.

We start by defining the symplectic transformation

$`\eta ={\displaystyle \frac{\eta _a+\eta _b}{2}},\kappa =\kappa _a+\kappa _b,`$
$`\upsilon =\eta _a\eta _b,\zeta ={\displaystyle \frac{\kappa _a\kappa _b}{2}}.`$ (90)
The monodromy matrix $`M_\gamma ^{}`$ describes the linearised motion of small perturbations around a DPPO $`\gamma ^{}`$ of the full phase space. In particular, if $`𝚼=(\kappa ,\eta ,\upsilon ,\zeta ,\xi _a,\xi _b)`$, then $`\delta 𝚼=M_\gamma ^{}\delta 𝚼_0`$. Consider an initial slight change in $`\eta `$ by the amount $`\delta \eta _0`$ while keeping all other coordinates constant. This implies that both $`\eta _a`$ and $`\eta _b`$ increase by $`\delta \eta _0`$. After time evolution for $`T_\gamma /2`$ and particle exchange, $`\delta \eta =\delta \eta _0`$ while all other coordinates are unchanged (in particular, the transverse coordinates are unaffected). Now consider an initial small change in $`\kappa `$ by the amount $`\delta \kappa _0`$. This implies that $`\kappa _a`$ and $`\kappa _b`$ change by $`\delta \kappa _0/2`$. After integrating for time $`T_\gamma /2`$ and interchanging the particles, we observe that $`\delta \kappa =\delta \kappa _0`$. However, this change of value in $`\kappa `$ does affect the $`\eta `$ coordinate. Under this change, the period of the orbit $`\gamma `$ also changes; let $`T_\gamma ^{}`$ denote the derivative of this period with respect to the single-particle energy. Since we are only integrating for half of the period, and the single-particle energies are changed by $`\delta \kappa _0/2`$, we find that $`\delta \eta =T_\gamma ^{}\delta \kappa _0/4`$. (The minus sign indicates that if the period increases and we integrate for the same amount of time as before, then the particles will fail to execute a complete loop, corresponding to a negative value of $`\eta `$.) Thus, the monodromy matrix of the pseudo-periodic orbit $`\gamma ^{}`$ in the full phase space has the form

$$M_\gamma ^{}=\left(\begin{array}{ccc}\begin{array}{cc}1& 0\\ \frac{T_\gamma ^{}}{4}& 1\end{array}& \mathrm{𝟎}& \\ \mathrm{𝟎}& \stackrel{~}{M}_\gamma ^{}& \end{array}\right).$$
(91)
We are interested in calculating $`det(\stackrel{~}{M}_\gamma ^{}I)`$ to evaluate the Gutzwiller amplitude. Note that the matrix $`\stackrel{~}{M}_\gamma ^{}`$ involves only the $`(4d2)`$ phase space coordinates other than $`\eta `$ and $`\kappa `$. We can understand the calculation up to now as follows. The transformation (B) can be thought of as a transformation to center-of-mass coordinates. We have removed the center-of-mass coordinates $`\eta `$ and $`\kappa `$ from consideration and are only left with the relative coordinates $`\upsilon `$ and $`\zeta `$ (as well as all the transverse coordinates $`\xi _a`$ and $`\xi _b`$.) It is reasonable that only the relative coordinates are important for determining the stability.

The next two coordinates we consider are $`\upsilon `$ and $`\zeta `$. Let us start with $`\upsilon `$. A small initial change in $`\upsilon `$ by the amount $`\delta \upsilon _0`$ implies that $`\eta _a`$ changes by $`\delta \upsilon _0/2`$ while $`\eta _b`$ changes by $`\delta \upsilon _0/2`$. After integrating for time $`T_\gamma /2`$, this remains unchanged, but after particle exchange, the final values of $`\delta \eta _a`$ and $`\delta \eta _b`$ are changed in sign so that the corresponding diagonal matrix element of $`\stackrel{~}{M}_\gamma ^{}`$ is $`1`$. Similarly, the diagonal matrix element corresponding to the $`\zeta `$ coordinate is also $`1`$. As before, an infinitesmal change in $`\zeta `$ implies an infinitesmal change in $`\upsilon `$. In this case, the corresponding matrix element is $`T_\gamma ^{}`$. Therefore, we can write

$$\stackrel{~}{M}_\gamma ^{}=\left(\begin{array}{ccc}\begin{array}{cc}1& 0\\ T_\gamma ^{}& 1\end{array}& \mathrm{𝟎}& \\ \mathrm{𝟎}& \stackrel{~}{N}& \end{array}\right).$$
(92)
Then, $`det(\stackrel{~}{M}_\gamma ^{}I)=4det(\stackrel{~}{N}I)`$, where we use appropriately dimensioned identity matrices on each side of the equality. It remains to calculate the determinant of the $`(4d4)\times (4d4)`$ matrix $`\stackrel{~}{N}I`$. The matrix $`\stackrel{~}{N}`$ involves only the coordinates $`\xi _a`$ and $`\xi _b`$. Since these two sets of coordinates live on different sections, we cannot immediately define a mapping between them. To do so, we note that we have defined coordinates on the two sections so that the exchange operation is a simple mapping of the form

$$\stackrel{~}{E}=\left(\begin{array}{cc}\mathrm{𝟎}& I\\ I& \mathrm{𝟎}\end{array}\right),$$
(93)
where $`I`$ is a $`(2d2)\times (2d2)`$ identity matrix. In terms of these coordinates, the one-particle stability matrices $`\stackrel{~}{N}_a`$ and $`\stackrel{~}{N}_b`$ are such that $`\stackrel{~}{N}_a\stackrel{~}{N}_b=\stackrel{~}{M}_\gamma `$, which is the stability matrix of the full periodic orbit for the one-particle dynamics. The combined operations of flow and exchange give

$$\stackrel{~}{N}=\left(\begin{array}{cc}\mathrm{𝟎}& \stackrel{~}{N}_b\\ \stackrel{~}{N}_a& \mathrm{𝟎}\end{array}\right)$$
(94)
and

$`det(\stackrel{~}{N}I)`$ $`=`$ $`det\left(\begin{array}{cc}I& \stackrel{~}{N}_b\\ \stackrel{~}{N}_a& I\end{array}\right)`$ (97)
$`=`$ $`det\left(\begin{array}{cc}\stackrel{~}{N}_a& I\\ I& \stackrel{~}{N}_b\end{array}\right)`$ (100)
$`=`$ $`det(\stackrel{~}{M}_\gamma I),`$ (101)
where after the second equality, we interchanged rows to put the matrix into a more useful form. (The matrix has even dimension so there is no sign introduced as a result of this interchange.) The final equality of Eq. (97) requires the following identity. If a matrix $`C`$ has the form

$$C=\left(\begin{array}{cc}A& I\\ I& B\end{array}\right),$$
(102)
then $`det(C)=det(ABI)`$. This can be shown by multiplying $`C`$ by the matrix

$$C^{}=\left(\begin{array}{cc}B& I\\ I& A\end{array}\right).$$
(103)
After multiplying them together, the product is block-diagonal with $`(ABI)`$ in one block and $`(BAI)`$ in the other. These have equal determinants. Since $`C^{}`$ has the same determinant as $`C`$, it follows that $`[det(C)]^2=[det(ABI)]^2`$ and thus we have identified the two determinants within a sign. The sign follows from the fact that the contribution to the determinant of the fully diagonal term $`A_{ii}B_{ii}`$ should be positive. Thus, we conclude that $`det(\stackrel{~}{M}_\gamma ^{}I)=4det(\stackrel{~}{M}_\gamma I)`$.

It is a straightforward extension to generalise this result to a cycle with $`n`$ particles on an orbit. We first have to find some appropriate set of variables so that we may isolate an upper-left block of the monodromy matrix in analogy to Eq. (92). This comes from the $`2n`$ coordinates $`\eta `$ and $`\kappa `$. As in the previous case, we separate the variables into center-of-mass coordinates (which subsequently play no role) and a set of relative (Jacobi) coordinates. Using similar arguments to the $`n=2`$ case, the determinant of the upper-left block is then $`n^2`$. The contribution from the rest of the matrix (i.e. the lower-right block) comes from the transverse coordinates. In terms of these transverse coordinates, the single-particle stability matrices $`\stackrel{~}{N}_a,\stackrel{~}{N}_b,\mathrm{},\stackrel{~}{N}_n`$ are such that $`\stackrel{~}{N}_a\stackrel{~}{N}_b\mathrm{}\stackrel{~}{N}_n=\stackrel{~}{M}_\gamma `$, which is the stability matrix of the full periodic orbit for the one-particle dynamics. Through a sequence of manipulations and transpositions similar to the $`n=2`$ case, the determinant of this lower-right block (i.e. $`det(\stackrel{~}{N}I)`$) can then be reduced to the form

$$det\left(\begin{array}{cccccc}\stackrel{~}{N}_a& I& 0& \mathrm{}& 0& 0\\ 0& \stackrel{~}{N}_b& I& \mathrm{}& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& \stackrel{~}{N}_{n1}& I\\ I& 0& 0& \mathrm{}& 0& \stackrel{~}{N}_n\end{array}\right),$$
(104)
which is a generalisation of the $`n=2`$ case. It can be shown that $`det(\stackrel{~}{N}I)=det(\stackrel{~}{M}_\gamma I)`$ and thus, $`det(\stackrel{~}{M}_\gamma ^{}I)=n^2det(\stackrel{~}{M}_\gamma I)`$. If the orbit is not primitive, but is a repetition of some simpler orbit, then we can absorb this into the definitions of the single-particle stability matrices and carry through all of the manipulations as before. The result is unchanged.

## Appendix C Monodromy matrix and phase indices of a harmonic oscillator

It is shown in Ref.  that the monodromy matrix for a primitive orbit along the $`x`$-axis is

$`M=\left(\begin{array}{cc}\mathrm{cos}(\omega _yT)& \frac{1}{\omega _y}\mathrm{sin}(\omega _yT)\\ \omega _y\mathrm{sin}(\omega _yT)& \mathrm{cos}(\omega _yT)\end{array}\right),`$ (107)
where $`\omega _y`$ is the frequency of the $`y`$-motion and $`T=2\pi /\omega _x`$ is the period of the $`x`$-motion. This is derived by integrating the harmonic oscillator equations of motion for time $`T`$. Then, $`det(MI)=2\mathrm{sin}(\omega _yT/2)`$. In higher dimensions, the monodromy matrix is simply block diagonal so that $`det(MI)=_j2\mathrm{sin}(\omega _jT/2)`$, where the product is over the directions other than $`x`$.

The only role played by the $`x`$ variable above was to specify the time of evolution in the determination of the arguments of the sinusoids. It was not important that it be a harmonic motion, it is enough that it be periodic. These results apply for any periodic orbit with period $`T`$ as long as the transverse motion is harmonic. This exactly describes a heterogeneous orbit. This justifies the amplitude factor in the denominator of Eq. (28) for $`d=1`$. In higher dimensions, the stability is given by both the single-particle monodromy matrix of the evolving particle and by the harmonic motion of the stationary particle about its potential minimum. But these motions are uncoupled and so are simply multiplicative in their combined contribution to the amplitude.

The phase factor can be determined in an analogous manner from the exact harmonic oscillator trace formula. For the primitive orbit along the $`x`$-axis, the phase index is 3. A factor of 2 arises from the two turning points experienced by the periodic orbit in traversing the $`x`$-motion independent of the harmonic motion transverse to the orbit. The remaining factor of 1 can be attributed to the harmonic $`y`$-motion and is related to the sign of the determinant of the monodromy matrix. For heterogeneous orbits, this means that we should simply include a phase factor of $`\pi /2`$ for the transverse harmonic motion in addition to any phase factors from the single-particle motion along the periodic orbit. In higher dimensions, each transverse direction is independent and the phase index is additive. This accounts for the phase factor of $`d\pi /2`$ in Eq. (28). The fact that each transverse direction is uncoupled from all the rest as well as from the single-particle dynamics transverse to the periodic orbit allows us to simply multiply the amplitudes and add the phase factors.

Finally, if the potential is a local maximum in one of the directions, this corresponds to the case of an unstable harmonic oscillator. It is straightforward to show that its contribution to $`det(MI)`$ is $`2\mathrm{sinh}(\omega _yT/2)`$. Furthermore, its phase index is trivially zero since an unstable periodic orbit running along a ridge does not fold back on itself and introduces no caustics in phase space. This fact is also consistent with the trace formula for an unstable harmonic oscillator as described in Ref. . In higher dimensions with a mixture of stable and unstable directions, we continue to multiply the amplitude factors and add the phase indices of the separate directions. This fully accounts for the modifications described immediately below Eq. (29).

## Appendix D Symmetrised $`N`$-particle Thomas-Fermi term

We now discuss the smooth contribution to $`\text{Tr}(\widehat{U}_\tau \widehat{G})`$ and we shall evaluate it using Wigner transforms as in Appendix A.1. We need to determine the smooth approximation to the symmetric/antisymmetric partition function

$$\begin{array}{cc}& Z_\pm (\beta )=\text{Tr}\left(\widehat{P}_\pm \mathrm{exp}(\beta \widehat{H})\right)\hfill \\ & =\frac{1}{N!}\underset{\tau }{}(\pm 1)^{n_\tau }\text{Tr}\left(\widehat{U}_\tau \mathrm{exp}(\beta \widehat{H})\right)\hfill \\ & =\frac{1}{N!}\underset{\tau }{}(\pm 1)^{n_\tau }\frac{\text{d}z}{(2\pi \mathrm{})^{Nd}}(\widehat{U}_\tau )_\text{W}(z)\left(\mathrm{exp}(\beta \widehat{H})\right)_\text{W}(z).\hfill \end{array}$$
(108)
Since each group element can be decomposed into independent cycles,

$$(\widehat{U}_\tau )_\text{W}(z)=\underset{k}{}f_k(z_k),$$
(109)
where $`k`$ indicates the different cycles comprising the group element $`\tau `$, $`f_k`$ is a function we discuss below and $`z_k`$ denotes the phase space coordinates of the $`n_k`$ particles being permuted by that cycle. (For each group element, the unique decomposition into cycles also provides a unique decomposition of the phase space into the subspaces corresponding to the cycles.) The function $`f_k`$ can be specified without loss of generality by choosing to label the particles being permuted by the cycle as $`1,2,\mathrm{},n_k`$ (i.e. $`12`$, $`23`$, $`\mathrm{},n_k1`$) and to leading order in $`\mathrm{}`$ ,

$$f_k(z_k)(2\pi \mathrm{})^{(n_k1)d}\delta (z_1z_2)\delta (z_2z_3)\mathrm{}\delta (z_{n_k1}z_{n_k}).$$
(110)
The first group element of the sum in Eq. (108) is the identity element for which the decomposition into cycles is the trivial one where each particle is in a cycle by itself so that all of the $`f_k`$ are identically unity. Integrating the smooth approximation to the Wigner function of $`e^{\beta \widehat{H}}`$ yields the smooth $`N`$-particle partition function. Using the generalisation of property (85), we observe that the leading-order term of $`\overline{Z}_\pm (\beta )`$ is just the $`N`$th power of the single-particle smooth partition function $`\overline{Z}_1(\beta )^N`$ and under the inverse Laplace transform, this is just the $`(N1)`$-fold convolution integral of the single-particle smooth density of states. The prefactor of $`1/N!`$ comes from the projection operator (41) and we conclude that the identity term is $`O(1/N!\mathrm{}^{Nd})`$. The first correction will come from group elements that consist of one $`2`$-cycle and $`(N2)`$ $`1`$-cycles. The contribution from this class will have the form $`\overline{Z}_1(\beta )^{N2}\overline{Z}_1(2\beta )`$. Compared to the leading-order term, this class contributes to the density of states with relative order $`O(N\mathrm{}^d/2)`$. The factor of $`N`$ is due to the fact that this class has $`N(N1)/2`$ members and they all contribute identically. The factor of $`2`$ comes from the inverse Laplace transform since the argument of one of the single-particle partition functions is $`2\beta `$. The general structure then emerges. For an arbitrary group element, the contribution to the smooth partition function is $`_k\overline{Z}_1(n_k\beta )`$. It contributes to the smooth density of states with relative order $`O(\mathrm{}^{(Nm_\tau )d}w_\tau )`$, where $`m_\tau `$ is the number of independent cycles in the decomposition of $`\tau `$. The factor $`w_\tau `$ is the size of the class (a combinatoric factor which can be found from Eq. (1-27) of Ref. ) divided by a factor arising from the inverse Laplace transform which equals $`_kn_k`$.

As a formal expansion in powers of $`\mathrm{}`$, this may be inconsistent since some of the neglected corrections from the first few group elements may be of more significant order than the leading-order contributions of later group elements. However, for large $`N`$, one could easily imagine that the combinatoric factor $`w_\tau `$ offsets this effect. Keeping the leading-order $`\mathrm{}`$ term of all the group elements then guarantees that one has a good approximation regardless of the relative size of $`1/\mathrm{}`$ and $`N`$.

## Appendix E Higher-order $`\mathrm{}`$ contributions from hetero-orbits in 2-D billiards

In this paper, we only discuss leading-order contributions to the oscillatory part of the density of states. For billiards, hetero-orbit families generally have higher dimensionality than dynamical orbit families and the corrections from the former can be quite significant. The formalism of Ref.  cannot be used to determine these corrections. We calculate correction terms arising from hetero-orbits using the convolution method as in Ref. . The first few terms of the asymptotic series can be determined by convolving the Weyl expansion (13) term-by-term with a two-particle trace formula. As a formal expansion in powers of $`\mathrm{}`$, this is inconsistent since we do not include corrections to the one-particle trace formula (15) itself. However, our numerics seem to indicate that the corrections to the Gutzwiller trace formula are negligible in this case, otherwise we could not reproduce the quantum results with such accuracy.

The contribution from the first type of hetero-orbit where one particle evolves while the others are stationary is calculated from

$`\stackrel{~}{\rho }_3^{\text{h1}}(E)`$ $`=`$ $`\overline{\rho }_1(E)(\overline{\rho }_1\stackrel{~}{\rho }_1)(E)`$ (111)
$`=`$ $`{\displaystyle _0^E}\overline{\rho }_1(\epsilon )(\overline{\rho }_1\stackrel{~}{\rho }_1)(E\epsilon )\text{d}\epsilon .`$
After convolving Eq. (13) with the oscillatory function $`(\overline{\rho }_1\stackrel{~}{\rho }_1)(E)`$ (which has been calculated in Ref. ), we find there are nine integrals to do, but three of these are trivial because of a delta function in the integrand. The remaining six integrals require careful analysis. As an example, we obtain the asymptotic expansion of one integral. The others are calculated in the same manner, but we forego the details. The first correction to the leading-order result (49) comes from two terms; the first (second) is the convolution of the area (perimeter) term of $`\overline{\rho }_1(E)`$ with the perimeter (area) term of $`(\overline{\rho }_1\stackrel{~}{\rho }_1)(E)`$. The integral involved in the first term is

$$I_𝒜(E)=_0^E(E\epsilon )^{1/4}\mathrm{cos}\left(a\sqrt{E\epsilon }+b\right)\text{d}\epsilon ,$$
(112)
where $`a=\sqrt{\alpha }L_\gamma ,b=\sigma _\gamma \pi /2\pi /4`$. We want to perform a local analysis about $`\epsilon =0`$. The reason for this is that small $`\epsilon `$ corresponds physically to the situation where the stationary particles have little energy and most of the energy belongs to the dynamically evolving particle. For $`\epsilon E`$, we have the opposite situation which does not make sense physically. We first make several changes of variable to simplify the calculation. A first change of variable $`u=(E\epsilon )^{1/2}`$ removes the square root in the argument of the sinusoid. Next, we wish to make a Taylor expansion about the point $`u=\sqrt{E}`$ and to facilitate this, we make a second change of variable $`x=\sqrt{E}u`$. The integral then becomes

$$\begin{array}{cc}& I_𝒜(E)=2_0^\sqrt{E}(\sqrt{E}x)^{1/2}\mathrm{cos}\left(a(\sqrt{E}x)+b\right)\text{d}x\hfill \\ & 2\text{Re}\left\{e^{ia\sqrt{E}+b}_0^{\mathrm{}}(\sqrt{E}x)^{1/2}e^{iax}\text{d}x\right\}.\hfill \end{array}$$
(113)
The integral $`_0^\sqrt{E}(\mathrm{})=_0^{\mathrm{}}(\mathrm{})_\sqrt{E}^{\mathrm{}}(\mathrm{})`$; the second integral is an endpoint correction, but asymptotic in $`E`$, this correction is negligible. Thus, we are justified in replacing $`\sqrt{E}`$ with $`\mathrm{}`$ in the second line above. At this point, we obtain the Taylor expansion of $`f(x)=(\sqrt{E}x)^{1/2}`$ about $`x=0`$:

$`f(x)=E^{1/4}{\displaystyle \frac{1}{2}}E^{1/4}x{\displaystyle \frac{1}{8}}E^{3/4}x^2+\mathrm{}.`$ (114)
Since the final correction \[$`I_𝒦(E)`$\] in the expansion (111) is $`O(E^{1/4})`$, it is only necessary to include terms in the Taylor series to $`O(E^{1/4})`$. Thus,

$$\begin{array}{cc}& _0^{\mathrm{}}(\sqrt{E}x)^{1/2}e^{iax}\text{d}x\hfill \\ & E^{1/4}_0^{\mathrm{}}e^{iax}\text{d}x\frac{1}{2}E^{1/4}_0^{\mathrm{}}xe^{iax}\text{d}x\hfill \\ & =E^{1/4}\left(\frac{i}{a}\right)\frac{1}{2}E^{1/4}\left(\frac{1}{a^2}\right)\hfill \end{array}$$
(115)
and asymptotically,

$$\begin{array}{cc}\hfill I_𝒜(E)& \frac{2}{a}\{E^{1/4}\mathrm{cos}(a\sqrt{E}+b\frac{\pi }{2})\hfill \\ & +\frac{E^{1/4}}{2a}\mathrm{cos}(a\sqrt{E}+b)\}.\hfill \end{array}$$
(116)
An equivalent approach is to evaluate the integral exactly and then replace the resulting functions with their asymptotic forms. Evaluating the integral (112) at the upper limit using this method then corresponds to the situation where one of the stationary particles has all of the energy while the dynamically evolving particle has no energy. Physically, this does not make sense and this is evident mathematically since the contribution that comes from evaluating the integral at the upper endpoint is a smooth function of $`E`$ and is therefore spurious in the sense discussed in Refs. . The result is only meaningful if we drop this smooth contribution. This is justified since we know that any smooth contribution to the density of states is already contained in the $`\overline{\rho }_3(E)`$ term. This is completely equivalent to what is done above (the spurious smooth contribution above is the endpoint correction that was dropped in the second line of Eq. (113)). All six convolution integrals can be analysed in this way. Collecting the contributions from all six integrals, the expansion up to $`O(1/\mathrm{}^{3/2})`$ is

$$\begin{array}{cc}\hfill \stackrel{~}{\rho }_3^{\text{h1}}(E)& =\underset{\gamma }{}\frac{L_\gamma ^0}{\sqrt{\left|det\left(\stackrel{~}{M}_\gamma I\right)\right|}}[\frac{\alpha ^{3/2}𝒜^2E^{1/2}}{8\pi ^3L_\gamma ^2}\mathrm{cos}\left(\mathrm{\Phi }\right)\frac{\alpha ^{5/4}𝒜E^{1/4}}{8\sqrt{2}\pi ^{5/2}L_\gamma ^{3/2}}\mathrm{cos}(\mathrm{\Phi }\frac{3\pi }{4})\hfill \\ & +\frac{\alpha }{2\pi ^2L_\gamma }(\frac{𝒜^2}{4\pi L_\gamma ^2}+\frac{^2}{32}+𝒜𝒦)\mathrm{cos}(\mathrm{\Phi }\frac{\pi }{2})\frac{\alpha ^{3/4}E^{1/4}}{8\sqrt{2}\pi ^{3/2}L_\gamma ^{1/2}}(\frac{3𝒜}{8\pi L_\gamma ^2}+𝒦)\mathrm{cos}(\mathrm{\Phi }\frac{\pi }{4})],\hfill \end{array}$$
(117)
where $`\mathrm{\Phi }=\sqrt{\alpha E}L_\gamma \sigma _\gamma \pi /2`$. The contribution from the second type of hetero-orbit where one particle is stationary while the others evolve is calculated from

$`\stackrel{~}{\rho }_3^{\text{h2}}(E)`$ $`=`$ $`\overline{\rho }_1(E)(\stackrel{~}{\rho }_1\stackrel{~}{\rho }_1)(E)`$ (118)
$`=`$ $`{\displaystyle _0^E}\overline{\rho }_1(\epsilon )(\stackrel{~}{\rho }_1\stackrel{~}{\rho }_1)(E\epsilon )\text{d}\epsilon .`$
After convolving Eq. (13) with the formula for $`(\stackrel{~}{\rho }_1\stackrel{~}{\rho }_1)(E)`$ (which has been calculated in Ref.  using stationary phase asymptotics), we find there are only two convolution integrals which require analysis. These are evaluated asymptotically using the same technique as above. The final result (up to $`O(1/\mathrm{}^{3/2})`$) is

$$\begin{array}{cc}\hfill \stackrel{~}{\rho }_3^{\text{h2}}(E)& =\underset{\gamma _1,\gamma _2}{}\frac{L_{\gamma _1}^0L_{\gamma _2}^0\left(L_{\gamma _1}^2+L_{\gamma _2}^2\right)^{1/4}}{\sqrt{\left|det(\stackrel{~}{M}_{\gamma _1}I)\right|}\sqrt{\left|det(\stackrel{~}{M}_{\gamma _2}I)\right|}}[\frac{\alpha ^{5/4}𝒜E^{1/4}}{(2\pi )^{5/2}L_{12}}\mathrm{cos}(\mathrm{\Phi }_{12}\frac{3\pi }{4})\hfill \\ & \frac{\alpha }{16\pi ^2L_{12}^{1/2}}\mathrm{cos}(\mathrm{\Phi }_{12}\frac{\pi }{2})+\frac{\alpha ^{3/4}E^{1/4}}{(2\pi )^{3/2}}(\frac{𝒜}{4\pi L_{12}^2}+𝒦)\mathrm{cos}(\mathrm{\Phi }_{12}\frac{\pi }{4})],\hfill \end{array}$$
(119)
where $`L_{12}=\sqrt{L_{\gamma _1}^2+L_{\gamma _2}^2}`$, $`\sigma _{12}=\sigma _{\gamma _1}+\sigma _{\gamma _2}`$ and $`\mathrm{\Phi }_{12}=\sqrt{\alpha E}L_{12}\sigma _{12}\pi /2`$.