text:
# The Giant X-Ray Flare of NGC 5905: Tidal Disruption of a Star, a Brown Dwarf, or a Planet?

## 1 Introduction

Tidal disruption of a star by a supermassive black hole in a galactic nucleus has been investigated by many authors. Early seminal investigations by Hills (1975); Lacy, Townes, & Hollenbach (1982); Rees (1988); Evans & Kochanek (1989) and Cannizzo, Lee, & Goodman (1990) were followed up with detailed analytical studies (Kochanek, 1994; Khokhlov, & Melia, 1996; Loeb & Ulmer, 1997; Ulmer, Paczyński, & Goodman, 1998; Ulmer, 1999) and numerical simulations (Khokhlov, Novikov, & Pethick, 1993; Frolov et al., 1994; Marck, Lioure, & Bonazzola, 1996; Diener et al., 1997; Kim, Park, & Lee, 1999; Ayal, Livio, & Piran, 2000; Ivanov & Novikov, 2001).

A tidal disruption event is expected to produce a luminous flare of electromagnetic radiation in the UV to X-ray band. Komossa (2001) argued that a convincing detection of a tidal disruption event must fulfill the following three criteria: (1) The event should be of short duration (a “flare”); (2) It should be very luminous (up to $`L_{\mathrm{max}}10^{45}`$ erg s-1 at maximum); (3) It should reside in a galaxy that is otherwise non-active (in order to exclude an upward fluctuation in the mass accretion rate of an active nucleus).

ROSAT observations have led to the discovery of four X-ray flaring “normal” galaxies that fulfill the above criteria \[summarized in Komossa (2001)\]: NGC 5905, RXJ1242-1119, RXJ1624+7554, RXJ1420+5334; plus a possible fifth candidate: RXJ1331-3243. All these galaxies have exhibited large X-ray flares (peak luminosity up to $`10^{44}`$ erg s-1), corresponding to a significant amplitude of variability (up to a factor $`200`$), an ultra-soft X-ray spectrum ($`k_BT_{\mathrm{bb}}0.040.1`$ keV, where $`T_{\mathrm{bb}}`$ is the temperature derived from a black-body fit to the spectral data), and no signs of Seyfert-like activity in the optical. These flares are promising candidates for tidal disruption events (Komossa, 2001). A UV flare has also been detected at the center of a mildly active elliptical galaxy, NGC 4552, with a bolometric luminosity $`10^{39}`$ erg s-1. Since the luminosity is several orders of magnitude less than would be expected from complete tidal disruption of a star, the event has been explained as the tidal stripping of a stellar atmosphere (Renzini et al., 1995).

Among the five galaxies listed above, long-term follow-up observations exist only for NGC 5905 (Komossa & Bade, 1999). In this paper, we model the light curve of NGC 5905 assuming that the flare was produced by the tidal disruption of a star by a supermassive black hole. We compare theoretical predictions with the observational data and constrain the properties of the progenitor star (e.g., its mass and radius) and estimate how much mass was accreted by the black hole.

## 2 The Light Curve of the Flare in NGC 5905

A giant X-ray flare in NGC 5905 was detected by ROSAT in July 1990 (Bade, Komossa, & Dahlem, 1996). The observed peak X-ray luminosity exceeded $`10^{42}`$ erg s-1, and the X-ray spectrum was very soft (photon index $`\mathrm{\Gamma }_X4`$, or a very cool black-body with $`k_BT_{\mathrm{bb}}0.06`$ keV). Follow-up observations were made in 1990 and 1993 with the PSPC, and in 1996 with the HRI (Komossa & Bade, 1999). The observational data are listed in Table 1. Assuming a Hubble constant $`H_0=50`$ km s-1 Mpc-1, the count rates may be converted to luminosities, which are plotted with error bars in Figure 1. The two downward arrows are upper limits. The light curve after the initial burst gradually declines with time. During the X-ray outburst as well as at later times, the optical brightness is constant to within observational errors ($`0.2^m`$).

Fitting the late decline of the NGC 5905 light curve (the four last data points in Table 1) with a $`(tt_D)^{5/3}`$ law that is predicted by the fallback model of tidal disruption (see §3), and treating the largest luminosity observed at $`t=1990.54`$ yr, viz. $`L_X=4.47\times 10^{42}\mathrm{erg}\mathrm{s}^1`$, as a lower limit (since this was still in the rising part of the light curve), we obtain

$`L_X=(0.30\pm 0.03)\times 10^{42}\mathrm{erg}\mathrm{s}^1\left[{\displaystyle \frac{t(1990.36\pm 0.02)\mathrm{yr}}{1\mathrm{yr}}}\right]^{5/3},`$ (1)
where the time $`t`$ is in units of years, and $`t_D=1990.36\pm 0.02`$ yr is the moment when tidal disruption happened.

During the first epoch of observations at time $`t=1990.54`$, the X-ray flux rose steadily (see Table 1), strongly suggesting that this epoch corresponded to the onset of the flare. For a start time of $`t_11990.54`$ yr, the total X-ray energy released during the entire outburst event is

$`\mathrm{\Delta }E_X={\displaystyle _{t_1}^{\mathrm{}}}L_X(t)𝑑t(4.5\pm 0.9)\times 10^{49}\mathrm{erg}.`$ (2)
A black-body fit to the X-ray spectrum during the outburst gives a temperature $`k_BT_{\mathrm{bb}}0.06`$ keV, and an intrinsic ($`0.12.4`$ keV) mean luminosity $`L_X3\times 10^{42}`$ erg s-1. Integrating the black-body spectrum over all energies (from zero to infinity), the bolometric luminosity is $`L1.1L_X`$. Thus, if the spectrum does indeed have a black-body form, the bolometric correction is not very significant. Then, the total mass associated with the released energy may be estimated to be

$`\mathrm{\Delta }M={\displaystyle \frac{\mathrm{\Delta }E}{ϵc^2}}{\displaystyle \frac{\mathrm{\Delta }E_X}{ϵc^2}}(2.5\pm 0.5)\times 10^4M_{}\left({\displaystyle \frac{ϵ}{0.1}}\right)^1,`$ (3)
where $`ϵ`$ is the efficiency of converting mass to radiated energy. This mass estimate is similar to that for the UV flare in NGC 4552 (Renzini et al., 1995). Note that our estimate of $`\mathrm{\Delta }E_X`$ depends only weakly on the assumed power law index $`5/3`$ in equation (1). This is because the observed light curve extends over a long enough time to capture most of the emission in the flare. The estimate of $`\mathrm{\Delta }M`$ suffers from the uncertainty in the bolometric correction (which we have taken to be unity).

In models of thin accretion disks around black holes, the color temperature $`T_{\mathrm{bb}}`$ of the emitted radiation is generally higher than the effective temperature $`T_{\mathrm{eff}}`$, defined by flux $`F=\sigma T_{\mathrm{eff}}^4`$. The deviation of $`T_{\mathrm{bb}}`$ from $`T_{\mathrm{eff}}`$ is caused by electron scattering and Comptonization (Ross, Fabian & Mineshige 1992; Shimura & Takahara 1993, 1995a, 1995b). Writing $`T_{\mathrm{bb}}=f_cT_{\mathrm{eff}}`$, the factor $`f_c`$ can be estimated from detailed self-consistent models of the radiative transfer in the disk. It is found that the value of $`f_c`$ depends on both the mass of the black hole (compare Shimura & Takahara 1993 and 1995b) and the mass accretion rate. Shimura & Takahara (1993) present results for a $`10^8M_{}`$ black hole emitting at a radius of 5 Schwarzschild radii, which is similar to the conditions we expect in our problem. Depending on the accretion rate, they obtain values of $`f_c`$ in the range $`1.33.2`$ (see their Fig. 4 and the associated discussion), with the highest value of $`f_c`$ being obtained for near-Eddington rate accretion. In what follows, we retain the factor $`f_c`$ in the equations, and we substitute a value of 3 whenever we require a numerical estimate.

Including the factor $`f_c`$ as defined above, we estimate the radius of the X-ray emitting region of the flare in NGC 5905 to be

$`R_X=\left({\displaystyle \frac{f_c^4L_X}{\pi \sigma T_{\mathrm{bb}}^4}}\right)^{1/2}2.4\times 10^{12}\mathrm{cm}\left({\displaystyle \frac{f_c}{3}}\right)^2,f_c{\displaystyle \frac{T_{\mathrm{bb}}}{T_{\mathrm{eff}}}},`$ (4)
where we have used $`L_X3\times 10^{42}`$ erg s-1 for the mean luminosity during the peak of the outburst (Komossa, 2001). The radius we derive is surprisingly small; even after setting $`f_c=3`$, it is barely equal to the Schwarzschild radius of a $`10^7M_{}`$ black hole, whereas an accretion disk is likely to have an effective radiating area $`10100`$ times larger. We return to this point in §5.

## 3 Tidal Disruption of a Star by a Supermassive Black Hole

The tidal disruption radius of a star in the vicinity of a Schwarzschild black hole is

$`r_T=\mu R_{}\left({\displaystyle \frac{M_H}{M_{}}}\right)^{1/3},`$ (5)
where $`M_H`$ is the mass of the black hole, $`M_{}`$ is the mass of the star, $`R_{}`$ is the radius of the star, and $`\mu `$ is a dimensionless coefficient. The exact value of $`\mu `$ is uncertain, though physically we expect it to be of order unity. We assume $`\mu =1`$ in our modeling.

The ratio of $`r_T`$ to $`r_H`$, where $`r_H=2GM_H/c^2`$ is the Schwarzschild radius of the black hole, is

$`{\displaystyle \frac{r_T}{r_H}}=5.08\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{2/3}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/3}\left({\displaystyle \frac{R_{}}{R_{}}}\right).`$ (6)
For a solar-type star with $`M_{}=M_{}`$ and $`R_{}=R_{}`$, we have $`r_T<r_H`$ when $`M_H>1.1\times 10^8M_{}`$. In such a case the star will be swallowed whole by the black hole, without a tidal disruption event (Hills, 1975).

After a star is tidally disrupted, there are two distinct stages in the dynamical evolution of the debris, each with its own time development and luminosity:

(1) The fallback stage (Rees, 1988): A fraction of the material in the disrupted star remains gravitationally bound to the black hole and falls back to the pericenter, giving rise to a mass accretion rate (and a luminosity) evolving with time as $`t^{5/3}`$ (Phinney, 1989). The returning material shocks and then circularizes to form an orbiting torus around the black hole at a radius $`r_c=2r_P`$ (by angular momentum conservation), where $`r_P`$ is the pericentric radius. The total energy radiated during the fallback is given by $`GM_H\mathrm{\Delta }M/2r_c`$, where $`\mathrm{\Delta }M`$ is the fallback mass. This energy can be less than the energy released in the subsequent viscous accretion stage discussed below. Nevertheless, since the fallback stage is relatively short-lived, it dominates the early luminosity of a disruption event.

(2) The viscous accretion stage (Cannizzo, Lee, & Goodman, 1990): The torus formed from returning debris gradually spreads inward and outward from $`r=r_c`$ under the action of viscosity, and gives rise to a mass accretion rate (and a luminosity) evolving with time approximately as $`t^{1.2}`$. The energy radiated during this stage is given by $`GM_\mathrm{H}\mathrm{\Delta }M[(1/2r_{ms})(1/2r_c)]`$, where $`r_{ms}`$ is the radius of the last stable circular orbit ($`=3`$ Schwarzschild radii for a non-spinning black hole). The accretion occurs on a viscous time scale, which is typically very long compared to the fallback time (Cannizo et al. 1990; Ulmer 1999; Appendix A of the present paper). Therefore, this stage is expected to last a long time, up to hundreds of years, depending on the magnitude of the viscosity in the disk, and the luminosity is significantly lower than that associated with the fallback stage (Appendix A).

In the case of NGC 5905, the X-ray luminosity was seen to rise rapidly during the first epoch of observations at $`t=1990.54`$ yr, indicating that the flare probably began around that time (see Table 1 and Fig. 1). The luminosity then dropped by more than a factor of $`100`$ between $`t=1990.54`$ yr and $`t=1996.89`$ yr. Such a rapid decline indicates that the observed outburst of NGC 5905 could not have been due to the accretion phase, but must have corresponded to the fallback stage (if it indeed was a tidal disruption event). Therefore, in this paper we assume that the flare in NGC 5905 corresponds to the fallback stage of tidal disruption.

Let us assume that the center of the star follows a nearly parabolic orbit (i.e., the binding energy of the star to the black hole is close to zero), and for definiteness let us assume that the pericenter of the orbit is at $`r_P=r_T`$, i.e. that the parameter $`\eta `$ defined below in equation (8) is equal to unity. Since the specific energy at the center of the star is $`E_c=v_P^2/2GM_H/r_T=0`$, the orbital velocity of the star at pericenter is $`v_P=(2GM_H/r_T)^{1/2}=c(r_H/r_T)^{1/2}`$, and the spread in the specific energy of the disrupted debris, $`2\mathrm{\Delta }E`$, is governed by the variation of the black hole gravitational potential across the star, and the spin-up of the star as a result of the tidal interaction (Rees, 1988). If $`M_HM_{}`$, we have

