text:
# Field Theory of Polymers: Escaping the Sign Problem

## I Introduction

Field theories have proved to be a useful theoretical framework for studying the equilibrium and time-dependent properties of a wide variety of complex fluid systems. In the 1960s, Edwards introduced field-theoretic methods to the field of polymer physics and applied these techniques fruitfully to examine the universal aspects of chain conformational statistics and thermodynamic properties of polymer solutions edwards65 . Perturbation methods, mean-field approximations, and variational schemes were introduced for tackling important problems such as the excluded volume effect and addressing collective phenomena, such as screening in semi-dilute solutions. These analytical methods were considerably enhanced in the 1970s and 1980s with the advent of renormalization group theory, the introduction of scaling concepts, and the discovery of the connection of polymer statistical mechanics with the $`n`$-vector field theory model for $`n0`$ degennes .

Unfortunately, most of these analytical techniques are difficult to apply in the important case of *inhomogeneous* polymers, which include multicomponent polymer alloys, block and graft copolymers, nano-composites, thin films, emulsions, and suspensions. Appropriate field theory models for such systems are straightforward to construct, but the analytical tools available for examining the properties of the models are rather limited. One powerful tool is self-consistent field theory (SCFT), which amounts to a saddle point approximation to the partition function of the field theory helfand75 ; schmid98 . However, the SCFT field equations determining the “mean fields” are nonlinear, nonlocal, and defy analytical solution. Certain limiting cases of SCFT, such as weak inhomogeneities or fields that vary slowly over the size of a polymer, and for simple geometries, are analytically tractable. In such special circumstances, however, the SCFT solutions still reflect the mean-field assumption intrinsic to the theory, which further restricts their applicability. In particular, such solutions are primarily limited to the description of concentrated solutions or melts of high molecular weight polymers. Even in the case of dense molten blends or copolymers, SCFT can fail in situations where soft modes appear, such as near critical points, unbinding transitions, or in micellar or microemulsion phases bates97 , since thermally-induced field fluctuations can then be significant.

In recent years, it has become apparent that field theory models of polymers and complex fluids, besides serving as the platform for analytical calculations, can also provide a flexible framework for numerical simulations ghfrev . Building on significant advances in numerical methods for solving the mean-field SCFT equations fraaije97 ; matsen94 , techniques have been developed for direct numerical sampling of field fluctuations ganesan01 ; ghfrev . These methods open many exciting possibilities for carrying out “field-theoretic simulations” of a wide variety of inhomogeneous polymer systems. Aside from the numerical errors associated with representing the fields and the usual statistical errors intrinsic to any computer simulation technique, field-theoretic simulations provide a route to studying the exact properties of a polymer fluid model in the absence of any simplifying approximations. Such field-based simulations also facilitate systematic coarse-graining, making it possible to equilibrate much larger systems with a broader range of structural scales than in conventional atomistic computer simulations. At the present time the methods are limited to *equilibrium*  structure and properties, although research is underway to extend them to nonequilibrium situations.

A difficulty with the current strategies for field-based simulations of polymeric fluids is that they rely on the introduction of Hubbard-Stratonovich “chemical potential” fields, first employed for polymers by Edwards edwards65 , to decouple the non-bonded interactions among monomers. This allows one to explicitly trace out the monomer coordinates in favor of the introduced fields, generating a field theory model for a polymer solution or melt. Unfortunately, this description leads to an effective Hamiltonian $`𝒢(\varphi )`$, expressing the energy of a chemical potential field configuration $`\varphi (r)`$, that is *complex*, rather than real. As a consequence, the Boltzmann factor $`\mathrm{exp}(𝒢)`$ is not positive definite, so conventional Monte Carlo sampling of the partition function is problematic lin86 . This difficulty in sampling field configurations is a special case of a more general “sign problem”, which involves the development of techniques for efficient numerical evaluation of high dimensional integrals with non-positive definite integrands. While promising new methods have appeared for battling the sign problem baeurle02 ; ghf03 , the problem is by no means “solved”. Oscillations in the phase factor $`\mathrm{exp}(i\mathrm{Im}𝒢)`$ severely limit the convergence of field-based simulations of polymeric fluids.

In the present paper, we describe an alternative formulation of polymer field theory that essentially side-steps the sign problem. Rather than retaining just a Hubbard-Stratonovich chemical potential field for each monomer species, we also retain a second density field. This two field approach is not new matsen94 ; hongnoolandi81 , but rather is a usual starting point for deriving SCFT. However, we show in the present paper that at a fixed value of the density field, the functional integral over the chemical potential field can be performed in the *saddle point approximation, and that this approximation is exact in the thermodynamic limit of large systems.* This we believe to be a new result and one that has profound consequences for field-theoretic computer simulations of polymeric fluids.

## II Basic Formalism: Polymer Solutions

### II.1 Transformation to a Field Theory

Consider a solution of $`M`$ homopolymer chains in a volume $`V`$, each chain composed of $`N`$ monomers of size $`a`$. The solvent is not considered explicitly, the temperature is denoted by $`T`$, and $`\beta =1/T`$. The formalism presented below can be easily generalized to blends of different chains and to melts or solutions of copolymers as will be seen later.

We assume that the monomers interact through a 2-body interaction potential (of mean force) denoted by $`v(r)`$. In the celebrated Edwards model, the interaction $`v(r)`$ is modelled as a Dirac delta function, $`\delta (r)`$, but in general, the interaction $`v`$ can represent both solvent effects, electrostatic effects or any other effective monomer-monomer interaction. The canonical partition function of the system reads

$$𝒵=\frac{1}{M!}\underset{k=1}{\overset{M}{}}𝒟r_k\mathrm{exp}\left(\frac{3}{2a^2}\underset{k=1}{\overset{M}{}}_0^N𝑑s\left(\frac{dr_k}{ds}\right)^2\frac{\beta }{2}\underset{k,l=1}{\overset{M}{}}_0^N𝑑s_0^N𝑑s^{}v(r_k(s)r_l(s^{}))\right)$$
(1)
where $`𝒟r_k`$ denotes a path integral over all conformations of the $`k`$th polymer chain.

Next, we introduce a monomer (number) density field $`\rho (r)`$ as a new integration field. The microscopic expression for the monomer density, $`\widehat{\rho }(r)`$, is given by

$$\widehat{\rho }(r)=\underset{k=1}{\overset{M}{}}_0^N𝑑s\delta (rr_k(s))$$
(2)
If we constrain $`\rho (r)`$ to $`\widehat{\rho }(r)`$ at each point $`r`$ of space by the use of a $`\delta `$-function, the partition function of the polymer system can be written as

$$𝒵=𝒟\rho \mathrm{exp}\left(\frac{\beta }{2}𝑑r𝑑r^{}\rho (r)v(rr^{})\rho (r^{})\right)Z(\rho )$$
(3)
where the function $`Z(\rho )`$ is defined by

$$Z(\rho )=\frac{1}{M!}\underset{k=1}{\overset{M}{}}𝒟r_k\mathrm{exp}\left(\frac{3}{2a^2}\underset{k=1}{\overset{M}{}}_0^N𝑑s\left(\frac{dr_k}{ds}\right)^2\right)\underset{r}{}\delta \left(\rho (r)\widehat{\rho }(r)\right)$$
(4)
and $`𝒟\rho `$ denotes a functional integral over the real, scalar density field $`\rho (r)`$. The object $`Z(\rho )`$ is the partition function of a set of $`M`$ continuous Gaussian chains of $`N`$ monomers, constrained in such a way that the total monomer density is fixed to $`\rho (r)`$ at any point $`r`$ of space. This partition function evidently counts the total number of chain conformations that are consistent with the density profile $`\rho (r)`$. We thus write

$$Z(\rho )=\mathrm{exp}[S(\rho )]$$
(5)
and identify $`S(\rho )`$ with the entropy of a system of $`M`$ chains constrained to have a density profile $`\rho (r)`$. From its very definition, it is clear that $`Z(\rho )`$ is real and positive for any real $`\rho (r)`$. In addition, we expect $`Z(\rho )`$ to scale with $`V`$, $`N`$, and $`M`$ as $`(Vg^N)^M/M!`$, where $`V`$ is the volume of the system and $`g`$ is some positive, real constant.

At this stage, it is clear that we have expressed the partition function of the polymeric system in terms of a functional integral over a real field with a positive definite Boltzmann weight. We have yet to show how to compute the entropic part of the weight, namely the factor $`Z(\rho )`$.

Using a Fourier representation of the delta function, the above expression for $`Z(\rho )`$ can be written as

$$Z(\rho )=\frac{1}{M!}𝒟\varphi \mathrm{exp}\left(i𝑑r\varphi (r)\rho (r)+M\mathrm{ln}Q(i\varphi )\right)$$
(6)
where

$$Q(i\varphi )=𝒟re^{\frac{3}{2a^2}_0^N𝑑s\left(\frac{dr}{ds}\right)^2i_0^N𝑑s\varphi (r(s))}$$
(7)
is the partition function of a *single* polymer in the presence of the potential $`i\varphi (r)`$. It is important to note that this is a *purely imaginary potential* because the functional integral in eq. (6) is over a real field $`\varphi (r)`$. We recognize in eq. (7) a Feynman path integral for a quantum mechanical problem with quantum Hamiltonian $`=\frac{a^2}{6}^2+i\varphi (r)`$. Using standard quantum mechanical notation, we may write

$$Q(i\varphi )=𝑑r𝑑r^{}<r|e^N|r^{}>$$
(8)
### II.2 Expansion about the Saddle Point

Of particular interest is an asymptotic expansion of $`Z(\rho )`$ for large $`M`$ to examine the thermodynamic limit of $`M\mathrm{},V\mathrm{},M/V\text{finite}`$. For this purpose, we perform a saddle point expansion of the $`\varphi `$ integral. The first step is to write an equation for the mean-field (saddle point) $`\varphi _0(r)`$ that produces the desired $`\rho (r)`$

$`\rho (r)`$ $`=`$ $`M{\displaystyle \frac{\delta \mathrm{ln}Q(i\varphi _0)}{\delta (i\varphi _0(r))}}`$ (9)
$`=`$ $`{\displaystyle \frac{M}{Q(i\varphi _0)}}{\displaystyle _0^N}𝑑s\mathrm{\Psi }(r,s)\mathrm{\Psi }(r,Ns)`$ (10)
Here, $`s`$ is denotes a contour location along the chain and the “propagator” $`\mathrm{\Psi }(r,s)`$ satisfies the diffusion equation (analogous to the Schrödinger equation)

$$(\frac{}{s}\frac{a^2}{6}^2+i\varphi _0(r))\mathrm{\Psi }(r,s)=0$$
(11)
subject to the initial condition $`\mathrm{\Psi }(r,0)=1`$. The single chain partition function is related to the field $`\mathrm{\Psi }`$ by

$$Q(i\varphi )=𝑑r\mathrm{\Psi }(r,N)$$
(12)
¿From the structure of the above equations, we can make a few general statements. First, a necessary condition for eq. (9) to have a solution $`\varphi _0(r)`$ is that $`\rho (r)`$ should be positive at any point $`r`$ and that its integral over space should be equal to $`MN`$. It is also clear from eqs. (9) and (10) that a self-consistent solution may exist with a real $`i\varphi _0(r)`$, i.e. purely imaginary $`\varphi _0(r)`$. From the structure of these equations, we see that if $`\varphi _0(r)`$ is a mean-field (saddle point) solution, one can add to it any *constant* and obtain a new solution. To lift this degeneracy of the saddle-point, we may specify the value of the field $`\varphi _0`$ at a specific point in space, or its spatially-averaged value. Having made such a choice, it can be proven that such a real solution $`i\varphi _0(r)`$ not only exists, but is a *unique* solution of the above equations.

