text:
# Simulation Software: Then, Now and Virtual Observatory

## 1. Introduction

The increased computing speed of both off-the-shelf and dedicated hardware such as the GRAPE series have made it possible to write increasingly complex simulation software for very large N-body systems. What started as simple few hundred line FORTRAN programs (for a review, see e.g. Aarseth 1999) with their individual analysis routines, have grown into mid-size packages of libraries and toolsets. Looking at other fields these will likely evolve into sophisticated large scale frameworks such as ROOT (Brun et al. 1999) and AIPS++ (Glendenning 1996) if the community combines their programming efforts and reuse code.

The plan of this paper is as follows. In Section 2 we will show some of the current techniques of simulation software, and where this could lead. In Section 3 the impact of a Virtual Observatory will be discussed, followed by conclusions and future developments.

## 2. Simulation Software

In the Dark Ages simulation software was developed in small, often one-person, teams for then available “supercomputers”. Possibly due to limited electronic communication, software also did not migrate very easily between researchers. The codes Aarseth developed have arguably been most widely spread and used (see also Binney & Tremaine, 1987). Another early example of shared code development was the OLYMPUS programming system, as described in Hockney & Eastwood (1981). These, and upcoming data reduction packages such as AIPS, GIPSY, IRAF and MIDAS in observational astronomy have led to a number of packages for particle simulations. Although these are still written by fairly small teams, they have now attracted a moderate number of users. But most notably, they have started to attract developers, as is very common in Open Source development these days. One can distinguish two types of packages: on the one hand there are NEMO and Starlab, which present themselves to the user as a collection of programs that can be called from a shell, or as a collection of subroutines and functions with which new tools can be built. On the other hand there are single programs such as tipsy and astroMD, which come with their own programmable interface. Our field is, of course, not that much different from that of for example High Energy Physics (cf. De Angelis 2002).

### 2.1. NEMO, ZENO

Initial work on NEMO started in 1986 by Joshua Barnes, Piet Hut and Peter Teuben (Barnes et al. 1987), and has been subsequently extended (Teuben 1995). This paper also serves as an update to document the recent upgrade from Release 2 to 3, of which many details can be found on the NEMO website111See also http://www.manybody.org which hosts a number of N-body resources.

#### Source Code

The source code consist of two source code trees: a “src” tree and a “usr” tree, which resp. hold the basic NEMO source code, and various public (mostly N-body) codes graciously supplied by their respective authors. Most codes in the “usr” tree are available “as-is”, some have enhancements for support within the NEMO environment. Currently the “usr” tree is already about 4 times larger than the “src” tree (860 KLOC vs. 193 KLOC).

Installation in Release 3 has been largely simplified by using current techniques like autoconf, and using a source code revision control system (CVS) to simplify shared development. This has also made it easier to create binary releases.

The source code is largely written in C, with some C++ and FORTRAN (and support to simplify linking the languages). Table 1 lists some of the public code now available in NEMO.

#### Packages

NEMO’s software is packaged and grouped around a number of data formats. In the nbody package, various programs exist to integrate N-body systems with a wide variety of types of integrators, codes to initialize N-body systems, and visualization and analysis programs. One of the versatile plotting programs within this group is called snapplot, with which any body variable can be plotted vs. any other body variable, using on-the-fly code generation (dynamic object loading) for fast and flexible analysis. One of its derived programs is snapgrid, which produces an image instead of a scatterplot, and includes effects such as optical depth, and can be more directly compared to observations. In addition to snapshots and images, two other data formats have a large set of analysis tools: orbits and tables. Associated with orbits are potential descriptors, which allow for user supplied potentials to be loaded into various orbit integrators without the need to recompile those programs.

In the following example from the nbody/image group of programs an exponential disk is created, and integrated through a few dynamical times such that a bar will form. The disk is then viewed from some angle and a first and second moment along the new Z axis is then used to compute a velocity field and velocity dispersion map on a grid in projected X-Y space. The resulting dataset can be converted to a FITS file and manipulated in external packages, such as saoimage. The resultant view of snapplot and ccdplot is shown in Figure 1.

```
% mkexpdisk - 20000 rcut=2 \                  #  make an exponential disk
  | hackcode - - tstop=4 \                    #  integrate to bar formation
  | snaprotate - - 60,45 xz \                 #  rotate around a bit
  | snapgrid - - zvar=-vz moment=-1 times=4 \ #  grid snapshot to velfield
  | tee ngc999vel.ccd \                       #  save a copy of the data
  | ccdplot - contour=-1:1:0.2 blankval=0     #  contour plot
% ccdfits ngc9999vel.ccd ngc9999vel.fits      #  convert image to FITS format
% ds9 ngc9999vel.fits &                       #  display with saoimage

```

#### Data Formats

NEMO’s data format is a structured binary format, where data elements are identified by name and type, and can be nested at an arbitrary level. I/O routines access these data in an associative manner, only retrieving the data needed at that time. Data is also interchangeable between machines with different data types.

With the introduction of more “foreign” integrators into NEMO, each with their own data format, it became necessary to be able to read and write a large variety of data formats. Because there is no such thing as a standard interchange format like FITS 222N-body data can often be described as a table, and thus the FITS BINTABLE format could very well serve as an an interchange format, see e.g. Teuben 1995, NEMO’s snapshot format has become the central format to interchange format X to Y.

#### ZENO

The ZENO software package is an evolutionary product of an earlier version of NEMO, written by Joshua Barnes, and is largely still source code compatible with NEMO. ZENO instead concentrates on N-body and SPH simulations, and particle representations are dynamically extendible instead.

### 2.2. Starlab

Starlab (Portegies Zwart et al. 2001) was loosely modeled after NEMO, but written completely from scratch to handle the more intricate physics of collisional dynamics. A new tree-based data structure was introduced to handle the more complex stellar interactions. In addition, data piped through the system would not lose information not known to the data-handler. The code is mostly written in C++, and as of this writing consists of about 236 KLOC in 911 files. The kira integrator (23 KLOC) can also be linked with a variety of GRAPE libraries to take advantage if this hardware is available. kira also contains the SeBa stellar evolution module.

Here is a simple example how to create a King model with a given IMF, binaries and stellar evolution, then integrated on a GRAPE-6 and output dumped in a file that can be processed using a number of tools available in Starlab.

```
% makeking -n 500 -w 5 -i -u \                # 500 particle king model
   | makemass -f 1 -x -2.0 -l 0.1 -u 20 \     # mass spectrum
   | makesecondary -f 0.1 -l 0.1 \            # make secondary stars
   | add_star -Q 0.5 -R 5 \                   # add stellar evolution
   | scale -M 1 -E -0.25 -Q 0.5 \             # scale to virial equilibrium
   | makebinary -f 1 -l 1 -u 1000 -o 2 \      # set orbitals for binaries
   | kira_grape6 -t 100 -d 1 -D 10 -f 0.3 -n 10 -q 0.5 -G 2 -S -B  > big.out

```

Visualization in Starlab is currently best done with partiview, which can take advantage of the hierarchical space-time nature of the data and also knows about double stars (Teuben et al. 2001).

### 2.3. Scripting

In the tradition of UNIX, the modular design of NEMO and Starlab (and their common general data format) programs can be combined in a large variety of ways. Prototyping can now be done in minutes, to produce complex and sophisticated analysis pipelines. However, we have found in practice that the resulting shell scripts (sh, csh, make) are often not very robust. Subsequently scripting has exhibited some fragility that more comprehensive control mechanisms could help prevent. Future versions or new software should make use of such control data in a more reliable way. Another approach is to use an embeddable scripting language, such as python or ruby, which can enforce a tighter connection between codes and data. In addition, these hybrid software environments often lend themselves to better GUI development. We are currently experimenting with this.

### 2.4. TIPSY

Another popular package to analyze particle simulations is TIPSY333See also http://www-hpcc.astro.washington.edu/tools/tipsy/tipsy.html (Quinn & Katz). This package has followed the philosophy of a single program, with a special command interpreter to operate on snapshots containing a combination of SPH, dark-matter, and pure Newtonian particles.

```
  % tipsy
  <yes, Peter>openascii run99.ascii
  <yes, Peter>readascii run99.bin
  read time 14.970800
  <yes, Peter>xall
  <yes, Peter>quit
  <I will miss you, master>
  % tipsysnap run99.ascii - | snapplot - xrange=-4:4 yrange=-4:4

```

The obvious advantage of this approach is speed, as the data always remains in memory. However, any modifications to the program means it will have to be abandoned in order to be recompiled. Plug-ins or dynamic objects can alleviate some of these problems.

### 2.5. IDL, and other

As the field of astrophysical particle simulations is a very specialized one, the majority of research is done with personal codes, and likewise their analysis. In recent years these codes have also incorporated more interesting physics, such as simple empirical sticky particle dynamics, SPH gas dynamics, stellar evolution, chemodynamical evolution, etc. These codes are mostly written in languages like FORTRAN, C or C++. In recent years commercial graphics packages like IDL and Open Source toolkits such as VTK (Visualization ToolKit) have also become popular to analyze and visualize such complex datasets. Generic visualization packages such as AVS, IRIS Explorer and IBM’s Data Visualizer are widely used, yet limitations in such generic packages continue to create a niche for programs like the recently developed AstroMD toolkit (Becciani et al. 2000). This program uses the VTK library to allow for very sophisticated multi-variate data analysis and visualization.

## 3. Virtual Observatory

The concept of a Virtual Observatory to federate various observational databases (e.g. Brunner et al. 2001) and make combined searching and analysis on these databases possible, has not gone un-noticed in the theoretical community. Teuben et al. 2002 argued that adding various types of theoretical data to the VO will benefit observers as well as theoreticians, and open up new and unexplored avenues for research.

For example, existence of a standard number of (benchmark) datasets will benefit authors of new codes to quickly compare and highlight differences between various codes. This is of course not something new, but still a relatively rare event in our community. After the first published code comparison by Lecar (1968) it still took nearly 30 years for the community to continue this effort when in 1997, Heggie (2002, this volume) reported on a comparison between different star cluster simulations and Sellwood (1997) published a comparison between five different N-body codes typically used in galaxy simulations. Setting up test problems is important (Heggie 1997, 2001), and has been a standard in many field of computational science (e.g. Stone & Norman 1992).

In a Virtual Observatory we can expect to select models, compare them to existing data and perform various types of fits to best describe the observations. In addition, new models can be compared to old models, and provide feedback to code development. In order to better understand the scope of the role of theory in a Virtual Observatory, we have started to construct a “toymodel”444See http://www.astro.umd.edu/nemo/tvo, which contains a growing collection of different types of theory data. The only condition for data to be added to this toymodel is that they must either be benchmark data, or datasets associated with published papers.

## 4. Conclusion

