text:
# Magnetic field topology of accreting white dwarfs

## 1. Introduction

Magnetic fields have been detected in $``$ 5 % of the 2300 known single white dwarfs with field strengths ranging from $``$ 0.1–1000 MG and peaking around 16 MG (McCook & Sion 1999, Wickramasinghe & Ferrario 2000). A similar fraction of magnetic white dwarfs has been found in accreting close binary systems (cataclysmic variables), possibly with a more restricted range of field strengths (7–230 MG, Beuermann 1998). The origin of the magnetic fields is not well understood. The magnetic fluxes, which should be conserved during stellar evolution, are similar to those of the magnetic main sequence stars such as the Ap and Bp stars which constitute $``$ 5 % of the normal main sequence stars. The decay times of the lowest multipole components are predicted to be long compared to the evolutionary ages of the white dwarfs. The magnetic field topologies, at least of isolated white dwarfs, are, therefore, likely to be relics of previous evolutionary phases. In accreting systems, the field structure in the outer layers of the white dwarf, however, may have been significantly changed if the accretion rate is high enough that accretion occurs more rapidly than ohmic diffusion (Cumming 2002).

Observational evidence suggests that the field topology of isolated and accreting white dwarfs can deviate drastically from that of a centered dipole (e.g. Wickramasinghe & Ferrario 2000, Schwope 1995). Information about the strength, orientation, and structure of the surface magnetic field of white dwarfs so far has been mainly derived from the analysis of photospheric Zeeman absorption lines and thermally broadened cyclotron harmonics from the polar regions of accreting white dwarfs. The interpretation of Zeeman and cyclotron intensity spectra alone is, however, often ambiguous and circular polarization spectra are required to investigate the magnetic field structure in detail.

## 2. Zeeman tomography

We have developed a systematic method, called Zeeman tomography, to derive the global magnetic field distribution on rotating white dwarfs (Euchner et al. 2002). At present, it utilizes a database of 46800 sets of flux and circular polarization Zeeman spectra of homogeneous magnetic white dwarf atmospheres calculated with a code developed by S. Jordan for a broad range of field strengths $`|\stackrel{}{B}|=`$ 1–400 MG, effective atmospheric temperatures $`T=`$ 8000–50000 K, and 9 field directions $`\psi `$ relative to the line of sight. All spectra are calculated for a surface gravity $`\mathrm{log}g=8`$ and use a simple limb darkening law which is independent of wavelength. The magnetic field structure is approximated by a multipole expansion of the scalar magnetic potential, using spherical harmonics with coefficients $`l=1,\mathrm{},5`$ and $`m=0`$ for the zonal and sectoral periodicity of the field. The database models are integrated for the $`\stackrel{}{B}`$ distribution of 900 elements of the visible white dwarf surface at a given rotational phase. A least-squares optimization code based on an evolutionary strategy is used to reconstruct the multipole parameters from a set of flux and polarization spectra obtained at different rotational phases of the white dwarf.

The Zeeman tomography method has been thoroughly tested with the reconstruction of various field geometries from synthetic spectra (Euchner et al. 2002). Here, we present first results of our application to real observational data.

## 3. Observations

We have obtained spin-phase resolved circular spectropolarimetry of a sample of isolated and accreting magnetic white dwarfs using FORS1 at the ESO VLT during four observing runs between May 1999 and December 2000. The instrument has been set up with a Wollaston prism and a rotatable quarter-wave retarder plate. Exposures have been alternately taken with retarder plate position angles $`\varphi =45^{}`$ and $`\varphi =+45^{}`$ in order to correct for instrumental polarization effects and linear-polarization cross-talk during the data reduction. Spectra of the target star and comparison stars in the field have been obtained simultaneously by using the multi-object spectroscopy mode of FORS. This allowed us to derive individual correction functions for the atmospheric absorption losses in the target spectra and to check the data for remnant instrumental polarization. The target spectra cover the wavelength range $`38007900`$ Å at a spectral resolution $`\lambda /\mathrm{\Delta }\lambda =440`$. A signal-to-noise ratio $`S/N100`$ has been reached for the individual flux spectra at a spin-phase resolution $`\mathrm{\Delta }\mathrm{\Phi }0.2`$. All raw data have been reduced using the context MOS of the ESO MIDAS package.

All our spectropolarimetric observations of magnetic CVs have been obtained while the systems were in a low-state of accretion. Hence, the Zeeman-split photospheric hydrogen absorption lines were clearly visible as the emission of the white dwarf was not diluted by emission from the accretion stream.

## 4. Results

### 4.1. HE 1045-0908

HE 1045-0908 is an isolated hydrogen-line white dwarf with a polar field strength $`B20`$ MG and a rotational period $`P_{\mathrm{rot}}24`$ hr (Reimers et al. 1994, Schmidt et al. 2001). Our observations cover $``$ 2 hr during which drastic variability in the flux and polarization spectra has been seen (Fig. 1). The field topology can be described by a combination of aligned dipole, quadrupole, and octopole terms with a common offset perpendicular to the magnetic axis (Fig. 2). An almost equally well fit has been obtained for a multipole expansion including the zonal components up to $`l=5`$. It is reassuring that both fits yield similar ranges for the observed field strengths and viewing angles.

### 4.2. EF Eri

Spectropolarimetric observations of EF Eri have been obtained during two epochs, in December 1999 and in November 2000 and cover in total 2 spin cycles of the white dwarf. Rotational phases have been calculated using the ephemeris given by Piirola et al. (1987), with $`\mathrm{\Phi }=0.0`$ corresponding to the narrow minimum in the infrared light curve of EF Eri, i.e. to the phase when the magnetic pole is closest to the observer. As our data show little variation in the Zeeman features modeling has been done for the phase averaged flux and polarization spectra (Fig. 3 and 4). The data are sufficiently well described by a multipole expansion including the zonal components up to $`l=5`$. We note, however, that there are still systematic residuals in our best fit solution which indicate that the field topology must be even more complex.

### 4.3. CP Tuc

CP Tuc (= AX J2315-592) has been discovered by Thomas & Reinsch (1996) who derived an upper limit of the magnetic field strength at the accreting pole $`B17`$ MG from the properties of the optically thin cyclotron emission. Our spectropolarimetric observations cover $``$0.5 spin periods of CP Tuc (Fig. 5 and 6). Rotational phases have been calculated using the ephemeris given by Ramsay et al. (1999) where $`\mathrm{\Phi }=0.0`$ is defined by the minimum of the X-ray light curve. Our best-fit solution for the field topology comprises a dipole, quadrupole plus octopole combination. Again, residual features indicate an even more complex field. The high polar field strength found here is not contradicting the much lower upper limit derived by Thomas & Reinsch (1996) for the field strength in the accretion region as the latter is probably offset from the magnetic poles.

## 5. Conclusions

We have shown that Zeeman tomography is a suitable systematic method to derive the global magnetic field distribution on rotating white dwarfs. With this method we have obtained for the first time detailed information about the range of field strengths and the field topology of a sample of isolated and accreting white dwarfs. Our results clearly demonstrate that a single value as obtained with hitherto available methods is not sufficient to quantify the field of magnetic white dwarfs. Overall our model fits are in excellent agreement with observations. Remaining differences indicate that the field topology is even more complex than described by an up to 5 component multipole expansion. A more detailed discussion of our results will be presented elsewhere (Euchner et al., in prep.).

#### Acknowledgments.

Based on observations collected at the European Southern Observatory, Chile under program numbers 63.P-0003(A), 64.P-0150(C), and 66.D-0128(A). This work was supported in part by BMBF/DLR grant 50 OR 9903 6.

## References

Beuermann, K. 1998, in High Energy Astronomy and Astrophysics, Tata Inst. of Fund. Res., 100

Cumming, A. 2002, MNRAS, 333, 589

Euchner, F., Jordan, S., Beuermann, K., Gänsicke, B.T., Hessman, F.V. 2002, A&A, 390, 633

McCook, G.P., Sion, E.M. 1999, ApJS, 121, 1

Piirola, V., Coye, G.V., Reiz, A. 1987, A&A186, 120

Ramsay, G., Potter, S.B., Buckley, D.A.H., Wheatley, P.J. 1999, MNRAS, 306, 809

Reimers, D., Jordan, S., Köhler, T., Wisotzki, L. 1994, A&A, 285, 995

Schmidt, G.D., Vennes, S., Wickramasinghe, D.T., Ferrario, L. 2001, MNRAS, 328, 203

Schwope, A.D. 1995, Rev. Mod. Astron., 8, 125

Thomas, H.-C., Reinsch, K. 1996, A&A, 315, L1

Wickramasinghe, D.T., Ferrario, L. 2000, PASP, 112, 873


format: text

text:
# Dynamics of Lieb-Liniger gases

## Abstract

It is proved that the Lieb-Liniger (LL) cusp condition implementing the delta function interaction in one-dimensional Bose gases is dynamically conserved under phase imprinting by pulses of arbitrary spatial form and the subsequent *many-body* dynamics in the thermodynamic limit is expressed approximately in terms of solutions of the time-dependent *single-particle* Schrödinger equation for a set of time-dependent orbitals evolving from an initial LL-Fermi sea. As an illustrative application, generation of gray solitons in a LL gas on a ring by a phase imprinting pulse is studied.

Recent advances and experiments on ultracold atomic vapors in atom waveguides Key ; Muller ; Dekker ; Thy ; Hinds ; Bongs ; Hansel ; Denschlag ; Schreck ; Greiner ; Arlt ; DePMcCWin99 ; GorVogLea01  and potential applicability to atom interferometry Berman  and integrated atom optics Dekker ; Sch  create a need for accurate theoretical modelling of such systems, and this is underscored by recent experimental achievement of optical atom interferometric structures of Mach-Zehnder and Michelson type Dumke . At sufficiently low temperatures and densities and high transverse frequencies $`\omega _0`$ the transverse degrees of freedom of a Bose-Einstein condensate (BEC) in an atom waveguide are frozen in the ground transverse mode and the dynamics is described by an effective 1D Hamiltonian with delta function interactions Olshanii ; PetShlWal00 . In the spatially uniform case the exact ground state was found in 1963 by Lieb and Liniger (LL) LL . The properties of this LL gas are determined by a dimensionless parameter $`\gamma =mg_{1D}/n\mathrm{}^2`$ where $`n`$ is the 1D density and the effective 1D coupling constant $`g_{1D}`$ is related to the 3D s-wave scattering length $`a`$ by Olshanii  $`g_{1\mathrm{D}}=2a\mathrm{}{}_{}{}^{2}/m\mathrm{}_0^2[1(𝒞/\sqrt{2})(a/\mathrm{}_0)]`$ with $`\mathrm{}_0=\sqrt{\mathrm{}/m\omega _0}`$ the transverse oscillator length and $`𝒞=1.4603\mathrm{}`$. If $`\gamma 1`$ the dynamics reduce to those of a 1D gas of hard core, or impenetrable, point bosons Olshanii ; PetShlWal00 , the “TG gas” Olsh01 ; Ohberg ; BMO02 ; OD02 ; GanShl02 ; AstGio02 . This is a model for which exact many-body energy eigensolutions and dynamics were found using a mapping from the Hilbert space of eigenstates of an ideal gas of spinless fermions to that of many-body eigenstates of hard core bosons map ; map2 ; 1dsho ; soliton ; breakdown ; ring1 ; ring2 ; double-X . The TG regime is only a tiny part of the LL regime, most of the latter having already been reached experimentally Bongs ; Schreck ; Greiner ; DePMcCWin99 ; GorVogLea01 , and it is desirable to extend the techniques applicable in the TG regime to the whole LL regime so as to reduce the *many-body* dynamics of free and trapped LL gases to solution of the time-dependent *single-particle* Schrödinger equation for a set of time-dependent orbitals evolving from initial orbitals comprising a LL-Fermi sea. It will be shown herein that although this cannot be done exactly for all values of $`\gamma `$ and all time intervals, it can be done in good approximation for limited time intervals for all $`\gamma `$, and exactly in some special cases.

LL cusp condition and energy eigenstates: The $`N`$-particle energy eigenstates satisfy $`\widehat{H}_0\psi _{B\alpha }=E_\alpha \psi _{B\alpha }`$ where $`\psi _{B\alpha }`$ is symmetric (Bose). LL considered a uniform system and used dimensionless units such that $`\mathrm{}=2m=1`$, yielding a Hamiltonian

