text:
# Third Neighbor Correlators of Spin-1/2 Heisenberg Antiferromagnet

(February 27, 2003)

## Abstract

We exactly evaluate the third neighbor correlator $`S_j^zS_{j+3}^z`$ and all the possible non-zero correlators $`S_j^\alpha S_{j+1}^\beta S_{j+2}^\gamma S_{j+3}^\delta `$ of the spin-1/2 Heisenberg $`XXX`$ antiferromagnet in the ground state without magnetic field. All the correlators are expressed in terms of certain combinations of logarithm $`\mathrm{ln}2`$, the Riemann zeta function $`\zeta (3)`$, $`\zeta (5)`$ with rational coefficients. The results accurately coincide with the numerical ones obtained by the density-matrix renormalization group method and the numerical diagonalization.

Since Bethe’s pioneering work of the spin-1/2 Heisenberg $`XXX`$magnet Bethe ,

$$H=J\underset{j=1}{\overset{L}{}}\left(S_j^xS_{j+1}^x+S_j^yS_{j+1}^y+S_j^zS_{j+1}^z\right),$$
(1)
exact evaluation of the correlation functions has been a long standing problem in mathematical physics. Especially significant are the spin-spin correlators $`S_j^zS_{j+k}^z`$ (or equivalently $`S_j^+S_{j+k}^{}/2`$; $`S_j^\pm =S_j^x\pm iS_j^y`$), for which only the first and second neighbor ($`k=1,2`$) have been calculated so far:

$`S_j^zS_{j+1}^z`$ $`={\displaystyle \frac{1}{12}}{\displaystyle \frac{1}{3}}\mathrm{ln}20.14771572685,`$ (2)
$`S_j^zS_{j+2}^z`$ $`={\displaystyle \frac{1}{12}}{\displaystyle \frac{4}{3}}\mathrm{ln}2+{\displaystyle \frac{3}{4}}\zeta (3)0.06067976996,`$ (3)
where $`\zeta (s)`$ is the Riemann zeta function and $`\mathrm{}`$ denotes the ground state expectation value of the antiferromagnetic model ($`J>0`$). Here we have taken the thermodynamic limit $`L\mathrm{}`$. In this letter, we would like to report our new results regarding the third neighbor correlators. Our main result is

$`S_j^z`$ $`S_{j+3}^z={\displaystyle \frac{1}{2}}S_j^+S_{j+3}^{}`$
$`={\displaystyle \frac{1}{12}}3\mathrm{ln}2+{\displaystyle \frac{37}{6}}\zeta (3){\displaystyle \frac{14}{3}}\zeta (3)\mathrm{ln}2{\displaystyle \frac{3}{2}}\zeta (3)^2`$
$`{\displaystyle \frac{125}{24}}\zeta (5)+{\displaystyle \frac{25}{3}}\zeta (5)\mathrm{ln}20.05024862726.`$ (4)
In addition, we obtain the third neighbor one-particle Green function $`c_j^{}c_{j+3}_\mathrm{f}`$

$`c_j^{}c_{j+3}_\mathrm{f}={\displaystyle \frac{1}{30}}2\mathrm{ln}2+{\displaystyle \frac{169}{30}}\zeta (3){\displaystyle \frac{10}{3}}\zeta (3)\mathrm{ln}2`$
$`{\displaystyle \frac{6}{5}}\zeta (3)^2{\displaystyle \frac{65}{12}}\zeta (5)+{\displaystyle \frac{20}{3}}\zeta (5)\mathrm{ln}20.08228771669,`$ (5)
for the isotropic spinless fermion model corresponding to (1) by the Jordan-Wigner transformation:

$$S_k^{}=\underset{j=1}{\overset{k1}{}}(12c_j^{}c_j)c_k^{},S_k^+=\underset{j=1}{\overset{k1}{}}(12c_j^{}c_j)c_k.$$
(6)
Here $`\mathrm{}_\mathrm{f}`$ denotes the expectation value in the half-filled state of the spinless fermion model. Moreover we exactly calculate all the possible non-zero correlators $`S_j^\alpha S_{j+1}^\beta S_{j+2}^\gamma S_{j+3}^\delta `$.

The result (2) comes from the ground state energy of (1) derived by Hulthén in 1938 hult . The result (3) was obtained by one of the authors in 1977 Taka ; Takabook  via the strong coupling expansion for the ground state energy of the half-filled Hubbard model. This result is also reproduced in the framework of the asymptotic Bethe ansatz for an integrable spin chain with variable range exchange DitIno . However, probably due to the complexity of the wave function for these models, no one has succeeded in generalizing the method to obtain the higher neighbor correlators.

On the other hand, utilizing the representation theory of the quantum affine algebra $`U_q(\widehat{sl_2})`$ and the associated vertex operators, in 1992, Jimbo et al. derived a universal multiple integral representation of arbitrary correlators for the massive $`XXZ`$ antiferromagnet JMMN ; JMbook . Their result has been extended to the $`XXX`$Naka ; KIEU , the massless $`XXZ`$ JM ; KMT  and the $`XYZ`$ Quano  antiferromagnets. However the explicit evaluation, even for the second neighbor correlator (3), was not achieved for a long time.

In this respect, it is quite remarkable that Boos and Korepin recently devised a general method to evaluate the multiple integral representation especially in the study of the Emptiness Formation Probability (EFP) for the $`XXX`$antiferromagnet BK1 ; BK2 . The EFP, $`P(n)`$ describes the probability of finding a ferromagnetic string of length $`n`$ in the antiferromagnetic ground state KIEU . Explicitly it reads

$$P(n)=\underset{j=1}{\overset{n}{}}\left(S_j^z+\frac{1}{2}\right).$$
(7)
By reducing the integrand of the multiple integral representation to certain canonical form, the EFP for $`n=3,4`$ BK1 ; BK2  and $`n=5`$ BKNS  were evaluated by Boos et al. (see also recent progress for $`n=6`$ BKS ). Note that $`P(2)`$ and $`P(3)`$ in (7) are related to the first and second neighbor correlators as $`P(2)=1/4+S_j^zS_{j+1}^z`$ and $`P(3)=1/8+S_j^zS_{j+1}^z+S_j^zS_{j+2}^z/2`$.

Here we quote the explicit form of $`P(4)`$ obtained in BK1 ; BK2 , which is closely related to the third neighbor correlator $`S_j^zS_{j+3}^z`$.

$`P(4)=`$ $`{\displaystyle \frac{1}{16}}+{\displaystyle \frac{3}{4}}S_j^zS_{j+1}^z+{\displaystyle \frac{1}{2}}S_j^zS_{j+2}^z+{\displaystyle \frac{1}{4}}S_j^zS_{j+3}^z`$
$`+S_j^zS_{j+1}^zS_{j+2}^zS_{j+3}^z`$
$`=`$ $`{\displaystyle \frac{1}{5}}2\mathrm{ln}2+{\displaystyle \frac{173}{60}}\zeta (3){\displaystyle \frac{11}{6}}\zeta (3)\mathrm{ln}2{\displaystyle \frac{51}{80}}\zeta (3)^2`$
$`{\displaystyle \frac{55}{24}}\zeta (5)+{\displaystyle \frac{85}{24}}\zeta (5)\mathrm{ln}2.`$ (8)
Note that on the antiferromagnetic ground state without magnetic field, all the correlators with an odd number of $`S^z`$ vanishes. Substituting (2) and (3) into (8), one finds the relation between the third neighbor correlator $`S_j^zS_{j+3}^z`$ and the four point correlator $`S_j^zS_{j+1}^zS_{j+2}^zS_{j+3}^z`$. However the exact value of $`S_j^zS_{j+3}^z`$ itself can not be determined solely from $`P(4)`$.

To determine $`S_j^zS_{j+3}^z`$, we consider the following auxiliary correlator:

$`P_{++}^{++}=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \frac{3}{4}}S_j^zS_{j+1}^z+{\displaystyle \frac{1}{2}}S_j^zS_{j+2}^z{\displaystyle \frac{1}{4}}S_j^zS_{j+3}^z`$
$`+S_j^zS_{j+1}^zS_{j+2}^zS_{j+3}^z.`$ (9)
Here and hereafter $`P_{\epsilon _1\epsilon _2\epsilon _3\epsilon _4}^{\stackrel{~}{\epsilon }_1\stackrel{~}{\epsilon }_2\stackrel{~}{\epsilon }_3\stackrel{~}{\epsilon }_4}`$ (also written as $`P_\epsilon ^{\stackrel{~}{\epsilon }}`$ for simplicity) denotes a correlator of the form

$$P_{\epsilon _1\epsilon _2\epsilon _3\epsilon _4}^{\stackrel{~}{\epsilon }_1\stackrel{~}{\epsilon }_2\stackrel{~}{\epsilon }_3\stackrel{~}{\epsilon }_4}=E_1^{\stackrel{~}{\epsilon }_1\epsilon _1}E_2^{\stackrel{~}{\epsilon }_2\epsilon _2}E_3^{\stackrel{~}{\epsilon }_3\epsilon _3}E_4^{\stackrel{~}{\epsilon }_4\epsilon _4},$$
(10)
where $`\epsilon _j`$, $`\stackrel{~}{\epsilon }_j=\{+,\}`$ and $`E_j^{\epsilon _j\stackrel{~}{\epsilon }_j}`$ is the 2 $`\times `$ 2 elementary matrix ($`E^{\pm \pm }=\pm S^z+1/2`$, $`E^+=S^+`$, $`E^+=S^{}`$) acting on $`j`$th site. In this notation $`P(4)=P_{++++}^{++++}`$. Note that because the Hamiltonian (1) has the symmetry under the group $`SU(2)`$, $`P_\epsilon ^{\stackrel{~}{\epsilon }}`$ possesses a property like

$$P_\epsilon ^{\stackrel{~}{\epsilon }}=P_{\stackrel{~}{\epsilon }}^\epsilon =P_\epsilon ^{\stackrel{~}{\epsilon }}=P_{\stackrel{~}{\epsilon }}^\epsilon .$$
(11)
As in the case of $`P(4)`$, the correlators (10) enjoy the multiple integral representation JMMN ; JMbook ; Naka ; JM :

$$P_\epsilon ^{\stackrel{~}{\epsilon }}=\underset{j=1}{\overset{4}{}}_C\frac{d\lambda _j}{2\pi i}U(\lambda _1,\mathrm{},\lambda _4)T(\lambda _1,\mathrm{},\lambda _4),$$
(12)
where the integration contour $`C`$ is taken to be a line $`[\mathrm{}i\alpha ,\mathrm{}i\alpha ]`$ $`(0<\alpha <1)`$. For convenience, we choose $`\alpha =1/2`$. The integrand $`U(\lambda _1,\mathrm{},\lambda _4)`$ is given by

$$U(\lambda _1,\mathrm{},\lambda _4)=\pi ^{10}\frac{_{1k<j4}\mathrm{sinh}\pi \lambda _{jk}}{_{j=1}^4\mathrm{sinh}^4\pi \lambda _j},$$
(13)
while $`T(\lambda _1,\mathrm{},\lambda _4)`$ depends on the selection of $`\epsilon `$ and $`\stackrel{~}{\epsilon }`$. Here and hereafter we use the notation $`\lambda _{jk}=\lambda _j\lambda _k`$ to save space. In particular for the correlator (9), $`T(\lambda _1,\mathrm{},\lambda _4)`$ is given by

$$T(\lambda _1,\mathrm{},\lambda _4)=\frac{\lambda _1(\lambda _1+i)^2\lambda _2^3\lambda _3(\lambda _3+i)^2\lambda _4^3}{(\lambda _{21}i)\lambda _{31}\lambda _{41}\lambda _{32}\lambda _{42}(\lambda _{43}i)}.$$
(14)
To calculate the multiple integral (12), we follow the method by Boos and Korepin BK2 . Roughly, their method is described as follows. First Taking carefully into account the property of $`U(\lambda _1,\mathrm{},\lambda _4)`$, we modify the integrand $`T(\lambda _1,\mathrm{},\lambda _4)`$ such that the integral gives the same result as the original one (“weak equivalence”). In this way it is likely that the integrand $`T(\lambda _1,\mathrm{},\lambda _4)`$ can always be reduced to the following form (we call it “canonical form”):

$`T_c=`$ $`P_0(\lambda _2,\lambda _3,\lambda _4)+{\displaystyle \frac{P_1(\lambda _1,\lambda _3,\lambda _4)}{\lambda _{21}}}+{\displaystyle \frac{P_2(\lambda _1,\lambda _3)}{\lambda _{21}\lambda _{43}}},`$ (15)
where $`P_0`$, $`P_1`$ and $`P_2`$ are certain polynomials. Once one derives the canonical form, one can perform the multiple integral (12) by using the Cauchy theorem BK2 . The result is written as combinations of the logarithm $`\mathrm{ln}2`$, the Riemann zeta function $`\zeta (3)`$ and $`\zeta (5)`$ and rational numbers. Consequently, the main part of the calculation for the multiple integral (12) reduces to finding the canonical form (15).

Now we consider the case (14) and show that the Boos–Korepin method is also applicable to our case. Let us introduce the following diagram.

$$\text{}:=\frac{1}{\lambda _{jk}},\text{}:=\frac{1}{\lambda _{jk}i},$$

where $`j>k`$. First we expand (14) through partial fractions. The result consists of 24 terms. Taking into account the antisymmetry of the function (13) under transposition of any two variables $`\lambda _j`$ and $`\lambda _k`$ and the symmetry of $`T`$ (14) under $`\lambda _1\lambda _3`$, $`\lambda _2\lambda _4`$, we can reduce the 24 partial fractions to eight ones. Diagrammatically, its denominators are written as

