text:
# The ΛΛ→𝑌𝑁 Weak Transition in Double-Λ Hypernuclei

## 1 Introduction

Recent analyses of hypernuclear nonmesonic decay reveal the importance of short-range interaction of baryonic weak transition ITO:NPA ; Par:PRC ; IOMI:NPA ; SIO01:NPA ; Jid:NPA ; Ito:PRC ; SIO02:NPA . It is, in fact, realized that the momentum transfer is large due to the $`\mathrm{\Lambda }N`$ mass difference, and therefore the short range interaction becomes important. We proposed to treat the short range part of $`\mathrm{\Lambda }NNN`$ interaction using the valence quark picture of the baryon and the effective four-quark weak Hamiltonian ITO:NPA ; IOMI:NPA ; SIO01:NPA ; SIO02:NPA . We found that the direct quark (DQ) process gives significantly large contribution and shows qualitatively different features from the meson exchange contribution, especially in its isospin structure. The kaon exchange contribution also largely enhance the $`\mathrm{\Gamma }_{nn}/\mathrm{\Gamma }_{pn}`$ ratio due to its spin structure Par:PRC ; SIO01:NPA ; Jid:NPA ; SIO02:NPA .

There remains, however, another problem, i.e., experimental values of the proton asymmetry from polarized hypernucleus are incompatible with the theoretical values, which is closely related to the $`J=0`$ transition amplitudes Nab:PRC . The $`J=0`$ amplitude is relatively small and therefore masked by the $`J=1`$ amplitude, but it is important to clarify the breaking of the $`\mathrm{\Delta }I=1/2`$ rule for the nonmesonic weak decay of hypernuclei. In fact, the quark model treatments of $`\mathrm{\Lambda }NNN`$ transition predict large breaking of the $`\mathrm{\Delta }I=1/2`$ rule in $`J=0`$ transition channel ITO:NPA ; IOMI:NPA ; SIO02:NPA ; Mal:PhL .

For the case of double-$`\mathrm{\Lambda }`$ hypernuclei, the nonmesonic decay is induced not only by the nucleon but also by the other $`\mathrm{\Lambda }`$ hyperon: $`\mathrm{\Lambda }NNN\mathrm{and}\mathrm{\Lambda }\mathrm{\Lambda }YN`$. The two-body part of the nonmesonic decay rate of double-$`\mathrm{\Lambda }`$ hypernuclei is classified as

$`\mathrm{\Gamma }_{nm}`$ $`=`$ $`\mathrm{\Gamma }_{NN}+\mathrm{\Gamma }_{YN}`$ (1)
$`=`$ $`\mathrm{\Gamma }_{nn}+\mathrm{\Gamma }_{pn}(+\mathrm{\Gamma }_{pp})+\mathrm{\Gamma }_{\mathrm{\Lambda }n}+\mathrm{\Gamma }_{\mathrm{\Sigma }^{}p}+\mathrm{\Gamma }_{\mathrm{\Sigma }^0n}`$
where $`\mathrm{\Gamma }_{pp}`$ stands for the $`\mathrm{\Sigma }^+ppp`$ transition which is available when the virtual $`\mathrm{\Sigma }`$ exists SIO02:NPA . The $`\mathrm{\Gamma }_{YN}`$ represents the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition, which is a characteristic transition for the multi-$`\mathrm{\Lambda }`$ system. Experimentally, the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition should be distinguishable from the $`\mathrm{\Lambda }NNN`$ decay mode because the hyperon with high momentum is emitted once it occurs in double-$`\mathrm{\Lambda }`$ hypernucleus. Thus the hyperon-induced decay mode provides us with valuable information on the $`|\mathrm{\Delta }S|=1`$ baryon-baryon interaction IUM:NPA ; PRB:PRC .

In sec.2, the hyperon-induced decay rate, $`\mathrm{\Gamma }_{YN}`$, is expressed in terms of two-body amplitudes corresponding to the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition. In sec.3, we construct the meson exchange and DQ potential for the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition. In sec.4, we show the result of the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition rate in $`{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{}{}_{}{}^{6}`$He. Finally, sec.5 is devoted for the conclusion of this work.

## 2 The $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition rate

For the nonmesonic weak decay of double-$`\mathrm{\Lambda }`$ hypernuclei, the hyperon induced decay mode, $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$, becomes possible besides the usual $`\mathrm{\Lambda }NNN`$ transition. Now we concentrate on the novel $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ decay. The $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition rate is given by

$$\mathrm{\Gamma }_{YN}=\frac{d^3p_1^{}}{(2\pi )^3}\frac{d^3p_2^{}}{(2\pi )^3}\frac{1}{2J_H+1}\underset{all}{}(2\pi )\delta (E_fE_i)|M_{fi}|^2$$
(2)
where $`M_{fi}`$ is the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition amplitude, $`J_H`$ is the total spin of initial hypernucleus, and $`p_{}^{}{}_{1}{}^{}`$ and $`p_{}^{}{}_{2}{}^{}`$ are momenta of emitted particles, i.e., hyperon and nucleon. The summation indicates a sum over all quantum numbers of the initial and final particle systems. The $`\delta `$-function for the energy conservation is given by

$$\delta (E_fE_i)=\delta \left(M_Y+M_N+\frac{p_{1}^{}{}_{}{}^{2}}{2M_Y}+\frac{p_{2}^{}{}_{}{}^{2}}{2M_N}2M_\mathrm{\Lambda }ϵ_{\mathrm{\Lambda }\mathrm{\Lambda }}\right)$$
(3)
where the hyperon binding energy is taken as $`ϵ_{\mathrm{\Lambda }\mathrm{\Lambda }}7.25`$MeV NAG:PRL . After the decomposition of angular momentum, the explicit form of $`M_{fi}`$ is

$$|M_{fi}|^2=(4\pi )^4\left|\mathrm{\Psi }_f^{L^{}S^{}J}(R,r^{})V_{SS^{}J}^{LL^{}}(r,r^{})\mathrm{\Psi }_i^{LSJ}(R,r)r^2𝑑rr^2𝑑r^{}R^2𝑑R\right|^2$$
(4)
where $`V_{SS^{}J}^{LL^{}}(r,r^{})`$ is the transition potential and $`\mathrm{\Psi }^{LSJ}(R,r)`$ is the wave function of the $`\mathrm{\Lambda }\mathrm{\Lambda }`$ or YN two-body system in the configuration space. The $`L`$, $`S`$, and $`J`$ indicate the angular momentum, spin, and total spin for two-body system, respectively.

In order to extract the $`\mathrm{\Lambda }\mathrm{\Lambda }`$ state from the double-$`\mathrm{\Lambda }`$ hypernucleus, we assume, for simplicity, that the $`\mathrm{\Lambda }\mathrm{\Lambda }`$ state with $`J=0`$ and $`I=0`$ is decoupled from the residual nuclei:

$$|{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{A}Z=|(\mathrm{\Lambda }\mathrm{\Lambda })_{J=0}^{I=0}|{}_{}{}^{A2}Z$$
(5)
where the double-$`\mathrm{\Lambda }`$ state is always in $`{}_{}{}^{1}S_{0}^{}`$ state in the ground state. The important point of the study the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ process is that the initial $`\mathrm{\Lambda }\mathrm{\Lambda }`$ pair is in the $`{}_{}{}^{1}S_{0}^{}`$ state and thus the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition occurs only in the $`J=0`$ channel. The $`J=0`$ transition amplitude is not well studied in the $`\mathrm{\Lambda }NNN`$ decay, because the $`J=1`$ amplitude is much larger. In contrast, the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ decay directly probes the $`J=0`$ transition amplitude and therefore plays a complementary role.

The initial double-$`\mathrm{\Lambda }`$ state is written in terms of the shell model wave function of $`0s`$ state :

$$\mathrm{\Psi }_{\mathrm{\Lambda }\mathrm{\Lambda }}^{000}(\stackrel{}{r}_{\mathrm{\Lambda }_1},\stackrel{}{r}_{\mathrm{\Lambda }_2})=\left[\varphi _{0s}^b(\stackrel{}{r}_{\mathrm{\Lambda }_1})\times \varphi _{0s}^b(\stackrel{}{r}_{\mathrm{\Lambda }_2})\right]^{J=0}=\left[\varphi _{0s}^{b_r}(\stackrel{}{r})\times \varphi _{0s}^{b_R}(\stackrel{}{R})\right]^{J=0}$$
(6)
where $`\varphi _{0s}^b`$ is assumed to be a Gaussian with size parameter $`b`$. The last expression of eq.(6) is the wave function in terms of the relative and center-of-mass coordinate with $`\stackrel{}{r}=\stackrel{}{r_{\mathrm{\Lambda }_1}}\stackrel{}{r_{\mathrm{\Lambda }_2}}`$ and $`\stackrel{}{R}=(\stackrel{}{r_{\mathrm{\Lambda }_1}}+\stackrel{}{r_{\mathrm{\Lambda }_2}})/2`$. The $`b_r`$ and $`b_R`$ are $`\sqrt{2}b`$ and $`b/\sqrt{2}`$, respectively. The wave function for the final $`YN`$ state is assumed as the plane wave with $`\stackrel{}{K}=\stackrel{}{r}_Y+\stackrel{}{r}_N`$ and $`\stackrel{}{k}=(m_N\stackrel{}{p}_{}^{}{}_{1}{}^{}m_Y\stackrel{}{p}_{}^{}{}_{2}{}^{})/(m_Y+m_N)`$. The $`b`$ parameter is taken as $`1.9`$fm which fit on the r.m.s. radius of $`\mathrm{\Lambda }`$ in $`{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{}{}_{}{}^{6}`$He AN:Pcom .

In addition, we introduce phenomenological correlation functions. For the initial $`\mathrm{\Lambda }\mathrm{\Lambda }`$ wave function, we multiply to the uncorrelated wave function

$$f_{\mathrm{\Lambda }\mathrm{\Lambda }}(r)=\left(1e^{\frac{r^2}{a^2}}\right)^n+br^2e^{\frac{r^2}{c^2}}$$
(7)
with $`a=0.8`$fm, $`b=0.12`$fm-2, $`c=1.28`$fm, $`n=1`$, according to ref.PRB:PRC . On the other hand, for the final $`YN`$ state, we multiply

$$f_{YN}(r)=1j_0(q_cr)$$
(8)
where $`j_0(x)`$ is a spherical Bessel function with $`q_c=3.93`$ fm-1.

## 3 Transition potential

In order to describe the weak $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition, we adopt the quark-meson hybrid model, which is focused on the quark degrees of freedom in a short-range region. At long and medium distances, the transition is induced mainly by the one pion exchange (OPE) and one kaon exchange (OKE) processes. The general form of $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition potential in the meson exchange picture is given by

$$V_{\mathrm{\Lambda }\mathrm{\Lambda }YN}(\stackrel{}{q})=g_s[\overline{u}_Y\gamma _5u_\mathrm{\Lambda }]\frac{1}{\stackrel{}{q}^2+\stackrel{~}{m}_i^2}\left(\frac{\mathrm{\Lambda }_i^2\stackrel{~}{m}_i^2}{\mathrm{\Lambda }_i^2+\stackrel{}{q}^2}\right)^2g_w[\overline{u}_N(A+B\gamma _5)u_\mathrm{\Lambda }]$$
(9)
where the coupling constants $`g_s`$, $`g_w`$, $`A`$ and $`B`$, shown in Table 1, are chosen properly for each transition. A monopole form factor with cutoff parameters, $`\mathrm{\Lambda }_\pi `$ = 800 MeV and $`\mathrm{\Lambda }_K=1300`$ MeV, are employed for each vertex. As the energy transfer is significantly large, we introduce the effective meson mass:

$$\stackrel{~}{m}=\sqrt{m^2(q^0)^2},q^0=\{\begin{array}{cc}\hfill 88.5\mathrm{MeV}& \mathrm{for}\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n\mathrm{transition}\hfill \\ \hfill 50.2\mathrm{MeV}& \mathrm{for}\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N\mathrm{transition}\hfill \end{array}.$$
(10)
The kaon exchange potential can be constructed similarly. Both the strong and weak coupling constants are evaluated under the assumption of the flavor SU(3) symmetry. Note that, for OKE, it involves the strangeness transfer and thus the strong and weak vertices should be exchanged.

For shorter distances, we employ the direct quark (DQ) transition potential based on the constituent quark picture of baryons. In this mechanism, a strangeness changing weak interaction between two constituent quarks induces the transition of two baryons.

The DQ potential is given as a nonlocal form as

$$V_{DQ}^{}{}_{SS^{}J}{}^{LL^{}}(r,r^{})=\frac{G_F}{\sqrt{2}}\times W\underset{i=1}{\overset{7}{}}\left\{V_i^ff_i(r,r^{})+V_i^gg_i(r,r^{})+V_i^hh_i(r,r^{})\right\}$$
(11)
where $`r`$ ($`r^{}`$) stands for the radial part of the relative coordinate in the initial (final) state. The explicit forms of $`f_i`$, $`g_i`$, and $`h_i`$ are given in ref. SIO01:NPA , and the coefficients, $`V_i^k`$, for the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transitions are listed in Table 2.

## 4 Nonmesonic Decay of Double-$`\mathrm{\Lambda }`$ Hypernuclei

Now we show the results of $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ decay rates of $`{}_{\mathrm{\Lambda }\mathrm{\Lambda }}{}^{}{}_{}{}^{6}`$He. In Table 3, the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition is estimated by extending the $`\pi +K+`$DQ model to the $`S=2`$ system.

For the case of $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n`$ transition, the OPE process is forbidden because of the isospin symmetry on the strong vertex. Although the short-range interaction becomes dominant, their contributions are drastically reduced by the short range correlation. The DQ contribution is much larger than the kaon exchange contribution for the PV transition channel, but it is about $`0.1\%`$ of the free $`\mathrm{\Lambda }`$ decay rate. In our full calculation, we find that the DQ reduces the total PC transition rate to almost zero. In contrast, the meson exchange and the DQ are added coherently in the PV channels, resulting in a relatively large PV transition rate, which is about $`0.2\%`$ of the free $`\mathrm{\Lambda }`$ decay.

The OPE contributions to the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transition is expected to be small in the $`J=0`$ transition. In fact, some previous analyses have revealed that the $`J=0`$ decay rate is only few $`\%`$ of the $`J=1`$ transition rate. Table 3 shows that the decay rate in OPE is about $`1\%`$ of the free $`\mathrm{\Lambda }`$ decay in $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }^{}p`$ transition. The kaon exchange process additively contributes to the OPE process, and the decay rate becomes about twice as large as the OPE alone in both the PC and PV channels. Thus the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ decay rate in meson exchange processes becomes much larger than the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n`$. The DQ transition rate for the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transition is also larger than the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n`$ case, and it behaves differently from the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n`$ transition. In our full calculation, the PV transition rate is largely reduced, while the PC is enhanced due to the DQ contribution. Therefore the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transition occurs mainly in the PC channel.

It is also interesting to check the breaking of the empirical $`\mathrm{\Delta }I=1/2`$ rule for the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition, which leads to

$$\mathrm{\Gamma }_{\mathrm{\Sigma }^{}p}=2\mathrm{\Gamma }_{\mathrm{\Sigma }^0n}.$$
(12)
This relation shows that the nonmesonic weak decay of double-$`\mathrm{\Lambda }`$ hypernuclei provides us with an opportunity of direct confirmation of the validity of the empirical $`\mathrm{\Delta }I=1/2`$ dominance. In our prediction, this relation is badly violated due to the DQ contribution, in which the $`\mathrm{\Delta }I=3/2`$ amplitude exists through the weak four-quark Hamiltonian. This situation is similar to the $`\mathrm{\Lambda }NNN`$ transition. Thus the test of the $`\mathrm{\Delta }I=1/2`$ rule in the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition may reveal the importance of the short-range part in the weak baryonic interaction, and may also help to understand the mechanism of the $`\mathrm{\Delta }I=1/2`$ enhancement observed in the nonleptonic weak decays of the free hyperons and kaons.

We compare our results with the previous calculations carried out by Itonaga et al. IUM:NPA  and Parreño et al. PRB:PRC . We all agree the bare pion-exchange amplitudes in the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transitions IP:Pdis . In Table 3, Ref. IUM:NPA  includes effects of two-pion exchange in the scalar ($`\sigma `$) channel as well as the vector ($`\rho `$) channel. Ref. PRB:PRC  contains all the pseudoscalar mesons and the vector mesons exchanged between $`\mathrm{\Lambda }\mathrm{\Lambda }`$. Ref.PRB:PRC  took two different approaches for the final state interaction. The full calculation (F) is based on the scattering $`T`$ matrix in the final states, while an effective treatment (E) evaluates the transition amplitudes by multiplying a correlation function to the relative wave function. They pointed out that the full calculation enhances $`\mathrm{\Gamma }_{\mathrm{\Lambda }n}`$ while it reduces $`\mathrm{\Gamma }_{\mathrm{\Sigma }N}`$. The effective treatment does not include the conversion from $`\mathrm{\Sigma }N`$ to $`\mathrm{\Lambda }N`$.

In the present work, we make an effective treatment of the final state interactions. Therefore, our results have to be compared to the ones in Ref.PRB:PRC (E). Our calculation is similar to their effective treatment and thus the $`\mathrm{\Gamma }_{\mathrm{\Lambda }n}`$ enhancement due to the final state interaction is not shown. The total decay rates in our results are reduced by a factor about two due to the larger $`b`$ parameter for the $`\mathrm{\Lambda }`$ wave function. We note that if the same $`b`$ value is used, our results agree with the (E) results in Ref.PRB:PRC . We also point out that the ratio $`\mathrm{\Gamma }_{\mathrm{\Sigma }^{}p}/\mathrm{\Gamma }_{\mathrm{\Sigma }^0n}`$ is fixed to $`2`$ by the $`\mathrm{\Delta }I=1/2`$ rule in the meson exchange calculations.

## 5 Conclusion

In the double $`\mathrm{\Lambda }`$ hypernuclear system, the novel decay mode, $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$, is allowed. The hyperon-induced decay rates are much smaller than the decay rate of free $`\mathrm{\Lambda }`$ ($`\mathrm{\Gamma }_\mathrm{\Lambda }`$). This is because the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition occurs only in the $`J=0`$ channel. Although the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition is not easy to detect experimentally, it is important to extract the information about the $`J=0`$ transition.

Our result shows that the DQ process gives significant contribution to the $`\mathrm{\Lambda }\mathrm{\Lambda }YN`$ transition. It strongly suppresses the PC amplitude for the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Lambda }n`$ transition and the PV amplitude for the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transition. Furthermore the DQ largely breaks the empirical $`\mathrm{\Delta }I=1/2`$ rule. Unlike the $`\mathrm{\Lambda }NNN`$ transition, we have an opportunity to check directly the breaking of $`\mathrm{\Delta }I=1/2`$ dominance for the $`\mathrm{\Lambda }\mathrm{\Lambda }\mathrm{\Sigma }N`$ transition.

## Acknowledgment

The authors acknowledge Profs. A. Parreño, K. Itonaga, H. Nemura, and Y. Akaishi for variable discussions. One of authors, K. S., thanks to JSPS Research Fellowship for financial support. This work is supported in part by Grants-in-Aid for Scientific Research (No.11640261) from Japan Society for the Promotion of Science (JSPS).


format: text

text:
# The puzzling dynamical status of the core of the globular cluster NGC 6752 1footnote 11footnote 1Based on observations with the NASA/ESA HST, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. Also based on WFI observations collected at the European Southern Observatory, La Silla, Chile, within the observing programmes 62.L-0354 and 64.L-0439.

## 1 Introduction

Although NGC 6752 is one of the nearest clusters, there is no consensus in the literature on its dynamical status. It was initially classified as a post-core collapse (PCC) by Djorgovski & King (1986, hereafter DK86) and Auriere & Ortolani (1989, hereafter AO89). Later Lugger, Cohn & Grindlay (1995, hereafter LCG95) argued that the radial profile was not inconsistent with a King model. As part of a project devoted to the study the characteristics of globulars harboring millisecond pulsars (MSPs), we present a new determination of the center of gravity (§3) and of the star density profile (§4) of NGC 6752 both based on new high-resolution and wide-field observations. By coupling the superior resolution of HST with the imaging capability of the wide field camera (WFI) mounted at the 2.2m ESO-MPI telescope, we have obtained the most accurate and extended radial profile ever published for this cluster. After modeling these data (§5) we examine the dynamical status of the cluster (§6).

Recently, observations of $`|\dot{P}/P|`$ and location of five millisecond pulsars in NGC 6752 (D’Amico et al. 2002) have suggested a surprisingly high mass-to-light ratio $`/_V`$ in its core and the occurrence of non thermal dynamics in the inner regions (see Colpi, Possenti & Gualandris 2002). The determination of $`/_V`$ by D’Amico et al. relied on results of pulsar timing and published optical data derived from medium resolution ground based observation only. Taking advantage from our new optical data, we re-examine in §7 the D’Amico et al. measurement, discussing the possibile origin and the consequences of the observed values of $`|\dot{P}/P|`$ at some length.

## 2 Observations and data analysis

The photometric data used here consist of two sets: (i)—High resolution set—a series of HST- WFPC2 images was obtained in March 2001, through the F555W ($`V`$) and F336W ($`U`$) filters as part of a long term project (GO-8709, PI: F. R. Ferraro) aimed to study the central stellar populations in a set of Galactic Globular Glusters (GGCs). In this dataset the planetary camera (PC, which has the highest resolution $`0\stackrel{}{\mathrm{.}}046/\mathrm{pixel}`$) was roughly centered on the cluster center while the Wide Field (WF) cameras (at lower resolution $`0\stackrel{}{\mathrm{.}}1/\mathrm{pixel}`$) sampled the surrounding outer regions; (ii)—Wide Field set—a complementary set of multi-filter ($`B`$, $`V`$, $`I`$) wide field images was secured during an observing run at the 2.2m ESO-MPI telescope at ESO (La Silla) in July 1999, using the Wide Field Imager (WFI). The WFI has exceptional imaging capabilities—the image consists of a mosaic of 8 CCD chips (each with a field of view of $`8^{}\times 16^{}`$) giving a global field of view of $`33^{}\times 34^{}`$. The cluster was roughly centered on chip $`\mathrm{\#}2`$.

### 2.1 Photometry

The raw WFI images were corrected for bias and flat field, and the over-scan region was trimmed using standard IRAF222IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Fundation. tools. The point spread function (PSF) fitting procedure was performed independently on each $`V`$ and $`B`$ images, using DAOPHOT II (Stetson 1994). A final catalog listing the instrumental $`B,V`$ magnitudes for all the stars in each field has been obtained by cross-correlating the $`B`$ and $`V`$ catalogs. The WFI catalog was finally calibrated by using the data-set of Buonanno et al. (1986).

The photometric reductions of the high resolution images were carried out using ROMAFOT (Buonanno et al. 1983), a package developed to perform accurate photometry in crowded fields and specifically optimized to handle under-sampled PSF (Buonanno & Iannicola 1989) as in the case of the HST-WF chips.

PSF-fitting instrumental magnitudes have been obtained using the standard procedure described in Ferraro et al. (1997a, 2001). The final catalog of the F555W and F336W magnitudes was calibrated by using the zero-points listed by Holtzman et al. (1995).

### 2.2 Astrometry

The recently released Guide Star Catalog ($`GSCII`$) was used to search for astrometric standards in the entire WFI image field of view. Several hundred astrometric $`GSCII`$ reference stars were found in each chip, allowing an accurate absolute positioning of the sources. An astrometric solution has been obtained for each of the 8 WFI chips independently, by using a procedure developed at the Bologna Observatory (Montegriffo et al. 2002, in preparation). At the end of the entire procedure, the rms residuals were of the order of $`0\stackrel{}{\mathrm{.}}2`$ both in RA and Dec.

The small field ($`2\stackrel{}{\mathrm{.}}5`$ on the side) of the high resolution WFPC2/HST images was entirely contained within the field of view of the WFI chip $`\mathrm{\#}2`$, which imaged the central part of the cluster. Thus we used more than 1200 bright stars in the WFI catalog lying in the WFPC2-FoV as secondary astrometric standards in order to properly find an astrometric solution for the WFPC2 catalog. We estimate that the global uncertainty in the astrometric procedure is less than $`0\stackrel{}{\mathrm{.}}4`$ both in RA and Dec. At the end of the procedure the two catalogs (WFPC2 and WFI) have a fully homogeneous absolute coordinate system.

The two lists were then matched together: stars in the overlapping area were used to recheck the homogeneity of the calibration of the $`V`$-magnitude (which is used to construct the surface profile) in the two catalogs. To do this the HST F555W band was transformed into the Johnson $`V`$ system and a homogeneous list of $`V`$ magnitudes and absolute coordinates for the sources in the HST and the WFI catalogs was produced. In order to avoid strong selection effects due to possible incompleteness of the samples (in particular the WFI catalog), we grouped the two samples as follows: (1) HST sample: all stars in the WFPC2 FoV with $`V<18.5`$ and $`r<96\mathrm{}`$ from the cluster center; (2) WFI sample: all the stars in the WFI FoV with $`V<18.5`$ and $`r>120\mathrm{}`$ from the cluster center. The Color-Magnitude Diagrams (CMDs) derived from these samples are shown in Figure 1.

## 3 The center of gravity

The first step toward the computation of the density profile is the determination of the center of gravity $`C_{\mathrm{grav}}`$ of the cluster. In doing this, we estimated the position of the geometrical center of the star distribution taking advantage from the knowledge of the exact star positions even in the innermost central region. We applied the procedure described in Montegriffo et al. (1995) which computed $`C_{\mathrm{grav}}`$ by simply averaging the $`\alpha `$ and $`\delta `$ coordinates of stars lying in the PC camera of the HST catalog. $`C_{\mathrm{grav}}`$ is located at $`\alpha _{\mathrm{J2000}}=19^\mathrm{h}\mathrm{\hspace{0.17em}10}^\mathrm{m}\mathrm{\hspace{0.17em}52}\stackrel{\mathrm{s}}{\mathrm{.}}04,\delta _{J2000}=59\mathrm{°}\mathrm{\hspace{0.17em}59}\mathrm{}\mathrm{\hspace{0.17em}04}\stackrel{}{\mathrm{.}}64`$ with a typical $`1\sigma `$ uncertainty of $`0\stackrel{}{\mathrm{.}}5`$ in both $`\alpha _{\mathrm{J2000}}`$ and $`\delta _{J2000}`$, corresponding to about $`10`$ pixels in the WFPC2/HST images.

AO89 derived the $`C_{\mathrm{grav}}`$ by using the barycenter of the bright ($`V<16`$) resolved stars from ground based observations. We transformed the center position in their Figure 3,333Note that the orientation in both the maps shown by AO89—Figures 1 and 3, respectively—is incorrectly reported in their captions since E is on the right and N is at the bottom. into our coordinate system. It well agrees with our determination, being located only $`2\mathrm{}`$ S.

The coordinates of the cluster center previously reported in the literature refer mostly to the center of luminosity $`C_{\mathrm{lum}}`$, which is usually determined by the so-called mirror autocorrelation (see Djorgovski 1988). $`C_{\mathrm{grav}}`$ turns out to be $`10\mathrm{}`$ S and $`2\mathrm{}`$ E of the $`C_{\mathrm{lum}}`$ reported in the Djorgovski (1993) compilation. Such a difference is not surprising—similar offsets have been found in other clusters (see for example 47 Tucanae—Calzetti et al. 1993, Montegriffo et al. 1995; M80—Ferraro et al. 1999a).

Figure 2 shows the central $`20\mathrm{}\times 20\mathrm{}`$ of the cluster with respect to the $`C_{\mathrm{grav}}`$ determined in this work (marked in Figure with a large cross at $`(0,0)`$). The center position by AO89 (small cross labeled with AO) and by Djorgovski (1993; small cross labeled with D) are also indicated.

Pooley et al. (2002) noticed that the highest concentration of X-ray sources detected by Chandra in a 30 ksec pointing toward NGC 6752 appeared surprisingly displaced with respect to the optical center by Djorgovski (1988). On the contrary, Figure 2 shows that our value of $`C_{\mathrm{grav}}`$ lies close to the barycenter of the nine more central objects (marked in Figure 2 with a large cross labeled with XRS at coordinates ($`\alpha _{\mathrm{J2000}}=19^\mathrm{h}\mathrm{\hspace{0.17em}10}^\mathrm{m}\mathrm{\hspace{0.17em}51}\stackrel{\mathrm{s}}{\mathrm{.}}82,\delta _{J2000}=59\mathrm{°}\mathrm{\hspace{0.17em}59}\mathrm{}\mathrm{\hspace{0.17em}03}\stackrel{}{\mathrm{.}}84`$), reconciling their projected positions with an almost spherical distribution.

## 4 The surface brightness and star density profiles

By using the combined data set shown in Figure 1 we computed star density and surface brightness profiles applying the standard procedure fully described in Calzetti et al. 1993 (see also Ferraro et al. 1999a, Paltrinieri et al. 2001). The entire photometric sample has been divided in 35 concentric annuli centered on $`C_{\mathrm{grav}}`$, spanning a spatial range from $`0\mathrm{}`$ to $`27\mathrm{}`$; each annulus has been in turn split in a number of sub-sectors (generally octants or quadrants depending on the shape of the FoV covered by the HST and WFI fields). The surface brightness of each sub-sector has been evaluated as the sum of the luminosity of all the stars lying in it, divided by the area (expressed in $`\mathrm{arcsec}^2`$). The average and the rms over the brightness of the sub-sectors determine the value and the uncertainty of the surface brightness of an annulus. However, evaluating the brightness over small regions can suffer by large fluctuations due to small number statistics of few bright giants. Hence, we computed three radial profiles, removing the stars brighter than $`V=12,V=13,V=14`$, respectively. Figure 3 shows that the overall structure of the profile does not change with the adopted magnitude limit. In contrast, the rms between sub-sectors significantly decreases once the brightest stars are excluded. The brightness profile by AO89 is shown (as open circles) for comparison in the lower panel of Figure 3. The agreement with the results of this work is excellent, although their profile is slightly fainter in the very central region, probably due to the incompleteness of their sample. This might be due to the fact that removal of the giants can lead to an unintentional subtraction of the contribution of close, fainter stars, which are often unresolved in ground-based observations.

We can take advantage from the fact that we resolve the stars even in the innermost region of the cluster and directly produce a star density profile. As noted by LCG95, a star density profile is not affected by the statistical fluctuations caused by the bright stars, and it is the most robust mean for determining the cluster structure parameters. Once the photometric sample has been divided in annuli and sub-sectors, the number of stars lying within each sub-sector has been counted and averaged. Star density has been eventually obtained by dividing the average star number by the area of the sub-sector (expressed in $`\mathrm{arcsec}^2`$). The star density profile values obtained from this procedure are listed in Table 1. This is the most complete and extended density profile ever published for this cluster, since it samples the cluster population from the very inner region out to $`r27\mathrm{}`$ from the cluster center.

## 5 Modeling the radial density profile

Usually GGCs are considered core-collapsed or not depending on how well their radial distribution of stars is fitted by King (1966) profiles (Trager et al. 1995). These models are characterized by two parameters, the core radius, $`r_c`$, and the tidal radius, $`r_t`$, or, alternatively, the concentration $`c=\mathrm{log}(r_t/r_c)`$. As stated by Meylan & Heggie (1997, hereafter MH97) a more general rule is that all clusters with a concentration parameter $`c>2`$–2.5 can be considered as collapsed, on the verge of collapsing or just beyond the collapse.

In order to reproduce the observed profile, we obtained the surface density by projecting the star density from a standard isotropic, single-mass King-model derived using the code described in Sigurdsson & Phinney (1995).444In our earlier work we have used an analytic representation of King profiles. Parameter values derived from the analytic fits differ somewhat from those derived using the current scheme. The result is shown in Figure 4. As can be seen, a King model which properly fits the outer region significantly overestimates the star density over an extended inner region ($`5`$$`20\mathrm{}`$). This fit closely resembles the result obtained by LCG95 who fitted the $`U`$-band surface brightness profile for $`r<120\mathrm{}`$ with a $`c=2`$ King model (see their Figure 2r). On the basis of that fit, they concluded that the cluster is not required to be in a post-core collapse state because the model did not differ from the data in a statistically significant way. However, there were aspects of the fit which led us to explore other possibilities. These included the bad fit in the region $`5<r<20\mathrm{}`$ and the fact that even in the outer parts there are regions where the data line systematically above or below the observed points. We suspected that unmodeled effects were degrading this fit.

For this reason, we searched for alternative solutions. As shown in Figure 5 the density profile can be well fit by a dual King model. The outer cluster ($`r>10\mathrm{}`$) is well fit by a model with $`r_c=28\mathrm{}`$ and $`c=1.9`$. The observed counts are significantly in excess of this model within the central $`8\mathrm{}`$. In this innermost region a King model with $`c=2.1`$ and $`r_c=5.7\mathrm{}`$ fits the data, but lies significantly below all of the points with $`r>10\mathrm{}`$. A dual King model does not represent a detailed equilibrium of cluster structure. However, it is consistent with the scenario in which the central regions have evolved away from a global King model that now only characterizes the outer regions.

This anomalous structure of the innermost radial profile of NGC 6752 has been evident since the very first studies (Da Costa 1979). The $`U`$-band data presented by DK86 showed a sharp shoulder in the radial profile and prompted the authors to classify NGC 6752 as a core-collapsed cluster. A few years later, AO89 modeled their detailed $`V`$-band profile with a power law in the region between $`3\mathrm{}`$ and $`60\mathrm{}`$ from the center, noting a significant flattening for $`r<3\mathrm{}`$.

LCG95 used also a modified power-law (see their eq. 3, also reported in Figure 6), in order to test the presence of resolved cores. We used the inner part of the star density profile of Figure 5 to repeat an analysis similar to that presented by LCG95. Besides a normalization factor, their expression is a function of two free parameters: the power-law index $`\alpha `$ and the scale-length of the core $`r_0`$, which can be easily converted in the usual $`r_c`$ (see their eq. 2). As can be seen in Figure 6 a proper fit (with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=0.28`$ for $`8`$ d.o.f) can be obtained, yielding $`r_0=3\stackrel{}{\mathrm{.}}1\pm 1\stackrel{}{\mathrm{.}}4`$ (corresponding to a core radius $`r_c=5\stackrel{}{\mathrm{.}}2\pm 2\stackrel{}{\mathrm{.}}4`$) and $`\alpha =1.05\pm 0.05`$ ($`1\sigma `$ errors). These values are consistent with the previous results by LCG95 ($`r_c=6\stackrel{}{\mathrm{.}}7\pm 1\stackrel{}{\mathrm{.}}9`$ and $`\alpha =0.97\pm 0.15`$) and by AO89, who estimated a radius of the central plateau of $`3\mathrm{}`$ (from the analysis of their images) and $`\alpha =0.95`$ (from the brightness density profile). At a distance of $`D=4.3\pm 0.4`$ kpc (Ferraro et al. 1999b), the core radius we inferred corresponds to a physical dimension of $`0.11\pm 0.05`$ pc.

## 6 The dynamical status of NGC 6752

The significant deviation of the star number density profile from a canonical King model is a clear indication that the innermost region NGC 6752 has experienced (or is experiencing) a collapse phase. Unfortunately, no similarly unambiguous signature is available for differentiating the in- and post-core collapse state (Meylan & Heggie 1997). However, in favorable cases the phase of the collapse can be evaluated from indirect evidence.

Our results solidify the earlier suggestions of post-core collapse bounce made by AO89. Post-core collapse clusters are expected to undergo large amplitude oscillations in core size due to the gravothermal instability of collisional systems (e.g. Cohn et al. 1991). The oscillating core spends most of the time at near maximum size and a radial extension of $`0.11`$ pc is consistent with the maximum radius of the core predicted by the models of post-core collapse bounce. The parameter most commonly used in theoretical studies is the ratio of core radius to half-mass radius, $`r_c/r_h`$. Using our value of $`r_c=5.2\mathrm{}`$ and the half-mass radius of $`115\mathrm{}`$ from Trager et al (1995), we find a ratio $`r_c/r_h=0.045`$. For comparison, the multi-mass Fokker-Planck models for the post-collapse evolution of M15 presented by Grabhorn et al (1992, see their Fig. 5) reach a comparable value of $`r_c/r_h=0.034`$ during the maximally expanded state of the most extreme core bounces. The inclusion of primordial binaries in Monte Carlo simulations by Fregeau et al (2003) results in an even wider range of predicted $`r_c/r_h`$ in the post-collapse phase.

Furthermore, we have found that there is an intermediate region in the composite density profile of NGC 6752 which is well represented by a power law profile with a slope $`\alpha 1,`$ compatible with the steepness predicted by single mass models of expanding bouncing cores (LCG95).

It is also worth noting that NGC 6752 has apparently retained a substantial primordial binary population (Rubenstein & Bailyn, 1997); these binaries may play an important role in supporting the core and delaying the core-collapse event. In this respect, Ferraro et al (1999a) have suggested that some species originated from binary evolution could be used as possible tracers of the cluster dynamical evolution. In particular, the large BSS population recently found in M80 by Ferraro et al (1999a) might be the signature of a transient dynamical state, during which stellar interactions are delaying the core-collapse process leading to an exceptionally large population of collisional-BSS. On the other hand the BSS population found in the central region of NGC 6752 is small (Ferraro et al 2003, in preparation) perhaps indicating that NGC 6752 is in a different dynamical evolutionary state than M80. Maybe the binary population in NGC 6752 has not been “burned out” producing collisional BSS while that in M80 has.

## 7 The interpretation of the MSPs accelerations

NGC 6752 hosts 5 known millisecond pulsars (D’Amico et al. 2001a; D’Amico et al. 2002). The positions in the plane of the sky of three of them (PSRs B, D and E, all isolated pulsars) are close to the cluster center, as expected on the basis of mass segregation in the cluster. In particular PSR-B and E show large negative values of $`\dot{P}`$, implying that the pulsars are experiencing an acceleration with a line-of-sight component $`a_l`$ directed toward the observer and a magnitude significantly larger than the positive component of $`\dot{P}`$ due to the intrinsic pulsar spin-down (see e.g. Phinney 1992). What is the origin of such acceleration? Given the location of NGC 6752 in the galactic halo and knowing its proper motion (Dinescu et al. 1999), it is possible to calculate the contributions to $`\dot{P}`$ due to centrifugal acceleration, differential galactic rotation and vertical acceleration in the Galactic potential, all of them resulting negligible (D’Amico et al. 2002). Hence, the remaining plausible explanations of the observed negative $`\dot{P}`$ are: the accelerating effect of the cluster gravitational potential well or the presence of some close perturbator(s) exerting a gravitational pull onto the pulsars. Taking advantage from the new results presented in this paper we discuss in the follow viability and implications of these two possibilities.

### 7.1 Case (i): Overall effect of the GC potential well

The hypothesis that the line-of-sight acceleration of the MSPs with negative $`\dot{P}`$ is dominated by the cluster gravitational potential has been routinely applied to many globulars. In particular from this assumption a lower limit to the mean projected mass-to-light ratio in the $`V`$-band $`/_V`$ in the central region of M15 (Phinney 1992) and 47 Tucanae (Freire et al. 2003), yielded $`/_V`$$`\stackrel{>}{_{}}3`$ and $`\stackrel{>}{_{}}0.7`$ respectively. Following Phinney (1992), a lower limit to $`/_𝒱`$ in the inner regions of NGC 6752 is given by

$`\left|{\displaystyle \frac{\dot{P}}{P}}(\theta _{})\right|<\left|{\displaystyle \frac{a_{l,max}(\theta _{})}{c}}\right|`$
$`1.1{\displaystyle \frac{G}{c}}{\displaystyle \frac{M_{cyl}(<\theta _{})}{\pi D^2\theta _{}^2}}=5.1\times 10^{18}{\displaystyle \frac{}{_V}}\left({\displaystyle \frac{\mathrm{\Sigma }_V(<\theta _{})}{10^4\mathrm{L}_\mathrm{V}\mathrm{pc}^2}}\right)\mathrm{s}^1,`$
where $`\mathrm{\Sigma }_V(<\theta _{})`$ is the mean surface brightness within a line of sight subtended by an angle $`\theta _{}`$ with respect to the cluster center and $`M_{cyl}(<\theta _{})`$ is the mass enclosed in the cylindrical volume of radius $`R_{}=D\theta _{}.`$ This equation holds to within $``$10% in all plausible cluster models and is independent of cluster distance, except for the effects of extinction. Since $`E(BV)`$ is very small for NGC 6752 (= 0.04 according to Harris 1996) the latter is a negligibe affect for this cluster.

Using the observed $`\dot{P}/P`$ of PSR-B and PSR-E (D’Amico et al. 2002), combined with our accurate determinations of $`C_{\mathrm{grav}}`$ (the cluster center of gravity) and $`\mathrm{\Sigma }_V(<r)`$ (the mean surface brightness radial profile), it turns out that $`/_V\stackrel{>}{_{}}6`$–7 (see Figure 7) for the case of NGC 6752. D’Amico et al. (2002) obtained a slightly larger $`/_V\stackrel{>}{_{}}10`$ using published values of $`C_{\mathrm{grav}}`$ and $`\mathrm{\Sigma }_V(<r)`$ derived from medium resolution ground based observations only. The difference between the two estimates is mainly due to our new position of the cluster center of gravity. Despite a residual uncertainty $`0\stackrel{}{\mathrm{.}}7`$ on $`C_{\mathrm{grav}},`$ under the hypothesis that the line-of-sight acceleration of PSR-B and PSR-E are entirely due to the cluster gravitational potential, a lower limit of $`/_V\stackrel{>}{_{}}5.5`$ can be firmly established. It is obtained assuming that the two millisecond pulsars were just symmetrically located (and hence at the minimum projected distance) with respect to the actual center of gravity.

The sample of the core collapsed clusters shows typical values of the projected central mass-to-light ratio in the interval 2–3.5 (Pryor & Meylan 1993), although larger $`/`$ ratios can be obtained when a Fokker–Planck model fit is used (e.g. the case of M15 Dull et al. 1997, 2003). If we take $`/_V\stackrel{>}{_{}}3`$ the expected total mass located within the inner $`r_{,B}=0.08`$ pc of NGC 6752 (equivalent to the projected displacement of PSR-B from $`C_{\mathrm{grav}}`$) would be $`1200`$$`2000M_{}.`$ On the other hand, the observed $`/_V`$$`6`$–7 implies the existence of further $`1500`$$`2000M_{}`$ of low-luminosity matter segregated inside the projected radius $`r_{,B}.`$ This extra amount of mass could be constituted by a relatively massive black hole, like the $`1700_{1700}^{+2700}M_{}`$ black hole in the center of the globular cluster M15 recently proposed by Gerssen et al. (2003). However, the results presented in §4 and §5 show that there are significant differences in the kinematics and mass distributions of the central regions of M15 and NGC 6752 (we note here that because the distance to NGC 6752 is less than half that to M15, its inner core is more easily studied): (1.) HST imaging of M15 (Guhathakurta et al. 1996—WFPC2; Sosin & King 1997—FOC) shows no evidence for a compact core, at variance with our observations of the core of NGC 6752; (2.) The derived stellar density profiles of M15 have power law slopes consistent with $`\alpha =0.75`$ (expected from single mass models with a dominant central black hole; LCG95). If a $`\stackrel{>}{_{}}10^3M_{}`$ black hole resides in the inner region of NGC 6752, its gravitational influence would extend more than $`2\mathrm{}`$ from the center of the cluster and probably produce a central power-law cusp, which we do not observe.

A very high $`/_V6`$–7 could be also due to central concentration of dark remnants of stellar evolution, like neutron stars (NSs) and heavy $`1.0M_{}`$ white dwarfs (WDs) (as also proposed for M15 by Dull et al. 1997, 2003 and by Baumgardt et al. 2003). In this case, one can constrain the initial mass function (IMF) and/or the neutron star retention fraction $`f_{\mathrm{ret}}`$ in NGC 6752. Based on the current population of turn-off stars (in the mass interval $`0.6`$$`0.8M_{}`$), the estimated number of upper main sequence initially present in NGC 6752 is $`4000`$ (D’Amico et al. 2002). This assumes a Salpeter-like IMF ($`\alpha _{\mathrm{IMF}}=2.35`$ which is consistent with that measured by Ferraro et al. 1997b). If the low-luminosity matter observed in the central 0.08 pc were entirely due to $`1300`$ NSs of $`1.4M_{}`$, then $`f_{\mathrm{ret}}30`$% (a reasonable value for collapsed clusters \[Drukier 1996\]). Alternatively, $`/_V6`$$`7`$ can be explained by a Salpeter IMF if $`\stackrel{>}{_{}}20`$% of the total population of heavy $`1.0M_{}`$-WDs sank into the NGC 6752 core during the cluster dynamical evolution. Either scenario must also be consistent with the observed shape of the star density profile.

According to Cohn (1985), during the core collapse phase the surface density slope for the most massive component is expected to be $`\alpha =1.23`$ rather than the projected isothermal slope of $`\alpha =1.0`$. In the central cusp that forms during core collapse, the surface density slope of a component of stellar mass $`m`$ is given approximately by:

$$\alpha =1.89\frac{m}{m_d}+0.65$$

where $`m_d`$ is the stellar mass of the dominant component. If the luminosity profile is dominated by turnoff stars of mass $`m=0.8M_{}`$ and has a slope of $`\alpha =1.05`$, then the implied mass of the dominant, nonluminous component should be $`m_d=0.89M_{}`$ , i.e. somewhat more massive than the adopted turnoff mass. This argument does suggest that the central gravitational potential is not likely to be dominated by a large number of neutron stars, but heavy white dwarfs still remain a possibility.

Velocity dispersions provide a further constraint on the nature of the cluster potential well. The stars dominating the dynamics in the inner part of the cluster should have (see eq. 3.5 of Phinney 1992) a 1-dimensional central velocity dispersion $`\sigma _{v,0}\stackrel{>}{_{}}9`$$`10\mathrm{km}\mathrm{s}^1`$. This is compatible both with the very wide published $`2\sigma `$ interval for the $`\sigma _{v,0}`$ of NGC 6752 ($`2.1`$$`9.7\mathrm{km}\mathrm{s}^1`$, Dubath, Meylan & Mayor 1997) and with preliminary proper motion measurements of stars in the central part of NGC 6752, which would suggest a significantly higher one-dimension velocity dispersion $`\sigma _{v,0}12\mathrm{km}\mathrm{s}^1`$ and the existence of strong velocity anisotropies (Drukier et al. 2003)555If confirmed, such anisotropies would support the hypothesis of a relatively massive, probably binary black hole. Recent Fabry-Perot spectroscopy of single stars in NGC 6752 shows a flat velocity dispersion profile with typical one-dimension velocity dispersion $`\sigma 6`$–7$`\mathrm{km}\mathrm{s}^1`$ within the central $`1\mathrm{}`$ (Xie et al. 2002). This is also marginally compatibile (given the $`10`$% uncertainty of the formula 3.5 of Phinney 1992).

### 7.2 Case (ii): Local perturbator(s)

We here explore the alternate possibility that the acceleration imparted to PSR-B and PSR-E are due to some local perturbator.

NGC 6752 is a highly concentrated cluster and its core could host $`n_{}\stackrel{>}{_{}}10^6`$ stars per cubic pc (assuming an average stellar mass of $`\overline{m}0.5M_{}`$). Hence close star-star encounters are a viable possibility. In order to produce the line-of-sight acceleration $`|a_l|=c|\dot{P}/P|=2.9\times 10^6\mathrm{cm}\mathrm{s}^2`$ seen in PSR-B (or PSR-E) a passing-by star of mass $`\overline{m}`$ must approach the pulsar at $`s(G\overline{m}|a_l|^1)^{1/2}=0.0015(\overline{m}/0.5M_{})^{1/2}`$ pc. An upper limit to the probability of occurrence of a suitable close encounter can be roughly estimated as $`s^3n_{}=3.7\times 10^3(\overline{m}/0.5M_{})^{3/2}.`$ Although this figure is not negliglible, the need of having two different canonical stars independently exerting their gravitational pull onto PSR-B and PSR-E, makes the joint probability of such configuration suspiciously low, $`\stackrel{<}{_{}}10^5`$. Rather than being randomly placed could the perturber be a binary companion to the MSPs? For such a companion not to have been aleady discovered by pulsar timing analysis would require an orbital period $`P_b\stackrel{>}{_{}}20`$ yr. Survival of such a wide binary is quite problematic in the core of a dense cluster like NGC 6752—indeed no binary pulsar with $`P_b>3.8`$ day has been detected in collapsed clusters to date. Seeing two such systems in NGC 6752 appears extremely unlikely.

One may wonder if a single object, significantly more massive than a typical star in the cluster, could simultaneously produce the accelerations detected both in PSR-B and PSR-E. Recently, Colpi, Possenti & Gualandris (2002) suggested the presence of a binary black-hole (BH) of moderate mass ($`M_{\mathrm{bh}+\mathrm{bh}}100`$$`200M_{}`$) in the center of NGC 6752 in order to explain the unprecedented position of PSR-A, a binary millisecond pulsar which is far away from the cluster center (D’Amico et al 2002). As shown in Figure 2, the projected separation of PSR-B and PSR-E is only $`d_{}=0.03`$ pc. A binary BH of total mass $`M_{\mathrm{bh}+\mathrm{bh}}`$, approximately located in front of them within a distance of the same order of $`d_{}`$, could be accelerating both the pulsars without leaving any observable signature on the photometric profile of the cluster. As the BH binarity ensures a large cross section to interaction with other stars, the recoil velocity $`v_{\mathrm{rec}}`$ due a recent dynamical encounter could explain the offset position (with respect to $`C_{\mathrm{grav}}`$) of the BH. However this scenario suffers of a probability at least as low as that of the previous two: placing the black hole at random within the core gives roughly a 1% chance that it would land in a location where it would produce the observed pulsar accelerations and the required $`v_{\mathrm{rec}}\stackrel{>}{_{}}4`$$`5\mathrm{km}\mathrm{s}^1`$ is also at the upper end of the expected distribution of $`v_{\mathrm{rec}}`$ for an intermediate mass black-hole in NGC 6752 (Colpi, Mapelli & Possenti 2003, in preparation).

In summary, on the basis of the available data, the accelerations shown by PSR-B and PSR-E can be easilier accounted by the usually adopted hypothesis (i) (which could explain also the large positive $`\dot{P}`$ of PSR-D if it were not instrinsic), although the nature of the required extra amount ($`1500`$$`2000M_{}`$) of low luminosity mass still remains puzzling. The existence of local perturber(s) of the pulsar dynamics (case (ii)) is a distinct possibility. While it seems extremely unlikely purely on the basis of the high value of $`\dot{P}/P`$ measured in PSR-B and PSR-E, there are other indications of a binary low-mass BH. In particular, it would explain: (a) the absence of any cusp in the radial density profile, (b) the flat velocity dispersion profile (Xie et al. 2002), (c) and (d) the ejection of PSR-A and PSR-C in the cluster outskirts (Colpi, Possenti & Gualandris 2002). However it is admittedly ad hoc hypothesis, requiring a fine tuned scenario.

Clearly additional information must be collected. A longer baseline for timing measurements will allow to better constrain the presence of companions in very large orbits around PSR-B and PSR-E. In particular it will permit to derive (at least) upper limits on the second derivative of the pulsars’ spin period, which is more influenced by by-passing stars rather than by the cluster potential well (Phinney 1992). Similarly, the (single massive or binary intermediate mass) black-hole hypotheses could be better investigated with spectroscopic determination of the dispersion velocity of the stars located in the pulsars’ neighborhood, likewise with star density counts reaching higher magnitudes (hence exploiting a larger sample of objects) and better spatial resolution (thus probing the innermost $`1\mathrm{}`$ of the cluster and the pulsars’ surroundings).

We thank Phyllis Lugger for useful comments and suggestions that significantly improved the presentation of the paper. We also thank P. Montegriffo for assistance with the astrometry procedure, E. Pancino for the pre-reduction of the WFI images and S. Sigurdsson for allow us to use his code to compute King model surface density profile. Bingrong Xie kindly provided results prior to publication. Financial support to this research is provided by the Agenzia Spaziale Italiana (ASI) and the Ministero dell’Istruzione, dell’ Università e della Ricerca (MIUR). The GSCII catalog was produced by the Space Telescope Science Institute and the Osservatorio Astronomico di Torino. RTR is partially supported by STScI grant GO-8709 and NASA LTSA grant NAG5-6403.


text:
# Depletion forces in non-equilibrium

## Abstract

The concept of effective depletion forces between two fixed big colloidal particles in a bath of small particles is generalized to a non-equilibrium situation where the bath of small Brownian particles is flowing around the big particles with a prescribed velocity. In striking contrast to the equilibrium case, the non-equilibrium forces violate Newton’s third law, are non-conservative and strongly anisotropic, featuring both strong attractive and repulsive domains.

Colloidal mixtures are excellent model systems to study fundamental questions of equilibrium and non-equilibrium phase transitions with important physical applications, such as protein crystallization Lekk  as well as gel and glass formation Fuchs . One of the reasons for the pivotal role of colloids as model systems lies in the fact that their interactions can be tailored Pusey . In the theoretical description of colloidal mixtures, the concept of the effective interaction has proved to be of paramount importance in advancing our understanding of how to steer the properties of these complex systems poon:jpcm:02 . It stems from integrating out the coordinates of the small particles, reducing thereby the number of statistical degrees of freedom considerably likos:physrep:01 ; Springer . This approach has been successfully applied to colloid/polymer mixtures AOmodel ; Lekkerkerker ; jusufi:jpcm:01 , to binary colloidal mixtures evans:pre:99 ; roth:pre:00  and to charged suspensions irene:prl:98 . The most widely known effect is the depletion attraction in strongly non-additive systems, such as the Asakura-Oosawa model AOmodel  pertaining to colloids and non-adsorbing ideal polymer coils. This attraction is entropic, i.e., it scales with the thermal energy $`k_BT`$ and results from a depletion zone of polymer between a pair of nearly touching colloidal particles, which causes an imbalance of the osmotic pressure of the polymers acting on the colloidal surface. While the range of the depletion attraction is typically associated with the radius of gyration of the polymers, its attractive depth is fixed by the polymer concentration, so that both the range and the depth are in principle tunable.

The crucial limitation of the concept of effective interactions is that it works, strictly speaking, only in thermodynamic equilibrium and can thereafter be applied for the calculations of static quantities, such as correlation functions and free energies (resulting into phase diagrams). In this Letter, we investigate whether and how the effective interaction can be applied also in non-equilibrium steady state situations. In order to be as transparent as possible, we study a minimal model of non-equilibrium for two fixed big colloidal particles in a quiescent solvent which are exposed to a flowing bath of small Brownian particles. The latter ones are ideal, in full analogy to to the equilibrium case of the Asakura-Oosawa model, and hydrodynamic interactions mediated by the solvent are neglected. Despite its simplicity, the model is realized in colloidal mixtures with weakly interacting depletants (which justifies the assumed ideality) which have a small physical volume fraction with well-separated excluded volumes of the big and small particles (which justifies neglecting hydrodynamic interactions). Examples include strongly asymmetric highly charged suspensions with a small physical volume fraction palberg:jcp:01  and sterically-stabilized colloidal particles in a bath of thin rods koenderink:prl:2003  or polymer coils verma:prl:98 ; rudhardt:prl:98 . Our situation is different to the classical hydrodynamic problem of two fixed spheres in a Stokes flow of a viscous incompressible fluid, felderhof:physica:1977 ; hinsen:jcp:94 ; jan:book  in which the force vanishes for zero flow. More importantly, the physical origin of the depletion forces is the high osmotic compressibility of the bath particles complementary to the incompressibility of the viscous fluid mediating hydrodynamic interactions.

Based on a standard Smoluchowski approach and on Brownian dynamics computer simulations in non-equilibrium, we show that the effective force can be generalized in a straightforward way to a stationary non-equilibrium state. The results for the effective force, however, differ from the corresponding equilibrium case both qualitatively and quantitatively. First of all, the non-equilibrium force field (as a function of the mutual distance between the particles) is non-conservative and violates Newton’s third law, implying that the concept of an effective interaction potential does not hold in non-equilibrium. Second, when the small particles are brought to flow, the equilibrium depletion attraction gets strongly anisotropic favoring an alignment of the colloidal pair along the flow direction. Conversely, there is mutual effective repulsion between the colloidal pair perpendicular to the drift direction of the bath due to a compression of the small particles between the colloids. Furthermore, the range of the depletion interaction in the drifting case is larger than that in equilibrium.

Equilibrium depletion forces can be measured in real colloidal samples with high accuracy koenderink:prl:2003 ; verma:prl:98 ; rudhardt:prl:98 ; bechinger:prl:99 ; crocker:prl:99 . Our predictions in non-equilibrium can be tested in experiments of binary colloidal mixtures as well. A drifting bath is generated either by gravity philipse:97 ; aarts:jpcm:03  or by an electric or magnetic field palberg:jcp:01  and the big particles can be fixed by optical tweezers crocker:prl:99 ; yodh:01 . The measurement of the relative displacement with respect to the tweezer’s position of the big particles provides direct access to the effective non-equilibrium forces.

In our theoretical model we consider two fixed big particles in a bath of small Brownian particles that drift with velocity $`𝐯=v\widehat{𝐳}`$. The bath has a bulk number density $`\rho _0`$. While there is no interaction between the bath particles, the interaction between the big and small particles is described with an exponential, spherically symmetric potential $`V(𝐫)`$ given by

$`V(𝐫)=V(r)=V_0\mathrm{exp}\left[(r/\sigma )^n\right],`$
where $`𝐫`$ is the vector separating the bath particle with the center of the colloid, $`r`$ is the magnitude of $`𝐫`$ and $`V_0=10k_BT`$. In this work we focus on an exponent $`n=6`$, on one hand to deal with a well-defined ‘hard’ interaction diameter of the colloidal sphere, $`\sigma `$, on the other hand to circumvent numerical inconveniences due to huge forces at $`r=\sigma `$ for high exponents $`n`$.

In what follows, we first discuss a theoretical Smoluchowski approach for a single driven colloid. Adopting a dynamical superposition approximation (DSA), we then address the case of two colloids and afterwards we compare the theoretical predictions arising from the DSA with Brownian dynamics (BD) computer simulations, finding good agreement between the two approaches. In detail, for the case of a single big colloidal particle, the Smoluchowski equation jan:book  for the steady state density field $`\rho (𝐫)`$ of the bath particles reads as $`\stackrel{}{}𝐣(𝐫)=0`$ with the current density field $`𝐣(𝐫)`$ composed of

$`𝐣(𝐫)={\displaystyle \frac{\mathrm{\Gamma }}{\beta }}\stackrel{}{}\rho (𝐫)+\mathrm{\Gamma }\rho (𝐫)\stackrel{}{}V(𝐫)+𝐯\rho (𝐫),`$
where $`\beta =1/k_BT`$ is the inverse thermal energy and $`\mathrm{\Gamma }`$ the mobility of the bath particles in the solvent. As the problem possesses azimuthal symmetry, these equations can be readily numerically solved in cylindrical coordinates. Let us define the resulting density profile of the bath particles around one big colloid as $`\rho ^{(1)}(𝐫𝐑;𝐯)`$ with $`𝐑`$ denoting the position of the colloid. Spatial homogeneity dictates that $`\rho ^{(1)}`$ depends only on the difference $`𝐫𝐑`$. The force $`𝐅_\mathrm{d}`$ acting on the big particle due to the drifting bath is then obtained as $`𝐅_\mathrm{d}=\mathrm{d}^3r\rho ^{(1)}(𝐫)V(|𝐫|)`$.

In the presence of two colloids, the corresponding quantity of interest is the bath density profile $`\rho ^{(2)}(𝐫,𝐑_1,𝐑_2;𝐯)`$ that depends on the locations $`𝐑_i`$ of the colloids, ($`i=1,2`$), and parametrically on the velocity $`𝐯`$. In this case, the numerical solution of the Smoluchowski equation is more difficult such that a direct particle-resolved computer simulations is more appropriate. However, let us motivate a simple superposition approximation. As the bath particles are noninteracting, in the case of thermodynamic equilibrium, $`v=0`$, $`\rho ^{(2)}(𝐫,𝐑_1,𝐑_2;0)`$ can be exactly factorized as the product of the density fields $`\rho ^{(1)}(𝐫𝐑_i;0)`$ divided by $`\rho _0`$. In what follows, we employ the same factorization ansatz also for the non-equilibrium situation, $`v0`$. Thereby, we introduce the dynamical analog of the superposition approximation that was first introduced by Attard attard:jcp:89  for equilibrium problems dealing with interacting particles, in order to calculate static depletion forces. Hence, in this dynamical superposition approximation (DSA), we write:

$`\rho ^{(2)}(𝐫,𝐑_1,𝐑_2;𝐯)\rho ^{(1)}(𝐫𝐑_1;𝐯)\rho ^{(1)}(𝐫𝐑_2;𝐯)/\rho _0.`$ (1)
Results from the DSA are shown in Fig. 1, where we plot the density distribution in the $`x`$-$`z`$ plane, on which the centers of both colloids are located. The profiles are plotted for four different relative positions of the two colloids and for a dimensionless drift velocity $`v^{}=5`$, defined as $`v^{}=\beta \sigma v/\mathrm{\Gamma }`$. The density profiles are strongly anisotropic and show long-ranged low-density regions in the wake region behind the spheres. The extension of the concept of the effective interaction to a steady-state non-equilibrium situation can now be performed via the definition of the effective force which is a statistical average over the forces of the small bath particles exerted onto the big ones likos:physrep:01 . Generalizing to a steady-state average in non-equilibrium, the effective force reads as footnote :

$`𝐅_i={\displaystyle \mathrm{d}^3r\rho ^{(2)}(𝐫,𝐑_1,𝐑_2;𝐯)\stackrel{}{}_{𝐑_i}V(|𝐑_i𝐫|)}.`$ (2)
The force has been calculated in theory by employing the DSA, Eq. (1), for various different configurations of the colloids and drift velocities, including the equilibrium case, $`v=0`$, in which the results are exact. In this way, a comparison between the static and dynamic depletion force $`𝐅_{\mathrm{depl},i}`$ acting on the $`i`$-th particle can be obtained. The latter is defined as the force $`𝐅_{\mathrm{depl},i}=𝐅_i𝐅_\mathrm{d}`$, relative to the constant the drift force $`𝐅_\mathrm{d}`$ acting on a single big particle ($`\beta F_\mathrm{d}\sigma =2.96`$ for $`v^{}=1`$ and $`\beta F_\mathrm{d}\sigma =11.76`$ for $`v^{}=5`$). We accompanied the DSA with standard Brownian Dynamic simulations Allen , in which the Langevin equations for each bath particle are solved iteratively. Here, the bath particles are driven by an external force $`𝐯/\mathrm{\Gamma }`$ for a prescribed flow velocity $`𝐯`$. The box lengths in all simulations were $`L_x=L_y=10\sigma `$ and $`L_z=20\sigma `$ and $`N=2000`$ particles were simulated, applying periodic boundary conditions in all directions.

Results for the averaged force field acting on a colloid placed at the position $`𝐑_2`$ away from another colloid at the origin ($`𝐑_1=0`$), are depicted in Fig. 2. A comparison between Figs. 2(c) and 2(d), (DSA vs. BD-simulation), shows that the former captures the features of the non-equilibrium depletion force quantitatively, even for drifting speeds as large as $`v^{}=5`$. Whereas for the case of equilibrium, Fig. 2(a), we obtain spherically symmetric attractive depletion forces, the dynamical depletion force shows new qualitative features: First, as is clear from inspection of Fig. 2, the dynamical depletion forces violate Newton’s third law. Second, the violation of the “$`\mathrm{action}=\mathrm{reaction}`$”-principle implies a non-conservative force. Indeed, if a putative depletion potential $`V_{\mathrm{depl}}`$ existed in the dynamical case at hand, translational invariance of space (which is not broken by the steady flow) would imply $`V_{\mathrm{depl}}=V_{\mathrm{depl}}(𝐑_1𝐑_2)`$. As $`𝐅_{\mathrm{depl},i}=\stackrel{}{}_{𝐑_i}V_{\mathrm{depl}}(𝐑_1𝐑_2)`$, $`i=1,2`$, this would have the consequence $`𝐅_{\mathrm{depl},1}=𝐅_{\mathrm{depl},2}`$. Third, the depletion force shows strong anisotropy, becoming attractive as the colloid starts positioning itself ‘behind’ the one at the origin and its magnitude (in the attractive case) and its range are vastly larger for $`v0`$ than in equilibrium.

In order to better quantify the features of the depletion force in non-equilibrium, we show in Fig. 3 plots of the projected depletion force at fixed $`X=0`$ (particles fully aligned along the flow direction) and at fixed $`Z=0`$ (particles fully aligned perpendicular to the flow direction). The depletion force exhibits a repulsive barrier in the $`X`$-direction due to compression of the small particles between the colloids, in qualitative difference to the purely attractive equilibrium case $`v=0`$. In the $`Z`$direction, on the other hand, there is much stronger attraction for $`Z>0`$ due to a long low density region in the wake of the colloids but a marked repulsive barrier for $`Z<0`$. Furthermore, the range of the depletion force is significantly larger than in equilibrium. Notice that we show forces also for the case of completely overlapping colloids footnote , $`2.0X/\sigma 2.0`$ and $`2.0Z/\sigma 2.0`$. For $`𝐑_1=𝐑_2`$, the depletion force on both colloids is obviously the same, $`𝐅_{\mathrm{depl},1}=𝐅_{\mathrm{depl},2}`$ but, in contrast to the equilibrium case, it does not vanish. We can finally deduce from Fig. 3 that the computer simulation results are in good quantitative agreement with the theory up to drifting speeds of $`v^{}=1`$, but for large drifting speeds as high as $`v^{}=5`$ deviations near the maxima and minima are getting more pronounced.

In conclusion, we have generalized the depletion concept to a simple non-equilibrium situation of two big particles in a Brownian flow of small ideal particles. The effective forces averaged in the steady state are strikingly different from those in equilibrium: they are strongly anisotropic with both attractive and repulsive domains and violate Newton’s third law, in strong resemblance to the so-called “social forces” that have been applied extensively in the modelling of pedestrian and automotive traffic flows helbing:book .

Although the present study has focused on noninteracting bath particles, the extension to interacting ones can be carried out within the framework of the recently proposed dynamical density functional theory (DDFT) marconi:jcp:99 . We expect that interaction effects will weaken the strength of the dynamical depletion forces, as they will tend to smoothen out regions of spatial inhomogeneity of the bath density that give rise to these forces in the first place. Finally, we remark that the strong attractive forces tending to align the colloids in the direction of the incoming flow can be regarded as precursors of the laning transition that has been analyzed recently for concentrated colloidal mixtures joe:laning1 .

We thank R. Finken, G. Nägele, and H. Stark for useful discussions. Financial support from the EPSRC and from the DFG through the SFB-TR6 are gratefully acknowledged. C.N.L. has been supported by the DFG through a Heisenberg Fellowship and wishes to thank the Department of Chemistry, University of Cambridge, and Prof. J.-P. Hansen for their hospitality.


text:
# 1 Introduction

## 1 Introduction

The fact that Planck mass, $`m_p=(\mathrm{}c/G)^{1/2}10^5`$ gms, is much larger than elementary particle masses, has sometimes led to the suggestion  that standard, linear quantum mechanics may not hold for idealized point particles whose mass approaches the Planck mass. In particular, it is possible that the Schrodinger equation for the free particle may have to be modified to include non-linear terms, when the particle’s mass becomes comparable to $`m_p`$. This is the possibility considered in this brief note.

From the experimental viewpoint, Planck mass is however, embarrasingly, in the macroscopic mass range, and its more likely that the fundamental mass scale $`m_c`$ where any non-linearities might set in is a few orders of magnitude smaller than $`m_p`$. The reason how this happens is not clear, though it is plausible that an underlying theory relates $`m_c`$ directly to $`m_p`$.

Even if $`m_c`$ is a few orders of magnitude smaller than Planck mass, there is still a very huge gap between elementary particle masses and $`m_c`$. It is significant that quantum mechanics has not been experimentally tested for a point particle with a mass $`m`$ somewhere in this range, say $`m=10^{15}`$ gms. The purpose of this note is to give one example of a mass-dependent non-linear modification of the Schrodinger equation, with a view to encouraging experimental searches of a possible mass-dependent non-linearity.

Considerations of non-linear generalizations of the Schrodinger equation have a long history, an early work being that of Bialynicki-Birula and Mycielski ; more recent discussions include those by Weinberg , Zeilinger , Pearle , Shimony  and Minic et al. . A more complete list of references can be found in the recent work of Adler . A mass dependent non-linearity does not seem to have been considered before, though.

## 2 A non-linear Schrodinger equation

Let us rewrite the free-particle Schrodinger equation

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}\frac{^2\psi }{q^2}$$
(1)
by defining $`\psi =e^{iS(q,t)/\mathrm{}}`$, where $`S(q,t)`$ is complex. This gives the following equation for $`S(q,t)`$:

$$\frac{S}{t}=\frac{1}{2m}\left(\frac{S}{q}\right)^2+\frac{i\mathrm{}}{2m}\frac{^2S}{q^2}.$$
(2)
In the so-called $`\mathrm{}0`$ limit the last term of this equation is dropped, and one recovers the classical Hamilton-Jacobi equation for the action function $`S(q,t)`$.

Let us now assume that Eqn. (2) has mass-dependent corrections - these corrections are assumed to be vanishingly small for masses which are much smaller or much larger than the fundamental mass scale $`m_c`$, but become significant when $`m`$ is comparable to $`m_c`$.

It is obvious that the equation (2) for the generalized action $`S`$ is non-linear, the non-linearity being introduced by the term quadratic in $`S/q`$. We assume that the proposed mass-dependent modification of this equation is such that precisely for $`m=m_c`$ equation (2) actually becomes linear! As $`m`$ moves away from $`m_c`$ (becoming either larger or smaller than $`m_c`$) the non-linear term in (2) is recovered. Thus, we propose the following modified equation, by introducing the modifying function $`A(m/m_c)`$:

$$\frac{S}{t}=\frac{1}{2m}\left(1+A\left(\frac{m}{m_c}\right)\right)\left(\frac{S}{q}\right)^2+\frac{i\mathrm{}}{2m}\frac{^2S}{q^2}.$$
(3)
The function $`A(m/m_c)`$ takes the value one at $`m=m_c`$ and goes to zero when $`m`$ is much larger or much smaller than $`m_c`$. Hence, for $`m=m_c`$ the modified equation for $`S(q,t)`$ is

$$i\mathrm{}\frac{S}{t}=\frac{\mathrm{}^2}{2m}\frac{^2S}{q^2}$$
(4)
which has the same form as the Schrodinger equation!

In Eqn. (3), after substituting for $`S`$ in terms of $`\psi `$, one gets the following non-linear Schrodinger equation for $`\psi (q,t)`$:

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}\frac{^2\psi }{q^2}+\frac{\mathrm{}^2}{2m}A(m/m_c)\left(\frac{\mathrm{ln}\psi }{q}\right)^2\psi .$$
(5)
Consider the implication of this non-linearity for the special case $`m=m_c`$. Let $`\psi _1`$ and $`\psi _2`$ be two solutions of the Schrodinger equation (1), so that the superposition $`\psi =a\psi _1+b\psi _2`$ also solves (1). Since (4) has the same structure as (1), one can choose solutions for the generalised action $`S`$ in (4) as $`S_1=\psi _1`$ and $`S_2=\psi _2`$. The non-linear Schrodinger equation thus has two solutions $`\psi _{1NL}=e^{i\psi _1/\mathrm{}}`$ and $`\psi _{2NL}=e^{i\psi _2/\mathrm{}}`$ which have the same information content as $`\psi _1`$ and $`\psi _2`$ but the superposition of $`\psi _{1NL}`$ and $`\psi _{2NL}`$ does not satisfy the non-linear equation (5). If one regards the quantum measurement process as a sudden jump in the mass of a system from a value $`mm_c`$ to a value $`mm_c`$, the breakdown in superposition induced by the non-linearity could possibly explain wave-function collapse. This is because during the measurement $`m`$ ‘passes‘ through the intermediate value $`m_c`$. Here, before the measurement, by system one means the quantum system (hence $`mm_c`$) and after the measurement by system one means \[measuring apparatus + quantum system\] (hence $`mm_c`$).

A plane wave $`\psi =\alpha \mathrm{exp}i(\omega tkq)`$ satisfies the non-linear equation (5) with the dispersion relation

$$\mathrm{}\omega =\frac{\mathrm{}^2}{2m}k^2(1A)$$
(6)
but a superposition of plane-waves does not satisfy this equation.

In Eqn. (3) let us make the substitutions

$$\stackrel{~}{q}=\frac{q}{\sqrt{1A(m/m_c)}},\mathrm{}_N=\frac{\mathrm{}}{1A(m/m_c)}.$$
(7)
This changes (3) to the form

$$\frac{S}{t}=\frac{1}{2m}\left(\frac{S}{\stackrel{~}{q}}\right)^2+\frac{i\mathrm{}_N}{2m}\frac{^2S}{\stackrel{~}{q}^2}.$$
(8)
This has the same form as (2), with $`\stackrel{~}{q}`$ replacing $`q`$, and $`\mathrm{}_N`$ replacing $`\mathrm{}`$. It is evident that if we now define $`\psi _N=e^{iS/\mathrm{}_N}`$ then (8) is transformed to the following linear equation

$$i\mathrm{}_N\frac{\psi _N}{t}=\frac{\mathrm{}_N^2}{2m}\frac{^2\psi _N}{\stackrel{~}{q}^2}$$
(9)
which has the same form as the Schrodinger equation, with $`\mathrm{}_N`$ replacing $`\mathrm{}`$, and $`\stackrel{~}{q}`$ replacing $`q`$. It is also easily seen that $`\psi _N`$ is related to the original wave-function $`\psi `$ as

$$\psi _N=e^{iS/\mathrm{}_N}=\left(e^{iS/\mathrm{}}\right)^{\mathrm{}/\mathrm{}_N}=\psi ^{1A(m/m_c)}.$$
(10)
Equivalently, one can verify that the non-linear equation (5) reduces to the linear equation (9) after the substitutions given in (7) and (10) are made. The wave-function $`\psi _N`$ satisfies the continuity equation, with the modified Planck’s constant $`\mathrm{}_N`$ replacing $`\mathrm{}`$, and with $`\stackrel{~}{q}`$ replacing $`q`$. In the limit that $`A`$ goes to zero, $`\psi _N`$, $`\mathrm{}_N`$ and $`\stackrel{~}{q}`$ reduce to $`\psi `$, $`\mathrm{}`$ and $`q`$ respectively.

There appears to be a certain kind of singularity at $`m=m_c`$ which is possibly of fundamental importance.

## 3 Corrections to energy levels of the harmonic oscillator

The Schrodinger equation for the harmonic oscillator

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}\frac{^2\psi }{q^2}+\frac{1}{2}m\omega ^2q^2\psi $$
(11)
is transformed, after the substitution $`\psi =e^{iS/\mathrm{}}`$, to

$$\frac{S}{t}=\frac{1}{2m}\left(\frac{S}{q}\right)^2+\frac{i\mathrm{}}{2m}\frac{^2S}{q^2}\frac{1}{2}m\omega ^2q^2.$$
(12)
Including the proposed mass-dependent correction, as in (3), changes this equation to

$$\frac{S}{t}=\frac{1}{2m}\left(1+A\left(\frac{m}{m_c}\right)\right)\left(\frac{S}{q}\right)^2+\frac{i\mathrm{}}{2m}\frac{^2S}{q^2}\frac{1}{2}m\omega ^2q^2.$$
(13)
The transformations (7) along with the re-definition $`\stackrel{~}{\omega }=\omega \sqrt{1A}`$ give

$$\frac{S}{t}=\frac{1}{2m}\left(\frac{S}{\stackrel{~}{q}}\right)^2+\frac{i\mathrm{}_N}{2m}\frac{^2S}{\stackrel{~}{q}^2}\frac{1}{2}m\stackrel{~}{\omega }^2\stackrel{~}{q}^2.$$
(14)
The substitution $`\psi _N=e^{iS/\mathrm{}_N}`$ then gives

$$i\mathrm{}_N\frac{\psi _N}{t}=\frac{\mathrm{}_N^2}{2m}\frac{^2\psi _N}{\stackrel{~}{q}^2}+\frac{1}{2}m\stackrel{~}{\omega }^2\stackrel{~}{q}^2\psi _N.$$
(15)
Noting that the time-dependent part of a stationary state $`\psi `$ goes as $`e^{iEt/\mathrm{}}`$ we can see that

$$i\mathrm{}_N\frac{\psi _N}{t}=E\psi _N$$
(16)
and hence that the non-linear correction modifies the energy levels to

$$E=(n+1/2)\mathrm{}_N\stackrel{~}{\omega }=(n+1/2)\frac{\mathrm{}\omega }{\sqrt{1A}}.$$
(17)
## 4 Discussion

The particular proposal for non-linearity given here is of course tentative; the main purpose of this note being to highlight the possibility of a mass-dependent non-linear correction, which should be probed by experiments.

The possibility of such corrections is also suggested by earlier work  on a spacetime-independent formulation of quantum mechanics. The structure of the formulation for a free particle suggests that if the mass of the particle is close to $`m_c`$ there should be non-linear corrections to the Schrodinger equation. The exact nature of these non-linear terms (for instance the functional form of $`A(m/m_c)`$) can only be known when a better understanding of spacetime-independent quantum mechanics is achieved. Furthermore, as discussed in , one expects the introduction of non-linearities to be accompanied by a change in the fundamental commutation relations.

It is known  that non-linearities in a deterministic modification of the Schrodinger equation can allow for superluminal communication, and this has sometimes been presented as an argument against such a modification. But one could well speculate that such a radical change in our basic ideas might actually be occurring in the experimentally unprobed domain around $`m=m_c`$. Physicists are not averse to the possibility of a modification in quantum mechanics at the Planck energy scale $`E_p10^{19}`$ GeV. By virtue of the equivalence between mass and energy, the possibility of such a modification should also be entertained around the Planck mass scale $`m_p=E_p/c^2`$.

I would like to thank C. S. Unnikrishnan for useful discussions which lead to a significant modification of an earlier version of the paper.


text:
# Spatiotemporal Behavior of Void Collapse in Shocked Solids

## Abstract

Molecular dynamics simulations on a three dimensional defective Lennard-Jones solid containing a void are performed in order to investigate detailed properties of hot spot generation. In addition to the temperature, I monitor the number of energetically colliding particles per unit volume which characterizes the intensity of shock-enhanced chemistry. The quantity is found to saturate for nanoscale voids and to be maximized after voids have completely collapsed. It makes an apparent comparison to the temperature which requires much larger void for the enhancement and becomes maximum during the early stage of the collapse. It is also found that the average velocity and the temperature of ejected molecules inside a cubic void are enhanced during the collapse because of the focusing of momentum and energy towards the center line of a void.

In liquid or solid explosives, it is known that heterogeneities such as bubbles or voids decrease a threshold shock strength above which a detonation takes place. When a shock passes through such defects, they form locally excited regions which are reffered to as “hot spots” dear . Since chemical reactions which turn shocks to detonations are believed to initiate inside such hot spots, the mechanism for the generation of hot spots is important for the understanding of detonations cheret . Since shock waves and void collapse are extremely fast microscopic phenomena, molecular dynamics simulation may be a powerful tool for the study of detailed properties of hot spot initiation. As well as the studies which treat perfect crystals using empirical reactive potentials brenner , there are a few works on the interaction between shocks and defects which can yield hot spots maffre ; phillips ; phillips2 ; mintmire . For example, Phillips speculated that chemistry would occur in a void below a detonation threshold phillips2 . However, it is yet to be understood what is essential to the initiation of hot spots. Holian et al. recently studied a two dimensional Lennard-Jones solid with a planar gap which might be regarded as a gap of infinite transverse width holian . They insisted that the essential ingredients to hot spot generation are: (i) Enough shock strength for the ejection (spallation) of solid molecules into a void. (ii) Compression of ejected particles by the far side of a void. They estimated the overshooting temperature via a straightforward model constructed on the basis of the above conjectures. However, since their estimation is not quantatively satisfactory (about $`3/2`$ of those obtained by the simulations), the validity of their proposal is not apparent. Moreover, since void collapse is so fast that particles ejected into a void may be far from local equilibrium, it is no longer obvious that the “temperature” could be a fair index for chemistry. We have to monitor other quantities to represent the degree of enhanced chemistry in this situation.

In this Letter, I focus on a dynamical aspect of hot spot generation to understand what is essential to shock-induced chemistry inside a void. Especially, various spatial profiles of the ejected particles are monitored in order to see how the void collapse results in energetic intermolecular collisions which are responsible for chemistry. As for the model, I adopt a three dimensional fcc crystal whose intermolecular potential is the Lennard-Jones one;

$$U(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right],$$
(1)
whose cut-off length is set to be $`4.0\sigma `$. Throughout this Letter, $`\sigma `$ and $`ϵ`$ are set to unity (the LJ unit). Although one might wonder about the use of nonreactive potential to this kind of issue, this may not be a serious limitation for the understanding of dynamical aspect of shock-induced chemistry; i.e. how the coherent (shock) energy relaxes into the thermal energy through the void collapse, which is the main aim of this study. The setup of the system is as follows. The axes of $`x`$, $`y`$, and $`z`$ are taken to coincide with $`100`$, $`010`$, and $`001`$ directions, respectively. Although there are various methods to realize shock waves in molecular dynamics simulations holian2 , the most popular method is adopted here that a moving piston of infinite mass hits the still target. The velocity of the piston is denoted by $`u_p`$. Throughout this study, $`u_p`$ is fixed to be $`3`$ which is strong enough to cause spallation. Shock propagates along the $`001`$ orientation; i.e. $`z`$ axis. I do not argue orientational dependence which is important to other shock-induced phenomena, such as plastic deformation germann  or phase transitions zhakhovskii . The whole system consists of $`200\times 200\times 200`$ unit cells. As usual, periodic boundary conditions are applied to $`x`$ and $`y`$ directions. A void consists of blank unit cells, the shape of which is cubic or rectangular-solid. Each face is parallel to $`\{100\}`$, $`\{010\}`$, and $`\{001\}`$, respectively. Its typical size is ranged from $`5`$ to $`40`$ unit cells for each direction. I fix the initial temperature to be $`0.1`$ and the corresponding density is about $`1.053`$.

Before proceeding to the results, I make a remark on the two kinds of averaging methods adopted in this study. One is the average regarding the particles whose initial distance from the upstream wall of the void is less than $`5`$ (about 3 unit cells). Note that all of these particles are eventually ejected into the void after the passage of shock front. I have crosschecked the simulations with different distances from $`2`$ to $`15`$, and the qualitative property obtained below is not influenced by a specific choice of distance. By this method, we can trace the most energetic part of the ejected particles. Most of the results except for Fig. 4 are based on this averaging method. The other method is to take spatial profiles inside a void. This method enables us to see where and when the hot spot generation takes place.

Now let us turn to the results of my simulations. First, time evolutions of the temperature during void collapse are monitored. However, note that the system we are considering is nonequilibrium, where the definitions of thermodynamic quantities including the temperature are precarious jou . Here I adopt the variance of $`z`$ velocity as the “temperature”, which is denoted by $`T_z`$. For comparison, I also monitor the variance of $`x`$ velocity which is denoted by $`T_x`$. Note that the variances can be different values for different directions; e.g. $`T_x`$ is smaller than $`T_z`$ in this situation. Temporal behaviors of $`T_z`$, $`T_x`$, and the average $`z`$ velocity denoted by $`v_z`$ are shown in Fig.1. Spallation occurs at $`t=5`$ where $`v_z`$ suddenly increases up to about $`2u_p=6`$ (the spall velocity). After the spallation, all of the quantities increase gradually until the ejected molecules collide with the far side of the void (approximately $`t=10`$). It should be stressed that $`v_z`$ increases up to about $`2.6u_p`$, which means that there is some kind of momentum transfer into the ejecta gas. (We see this acceleration by observing the velocity distribution functions in the next paragraph.) Also it should be remarked that the maxima of $`T_z`$ and $`v_z`$ are realized simulataneously, which is just before the collision of ejecta gas and the downstream wall of the void. It is not only the compression by the void wall that is responsible for heating of the ejecta gas but also some kind of focusing of momentum and energy. However, recall that the above results only concern with the sample molecules which constitute the proceeding part of ejecta gas. As for the whole gas, observation of the spatial distributions inside the void is needed, which we will see later in detail.

To see the above circumstances regarding the acceleration of $`v_z`$ in detail, velocity distribution function is shown in Fig. 2. Strong displacement from the Maxwellian are seen after the passage of the shock front. Note that a distribution function of $`v_x`$ is also observed to show an exponential distribution just after the spallation, which then eventually relaxes to the Maxwellian. On the other hand, we can see that the distribution function of $`v_z`$ has a widely stretched tail for the positive direction, which develops until the collision with the far side of the void. Note that the tail by far exceeds the spall velocity. It is obviously responsible for the increase of $`v_z`$, since the tail develops only towards the positive direction. The developing tail means that there is a momentum transfer into the ejecta gas, which may be attributed to the focusing. As is visualized in Fig. 3, since the shock pressure also pushes the side wall of the void to form a parabolic surface, the momenta of the collapsing wall can focus into the center line of a void, which eventually causes jetting and accelerates $`v_z`$. Note that the acceleration during the flight does not happen in a rectangular-solid void with sufficiently large longitudinal length where the parabolic surface is not formed. It is obvious that larger $`v_z`$ yields more energetic collision, which eventually leads to higher overshooting temperature inside a hot spot.

Phillips et al. phillips  also observed deviations from equilibrium distributions in a collapsing void using a two dimensional LJ crystal. However, their observation focused only on a relaxation to equilibrium, whereas the transition from equilibrium to nonequilibrium is traced here. Also, the positive tail which prevails spall velocity are not seen there, which is the original contribution of this work.

Apart from the temperature and the spall velocity, we then turn to the enhancement of chemical reactions which is our original concern regarding hot spots. Since ejected particles are strongly nonequilibrium as we have seen above, it is no longer expected that the “temperature” $`T_z`$ can be a fair index of chemistry. Although chemistry involves many intramolecular degrees of freedom, strong intermolecular collision is essential to the excitation of those degrees of freedom. Therefore it is quite natural and direct to count intermolecular collisions for the estimation of chemistry. In this study, I count the couple of particles whose distance is smaller than a certain length, $`r_c`$. Then the number of the colliding particles is normalized by an initial volume of a void. Hereafter this intensive quantity is denoted by $`N_c`$. The value of $`r_c`$ can be various according to the molecules we are to consider; in this case, I choose $`r_c=0.88`$ which corresponds to $`U(r_c)=2.5`$. Note that $`U(r_c)`$ gives an energy threshold for a reaction to occur. However, since LJ potential is adopted here, we do not have the quantitative validity for the choice of $`r_c`$. Therefore, I crosscheck the behaviors of $`N_c`$ for several values of $`r_c`$ and no qualitative difference is seen regarding the results obtained in this Letter. Typical temporal behavior of $`N_c`$ is shown in Fig. 1. Note that $`N_c`$ is zero until the void collapse, which means that strong intermolecular collisions take place only in the void. Also it should be remarked that the peaks of $`N_c`$ and of $`T_z`$ are not simultaneous; It takes a certain amount of time that $`N_c`$ becomes maximum after the peak of $`T_z`$. Indeed, this delay means that enhanced chemistry needs higher density than overheating; Namely, the delay is a relaxation time for the density to increase in the collapsing void. To confirm this idea, it is observed that the delay becomes longer for larger voids. This situation becomes clear by monitoring spatiotemporal behavior of the density profile inside the void. Fig .4 shows the spatial profiles of $`N_c`$, $`T_z`$, and the density, through which we can see in datail where and how the energetic molecular collisions and the overheating take place. Note that these two snapshots are taken at the moments when $`T_z`$ and $`N_c`$ become maxima, respectively. The maximum $`N_c`$ is realized after the whole region of the void has collapsed, while the preceeding maximum temperature takes place in the relatively dilute region on the far side of the collapsing void. It is important to notice that $`N_c`$ becomes maximum after the void has totally collapsed, whereas $`T_z`$ is maximized at the beginning of collapse. The total collapse yields high density to enhance energetic intermolecular collisions. Therefore the time lag between the maxima of $`T_z`$ and $`N_c`$ are the time needed for the total collapse of the void.

The discussions so far finish spatiotemporal behaviors of the ejecta gas and enhanced chemistry. As the final part of my results, I show how the size and the shape of voids affect chemistry. Although voids in real materials can take various forms, it may be decomposed into two essential ingredients: the longitudinal length and the normal cross section. In order to specify which component is important to the initiation of chemical reactions, I check rectangular-solid voids changing the transverse and the longitudinal lengths independently. First I show the longitudinal length dependence with the transverse length (and hence the cross section) fixed. We can see that $`N_c`$ decreases for larger longitudinal length, whereas $`T_z`$ increases. As for $`T_z`$, this result is consistent with holian  where a considerable longitudinal length was needed for the enhancement of temperature. However, the decrease of $`N_c`$ again shows that the temperature cannot be a reliable index for enhanced chemistry in this nonequilibrium situation.

On the other hand, the cross section dependence of $`T_z`$ and $`N_c`$ seems quite different from the longitudinal length dependence. The temperature $`T_z`$ is maximized when the void is cubic, whereas $`N_c`$ increases to saturate for the larger cross section. It should be stressed that the cross section of $`45\times 45`$ (approximately $`30\times 30`$ unit cells) is enough for the saturation. The maximization of $`T_z`$ with a cubic void may be a result of focusing discussed in the context of the acceleration of $`v_z`$. Also note that a larger cubic void can yield higher temperature and $`v_z`$. The saturation of $`N_c`$ for larger cross section means that the focusing in a cubic void does not seriously influence the number of energetic intermolecular collisions. Recalling that higher density (the total collapse) is also required for the increase of $`N_c`$, it is plausible that the intensity of void collapse is not affected by the focusing, whereas the foregoing part of ejecta gas is accelerated to enhance the temperature. Hence it is concluded that $`N_c`$ is rather insensitive to the shapes of voids, whereas $`T_z`$ has the apparent tendencies. In any case, nanoscale voids (approximately $`30`$ unit cells in each dimension) is enough for the enhancement of energetic intermolecular collisions.

To summarize, through the observations of spatiotemporal evolutions of temperature and density inside a void, we have seen that energetic intermolecular collisions are most likely to occur when a void is totally collapsed. This is not simultaneous with the maximization of temperature which is realized at the beginning of the collapse. This time lag is due to the difference between the optimum densities where the each quantity is maximized. The shape of void affects the overshooting temperature, but cannot seriously affect the number of energetically colliding particles per volume. Especially, the intensity of enhanced chemistry is found to saturate for nanoscale void whereas the enhancement of temperature needs much larger void. Since chemical reactions originate from intermolecular collisions, monitoring only the temperature is not enough to explore the various properties of hot spot generation, such as the void-shape dependence on the enhanced chemistry. Also, velocity distribution functions are found to deviate from the Maxwellian. These strong nonequilibrium effects must be taken into account for the analysis of spallation phenomena.

The author is grateful to H. Kaburaki and F. Shimizu for helpful discussions and encouragements.