To perform the asymptotic expansion about the saddle point $`\varphi _0`$ at higher orders, we shift the integration variable $`\varphi `$ by writing

$$\varphi (r)=\varphi _0(r)+\chi (r)$$

The partition function then reads

$$Z(\rho )=Z_0(\rho )𝒟\chi \mathrm{exp}\left(i𝑑r\chi (r)\rho (r)+M\mathrm{ln}𝒬(i\chi )\right)$$
(13)
where

$$𝒬(i\chi )=\frac{𝒟re^{\frac{3}{2a^2}_0^N𝑑s\left(\frac{dr}{ds}\right)^2i_0^N𝑑s\varphi _0(r(s))i_0^N𝑑s\chi (r(s))}}{𝒟re^{\frac{3}{2a^2}_0^N𝑑s\left(\frac{dr}{ds}\right)^2i_0^N𝑑s\varphi _0(r(s))}}$$
(14)
and $`Z_0`$ denotes the mean-field value of the partition function

$$Z_0(\rho )=\frac{1}{M!}\mathrm{exp}\left(i𝑑r\varphi _0(r)\rho (r)+M\mathrm{ln}Q(i\varphi _0)\right)$$
(15)
Next, we expand the exponent in the integrand of eq. (13) in powers of $`\chi `$. For this purpose, it is convenient to reexpress eq. (14) as

$$𝒬(i\chi )=<e^{i_0^N𝑑s\chi (r(s))}>_0$$
(16)
where the notation $`<\mathrm{}>_0`$ describes the expectation value over the conformations of a single polymer subjected to the mean-field potential $`i\varphi _0`$. Expanding in powers of $`\chi `$ and rescaling the field $`\chi `$ by a factor of $`\sqrt{M}`$, we obtain

$$Z(\rho )=Z_0(\rho )𝒟\chi \mathrm{exp}\left(\frac{1}{2}𝑑r𝑑r^{}\chi (r)G_2(r,r^{})\chi (r^{})+W(\chi )\right)$$
(17)
where

$$W(\chi )=\underset{p=3}{\overset{\mathrm{}}{}}\frac{(i)^p}{p!M^{p/21}}𝑑r_1\mathrm{}𝑑r_p\chi (r_1)\mathrm{}\chi (r_p)G_p(r_1,\mathrm{},r_p)$$
(18)
and $`G_p(r_1,\mathrm{},r_p)`$ denotes the $`p`$-point *connected* (cumulant) correlation function for the density of a *single chain*, defined by

$$G_p(r_1,\mathrm{},r_p)=<\overline{\rho }(r_1)\mathrm{}\overline{\rho }(r_p)>_c$$
(19)
where $`\overline{\rho }(r)=_0^N𝑑s\delta (rr(s))`$ is the single chain microscopic density. This expectation value is taken with respect to the conformations of a single chain in the mean-field according to eq. (16) with the usual definition of connected averages amit84 . Note that there is no linear term in $`\chi `$ appearing in the exponent of eq. (17) since the action is stationary with respect to $`\varphi _0`$.

The partition function of eq. (17) can be calculated in a perturbation expansion in powers of $`1/M`$. This is most conveniently done by the use of Feynman diagrams in which propagator lines are associated with factors of $`S_2(r,r^{})=G_2^1(r,r^{})`$ and $`p`$th-order vertices are assigned factors of $`[(i)^p/M^{p/21}]G_p(r_1,\mathrm{},r_p)`$. The result is that

$$Z(\rho )=Z_0(\rho )\mathrm{exp}\left(\frac{1}{2}\mathrm{Tr}\mathrm{ln}G_2(r,r^{})+\mathrm{Sum}\mathrm{of}\mathrm{all}\mathrm{Connected}\mathrm{Diagrams}\right)$$
(20)
where the connected diagrams have one or more vertices of degree $`p3`$. Since the vertices in eq. (19) are connected correlation functions, all terms in the exponent of eq. (20) are extensive, that is they are proportional to the volume $`V`$ of the system. The general expansion of $`Z(\rho )`$ thus takes the form

$$Z(\rho )=\mathrm{exp}\left[\beta V(f_0+\frac{f_1}{M}+\frac{f_2}{M^2}+\mathrm{})\right]$$
(21)
where the $`f_j`$ are contributions to the free energy density that are *independent* of $`M`$ and $`V`$ in the thermodynamic limit. Since all vertices carry powers of $`1/M`$, the sum of all connected graphs is negligible in the thermodynamic limit, and eq.(20) yields

$$Z(\rho )Z_0(\rho )\mathrm{exp}\left(\frac{1}{2}\mathrm{Tr}\mathrm{ln}G_2(r,r^{})\right)$$
(22)
in the limit $`M\mathrm{}`$. We note that $`G_2(r,r^{})`$ is a functional of $`\rho `$, because it represents the 2-point connected correlation function of a chain subjected to the mean-field potential $`i\varphi _0`$ and $`\varphi _0`$ is in turn related to $`\rho `$ by the saddle point condition eqs. (9)-(10).

¿From the form of eq. (22), it would appear that there is a correction to the mean-field value $`Z_0(\rho )`$ of the restricted density partition function $`Z(\rho )`$ that is associated with Gaussian fluctuations in $`\chi (r)`$ about the mean chemical potential field $`\varphi _0(r)`$. Interestingly, these corrections from Gaussian field fluctuations play *no role* in the statistical properties of the field $`\rho `$ for $`M\mathrm{}`$. This is demonstrated by two independent routes in the next section.

### II.3 Significance of the Fluctuation Term

To examine the significance of the fluctuation term $`\frac{1}{2}\mathrm{Tr}\mathrm{ln}G_2(r,r^{})`$, we return to the original formulation of the field theory prior to the asymptotic expansion in $`M`$. From eq. (3) it is clear that we can write an effective, real field theory for the $`\rho `$ field as

$$𝒵=𝒟\rho \mathrm{exp}[H(\rho )]$$
(23)
where the effective Hamiltonian $`H(\rho )`$ is given by

$$H(\rho )=\frac{\beta }{2}𝑑r𝑑r^{}\rho (r)v(rr^{})\rho (r^{})S(\rho )$$
(24)
where $`S(\rho )=\mathrm{ln}Z(\rho )`$ is given in eq. (6). One way to simulate such a field theory is through a Langevin stochastic dynamics parisi88  for $`\rho `$. Note that the density field $`\rho (r)`$ is constrained to be positive at each point $`r`$. A simple way to enforce this condition in the functional integral (23) is to add an external potential which is infinite when $`\rho (r)`$ is negative and $`0`$ when $`\rho (r)`$ is positive. A practical realization of such a potential can be obtained by adding an exponential function to the Hamiltonian, and define

$$H^{}(\rho )=H(\rho )+𝑑re^{\lambda \rho (r)}$$

where $`\lambda `$ is a constant, chosen large enough to ensure the positivity constraint.

If we keep in mind that the total number of monomers is fixed, equal to $`MN`$, we may use the conserved dynamics (model B) halperin\_hohenberg  and write

$`{\displaystyle \frac{}{t}}\rho (r,t)`$ $`=`$ $`\mathrm{\Gamma }^2{\displaystyle \frac{\delta H^{}}{\delta \rho (r,t)}}+\eta (r,t)`$ (25)
$`=`$ $`\mathrm{\Gamma }^2\left(\beta {\displaystyle 𝑑r^{}v(rr^{})\rho (r^{},t)}<i\varphi (r,t)>_\rho \lambda e^{\lambda \rho (r)}\right)+\eta (r,t)`$
where we have used

$$\frac{\delta S(\rho )}{\delta \rho (r,t)}=\frac{\delta \mathrm{log}Z(\rho )}{\delta \rho (r,t)}=<i\varphi (r,t)>_\rho $$

and $`\eta (r,t)`$ is a Gaussian, white noise with vanishing first moment, $`\eta (r,t)=0`$, and second moment satisfying the usual fluctuation-dissipation theorem for a conserved noise,

$$<\eta (r,t)\eta (r^{},t^{})>=2\mathrm{\Gamma }^2\delta (rr^{})\delta (tt^{})$$
(26)
The Laplacian in front of the r.h.s. of (25) ensures that the space integral of $`\rho (r,t)`$, which is equal to the total number of monomers $`MN`$, remains constant in time. In the above equation, the expectation value $`<i\varphi (r)>_\rho `$ denotes an equilibrium average of the field $`i\varphi (r)`$ at constrained $`\rho (r)`$ using the statistical weight contained in eq. (6) for $`Z(\rho )`$.

With the use of eq. (17), it is easily seen that

$$<i\varphi (r)>_\rho =i\varphi _0(r)+\frac{i}{\sqrt{M}}<\chi (r)>_\rho $$

where

$$<\chi (r)>_\rho =\frac{𝒟\chi \chi (r)\mathrm{exp}\left(\frac{1}{2}𝑑r𝑑r^{}\chi (r)G_2(r,r^{})\chi (r^{})+W(\chi )\right)}{𝒟\chi \mathrm{exp}\left(\frac{1}{2}𝑑r𝑑r^{}\chi (r)G_2(r,r^{})\chi (r^{})+W(\chi )\right)}$$

and $`W(\chi )`$ is given by eq. (18). From this expression it follows immediately that

$$<i\varphi (r)>_\rho =i\varphi _0(r)+O(\frac{1}{M})$$
(27)
and thus, in the Langevin equation formalism, it is sufficient to evaluate $`<i\varphi (r)>_\rho `$ at the mean-field level to capture the leading asymptotic behavior in the thermodynamic limit. The evaluation of the Langevin forces on the right hand side of eq. (25) therefore requires only the computation of $`\varphi _0(r,t)`$, which follows from the solution of the saddle point equations at prescribed $`\rho (r,t)`$. The forces entering the Langevin equation are real (since $`i\varphi _0`$ is real), and thus we have a conventional real Langevin dynamics. The key result is that in the thermodynamic limit, there is no signature of the fluctuation term $`(1/2)\mathrm{Tr}\mathrm{ln}G_2(r,r^{})`$ that appears in eq. (22) for $`Z(\rho )`$.

The same result can be directly obtained from eq.(22) by means of some straightforward manipulations. Indeed, from that equation and eq. (15), we have

$`{\displaystyle \frac{\delta \mathrm{ln}Z(\rho )}{\delta \rho (r_0)}}`$ $`=`$ $`<i\varphi (r_0)>_\rho `$ (28)
$`=`$ $`i\varphi _0(r_0){\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta \mathrm{Tr}\mathrm{ln}G_2(r,r^{})}{\delta \rho (r_0)}}`$ (29)
The $`\mathrm{Tr}\mathrm{ln}`$ term can be further simplified. Using some simple identities, we get