$$\widehat{H}_0=\underset{j=1}{\overset{N}{}}\frac{^2}{x_j^2}+2c\underset{1j<kN}{}\delta (x_jx_k)$$
(1)
where $`2c=g_{1D}=2\gamma n`$. They used periodic boundary conditions $`\psi _{B\alpha }(\mathrm{}x_j\pm L\mathrm{})=\psi _{B\alpha }(\mathrm{}x_j\mathrm{})`$. The $`N`$-particle configuration space decomposes into $`N!`$ permutation sectors, each defined by one particular ordering of numerical values of the $`x_1,\mathrm{},x_N`$. In view of the Bose symmetry of the $`\psi _{B\alpha }`$, it is sufficient to define them only in the fundamental permutation sector $`_1`$ wherein $`x_1<x_2<\mathrm{}<x_N`$, then extending to all $`N!`$ sectors by symmetry. The cusp condition following from the $`\delta (x_jx_k)`$ interactions is LL

$$\left(\frac{}{x_j}\frac{}{x_k}\right)\psi _{B\alpha }|_{x_j=x_k+}=c\psi _{B\alpha }|_{x_j=x_k+}.$$
(2)
Define $`\widehat{B}_j=/x_{j+1}/x_jc\widehat{1}`$ where $`\widehat{1}`$ is the unit operator. Then in $`_1`$

$$\widehat{B}_j\psi _{B\alpha }|_{x_{j+1}=x_j+}=0.$$
(3)
In $`_1`$ the energy eigenstates satisfying this are

$$\psi _{B\alpha }(x_1,\mathrm{},x_N)=\frac{1}{\sqrt{N!L^N}}\underset{P}{}a_P(K_\alpha )e^{i(PK_\alpha )X}$$
(4)
and their energies are $`E_\alpha =_jk_j^2`$. Here $`K_\alpha `$ is an N-dimensional vector $`(k_1,k_2,\mathrm{},k_N)`$ with $`k_1<k_2<\mathrm{}<k_N`$, the sum runs over the $`N!`$ permutations $`P`$ of this set, and $`X`$ is the $`N`$-dimensional position vector $`X=(x_1,x_2,\mathrm{},x_N)`$, assuming $`x_1<x_2<\mathrm{}<x_N`$. There are infinitely many such choices $`k_1,\mathrm{},k_N`$, and $`\alpha `$ labels each particular choice. The coefficients are GanShl02 ; Korepin

$$a_P(K_\alpha )=ϵ_P\underset{1j<\mathrm{}N}{}\frac{\gamma n+i(Pk_{\mathrm{}}Pk_j)}{\sqrt{(\gamma n)^2+(Pk_{\mathrm{}}Pk_j)^2}}$$
(5)
where $`Pk_j`$ is the image of $`k_j`$ under $`P`$ and $`ϵ_P`$ is the parity $`\pm 1`$ of $`P`$. This differs from the $`a_P`$ of LL Eqs. (2.12) ff. only by a constant phase.

Although the symmetrical extension of $`\psi _{B\alpha }`$ to all permutation sectors is bosonic (totally symmetric under permutations of $`x_1,\mathrm{},x_N`$), it can be shown Korepin  that it is nevertheless *anti*symmetric under permutations of $`k_1,\mathrm{},k_N`$. Thus $`\psi _{B\alpha }`$ must be built from $`N`$ *singly-occupied* orbitals. The ground state $`\psi _{B0}`$ is obtained by a particular choice $`K_0`$ consisting of the smallest $`N`$ different $`k_j`$ consistent with the periodic boundary condition on the $`\psi _{B\alpha }`$ LL , a “LL-Fermi sea”.

Thermodynamic limit: The allowed $`k_j`$ were determined in the LL paper by requiring that $`\psi _{B\alpha }(x_1,\mathrm{},x_N)`$ be periodic in each $`x_j`$ with periodicity length $`L`$, leading to $`(k_{j+1}k_j)L=_{s=1}^N(\theta _{sj}\theta _{s,j+1})+2\pi n_j`$ where $`\theta _{sj}=2\mathrm{tan}^1[(k_sk_j)/\gamma n]`$ and the $`n_j`$ are any integers $`1`$, whose values distinguish between all the $`N`$-particle energy eigenstates LL ; LL2 . The ground state is given by $`n_j=1`$ for all $`j`$. It follows that the individual orbitals $`e^{ik_jx}`$ in (4) are *not* L-periodic; the $`k_j`$ are not actual particle momenta but rather nonuniformly spaced quasimomenta LL . Since our goal is reduction of the many-particle time-dependent Schrödinger equation (MBTDSE) to the single-particle time-dependent Schrödinger equation (SPTDSE), this leads to difficulties related to nonhermiticity of the single-particle Hamiltonian and an ill-posed initial value problem for the orbital evolutes $`\varphi _j(x,t)`$ on the finite interval $`\frac{L}{2}<x<\frac{L}{2}`$. These can be mitigated in the thermodynamic limit $`N\mathrm{}`$, $`L\mathrm{}`$, $`N/Ln=\text{const.}`$ wherein both the quasimomenta $`k_j`$ and the actual particle momenta $`k_\nu `$ become dense and every quasimomentum can be approximated with arbitrary accuracy by the $`k_\nu =\nu 2\pi /L`$ with $`\nu =0,\pm 1,\pm 2,\mathrm{}`$, appropriate to periodic boundary conditions. The initial orbitals $`e^{ik_\nu x}/\sqrt{L}`$ and their evolutes $`\varphi _\nu (x,t)`$ are then orthonormal on $`\frac{L}{2}<x<\frac{L}{2}`$. In this limit any thermodynamically intensive or extensive quantity expressible as a sum of contributions of “quasiorbitals” evolving from the quasimomentum plane waves $`e^{ik_jx}`$ approaches a similar expression expressed in terms of contributions of L-periodic orbitals evolving from L-periodic plane waves $`e^{ik_\nu x}`$. In doing this it is necessary to take into account the nontrivial density of states of the allowed quasimomenta. For example, the ground state energy in the thermodynamic limit is $`_{\nu =\nu _0}^{\nu _0}\beta (k_\nu )k_\nu ^2`$ where $`\beta (k_\nu )`$ is the ratio of the LL ground state density of states to its limiting value $`1/2\pi `$ in the TG limit, and $`\nu _0=Lk_{LLF}/2\pi `$ where $`k_{LLF}`$, the wave vector at the top of the LL-Fermi sea, is determined by the condition $`_{k_{LLF}}^{k_{LLF}}\beta (k)𝑑k=2\pi n`$ LL .

Response to phase imprinting: Solitons have been generated in experiments on cigar-shaped BECs by application of phase-imprinting pulses Burger ; Denschlag2 . We have recently studied soliton  the exact dynamics of such a process in the TG limit, as well as phase-imprinting pulse generated interference effects in a TG gas on a ring ring1 ; ring2  and in a double-X atom interferometer double-X . A similar LL model will be studied here. Consider a LL gas on a ring of circumference $`L`$ initially in its ground state, and represent the positions of the $`N`$ atoms by a 1D coordinate $`\frac{L}{2}<x<\frac{L}{2}`$ measured around the circumference. Suppose that this system is subjected at $`t=0`$ to a phase-imprinting pulse $`v(x,t)=\delta (t)S(x)`$. Denoting the orbitals of the LL solution by $`e^{ik_jx}`$, the orbitals just after the pulse are Rojo  $`\varphi _j(x,t=0+)=e^{iS(x)}e^{ik_jx}`$, and by (4) just after the pulse $`\psi _B(x_1,\mathrm{},x_N;t=0+)=\psi _{B\alpha }(x_1,\mathrm{},x_N)\mathrm{exp}(i_{j=1}^NS(x_j))`$. Supposing that $`x_1<x_2<\mathrm{}<x_N`$, and using the $`\widehat{B}_j`$ of Eqs. (2) and (3), one finds

$$\widehat{B}_j\psi _B(0+)=i\psi _B(0+)[S^{^{}}(x_{j+1})S^{^{}}(x_j)]+e^{iS}\widehat{B}_j\psi _{B\alpha }$$
(6)
with $`S=_jS(x_j)`$. This vanishes as $`x_{j+1}x_j+`$. It follows that the cusp condition need not be imposed as a constraint during the pulse, being exactly valid at $`t=0+`$ independently of the spatial form and strength of the pulse.

Dynamics: Consider next free propagation according to $`\widehat{H}_0`$ following such a pulse. This is generated in $`_1`$ by the kinetic energy alone, i.e., by Eq. (1) with the delta function interactions omitted. Let $`\psi _B(x_1,\mathrm{},x_N;t)`$ be a wave function differing from (4) only by replacement of the orbitals $`e^{ik_jx}`$ by the solutions $`\varphi _j(x,t)`$ of the *single*-particle TDSE $`[i/t+^2/x^2]\varphi _j(x,t)=0`$ reducing to $`\varphi _j(x,t=0+)`$ as $`t0+`$. Then $`\psi _B(x_1,\mathrm{},x_N;t)`$ satisfies the MBTDSE $`[i\mathrm{}/t\widehat{H_0}(t)]\psi _B=0`$ in the *interior* of $`_1`$ for $`t>0`$ and reduces to the previously discussed wave function $`\psi _B(x_1,\mathrm{},x_N;t=0+)`$ as $`t0+`$. However, a subtle complication enters here: Although the derivative operators in $`\widehat{B}_j`$ of Eq. (2) commute with the kinetic energy, the ratio of third to second derivative of the wave function is different from the ratio of first to zeroth derivative required by the cusp condition LL ; OD02 , so a wave function propagated as above does not in general satisfy the cusp condition for $`t>0`$. However, it is satisfied trivially for $`\gamma 0`$, and also for $`\gamma \mathrm{}`$ (TG limit). Furthermore, by temporal continuity of the solution of the TDSE it is satisfied in good approximation for a nonzero range of times $`t>0`$ for all $`\gamma `$, assuming as in (6) that $`\psi _B(\mathrm{};t=0)`$ is an energy eigenstate (4), (5). By a temporal phase mixing argument the length of this range should increase with decreasing values of $`\mathrm{\Delta }E/E`$, suggesting that constraint-free propagation as above is a reasonable approximation for $`\mathrm{\Delta }E/E1`$ where $`E`$ is the energy and $`\mathrm{\Delta }E`$ is its dispersion.

Single particle density and soliton propagation: Suppose now that the initial $`N`$-particle wave function $`\psi _B(\mathrm{};t=0+)`$ differs from the LL ground state $`\psi _{B0}`$ at most by multiplication of the LL orbitals $`e^{ik_jx}/\sqrt{L}`$ by $`e^{iS(x)}`$ as in the previously discussed case of phase imprinting. Then $`\psi _B(\mathrm{};t=0+)`$ will satisfy the LL cusp condition exactly. Suppose furthermore that $`\mathrm{\Delta }E/E1`$ where $`\mathrm{\Delta }E`$ is the energy uncertainty of $`\psi _B(\mathrm{};t=0+)`$. Then to a good approximation one has in $`_1`$

