text:
# Stability of the trapped nonconservative Gross-Pitaevskii equation with attractive two-body interaction

## I Introduction

The stability of the condensed state is governed by the nature of the effective atom-atom interaction: the two-body pseudopotential is repulsive for a positive $`s`$wave atom-atom scattering length and it is attractive for a negative scattering length . The ultra-cold trapped atoms with repulsive two-body interaction undergoes a phase-transition to a stable condensed state, in several cases found experimentally, as for 87Rb , 23Na  and 1. However, a condensed state of atoms with negative $`s`$wave atom-atom scattering length (as in case of 7Li ) would be unstable, unless the number of atoms $`N`$ is small enough such that the stabilizing force provided by the zero-point motion and the harmonic trap overcomes the attractive interaction, as found on theoretical grounds . Particularly, in the case of 7Li gas , for which the $`s`$wave scattering length is $`a=14.5\pm 0.4`$ Å, it was experimentally observed that the number of allowed atoms in the Bose condensed state was limited to a maximum value between 650 and 1300, a result consistent with the mean-field prediction , where the term proportional to the two-body scattering length (negative) dominates the nonlinear part of the interaction.

More recently, the maximum critical number of atoms for Bose-Einstein condensates with two-body attractive interactions have been deeply investigated by the JILA group, considering experiments with 85Rb . They have considered a wide tunning of the scattering length, $`a`$, from negative to positive, by means of Feshbach resonance , and observed that the system collapses for a number of atoms smaller than the theoretically predicted number. Their experimental results, when compared with theoretical predictions for spherical traps, show a deviation of up to 20% in the critical number. More precisely, it was shown in Ref.  that part of this discrepancy is due to the non-spherical symmetry that was considered in Ref.. Such a deviation can also be an indication of higher order non-linear effects that one should take into account into the mean-field description. In Ref. , it was considered the possibility of a real and positive quintic term, due to three-body effects, in the Gross-Pitaevskii formalism. A negative quintic term would favor the collapse of the system for a smaller critical number of atoms, as verified in the JILA’s experiments. However, the real significance of a quintic term in the formalism is still an open question.

Our main motivation in the present work is to analyze the dynamics represented by an extension of the mean-field or Gross-Pitaevskii approximation, with non-conservative imaginary terms that are added to the real part of the effective interaction, the two-body nonlinear term with a spherically symmetric harmonic trap. For the imaginary part, the interaction is a combination of a linear term, related to atomic feeding, and a quintic term, due to three-body recombination, that is responsible for the atomic dissipation. This is an approximation that is commonly used to study the properties of Bose-Einstein condensed systems. We consider a wide variation of the nonconservative parameters, in particular motivated by the actual realistic scenario, that already exists, of altering experimentally the two-body scattering length . As it will be clear in the following, this possibility will lead effectively to a modification of the dissipation parameter. By changing the absolute value of the scattering length, from zero to very large absolute values, one can change in an essential way the behavior of the mean-field description. As it will be shown from the present numerical approach, the results for the dynamical observables of the system can be very stable (solitonic-type) or very unstable (chaotic-type); the characteristic of the results will depend essentially on the ratio between the nonconservative parameters related to the atomic feeding and dissipation.

In the next section, we review the formalism. The main results are presented in section III, followed by our conclusions in section IV.

## II Mean field approximation

The mean field approximation has shown to be appropriate to describe atomic Bose-Einstein condensation of a dilute gas of atoms confined by a magnetic trap . In the case of positive scattering length, we have a very good agreement with experimental data, as the thermal cloud is practically absent (removed by cooling evaporation) and almost all the particles are in the condensed state. In this case, the mean field approximation results in a nonlinear Schrödinger equation (NLSE) known as Gross-Pitaevskii equation (GPE). If we have $`N`$ particles trapped in a spherical harmonic potential, this equation is given by

$`i\mathrm{}{\displaystyle \frac{\mathrm{\Psi }}{t}}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}\stackrel{}{}^2+{\displaystyle \frac{1}{2}}m\omega ^2r^2+{\displaystyle \frac{4\pi \mathrm{}^2a}{m}}\mathrm{\Psi }^2\right)\mathrm{\Psi },`$ (1)
where $`\mathrm{\Psi }\mathrm{\Psi }(\stackrel{}{r},t)`$, the wave-function of the condensate, is normalized to the number $`N`$, $`m`$ is the mass of a single atom, $`\omega `$ is the angular frequency of the trap, and $`a`$ is the two-body scattering length.

In this work, we have concentrated our study in the interesting dynamics that occurs when the scattering length is negative ($`a=|a|`$). In this case, it is well known that the system is unstable without the harmonic trap; and the trapped system has a critical limit $`N_c`$ in the number of condensed atoms. The mean-field approximation has also shown to be reliable in determining the critical number of particles and even collapse cycles in the condensate . Actually, systems with attractive two-body interaction are being intensively investigated experimentally , by using the so-called Feshbach resonance . The scattering length can be tuned over a large range by adjusting an external magnetic field (for more details, see Ref. ). Here, we are interested in the dynamics of a realistic system, where we add two non-conservative terms: one (linear) related to the atomic feeding from the non-equilibrium thermal cloud; and another, dissipative due to three-body recombination processes (quintic). It is true that other dissipative terms can also be relevant for an arbitrary trapped atomic system, as a cubic one, that can be related with dipolar relaxation or with an imaginary part of the two-body scattering length. However, in order to simplify the study and better analyze the results, we restrict our considerations to the case that we have just one parameter related with the feeding and another related with dissipation. We have considered only the three-body recombination parameter for dissipation also motivated by the observation that, for higher densities, this term dominates the two-body loss . So, for the generalization of Eq. (1), we add the imaginary terms in the interaction, such that

$`\mathrm{i}\mathrm{}{\displaystyle \frac{\mathrm{\Psi }}{t}}=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\stackrel{}{}^2\mathrm{\Psi }+{\displaystyle \frac{1}{2}}m\omega ^2r^2\mathrm{\Psi }+{\displaystyle \frac{4\pi \mathrm{}^2a}{m}}\mathrm{\Psi }^2\mathrm{\Psi }`$ (3)
$`+\mathrm{i}G_\gamma \mathrm{\Psi }\mathrm{i}G_\xi \mathrm{\Psi }^4\mathrm{\Psi },`$
where $`G_\xi `$ is the dissipation parameter, due to three-body collisions, and $`G_\gamma `$ is a parameter related to the feeding of the condensate from the thermal cloud. The Eq. (3) was first suggested in Ref.  to simulate the condensation of 7Li.

In order to recognize easily the physical scales in Eqs. (1) and (3), it is convenient to work with dimensionless units. By making the transformations

$`\stackrel{}{r}\sqrt{{\displaystyle \frac{\mathrm{}}{2m\omega }}}\stackrel{}{x},t{\displaystyle \frac{\tau }{\omega }},`$ (4)
$`G_\gamma {\displaystyle \frac{\gamma }{2}}\mathrm{}\omega ,G_\xi 2\xi \left({\displaystyle \frac{4\pi |a|\mathrm{}}{m\omega }}\right)^2\mathrm{}\omega \text{and}`$ (5)
$`\mathrm{\Phi }\mathrm{\Phi }(x,\tau )\sqrt{8\pi |a|}|\stackrel{}{r}|\mathrm{\Psi }(\stackrel{}{r},t),`$ (6)
we obtain the radial dimensionless $`s`$wave equation:

$`\mathrm{i}{\displaystyle \frac{\mathrm{\Phi }}{\tau }}`$ $`=`$ $`\left[{\displaystyle \frac{d^2}{dx^2}}+{\displaystyle \frac{x^2}{4}}{\displaystyle \frac{|\mathrm{\Phi }|^2}{x^2}}2\mathrm{i}\xi {\displaystyle \frac{|\mathrm{\Phi }|^4}{x^4}}+\mathrm{i}{\displaystyle \frac{\gamma }{2}}\right]\mathrm{\Phi }.`$ (7)
As $`\mathrm{\Psi }(\stackrel{}{r},t)`$ is normalized to the number of atoms $`N(t)`$ in Eq. (3), the corresponding time-dependent normalization of $`\mathrm{\Phi }(x,\tau )`$ is given by the reduced number $`n(\tau )`$:

$$_0^{\mathrm{}}𝑑x|\mathrm{\Phi }(x,\tau )|^2=n(\tau )2N(t)|a|\sqrt{\frac{2m\omega }{\mathrm{}}}.$$
(8)
The nonconservative GPE (7) is valid in the mean-field approximation of the quantum many-body problem of a dilute gas, when the average inter-particle distances are much larger than the absolute value of the scattering length; and also when the wave-lengths are much larger than the average inter-particle distance. The nonconservative terms are important when the condensate oscillates, fed by the thermal cloud, while losing atoms due to three-body inelastic collisions, which happen mainly in the high density regions.

In order to verify the stability and the time evolution of the condensate, as observed in Refs. , two possible relevant observables are the number of particles normalized by the critical number of atoms of the static case ($`N(t)/N_c`$) and the mean square radius (msr),