$`{\displaystyle \frac{\delta \mathrm{Tr}\mathrm{ln}G_2(r,r^{})}{\delta \rho (r_0)}}`$ $`=`$ $`{\displaystyle 𝑑r𝑑r^{}\frac{\delta G_2(r,r^{})}{\delta \rho (r_0)}G_2^1(r^{},r)}`$ (30)
$`=`$ $`{\displaystyle 𝑑r𝑑r^{}𝑑r_1\frac{\delta (i\varphi _0(r_1))}{\delta \rho (r_0)}\frac{\delta G_2(r,r^{})}{\delta (i\varphi _0(r_1))}G_2^1(r^{},r)}`$
¿From our previous definitions, we see that

$$\frac{\delta G_2(r,r^{})}{\delta (i\varphi _0(r_1))}=G_3(r,r^{},r_1)$$
(31)
and the saddle point condition implies

$$\frac{\delta (i\varphi _0(r_1))}{\delta \rho (r_0)}=\frac{1}{M}G_2^1(r_1,r_0)$$
(32)
Substitution of these results into eq. (30) leads to a final expression

$$\frac{\delta \mathrm{ln}Z(\rho )}{\delta \rho (r_0)}=i\varphi (r_0)\frac{1}{2M}𝑑r𝑑r^{}𝑑r_1G_2^1(r_1,r_0)G_3(r,r^{},r_1)G_2^1(r^{},r)$$
(33)
Again, since $`G_2`$ and $`G_3`$ are connected Green’s functions, the triple integral in eq. (33) is finite in the thermodynamic limit (non extensive) and this shows that to leading order in $`1/M`$, the expectation value of $`\varphi `$ is given by its mean-field value $`\varphi _0`$, in agreement with eq. (27). In addition, it is easy to see that the correction term calculated here is exactly the same as the one that would have been obtained as the leading term of the perturbation expansion using eq. (17). We conclude that the fluctuation term $`(1/2)\mathrm{Tr}\mathrm{ln}G_2(r,r^{})`$ can be safely omitted from the effective Hamiltonian $`H(\rho )`$ when computing the Langevin force $`\delta H/\delta \rho `$ for the purpose of simulating large systems. Moreover, it follows that the fluctuation correction can also be dropped from the Hamiltonian itself, if Monte Carlo simulations are to be employed.

We summarize with the statement that in the thermodynamic limit of large systems, the field theory expressed in density variables can be written as

$$𝒵=𝒟\rho \mathrm{exp}[H^{}(\rho )]$$
(34)
with

$$H^{}(\rho )=\frac{\beta }{2}𝑑r𝑑r^{}\rho (r)v(rr^{})\rho (r^{})+𝑑re^{\lambda \rho (r)}𝑑ri\varphi _0(r)\rho (r)M\mathrm{ln}Q(i\varphi _0)$$
(35)
and where the real field $`i\varphi _0(r)`$ is obtained at a particular $`\rho (r)`$ by the solution of eq. (9).

### II.4 Langevin Dynamics and Monte Carlo: Implementation

The Langevin equation (25) can be discretized both in time and space. A variety of time integration schemes for such nonlinear stochastic differential equations are available kloeden .

It is of interest to compare this *real* Langevin scheme with the *complex* Langevin scheme that has been used to numerically simulate field theories of polymer solutions and melts expressed in chemical potential \[i.e. $`\varphi (r)`$\] variables ganesan01 ; katz02 . Because the effective Hamiltonian $`𝒢(\varphi )`$ in such field theories is complex, one must address the sign problem. A useful strategy, originally devised by Parisi and Klauder parisi83 ; klauder84 , has been to extend the field variables to the complex plane ($`\varphi =\varphi _R+i\varphi _I`$) and write a Langevin dynamics

$$\frac{}{t}\varphi _R(r,t)=\mathrm{\Gamma }\mathrm{Re}\frac{\delta 𝒢}{\delta \varphi (r,t)}+\eta (r,t)$$
(36)
$$\frac{}{t}\varphi _I(r,t)=\mathrm{\Gamma }\mathrm{Im}\frac{\delta 𝒢}{\delta \varphi (r,t)}$$
(37)
where $`\eta (r,t)`$ is a real Gaussian, white noise with the following covariance

$$<\eta (r,t)\eta (r^{},t^{})>=2\mathrm{\Gamma }\delta (rr^{})\delta (tt^{})$$
(38)
This procedure, if it converges gausterer98 ; lee94 ; schoenmaker87 ; ghf03 , will properly describe the stationary fluctuation spectrum of a model with complex $`𝒢(\varphi )`$. We have found it to be reasonably effective at suppressing phase oscillations in simulations of polymer solutions and melts, although sign problem is still very pronounced in dilute and semidilute systems. An advantage of eqs. (36)-(37) is that the computation of the complex force $`\delta 𝒢/\delta \varphi `$ requires only a *single pass* at solving the diffusion equation (11) (with $`\varphi _0\varphi `$), while computation of $`\delta ^{}/\delta \rho `$ in the real Langevin scheme of eq. (25) requires multiple passes at solving the diffusion equation in order to establish the field $`\varphi _0`$ consistent with the current $`\rho `$. Whether the extra computational burden of the force evaluation in the present density-based real Langevin scheme is offset by avoidance of the sign problem remains to be seen and, indeed, may prove to be system and problem dependent.

Note that in a Monte Carlo implementation, the minor complications due to the positivity of $`\rho (r)`$ and the conservation of the integral $`𝑑r\rho (r)`$ are easy to overcome. First, one discretizes space. To implement the conservation of the total number of monomers, one starts with a configuration of $`\rho (r)`$ that has the right number of monomers, and then one chooses randomly two lattice sites, say $`r_1`$ and $`r_2`$ . The elementary Monte Carlo variation of $`\rho (r)`$ consists in modifying simultaneously $`\rho `$ at two points $`r_1`$ and $`r_2`$ according to $`\rho (r_1)\rho (r_1)+\delta \rho `$ and $`\rho (r_2)\rho (r_2)\delta \rho `$ . This procedure obviously conserves the total integral of $`\rho (r)`$ in space. Also, to enforce the positivity of $`\rho (r)`$, one just has to reject any variation of $`\rho `$ that produces a negative $`\rho `$.

## III More complicated systems

The concepts and methods presented above can be extended in a straightforward way to virtually any type of inhomogeneous polymer system. We illustrate in the context of a solution of AB diblock copolymers.

Let $`f`$ be the volume fraction of A monomers and $`1f`$ the volume fraction of B monomers on each diblock copolymer chain with a total of $`N`$ monomers. Denoting by $`v_{AA}(r)`$, $`v_{BB}(r)`$, and $`v_{AB}(r)`$ the respective potentials of mean force (mediated by the solvent) between the different species of monomers, and using the same notations as in the previous section, we may write

$`𝒵`$ $`=`$ $`{\displaystyle \frac{1}{M!}}{\displaystyle 𝒟\rho _A𝒟\rho _B\mathrm{exp}\left(\frac{\beta }{2}𝑑r𝑑r^{}\rho _A(r)v_{AA}(rr^{})\rho _A(r^{})\right)}`$ (39)
$`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{\beta }{2}}{\displaystyle 𝑑r𝑑r^{}\rho _B(r)v_{BB}(rr^{})\rho _B(r^{})}\right)`$
$`\times `$ $`\mathrm{exp}\left(\beta {\displaystyle 𝑑r𝑑r^{}\rho _A(r)v_{AB}(rr^{})\rho _B(r^{})}\right)Z_A(\rho _A)Z_B(\rho _B)`$
1where the functions $`Z_A(\rho _A)`$ and $`Z_B(\rho _B)`$ are defined by

$`Z_A(\rho _A)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{M}{}}𝒟r_k\mathrm{exp}\left(\frac{3}{2a^2}\underset{k=1}{\overset{M}{}}_0^{fN}𝑑s\left(\frac{dr_k}{ds}\right)^2\right)}`$ (40)
$`\times `$ $`{\displaystyle \underset{r}{}}\delta \left(\rho _A(r){\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle _0^{fN}}𝑑s\delta (rr_k(s))\right)`$
and

$`Z_B(\rho _B)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{M}{}}𝒟r_k\mathrm{exp}\left(\frac{3}{2a^2}\underset{k=1}{\overset{M}{}}_{fN}^N𝑑s\left(\frac{dr_k}{ds}\right)^2\right)}`$ (41)
$`\times `$ $`{\displaystyle \underset{r}{}}\delta \left(\rho _B(r){\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle _{fN}^N}𝑑s\delta (rr_k(s))\right)`$
It is clear from their definitions that both $`Z_A`$ and $`Z_B`$ are positive definite, and using a Fourier representation for the delta functionals, one can again show that in the limit of an infinite number of chains (with infinite volume and finite concentration), the corresponding “chemical potential” field integrals can be asymptotically evaluated by the saddle point method. The resulting expression for the entropy functional

$$S(\rho _A,\rho _B)=\mathrm{ln}[Z_A(\rho _A)Z_B(\rho _B)]$$
(42)
is

$$S(\rho _A,\rho _B)=𝑑r(i\varphi _A\rho _A+i\varphi _B\rho _B)+M\mathrm{ln}Q(i\varphi _A,i\varphi _B)$$
(43)
where $`Q(i\varphi _A,i\varphi _B)=𝑑r\mathrm{\Psi }(r,N)`$. The propagator $`\mathrm{\Psi }(r,s)`$ satisfies (again, subject to $`\mathrm{\Psi }(r,0)=1`$)

$$(\frac{}{s}\frac{a^2}{6}^2+w(r,s))\mathrm{\Psi }(r,s)=0$$
(44)
where $`w(r,s)=i\varphi _A(r)`$ for $`s(0,Nf)`$ and $`w(r,s)=i\varphi _B(r)`$ for $`s(Nf,N)`$. The density fields are related to the saddle point chemical potentials, $`\varphi _A`$ and $`\varphi _B`$, by equations analogous to eq. (10)

$$\rho _A(r)=\frac{M}{Q(i\varphi _A,i\varphi _B)}_0^{Nf}𝑑s\mathrm{\Psi }(r,s)\mathrm{\Psi }^{}(r,Ns)$$
(45)
$$\rho _B(r)=\frac{M}{Q(i\varphi _A,i\varphi _B)}_{Nf}^N𝑑s\mathrm{\Psi }(r,s)\mathrm{\Psi }^{}(r,Ns)$$
(46)
A second propagator $`\mathrm{\Psi }^{}(r,s)`$ enters these expressions due to the head-to-tail asymmetry of diblock copolymers. It satisfies

$$(\frac{}{s}\frac{a^2}{6}^2+w^{}(r,s))\mathrm{\Psi }^{}(r,s)=0$$
(47)
subject to $`\mathrm{\Psi }^{}(r,0)=1`$, where $`w^{}(r,s)=i\varphi _B(r)`$ for $`s(0,N(1f))`$ and $`w^{}(r,s)=i\varphi _A(r)`$ for $`s(N(1f),N)`$.

A formalism very similar to this was employed recently, without formal justification, for a real Langevin dynamics study of spinodal decomposition in binary polymer blends reister  and in studying microemulsion formation in ternary copolymer alloys dominik . The present analysis shows that this is indeed a rigorous approach to simulating fluctuating systems of macroscopic size.

## IV Summary and Discussion

In the present paper we have examined statistical field theories of polymeric fluids with the objective of identifying field-based models that are well-suited to numerical simulation. It has been known for some time that essentially any microscopic (particle-based) model of a polymeric fluid can be converted to a field theory in density and chemical potential variables matsen94 ; hongnoolandi81 , e.g. for a one-component system