$`\mathrm{\Delta }E=k{\displaystyle \frac{GM_{}}{R_{}}}\left({\displaystyle \frac{M_H}{M_{}}}\right)^{1/3},`$ (7)
where $`k`$ depends on the spin-up state of the star. If the star is spun up to the break-up spin angular velocity, we have $`k3`$. On the other hand, if the spin-up effect is negligible, then we have $`k1`$ (Rees, 1988; Lacy, Townes, & Hollenbach, 1982; Evans & Kochanek, 1989; Ayal, Livio, & Piran, 2000). Since the spin-up velocity is always $`v_P`$, the energy used to spin up the star is much less than the orbital kinetic energy of the star, and the stellar orbit remains parabolic to a very good approximation.

The spin-up of a star via tidal interaction is a complex process. In linear perturbation theory, the spin-up angular velocity is given by (Press & Teukolsky, 1977; Alexander & Kumar, 2001; Alexander & Livio, 2001)

$`\omega _s{\displaystyle \frac{T_2(\eta )}{2\alpha \eta ^2}}\omega _P,\eta \left({\displaystyle \frac{r_P}{r_T}}\right)^{3/2},`$ (8)
where $`\eta 1`$ in our problem, $`\omega _Pv_P/r_P`$ is the orbital angular velocity of the star at the pericenter, $`\alpha `$ is the stellar momentum of inertia in units of $`M_{}R_{}^2`$, and $`T_2`$ is the second tidal coupling coefficient which is a function of $`\eta `$. In equation (8) we have assumed that $`M_HM_{}`$. For a main sequence star of mass $`0.76M_{}`$ and radius $`0.75R_{}`$, Alexander & Kumar (2001) give $`\alpha 0.07`$ and $`T_2(1)0.06`$, which implies that $`\omega _s0.43\omega _P`$ if $`\eta =1`$. For an $`n=1.5`$ polytrope star, they give $`\alpha 0.21`$ and $`T_2(1)0.36`$, which corresponds to $`\omega _s0.86\omega _P`$. Furthermore, they show that numerical simulations that include non-linear effects lead to a larger energy transfer from the orbit to the star and a larger spin-up than predicted by linear theory. For these reasons, we feel that it is likely a tidally disrupted star will be spun up to nearly the break-up limit, with $`\omega _s\omega _P`$ and $`k3`$. We therefore assume $`k=3`$ for most of our analysis; however, for completeness we also present a few results for the non-spinning case $`k=1`$.

Part of the disrupted star has positive energy and escapes from the system, and part remains bound. The bound material follows a highly eccentric orbit and returns to the central black hole after completing one orbit. The flare begins when the most bound material, with specific energy $`E_{\mathrm{min}}=\mathrm{\Delta }E`$, returns. The time since disruption for this to happen is given by

$`\mathrm{\Delta }t_1=2\pi GM_H(2\mathrm{\Delta }E)^{3/2}0.068\mathrm{yr}\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{1/2}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^1\left({\displaystyle \frac{R_{}}{R_{}}}\right)^{3/2},`$ (9)
where we have used $`r_P=r_T`$ and $`k=3`$ 111If the star is on a parabolic orbit with $`r_P<r_T`$, then the spin-up effect will be less important and we will have a smaller $`k`$ (but $`k`$ is still $`1`$), which will give rise to a $`\mathrm{\Delta }t_1`$ that is somewhat larger than that given by eq. (9).. Material that is less bound takes a progressively longer time to return, and the overall light curve is determined by the distribution of debris mass as a function of binding energy.

Rees (1988) assumed that the debris is uniformly distributed in mass between $`\mathrm{\Delta }E`$ and $`+\mathrm{\Delta }E`$. Numerical simulations (Evans & Kochanek, 1989; Ayal, Livio, & Piran, 2000) have shown that this is a reasonable approximation. The bound material then returns to pericenter at the rate (Phinney, 1989; Evans & Kochanek, 1989)

$`\dot{M}={\displaystyle \frac{2\mathrm{\Delta }M}{3\mathrm{\Delta }t_1}}\left({\displaystyle \frac{tt_D}{\mathrm{\Delta }t_1}}\right)^{5/3}=A\left({\displaystyle \frac{tt_D}{1\mathrm{y}\mathrm{r}}}\right)^{5/3}\mathrm{\Delta }MfM_{},`$ (10)
where $`t_D`$ is the time of the initial tidal disruption, $`\mathrm{\Delta }M`$ is the actual mass that falls back to pericenter, which is a fraction $`f`$ of the original mass of the star, and

$`A7.0\times 10^{23}\mathrm{g}\mathrm{s}^1\left({\displaystyle \frac{f}{0.1}}\right)\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{1/3}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/3}\left({\displaystyle \frac{R_{}}{R_{}}}\right).`$ (11)
In Rees’ original model, the whole star is disrupted and half the debris is on bound orbits, so that $`f=0.5`$. However, a recent numerical simulation shows that not all the returning material is captured by the black hole: about $`75\%`$ of the returning mass becomes unbound following the large compression it experiences on the way back (Ayal, Livio, & Piran, 2000). This gives rise to a smaller $`f0.1`$. Another possibility giving rise to a small $`f`$ is that the star is only partially disrupted: its envelope could be stripped by the black hole, leaving most of its core nearly intact (Renzini et al. 1995; see also §5 below). Here we treat $`f`$ as a free parameter.

The gravitational potential energy available from fallback is determined by the difference between the specific binding energy of the circularization orbit at $`r=2r_P=2r_T`$ and the specific binding energy of the incoming material. Since $`M_HM_{}`$, all the bound debris is on highly eccentric orbits with a specific binding energy much smaller than the binding energy of the final circular orbit. Thus, assuming that the fallback material radiates the energy release promptly, the radiation efficiency $`ϵ`$ is independent of time during fallback:

$`ϵ{\displaystyle \frac{GM_{}}{4R_{}c^2}}\left({\displaystyle \frac{M_H}{M_{}}}\right)^{2/3}0.025\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{2/3}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/3}\left({\displaystyle \frac{R_{}}{R_{}}}\right)^1.`$ (12)
The luminosity of the fallback process is then given by

$`L=ϵ\dot{M}c^21.55\times 10^{43}\mathrm{erg}\mathrm{s}^1\left({\displaystyle \frac{f}{0.1}}\right)\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{2/3}\left({\displaystyle \frac{tt_D}{1\mathrm{y}\mathrm{r}}}\right)^{5/3}.`$ (13)
The luminosity peaks at $`t=t_D+\mathrm{\Delta }t_1`$, i.e. when the most bound debris falls back to the pericenter, so we have

$`L_{\mathrm{peak}}1.36\times 10^{45}\mathrm{erg}\mathrm{s}^1\left({\displaystyle \frac{f}{0.1}}\right)\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{1/6}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{7/3}\left({\displaystyle \frac{R_{}}{R_{}}}\right)^{5/2}.`$ (14)
## 4 Comparison with the Flare in NGC 5905

In §2 we have fitted the decline of the NGC 5905 light curve to a $`(tt_D)^{5/3}`$ law, and obtained the result given in equation (1) with a disruption time of $`t_D=1990.36`$ yr. The fit is acceptable (see Fig. 1), and the time at which the light curve reached its peak is fairly well constrained. Since the initial measurements at $`t1990.54`$ yr showed flux rising with time, the maximum in the light curve was clearly later than this time. At the same time, the maximum could not have been later than 1990.56 yr, since the fitted light curve in (1) falls below $`L4.47\times 10^{42}`$ erg s-1 at that time. Thus, the light curve fit gives $`\mathrm{\Delta }t_1`$ within the range 0.18 to 0.20 yr, where $`\mathrm{\Delta }t_1`$ is the time lag between the peak of the light curve and the moment of tidal disruption. However, there is an error of $`\pm 0.03`$ yr on each limit. Therefore, we use conservative bounds on $`\mathrm{\Delta }t_1`$

$`0.15\mathrm{yr}<\mathrm{\Delta }t_1<0.23\mathrm{yr}.`$ (15)
¿From the discussion in §2, the mass involved in the outburst of NGC 5905 is only $`10^4M_{}`$. We therefore consider the following two possibilities for the candidate stellar object disrupted in NGC 5905:

(1) A low mass star (LMS), whose radius and mass are related by $`R_{}=R_{}(M_{}/M_{})`$ $`(0.08M_{}<M_{}<1M_{}`$).

(2) A substellar objects (SSOs, including brown dwarfs and planets), whose radius and mass are related by $`R_{}=0.06R_{}(M_{}/M_{})^{1/8}`$ ($`0.01M_{}<M_{}<0.08M_{}`$) (Chabrier & Baraffe, 2000)222Chabrier & Baraffe (2000) state that their fit is good down to $`M_{}=0.01M_{}`$, but a comparison with their plot shows that the fit is acceptable even for lower masses. We apply their scaling down to $`M_{}0.006M_{}`$..

Substituting these two radius-mass relations into equation (9) and noting the limits on $`\mathrm{\Delta }t_1`$ in equation (15), we find that the black hole mass in NGC 5905 and the mass of the disrupted star are restricted to the two relatively narrow shaded bands shown in Figure 2. The figure makes use of an additional constraint, namely a lower limit on the tidal radius $`r_T`$. If the black hole is not rotating, then the smallest distance to which a star on a parabolic orbit can approach and yet not be captured by the black hole is $`r=4r_g`$, where $`r_gGM_H/c^2=r_H/2`$. If we have a maximally rotating black hole and if the angular momentum of the stellar orbit is parallel to the spin of the black hole, then the closest approach possible is $`r=r_g`$. The real situation is probably in between these two limits, so we use a fiducial limiting radius of $`r=2r_g`$. Using equation (6), we have calculated relations connecting $`M_H`$ and $`M_{}`$ for all three limiting radii and plotted them as dotted lines in Figure 2.

The mass of the black hole in NGC 5905 is not very well constrained, but it is believed to be in the range $`M_H10^710^8M_{}`$ (Komossa, 2001), based on the correlation between bulge blue luminosity and black hole mass for spiral galaxies (Salucci et al., 2000) and on the correlation between bulge velocity dispersion and black hole mass for ellipticals and spirals (Gebhardt et al., 2000; Merritt & Ferrarese, 2001). This information is sufficient to constrain the mass of the disrupted star fairly tightly: the mass must be either in the range $`0.010.02M_{}`$ in the case of a SSO or $`0.61M_{}`$ in the case of a LMS.

The result, however, depends sensitively on the choice we make for the tidal spin-up parameter $`k`$ discussed in §2. Figure 2 has been calculated for $`k=3`$, a reasonable and possibly likely value, but in principle $`k`$ could be as small as 1. Figure 3 shows how the mass constraints change when we use the latter value. We see that the allowed range of models is limited to somewhat lower black hole masses, unless the black hole spins very rapidly and allows tidal disruptions down to $`r_Tr_g`$. We do not discuss the $`k=1`$ case further.

¿From equations (13) and (1) we infer that

$`{\displaystyle \frac{M_H}{10^7M_{}}}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{2/3}={\displaystyle \frac{0.0019}{f}},`$ (16)
where $`f`$ is the fraction of the mass of the star that returns as fallback debris. The dashed lines in Figure 2 correspond to $`f=0.0005`$, $`0.012`$, and $`0.12`$ respectively. We see that, for the case of a LMS, $`f`$ must be $`<0.0005`$, which means that an unusually small fraction of the star must have participated in the fallback. For the case of a SSO, somewhat larger values of $`f`$ are obtained (but still rather small).

In Figures 4 and 5 we plot the luminosities given by equation (13) for different values of $`M_H`$, $`M_{}`$, and $`f`$, where the disrupted object is assumed to be a LMS (Fig. 4) and a SSO (Fig. 5), respectively, and the disruption time is taken to be $`t_D=1990.36`$ yr. We have chosen values of $`(M_H,M_{})`$ from the shaded regions in Figure 2 so that equation (15) and the condition $`r_T>2r_g`$ are satisfied, and chosen $`f`$ such that $`L_{\mathrm{peak}}4.47\times 10^{42}`$ erg s-1 where $`L_{\mathrm{peak}}L(t=t_D+\mathrm{\Delta }t_1)`$ is the peak luminosity.

We confirm that, in order for the models to fit the observational data, a very small $`f`$ ($`<0.0003`$) is required for the case of a LMS (Fig. 4). For the case of a SSO (Fig. 5) a fairly small $`f`$ ($`0.04`$) is again required if $`M_H2\times 10^7M_{}`$. A largish $`f`$ ($`>0.1`$) can be obtained only if the black hole has a very small mass: $`M_H<4\times 10^6M_{}`$.

In Figures 4 and 5 we also indicate the ratio of the peak flux, $`F_{\mathrm{peak}}`$, in the fitted light curve to the Eddington flux, $`F_{\mathrm{Edd}}`$, where $`F_{\mathrm{peak}}L_{\mathrm{peak}}/(\pi R_X^2)`$ with $`R_X`$ given by equation (4), and $`F_{\mathrm{Edd}}L_{\mathrm{Edd}}/(4\pi r_T^2)`$, with $`L_{\mathrm{Edd}}=1.3\times 10^{45}\mathrm{erg}\mathrm{s}^1(M_H/10^7M_{})`$ and $`r_T`$ given by equation (6). We see that the ratio $`F_{\mathrm{peak}}/F_{\mathrm{Edd}}`$ is significantly larger than unity if the color temperature factor $`f_c1`$, and it is not much less than unity even if $`f_c3`$ (the likely value). We discuss the implications of this result in the next section.

