text:
# Quantum vortices in systems obeying a generalized exclusion principle

## Abstract

The paper deals with a planar particle system obeying a generalized exclusion principle (EP) and governed, in the mean field approximation, by a nonlinear Schrödinger equation. We show that the EP involves a mathematically simple and physically transparent mechanism, which allows the genesis of quantum vortices in the system. We obtain in a closed form the shape of the vortices and investigate its main physical properties.

PACS numbers: 03.65.-w, 03.65.Ge, 05.45.Yv

In literature experimental evidences has been presented of formation of quantized vortex structures in different macroscopic quantum systems. Vortices have been observed both in Bose (Bose-Einstein condensate of alkali atom clouds) and in Fermi (3He-A superfluids, heavy fermion superconductors UPt3 and U0.97Th0.03Be13) systems . Actually, there exist different theoretical models admitting stationary solutions of vortex type . It is apparent that the vortex type solution is imposed by the non linearities introduced in the particular model adopted. For instance, by considering a non relativistic matter system, one can construct different models starting from Non Linear Schrödinger Equations (NLSEs) or from the gauged NLSEs in the frame of a Maxwell or Chern-Simons theory.

Principal goal of the present work is to show that there exists a mathematically simple and physically transparent mechanism, imposed by the generalized Exclusion Principle (EP), which allows the genesis of quantum vortices in planar non relativistic particle systems.

Since 1932, it was clear that the effects due to the statistics and imposed to a system of free fermions by the Pauli exclusion principle can be simulated by a repulsive potential in the coordinate space . In the same way for bosons, an attractive potential can simulate the statistical behaviour of the system. After 1940, in different works, particle systems obeying to statistics, which are different from the standard Bose-Einstein and Fermi-Dirac ones, have been considered .

Recently, it has been studied a many body quantum system obeying a generalized Exclusion-Inclusion Principle (EIP) . The peculiar property of this system is that the EIP introduces an attractive or repulsive potential in the coordinate space and can simulate the intermediate (between bosonic and fermionic) behaviour of the system. We consider some examples of real physical systems where EIP can be usefully applied. The Bose-Einstein condensation originates from an attraction of statistical nature (Bose-Einstein statistics) among the particles. In several systems the Bose-Einstein condensation is studied by means of a cubic NLSE which describes in mean field approximation an attractive interaction between two bodies. In place of the cubic and simplest interaction, other interactions can be considered as, for instance, the one introduced by EIP to simulate an attraction among the particles. Analogously, in the case of superfluids or semiconductors Fermionic systems, the repulsive interaction among the particles originated from the Fermi-Dirac statistics can be simulated by the EIP with $`\kappa <0`$. In condensed matter the motion of the couple electron-hole in a semiconductor can be described again by EIP. In fact, while electrons and holes are fermions, together can be considered excited states behaving differently from a fermion or a boson. Moreover, in nuclear physics the interaction among the fermionic valence nucleons outside the core produces pairs of correlated nucleons that can be approximated as particles with a behaviour intermediate between fermionic and bosonic ones. This nuclear state (quasideuteron state) can be viewed as a particle system that obeys to EIP.
The dynamics of the canonical quantum system obeying the EIP is governed by the following NLSE:

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}={\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }\psi +W(\rho ,𝒋)\psi +i𝒲(\rho ,𝒋)\psi +V\psi ,`$ (1)
where the real and imaginary parts of the nonlinearity are given respectively by:

$`W(\rho ,𝒋)=\kappa {\displaystyle \frac{m}{\rho }}\left({\displaystyle \frac{𝒋}{1+\kappa \rho }}\right)^2,`$ (2)
$`𝒲(\rho ,𝒋)=\kappa {\displaystyle \frac{\mathrm{}}{2\rho }}\mathbf{}\left({\displaystyle \frac{𝒋\rho }{1+\kappa \rho }}\right).`$ (3)
The free parameter $`\kappa `$ is a constant which takes into account the intensity of the exclusion-inclusion statistical effects. It is easy to verify that the system described by Eq. (1) obeys the continuity equation

$`{\displaystyle \frac{\rho }{t}}+\mathbf{}𝒋=0,`$ (4)
where $`\rho =|\psi |^2`$, while the quantum current $`𝒋`$ is given by:

$`𝒋={\displaystyle \frac{i\mathrm{}}{2m}}(1+\kappa \rho )(\psi ^{}\mathbf{}\psi \psi \mathbf{}\psi ^{}).`$ (5)
Eq. (4) assures the conservation of the particle number $`N=\rho d^Dx`$ of the system.

We discuss now briefly the origin of the model described by Eq. (1). We start by considering the following nonlinear Fokker-Planck equation :

$`{\displaystyle \frac{\rho }{t}}+\mathbf{}\left[{\displaystyle \frac{\mathbf{}S}{m}}\rho (1+\kappa \rho )+D\mathbf{}\rho \right]=0,`$ (6)
being $`\mathbf{}S/m`$ the drift velocity and $`D`$ the diffusion coefficient. In the case $`\kappa =0`$ and $`D=0`$, the above Fokker-Planck equation reduces to the well-known continuity equation for $`\rho =|\psi |^2`$ and the linear quantum mechanics can be obtained if we use the ansatz $`\psi =\rho ^{1/2}\mathrm{exp}(iS/\mathrm{})`$. In Ref.  has been considered the case $`\kappa =0,D0`$ and starting from Eq. (6) a new NLSE was obtained. Differently, starting from Eq. (6), after posing $`\kappa 0`$ and $`D=0`$, Eq. (1) can be obtained.

The introduction in Eq. (1) of the factor $`1+\kappa \rho `$ originates from the presence of the EIP and allows us to take into account many particle quantum effects. In fact, the transition probability from the site $`𝒙`$ to $`𝒙^{}`$ is defined as $`\pi (t,𝒙𝒙^{})=r(t,𝒙,𝒙^{})\rho (t,𝒙)[1+\kappa \rho (t,𝒙^{})]`$ with $`r(t,𝒙,𝒙^{})`$ the transition rate. The transition probability depends on the particle population $`\rho (t,𝒙)`$ of the starting point $`𝒙`$, and also on the population $`\rho (t,𝒙^{})`$ of the arrival point $`𝒙^{}`$. For $`\kappa 0`$ the EIP holds and the parameter $`\kappa `$ quantifies how much the particle kinetics is affected by the particle population of the arrival point. If $`\kappa >0`$ the $`\pi (t,𝒙𝒙^{})`$ introduces an inclusion principle. In fact the population at the arrival point $`𝒙^{}`$ stimulates the transition, and the transition probability increases linearly with $`\rho (t,𝒙^{})`$. In the case $`\kappa <0`$ the $`\pi (t,𝒙𝒙^{})`$ takes into account the Pauli exclusion principle. If the arrival point $`𝒙^{}`$ is empty, $`\rho (t,𝒙^{})=0`$, and the $`\pi (t,𝒙𝒙^{})`$ depends only on the population of the starting point. If the arrival site is populated $`0<\rho (t,𝒙^{})\rho _{\mathrm{max}}`$, and the transition is inhibited. The range of values the parameter $`\kappa `$ can assume is bounded by the condition that $`\pi (t,𝒙𝒙^{})`$ be real and positive, as $`r(t,𝒙,𝒙^{})`$ is. Thus we conclude that $`\kappa 1/\rho _{\mathrm{max}}`$.

The form of the nonlinearity $`𝒲(\rho ,𝒋)`$ in (1) is imposed by the continuity equation (4), while the form of $`W(\rho ,𝒋)`$ is imposed by the requirement of the canonicity of the system. The Hamiltonian density of the system (1) is given by:

$`={\displaystyle \frac{\mathrm{}^2}{2m}}\left|\mathbf{}\psi \right|^2+U_{_{\mathrm{EIP}}}+V\rho ,`$ (7)
being

$`U_{_{\mathrm{EIP}}}=\kappa \rho ^2{\displaystyle \frac{(\mathbf{}S)^2}{2m}},`$ (8)
the nonlinear potential introduced by the EIP. In Ref.  it has been shown that the system (1) admits 1-D solitons.

Following the standard procedure of the linear quantum mechanics, we define the quantum velocity $`𝒗`$ through $`𝒋=𝒗\rho `$. Taking into account Eq. (5) and writing $`\psi `$ in terms of the hydrodynamic variables $`\psi =\rho ^{1/2}\mathrm{exp}(iS/\mathrm{})`$, we have:

$`𝒗=(1+\kappa \rho ){\displaystyle \frac{\mathbf{}S}{m}}.`$ (9)
The quantum velocity can be expressed in terms of the Clebsch potentials $`S,\lambda ,\mu `$ through $`m𝒗=\mathbf{}S+\lambda \mathbf{}\mu `$. It results that the EIP imposes the choice $`\lambda =\kappa \rho `$ and $`\mu =S`$ for the Clebsch potentials. Now, we define the vorticity $`𝝎`$ of the system, through $`𝝎=\mathbf{}\times 𝒗`$ and obtain

$`𝝎={\displaystyle \frac{\kappa }{m}}\mathbf{}\rho \times \mathbf{}S.`$ (10)
From (10), it results that Eq. (1) describes a vorticose system. In the present contribution, we will consider the planar, $`D=2`$, static vortex solutions of Eq. (1) in the case $`\kappa =\xi <0`$, when EIP is reduced to an exclusion principle (EP). By introducing the plane polar coordinates $`r=\sqrt{x^2+y^2}`$ and $`\theta =\mathrm{arctan}(y/x)`$ in the plane, we search for solutions of Eq. (1) in which $`S=S(\theta )`$ and $`\rho =\rho (r)`$. The continuity equation imposes $`S=\mathrm{}n\theta `$ and consequently we write $`\psi `$ as:

$`\psi (r,\theta )=\rho (r)^{1/2}e^{in\theta }.`$ (11)
The parameter $`n`$ must be integer in order to make $`\psi `$, given by Eq. (11), a single value function. In the following, we impose $`\rho d^2x=N`$, and therefore $`\rho (\mathrm{})=0`$. Let us note that Eq. (1) is not Galilei invariant  so that traveling solutions can not be obtained by boosting static solutions of Eq.(1).

The quantum velocity $`𝒗`$ for the vortices (11) becomes:

$`𝒗={\displaystyle \frac{\mathrm{}}{m}}{\displaystyle \frac{n}{r}}(1\xi \rho )\widehat{𝒆}__\theta ,`$ (12)
where $`\widehat{𝒆}__\theta `$ is the unitary vector orthogonal to the vector $`𝒓=(x,y)`$. After integration of Eq. (12) on the circle $`\gamma _{_{\mathrm{}}}`$ with center at the vortex core and with a radius $`R\mathrm{}`$, we obtain the following relevant property:

$`m{\displaystyle _\gamma _{_{\mathrm{}}}}𝒗𝑑\mathbf{}=2\pi \mathrm{}n,`$ (13)
which justifies the name of vorticity index $`n`$. The vorticity $`𝝎=\omega \widehat{𝒆}__z`$ of the system is given by:

$`\omega ={\displaystyle \frac{2\pi \mathrm{}n}{m}}\delta ^2(𝒓)\xi {\displaystyle \frac{\mathrm{}}{m}}{\displaystyle \frac{n}{r}}{\displaystyle \frac{d\rho }{dr}},`$ (14)
and taking into account $`\rho (\mathrm{})=0`$, it is easy to verify that the total vorticity depends exclusively on the behaviour of the vortex core

$`{\displaystyle \omega d^2x}=2\pi n{\displaystyle \frac{\mathrm{}}{m}}.`$ (15)
The definition of the vorticity $`𝝎=\mathbf{}\times 𝒗`$, imposes the relation:

$`{\displaystyle \omega d^2x}={\displaystyle _\gamma _{_{\mathrm{}}}}𝒗𝑑\mathbf{},`$ (16)
which can be easily verified by comparing Eqs (13) and (15).

The angular momentum of the vortex, which is a vector orthogonal to the vortex plane, can be calculated as mean value of the operator $`\widehat{L}_z=i\mathrm{}/\theta `$ and assumes the following quantized value:

$`L_z=nN\mathrm{},`$ (17)
which does not depend on the parameter $`\xi `$.

We can calculate now the complex non linearity $`W+i𝒲`$ in Eq. (1) for the vortex given by Eq. (11). The non linearity becomes a real one and Eq. (1) reduces to:

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}={\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }\psi {\displaystyle \frac{\mathrm{}^2}{m}}{\displaystyle \frac{\xi n^2}{r^2}}\rho \psi +V\psi .`$ (18)
Eq. (18) contains a nonlinearity very close to the one of cubic NLSE and describes a canonical quantum system with a Hamiltonian density given by Eq. (7), where the nonlinear potential introduced by the EP assumes the form:

$`U_{_{\mathrm{EP}}}={\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{\xi n^2}{r^2}}\rho ^2.`$ (19)
As regard the spatial shape of the vortex $`\rho =\rho (r)`$, we insert Eq. (11) into Eq. (18), and obtain the following second order ordinary differential equation:

$`{\displaystyle \frac{1}{r\rho }}{\displaystyle \frac{d}{dr}}\left(r{\displaystyle \frac{d\rho }{dr}}\right){\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d\rho }{dr}}\right)^2`$ (20)
$`{\displaystyle \frac{2n^2}{r^2}}(12\xi \rho ){\displaystyle \frac{4m}{\mathrm{}^2}}V(r)=0.`$ (21)
In the following, we are interested to study the free vortices $`V=0`$. After introducing the dimensionless variable

$`z=2n\mathrm{log}{\displaystyle \frac{r}{r_n}},`$ (22)
where $`r_n`$ is an arbitrary constant, Eq. (21) becomes:

$`{\displaystyle \frac{2}{\rho }}{\displaystyle \frac{d^2\rho }{dz^2}}\left({\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d\rho }{dz}}\right)^2+2\xi \rho 1=0.`$ (23)
Now we consider the auxiliary function $`y(\rho )`$

$`y(\rho )=\left({\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d\rho }{dz}}\right)^2,`$ (24)
which, taking into account Eq. (23), obeys the differential equation

$`{\displaystyle \frac{dy}{d\rho }}+{\displaystyle \frac{y}{\rho }}{\displaystyle \frac{1}{\rho }}(12\xi \rho )=0.`$ (25)
The function $`y(\rho )`$ after integration of (25) assumes the form

$`y(\rho )={\displaystyle \frac{\alpha }{\rho }}+1\xi \rho ,`$ (26)
being $`\alpha `$ an integration constant. Combining Eq. (24) and Eq. (26), we arrive to the following first order ordinary differential equation for the shape of the vortex:

$`\left({\displaystyle \frac{d\rho }{dz}}\right)^2=\alpha \rho +\rho ^2(1\xi \rho ).`$ (27)
All the physical solutions of Eq. (27), whatever the value of the parameter $`\alpha IR`$, must be non singular, non negative, and normalizable.

The Hamiltonian $`H=d^2x`$ of the system can be calculated starting from (7) and the expression (19) for $`U_{_{\mathrm{EP}}}`$:

$`H={\displaystyle \frac{\mathrm{}^2n^2}{m}}{\displaystyle \left[\frac{\alpha }{2}\frac{1}{r^2}+\frac{1}{r^2}\rho (1\xi \rho )\right]d^2x}.`$ (28)
We must choose $`\alpha =0`$ in order to have a finite value of the Hamiltonian. After integration on the variable $`\theta `$, we obtain the following simple expression for the Hamiltonian

$`H=2\pi n\mathrm{}{\displaystyle _0^{\mathrm{}}}j𝑑r,`$ (29)
being $`j=(\mathrm{}n/mr)\rho (1\xi \rho )`$ the value of the current $`𝒋=j\widehat{𝒆}__\theta `$ given by Eq. (5). Finally, also the differential equation for the vortex shape assumes a very simple form

$`\left({\displaystyle \frac{d\rho }{dz}}\right)^2=\rho ^2(1\xi \rho ),`$ (30)
and can be easily integrated, obtaining the vortex profile in an explicit form

$`\rho (r)={\displaystyle \frac{4}{\xi }}\left[\left({\displaystyle \frac{r}{r_n}}\right)^n+\left({\displaystyle \frac{r_n}{r}}\right)^n\right]^2.`$ (31)
Therefore, the wave function of the vortex becomes

$`\psi (r,\theta )={\displaystyle \frac{2}{\sqrt{\xi }}}\left[\left({\displaystyle \frac{r}{r_n}}\right)^n+\left({\displaystyle \frac{r_n}{r}}\right)^n\right]^1\mathrm{exp}(in\theta ).`$ (32)
The free parameter $`r_n`$, which absorbs the integration constant related to Eq. (30), can be calculated from the normalization condition $`2\pi _0^{\mathrm{}}\rho (r)r𝑑r=N`$ and assumes the value