$$𝒵=𝒟\rho 𝒟\varphi \mathrm{exp}[(\rho ,\varphi )]$$
(48)
where the effective Hamiltonian $`(\rho ,\varphi )`$ is a *complex* functional of the two fields. In the case of pair interactions among monomers, $`(\rho ,\varphi )`$ is a quadratic form in $`\rho `$, so the density field can be integrated out exactly (for potentials with an inverse) to obtain a statistical field theory in $`\varphi `$

$$𝒵=𝒟\varphi \mathrm{exp}[𝒢(\varphi )]$$
(49)
In this case, the effective Hamiltonian $`𝒢(\varphi )`$ is complex. This type of field theory has been studied using approximate analytical methods for many years edwards65  and has recently been examined with complex Langevin simulation techniques ganesan01 ; ghfrev ; katz02 . Unfortunately, the non-positive definite character of the Boltzmann factor $`\mathrm{exp}(𝒢)`$ gives the integrand an oscillatory character (“sign problem”), which can dramatically hinder the convergence of numerical simulations.

Rather than transforming to a field theory in the single field $`\varphi `$, a more typical strategy for studying inhomogeneous polymers is to invoke a mean-field approximation. This amounts to evaluating both integrals in eq. (48) at leading order by the saddle point method. This reproduces the familiar self-consistent field theory (SCFT) helfand75  and amounts to

$$𝒵\mathrm{exp}[(\rho ^{},\varphi ^{})]$$
(50)
where $`\rho ^{}`$ and $`\varphi ^{}`$ are the “self-consistent” density and chemical potential fields that satisfy $`\delta /\delta \rho ^{}=0`$, $`\delta /\delta \varphi ^{}=0`$. This is a powerful approach, but it neglects all field fluctuations, so the utility of the SCFT is limited to highly concentrated systems of high molecular weight polymers, far from critical points and phase transitions where soft fluctuation modes occur.

In the present paper, we have discussed a third approach in which the chemical potential field $`\varphi `$ is integrated out of eq. (48), leading to a field theory solely in the monomer density

$$𝒵=𝒟\rho \mathrm{exp}[H(\rho )]$$
(51)
This approach has the benefit of producing a *real* field theory in which $`\mathrm{exp}(H)`$ is positive definite so that the sign problem is avoided and standard Monte Carlo and real Langevin simulation methods can be applied. In the past it was assumed that the $`\varphi `$ integral necessary to reduce eq. (48) to eq. (51) was intractable, leading various authors to impose a variety of approximations in order to construct the density functional $`H(\rho )`$ leibler80 ; tang91 . We have demonstrated here, however, that the $`\varphi `$ integral in eq. (48) can be *exactly* evaluated by the saddle point method provided that the thermodynamic limit is taken. The resulting real energy functional $`H(\rho )`$ is given (for a homopolymer solution) by eq. (35) and involves the constrained (real) saddle point $`i\varphi _0(r)`$, which is the mean chemical potential field required to generate a particular density pattern $`\rho (r)`$.

We believe that our discovery will have a number of significant implications. First of all, it provides theoretical justification for ad hoc approaches that have been used to include thermal fluctuations in density field-based simulations of polymers reister ; fraaije96 . Secondly, by circumventing the sign problem, it avoids the difficulties that have been encountered in field-theoretic polymer simulations based on eq. (49). Moreover, with a real field theory, a much wider variety of options exist for carrying out stochastic computer simulations. Finally, we note that our proof that the $`\varphi `$ integral in eq. (48) is dominated for large systems by the saddle point carries over to a much broader class of classical fluids, both simple and complex. For example, real density functionals can be constructed using the present methods for alternative polymer chain models, such as the worm-like chain used to describe semiflexible polymers. Real density-based field theories for a wide variety of polymers, copolymers, and alloys can be derived, including solutions of charged polyelectrolyte systems. In the case of simple fluids, we believe that the present formalism will provide new insights and numerical strategies for improving density functional methods, including the treatment of crystallization and melting.

A clear disadvantage of the density field approach, at least for polymers, is that force or energy evaluations require the computation of the mean-field $`i\varphi _0(r)`$ from equations such as (10) at constrained monomer density $`\rho (r)`$. This must be done iteratively and each iteration requires an independent solution of the diffusion equation (11). In constrast, evaluation of the energy $`𝒢(\varphi )`$ in the potential field formulation of eq. (49) requires only a single solution of the diffusion equation. Since the vast majority of the computational effort in a field-theoretic computer simulation of a polymeric fluid is spent in solving the diffusion equation, it is essential that a very efficient scheme be devised for computing $`i\varphi _0`$ given $`\rho `$. We are cautiously optimistic, however, that the extra computational burden associated with simulating in the density variables will be more than offset by the avoidance of the sign problem.

###### Acknowledgements.

This work was supported by the National Science Foundation under award DMR 03-12097. One of the authors (H.O.) wishes to thank the MRL at UCSB for its kind hospitality during the course of this work.


format: text

text:
# Optical, Ultraviolet, and Infrared Observations of SN 1993J

## 1 Introduction

Supernova (SN) 1993J was visually discovered in the nearby galaxy M81 (NGC 3031; $`d`$ = 3.6 Mpc; Freedman et al. 1994) by Francisco Garcia on 1993 March 28.906 UT (Ripero, Garcia, & Rodriguez 1993). Ten years later, we still vividly remember our excitement after receiving by fax a rough finding chart from the discoverers, via the AAVSO. Ever since SN 1987A, those of us in the northern hemisphere had been eagerly waiting for a supernova in M31, or (more likely) in a somewhat more distant galaxy. What a golden opportunity! Unfortunately, our own robotic SN search, at that time being conducted with a 0.76-m telescope at Leuschner Observatory, missed discovery of SN 1993J because its field of view was somewhat too small to include the position of the supernova when centered on the galaxy nucleus, and because of bad weather (Filippenko 1993). It was a pleasure for both of us to finally meet Mr. Garcia, at this conference.

SN 1993J reached $`V=10.8`$ mag (e.g., Richmond et al. 1994), becoming the brightest SN in the northern hemisphere since SN 1954A ($`m_{\mathrm{pg}}=9.95`$; Wild 1960; Barbon, Ciatti, & Rosino 1973). In terms of observational coverage, both in temporal consistency (almost nightly observations at early times) and in the details of individual observations (including observations with signal-to-noise ratios, spectral resolutions, and wavelength regions not typically found in studies of supernovae), SN 1993J is surpassed only by SN 1987A. Early spectra showed an almost featureless blue continuum, possibly with broad, but weak, H$`\alpha `$ and He I $`\lambda `$5876 lines. This led to a Type II classification (Filippenko et al. 1993; Garnavich & Ann 1993; see Filippenko 1997 for a general discussion of SN types). Wheeler & Filippenko (1996) present a thorough review of the early work on SN 1993J.

Both the spectra and the light curves of SN 1993J quickly began to indicate that this was not a typical Type II SN. Indeed, the initially unusual light curves and the appearance of He I lines in the spectra were interpreted as evidence that SN 1993J was similar to a SN Ib, with a low-mass outer layer of hydrogen that gave the early impression of a SN II (see discussion and references below). Following Woosley et al. (1987), it was described as a “Type IIb” SN. This transformation from SN II to nearly SN Ib indicates a common mechanism (core collapse) for these two observationally defined subclasses. SN 1993J is thus one of the most significant SNe ever studied, not only for its role in linking Types II and Ib (and possibly Ic), but also because it was observed with such great detail.

Here we review observations of SN 1993J at optical, ultraviolet (UV), and infrared (IR) wavelengths, concentrating on optical because of the vast amount of data obtained in this wavelength range. Because of time and space constraints, we exclude light echoes and interstellar absorption lines, which have been used to probe the interstellar medium near SN 1993J. A recent optical spectrum of SN 1993J, obtained with the Keck-I 10 m telescope 9 years and 11 months after the explosion, is shown in Figure 1. This is very close to 10 years, so we shall dub this the official “Valencia spectrum of SN 1993J.”

## 2 Photometry of SN 1993J

The evolution of the light curve of SN 1993J did not follow either of the two typical paths for SNe II. SN 1993J did not remain at a relatively constant brightness after a slight decline from maximum, as a normal Type II plateau SN would, nor did the brightness decline in the pattern of a Type II linear SN. (For representatives of these Type II light curves, see, for example, Doggett & Branch 1985.) Instead, SN 1993J rose quickly, then rapidly declined for $``$1 week, only to brighten a second time over the next two weeks. This led to another rapid decrease in brightness for $``$3 weeks, followed by an approximately exponential decline. For a complete discussion of the photometry of SN 1993J, see Okyudo et al. (1993), Schmidt et al. (1993), van Driel et al. (1993), Wheeler et al. (1993), Benson et al. (1994), Lewis et al. (1994), Richmond et al. (1994, 1996), Barbon et al. (1995), Doroshenko et al. (1995), and Prabhu et al. (1995).

The unusual initial behavior of the light curve rapidly led many SN modelers to conclude that SN 1993J was the result of a core-collapse explosion in a progenitor that had lost a significant fraction of its hydrogen envelope, leaving only $``$0.1–0.5 $`M_{}`$ of hydrogen. The original envelope could have been lost through winds (Höflich, Langer, & Duschinger 1993) from a fairly massive star (25–30 $`M_{}`$). Another possibility explored by Hashimoto, Iwamoto, & Nomoto (1993; see also Nomoto et al. 1993) is that SN 1993J was the result of the explosion of an asymptotic giant branch star having main-sequence mass $`M_{ms}7`$–10 $`M_{}`$, with a helium-rich envelope.

A more likely solution is that the progenitor of SN 1993J was a member of a binary system and the companion had stripped away a considerable amount of hydrogen. The progenitor was observed during prior studies of M81. Aldering, Humphreys, & Richmond (1994) analyzed several sets of pre-existing images and deduced that the photometry was inconsistent with a single star at the position of SN 1993J. They found that the best fit for the progenitor itself was a K0 I star with $`M_{\mathrm{bol}}7.8`$ mag and $`VR0.7`$ mag. Cohen, Darling, & Porter (1995) derived a similar color from a five-month series of images of M81 from 1984; there was no apparent variability.

Using the scenario of a star that had been stripped of most of its hydrogen envelope, Nomoto et al. (1993) and Shigeyama et al. (1994) found a best fit to the light curve from their model of a $`4M_{}`$ helium core, although a range of 3–6 $`M_{}`$ for the core is reasonable. The main-sequence mass of the star would have been 15 $`M_{}`$, while the residual hydrogen envelope is less than $`0.9M_{}`$. Starting with a star of initial mass of 13–16 $`M_{}`$, Woosley et al. (1994) could reproduce the light curve from the explosion of a remaining helium core with mass $`4.0\pm 0.5M_{}`$ and hydrogen envelope with mass $`0.20\pm 0.05M_{}`$. A similar model by Podsiadlowski et al. (1993) had $`0.2M_{}`$ of hydrogen remaining on a star of initial mass $`M_i15M_{}`$. Ray, Singh, & Sutaria (1993) also invoked a binary system for SN 1993J with a residual hydrogen envelope mass of $`0.2M_{}`$. Utrobin (1994) used an envelope mass of $`0.1M_{}`$ remaining on a $`3M_{}`$ helium core from an initial mass of $`12M_{}`$. Bartunov et al. (1994) achieved a good fit to the light curve with a helium core mass of 3.5 $`M_{}`$, but a larger hydrogen envelope ($`M_{env}0.9M_{}`$). Later studies continued to conclude that a low-mass envelope of hydrogen on a helium core was the most likely scenario for the progenitor (Young, Baron, & Branch 1995; Utrobin 1996). Intercomparison of two methods also indicated that the results were robust (Blinnikov et al. 1998). Houck & Fransson (1996) used a non-local thermodynamic equilibrium (NLTE) synthetic spectrum code to fit nebular spectra and found that the Nomoto et al. (1993) models could explain the late-time spectra. They found a best fit with a 3.2 $`M_{}`$ helium core with a 0.2–0.4 $`M_{}`$ hydrogen envelope. Patat, Chugai, & Mazzali (1995) also used the late-time spectra, specifically the H$`\alpha `$ line, to derive an ionized hydrogen mass of 0.05–0.2 $`M_{}`$; this is a lower limit to the envelope mass.