We have reviewed the evolution of particle simulation codes, and seen them adapt to the ever growing speeds of hardware. This will have to include the rapid development of PC cluster hardware, which will require enhancing the scalability of parallel algorithms. The simulation software itself has more slowly matured (e.g. compare different software engineering practices) with promises of code reuse and extendibility. The future in software development is likely going to be in frameworks such as ROOT and AIPS++ with plenty of room for niche applications, as long as they can easily share their data!

A Virtual Observatory framework will allow for a more seamless integration of observations and simulations, and allow astronomers to compare observations and obtain best fit models. It should also enable theory to develop more in pace with simulations and encourage code reuse and data format sharing.

It is expected that the N-body simulation community will continue to contribute, as in the past, a diversity of insights that will continue to make the field exciting and productive for years to come.

##### Acknowledgments.

I acknowledge support from the Alfred P. Sloan foundation, through a grant to Piet Hut for research at the Hayden Planetarium of the American Museum of Natural History in New York. I also wish to thank the AAS International Travel Grant for support; Joshua Barnes and Piet Hut for our collaborations in NEMO; the numerous users of NEMO in (requests for) enhancing the software; Vicky Johnson for a careful reading; George Lake and Tom Quinn for long discussions around the future of particle pushing analysis frameworks, a.k.a. the N-Chilada.

## References

Aarseth, S. 1999, Publications of the Astronomical Society of the Pacific, 111, 1333

Angelis de, A. 2002 astro-ph/0112371

Barnes, J.E., & Hut, P. 1986, Nature, 324, 446

Barnes, J.E., Hernquist, L, Hut, P. & Teuben, P. 1988, B.A.A.S. 20, 706

Becciani, U., Antonuccio-Delogu, V., Gheller, C., Calori, L., Buonomo, F., & Imboden, S. 2001, astro-ph/0006402. (IEEE conference proceeding)

Binney, J.J. & Tremaine, S. 1987, “Galactic Dynamics”, Princeton University Press

Brun, R. & Rademakers, F. 1997 Proceedings AIHENP’96 Workshop, Lausanne, Sep. 1996, Nucl. Inst. & Meth. in Phys. Res. A 389 (1997) 81-86.

Brunner, R.J., Djorgovski, S.G. & Szalay, A.S. (eds) ASP Conf. Ser., Vol. 225, Virtual Observatories of the Future, p.353

Dehnen, W. 2000, ApJ, 536, L39

Glendenning, B. 1996 in ASP Conf. Ser., Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco: ASP), p.271

Heggie, D.C 2002 - Kyoto-I and II experiment, this volume; astro-ph/0110021

Hockney, R.W., Eastwood, J.W., 1981, “Computer Simulations using Particles”, McGraw-Hill International, New York.

Hut, P., and Makino, J. 1999, Science 283, 64

Lecar, M. 1968, (in: Colloque sur le problem gravitationnel des N corps) Bull.Astr. 3, 91

Portegies Zwart, S.F., McMillan, S.L.W., Hut, P. & Makino, J., 2001, MNRAS 321, 199

Sellwood, J. A. 1997, in Computational Astrophysics, eds. D. A. Clarke & M. J. West (San Francisco: ASP Conference Series 123), 215

Stone, J.M., & Norman, M.L. 1992 ApJS, 80, 753

Teuben, P.J. 1995, in ASP Conf. Ser., Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco: ASP), p398

Teuben, P.J. et al. 2001, in ASP Conf. Ser., Vol. 238, Astronomical Data Analysis Software and Systems X, ed. F.R. Harnden, Jr., F.A. Primini, & H.E. Payne (San Francisco: ASP), p.499

Teuben, P.J. et al. 2002, in ASP Conf. Ser., Vol. xxx, Astronomical Data Analysis Software and Systems XI, ed. D.A. Bohlender, D. Durand, T.H. Handley (San Francisco: ASP), p.??? (see also astro-ph/0111478)

## Appendix: Data Usage Survey

Before the final “Open Discussion” (chaired by Prof. Sugimoto) a small survey was handed out to all participants in order get some more insight in the types of data that are produced in our simulations, and the current habits of its practitioners. Although no head-count was made of the number of people present at this session, 53 (exactly 50% of the 106 officially registered participants) forms were returned and analyzed.

A similar (but unpublished) survey was held in 1994 amongst a much smaller but similar focus group of astronomers at the Lake Tahoe meeting to celebrate Sverre Aarseth’s 60th birthday. The current ”survey” was intended to aid the discussion on typical current ”N-body” data usage and future Virtual Observatories. A quick overview of the results was given during the summary session on Friday afternoon. Here we reproduce some of that discussion.

### The Survey

```
1. Your name or alias:

2. Do you:  (in all of the items below you are allowed to
Ψ    mark multiple items, perhaps you can indicate
Ψ    the percentage of relevance in each)
      a  develop your own code(s)Ψ[name(s):
      bΨ use existing code(s)ΨΨ[name(s):
      c  other
      d  N/A

3. What is your typical range in N:

4. Do you
      aΨ save particle/grid data (see also 5)
      bΨ discard data (e.g. analysis within simulation code)
      cΨ other:
      d  N/A

5. Your data format is mostly:
      aΨ table of particle attributes (m, x, y, z, vx, ..)
      bΨ grid of cell attributes (den(i,j,k), vx(i,j,k), ..)
      cΨ hydrid of 1 and 2 (e.g. P3M)
      d  a tree structure
      e  other: (try and describe)
      f  N/A

6. What analysis software do you use:
      a  my ownΨΨΨ[name(s):    / language(s):
      bΨ sm
      cΨ IDL
      dΨ NEMO
      eΨ ZENO
      fΨ Starlab
      gΨ Tipsy
      hΨ other:
      i  N/A

7. What kind of ancillary data do you also store:
      aΨ minor diagnostics (E, Lx, Ly, Lz, ...)
      bΨ neighrest neighbor list
      c  detailed grouping info
      dΨ a tree structure representing particles
      eΨ a tree structure representing grid
      fΨ other:

8. Do you use models
      a  to compare to your own other models
      bΨ to compare to other people models
      cΨ to compare/fit to observational data
      dΨ other:
      e  N/A

9. If you do any of those, do you
      a  compare in configuration (pos, vel) space (e.g. velocity field)
      b  compare in derived quantities  (e.g. power spectrum, luminosity function)
      cΨ other:
      d  N/A

10.Any suggestions for Virtual Observatories?
      a  a waste of time, because ....
      b  a great idea, but ....
      c  other:
      d  N/A

11.Do you have any data yourself you have available?
      a  no
      b  yes, but I can’t make it available
      c  yes, and I might be able to make it available
      d  yes, they are available on the web already

12. Any remaining comments?

```

### Results

It quickly became clear that making a good survey is harder than it looks. Many questions had multiple answers, thus the numbers quotes will be percentages and will generally add up to over 100%.

1. Nobody was using an alias, everybody choose to use their known name (one person choose his/her Japanese name).
2. a) 72% b) 68% c) 4% d) 2
Codes mentioned (quite a few respondents actually missed the request to mention names of their code(s)): AMR using boxlib, ap3m, asph, c++tree, chameleon, COSMIC, eurostar, gadget, gasoline, gizmo, grapesph, hydra, kira, nemo, nbodyN, p3m, p3msph, pg, pkdgrav, pmtreecode, ptreecode, scf, superbox. tipsy, treescf, treesph.
3. 2 - 1,000,000,000, with one respondent quoting $`32768`$!
4. a) 90% b) 6% c) - d) -
data summaries, expansion coefficients, test vector.
5. a) 84% b) 14% c) 14% d) 14% e) 6% f) 4%
series of nested adaptive grids, distribution function,
6. a) 60% b) 36% c) 32% d) 10% e) - f) 10% g) 20% h) - i) :
xmgr(2), gnuplot(8), mongo(4), dx(2), amrvis(1), pgplot(1), pdl (1), midas(1), sm(1)
Languages: Fortran (8), C (6), C++ (4), Perl (2), Awk (1)
7. a) 74% b) 22% c) 12% d) 12% e) - f) -
special purpose diagnostic files, likelyhoods, never store data.
8. a) 74% b) 66% c) 70% d) - e) 4%
9. a) 64% b) 64% c) - d) 4%
10. ambitious, good luck, hope it takes off, would not use it  (2x)
11. a) 18% b) 6% c) 50% e) 16%
12. A nice variety of comments: The VO might be a good place to offer their codes to the public. There needs to be a common data format (CDF, HDF). The project is very ambitious. Worry that the credit to a persons work is lost. VO should be looked at as an internationally funded observatory. Theory should also play their parts in making images available for public outreach. Programming resources are needed to support the scientists who are supposed to this work. Reliability and possible refereeing needed for the submitted data. Good luck! Useful to chain simulations: the end of one simulation can be used as input to the next.

### Discussion and Conclusions

Although a very large programming effort is shared amongst this community, at least an equal amount is using simulation software from colleagues (and probably expanding it). The survey unfortunately did not address the question how data is interchanged between codes, if any, despite that a good fraction of people use more than one code.

Data analysis fills an almost equally wide spectrum, from a variety of special purpose software (TIPSY, NEMO) to utilizing generic, and sometimes even commercial, software (IDL, sm, perl, dx). Part of this is no doubt sociological and the familiarity of the user with that specific software.

Although a large amount of software is still written in FORTRAN, the majority is now in C/C++. Perhaps notable in this survey was the absence of Java. perl and awk are popular scripting languages, although nobody mentioned the basic shells sh and csh. Also, no explicit operating systems were mentioned. The large amounts of code, and interesting genealogy between them, shows that most of them fill a specific niche in the market.


format: text

text:
# Anisotropic generalizations of de Sitter spacetime

## 1. Introduction

As is well known, the energy-momentum tensor describing ordinary matter is of the Plebański type $`[TS_1S_2S_3]_{(1111)}`$, i.e. it admits one timelike and three spacelike eigenvectors ($`u^\mu `$ and $`(m^a)^\mu `$ with $`a=1,2,3`$, say). The eigenvalue corresponding to the timelike vector field $`u`$ is (minus) the energy density of the matter $`ϵ`$, while the three spacelike eigenvectors are the principal stresses $`\lambda _a`$ ($`a=1,2,3`$) (obviously, the perfect fluid is space-isotropic, so that for a fluid the spacelike eigenvalues coincide). Using as tetrad basis that defined by the eigenvectors, the energy momentum tensor can be written in the form

(1.1)
$$T_\nu ^\mu =ϵu^\mu u_\nu +\underset{a}{}\lambda _a(m^a)^\mu (m^a)_\nu $$

From now on we work in spherical symmetry (see ), where only two spacelike eigenvalues can be distinct, so that the canonical form is

(1.2)
$$T_\nu ^\mu =ϵu^\mu u_\nu +p_rm^\mu m_\nu +p_t\mathrm{\Delta }_\nu ^\mu ,$$