$`\text{}`$ $`\text{}\text{}+i\left(\text{}+\text{}\right)`$
$`+2i`$ $`\left(\text{}\text{}\text{}\text{}\right).`$
From a symmetry of the denominator, the second term can further be simplified as

$`g\text{}{\displaystyle \frac{1}{4}}(g(\lambda _1,\lambda _3,\lambda _2,\lambda _4)g(\lambda _4,\lambda _1,\lambda _3,\lambda _2)`$
$`+g(\lambda _2,\lambda _4,\lambda _1,\lambda _3)g(\lambda _3,\lambda _2,\lambda _4,\lambda _1))\text{}`$
$`{\displaystyle \frac{1}{2}}\lambda _1\lambda _2\lambda _3\lambda _4(i\lambda _1+i\lambda _3+2\lambda _1\lambda _3)(i\lambda _2+i\lambda _4+2\lambda _2\lambda _4)`$
$`\times \text{}\times \text{},`$
where $`g`$ denotes the numerator of (14). Next using the Cauchy theorem, we shift variables $`\lambda _j\lambda _j\pm i`$ such that the denominators do not contain $`i`$. For instance, we have

$`g\text{}`$ $`=g^{(1)}\text{}+g^{(2)}\text{}`$
$`g^{(1)}\text{}+g_{\lambda _2\lambda _2+i}^{(2)}\text{},`$
where

$`g^{(1)}`$ $`=\lambda _1(\lambda _1+i)\lambda _2^3\lambda _3(\lambda _3+i)^2\lambda _4^3`$
$`g^{(2)}`$ $`=\lambda _1(\lambda _1+i)\lambda _2^4\lambda _3(\lambda _3+i)^2\lambda _4^3.`$
Finally again using the antisymmetric property of (13), we eliminate the symmetric part with respect to $`\lambda _j\lambda _k`$. In this way we have obtain the canonical form as

$`P_0=`$ $`{\displaystyle \frac{56}{5}}\lambda _2\lambda _3^2\lambda _4^3,`$
$`P_1=`$ $`{\displaystyle \frac{27}{10}}\lambda _4i\lambda _4^2+{\displaystyle \frac{33}{5}}\lambda _3\lambda _4^2+{\displaystyle \frac{4}{5}}\lambda _4^3+2i\lambda _3\lambda _4^3+4\lambda _3^2\lambda _4^3`$
$`+\lambda _1(4i\lambda _4+7\lambda _4^232i\lambda _3\lambda _4^210i\lambda _4^312i\lambda _3^2\lambda _4^3)`$
$`+\lambda _1^2(4\lambda _419i\lambda _4^228\lambda _3\lambda _4^210\lambda _4^3`$
$`28i\lambda _3\lambda _4^332\lambda _3^2\lambda _4^3),`$
$`P_2=`$ $`{\displaystyle \frac{3}{10}}+{\displaystyle \frac{3}{2}}i\lambda _3+{\displaystyle \frac{3}{2}}\lambda _1\lambda _3{\displaystyle \frac{1}{2}}\lambda _3^2+4i\lambda _1\lambda _3^2+6\lambda _1^2\lambda _3^2.`$
Subsequently, applying the method developed in BK2 , we calculate the multiple integral by substituting the above canonical form into (12). Explicitly,

$`P_{++}^{++}=J_0+J_1+J_2`$
$`J_0={\displaystyle \frac{7}{10}},J_1={\displaystyle \frac{2}{3}}+{\displaystyle \frac{3}{10}}\zeta (3)+{\displaystyle \frac{35}{32}}\zeta (5),`$
$`J_2={\displaystyle \frac{1}{2}}\zeta (3)+{\displaystyle \frac{1}{2}}\zeta (3)\mathrm{ln}2+{\displaystyle \frac{9}{80}}\zeta (3)^2{\displaystyle \frac{25}{32}}\zeta (5)`$
$`{\displaystyle \frac{5}{8}}\zeta (5)\mathrm{ln}2,`$ (16)
where $`J_k`$ denotes the result of the integration regarding the term $`P_k`$. We remark that the canonical form is not unique due to the non-uniqueness of partial fraction expansions. Accordingly, the explicit value of each $`J_k`$ depends on the choice of the canonical form. The final result $`J_0+J_1+J_2`$, however, is always unique as a matter of course.

Combining the result (16) with (8) and (9), we obtain the third neighbor correlator $`S_j^zS_{j+3}^z`$ (4) and at the same time the correlator $`S_j^zS_{j+1}^zS_{j+2}^zS_{j+3}^z`$ as

$`S_j^zS_{j+1}^zS_{j+2}^zS_{j+3}^z`$ $`={\displaystyle \frac{1}{80}}{\displaystyle \frac{1}{3}}\mathrm{ln}2+{\displaystyle \frac{29}{30}}\zeta (3){\displaystyle \frac{2}{3}}\zeta (3)\mathrm{ln}2`$
$`{\displaystyle \frac{21}{80}}`$ $`\zeta (3)^2{\displaystyle \frac{95}{96}}\zeta (5)+{\displaystyle \frac{35}{24}}\zeta (5)\mathrm{ln}2.`$ (17)
Now let us consider the four-point correlators of the form $`S_j^\alpha S_{j+1}^\beta S_{j+2}^\gamma S_{j+3}^\delta `$ ($`\{\alpha ,\beta ,\gamma ,\delta \}\{x,y,z,0\};S_j^0=𝟙`$). Because the correlators with an odd number of $`S^{\{\alpha \}\mathrm{}0}`$ vanish, the possible non-zero correlators are restricted to the following three types: $`S_j^\alpha S_{j+1}^\alpha S_{j+2}^\beta S_{j+3}^\beta `$, $`S_j^\alpha S_{j+1}^\beta S_{j+2}^\alpha S_{j+3}^\beta `$ and $`S_j^\alpha S_{j+1}^\beta S_{j+2}^\beta S_{j+3}^\alpha `$. Further due to the isotropy of the Hamiltonian (1), one can find that the independent correlators are written as the following seven ones: $`S_j^xS_{j+1}^zS_{j+2}^zS_{j+3}^x`$, $`S_j^zS_{j+1}^xS_{j+2}^zS_{j+3}^x`$, $`S_j^zS_{j+1}^zS_{j+2}^xS_{j+3}^x`$ and already obtained ones in (2)–(4) and (17). Then we shall calculate the remaining three correlators here. For convenience we use the operator $`S^\pm (=S^x\pm iS^y)`$ instead of $`S^x`$. First we consider the correlator $`S_j^+S_{j+1}^zS_{j+2}^zS_{j+3}^{}(=2S_j^xS_{j+1}^zS_{j+2}^zS_{j+3}^x)`$. From (10) and the property (11), this correlator is expressed as $`S_j^+S_{j+1}^zS_{j+2}^zS_{j+3}^{}=\left(P_{+++}^{+++}P_{++}^{++}\right)/2`$. Using the relation $`S_j^+S_{j+3}^{}=2(P_{+++}^{+++}+P_{++}^{++})`$, we obtain

$$S_j^+S_{j+1}^zS_{j+2}^zS_{j+3}^{}=P_{+++}^{+++}\frac{1}{2}S_j^zS_{j+3}^z.$$
(18)
Here we have used the relation $`S_j^+S_{j+3}^{}=2S_j^zS_{j+3}^z`$. Similarly the other correlators are given by

$`S_j^zS_{j+1}^+S_{j+2}^zS_{j+3}^{}`$ $`=P_{+++}^{+++}{\displaystyle \frac{1}{2}}S_j^zS_{j+2}^z,`$
$`S_j^zS_{j+1}^zS_{j+2}^+S_{j+3}^{}`$ $`=P_{+++}^{+++}{\displaystyle \frac{1}{2}}S_j^zS_{j+1}^z.`$ (19)
Therefore our goal is to evaluate the auxiliary correlators $`P_{+++}^{+++}`$, $`P_{+++}^{+++}`$ and $`P_{+++}^{+++}`$. They are given if we replace the integrand $`T(\lambda _1,\mathrm{},\lambda _4)`$ by

$$T^{(l)}=\frac{(\lambda _1+i)^3\lambda _2(\lambda _2+i)^2(\lambda _3+i)\lambda _3^2\lambda _4^{4l}(\lambda _4+i)^{l1}}{(\lambda _{21}i)(\lambda _{31}i)(\lambda _{32}i)\lambda _{41}\lambda _{42}\lambda _{43}}.$$

in the multiple integral representation. Here the correlator $`P_{+++}^{+++}`$, $`P_{+++}^{+++}`$ and $`P_{+++}^{+++}`$ correspond to $`l=1`$, 2 and 3, respectively. Using the procedure similar to the case of $`P_{++}^{++}`$, one obtains the explicit values of the above auxiliary correlators. As a result, combining the identity (18) and (19) with (2)–(4), we arrive at

$`S_j^+S_{j+1}^zS_{j+2}^zS_{j+3}^{}`$ $`={\displaystyle \frac{1}{120}}{\displaystyle \frac{1}{2}}\mathrm{ln}2+{\displaystyle \frac{169}{120}}\zeta (3)`$
$`{\displaystyle \frac{5}{6}}\zeta (3)\mathrm{ln}2`$ $`{\displaystyle \frac{3}{10}}\zeta (3)^2{\displaystyle \frac{65}{48}}\zeta (5)+{\displaystyle \frac{5}{3}}\zeta (5)\mathrm{ln}2,`$ (20)
$`S_j^+S_{j+1}^zS_{j+2}^{}S_{j+3}^z`$ $`={\displaystyle \frac{1}{120}}{\displaystyle \frac{1}{3}}\mathrm{ln}2+{\displaystyle \frac{77}{60}}\zeta (3)`$
$`{\displaystyle \frac{5}{6}}\zeta (3)\mathrm{ln}2`$ $`{\displaystyle \frac{3}{10}}\zeta (3)^2{\displaystyle \frac{65}{48}}\zeta (5)+{\displaystyle \frac{5}{3}}\zeta (5)\mathrm{ln}2,`$ (21)
$`S_j^+S_{j+1}^{}S_{j+2}^zS_{j+3}^z`$ $`={\displaystyle \frac{1}{120}}+{\displaystyle \frac{1}{6}}\mathrm{ln}2{\displaystyle \frac{91}{120}}\zeta (3)`$
$`+{\displaystyle \frac{1}{3}}\zeta (3)\mathrm{ln}2+{\displaystyle \frac{3}{40}}`$ $`\zeta (3)^2+{\displaystyle \frac{35}{48}}\zeta (5){\displaystyle \frac{5}{12}}\zeta (5)\mathrm{ln}2.`$ (22)
We mention a few remarks of our results. (i) All the above correlators are written as the logarithm $`\mathrm{ln}2`$, the Riemann zeta function $`\zeta (3)`$ and $`\zeta (5)`$. This agrees with the general conjecture by Boos and Korepin: arbitrary correlators of the $`XXX`$antiferromagnet are described as certain combinations of logarithm $`\mathrm{ln}2`$, the Riemann zeta function with odd arguments and rational coefficients. Especially intriguing is the existence of the non-linear terms such as $`\zeta (3)^2`$, $`\zeta (3)\mathrm{ln}2`$ and $`\zeta (5)\mathrm{ln}2`$. (ii) The correlator (20) is interpreted as the third neighbor one-particle Green functions $`c_j^{}c_{j+3}_\mathrm{f}/4`$ via the Jordan-Wigner transformation (6). Obviously, the first neighbor one-particle Green function is expressed as

$$c_j^{}c_{j+1}_\mathrm{f}=\frac{1}{6}\frac{2}{3}\mathrm{ln}20.295431453707,$$

which coincides with $`S_j^+S_{j+1}^{}`$. Due to the characteristic $`\pi /2`$-oscillation: $`c_j^{}c_k_\mathrm{f}A_{jk}\mathrm{cos}(\pi (kj+1)/2)`$, one finds $`c_j^{}c_{j+2k}_\mathrm{f}=0`$. Therefore the quantity (5) is the first non-trivial exact result of the correlators containing the fermionic nature. (iii) The difference between (20) and (21) gives the nearest chiral correlator

$`(𝑺_j\times `$ $`𝑺_{j+1})(𝑺_{j+2}\times 𝑺_{j+3})=3(S_j^+S_{j+1}^zS_{j+2}^{}S_{j+3}^z`$
$`S_j^+S_{j+1}^zS_{j+2}^zS_{j+3}^{})={\displaystyle \frac{1}{2}}\mathrm{ln}2{\displaystyle \frac{3}{8}}\zeta (3),`$
which exactly agrees with the one derived from the ground state energy of an integrable two-chain model with four-body interactions MT .

To confirm the validity of our formulae, we performed numerical calculations by using the density-matrix renormalization group (DMRG) White92 ; White93  and numerical diagonalization. As for the DMRG, we followed standard algorithm Peschel99 . We have repeated renormalization 500-times. At each renormalization, we kept, at most, 200 relevant states for a (new) block. The numerical diagonalization was performed for the system size $`L=24`$, 28 and $`32`$. We extrapolate the data from a fitting function $`a_0+a_1/L^2+a_2/L^4`$. All our analytical results coincide quite accurately with both numerical ones (TABLE I).

In closing we would like to comment on generalizations of the present results. The extension to the calculation of higher neighbor correlators $`S_j^zS_{j+k}^z_{k4}`$ is of great interest. The fourth neighbor one $`S_j^zS_{j+4}^z`$, for example, will be calculated by combination of the EFP, $`P(5)`$ and two independent auxiliary correlators, which can in principle be evaluated. In fact $`P(5)`$ has been already obtained in BKNS . The computation, however, will be much more complicated. Alternatively, extending the present result to the inhomogeneous case as in BKS  and taking into account the property of the quantum Knizhnik-Zamolodchikow equation, we may derive higher neighbor correlators. Using this, eventually we hope to extract the long-distance asymptotics $`S_j^zS_{j+k}^z_{k1}`$, which is a crucial problem in conformal field theory Affleck ; LT .