## 3 Spectroscopy of SN 1993J

Woosley et al. (1987) had already considered the above possibility for core-collapse SNe, giving them a new name: SNe IIb. The low-mass outer layer of hydrogen would give the initial appearance of a SN II, but the spectrum would slowly change to one more similar to that of a SN Ib, dominated by helium lines with the hydrogen either appearing weakly or completely gone. Indeed, Nomoto et al. (1993) predicted that the spectrum of SN 1993J would show this behavior. This was first confirmed by Filippenko & Matheson (1993), followed rapidly by Schmidt et al. (1993) and Swartz et al. (1993). Studies of the early optical spectra include those of Filippenko, Matheson, & Ho (1993), Wheeler et al. (1993), Taniguchi et al. (1993), Garnavich & Ann (1994), Ohta et al. (1994), Prabhu et al. (1995), and Metlova et al. (1995).

Jeffery et al. (1994) present an early UV spectrum of SN 1993J taken with the Hubble Space Telescope (HST) on 1993 April 15 UT. The other core-collapse SNe that had been observed in the UV to that point were compared with SN 1993J, and there were striking differences. SN 1993J had a relatively smooth UV spectrum and was more similar to SN 1979C and SN 1980K, both of which are radio sources and thus likely to have thick circumstellar envelopes (e.g., Weiler et al. 1986). The UV spectra of SN 1987A, in contrast, showed broad absorption features. The illumination from circumstellar interaction may reduce the relative strengths of line features compared to the continuum and thus produce the featureless UV spectra of SNe 1979C, 1980K, and 1993J (Branch et al. 2000).

SN 1993J then evolved fairly rapidly into the nebular phase. The nebular-phase spectra were similar to those of a typical SN Ib, but the hydrogen lines never faded completely. In fact, H$`\alpha `$ began to dominate the spectrum at late times, certainly the result of circumstellar interaction. There were several other papers that considered the nebular-phase spectra (and some relatively late-time spectra). Filippenko, Matheson, & Barth (1994) show the transition to the nebular phase with spectra from Lick Observatory. Lewis et al. (1994) present the complete La Palma archive covering days 2 through 125. Li et al. (1994) discuss the nebular-phase spectra observed from the Beijing Astronomical Observatory. Barbon et al. (1995) show the first year of observations from Asiago; the transformation of the SN from Type II to IIb is evident, as is the return of H$`\alpha `$ at late times (by $``$ 200 days). A longer baseline ($`500`$ days) for the spectra is found in the work of Finn et al. (1995).

There were optical spectropolarimetric observations of SN 1993J. Trammell, Hines, & Wheeler (1993) found a continuum polarization of $`P=1.6\%\pm 0.1\%`$ on day 24 (assuming 1993 March 27.5 UT as the explosion date). Trammell et al. (1993), as well as later considerations of the same data (Höflich 1995; Höflich et al. 1996), argued that this polarization implied an overall asymmetry, but the source of this asymmetry was undetermined. The presence of SN 1993J in a binary system was implicated as a potential source for the asymmetry. With more epochs of observation, Tran et al. (1997) also found a polarization in the continuum of $``$ 1%, but a different level for the interstellar polarization. Nevertheless, they also concluded that SN 1993J was asymmetric. It is interesting to note that a subsequent SN IIb, SN 1996cb, showed substantially similar polarization of its spectra (Wang et al. 2000).

The analysis of individual aspects of the spectra has yielded some interesting results. Both Wang & Hu (1994) and Spyromilio (1994) found evidence for clumpy ejecta with blueshifted emission lines. Houck & Fransson (1996) argue that the lines are not actually blueshifted, but that contamination from other lines appears to shift them. Nonetheless, the lines do show substructure that indicates clumpy ejecta.

Models of the early spectra could reproduce their overall spectral shape, but the line strengths were problematic. Baron et al. (1993) found a photospheric temperature of $``$8000 K for day 10, but the predicted hydrogen and helium lines were too weak, possibly indicating unusual abundances or non-thermal effects. A later analysis including the HST UV spectrum was fit well by including enhanced helium abundance and NLTE effects (Baron, Hauschildt, & Branch 1994). Jeffery et al. (1994) also had difficulties fitting line strengths for transitions that are susceptible to NLTE effects. Clocchiatti et al. (1995) studied the early spectra to follow the evolution of color temperature and to calculate a distance to M81 ($``$ 3.5 Mpc) using the expanding photosphere method (e.g., Eastman, Schmidt, & Kirshner 1996). The NLTE treatment of calcium is explored by Zhang & Wang (1996), who found a best fit with a reduced calcium abundance.

## 4 Recent Studies of SN 1993J

Most of the above papers were published before year 2000, but in this section we discuss several more recent studies.

Matheson et al. (2000a) present a series of 42 Lick and Keck low-resolution optical spectra of SN 1993J from day 3 after explosion to day 2454, as well as one Keck high-dispersion spectrum from day 383. The spectra are studied in detail by Matheson et al. (2000b). Spectra during the nebular phase, but within the first two years after explosion, exhibit small-scale structure in the emission lines of some species, notably oxygen and magnesium, showing that the ejecta of SN 1993J are clumpy. On the other hand, a lack of structure in emission lines of calcium implies that the source of calcium emission is uniformly distributed throughout the ejecta. These results are interpreted as evidence that oxygen emission originates in clumpy, newly synthesized material, while calcium emission arises from material pre-existing in the atmosphere of the progenitor (Li & McCray 1992, 1993). Spectra spanning the range 433–2454 days after the explosion show box-like profiles for the emission lines, clearly indicating circumstellar interaction in a roughly spherical shell. This is interpreted within the Chevalier & Fransson (1994) model for SNe interacting with mass lost during prior stellar winds. At very late times, the emission lines have a two-horned profile, implying the formation of a somewhat flattened or disk-like structure that is a significant source of emission.

Matthews et al. (2002) conducted IR photometry (windows in the interval 1.25–3.7 $`\mu `$m) and IR spectroscopy (windows in the interval 1.2–2.4 $`\mu `$m) of SN 1993J at early times, through about day 250. As in the case of the optical bands, the IR brightness rose to a secondary maximum and then dropped exponentially. However, the $`L^{}`$ (3.7 $`\mu `$m) light curve exhibited an excess, beginning at day 130, which Matthews et al. interpret as thermal emission from dust, as in the SNe II 1987A and 1998S. At early times, during the rise to the secondary maximum, the spectral energy distribution (SED) of SN 1993J could be fitted with black bodies, but such fits were too broad for the observed SEDs during the exponential decline. Though initially featureless, the IR spectra subsequently (during the exponential decline) became dominated by H, He, and probably Fe line emission.

The HST SINS collaboration (Supernova INtensive Study) obtained a series of UV spectra of SN 1993J over the course of about 7 years (C. Fransson et al. 2004, in preparation). It is quite clear that at late times, the emission is almost entirely coming from the shock interaction between the ejecta and circumstellar gas. The model fit to the observed spectrum is quite good, and its total luminosity is fixed by the observed X-ray emission from the shock – it is not a free parameter. Moreover, the detailed line profiles show that the emission is coming from a finite shell rather than a centrally peaked distribution of gas. There are strong, broad emission lines of N II\], N III\], and N IV\], but only weak lines of C III\] and C IV — characteristic of CNO processing. Indeed, a model spectrum having C:N = 1:13 by number in the circumstellar gas (with which the ejecta collide) gives a relatively good fit to the observations. CNO abundance ratios have previously been measured for a number of SNe II, with similar (but perhaps less extreme) results.

¿From HST images with $`0.^{\prime \prime }05`$ resolution, Van Dyk et al. (2002) identify four stars brighter than $`V=25`$ mag within $`2.^{\prime \prime }5`$ of SN 1993J that contaminated previous ground-based brightness estimates for the supernova progenitor. Correcting for the contamination, they find that the energy distribution of the progenitor is consistent with that of an early K-type supergiant star with $`M_V7.0\pm 0.4`$ mag and an initial mass of 13–22 $`M_{}`$. The brightnesses of the nearby stars are sufficient to account for the excess blue light seen from the ground in pre-explosion observations (Aldering et al. 1994). Therefore, the SN 1993J progenitor did not necessarily have a blue companion, although by 2001, fainter blue stars are seen in close proximity to the supernova. These observations do not strongly limit the mass of a hypothetical companion. A blue dwarf star with a mass up to 30 $`M_{}`$ could have been orbiting the progenitor without being detected in the ground-based images.

Explosion models and observations show that the SN 1993J progenitor had a helium-rich envelope. To test whether the helium abundance could influence the energy distribution of the progenitor, Van Dyk et al. (2002) calculated model supergiant atmospheres with a range of plausible helium abundances. The models show that the pre-supernova colors are not strongly affected by the helium abundance longward of 4000 Å, and abundances ranging between solar and 90% helium (by number) are all consistent with the observations.

Recent optical spectra of SN 1993J (e.g., Figure 1) show that the optical continuum is still very bright ($`BV21`$ mag), certainly from the circumstellar interaction. Van Dyk et al. (2002) suggest that the putative, blue companion star has $`B>23`$ mag and $`U>22`$ mag, so perhaps the most likely spectral region in which it might be observed is in the $`U`$ band, where the actual observed continuum has $`U22`$ mag (Figure 1). At this conference, Stephen Smartt and collaborators reported the possible detection of hydrogen Balmer lines from the putative companion. Given the complexity of the SN spectrum (blends of emission lines, uncertain continuum level), however, this must be verified with future observations; most or all of the claimed Balmer lines may have other explanations. It is also not clear whether the Balmer lines, if real, come from a physical companion star or from an unrelated star along the same line of sight; the HST images shown by Van Dyk et al. (2002) reveal several possible contaminants.

## 5 Final Remarks

Although SN 1993J has provided the best observational evidence for the transformation of a SN from one type to another, there have been other examples. The early spectra of SN 1987K showed hydrogen lines, but the late-time spectra more closely resembled those of SNe Ib (Filippenko 1988). The transition itself was not observed, occurring while SN 1987K was in conjunction with the Sun. SN 1996cb underwent a very similar metamorphosis from SN II to SN Ib; Qiu et al. (1999) present a complete spectroscopic record of the transformation. A number of other examples of genuine SNe IIb have been found in recent years. In addition, there were some suggestions of hydrogen in spectra of the Type Ic SN 1987M (Jeffery et al. 1991; Filippenko 1992) and the SN Ic 1991A (and perhaps SN Ic 1990aa; Filippenko 1992). SN 1993J is clearly a significant object in the study of SNe. By providing a link between SNe II and SNe Ib, it has strengthened the argument that SNe Ib (and, by extension, SNe Ic) are also core-collapse events.