where $`m^\mu `$ denotes the unit radial vector and $`\mathrm{\Delta }`$ denotes the projector onto the two-dimensional spacelike subspace orthogonal to $`m`$. The radial and tangential stresses are respectively denoted by $`p_r`$ and $`p_t`$.

The de Sitter spacetime can be viewed as a vacuum solution with nonzero ”lambda” term but also, as is well known, as the spacetime originated by a lambda term ”source”, i.e. as the solution of the Einstein field equation in matter having $`T_\nu ^\mu =ϵ_0\delta _\nu ^\mu `$ where $`ϵ_0`$ is now constant due to the field equations for matter, which obviously imply $`_\nu ϵ_0=0`$. In the formula (1.2), one thus has $`ϵ=p_t=p_r=ϵ_0`$ so that the energy tensor is completely degenerate.

One may now ask if there are ways, in which it is possible to weaken the hypothesis of complete degeneracy in a ”minimal” way, i.e. retaining two degenerate subspaces, thereby obtaining ”anisotropic generalizations” of de Sitter spacetime. It is obvious, that perfect fluids cannot achieve this goal. However, if the stress is anisotropic then one can search for solutions described by

(1.3)
$$T_\nu ^\mu =A\mathrm{\Gamma }_\nu ^\mu +p_t\mathrm{\Delta }_\nu ^\mu $$

where $`A`$ is some function and $`\mathrm{\Gamma }`$ is the unit tensor living in the two-dimensional subspace spanned by the timelike and the radial spacelike eigenvectors (i.e. $`\mathrm{\Gamma }_\nu ^\mu =u^\mu u_\nu +m_\mu m^\nu `$). Interestingly enough, the equations of motion for the matter now shows that $`ϵ`$ is a function depending on $`R`$ only, where $`R`$ is the comoving area radius.

In the present paper this class of solutions is derived, together with minimal physical requirements that have to be imposed. It also includes models used to build up regular black hole interiors , where a regular solution is used to replace a singular core (e.g. Schwarzschild). It may be worthwhile noticing that here the starting point is given by a condition on the constitutive equation (2.4), (see (2.7) below) recovering a solution which in principle is neither static nor regular.

## 2. The solution

Consider a spherically symmetric object, whose general line element in comoving coordinates may be written as

(2.1)
$$\text{d}s^2=e^{2\nu }\text{d}t^2+e^{2\lambda }\text{d}r^2+R^2(\text{d}\theta ^2+\mathrm{sin}^2\theta \text{d}\varphi ^2)$$

(where $`\nu ,\lambda `$ and $`R`$ are function of $`r`$ and $`t`$). With the energy–momentum tensor describing the matter given by (1.2), Einstein field equations for this model reads

(2.2a) $`m^{}=4\pi ϵR^2R^{},\dot{m}=4\pi p_rR^2\dot{R},`$
(2.2b) $`\dot{R}^{}=\dot{\lambda }R^{}+\nu ^{}\dot{R},`$
(2.2c) $`p_r^{}=(ϵ+p_r)\nu ^{}2{\displaystyle \frac{R^{}}{R}}(p_rp_t),`$
where a prime and a dot denote partial derivative with respect to $`r`$ and $`t`$ respectively, and $`m`$ is the Misner–Sharp mass, defined as

(2.3)
$$m(r,t)=\frac{R}{2}\left[1(R^{})^2e^{2\lambda }+\dot{R}^2e^{2\nu }\right].$$

The equation of state for a general material in spherical symmetry can be given in terms of a state function (see e.g. )

(2.4)
$$ϵ=ϵ(r,R,\eta ),$$

where $`\eta =e^{2\lambda }`$, in such a way that the stresses, which are in general anisotropic, are given by the following relations:

(2.5)
$$p_r=ϵ+2\eta \frac{ϵ(r,R,\eta )}{\eta },p_t=ϵ\frac{R}{2}\frac{ϵ(r,R,\eta )}{R}.$$

It follows, that since, from (1.3), the solutions we are searching for are uniquely characterized by the condition that

(2.6)
$$ϵ+p_r=0,$$

then $`ϵ`$ has to be independent from $`\eta `$:

(2.7)
$$\frac{ϵ}{\eta }=0.$$

Substituting (2.5) in (2.2c) and using $`ϵ^{}(r,t)=\frac{ϵ(r,R)}{r}+\frac{ϵ(r,R)}{R}R^{}`$ coming from (2.7), we get also $`\frac{ϵ}{r}=0`$, that is $`ϵ`$ given by (2.4) must be a function of $`R`$ only. We are thus left with the following energy-momentum tensor:

(2.8)
$$T_\nu ^\mu =ϵ\mathrm{\Gamma }_\nu ^\mu (ϵ+\frac{R}{2}\frac{\text{d}ϵ}{\text{d}R})\mathrm{\Delta }_\nu ^\mu ,$$

and since equations (2.2a) now read

$$m^{}=4\pi ϵ(R)R^2R^{},\dot{m}=4\pi ϵ(R)R^2\dot{R},$$

it must be

(2.9)
$$m(R)=4\pi _0^Rϵ(\sigma )\sigma ^2\text{d}\sigma +m_0.$$

In order for this solution to satisfy minimal requirements of acceptability, the weak energy condition (w.e.c.) is imposed, that reads

(2.10)
$$ϵ0,ϵ+p_r0,ϵ+p_t0$$

(a basic reference for a discussion of energy conditions is ). In this case (2.5) implies that w.e.c. is satisfied if $`ϵ(R)`$ is a non negative and not increasing function of $`R`$.

We are left with the system (2.2b), (2.3), and (2.9) in the unknown $`R,\lambda ,\nu ,ϵ`$ and $`m`$. In principle (2.2b) should be integrated in order to obtain exact solutions, but recalling, from (2.3) and (2.9), that $`(R^{}e^\lambda )^2(\dot{R}e^\nu )^2`$ is a function of $`R`$ only, we will limit ourselves to the case when the two addenda are separately function of $`R`$ only:

(2.11)
$$R^{}e^\lambda a(R),\dot{R}e^\nu b(R),$$

where $`a`$ and $`b`$ are two prescribed functions of $`R`$.

Using condition (2.2b) we find, up to time reparameterization, the following expressions:

$`\lambda =\mathrm{log}(f(r)b(R)),\nu =\mathrm{log}a(R),`$
$`m(R)={\displaystyle \frac{R}{2}}[1a^2(R)+b^2(R)],ϵ(R)={\displaystyle \frac{1}{4\pi R^2}}{\displaystyle \frac{\text{d}m(R)}{\text{d}R}},`$
and the two (compatible) equations for $`R(r,t)`$:

(2.12)
$$\dot{R}=a(R)b(R),R^{}=f(r)a(R)b(R).$$

Here the function $`f(r)`$ arises as an integration term, and since the curve $`R(r,0)`$ is the initial data that will be conveniently taken equal to $`r`$, it is $`R^{}(r,0)=1`$ and then $`f(r)=(a(r)b(r))^1`$.

With the position

$$v(R):=a(R)b(R),$$

equations (2.12) reads

(2.13)
$$\dot{R}=v(R),R^{}=\frac{v(R)}{v(r)},$$

that integrate to give

(2.14)
$$R(r,t)=V^1(t+V(r)),$$

where $`V(\sigma )`$ is a primitive of $`\frac{1}{v(\sigma )}`$ (note that $`V`$ is invertible). The line element then takes the form

(2.15)
$$\text{d}s^2=a^2(R)\text{d}t^2+\frac{v^2(R)}{v^2(r)a^2(R)}\text{d}r^2+R^2\text{d}\mathrm{\Omega }^2,$$

with $`R=R(r,t)`$ given by (2.14).

## 3. Physical interpretation and remarks

Defining the function

(3.1)
$$\chi (R)=1\frac{2m(R)}{R}=a^2(R)b^2(R)=a^2(R)\frac{v^2(R)}{a^2(R)},$$

it can be seen that two suitable coordinate changes map the regions $`\{\chi >0\}`$ and $`\{\chi <0\}`$ respectively into a variation of mass of Schwarzschild solution. The sign of $`\chi `$ defines the character of $`R`$ viewed as a coordinate: if $`\chi `$ is positive, $`R`$ may be regarded as a ”length”, if $`\chi `$ is negative $`R`$ can be seen as a ”time” instead.

It can be noticed that these solutions can arise from Kerr-Schild geometry, since another suitable coordinate change may be applied, to recover the form

(3.2)
$$\text{d}s^2=(\text{d}\overline{t}^2+\text{d}\overline{r}^2+\overline{r}^2\text{d}\mathrm{\Omega }^2)+(1\chi (\overline{r}))(\text{d}\overline{t}+\text{d}\overline{r})^2,$$

that shows, by the way, that $`\chi =0`$ is a removable singularity. The family outlined here in fact contains Minkowski, Schwarzschild and de Sitter spacetimes as particular cases, corresponding to choosing the function $`m(R)`$ respectively equal to 0, to $`m_0`$ (constant), and to $`\frac{4}{3}ϵ_0\pi R^3`$.

As sketched before, the line element of these metrics are formally analogue to Schwarzschild exterior and black hole (coinciding with them if $`m(R)`$ is constant), so it is not surprising that $`\chi =0`$ is not a true singularity. Computation of Kretschmann scalar $`K=R_{abcd}R^{abcd}`$ yields

(3.3)
$$K=4\frac{(\chi 1)^2}{R^4}+4\frac{(\chi ^{})^2}{R^2}+(\chi ^{\prime \prime })^2,$$

and this suggest that the only singularity can occur at $`R=0`$. Let us assume $`m`$ (and therefore $`\chi `$) analytic in a right neighborhood of 0, and let $`\alpha 0`$ such that $`m(R)R^\alpha `$. It is easily seen that metric regularity at $`R=0`$ occurs if and only if $`\alpha 3`$. On the other side, since w.e.c. (2.10) in terms of $`m`$ reads

(3.4)
$$m^{}(R)0,m^{\prime \prime }(R)\frac{2}{R}m^{}(R)0,$$

in order for (3.4) to be satisfied near $`R=0`$ it must be $`\alpha 3`$. Therefore, the only case for a physically acceptable metric not to be singular at $`R=0`$ is when it behaves asymptotically as de Sitter spacetime as $`R`$ approaches 0 (that is, $`m(R)R^3`$). In all other cases (allowed by w.e.c.) $`R=0`$ is a true singularity.

The solutions here found can be matched to Schwarzschild spacetime if and only if $`ϵ`$ vanishes on a $`R=\text{const.}`$ surface, that is there exists $`R_b>0`$ such that $`\frac{\text{d}m}{\text{d}R}(R_b)=0`$. This can be performed either in the static region (i.e. where $`\chi >0`$), or in the non static one. In the first case, we obtain a globally static object.