$`\psi _B(x_1,\mathrm{},x_N;t)`$
$`={\displaystyle \frac{1}{\sqrt{N!}}}{\displaystyle \underset{P}{}}a_P(K_0)\varphi _{P1}(x_1,t)\mathrm{}\varphi _{PN}(x_N,t)`$ (7)
where the $`\varphi _j(x,t)`$ are evolved by the single-particle TDSE starting from $`\varphi _j(x,0+)=e^{ik_jx}e^{iS(x)}/\sqrt{L}`$. The single particle density at time $`t`$ is $`n(x,t)=N\mathrm{}|\psi _B(x,x_2,\mathrm{},x_N;t)|^2𝑑x_2\mathrm{}𝑑x_N`$ where the integrals run over one periodicity cell $`L/2x_jL/2`$ and the thermodynamic limit is assumed so that periodic exponential orbitals are used. Upon substitution from (Dynamics of Lieb-Liniger gases) one obtains a double permutation sum $`_{PP^{}}`$. Consider first the terms $`PP^{}`$. Recall that (Dynamics of Lieb-Liniger gases) holds only in $`_1`$ and is to be extended to other sectors by Bose symmetry. Then since $`P^{}`$ differs from $`P`$ by at least one interchange, and in most cases by many, integrals of the form $`_{L/2}^{x_p}\varphi _j^{}(x_m,t)\varphi _{\mathrm{}}(x_m,t)𝑑x_m`$ and $`_{x_p}^{L/2}\varphi _j^{}(x_m,t)\varphi _{\mathrm{}}(x_m,t)𝑑x_m`$ arise, where $`j`$, $`l`$, $`m`$, and $`p`$ are all different. The phases of these integrals have denominators $`L`$ following from the normalization of the $`\varphi _j`$ and their phases oscillate as a function of $`j\mathrm{}`$, and this phase variation is compounded by oscillation of the phases of the factors $`a_P^{}(K_0)a_P^{}(K_0)`$ from Eq. (5). One therefore expects that in the thermodynamic limit the net off-diagonal contribution will be small compared with the diagonal ones and will be dropped, an approximation in the same spirit as the random phase approximation. The diagonal contributions $`P=P^{}`$ behave very differently. Since (Dynamics of Lieb-Liniger gases) holds only in $`_1`$ and is to be extended to other permutation sectors by Bose symmetry, the integrations over $`x_2,\mathrm{},x_N`$ have different limits for each different ordering of $`x,x_2,\mathrm{},x_N`$. However, summing over all permutations, collecting coefficients of $`|\varphi _1(x)|^2`$, and performing appropriate permutations of the names of integration variables, one finds that the integrations over $`x_2,\mathrm{},x_N`$ sum to a multiple of $`_{L/2}^{L/2}𝑑x_2\mathrm{}_{L/2}^{L/2}𝑑x_N|\varphi _2(x_2)|^2\mathrm{}|\varphi _N(x_N)|^2=1`$, and the same holds for the coefficients of the other $`|\varphi _j(x)|^2`$. Noting also that $`|a_P(K_0)|^2=1`$ by (5) and that there are $`(N1)!`$ permutations $`P`$ such that $`P1=j`$ for each $`j=1,\mathrm{},N`$, one finds that the final expression collapses simply to the sum of absolute squares of all $`N`$ orbitals. Recalling that approximation of the nonperiodic LL exponentials $`e^{ik_jx}`$ by the periodic exponentials $`e^{ik_\nu x}`$ in the thermodynamic limit requires inclusion of a density of states factor, one obtains

$$n(x,t)=\underset{\nu =\nu _0}{\overset{\nu _0}{}}\beta (k_\nu )|\varphi _\nu (x,t)|^2$$
(8)
where $`\nu _0`$ and $`\beta (k_\nu )`$ were defined in the previous discussion of the thermodynamic limit.

To apply this to gray soliton generation consider a ring of circumference $`L`$ with circumferential coordinate $`x`$ with $`L/2xL/2`$ and impose a phase shift $`\pi `$ on the left half $`L/2x0`$. The density of states was determined by solution of Eqs. (3.18) and (3.20) of LL LL ; Atkinson , and the TDSE propagation was effected by FFT of the initial orbitals $`e^{ik_\nu x}e^{iS(x)}`$ after which the resultant Fourier components have trivial time dependence $`e^{ik^2t}`$. The density functions showing soliton propagation are shown in Fig. 1 for the cases $`\gamma =\mathrm{}`$ (TG limit), $`\gamma =1`$ (intermediate regime), and $`\gamma =0.05`$ (Bogoliubov regime) for $`\nu _0=32`$, corresponding to a LL-Fermi sea of $`65`$ orbitals. The periodicity cell lengths $`L`$ were adjusted to keep the mean density ($`65`$ in the units used) and hence the LL-Fermi momentum fixed, using Eqs. (3.14)-(3.21) of LL LL . The shape of the density profile is quite insensitive to $`\gamma `$ apart from the indicated scaling of the periodicity length $`L`$ and observation time $`t`$, in spite of strong peaking of the density of states about low momenta as $`\gamma `$ varies from the TG regime $`\gamma 1`$ to the Bogoliubov regime $`\gamma 1`$. This insensitivity is consistent with scaling of the allowed momenta $`k_\nu =2\pi /L`$ inversely with $`L`$ and scaling of observation times $`t`$ in their propagation factors $`e^{ik_\nu ^2t}`$ directly with $`L^2`$.

One can define a LL-Fermi sea quasimomentum distribution function $`f(k)`$ as the absolute square of the Fourier transform of each orbital, weighted by the LL quasimomentum density of states and summed over all orbitals. This was evaluated with the aid of an FFT code and is shown in Fig. 2 for the case $`\gamma =1`$ before and just after the $`\pi `$ phase imprint.The perturbation promotes particles just under the LL-Fermi surface to states just above the surface, indicating that Pauli blocking is still quite effective at $`\gamma =1`$, to be expected since the density of states is still almost flat. The quasimomentum distribution is not the same as the true momentum distribution, which is peaked about $`k=0`$ even in the TG limit Olshanii  and becomes increasingly peaked as $`\gamma `$ decreases AstGio02 . Nevertheless, the energy is given exactly by $`E=_kk^2f(k)`$ both for the LL ground state LL  and the LL excited states LL2 . It is thus reasonable to approximate the energy of the phase-imprinted state by the same simple expression. Furthermore, since one starts from the ground state which has zero dispersion, the energy dispersion after the pulse should be of the same order as the mean excitation energy $`E_{after}E_{before}`$. The FFT evaluation yields $`(E_{after}E_{before})/E_{before}=0.0831`$, suggesting that the soliton density profiles of Fig. 1 are good approximations to exact profiles which are not currently calculable.

*Discussion:* The LL solitons exhibit density ripples absent from Gross-Pitaevskii solitons even for $`\gamma `$ as small as $`.05`$, being many-orbital features stemming from fermionization which is present for any nonzero $`\gamma `$, no matter how small LL ; Korepin . Even the ground state of a trapped LL gas is currently unknown, so it is not clear whether or not such ripples will be present in soliton-like dynamics of a trapped LL-gas, an important topic for future investigations. The approach used for $`\delta (t)`$ pulses could be extended to arbitrary potentials $`v(x,t)`$ by split-operator propagation of orbitals. Could it be extended to the trapped LL gas?

###### Acknowledgements.

I am grateful to Maxim Olshanii and Dimitri Gangardt for helpful comments. I thank Ewan Wright for such comments and for copies of his FFT TDSE solver and quasimomentum distribution codes. This work was begun at the Institute for Theoretical Atomic and Molecular Physics, Harvard-Smithsonian Center for Astrophysics, where it was partially supported by the National Science Foundation, and at the Benasque, Spain 2002 workshop Physics of Ultracold Dilute Atomic Gases. I thank Kate Kirby for the hospitality of ITAMP and Anthony Leggett and Fernando Sols for that of the Benasque workshop. This work has been supported at the University of Arizona by Office of Naval Research grant N00014-99-1-0806.


text:
# Signatures of hadronic cosmic rays in starbursts? High-energy photons and neutrinos from NGC253

## 1 Introduction

Very recently, the CANGAROO collaboration reported the detection of the nearby starburst galaxy NGC 253 at TeV $`\gamma `$-ray energies (Itoh et al. 2002). The purpose of this Letter is to show that it is possible for starburst galaxies, like NGC 253, to produce TeV photons and neutrinos through hadronic processes.

## 2 A plausible scenario

Collective effects of strong stellar winds and supernova explosions in star-forming regions are expected to result in particle acceleration up to multi-TeV energies, or even beyond (e.g. Montmerle 1979; Cassé & Paul 1980; Bykov & Fleishman 1992a,b; Anchordoqui et al. 1999; Bykov 2001). Recently, the region co-spatial with the Cygnus OB2 association has been detected as a TeV $`\gamma `$-ray source by the HEGRA telescope array (Aharonian et al. 2002). TeV cosmic rays (CRs) accelerated in the association might be responsible for the high-energy $`\gamma `$-ray emission through the hadronic illumination of some suitable target (Butt et al. 2003). A nearby EGRET source (3EG J2033+4118) has also a likely stellar origin (Chen et al. 1996, Romero et al. 1999, Benaglia et al. 2001).

In a starburst, where the ambient density of CRs is enhanced by a local high rate of supernova explosions and the collective effects of strong stellar winds (e.g. Bykov 2001), $`\gamma `$-rays can be produced by interactions of relativistic protons with the rich interstellar medium. This interaction produces neutral pions that quickly decay into $`\gamma `$-rays. The latter can also be produced by leptonic processes, like inverse Compton interactions with the strong FIR field, and relativistic bremsstrahlung in the ambient gas. All these processes have been modelled with a set of reasonable parameters for NGC 253 by Paglione et al. (1996). The expected total $`\gamma `$-ray flux, however, is below the EGRET sensitivity limit, in accordance with the non-detection of the galaxy in the MeV-GeV range (Bhattacharya et al. 1994, Blom et al. 1999).

A source of $`\gamma `$-rays not previously considered in the mentioned analysis is the production at the base of the strong stellar winds of early-type stars. When the relativistic particles are accelerated in the stellar wind itself, for instance through multiple shocks produced by line-driven instabilities (e.g. White & Chen 1992), the $`\gamma `$-ray flux from pion decays is limited to MeV energies. However, if the star is immersed in an external bath of relativistic particles accelerated at larger scales, e.g. in the core of the starburst, TeV protons can penetrate the base of the wind to produce TeV $`\gamma `$-rays. The injection of MeV—GeV particles into the wind of the stars will be strongly attenuated by the modulation effects of the wind itself.

The dense medium in which these $`\gamma `$-rays are produced is, nonetheless, transparent to $`\gamma `$-ray propagation. The optical depth to pair production is $`\tau n\sigma R`$, where $`n`$ is the photon number density, $`\sigma `$ is the photon-photon cross section, and $`R`$ is the distance that the photon must travel to escape. In order to compute $`n`$, we need the photon energy density at the base of the stellar wind, $`U`$. We shall assume that target photons in the wind have typical energies of $``$ 1 eV. Then, $`U=3L/(4\pi R^2c)2\times 10^{10}`$ eV cm-3, where we have taken, for order-of-magnitude estimates, the following typical values: $`R5\times 10^{14}`$ cm and $`L10^{39}`$ erg s-1. Then, $`n=U/(1\mathrm{eV})`$, and the optical depth is very small, $`\tau 2\times 10^5`$. The computation of $`\tau `$ makes use of the $`\gamma \gamma `$ (1 TeV—1eV) cross section, $`\sigma 2\times 10^{30}`$ cm2 (Lang 1999). Once TeV photons are produced and escape from the wind region, they must yet leave the galaxy to reach the Earth. The starburst has a FIR field with $`U200`$ eV cm-3 (e.g. Paglione et al. 1996), and photons have to travel distances of about 70-100 pc to leave the active region of the galaxy. The optical depth is, in this case, $`\tau =(200`$ photons cm$`{}_{}{}^{3})\times (2\times 10^{30}`$ cm$`{}_{}{}^{2})\times (100`$ pc) $``$ 1.2 $`\times 10^7`$, which is also small enough to secure that photons can escape from the galaxy.

## 3 NGC 253

NGC 253 has been described as an archetypal starburst galaxy by Rieke et al. (1980), and it has been extensively studied from radio to $`\gamma `$-rays (e.g. Beck et al. 1994, Paglione et al. 1996, Ptak 1997). The supernova rate is estimated to be as high as $`0.20.3`$ yr-1, comparable to the massive star formation rate, $`0.1`$M yr-1 (Ulvestad et al. 1999, Forbes et al. 1993). The central region of this starburst is packed with massive stars. Watson et al. (1996) have discovered four young globular clusters near the center of NGC 253; they alone can account for a mass well in excess of 1.5$`\times 10^6M_{}`$ (see also Keto et al. 1999). Assuming that the star formation rate has been continuous in the central region for the last 109 yrs, and a Salpeter IMF for 0.08-100 $`M_{}`$, Watson et al. (1996) find that the bolometric luminosity of NGC 253 is consistent with 1.5 $`\times 10^8M_{}`$ of young stars. Based on this evidence, it appears likely that there are at least tens of millions of young stars in the central region of the starburst. Many stars might be obscured by the large amount of molecular material, as in the case of the massive Galactic cluster Westerlund 1 (Clark and Negueruela 2002).