## 6 Acknowledgements

A.V.F.’s research on SN 1993J is currently supported by NSF grant AST-0307894, as well as by NASA grants AR-9529 and GO-9114 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. We thank the conference organizers for partial travel funds. We are grateful to the staffs of the Lick and Keck Observatories for help with the observations. The W. M. Keck Observatory is operated as a scientific partnership among Caltech, the University of California, and NASA; the Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

## Index


text:
# NGC 4261 and NGC 4697: rejuvenated elliptical galaxies

## 1 Introduction

Numerical simulations of galaxy interactions suggest that mergers of spiral galaxies can lead to the formation of elliptical galaxies (e.g. Toomre & Toomre 1972; Barnes 1988, 1992; Hernquist 1992, 1993). These expectations are supported by optical imaging observations which show that several elliptical galaxies exhibit shells, ripples, arcs, counter-rotating cores, or faint tidal tails; structures which are interpreted as evidence for galaxy interactions or mergers (e.g. Schweizer et al. 1990).

The availability of high spatial resolution X-ray data allows us to investigate these indications using X-ray binary populations as probes of the star-formation histories of galaxies: High-Mass X-ray binaries (HMXBs) form much more efficiently than Low Mass X-ray binaries (LMXBs) (Kalogera & Webbink, 1998; Portegies Zwart & Verbunt, 1996), leading to a higher number of X-ray sources per star in young stellar populations. Thus, using number counts of discrete X-ray sources, we can identify regions of recent or enhanced star-formation within a galaxy. For example, the spatial distribution of LMXBs (which trace the old stellar populations) in early-type galaxies is expected to be smooth, generally following the distribution of optical star-light. On the other hand, if there are any sites of recent star-formation (e.g. triggered by galaxy interactions) they are expected to host HMXBs, the larger numbers of which (compared to the LMXBs) may result in an overall non-uniform X-ray source spatial distribution.

In this paper we present evidence for a non-uniform spatial distribution of X-ray sources in NGC 4697, an elliptical galaxy without any indication for merging activity from optical data, and we compare these results with similar observations of other early-type galaxies, finding a second example in NGC 4261. A more detailed investigation of the X-ray source populations, their degree of non-uniformity and links to other merger activity diagnostics will be presented in a forthcoming paper.

## 2 NGC 4261

NGC 4261 is a nearby (35.1 Mpc for $`\mathrm{H}_\mathrm{o}`$=75$`\mathrm{km}/\mathrm{s}`$) E2 type galaxy, belonging to a poor group of galaxies (Garcia, 1993) but without any evidence for interaction with other members of the group. Optical observations show a 20″($`3.4`$ pc) long dust lane along the North-South axis of the galaxy (Martel et al. 2000 and references therein) and a spectacular nuclear dust disk ($`2\mathrm{}`$; Jaffe et al. 1996, Ferrarese et al. 1996). Its optical isophotes in larger scales have a very strong boxy morphology (e.g. Peletier et al. 1990). Optical spectroscopy shows that its dominant stellar population has an age of $`15`$ Gyr (Trager et al. 2000). This galaxy also features two prominent radio jets (e.g. Birkinshaw & Davies 1985, Jones et al. 2000) emanating from an active nucleus. The mass of the supermassive nuclear black-hole is estimated to be $`(4.9\pm 1.0)\times 10^8M_{}`$ (Ferrarese et al. 1996).

### 2.1 Observations and Data Analysis

NGC4261 was observed for 35 ks (OBSID 834; PI M. Birkinshaw) with the ACIS-S3 CCD on-board the Chandra X-ray Observatory. The spatial resolution of Chandra is 0.5′′ (van Speybroeck et al. 1997) corresponding to a physical resolution of $`80`$ pc at the distance on the galaxy. Although the observation was performed in 1/2 subarray mode in order to mitigate pile-up on the nuclear source, the active field of view ($`4.1^{}\times 8.3^{}`$) was large enough to cover the whole galaxy ($`\mathrm{D}_{25}=4.07^{}`$; de Vaucouleurs et al., 1991). The initial processing of the data and the results from the analysis of the nuclear spectrum are presented in Zezas et al. (2003). In this paper we concentrate on the spatial distribution of the discrete sources.

We used the wavdetect algorithm, within the CIAO v2.2 data analysis suite in order to search for discrete sources in four different energy bands: full band (0.3-7.0 keV), soft band (0.3-2.0 keV) and hard band (2.0-7.0 keV). The limit of 7.0 keV was set in order to minimize contamination by the particle background, while the boundary of 2.0 keV was chosen in order to limit the contribution of the diffuse emission (which is mainly thermal with $`\mathrm{kT}=0.6\mathrm{keV}`$) only in the soft band. The source detection was performed in scales of 1, 2, 8 and 16 pixels and a limiting probability of $`1`$ chance detection over this field. We found a total of 62 sources within the S3 CCD. Of those, 45 sources are within the D25 radius of NGC 4261, 40 of which have significance higher than 3$`\sigma `$ above the local background. The absorption corrected luminosities of the latter are above the Eddington limit for an 1.4 Mneutron star, ranging from $`2.5\times 10^{38}`$ erg s-1 up to $`4.4\times 10^{39}`$ erg s-1 (with two sources above $`10^{39}`$ erg s-1) assuming a power-law model with $`\mathrm{\Gamma }=1.7`$ and Galactic line-of-sight $`\mathrm{H}_\mathrm{I}`$ column density ($`5.8\times 10^{20}`$ cm-2 ; Stark et al. 1992).

### 2.2 Spatial distribution of X-ray sources

Figure 1a shows an adaptively smoothed image of NGC 4261 in the full band, together with the D25 ellipse of the galaxy and the active field of the S3 ACIS-CCD. For reference in Figure 1b we present a DSS2 R-band image of the galaxy in the same scale, with the discrete sources marked by the best fit $`3\sigma `$ ellipses to the spatial distribution of their photons111http://asc.harvard.edu/ciao/download/doc/detect\_html\_manual/. From the X-ray image is clear that the spatial distribution of the discrete sources does not follow the distribution of optical star-light, but shows instead a distinct asymmetric pattern. In the following discussion we focus on the sources within the D25 ellipse.

In order to assess the probability that the non-uniform distribution may be observed by chance, we performed a Kolmogorov-Smirnov (KS) test to compare the distribution of the position angles (PA) of the sources (measured clockwise from the North-South direction), with (a) a uniform distribution and (b) the azimuthal profile of the optical light of the galaxy. For the latter we created a sample of “optical sources” based on the star-light distribution in the DSS2 R-band image, by randomly drawing position angles in 2 sectors from a flat parent distribution (equivalent to assuming constant surface brightness within each 2 sector). The number of draws in each sector was proportional to the fraction of the total galactic light included in it, after removing any strong point-like optical sources (most probably associated with foreground stars or background galaxies/AGNs). We found that both comparison distributions are inconsistent with the observed X-ray source azimuthal distribution above the 99.9% confidence level (Fig. 2a). This result holds even for model distributions with as few as 10 sources.