$`\omega ={\displaystyle \frac{2\pi \mathrm{}n}{m}}\delta ^2(𝒓){\displaystyle \frac{8n^2\mathrm{}}{mr_n^2}}\left({\displaystyle \frac{r}{r_n}}\right)^{2(n1)}\left[1\left({\displaystyle \frac{r}{r_n}}\right)^{2n}\right]\left[1+\left({\displaystyle \frac{r}{r_n}}\right)^{2n}\right]^3.`$ (33)
Figure 1 shows the plot of the dimensionless quantum velocity $`v^{}=vmr_2/\mathrm{}`$ with $`r_2=\sqrt{\xi N}/\pi `$, obtained by combining Eqs. (12) and (31), as a function of the dimensionless distance from the vortex core $`r^{}=r/r_2`$ for the vortices with $`n=2,4,6,8,10`$. We observe that $`𝒗`$ is equal to zero for $`r=r_n`$ where $`\rho `$ reaches its maximum value. Figure 2 reports the behaviour of the non singular part of the dimensionless vorticity $`\omega ^{}=\omega mr_2^2/\mathrm{}`$, given by (33) for the same vortices of figure 1. Figure 3 shows the dimensionless profile $`\xi \rho `$ as a function of $`r^{}`$ for the same vortices of the previous figures. Finally, a 3D representation of the profile $`\xi \rho `$ for the ground vortex $`n=2`$ in the 2D-dimensionless space: ($`x^{}=x/r_2,y^{}=y/r_2`$) is reported in Figure 4.

The energy $`E=H`$ of the n-vortex can be calculated easily by substituting (31) into (29) and performing the integration

$`E=|n|{\displaystyle \frac{\mathrm{}^2k^2}{2m}};k^2={\displaystyle \frac{8\pi }{3\xi }}.`$ (34)
The energy of the system results to be quantized and the energy spectrum lower bounded.

We recall that the family of the vortices (31) is related to the vortices $`\rho _{_{JP}}(r)`$ of ref.  through $`\rho (r)r^2\rho _{_{JP}}(r)`$. We remark that the vortices $`\rho _{_{JP}}(r)`$ are obtained as self-dual static solutions of a Chern-Simons model and correspond to the same energy state with energy equal to zero. On the contrary, for the vortex family (31), we have that any vortex corresponds to a different state of the system whose energy is given by Eq. (34).

We conclude by noting that very recently , it has been considered a generalization of the present model (describing neutral particles which obey the EIP) in the cases of non relativistic charged particles. In this modified model the matter field obeying the EIP is minimally coupled to a gauge field whose dynamics are described within the frame of the Chern-Simons picture. The model is a canonical one and admits self-dual static non topological vortex solutions with zero energy and linear momentum. The expressions of the main physical quantities associated to these solutions are obtained. The electric charge and the angular momentum are derived analytically, while the shape, together with the electric and magnetic fields of the vortex, are obtained numerically. This model can be considered as a continuous deformation of the Jackiw and Pi one , performed by the parameter $`\kappa `$ which takes into account the EIP.


format: text

text:
# Extraordinary Late-Time Infrared Emission of Type IIn Supernovae 1footnote 11footnote 1Based in part on data collected at Subaru Telescope, which is operated by the National Astronomy Observatory of Japan

## 1 Introduction

### 1.1 Background

On 3 January 1980, Merrill (1980) detected infrared emission from SN 1979C that was stronger toward longer wavelengths. Since then similar “infrared excesses” have been observed in a number of other supernovae including SN 1980K (Telesco et al., 1981), SN 1982E (Graham et al., 1983; Graham & Meikle, 1986), SN 1982L, SN 1982R (Graham & Meikle, 1986), SN 1985L (Elias et al., 1986), SN 1993J (Lewis et al., 1994), SN 1994Y (Garnavich, Noriega-Crespo & Moro-Martin, 1996), SN 1998S (Gerardy et al., 2000; Fassia et al., 2000) and SN 1997ab (Gerardy et al., 2000). These infrared excesses have generally been interpreted as thermal emission from dust, either forming in the ejecta (Merrill, 1980; Telesco et al., 1981; Dwek et al., 1983; Elias et al., 1986; Kozasa, Hasegawa, & Nomoto, 1989a, b, 1991; Lucy et al., 1991) or lying in a pre-existing circumstellar medium and heated by the supernova flash as an “infrared echo” (Bode & Evans, 1980; Dwek, 1983; Graham et al., 1983; Graham & Meikle, 1986; Lewis et al., 1994; Garnavich, Noriega-Crespo & Moro-Martin, 1996).

In the first scenario, dust forms in the cooling ejecta after it reaches a critical temperature. Late time thermal emission from such dust has been invoked many times to explain NIR excesses observed in nova outbursts (Gerhz, 1988), where the critical formation temperature is typically 1000–1200 K. While such a picture is perhaps the simplest explanation for the dust emission (e. g., SN 1987A: Kozasa, Hasegawa, & Nomoto 1989a, b, 1991; Danziger et al. 1991; Lucy et al. 1991), it can often be ruled out by energy considerations. Since the dust grain cooling time is short and dust is localized in the ejecta where light travel times are small, the infrared luminosity must be limited by the energy available in the supernova at the time of observation. These sources include the thermal energy stored in the dusty gas, UV/X-ray radiation from the supernova, and gamma radiation from radioactive heavy elements in the ejecta. At late times ($`t200`$ d) the luminosity of observed IR excesses becomes increasingly difficult to explain within the dust formation scenario (Dwek, 1983; Gerardy et al., 2000).

In light of this problem, a number of authors have suggested that the observed IR excesses are due to infrared dust echos. In this scenario, the emitting dust is not localized to the ejecta, but is distributed in an extended circumstellar nebula. As the supernova flash propagates outward into the nebula, the dust is heated and then radiates thermally in the infrared. Thus the late-time IR luminosity is powered by the maximum-light flash which, even after geometric dilution over a large sphere, is potentially a larger energy source than that available to dust in the SN ejecta.

Such energy considerations tended to make IR echos a more attractive explanation for the observed infrared excesses. For example, Dwek (1983) found that the IR emission from the Type IIL events SN 1979C and SN 1980K was too bright, and that there was not enough energy in these objects at late times for a dust formation scenario. They were, however, able to fit the observed NIR excesses with an echo model. Similarly, Graham & Meikle (1986) were able to use a IR echo model to reproduce the observed IR light curve of the untyped object SN 1982E.

On the other hand, Elias et al. (1986) observed IR emission in SN 1985L that became brighter before fading, in contradiction to IR echo models that predict a plateau followed by fading and only rise at very early times when the supernova is bright. Elias et al. (1986) argued that this was therefore a strong case for dust formation rather than an infrared echo. However, Emmering & Chevalier (1988) published models of echos in aspherical dust clouds which showed that IR echos from highly flattened dust distributions could exhibit more complicated rising and falling IR light curves.

Much stronger evidence for dust formation was seen in SN 1987A where increased infrared emission was accompanied by a corresponding decrease in optical emission, while at the same time emission line profiles were observed to shift toward the blue (McCray, 1993, and references therein). Both effects argue in favor of dust formed with the supernova ejecta. Similarly, in SN 1998S, the evolution of the hydrogen and helium line profiles also suggested dust formed in the ejecta as the red-shifted side of the profile steadily faded while the blue-shifted side remained nearly constant (Gerardy et al., 2000). In this case, however, the thermal emission in the NIR was too bright to have come solely from the ejecta (Fassia et al., 2000) suggesting that there may be two different dust components in SN 1998S: newly formed dust in the ejecta and brightly emitting dust in the dense CSM.

### 1.2 Recent Results On Type IIn Supernovae

The NIR excesses observed in SN 1998S and SN 1997ab have added new dimensions to the problem (Gerardy et al., 2000; Fassia et al., 2000). Both objects are classified as Type IIn events, a sub-class of Type II supernovae that exhibit multi-component “Seyfert-like” hydrogen and helium emission lines (Schlegel, 1990). This class is often associated with very slowly decaying light curves and strong X-ray emission (Filippenko, 1997, and references therein). They are thought to be interacting with a dense and perhaps clumpy circumstellar wind resulting from significant mass-loss prior to core-collapse (Chugai & Danziger, 1994).

The IR emission in SN 1997ab and SN 1998S was quite luminous ($`10^{4041}`$ erg s-1) at late epochs. For comparison, the excesses seen in SN 1979C and SN 1980K were first observed around day 200 and decayed exponentially with e-folding times $`100`$ days. Yet the strong NIR excess in SN 1997ab was observed more than three years after the explosion.