Physical, morphological, and kinematic evidence for the existence of a galactic superwind has been found for NGC 253 (e.g. McCarthy et al. 1987, Heckman et al. 1990). This superwind creates a cavity of hot ($`10^8`$ K) gas, with cooling times longer than the typical expansion time scales. As the cavity expands, a strong shock front is formed on the contact surface with the cool interstellar medium. The shock velocity can reach thousands of kilometers per second, and iron nuclei can be efficiently accelerated up to high energies $`10^{20}`$ eV by diffusive mechanisms operating at large scales. Indeed, NGC 253 has been proposed as the origin for some of the observed ultra-high energy CRs (Anchordoqui et al. 1999, 2002).

The CANGAROO collaboration did not give a precise photon spectrum for NGC 253. The reported data were only the integrated flux above 0.5 TeV (securing an 11$`\sigma `$ confidence detection) and the possible range of power law indices $`\mathrm{\Gamma }=`$(1.85, 3.75) (Itoh et al. 2002). We can summarize these results expressing the differential photon spectrum as: $`F_\gamma (E_\gamma )=B\left(E_\gamma /1\mathrm{TeV}\right)^\mathrm{\Gamma }\mathrm{cm}^2\mathrm{s}^1\mathrm{TeV}^1,`$ where $`B=(3.8\pm 1.2)\times 10^{12}`$ for $`\mathrm{\Gamma }=1.8`$, or $`B=(3.2\pm 1.0)\times 10^{12}`$ for $`\mathrm{\Gamma }=3.75`$. It is worth noticing that the direct extrapolation of the spectrum at TeV energies down to the GeV band would yield an integrated value compatible with the EGRET upper limits for NGC 253 only for power law indices $`\mathrm{\Gamma }2`$. Spectra softer than those with $`\mathrm{\Gamma }=2`$ at TeV energies should present a break at a few hundred GeV.

## 4 Hadronic origin of the TeV emission from NGC 253?

The differential $`\gamma `$-ray photon number distribution from neutral pion decays at the source is given by (e.g. Gaisser 1990) $`N_\gamma (E_\gamma )=N_p(E_\gamma )\xi _A\mathrm{\hspace{0.17em}4}\pi \sigma _{pp}(E_\gamma )N\mathrm{\hspace{0.17em}2}Z_{p\pi ^0}^{(\alpha )}/\alpha .`$ Here, $`E_\gamma `$ is the photon energy, $`N=\mathrm{\Sigma }/m_p`$ is the column density, $`m_p`$ is the proton mass, and $`\alpha `$ is the proton spectral index, such that $`N_p=K_pE^\alpha `$. In addition, $`\xi _A1.5`$ is a correction factor that takes into account possible effects of heavier nuclei, $`\sigma _{pp}(E_\gamma )35`$ mb is the hadronic interaction cross section, and $`Z_{p\pi ^0}^{(\alpha )}0.17`$ stands for the fraction of kinetic energy of the parent proton that is transferred to the neutral pion during the collision. For a photon/proton spectral index $`\alpha 2`$ we get $`N_\gamma (E_\gamma )=1.12\times 10^{25}K_pE_\gamma ^2`$ cm-3 TeV-1 (notice that 1 erg $``$ 1 TeV). Assuming typical grammages for the stellar winds of OB stars in the range $`\mathrm{\Sigma }50150`$ g cm-2 (e.g. White 1985), we get $`N_\gamma (E_\gamma )=(3.410.1)K_pE_\gamma ^2`$ cm-3 TeV-1 . The luminosity in the CANGAROO observing range is

$`L`$ $`=`$ $`4\pi R^2c{\displaystyle _{0.5\mathrm{TeV}}^{20\mathrm{T}\mathrm{e}\mathrm{V}}}N_\gamma (E_\gamma )E_\gamma 𝑑E_\gamma `$
$`=`$ $`(4.713.9)\times 10^{12}R^2K_p\mathrm{erg}\mathrm{s}^1,`$
where $`R`$ is the size of the base of the wind and $`c`$ is the speed of light (all numerical coefficients are for cgs units). The CR energy density is, on the other hand,

$$\omega _{\mathrm{CR}}=_{1\mathrm{G}\mathrm{e}\mathrm{V}}^{20\mathrm{T}\mathrm{e}\mathrm{V}}N_p(E_p)E_p𝑑E_p=9.9K_p\mathrm{erg}\mathrm{s}^1\varsigma \omega _{\mathrm{CR},},$$
(2)
where $`\varsigma `$ is the amplification factor of the CR energy density with respect to the local value, $`\omega _{\mathrm{CR},}10^{12}`$ erg cm-3. This implies that the TeV luminosity can be written as $`L=(0.471.4)R^2\varsigma `$ erg s-1. The latter corresponds to the luminosity of a single star illuminated by TeV CR protons. In the central region of the starburst, millions of stars could be illuminated at the same time. The total luminosity would then be $`L_{\mathrm{total}}=_{i=1}^N_{}(0.471.4)R_i^2\varsigma _i\mathrm{erg}\mathrm{s}^1,`$ where $`N_{}`$ is the number of contributing stars, and we have allowed for the possibility that not all of them will have the same $`R`$, nor will they be located within regions having the same CR enhancement factor.

If we now assume the following ranges for the different parameters involved: $`1\times 10^{12}<R[\mathrm{cm}]<1\times 10^{15}`$ (Lamers 1999), $`300<\varsigma <3000`$ (Suchkov et al. 1993), and $`10^6<N_{}<5\times 10^7`$ (see above), there is ample room for the total luminosity to reach values above 1039 erg s-1, i.e. the observed CANGAROO luminosity above 0.5 TeV, given the distance to NGC 253 (2.5 Mpc, de Vacoulers 1978). To show this, we have carried out Monte Carlo numerical simulations where the number of stars illuminated by the accelerated protons, the grammage, the size of the base of the wind for each star, and the ambient enhancement of CRs are independently allowed to take values within their assumed ranges. For each of these cases, the total luminosity was computed. We show the results in Fig. 1a, for a sample of 1000 trials. Only for less than 40% of the trials, the total luminosity is below 1039 erg s-1. This scenario then appears to be a plausible explanation for the TeV observations of NGC 253.