## 5 Summary and Discussion

The X-ray light curve of the flare in NGC 5905 agrees with the $`(tt_D)^{5/3}`$ dependence predicted for fallback of debris after tidal disruption of a star by a supermassive black hole (Fig. 1). If the black-body fit presented in Bade et al. (1996) is a good representation of the spectrum of the emission, then the bolometric correction to the X-ray flux is not very large and the observed light curve gives a meaningful measure of the energetics and luminosity of the flare. We note that, during and after the X-ray outburst, the optical brightness did not show any significant variation (Komossa & Bade, 1999). This confirms that the bulk of the fallback emission probably did occur in X-rays. (If the flare had substantial hidden emission in the far ultraviolet band, it is hard to see how there could have been no variations in the optical.)

In addition to fallback, a second independent stage of evolution is expected when the circularized fallback material viscously accretes onto the black hole. This stage is expected to develop much more slowly and to be very much dimmer than the fallback stage (Appendix A). It is not relevant for the particular flare in NGC 5905 considered here.

The time lag between the estimated disruption time $`t_D=1990.36`$ yr of the event and the time $`t_1`$ corresponding to the peak of the light curve gives the orbital period $`\mathrm{\Delta }t_1`$ of the most tightly bound debris. The data constrain $`\mathrm{\Delta }t_1`$ fairly tightly (eq. ). By the arguments described in §3, this translates into constraints on the mass $`M_H`$ of the black hole and the mass $`M_{}`$ of the disrupted star. The constraints are shown in Figure 2.

In deriving equation (9), we have estimated the binding energy of the most tightly bound debris assuming that the star is fully spun up by tidal interaction before being disrupted. This corresponds to a spin-up parameter $`k=3`$ (§3); this value is supported both by linear analysis (Press & Teukolsky, 1977) and numerical simulations (Alexander & Kumar, 2001). If on the other hand there is no tidal spin-up, then we expect $`k=1`$, and the results change significantly (Fig. 3). We consider $`k=1`$ to be somewhat unlikely, and have focused on the $`k=3`$ results.

There is also ambiguity in the value of the parameter $`\mu `$ in the definition of $`r_T`$ \[eq. (5)\]. Numerical simulations show that $`\mu `$ is usually of order unity, the value we use, but that it can be as large as $`1.7`$ in specific cases (Frolov et al., 1994; Diener et al., 1997). The time interval $`\mathrm{\Delta }t_1`$ \[eq. (9)\] is very sensitive to the value of $`\mu `$: $`\mathrm{\Delta }t_1\mu ^3`$. Because of the tight observational constraint on $`\mathrm{\Delta }t_1`$ (Figs. 2 and 3), $`\mu `$ cannot be much larger than $`1`$ since otherwise it would imply too small a black hole mass.

Another constraint on $`M_H`$ and $`M_{}`$ is obtained from the requirement that the radius $`r_T`$ at which the star is tidally disrupted must lie outside the radius of the marginally bound circular orbit. This requires $`r_T>4r_g`$ for a Schwarzschild black hole and $`r_T>r_g`$ for the most favorably oriented orbit around a maximally spinning Kerr hole. The uncertainty in the limiting $`r_T`$ is indicated by the dotted lines in Figure 2. For concreteness, we have assumed that the limiting radius is intermediate between the two extreme cases mentioned above; we have used the limit $`r_T>2r_g`$. Yet another source of uncertainty derives from the fact that we have used Newtonian dynamics in our analysis, whereas it is clear that relativistic effects should be important (Luminet & Marck, 1985; Laguna et al., 1993; Frolov et al., 1994; Diener et al., 1997).

Ignoring these caveats, we obtain a tight correlation between the mass of the disrupted star in NGC 5905 and the mass of the central black hole. As Figure 2 shows, the data are consistent with two distinct solution branches.

One solution branch corresponds to a relatively massive black hole, $`M_H`$ close to $`10^8M_{}`$, tidally disrupting a main sequence star with $`M_{}M_{}`$. The other branch consists of a lower mass black hole, $`M_H<10^{7.5}M_{}`$, disrupting a brown dwarf or a giant planet with $`M_{}0.02M_{}`$. Assuming that NGC 5905 satisfies the well-known correlations between black hole mass and host galaxy properties (Salucci et al., 2000; Gebhardt et al., 2000; Merritt & Ferrarese, 2001), the black hole mass is estimated to be $`M_H10^710^8M_{}`$ (Komossa, 2001). This mass range permits both solution branches.

Figures 4 and 5 show sample theoretical light curves corresponding to the two solution branches. In each case, there is one free parameter, namely the fraction $`f`$ of the mass of the disrupted star that returns as fallback debris to produce the observed radiation. We find that $`f`$ is unusually small if the disrupted star is a main sequence star; for the particular example shown in Figure 4, we require $`f<0.0003`$, which is an extremely small value compared to the theoretically expected value of $`f0.10.5`$. Brown dwarf models give somewhat larger estimates for $`f`$ (Fig. 5), though still uncomfortably small, unless the black hole has an unusually low mass ($`<10^{6.5}M_{}`$). Though the estimated $`f`$ can be increased by adopting a larger ambiguity factor $`\mu `$ in the definition of $`r_T`$, we do not consider this likely since a larger $`\mu `$ will lead to an unacceptably small value of $`M_H.`$

The above results on $`f`$ are not surprising. As we showed in §2, the total radiative energy in the flare is only $`4\times 10^{49}`$ erg, which for any reasonable radiative efficiency corresponds to an unexpectedly small mass $`\mathrm{few}\times 10^4M_{}`$. Figures 4 and 5 merely confirm this result using a detailed estimate of the radiative efficiency (eq. 12).

Can $`f`$ be so small for a tidally disrupted solar mass star? Rees (1988) and Renzini et al. (1995) have argued that a star could lose just its outer layers if the minimum distance from the star to the black hole is somewhat larger than $`r_T`$. Here we point out a related possibility. The usual formula for the tidal disruption radius, $`r_T=R_{}(M_H/M_{})^{1/3}`$, is correct only for a homogeneous star with a uniform mass density. For a real star, we could define a variable tidal disruption radius,

$`r_T(R)=R\left[{\displaystyle \frac{M_H}{M(R)}}\right]^{1/3}=\left[{\displaystyle \frac{3M_H}{4\pi \overline{\rho }(R)}}\right]^{1/3},`$ (17)
where $`RR_{}`$ is a radius inside the star, $`M(R)`$ is the mass contained within $`R`$, and $`\overline{\rho }(R)`$ is the corresponding mean density. The meaning of $`r_T(R)`$ is as follows: if the star is located at a distance $`r=r_T(R)`$ from the black hole, then the part of the star outside radius $`R`$ (i.e., the spherical envelope between $`R`$ and $`R_{}`$) is tidally stripped off, and the core inside $`R`$ remains intact. Obviously, $`r_{T,\mathrm{max}}=r_T(R_{})`$ and $`r_{T,\mathrm{min}}=(3M_H/4\pi \rho _c)^{1/3}`$, where $`\rho _c`$ is the central density of the star. If the star is at a distance $`r`$ from the black hole such that $`r_{T,\mathrm{min}}<r<r_{T,\mathrm{max}}`$, then only a part of the envelope of the star is peeled off.

One important question is how the core of the star responds when it loses its envelope. What is relevant is the response on a dynamical time. Webbink (1985) has calculated the logarithmic derivative of the radius of a main sequence star with respect to mass, $`\zeta _s(\mathrm{ln}R_{}/\mathrm{ln}M_{})_s`$, where the subscript $`s`$ indicates that the entropy profile of the star is held fixed (as appropriate for a dynamical process). Consider a star that loses a fraction of its mass as a result of a tidal interaction and assume that it dynamically adjusts to a new equilibrium corresponding to its reduced mass. If $`\zeta _s>1/3`$, then the mean density of the new configuration will be higher than the original mean density of the star. The leftover mass is then stable to further stripping. However, if $`\zeta _s<1/3`$, the new configuration will be prone to further tidal mass loss, and one has a runaway situation in which the entire star is likely to be disrupted.

Webbink (1985) shows that $`\zeta _s>1/3`$ for $`\mathrm{log}(M_{}/M_{})0.2`$. Since this is precisely the mass range of interest if the disrupted star in NGC 5905 was a main sequence star (Fig. 2), we conclude that partial stripping is a viable possibility. To estimate what fraction of its mass the star might lose, consider a star like the Sun as an example. For the Sun, the mass inside a radius $`R_1=0.7R_{}`$ is $`0.975M_{}`$, and correspondingly, $`r_T(R_1)=0.706r_{T,\mathrm{max}}`$. Thus, for a range of distances between $`r_{T,\mathrm{max}}`$ and $`0.7r_{T,\mathrm{max}}`$ (a non-negligible range), such a star will lose only a couple of percent of its mass. If we assume that 10% of this mass returns via fallback, then we can have $`f`$ as small as $`0.002`$. This is still much larger than the value needed to explain the flare in NGC 5905, but it is at least within striking distance.

Fully convective lower main sequence stars and degenerate brown dwarfs have $`\zeta _s1/3`$. These stars are likely to be fully disrupted in tidal encounters since they expand to a larger radius as they lose mass.

Apart from the small value of $`f`$, another unexpected result of our analysis is that the peak radiative flux emitted during the flare, $`F_{\mathrm{peak}}=\sigma T_{\mathrm{eff}}^4=\sigma (T_{\mathrm{bb}}/f_c)^4`$, is close to the Eddington flux. (Here $`T_{\mathrm{bb}}`$ is the black-body temperature of the radiation and $`f_c`$ is a color correction factor.) The results are shown for our models in Figures 4 and 5. Even with $`f_c3`$ (§2), we find $`F_{\mathrm{peak}}/F_{\mathrm{Edd}}0.10.3`$; the value increases rapidly with decreasing $`f_c`$ (as $`f_c^4`$). What makes this surprising is that the luminosity of the source, even at the peak of the flare, is much below the Eddington luminosity: compare $`L_{\mathrm{peak}}4.5\times 10^{42}\mathrm{erg}\mathrm{s}^1`$ with $`L_{\mathrm{Edd}}10^{45}(M_H/10^7M_{})\mathrm{erg}\mathrm{s}^1`$. Clearly, the radiation in the flare must have been emitted from a small patch on the source, and the patch must have projected a small solid angle at the central black hole. This result was anticipated in the model-independent analysis of §2, where we showed that the effective size of the emitting patch is small compared to the area we might expect if the entire accretion disk out to a few $`r_g`$ is involved in the emission.

The small size of the emitting region can be understood in analogy with mass transfer binary stars. In the latter systems, it is known that the cold gas stream from the secondary star impacts the accretion disk at a compact “bright spot” on the outer edge of the disk. Most of the radial kinetic energy of the stream is converted to thermal energy near this spot and radiated locally. The analogy with the tidal disruption problem is close because the returning material after a disruption event comes in similarly as a narrow cold stream (Kochanek, 1994; Kim, Park, & Lee, 1999). It is therefore natural that the radiation during the fallback stage should be emitted from a relatively small patch on the source. The effect is probably even more severe in the tidal disruption problem. Lense-Thirring precession may well cause the incoming stream to be misaligned with respect to the orbital plane of the accretion disk (Kochanek, 1994), causing the geometric area of the impact spot to become smaller. (We note that in binaries the stream appears to skim over the surface of the disk after impacting at the bright spot, thereby increasing the effective area over which the kinetic energy of the stream is dissipated. This effect would be reduced if the orbit of the stream is not in the plane of the disk.)

Having decided that emission from a small patch is natural for this problem, we should next ask whether it is a pure accident that the flux from the patch in the NGC 5905 flare happens to be nearly equal to the Eddington flux. We would like to suggest that it is perhaps not an accident.

Imagine a narrow stream with a mass accretion rate $`\dot{M}`$ impacting on the accretion torus at a radius $`r_T`$ and over a small patch that subtends a solid angle $`\mathrm{\Delta }\mathrm{\Omega }_{\mathrm{stream}}`$ at the black hole. Let us suppose that the kinetic energy in the stream is thermalized and radiated immediately from the patch. The radiated flux is locally given by $`ϵ\dot{M}c^2/\mathrm{\Delta }\mathrm{\Omega }_{\mathrm{stream}}r_T^2`$. What would happen if this flux is greater than the local $`F_{\mathrm{Edd}}`$? Jeremy Goodman (private communication) suggests that the super-Eddington flux would cause the local radiating atmosphere to puff up and that the radiation would be emitted from a larger patch than the original impact region. Moreover, it would be natural for the radiating patch to self-adjust so that the flux is close to the Eddington limit. According to this picture then, the radius $`R_X`$ estimated in equation (4) refers to the size of the puffed up atmosphere, while the actual transverse size of the stream is smaller. If the stream is extremely narrow, then the flux would be close to Eddington during a significant fraction of the light curve. Equivalently, the black-body temperature of the radiation would remain constant, even while the luminosity decreases.