Large IR luminosities at late epochs strained both the dust formation and IR echo models, especially for SN 1997ab. Furthermore, the NIR spectra of SN 1997ab and SN 1998S were nearly identical even though SN 1997ab was approximately three times older. This seemed to imply either that these events were coincidentally observed when they would show the same spectrum, or that the evolutionary time scale of these objects was quite long. These issues, plus the Type IIn classification of SN 1997ab and SN 1998S, led Gerardy et al. (2000) to propose that the IR emission was powered by the circumstellar interaction either directly through shock heating or via UV/X-ray emission from the shocked gas. They were unable to distinguish between circumstellar dust and newly formed ejecta dust as the emission source, but Fassia et al. (2000) used K and L photometry from an earlier epoch to rule out ejecta dust as the principle NIR emission source in SN 1998S.

In this paper, we present further near-infrared observations following the evolution of SN 1998S and SN 1997ab. We also present NIR observations of three other Type IIn supernovae, SN 1995N, SN 1999Z, and SN 1999el, which also exhibit large late-time IR excesses. In §2, we discuss the observations and data reduction procedures. Near-IR spectral energy distributions are presented for all objects in §3, along with NIR spectra of SN 1995N at four epochs covering nearly five years of evolution. We then consider several emission scenarios, including interaction shock precursor, shock heated dust, simple IR echo, and interaction powered echo models. We summarize our conclusions in §4.

## 2 Observations

### 2.1 Supernova Dates and Distances

SN 1995N was discovered on 5 May 1995 by Pollas (1995) in MCG-02-38-017. Based on optical spectra Benetti, Bouchet, & Schwarz (1995) estimated that it was at least 10 months old and, following Fransson et al. (2001), we will adopt 4 July 1994 as the date of zero epoch.

SN 1998S was discovered in NGC 3877 by Z. Wan (Li, Li & Wan, 1998; Qiu, 1998) on 2 March 1998. Following Gerardy et al. (2000), we chose the date of maximum light of 18 March 1998 (Garnavich 1998, private communication) as the date of zero epoch.

The other three supernovae were discovered either near or after maximum light and we adopt the discovery date as zero epoch for these objects. SN 1997ab was discovered in an objective prism image taken on 11 April 1996 for the Hamburg Quasar Survey, (Hagen, Engels, & Reimers, 1997). The host galaxy is the faint dwarf irregular galaxy HS 0948+2018.

SN 1999Z was discovered on 8 Feb 1999 by Schwartz (1999) in UGC 5608, and SN 1999el was discovered in NGC 6951 by the BAO Supernova Survey on 20 October 1999 (Cao et al., 1999).