The authors are grateful to H.E. Boos and V.E. Korepin for many valuable discussions. KS is supported by the JSPS research fellowships for young scientists. MS and YN are supported by Grant-in-Aid for young scientists No. 14740228 and No. 13740240, respectively. This work is in part supported by Grant-in Aid for the Scientific Research (B) No. 14340099 from the Ministry of Education, Culture, Sports, Science and Technology, Japan.


format: text

text:
# The Custodial Symmetry

## 1 Introduction

The Standard Model (SM) of the strong and electroweak interactions which has the $`SU(3)_CSU(2)_LU(1)_Y`$ gauge group symmetry has been quite succesful. Its predictions are in excellent agreement with the experimental data. Moreover, the model has predicted new particles like the $`Z`$ gauge boson which is the quantum of the weak interaction with zero electric charge. It has also predicted the charm quark which was introduced to avoid the flavor changing neutral currents at tree level, and top quark. However, the SM has some problems. It does not unify the strong and electroweak interactions and has three coupling constants associated to each gauge group. Currently there are some theories which unify these interactions and they only contain one coupling constant and a simple Lie group; however, they predict a proton decay faster than the experimental limit. Some of these groups are: $`SU(5)`$, $`SU(6)`$, $`SO(10)`$, $`E_6`$ etc.

One of the most interesting problems of the SM is that there is no explanation of the origin of masses, mixing of the fermions and the number of families appears as an identical pattern, although there are some models which introduce horizontal symmetries to explain the relation of the fermion masses by spontaneous symmetry breaking (SSB). The physicists expect to find a symmetry which relates the fermion generations; the general aim is to distinguish among them by breaking this symmetry and generating hierarchical mass scales that give different radiative corrections to the fermions. This symmetry would commute with the SM or a simple group that contains it. Some of the simple gauge groups that unify the SM with horizontal symmetries are: $`SU(11)`$, $`SO(14)`$, $`SO(18)`$, $`E_8`$, etc. At the present time there is not a realistic model that can explain satisfactorily the problem of the fermion masses.

To give mass to the fermions and the gauge bosons in the SM, it is necessary to introduce a multiplet of scalar fields which has some components that are ‘eaten‘ by the gauge bosons, the would-be Goldstone bosons, and another component which is called the Higgs particle. This is a physical degree of freedom, but the theory cannot predict its mass. There are some energy regions to look for the Higgs according to the production and decay channels. Up to now there is not evidence of it.

The left (right) quirality fermions in the SM transform according to the global symmetry $`SU(2)_{L(R)}`$ and the Higgs doublet, four real fields, is a bidoublet under this global symmetry. Before the breaking of $`SU(2)_LU(1)_Y`$ , the Higgs potential has a global $`SU(2)_LSU(2)_R`$ symmetry which reduces to $`SU(2)_V`$when the symmetry is broken. This remanent global symmetry is called the ”custodial symmetry”. Additionally, it is possible to write down the Yukawa Lagrangian to see the same $`SU(2)_L\times SU(2)_R`$ global symmetry explicitly.

In the present work, we explain the meaning of the custodial symmetry specially when dealing with Higgs multiplets appearing in extensions of the SM. We also explain that if the gauge symmetry of the electroweak model is broken by the Higgs doublet, there is a custodial symmetry which protects the mass relation of the $`W`$ and $`Z`$ gauge fields, i. e.,

$$\rho \frac{M_W^2}{M_Z^2\mathrm{cos}^2\theta }=1.$$
(1)
This relation is also valid after the radiative corrections if the custodial symmetry is not broken. Experimentally this relation is satisfied at the $`1\%`$ level, which restricts the new physics beyond the SM and allows us to distinguish among different models.

## 2 The gauge field masses

The kinetic energy parts of the Lagrangian of the scalar fields produce the mass of the gauge fields after the SSB. If we have an arbitrary scalar multiplet of $`SU(2)_LU(1)_Y`$, the gauge field masses arise from

$$D_\mu <\varphi >_oig\stackrel{}{A_\mu }[\stackrel{}{T},<\varphi >_o]+i\frac{g^{}}{2}B_\mu [Y,<\varphi >_o]$$

where $`<\varphi >_o`$ is the vacuum expectation value (v.e.v.) of the scalar multiplet, i.e., $`<\varphi >_o=<0|\varphi |0>`$. The covariant derivative for a general representation of $`SU(2)_LU(1)_Y`$ can be written as

$`D_\mu <\varphi >_o`$ $`i`$ $`{\displaystyle \frac{g}{\sqrt{2}}}W_\mu ^+[T^+,<\varphi >_o]+i{\displaystyle \frac{g}{\sqrt{2}}}W_\mu ^{}[T^{},<\varphi >_o]`$
$`+`$ $`ig`$ $`A_{3\mu }[T_{3L},<\varphi >_o]+i{\displaystyle \frac{g^{}}{2}}B_\mu [Y,<\varphi >_o].`$ (2)
where we have defined $`W_\mu ^\pm =(A_\mu ^1iA_\mu ^2)/\sqrt{2}`$.

The SM has a $`U(1)_Q`$ remanent symmetry and $`<\varphi >_o`$ is different from zero because of the neutral component, so that we get

$$T_{3L}<\varphi >_o=\frac{1}{2}Y<\varphi >_o.$$
(3)
Using the Lie algebra, the mass terms of the gauge fields can be written as

$$\left(\frac{g^2}{2}\left[t(t+1)t_{3L}^2\right]\left(A_\mu ^1A^{1\mu }+A_\mu ^2A^{2\mu }\right)+t_{3L}^2\left(gA_{3\mu }g^{}B_\mu \right)^2\right)<\varphi >_o^2,$$
(4)
where $`t(t+1)`$ and $`t_{3L}`$ are the eigenvalues of the $`_{i=1}^3T_{iL}^2`$ and $`T_{3L}`$ operators, respectively. Under the assumption that $`g^{}=0`$ or equivalently $`\mathrm{sin}\theta _W=0,`$ and that the $`SU(2)_L`$ gauge fields transform as a triplet of the $`SU(2)`$ global symmetry, and the kinetic energy Lagrangian is invariant under these transformation, the following relation is satisfied

$$t(t+1)t_{3L}^2=2t_{3L}^2.$$
(5)
This symmetry is only exact when the $`B_\mu `$ fields is zero because it breaks the custodial symmetry.

From the above relations and the definition of the mass terms

$`M_W^2`$ $`=`$ $`g^2\left(t(t+1)t_{3L}^2\right)<\varphi >_o^2,`$ (6)
$`M_Z^2`$ $`=`$ $`2(g^2+g^{}_{}{}^{}2)t_{3L}^2<\varphi >_o^2,`$
for the $`\rho `$ parameter we get

$`\rho `$ $`=`$ $`{\displaystyle \frac{M_W^2}{M_Z^2\mathrm{cos}^2\theta _W}}`$ (7)
$`=`$ $`{\displaystyle \frac{t(t+1)t_{3L}^2}{2t_{3L}^2}}{\displaystyle \frac{g^2}{(g^2+g^{}{}_{}{}^{2})}}{\displaystyle \frac{1}{\mathrm{cos}^2\theta _W}}`$
$`=`$ $`1`$
where $`\mathrm{tan}\theta _W=g^{^{}}/g`$ and $`Z_\mu =\mathrm{cos}\theta _WA_{3\mu }\mathrm{sin}\theta _WB_\mu `$. Taking into account the eq.(5), we can get different values for $`t`$ when $`\rho =1`$

$`t`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1+\sqrt{1+12t_{3L}^2}\right)`$ (8)
$`=`$ $`0,{\displaystyle \frac{1}{2}},\mathrm{\hspace{0.33em}3},\mathrm{}.`$
The first and second solutions correspond to singlets and doublets, respectively. The other solutions are not important in the SM because they do not couple to the fermions in the Yukawa Lagrangian and they do not produce fermion mass terms. However, for other fermion representations, physics beyond the SM, Higgs fields with other values of $`t`$ are important to construct Yukawa terms. One interesting example is the see-saw mechanism implemented with new fermions and triplet Higgs fields. In the minimal SM with only left handed neutrino, it is neccesary a Higgs triplet to get a massive neutrino; however the vev of the new Higgs representation is restricted by the experimental value of the $`\rho `$ parameter. With the above arguments we can say that if we use doublets to break the SM, the $`SU(2)_L`$ gauge fields will be a triplet of the $`SU(2)`$ global symmetry or the custodial symmetry and $`\rho `$ is equal to one. Radiative corrections to the $`\rho `$ parameter by gauge fields violate softly this symmetry.

## 3 The Higgs and fermions Lagrangian

The most general Higgs potential which is renormalizable and invariant under $`SU(2)_LU(1)_Y`$ gauge transformations have the form

$$V=\lambda \left(\varphi ^{}\varphi \mu ^2\right)^2,$$
(9)
where

$$\varphi =\frac{1}{\sqrt{2}}\left(\begin{array}{c}\varphi _1+i\varphi _2\hfill \\ \varphi _4+i\varphi _3\hfill \end{array}\right).$$

The potential as function of the $`\varphi _i`$ scalar fields can be written as

$$V=\frac{\lambda }{4}\left(\varphi _1^2+\varphi _2^2+\varphi _3^2+\varphi _4^22\mu ^2\right)^2.$$
(10)
It is invariant under rotations of the four fields which lead to $`SO(4)`$ as the global symmetry group. This group is isomorphic to $`SU(2)_LSU(2)_R`$ because both have the same Lie algebra. This symmetry is global and it does not necessary to introduce gauge fields.

When the symmetry is broken, the scalar field $`\varphi _4`$ get a v.e.v different from zero, and it can be redefined as follows

$$\varphi _4=H+v$$
(11)
where $`H`$ gets its mass and is called the Higgs particle. Moreover, it has a v.e.v equal to zero. This field is a physical degree of freedom and its mass is proportional to the $`\lambda `$ parameter which is unknown in the model. The other scalar fields remain massless. They are the would-be Goldstone bosons and correspond to the degrees that the gauge fields ‘eat‘ in order to get mass or longitudinal component.

The Higgs potential can be written as a function of the new fields as follows

$$V=\frac{\lambda }{4}\left(\varphi _1^2+\varphi _2^2+\varphi _3^2+H^2+2Hv+v^22\mu ^2\right)^2.$$
(12)
In this new potential the global symmetry is broken to $`SO(3)`$, which only rotates three scalar fields. It is isomorphic to $`SU(2)_V`$, the diagonal part of $`SU(2)_LSU(2)_R`$. We can say that the symmetry was broken according to the following scheme

$$SU(2)_LSU(2)_RSU(2)_V.$$
(13)
To evidentiate the isomorphism between $`SO(4)`$ and $`SU(2)_LSU(2)_R`$, we will express the fields according to

$`\varphi ^{^{}}`$ $`=`$ $`U_L\varphi U_R^{},`$
$`\varphi `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}\hfill \varphi _4i\varphi _3& \hfill \varphi _1+i\varphi _2\\ \hfill (\varphi _1i\varphi _2)& \hfill \varphi _4+i\varphi _3\end{array}\right],`$ (16)
where $`U_{L(R)}SU(2)_{L(R)}`$. The elements of the $`SU(2)_L`$ group act over the columns and the elements of the $`SU(2)_R`$ group act over the rows. These transformations leave invariant the following quadratic form

$$Tr(\varphi ^{}\varphi ),$$
(17)
and, obviously, leave invariant the Higgs potential.

If $`\theta _L=\theta _R`$ the group of transformations becomes $`SU(2)_V`$, and if $`\theta _L=\theta _R`$ the group of transformations reduces to $`SU(2)_A`$. Under the $`SU(2)_V`$ infinitesimal transformations the fields of the multiplet transform as follows

$`\varphi ^{^{}}`$ $`=`$ $`e^{i\stackrel{}{\theta }\stackrel{}{T}}\varphi e^{i\stackrel{}{\theta }\stackrel{}{T}},`$
$`\delta \varphi _i`$ $``$ $`ϵ_{ijk}\theta _jT_k,i,j,k=1,2,3,`$ (18)
$`\delta \varphi _4`$ $`=`$ $`\delta H=0.`$
After the symmetry breaking, the quadratic form still invariant under $`SU(2)_V`$ is

$$\varphi _1^2+\varphi _2^2+\varphi _3^2=\varphi ^+\varphi ^{}+\varphi ^{}\varphi ^++\varphi _3^2,$$
(19)
The interactions between fermions and the Higgs fields are given by the Yukawa Lagrangian. It also gives mass to the fermions when the vacuum is aligned. The renormalizable Yukawa Lagrangian is giving by

$$L_Y=h_uQ_L\stackrel{~}{\varphi }u_R+h_dQ_L\varphi d_R+h.c.,$$
(20)
where $`Q_L^T=(u_Ld_L)`$, $`m_u=h_u<\varphi _o>`$ and $`m_d=h_d<\varphi _o>`$. If we assume that the quarks have the same masses, $`h_u=h_d`$, then the Yukawa Lagrangian can be written as

$$h(\begin{array}{cc}\hfill \overline{u}_L& \hfill \overline{d}_L\end{array})\left[\begin{array}{cc}\hfill \varphi _4i\varphi _3& \hfill \varphi _1+i\varphi _2\\ \hfill (\varphi _1i\varphi _2)& \hfill \varphi _4+i\varphi _3\end{array}\right]\left(\begin{array}{c}\hfill u_R\\ \hfill d_R\end{array}\right)+h.c..$$