This picture (which is predicated on the assumption of a narrow fallback stream) suggests that most tidal disruption events should be visible in soft X-rays; the condition $`F=F_{\mathrm{Edd}}`$, coupled with $`f_c3`$ (see §2), ensures that the blackbody temperature of the radiation will be $`50`$ eV for a $`10^7M_{}`$ black hole, which is well-matched to X-ray observations. Soft X-ray surveys should thus be well-suited for an investigation of both the frequency of tidal disruption events and the amount of mass accreted in each event. In this connection, it is interesting that, with the exception of the UV flare of Renzini et al. (1995), all known candidate flares for tidal disruption events have been identified in the X-ray band (see Komossa 2001 for five candidates and Choi et al. 2002 for a possible sixth candidate).

Menou & Quataert (2001) have emphasized that, given current estimates of the rate of tidal disruptions, the nuclei of many nearby normal galaxies should be much brighter than observed if a fair fraction of a stellar mass is accreted by the black hole after each tidal disruption event. This is not relevant for the late accretion phase of the particular NGC 5905 flare discussed here because, as we have discussed (see also the Appendix below), very little mass ($`1M_{}`$) returned to the BH in this case. However, whether we attribute the small amount of mass in the NGC 5905 flare to partial stripping of a star or to the disruption of a substellar object, it does not help to explain the general problem of why nearby galactic nuclei are so faint. In addition to partial disruptions of the sort that might have occurred in NGC 5905, we expect total stellar disruptions in significant numbers and these should provide a significant mean mass accretion rate in galactic nuclei. Disruptions of substellar objects only add to the purely stellar disruptions considered to date. Thus, alternative ways of explaining the dimness of nearby galactic nuclei are required (Menou & Quataert 2001; Narayan 2002), beyond the specific scenario outlined here for the giant flare of NGC 5905. Either the mass must accrete in a radiatively inefficient mode or the accretion must be suppressed or slowed down significantly.

The authors gratefully thank Stefanie Komossa for kindly providing the ROSAT data listed in Table 1 along with the conversion factors from count rate to X-ray luminosity. The authors also thank Jeremy Goodman, Bohdan Paczyński and the referee for helpful comments. LXL and RN thank the Institute for Advanced Study and the Department of Astrophysical Sciences, Princeton, for hospitality while this work was being done. LXL’s and KM’s research was supported by NASA through Chandra Postdoctoral Fellowship grant numbers PF1-20018 (LXL) and PF9-10006 (KM) awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073. RN’s research was supported in part by NSF grant AST-9820686.

## Appendix A Viscous Evolution of the Thin Disk Formed from Tidal Disruption

After circularization, the fallback material from a tidal disruption event will form a torus at a radius $`r2r_T`$. The torus will then spread out viscously and accrete onto the black hole.

The time-scale for this stage of the evolution is determined by the viscous time-scale (Cannizzo, Lee, & Goodman, 1990; Ulmer, 1999)

$`t_{\mathrm{vis}}={\displaystyle \frac{t_{\mathrm{Kep}}(2r_T)}{\pi \alpha h^2}}2.9\times 10^3\mathrm{yr}\left({\displaystyle \frac{\alpha }{0.1}}\right)^1\left({\displaystyle \frac{h}{10^3}}\right)^2\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/2}\left({\displaystyle \frac{R_{}}{R_{}}}\right)^{3/2},`$ (A1)
where $`t_{\mathrm{Kep}}`$ is the Keplerian orbital period, $`\alpha `$ is the standard viscosity parameter (Shakura & Sunyaev, 1973), and $`h`$ is the ratio of the disk height to radius. It is easy to check that for the case of the flare in NGC 5905, where $`L_{\mathrm{peak}}10^3L_{\mathrm{Edd}}(M_H/10^7M_{})^1`$, the accretion disk would correspond to the regime of the “middle disk” in which electron scattering opacity dominates over the Kramers opacity and gas pressure dominates over radiation pressure (Shakura & Sunyaev, 1973; Novikov & Thorne, 1973; Frank, King & Raine, 1992). Then, $`h`$ is estimated to be

$`h1\times 10^3\left({\displaystyle \frac{\alpha }{0.1}}\right)^{1/10}\left({\displaystyle \frac{\dot{M}_2}{10^3\dot{M}_{\mathrm{Edd}}}}\right)^{1/5}\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{1/10}\left({\displaystyle \frac{r}{r_H}}\right)^{1/20},`$ (A2)
where $`\dot{M}_2`$ is the mass accretion rate, $`\dot{M}_{\mathrm{Edd}}L_{\mathrm{Edd}}/(0.1c^2)`$, and $`r`$ is the radius in the disk. In equation (A2) we have neglected a dimensionless factor $`1`$.

Inserting equation (A2) into equation (A1) and setting $`r=2r_T`$, we obtain

$`t_{\mathrm{vis}}`$ $``$ $`2.3\times 10^3\mathrm{yr}\left({\displaystyle \frac{\alpha }{0.1}}\right)^{4/5}\left({\displaystyle \frac{\dot{M}_2}{10^3\dot{M}_{\mathrm{Edd}}}}\right)^{2/5}`$ (A3)
$`\times \left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{4/15}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{7/15}\left({\displaystyle \frac{R_{}}{R_{}}}\right)^{7/5}.`$
The large value of $`t_{\mathrm{vis}}`$ confirms that the flare in NGC 5905 corresponded to the fallback stage, not the viscous accretion stage.

From the discussion in §3 on the fallback and viscous accretion stages associated with a tidal disruption event, the ratio of the total energy radiated during the accretion stage to the total energy radiated during the fallback stage is $`2r_P/r_{ms}1`$. Therefore, we can estimate the ratio of the peak of the viscous luminosity to the peak of the fallback luminosity by

$`{\displaystyle \frac{L_{\mathrm{vis}}}{L_{\mathrm{fb}}}}`$ $``$ $`\left({\displaystyle \frac{2r_P}{r_{ms}}}1\right){\displaystyle \frac{\mathrm{\Delta }t_1}{t_{\mathrm{vis}}}}`$ (A4)
$``$ $`2.4\times 10^5\left({\displaystyle \frac{\alpha }{0.1}}\right)\left({\displaystyle \frac{h}{10^3}}\right)^2\left({\displaystyle \frac{M_H}{10^7M_{}}}\right)^{1/2}\left({\displaystyle \frac{M_{}}{M_{}}}\right)^{1/2}\left({\displaystyle \frac{2r_P}{r_{ms}}}1\right),`$
where equations (9) and (A1) have been used.

The extremely large value of $`t_{\mathrm{vis}}`$ and the extremely small value of $`L_{\mathrm{vis}}/L_{\mathrm{fb}}`$ imply that the accretion stage will never be observed in practice for the particular case of NGC 5905. \[In fact, according the above estimates, one may expect a new tidal disruption before the disk evolves significantly (Syer & Ulmer, 1999; Magorrian & Tremaine, 1999).\]


format: text

text:
# Statistical Fluctuations of Electromagnetic Transition Intensities and Electromagnetic Moments in 𝑝𝑓-Shell Nuclei

## Abstract

We study the fluctuation properties of $`\mathrm{\Delta }T=0`$ electromagnetic transition intensities and electromagnetic moments in $`A60`$ nuclei within the framework of the interacting shell model, using a realistic effective interaction for $`pf`$-shell nuclei with a 56Ni core. The distributions of the transition intensities and of the electromagnetic moments are well described by the Gaussian orthogonal ensemble of random matrices. In particular, the transition intensity distributions follow a Porter-Thomas distribution. When diagonal matrix elements (i.e., moments) are included in the analysis of transition intensities, we find that the distributions remain Porter-Thomas except for the isoscalar $`M1`$. The latter deviation is explained in terms of the structure of the isoscalar $`M1`$ operator.

Random matrix theory (RMT)  was originally introduced to explain the statistical fluctuations of neutron resonances in compound nuclei . The theory assumes that the nuclear Hamiltonian belongs to an ensemble of random matrices that are consistent with the fundamental symmetries of the system. In particular, since the nuclear interaction preserves time-reversal symmetry, the relevant ensemble is the Gaussian orthogonal ensemble (GOE). The use of RMT in the compound nucleus was justified by the complexity of the nuclear system. Bohigas et al  conjectured that RMT describes the statistical fluctuations of a quantum system whose associated classical dynamics is chaotic. RMT has become a standard tool for analyzing the universal statistical fluctuations in chaotic systems .

The chaotic nature of the single-particle dynamics in the nucleus can be studied in the mean-field approximation. The interplay between shell structure and fluctuations in the single-particle spectrum has been understood in terms of the classical dynamics of the nucleons in the corresponding deformed mean-field potential . However, the residual nuclear interaction mixes different mean-field configurations and affects the statistical fluctuations of the many-particle spectrum and wave functions. These fluctuations can be studied using various nuclear structure models. The statistics of the low-lying collective part of the nuclear spectrum have been studied in the framework of the interacting boson model , in which the nuclear fermionic space is mapped onto a much smaller space of bosonic degrees of freedom. Because of the relatively small number of degrees of freedom in this model, it was also possible to relate the statistics to the underlying mean-field collective dynamics. At higher excitations, additional degrees of freedom (such as broken pairs) become important , and the effects of interactions on the statistics must be studied in larger model spaces. The interacting shell model offers an attractive framework for such studies, in which realistic effective interactions are available and where the basis states are labeled by exact quantum numbers of angular momentum, isospin and parity .

RMT makes definite predictions for the statistical fluctuations of both the eigenfunctions and the spectrum. The electromagnetic transition intensities in a nucleus are observables that are sensitive to the wave functions, and the study of their statistical distributions should complement  the more common spectral analysis. Using the interacting shell model, $`B(M1)`$ and $`B(E2)`$ transitions were recently analyzed in 22Na  and found to follow the Porter-Thomas distribution , in agreement with RMT and with the previous finding of a Gaussian distribution for the eigenvector components -.

Most studies of statistical fluctuations in the shell model have been restricted to lighter nuclei ($`A40`$) where complete $`0\mathrm{}\omega `$ calculations are feasible (e.g., $`sd`$-shell nuclei). Here we study the fluctuation properties of $`\mathrm{\Delta }T=0`$ electromagnetic transition intensities in nuclei with $`A60`$. We also study the statistics of the diagonal matrix elements, which describe the corresponding electric or magnetic moments. We find that the statistics of both the transition intensities and the electromagnetic moments are well described by RMT. In RMT, off-diagonal matrix elements have a Gaussian distribution with zero mean, leading to a Porter-Thomas distribution of the transition intensities. Diagonal elements also have a Gaussian distribution but with non-zero mean. Thus the squares of the diagonal matrix elements should also follow a Porter-Thomas distribution but only after the mean value is subtracted from the diagonal elements. Experimentally it is difficult to measure the moments of excited states, while transition intensities are more readily accessible.

When the good quantum numbers (isospin, spin and parity) of the initial and final states are the same, we can include the diagonal matrix elements in the analysis of the transition intensity distribution. The number of diagonal elements is much smaller than the number of off-diagonal elements, and generically we do not expect the diagonal elements to modify the Porter-Thomas statistics of the transition intensities. This is confirmed in our shell model calculations with the exception of the $`M1`$ operator in self-conjugate nuclei ($`T_z=0`$), where a significant deviation from a Porter-Thomas distribution is found.***We find a similar deviation in $`sd`$-shell nuclei. This deviation is explained in terms of the special structure of the isoscalar $`M1`$ operator, which leads to diagonal matrix elements that are much larger than its off-diagonal elements.

Full $`pf`$-shell calculations are currently not feasible for $`A60`$ nuclei. We therefore perform $`pf`$-shell calculations with 56Ni as a core, i.e., we assume fully occupied $`f_{7/2}`$ orbits and consider all possible many-nucleon configurations defined by the $`0f_{5/2}`$, $`1p_{3/2}`$ and $`1p_{1/2}`$ orbitals. The effective interaction is chosen to be the isospin-conserving F5P interaction . This interaction is successful in describing the mass range $`A5768`$. The calculations were performed using the shell model program OXBASH .

To test the validity of RMT in the above model space, we first study the spectral fluctuations for states with good spin and isospin (all states in the $`pf`$-shell model space have the same parity). Fig. 4(a) shows the nearest-neighbors level spacing distribution $`P(s)`$ for the unfolded $`T=0`$ levels. To improve the statistics we combine the spacings of the $`J=2`$ states and the $`J=4`$ states. The calculated distribution (histograms) is in agreement with the Wigner-Dyson distribution (solid line), and level repulsion is clearly observed at small spacings. The Poisson distribution, which corresponds to a random sequence of levels (and characterizes regular systems), is shown by the dashed line for comparison. Another measure of fluctuations is the spectral rigidity described by Dyson’s $`\mathrm{\Delta }_3`$ statistic . The bottom panels of Fig. 4 show $`\mathrm{\Delta }_3(L)`$ versus $`L`$ for $`J=2`$ states with $`T=0`$ \[panel (b)\] and $`T=1`$ \[panel (c)\]. The results agree well with the GOE limit (solid lines).

We now consider the electromagnetic transition intensities. Denoting by $`B(\overline{\omega }L;if)`$ the reduced transition probability from an initial state $`|i`$ to a final state $`|f`$, with $`\overline{\omega }`$ indicating the electric ($`E`$) or magnetic ($`M`$) character of the transition, and $`2^L`$ the multipolarity, we have