Approximate distances to the supernovae were calculated from the recession velocities of the host galaxies, using H$`{}_{0}{}^{}=72`$ km s-1 Mpc-1 (Freedman et al., 2001). For the three nearest objects, SN 1995N, SN1998S, and SN 1999el, recession velocities corrected for Virgo-centric infall (taken from the Lyon-Meuden Extragalactic Database222http://leda.univ-lyon1.fr/) were used. The resulting distances are 26 Mpc (1857 km s-1), 15.5 Mpc (1114 km s-1), and 24 Mpc (1705 km s-1) respectively. For the more distant objects SN 1997ab and SN 1999Z, raw observed redshifts of 3600 km s-1 (Hagen & Reimers, 1997) and 15100 km s-1 (Jha et al., 1999) were used yielding distances of 50 Mpc and 210 Mpc respectively.

### 2.2 Near-Infrared Photometry

Near-infrared J-, H-, and K-band images of SN 1998S, SN 1997ab, SN 1999Z, and SN 1999el were obtained from June 1998 to November 2000 using TIFKAM on the 2.4 m Hiltner telescope at MDM Observatory. TIFKAM is a high-throughput infrared imager and spectrograph with a $`512\times 1024`$ ALLADIN InSb detector. Dithered sets of on target images were bracketed with images of nearby empty fields for sky subtraction. Most nights were photometric and photometric solutions were achieved by observing HST/2MASS photometric standards (Persson et al., 1998) at a variety of airmasses throughout the night. The data from non-photometric nights was calibrated by matching the magnitudes of field stars to images from photometric nights.

J-, H-, and K-band images of SN 1995N were taken during the period from July 1996 to May 2001. The majority of observations were performed on the 1.2 m telescope at F. L. Whipple Observatory (FLWO) and with TIFKAM on the 2.4 m Hiltner telescope at MDM Observatory. Additional observations were obtained using the NASA Infrared Telescope Facility (IRTF), Apache Point Observatory (APO), and IRCS on the Subaru telescope. In order to avoid mixing data fluxed with different photometric systems, all SN 1995N observations were fluxed relative to several nearby field stars. A photometric solution for the field stars was obtained by observing Persson et al. (1998) standards from MDM. Thus, the SN 1995N data is on the same photometric system as the data from the other four supernovae.

### 2.3 Near-Infrared Spectroscopy

Near-infrared spectra of SN 1995N were observed from April 1996 to March 2001. H\- and K-band spectra were obtained on April 4 and 5 1996 using FSPEC on the Multiple-Mirror Telescope (MMT), and using the Infrared Cryogenic Spectrometer (CRSP) on the 4m Mayall telescope at Kitt Peak National Observatory (KPNO) on 16 June 1997 and 12 February 1998. Each night of the MMT data was reduced separately and then combined into a single spectrum. An additional K-band spectrum of SN 1995N was observed on 5 March 2001 using the Infrared Camera and Spectrograph (IRCS; Kobayashi et al. 2000) on the Subaru Telescope. A log of all the near-infrared spectroscopy is presented in Table 1.

The spectroscopic observations were taken in many short exposures. Between exposures the target was dithered across the slit to sample the chip at multiple locations and provide first-order background subtraction. Then 1D spectra were extracted using standard IRAF tasks. Wavelength calibration was achieved using OH lines in the night sky background for the MMT and KPNO spectra, and using arc-lamps for the Subaru data.

A nearby F5V star from the Bright Star Catalog (Hoffleit & Jaschek, 1982) was used as a telluric standard for the MMT and KPNO data. Subaru data was corrected for telluric absorption using nearby AV and GV stars selected from the Gemini Spectroscopic Standard Star Catalogues333http://www.us-gemini.noao.edu/sciops/instruments/niri/NIRIIndex.html. The G-dwarf was divided by a solar spectrum (Livingston & Wallace, 1991; Wallace, Hinkle, & Livingston, 1993)444NSO Kitt Peak FTS data used here were produced by NSF/NOAO to remove stellar features. This was used to correct for telluric absorption in the A dwarf spectrum. The stellar features were fit from the corrected A-star spectrum and removed from the raw A-star spectrum. The raw A-star spectrum with the stellar features removed was then used to correct the target data for telluric absorption. (See Hanson, Conti & Rieke 1996; Hanson, Rieke, & Luhman 1998 for a more detailed discussion of this technique.) The instrumental response for all the runs was removed by matching the observed spectra of the telluric standards to model spectra from the stellar atmosphere calculations of Kurucz (1994)

The absolute flux levels of the spectra were set by matching the average flux in each pass band to the NIR photometric data. For two of the epochs (days 1076 and 2435), photometric images were obtained within a few days and were used to set the flux. For the other two epochs, the NIR photometry was extrapolated using the two nearest photometric points for each band. In all cases, the pass bands of Persson et al. (1998) were assumed. We estimate the flux calibration is accurate to about $`20\%`$.

There are a few issues that we have not accounted for which could slightly alter the results of our analysis. For example, the local reddening of the supernovae is not known for many of these objects and thus no correction has been applied. We have also ignored the potential problems of exact photometric calibration of objects which have spectral energy distributions (SEDs) quite different from those of the photometric standards; in this case very red continua in the H\- and K-bands and strong line emission in J-band. In light of these issues, we have emphasized the qualitative results from our analysis and treat the quantitative results as approximations rather than exact measurements.

## 3 Analysis

Table 2 lists our infrared photometric results. The observed magnitudes have been transformed into fluxes using the zero points of Beckwith et al. (1976) and Schultz et al. (2001). Near-infrared spectral energy distributions are plotted for each of the supernovae in Figure 1. The SEDs have been corrected for galactic extinction using the measurements of Schlegel, Finkbeiner, & Davis (1998) and assuming an $`R_V=3.1`$ extinction curve. Not all of the data for SN 1995N is plotted. Two epochs (days 1375 and 1785) had very poor H-band detections and are omitted in Figure 1, as are all epochs with only a K-band detection.

The overall trend is a slow evolution toward redder and fainter near-infrared continua. In three of the objects (SN 1998S, SN 1999Z, and SN 1999el) we observe the onset of the NIR excess. For example, in SN 1998S, the K-band flux remains essentially constant from day 94 to day 355 while the J\- and H-bands fade as if the blue photospere fades to reveal a much more slowly evolving red continuum underneath. In contrast, the K-band flux increases during the blue-to-red transition in SN 1999Z suggesting that the IR continuum might still have been rising during this period. SN 1995N and SN 1997ab already exhibited strong NIR excesses when they were first observed. During the observed period, SN 1997ab slowly evolved to redder colors while SN 1995N exhibited virtually no color evolution over nearly 1800 days.

### 3.1 Single-Temperature Fits

For data showing clear infrared excesses, we fit $`HK`$ color temperatures with dust emission/absorption efficiencies $`Q_\lambda =`$ constant, $`Q_\lambda \lambda ^1`$, and $`Q_\lambda \lambda ^2`$. The results are listed in Table 3. We note that single-temperature dust emission is probably not a good description for the IR excesses, which would generally be a superposition of many dust temperatures. (It would be a valid description if the emitting dust was confined to a thin shell or very far away from the SN light source so that $`r^2`$ variations would be small.) However, since we only have NIR data it is not possible to constrain the long-wavelength dust emission spectrum. Thus, our data is really only sensitive to the hottest components of the dust emission.

The best-fit curves for $`Q_\lambda \lambda ^1`$ are shown in Figure 1. We note that the J-band flux always lies significantly above the continuum fit by the $`HK`$ color. IR dust emission is only a very small component of the J-band flux which is likely dominated by the Pa$`\beta `$ emission line. As a result, we ignored the J-band flux in the analysis of the dust emission. Also, while we have plotted the $`Q_\lambda \lambda ^1`$ curves, the $`Q_\lambda =`$ constant and $`Q_\lambda \lambda ^2`$ curves also fit the H\- and K-band spectra of SN 1995N. Without longer wavelength observations we cannot constrain the dust absorption/emission efficiency. For the rest of the analysis we will assume $`Q_\lambda \lambda ^1`$ unless otherwise noted.

The observed $`HK`$ temperature evolution of SN 1995N, SN 1997ab, and SN 1998S are plotted in the upper panels of Figure 2. SN 1995N exhibits almost no temperature evolution, cooling only from about 900 K to about 800 K between the first and last epochs. In contrast, the other two supernovae show significant color evolution. SN 1997ab appears to cool relatively linearly from 930 K to 665 K over 600 days, while SN 1998S cooled quickly from 1140 K to 900 K over 128 days (between day 227 and day 355), and then more slowly down to 810 K over the next year. Apparently the observation of nearly identical spectra/temperatures in SN 1998S and SN 1997ab noted by Gerardy et al. (2000) was merely a coincidence as these objects are clearly evolving, albeit over a longer time period than the previously observed NIR excesses in supernovae.

The total infrared fluxes, calculated by integrating the $`HK`$ temperature fits over all wavelengths, are also listed in Table 3 and plotted in the lower panels of Figure 2. These luminosities should be regarded as lower limits on the total infrared luminosities, as the near-IR is only sensitive to the hottest dust and consequently a significant flux from much colder dust could be missed in our data.

Although SN 1995N showed little temperature evolution, there is significant evolution of its infrared luminosity. For the first four epochs, the luminosity seems relatively constant at about $`7\times 10^{41}`$ erg s-1. Then after day 1000 the luminosity appears to decline exponentially with an e-folding time $`600`$ d. While no plateau phase was observed in SN 1997ab, the infrared luminosity also appears to decay exponentially, in this case with an e-folding time of $`450`$ d. However, the NIR behavior of SN 1998S was somewhat different. The inferred IR luminosity of SN 1998S increased $`50\%`$ between day 260 and day 355, and then faded between day 355 and day 796 at a rate similar to that of SN 1997ab.

While we interpret the IR continuum as thermal emission from dust, it is important to estimate the extent that free-free emission might contribute to the NIR continuum as there are strong recombination lines observed in these supernovae. Fransson et al. (2001) found that in SN 1995N, H$`\alpha `$ was $`8`$ times stronger than expected for standard H I recombination due to collisional excitation and a high optical depth for H$`\alpha `$. However, the higher order Balmer lines were consistent with Case B recombination. Assuming standard nebular conditions and Case B recombination, the H$`\beta `$ emission per unit volume is $`4\pi j_{H\beta }1.2\times 10^{25}n^2`$ erg cm-3 s-1 (Osterbrock, 1989). Using an abundance ratio of He to H of 0.35 by number (Fransson et al., 2001) the free-free emission per unit volume $`ϵ_{ff}3.2\times 10^{27}T^{1/2}n^2`$ erg cm-3 sec-1 (Rybicki & Lightman, 1979), where we have assumed a velocity averaged Gaunt factor $`\overline{g}_B1.2`$. For $`T=10^4`$ K, this implies that $`ϵ_{ff}3ϵ_{H\beta }`$. Fransson et al. (2001) report an H$`\beta `$ flux of $`8.8\times 10^{15}`$ erg s-1 cm-2 from SN 1995N on day 718, which translates to a free-free luminosity of $`L_{ff}2\times 10^{39}`$ erg s-1. Compared to the IR luminosity on day 730 of $`8\times 10^{41}`$ erg s-1, the free-free continuum is down by two orders of magnitude and thus does not contribute significantly to the observed NIR continuum.

### 3.2 Comparison to X-Ray Emission

Comparison of the observed X-ray and infrared light curves can be used to provide clues into the origin of the IR emission. Both SN 1995N (Fox et al., 2000) and SN 1998S (Pooley et al., 2001) have been detected in X-rays and the observed X-ray luminosities are plotted as stars in Figure 2.

One possible dust heating scenario is that the IR emission is from a precursor region in the circumstellar medium (CSM) just ahead of the shock front that is being heated by X-rays from the post shock gas. In such a situation, where the light-travel time is small, we might expect the X-ray luminosity and the IR luminosity to be well correlated. What we find, however, is that the IR evolution does not follow the X-ray light curve in SN 1995N, as the X-ray flux shows no signs of fading nearly 200 days after the IR begins to fall. Furthermore, the IR flux is significantly greater (as much as a factor of five) than the X-rays. This would require that the actual emitted X-ray flux is much larger and that nearly all of the X-rays are being absorbed by the circumstellar dust and re-emitted in the IR. However, Fransson et al. (2001) found no indication of strong reddening in the optical spectra of SN 1995N, suggesting that the dust extinction in the CSM is low. These problems argue against an IR precursor interpretation.

We note that large IR/X-ray flux ratios are seen in supernova remnants with observed ratios varying from $`2`$ to $`100`$ (Dwek et al., 1987). In remnants the IR emission is from hot, dusty, shocked gas where thermal emission from collisionally heated dust can be the dominant cooling mechanism. However, the IR emission from supernova remnants is typically much cooler and it seems unlikely that a significant amount of dust could be collisionally heated up to 800 – 1000 K without destroying the dust grains. Collisional heating of dust to $`1000`$ K would also require prohibitively large ion densities ($`>10^{21}`$ cm-3 for $`T_e=10^7`$ K; N. Evans 2002, private communication). Furthermore, while there could be some time lag between changes in the X-ray flux and changes in the IR flux, we would not expect to see the luminosity drop in the IR before it does in the X-rays as seen in SN 1995N.

### 3.3 Infrared Echo Models

On the other hand, the shape of the SN 1995N IR light curve is similar to the luminosity evolution of simple infrared echo models. A number of authors have presented models describing the infrared echos of supernovae in spherically symmetric dust distributions (Bode & Evans, 1980; Dwek, 1983; Graham et al., 1983; Graham & Meikle, 1986). In this paper, we shall adopt the formalism of Emmering & Chevalier (1988, EC88 hereafter) which adds the generality of non-spherical flattened dust distributions. A common feature of spherical models is that they predict an infrared light curve that rises quickly (only while the supernova is at its brightest) and then levels off in a plateau before finally decaying.

The light curve plateau results because the supernova flash clears out a large dust-free cavity in the center of the circumstellar nebula by vaporizing all the dust. At the edge of this cavity is the hottest and brightest dust emission. Outside the cavity, the dust is cooler and fainter as the supernova flash is spread out over a larger surface. For a pulse (delta function) supernova light curve that explodes at time $`t^{}`$, the observed flux at any time $`t`$ will come from dust lying on a paraboloid described by $`t=t^{}+\frac{r}{c}(1\mathrm{cos}\theta )`$. As time proceeds, the paraboloids become wider and the vertex recedes from the supernova. At the time $`t_v=2r_v/c`$, the vertex of the paraboloid will reach the edge of the edge of the cavity. Before this time, the IR luminosity, dominated by emission from the cavity’s edge, will be nearly constant forming the plateau in the IR light curve. After the paraboloid leaves the cavity, the emission will come from the fainter and cooler dust and thus the IR luminosity will fade. A plateau is also expected in models with more realistic SN light curves as long as the SN fades in significantly less time than it takes for the echo paraboloid to reach the cavity wall.

The observed IR light curve for SN 1995N is at least qualitatively similar to the light curves of spherically symmetric infrared echos. However, the evolution timescale for SN 1995N is much longer than for the models that were used to fit infrared echos to previous IR excesses in supernovae. In SN 1979C, SN 1980K and SN 1982E, the IR excesses were detected a couple hundred days after maximum light, and the IR light curve was already declining. The time $`t_v`$ at which the re-emission paraboloids reached the cavity wall in the model fits was typically $``$ 50 – 100 days. For SN 1995N this would have to be more like 1000 days where the observed IR light curve begins to decay rapidly. This time scale is set by the radius of the dust-free cavity, which scales as the square root of the supernova luminosity. Thus, the much longer time scale would imply either that the dust around SN 1995N is much easier to vaporize or that SN 1995N is much brighter ($`100\times `$) than those other objects. A more likely solution might be that some event prior to core collapse cleared out the cavity.

#### 3.3.1 Numerical Calculations

To test the infrared echo scenario, we have calculated the total IR luminosity evolution and near-IR spectra for echo models based on EC88. These models assume an axis-symmetric dust distribution where the isodensity curves are ellipsoids of revolution where the dust distribution has a $`b^2`$ radial dependence. The circumstellar dust is assumed to have an emission/absorption efficiency with a power-law wavelength dependence $`Q_\lambda \lambda ^n`$. The model dust cloud has an inner dust-free zone, as discussed above, and has no sharp outer edge. The models implicitly assume that the optical depth of the dust cloud is low and do not account for absorption inside the circumstellar dust nebula, either between the SN and the heated dust or after re-emission. Our calculations for the total infrared luminosity $`F_{\mathrm{IR}}(t)`$, and the infrared spectrum $`F_\lambda (t)`$ were performed by numerically integrating Equations 8 and 15 of EC88.

The IR dust echo models have a number of parameters that can be varied. These include the dust emission/absorption efficiency exponent $`n`$, the size of the inner cavity $`r_v`$, the dust temperature at the inner surface of the cavity $`T_v`$, the SN light curve, the degree to which the dust distribution is flattened, and for non-spherical dust distributions the angle between the axis-of-symmetry and the line of sight. Since the observed plateau in the IR light curve was suggestive of a spherically symmetric dust distribution, we used spherically symmetric echo models for SN 1995N. This eliminates two of the model parameters and also makes the numerical calculations easier.

In principle, if the dust-free cavity were created by the supernova flash vaporizing the innermost dust, the cavity radius $`r_v`$, the peak luminosity of the supernova, and the dust temperature at the inner edge of the cavity $`T_v`$ would all be related. If the dust has a vaporization temperature $`T_v`$, then $`r_v`$ would be the radius at which the maximum luminosity of the supernova would heat the dust to $`T_v`$. In practice, however, a number of other details about the dust distribution need to be known to relate these parameters (e. g., grain size, composition, and absolute values of the absorption/emission efficiency and the dust density). Furthermore, if the cavity is the result of pre-explosion evolution and not vaporization due to the SN flash, then no relationship would exist. Therefore, we treat $`T_v`$ and $`r_v`$ as independent parameters.

It is convenient to parameterize the cavity size by $`t_v=2r_v/c`$, the epoch at which we first observe light re-emitted by dust on the cavity surface directly behind the supernova along the line of sight. This parameter sets the timescale for the infrared evolution. For models where the supernova light curve decays quickly compared to this timescale, $`t_v`$ is near the “kink” in the IR light curve, where the plateau ends and the dust emission begins to decay. For SN 1995N, the kink in the IR light curve occurs near $`t=1000`$ d (see Fig. 2) so we consider models with $`t_v`$ near this epoch.

#### 3.3.2 Short-Flash Echo Model

Since we have more data on SN 1995N than on the other four Type IIn SNe, we use this event to test the IR echo scenario. The solid line in the lower left panel of Figure 2 shows the IR light curve from a model where the input SN light curve is a short flash compared to the echo timescale $`t_v`$. For the model shown, the SN light curve was exponentially decaying with an e-folding time $`\mathrm{\Gamma }=t_v/10`$ and dust emission/absorption exponent $`n=1`$. For short-flash models where the supernova light curve fades much faster than $`t_v`$, the results are quite insensitive to the exact details of the input light curve and a short “top-hat” light curve would have produced essentially the same results.

The model IR luminosity was matched to the observed plateau luminosity with $`t_v`$ scaled to 950 d to match the observed kink in the IR light curve. As can be seen in Figure 3, the resulting light curve is qualitatively similar to the observed IR luminosity evolution of SN 1995N but it fades slower after the kink. However, this is not a strong test of the dust echo model. The “observed” IR luminosity is an extrapolation of the single temperature continuum fit to the $`HK`$ color and thus is only sensitive to the hottest dust. In contrast, the model light curve shows emission from a superposition of dust over a wide range of temperatures. As the echo ages and the emission comes from farther out in the nebula, the IR emission becomes more and more dominated by cooler dust that is missed by our near-IR photometry.

A better test is to directly compare the observed flux-calibrated spectra to the evolution of the model spectrum. The echo model calculates spectra for each epoch in dimensionless wavelength units $`\lambda `$/$`\lambda _v`$ where $`\lambda _vkT_v/h`$. To make a fit for a given set of parameters, the model spectrum is matched to a single epoch of the data. The model is shifted vertically to match the flux levels and horizontally match the shape of the observed spectrum, essentially fitting $`T_v`$. The same shifts are then applied to the other epochs and the resulting spectral evolution can be compared to the data.

The results of such a procedure are plotted in Figure 3. In the upper panel we plot the observed spectra of SN 1995N with spectra from the short-flash echo model for the corresponding epochs. The model luminosity and $`T_v`$ have been fit to the second spectrum on day 1078. The cavity time $`t_v`$ for this model is 950 d and $`T_v1100`$ K. The resulting spectral evolution compares well with the observations, although the model luminosity fades more quickly. This is in contrast to the plot in the lower left panel of Figure 2 where the model fades more slowly than the “observed” IR luminosity. Again, the reason for the discrepancy is that the “observed” IR luminosities are only extrapolations based on the $`HK`$ temperature fits and are therefore missing a significant component of emission from cool dust. Since the spectral evolution provides a much more direct test, we use this as our primary model discriminant.

While the short-flash echo model is successful enough to suggest that the IR echo scenario may be correct, it involves an unrealistic SN light curve. This model describes an echo where most of the SN energy is released quickly, in this case within the first 100 days. This timescale, while typical for most supernovae, is not really a good description of a Type IIn supernova like SN 1995N. SN 1995N was observed to have a very slowly decaying light curve (Baird et al., 1998), a feature seen in several other Type IIn events (SN 1987F, SN 1978K, SN 1988Z: Filippenko 1997; SN 1997ab: Hagen, Engels, & Reimers 1997). Also, bright X-ray emission was observed in SN 1995N at very late times, even beyond the beginning of the observed IR decline. Both of these facts suggest that a much longer timescale should be used for the light source of the echo model.

A longer source light curve is also easier to reconcile with the very large IR luminosity at late time. A simple interpolation of the “observed” infrared luminosities shown in Figure 2 yields a total output of about $`5\times 10^{49}`$ erg in between days 730 and 2435. Integration of the model light curve yields a value closer to $`10^{50}`$ erg. Thus the energy released in the infrared is significantly larger than a typical $`100`$ day supernova flash, but is perhaps consistent with the long lasting interaction powered emission seen in other bright Type IIn supernovae. Aretxaga et al. (1999) found a lower limit of $`2\times 10^{51}`$ erg for the emission from SN 1988Z, during 8.5 years of observation. Similar emission from SN 1995N would certainly provide enough energy to power the observed dust emission, but would again call for a much longer light curve for the echo model.

#### 3.3.3 CSM Interaction Echo Model

While a short-flash model is relatively successful at matching the observed NIR spectral evolution, the input SN light curve was not a satisfactory description of SN 1995N and we therefore examined echo models with longer supernova light curves. The late-time X-ray flux appears to be relatively constant thereby making a long “top-hat” light curve an attractive possibility. For a long top-hat model, the IR continuum continues to rise until after the SN light source shuts off. If the cutoff occurs prior to $`t_v`$ then the IR luminosity levels off until after $`t_v`$, and then begins to fade. Conversely, if the cutoff epoch is after $`t_v`$, then the IR flux continues to rise until the light source shuts off and then the IR flux quickly decays.

However, what we observe in SN 1995N is not consistent with such a model. The NIR flux faded very slowly prior to the $`t1000`$ d kink and fairly rapidly thereafter with no sign of an X-ray cutoff $`200`$ days after the IR kink. Furthermore, the qualitative results from the long top-hat model hold true for any fairly flat input light-curve. Thus the rapid decline of the IR flux after the kink forces us to the conclusion that most of the total emission from SN 1995N was released well before $`t=1000`$ d even though the late-time X-ray emission was still quite bright.

This suggests a two-component light curve which we model as a top-hat plus constant late-time emission. The results of such a model are plotted in the bottom half of Figure 3. For the model shown, the absorption/emission exponent $`n=2`$, and the echo timescale parameter $`t_v=800`$ d which implies a dust-free cavity radius $`r_v0.3`$ pc. The model supernova light curve consists of a bright “top-hat” plateau that lasts for 400 days, followed by constant luminosity emission 20 times fainter than the bright plateau. The fit value of $`T_v900`$ K for this plot.

Figure 3 shows that this model reproduces the observed IR spectral evolution at least as well as the short-flash echo model. However, where the short-flash model was based on an unsatisfying description of the SN light curve, the two-component model is based on a more reasonable approximation for the UV/X-ray light curve of SN 1995N. What this shows is that the IR excess observed in SN 1995N can plausibly be explained as an infrared echo. In this case the echo is not powered by the flash of the supernova at maximum light but rather by the much more slowly evolving emission from the circumstellar interaction. We stress, however, that the particular model fit shown in Figure 3 is probably not unique. There are enough parameters that one could likely fit a number of echo models to the data.

Because we did not have nearly as much data from the other Type IIn supernovae, we did not attempt to fit specific echo models for them. But the infrared behavior of all of these objects can probably be explained by an interaction echo model. The very large IR luminosities of SN 1997ab and SN 1999Z are easiest to explain with such a model. In SN 1999Z, the rising IR flux seen between day 296 and 465 would suggest that the light source was still emitting a significant amount of energy in this time period. Since we only observed a declining IR flux in SN 1997ab, this would suggest that we caught the IR emission in SN 1997ab after $`t_v`$ in the fading era of the echo.

The behavior of SN 1998S is not suggestive of a spherically symmetric dust echo, but is more consistent with an echo from a highly flattened dust distribution. Echo models for flattened dust nebulae with the line of sight near the equatorial plane do predict a bump in the IR light curve near $`t_v`$ (EC88). This would be in line with other evidence for deviations from spherical symmetry in SN 1998S. Spectropolarimetry indicated that both the SN ejecta and the circumstellar environment were aspherical (Leonard et al., 2000; Wang et al., 2001). Also a well-resolved triple-peak line profile observed for the hydrogen and helium lines of SN 1998S was interpreted by Gerardy et al. (2000) as emission from interaction with dense, clumpy circumstellar gas in a ring or disk surrounding the supernova and probably observed from a fairly high inclination.

## 4 Conclusions

We have presented near-infrared photometry and spectra from five recent Type IIn supernovae that exhibited large late-time infrared excesses. The infrared excesses in these objects seem to have a much slower evolution than those that have been previously presented in the literature, remaining bright many years after maximum light. The resulting energy release via IR dust emission is quite large, approximately $`0.5`$$`1\times 10^{50}`$ erg from SN 1995N. The long IR light curve and large luminosity cannot be reconciled with emission from newly formed dust in the ejecta. In addition, for SN 1995N and SN 1998S, the IR luminosity is significantly larger than the late-time X-ray emission, ruling out emission from a shock precursor. IR emission from collisionally excited dust in shocked gas also seems unlikely as the observed temperatures are fairly high ($`1000`$ K). Moreover, in SN 1995N the IR flux declines rapidly after day 1000, while the X-ray flux shows little sign of fading even $`200`$ days later.

On the other hand, the data do appear to be consistent with an infrared echo of the bright, slowly-fading emission from the circumstellar interaction. Simple IR echo models are able to reproduce the observed near-infrared spectral evolution of SN 1995N and such a scenario helps explain the slow evolution and large IR luminosities seen in these objects. In principle, detailed studies of the IR echos of these objects could provide important information about the circumstellar environment of Type IIn supernovae, and might lend some insight into their progenitors. In practice, however, to provide significant constraints the supernova will have to be very well observed, including early time UV and X-ray observations, as well as good near- and mid-infrared observations of the echo.

We would like to thank the observatory staffs at MDM, FLWO, KPNO, APO, MMT, IRTF, and Subaru for their excellent support, especially Dr. Hiroshi Terada who was enormously helpful during the Subaru observations. C. L. G. and R. A. F.’s research is supported by NSF grant 98-76703. K. N. has been supported in part by the Grant-in-Aid for Scientific Research (07CE2002, 12640233) of the Ministry of Education, Culture, Sports, and Technology in Japan.


text:
# Numerical Implementation of Three-Body Forces in Bound State Faddeev Calculations in Three Dimensions

## 1 Introduction

Three nucleon (3N) bound state calculations are traditionally carried out by solving Faddeev equations on a partial wave basis. After truncation this leads to a set of a finite number of coupled equations in two variables for the amplitude. The calculations are performed either in momentum space , in configuration space , or in a hybrid fashion using both spaces . Though a few partial waves often provide qualitative insight, modern three nucleon bound state calculations need 34 or more different isospin, spin and orbital angular momentum combinations . It appears therefore natural to avoid a partial wave representation completely and work directly with vector variables . This is especially the case, if one wants to consider genuine three-nucleon force (3NF) effects. The true advantage stems from the fact that a 3N calculation is carried out in Jacobi variables, whereas a typical 3NF has the form of two consecutive meson-exchange propagators between e.g. particles 1 and 2, then particles 2 and 3. Using vector variables, the necessary coordinate transformations are numerically realized by interpolations, whereas in a partial wave based calculation a large amount of coupling coefficients has to be evaluated and stored , requiring large storage and memory capabilities of a computer architecture.

## 2 Bound State Equation with a Three-Body Force

The Faddeev component describing the bound state of three identical particles interacting via pairwise forces as well as genuine three-body forces can be written as

$$\psi _i=G_0t_iP\psi _i+(1+G_0t_i)G_0W^{(i)}(1+P)\psi _i,$$
(1)
where $`i=1,2,3`$. In the following, we choose $`i=1`$ without loss of generality. The operator $`t_i`$ stands for the two-body t-matrix in the subsystem $`(jk)`$ summing up the pair interaction in this system. The quantity $`W^{(i)}`$, shown diagrammatically in Fig. 1, is defined via

$$V_{123}=W^{(1)}+W^{(2)}+W^{(3)},$$
(2)
where $`W^{(i)}`$ is that part of the full 3N force $`V_{123}`$ which is symmetric for the exchange of nucleons $`j`$ and $`k`$ ($`jik`$). The decomposition of Eq. (2) is natural, e.g. for the $`\pi \pi `$ exchange 3N force which is present in all currently available 3N forces. The free 3N propagator is given by

$$G_0^1=E\frac{p^2}{m}\frac{3}{4m}q^2,$$
(3)
where p and q are the standard Jacobi momenta

$`𝐩_i`$ $`=`$ $`\frac{1}{2}(𝐤_j𝐤_k)`$ (4)
$`𝐪_i`$ $`=`$ $`\frac{2}{3}(𝐤_i\frac{1}{2}(𝐤_j+𝐤_k)),`$
where $`ijk=123`$ and cyclic permutations thereof. The permutation operator $`P`$ is given as $`P=P_{12}P_{23}+P_{13}P_{23}`$, and the full 3N wave function is related to the Faddeev component by

$$|\mathrm{\Psi }=(1+P)|\psi .$$
(5)
The solution of Eq. (1) with $`W^{(1)}=0`$ is described in detail in Ref.  and shall not be discussed here. For our model calculations we also use Yukawa interactions of the Malfliet-Tjon type here, however, we modify the interaction with a cutoff function of a dipole type .

In order to develop the algorithm, we concentrate on a model 3N force of scalar meson exchange character, which can be written as

$$W^{(1)}=\frac{1}{(2\pi )^6}\frac{a_s}{m_s}g_s^2\frac{F(Q^2)}{Q^2+m_s^2}\frac{F(Q^2)}{Q^2+m_s^2},$$
(6)
where $`m_s`$ represents the mass of the exchanged meson. The momenta $`𝐐`$ and $`𝐐^{}`$ are defined in Fig. 1 and given by

$`𝐐`$ $`=`$ $`𝐩𝐩^{}{\displaystyle \frac{1}{2}}(𝐪𝐪^{})`$
$`𝐐^{}`$ $`=`$ $`𝐩𝐩^{}+{\displaystyle \frac{1}{2}}(𝐪𝐪^{}).`$ (7)
The function $`F(Q^2)=\left[(\mathrm{\Lambda }_s^2m_s^2)/(\mathrm{\Lambda }_s^2+Q^2)\right]^2`$ represents a cutoff function for large momenta $`Q^2`$. The parameters of the 3NF used in the presented calculation are given in Table 1.

For evaluating Eq. (1) we need to calculate matrix elements like

$$\mathrm{𝐩𝐪}|G_0W^{(1)}|\mathrm{\Psi }.$$
(8)
Considering the momentum dependence in the meson propagators of $`W^{(1)}`$, one can easily imagine that the efficiency of any algorithm will crucially depend on the choice of coordinate system(s) when carrying out the integrations over the intermediate momenta. Our numerical evaluation of Eq. (8) is based on the realization that $`W^{(1)}`$ can be interpreted as two independent interactions, first in the subsystem $`(12)`$, then in $`(31)`$. Explicitly, we write

$`\mathrm{𝐩𝐪}`$ $`|W^{(1)}|\mathrm{\Psi }={\displaystyle }d^3p^{}{\displaystyle }d^3q^{}_1\mathrm{𝐩𝐪}|𝐩^{}𝐪^{}_2`$ (9)
$`{\displaystyle d^3p^{\prime \prime }\frac{F((𝐩^{}𝐩^{\prime \prime })^2)}{(𝐩^{}𝐩^{\prime \prime })^2+m_s^2}}`$
$`{\displaystyle d^3p^{\prime \prime \prime }d^3q_2^{\prime \prime \prime }𝐩^{\prime \prime }𝐪^{}|𝐩^{\prime \prime \prime }𝐪^{\prime \prime \prime }_3}`$
$`{\displaystyle d^3p^{\prime \prime \prime \prime }\frac{F((𝐩^{\prime \prime \prime }𝐩^{\prime \prime \prime \prime })^2)}{(𝐩^{\prime \prime \prime }𝐩^{\prime \prime \prime \prime })^2+m_s^2}_3𝐩^{\prime \prime \prime \prime }𝐪^{\prime \prime \prime }|\mathrm{\Psi }}.`$
Here the subscripts $`1,2,3`$ of the bra and ket vectors define the meaning of the related vectors $`𝐩`$ and $`𝐪`$, namely the particle number $`i`$ singled out by $`𝐪_𝐢`$ as indicated in Eq. (4). We would like to point out that in each integration over a piece of the 3NF we only integrate over three variables, namely the magnitude of a momentum and two angle variables. It turns out that it is most favorable to choose momenta $`𝐪^{}`$ ($`𝐪^{\prime \prime \prime }`$) parallel to the z-axis, so that no interpolation on the unknown function $`\mathrm{\Psi }(p,q^{},\widehat{p}\widehat{q}^{})`$ needs to be performed. The calculation of the transformation from one subsystem to another, e.g. $`{}_{2}{}^{}𝐩^{\prime \prime }𝐪^{}|𝐩^{\prime \prime \prime }𝐪^{\prime \prime \prime }_{3}^{}`$, requires three-dimensional interpolations. We employ cubic Hermite splines , which prove to be both accurate and fast.

## 3 Computational Approach

The discretized Faddeev equation for a bound state (neglecting spin degrees of freedom) is an integral equation in three variables on a typical grid size of 65\*65\*42 (momentum magnitudes, $`p`$, $`q`$, and angle between the momentum vectors). A priori the multidimensionality of the integral equation to be solved requires more memory. However, on an MPP system the number of variables and thus the memory do not pose a computational problem, since a variable defining a specific dimension of the grid can be distributed over a number of processors, leaving a lower dimensional grid on each processor. As such, our three-dimensional approach is ideally suited as MPP application, and we can achieve an almost perfect load balance in our runs.

The eigenvalue equation for the bound state is solved iteratively by using Lanczo’s type techniques, here the method of iterated orthogonal eigenvectors . For a typical run ten orthogonal eigenvectors are calculated per energy. We need about 5 to 7 energy iterations to find the ground state energy.

The calculation of the kernel of the integral equation means evaluating the matrix elements given on the right-hand side of Eq. (1) on a fixed grid $`p`$, $`q`$, and angle $`x=\mathrm{cos}(\widehat{p}\widehat{q})`$. The two-body t-matrix (with the two-nucleon pair interaction as the driving term) is obtained by solving a system of linear equations of the form $`Ax=b`$, where A is typically a 2500\*2500 matrix. This system is solved for about 60 different vectors $`b`$, distributed over correspondingly many processors. The integrations over the 3NF terms are also distributed over the same number of processors, which depends on the size of the q-grid. Details about the different grid sizes and the dependence of the numerical accuracy on them can be found in Ref. .

## 4 Results and Discussion

For our model calculation we use Yukawa interactions of the Malfliet-Tjon  type

$`V(𝐩^{},𝐩)`$ $`=`$
$`{\displaystyle \underset{\alpha =a,r}{}}`$ $`{\displaystyle \frac{g_\alpha ^2}{(𝐩^{}𝐩)^2+m_\alpha ^2}}\left({\displaystyle \frac{\mathrm{\Lambda }_\alpha ^2m_\alpha ^2}{\mathrm{\Lambda }_\alpha ^2+(𝐩^{}𝐩)^2}}\right)^2`$ (10)
with a short-ranged repulsive and a long-ranged attractive piece. The parameters used in our calculation are given in Table 2. A calculation of the three-body binding energy with this pair force gives $`E_t`$=7.699 MeV. The parameters of the 3NF are then adjusted such that its inclusion in the calculation gives a binding energy $`E_t`$=8.590 MeV, a value close to the measured one. The full 3N wave function is calculated from the solution of Eq. (1) using Eq. (5) and shown in Fig. 2 as a function of the momenta $`p`$ and $`q`$ at the fixed angle $`x=1`$. A coparison of the 3N wave functions calculated with the 2N pair forces alone and with the inclusion of the 3NF is shown in Figs. 3 and 4. Both figures depict slices through the wave functions at specific values.

The 3NF in our model calculation is attractive, as expressed by the larger binding energy. One would expect that the system becomes slightly smaller and acquires more high momentum components. Since the wave functions are both normalized to one, an increase of high momentum components will be seen as a decrease in low momentum components, as seen in Fig. 3. Once the wave function decreases by some orders of magnitude, the differences disappear.

In summary, an alternative approach to state-of-the-art three-nucleon bound state calculations, which are based on solving the Faddeev equations on a partial wave basis, is to work directly with momentum vectors. We formulate the Faddeev equations for identical particles as a function of vector Jacobi momenta, specifically the magnitudes of the momenta and the angle between them, for the case when not only pair forces act between the particles but also genuine three-body forces. As model forces we concentrate on scalar forces, a superposition of an attractive and repulsive Yukawa interaction for the pair force and an attractive Yukawa interaction for the three-body force. We demonstrate the numerical feasibility and accuracy of the solution. We want to point out that the incorporation of the 3N forces is less cumbersome in a three-dimensional approach, and the algorithm can be made quite efficient on parallel architectures.

Acknowledgments This work was performed in part under the auspices of the U. S. Department of Energy under contract No. DE-FG02-93ER40756 with Ohio University. We thank the the Neumann Institute for Computing (NIC) and the National Energy Research Supercomputer Center (NERSC) for the use of their Cray-T3E computers.