To get $`L_Y`$ invariant under the global $`SU(2)_LSU(2)_R`$ symmetry, the quarks have to transform according to

$$Q_{L(R)}^{^{}}=U_{L(R)}Q_{L(R)}.$$
(21)
After SSB the quarks get mass and the global symmetry is broken down to $`SU(2)_V`$. The mass terms are given by

$$L_m=h(\overline{u}_Ru_L+\overline{d}_Rd_L)+h.c.,$$
(22)
which are invariant under $`U_L=U_R=SU(2)_V`$, i.e., we have the isospin symmetry. It is the same custodial symmetry.

## 4 Custodial symmetry and radiative corrections

Considering the contribution of the first doublet of fermions $`(u,d)`$ to the radiative corrections of the $`\rho `$ parameter, we have

$`\mathrm{\Delta }\rho `$ $`=`$ $`\mathrm{\Pi }_{ZZ}(0)\mathrm{\Pi }_{WW}(0),`$
$`\mathrm{\Pi }_{ZZ}(0)`$ $`=`$ $`{\displaystyle \frac{3G_F}{8\sqrt{2}\pi ^2}}(2m_u^2\mathrm{log}m_u^22m_d^2\mathrm{log}m_d^2),`$ (23)
$`\mathrm{\Pi }_{WW}(0)`$ $`=`$ $`{\displaystyle \frac{3G_F}{8\sqrt{2}\pi ^2}}\left(m_u^2+m_d^2+{\displaystyle \frac{2m_u^4\mathrm{log}m_u^22m_d^4\mathrm{log}m_d^2}{m_u^2m_d^2}}\right)`$
where $`\mathrm{\Pi }_{ii}(0)=_{ii}(0)/M_i^2`$ and $`_{ii}(0)`$ is the diagonal contribution to the self energy. Summarizing all terms yielded  we get

$`\mathrm{\Delta }\rho `$ $`=`$ $`{\displaystyle \frac{3G_F}{8\sqrt{2}\pi ^2}}F(m_u,m_d),`$
$`F(m_u,m_d)`$ $`=`$ $`m_u^2+m_d^2{\displaystyle \frac{2m_u^2m_d^2}{m_u^2m_d^2}}\mathrm{log}{\displaystyle \frac{m_u^2}{m_d^2}}.`$ (24)
In the limit $`m_u=m_d`$, i. e., when the isospin symmetry or custodial symmetry is valid the function $`F(m_u,m_d)`$ is equal to zero. The radiative corrections to $`\rho `$ are different from zero when the custodial symmetry is broken. For example, the third family of the SM where the up and down quarks are replaced by the top and the bottom quarks, respectively, the global symmetry is broken because $`m_tm_b`$ and the correction to $`\rho `$ is proportional to $`m_t^2`$

$$\mathrm{\Delta }\rho =\frac{3G_Fm_t^2}{8\sqrt{2}\pi ^2}.$$
(25)
## 5 Conclusions

When the SM symmetry is broken by Higgs doublets, there is a global symmetry that protects the mass relation of the $`W`$ and $`Z`$ fields, which transforms as the components of a triplet. This relation is violated when this symmetry is not exact. If the isospin symmetry or custodial symmetry is broken, then $`\rho `$ parameter may get radiative corrections different from zero. To look for new physics beyond the SM it is important to know if this symmetry is still exact.

We thank to COLCIENCIAS for financial support.


text:
# Exact SO(5) Symmetry in spin 3/2 fermionic system

## Abstract

The spin 3/2 fermion models with contact interactions have a generic SO(5) symmetry without any fine-tuning of parameters. Its physical consequences are discussed in both the continuum and lattice models. A Monte-Carlo algorithm free of the sign problem at any doping and lattice topology is designed when the singlet and quintet interactions satisfy $`U_0U_2\frac{3}{5}U_0(U_00)`$, thus making it possible to study different competing orders with high numerical accuracy. This model can be accurately realized in ultra-cold atomic systems.

With the rapid progress in ultra-cold atomic systems, many alkali fermions have been cooled below Fermi temperatures demarco1999a ; truscott2001a ; schreck2001a . All of them except 6Li have spins higher than 1/2 in the lowest hyper-fine multiplets. The spin degrees of freedom become free in the optical traps, which has attracted interest in their effects on Cooper pair structures and collective modes ho1999a ; yip1999a . The proposal of the optical lattice jacksh1998a  has led to a tremendous progress in studying the strongly correlated bosonic lattice systems anderson1998a ; cataliotti2001a ; morsch2001a ; greiner2002a ; greiner2002b . Recently, fermionic lattice systems are also exciting. For example, the degenerate 40K gas has been prepared in a one-dimensional optical lattice modugno2003a .

Comprehensive analysis of symmetries is helpful in understanding the physics in strongly correlated systems. For example, the SO(5) theory zhang1997a  of high Tc cuprates unifies the d-wave superconductivity (SC) and antiferromagnetism (AF) orders, leading to many experimental consequences. The sharp neutron scattering mode can be interpreted as the pseudo-Goldstone modezhang1997a . The prediction of the antiferromagnetic vortex core arovas1997a  has also been verified in recent experiments levi2002a .

In this article, we focus on the symmetry properties and corresponding consequences in the spin 3/2 system with contact interactions, including both the continuum model with $`s`$-wave scattering and the generalized lattice Hubbard model with on-site interactions. For neutral atoms, these interactions are generally described by two parameters in the total spin $`S_T=0,2`$ channels as $`g_{0,2}=4\pi \mathrm{}^2a_{0,2}/M`$ in the continuum model with $`a_{0,2}`$ the corresponding $`s`$-wave scattering lengths and $`M`$ the atom mass; or $`U_{0,2}`$ in the lattice model. Interactions in the odd total spin ($`S_T=1,3`$) channels are forbidden by Pauli’s exclusion principle. Remarkably, in addition to the explicit spin SU(2) symmetry, an enlarged SO(5) symmetry is present without any fine tuning of parameters. In the continuum model, this symmetry has direct consequences on the collective modes and pairing structures. In the lattice model, exact phase boundaries of various competing phases can be determined directly from symmetries. Because of the time-reversal symmetry of the Kramers doublets, a Monte-Carlo algorithm free of the notorious sign problem is designed when $`U_0U_23/5U_0(U_00)`$ at any filling level and lattice topology.

We start with the standard form of the spin $`3/2`$ Hamiltonian of the continuum model ho1999a ; yip1999a

$`H`$ $`=`$ $`{\displaystyle }d^d𝐫\{{\displaystyle \underset{\alpha =\pm 3/2,\pm 1/2}{}}\psi _\alpha ^{}(𝐫)({\displaystyle \frac{\mathrm{}^2}{2m}}^2\mu )\psi _\alpha (𝐫)`$ (1)
$`+`$ $`g_0P_{0,0}^{}(𝐫)P_{0,0}(𝐫)+g_2{\displaystyle \underset{m=\pm 2,\pm 1,0}{}}P_{2,m}^{}(𝐫)P_{2,m}(𝐫)\},`$
with $`d`$ the space dimension, $`\mu `$ the chemical potential and $`P_{0,0}^{},P_{2,m}^{}`$ the singlet ($`S_T=0`$) and quintet ($`S_T=2`$) pairing operators defined through the Clebsh-Gordan coefficient for two indistinguishable particles as $`P_{F,m}^{}(𝐫)=_{\alpha \beta }\frac{\mathrm{𝟑}}{\mathrm{𝟐}}\frac{\mathrm{𝟑}}{\mathrm{𝟐}};𝐅,𝐦|\frac{\mathrm{𝟑}}{\mathrm{𝟐}}\frac{\mathrm{𝟑}}{\mathrm{𝟐}}\alpha \beta \psi _\alpha ^{}(𝐫)\psi _\beta ^{}(𝐫),`$ where $`F=0,2`$ and $`m=F,F+1,\mathrm{},F`$.

We first construct the SO(5) algebra by introducing the five Dirac $`\mathrm{\Gamma }^a`$ $`(1a5)`$ matrices

$`\mathrm{\Gamma }^1=\left(\begin{array}{cc}0& iI\\ iI& 0\end{array}\right),\mathrm{\Gamma }^{2,3,4}=\left(\begin{array}{cc}\stackrel{}{\sigma }& 0\\ 0& \stackrel{}{\sigma }\end{array}\right),\mathrm{\Gamma }^5=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right).`$ (8)
Then the ten SO(5) generators are defined as $`\mathrm{\Gamma }^{ab}=\frac{i}{2}[\mathrm{\Gamma }^a,\mathrm{\Gamma }^b](1a,b5)`$, where $`I`$ and $`\stackrel{}{\sigma }`$ are the 2$`\times `$ 2 unit and Pauli matrices. The four-component spinor can be defined by $`\psi (𝐫)=(\psi _{\frac{3}{2}}(𝐫),\psi _{\frac{1}{2}}(𝐫),\psi _{\frac{1}{2}}(𝐫),\psi _{\frac{3}{2}}(𝐫))^T`$. Furthermore, the bilinear operators can be classified according to their properties under the SO(5) transformations. The 16 bilinear operators in the particle-hole (p-h) channel can be classified as SO(5)’s scalar, vector, and anti-symmetric tensors (generators) as

$`n(𝐫)=\psi _\alpha ^{}(𝐫)\psi _\alpha (𝐫),n_a(𝐫)={\displaystyle \frac{1}{2}}\psi _\alpha ^{}(𝐫)\mathrm{\Gamma }_{\alpha \beta }^a\psi _\beta (𝐫),`$
$`L_{ab}(𝐫)={\displaystyle \frac{1}{2}}\psi _\alpha ^{}(𝐫)\mathrm{\Gamma }_{\alpha \beta }^{ab}\psi _\beta (𝐫).`$ (9)
$`L_{ab}`$ and $`n_a`$ together form the SU(4), or isomorphically, the SO(6) generators. The spin SU(2) generators $`J_{x,y,z}`$ are expressed as $`J_+=J_x+iJ_y=\sqrt{3}(L_{34}+iL_{24})+(L_{12}+iL_{25})i(L_{13}+iL_{35})`$, $`J_{}=J_+^{}`$, and $`J_z=L_{23}+2L_{15}`$. $`n`$ and $`n_a`$ have spin 0 and 2, and $`L_{ab}`$ contains both the spin 1 and 3 parts. Pairing operators can also be organized as SO(5) scalar and vectors through the matrix $`R=\mathrm{\Gamma }_1\mathrm{\Gamma }_3`$

$`\eta ^{}(𝐫)`$ $`=`$ $`\text{Re}\eta +i\text{Im}\eta ={\displaystyle \frac{1}{2}}\psi _\alpha ^{}(𝐫)R_{\alpha \beta }\psi _\beta ^{}(𝐫),`$
$`\chi _a^{}(𝐫)`$ $`=`$ $`\text{Re}\chi _a+i\text{Im}\chi _a={\displaystyle \frac{i}{2}}\psi _\alpha ^{}(𝐫)(\mathrm{\Gamma }^aR)_{\alpha \beta }\psi _\beta ^{}(𝐫)`$ (10)
where $`P_{0,0}^{}=\eta ^{}/\sqrt{2}`$, and $`P_{2,\pm 2}^{}=(\chi _1^{}+i\chi _5^{})/2`$, $`P_{2,\pm 1}^{}=(\chi _3^{}\pm i\chi _2^{})/2`$, $`P_{2,0}^{}=i\chi _4^{}/\sqrt{2}`$. That is, $`\chi _a^{}`$ are polar combinations of $`J_z`$’s eigenoperators $`P_{2,m}^{}`$. The existence of the $`R`$ matrix is related to the pseudoreality of SO(5)’s spinor representation. It satisfies $`R^2=1,R^{}=R^1={}_{}{}^{t}R=R`$ and $`R\mathrm{\Gamma }^aR=^t\mathrm{\Gamma }^a,R\mathrm{\Gamma }^{ab}R={}_{}{}^{t}\mathrm{\Gamma }_{}^{ab}`$ scalapino1998a . The anti-unitary time-reversal transformation can be expressed as $`T=RC`$, where $`C`$ denotes complex conjugation and $`T^2=1`$. $`N`$, $`n_a`$, and $`L_{ab}`$ transform differently under the $`T`$ transformation

$`TnT^1=n,Tn_aT^1=n_a,TL_{ab}T^1=L_{ab}.`$ (11)
With the above preparation, the hidden SO(5) symmetry becomes manifest. The kinetic energy part has an explicit SU(4) symmetry which is the unitary transformation among four spin components. The singlet and quintet interactions are proportional to $`\eta ^{}(𝐫)\eta (𝐫)`$ and $`\chi _a^{}(𝐫)\chi _a(𝐫)`$ respectively, thus reducing the symmetry group from SU(4) to SO(5). When $`g_0=g_2`$, the SU(4) symmetry is restored because $`\chi _a^{},\eta ^{}`$ together form its 6 dimensional antisymmetrical tensor representation. In the continuum model, interactions in other even partial wave channels also keep the SO(5) symmetry. The odd partial wave scattering include spin 1 and 3 channel interactions $`g_1`$ and $`g_3`$, which together could form the 10-d adjoint representation. of SO(5), if and only $`g_1=g_3`$. However, to the leading order, $`p`$-wave scattering is much weaker than the $`s`$-wave one for neutral atoms, and can thus be safely neglected.

The SO(5) symmetry implies more degeneracies in the collective excitations in the spin 3/2 Landau fermi liquid theory, which generally requires four fermi liquid functions in total spin $`S_T=0,1,2,3`$ channels. The SO(5) symmetry of the microscopic Hamiltonian reduces these to three independent sets, classified according to the SO(5) scalar, vector, and tensor channels as