$$B(\overline{\omega }L;J_iT_iT_zJ_fT_fT_z)=\frac{|\delta _{T_iT_f}M_{is}(\overline{\omega }L)(T_iT_z10|T_fT_z)M_{iv}(\overline{\omega }L)|^2}{(2J_i+1)(2T_i+1)}.$$
(1)
Here $`M_{is}(\overline{\omega }L)`$ and $`M_{iv}(\overline{\omega }L)`$ are the triply reduced matrix elements for the isoscalar and isovector components of the transition operator, respectively. Note that these matrix elements depend on $`J_i,T_i`$ and $`J_f,T_f`$ but not on $`T_z`$. For $`\mathrm{\Delta }T=0`$ transitions ($`T_i=T_f=T`$) the isospin Clebsch-Gordan coefficient in Eq. (1) is simply given by

$$(TT_z10|TT_z)=T_z/\sqrt{T(T+1)}.$$
(2)
It follows that the isovector component in Eq. (1) is absent for self-conjugate (i.e. $`T_z=0`$) nuclei. Thus the statistics of the isoscalar component of an electromagnetic transition operator can be inferred directly from the study of $`\mathrm{\Delta }T=0`$ transitions in $`T_z=0`$ nuclei. Consequently, we are able to test the sensitivity of the statistics to the isovector and the isoscalar contributions. Below we present results for E2 and M1 transitions. The $`E2`$ transitions were calculated using standard effective charges of $`e_p=1.5`$ and $`e_n=0.5`$. We note, however, that the isoscalar and isovector components of $`E2`$ are (up to a proportionality constant) independent of the effective charges, and their corresponding statistics are thus also independent of the particular choice of effective charges.

To study the fluctuation properties of the transition rates, it is necessary to divide out any secular variation of the average strength function versus the initial and final energies. We do this by applying the method of Refs.  and . We calculate an average transition strength at an initial energy $`E`$ and final energy $`E^{}`$ from

$`B(\overline{\omega }L;E,E^{})={\displaystyle \frac{_{i,f}B(\overline{\omega }L;if)e^{(EE_i)^2/2\gamma ^2}e^{(E^{}E_f)^2/2\gamma ^2}}{_{i,f}e^{(EE_i)^2/2\gamma ^2}e^{(E^{}E_f)^2/2\gamma ^2}}},`$ (3)
where $`\gamma `$ is a parameter chosen as described below. For fixed values of the initial ($`J_i,T`$) and final ($`J_f,T`$), we calculate from Eq. (1) the intensities $`B(\overline{\omega }L;if)`$. All transitions of a given operator (e.g., $`M1`$ or $`E2`$) between the initial and final states of the given spin and isospin classes have been included in the statistics. We remark that the energy levels used in (3) are the unfolded energy levels characterized by a constant mean spacing. The value of $`\gamma `$ in (3) has been chosen to be large enough to minimize effects arising from the local fluctuations in the transition strength but not so large as to wash away the secular energy variation of the average intensity. In the present calculations we use $`\gamma 2.5`$. We renormalize the actual intensities by dividing out their smooth part

$$y_{fi}=\frac{B(\overline{\omega }L;J_iTT_zJ_fTT_z)}{B(\overline{\omega }L;E,E^{})},$$
(4)
and construct their distribution using bins that are equally spaced in $`\mathrm{log}_{10}y`$. The choice of $`\mathrm{log}_{10}y`$ as the variable allows us to display the distribution of the weak transitions over several orders of magnitude .

In RMT we expect a Porter-Thomas distribution for $`P(y)`$, i.e., a $`\chi ^2`$ distribution distribution in $`\nu =1`$ degrees of freedom . A $`\chi ^2`$ distribution in $`\nu `$ degrees of freedom is given by

$$P_\nu (y)=(\nu /2<y>)^{\nu /2}y^{\nu /21}e^{\nu y/2<y>}/\mathrm{\Gamma }(\nu /2).$$
(5)
Indeed, consider a transition operator $`𝒯`$ and its matrix elements $`f|𝒯|i`$ between an initial state $`|i`$ and a final state $`|f`$. Using a complete (and fixed) set of states $`|m`$, we have

$$f|𝒯|i=\underset{m,n}{}𝒯_{m,n}\psi _f^{}(m)\psi _i(n),$$
(6)
where $`\psi _i(n)n|i`$ represents the eigenfunction $`|i`$ in the fixed chosen basis and $`T_{m,n}=m|T|n`$. In RMT, $`\psi _i(n)`$ are Gaussian random variables, and by the central limit theorem it follows from (6) that $`f|𝒯|i`$ is also Gaussian. In RMT, two different eigenfunctions are uncorrelated $`\overline{\psi _f^{}(m)\psi _i(n)}=0`$, and consequently $`\overline{f|𝒯|i}=0`$. Similarly, using Wick’s theorem for Gaussian variables, the variance is given by

$$\sigma ^2(f|𝒯|i)=\underset{m,n,m^{},n^{}}{}𝒯_{m,n}^{}𝒯_{m^{},n^{}}\overline{\psi _f(m)\psi _f^{}(m^{})}\overline{\psi _i(n)\psi _i^{}(n^{})}=\frac{1}{N_iN_f}\mathrm{Tr}\left(P_i𝒯^{}P_f𝒯\right),$$
(7)
where $`P_i`$ is the projection operator on the $`N_i`$-dimensional subspace of eigenstates with quantum numbers $`J_i^{\pi _i},T_i`$ ($`P_f`$ is similarly the projection operator on the corresponding subspace of final states). The transition intensity is proportional to the square of the matrix element, and thus has a Porter-Thomas distribution (i.e., Eq. (5) with $`\nu =1`$).

We examined $`E2`$ and $`M1`$ $`2^+,T=12^+,T=1`$ transitions in $`A=60`$ nuclei. The $`T=1`$ states form isobaric multiplets with $`T_z=0,\pm 1`$ (i.e., 60Co, 60Zn and 60Cu). We studied the statistics of the transition intensities in both $`T_z=0`$ and $`T_z=1`$ nuclei. Because of the vanishing of the isovector Clebsch-Gordan coefficient in Eq. (1), the transitions in the self-conjugate nucleus 60Zn ($`T_z=0`$) are purely isoscalar. For each transition operator we sampled $`66^266=4290`$ transition matrix elements. The calculated distributions (histograms) of the $`B(E2)`$ (left panels) and $`B(M1)`$ (right panels) $`2^+,T=12^+,T=1`$ transitions are shown in Fig. 4. The top and bottom panels correspond to $`T_z=1`$ and $`T_z=0`$ nuclei, respectively. All distributions are well described by the Porter-Thomas distribution. Similar results were found by Adams et al  for the $`N=Z`$ ($`T_z=0`$) nucleus 22Na in the $`sd`$-shell.

When the quantum numbers of the initial and final states are identical, we can also examine the statistics of the diagonal matrix elements $`i|𝒯|i`$. Using (6) with $`i=f`$ and the central limit theorem we still find that $`i|𝒯|i`$ is a Gaussian variable but with a non-zero average

$$\overline{i|𝒯|i}=\underset{m,n}{}𝒯_{m,n}\overline{\psi _i^{}(m)\psi _i(n)}=\frac{1}{N_i}\mathrm{Tr}\left(P_i𝒯\right).$$
(8)
Using Wick’s theorem and the fact that the GOE wave functions $`\psi _i(m)`$ are real, the variance of a diagonal element is

$$\sigma ^2(i|𝒯|i)=\underset{\genfrac{}{}{0pt}{}{m,n}{m^{},n^{}}}{}𝒯_{m,n}^{}𝒯_{m^{},n^{}}\left[\overline{\psi _i(m)\psi _i^{}(m^{})}\overline{\psi _i(n)\psi _i^{}(n^{})}+\overline{\psi _i(m)\psi _i^{}(n^{})}\overline{\psi _i(n)\psi _i^{}(m^{})}\right].$$
(9)
For a real and hermitean operator $`𝒯`$, we find

$$\sigma ^2(i|𝒯|i)=\frac{2}{N_i^2}\mathrm{Tr}\left(P_i𝒯P_i𝒯\right),$$
(10)
which is a factor of two larger than the variance of an off-diagonal element (see Eq. (7) for the case when the initial and final states span the same subspace of states, i.e., $`P_i=P_f`$).

Since the mean value of $`i|𝒯|i`$ is non-zero, its square does not follow a Porter-Thomas distribution. However, once the mean value (8) has been subtracted from the diagonal matrix elements, their squares $`z`$ are predicted to have a Porter-Thomas distribution. The distributions of the squares of these shifted reduced diagonal matrix elements are shown in Fig. 4 for the $`M1`$ magnetic moments of the $`J=2^+,T=1`$ states in $`A=60`$ nuclei. The values of $`z`$ are renormalized by dividing out the secular variation with energy of the square of the reduced diagonal matrix element, defined by an equation similar to (3) but with a single Gaussian (since $`i=f`$). Considering the small number of data points used to compute the distribution $`P(\mathrm{log}_{10}z)`$ (there are only 66 diagonal matrix elements), the agreement with a Porter-Thomas distribution (solid lines) is reasonable.

It is of interest to see the effect that the inclusion of diagonal matrix elements has on the distribution of transition intensities. Since the number of diagonal elements is much smaller than the number of transition intensities, we do not expect them to modify the Porter-Thomas statistics of the transition intensities (even though the squares of the unshifted diagonal reduced elements do not follow a Porter-Thomas distribution). Fig. 4 shows the same distributions as in Fig. 4 but when the squares of the reduced diagonal elements of the corresponding transition operator are included in the analysis. The $`B(E2)`$ distributions as well as the $`B(M1)`$ distribution in the $`T_z=1`$ nucleus are all well described by a Porter-Thomas distribution, but the $`B(M1)`$ distribution in the self-conjugate nucleus ($`T_z=0`$) deviates significantly from the Porter-Thomas distribution (right bottom panel of Fig. 4). The $`M1`$ transitions in self-conjugate nuclei ($`T_z=0`$) have the property that they are purely isoscalar. In $`T_z=1`$ nuclei both isoscalar and isovector components contribute to the $`M1`$ transitions, but the isoscalar $`M1`$ matrix elements are much weaker than the isovector $`M1`$ , and the latter dominate the distributions. To check that our conclusion does not depend on the size of the model space, we performed similar calculations in a model space with 6 active particles, i.e., for $`A=62`$ nuclei. Again, we observed that the distributions of the $`M1`$ transitions in the self-conjugate nucleus $`T_z=0`$ (i.e., 62Ga) deviate from a Porter-Thomas distribution.

The deviation from a Porter-Thomas distribution for the isoscalar $`M1`$ can be traced to the diagonal matrix elements. Although their number is much smaller than the number of off-diagonal elements, they are in magnitude much larger than the off-diagonal elements. Indeed, the isoscalar $`M1`$ has a particularly simple structure and can be written (in units of nuclear magneton $`\mu _N`$) as

$$\frac{g_p^{\mathrm{}}+g_n^{\mathrm{}}}{2}𝐋+\frac{g_p^s+g_n^s}{2}𝐒,$$
(11)
where $`𝐋`$ and $`𝐒`$ are the total orbital angular momentum and spin, respectively, and $`g_{p,n}^{\mathrm{}}`$, $`g_{p,n}^s`$ are the proton and neutron orbital and spin $`g`$ factors. Since the transition matrix elements of the total angular momentum $`𝐉=𝐋+𝐒`$ between different states vanish, the orbital and spin contributions to the isoscalar $`M1`$ matrix elements always have opposite signs. Consequently, the off-diagonal elements are small. However the interference between the orbital and spin contributions to the diagonal matrix elements can be constructive, leading to larger diagonal matrix elements. These diagonal matrix elements are related to the magnetic moments, which in $`T_z=0`$ nuclei are given by $`J/2+0.38S_z`$ . The mean value (8) of the isoscalar magnetic moment is large (relative to a typical off-diagonal element) and leads to the observed deviation from a Porter-Thomas distribution. In the case of the isovector $`M1`$ matrix elements, the relative phase of the orbital and spin contributions is more random , and the corresponding statistics are Porter-Thomas even when diagonal elements are included.

In conclusion, we have studied $`\mathrm{\Delta }T=0`$ $`B(E2)`$ and $`B(M1)`$ transition strength distributions in $`pf`$-shell nuclei with $`A60`$ and found that they agree well with Porter-Thomas statistics. We have also studied the statistics of diagonal elements (e.g., $`M1`$ magnetic moments) in these nuclei. The distribution of the squares of these moments is found to be in fair agreement with a Porter-Thomas distribution, but only when the mean value is subtracted from the moments. In general, the inclusion of off-diagonal elements in the statistics of transition intensities does not affect the observed Porter-Thomas statistics. An exception is the isoscalar $`M1`$ operator whose diagonal matrix elements are much larger than its off-diagonal elements.

This work was supported in part by the Department of Energy grant No. DE-FG-0291-ER-40608 and by the Russian Foundation for Basic Research, project No. 00-02-17194.


text:
# Parametric Resonance versus Forced Oscillation in Time-Evolution of Quantum Meson Fields