$`r^2(t)`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}}{2m\omega }}\right){\displaystyle \frac{1}{n(\tau )}}{\displaystyle _0^{\mathrm{}}}𝑑xx^2|\mathrm{\Phi }(x,\tau )|^2`$ (9)
$``$ $`\left({\displaystyle \frac{\mathrm{}}{2m\omega }}\right)x^2(\tau )\left({\displaystyle \frac{\mathrm{}}{2m\omega }}\right)X^2,`$ (10)
where $`X\sqrt{x^2(\tau )}`$ is the dimensionless root mean square radius (rmsr). In our analysis of stability, we calculate the time evolution of these quantities. We explore several combinations of the dimensionless nonconservative parameters $`\xi `$ and $`\gamma `$. We first consider the case in which the atomic feeding is absent or when its parameter is smaller than the atomic dissipation parameter. Next, we explore variations of both parameters of about five orders of magnitude, from $`10^5`$ to $`10^1`$. This wide spectrum includes the parameters considered by Kagan et al. , as well as other combinations that can be considered more realistic due to recent experimental results .

Actually, the relevance of a wider relative variation of the nonconservative parameters $`\gamma `$ and $`\xi `$, presented in Eq. (7), can be better appreciated in face of the experimental possibilities that exist to alter the two-body scattering length . As one should note from Eq. (5), any variation of the scattering length will also affect the effective dissipation parameter $`\xi `$ and, consequently, its relation with the feeding parameter $`\gamma `$. This implies that, by changing the value of the scattering length, from positive to negative, and from zero to very large absolute values, one can change in an essential way the behavior of the mean-field description. In the present work, we are concerned with negative two-body scattering length, where the collapsing behavior of the Eq. (7) shows a very interesting dynamical structure. Even considering the possible limitations on the validity of the mean-field approach after the first collapse (in cases of parametrization where it can occur), it is worthwhile to verify experimentally the behavior of a system in such a situation, by varying $`|a|`$. At least, one can verify how far the theoretical description can be qualitatively acceptable.

As already verified for systems with attractive interaction, as the 7Li, it has been possible, via the mean-field approach, to describe properties like the critical number of atoms in the condensate and growth and collapse cycles ; besides, in the long time evolution, for certain sets of parameters, the calculations have also shown the presence of strong instabilities of the condensate, with signals of spatio-temporal chaotic behavior.

In order to characterize a chaotic behavior, it is necessary to show that the largest Lyapunov exponent related with the solutions of the equation is positive. We follow the criterion used by Deissler and Kaneko  to characterize spatio-temporal chaos. This criterion prescribes that the largest Lyapunov exponent for the system, in an arbitrary time interval, is obtained by plotting the logarithm of a function $`\zeta `$, which is defined by

$`\zeta (\tau )`$ $``$ $`\left({\displaystyle _0^{\mathrm{}}}|\delta \mathrm{\Phi }(x,\tau )|^2𝑑x\right)^{1/2}.`$ (11)
$`\delta \mathrm{\Phi }(x,\tau )`$ will give us the separation between two nearby trajectories; it is obtained in the following form: we numerically evolve in time an initial $`\mathrm{\Phi }_0(x)`$, obtaining $`\mathrm{\Phi }(x,\tau )`$. Independently, we evolve $`\mathrm{\Phi }_0(x)+ϵ(x)`$, and get $`\mathrm{\Phi }^{}(x,\tau )`$, where $`ϵ(x)`$ is a very small random perturbation. $`\delta \mathrm{\Phi }(x,\tau )`$ is given by $`\mathrm{\Phi }^{}(x,\tau )\mathrm{\Phi }(x,\tau )`$. The chaotic behavior is characterized by a positive slope of $`\mathrm{ln}\zeta (\tau )`$, which gives the largest Lyapunov exponent .

## III Numerical Results

In the next, we present the most significant results that characterize the time evolution of the normalized number of particles ( $`N(t)/N_c`$ ), the dimensionless mean square radius $`x^2`$, and, in order to characterize the stability of the system, the function related to the largest Lyapunov exponent. Further, we present a representative case of the phase-space for the root mean square radius. We have studied a wide region of parameters $`\gamma `$ and $`\xi `$, covering about five orders of magnitude, from $`10^5`$ to $`10^1`$, including the case with no feeding ($`\gamma =0`$).

In order to have a clear and useful map of the regions where one should expect stable results, as well as regions with instabilities or chaos, we summarize the present numerical results in a diagrammatic picture that relates these two nonconservative parameters. In general, it is expected that the system is more stable when the parameter related to the feeding of atoms ($`\gamma `$) from the thermal cloud is significantly smaller than the parameter related to the dissipation ($`\xi `$). However, it is interesting to find out the region of parameters where this transition (from stable to unstable results) occurs. Analysis of experimental results can provide a test to the present mean-field description in case of negative two-body interaction. As previously observed, we are considering dimensionless observables and parameters. For any realistic comparison with experimental parameters, one should convert $`\gamma `$ and $`\xi `$ to the parameters $`G_\gamma `$ and $`G_\xi `$, as given in Eq. (5).

The numerical solutions of Eq. (7) were obtained by applying the semi-implicit Crank-Nicolson algorithm for nonlinear problems, as implemented in Ref. . This method is stable and, therefore, very convenient and reliable to treat time-dependent non-linear partial differential equations. The initial condition for the number of atoms $`N`$ in the condensate was such that $`N(0)/N_c=n(0)/n_c=0.75`$. The evolution of the observables have been extended up to $`\tau =\omega t=500`$.

In general, as expected, the smaller is the dissipation parameter, the longer is the life of the condensate. The mean square radius presents an oscillatory behavior while one increases $`\xi `$. One observes that, in the regime of small feeding ($`\gamma 10^4`$), the extended Lyapunov presents no positive slope. For larger values of $`\gamma `$, from $`10^3`$ and $`10^2`$, we have studied a few cases where the interplay between the nonconservative behaviors are significant.

In Fig. 1, we show the dynamical behavior of the number of atoms for $`\gamma =10^2`$ and several values of $`\xi `$; and, in the Fig. 2, the corresponding time evolution of $`x(\tau )^2`$. We realize an interesting behavior, that occurs when the dissipation is larger than the feeding process: there are solutions of stability or dynamical equilibrium between both nonconservative processes. This phenomenon was already discussed in Ref. , for a few values of the dissipation and feeding parameters, using the time dependent variational approach and also the Crank-Nicolson method. In the present work, we observe a wide region of parameters where it is possible the formation of autosolitons . However, when the feeding process is much larger than the dissipation, of about one or more orders of magnitude, we can also observe chaotic behaviors. See, for example, the case with $`\xi =10^3`$.

The time evolution of the number of particles, represented in Fig. 1, shows a collapse for $`\tau `$30, followed by several other collapses, with the number of particles going above the critical limit $`N_c`$. So, after a sequence of collapses, the critical limit for the number of particles is no more followed, as already shown in Ref. .

The corresponding time evolution of $`x(\tau )^2`$ is shown in the upper frame of Fig. 2. We observe that, following each collapse, after the shrinking of the system, the radius is multiplied by a large factor, with indication of being populated by radial excited states. In the lower frame of Fig. 2, we can observe the corresponding transition from the stable region (where the system finds the equilibrium at a fixed value of the radius, corresponding to autosoliton formation) to the unstable region. As shown, the instability starts to occur when $`\xi =2\times 10^3`$, and it can be developed to a spatio-temporal chaos. The chaotic behavior can be verified through the Deissler and Kaneko criterion .

In Fig. 3 we illustrate the application of the Deissler and Kaneko criterion to the system given by Eq. (7), for a fixed value of the feeding parameter $`\gamma =0.01`$, and a set of values of the dissipation parameter $`\xi `$. It was plotted the time evolution of the function ln$`(\zeta )`$, where $`\zeta `$ is given by Eq. (11), following the prescription given in Ref.  to obtain the largest Lyapunov exponents for the system. Within this prescription, the system becomes chaotic when ln$`(\zeta )`$ has a positive slope. As shown in Fig. 3, this clearly occurs, for example, when $`\xi =10^4`$. In case of $`\xi =10^5`$ we note a much faster increasing in ln$`(\zeta )`$, with an observed saturation that happens due to the fact that such function has reached the maximum separation between the trajectories. The saturation properties is also verified when studying chaotic behaviors in ordinary differential equations . The plot of ln$`(\zeta )`$ corresponds to the same value of $`\gamma `$ ($`=10^2`$) used in Figs. 1 and 2. As shown, a clear characterization of chaotic behaviors starts to occur only for values of the dissipation parameter $`\xi `$ much smaller than $`\gamma `$. In the cases presented in Fig. 3, for $`\xi 10^3`$.

In Fig. 4, we present another significative illustration of chaotic behavior, through the phase-space behavior of the mean-square-radius, considering one case that was characterized as chaotic by using the Deissler-Kaneko criterion. We have plotted in this figure the root mean-square-radius phase space for the case with $`\gamma =0.01`$ and $`\xi =10^4`$. The irregular behavior of the trajectories, observed in Fig. 4, with the classical strange attractors being observed, clearly resembles chaos. This behavior is similar to the chaotic behavior observed in ordinary cases .

As a general remark that one can make from the presented results, we should note that, in order to observe unstable chaotic behaviors, the dissipation must be much smaller than the feeding parameter.

In a diagrammatic picture, given in Fig. 5, we resume our results. We show the relation between the two nonconservative parameters, $`\xi `$ and $`\gamma `$, in order to characterize the parametric regions where one should expect stability or instability in the solutions for the Eq. (7). The stable results of the Eq. (7) are represented by bullets; the nonstable results that clearly present positive slope for $`\mathrm{ln}\zeta (\tau )`$ (chaotic behavior) are represented by empty squares; with $`\times `$, we show other intermediate nonstable results, in which the characterization of chaotic behavior was not so clear, through the Deissler-Kaneko criterion. In this figure, in order to observe the approximate consistency of the numerical results, we also include the variational analysis presented in Fig.1 of Ref. , represented by the dashed-line. It is separating the stable region (upper part) from the unstable one (lower part).

We should note that, in section V of Ref., it was also considered the dynamics of growth and collapse, with non-conservative terms related to feeding ($`\gamma _0`$) and dissipation ($`\gamma _1`$ and $`\gamma _2`$) in a specific example. For the dissipation they have also considered a term related to dipolar relaxation, given by $`\gamma _2`$. Here, in our systematic study of the regions of instability, we took into account previous experimental  observations that the dominant process for the dissipation is the three-body recombination. By comparing the parameters of Ref.  with the parameters that we have used, and observing that our parameter $`\xi `$ should be related to both dissipation parameters used in Ref.  ($`\gamma =\gamma _0=2.6\times 10^3`$, $`\xi 10^5`$) one can verify from the results given in Fig. 5 that the model of Ref.  is inside the intermediate region, where the system is unstable, without a clear signature of chaos.

## IV Conclusions

In summary, we have studied the dynamics associated with the extended nonconservative Gross-Pitaevskii equation for a wide region of the dimensionless nonconservative parameters, $`\xi `$ and $`\gamma `$, that, respectively, are related to atomic dissipation and feeding in a trapped atomic condensed system. We consider systems with attractive two-body interaction in a spherically symmetric harmonic trap. In Fig. 5, we resume our results, by mapping the space of $`\gamma `$ versus $`\xi `$, showing the regions of equilibrium and the regions of instability, as well as the regions where we are able to characterize chaotic behaviors, using a criterion given in Ref. . It was also confirmed that chaotic behaviors occur mainly when $`\gamma `$ is big enough and $`\gamma /\xi `$ is large (at least, when $`\gamma `$ is one or two orders of magnitude larger than $`\xi `$). A wide variation of the nonconservative parameters was analyzed, in particular motivated by the actual realistic scenario, that already exists, of altering experimentally the two-body scattering length . By changing the absolute value of the scattering length, one can change in an essential way the behavior of the mean-field description.

## Acknowledgements

We thank Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) for partial support. LT and TF also thank partial support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).


format: text

text:
# Statistical properties of one dimensional “turbulence”.

## 1 Introduction

Although studies of turbulence have been carried out for decades, they go on now-days because the problem is very hard and challenging but new ideas have emerged in the last few years from the investigations of statistical properties of turbulence, both at a local or a global scale. First, studies of local longitudinal velocity differences measured at distance $`r`$ showed that their probability densities are very well described by a theoretical model based on nonextensive statistical mechanics , which, in the context of turbulence, can be justified by a phenomenological model based on dynamical foundations . Second, interesting properties showed up in the analysis of the statistical properties of global quantities such as the power consumption of the confined turbulent Von Kármán flow  and they culminated with the observation of the universality of rare fluctuations in turbulence and critical phenomena .

These recent developments are based on an analysis of experimental results . Indeed only such data contain the full features of turbulence, but for an understanding of the results it would be desirable to have a simple model, showing the basic properties of turbulence, that would be amenable to very detailed studies. The search for simple models of turbulence lead to the introduction of the so-called shell models  which describe the energy exchange mechanism between various scales. However shell models are not suitable to study properties that depend on space or to analyze how the dynamics can lead to the observed statistics.

This is why it might be interesting to introduce models of turbulence that are sufficiently simple to allow a detailed study, and nevertheless preserve a spatial structure. In this letter we introduce such a model derived from studies on one-dimensional nonlinear lattices. These systems show various features which are also found in turbulence, and, in particular an exchange of energy between various spatial scales which tends to cause a flow of energy from large scales (long wavelengths) to smaller scales, through the mechanism of modulational instability . The idea to use one-dimensional systems as simple models for turbulence is not new and continuous one-dimensional media have been introduced to study wave turbulence , but these studies only lead to weak turbulence. Discrete lattice are attractive for two reasons. First they are simpler since they contain only a limited number of degrees of freedom (and therefore modes), but, more importantly because they have exact localized solutions which could play the role of vortices, the discrete breathers , opening the possibility to study an analogous of fully developed turbulence.

Many interesting experimental results on turbulence have been obtained with a simple model flow, the Von Kármán flow produced in the gap between coaxial disks . We investigate here a nonlinear one-dimensional lattice excited by both ends, in a configuration which drives a large scale structure in the system, similar to the vortex of the Von Kármán flow and we show that it exhibits remarkable similarities with turbulence, both for local and global properties.

## 2 The one-dimensional lattice model.

The system that we study is a one-dimensional chain of anharmonic oscillators, coupled to each other by harmonic interactions. We have chosen Morse oscillators because, the presence of a cubic term in its expansion leads to a suitable nonlinear coupling between the modes, and moreover the lattice of Morse oscillators has localized oscillatory solutions which have interesting properties . In analogy with the Von Kármán flow, we want to drive a finite lattice from both ends, i.e. to inject energy in the system through an external driving. As in actual turbulence experiment where the viscosity of the fluids contributes do dissipate energy, we also need to introduce a dissipation mechanism so that the system can reach a steady state. In analogy with the dissipation term $`\nu ^2v`$ of the Navier Stokes equations, where $`\nu `$ is the kinematic viscosity, we introduce a dissipation term proportional to the second order finite difference of the velocity of the oscillators, so that the equations of motion of the nonlinear lattice are

$$\ddot{y}_nK(y_{n+1}+y_{n1}2y_n)2\omega _0^2e^{y_n}\left(e^{y_n}1\right)\nu (\dot{y}_{n+1}+\dot{y}_{n1}2\dot{y}_n)=0,$$
(1)
where $`y_n`$ denotes the position of the $`n^{\mathrm{th}}`$ oscillator, $`K`$ is the coupling constant between neighboring sites, $`\omega _0^2`$ determines the depth of the Morse on-site potential, and $`\nu `$ is a coefficient which plays the role of the kinematic viscosity. As usual, the dot designates a time derivative.

Similarly to the Von Kármán flow where the medium is excited from both ends by a set of rotating disks which create a large scale structure in the system, we consider a finite lattice of $`N+2`$ particles and we excite it by imposing the displacements of the two end particles, particle $`0`$ and particle $`N+1`$. Their motion is chosen according to a linear mode of the lattice, with an amplitude $`A`$, a wavevector $`\kappa =p\pi /(N+3)`$, where $`1<p<N+1`$ is an integer, and a frequency $`\omega =\sqrt{2\omega _0^2+4K\mathrm{sin}\kappa /2}`$ (assuming that particles $`1`$ and $`N+2`$ are fixed). Therefore the driving of the system is defined by two parameters, $`p`$ which fixes its wavevector and frequency and $`A`$ which determines its amplitude. All simulations start from a lattice at rest and the equations of motion are integrated with a fourth-order Runge-Kutta scheme.

Indeed, for such a nonlinear lattice, the fact that the driving corresponds to a single mode of the system does not imply that its response will be merely a single mode. Similarly to the Von Kármán flow, driven by two co-rotating disks, which has a very complex response, the energy density nonlinear lattice also shows a complex spatio-temporal pattern, as displayed in Fig. 1. Localized modes that move in the system, interact with each other, and die out due to dissipation are observed, and they appear to play a role very similar to the vortex filaments in fluids. The point of this study is to observe and analyze the response of the nonlinear lattice to a simple excitation. More complex cases could easily be imagined, such as an excitation with two different frequencies at both ends, which would be similar to some fluid experiments where the two disks have different velocities, but, as shown below the one-frequency driving is sufficient to yield interesting and puzzling results.

## 3 Local analysis.

Numerical investigations of spatial the Fourier spectrum of the energy density $`e(\kappa ,t)`$ in this simple model have exhibited properties analogous to turbulence, such as a power law spectrum $`e(\kappa )=E_0\kappa ^\alpha `$ with $`\alpha 2.2`$ which is remarkably insensitive to the parameters of the model or of the excitation in a very broad range of parameters . When the excitation amplitude increases even more, a transition to a spectrum decaying exponentially versus $`\kappa `$ is found when the large scale structure created by the excitation breaks down into many smaller scale localized modes, as shown in Fig. 1. This is the regime that we study in this letter. In actual turbulence experiments the energy density cannot be directly measured, but the local velocity is easily accessible. A map of the velocity field would require an array of probes that perturb the flow, and this is why the investigations only evaluate velocity differences measured in 2 points separated by a distance $`r`$, which provide some information on the spatial structure of the flow, and are also related to the viscous dissipation. For the nonlinear lattice it is easy to calculate an equivalent quantity in the numerical simulations $`u(\delta n,t)=\dot{y}_{n_0+\delta n}(t)y_{n_0}(t)`$, where $`n_0`$ is a fixed reference site and $`\delta n`$ is an index that measures the distance between the two points, in units of the lattice spacing. It should be stressed that this quantity is local because there is no spatial averaging over the reference site position $`n_0`$.

In a typical numerical experiment, we simulate the dynamics of the system for a time interval of $`210^7`$ time units (t.u.). The first $`510^5`$ t.u. are discarded in order to allow the system to reach a steady state before starting the recording of $`u(\delta n,t)`$. Then we store $`u(\delta n,t_i)`$ at times $`t_i`$ separated by an interval $`t_{i+1}t_i=40`$ t.u. which is large enough to ensure that successive values are statistically independent. This gives us a set of 375000 values $`u_i`$ for each $`\delta n`$. We evaluate the reduced probability density of this distribution by first calculating $`\stackrel{~}{u}_i=(u_iu_i)/\sigma `$ where $``$ designates a time average and $`\sigma =\sqrt{u_i^2u_i^2}`$ is the variance of the distribution. Then we make a normalized histogram of the $`\stackrel{~}{u}_i`$ values in 100 bins to get the normalized probability density $`\stackrel{~}{P}(\stackrel{~}{u})`$.

Figure 2a shows such histograms for $`\delta n`$ varying from $`2`$ to $`40`$. The analogy with similar figures deduced form experiments in fluid turbulence is striking (see for instance Fig. 1 in Ref. ). The same scenario from almost Gaussian distribution at large distances to stretched distributions at small distances is found. The analogy can be made more quantitative by analyzing the shape of the histograms. The curves of Fig. 2a are very well fitted by the probability distribution deduced from nonextensive statistical mechanics

$$\stackrel{~}{P}(\stackrel{~}{u})=\frac{1}{Z_q}\frac{1}{\left[1+\frac{1}{2}\beta (q1)\stackrel{~}{u}^2\right]^{\frac{1}{q1}}},$$
(2)
with $`\beta =2/(53q)`$, determined by the condition that the reduced probability distribution should have variance 1. As $`Z_q`$ is a normalization constant, it should be noticed that each curve is fitted with only one free parameter, the parameter $`q1`$, which, in the formalism of nonextensive statistical mechanics measures the degree of nonextensivity. With $`q1`$ ordinary statistical mechanics is recovered and the probability distribution tends to a Gaussian. In the case of turbulence, a formula that generalizes (2) has been derived to take into account the slight skewness of the distributions . The distributions obtained for the nonlinear lattice also have a slight skewness so that the extended formula improves the fit. However, as it is already very good with Eq. (2) we sticked to the expression that minimizes the number of parameters.

The variation of $`q`$ versus $`\delta n`$ is shown in Fig. 2b for two values of the viscosity $`\nu `$ and a mean amplitude $`y1.110^2`$ in both cases. As for the shape of the probability densities, the similarity with the experimental results on turbulence is strong. At small distances $`q`$ is significantly larger than 1 ($`q1.3`$ in our case) and it decays toward a value close to 1, but above 1, at large distances ($`q1.08`$ in our case). In turbulence experiments the variation of $`q`$ versus distance is essentially determined by the rescaled distance $`r/\eta `$, where $`\eta `$ is the Kolmogorov scale which is the typical size of the smallest vortices. In the discrete nonlinear lattice, the size of the smallest breathers is of the order of the lattice spacing, and, in the low viscosity range (high “Reynolds number”) that we are investigating here, it does not vary significantly with the viscosity. This is why $`q`$ only depends on $`\delta n`$ for the lattice as shown in Fig. 2b.

## 4 Global properties of the dynamics.

Our aim in this section is to examine to what extend the one-dimensional nonlinear lattice that we have investigated also shows the same “universal” statistical properties that where observed in critical magnetic systems, self-organized criticality or turbulent flows when one investigate global quantities.

In actual turbulence experiments, a global quantity which can be easily evaluated is the power injected by the engines driving the flow, for instance through the disks in the Von Kármán experiment. In the numerical simulations we have however the freedom to evaluate other global quantities directly in the bulk of the system. The power dissipated in the bulk

$$P_{\mathrm{d}iss}=\underset{n=N/4+1}{\overset{n=3N/4}{}}F_v(n)\dot{y}_n=\nu (\dot{y}_{n+1}+\dot{y}_{n1}2\dot{y}_n)\dot{y}_n$$
(3)
(where $`F_v`$ is the viscous force) is particularly interesting because it probes the dissipation, an essential physical process for turbulence, at a global scale without being affected by the boundary effects. To avoid any perturbation due to boundaries, we evaluate $`P_{\mathrm{d}iss}`$ by summing all the dissipative contribution from $`n=N/4+1`$ to $`3N/4`$, i.e. in the central half of the system.

As $`P_{\mathrm{d}iss}`$ is a global quantity, one could think that it would be smooth and featureless due to averaging. This is not the case. In addition to small amplitude fluctuations, the signal shows sharp peaks corresponding to strong bursts of dissipation. This feature shows up on the normalized probability distribution function (PDF) obtained from an histogram of the recorded values of $`P_{\mathrm{d}iss}`$ in numerical simulations. Figure 3 show such histograms plotted in logarithmic scales as a function of $`\theta =(P_{\mathrm{d}iss}P_{\mathrm{d}iss})/\sigma `$ where $`P_{\mathrm{d}iss}`$ is the time average of the dissipated power and $`\sigma _{\mathrm{d}iss}`$ its standard deviation. The dissipation bursts give rise to a highly asymmetric PDF with a large tail on the negative side. The dissipation bursts are associated to collisions of the localized excitations, the discrete breathers, that are visible on Fig. 1. Using an approximate expression of the breathers we can evaluate, in the continuum limit, the time average of the power dissipated by a breather. It shows that the dissipation grows with the third power of the breather amplitude, so that the large peaks observed during breather collisions generate strong bursts of dissipation.

Figure 1 showing the space-time evolution of the energy density illustrates a typical situation where one observes localized structures having a long life-time and a complex pattern of interactions. They travel on long distances introducing long range correlations in the dynamics of the lattice. In the context of the general studies of correlated systems which have been performed recently , systems showing such correlations have been called “inertial systems” and they display a strongly non Gaussian probability distribution function of some global quantities. As shown by Fig. 3, this is also true for the nonlinear lattice that we investigate here. The strong asymmetry of the PDFs is similar to the observations made on turbulence and critical phenomena . In order to evaluate this quantitatively, we have fitted the numerical PDFs with a generalized Gumbell distribution

$$\mathrm{\Pi }(\theta )=w\mathrm{exp}\left[a\left(b\theta e^{b\theta }\right)\right],$$
(4)
where $`w`$ is the normalization factor and $`a`$, $`b`$ are parameters of the fit. Figure 3 shows that this expression gives a very good fit of the observed distributions, and is able to describe properly their evolution as a function of the system size $`N`$. In the limit $`b0`$ the distribution (4) tends to a Gaussian distribution $`\mathrm{\Pi }(\theta )=w^{}\mathrm{exp}[\theta ^2/(2\sigma _0^2)]`$ with a standard deviation $`\sigma _0=1/(ab^2)`$. The value $`b=1`$ corresponds to the standard Gumbell distribution observed in a turbulent flow , and larger values of $`b`$ correspond to even larger deviations from Gaussianity. Figure 3 shows that, even for large systems ($`N=4096`$) the fluctuations of $`P_{\mathrm{d}iss}`$ stay highly non Gaussian.

The origin of this non-Gaussian behavior can be understood from the same argument as for the case of turbulence or critical phenomena : there are long range correlations in the system which cannot be divided into statistically independent regions. However the origin of these long range correlations appears to be different. While in the vicinity of a phase transition correlations are due to the existence of structures at all scales, in the present case, we observe small scale structures (localized modes) which can move on long distances and preserve their identities in collisions (see Fig. 1). However the localized modes are not solitons and are therefore perturbed by the collisions, which limits the correlations to a finite range. As a result the PDF of the dissipated power depends on the size of the system and on the density of localized modes, both of which determine the average number of collisions that a localized excitation experiences when it travels inside the part of the lattice in which we compute the dissipated power (i.e. along $`N/2`$ sites).

The effect of the size of the system is shown in Fig. 3b which shows that $`b`$ decreases with $`N`$ approximately with a power law $`bN^{0.67}`$ for $`128N4096`$ so that the approach of a Gaussian distribution for large systems is very slow. The study of the effect of the mean amplitude of the displacement, which is correlated to the density of localized modes in the system, shows that a strong non-Gaussianity of the power dissipated in the system persists for variations of $`y`$ that extend on more than two orders of magnitude.

Consequently, contrary to the case of turbulence where a universal form of the fluctuations of a global quantity, independent of the Reynolds number, was found , the shape of the PDF of the power dissipated in the nonlinear lattice changes with the size of the system or amplitude of the excitation. However it stays strongly non-Gaussian and well described by a generalized Gumbell distribution for all the cases that we could investigate. One may wonder what are the differences between the turbulent flow and our model system. One possibility is the existence in the turbulent flow of structures at all scales, including the integral scale, i.e. flow patterns that have a characteristic size equal to the size of the system, which is not true for the nonlinear lattice. One may ask however whether larger Reynolds number would not break the large vortex into smaller filaments, leading to a gradual change in the shape of the PDF. Such a regime, if it exists, is probably beyond experimental possibilities, while, for the one-dimensional nonlinear lattice, the simulation can examine a very broad range of dynamical behaviors.

## 5 Conclusion

In this work we examined a simple model related to turbulence. The statistics of its local and global properties are remarkably similar to the experimental results of actual turbulence so that it could be useful for its understanding because it is sufficiently simple to allow detailed studies and test ideas. The remarkable features of the local and global properties could perhaps be understood in a single framework. Global properties of turbulence have been analyzed by Bramwell et al.  by pointing out that the flow cannot be divided into mesoscopic regions which are statistically independent. Such a “long range interaction” is precisely what leads to non-extensive thermodynamics  which models so well the local properties.

Acknowledgements.

Part of this work has been supported by the EU Contract No. HPRN-CT-1999-00163 (LOCNET network). M.P. wants to thank T. Dauxois and B. Castaing (ENS Lyon) for helpful discussions and E.D.G. Cohen (Rockefeller University) for an illuminating lecture that motivated part of this work.


text:
# Correlation functions of the BC Calogero–Sutherland model

## I Introduction

It is a well known fact that there is an intimate relation between the one–dimensional quantum problem with the inverse–square interaction potential, i.e. Calogero–Sutherland model (CSM) Cal ; Sut , and Dyson random matrix ensembles (RMEs) Meh . On the most elementary level the correspondence goes as follows: for the three particular values of the coupling constant ($`\lambda =1/2,1,2`$) square of the many–body ground–state wave function of the CSM coincides with the joint probability distribution (JPD) of RMEs with Dyson index $`\beta =2\lambda `$. Consequently the knowledge of the RME correlation functions may be immediately translated to the information about the CSM. Historically this correspondence proved to be very fruitful for advancing the understanding of the CSM.

It was later realized that Dyson’s classification of the RMEs was not exhaustive. Studies of the two–sublattice model, mesoscopic transport, and QCD Dirac spectra initiated the introduction of their “chiral” counterparts Gad ; SN ; SV ; ForNPB ; Ver . Subsequently, based on Cartan’s classification of Riemannian symmetric spaces, Altland and Zirnbauer AZ ; Zir2 ; AZ97  have further added six “superconducting” chiral symmetry classes. These chiral classes are characterized by the special role played by the zero energy. Namely, there is a certain number of eigenvalues that must have zero energy, while all other eigenvalues occur in symmetric pairs (mirror images) around zero. As a result, the mean density of states (DoS) exhibits either a hollow or a bump around the zero energy followed by decaying oscillations at larger distances. Such a structure on the level of the mean DoS does not show up in the Dyson classes and, as we shall explain below, may be called Friedel oscillations.

The question is whether one can find an appropriate generalization of the CSM whose ground–state wave function possess the same reflection symmetry property as the JPD of chiral ensembles. The answer is known to be affirmative. Indeed, one can write down an integrable one–dimensional model with inverse–square interaction, reflection symmetry and special single–particle potential centered at zero having the required ground state. It is known as the BC–type CSM in the literature OP . Due to the presence of the mirror boundary and localized single–particle potential (an impurity) the model lacks translational invariance and the resulting ground–state density is not uniform. In particular the density profile develops the Friedel oscillations far enough from the impurity. Accordingly one has a unique example of integrable strongly interacting models that exhibit Friedel phenomena. That gives one a possibility to gain the exact information on the amplitude decay rate, spectral characteristics and phase shifts of the Friedel oscillations in the interacting system. Despite the proven integrability, explicit form of the correlation functions of the BC–CSM with generic coupling constant were not established, except for some partial results ForJMP ; YKY2 ; AZ97 ; FM .

The purpose of this paper it to fill this gap for rational values of the coupling constant. To this end we employ the recently developed approach based on the replica trick KM1 ; YL ; NK . It was previously tested on the pair correlation function of the ordinary A–CSM GK , where it perfectly agrees with the exact results of Haldane Hal  and Ha Ha  for any rational coupling constant. The idea is to explore the relation with the RMEs, where the replica trick was found to be accurate in the asymptotic regime. We thus develop first the fermionic replica approach to the chiral symmetry classes of RMEs. Not surprisingly, we are able to reproduce the asymptotic behavior of the known DoS profiles for chiral RMEs. We then extend the treatment away from the RMEs values of the coupling constant and obtain closed analytic results for any rational two–body coupling constant $`\lambda `$ and any impurity phase shift.

We found that for rational values of $`\lambda =p/q`$ ($`p,q`$ coprime), the spectrum of the Friedel oscillations contains exactly $`p`$ harmonics, corresponding to $`2k_F,4k_F,\mathrm{},2pk_F`$ density oscillations. The $`l`$–th harmonic ($`l=1,\mathrm{},p`$) decays algebraically as $`\theta ^{l^2/\lambda }`$, where $`\theta `$ is the distance from the impurity. The amplitudes of the harmonics depend on the number, $`l`$, and the coupling constant, $`\lambda `$, but are not sensitive to the strength and details of the impurity expressed through a phase shift, $`\nu `$. Moreover, the harmonics amplitudes are closely related to those of the two–point correlation function of the homogenous A–CSM. We provide a conformal field theory account for this fact. The impurity phase shift, $`\nu `$, affects only phases of the harmonics in the asymptotic regime. We explicitly show that the relation between the total charge expelled by the impurity and the asymptotic phase shifts, known as the Friedel sum rule, holds for the interacting system. We also compute the distribution of the locus of the particle nearest to the mirror boundary point for the BC–CSM at generic values of the coupling constants.

The paper is organized as follows: in section II we briefly introduce chiral ensembles, and develop the appropriate fermionic replica $`\sigma `$–models. Section III is devoted to the introduction and replica treatment of the BC–CSM. In section IV we perform the analytical continuation and the replica limit, and extract the density profile. In section V we compare this result with the effective conformal field theory description. In section VI we compute the nearest particle distribution. A summary and discussions are provided in section VII. Technicalities of the calculations are relegated to the two appendices.

## II Chiral Circular Ensembles

A circular RME is defined as an ensemble of unitary matrices $`U`$ representing a Riemannian symmetric space $`𝒟`$, stochastically distributed according to the Haar measure $`dU`$ of $`𝒟`$. Expressing a symmetric space as a coset $`𝒟=G/H`$ with a compact Lie group $`G`$ and $`HG`$, the Cartan mapping $`GG/H`$, $`gU(g)`$ is as shown in the table below:

$$\begin{array}{cccc}& & & \\ \text{class}\hfill & G\hfill & H\hfill & U(g)\hfill \\ & & & \\ \mathrm{A}(\mathrm{CUE})\hfill & \mathrm{U}(N)\hfill & 1\hfill & g\hfill \\ \mathrm{AI}(\mathrm{COE})\hfill & \mathrm{U}(N)\hfill & \mathrm{O}(N)\hfill & g^Tg\hfill \\ \mathrm{AII}(\mathrm{CSE})\hfill & \mathrm{U}(2N)\hfill & \mathrm{Sp}(2N)\hfill & g^TJg\hfill \\ \mathrm{AIII}(\mathrm{chCUE})\hfill & \mathrm{U}(N+N^{})\hfill & \mathrm{U}(N)\times \mathrm{U}(N^{})\hfill & Ig^{}Ig\hfill \\ \mathrm{BDI}(\mathrm{chCOE})\hfill & \mathrm{SO}(N+N^{})\hfill & \mathrm{SO}(N)\times \mathrm{SO}(N^{})\hfill & Ig^TIg\hfill \\ \mathrm{CII}(\mathrm{chCSE})\hfill & \mathrm{Sp}(2(N+N^{}))\hfill & \mathrm{Sp}(2N)\times \mathrm{Sp}(2N^{})\hfill & Ig^DIg\hfill \\ \mathrm{D},\mathrm{B}\hfill & \mathrm{SO}(2N),\mathrm{SO}(2N+1)\hfill & 1\hfill & g\hfill \\ \mathrm{C}\hfill & \mathrm{Sp}(2N)\hfill & 1\hfill & g\hfill \\ \mathrm{CI}\hfill & \mathrm{Sp}(4N)\hfill & \mathrm{U}(2N)\hfill & Ig^{}Ig\hfill \\ \mathrm{DIII}_{e,o}\hfill & \mathrm{SO}(4N),\mathrm{SO}(4N+2)\hfill & \mathrm{U}(2N),\mathrm{U}(2N+1)\hfill & g^Dg\hfill \end{array}$$
(1)
Here

$$J=\left[\begin{array}{cc}0& 𝟙_{}\\ 𝟙_{}& 0\end{array}\right],I=\left[\begin{array}{cc}𝟙_{}& 0\\ 0& 𝟙_{^{}}\end{array}\right]$$
(2)
($`N2N,N^{}2N^{}`$ for CII and $`N,N^{}2N`$ for CI), and $`g^D=Jg^TJ^1`$ denotes the quaternion dual of the matrix $`g`$. Out of these twelve classes, last non–classical nine possess chirality Zir1 , i.e. non–zero eigenphases appear in complex conjugate pairs. The JPD of these non–zero eigenphases

$$U=V\mathrm{diag}(\mathrm{e}^{i\theta _1},\mathrm{},\mathrm{e}^{i\theta _N},\mathrm{e}^{i\theta _1},\mathrm{},\mathrm{e}^{i\theta _N},1,\mathrm{},1)V^{},$$
(3)
for the circular ensemble is given by ForJE  $`(0\theta \pi )`$

$`P(\theta _1,\mathrm{},\theta _N)d\theta _1\mathrm{}d\theta _N`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left(d\theta _i\mathrm{sin}^{c_1}{\displaystyle \frac{\theta _i}{2}}\mathrm{cos}^{c_2}{\displaystyle \frac{\theta _i}{2}}\right)|\mathrm{\Delta }_N(\mathrm{cos}\theta )|^\beta `$ (4)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left(dy_iy_i^{(c_11)/2}(1y_i)^{(c_21)/2}\right)|\mathrm{\Delta }_N(y)|^\beta ,y_i\mathrm{sin}^2(\theta _i/2),`$ (5)
where $`\mathrm{\Delta }_N`$ is the Vandermonde determinant of the rank $`N`$ and the constants $`c_1,c_2`$ and $`\beta `$ for all nine chiral ensembles are given in the table below.

$$\begin{array}{ccccccc}& & & & & & \\ \text{class}\hfill & c_1& c_2& \beta & \nu & 𝒟\text{ (circular RME)}\hfill & \text{ (F-replica NL}\sigma \text{M)}\hfill \\ & & & & & & \\ \text{AIII}\hfill & 2(N^{}N)+1& 1& 2& N^{}N& \mathrm{U}(N+N^{})/\left(\mathrm{U}(N)\times \mathrm{U}(N^{})\right)\hfill & \mathrm{U}(n)\hfill \\ \text{BDI}\hfill & N^{}N& 0& 1& N^{}N& \mathrm{SO}(N+N^{})/\left(\mathrm{SO}(N)\times \mathrm{SO}(N^{})\right)\hfill & \mathrm{U}(2n)/\mathrm{Sp}(2n)\hfill \\ \text{CII}\hfill & 4(N^{}N)+3& 3& 4& N^{}N& \mathrm{Sp}(2(N+N^{}))/\left(\mathrm{Sp}(2N)\times \mathrm{Sp}(2N^{})\right)\hfill & \mathrm{U}(2n)/\mathrm{O}(2n)\hfill \\ \text{D}\hfill & 0& 0& 2& 1/2& \mathrm{SO}(2N)\hfill & \mathrm{SO}(2n)/\mathrm{U}(n)\hfill \\ \text{B}\hfill & 2& 0& 2& 1/2& \mathrm{SO}(2N+1)\hfill & \mathrm{SO}(2n)/\mathrm{U}(n)\hfill \\ \text{C}\hfill & 2& 2& 2& 1/2& \mathrm{Sp}(2N)\hfill & \mathrm{Sp}(2n)/\mathrm{U}(n)\hfill \\ \text{CI}\hfill & 1& 1& 1& 1& \mathrm{Sp}(2N)/\mathrm{U}(N)\hfill & \mathrm{Sp}(2n)\hfill \\ \text{DIII}_e\hfill & 1& 1& 4& 1/2& \mathrm{SO}(4N)/\mathrm{U}(2N)\hfill & \mathrm{SO}(2n)\hfill \\ \text{DIII}_o\hfill & 5& 1& 4& 1/2& \mathrm{SO}(4N+2)/\mathrm{U}(2N+1)\hfill & \mathrm{SO}(2n)\hfill \end{array}$$
(6)
The JPD Eq. (5) defined over $`y[0,1]`$ may be comprehensively called Jacobi ensemble. The constant $`\nu `$, defined as

$$\nu =\frac{c_1+1}{\beta }1,$$
(7)
may be called “topological charge”, when the tangent element of a random matrix $`U`$ in the first three cases is interpreted as the QCD Dirac operator in even dimensions SV ; Ver .

The real characteristic polynomial, or so called the fermionic replicated partition function, for these circular ensembles is defined as

$$Z_{n,N}(\theta )=\underset{𝒟}{}𝑑Udet(\mathrm{e}^{i\frac{\theta }{2}}\mathrm{e}^{i\frac{\theta }{2}}U)^n\underset{G}{}𝑑gdet(\mathrm{e}^{i\frac{\theta }{2}}\mathrm{e}^{i\frac{\theta }{2}}U(g))^n.$$
(8)
After the color-flavor transformation Zir1  and the thermodynamic limit where $`N\mathrm{},\theta 0`$ with their product (denoted by the same $`\theta `$ for the sake of simplicity) fixed finite, it takes the form,

$$Z_n(\theta )\underset{N\mathrm{}}{lim}Z_{n,N}\left(\frac{\theta }{N}\right)=\theta ^{n\nu }\underset{}{}𝑑u(detu)^\nu \mathrm{e}^{i\frac{\theta }{2}\mathrm{tr}(u+u^{})},$$
(9)
where $``$ are the ‘dual’ symmetric spaces of the fermionic nonlinear $`\sigma `$–models which are listed in the table (6). Hereafter we suppress irrelevant normalization constants, that goes to unity in the replica limit, $`n0`$. The derivation of Eq. (9) from Eq. (8) for the AIII, BDI, and CII classes is summarized in Appendix A. In order to derive Eq. (9) for symmetry classes whose pertinent color–flavor transformations are not immediately available in the literature, one could as well adopt an alternative method of magnifying the origin of the circular ensembles first (i.e. employing Gaussian ensembles), performing Hubbard-Stratonovich transformation and taking the thermodynamic limit. For a Gaussian ensemble treatment of the BDI and CII classes, see Ref. HV .

Performing the integration over angular degrees of freedom $`v`$ of $`u=v\mathrm{diag}(\mathrm{e}^{i\varphi _a})v^{}`$, one obtains

$$Z_n(\theta )=\theta ^{n\nu }\underset{0}{\overset{2\pi }{}}\underset{a=1}{\overset{n}{}}\left(d\varphi _a\mathrm{e}^{i\nu \varphi _a+i\theta \mathrm{cos}\varphi _a}\right)\left|\underset{a>b}{\overset{n}{}}\mathrm{sin}\left(\frac{\varphi _a\varphi _b}{2}\right)\right|^{4/\beta }.$$
(10)
This expression is valid for all nine chiral symmetry classes. Note that Eq. (10) depends only on $`\beta `$ and $`c_1`$ but not on $`c_2`$, because we have magnified the vicinity of the origin, $`\theta =0`$. Accordingly the class C gives the same $`Z_n(\theta )`$ as B, and the class CI the same as BDI at $`c_1=1`$, reducing nine symmetry classes of chiral RMEs to seven universality classes. The number could be further reduced by introducing three universality classes of Laguerre ensembles, having continuous $`\nu `$ and discrete $`\beta =2,1,4`$.

Consequently we succeeded in expressing the integration over the initial $`N`$–variable JPD, Eq. (4), through the $`n`$–fold integral, Eq. (10). We shall discuss its evaluation, analytical continuation and the replica limit after we introduce the BC–CSM. This replica treatment was previously performed for the AIII class in Ref. DV .

## III BC Calogero-Sutherland Model

The generalized Calogero-Sutherland Hamiltonian OP  is defined as

$$H=\underset{i=1}{\overset{N}{}}\frac{^2}{\theta _i^2}+\underset{\alpha \mathrm{\Delta }}{}\frac{g_{|\alpha |}}{\mathrm{sin}^2(\theta \alpha /2)}.$$
(11)
Here $`\mathrm{\Delta }`$ is a root system of a Lie algebra in $`N`$–dimensional vector space, $`\theta `$ is the vector $`(\theta _1,\mathrm{},\theta _N)`$, and $`g_{|\alpha |}`$ is a coupling constant depending only on the root length. Quantum integrability is ensured by these conditions. The ordinary, translationally invariant model corresponds to the $`\mathrm{A}_{N1}`$ root system.

The quantum one–dimensional model of $`N`$ interacting particles on a semicircle, $`0\theta _i\pi `$, where $`i=1,\mathrm{},N`$ with the Hamiltonian

$$H=\underset{i=1}{\overset{N}{}}\frac{^2}{\theta _i^2}+\frac{\lambda (\lambda 1)}{2}\underset{i>j}{\overset{N}{}}\left[\frac{1}{\mathrm{sin}^2\frac{\theta _i\theta _j}{2}}+\frac{1}{\mathrm{sin}^2\frac{\theta _i+\theta _j}{2}}\right]+\frac{\lambda _1(\lambda _11)}{4}\underset{i=1}{\overset{N}{}}\frac{1}{\mathrm{sin}^2\frac{\theta _i}{2}}+\frac{\lambda _2(\lambda _21)}{4}\underset{i=1}{\overset{N}{}}\frac{1}{\mathrm{cos}^2\frac{\theta _i}{2}}$$
(12)
corresponds to the $`\mathrm{BC}_N`$ root system, hence called the BC–CSM. As the $`\mathrm{BC}_N`$ root system contains roots of length $`1,\sqrt{2},2`$, there are three independent coupling constants, $`\lambda ,\lambda _1,\lambda _2`$. This family contains CSMs corresponding to $`\mathrm{B}_N`$ ($`\lambda _2=0`$), $`\mathrm{C}_N`$ ($`\lambda _1=\lambda _2`$), and $`\mathrm{D}_N`$ ($`\lambda _1=\lambda _2=0`$) root systems as its subfamilies. In addition to the pairwise inverse–square interactions, the particles interact with their own mirror images (occupying the other semicircle, $`\pi \theta 2\pi `$) and with the two single–particle impurity potentials placed at $`\theta =0`$ and $`\theta =\pi `$. In what follows we shall assume the thermodynamic limit, $`N\mathrm{}`$, and shall focus on the vicinity of the impurity at $`\theta =0`$. The other impurity may be treated in exactly the same manner. The model is known to have the following ground–state energy OP

$$H\mathrm{\Psi }_0=E_0\mathrm{\Psi }_0,E_0=\underset{i=1}{\overset{N}{}}\left((Ni)\lambda +\frac{\lambda _1+\lambda _2}{2}\right)^2,$$
(13)
and the ground–state wave function

$$\mathrm{\Psi }_0(\theta _1,\mathrm{},\theta _N)=\underset{i=1}{\overset{N}{}}\left(\mathrm{sin}^{\lambda _1}\frac{\theta _i}{2}\mathrm{cos}^{\lambda _2}\frac{\theta _i}{2}\right)\mathrm{\Delta }_N(\mathrm{cos}\theta )^\lambda .$$
(14)
The absolute square of the ground–state wave function coincides with the JPD of the chiral RMEs (4) YKY , through a change of the coupling constants

$$\lambda =\beta /2,\lambda _1=c_1/2,\lambda _2=c_2/2.$$
(15)
Notice, that Eq. (14) is not restricted to the special values of $`\beta ,c_1,c_2`$ listed in the table (6). The energy and wave functions of the excited states were studied in Refs. Yam ; BPS ; Ser ; KO .

The particle density is defined as,

$$\rho (\theta )\underset{j=1}{\overset{N}{}}\delta (\theta \theta _j)=\sqrt{y(1y)}\underset{j=1}{\overset{N}{}}\delta (yy_j),$$
(16)
where $`y\mathrm{sin}^2(\theta /2)`$. The angular brackets denote ground–state expectation values, or equivalently averaging over the normalized JPD, Eqs. (4) and (5), for the first and second equalities correspondingly. One may then employ the replica trick to write

$$\underset{j=1}{\overset{N}{}}\delta (yy_j)=\underset{n0}{lim}\frac{1}{n\pi }\mathrm{}m\frac{d}{dy}\underset{j=1}{\overset{N}{}}(yy_jiϵ)^n.$$
(17)
As a result, one obtains

$$\rho (\theta )=\sqrt{y(1y)}\underset{n0}{lim}\frac{1}{n\pi }\mathrm{}m\frac{d}{dy}Z_{n,N}(yiϵ)|_{y=\mathrm{sin}^2(\theta /2)},$$
(18)
where the “replicated partition function” is defined as

$$Z_{n,N}(y)=\underset{0}{\overset{1}{}}\underset{i=1}{\overset{N}{}}\left(dy_iy_i^{\lambda _11/2}(1y_i)^{\lambda _21/2}(yy_i)^n\right)\left(\mathrm{\Delta }_N(y)^2\right)^\lambda .$$
(19)
Baker and Forrester BF  have noticed the integral equality due to Kaneko Kan  and Yan Yan , which we suggestively call the “generalized color-flavor transformation”. With its help one may express the partition function in the following way

$$Z_{n,N}(y)=\underset{𝒞}{}\underset{a=1}{\overset{n}{}}\left(dx_ax_a^{\frac{\lambda _1+\lambda _2+1}{\lambda }2}(1x_a)^{\frac{\lambda _2+n1/2}{\lambda }}\left[\frac{x_a(1yx_a)}{1x_a}\right]^N\right)\left(\mathrm{\Delta }_n(x)^2\right)^{1/\lambda }.$$
(20)
The integration contour $`𝒞`$ encircles the cut between $`x_a=0`$ and $`x_a=1`$. The general form of the integral identity is given in Appendix B. So far no approximation has been made. Now we pass to the thermodynamic limit, $`N\mathrm{}`$, and magnify the vicinity of the $`\theta =0`$ impurity. To this end we rescale the variable as $`\theta \theta /N`$ and correspondingly $`y\theta ^2/(4N^2)`$. By redefining the integration variables as $`x_a=12iN\theta ^1\mathrm{e}^{i\varphi _a}`$ and taking the thermodynamic limit, $`N\mathrm{}`$, one finds (see Appendix B for details)

$$Z_n(\theta )\underset{N\mathrm{}}{lim}Z_{n,N}\left(\frac{\theta }{N}\right)=\theta ^{n(\frac{\lambda _1}{\lambda }+\frac{1}{2\lambda }1)}\underset{0}{\overset{2\pi }{}}\underset{a=1}{\overset{n}{}}\left(d\varphi _a\mathrm{e}^{i(\frac{\lambda _1}{\lambda }+\frac{1}{2\lambda }1)\varphi _a+i\theta \mathrm{cos}\varphi _a}\right)\left[\underset{a>b}{\overset{n}{}}\mathrm{sin}^2\left(\frac{\varphi _a\varphi _b}{2}\right)\right]^{1/\lambda }.$$
(21)
One may introduce notation

$$\nu =\frac{\lambda _1}{\lambda }+\frac{1}{2\lambda }1,$$
(22)
to notice the exact coincidence with the $`\sigma `$–model representation of the chiral RMEs, Eq. (10), provided that the coupling constants are related via Eq. (15). The important difference is that the BC–CSM representation in the form of Eq. (21) is not restricted to the RMEs values $`\lambda =1/2,1,2`$ and special values of the topological charge, $`\nu `$.

## IV Analytic Continuation and Replica Limit

Consider the $`n`$–fold integral

$$Z_n(\theta )=\theta ^{n\nu }\underset{0}{\overset{2\pi }{}}\underset{a=1}{\overset{n}{}}\left(d\varphi _a\mathrm{e}^{i\nu \varphi _a+i\theta \mathrm{cos}\varphi _a}\right)\left[\underset{a>b}{\overset{n}{}}\mathrm{sin}^2\left(\frac{\varphi _a\varphi _b}{2}\right)\right]^{1/\lambda }.$$
(23)
One may stretch the integration contour from the unit circle in the complex plane of $`z_a=\mathrm{e}^{i\varphi _a}`$ into two lines parallel to the imaginary axis with $`\mathrm{}ez_a=\pm 1`$. The original integral, Eq. (23), splits into the sum of $`n`$ terms with $`l`$ integrals having $`\mathrm{}ez_a=1`$ and remaining $`nl`$ ones $`\mathrm{}ez_a=1`$; here $`l=1,\mathrm{},n`$. The further progress is made possible in the asymptotic limit, $`\theta 1`$. In this case the integrals are dominated by the vicinities of the saddle points $`z_a=\pm 1`$ and thus may be evaluated employing the Selberg integral. This strategy was described in details in Refs. KM1 ; YL ; NK ; GK . Proceeding this way, one finds for the replicated partition function

$$Z_n(\theta )=\theta ^{n\nu }\mathrm{e}^{in\theta }\underset{l=0}{\overset{n}{}}F_n^l(\lambda )\mathrm{e}^{i\pi \nu l2il\theta }\mathrm{\hspace{0.17em}2}^{\frac{2l^2}{\lambda }}(i\theta )^{\frac{nl}{2}\left(1+\frac{nl1}{\lambda }\right)}(i\theta )^{\frac{l}{2}\left(1+\frac{l1}{\lambda }\right)},$$
(24)
where we have omitted a normalization constant that goes to unity in the replica limit, $`n0`$ and, following Refs. KM1 ; YL ; NK ; GK , introduced the notation

$$F_n^l(\lambda )\left(\genfrac{}{}{0pt}{}{n}{l}\right)\underset{a=1}{\overset{l}{}}\frac{\mathrm{\Gamma }(1+a/\lambda )}{\mathrm{\Gamma }(1+(na+1)/\lambda )}.$$
(25)
Employing the observation that $`F_n^l(\lambda )0`$ for $`l>n`$, one may extend summation over $`l`$ in Eq. (24) to infinity and then perform the analytic continuation $`n0`$. As a result, one finds via $`\rho (\theta )=lim_{n0}(\pi n)^1\mathrm{}m_\theta Z_n(\theta )`$ (cf. Eq. (18)),

$$\rho (\theta )=\frac{1}{\pi }\left[1+2\underset{l=1}{\overset{\mathrm{}}{}}\frac{d_l(\lambda )}{(2\theta )^{l^2/\lambda }}\mathrm{cos}\left(2l\theta l\pi \left(\nu +\frac{1}{2}\frac{1}{2\lambda }\right)\right)\right],$$
(26)
where

$$d_l(\lambda )\frac{(1)^l}{2^{l^2/\lambda }}\underset{a=1}{\overset{l}{}}\frac{\mathrm{\Gamma }(1+a/\lambda )}{\mathrm{\Gamma }(1(a1)/\lambda )}.$$
(27)
Eq. (26) for the asymptotic of the ground–state density of the BC–CSM is the central result of this paper. For the integer values of the coupling constant, $`\lambda `$, it may be extracted from the expressions derived by Baker and Forrester BF , see also Ref. FM .

We remark that each term in Eq. (26) is the one with the lowest power in $`\theta ^1`$ among all terms carrying the same frequency $`2l`$. This corresponds to truncating all contributions from the descendent fields in the conformal description, see the next section. One could compute these secondary terms, subleading in $`\theta `$, by performing a perturbative expansion around each saddle point with the help of the loop equations NK .

## V Conformal Field Theory Description

By comparing the one–particle correlation function of the BC–CSM (26) with the equal–time two–particle correlation function of the ordinary A–CSM Ha ; GK

$$\rho (\theta )\rho (0)=\frac{1}{(2\pi )^2}\left[1\frac{1}{2\lambda \theta ^2}+2\underset{l=1}{\overset{\mathrm{}}{}}\frac{d_l(\lambda )^2}{\theta ^{2l^2/\lambda }}\mathrm{cos}(2l\theta )\right],$$
(28)
one notices that it is expressed through the very same coefficients $`d_l(\lambda )`$, Eq. (27). It is thus clear that the harmonics amplitudes, $`d_l(\lambda )`$, are properties of the homogenous interacting system and not of the localized impurity. As explained below, this result could be anticipated from the effective conformal field theory description. The latter is capable of predicting the low energy properties of the system apart from numerical values of the coefficients.

Based on the finite–size scaling analysis, Kawakami and Yang KY  identified the low–energy effective theory of the CSM in the thermodynamic limit to be the $`c=1`$ Gaussian conformal field theory at radius $`R=\sqrt{\lambda }/2`$, either non–chiral (A–CSM)

$$=\frac{1}{2\pi }_z\mathrm{\Phi }_{\overline{z}}\mathrm{\Phi },\mathrm{\Phi }(z,\overline{z})=\varphi (z)+\overline{\varphi }(\overline{z}),\mathrm{\Phi }\mathrm{\Phi }+2\pi R,$$
(29)
or chiral (BC–CSM) YKY ; YKY2 . Namely, they have found the density operator should have zero winding number. It does not have a definite conformal dimension, and therefore is expanded in terms of primary and secondary operators whose left– and right–moving vertex momenta are equal. Here we shall consider only contributions from the $`\mathrm{U}(1)`$ current and the primary fields (vertex operators with charges allowed by the compactification (29)),

$$\rho (z,\overline{z})=\rho _0\left[b\left(_z\varphi (z)+_{\overline{z}}\overline{\varphi }(\overline{z})\right)+\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}d_l\mathrm{e}^{il(z+\overline{z})}\mathrm{e}^{il\varphi (z)/R}\mathrm{e}^{il\overline{\varphi }(\overline{z})/R}\right],$$
(30)
where the expansion coefficients, $`b`$ and $`d_l(=d_l)`$, are not determined from the conformal field theory. Neither are the oscillation factors $`\mathrm{e}^{il(z+\overline{z})}`$, which are set by hand to describe transport of $`l`$ pseudo–particles Ha  from the left Fermi point ($`k_F=1`$ by normalization) to the right Fermi point ($`k_F=1`$) KY . We have factored out $`\rho _0`$ such that the constant term carries $`d_0=1`$. The propagator and the vertex correlator are given by

$$\varphi (z)\varphi (z^{})=\frac{1}{4}\mathrm{log}(zz^{});\mathrm{e}^{il\varphi (z)/R}\mathrm{e}^{il^{}\varphi (z^{})/R}=\frac{\delta _{l,l^{}}}{(zz^{})^{l^2/(4R^2)}}.$$
(31)
The coefficient in the second equation is a matter of convention and reflects a particular choice of the ultraviolet regularization. Another choice of the regularization would change coefficients $`d_l`$, but not the final result. Employing Eqs. (30), (31), one obtains for the equal–time two–particle correlation function of the A–CSM (we denote $`z=\theta +i\tau `$)

$$\rho (\theta +i0)\rho (\theta ^{}+i0)=\rho _0^2\left[\frac{b^2}{2(\theta \theta ^{})^2}+\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}d_l^2\frac{\mathrm{e}^{2il(\theta \theta ^{})}}{(\theta \theta ^{})^{l^2/(2R^2)}}\right].$$
(32)
Comparing this expression with Eq. (28), one finds that the $`b=1/\sqrt{\lambda }`$ Ha , while $`d_l=d_l(\lambda )`$, cf. Eq. (27). Let us consider now the BC–CSM and concentrate on the case without the phase shift for simplicity. The Dirichlet boundary condition at $`\mathrm{}ez=0`$ is translated into the open boundary bosonization rule GNT

$$\varphi (z)=\overline{\varphi }(z)\text{for}\mathrm{}mz<0.$$
(33)
As the right mover is identified as $`(1)`$ times the left mover at the mirror–imaged point, there exists a nonzero matrix element between the left– and right–moving vertex operators:

$$\mathrm{e}^{il\varphi (z)/R}\mathrm{e}^{il^{}\overline{\varphi }(\overline{z}^{})/R}=\frac{\delta _{l,l^{}}}{(z+\overline{z}^{})^{l^2/(4R^2)}}.$$
(34)
As a result, the mean density of the BC–CSM becomes nontrivial,

$$\rho (\theta +i0)=\rho _0\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}d_l\frac{\mathrm{e}^{2il\theta }}{(2\theta )^{l^2/(4R^2)}}.$$
(35)
It is in exact agreement with Eq. (26) if one disregards the phase shift. The phase shift may also be included in the conformal description by shifting the identification of the right and left movers, Eq. (33), by a constant factor: $`\pi \left(\nu +1/21/(2\lambda )\right)`$.

We found that the knowledge of the asymptotic behavior of the correlator of the homogenous A–CSM, Eq. (28), supplemented by the conformal field theory description is, in principle, sufficient to predict the BC–CSM correlation function, Eq. (26). This agreement indicates the fact that the ultraviolet property of the field that is responsible for the normalization of vertex operators are not affected by the presence or absence of the boundary. Well anticipated as it is, we nevertheless consider this fact to be worth verifying, as done in this paper. This fact can be put on a further test by computing the asymptotics of e.g. two–particle correlation function for the BC–CSM, though it is technically more challenging. We could as well reverse the logic and conjecture the asymptotically expanded form of any $`p`$-point correlation function of density operators for the A–CSM ($`p3`$) or for the BC–CSM ($`p2`$), by using Eqs. (30), (31), and (27).

## VI Nearest Particle Distribution

Another quantity of interest in the theory of interacting electrons is the probability $`E[s^{},s]`$ of finding no particle within an interval $`[s^{},s]`$ or the particle spacing distribution $`p(s)`$ that is a derivative of the former. In the context of spin chains a similar quantity was recently discussed in Ref. Abanov . For nonchiral as well as chiral RMEs, Tracy and Widom TW  has developed Mehta’s computation Meh  of $`E[s^{},s]`$ as a Fredholm determinant into a systematic and powerful method. As their method determines $`E[s^{},s]`$ as a solution ($`\tau `$ function) to a transcendental equation of Painlevé type relies upon the orthogonal polynomials, its validity is necessarily limited to $`\lambda =1/2,1,2`$. On the other hand, $`E(s)E[0,s]`$ for the chiral RMEs has been computed by an alternative and far simpler “shifting” method SN ; ForNPB ; ForJMP ; FH ; DN  as explained below. We show that with a help of the generalized color–flavor transformation, this method is applicable also to the BC–CSM at generic values of the coupling constants. Consider for a moment $`\lambda _1=n+1/2`$, where $`n=0,1,2,\mathrm{}`$ is an integer and $`\lambda _2=1/2`$. As we are interested in the universal behavior in the vicinity of the reflection point $`\theta =0`$, the restriction on $`\lambda _2`$ is irrelevant. The probability of having no particle within an interval $`0\theta s`$, or $`0yY`$ with $`Y=\mathrm{sin}^2(s/2)`$, is defined as

$$E_N(s)=\mathrm{const}.\underset{Y}{\overset{1}{}}\underset{i=1}{\overset{N}{}}\left(dy_iy_i^n\right)\left(\mathrm{\Delta }_N(y)^2\right)^\lambda .$$
(36)
The constant should be chosen to ensure $`E_N(0)=1`$. By shifting and rescaling the integration variable as $`y(1Y)y+Y`$, one obtains

$$E_N(s)=\mathrm{const}.(1Y)^{(1+n)N+\lambda N(N1)}\underset{0}{\overset{1}{}}\underset{i=1}{\overset{N}{}}\left[dy_i\left(y_i+\frac{Y}{1Y}\right)^n\right]\left(\mathrm{\Delta }_N(y)^2\right)^\lambda .$$
(37)
Now we apply the generalized color–flavor transformation (20) to reexpress $`E_N(s)`$ as an $`n`$–fold integral,

$$E_N(s)=\mathrm{const}.(1Y)^{(1+n)N+\lambda N(N1)}\underset{𝒞}{}\underset{a=1}{\overset{n}{}}\left[dx_ax_a^{\frac{2}{\lambda }2}(1x_a)^{\frac{n}{\lambda }}\left[\frac{x_a(1+\frac{Y}{1Y}x_a)}{1x_a}\right]^N\right](\mathrm{\Delta }_n(x)^2)^{1/\lambda }.$$
(38)
We finally rescale $`ss/N;Ys^2/(4N^2)`$ and take the thermodynamic limit. Following the same procedure that lead from Eq. (20) to Eq. (21), one finds

$`E(s)`$ $``$ $`\underset{N\mathrm{}}{lim}E_N\left({\displaystyle \frac{s}{N}}\right)`$ (39)
$`=`$ $`\mathrm{const}.\mathrm{e}^{\frac{\lambda }{4}s^2}s^{n(1\frac{1}{\lambda })}{\displaystyle \underset{0}{\overset{2\pi }{}}}{\displaystyle \underset{a=1}{\overset{n}{}}}\left(d\varphi _a\mathrm{e}^{i(\frac{1}{\lambda }1)\varphi _a+s\mathrm{cos}\varphi _a}\right)\left[{\displaystyle \underset{a>b}{\overset{n}{}}}\mathrm{sin}^2\left({\displaystyle \frac{\varphi _a\varphi _b}{2}}\right)\right]^{1/\lambda }.`$
This exact result has previously been derived from Laguerre (chiral Gaussian) ensembles at arbirtary $`\lambda `$ and at $`\lambda _1=n+1/2`$ ForJMP .

To compute the asymptotics of $`E(s)`$ for $`s1`$, one can relax the restriction on $`\lambda _1`$, by first evaluating the $`n`$–fold integral (39) by the saddle point method and then performing the analytical continuation $`n\lambda _11/2`$. In the large–$`s`$ limit we pick only the contribution of the replica-symmetric saddle point $`\varphi _a=0`$, since the contribution of all other saddle points is exponentially smaller (the difference with the previous computation is that $`s`$ enters Eq. (39) without imaginary unit). This way Forrester ForJMP  has derived the asymptotic $`(s1)`$ result

$`E(s)`$ $`=`$ $`\mathrm{const}.s^{\frac{n}{2}\frac{n(n+1)}{2\lambda }}\mathrm{e}^{\frac{\lambda }{4}s^2}\mathrm{e}^{ns}`$ (40)
$`=`$ $`\mathrm{const}.s^{\nu \left(\lambda _1\frac{1}{2}\right)}\mathrm{e}^{\frac{\lambda }{4}s^2+\left(\lambda _1\frac{1}{2}\right)s},`$
where $`\nu `$ is the topological charge defined by Eq. (22). The Gaussian factor $`\mathrm{exp}(\lambda s^2/4)`$ could be anticipated from the mean–field treatment of the classical logarithmic gas Abanov . The other factors in the asymptotic expression Eq. (40) could not be found in any simpler way, to the best of our knowledge.

The distribution $`p(s)`$ of the locus $`s`$ of the particle nearest to the reflection point is given by $`p(s)=_sE(s)`$. In the other limiting case, $`s1`$, $`p(s)`$ is determined by the interaction of the particle closest to the reflection point to its own mirror image. Inspecting Eq. (5) one immediately finds $`p(s)s^{2\lambda _1}`$.

## VII Discussions

Let us now take a closer look at our main result, Eq. (26). The constant term on its r.h.s., $`\rho _0=1/\pi `$, represents the uniform density of particles ($`N`$ particles within $`[0,\pi ]`$) far away from the impurity. It may be traced back to the replica symmetric contribution (all $`n`$ integrals are taken at $`z_a=1`$ saddle point) to the partition function. The decaying (as $`\theta \mathrm{}`$) oscillatory terms on the r.h.s. of Eq. (26) are the Friedel oscillations of the particle density induced by the mirror boundary and the impurity potential. These terms may be identified as the replica symmetry broken contributions to the partition functions ($`l`$ integrals are taken at the “wrong” saddle point, $`z_a=1`$). In general, there is an infinite number of harmonics (unlike a single “$`2k_F`$” harmonic in the non–interacting system!) in the oscillation spectrum. Since $`d_l\mathrm{exp}(\lambda ^1l^2\mathrm{ln}l)`$ for $`l1`$, the sum over harmonics on the r.h.s. of Eq. (26) is, in general, not convergent. It is not clear to us at the moment, whether there is a consistent resummation scheme. There is, however, an important class of the parameters, where Eq. (26) is mathematically rigorous.

For any rational coupling constant $`\lambda =p/q`$ the coefficient $`d_l(p/q)0`$ for $`l>p`$ and therefore the sum terminates after exactly $`p`$ oscillatory components. One finds that, in addition to the usual $`2k_F`$ Friedel oscillation, the system possess $`4k_F,\mathrm{},2pk_F`$ oscillatory components of the density (in the unit $`k_F=1`$ accepted here). This fact might be expected from the form of the density–density response function of the homogenous A–CSM Hal ; Ha . However, the algebraic decay rate of the harmonics could not be determined employing linear response of the A–CSM. Indeed, the latter predicts that the $`l`$–th harmonic decays as $`\theta ^{2l^2/\lambda +1}`$, while the correct decay rate is $`\theta ^{l^2/\lambda }`$. Notice, that for non–interacting particles, $`\lambda =1`$ and thus $`l=1`$, both ways give the correct one–dimensional decay of the $`2k_F`$ Friedel oscillations: $`\theta ^1`$. For any interacting system, $`\lambda 1`$, the linear response is bound to fail in the asymptotic regime. These observations was already made in the Tomonaga–Luttinger liquid literature KF ; EG . Now we can confirm them having the exactly solvable model system.

In the BC–CSM one finds not only the decay law, but also the relative amplitudes of the harmonics: coefficients $`d_l(\lambda )`$, Eq. (27). Notice that these amplitudes are determined by the interaction strength, $`\lambda `$, only and are independent on the impurity strength, $`\lambda _1`$. This is due to the fact that the mirror boundary condition induces oscillations of the maximal possible amplitude. The additional single–particle potential centered at $`\theta =0`$ and characterized by $`\lambda _1`$ changes the phase of the oscillations only. The entire information about the impurity strength, $`\lambda _1`$, is incorporated in the parameter $`\nu `$, Eq. (22). In the asymptotic regime the latter affects the phase of the Friedel oscillations only and therefore may be associated with the impurity phase shift. (Unlike the leading order, the amplitudes of sub–leading perturbative corrections in negative powers of $`\theta `$ do depend on the phase shift, $`\nu `$.)

To verify Eqs. (26), (27) one may compare them with the available exact DoS of the chiral RMEs (see Ref. Iva  and references therein). Employing Hankel’s asymptotic expansion of the Bessel function, one may check that these asymptotic perfectly agree with Eqs. (26), (27). One may also notice that for D, B, and C symmetry classes we have obtained the exact rather than the asymptotic results. This coincidence is due to the Duistermaat–Heckman localization theorem DH ; Zir99 . We see that having unitary symmetry class, $`\beta =2`$, is not sufficient to satisfy the Duistermaat–Heckman theorem. One should also have the special value of the topological charge, $`\nu =\pm 1/2`$, to secure cancellation of all higher–order perturbative corrections DV .

One may define the total charge attracted (expelled) by the impurity to (from) the region near $`\theta =0`$ as:

$$Q\underset{0}{\overset{\mathrm{}}{}}𝑑\theta (\rho (\theta )\rho _0),$$
(41)
where $`\rho _0=1/\pi `$ is the uniform asymptotic density. For the chiral RMEs, where the exact expressions including small $`\theta `$ region are available, the result is (cf. Eq. (26) ):

$$Q=\frac{1}{2}\left(\nu +\frac{1}{2}\frac{1}{2\lambda }\right)=\frac{1}{4}\frac{\lambda _1}{2\lambda }.$$
(42)
Notice, that the attracted charge depends both on the impurity amplitude, $`\lambda _1`$ and the interaction strength, $`\lambda `$, reflecting the fact that the impurity is screened due to the interactions. The pure mirror boundary, without the single–particle potential attracts quarter of a particle irrespective to the interaction strength. Equation (42) is a manifestation of the famous Friedel sum rule: the total expelled charge is equal to the impurity phase shift (divided by $`2\pi `$); the latter also determines the phase of the density oscillations far from the impurity. We conjecture, in accordance with the earlier works Langer , that Eq. (42) is valid for any values of $`\lambda `$ and $`\lambda _1`$.

###### Acknowledgements.

We are grateful to the hospitality of Department of Physics, Technion–Israel Institute of Technology, Haifa, Israel, where our collaborated work has been initiated. We also acknowledge A. Abanov, A. Altland, G. Dunne, A. García-García, P. Forrester, J. Verbaarschot, and M. Zirnbauer for valuable discussions and correspondences. This work was supported in part (SMN) by the DOE grant no. DE-FR02-92ER40716, and (AK) by the BSF grant no. 9800338. Le LKB est UMR 8552 du CNRS, de l’ENS et de Université P. et M. Curie.

## Appendix A $`\sigma `$–Model Derivation via Color–Flavor Transformation

We first consider the simplest case AIII with $`\nu =0`$. The Grassmannian $`\mathrm{U}(2N)/\left(\mathrm{U}(N)\times \mathrm{U}(N)\right)`$ is a complex Kähler manifold, and its unitary matrix representative $`U`$ in Table (1) is conveniently parameterized by the complex stereographic coordinate $`Z_{ij}`$, $`i,j=1,\mathrm{},N`$, as

$$U=Ig^{}Ig=I\gamma I\gamma ^1,\gamma =\left[\begin{array}{cc}1& Z\\ Z^{}& 1\end{array}\right],Z𝐂^{N\times N},I=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$
(43)
The Kähler potential $`k(Z,\overline{Z})=\mathrm{tr}\mathrm{log}(1+ZZ^{})`$ leads to the Haar measure

$$dU=\frac{_{i,j=1}^Nd^2Z_{ij}}{det(1+ZZ^{})^{2N}}.$$
(44)
The replicated partition function then reads

$`Z_{n,N}(\theta )`$ $`=`$ $`{\displaystyle _{\mathrm{U}(2N)/\left(\mathrm{U}(N)\times \mathrm{U}(N)\right)}}𝑑Udet(\mathrm{e}^{i\frac{\theta }{2}}\mathrm{e}^{i\frac{\theta }{2}}U)^n`$ (45)
$`=`$ $`{\displaystyle _{𝐂^{N\times N}}}{\displaystyle \frac{d^2Z}{det(1+ZZ^{})^{2N}}}det\left(\mathrm{e}^{i\frac{\theta }{2}}\mathrm{e}^{i\frac{\theta }{2}}I\gamma I\gamma ^1\right)^n`$
$`=`$ $`{\displaystyle _{𝐂^{N\times N}}}{\displaystyle \frac{d^2Z}{det(1+ZZ^{})^{2N+n}}}det\left(\mathrm{e}^{i\frac{\theta }{2}}\gamma \mathrm{e}^{i\frac{\theta }{2}}I\gamma I\right)^n.`$
We introduce a set of $`(N\times n)`$–component independent Grassmannian numbers $`\psi _i^a,\chi _i^a,\overline{\psi }_i^a,\overline{\chi }_i^a,i=1,\mathrm{},N,a=1,\mathrm{},n`$ to exponentiate the determinant, ($`I\gamma I=\gamma ^{}`$)

$$Z_{n,N}(\theta )=𝑑\psi 𝑑\overline{\psi }𝑑\chi 𝑑\overline{\chi }_{𝐂^{N\times N}}\frac{d^2Z}{det(1+ZZ^{})^{2N+n}}\mathrm{exp}\left([\overline{\psi }\overline{\chi }]\left(\mathrm{e}^{i\frac{\theta }{2}}\gamma \mathrm{e}^{i\frac{\theta }{2}}\gamma ^{}\right)\left[\genfrac{}{}{0pt}{}{\psi }{\chi }\right]\right).$$
(46)
Now we employ Zirnbauer’s color–flavor transformation Zir1

$$_{𝐂^{N\times N}}\frac{d^2Z}{det(1+ZZ^{})^{2N+n}}\mathrm{exp}\left(\overline{\psi }_i^aZ_{ij}\chi _j^a\overline{\chi }_i^aZ_{ij}^{}\psi _j^a\right)=_{\mathrm{U}(n)}𝑑u\mathrm{exp}\left(\overline{\psi }_i^au^{ab}\psi _i^b+\overline{\chi }_i^au^{}{}_{}{}^{ab}\chi _{i}^{b}\right),$$
(47)
to obtain

$`Z_{n,N}(\theta )`$ $`=`$ $`{\displaystyle 𝑑\psi 𝑑\overline{\psi }𝑑\chi 𝑑\overline{\chi }_{\mathrm{U}(n)}𝑑u\mathrm{exp}\left((\mathrm{e}^{i\frac{\theta }{2}}\mathrm{e}^{i\frac{\theta }{2}})(\overline{\psi }\psi +\overline{\chi }\chi )+(\mathrm{e}^{i\frac{\theta }{2}}+\mathrm{e}^{i\frac{\theta }{2}})(\overline{\psi }u\psi +\overline{\chi }u^{}\chi )\right)}`$ (48)
$`=`$ $`{\displaystyle _{\mathrm{U}(n)}}𝑑udet\left(\mathrm{cos}\theta +i\mathrm{sin}\theta {\displaystyle \frac{u+u^{}}{2}}\right)^N.`$
We stress that no approximation has been made in the above procedure.

In the thermodynamic limit, $`N\mathrm{}`$ and $`\theta 0`$, one has $`\mathrm{cos}\theta 1`$ and $`\mathrm{sin}\theta \theta `$. As a result, the determinant in Eq. (48) may be exponentiated:

$$Z_n(\theta )\underset{N\mathrm{}}{lim}Z_{n,N}\left(\frac{\theta }{N}\right)=_{\mathrm{U}(n)}𝑑u\mathrm{e}^{i\frac{\theta }{2}\mathrm{tr}(u+u^{})}.$$
(49)
One can repeat the above procedure for the BDI and CII classes. The parametrization (43) of the real and quaternionic Grassmannian manifolds, $`\mathrm{SO}(2N)/(\mathrm{SO}(N)\times \mathrm{SO}(N))`$ and $`\mathrm{Sp}(4N)/(\mathrm{Sp}(2N)\times \mathrm{Sp}(2N))`$, involves $`N\times N`$ real and quaternion–real matrices $`Z`$, instead of complex. The color–flavor transformation trades the integrations over these ‘colored’ variables $`Z`$ with the ones over ‘flavored’ variables $`u`$ that are antisymmetric and symmetric unitary matrices, respectively (straightforward as they are, such types of color–flavor transformation have yet to be exhibited explicitly in the literature, to the best of our knowledge). Accordingly, the integration domain of the transformed partition function (48) becomes antisymmetric unitary matrices ($`\mathrm{U}(2n)/\mathrm{Sp}(2n)`$) and symmetric unitary matrices ($`\mathrm{U}(2n)/\mathrm{O}(2n)`$), with the rest being unaltered. Inclusion of non-zero $`\nu `$ is straightforward by considering a rectangular $`Z`$. It merely shifts the the power of the determinant in the measure (44) by $`\nu `$, and modifies Eq. (49) into Eq. (9) in the thermodynamic limit.

For the B–D and C or DIII and CI classes, one utilizes the color–flavor transformation between the orthogonal or symplectic group and the associated dual symmetric space parametrized by antisymmetric or symmetric complex matrices, respectively Zir1 ; Zir98 ; NN .

## Appendix B Generalized Color–Flavor Transformation

Kaneko Kan  (see also Yan Yan ) has derived the following remarkable integral identity:

$`Z_{n,N}(t)`$ $`=`$ $`{\displaystyle \frac{1}{S_N(\mathrm{\Lambda }_1+n,\mathrm{\Lambda }_2,\lambda )}}{\displaystyle \underset{0}{\overset{1}{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(dy_iy_i^{\mathrm{\Lambda }_1}(1y_i)^{\mathrm{\Lambda }_2}(y_it)^n\right)\left(\mathrm{\Delta }_N(y)^2\right)^\lambda `$ (50)
$`=`$ $`{\displaystyle \frac{1}{S_n(V_1,V_2,1/\lambda )}}{\displaystyle \underset{𝒞}{}}{\displaystyle \underset{a=1}{\overset{n}{}}}\left(dx_ax_a^{V_1}(1x_a)^{V_2}(1tx_a)^N\right)\left(\mathrm{\Delta }_n(x)^2\right)^{1/\lambda },`$
where the constants are related as

$$V_1=\frac{\mathrm{\Lambda }_1+\mathrm{\Lambda }_2+2}{\lambda }+N2;V_2=\frac{\mathrm{\Lambda }_2+n}{\lambda }N.$$
(51)
The normalization constant $`S_k(a,b,c)`$ is given by the Selberg integral

$$S_k(a,b,c)=\underset{j=0}{\overset{k1}{}}\frac{\mathrm{\Gamma }(a+1+cj)\mathrm{\Gamma }(b+1+cj)\mathrm{\Gamma }(1+c(j+1))}{\mathrm{\Gamma }(a+b+2+c(k+j1))\mathrm{\Gamma }(1+c)}.$$
(52)
The integration contour $`𝒞`$ encircles the cut between $`x_a=0`$ and $`x_a=1`$. In our case (cf. Eq. (19))

$$\mathrm{\Lambda }_{1,2}=\lambda _{1,2}\frac{1}{2},V_1=\frac{\lambda _1+\lambda _2+1}{\lambda }+N2;V_2=\frac{\lambda _2+n1/2}{\lambda }N.$$
(53)
We call the identity (50) fermionic replica “generalized color–flavor transformation”. Indeed, for the RMEs values of the parameters ($`\lambda =1/2,1,2`$ and special values of $`\lambda _{1,2}`$) Eq. (50) essentially coincides with the fermionic replica version of Zirnbauer’s color–flavor transformation Zir1 ; Zir2 , after proper parametrization of the symmetric space elements and integration out of the irrelevant angles. The questions whether there is a geometrical “dual pair” interpretation of Eq. (50) for arbitrary parameters, and whether there is a supersymmetric (say for rational $`\lambda `$) or bosonic replica analogue, are currently open.

In the large–$`N`$ limit we collect all the terms having the $`N`$–th power into $`\mathrm{exp}[N_{a=1}^nS(x_a,t)]`$, where

$$S(x,t)=\mathrm{log}(1tx)\mathrm{log}x+\mathrm{log}(1x).$$
(54)
We then look for the stationary points of the “action” $`S(x)`$ given by solutions of $`_xS=0`$. A simple algebra gives for the stationary points

$$x_\pm =1\pm i\sqrt{\frac{1t}{t}}.$$
(55)
We then magnify the vicinity of $`t=0`$ by introducing $`\theta `$ as $`t=\mathrm{sin}^2(\theta /(2N))\theta ^2/(2N)^2`$ and changing the integration variable $`x_a`$ to $`\varphi _a`$ as

$$x_a=1i\sqrt{\frac{1t}{t}}\mathrm{e}^{i\varphi _a}.$$
(56)
The two saddle points (55) are at $`\varphi _a=0,\pi `$. Taking the limit $`N\mathrm{}`$, one obtains for the action $`NS(x_a,t)i\theta \mathrm{cos}\varphi _a`$. A straightforward algebra yields

$$\left(\mathrm{\Delta }_n(x)^2\right)^{1/\lambda }\theta ^{\frac{n(n1)}{\lambda }}\underset{a=1}{\overset{n}{}}\mathrm{e}^{i\frac{n1}{\lambda }\varphi _a}\left[\underset{a>b}{\overset{n}{}}\mathrm{sin}^2\left(\frac{\varphi _a\varphi _b}{2}\right)\right]^{1/\lambda }.$$
(57)
The closed contour $`𝒞`$ in the $`x_a`$–plane can be taken to be a circle, so that $`\varphi _a[0,2\pi ]`$. As a result, one obtains Eq. (21) of the main text.