$`f_{\alpha \beta ,\gamma \delta }(p,p^{})`$ $`=`$ $`f_s(p,p^{})+f_v(p,p^{})(\mathrm{\Gamma }^a/2)_{\alpha \beta }(\mathrm{\Gamma }^a/2)_{\gamma \delta }`$ (12)
$`+`$ $`f_t(p,p^{})(\mathrm{\Gamma }^{ab}/2)_{\alpha \beta }(\mathrm{\Gamma }^{ab}/2)_{\gamma \delta }.`$
In other words, the effective interaction functions in the $`S_T=1,3`$ channels are exactly identical in all orders in perturbation theory. Furthermore, within the $`s`$-wave scattering approximation, the interaction functions become constants, and are given as $`f_s=(g_0+5g_2)/16`$, $`f_v=(g_03g_2)/4`$, $`f_t=(g_0+g_2)/4`$. Experiments in the fermi liquid regime can determine the four fermi liquid constants in the $`S_T=0,1,2,3`$ channels separately and verify the degeneracy between spin 1 and 3 channels. This degeneracy appears to be accidental in Ref. yip1999a . It is in fact exact and protected by the generic SO(5) symmetry.

The SO(5) symmetry also enriches the Cooper pair structures. Ref. ho1999a  showed that, in addition to the singlet pairing when $`g_0<0`$, the spin 3/2 system energetically favors the polar pairing state in the quintet channel when $`g_2<0`$ with the order parameter $`\mathrm{\Delta }_{real}^{}=\xi _1P_{2,0}^{}+\xi _2(P_{2,2}^{}+P_{2,2}^{})\xi _1\chi _4^{}+\xi _2\chi _5^{}`$ ($`\xi _{1,2}`$ are real). We understand that this is only a special case of the general pairing structures spanned by all the $`\chi _{15}^{}`$. The polar pairing states break the SO(5)$``$ U(1) (charge) symmetry to $`SO(4)Z_2`$, and thus the Goldstone (GS) manifold is the quotient space $`[SO(5)U(1)]/[SO(4)Z_2]=[S^4U(1)]/Z_2`$. Its dimension, 5, is the number of GS modes. When both $`g_{0,2}`$ are positive, s-wave pairing is not favorable. However, similarly to the spin fluctuation exchange mediated p-wave pairing in 3He anderson1973a ; leggett1975a , the spin fluctuations in the tensor channel can induce effective attractions between two atoms with total spin 1 and 3. This may lead to p-wave pairing where the spin part forms the 10-d adjoint representation of SO(5).

Now we consider the more interesting case of spin 3/2 fermions in the optical lattice. The periodic potential is $`V(x,y,z)=V_0(\mathrm{sin}^2(kx)+\mathrm{sin}^2(ky)+\mathrm{sin}^2(kz))`$ with $`V_0`$ the potential depth , $`k=\pi /l_0`$ the wavevector, and $`l_0`$ the lattice constant. The hopping integral $`t`$ between neighboring sites decreases exponentially with increasing $`V_0`$. Within the harmonic approximation, the parameter $`U/\mathrm{\Delta }E(\pi ^2/2)(a_s/l_0)(V_0/E_r)^{1/4}`$, with $`U`$ the repulsion of two fermions on one site, $`\mathrm{\Delta }E`$ the gap between the lowest and first excited single particle state in one site, $`a_s`$ the $`s`$-wave scattering length in the corresponding channel, and $`E_r=\mathrm{}^2k^2/2M`$ the recoil energy. With the typical estimate of $`a_s100a_B`$ ( $`a_B`$ the Bohr radius), $`l_05000`$A, and $`(V_0/E_r)^{1/4}12`$, we arrive at $`U/\mathrm{\Delta }E<0.1`$. Thus this system can be approximated by the one-band Hubbard model

$`H=t{\displaystyle \underset{ij,\sigma }{}}\{c_{i\sigma }^{}c_{j\sigma }+h.c.\}\mu {\displaystyle \underset{i\sigma }{}}c_{i\sigma }^{}c_{i\sigma }`$
$`+U_0{\displaystyle \underset{i}{}}P_{0,0}^{}(i)P_{0,0}(i)+U_2{\displaystyle \underset{i,m=\pm 2,\pm 1,0}{}}P_{2,m}^{}(i)P_{2,m}(i).`$ (13)
for particle density $`n4`$. At half-filling on a bipartite lattice, $`\mu `$ is given by $`\mu _0=(U_0+5U_2)/4`$ to ensure the particle-hole (p-h) symmetry under the transformation $`c_{i,\sigma }()^ic_{i,\sigma }^{}`$. The lattice fermion operators and their continuum counterparts are related by $`\psi _\alpha (𝐫)=c_\alpha (i)/(l_0)^{d/2}`$. We use the same symbols for bilinear fermion operators as in the continuum model.

The proof of SO(5) invariance in the continuum model applies equally well in the lattice model at any lattice topology and at any filling level. Eq. Exact SO(5) Symmetry in spin 3/2 fermionic system can be conveniently rewritten in another manifestly SO(5) invariant form as

$`H_0`$ $`=`$ $`t{\displaystyle \underset{i,j}{}}\{\psi ^{}(i)\psi (j)+h.c.\}`$ (14)
$`H_I`$ $`=`$ $`{\displaystyle \underset{i,1a5}{}}\left\{{\displaystyle \frac{3U_0+5U_2}{16}}(n(i)2)^2{\displaystyle \frac{U_2U_0}{4}}n_a^2(i)\right\}`$ (15)
$``$ $`(\mu \mu _0){\displaystyle \underset{i}{}}n(i),`$
where the SU(4) symmetry appears at $`U_0=U_2`$ as before.

The lattice Hamiltonian, Eq. Exact SO(5) Symmetry in spin 3/2 fermionic system, contains even higher symmetries under certain conditions. One can construct the largest SO(8)lin1998a  algebra using all the independent fermionic bilinear operators. Its generators $`M_{ab}(0a<b7)`$ including $`L_{ab}(1a<b5)`$ as its SO(5) sub-algebra, are denoted as

$`M_{ab}=\left(\begin{array}{cccc}0& \text{Re}\chi _1\text{Re}\chi _5& N& \text{Re}\eta \\ & & \text{Im}\chi _1& n_1\\ & L_{ab}& & \\ & & \text{Im}\chi _5& n_5\\ & & 0& \text{Im}\eta \\ & & & 0\end{array}\right),`$ (22)
with $`N=(n2)/2`$. Its Casimir is a constant $`C_{so8}=_{0a<b7}M_{ab}^2(i)=7`$. The global SO(8) generators are defined to be uniform in the p-h channel as $`M_{ab}=_iM_{ab}(i)`$ and staggered in the p-p channel as $`M_{ab}=_i()^iM_{ab}(i)`$ on the bipartite lattice. These global generators commute with the hopping term $`[M_{ab},H_0]=0`$. On the other hand, order parameters transformed under the SO(8) group should be staggered (uniform) in the p-h (p-p) channel respectively. The SO(8) symmetry is always broken by the interaction, but its subgroup symmetry, the SO(5)$``$ SU(2) and SO(7) symmetries, appear under special conditions as shown below.

At $`U_0=5U_2`$, $`H_I`$ can be rewritten as $`H_I=_{i,1a,b5}\{U_2L_{ab}^2(i)(\mu \mu _0)n(i)\}`$, using the Fierz identity $`_{1a5}L_{ab}^2(i)+_{1a5}n_a^2(i)+5N^2(i)=5`$. As a generalization of the pseudospin algebra in the usual Hubbard model yang1989a , we construct them as $`\eta ^{},\eta ,N`$. The symmetry at half-filling is SO(5) $``$ SU(2), which unifies the charge density wave (CDW) and the singlet pairing (SP) order parameters. Away from half-filling, this symmetry is broken but $`\eta ,\eta ^{}`$ are still eigen-operators since $`[H,\eta ^{}]=(\mu \mu _0)\eta ^{}`$, and $`[H,\eta ]=(\mu \mu _0)\eta `$.

The p-h channel SO(5) $``$ U(1) symmetry can also be extended to SO(7) at $`U_0=3U_2`$ where $`H_I`$ can be rewritten as $`H_I=_{i,0a<b6}\left\{\frac{2}{3}U_2M_{ab}(i)^2(\mu \mu _0)n(i)\right\}`$. The SO(7) symmetry is exact at half-filling. Its 7-d vector Rep unifies the staggered 5-vector and SP order parameters. Its 21-d adjoint Rep unifies the staggered SO(5) adjoint Rep order parameters, CDW, and quintet pairing (QP) order parameters. Away from half-filling, QP operators are spin 2 quasi-Goldstone operators $`[H,\chi _a^{}]=(\mu \mu _0)\chi _a^{}`$ and $`[H,\chi _a]=(\mu \mu _0)\chi _a`$. These $`\chi `$ modes are just the analogs of the $`\pi `$ modes in the high Tc contextzhang1997a .

The above symmetry structures guide the mean field (MF) analysis of the phase diagram. In the weak coupling limit, the complete MF decoupling is performed in the direct, exchange, and pairing channels. We take the MF ansatz on the 2d square lattice

$`n_a(i)=()^i\overline{n}_a,N(i)=()^i\overline{N},`$
$`L_{ab}(i)=()^i\overline{L}_{ab}\eta (i)=\overline{\eta },\chi _a(i)=\overline{\chi }_a,`$ (23)
and solve it self-consistently at half-filling to obtain the phase diagram shown in Fig. 1. Higher symmetry lines E, F, G, H separate phases A, B, C, D as first order phase transition boundaries where order parameters smoothly rotate from one phase to another. Symmetries on lines E, F, G, H and the order parameters are SU(4)(adjoint Rep), SO(7)(vector Rep), SO(5)$``$ SU(2) (scalar $``$ vector Rep), SO(7) (adjoint Rep) as discussed before. Phases A and B spontaneously break the SO(5) symmetry in the adjoint and vector Rep channels respectively. Phases C and D have singlet pairing SC and CDW as order parameters, respectively. Order parameters in each phase and corresponding GS modes are summarized in Table 1. The effective theory is generally given by a quantum non-linear $`\sigma `$ model defined on the GS manifold.

One major difficulty of Monte-Carlo simulations in fermionic systems, the sign problemblankenbecler1981a , is absent in the spin 3/2 model when $`U_0U_23/5U_0`$. By the Hubbard-Stratonovich (HS) transformation, the partition function can be written as below when $`V=(3U_0+5U_2)/8>0`$ and $`W=(U_2U_0)/2>0`$, or equivalently $`U_0U_23/5U_0`$,

$`Z`$ $`=`$ $`{\displaystyle }Dn{\displaystyle }Dn^a\mathrm{exp}\{{\displaystyle \frac{V}{2}}{\displaystyle _0^\beta }d\tau {\displaystyle \underset{i}{}}n(i,\tau )^2`$
$``$ $`{\displaystyle \frac{W}{2}}{\displaystyle _0^\beta }d\tau {\displaystyle \underset{i,a}{}}n_a^2(i,\tau )\}det\{I+B\},`$
where $`B=𝒯e^{_0^\beta 𝑑\tau H_0+H_I(\tau )}`$ and $`𝒯`$ is the time order operator. Its discrete version is

$`B=e^{\mathrm{\Delta }\tau H_0}e^{\mathrm{\Delta }\tau H_I(\tau _L)}\mathrm{}e^{\mathrm{\Delta }\tau H_0}e^{\mathrm{\Delta }\tau H_I(\tau _2)}e^{\mathrm{\Delta }\tau H_0}e^{\mathrm{\Delta }\tau H_I(\tau _1)},`$
$`H_I(\tau )={\displaystyle \underset{i}{}}\psi _\alpha ^{}(i)\psi _\alpha (i)\left\{V(n(i,\tau )2)+(\mu \mu _0)\right\}`$
$`W{\displaystyle \underset{i,a}{}}\psi _\alpha ^{}(i)\mathrm{\Gamma }_{\alpha \beta }^a\psi _\beta (i)n^a(i,\tau ),`$ (24)
where $`\mathrm{\Delta }\tau =\beta /L`$. $`I+B`$ is invariant under the time-reversal transformation: $`T(I+B)T^1=I+B`$. If $`\lambda `$ is an eigenvalue of $`I+B`$ with the eigenvector $`|\varphi `$, then $`\lambda ^{}`$ is also an eigenvalue with the eigenvector $`T|\varphi `$. From $`T^2=1`$, it follows that $`\varphi |T\varphi =T^2\varphi |T\varphi =0`$, i.e. $`|\varphi `$ and $`T|\varphi `$ are orthogonal. Thus although $`I+B`$ may not be Hermitian because of the $`𝒯`$ operator, its determinant, a product of $`\lambda ^{}\lambda `$, is always positive semi-definite. Our proof is equally valid in the practical sampling with the discrete HS transformation as in Ref. motome1997a , and has been confirmed numericallypreparation . We emphasize that this proof is valid for any filling and lattice topology. A similar model has recently been introduced in Ref. assaad2002a , where the sign problem is also absent. However, their model keeps only the diagonal $`n_4^2`$ interaction and is not spin rotationally invariant. The valid region for the above algorithm (see Fig. 1) includes the 5-vector phases B, SP phase C and their SO(7) boundary, which are analogs of the competitions between antiferromagnetism and superconductivity in the high Tc context. It would be interesting to study the doping effect, the frustration on the triangular lattice, etc, which are difficult at low temperatures for previous Monte-Carlo works. Extensive numerical simulations are currently being carried outpreparation .