A search for irregularities in the azimuthal distribution of the PA of the sources, using the Bayesian blocks method (Scargle 1998; Scargle et al. 2003, in prep) also shows that there is an excess of sources between PA=140 and PA=190 at a confidence level greater than 99% (estimated from the prior of $`\gamma =4.5`$ used in the decomposition; however see Scargle et al. 2003, in prep222see also http://space.mit.edu/CXC/analysis/SITAR/functions.html, for a caveat on this interpretation of the value of the prior). Therefore, we conclude that the spatial distribution of the X-ray sources is not the same as that of the optical light tracing the stellar population.

## 3 Other examples: NGC 4697

Indication for a different distribution of X-ray sources and optical star-light has been also seen in NGC 720 (Jeltema et al. 2003), which has relatively young stellar populations ($`5`$ Gyrs; Trager et al. 2000). The X-ray sources in NGC 720 appear to trace out spiral arms, but performing the same analysis as above we confirm that this result is not statistically significant (Jeltema et al. 2003). Based on the KS test we find that the X-ray source distribution is different from the optical distribution or a constant at the 50% and 70% confidence levels respectively. Similarly, the Bayesian blocks method showed that any local enhancement of sources is significant only at the 63% confidence level.

We searched the Chandra archive for deep observations of nearby early-type galaxies, which also did not show any evidence for fine structure (Schweizer & Seitzer 1992) or recent merging activity, in order to construct a comparison sample. We obtained data for NGC3379, NGC4636 and NGC4697, the general properties of which are presented in Table 1. We analyzed the data for these three galaxies in the same way as for NGC 4261 (§2.1). A KS test between the X-ray source distribution and the comparison distributions (calculated as described in §2.2) showed that only in the case of NGC 4697 we can rule out at the 98% confidence level (Fig. 2b) that they belong to the same population ($`95\%`$ probability that there is a local enhancement of sources based on Bayesian blocks). In the other two cases (NGC 3379 and NGC 4636) these confidence levels are $`30`$% and $`25`$% respectively, not allowing us to draw any definitive conclusions.

## 4 Discussion and Conclusions

Our results show evidence for a non-uniform spatial distribution of the X-ray source population in two nearby apparently normal elliptical galaxies. These galaxies have very low fine-structure parameters ($`\mathrm{\Sigma }=1`$ for NGC 4261 and $`\mathrm{\Sigma }=0`$ NGC4697) indicating that if they are merger products the merger event took place at least a few Gyrs ago. This is because merger simulations indicate that ripples, shells or strong tidal tails are not easily observable for much longer than $`10^9`$ yr after the relaxation of the merged system (e.g. Quinn 1984; Hernquist & Quinn 1988, 1989). Moreover, faint traces of tidal tails surviving for much longer might not be easily detectable in optical observations because of their very low surface brightness (Mihos 1995). In fact, optical observations indicate stellar populations as old as 9 and 15 Gyrs in the nuclear regions of NGC 4697 and NGC 4261 respectively (Trager et al. 2000).

We propose that one or more localized star-formation events which occurred at most a few hundred Myr ago, are responsible for the discrepancy between the optical and X-ray morphology of these two galaxies. According to X-ray binary formation models, LMXBs are more susceptible to effects such as supernova kicks and common envelope phases (Kalogera & Webbink, 1998) than HMXBs, leading to HMXB formation efficiencies 10 to 100 times higher than for LMXBs (Portegies-Zwart & Verbunt 1996). In fig. 3 we plot the $`\mathrm{L}_\mathrm{X}/\mathrm{L}_\mathrm{B}`$ ratios of star-forming galaxies from the Einstein survey of Shapley, Fabbiano & Eskridge (2002) together with the mean $`\mathrm{L}_\mathrm{X}/\mathrm{L}_\mathrm{B}`$ ratio for the discrete X-ray sources in elliptical galaxies ($`\mathrm{log}(\mathrm{L}_\mathrm{X}/\mathrm{L}_\mathrm{B})=3.05\pm 0.4`$333where $`\mathrm{L}_\mathrm{X}`$ is in erg s-1, in the 0.35-10.0 keV band, $`\mathrm{L}_\mathrm{B}`$ is the total B band luminosity in erg/s; Athey 2003). From this plot is clear that even after accounting for a 50% contribution from the diffuse emission in star-forming galaxies (e.g. Zezas et al. 2001), the latter have systematically higher $`\mathrm{L}_\mathrm{X}/\mathrm{L}_\mathrm{B}`$ ratios than the X-ray binary component of early type galaxies. This together with the fact that individual HMXBs have significantly lower X-ray to B-band flux ratios than LMXBs (van Paradijs & McClintock 1995) supports the picture that HMXBs are forming much more efficiently than LMXBs.

Recently Barnes (2003) suggested that shock-induced star-formation can explain the star-forming activity observed along the tidal tails of “The Mice”(Stockton 1974; de Grijs et al. 2003). According to this model shock waves developing along the tidal tails can compress the neutral gas and trigger star-formation. This picture is consistent with the tail-like distribution of the observed X-ray sources in NGC 4261, if the tidal tails are projected against the body of the merger remnant. If the age of this young stellar population is less than a few hundred Myr, it forms X-ray binaries much more efficiently than the populations in the relaxed merger remnant. On the other hand its optical emission is diluted by the optical light of the much stronger old population making its detection in the optical band very difficult. Depending on the strength of the star-formation, this may result in an overall projected spatial distribution of X-ray sources which is inconsistent with that of the star-light. This star-formation event can be well-approximated by an instantaneous burst, in which case the young X-ray source populations are expected to decay in a few hundred Myr after the passage of the shock. Therefore, we estimate that the shock should propagate with a speed of $`100`$ km s-1in order to cover the length of the region we observe X-ray sources (20.5 kpc), within the timescale of HMXB formation ($`200`$ Myr). This speed is realistic for shock waves developing in the interstellar medium.

Alternatively, N-body simulations show that structures resembling dwarf galaxies can form in the tidal tails as a merger completes (Barnes & Hernquist 1992, 1996). These structures remain bound to the body of the remnant, orbiting it. Because the objects formed in this manner have a range of binding energies in the potential well of the remnant, the characteristic time for them to fall back onto the remnant can be much longer than the time for the body of the remnant to relax (Hernquist & Spergel 1992; Mihos & Hernquist 1996). Thus, depending on the orbital distribution, the remnant will eventually lose evidence for a merger, while tidal dwarfs continue to orbit at large radii. The interaction between the body of the merger and the condensations in the tidal tails can trigger small scale star-formation events and locally enrich the galaxy with a young stellar population.

As was mentioned earlier, a young stellar population can form X-ray binaries very efficiently, while it is very difficult to detect it against the much brighter stellar populations of the merger remnant. Therefore, depending on the strength of the star-formation, this may result in an overall spatial distribution of X-ray sources which is inconsistent with that of the star-light. Since fall-back is estimated to occur for several Gyrs after relaxation, we expect that even elliptical galaxies with stellar populations of $`10`$ Gyrs may exhibit non-uniform X-ray source distribution, if they are the end-points of galaxy mergers.

Although the spatial distribution of the sources in NGC 4261 indicates that they are associated with the tidal tails, probably both mechanisms can produce populations of numerous, young X-ray binaries in merger remnants. Given the timescales of the orbits of dwarf galaxies (up to 1 Gyr) or the velocities of the shocks along the tidal tails, and the lifetime of HMXBs (up to $`100`$ Myr), we would expect long periods during which these galaxies have uniform X-ray source distributions. These periods are expected to be longer in the more evolved systems since the dwarf companions with short orbits are accreted first, which is consistent with our finding that this phenomenon appears in only a few of the most evolved galaxies examined. Studies of a more complete sample of elliptical galaxies and merging galaxies in their latest stages of merging will allow us to further investigate this hypothesis.

We thank the referee (J. Rose) for a very helpful report. We also thank J. Barnes for providing results prior to publication, and E. O’Sullivan and A. Athey for useful discussions on the X-ray emission of X-ray binaries in elliptical galaxies. We thank P. Ratzlaff and J. J. Drake for making available their code for Bayesian blocks analysis. This work has been partly supported by NASA Grants G01-2116X, G01-3150X, and ATP NAG5-12140 and NSF grants ACI 96-19019, AST 98-02568, AST 99-00877, and AST 00-71019.


text:
# Brownian motion approach to the ideal gas of relativistic particles

## Abstract

The relativistic generalization of a free Brownian motion theory is presented. The global characteristics of the relaxation are explicitly found for the velocity and momentum (stochastic) kinetics. It is shown that the thermal corrections, to the both relaxation times $`T`$ (of stationary autocorrelations) and transient relaxation time of momentum, appear slowing down the processes. The transient relaxation time of the velocity does not depend explicitly on temperature, $`T(v_0)=m(v_0)/\gamma \epsilon _0/\gamma c^2`$, and it is proportional to the initial energy of a relativistic Brownian particle.

During the last twenty years the Brownian motion approach has become the most important tool for discovering and studying the spectacular phenomena, the critical slowing down , noise–induced transitions , stochastic resonance , ratchet dynamics , and resonant activation , in (nonlinear) nonequilibrium systems. In the most of thermodynamical applications the nonlinear Langevin equation appears as an overdamped limit of a dissipative dynamic in a certain external potential . The important conclusion of the present paper is that, within the Markovian diffusion theory, the nonlinear stochastic kinetic results already on the basic level of (spatially homogeneous) ideal gases description if the relativistic particles are considered. Particularly the relativistic free Brownian motion theory based on identical general assumption as the original one \[i.e., including the interaction with the thermal bath by a systematic force, the dissipative kinematic friction proportional to the velocity $`\gamma v`$, and an additive random force, the thermal fluctuation, given by Gaussian white noise $`\xi _t`$\]

$$\dot{p}_t=\gamma v+\xi _t,$$
(1)
where \[$`c=1`$ is used throughout the paper\] $`p=mv/\sqrt{1v^2}`$ is a relativistic momentum of a Brownian particle, $`v=p/\sqrt{p^2+m^2}=p/\epsilon `$, $`\epsilon =\sqrt{p^2+m^2}`$, and $`\xi _t\xi _s=2D\delta (ts)`$, is consistent both with the thermodynamical requirements of the equilibrium Gibbs ensemble theory and with the special relativity requirements of transformational properties of the momentum probability density distribution.

First note that the general Kramers–Fokker–Planck kinetic equation (see, e.g., the equation (26.26) of the Ref. ) of a Markovian diffusion description simplifies to the form

$$_tP=_pL_{22}[(_pH_1)P+k_BT_pP]$$
(2)
for the spatially homogeneous system of free noninteracting particles \[of the Hamiltonian $`H_1(p)`$\]. The kinetic coefficient $`L_{22}=L_{22}(p,T)`$ is in principle determined by the properties of the interaction with the thermal bath and do not depend on $`p`$ for the potential interaction. For the latter case (of the state independent diffusion) the Langevin description Eq. (1), and the corresponding Fokker–Planck one , Eq. (2), are equivalent in view of the identification $`L_{22}=\gamma =\gamma (T)`$, $`D=\gamma k_BT`$ (the Einstein relation, or the fluctuation–dissipation theorem), and $`_pH_1=v`$. The Eq. (2) is the most general form of the diffusion in the momentum space leading to the Gibbs–Boltzmann stationary state $`W(p)e^{\beta H_1(p)}`$ of thermal equilibrium ($`\beta =1/k_BT`$). The corresponding (nonlinearly coupled) Langevin equation (in the Stratonovich interpretation ) is

$$\dot{p}_t=\gamma g^2v+Dgg^{}+g\xi _t,$$
(3)
where $`L_{22}=\gamma g^2`$ and $`g^{}=g/p`$.

The (any) corresponding probability density distributions $`W(p)`$ and $`W^{}(p^{})`$ in a resting and moving (with a constant velocity $`V`$) reference frame are related by non-covariant formula $`\epsilon W(p)=\epsilon ^{}W^{}(p^{})`$, where $`\epsilon `$ and $`p`$ should be expressed as functions of $`\epsilon ^{}`$ and $`p^{}`$, according to the Lorentz transformation . Because the time increments (counted from the initial preparation of the system) are given by $`t^{}`$ and $`t=t^{}/\sqrt{1V^2}`$, respectively, the relation

$$P^{}(p^{},t^{})=\frac{1+Vv^{}}{\sqrt{1V^2}}P(\frac{p^{}+V\epsilon ^{}}{\sqrt{1V^2}},\frac{t^{}}{\sqrt{1V^2}})$$
(4)
is required. One has $`v=(v^{}+V)/(1+Vv^{})`$, $`\xi _t=(1V^2)^{1/4}\xi _t^{}`$, and, using the Stratonovich calculus, $`dp/dt=(1+Vv^{})dp^{}/dt^{}`$, thus the transformed stochastic kinetic (1) takes the form

$$\frac{dp^{}}{dt^{}}=\gamma \frac{v^{}+V}{(1+Vv^{})^2}+\frac{(1V^2)^{1/4}}{1+Vv^{}}\xi _t^{},$$
(5)
where $`v^{}=p^{}/\sqrt{p_{}^{}{}_{}{}^{2}+m^2}`$. It is easy to verify that the r.h.s. of (4) solves the Fokker–Planck equation corresponding to the Eq. (5) if (and only if) $`P(p,t)`$ satisfies the proper Eq. (2) related to the kinetic (1). The same conclusion applies also to the general case (3).

The Eq. (1), $`\dot{p}_t=\gamma p/\sqrt{p^2+m^2}+\xi _t`$, for the momentum or, equivalently, the (Stratonovich) equation for the velocity,

$$\dot{v}_t=\gamma v(1v^2)^{3/2}/m+\xi _t(1v^2)^{3/2}/m,$$
(6)
thus provides the relativistic generalization of the diffusion Markovian description of equilibration for noninteracting free particles system . The normalized stationary distributions are

$$W(p)=\frac{1}{2mK_1(m\gamma /D)}e^{(\gamma /D)\sqrt{p^2+m^2}}$$
(7)
$$\mathrm{and}W(v)=\frac{(1v^2)^{3/2}}{2K_1(m\gamma /D)}e^{(\gamma /D)m/\sqrt{1v^2}},$$
(8)
where $`\gamma =D\beta `$. $`K_n`$ is a modified Bessel function. The second moments read

$$p^2=\frac{m}{\beta }\frac{K_2(m\beta )}{K_1(m\beta )},v^2=1_{m\beta }^{\mathrm{}}𝑑z\frac{K_0(z)}{K_1(m\beta )},$$
(9)
so the (even) moments of the velocity, as incomplete integrals of Bessel functions, are not given in closed analytical form.

The nonlinear kinetic (1) or (6) does not belong to the class of the known solvable models and the time–dependent solution is not known. Nevertheless, there are still the quantities, useful to characterize the relaxation of the nonlinear systems, which can be expressed by quadratures of the (necessarily autonomous) drift and diffusion coefficients, without solving the kinetic equation. These are the complete (over time) integrals of (various) transient (single-time) moments  and stationary (two-times) correlation functions . In fact, making use of $`P(x,t|x_0)=e^{t\widehat{L}}\delta (xx_0)`$, where $`\widehat{L}(x_0)`$ is a backward Fokker–Planck operator , and computing $`X(x_0)=_0^{\mathrm{}}𝑑tx_t_{x_0}`$, where the integrand is a conditional moment equal $`x_0`$ at the initial time $`t=0`$, by a formal integration $`_0^{\mathrm{}}𝑑te^{t\widehat{L}}=\widehat{L}^1`$ one obtains

$$\widehat{L}X=x_0,$$
(10)
which is a nonhomogeneous linear equation of the first order for the function $`dX/dx_0`$. The $`x_{\mathrm{st}}=0`$ is required for the convergence of $`X`$ and the condition $`X(0)=0`$ determines a first constant of integration for Eq. (10). The second one should be chosen in such way that the correct (deterministic) result is recovered when the noise strength $`D0`$. The Eq. (10) differs from the Pontryagin equation for the mean first passage time  in nonhomogeneous term ($`x_0`$ instead of $`1`$). The integrals of stationary autocorrelations are then given by $`S_x=_0^{\mathrm{}}𝑑tx_tx_0_{\mathrm{st}}=𝑑xW(x)X(x)x`$. The last quantity (i.e., the power spectrum function taken at zero frequency $`S_x=S_x(0)`$), is called, after some normalization (usually divided by the stationary variance) resulting in the proper (time) dimension, the relaxation time $`T`$ . The analytical expression for this quantity has been first obtained by Jung and Risken . Then, by a similar consideration, the formulae for the quantities of the former type, called generally nonlinear relaxation times, have been found .

Identifying the appropriate coefficients from Eqs. (1) and (6), and carrying out the integrals of the Jung and Risken formula  we obtain the exact results

$$S_v=D/\gamma ^2,$$
(11)
$$S_p=Dm^2/\gamma ^2+\frac{3mD^2}{\gamma ^3}\frac{K_2(m\gamma /D)}{K_1(m\gamma /D)}D^3/\gamma ^4.$$
(12)
The same Eq. (11) is obtained for the ordinary Brownian motion, $`S_p^o=m^2S_v^o=Dm^2/\gamma ^2`$. Introducing $`T_x=S_x/x^2`$ and using $`m`$ and $`m/\gamma `$ as units of momentum and time \[which is equivalent with considering the generic forms of Eqs. (1), (6) with a single (implicit) parameter $`D=\tau `$ and $`m=\gamma =1`$\] one gets for dimensionless quantities

$$T_v=\left\{_0^1𝑑v\frac{\mathrm{exp}[(\tau \sqrt{1v^2})^1]}{K_1(\tau ^1)}\right\}^1,$$
(13)
$$T_p=\frac{1\tau ^2+3\tau K_2(\tau ^1)/K_1(\tau ^1)}{K_2(\tau ^1)/K_1(\tau ^1)},$$
(14)
where $`\tau `$ is a dimensionless temperature, $`\tau =D/m\gamma =k_BT/mc^2`$. Note that (13) is the alternative form of $`T_v=\tau /v^2`$, where $`v^2`$ is given by Eq. (9), which is more suitable for numerical computations. The respective quantities of the ordinary Brownian motion are simply equal to the inverse of the relaxation rate of a linear model (1), $`T_v^o=T_p^o=1`$, and do not depend explicitly on temperature . In contrast both $`T_v`$ or $`T_p`$ increases with temperature, see Fig. 1, approaching the asymptotic $`T_v(\tau )\tau +\pi /2`$ or $`T_p(\tau )(5/2)\tau `$, respectively. The former asymptotic law follows from representation (9) \[the complete integral of $`K_0`$ is equal to $`\pi /2`$\], the latter is obtained by the use of the asymptotic formula $`K_n(1/\tau )(n1)!(2\tau )^n/2`$ (for $`\tau \mathrm{}`$)  in Eq. (14). For small $`\tau 0`$ $`T_v(\tau )1+(3/2)\tau `$ and $`T_p(\tau )1+3\tau `$. So, the decay of both (stationary) autocorrelations in a relativistic case proceed slower in higher temperatures, and slower than for the ordinary Brownian motion. We want to stress however that, in the case of nonlinear system, the relaxation time $`T_x`$ should not be interpreted as a characteristic of the equilibration process from a certain nonstationary state. One has in particular $`T_x(\tau =0)=1`$, which would be misunderstood that the (deterministic) evolution of the relativistic and classical system proceeds with the same (or at least similar) speed. Meanwhile it may be generally proved that at the limit $`\tau 0`$ only the lowest order terms of the drift and diffusion coefficients contribute to $`T_x`$. Thus, such a conclusion applies at most to the relaxation from the close to the equilibrium states.

In order to compare the nonstationary behavior let us define the transient relaxation time

$$T_\tau (x_0)=x_0^1_0^{\mathrm{}}𝑑tx_t_{x_0},$$
(15)
which is related to the quantity $`X`$, Eq. (10), via $`x_0T_\tau (x_0)=X_\tau (x_0)`$. Then $`P_\tau (p_0)`$ satisfies

$$p_0P_\tau ^{}/\sqrt{1+p_0^2}+\tau P_\tau ^{\prime \prime }=p_0.$$

The particular solution of the nonhomogeneous equation is $`P_\tau ^{}=\sqrt{1+p_0^2}+\tau `$, whereas the (general) solution of the homogeneous equation is singular at $`\tau =0`$. Thus,

$$T_\tau (p_0)=[p_0\sqrt{1+p_0^2}+\mathrm{log}(p_0+\sqrt{1+p_0^2})]/2p_0+\tau .$$
(16)
Similarly,

$$[v_0(1v_0^2)^{3/2}3\tau v_0(1v_0^2)^2]V_\tau ^{}+\tau (1v_0^2)^3V_\tau ^{\prime \prime }=v_0$$

has a (regular at $`\tau =0`$) particular solution $`V_\tau ^{}=(1v_0^2)^{3/2}`$, so

$$T_\tau (v_0)=1/\sqrt{1v_0^2}\epsilon _0.$$
(17)
The deterministic (generic) Stokes equations, $`\dot{p}=p/\sqrt{p^2+1}`$ and $`\dot{v}=v(1v^2)^{3/2}`$, are solvable by quadratures. The solution, with the initial condition $`x_0`$, is

$$t(x)=t_x(x)t_x(x_0),$$
(18)
where $`t_p(p)=\sqrt{1+p^2}+\mathrm{log}|(1+\sqrt{1+p^2})/p|`$, and $`t_v(v)=(1v^2)^{1/2}+\mathrm{log}|(1+\sqrt{1v^2})/v|`$, respectively. One can verify that the direct calculation of $`_0^{\mathrm{}}𝑑tx(t)/x_0=_{x_0}^0𝑑xt(x)/x_0`$ yields $`T_0(x_0)`$ as given by Eq. (16) or (17). The previous results (11) and (12) are also recovered when $`T_\tau (x_0)`$ is integrated with $`x_0^2W(x_0)`$, where $`W`$ is a stationary probability density distribution (8) or (7), respectively. Both (dimensionless) transient relaxation times for nonrelativistic Brownian motion are equal to unity. The relativistic relaxation of momentum depends both on the temperature and the initial state. $`T_\tau (p_0)`$, Eq. (16), is an increasing function of both arguments and has the following asymptotic $`T_\tau (p_0)1+\tau +p_0^2/6`$, for $`p_00`$, and $`T_\tau (p_0)|p_0|/2+\tau `$, for $`p_0\mathrm{}`$. The result for the velocity is somehow unexpected. The Eq. (17) shows that indeed the relativistic relaxation proceeds slower, however does not depend on temperature. It turns out that with the increase of temperature the early stage of evolution becomes faster in such a way that compensates the subsequent slower long-time behavior, see Fig. 3. The transient relaxation time $`T_\tau (v_0)`$ is equal to the (dimensionless) initial energy of the relativistic Brownian particle $`\epsilon _0`$. The plots of $`T_0(x_0)`$ are shown in Fig. 2.

The equation (1) (in a generic form) has been solved numerically by the Runge–Kutta method for $`10^4`$ sample realizations of the white Gaussian noise, for different initial conditions and temperatures. $`p_t`$ and $`v_t`$ have then been computed as the arithmetical average of appropriate values obtained for different trajectories of the noise. The typical curves are plotted in Figs. 3 and 4, together with the corresponding deterministic results (18). The solution of the relativistic Stokes equation for velocity changes qualitatively, when $`v_0`$ exceeds $`1/2`$. In fact below this value of the initial state the function $`v(t)`$ is convex down for all $`t`$ \[as well as (always) $`p(t)`$\]. However for $`v_0>1/2`$ at the early stage of evolution the $`v(t)`$ appears convex up \[in order to achieve the required asymptote $`v(t)=1=c`$ for $`v_01`$\]. The former subrelativistic case, with particular $`v_0=2/5`$ ($`p_00.436`$), and the latter, with $`p_0=4`$ ($`v_00.970`$), are presented.

The results of the paper may be summarized as follows. We have shown that using the relations of the relativistic dynamics we can generalize the Brownian motion theory in a way consistent with thermodynamical requirements. The general solvable equation (10) for the quantities simply related to the transient relaxation times (15) has been found. These, as well the stationary relaxation times $`T`$, have been found analytically for the relativistic momentum and velocity (stochastic) kinetics, Eqs. (13, 14, 16, 17). All these quantities (considered as dimensionless) are equal unity for the ordinary Brownian motion. The relaxation times $`T_x(\tau )`$, by the normalization equal unity in the deterministic limit $`\tau 0`$, exhibit positive first order corrections for small temperatures, becoming just proportional to $`\tau `$ in a high temperature asymptotic. Thus the decay of correlations in a stationary state proceed slower than for the nonrelativistic Brownian motion. The transient relaxation time \[describing the relaxation from a given initial state $`x_0`$, and by the normalization equal unity for the deterministic system set initially at the equilibrium, $`T_0(0)=1`$\] of momentum $`T_\tau (p_0)`$ is an affinic function of temperature. The asymptotic behavior for the large argument(s) reads $`T_\tau (p_0)|p_0|/2+\tau `$. The transient relaxation time for the velocity (17) does not depend on temperature, and it is equal to the dimensionless initial energy of the relativistic Brownian particle. The thermal corrections to the considered global characteristics of the relaxation are linear with the slope of the unity order and appear (even) for the noninteracting free relativistic Brownian particles system. The corrections can be observed for sufficiently high physical temperatures. For, e.g., hydrogen particle (proton) the dimensionless quantity, which characterize the deviation from the ordinary nonrelativistic behavior, $`\tau =k_BT/mc^210^{13}T`$\[K\]. So for the temperature $`T1\mathrm{M}\mathrm{e}\mathrm{V}10^{10}`$K the relative correction is of the order $`0.1`$%. The correction to transient relaxation time of the velocity has the kinematic origin only and for $`T_\tau (p_0)`$ both kinds of corrections are decoupled. Note that the general dependence on temperature in physical units, introduced (multiplicatively) by the coefficient $`\gamma (T)`$, should be counted separately .

The three-dimensional generalization of the proposed diffusion description is obvious, however in contrast to the ordinary Brownian motion theory it leads to the set of the three coupled ($`\stackrel{}{v}=\stackrel{}{p}/\epsilon =\stackrel{}{p}/\sqrt{m^2+\stackrel{}{p}^2}`$) Langevin equations (with statistically independent noises in each direction) or equivalently to the (forward) Fokker–Planck description  \[$`=/\stackrel{}{p}`$\]

$$_tP(\stackrel{}{p},t|\stackrel{}{p_0})=\gamma (\stackrel{}{p}P/\epsilon )+D^2P.$$
(19)
The required form of the (Gibbs-Boltzmann) stationary distribution is still recovered, however the calculations of relaxation times cannot be analytically completed in this case.