Note that multi-TeV electrons in the starburst central region will experience inverse Compton losses with the FIR fields (typical energy density $`200`$ eV cm-3) in the Klein-Nishina regime, which result in a very steep $`\gamma `$-ray spectrum (that can be approximated by power-law indices $`3.2`$ (although notice that the actual spectrum is not a pure power law, see e.g. Georganopoulos et al. 2001). The total leptonic contribution from the central starburst region at TeV energies is then negligible in comparison to the hadronic one.

MeV—GeV particles cannot easily penetrate the stellar winds, and hence these winds will glow mainly at TeV energies. The modulation of CRs by stellar winds is a complex subject that has been studied in detail only in the context of the relatively weak solar wind. In such case, three-dimensional models including diffusion and the effects of the terminal shocks have been developed (e.g. Kóta & Jokipii 1983, Jokipii et al. 1993). These models can explain fairly well the basic observed features of the nucleonic CR component in the solar system, where strong gradients in the CR density are known to exist from interplanetary spacecraft measurements (e.g. Quenby et al. 1990).333Although it is highly dependent of the solar activity and the inclination angle of the incident particles, the proton flux decreases to less than $`30`$% of its initial value at 80 A.U. for 2.5 GeV particles, and to less than $`10`$% for 250 MeV particles. In OB stars, with extremely supersonic winds and mass loss rates orders of magnitude higher than the solar values, these effects should be far more pronounced, leading to an almost complete rejection of externally injected MeV—GeV particles at the base of the winds and thus to an absence of a strong (above EGRET sensibility) $`\gamma `$-ray counterpart at these energies. The modulation effect depends on the parameter $`ϵuR_\mathrm{w}/D`$, where $`u`$ is the stellar wind velocity, $`R_\mathrm{w}`$ is the radius of the wind shell, and $`D`$ is the diffusion coefficient. $`ϵ`$ measures the ratio between the diffusive timescale of the particle to their convective timescale. The former should be shorter than the latter for a TeV particle to enter into the wind cavity. For a typical O star, $`R_\mathrm{w}>10`$ pc and $`u2500`$ km/s. Hence, the diffusion coefficient of the TeV particles should be at least $`10^{28}`$ cm2 s-1, similar to what is estimated for our Galaxy (e.g. Berezinskii et al. 1990).

## 5 Neutrino signal

If the $`\gamma `$-ray emission is explained by the decay of neutral pions, their charged counterparts must produce a neutrino flux. Following Alvarez-Muñiz & Halzen (2002), the $`\nu _\mu +\overline{\nu }_\mu `$ neutrino flux $`F_\nu (E_\nu )`$ produced by the decay of charged pions in the source can be derived from the observed $`\gamma `$-ray flux $`F{}_{\gamma }{}^{}(E_\gamma )`$ by imposing energy conservation (see also Stecker 1979):

$$_{E_\gamma ^{\mathrm{min}}}^{E_\gamma ^{\mathrm{max}}}E_\gamma F{}_{\gamma }{}^{}(E_\gamma )dE_\gamma =K_{E_\nu ^{\mathrm{min}}}^{E_\nu ^{\mathrm{max}}}E_\nu F_\nu (E_\nu )dE_\nu .$$
(3)
Here, $`E_{\gamma [\nu ]}^{\mathrm{min}}`$ ($`E_{\gamma [\nu ]}^{\mathrm{max}}`$) is the minimum (maximum) energy of the photons (neutrinos) and the pre-factor $`K=1`$.444Note that there are two muon-neutrinos (out of the pion-to-muon-to-electron decay chain) for every photon produced. Each of the latter carries, on average, an energy $`E_\pi /4`$, whereas each photon carries on average half the energy of the pion. Therefore, there is as much energy in photons as there is in neutrinos, and $`K=1`$. The maximum neutrino energy is fixed by the maximum energy of the accelerated protons ($`E_p^{\mathrm{max}}`$) which can in turn be obtained from the maximum observed $`\gamma `$-ray energy ($`E_\gamma ^{\mathrm{max}}`$) as: $`E_p^{\mathrm{max}}6E_\gamma ^{\mathrm{max}},E_\nu ^{\mathrm{max}}1/12E_p^{\mathrm{max}}`$ (Alvarez-Muñiz & Halzen 2002). The minimum neutrino energy is fixed by the threshold for pion production, $`E_p^{\mathrm{min}}=\delta \frac{(2m_p+m_\pi )^22m_p^2}{2m_p},`$ where $`\delta `$ is the Lorentz factor of the accelerator relative to the observer. Assuming a Lorentz factor of order 1, the $`\nu `$-spectrum results in $`F_\nu (E_\nu )4.2\times 10^{12}(E/\mathrm{TeV})^{1.85}`$ TeV-1 cm-2 s-1 for the harder photon index, and $`F_\nu (E_\nu )9.5\times 10^{13}(E/\mathrm{TeV})^{3.75}`$ TeV-1 cm-2 s-1 for the softer photon index.

The event rate of atmospheric $`\nu `$-background is (see Anchordoqui et al. 2002b)

$$\frac{dN}{dt}|_\mathrm{B}=A_{\mathrm{eff}}𝑑E_\nu \frac{d\mathrm{\Phi }_\mathrm{B}}{dE_\nu }P_{\nu \mu }(E_\nu )\mathrm{\Delta }\mathrm{\Omega }_{1^{}\times 1^{}},$$
(4)
where $`A_{\mathrm{eff}}`$ is the effective area of the detector, $`\mathrm{\Delta }\mathrm{\Omega }_{1^{}\times 1^{}}3\times 10^4`$ sr is an assumed bin size for the observation (ICECUBE-generation detector, Karle et al. 2002), and $`d\mathrm{\Phi }_\mathrm{B}/dE_\nu 0.2(E_\nu /\mathrm{GeV})^{3.21}`$ GeV-1 cm-2 s-1 sr-1 is the $`\nu _\mu +\overline{\nu }_\mu `$ atmospheric $`\nu `$-flux (Volkova 1980, Lipari 1993). Here, $`P_{\nu \mu }(E_\nu )1.3\times 10^6(E_\nu /\mathrm{TeV})^{0.8}`$ denotes the probability that a $`\nu `$ of energy $`E_\nu >1\mathrm{TeV}`$, on a trajectory through the detector, produces a muon (Gaisser et al. 1995). On the other hand, the $`\nu `$-signal is

$$\frac{dN}{dt}|_\mathrm{S}=A_{\mathrm{eff}}𝑑E_\nu F_\nu (E_\nu )P_{\nu \mu }(E_\nu ).$$
(5)
The signal-to-noise (S/N) ratio in the 1–20 TeV band is then $`(dN/dt|_\mathrm{S})/(\sqrt{dN/dt|_\mathrm{B}})=0.12.5`$ for one year of operation, depending on the photon spectral index. A NEMO-generation detector, with a bin size $`\mathrm{\Delta }\mathrm{\Omega }0.3^{}\times 0.3^{}`$ (Riccobene et al. 2002), would have a S/N ratio roughly three times larger for the same integration time. The effect of neutrino oscillations may reduce the flux by a factor of three at Earth (Bilenky et al. 1999). However, for some possible values of the photon spectral index, this signal could be detected in a reasonable time span. The possible starburst neutrino spectra producing the previously quoted S/N ratio are plotted, together with an upper limit for the atmospheric neutrino background with vertical incidence, in Fig. 1b.555Regrettably, the south pole location of the ICECUBE telescope would render its observation of NGC 253 nearly impossible; only neutrinos going through Earth can be identified. The smaller northern telescope ANTARES (0.1 km2, Becherini et al. 2002) would need much longer integration times to be able to distinguish NGC 253 from the background, if at all possible. Nevertheless, the analysis presented here can be directly applied to M82, a northern-sky starburst at the same distance as NGC 253.

## 6 Concluding remarks

Through hadronic interactions in their central regions, nearby starburst galaxies can generate high-energy $`\gamma `$-ray emission detectable at Earth. Proton illumination of the inner winds of massive stars in a high density CR environment may result in TeV sources without strong counterparts at MeV—GeV energies. The high number of early-type stars with strong stellar winds in this kind of galaxies might produce a strong collective effect, providing the setting for accelerating protons up to multi-TeV energies. In addition, we have found that starbursts might also be potential sources for km-scale neutrino telescopes within reasonable integration times, depending on the photon spectral index. The detection —or non-detection— of the neutrino signal may be an essential piece of evidence in the determination of the origin of the high-energy radiation.

Gamma-rays from neutral pion decays will be produced not only in the stellar winds at the core of the starburst, but also in the diffuse matter among the stars. If we model the starburst region as a disk with a radius of $`300`$ pc and a thickness of $`70`$ pc, where the average density of H atoms in the diffuse ISM is $`300`$ cm-3 (Paglione et al. 1996), we get that the expected flux above 1 TeV is $`F(E>1\mathrm{TeV})3.9\times 10^{13}`$ ph cm-2 s-1, for an average CR enhancement of $`\varsigma 1000`$ (see Torres et al. 2002 for details of calculation). This means an additional diffuse contribution of $`3\times 10^{38}`$ erg s-1 to the total TeV $`\gamma `$-ray luminosity of the starburst. In most of our models this diffuse contribution is only a small fraction of the total $`\pi ^0`$-decay $`\gamma `$-ray emission, which is dominated by interactions occurring in the stellar winds.

Finally, it is important to remind that the CANGAROO observations are consistent with an extended source, whose size might be even bigger than the galaxy itself. If this is true, it would imply that there might be another component generating TeV radiation, related with the galactic halo. This component might probably be the result of Inverse Compton scattering of CMB photons by TeV electrons accelerated at the superwind. In such case, since the interactions occur in the Thomson regime, the spectrum should be harder than that produced by hadrons in the center of the starburst.666After the submission of this paper we were informed that Itoh et al. (2003) discussed in detail this point. We thank Prof. Yanagita for calling our attention to this.

Comments by L. Anchordoqui, P. Benaglia, Y.M. Butt, C. Mauche, M.M. Kaufman-Bernadó, M. Pohl, and an anonymous referee are gratefully acknowledged. G.E.R.’s research is mainly supported by Fundación Antorchas. Additional support was provided by CONICET (under grant PIP No 0430/98) and ANPCT (PICT 03-04881). The work of D.F.T. was performed under the auspices of the U.S. DOE (NNSA), by UC’s LLNL under contract No. W-7405-Eng-48. His research was done while he was a Visiting Scientist at SISSA (Italy) and a Visiting Professor at the University of Barcelona (Spain), to which he gratefully acknowledge their hospitality.


text:
# Clausius’ Virial vs. Total Potential Energy in the dynamics of a two-component system.

## 1 Introduction

In order to understand the dynamics of a two-component system which a galaxy is, during its virial evolution, we need to explain what the trends of the Clausius’ Virial and of the total potential energy of each subsystem, together with that of the whole system, are, as soon as the inner (baryonic) component contracts inside the potential well of the outer dark matter halo. In spite of the limits of this model which, of course, is able to reproduce only some essential features of a real galaxy, its relative analytical simplicity, due to the use of tensor virial theorem, provides us with a powerful tool for understanding the role of the physical parameters involved in the dynamic evolution of a real system. For the sake of simplicity we will assume that the outer component is at fixed size and shape without considering this constraint to be too essential in order to determine the main features of the dynamic evolution we are dealing with ( see, (Secco, 2001) , hereafter LS1). Indeed the masses of the two components are not equal, the outer one being about ten times the inner one. As a consequence, tidal influences between the subsystems are not symmetric; that acting from inner to the outer is actually weaker than the reverse (Caimmi & Secco 1993). Even if a contraction effect is induced from inner density distribution on the inner regions of outer halo, during the dynamical evolution, as already underlined by Barnes & White (1984) , that effect seems, indeed, does not cause a dramatic modification of the outer mass distribution, if there is a supernovae- driven outflow, according to some recent N-body simulations (Lia et al. 2000). On the other side, models in which this constraint has been changed with some less stringent additional conditions (Caimmi 1994), seem to prove the same conclusion. Then, in order to underline the main effect of inducing a tidal scale length from the dark halo over the luminous component, we will neglect this possible modification. It seems to us that this may correspond to change from a given mass distribution to an other with a new exponent $`d`$ (see, sect.3) not too different from the previous one; that might be included in the different cases considered without cancelling the effect we are dealing with.

The two subsystems will be modelled by two penetrating similar-strata spheroids with two power-law density profiles. The advantage of this model is to be able to handle all the energetics of a two-component system during its dynamic evolution in an analytical way. Nevertheless to obtain some physical meaningful relationships, we prefer to sacrifice the more correct, but very boring, exact analytical expressions and to use some suitable mathematical approximations. Through them it will be possible to describe the whole virial evolution by means of relatively simple equations which involve only the two exponents ($`b`$ and $`d`$ of sect.3) which characterize the two subsystem power-law density distributions ($`B`$, the inner and $`D`$ the outer, respectively), as soon as the common form factor $`F`$ of sect.3 is fixed (=2, in the spherical case). The appearance of a non-monotonic trend for the inner Clausius’ virial energy, the contrary of that which refers to the inner total potential energy, is the most important result for the dynamics of a two-component system during its virial phase. Its relevance has already been stressed by demonstrating how some of the main features of elliptical galaxy Fundamental Plane may follow as a consequence of it (LS1). Its occurence seems to be a general feature when one component is completely embedded in another. This has also been tested by considering different models from those taken into account here as such the case of heterogeneous density profiles with different eccentricities (Caimmi & Secco, 2001).

Moreover, how much the extensions of the cores, at the center of the two substructures, have to be in order to obtain a maximum of the virial energy in the inner component, will be investigated when the luminous matter distribution becomes particularly steep. The present scheme may be extended to other two-component systems built up of baryonic + dark matter as, e.g., the galaxy clusters are.

## 2 Virial and total potential energy

By referring to Caimmi & Secco (1992) and the references therein, the Clausius’ virial, or virial potential energy, $`V_u`$ ($`u=B,D`$), appears with the kinetic energy, $`T_u`$, in the scalar virial equations related to the two components, the inner, $`B`$, embedded in the outer, $`D`$, of a virialized system:

$$2T_u+V_u=0.$$
(1)
It refers to the energy contribution due to the two active forces on the respective subsystem: the self-gravity and the tidal-gravity due to the other component.

Generally speaking, it turns out to be:

$$V_u=\mathrm{\Omega }_u+V_{uv},$$
(2)
where,

$$\mathrm{\Omega }_u=\rho _u\underset{r=1}{\overset{3}{}}x_r\frac{\mathrm{\Phi }_u}{x_r}d\stackrel{}{x}_u,$$
(3)
and

$$(V_{uv})=\rho _u\underset{r=1}{\overset{3}{}}x_r\frac{\mathrm{\Phi }_v}{x_r}d\stackrel{}{x}_u,$$
(4)
$`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_v`$ being the gravitational potentials due to $`u`$-matter an $`v`$-matter distributions, respectively.

In order to obtain the total potential energy of the $`u`$-component, we have to add at the self-potential energy the interaction energy $`W_{uv}`$, i.e.:

$$(E_{pot})_u=\mathrm{\Omega }_u+W_{uv},$$
(5)
The difference between the total and the virial $`u`$-energy is:

$$V_u(E_{pot})_u=Q_{uv},$$
(6)
$`Q_{uv}`$ being the residual energy, which is anti-symmetric with respect to the exchange of one component with the other. The following equations also hold:

$`V_{uv}`$ $`=`$ $`W_{uv}+Q_{uv},`$ (7)
$`W_{uv}`$ $`=`$ $`W_{vu};`$ (8)
the latter shows the symmetry of the interaction energy tensor.

As consequence, the total potential energy tensor trace of the whole system, in the principal inertia axes frame of reference, is expressed by:

$$E_{pt}=(E_{pot})_u+(E_{pot})_v=\mathrm{\Omega }_u+\mathrm{\Omega }_v+W_{uv}+W_{vu}=V_u+V_v.$$
(9)
It should be underlined that the residual energy may become zero, and then the total potential energy of each subsystem is equal to the Clausius’ virial energy, only for some specific configurations which depend on the mass distributions of the two subsystems (see, sect.7).

## 3 Similar and similar-strata spheroids

According to Caimmi (1993) (hereafter RC) and to Secco (2000) (hereafter LS), we consider, now, the case of two similar spheroids (that means, with the same axis ratio $`ϵ_B=ϵ_D=ϵ`$ (then the same form factor $`F=(2\frac{\alpha (ϵ)}{ϵ}+\frac{ϵ^2\gamma (ϵ)}{ϵ})`$), with similar and coaxial strata mass distributions. We set the radial density profiles according to two power-laws as follows:

$`\rho _B`$ $``$ $`{\displaystyle \frac{1}{r^b}},`$ (10)
$`\rho _D`$ $``$ $`{\displaystyle \frac{1}{r^d}};`$ (11)
and to avoid the central divergence we add an inner homogeneous core to both mass distributions. For the sake of simplicity we assume here that these two cores have the same extension ($`\xi _c`$) in the a-dimensional coordinates (see, RC). Then, the traces of B and D self-potential energy tensors are, respectively:

$`\mathrm{\Omega }_B`$ $`=`$ $`\nu _{\mathrm{\Omega }B}{\displaystyle \frac{GM_B^2}{a_B}}F,`$ (12)
$`\mathrm{\Omega }_D`$ $`=`$ $`\nu _{\mathrm{\Omega }D}{\displaystyle \frac{GM_D^2}{a_D}}F.`$ (13)
where $`M_u`$ is the mass and $`a_u`$ the semi-major axis of the considered subsystem. The coefficient $`\nu _{\mathrm{\Omega }u}`$ depends on the mass distribution via the functions $`h`$ and $`g_{3u}`$, from Table 1 of the Appendix A in RC.

The traces of B and D tidal-potential energy tensors are:

$`V_{BD}`$ $`=`$ $`\nu _V{\displaystyle \frac{GM_B^2}{a_B}}F,`$ (14)
$`V_{DB}`$ $`=`$ $`\nu _{VD}{\displaystyle \frac{GM_D^2}{a_D}}F.`$ (15)
The symmetric interaction energy tensor trace is

$$W_{BD}=\nu _W\frac{GM_B^2}{a_B}F.$$
(16)
The following expressions define the $`\nu `$-coefficients of the previous energies:

$`\nu _V`$ $`=`$ $`{\displaystyle \frac{9}{8}}\left[(\nu _B)_M(\nu _D)_M\right]^1mw^{ext}(x),`$ (17)
$`\nu _{VD}`$ $`=`$ $`{\displaystyle \frac{9}{8}}\left[(\nu _B)_M(\nu _D)_M\right]^1{\displaystyle \frac{1}{mx}}w^{int}(x),`$ (18)
$`\nu _W`$ $`=`$ $`{\displaystyle \frac{9}{16}}\left[(\nu _B)_M(\nu _D)_M\right]^1m[w^{int}(x)+w^{ext}(x)],`$ (19)
where $`x=a_B/a_D`$ e $`m=M_D/M_B`$. The $`(\nu _u)_M`$ coefficients are defined in Tab. 1 of RC. Due to the eqs.(17, 19) the following relation also holds:

$$W_{BD}=\frac{9}{16}[(\nu _B)_M(\nu _D)_M]^1mw^{int}(x)\frac{GM_B^2}{a_B}F+\frac{1}{2}V_{BD}$$
(20)
The definitions of $`w^{int}(x)`$ and $`w^{ext}(x)`$ and their approximations are given in the Appendix.

It should be stressed that $`w^{ext}`$ contains the expression of the $`D`$-mass fraction which has dynamical effect over the $`B`$-component at every value of $`x`$ (see Appendix).

## 4 Dynamical quantities

In order to gain some physical meaningful relationships, we will adopt here some mathematical simplifications. It should be underlined that all the approximations we consider in the following section and in the Appendix, refer to the case in which the two cores have the same a-dimensional extension $`\xi _c`$ and under the limitation that $`x>\xi _c`$. In a further paper we will relax these constraints.

Under the approximations we have considered in the Appendix, eq.(14) becomes:

$$V_{BD}\nu _V^{}G\frac{M_B\stackrel{~}{M_D}}{a_B}F$$
(21)
$`B`$-Clausius’ virial may be defined as:

$$V_BGM_B^2F\left[\nu _{\mathrm{\Omega }B}\frac{1}{a_B}\nu _V^{}\frac{M_D}{M_B}\frac{1}{a_D}(\frac{a_B}{a_D})^{2d}\right]$$
(22)
or in normalized form:

$$\stackrel{~}{V_B}=\frac{V_Ba_D}{GM_B^2F}\frac{\nu _{\mathrm{\Omega }B}}{x}\nu _V^{}\left(\frac{M_D}{M_B}\right)x^{2d};x=\frac{a_B}{a_D}$$
(23)
By using the eq.(20) and the approximation (71) ($`ϵ0`$) we obtain:

$$W_{BD}\frac{9}{4}[(\nu _B)_M(\nu _D)_M]^1\frac{GM_BM_D}{a_D}\frac{\xi _c^{b+d}}{(2d)(3b)}F\frac{1}{2(2d)}V_{BD}$$
(24)
and from eq.(7):

$$Q_{BD}\frac{9}{4}[(\nu _B)_M(\nu _D)_M]^1\frac{GM_BM_D}{a_D}\frac{\xi _c^{b+d}}{(2d)(3b)}F+\frac{52d}{2(2d)}V_{BD}$$
(25)
Finally from the definition of $`(E_{pot})_B`$ (eq.(5)) and by the approximations of eqs.(24, 21), we have:

$`(E_{pot})_B`$ $``$ $`\mathrm{\Omega }_B{\displaystyle \frac{9}{4}}[(\nu _B)_M(\nu _D)_M]^1{\displaystyle \frac{GM_BM_D}{a_D}}{\displaystyle \frac{\xi _c^{b+d}}{(2d)(3b)}}F`$ (26)
$`+{\displaystyle \frac{\nu _V^{}}{2(2d)}}G{\displaystyle \frac{M_B^2m}{a_B}}({\displaystyle \frac{a_B}{a_D}})^{3d}F`$
By normalization to the same factor of $`\stackrel{~}{V_B}`$ (eq.(23)), it follows:

$$(\stackrel{~}{E}_{pot})_B\frac{\nu _{\mathrm{\Omega }_B}}{x}\frac{9}{4}[(\nu _B)_M(\nu _D)_M]^1m\frac{\xi _c^{b+d}}{(2d)(3b)}+\frac{\nu _V^{}}{2(2d)}mx^{2d}$$
(27)
or:

$$(\stackrel{~}{E}_{pot})_B\frac{\nu _{\mathrm{\Omega }_B}}{x}\frac{1}{2}\nu _V^{}m\frac{(3d)[5(b+d)]}{(2d)(3b)}+\frac{\nu _V^{}}{2(2d)}mx^{2d}$$
(28)
where is:

$$\stackrel{~}{W}_{BD}=\stackrel{~}{W}_{DB}\frac{1}{2}\nu _V^{}m\frac{(3d)[5(b+d)]}{(2d)(3b)}+\frac{\nu _V^{}}{2(2d)}mx^{2d}$$
(29)
Moreover the residual energy of eq.(25) becomes:

$$\stackrel{~}{Q}_{BD}=\stackrel{~}{Q}_{DB}\frac{1}{2}\nu _V^{}m\frac{(3d)[5(b+d)]}{(2d)(3b)}\frac{(52d)}{2(2d)}\nu _V^{}mx^{2d}$$
(30)
The Clausius’ virial related to the $`D`$ subsystem, eq.(2), is:

$$V_D=\mathrm{\Omega }_D+V_{DB}$$

where $`\mathrm{\Omega }_D`$ and $`V_{DB}`$ correspond to the eqs.(13, 15), respectively. In the limitations under which the approximation (72) holds, the approximative tensor trace of Clausius’ virial, related to $`D`$ component, becomes:

$$V_D\nu _{\mathrm{\Omega }D}G\frac{M_D^2}{a_D}F\nu _V^{}G\frac{M_DM_B}{a_D}F\left\{\frac{(3d)[5(b+d)]}{(2d)(3b)}\frac{3d}{2d}\left(\frac{a_B}{a_D}\right)^{2d}\right\}$$
(31)
By normalization at the usual factor $`\frac{GM_B^2F}{a_D}`$, we obtain:

$$\stackrel{~}{V_D}\nu _{\mathrm{\Omega }D}m^2\nu _V^{}m\frac{3d}{2d}\left\{\frac{[5(b+d)]}{(3b)}x^{2d}\right\}$$
(32)
By using the definition of the total potential energy of $`D`$-component of eq.(5) and passing through the eq.(29), we obtain:

$$(\stackrel{~}{E}_{pot})_D\nu _{\mathrm{\Omega }D}m^2\frac{1}{2}\nu _V^{}m\frac{(3d)[5(b+d)]}{(2d)(3b)}+\frac{\nu _V^{}}{2(2d)}mx^{2d}$$
(33)
The definition of total potential energy of the whole system, eq.(9), yields the following normalized and approximate expression:

$$\stackrel{~}{E}_{pt}\frac{\nu _{\mathrm{\Omega }B}}{x}\nu _{\mathrm{\Omega }D}m^2\nu _V^{}m\frac{(3d)[5(b+d)]}{(2d)(3b)}+\frac{\nu _V^{}}{(2d)}mx^{2d}$$
(34)
As a concluding remark of this section we wish open the question of which is the weight of the approximations on the trends of the energies if they are compared with those coming from the rigorous analytical expressions. This is remarkable (see the next section) when it refers to the non-monotonic character of the luminous Clausius’ virial energy which turns out to be of extremle importance for the dynamics of the two-component system. Generally speaking, the approximations considered ($`b<2`$) are as good as $`\xi _c0`$ (e.g., the maximum score for the main physical quantities of Tab.1 is of $`0.5\%`$ if $`\xi _c=0.01`$ and $`b=1.5`$), and that will constraint us to explore what happens on this limit.

## 5 Luminous component energy trends

We will analyze the trends of the two energies, the Clausius’ virial of eq.(2) and the total potential energy of eq.(5) normalized by the factor $`\frac{FGM_B^2}{a_D}`$, which correspond to the $`B`$ component mass distribution (here assumed to be the luminous one). They are approximated by the eq.(23) and the eq.(28) respectively. The related rigorous functions are:

$$\stackrel{~}{V}_B=\frac{\nu _{\mathrm{\Omega }B}}{x}\frac{\nu _V}{x}$$
(35)
and:

$$(\stackrel{~}{E}_{pot})_B=\frac{\nu _{\mathrm{\Omega }B}}{x}\frac{\nu _W}{x}$$
(36)
Their corresponding trends are given in the Fig.1 for the cases 2) and 4) of Table 1, by assuming $`\xi _c=0.01`$ and $`m=8.5`$. From the comparison between the values of the main physical quantities of Tab.1, which are evaluated by using the rigorous trends shown in the figures, and those yielded by the approximate formulae, the conclusion is that the errors are less than $`0.5\%`$. This maximum error is reached only in the case 4).