Besides the alkali atoms, the trapping and cooling of the alkaline-earth atoms are also exciting recentlymachholm2001 ; maruyama2003 . Among these two families, 132Cs, 9Be, 135Ba and 137Ba are spin 3/2 atoms. 132Cs is unstable and the $`2s^22s^12p^1`$ resonance of Be lies in the ultraviolet region, making them difficult for experimental use. The resonances of the last two Ba atoms are $`6s^26s^16p^1`$ at 553.7 nm he1991a , thus they are possible candidates. Their scattering lengths are not available now, but that of 138Ba (spin 0) was estimated as $`41a_B`$ machholm2001 , Because the $`6s`$ shell of Ba is full-filled, both the $`a_0`$, $`a_2`$ of 135Ba and 137Ba should have the similar value. Considering the rapid development in this field, we expect more and more spin 3/2 systems will be realized experimentally, allowing us to explore the full phase diagram.

In summary, we found an exact and generic SO(5) symmetry in spin 3/2 models with local interactions. This model can be accurately realized in cold atom systems and the theoretical predictions can be tested experimentally by the exact relationship among the Landau fermi liquid parameters, spin wave functions of the Cooper pairs, exact boundaries of quantum phase transitions among various competing states and the number of the collective modes. In the regime where accurate Monte Carlo simulations can be carried out without the sign problem, detailed quantitative comparisons with experiments are possible, including the quantum phase transition from the AF to the SC phases as a function of doping.

We thank D. P. Arovas, A. Auerbach, B. A. Bernevig, S. Capponi, H. D. Chen, C. Chin, A. L. Fetter, F. Kasevich, T. K. Ng for helpful discussions. We especially thank E. Mukamel for carefully reading the manuscript. This work is supported by the NSF under grant numbers DMR-9814289, and the US Department of Energy, Office of Basic Energy Sciences under contract DE-AC03-76SF00515. CW is also supported by the Stanford Graduate Fellowship program.


text:
# Sub-Poissonian statistics in order-to-chaos transition

## I Introduction

The subject of this article is at the boundary of two basic phenomena in quantum physics nowadays attracting significant interest of a broad readership. One of them is quantum chaos, which has been intensively investigated in last two decades and has been found in many nonlinear systems. Quantum effects inherent to the chaotic behavior of classical systems are an interesting field pertaining to many problems of fundamental interest . Another phenomenon is the generation of nonclassical states, including nonclassical states of light with a range of unique properties, e.g. sub-Poissonian statistics and squeezing. At present the problem of nonclassical states generation draws a lot of attention, both from the standpoints of pure knowledge and possible applications . The principal purpose of the present paper is to study the interesting crossovers between these phenomena for an open time-dependent quantum system.

The majority of studies of quantum chaos for isolated or so called Hamiltonian systems, the classical counterparts of which are chaotic, focus on static properties such as spectral statistics of energy levels and transition probabilities between eigenstates of the system. A variety of studies have also been carried out to understand the features of time-dependent chaotic systems. By contrast to that, very little work has been done to investigate the quantum chaos for open nonlinear systems. The beginning of study of an open chaotic system can be dated back to the papers , where the authors have analyzed the kicked rotor and similar systems with discrete time interacting with a heat bath. Quite generally, chaos in classical conservative and dissipative systems with noise has completely different properties, e.g., strange attractors can appear only in dissipative systems. Noise can induce a number of interesting phenomena in transition from regular to chaotic behavior. In particular, it was shown that under definite conditions noise can induce chaos as well as destroy chaotic behavior (see for example ). In recent years much effort has been expended, both theoretically and experimentally, to explore the role of quantum fluctuations and noise in order-to-chaos transition for open systems. It is obvious, that the investigations in this area are connected with the quantum-classical correspondence problem in general and with environment induced decoherence and dissipation in particular. Recently this topic has been the focus of theoretical investigations. As a part of these studies it has been recognized  that the decoherence has rather unique properties for systems classical analogues of which are chaotic. In particular, the formation of sub-Plank structure in phase space has been discussed for chaotic system . The connection between quantum and classical treatments of chaos was also realized by means of comparison between strange attractors on the classical Poincaré section and the contour plots of the Wigner functions . Some studies have explored the connection between quantum dynamical manifestation of chaos and quantum entanglement which is possibly the most typical property of composite quantum systems . From the experimental viewpoint, observation of dissipative effects and environment induced decoherence of dynamically localized states in the quantum delta-kicked rotor is carried out with the gas of ultracold caesium atoms in a magneto-optical trap subjected to a pulsed standing wave ,. Recently, new problems of chaotic motion have been studied in the experimental scheme with ultra-cold atoms in magneto-optical double-well potential . Fluctuation and decoherence in three-level model of chaos-assisted tunneling have been experimentally studied for samples of caesium atoms in an amplitude-modulated standing wave of light .

The other problem of interest in this area relates to the production of various nonclassical states in systems with chaotic dynamics. Generation of non-classical states of light has been the subject of an intense theoretical and experimental activity since the first observations of squeezed states of light. Generation of squeezed states in chaotic systems was studied numerically and analytically by using a specific form of perturbation theory for the model of quantum kicked rotator in Refs. . The result of numerical experiments on quadrature squeezing in simple quantum models that allow the transition to chaos in the classical limit are presented in . In a recent paper sub-Poissonian statistics of oscillatory excitations numbers is established for chaotic dynamics of nonlinear oscillator . The study of these phenomena provides a fundamental understanding of quantum fluctuations in quantum chaos and opens a way for new experimental studies of the quantum dissipative chaos in the field of Quantum Optics.

In addition to these important developments in investigations of quantum fluctuations and noise assisted to chaotic motion, in the present paper we examine mainly quantum statistics of elementary excitations for the model of nonlinear dissipative oscillator which shows the order-to-chaos transition. Our goal is twofold: first we study in detail what kind of statistics of oscillatory number states takes place in the quantum chaos, for both cases of vacuum- and temperature-reservoir, and second we analyze Wigner functions in a chaotic regime and strange attractors from the point of view of nonclassical statistics. It should be noted that sub-Poissonian statistics of oscillatory excitation numbers was established for chaotic dissipative dynamics and for the case of vacuum-reservoir in our previous paper . In the part of the present paper we essentially expand this study particularly considering the important details as well as the case of coupling of the nonlinear oscillator with nonzero temperature reservoir. The other part of the paper is devoted to the problems of generation of sub-Poissonian light, parameter scaling in chaotic dynamics and to the external noise-induced effects.

The requirement in realization of this study is to have a proper quantum model showing both regular and chaotic dynamics in the classical limit. We propose a nonlinear oscillator driven by two forces at different frequencies for this goal. This model was proposed to study the quantum chaos in one of the our previous papers , where it was shown that the model is apt to verification in experiments. This model also is attractive for the following reasons: (i) it is different from a single driven nonlinear oscillator and kicked rotor, where a pulsed pump field could be used for realization of a chaotic regime, and may be proposed for experimental studies of quantum chaos in the area of quantum optics with a cw laser; (ii) dynamics of this model exhibits a rich phase-space structure including regimes of regular and chaotic motion depending on the system’s parameters as well as displays remarkable quantum properties, in particular - nonclassical sub-Poissonian statistics of elementary excitations.

Open quantum systems are usually studied in the framework of reduced density matrix obtained by tracing over the degrees of freedom of environment. There have been suggestions (see, for example,) that chaotic dissipative dynamics should be described on individual quantum trajectories in the framework of quantum state diffusion method (QSD) . In addition to these studies, here we describe the quantum dissipative chaos using a statistical ensemble of trajectories which is usually realized in nature. Such description allows us to consider quantum manifestation of chaos in measurable quantities.

The plan of this paper is as follows. In Sec.II nonlinear oscillator driven by two forces at different frequencies is described and relevant regular and chaotic operational regimes are discussed. In Sec. III the formation of sub-Poissonian statistics of oscillatory excitations number for regular dynamics is analyzed and application to the generation of strong sub-Poissonian light is considered. The important novelty here is that the radial squeezing effect is much stronger, than an analogous one for the model of a single driven anharmonic oscillator. In Sec. IV we identify the kind of statistics of oscillatory number states which takes place for quantum chaos. Our central result here is that nonclassical, sub-Poissonian statistics can be realized for chaotic dynamics of the system under consideration. We also discuss the role of the dissipation in scaling relationships for the regular and chaotic motion as well as the dependence of results on the level of external noise. We summarize our results in Sec. V.

## II Order to chaos transition in doubly driven nonlinear oscillator

In this section we briefly review the results on the model . We treat dissipation and decoherence microscopically using a master equation which is solved numerically in the framework of QSD approach. This model is a nonlinear oscillator driven by two periodic monochromatic forces at different frequencies $`\omega _1`$ and $`\omega _2`$ for which the evolution of the reduced density operator $`\rho `$ is governed by the master equation

$$\frac{\rho }{t}=\frac{i}{\mathrm{}}\left[H_0+H_{int,}\rho \right]+\underset{i=1,2}{}\left(L_i\rho L_i^+\frac{1}{2}L_i^+L_i\rho \frac{1}{2}\rho L_i^+L_i\right).$$
(1)
We use the interaction picture corresponding to the transformation $`\rho e^{i\omega _1a^+at}\rho e^{i\omega _1a^+at}`$ and the Hamiltonians are:

$`H_0`$ $`=`$ $`\mathrm{}\mathrm{\Delta }a^+a,`$ (2)
$`H_{int}`$ $`=`$ $`\mathrm{}\left[\left(\mathrm{\Omega }_1+\mathrm{\Omega }_2\mathrm{exp}\left(i\delta t\right)\right)a^++\left(\mathrm{\Omega }_1^{}+\mathrm{\Omega }_2^{}\mathrm{exp}\left(i\delta t\right)\right)a\right]+\mathrm{}\chi (a^+a)^2.`$ (3)
Here$`\mathrm{\Delta }=\omega _0`$ $`\omega _1`$ is the detuning, and $`\delta =\omega _2`$ $`\omega _1`$ is the difference between driving frequencies, which works as modulation frequency. $`a,a^+`$ are boson annihilation and creation operators and $`L_i`$ are the Lindblad operators:

$$L_1=\sqrt{\left(N+1\right)\gamma }a,L_2=\sqrt{N\gamma }a^+,$$
(4)
where $`\gamma `$ is the spontaneous decay rate of the dissipation process, and $`N`$ denotes the mean number of quanta of a heat bath. The couplings with two driving forces are given by Rabi frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, and $`\chi `$ is the strength of anharmonicity. The last terms in Eq. (1) concern the influence of the environment induced diffusion. We have followed the standard approach ,  to dissipative quantum dynamics in the range of weak coupling of the system with the reservoir under the condition: $`\gamma k_BT/\mathrm{},`$ where $`k_BT`$ is the Boltzman’s constant times temperature. Equation (1) is obtained in both the rotating wave and Markov approximations, without regard to the driving-induced noise effects .

To study the pure quantum effects we focus on the cases of very low reservoir temperatures which, however, ought to be still larger than the characteristic temperature $`TT_{ch}=\mathrm{}\gamma /k_B.`$ This restriction implies that dissipative effects can be described selfconsistently in the frame of the Lindblad equation (1). Note here that for $`N1`$ and, even for vacuum reservoir, for $`N0`$ the above restriction is valid for the majority of problems of quantum optics.

For $`\mathrm{\Omega }_2=0`$ this equation describes the single driven, dissipative anharmonic oscillator, which is a well-known and archetypal model in nonlinear physics . In case of doubly driven oscillator ($`\mathrm{\Omega }_20`$), the interaction Hamiltonian (2) is explicitly time-dependent and the system exhibits regions of regular and chaotic motion. It should be noted that this model corresponds to a modified model of Duffing oscillator, i.e. an oscillator where the oscillatory frequency depends on the amplitude. To illustrate it, we first pay attention to the classical description.

### A Classical equation and connection to the Duffing oscillator

In the classical limit, the corresponding equation of motion for the dimensionless amplitude $`\alpha (t)=a(t)=Tr(a\rho (t))`$ has the form

$$\frac{d}{dt}\alpha =\frac{1}{2}\gamma \alpha i\left(\mathrm{\Delta }+\chi (1+2\left|\alpha \right|^2)\right)\alpha i\left(\mathrm{\Omega }_1+\mathrm{\Omega }_2\mathrm{exp}\left(i\delta t\right)\right).$$
(5)
It is easy to demonstrate that this equation is the rotating-wave approximation of the equation for Duffing oscillator driven by two periodic forces. Indeed, let us start with the following equation for a time-dependent real amplitude

$$\stackrel{}{E}+\gamma \stackrel{}{E}+\omega _0^2[1+\frac{2\chi }{\omega _0}(1+E^2/2)]E=4\omega _0(\mathrm{\Omega }_1\mathrm{cos}\omega _1t+\mathrm{\Omega }_2\mathrm{cos}\omega _2t).$$
(6)
which describes such oscillator. This equation is a generalization of the well-known Duffing oscillator to the case of two driving forces at different frequencies $`\omega _1`$ and $`\omega _2`$. To perform the rotating-wave approximation of the equation (6) we introduce the dimensionless complex amplitude $`\alpha (t)`$ as

$$E(t)=\alpha (t)\mathrm{exp}(i\omega _1t)+c.c.$$
(7)
Neglecting also the second-order derivative of the slowly varying amplitude $`\alpha (t)`$ and using the inequality $`\gamma \left|\dot{\alpha }\right|\omega _0^2\left|\alpha \right|`$ one obtains (5).

### B Operational regimes and measurable quantities