It is interesting to study the dynamics of the chiral phase transition in the context of relativistic heavy-ion collisions, such as the problem of the formation of a disoriented chiral condensate. One of the important theoretical aspects is to investigate the time-evolution of the order parameter of the chiral phase transition. As for the fluctuation modes around the order parameter, the amplification of the amplitudes must occur accompanying with the relaxation of the chiral order parameter. Recently, the present author has investigated the time-evolution of a collective meson field in the context of the dynamical chiral phase transition. It has been seen that the amplitudes of quantum fluctuations with low momenta around the mean field configuration have been amplified with the direction of pion modes in the $`O(4)`$-linear sigma model. As is indicated by many authors, this phenomena may be understood in terms of a parametric amplification. However, there may be another mechanism to amplify the fluctuation modes. Actually, it is seen that the amplitude of the pion mode with higher momentum without $`𝒌=0`$ can be amplified in the previous paper. In this paper, it is pointed out that there is another mechanism to amplify the quantum fluctuation modes through the chiral phase transition, namely, that a forced oscillation works as well as a parametric resonance.

Before treating the $`O(4)`$-linear sigma model, it is instructive to recapitulate the ingredients of the parametric resonance and the forced oscillation. Let us consider the following equation of motion for $`x(t)`$ :

$$\ddot{x}+\omega _0^2(1h\mathrm{cos}\gamma t)x=0.$$
(1)
The parametric resonance occurs around $`\gamma =2\omega _0/m(m=1,2,\mathrm{})`$. For $`m=1`$, the resonance region of the frequency $`\gamma `$ is given as

$`h\omega /2<ϵ<h\omega /2,\text{for}\gamma =2\omega _0+ϵ,ϵ\omega _0,`$ (2)
where we treated $`h`$ as a small parameter, $`h1`$. Similarly, for $`m=2`$, we obtain

$`(5/24)h^2\omega _0<ϵ<(1/24)h^2\omega _0,\text{for}\gamma =\omega _0+ϵ,ϵ\omega _0.`$ (3)
In general, for given integer $`m`$, the unstable solution of $`x(t)`$ is obtained as the form $`x(t)e^{st}\times `$(oscillation and constant parts), where $`s`$ is of the order of $`h^m`$. Namely, the amplification becomes slower in time as the integer $`m`$ becomes larger. The resonance region also becomes narrow with the order of $`h^m`$, that is, $`O(h^m)<ϵ<O(h^m)`$. If $`h`$ is not so small, we have to deal with the equation (1) directly. The equation (1) is known as Mathieu’s equation and the properties of this equation are investigated in detail. In the following investigation, $`h`$ is actually small fortunately. Thus, it is not necessary to take care of the full treatment of this Mathieu’s equation. Next, let us consider the equation of motion which describes the forced oscillation :

$$\ddot{x}+\omega _0^2x=F\mathrm{cos}\mathrm{\Omega }t.$$
(4)
The solution with initial conditions $`x(0)=\dot{x}(0)=0`$ is obtained as

$$x(t)=F\frac{\mathrm{cos}\mathrm{\Omega }t\mathrm{cos}\omega _0t}{\omega _0^2\mathrm{\Omega }^2}=\frac{2F}{(\omega _0\mathrm{\Omega })(\omega _0+\mathrm{\Omega })}\mathrm{sin}\frac{\omega _0\mathrm{\Omega }}{2}t\mathrm{sin}\frac{\omega _0+\mathrm{\Omega }}{2}t.$$
(5)
As is well known, when $`\mathrm{\Omega }`$ has a value near $`\omega _0`$, the behavior of the above solution shows the beat. Especially, $`\mathrm{\Omega }\omega _0`$, the amplitude of the solution becomes large linearly in time :

$$x(t)(F/2\omega _0)t\mathrm{sin}\omega _0t$$
(6)
Now let us return to the main subject, namely, let us treat the $`O(4)`$-linear sigma model in the squeezed state approach. We start from the following Hamiltonian density of $`O(4)`$-linear sigma model:

$$=\frac{1}{2}\pi _a(𝒙)^2+\frac{1}{2}\varphi _a(𝒙)\varphi _a(𝒙)+\lambda \left(\varphi _a(𝒙)^2m^2/4\lambda \right)^2H\varphi _0(𝒙).$$
(7)
The squeezed state is adopted as a trial state in the time-dependent variational principle :

$`|\mathrm{\Phi }(t)=N\mathrm{exp}\left(i(\overline{\pi }\varphi \overline{\phi }\pi )\right)\mathrm{exp}(\varphi [(G^1G^{(0)}{}_{}{}^{1})/4+i\mathrm{\Sigma }]\varphi )|0,`$ (8)
where $`N`$ represents a normalization factor and $`\overline{\pi }\varphi =_{a=0}^3d^3𝒙\overline{\pi }_a(𝒙,t)\varphi _a(𝒙)`$ and $`\varphi G^1\varphi =_{a=0}^3d^3𝒙d^3𝒚\varphi _a(𝒙)G_a^1(𝒙,𝒚,t)\varphi _a(𝒚)`$ and so on. Here, $`G_a^{(0)}=0|\varphi _a(𝒙)\varphi _a(𝒚)|0`$. The expectation values of the field operators are obtained, for example, as

$`\mathrm{\Phi }(t)|\varphi _a(𝒙)|\mathrm{\Phi }(t)=\overline{\phi }_a(𝒙,t),`$
$`\mathrm{\Phi }(t)|\varphi _a(𝒙)\varphi _a(𝒚)|\mathrm{\Phi }(t)=\overline{\phi }_a(𝒙,t)\overline{\phi }_a(𝒚,t)+G_a(𝒙,𝒚,t).`$ (9)
$`\overline{\phi }_a(𝒙,y)`$ and $`G_a(𝒙,𝒚,t)`$ represent the mean field and the fluctuation around it, respectively. The time-dependence of these functions, together with $`\overline{\pi }_a`$ and $`\mathrm{\Sigma }_a`$, are determined in the time-dependent variational principle :

$$\delta _{t_1}^{t_2}𝑑t\mathrm{\Phi }(t)|i\frac{}{t}\widehat{H}|\mathrm{\Phi }(t)=0,$$
(10)
where $`\widehat{H}=d^3𝒙`$. The two-point function $`G_a(𝒙,𝒚,t)`$ can be further expanded in terms of the mode functions under the assumption of the translational invariance :

$$G_a(𝒙,𝒚,t)=_𝐤e^{i𝐤(𝐱𝐲)}\eta _𝐤^a(t)^2,$$
(11)
where $`_𝐤=d^3𝐤/(2\pi )^3`$. The equations of motion for the condensate $`\overline{\phi }_a(t)`$ and the fluctuation modes $`\eta _𝐤^a(t)`$ under the translational invariance can be expressed as

$`\ddot{\overline{\phi }}_a(t)m^2\overline{\phi }_a(t)+4\lambda \overline{\phi }_a(t)^3+12\lambda {\displaystyle _𝐤}\eta _𝐤^a(t)^2\overline{\phi }_a(t)`$
$`+4\lambda {\displaystyle \underset{ba}{}}\left(\overline{\phi }_b(t)^2+{\displaystyle _𝐤}\eta _𝐤^b(t)^2\right)\overline{\phi }_a(t)H\delta _{a0}=0,`$
$`\ddot{\eta }_𝐤^a(t)+[𝒌^2m^2+12\lambda \overline{\phi }_a(t)^2+12\lambda {\displaystyle _𝐤^{}}\eta _𝐤^{}^a(t)^2`$
$`+4\lambda {\displaystyle \underset{ba}{}}(\overline{\phi }_b(t)^2+{\displaystyle _𝐤^{}}\eta _𝐤^{}^b(t)^2)]\eta _𝐤^a(t)1/4\eta _𝐤^a(t)^3=0.`$ (12)
The derivation and the physical viewpoint in this approach are given in detail in Ref.?. Hereafter, let us assume that the chiral condensate points in the sigma direction, that is, $`\overline{\phi }_00`$ and $`\overline{\phi }_i=0`$ for $`i=13`$. Also, fluctuation modes of the pi directions are identical one another, which are denoted as $`\eta _𝐤^\pi (t)`$, that is, $`\eta _𝐤^1=\eta _𝐤^2=\eta _𝐤^3\eta _𝐤^\pi `$. When the explicit chiral symmetry breaking term, $`H`$, is small, the static solutions in Eq.(Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields) are given as

$$\overline{\phi }_0\phi _0=\sqrt{M_\sigma ^2/8\lambda }H/2M_\sigma ^2,\eta _{𝐤}^{0}{}_{}{}^{2}=1/2\sqrt{𝒌^2+M_\sigma ^2},\eta _{𝐤}^{\pi }{}_{}{}^{2}=1/2\sqrt{𝒌^2+M_\pi ^2},$$
(13)
where the sigma meson mass $`M_\sigma `$ and the pion mass $`M_\pi `$ are defined as

$`M_\sigma ^2=m^2+12\lambda \phi _0^2+12\lambda {\displaystyle _𝐤^{}}\eta _{𝐤^{}}^{0}{}_{}{}^{2}+12\lambda {\displaystyle _𝐤^{}}\eta _{𝐤^{}}^{\pi }{}_{}{}^{2},`$
$`M_\pi ^2=m^2+4\lambda \phi _0^2+4\lambda {\displaystyle _𝐤^{}}\eta _{𝐤^{}}^{0}{}_{}{}^{2}+20\lambda {\displaystyle _𝐤^{}}\eta _{𝐤^{}}^{\pi }{}_{}{}^{2}.`$ (14)
Let us examine the time-dependent solutions around the static configurations. In Eq.(Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields), the time-dependent variables can be expanded around the static solutions of Eq.(13) :

$`\overline{\phi }_0(t)=\phi _0+\delta \phi (t),\eta _𝐤^a(t)=\eta _𝐤^a+\delta \eta _𝐤^a(t).`$ (15)
Here, we consider the late time of the chiral phase transition. Thus, $`\delta \phi (t)`$ is small compared with $`\phi _0`$. Further, we assume that $`|\delta \eta _𝐤^a(t)||\eta _𝐤^a|`$ :

$`|\delta \phi (t)|/\phi _01,|\delta \eta _𝐤^a(t)/\eta _𝐤^a|1.`$ (16)
From Eq.(Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields), it is seen that the two-point function $`G_a`$ represents quantum fluctuations around the mean field. Namely, $`G_a`$ is small compared with $`\phi _0^2`$. Thus, from Eq.(11), we conclude that the following relation is satisfied :

$`\phi _0^2{\displaystyle _𝐤}\eta _𝐤^{a2},\mathrm{then},\left|\phi _0\delta \phi (t)\right|\left|{\displaystyle _𝐤}\eta _𝐤^a\delta \eta _𝐤^a(t)\right|.`$ (17)
Under the above approximations and the small explicit chiral symmetry breaking, the equations of motion for $`\delta \phi (t)`$ and $`\delta \eta _𝐤^a(t)`$ are obtained from (Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields) and (15) as

$`\delta \ddot{\phi }(t)+M_\sigma ^2\delta \phi (t)=0,`$ (18)
$`\delta \ddot{\eta }_𝐤^\sigma (t)+\left[4(𝒌^2+M_\sigma ^2)+24\lambda \phi _0\delta \phi (t)\right]\delta \eta _𝐤^\sigma (t)=24\lambda \eta _𝐤^\sigma \phi _0\delta \phi (t),`$
$`\delta \ddot{\eta }_𝐤^\pi (t)+\left[4(𝒌^2+M_\pi ^2)+8\lambda \phi _0\delta \phi (t)\right]\delta \eta _𝐤^\pi (t)=8\lambda \eta _𝐤^\pi \phi _0\delta \phi (t),`$ (19)
where $`\eta _𝐤^0`$ and $`\delta \eta _𝐤^0(t)`$ have been rewritten as $`\eta _𝐤^\sigma `$ and $`\delta \eta _𝐤^\sigma (t)`$, respectively. The equation (18) is easily solved by using the sigma meson mass :

$$\delta \phi (t)=\delta \sigma \mathrm{cos}(M_\sigma t+\alpha ),$$
(20)
where $`\delta \sigma `$ and $`\alpha `$ are determined by the initial conditions. Here, we adopt $`\alpha =\pi `$ for simplicity. We further define the dimensionless variables :

$`\delta \stackrel{~}{\eta }_𝐤^a(\tau )=\sqrt{2M_a}\delta \eta _𝐤^a(t),\tau =M_\sigma t/2,`$ (21)
$`\omega _\sigma =4\sqrt{𝐤^2+M_\sigma ^2}/M_\sigma ,\omega _\pi =4\sqrt{𝐤^2+M_\pi ^2}/M_\sigma ,`$
$`h_\sigma =6\lambda \phi _0^2/(𝐤^2+M_\sigma ^2)\delta \sigma /\phi _0,h_\pi =2\lambda \phi _0^2/(𝐤^2+M_\pi ^2)\delta \sigma /\phi _0,`$
$`F_\sigma =96{\displaystyle \frac{\lambda \phi _0^2}{M_\sigma ^2}}\left({\displaystyle \frac{M_\sigma ^2}{𝐤^2+M_\sigma ^2}}\right)^{1/4}{\displaystyle \frac{\delta \sigma }{\phi _0}},F_\pi =32{\displaystyle \frac{\lambda \phi _0^2}{M_\sigma ^2}}\left({\displaystyle \frac{M_\pi ^2}{𝐤^2+M_\pi ^2}}\right)^{1/4}{\displaystyle \frac{\delta \sigma }{\phi _0}}.`$ (22)
By introducing the above dimensionless variables, the solution (20) and the equations of motion (19) are recast into the simple forms as