If, instead, we allow $`\chi `$ to vanish for some $`R_0>0`$ (Cauchy horizon), the solution may enter the non static region, where the matching can be performed obtaining physically valid black hole interior models. This is not in contradiction with the result proved by Baumgarte and Rendall in , stating that the inequality $`\chi (R)>0`$ remains true until radial pressure vanishes. Indeed, the hypothesis $`ϵ+p_r>0`$ is crucial for the argument used in , whereas in our case the w.e.c. is satisfied at its borderline, that is $`ϵ+p_r0`$. As pointed out before, these solutions are regular only if they have a de Sitter–like behavior as $`R0^+`$, other regular choices being forbidden by w.e.c.. The junction can even be done between a true de Sitter core and a Schwarzschild spacetime, but this time an intermediate thick shell is needed in the middle, to ensure continuity of radial pressure. An example is given by

(3.5a)
$$m(R)=\{\begin{array}{cc}\frac{4}{3}\pi ϵ_0R^3,\hfill & R[0,kM[\hfill \\ \frac{4}{3}\pi ϵ_0(R2kM)^3+2M,\hfill & R[kM,2kM[\hfill \\ 2M,\hfill & R[2kM,+\mathrm{})\hfill \end{array}$$

where $`k`$ is a constant such that, to ensure continuity, the Schwarzschild mass $`M`$ equals twice the mass of de Sitter $`\frac{4}{3}\pi ϵ_0(kM)^3`$. This choice yields the following energy density function:
(3.5b)
$$ϵ(R)=\{\begin{array}{cc}ϵ_0,\hfill & R[0,kM[\hfill \\ ϵ_0\left(1\frac{2kM}{R}\right)^2,\hfill & R[kM,2kM[\hfill \\ 0,\hfill & R[2kM,+\mathrm{}).\hfill \end{array}$$

The limitation $`k1`$ gives a relation between the parameters of de Sitter and Schwarzschild solutions (i.e. $`ϵ_0M^2\frac{3}{8\pi }`$) to be satisfied in order to obtain a regular black hole interior model with a single Cauchy horizon.

Regular Schwarzschild black holes have been studied by several authors , and explicit solutions have been build up, with a number of Cauchy horizons even greater than one before energy density vanishes (see  for an example).

Of course, matching with Schwarzschild spacetime may be also performed if $`m(R)`$ $`R^\alpha `$ with $`\alpha <3`$, this operation resulting in replacing Schwarzschild black hole with a non flat core still possessing a singularity at $`R=0`$. As an example, if $`m(R)`$ is a second order polynomial in $`R`$,

(3.6)
$$m(R)=kR^2+4kMR,R[0,2M],$$

the junction is done at $`R=2M`$, and the condition $`k\frac{1}{4M}`$ ensures that this matching is made in the non static region. Therefore, in this case there are no Cauchy horizons.

We also notice that the above constructions are made supposing continuity of both the metric and radial pressure $`p_r`$. For instance, if the core is a Schwarzschild solution, we cannot perform a matching with one of our solution at $`R_0>0`$, since continuity of $`p_r`$ is lost at $`R_0`$. Anyway, if discontinuities of $`p_r`$ at junctions are allowed, other spacetimes can be built. Indeed, matching between two variations of mass of Schwarzschild spacetime – generated by an ”inside–shell” mass $`m(R)`$ and an ”outside–shell” mass $`M(R)`$ – can be performed at $`R_0`$ such that, for instance, both $`\chi _m(R_0)`$ and $`\chi _M(R_0)`$ are positive. In this case it must be checked that the surface stress–energy tensor (i.e. related to the surface $`𝒮=\{R=R_0\}`$) obeys weak energy condition. Following the method in  the surface energy density and pressure are found to be

(3.7a) $`ϵ_𝒮={\displaystyle \frac{1}{4\pi R_0}}\left(\sqrt{\chi _m(R_0)}\sqrt{\chi _M(R_0)}\right),`$
(3.7b) $`p_𝒮={\displaystyle \frac{1}{8\pi R_0^2}}\left[{\displaystyle \frac{(1M^{}(R_0))R_0M(R_0)}{\sqrt{\chi _M(R_0)}}}{\displaystyle \frac{(1m^{}(R_0))R_0m(R_0)}{\sqrt{\chi _m(R_0)}}}\right],`$
and for the solution to be physically acceptable it must be $`ϵ_𝒮0`$ (that is, $`M(R_0)m(R_0)`$) and $`ϵ_𝒮+p_𝒮0`$. A particular case of a Schwarzschild region surrounded by a massive shell is studied in , where the region outside shell can also be viewed as a Schwarzschild solution, under a suitable coordinate change.


text:
# Detailed Spectroscopic Analysis of SN 1987A: The Distance to the LMC using the SEAM method

## 1 Introduction

Observations of SN 1987A produced extensive broad-band photometry and excellent spectroscopic coverage. This motivated the development of models to predict its spectral and photometric evolution. In our analysis here we make use of the optical data obtained at the Cerro Tololo Inter-American Observatory (CTIO) (Phillips et al., 1988) and the extensive International Ultraviolet Explorer (IUE) data which were obtained, re-reduced and analyzed by Pun et al. (1995).

While SN 1987A led to the growth of observational data, so too did it lead to an improvement in theoretical models. Höflich (1991) modeled the early phase of SN 1987A, using a pure hydrogen atmosphere and accounting for Non-Local Thermodynamic Equilibrium (NLTE). Eastman & Kirshner (1989) modeled the first ten days of the SN 1987A explosion. H I and He I were modeled in NLTE, while metal lines were treated as pure scattering LTE lines. The pure H/He non-LTE model of Takeda (1991) resulted in weak Balmer lines, as did a model by Hauschildt & Ensman (1994), which used an earlier version of PHOENIX with H I, He I, Mg II, and Ca II in NLTE and metal line blanketing in LTE. These two groups suggested that differences between the theoretical and actual density distributions might be responsible for the discrepancy between the synthetic and observed Balmer lines.

We use the light curve model of Blinnikov and collaborators calculated using the STELLA software package (Blinnikov et al., 1998; Blinnikov, 1999; Blinnikov et al., 2000). This model simulates the light curve for up to six months after the explosion, and includes allowance for time-dependent, multi-group radiation hydrodynamics, monochromatic scattering effects, and the effects of spectral lines on the opacity (Blinnikov et al., 2000). The light curve model was based on the stellar evolution models of Saio, Nomoto, & Kato (1988), Saio, Kato, & Nomoto (1988), and Nomoto & Hashimoto (1988). The procedure for synthesizing SN 1987A spectra from Blinnikov et al.’s model and the results are described in § 2. In § 3 we compare our synthetic spectra to observations. § 4 compares our synthetic light curves with observations and we derive a distance to SN 1987A. In a final section we discuss the implications of our results for the use of SNe II as distance indicators.

## 2 Models

### 2.1 PHOENIX

PHOENIX is a general radiation transfer code that computes temperature, opacity, and level populations for each of typically 50–100 radial zones in a moving stellar envelope. PHOENIX solves the full NLTE rate equations and calculates level populations for a multitude of different atomic species in LTE or NLTE (Hauschildt & Baron, 1999; Short et al., 1999). For this study, the following species were calculated in NLTE: H I(30/435), He I(19/37), He II(10/45), O I(36/66), Ne I(26/37), Na I(53/142), Mg II(72/340), Si II(93/436), S II(84/444), Ca II(87/455), Fe I(494/6903), Fe II(617/13675), and Fe III(566/9721), where e.g. H I(30/435) indicates that our model atom for H I contains 30 levels and 435 permitted transitions treated in full NLTE.

The Blinnikov et al. light curve model consists of an ejecta envelope divided into 300 zones, with compositions defined for each layer (Blinnikov, 1999; Blinnikov et al., 2000). The outer layers of this model were re-zoned onto the 75–100 zone grid that is used by PHOENIX, preserving the compositions in velocity-space. The radius, density, and compositions were taken from the Blinnikov et al. model and interpolated onto the new grid. The inner boundary condition was taken to be diffusive, which is well justified in these calculations. The total bolometric luminosity provides the outer boundary condition and we generally take the value from the results of Blinnikov et al. (2000). However since the bolometric luminosity is a parameter, we study the effects of varying it to improve the fit to the observed spectra. After about 18 days from the explosion we find that the values predicted by the light curve calculations need to be adjusted in order for the synthetic spectra to match the observations. We adjust the bolometric luminosity to find the best fit to the observed spectra, we do not attempt to use the photometry as a fitting parameter in any way. All synthetic photometry is an output of our calculations once we have calculated the synthetic spectrum. The temperature is determined self-consistently from the generalized condition of radiative equilibrium. We use the initial temperature from Blinnikov et al. as our initial guess. The gamma-ray energy deposition rate, which is the energy in gamma-rays deposited in the matter in the envelope per unit mass per unit time, is calculated for each new zone. We calculate the gamma-ray deposition function using our self-consistent spherically symmetric radiative transfer code, assuming a gray opacity for gamma-rays of $`0.06`$ cm2 g-1 (Swartz et al., 1995; Colgate et al., 1980). Our calculational grid was chosen so as to maintain both composition and density profiles with adequate resolution.

## 3 Synthetic Spectra

Figures 13 show the spectra resulting from PHOENIX computations on days 1.36, 2.67, and 3.59, respectively, of SN 1987A. The fits between the IUE spectra and the synthetic UV spectra are very good, an indicator of the importance of NLTE effects in the supernova envelope (Baron et al., 1996b). Blinnikov et al. (2000) noted that the predicted UV flux 5–10 days after shock breakout was much greater than the actual flux when computed in LTE, but this also may have been due to not enough line blanketing in these initial calculations (E. Sorokina et al., in preparation). We included over 1 million lines in all our calculations and they are dynamically selected to be the most important lines for the conditions, thus line blanketing is well accounted for in these calculations. For days 1.36 and 2.67, the agreement between the synthetic and observed optical spectra is also reasonably good. In day 1.36 the blueshifted absorption of the Mg h+k line does not extend far enough to the blue. The observed feature extends to 43,000 km s-1, however the hydro model only extends to 33,500 km s-1. Since Mg h+k is a resonance line it is optically thick all the way to the surface and the blue absorption is sensitive to the outermost layers of the model. Similar results were found by Hauschildt & Ensman (1994) using a different hydro model and an earlier version of PHOENIX.

By day 3.59 a discrepancy in the strengths of the Balmer lines begins to appear in the synthetic spectra, most notably H$`\alpha `$, even when NLTE effects are included. Figure 4 shows the synthetic optical/near-UV spectrum for day 4.52 compared with the observations (Phillips et al., 1988; Pun et al., 1995), and here the Balmer lines are clearly much weaker in the synthetic spectrum than they are in the observed spectrum. As we showed in Mitchell et al. (2001), we can improve the fit to the Balmer lines significantly by increasing the gamma-ray deposition in the outer layers. The result is displayed in Figure 5 where we have replaced the self-consistent gamma-ray deposition function with one that would be obtained from a uniform nickel mass fraction of $`\mathrm{X}_{\mathrm{Ni}}=1.0\times 10^3`$ in the envelope. The fit between the predicted lines and the actual lines is much better in the optical, but the UV is no longer fit as well. This shows that the true nickel deposition is more complicated and we have overestimated the nickel deposition in the very outer layers where the line-blanketed UV is formed. We have chosen a simple parameterization to enhance the gamma-ray deposition, i.e., the gamma-ray deposition follows the density of the model. It is certainly possible (and even plausible) that there exists a spherically symmetric gamma-ray deposition function that fits both the optical and the UV, however the parameter space is large and it is beyond the scope of this work to find “the” correct gamma-ray deposition for each epoch. Nevertheless our results are a strong indication for enhanced nickel mixing in SN 1987A. As noted by Mitchell et al. (2001) the enhanced nickel mixing improves both the strength of H$`\alpha `$ and the velocity (as can be seen by comparing Figs. 5 and 4). Utrobin & Chugai (2002) find that the velocity and strength of H$`\alpha `$ can be well fit by taking into account the effect time dependence of the hydrogen non-thermal excitation and ionization (the ionization freeze-out) up to about 30 days, without invoking large 56Ni mixing.

Figure 6 shows the synthetic spectrum compared to the optical and IUE spectra on day 8. We have assumed uniform nickel mixing of $`\mathrm{X}_{\mathrm{Ni}}=1.0\times 10^4`$. Comparing this to the day 4 spectrum, we can see that there is not quite enough nickel mixing to produce the Balmer lines, but the UV won’t allow for more mixing in the simple parameterization we employ. Figure 7 shows the day 10 spectrum with no additional nickel mixing and again the Balmer lines are far too weak, and only the overall color in the optical is reproduced. Figure 8 shows that nickel mixing improves the blue and the strength of the Balmer lines, however the overall fit is not better. Figure 9 shows the day 14 spectrum, with nickel mixing assumed. The overall fit is not bad in the optical, but again with the too simple deposition function the UV fit is poor. The synthetic spectrum shows a weak feature that has been identified as Ba II (see Williams, 1987; Jeffery & Branch, 1990), but since s-process enhancements have not been included in the stellar evolution model we use we don’t expect to reproduce the observed strength of this or other s-process elements. The Na I feature is far too weak in the synthetic spectrum.

Figure 10 shows the synthetic and observed spectra of day 19, where we have both assumed uniform mixing with $`\mathrm{X}_{\mathrm{Ni}}=1.0\times 10^4`$, and we have increased the bolometric luminosity (a boundary condition for our modeling) of the model by 26% over that obtained from the light curve calculations. We find that if we simply use the bolometric luminosity of Blinnikov et al. (2000) we can not reproduce the observed spectra, but the variations we employ are well within the errors of Blinnikov et al. (2000) (see their figure 15). We have performed a series of calculations where the bolometric luminosity is varied and only present the best fits here. The best fits are obtained using the spectra only, no attempt was made to use photometric colors or any photometric data. The synthetic photometry is an *output* of our calculations. Except for a significant over-prediction of the flux just redward of Ca H+K, the fit is not bad. The Balmer lines are mostly too weak, but the Na D feature is significantly stronger than just 4 days earlier, even though we have increased the target luminosity of the model over that obtained from the light-curve calculation. Figure 11 shows the spectra at day 24, where again we assume constant nickel mixing and have increased the bolometric luminosity by 31%. Figure 12 shows the spectra at day 31 and the fit from the near UV to the optical is quite good. The flux is too low compared to that in the *IUE* SWP band. Even without additional nickel mixing the H$`\alpha `$ line is reasonably reproduced. Figure 13 shows the optical spectrum at day 58 and again the H$`\alpha `$ line is well reproduced. The structure of the Balmer lines is in fact dependent not only on the amount of nickel mixing, but also on the details of the density structure as the pseudo-photosphere recedes in mass. We have performed a series of test calculations to verify this (P. Nugent et al., in preparation). In the observations H$`\beta `$ is blended with an Fe II line just to the red, in the models the lines are unblended and H$`\beta `$ appears too strong. This may be due to the nickel mixing in the supernova, which leads to stronger blending of the iron than is present in the model.

Figure 14 shows that the LWP portion of the *IUE* spectrum is well-fit, however the SWP is not. Figure 15 shows the day 81 spectra with an increase in the bolometric luminosity of 17.5% over the light-curve value.

We have found that the calculated bolometric luminosity in the lightcurve models is too low at some epochs. We also find that the assumed mixing of nickel is also too small, particularly at early times. Our simple enhanced gamma-ray deposition (that the gamma-ray deposition follows the density) does not contain enough parameters to fit both the blue and red portions of the observed spectra. We note that the synthetic spectra are sensitive to the *rate* of gamma-ray deposition as a function of position in the ejecta. At early times almost the entire envelope is optically thick to gamma-rays and therefore we can estimate the total mass of nickel in the outer layers (Mitchell et al., 2001). At later times gamma-rays are transported from the inner layers into the outer shell, and gamma-rays produced in the outermost layers escape, therefore our simple gamma-ray deposition prescription cannot be used to infer a nickel mass fraction as a function of velocity, since gamma-ray deposition is non-local.

Even though it is likely that there exists a spherically symmetric deposition function that will better reproduce the observed spectra, it is clear that the distribution of nickel in SN 1987A is non-spherical. Chugai (1991) found evidence from the He I $`\lambda 10830`$ line for a “jetlike” ejection of 56Ni. He found that he could fit the profile with $`10^4`$ $`\mathrm{M}_{}`$ of 56Ni above 5000 km s-1 and that it was consistent with $`10^310^5`$ $`\mathrm{M}_{}`$ of nickel needed to fit the emission spectrum of the hydrogen lines. This is also consistent with our results. Haas et al. (1990) observed the \[Fe II\] 17.94 and 25.99 $`\mu `$m lines and found clear asymmetry extending to about 3500 km s-1, these results were modeled by Nagataki (2000) using 2-D hydro and confirm the existence of an asymmetric distribution of nickel.

## 4 Synthetic Light Curves

### 4.1 SEAM Distance

We have developed an improved version of the “Expanding Photosphere Method” (EPM): the Spectral-fitting Expanding Atmosphere Method (SEAM). The EPM method (Baade, 1926; Branch & Patchett, 1973; Kirshner & Kwan, 1974) calculates distances to supernovae by assuming that the bolometric luminosity can be approximated as

$$L\zeta ^24\pi R^2\sigma T^4,$$

where $`\zeta `$ is the dilution factor (Hershkowitz et al., 1986a, b; Hershkowitz & Wagoner, 1987), and represents the fact that flux is diluted since the effective temperature differs from that of the photospheric temperature (By energy conservation, $`\zeta `$ can be both less than one or greater than one depending on how T is determined). In supernovae, homology is valid, except during the first instants of the explosion and hence $`R=vt`$. The velocity can be measured from the spectral lines, and $`T`$ is obtained from photometric colors. The time since explosion can be determined by demanding self-consistency among epochs or from the light curve (Eastman & Kirshner, 1989; Schmidt et al., 1992; Eastman et al., 1996; Hamuy et al., 2001; Leonard et al., 2002). EPM clearly suffers from the limitation that supernovae are not blackbodies and spectral line features affect the value of the colors determined observationally and the theoretical colors that must be determined from the Planck function to obtain a distance. It is also unclear just which spectral lines should be used to determine the velocity that is needed to ascertain the radius (Hamuy et al., 2001; Leonard et al., 2002). The other obvious limitation is that the method depends sensitively on the value of the dilution factors which are assumed to be a monotonic function of the color temperature.

The SEAM method (Baron et al., 1993, 1994, 1995a, 1995b, 1996a) retains the assumption of spherical symmetry, and homology $`R=vt`$, but uses detailed spectral models to calculate a synthetic spectrum. Then (modulo extinction, see below) since we know the entire spectral energy distribution, synthetic photometry can be used to calculate a distance modulus in every photometric band. No blackbody assumption is required, no dilution factors are needed, and the velocity is determined by the fit to all the spectral lines, not just some pre-selected set whose identities may not be certain. A time series constrains the time of explosion.

Here, we obtain the radius using the results of the hydrodynamical model and each of our models predicts the total emitted spectral energy distribution (SED) at a particular time. We can convolve our theoretical SED with any filter function and calculate the absolute magnitude in that band. If we denote the absolute magnitude in band $`X`$ as $`M_X`$ then the distance modulus in a particular band is $`\mu _X`$. For every band at every time we can calculate a value of $`\mu _X`$ by the relation

$`\mu _X`$ $`=`$ $`m_XA_XM_X`$
$`=`$ $`m_X(M_X+A_X),`$
where $`m_X`$ is the apparent magnitude in band $`X`$, and $`A_X`$ is the extinction in band $`X`$. We obtain $`A_X`$ by applying the reddening law of Cardelli et al. (1989) with $`E(BV)=0.16`$ (Lundqvist & Fransson, 1996) and $`R=3.1`$ to the theoretical SED. This procedure defines the SEAM method. No assumption about the shape of the SED is made since it is determined by the fit to the observed spectrum.

The distance to the LMC is of great importance in establishing the cosmological distance ladder. The value of the Hubble constant depends crucially on the calibration of Cepheids (Tonry, 2001), which is dependent on the distance to the LMC. Recent reviews of the distance to the LMC are in Gibson (2002) and Feast (2001). We focus here specifically on the distance to SN 1987A, but note that the distance to the center of the LMC is about $`\mu =0.03\pm 0.03`$ nearer than to the supernova (McCall, 1993) and that a recent determination to the eclipsing binary HV 2274 gives $`\mu =18.46\pm 0.06`$ (Groenewegen & Salaris, 2001).

Figure 16 shows the results for the bands $`UBVRI`$. The observed photometry was taken from Catchpole et al. (1987) and Menzies et al. (1987) for the earliest times ($`t<4`$ days), and from Hamuy & Suntzeff (1990) for all other epochs. The $`U`$ band as expected has a lot of scatter since it is quite sensitive to the assumed nickel distribution. The mean distance is $`\mu =18.6\pm 0.3`$ using all the bands. The quoted error here and below is just the 1-$`\sigma `$ error, assuming the errors are Gaussian random errors. If we limit ourselves to $`VRI`$ then we obtain $`\mu =18.5\pm 0.2`$.

Table 1 reveals the surprising result that the large scatter is not limited to $`U`$ and $`B`$ but that the $`I`$ band shows more scatter than does $`R`$. Conventional wisdom dictates that $`R`$ should not be used for distance determinations because it includes H$`\alpha `$, however our results don’t bear this out. Hamuy et al. (2001) found that synthetic $`I`$ magnitudes on *observed* spectra had errors at the level of $`0.1`$ mag, which may explain, e.g, the large discrepancy in the value of $`\mu _I`$ on day 81 in Figure 16.

So far we have used all the model spectra presented in this paper in determining a distance. The purpose of this is to scrupulously convince the reader that we have not simply chosen models which provide us with a desired answer. However, the spirit of the SEAM method is to use *only* the best fit models, i.e. to use the full information of the observed SED to determine which models to include in a distance determination. If we restrict our distances to just the earliest times ($`t<5`$ days) where the spectra are well fit, we find $`\mu =18.44\pm 0.08`$ and if we limit our distances to just the best fit spectra ($`t<5`$ days, 31 days, 58 days) we find $`\mu =18.5\pm 0.2`$ using $`UBVRI`$, and $`\mu =18.46\pm 0.12`$ if we restrict ourselves to $`BVRI`$. Figure 17 shows the results for the SEAM distance, where we require a good fit. This result shows the advantage of the SEAM method when the quality of the fit SED to the data is evaluated to calculate distances. This should tend to reduce systematic effects when applying SEAM since the sign of the error varies with both epoch and photometric band. As our final value for $`\mu `$ we shall adopt the conservative $`\mu =18.5\pm 0.2`$, but the least scatter over the widest wavelength range leads to the slightly lower value of $`\mu =18.46\pm 0.12`$. Table 2 compares determinations of the distance to SN 1987A using both EPM and the ring. Our result is consistent and competitive with even the error bars of the nearly geometrical ring determination.

While the neutrino signal determines the time of the explosion of SN 1987A very accurately, it is a useful exercise to examine how the assumption of the neutrino derived time of explosion affects our models. From a $`\chi ^2`$ analysis, the neutrino signal time is the “best” explosion time in that it minimizes the variance in the distance modulus, and from the SEAM analysis alone the time is known to 68% confidence ($`\mathrm{\Delta }\chi ^2=1`$) within $`\delta t=\pm 0.9`$ days, where the confidence intervals are symmetric about the time of the neutrino signal. Thus, the SEAM method is able to determine both the explosion date and the distance with good accuracy, although our quantitative uncertainty in the explosion time for SN 1987A likely depends on having an extremely well sampled early light curve.

### 4.2 Light Curves

Using our derived value of the distance modulus $`\mu =18.5`$, we can transform our synthetic (reddened) absolute magnitudes into (reddened) apparent magnitudes. Figure 18 compares our bolometric light curve to the results of Suntzeff & Bouchet (1990) and Bouchet et al. (1991). The early times agree quite well, although the agreement is poor for days 8-24 and then improves significantly for days $`>24`$. This just reflects the fact that the hydrodynamical model that we have used is not *the* model for SN 1987A. It is clear that days 8-24 should have a higher bolometric luminosity in order to match the observations, however we have already increased the bolometric luminosity at these times above that obtained by Blinnikov et al. (2000) and larger luminosities than we have used lead to a significantly poorer fit. Thus, there is not a one-parameter correspondence between the fit of the observed spectra to the luminosity.

Figures 1923 show the optical $`UBVRI`$ light curves compared to the observations (Hamuy & Suntzeff, 1990). As we found in §4.1, $`V`$ and $`R`$ are most reliable and $`UBI`$ are less reliable, in fact $`B`$ fails completely. The failure in $`B`$ is disheartening, but shows that the $`B`$ band is highly sensitive to the assumed nickel distribution. Our preliminary calculations show that early $`B`$ band is extremely well fit in normal SNe IIP, which we believe to be better candidates for the SEAM method than are 87A-like events. The fact that $`I`$ is not consistent is largely due to the lack of observed spectra in the range 7000–9000 Å. If we had had spectra to fit this region we would have altered our outer boundary condition and probably achieved higher accuracy. The fact that $`R`$ (which includes H$`\alpha `$) is reasonably well reproduced is interesting. Figure 24 shows that the $`I`$-band only covers the absorption component of the Ca IR-triplet and is dropping in sensitivity in the emission component. The two synthetic spectra differ only in the total amount of nickel which clearly affects the shape of the continuum in addition to that of the lines. The synthetic $`R`$ and $`I`$ both differ by 0.3 mag between the two spectra, but $`RI`$ differs only by 0.15 mag, so that EPM analyses that use $`VI`$ to determine the color temperature might be less reliable than using $`RI`$.

As we have already noted the bolometric light curve would be improved if the models were brighter in the 8-24 day epochs in particular. However, the $`U`$ and $`B`$ band light curves are too bright in this region. Thus, the density structure of the models must also be altered in order to match the observed spectra. In particular, since the blue part of the spectrum forms farthest out (due to line-blanketing from Fe II lines) it is possible that the poor fit in $`U`$ and $`B`$ is due to the atmosphere not extending to high enough velocities, thus there is not enough line blanketing. Line blanketing leads to both back warming and line cooling and so it may be possible to reduce the $`U`$ and $`B`$ magnitudes, while at the same time increasing the total bolometric luminosity. It is also possible that asymmetry is playing a role at these times, but we would expect the effects of asymmetry to increase with time and so we think that asymmetry is unlikely to be the cause here. However, the proper gamma-ray deposition function is certainly needed to reproduce the blue part of the spectrum.

Figures 2527 show the IR ($`JHK`$) light curves where even without guidance from observed spectra we obtain very good agreement. The observations are taken from Bouchet et al. (1989). This result gives us confidence in the value obtained for $`\mu `$ and it indicates that only small changes are needed in order to obtain consistent values of $`I`$. On the other hand it shows that spectroscopy in the $`I`$ band is required in order to obtain consistent results. Figure 28 shows the IR spectrum of the 30.9 day model with the filter and atmospheric transmission functions. The SED is nearly a blackbody with only the hydrogen lines P$`\alpha `$, P$`\beta `$, and P$`\gamma `$ displaying strong P-Cygni profiles.

## 5 Discussion

Our results indicate that the SEAM method has promise as an accurate, competitive distance indicator. It has the advantage over other methods that its results can be directly compared for consistency and that the SED includes important effects due to differential expansion and line blanketing. On the other hand, one would like to see perfect consistency in every band at every epoch. While perfect consistency is a unreachable goal, we expect improved results with early time normal Type IIP supernovae. We have used here only a single hydrodynamical model which is clearly not *the* model of SN 1987A. We believe that the errors can be significantly reduced with a large grid of models at each epoch which we will pursue in future work. It is particularly important to have a large range of density profiles for the lower velocity ejecta in order to correctly reproduce the Balmer lines consistently over time. It is also important for the model to extend to the highest velocity seen in the observed spectra.

On the face it may appear that the SEAM method is less consistent than the EPM method as applied to SN 1999em (Hamuy et al., 2001; Leonard et al., 2002). However, the EPM is more prone to systematic effects than is SEAM, where the systematic errors are more likely to average out. In EPM the “temperature” is obtained from a color temperature from different filter bands. Hamuy et al. (2001) found that the EPM leads to a deviation in the distance obtained with epoch (and they did not use either the $`U`$-band or $`I`$-band data). They also found that $`BV`$ were less reliable than the redder bands, which is consistent with our results. Leonard et al. (2002) found that they obtained different preferred distances when using different filter bands and that determining the velocities that are required in EPM is non-trivial. Leonard et al. (2002) found that even obtaining the rest wavelength of multiplets such as Sc II is subject to uncertainty. In reality, there is probably no single unblended line in the entire spectrum, and the problems of blending become more important for weak lines that are desired to obtain the photospheric radius. The need for detailed synthetic spectral modeling is clear. In fact, the velocity is needed in EPM not only to determine the radius of the supernova, but also to compare with models in order to obtain the dilution factor, $`\zeta `$. Figure 29 shows the results of our values of $`\zeta `$ (obtained using $`BV`$) with that of Eastman et al. (1996). The filled circles are our result, the triangles the power-law model, p6, of Eastman et al. (1996) and the squares the 15 $`\mathrm{M}_{}`$ model, s15, of Eastman et al. (1996). There does appear to be some difference between our values and that of Eastman et al. (1996), particularly for the important color temperature range of $`50008000`$ K. Although Hamuy et al. (2001) found that the power-law model gave more consistent results for SN 1999em than s15, this cannot be physically correct since real models must flatten out at low velocity (or they would be singular).

The fact that the scatter is minimized and the fits are quite good at early times, leads us to believe that early time normal type IIP supernovae are likely the best candidates for SEAM. Since polarization is lowest at early times (Leonard et al., 2001) it is also the time that the assumption of spherical symmetry is likely to be most valid. Thus we believe that future SEAM (and EPM) analyses should concentrate on early time supernovae. It is now up to dedicated supernova searches to find supernovae at early times and *to follow them both photometrically and spectroscopically*.

We thank Mario Hamuy for providing us with IR-filter and transmission functions and for helpful discussions on synthetic photometry and we thank Claes Fransson for helpful discussions. This work was supported in part by NSF grants AST-9731450, NASA grant NAG5-3505, and an IBM SUR grant to the University of Oklahoma; and by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619 to the University of Georgia. SB is supported in Russia by the grant RFBR 02-02-16500. PHH was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. PL acknowledges support from the Royal Swedish Academy and the Swedish Research Council. Some of the calculations presented here were performed at the San Diego Supercomputer Center (SDSC), supported by the NSF, and at the National Energy Research Supercomputer Center (NERSC), supported by the U.S. DOE. We thank both these institutions for a generous allocation of computer time.


text:
# On states, channels and purificationWork partially supported by INTAS grant 00-738.

## Abstract

In this note we introduce purification for a pair $`(\rho ,\mathrm{\Phi }),`$ where $`\rho `$ is a quantum state and $`\mathrm{\Phi }`$ is a channel, which allows in particular a natural extension of the properties of related information quantities (mutual and coherent informations) to the channels with arbitrary input and output spaces.

Given a state $`\rho `$ and a channel $`\mathrm{\Phi }`$ there are three fundamental entropy quantities: the input entropy $`H(\rho )`$, the output entropy $`H(\mathrm{\Phi }[\rho ]),`$ and the entropy exchange $`H(\rho ,\mathrm{\Phi }),`$ the last being defined as $`H((\mathrm{\Phi }\mathrm{Id})[|\psi _\rho ><\psi _\rho ]),`$ where $`|\psi _\rho >`$ is a purification of $`\rho .`$ From these entropies one derives several information quantities of which the mutual information $`I(\rho ,\mathrm{\Phi })=H(\rho )+H(\mathrm{\Phi }[\rho ])H(\rho ,\mathrm{\Phi })`$ and the coherent information $`I_c(\rho ,\mathrm{\Phi })=H(\mathrm{\Phi }[\rho ])H(\rho ,\mathrm{\Phi })`$ are of the main importance. The quantities $`I,I_c`$ have a number of useful properties, the standard proof of which uses a dynamical representation of the channel $`\mathrm{\Phi }`$ via interaction of the system in question with an environment (see e. g. , ). Such a representation presupposes that the system preserves its identity during the interaction and hence the input and the output spaces of the channel are necessarily identical too. However from the viewpoint of information theory this is quite an unnatural restriction. One of the properties which we use as an illustration is the data processing inequality for two channels $`\mathrm{\Phi }_1,\mathrm{\Phi }_2:`$

$$I_c(\rho ,\mathrm{\Phi }_2\mathrm{\Phi }_1)I_c(\rho ,\mathrm{\Phi }_1).$$

It is natural to expect that it holds for arbitrary channels satisfying the only restriction that the output of the first channel is equal to the input of the second. But using dynamical representations restrict us to the case where both channels act in the same space. The present note removes this restriction by introducing a new tool: a purification of the couple $`(\rho ,\mathrm{\Phi }),`$ from which all the three entropy quantities emerge with equal status.

Consider the spectral decomposition

$$\rho =\underset{j=1}{\overset{d}{}}\lambda _j|e_j><e_j|,\lambda _j0,\underset{j=1}{\overset{d}{}}\lambda _j=1$$

in the input Hilbert space $`_{in},`$ and the Kraus decomposition of the channel

$$\mathrm{\Phi }[\sigma ]=\underset{\alpha =1}{\overset{N}{}}A_\alpha \sigma A_\alpha ^{},\underset{\alpha =1}{\overset{N}{}}A_\alpha ^{}A_\alpha =I,$$

where $`A_\alpha :_{in}_{out}.`$ We introduce the three systems $`R,Q,E`$ with the Hilbert spaces defined as follows:

$$_R=l_d^2_{in};_Q=_{out};_E=l_N^2,$$

where $`l_n^2`$ is the standard Hilbert space of $`n`$-dimensional vectors and $``$ denotes unitary equivalence. Consider the tensor product $`=_R_Q_E`$ which can be realized as the space of vectors $`|\psi >=\left[\psi _{j\alpha }\right]`$ with the components $`\psi _{j\alpha }_{out}.`$ An operator $`X`$ in $``$ is represented by a square matrix $`\left[X_{k\beta }^{j\alpha }\right]`$, the elements of which are operators in $`_{out}.`$ Partial traces in $``$ are computed according to the formulas

$$\mathrm{Tr}_RX=\left[\underset{j=1}{\overset{d}{}}X_{j\beta }^{j\alpha }\right];\mathrm{Tr}_QX=\left[\mathrm{Tr}X_{j\beta }^{j\alpha }\right];\mathrm{Tr}_EX=\left[\underset{\alpha =1}{\overset{N}{}}X_{k\alpha }^{j\alpha }\right].$$

Take the unit vector

$$|\psi _{(\rho ,\mathrm{\Phi })}>=\left[\sqrt{\lambda _j}A_\alpha |e_j>\right],$$

and the corresponding pure state $`\mathrm{\Omega }=|\psi _{(\rho ,\mathrm{\Phi })}><\psi _{(\rho ,\mathrm{\Phi })}|,`$ then the partial states are

$$\mathrm{\Omega }_R=\left[\lambda _j\delta _k^j\right]\rho ;\mathrm{\Omega }_Q=\mathrm{\Phi }[\rho ];\mathrm{\Omega }_E=\left[\mathrm{Tr}A_\alpha \rho A_\beta ^{}\right],$$

therefore

$$H(\mathrm{\Omega }_R)=H(\rho );H(\mathrm{\Omega }_Q)=H(\mathrm{\Phi }[\rho ]);H(\mathrm{\Omega }_E)=H(\rho ,\mathrm{\Phi }),$$

where the last equality follows from the facts that $`H(\mathrm{\Omega }_E)=H(\mathrm{\Omega }_{RQ})`$ and

$$\mathrm{\Omega }_{RQ}=[\sqrt{\lambda _j}\sqrt{\lambda _k}\mathrm{\Phi }[|e_j><e_k|]]=(\mathrm{Id}\mathrm{\Phi })[|\psi _\rho ><\psi _\rho |],|\psi _\rho >=[\sqrt{\lambda _j}|e_j>].$$

It also coinsides with the well-known expression for the entropy exchange in the dynamical case, see e. g. .

Let us see how this can be used for the proof of the general data processing inequality. Let $`\left\{A_\alpha \right\},\left\{B_\mu \right\}`$ be the components of the Kraus decompositions for the channels $`\mathrm{\Phi }_1,\mathrm{\Phi }_2.`$ By using the purification $`\mathrm{\Omega }^1=|\psi _{(\rho ,\mathrm{\Phi }_1)}><\psi _{(\rho ,\mathrm{\Phi }_1)}|`$ in $`_R_{Q_1}_{E_1}`$ we obtain

$$I_c(\rho ,\mathrm{\Phi }_1)=H(\mathrm{\Omega }_{RE_1}^1)H(\mathrm{\Omega }_{E_1}^1).$$

For the superposition $`\mathrm{\Phi }_2\mathrm{\Phi }_1`$ we use the purification

$$|\psi _{(\rho ,\mathrm{\Phi }_2\mathrm{\Phi }_1)}>=\left[\sqrt{\lambda _j}B_\mu A_\alpha |e_j>\right]_R_{Q_2}_{E_1}_{E_2},$$

where $`_{E_2}=l_M^2,M`$ being the number of the components in the Kraus decomposition of $`\mathrm{\Phi }_2.`$ If $`\mathrm{\Omega }^{12}`$ is the corresponding pure state, then

$$I_c(\rho ,\mathrm{\Phi }_2\mathrm{\Phi }_1)=H(\mathrm{\Omega }_{RE_1E_2}^{12})H(\mathrm{\Omega }_{E_1E_2}^{12}).$$

The data processing inequality will follow from strong subadditivity if we show that $`\mathrm{\Omega }_{RE_1}^1=\mathrm{\Omega }_{RE_1}^{12}.`$ But

$$\mathrm{\Omega }_{RE_1}^1=\mathrm{Tr}_{Q_1}\mathrm{\Omega }^1=\left[\sqrt{\lambda _j\lambda _k}<\psi _k|A_\beta ^{}A_\alpha |\psi _j>\right],$$

while

$$\mathrm{\Omega }_{RE_1}^{12}=\mathrm{Tr}_{Q_2E_2}\mathrm{\Omega }^{12}=\left[\sqrt{\lambda _j\lambda _k}\underset{\mu =1}{\overset{M}{}}<\psi _k|A_\beta ^{}B_\mu ^{}B_\mu A_\alpha |\psi _j>\right],$$

which is indeed the same.

To see that the data processing inequality is not special, and other properties can be treated in a similar way, let us establish subadditivity of quantum mutual information

$$I(\rho _{12},\mathrm{\Phi }_1\mathrm{\Phi }_2)I(\rho _1,\mathrm{\Phi }_1)+I(\rho _2,\mathrm{\Phi }_2).$$

Let $`\mathrm{\Omega }^{12},\mathrm{\Omega }^1,\mathrm{\Omega }^2`$ be the purifying states for the corresponding systems, constructed according to the recepee above. Then in obviuos notations the inequality becomes

$`H(\mathrm{\Omega }_{R_{12}}^{12})+H(\mathrm{\Omega }_{Q_1Q_2}^{12})H(\mathrm{\Omega }_{E_1E_2}^{12})`$
$`H(\mathrm{\Omega }_{R_1}^1)+H(\mathrm{\Omega }_{Q_1}^1)H(\mathrm{\Omega }_{E_1}^1)`$
$`+H(\mathrm{\Omega }_{R_2}^2)+H(\mathrm{\Omega }_{Q_2}^2)H(\mathrm{\Omega }_{E_2}^2),`$
or, taking into account purifications,

$`H(\mathrm{\Omega }_{Q_1Q_2E_1E_2}^{12})+H(\mathrm{\Omega }_{Q_1Q_2}^{12})H(\mathrm{\Omega }_{E_1E_2}^{12})`$
$`H(\mathrm{\Omega }_{Q_1E_1}^1)+H(\mathrm{\Omega }_{Q_1}^1)H(\mathrm{\Omega }_{E_1}^1)`$
$`+H(\mathrm{\Omega }_{Q_2E_2}^2)+H(\mathrm{\Omega }_{Q_2}^2)H(\mathrm{\Omega }_{E_2}^2).`$
This again will follow from repeated use of strong subadditivity and from subadditivity of $`H(\mathrm{\Omega }_{Q_1Q_2}^{12})`$ if we show that

$$\mathrm{\Omega }_{Q_1E_1}^1=\mathrm{\Omega }_{Q_1E_1}^{12},\mathrm{\Omega }_{Q_2E_2}^2=\mathrm{\Omega }_{Q_2E_2}^{12}.$$

But the first equality just means that

$$\left[\underset{\mu =1}{\overset{M}{}}\mathrm{Tr}_{Q_2}(A_\alpha B_\mu )\rho _{12}(A_\beta B_\mu )^{}\right]=\left[A_\alpha \rho _1A_\beta ^{}\right],$$

and similarly the second.

Let us make some general remarks. Unlike the case of dynamical representation of the channels, where one considers the “in” systems $`R,Q,E`$ and the “out” systems $`R^{},Q^{},E^{}`$ etc., here we have only the “out” systems with $`R`$ being identical to the input of $`Q`$ so there is no sense in “primed” notations. As in other cases, we observe a kind of complementarity between statistical and dynamical aspects of quantum description: our purification is “less physical” in that it does not admit a dynamical interpretation, but it is more relevant to the statistical structure related to the pair $`(\rho ,\mathrm{\Phi })`$ since it does not involve any arbitrariness from outside.


text:
# Chaotic diffusion of particles with finite mass in oscillating convection flows

## Abstract

Deterministic diffusion in temporally oscillating convection is studied for particles with finite mass. The particles are assumed to obey a simple dissipative dynamical system and the particle diffusion is induced by the strange attractor. The diffusion constants are numerically calculated for convection models with free and rigid boundary conditions.
PACS numbers:05.40.-a, 05.45.-a, 47.52.+j

The motion of a fluid particle in the velocity field $`𝐯(𝐱)`$ is determined by the differential equation.

$$\frac{d𝐱}{dt}=𝐯(𝐱,t)$$
(1)
with the initial condition $`𝐱(0)`$. Even for regular velocity fields, the particle motion can become chaotic and it is called the Lagrangian chaos. For an incompressible fluid, the dynamical system (1) is a conservative system. In two dimensions, the evolution equations become

$$\frac{dx}{dt}=\frac{\psi }{z},\frac{dz}{dt}=\frac{\psi }{x},$$
(2)
where $`\psi `$ is the stream function and $`𝐱=(x,z)`$. Equation (2) is a Hamiltonian system. Solomon and Gollub studied experimentally and numerically the particle diffusion in Rayleigh-Bénard convection which oscillates temporally . In the oscillating convection, the velocity field is derived from a stream function $`\psi =A/k\mathrm{sin}\{kx+B\mathrm{sin}(\omega t)\}W(z)`$, where $`k`$ is the wavenumber, $`A`$ is the maximum vertical velocity, and $`B`$ and $`\omega `$ are the amplitude and the frequency of the oscillation and $`W(z)`$ is a function which satisfies the boundary conditions at the top and bottom surfaces ($`z=1`$ and 0).

In this report, we assume that a particle embedded in the fluid has a finite mass. The equation of motion is in general rather complicated . We assume a simple model as

$$\frac{dv_x}{dt}=\gamma \left(\frac{\psi }{z}v_x\right),\frac{dv_z}{dt}=\gamma \left(\frac{\psi }{x}v_z\right),$$
(3)
where $`v_x=dx/dt,v_z=dz/dt`$, and $`\gamma ^1`$ denotes a response time and is proportial to the mass of the particle. For the massless case, $`\gamma `$ becomes infinity and Eq. (2) is recovered. The volume contraction rate of the four dimensional phase space $`(x,v_x,z,v_z)`$ is $`2\gamma `$, and therefore the dynamical system (3) is a dissipative system. We assume the same form of stream function $`\psi =A/\pi \mathrm{sin}\{\pi x+B\mathrm{sin}(2\pi t)\}W(z)`$. For simplicity, we assume $`W(z)=\mathrm{sin}(\pi z)`$ for the convection with free boundary conditions, and $`W(z)=z^2(1z)^2`$ for the convection with rigid boundary conditions. More complicated function $`W(z)`$ is used in the exact linear stability analysis of the Bénard convection with rigid boundary conditions .

The fixed points of the dynamical system (3) at $`B=0`$ are $`(x_0,v_{x0},z_0,v_{z0})=(1/2,0,1/2,0)`$, and $`(x_1,v_{x1},z_1,v_{z1})=(0,0,0,0)`$ and $`(x_2,v_{x2},z_2,v_{z2})=(0,0,1,0)`$. The fixed point $`(x_0,z_0)`$ is an elliptic point for the Hamiltonian system (2) and the orbit near the elliptic point is stable for small perturbation in the Hamiltonian system. It is called the KAM theorem and the chaotic behaviors cannot be expected near the elliptic point. The fixed points $`(x_1,z_1)`$ and $`(x_2,z_2)`$ are saddle points and the chaotic orbits appear near the saddle points for nonzero $`B`$ in the Hamiltonian system.

On the other hand, the fixed point $`(x_0,v_{x0},z_0,v_{z0})`$ for Eq. (3) at $`B=0`$ is an unstable focus for finite $`\gamma `$ and the orbit approaches the heteroclinic orbits which connects the different saddle points. Figure 1(a) displays a trajectory starting from $`(x,v_x,z,v_z)=(0.5001,0,0.5,0)`$ for $`A=1,B=0,\gamma =20`$ and $`W(z)=\mathrm{sin}\pi z`$. As the time evolves, all trajectories approach the saddle connections which connect four saddle points $`(0,0),(1,0),(1,1)`$ and $`(1,0)`$. The numerical simulation was performed with the Runge-Kutta method of timestep $`\mathrm{\Delta }t=0.0005`$. The orbits seem to be attracted rapidly to the saddle points.

For the rigid boundary model with $`W(z)=z^2(1z)^2`$, there are fixed lines at $`z=0`$ and 1. The eigenvalues around the points $`(x,z)=(x,0)`$ and $`(x,1)`$ on the fixed lines are calculated from

$$\frac{d\delta x}{dt}=\delta v_x,\frac{d\delta v_x}{dt}=\gamma (2A/\pi \mathrm{sin}(\pi x)\delta z\delta v_x),\frac{d\delta z}{dt}=\delta v_z,\frac{d\delta v_z}{dt}=\gamma \delta v_z,$$

where $`\delta x,\delta z,\delta v_x,\delta v_z`$ are small deviations from the points on the fixed lines. The eigenvalues are $`0,0,\gamma `$ and $`\gamma `$. It is a singular situation in a viewpoint of dynamical systems. The trajectories are attracted to the saddle connection, however, the flow near the fixed lines is very slow. Figure 1(b) displays the time evolution of $`x(t)`$ starting from $`(x,v_x,z,v_z)=(0.1,0,0.1,0)`$ for $`A=4\pi ,B=0,\gamma =20`$ and $`W(z)=z^2(1z)^2`$. The time evolution of $`x(t)`$ may seem to be like a limit cycle, however, the time axis is plotted with a logarithmic scale. It implies that it takes exponentially longer time for the particle to circulate one convection cell as the particle is close to the saddle connection.

The saddle connection is unstable for general perturbations. The unstable manifolds and stable manifolds starting from saddle points intersect and a strange attractor appears near the saddle points for nonzero $`B`$ . All trajectories are attracted to the strange attractor near the saddle points and the chaotic diffusion is induced.

Figure 2(a) is a projection of the strange attractor into $`([x+1/2],[z+1/2])`$ plane at $`t=n=`$integer for $`A=1,B=0.001,\gamma =20`$ and $`W(z)=\mathrm{sin}\pi z`$, where $`[y]`$ implies the fractional part of $`y`$. The strange attractor around the saddle points $`(x,z)=(0,0),(1,0),(0,1)`$ and $`(1,1)`$ is shifted around $`(1/2,1/2)`$ to clarify the attractor. We have performed several numerical simulations using different initial conditions, however, we have obtained the same strange attractor. The stationary distribution of the particle position is uniform for the Hamiltonian system (2) which corresponds to the massless case, however, the stationary distribution of the particle position is not uniform in the dissipative system (3). Figure 2(b) displays a double logarithmic plot of the distribution $`P(z)`$ of $`z`$ near $`z=0`$ for $`A=1,B=0.001,\gamma =20`$ and $`W(z)=\mathrm{sin}\pi z`$. The particle is attracted to the strange attractor near the saddle points and the stationary distribution seems to obey a power law with exponent about 0.93. Figure 2(c) displays the time evolution of $`x(t)`$ for the same parameters from the initial condition $`(x,v_x,z,v_z)=(0.5001,0,0.5,0)`$. Chaotic diffusion is induced even for the very small perturbation.

There are many ”windows” of stable limit cycle attractors when the parameter $`B`$ is changed. If a limit cycle which runs across convection cells exists and it is stable, $`x(t)`$ increases or decreases with a constant average velocity. The limit cycle corresponds to the accelerator mode in the Hamiltonian system (2) . As $`B`$ is changed, the stable limit cycle disappears at a critical point. However, the velocity correlation is very long near the bifurcation point of the limit cycle. The velocity $`v_x`$ keeps positive or negative values for a long time and the direction of motion changes intermittently. Figure 3 displays such a constantly increasing time evolution of $`x(t)`$ at $`B=0.002`$ and an intermittent time evolution of $`x(t)`$ at $`B=0.002165`$ for $`A=1,\gamma =20`$ and $`W(z)=\mathrm{sin}\pi z`$.

The diffusion constant $`D=lim_\tau \mathrm{}(x(t+\tau )x(t))^2/2\tau `$ is enhanced near the bifurcation point. Figure 3(b) displays the diffusion constant as a function of $`B`$. The average $`\mathrm{}`$ to evaluate the diffusion constant is numerically calculated as double averages of a long time average and the average with respect to 100 different initial conditions. The diffusion constant increases as $`1/(BB_c)^\alpha `$, where $`B_c0.00216`$ and $`\alpha =0.50.6`$.

The chaotic diffusion occurs both for free boundary models and rigid boundary models. The difference between the two boundary conditions is clearly seen, when the perturbation amplitude $`B`$ is very small. Figure 4 displays the averge period for a particle to circulate one convection cell as a function of $`B`$. The period is numerically obtained as the average interval between neighboring two times where $`z(t)=0.5`$ and $`dz/dt>0`$. Figure 4(a) shows the results for $`A=1,\gamma =20,W(z)=\mathrm{sin}\pi z`$. The average period $`T`$ is approximately $`T4.1\mathrm{log}_{10}B`$. The logarithmic dependence is characteristic of the saddle connection. The period increases very slowly as $`B`$ is decreased. Figure 4(b) shows the results for $`A=4\pi ,\gamma =20`$ and $`W(z)=z^2(1z)^2`$.

The average period increases approximately by a power law: $`TB^{1/2}`$ as $`B`$ is decreased for the rigid boundary model. It may be due to the existence of the fixed lines for the rigid boundary conditions. The average period $`T`$ plays a role of a unit time for the chaotic diffusion. For example, if the direction of motion is assumed to be randomly changed after each circulation of convection cells, the total steps $`N`$ of random changes of direction during the time interval $`\tau `$ is estimated as $`N\tau /T`$, therefore, the diffusion constant is estimated as $`D1/T`$. Figure 5(a) displays numerically obtained diffusion constants at several $`B`$’s for $`A=1,\gamma =20,W(z)=\mathrm{sin}\pi z`$. The diffusion constants depend on $`B`$ in a complicated manner. It is partly due to the ”windows” structures by the stable limit cycles. For example, the diffusion constant $`D`$ is zero at $`B=3\times 10^5`$, where there exists a stable limit cycle which circulates in a single convection cell. However, the diffusion constant has a tendency to decrease very slowly as $`B`$ is decreased. Figure 5(b) displays numerically obtained diffusion constants at several $`B`$’s for $`A=4\pi ,\gamma =20`$ and $`W(z)=z^2(1z)^2`$. The diffusion constant decreases with a power law: $`DB^{1/2}`$ as $`B`$ is decreased. The numerical results are consistent with the above rough argument.

To summarize, we have numerically studied the chaotic diffusion of particles with finite mass in oscillating convection. The diffusion is induced by the strange attractor near the saddle points. In a sence, the chaotic diffusion is enhanced in such dissipative systems, since there exist no stable oribits near the center of the convection cell and the orbits are attracted toward the saddle points in contrast to the Hamiltonian system. The strange attractor exists even for very small perturbations, since the attractor of the nonperturbed system is a singular one, that is, the saddle-connection. Special enhancement of the diffusion constant occurs near the bifurcation points of limit cycles corresponding to the accelerator modes in the Hamiltonian system. The behaviors of the chaotic diffusion at very weak $`B`$ depend on the boundary conditions. For the rigid boundary conditions, the motion of the particle becomes very slow near the top and bottom boundaries, since the flow velocity is zero at the boundaries. Therefore, the average period circulating around one cell becomes very large for small $`B`$ and the diffusion constant decreases more rapidly than the case of the free boundary conditions. However, it is left to future study to explain the $`B`$ dependence theoretically.

We would like to thank Professor S.Kai, Dr.Y.Hidaka and Mr.K.Tamura for valuable discussions.