The dynamics of doubly driven anharmonic oscillator exhibits a rich phase-space structure, including regimes of regular and chaotic motion, with two Rabi frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, and the difference $`\delta `$ between driving frequencies and detuning $`\mathrm{\Delta }`$ being the control parameters. In study of the order-to-chaos transition in classical systems a useful tool is the examination of a constant phase map in the phase-space. In this context, our numerical analysis of the classical equation of motion in the $`(X,Y)`$ plane ($`X=Re\alpha `$, $`Y=Im\alpha )`$ shows that the classical dynamics of the system is regular in domains of small and large values of modulation frequency, i.e. $`\delta \gamma `$ and $`\delta \gamma `$, and also when one of the perturbation forces is much greater than the other: $`\mathrm{\Omega }_1\mathrm{\Omega }_2`$ or $`\mathrm{\Omega }_2\mathrm{\Omega }_1.`$ The dynamics is chaotic in the range of parameters $`\delta \gamma `$ and $`\mathrm{\Omega }_1\mathrm{\Omega }_2`$, where the classical results for $`X=Re\alpha ,Y=Im\alpha `$ obtained from Eq.(5) show that the classical strange attractors for Poincaré section are realized . Thus, two ways of producing chaos in a controlled manner can be considered by monitoring the system through varying either the strength of the driving force, or the difference frequency $`\delta `$.

We will examine macroscopic quantum effects contributing to the chaotic behavior by considering both the Mandel $`Q`$ parameter  which describes the deviation of excitation number uncertainty from the Poissonian variance, i. e. $`Q=\left(\left(\mathrm{\Delta }n\right)^2n\right)/n,`$ $`\left(\mathrm{\Delta }n\right)^2=\left(a^+a\right)^2a^+a^2`$, and the Wigner function. Thus, we suggest that the distinction between regular and chaotic dynamics can be most easily understood by studying the dynamics of essentially quantum properties. In fact, our numerical analysis has shown that the quantum dynamical manifestation of chaotic behavior does not appear on ensemble averaged oscillatory excitation numbers, but is clearly seen on the probability distributions . The connection between quantum and classical treatments of chaos can be realized by means of comparison between strange attractors in the classical Poincaré section, and the contour plots of the Wigner functions. Indeed, it was demonstrated in details  that for comparatively small values of the ratio $`\chi /\gamma 0.1`$ the contour plots of Wigner functions are relatively close to the strange attractors. Such likeness of quantum and classical distribution vanishes in the deep quantum regime.

## III Sub-Poissonian statistics for regular dynamics

Our numerical analysis of quantum effects assisting to chaotic dynamics is based on QSD approach that represents the reduced density operator averaging over the projectors onto the stochastic states $`|\mathrm{\Psi }_\xi `$ of the ensemble: $`\rho (t)=M\left(|\mathrm{\Psi }_\xi \mathrm{\Psi }_\xi |\right)`$, where $`M`$ denotes the ensemble averaging. The corresponding equation of motion is

$`|d\mathrm{\Psi }_\xi `$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}H|\mathrm{\Psi }_\xi dt`$ (9)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1,2}{}}\left(L_i^+L_i2L_i^+L_i+L_iL_i^+\right)|\mathrm{\Psi }_\xi dt+{\displaystyle \underset{i=1,2}{}}\left(L_iL_i\right)|\mathrm{\Psi }_\xi d\xi _i,`$
where $`\xi `$ indicates the dependence on the stochastic process, the complex Wiener variables $`d\xi _i`$ satisfy the fundamental properties $`M\left(d\xi _i\right)=0,M\left(d\xi _id\xi _j\right)=0,M\left(d\xi _id\xi _j^{}\right)=\delta _{ij}dt,`$ and the expectation value $`L_i=\mathrm{\Psi }_\xi \left|L_i\right|\mathrm{\Psi }_\xi `$. For the oscillatory mean excitations number and the variance this method gives:

$$n=\underset{m\mathrm{}}{lim}\frac{1}{m}\underset{\xi }{\overset{m}{}}\mathrm{\Psi }_\xi \left|a^+a\right|\mathrm{\Psi }_\xi ,$$
(10)
$$\left(\mathrm{\Delta }n\right)^2=\underset{m\mathrm{}}{lim}\frac{1}{m}\underset{\xi }{\overset{m}{}}\left(\mathrm{\Psi }_\xi \left|(a^+a)^2\right|\mathrm{\Psi }_\xi \left|\mathrm{\Psi }_\xi \left|(a^+a)\right|\mathrm{\Psi }_\xi \right|^2\right),$$
(11)
while the Wigner function is calculated with the formulas

$$W(r,\theta )=\underset{m,n}{}\rho _{nm}W_{mn}(r,\theta ),$$
(12)
in terms of matrix elements $`\rho _{nm}=M\left(n|\mathrm{\Psi }_\xi \mathrm{\Psi }_\xi |m\right)`$ of density operator in the Fock state representation, where $`(r,\theta )`$ are the polar coordinates in the complex phase-space plane $`X=r\mathrm{cos}\theta `$, $`Y=r\mathrm{sin}\theta .`$ The coefficients $`W_{mn}(r,\theta )`$ are Fourier transform of matrix elements of the Wigner characterization function

$$W_{mn}(r,\theta )=\left\{\genfrac{}{}{0pt}{}{\frac{2}{\pi }(1)^n\sqrt{\frac{n!}{m!}}e^{i(mn)\theta }(2r)^{mn}e^{2r^2}L_n^{mn}(4r^2),mn}{\frac{2}{\pi }(1)^m\sqrt{\frac{m!}{n!}}e^{i(mn)\theta }(2r)^{nm}e^{2r^2}L_m^{nm}(4r^2),nm}\right\},$$
(13)
where $`L_p^q`$ are the Laguerre polynomials.

### A Controlled strong sub-Poissonian statistics

First we discuss the case of classically regular behavior of the considered model assuming the interaction of the system with vacuum reservoir, $`N=0`$. In our numerical analysis we focus, for simplicity, on the case of strong anharmonicity considering the parameters from $`\chi /\gamma =0.05`$ to $`\chi /\gamma =0.7`$. In the latter case the maximum mean number of oscillatory excitations $`n=15`$, while in the former case it is equal to $`n=250.`$

The results for time evolution of the mean oscillatory excitation number as well as the Mandel $`Q`$ parameter averaged over quantum trajectories, are depicted in Fig.1 and 2 for the parameters: (curve 1) $`\chi /\gamma =0.7`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =10.2`$, $`\mathrm{\Omega }_2/\gamma =13.5`$, and $`\delta /\gamma =5`$; (curve 2) $`\chi /\gamma =0.05`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =38.18`$, $`\mathrm{\Omega }_2/\gamma =49.5`$, and $`\delta /\gamma =5`$. As we see the $`Q`$ parameter shows a time-dependent modulation with a period $`2\pi /\delta `$ and formation of sub-Poissonian statistics, ( $`Q<0`$) in which the oscillatory excitation number fluctuations are squeezed below the coherent level $`\left(\mathrm{\Delta }n\right)^2<n),`$ for time intervals exceeding the transient time, $`t\gamma 5.`$ The Mandel parameter reaches its minimum $`Q_{\mathrm{min}}0.9`$ and maximum $`Q_{\mathrm{max}}0.38`$ values for definite time intervals, for the case illustrated in Fig. 2 (curve 2). Moreover, the degree of sub-Poissonian statistics increases with increasing of mean oscillatory excitation numbers at definite time intervals. Below we demonstrate the formation of sub-Poissonian statistics by analyzing the Wigner function. Figs. 3 show the Wigner function and its contour-plot at the fixed moments of time $`t_k=[7.13+(2\pi /\delta )k]\gamma ^1,`$ $`(k=0,1,2,\mathrm{})`$ exceeding transient time. We find that the Wigner function is located around the point $`X=0,`$ $`Y=10`$ and its contour-plot has narrow crescent form with the origin of phase space as its centrum. Another peculiarity is that the Wigner function is nonstationary. As calculations show, during the period of modulation $`2\pi /\delta `$ it revolves around the origin of the phase-space. The radial squeezing that represents the known property of the anharmonic oscillator to produce the excitation number squeezing is also clearly seen in the figure. The important novelty here is that the radial squeezing effect in this model is much stronger, than an analogous one for the model of single driven anharmonic oscillator. Indeed, for $`\mathrm{\Omega }_2=0`$, when we turn to the single driven anharmonic oscillator, the Mandel parameter is $`Q=0.72`$ for the same parameters as in Fig.2 (curve 2).

### B Application to the generation of nonclassical light

This model can be proposed for the realization of optical schemes generating light with strong sub-Poissonian statistics. One of the candidates may be a scheme involving self-phase modulation in an optical cavity pumped by two coherent laser fields. In this scheme the anharmonicity of the oscillatory mode comes from photon-photon interactions in the $`\chi ^{(3)}`$ medium inserted in the cavity, while dissipative effects arise from leakage of photons through the cavity mirrors, which damps the mode. In the case of a ring cavity the evolution of the mode at frequency $`\omega _0`$ is described by the Hamiltonian (2) with operators $`a,`$ $`a^+`$ being the annihilation and creation operators. The coupling constant $`\chi `$ of the Kerr interaction is proportional to the third-order nonlinear susceptibility $`\chi ^{(3)},`$ and $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are the Rabi frequencies corresponding to two classical coherent fields at the frequencies $`\omega _1`$ and $`\omega _2,`$ respectively. The cavity damping and noise effect are described by the Lindblad operator (4), where $`\gamma `$ is the damping rate and $`N`$ is the mean number of thermal reservoir photons.

For the concrete scheme of light generation in a ring cavity, when the coupling of ”in” and ”out” fields occurs at the different mirrors, the cavity-output field operator $`a_{out}=\sqrt{2\gamma }a`$, while the cavity-output intensity in photon number units per unit time is determined by $`n_{out}=a_{out}^+a_{out}=2\gamma a^+a.`$

It is easy to verify that in this case the results depicted on the Fig.1 relate also to the normalized intensity $`n_{out}/2\gamma `$ of output light. For this scheme the measure of photocount statistics of detected output photons can also be directly expressed in terms of itracavity photons, namely Mandel factor $`Q.`$ Indeed, in sub-Poissonian light experiments the observable quantity is the deviation of the variance of photocount number of output photons defined as $`Q_i=\left(\left(\delta n_i\right)^2n_i\right)/n_i`$, where $`n_i`$ is the mean number of photocounts . We consider a short counting time $`T2\pi /\delta `$ to be accurate in detection of the time-dependent output field modulated with period $`2\pi /\delta `$. In this approximation $`Q_i=\alpha Tn_{out}\left(g_{out}^{(2)}1\right)`$, where $`\alpha `$ is a dimensionless quantum output of a detector, and $`g_{out}^{(2)}=a_{out}^+a_{out}^+a_{out}a_{out}/n_{out}^2`$ is the second order correlation function of the output field. Further, in terms of intracavity photon numbers $`Q_i=2\gamma \alpha TQ`$, and hence the normalized quantity $`Q_i/2\gamma \alpha T`$ reduces to the Mandel parameter for intracavity photons. Thus, the results of Fig.2 are applicable to quantum statistics of light and indicate the time-dependence of the normalized quantity $`Q_i/2\gamma \alpha T`$. In particular, the most strong sub-Poissonian statistics of generated light takes place for time intervals at which the output intensity reaches its maximal values.

It is natural to explain such strong sub-Poissonian statistics of light by a nontrivial interference effect or by the mutual interference of two coherent fields in nonlinear medium. However such analysis is difficult to perform because the system cannot be analyzed analytically in details. In this context we restrict ourselves by the remark that the system under consideration is an example of a quantum systems in contact with environment, where the controlling and ordering of dissipative dynamics as well as quantum statistics are realized through an external time-dependent force. Some examples have been proposed in Refs. .

## IV Nonclassical statistics for chaotic dynamics

We now study the emergence of quantum chaos, which is expected to manifest itself as crucial changes in the above results in the classically chaotic operational regime. One of the ways to realize the controlling transition from regular to chaotic dynamics is to vary the strength $`\mathrm{\Omega }_2`$ of the second force in the range from $`\mathrm{\Omega }_2\mathrm{\Omega }_1`$ to $`\mathrm{\Omega }_2\mathrm{\Omega }_1`$. In the limit $`\mathrm{\Omega }_2\mathrm{\Omega }_1`$ the system is reduced effectively to the model of single driven anharmonic oscillator, which exhibits bistability for a definite range of parameters $`\chi ,\mathrm{\Delta },\mathrm{\Omega }_1,`$ and $`\gamma `$  On increasing the amplitude $`\mathrm{\Omega }_2`$ the system transits from the regular to the chaotic regimes at $`\mathrm{\Omega }_1\mathrm{\Omega }_2.`$ In the limit $`\mathrm{\Omega }_2>\mathrm{\Omega }_{1\text{ }}`$the regular dynamics is restored. The results of numerical simulation of the Mandel parameter in the regime of chaos are shown on Figs.4 for the parameters: (curve 1) $`\chi /\gamma =0.7`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =\mathrm{\Omega }_2/\gamma =10.2`$, and $`\delta /\gamma =5`$; (curve 2) $`\chi /\gamma =0.05`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =\mathrm{\Omega }_2/\gamma =38.18`$, and $`\delta /\gamma =5.`$ Surprisingly, the excitation-number fluctuations are also squeezed below the coherent level for the considered chaotic regime. However, in contrast to the previous regular regime, the excitation number exhibits both sub-Poissonian ($`Q<0`$) and super-Poissonian ($`Q>0`$) statistics, that alternate in definite time intervals. The minimum and maximum values of $`Q`$ in time intervals during one modulation period are equal to $`Q_{\mathrm{min}}0.7`$ and $`Q_{\mathrm{max}}0.6`$ \[for case (b)\]. Thus, Figs.2 and 4 show the drastic difference between the behavior of Mandel parameter for regular and chaotic dynamics. It means that the variance of oscillatory number fluctuations may be used for testing of quantum chaos.