$`\delta \phi (\tau )=\delta \sigma \mathrm{cos}\gamma \tau ,\gamma =2,`$ (23)
$`\left({\displaystyle \frac{d^2}{d\tau ^2}}+\omega _\sigma ^2[1h_\sigma \mathrm{cos}\gamma \tau ]\right)\delta \stackrel{~}{\eta }_𝐤^\sigma (\tau )=F_\sigma \mathrm{cos}\gamma \tau ,`$ (24)
$`\left({\displaystyle \frac{d^2}{d\tau ^2}}+\omega _\pi ^2[1h_\pi \mathrm{cos}\gamma \tau ]\right)\delta \stackrel{~}{\eta }_𝐤^\pi (\tau )=F_\pi \mathrm{cos}\gamma \tau .`$ (25)
If $`F_\sigma `$ ($`F_\pi `$) is negligible, then the equation (24) ((25)) is reduced to Eq.(1). In this case, the existence of the unstable solutions for $`\delta \stackrel{~}{\eta }_𝐤^a(\tau )`$ may be expected. On the other hand, if $`h_\sigma `$ ($`h_\pi `$) is negligible or even if the parameters do not offer the unstable regions, the forced oscillation may be realized by the effect of $`F_\sigma `$ ($`F_\pi `$).

In the previous paper, the behavior of the time-dependent solutions in Eq.(Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields) was investigated numerically. It was shown that the amplification of the amplitude of the pion modes with low momenta was realized even in the late time. However, it was seen that the amplitude of the second excited mode was also amplified (See, Fig.3 in Ref.?). We can now understand their behaviors by using Eqs.(24) and (25) approximately. The model parameters $`m`$, $`\lambda `$ and $`H`$ are adopted so as to reproduce the pion mass, the sigma meson mass and the pion decay constant in the static case in Eqs.(13) and (Parametric Resonance versus Forced Oscillation in  Time-Evolution of Quantum Meson Fields). The box normalization with the spatial length $`L`$ ($`=10`$ fm) is applied, in which the periodic boundary conditions for the fluctuation modes are imposed. Thus, the allowed values of momenta are restricted as follows : $`k_x=(2\pi /L)n_x`$ and so on, where $`n_x`$ is integer. The three momentum cutoff is also applied to regularize the momentum integration. The fluctuation modes are taken into account up to $`n^2=n_x^2+n_y^2+n_z^28^2`$ corresponding to the three momentum cutoff about 1 GeV (990 MeV). The values of variables given in Eq.(22) are collected up to $`n=4`$ in Table I., where $`\delta =\delta \sigma /\phi _0O(10^1)`$ in the late time of the time-evolution of the chiral condensate.

First, let us consider the pion modes. It should be noted that $`\gamma =2`$. For the lowest pion mode with $`n=0`$, $`\omega _\pi =0.92`$. Therefore, it would not expected that the amplification due to the forced oscillation occurs. We thus investigate the possibility of the parametric resonance by setting $`F_\pi 0`$ safely. By noting $`2\omega _\pi 2`$, let us examine the parametric resonance with $`m=1`$ in the case Eq.(2). The unstable region is given as

$$h_\pi \omega _\pi /2<ϵ<h_\pi \omega _\pi /2\text{with}\gamma =2=2\omega _\pi +ϵ,(ϵ0.16).$$
(26)
The condition (26) reads together with Table I. as

$$\delta \delta \sigma /\phi _0>0.072.$$
(27)
Thus, in the above region for $`\delta \sigma `$, the parametric resonance occurs and the amplitude of the lowest pion mode is amplified. This behavior appears in Fig.3 in Ref.?. Similarly, let us consider the pion mode with $`n=1`$, the first excited mode. As is similar to the case $`n=0`$, it is not necessary to take care of the amplification due to the forced oscillation because $`\omega _\pi =1.24`$. For $`m=1`$, the same condition as (26) should be satisfied to realize the parametric resonance, except for the value of $`ϵ`$ : $`ϵ=2(1\omega _\pi )0.48`$. Thus, the unstable region is rewritten as

$$\delta \delta \sigma /\phi _0>0.29.$$
(28)
The parametric resonance is then realized. For $`m=2`$, the unstable region is given as

$$(5/24)h_\pi ^2\omega _\pi <ϵ<(1/24)h_\pi ^2\omega _\pi \mathrm{with}\gamma =2=\omega _\pi +ϵ(ϵ0.76).$$
(29)
However, $`(1/24)h_\pi ^2\omega _\pi 0.37(\delta \sigma /\phi _0)^2`$ is obtained in this parameterization. Therefore, the parametric resonance with $`m=2`$ does not occur.

Secondly, let us investigate the pion mode with $`n=2`$, the second excited mode. If the parametric amplification occur, the case $`m=2`$ is realized because $`\omega _\pi =1.89`$. This condition is given in Eq.(29) except for $`ϵ`$ : $`ϵ0.11`$. However, the upper limit $`(1/24)h_\pi ^2\omega _\pi `$ gives a severe limit : $`0.10(\delta \sigma /\phi _0)^2`$. Thus, the condition (29) can not be satisfied. On the other hand, it is possible that the amplitude of the pion mode with $`n=2`$ is amplified by the mechanism of the forced oscillation. The value of $`\omega _\pi (=1.89)`$ is near $`\gamma (=2)`$. If we neglect $`h_\pi `$, the equation (25) is same as (4) and the solution is given in Eq.(6). Thus, we obtain

$$\delta \stackrel{~}{\eta }_{n=2}^\pi (\tau )0.71(\delta \sigma /\phi _0)\tau \mathrm{sin}\omega _\pi \tau .$$
(30)
Thus, the amplitude increases linearly in time. This behavior is seen in Fig.3 in Ref.?.

Finally, let us consider the sigma mode. The amplification due to the forced oscillation is impossible because the frequency $`\omega _\sigma `$ is always greater than 2 $`(=\gamma )`$. Here, let us investigate the possibility of the parametric amplification. The lowest mode with $`n=0`$ is the most powerful candidate. Since $`\omega _\sigma =4`$, the case $`m=4`$ is applicable : $`\gamma =2=\omega _\sigma /2+ϵ`$, where $`ϵ=0`$. The unstable region is given as $`O(h_\sigma ^4)<ϵ<O(h_\sigma ^4)`$ for $`m=4`$ case. Since $`ϵ=0`$, the parametric amplification occurs. However, the unstable solution is proportional to $`e^{s\tau }`$, in which $`s`$ is of the order of $`h_\sigma ^4`$ for $`m=4`$ case. From Table I., $`h_\sigma ^4`$ is about $`3.4\times 10^5`$ with $`\delta \sigma /\phi _0=10^1`$. This is very small value. Thus, the amplification occurs slowly in time. In the time region under consideration, the amplification may not be seen explicitly. In the sigma modes with higher momenta, the parametric resonance would not occur because the unstable window is very narrow.

It should be noted here that the mechanism amplifying the amplitudes of fluctuation modes, i.e., the parametric resonance or the forced oscillation, depends on the model parameters. However, in conclusion, it should be pointed out that there is a possibility of the amplification due to the forced oscillation as well as the parametric resonance. In this paper, the parametric resonance and the forced oscillation are treated separately in Eqs.(24) and (25). It is interesting to investigate both the mechanisms in the unified way. The important ingredients may be the properties of Mathieu’s equation and the Mathieu function.

This work was supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology (No.13740159).


text:
# Generalised Wigner surmise for (2×2) random matrices

## 1 Introduction

The level statistics of a quantum system represents the most significant, although not the only , signature of quantum chaos. The Poisson and Wigner distributions of dimensionless nearest-neighbour spacing, $`s`$,

$$\stackrel{~}{p}_P(s)=\mathrm{exp}(s),$$
(1)
$$\stackrel{~}{p}_W(s)=\frac{\pi s}{2}\mathrm{exp}\left(\pi s^2/4\right),$$
(2)
are known in quantum chaos theory as two universalities that correspond to two extreme cases of classical dynamics, namely purely regular and completely chaotic (see , e.g., ). The majority of many-body systems such as nuclei, molecules, atoms or solids (see  and references therein) have been found to be chaotic although for such complex systems no classical limit can be constructed.

As was first recognised by Wigner , the nearest-neighbour spacing distribution (NNSD) (2) well corresponds to the eigenvalue distributions of random matrices, and this explains the importance of the Random Matrix Theory  for studying statistical properties of many-body systems. Particular attention was paid to Gaussian ensembles. In fact, assuming that (i) the elements of the Hamiltonian matrix are independent real variables and (ii) the matrix distribution is invariant under an orthogonal transformation of the basis states (see, e.g., chapter 3 in ), one finds that the matrix elements are independent Gaussian variables with zero mean and with variance satisfying the conditions $`\sigma _{ij}^2=(1+\delta _{ij})\sigma ^2`$. Imposing particular symmetries on the Hamiltonian, one gets  Gaussian Orthogonal, Unitary or Symplectic Ensembles (GOE, GUE or GSE, respectively), which are widely and successfully applied in many fields of physics (see e.g. reviews ).

Up to now, excited atomic nuclei are considered to be the best examples of chaotic quantum systems . Starting from slow-neutron scattering experiments, which were first described in terms of random matrices , a lot of nuclear structure data has been analysed in the context of chaos. However, it was repeatedly noticed that the experimental data do not exactly match the distribution (2) and exhibit slight deviations . These deviations are thought to be caused by the fact that the real system is not purely chaotic, but can be the quantum analog of a classical system that is transitional between chaotic and integrable. In this context, a few phenomenological formulae were proposed and analysed. The most famous are the Brody  distribution,

$$\stackrel{~}{p}_\omega (s)=(\omega +1)\alpha s^\omega \mathrm{exp}\left(\alpha s^{\omega +1}\right),\alpha =\left[\mathrm{\Gamma }\left(\frac{\omega +2}{\omega +1}\right)\right]^{\omega +1},$$
(3)
which was shown to match better the experimental data  on both high- and low-energy nuclear spectra, and the Berry-Robnik distribution

$$\stackrel{~}{p}_{BR}(s)=e^{(q1)s}\{(1q)^2\text{erfc}(\sqrt{\pi }qs/2)+[2q(1q)+(\pi /2)q^3s]e^{(\pi /4)q^2s^2}\}.$$
(4)
Although the Brody distribution takes the form of Poisson for $`\omega =0`$ and Wigner for $`\omega =1`$, it has the unrealistic property that its derivative goes to infinity at $`\epsilon =0`$ . Moreover, the parameter $`\omega `$ has no clear physical meaning. On the other hand, the distribution (4) does not give a level-repulsion for non-integrable systems. Caurier et al  considered a model which allowed to simulate the transition from integrability to chaos and succeeded to derive the asymptotic limits for small and large neighbour spacings in the near-integrable limit.

The idea to derive a third universality class corresponding to intermediate statistics has been actively pursued in recent years. In particular, in Refs. , the distribution of $`n`$ nearest neighbours,

$$\stackrel{~}{p}^{(\beta )}(n,s)=\frac{(\beta +1)^{n(\beta +1)}}{\mathrm{\Gamma }[n(\beta +1)]}s^{[n(\beta +1)1]}\mathrm{exp}\left[(\beta +1)s\right],$$
(5)
was shown to be relevant to a certain class of exactly solvable models with nearest and next-to-nearest neighbour interactions, generalizing the results obtained earlier  on pseudointegrable billiards and the short-range Dyson models. Remarkably, the NNSD given by (5) with $`n=1`$ exhibits a level repulsion $`s^\beta `$ and falls to zero at $`s`$ as $`\mathrm{exp}[(\beta +1)s]`$. For $`\beta =1`$ it is referred to as the semi-Poisson distribution .

In this letter we argue that one might search for an explanation of the discrepancy between data and random matrix theory simply by generalising the random matrix ensemble used. In fact, the Gaussian ensemble is defined by the two assumptions (i) and (ii) mentioned above, and their applicability should be carefully checked for a given physical system. As far as nuclear physics is concerned, the deviations of the experimental data on slow-neutron or $`(p,p^{})`$ resonances from a random matrix description could simply arise from the non-invariance of the random matrix ensemble under an orthogonal transformation of a basis, i.e. the assumption (ii) is violated.

In addition, it is well known that realistic interactions in many-body nuclear, molecular or atomic systems are predominantly of one- and two-body nature, implying that the distribution of the matrix is not only not invariant under an orthogonal (unitary) transformation of the basis, but also that the elements of the Hamiltonian matrix are not independent (if the number of particles is more than two), i.e. both assumptions for a Gaussian ensemble do not hold any longer. In this context, French and Wong  and Bohigas and Flores  independently introduced the Two-Body Random Ensemble (TBRE), which is characterised by a Gaussian level density distribution, rather than a semi-circle provided by a GOE . The level distributions of experimental nuclear spectra did confirm this result. The NNSD relevant for a TBRE was found numerically to be fairly represented by the Wigner surmise , although recently it has been pointed out  that NNSD given by (5) with $`n=1`$ and a certain real value of $`\beta `$ fits better the shell–model spectrum obtained with the $`sd`$-interaction of Wildenthal .