It appears the relevant non-monotonic character of the Clausius’ virial energy (curves a, a’) in respect to the monotonic one of the total potential energy of the same component (curves b, b’).

As already shown (see, LS and LS1), when we describe a gravitational structure as in section 3, a configuration which minimizes the Clausius’ virial of inner component exists, if the outer density distribution is described by a power law with an exponent less then $`2`$. This special configuration of the embedded component is also analytically defined by the d’Alembert’s Principle of virtual works. The condition of virial equilibrium means that, if the virial energy has a minimum (in absolute value), the total kinetic energy also exhibits a minimum. This implies, in turn, that if the kinetic energy of the inner bright component is dominated by the random velocity dispersion (as in ellipticals) the corresponding macroscopic pressure support has a minimum; if, on the contrary, it is the rotational kinetic energy which dominates (as in the spirals), the rotation velocity distribution has to be that which is able to minimize the corresponding total ordered kinetic energy. The semi-major axis $`a_t`$, which characterizes this configuration, is a scale length acquired by the gravitational field of the inner component, typically without any. It is expressed, by deriving the approximated eq.(23) in respect to $`x`$:

$$a_t=\left(\frac{\nu _{\mathrm{\Omega }B}}{\nu _V^{}}\frac{1}{(2d)}\frac{M_B}{M_D}\right)^{\frac{1}{3d}}a_D,$$
(37)
The result is the same as Hoerner (1958) found for the tidal radius of a spherical star cluster embedded in the galaxy tidal potential and represents a generalization of it. This is the reason the dimension given by the eq.(37) is named as tidal radius. By assuming $`a_D=const.`$, $`M_D/M_B10`$ (on the basis of cosmological arguments) and by starting with the dynamical evolution at $`a_B=a_{Bo}=a_D`$ ($`x`$=1), we can see from eq.(23) that the ratio between the tidal-energy term ($`\stackrel{~}{V}_{BD}`$, the second one) and the self-energy term ($`\stackrel{~}{\mathrm{\Omega }}_B`$, the first one) is:

$$\zeta _1=\left(\frac{\stackrel{~}{V}_{BD}}{\stackrel{~}{\mathrm{\Omega }}_B}\right)_{x=1}\frac{\nu _V^{}}{\nu _{\mathrm{\Omega }B}}m$$
(38)
For all the realistic values of the couple ($`b`$,$`d`$) we have considered (see, Tab.1 and Fig.1 in LS) the ratio $`\zeta _1`$ is greater than one and then the tidal-energy term is dominant with respect to the self-energy one. As $`a_B`$ decreases at $`a_D=const.`$, the tidal-energy decreases, in absolute value, because of the decrease in $`D`$-mass inside the $`B`$-boundary, which is $`M_Dx^{3d}`$. Then the integral of force times the the position, $`V_{BD}`$, (21) decreases as $`x^{2d}`$, meanwhile the self-energy increases in absolute value as $`x^1`$. The balance between the self- and tidal-energy is reached at about $`a_t`$ and then, when $`x`$ becomes less than $`x_t`$ ($`=a_t/a_D)`$, the self-energy overwhelms the tidal-energy. The $`\zeta _t`$ ratio between the two energies at this minimum becomes:

$$\zeta _t\frac{1}{2d}$$
(39)
The approximated eq.(23) allows us to connect analytically the $`x_e`$ value at which the two energies, self and tidal, are equal with $`x_t`$. Indeed, it turns to be:

$$x_ex_t(2d)^{\frac{1}{3d}}$$
(40)
If $`d=1`$ the equipartition between the two energies properly occurs at $`x_t`$.

Neverthless the equation (37) is not the most general expression: the two components considered here are, indeed, always homothetic as $`a_B`$ changes. In the homogeneous ($`b=d=0`$), non similar case, when the axis ratio depends weakly on $`a_B`$, one approximative analytical expression of $`a_t`$ exists (see, eq.(53) in LS). In order to look for a more general extension, as in the case of two different eccentricities with two heterogeneous density profiles, only a first order approach has been taken into account (Caimmi & Secco, 2001), because in this case the tensor, $`(V_{uv})_{ij}(u,v=B,D;i,j=x,y,z)`$, in an analytical form, is not yet available.

It is interesting to consider where the tidal radius is located on the energies curves , in order to gain more insight into its physical meaning. $`x_t`$ (the vertical mark on curves of Fig.1,) corresponds to the beginning of the self-energy domain in the Clausius’ energy trend. Indeed, as soon as $`x<x_t`$ both the virial and the total potential energy of the luminous component begin to increase (in absolute value) very steeply as the size decreases, with approximately the same trend of a single $`B`$-component self-energy (curve c in Fig.1). At dimension greater than $`x_t`$, the total potential energy curve of $`B`$-subsystem exibits slopes smaller than the previous ones at $`x<x_t`$ (see curve b of Fig.2), by changing its concavity a little bit before $`x=1`$ (the zero of the c-curve(= $`\frac{d^2(\stackrel{~}{E}_{pot})_B}{dx^2}`$) in Fig.2). If one looks at the derivative $`\frac{d(\stackrel{~}{E}_{pot})_B}{dx}`$ , curve b of Fig.2, it is manifest that $`a_t`$ gives the scale length at which the slope of the $`B`$ total potential energy changes. Indeed, if for $`x<x_t`$ the slope is about the same of that of a single contracting component (curve b’), as soon as $`x>x_t`$ the slope changes step by step in order to stabilize itself on one, which characterizes the tidal-energy domain, at about $`x3x_t1`$. It should be noted that the flex point is included inside the interval $`x_t3x_t`$ (see Tab.1, even if in the case $`4)`$, strictly speaking, the mathematical flex point does not correspond to a physical configuration in our model because its coordinate is a little bit greater than $`1`$). For $`x>x_t`$ the rate at which the slope of the $`B`$ total potential energy (curve c) varies is smaller than that at which the slope of same component without the potential well of dark matter (curve c’) varies.

In order to understand under which conditions the flex point does appear and where it is located, we will take into account the approximation given by eq.(28). Indeed, it allows us to handle the following analytical expressions. At fixed outer component, the first and second derivative of $`(\stackrel{~}{E}_{pot})_B`$ are, respectively:

$`{\displaystyle \frac{d(\stackrel{~}{E}_{pot})_B}{dx}}{\displaystyle \frac{\nu _{\mathrm{\Omega }_B}}{x^2}}+{\displaystyle \frac{1}{2}}\nu _V^{}mx^{1d}`$ (41)
$`{\displaystyle \frac{d^2(\stackrel{~}{E}_{pot})_B}{dx^2}}2{\displaystyle \frac{\nu _{\mathrm{\Omega }_B}}{x^3}}+{\displaystyle \frac{1}{2}}\nu _V^{}m{\displaystyle \frac{1}{x^d}}(1d)`$ (42)
Then the flex point will occur at:

$$x_f\left(4\frac{2d}{1d}\right)^{\frac{1}{3d}}x_t$$
(43)
It should be noted that the flex point may exist only if $`d<1`$ and that it occurs at $`x_t`$ times a factor $`f=\left(4\frac{2d}{1d}\right)^{\frac{1}{3d}}`$ (see, Tab.1) which depends only on the exponent $`d`$ and which changes in a monotonic way from $`2`$ to $`\mathrm{}`$ as $`d`$ increases from $`0`$ to $`1`$.

We now turn back on the question, opened at the end of the previous section, related to the dependence of the non-monotonic character of the luminous Clausius’ virial energy, and of the location of its maximum, on the approximations used. As we have seen the approximated expressions of energies tend to the rigorous ones as soon as the $`\xi _c0`$. Then we have to explore what happens when the core dimensions in the two components decrease. In accordance with our previous discussion on this point (see, sect.A.2 in LS), we again conclude that, if the exponent $`b`$ is not too high ($`b<2`$) and the two power-law exponents satisfy the constraints given by (65), the maximum on the $`V_B`$ trend appears and its coordinate $`x_t`$ changes very little as $`\xi _c0`$. As soon as $`2\stackrel{<}{}b<2.5`$ the maximum also exists under the limitations (65), but its location moves towards $`x=1`$ with $`\xi _c0`$, by reaching a limit coordinate which is always less than $`1`$. When $`2.5b<3`$ and again the limitations (65) are satisfied, the maximum appears if $`\xi _c`$ is not too small. Its coordinate $`x_t`$ moves towards $`1`$ as soon as $`\xi _c`$ decreases and then it disappears as $`\xi _c0`$ (in the case of, b=2.5, which happens at $`\xi _c=0.00038)`$. These last two results are something new to add to the previous analysis. Then the following warning has to be pointed out: the non-monotonic character of $`V_B`$, may be lost<sup>1</sup><sup>1</sup>1It may be conserved from a mathemathical point of view, but the maximum corresponds to a non-physical configuration in our model, because it occurs at $`x_t>1`$. if the dimensions of the homogeneous cores become very small with $`2.5b`$, even if the conditions (65) are fulfilled. To have or not to have the divergency at the center of the two configurations, seems to become an essential feature in order to obtain or not an induced scale length on the inner component. This new message again requires, with the observations made at the beginning of the previous section, a further paper on this matter.

## 6 Thermodynamical arguments

To underline as special is the tidal radius configuration we will consider here some aspects of the thermodynamics related to the luminous component. As well kwown (Chandrasekhar 1939) when a single virialized component evolves through a quasi-static sequence of contractions, it has to lose an heat amount, $`\mathrm{\Delta }Q`$, which is one half of the variation of its self potential energy due to the contraction, in the meanwhile the other half goes to heat the structure. That means we need of a dissipation mechanism (e.g., gas cloud collisions) and of cooling processes. The consequence will be a local decreasing of the structure entropy and a corresponding increase of the universe entropy according to the $`II^o`$ thermodynamic principle (Secco 1999). Now, the question is: what happens if we take into account the contraction process of the same component ($`B`$) if it is embedded inside an other non-dissipative one ($`D`$)? In order to answer the question we have to re-consider the two fundamental equations: 1) the first thermodynamic Principle applied to $`B`$ system and ii) the variation of the virial quantities of the same component during a quasi-static transition. In the $`I^o`$ Principle equation it has to appear now the work done by both the forces on the system: the self gravity and the gravity which the dark matter distribution exerts on it. The corresponding potential energy variation for a small contraction will be $`\mathrm{\Delta }V_B`$. Combining the two equations and considering that the internal energy is the total kinetic energy of the system, $`T_B`$, (Schwarzshild 1958), we obtain:

$`\mathrm{\Delta }Q=\mathrm{\Delta }V_B/2`$ (44)
$`\mathrm{\Delta }T_B=\mathrm{\Delta }V_B/2`$ (45)
where $`\mathrm{\Delta }Q`$ is the amount of heat exchanged between the system and the universe. Again the two requests that the system settles again in virial equilibrium and obeys to the $`I^o`$ Principle yields the equipartition of the Clausius’ virial energy variation that is one half of it has to be exchanged with the Univese, the other half has to contribute to change the ”temperature” of the system which is proportional to the total kinetic energy. The result is formally the same found for the single component for which the Clausius’ virial energy coincides with the self potential energy, but the difference for the $`B`$ component inside the $`D`$ one is that the potential energy $`V_B`$ now is a non-monotonic function of $`x`$ (Fig.1). The consequence is that the differential eq.(44) becomes the following:

$$\mathrm{\Delta }\stackrel{~}{Q}=\frac{1}{2}\left(\frac{\nu _{\mathrm{\Omega }B}}{x^2}\nu _V^{}m(2d)x^{1d}\right)\mathrm{\Delta }x$$
(46)
where $`\mathrm{\Delta }\stackrel{~}{Q}`$ is the exchanged amount of heat in normalized units (see, eq.23) and $`\mathrm{\Delta }x`$ the normalized variation of the $`B`$ structure semimajor axis. Moreover, by defining the entropy variation of the $`B`$ system, in arbitrary units, as:

$$\mathrm{\Delta }\stackrel{~}{S}=\frac{1}{\stackrel{~}{T_B}}\mathrm{\Delta }\stackrel{~}{Q}$$
(47)
The eqs.(47) and (23) yield:

$$\mathrm{\Delta }\stackrel{~}{S}=\left(\frac{\frac{\nu _{\mathrm{\Omega }B}}{x^2}\nu _V^{}m(2d)x^{1d}}{\frac{\nu _{\mathrm{\Omega }B}}{x}+\nu _V^{}mx^{2d}}\right)\mathrm{\Delta }x$$
(48)
where $`\stackrel{~}{T_B}=\frac{1}{2}\stackrel{~}{V_B}`$, according to eq.(1). It is easy to see that the derivative of heat is stationary at $`x=x_t`$, and as consequence also the derivative of the $`B`$ entropy. Then we have to wait a minimum or a maximum of the entropy of the luminous component at the tidal radius configuration. Indeed, the integration of eq.(48) yields always a non-monotonic trend for the function $`\stackrel{~}{S}(x)`$ with a maximun at the tidal radius, as soon as the the maximun in Clausius’ virial energy does exist. In the special case of Tab.1, $`b=0;d=0.5`$, we obtain:

$$\stackrel{~}{S}(x)\stackrel{~}{S}(1)=ln(12539)ln\left(\frac{11339x^{5/2}+1200}{x}\right)$$
(49)
Its trend is shown in Fig.3.
The existence of this entropy maximum is able to stress the relevance of the Clausius’ energy maximum. Indeed, according to Layzer (1976), the definition of thermodynamic information is :

$$I=S_{max}S$$
(50)
where $`S_{max}`$ means the maximum value the entropy of the system may have as soon as the constraints on it, which fix the actual value of its entropy to $`S`$, are relaxed. Then, in order to increase its information a system has to decrease its entropy in respect to that of the universe. For a luminous component as $`B`$, which is embedded in an other $`D`$, a contraction decreases the entropy of $`B`$ only if its dimension is smaller or equal to $`a_t`$. It means that only from this dimension downwards the component $`B`$ may gain information during its evolution due to dissipation until to become more structured by changing its gas amount into stars (Secco 1999). Even if the argument has to make wider and deeper in a future work it seems to us it may become the ground for the interpretation of the oserved cut-off radius observed in many edge-on spiral galaxies (see, e.g., van der Kruit, 1979; Pohlen et al. 2000a, 2000b and Kregel et al., 2002). From a preliminary analysis the observed cut-off adius correlates with the $`B`$ mass as the tidal radius does (Guarise et al. 2001; Secco & Guarise 2001).

## 7 DM energy trends

By looking at Fig.4 which shows the energy trends of the $`D`$-component, the non-specularity of the virial energies behavior between the two components is manifest. Indeed, $`\stackrel{~}{V}_D`$ is monotonic (curve a), moreover $`(\stackrel{~}{E}_{pot})_D`$ increases very slowly at decreasing $`x`$.

It should be underlined that the two energies become equal as soon as the residual energies, given by the eq.(30), are zero. That occurs at the configuration characterized by the size ratio:

$$x_o\left[\frac{(3d)}{(3b)}\frac{[5(b+d)]}{(52d)}\right]^{\frac{1}{2d}}$$
(51)
If the mass distributions of the two components are the same ($`b=d`$), the equality appears only when the two configurations are coinciding ($`x_o=1`$). In the cases of different mass distributions, it may be possible to find that the virial energy becomes equal to the potential energy only at one configuration which satisfies the physical condition: $`0<x<1`$; in turn, that means :

$$(2d)(bd)<0$$

i.e., $`b<d`$ (the case of the Fig.4), if we have taken into account the dynamics of two-component system at the presence of a tidal radius ($`d<2`$). Then, the difference between the energies, at $`x=1`$, becomes:

$$(\stackrel{~}{Q}_{BD})_{x=1}=(\stackrel{~}{Q}_{DB})_{x=1}\frac{1}{2}\nu _V^{}m\frac{(bd)}{(3b)}$$
(52)
All that is true not only for the $`D`$-component but, of course, also for the $`B`$-one, due to the anti-symmetric nature of the residual terms.

## 8 The energies of the whole system

In the last step we will investigate the potential energy trend of the whole system, which means the trend of the function (34) at decreasing $`x`$. It is shown in Fig.5. Again a flex point appears. Indeed, the first and second derivative turn out to be, respectively:

$`{\displaystyle \frac{d\stackrel{~}{E}_{pt}}{dx}}{\displaystyle \frac{\nu _{\mathrm{\Omega }_B}}{x^2}}+\nu _V^{}mx^{1d}`$ (53)
$`{\displaystyle \frac{d^2\stackrel{~}{E}_{pt}}{dx^2}}2{\displaystyle \frac{\nu _{\mathrm{\Omega }_B}}{x^3}}+\nu _V^{}m{\displaystyle \frac{1}{x^d}}(1d)`$ (54)
Then the flex point will occur for:

$$x_f^{}\left(2\frac{2d}{1d}\right)^{\frac{1}{3d}}x_t=x_f/2^{1/(3d)}$$
(55)
It should be noted that again the flex point may exist only if $`d<1`$ and it is located a factor $`1/2^{\frac{1}{3d}}`$ before the flex point of $`(\stackrel{~}{E}_{pot})_B`$.

The total potential energy trend of Fig.5 for the whole system is also that of its total mechanical energy when it is in virial equilibrium. Indeed, according to eq.(9), we have:

$$E_{mec}=E_{pt}+T_B+T_D=V_B+V_D\frac{1}{2}V_B\frac{1}{2}V_D=\frac{1}{2}E_{pt}$$
(56)
The result is different for the single u-subsystem. From eqs.(5,6) we obtain:

$$(E_{mec})_u=(E_{pot})_u+T_u=\mathrm{\Omega }_u+W_u\frac{1}{2}V_u=\frac{1}{2}V_uQ_{uv}$$
(57)
The trends of the u-mechanical energies are given in Fig.6. It is difficult to look for a physical meaning related to the small increase of the mechanical energy of the $`D`$ component as $`x`$ decreases (c curve). It is probably due to the artificial constraint, we have considered, that the dark matter system has a size and a mass distribution insensitive to the contraction of the baryonic component which occurs inside. As already pointed out by Barnes & White (1984), there is an iduced effect from $`B`$ density distribution on the inner regions of dark halo which, even if non-dramatic, might change the coefficient $`\nu _{\mathrm{\Omega }D},(\nu _D)_M`$ and then $`\nu _V^{}`$.

## 9 Discussion and Conclusion

The energetics of a two-component virial system has been analized during its dynamic evolution due to the contraction of a inner (baryonic) component inside an outer of dark matter, which is assumed, for sake of simplicity, to be fixed. The results are also compared with those we obtain without the presence of the dark potential well. The analysis was done with the help of some simple, good approximate expressions for all the energies (in the cases of $`b<2`$ and with $`\xi _c=0.01`$) which, once the form factor is assigned, are only functions of the $`b`$ and $`d`$ exponents in the two power-law density distributions. The different character between the energy trends of Clausius’ virial and of the total potential energy for the baryonic (or bright component $`B`$), has been underlined and the physical main consequences has been derived for the minimum of the virial energy. The particular location of the special $`B`$ size corresponding to the minimum (tidal radius) on the total potential energy curve and the discover of one flex point on it, has stressed the character of the tidal radius as a scale length induced by the dark matter halo component on the gravitational field of the inner baryonic one. This scale length marks the beginning of the two different regimes as follows: the self-energy regime, when $`x<x_t`$, and the tidal-energy regime at about $`3x_t`$.

We claim as observable counter-parties of these theoretical reasons for the real existence of this induced scale length, at least the two followings: the first one is more direct: it is the cut-off radius observed in many edge-on spiral galaxies. Thermodynamical arguments related to the presence of an entropy maximum at this scale length and some preliminary comparison between correlations from observables and theoretical expectations seem to prove the connection. The other one is more indirect but it has an important, observable consequence: if there is this length, we are able to reproduce some of the main features of the Fundamental Plane (FP) for the elliptical galaxies (LS1). From a general point of view, the two arguments are connected by the thermodynamics in the following way: the elliptical galaxies stay on the maximum entropy stage without a significant structure evolution, the cut-off spirals, practically, begins their evolution from this stage forwards.

How the main physical results depend on the approximations used has also been considered. The main result is that: a too small homogeneous cores in both the two subsystems might cancel this induced length as soon as $`b2.5`$. Avoiding the divergency at the center of the two configurations, seems to become an essential condition in order to obtain the FP features. Here we are constrained by the restriction of handling together the two cores and then we need of a further paper on this matter in order to disentangle the probable different core-dimension play.

Nevertheless, from this point of view an appealing application to the galaxy clusters may be thought of. If we take into account the values of the exponent $`\alpha `$, used by Girardi et al.(1995) in the power-laws which give the best fit for the galaxy distributions, at a distance $`r>>R_c`$ ($`R_c`$= core radius), in the clusters, we can see that the corresponding range of our exponent $`b`$ turns out to be approximately the critical one we have found in the section (5), $`2.5b<3`$, with ratios $`R_c/R_{vir}`$ (our $`\xi _c`$) which are big enough for obtaining the maximum virial energy at $`x<1`$.

The energies of the $`D`$-subsystem were also considered. Its total potential energy is practically unaffected by the contraction of the $`B`$-subsystem whereas the virial energy slowly increases, in absolute value, as $`x`$ decreases. The conditions to obtain the equality between these two energies was with the main result that only one configuration may enjoy this property. It corresponds to the initial one, when the two substructures fill the same volume, only if the two mass distributions are the same.

Finally the trend of the potential energy of the whole system tells us that a configuration exists of minimum variation of the total potential energy (and then of the total mechanical energy), a flex point, inside its monotonic behavior. The conditions required in order to obtain either the flex on the $`B`$ total potential energy , or on the total potential energy of the whole system were investigated and the result again is a constrain only on the dark matter distribution (see the tilt of the FP). The mechanical energies either of the whole system or of both the subsystems were also taken into account.

Even if we are aware that the next step we have to do is to investigate the energetics of a two-component system with or without central divergencies and characterized by two more realistic density distributions, e.g., of Hernquist (1990) kind for the inner component and of Zhao (1996) general kind for the dark matter halo (included the NFW profile (Navarro et al., 1997)), as already considered, e.g., by Caimmi & Marmo (2002), the present results seem to be very attractive for the future work.

## Acknowledgements

We thank our friend, Prof. R. Caimmi for helpful discussions and mathematical support.

## Appendix A Appendix

According to RC we define:

$`w^{int}(x)=L(\xi _c,x)+M(\xi _c,x)`$ (58)
$`w^{ext}(x)=G(\xi _c,x)+H(\xi _c,x)`$ (59)
$`L(\xi _c,x)={\displaystyle _0^{\xi _c}}F_D{\displaystyle \frac{dF_B}{d\xi _D}}\xi _D𝑑\xi _D`$ (60)
$`M(\xi _c,x)={\displaystyle _{\xi _c}^x}F_D{\displaystyle \frac{dF_B}{d\xi _D}}\xi _D𝑑\xi _D`$ (61)
$`G(\xi _c,x)={\displaystyle _0^{\xi _c}}F_B{\displaystyle \frac{dF_D}{d\xi _D}}\xi _D𝑑\xi _D`$ (62)
$`H(\xi _c,x)={\displaystyle _{\xi _c}^x}F_B{\displaystyle \frac{dF_D}{d\xi _D}}\xi _D𝑑\xi _D`$ (63)
The two functions $`F_B(\xi _B=\frac{a_D}{a_B}\xi _D)`$ and $`F_D(\xi _D)`$ in the ranges $`(0,\xi _c)`$ and $`(\xi _c,x)`$, for different values of exponents $`b`$ and $`d`$, are given in Tab.3 of Appendix A in LS, with the assumption that $`(\xi _c<\frac{a_B}{a_D})`$.

In Table 1 of Appendix A (LS1), the corresponding functions $`L(\xi _c,x)`$ and $`M(\xi _c,x)`$ are shown with the approximate analytical expressions of $`w^{int}(x)`$ we will use in the next sect.. In Table 4 of Appendix A (LS), the corresponding functions $`G(\xi _c,x)`$ and $`H(\xi _c,x)`$ are given with the approximate analytical expressions of $`w^{ext}(x)`$ which we will also need. We take into account the approximation of $`w^{ext}`$, towards which all the exact analytic expressions of it converge, by disregarding terms of higher order in respect to $`\xi _c`$ at a different degree (see Table 4 in Appendix A of LS), which is:

$$w^{ext}(\frac{a_B}{a_D})w_{prox}^{ext}=\frac{4\xi _c^{b+d}}{(3d)[5(b+d)]}(\frac{a_B}{a_D})^{3d}$$
(64)
despite the $`b`$ and $`d`$ values, provided the following limitations hold:

$$0b<3;0d<2(b+d)<5$$
(65)
These limitations also ensure that the tidal and the self potential-energy tensor have the same sign, in all the cases we take into account. In order to highlight the $`D`$ mass fraction, $`\stackrel{~}{M}_D=M_D(\frac{a_B}{a_D})^{3d}`$, which exerts dynamic effects on $`B`$ according to Newton’s theorem (the equivalent of $`M_D^{}`$ defined for the homogeneous case, see LS), we define the factor:

$$\nu _V^{}=\frac{9}{2}[(\nu _B)_M(\nu _D)_M]^1\frac{\xi _c^{b+d}}{(3d)[5(b+d)]}$$
(66)
In this way it follows that:

$$\nu _V\nu _V^{}\frac{M_D}{M_B}(\frac{a_B}{a_D})^{3d}$$
(67)
Moreover under the limitations given by the (65), the following approximation also holds (see LS1):

$`w^{int}({\displaystyle \frac{a_B}{a_D}}){\displaystyle \frac{4\xi _c^{b+d}}{(2d)[5(b+d)]}}({\displaystyle \frac{a_B}{a_D}})^{3d}{\displaystyle \frac{4\xi _c^{b+d}}{(2d)(3b)}}({\displaystyle \frac{a_B}{a_D}})+ϵ`$ (68)
$`ϵ=4{\displaystyle \frac{\xi _c^{d+3}}{(2d)3}}({\displaystyle \frac{a_D}{a_B}})^2\delta `$ (69)
$`\delta =\left[1{\displaystyle \frac{3}{3b}}\left({\displaystyle \frac{a_B}{a_D}}\right)^b\right]`$ (70)
As we have already underlined (LS1) in case $`b=0`$, $`\delta `$ and then $`ϵ`$ of eq.(68) become $`0`$ , whatever the value of $`d`$ may be, and that the relevance of $`ϵ`$, in respect to the other terms of the same equation, increases quickly, at fixed $`d`$, as soon as $`b`$ increases. Moreover, at fixed $`b`$, the weight of $`ϵ`$ grows as $`d`$ increases, even if less than before; at fixed $`b`$ and $`d`$ the contribution of $`ϵ`$ increases as the ratio $`a_B/a_D`$ becomes smaller. If we limit ourselves to choosing $`b1.5`$ and $`d0.5`$, we may disregard $`ϵ`$ in the approximation of $`w^{int}`$ with a maximum error of $`12\%`$ when $`b=1.5`$ and $`a_B/a_D`$ decreases to $`0.4`$. That occurs if we choose $`\xi _c=0.1`$, but the error decreases as $`\xi _c`$ decreases.

It should be noted that the approximation of $`w^{int}`$ is different from that of $`w^{ext}`$, even if the core dimensions $`\xi _c`$, are the same, because of the $`b`$ value. Indeed, for $`2\stackrel{<}{}b`$ one needs to introduce in the approximated $`w^{int}`$ the term $`\delta `$ in which the ratio $`a_B/a_D`$ has $`b`$ as exponent. That does not appear in the $`w^{ext}`$ approximation.

By also taking into account the approximation of eq.(64) the following relation holds, as soon as $`ϵ`$ is negligible:

$$w^{int}\frac{4\xi _c^{b+d}}{(2d)(3b)}(\frac{a_B}{a_D})w_{prox}^{ext}\frac{(3d)}{(2d)}$$
(71)
By using eq.(71), ($`ϵ0`$), we obtain the following approximation for the coefficient of eq.(18), as a function of $`\nu _V^{}`$:

$$\nu _{VD}\nu _V^{}\frac{1}{m}\frac{3d}{2d}\left\{\frac{5(b+d)}{3b}(\frac{a_B}{a_D})^{2d}\right\}$$
(72)