Now we illustrate the emergence of nonclassical sub-Poissonian statistics in the doubly driven nonlinear oscillator in its transition from regular to chaotic dynamics using the phase-space symmetry properties of the Wigner function. One can make sure of that by comparing the contour-plots of Wigner function for sub-Poissonian and super-Poissonian statistics. The results of ensemble-averaged numerical calculations of both the Wigner function and its contour-plot at fixed time intervals $`t_k=[6.96+(2\pi /\delta )k]\gamma ^1,`$ $`(k=0,1,2,\mathrm{})`$ are shown in Figs. 5 (a, b) respectively. As we see the contour-plot for chaotic motion still has the radial squeezed form \[see Fig. 5 (b)\]. This result takes place for $`t_k=[6.96+(2\pi /\delta )k]\gamma ^1,`$ $`(k=0,1,2,\mathrm{})`$ at which the Mandel parameter reaches its minimum value $`Q_{\mathrm{min}}0.7`$. In subsequent time intervals during the period of modulation the level of excitation number fluctuations increases, and as a result the radial squeezing in contour-plot decreases. This result is depicted in Figs. 6 (a) for the same parameters as in Fig.5(a), but for time-intervals $`t_k=[7.325+(2\pi /\delta )k]\gamma ^1,`$ $`(k=0,1,2,\mathrm{})`$. The Poincaré section corresponding to these time moments are shown in Fig. 6 (b).

As we see there is a possibility to control the statistics of excitation numbers of the oscillator by means of variation of the field strengths $`\mathrm{\Omega }_2`$. One can obtain the following regime of generation: purely sub-Poissonian statistics (Fig. 2 (curves 1, 2)), oscillating between super- and sub-Poissonian (Fig. 4 (curves 1, 2)). Moreover, for certain values of the parameters we achieve an essential improvement of sub-Poissonian statistics in comparison with a single driven anharmonic oscillator, $`\mathrm{\Omega }_2=0`$. In this context it seems important to analyze in more detail the Mandel parameter in order-to-chaos transition. To this end we have also studied the behavior of the $`Q`$ factor versus the controlling parameter $`\mathrm{\Omega }_2/\gamma `$. The results of numerical calculations at a definite time moments exceeding the transient time, when the minimum of Mandel parameter was realized, are displayed in Fig.7. As is seen the result is displayed in the complicated form and shows the essential improvement of sub-Poissonian statistics in the chaos-to-order transition taking place for high values of $`\mathrm{\Omega }_2/\gamma `$. The sub-Poissonian statistics goes bad for small values of Rabi frequency $`\mathrm{\Omega }_2/\gamma `$ as well as for its intermediate values corresponding to the chaotic regime.

There is a particularly interesting feature to notice in the comparative analysis of strange attractors and contour plots of the Wigner function \[Figs.5 (b,c) and Fig. 6 (a,b)\]. We can conclude that it is possible to predict the shape of Wigner function by the knowledge of Poincaré section. Therefore, it is possible to predict the squeezing in quantum system by the knowledge of Poincaré section of its chaotic classical counterpart, that occurs by using Fig. 5 (c). This opens new possibilities to guess intuitively what kind of states will appear in quantum chaotic system by treating its classical chaotic map.

### A Parameter scaling and dependence on external noise

Further investigation of this model allows to establish other properties of dissipative quantum chaos. It is easy to verify that Eq.(5) is invariant with respect to the following scaling transformation of the complex amplitude $`\alpha \alpha ^{}=\lambda \alpha ,`$ where $`\lambda `$ is a real positive dimensionless coefficient, if the parameters $`\chi ,`$ $`\mathrm{\Delta }`$, $`\mathrm{\Omega }_1,`$ $`\mathrm{\Omega }_2`$ are correspondingly transformed as: $`\chi \chi ^{}=\chi \mathrm{}\lambda ^2,`$ $`\mathrm{\Delta }\mathrm{\Delta }^{}=\mathrm{\Delta }+\chi \left(11\mathrm{}\lambda ^2\right),`$ $`\mathrm{\Omega }_{1,2}\mathrm{\Omega }_{1,2}^{}=\lambda \mathrm{\Omega }_{1,2}.`$ As we see the scaling transformations do not include the rate of the dissipative process $`\gamma `$ and hence can be written by using of dimensionless parameters as the follows: $`\mathrm{\Delta }/\gamma \mathrm{\Delta }^{}/\gamma ^{}=\mathrm{\Delta }/\gamma +(\chi /\gamma )\left(11\mathrm{}\lambda ^2\right)`$, $`\chi /\gamma \chi ^{}/\gamma ^{}=\chi \mathrm{}\gamma \lambda ^2`$, $`\mathrm{\Omega }_{1,2}/\gamma \mathrm{\Omega }_{1,2}^{}/\gamma ^{}=\lambda \mathrm{\Omega }_{1,2}/\gamma `$ . This scaling property of the classical equation for the chaotic dynamics leads to the symmetry of strange attractors which, for definite different sets of parameters, have the same form in the phase-space and differ from each other only in scale. Some examples of such behavior have been presented in our previous paper . So, we conclude that the parameter scaling happens for the classical trajectory dynamics for both regular and chaotic regimes. It is interesting to analyze such scaling from the point of view of quantum-statistical theory and its relevance for decoherence and dissipation. It should be noted, that recently several studies have been devoted to the scaling relationships in the decoherent quantum-classical transition for chaotic systems (see  and reference therein). The quantum-classical transition is now understood to be affected by the relative size of $`\mathrm{}`$ for a given system as well as by the parameter $`D`$, a measure of decoherence of quantum system of interest. Further, in systems where the classical equation is chaotic, the transition is also affected by the chaos in the system, and thus by $`\theta `$, the Lyapunov exponent of the classical trajectory dynamics. Indeed, above mentioned study  was devoted to scaling relationships involving $`\mathrm{},`$ $`D,`$ $`\theta `$. Here we hold on another approach and present the results on scaling relationships, involving the system parameters $`\chi /\gamma `$, $`\mathrm{\Delta }/\gamma `$, and $`\mathrm{\Omega }_{1,2}/\gamma `$, in the framework of oscillatory numbers and its statistics. The goal is to analyze what kind of parameter scaling we can achieve in a quantum ensemble theory in the presence of quantum noise. We show below that such parameter scaling occurs for wider ranges of the parameters, but for not large values of the parameter $`\chi /\gamma `$. As our numerical calculations show $`Q`$ Mandel factor is approximately invariant under scaling transformations, i.e. $`Q(\chi /\gamma ,\mathrm{\Delta }/\gamma ,\mathrm{\Omega }_{1,2}/\gamma )Q(\chi /\lambda ^2\gamma ,\mathrm{\Delta }/\gamma +\chi /\gamma (11/\lambda ^2),\lambda \mathrm{\Omega }_{1,2}/\gamma )`$. This is illustrated in Fig.8 (a) for a time interval during the period of modulation and for the regular dynamics. Two sets of the parameters are: (curve 1) $`\chi /\gamma =0.1`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =27,`$ $`\mathrm{\Omega }_2/\gamma =35`$, and $`\delta /\gamma =5`$; (curve 2) $`\chi /\gamma =0.05`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =38.18,`$ $`\mathrm{\Omega }_2/\gamma =49.5`$, and $`\delta /\gamma =5`$ with scaling parameter $`\lambda =\sqrt{2}`$. As we see, the scaling invariance is valid within $`0.02`$. We also easily observe that the degree of violation of the scaling increases with increasing of the ratio $`\chi /\gamma `$ due to the moving of the system to a deeper quantum regime. See, in particular Fig. 2 and Fig. 4 where we have shown the results for two sets of parameters which are connected with a scaling transformation with $`\lambda =\sqrt{14}`$. Analogous results take place for the ranges of the parameters leading to chaotic dynamics. In Fig. 8(b) we illustrate the scaling of the $`Q`$ factor for the parameters: (curve 1) $`\chi /\gamma =0.1`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =`$ $`\mathrm{\Omega }_2/\gamma =27,`$ and $`\delta /\gamma =5`$; (curve 2) $`\chi /\gamma =0.05`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =`$ $`\mathrm{\Omega }_2/\gamma =38.18,`$ and $`\delta /\gamma =5`$ with scaling parameter $`\lambda =\sqrt{2}.`$ The breaking of scaling invariance in deep quantum chaotic regime is shown in Fig. 4 with scaling parameter $`\lambda =\sqrt{14}.`$We emphasize that even though we have motivated the study of the parameter scaling with connection to dissipation and decoherence, the results obtained can be useful for other problems.

In the preceding sections we have studied statistics of excitation number in the pure quantum regime assuming coupling of nonlinear oscillator with a reservoir at zero temperature. In this part of the section we address the question of how external noise affects statistics of excitation numbers. We calculate the Mandel parameter with formulas (10), (11) as in the previous case of vacuum reservoir. Calculating quantum trajectories, we set the system initially to the vacuum oscillatory state and integrate Eq. (9) with Lindblad operators (4) for series of definite values of $`N`$. The results of numerical calculations at definite time moments exceeding the transient time are depicted in Figs.9 , where minimum values of the $`Q`$ factor in time during the period $`2\pi /\delta `$ versus external noise strength are presented. The curve (1) describes the case of a single driven anharmonic oscillator ($`\mathrm{\Omega }_2=0`$) for the same parameters as for the curve (2). As we see for the regular regime (curve(2)) $`Q`$ increases linearly with $`N`$. The Mandel parameter is $`Q=0.88`$ for $`T=0`$ K ($`N=0`$), but is 30 % larger for $`N=0.5`$ (curve 2). An analogous linear dependence takes place for the single driven anharmonic oscillator . However, for noisy chaotic regime the $`Q`$ factor is approximately independent of noise and even slightly decreases for comparatively large ranges of $`N`$ \[Fig. 9 (curve 3)\]. Such behavior can be explained by additional noise-induced mechanisms in chaotic evolution, namely, by chaos-to-order transition under influence of external noise.

## V Conclusion

Thus, in this work we have numerically studied the phenomena at the overlap of chaos, dissipation and nonclassical statistics for the time-dependent nonlinear model, namely Duffing oscillator driven by two periodic forces, showing order-to-chaos transition. We have shown that physical systems based on this model have a potential for generation of high-degree sub-Poissonian light as well as for observation of quantum-statistical effects that accompany chaotic dynamics. We have analyzed the mean oscillatory excitation number, Mandel $`Q`$ factor and the Wigner function for both regular and chaotic dynamics. We have shown in detail that the statistics of oscillatory excitation number is drastically changed in the order-to-chaos transition. Particularly interesting is the essential improvement of sub-Poissonian statistics in the regular operational regime of the scheme outlined above in comparison with an analogous one for the standard model of a single driven anharmonic oscillator. We have observed that in the chaotic regime the system exhibits ranges of different statistics depending on the modulation time, which alternate between sub- and super-Poissonian statistics. Another important conclusion concerns the conformity between strange attractors and contour plots of the Wigner function, that may be used for qualitative testing of quantum states by treating its classical chaotic map. Considering noise-induced effects we conclude that $`Q`$ parameter increases with increasing of the noise level $`N`$ for regular dynamics. This result is in good agreement with the results obtained for the single driven anharmonic oscillator . However, $`Q`$ parameter is almost insensitive to the level of external noise for the case of chaos.

Acknowledgments

G. Yu. K. gratefully acknowledges useful discussions with H.K. Avetisian and support by the ISTC grant No. A-353. S. B. M. thanks Zh. Haroutunyan for technical support.

FIGURE CAPTIONS

Fig. 1. Excitation numbers in regular regime for two sets of parameters (see text). Excitation numbers for curve 1 is multiplied on 14. The averaging is over 2000 trajectories.

Fig. 2. The Mandel factor for the same as in Fig.1 parameters. The averaging is over 2000 trajectories.

Fig. 3. The Wigner function (a) and its contour-plot (b) for the regular regime averaged over 3000 trajectories.

Fig. 4. The Mandel factor in chaotic regime for two sets of parameters (see text).

Fig. 5. The Wigner function (a) and its contour-plot (b) averaged over 3000 trajectories. The Poincaré section (c) (approximately 20000 points) plotted at times of the constant phase, when the maximal sub-Poissonian statistics in quantum regime is reaslized. In both cases the dimensionless parameters are in the range of chaos:

Fig. 6. The contour-plot of Wigner function (a) and corresponding Poincaré section (b) (approximately 10000 points) constructed at time moments when the maximal super-Poissonian statistics is reaslized. The parameters and averagind for case (a) are the same as in Fig. 5.

Fig. 7. The behavior of the minimums of the Mandel factor $`Q_{\mathrm{min}}`$ versus the controlling parameter $`\mathrm{\Omega }_2/\gamma `$. The parameters are: $`\chi /\gamma =0.1,`$ $`\mathrm{\Delta }/\gamma =15,`$ $`\mathrm{\Omega }_1/\gamma =27,`$ $`\delta /\gamma =5.`$

Fig. 8. The scaling of the $`Q`$ factor for two sets of parameters (see text) in regular (a) and chaotic (b) regimes. The averaging is over 1000 trajectories.

Fig. 9. The dependence of $`Q_{\mathrm{min}}`$ on the noise intensity for: (curve 1) single driven nonlinear oscillator ($`\chi /\gamma =0.1`$, $`\mathrm{\Delta }/\gamma =15`$, $`\mathrm{\Omega }_1/\gamma =`$ $`27,`$ and $`\delta /\gamma =5`$), (curve 2) doubly driven nonlinear oscillator in regular ($`\mathrm{\Omega }_2/\gamma =35`$), and (curve 3) in chaotic ($`\mathrm{\Omega }_2/\gamma =`$ $`\mathrm{\Omega }_1/\gamma =27`$) regimes.