Given the importance and actuality of these investigations, we present in this letter some analytical results concerning the properties of (2$`\times `$2) random symmetric matrices for which the assumption (ii) mentioned above is not satisfied. First, we derive the Hamiltonian distribution as a function of its eigenvalues and we calculate the NNSD which generalises the well-known Wigner surmise . We show that the model allows to describe the transition from purely chaotic to asymptotically nearly integrable limits, being different from the intermediate statistics mentioned above. Then, for a particular case, we give the analytical expressions for the lowest moments of this distribution. Finally, we propose a method to derive the moments of the eigenvalue distribution without knowledge of an explicit expression for the distribution.

## 2 Generalised Wigner surmise

Let us consider a (2$`\times `$2) real symmetric matrix

$$\begin{array}{cc}H=\left(\begin{array}{cc}H_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right),& \end{array}$$
(6)
whose elements are independent Gaussian variables with zero mean and variance $`\sigma _{ij}^2`$, and $`H_{12}=H_{21}`$. The probability density of the matrix $`H`$ is then given by

$$p(H)=\frac{1}{(2\pi )^{3/2}\sqrt{\sigma _{11}^2\sigma _{12}^2\sigma _{22}^2}}\mathrm{exp}\left[\left(\frac{H_{11}^2}{2\sigma _{11}^2}+\frac{H_{12}^2}{2\sigma _{12}^2}+\frac{H_{22}^2}{2\sigma _{22}^2}\right)\right].$$
(7)
Each matrix $`H`$ can be diagonalised in an orthogonal basis and therefore $`H=O^tDO`$, with

$$O=\begin{array}{cc}\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right),& \end{array}$$

and

$$D=\begin{array}{cc}\left(\begin{array}{cc}E_\alpha & 0\\ 0& E_\beta \end{array}\right).& \end{array}$$

Similar to the case of GOE , we find that in the general case

$$\begin{array}{cc}\begin{array}{ccc}H_{11}& \text{=}& E_\alpha \mathrm{cos}^2\theta +E_\beta \mathrm{sin}^2\theta \\ H_{12}& \text{=}& (E_\alpha E_\beta )\mathrm{cos}\theta \mathrm{sin}\theta \\ H_{22}& \text{=}& E_\alpha \mathrm{sin}^2\theta +E_\beta \mathrm{cos}^2\theta \end{array}& \end{array}.$$

We deduce that the probability density expressed in terms of the eigenvalues and the angle $`\theta `$ is

$`p(E_\alpha ,E_\beta ,\theta )`$ $`=`$ $`{\displaystyle \frac{E_\alpha E_\beta }{(2\pi )^{3/2}\sqrt{\sigma _{11}^2\sigma _{22}^2\sigma _{12}^2}}}\mathrm{exp}\left\{{\displaystyle \frac{\left[E_\alpha \mathrm{\Sigma }^2(E_\alpha E_\beta )\left(\sigma _{11}^2\mathrm{cos}^2\theta +\sigma _{22}^2\mathrm{sin}^2\theta \right)\right]^2}{2\sigma _{11}^2\sigma _{22}^2\mathrm{\Sigma }^2}}\right\}`$ (8)
$`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}(E_\alpha E_\beta )^2\left({\displaystyle \frac{\mathrm{cos}^2(2\theta )}{\mathrm{\Sigma }^2}}+{\displaystyle \frac{\mathrm{sin}^2(2\theta )}{4\sigma _{12}^2}}\right)\right]`$
where $`\mathrm{\Sigma }^2=\sigma _{11}^2+\sigma _{22}^2`$ and $`E_\alpha E_\beta 0`$.

The nearest-neighbour spacing distribution for the variable $`\epsilon =E_\alpha E_\beta =Ds`$, with the mean spacing $`D=\epsilon \stackrel{~}{p}\left(\epsilon \right)𝑑\epsilon `$, is given by the following integral

$$\stackrel{~}{p}\left(\epsilon \right)=_{\pi /2}^{\pi /2}𝑑\theta _{\mathrm{}}^{\mathrm{}}𝑑E_\alpha _{\mathrm{}}^{E_\alpha }𝑑E_\beta p(E_\alpha ,E_\beta ,\theta )\delta \left(\epsilon E_\alpha +E_\beta \right),$$
(9)
from which we obtain

$$\stackrel{~}{p}\left(\epsilon \right)=\frac{\epsilon }{2\sqrt{\mathrm{\Sigma }^2\sigma _{12}^2}}\mathrm{exp}\left[\frac{\epsilon ^2\left(\mathrm{\Sigma }^2+4\sigma _{12}^2\right)}{16\mathrm{\Sigma }^2\sigma _{12}^2}\right]I_0\left(\frac{\epsilon ^2\left(\mathrm{\Sigma }^24\sigma _{12}^2\right)}{16\mathrm{\Sigma }^2\sigma _{12}^2}\right)$$
(10)
where $`I_0`$ is a modified Bessel function of the first kind.

The expression (10) looks like a Rayleigh-Rice distribution, well known in signal theory , except for the argument of $`I_0`$, which is not linear as in the usual Rayleigh-Rice distribution but quadratic. This is why we will refer to $`\stackrel{~}{p}\left(\epsilon \right)`$ as to a quadratic Rayleigh-Rice distribution.

Let us consider a particular case when the diagonal matrix elements have the same variance $`\sigma _{11}^2=\sigma _{22}^2`$, which is $`\chi `$ times larger than the variance, $`\sigma _{12}^2=\sigma ^2`$, of the non-diagonal matrix elements, i.e. $`\chi =\sigma _{11}^2/\sigma _{12}^2`$. Then the eigenvalue distribution (8) reduces to

$$p_\chi (E_\alpha ,E_\beta ,\theta )=\frac{E_\alpha E_\beta }{(2\pi \sigma ^2)^{3/2}\chi }\mathrm{exp}\left[\frac{E_\alpha ^2+E_\beta ^2+\frac{1}{4}(E_\alpha E_\beta )^2(\chi 2)\mathrm{sin}^2(2\theta )}{2\chi \sigma ^2}\right],$$
(11)
while for the nearest-neighbour spacing we get

$$\stackrel{~}{p}_\chi (\epsilon )=\frac{\epsilon }{\sqrt{2\chi }2\sigma ^2}\mathrm{exp}\left(\frac{\left(\chi +2\right)\epsilon ^2}{16\chi \sigma ^2}\right)I_0\left(\frac{(\chi 2)\epsilon ^2}{16\chi \sigma ^2}\right).$$
(12)
For $`\chi =2`$, expression (12) reduces to the Wigner surmise. The distributions $`\stackrel{~}{p}_\chi (\epsilon )`$ are plotted in figure 1 for $`\chi =1`$, $`\chi =2`$ and $`\chi =5`$.

For small values of $`\epsilon `$, the distribution (12) goes linearly to zero,

$$\stackrel{~}{p}_\chi (\epsilon )\frac{1}{\sqrt{2\chi }2\sigma ^2}\epsilon ,$$
(13)
while for large $`\epsilon `$,

$$\stackrel{~}{p}_\chi (\epsilon )\{\begin{array}{cc}\sqrt{\frac{2}{\chi 2\sigma ^2}}\mathrm{exp}\left(\frac{\epsilon ^2}{4\chi \sigma ^2}\right)\hfill & \chi 2\hfill \\ \frac{\epsilon }{4\sigma ^2}\mathrm{exp}\left(\frac{\epsilon ^2}{8\sigma ^2}\right)\hfill & \chi =2\hfill \end{array}$$
(14)
The asymptotic behaviour has a functional dependence on $`\epsilon `$ similar to that of the Wigner surmise, i.e. the NNSD (11) goes linearly to zero for $`\epsilon 0`$ and it falls down according to $`\mathrm{exp}(\epsilon ^2)`$ for $`\epsilon \mathrm{}`$. However, as seen from (13)–(14), the natural dependence on $`\chi `$ provides a certain scaling.

To calculate various statistical characteristics, it is often required to know certain moments of the distribution. Thus, we have derived analytical expressions for some moments of the distribution (12), and the lowest are given in table 1.

From (12) and the expression of the mean spacing, $`D`$ (cf table 1), we have derived the distribution of the dimensionless nearest-neighbour spacing, $`s=\epsilon /D`$ :

$$q_\chi (s)=D\stackrel{~}{p}_\chi (sD).$$
(15)
For $`\chi =2`$, one finds the Wigner distribution (2). It can be shown that for $`\chi `$ and $`\chi ^{}=4/\chi `$, the two functions $`q_\chi `$ and $`q_\chi ^{}`$ are equal.

The distributions $`q_\chi `$ for $`\chi =4`$, $`\chi =30`$, $`\chi =500`$ and $`\chi =1000`$ are plotted in figure 2(a-d) and compared to the Poisson (dashed), Wigner (dotted) and semi-Poisson (dashed-dotted) distributions. Note that the semi-Poisson distribution can be fairly approximated for $`\chi =30`$ although its asymptotic behaviour for large $`s`$ is different. As follows from figure 2, the parameter $`\chi `$ allows to describe a transition of a quantum system from a completely chaotic limit ($`\chi =2`$) to a nearly integrable one ($`\chi \mathrm{}`$ or $`\chi 0`$). However, the integrability characterised by the Poisson distribution (1) can never be reached (compare with the model of Caurier et al ).

Indeed, integrating $`p_\chi `$ in (11) over $`E_\beta `$ from $`\mathrm{}`$ to $`E_\alpha `$ and over $`E_\alpha `$ from $`\mathrm{}`$ to $`+\mathrm{}`$, we obtain the angular distribution

$$r_\chi \left(\theta \right)=\frac{1}{\pi }\sqrt{\frac{\chi }{2}}\frac{1}{\left[1+\frac{1}{2}\left(\chi 2\right)\mathrm{sin}^2\left(2\theta \right)\right]}.$$
(16)
This distribution is represented in figure 3 for different values of $`\chi `$. For $`\chi =2`$ it is an exactly uniform distribution, which means that there is no privileged basis (the orthogonal invariance holds). For high values of $`\chi `$, the initial basis is nearly the eigenbasis (the diagonal elements are much larger than the non-diagonal ones), thus $`r_\chi `$ takes its maximum absolute values for $`\theta =0`$ and $`\theta =\pi /2`$, whereas for small values of $`\chi `$, the eigenstates are more likely obtained after a rotation of $`\pi /4`$ of the initial basis and $`r_\chi `$ is maximum for $`\theta =\pi /4`$. For $`\chi `$ and $`\chi ^{}=4/\chi `$ the two curves are in quadrature.

We can re-express $`p_\chi (E_\alpha ,E_\beta )`$ as a function of $`\epsilon `$ and $`S=E_\alpha +E_\beta `$. Then $`p_\chi `$ can be factorised into a function depending on $`\epsilon `$ times a function depending on $`S`$, i.e. these variables are independent. Moreover since $`S`$ is the trace of the matrix it is a Gaussian variable with zero mean and all its odd moments are zero. From the independence of $`\epsilon `$ and $`S`$ we deduce that the moments of the eigenvalues fulfill

$`E_\alpha ^n`$ $`=`$ $`(1)^nE_\beta ^n`$ (17)
$`E_\alpha ^n`$ $`=`$ $`{\displaystyle \frac{1}{2^n}}{\displaystyle \underset{p=0}{\overset{n/2}{}}}({}_{2p}{}^{n})S^{2p}\epsilon ^{n2p}`$ (18)
From (17) and (18), we can derive the moments of the eigenvalues whose distributions are difficult to compute. Note that we deduce from (17) that the highest and the lowest eigenvalues have opposite mean values and the same variance.

## 3 Conclusion

The study of the statistical properties of spectra of realistic Hamiltonians requires the consideration of a random matrix ensemble whose elements are not independent or whose distribution is not invariant under orthogonal transformation of a chosen basis. In this letter we have concentrated on the properties of (2$`\times `$2) real symmetric matrices whose elements are independent Gaussian variables with zero means but do not belong to the GOE. We have derived the distribution of eigenvalues for such a matrix, the NNSD which generalises the Wigner surmise and we have calculated some important moments. The asymptotic properties of the distribution obtained are functionally identical to those of the ordinary Wigner surmise. For finite $`\chi `$, the model considered here allows to describe the transition from chaos to near integrability (the exact integrable limit is never realized). Thus it represents a chaotic system although with a degree of disorder less important than in the Wigner surmise. The derivation of similar analytical expressions for matrices of larger dimensions is technically difficult. However, we believe that the present results already justify the use of NNND of type (12) to fit the data as an alternative to the Brody distribution. We also think that these results can be extended to hermitian matrices.

## Acknowledgements

We would like to thank Olivier Juillet for his help to calculate some integrals and Alejandro Frank for his helpful and enthusiastic discussions. PCH-T would like to thank Sylvie Hudan for her help to draw the figures. NAS is a post-doctoral researcher of the FWO-Vlaanderen, Belgium.

## References