text:
# Kinematics of stock prices

(4 September 2002, subm. Phys. Rev. Lett.)

## Abstract

We investigate the general problem of how to model the kinematics of stock prices without considering the dynamical causes of motion. We propose a stochastic process with long-range correlated absolute returns. We find that the model is able to reproduce the experimentally observed clustering, power law memory, fat tails and multifractality of real financial time series. We find that the distribution of stock returns is approximated by a Gaussian with log-normally distributed local variance and shows excellent agreement with the behavior of the NYSE index for a range of time scales.

Remarkable progress has been made in quantitatively describing nonstationary and non-Gaussian phenomena, including those observed in economic econophysicsbook  and social systems. The behavior of financial markets has recently econophysicsbook ; prl1 ; plamen1 ; lrc1 ; lrc-old0 ; gopi ; lrc2 ; mantegna2 ; mantegna  become a focus of interest to physicists as well as an area of active research because of its rich and complex bacchelier ; ding ; crato ; bouchaud ; levy ; baillie ; galluccio ; vandewall ; Liu ; volatility-stanley  dynamics. The daily returns $`r_t`$ for a given stock can be defined as

$$r_t=\mathrm{log}_{10}\frac{S_{t+1}}{S_t},$$

at time $`t`$ (Fig. 1), where $`S_t`$ is the price of a an asset or the value of an index or a currency exchange rate. In this paper we address the important yet unsolved problem of how to model the kinematics of stock prices without handling the problem of the causes of motion. Our point of departure is the family of ARCH-GARCH models. These families of stochastic kinematic models neglect the causes underlying the price variations and focus only on the equations of motion governing the fluctuating returns. Specifically, the interacting buying and selling processes are not considered. ARCH models, introduced by Engle in 1982 engle , are stochastic Auto-Regressive Conditional Heteroscedasticity models, characterized by zero mean and nonconstant variances dependent on the past. Such models can simultaneously have global stationarity as well as local nonstationary behavior. The simplest ARCH model of stock returns can be defined by

$$\sigma _t^2=\alpha _0+\alpha _1r_{t1}^2,$$

with returns defined as

$$r_t=\sigma _t\omega _t,$$

where the random Gaussian variable $`\omega _t`$ has unitary variance and zero mean, while $`\alpha _0`$ and $`\alpha _1`$ are tunable parameters. Hence, the returns $`r_t`$ are Gaussian distributed with zero mean and local variance $`\sigma _t^2.`$ Only very recently has the physics community begun to study such models, although they are common in the economics literature. In more general ARCH processes, the local variance can depend not only on the previous value $`r_{t1}`$, but on any finite number $`n`$ of previous values $`r_{t1}\mathrm{}r_{tn}`$. GARCH models are a further generalization in which the local variance can depend not only on previous values of the returns but also on the previous values of the locally measured variances.

ARCH-GARCH models have succeeded in capturing important features (such as volatility clustering) of real financial data. They are, therefore, widely used, both in the financial community and in econophysics to model the kinematics of price variations. Nevertheless, several characteristic features of the real data are not well described. The three most important (and difficult to model) features not well accounted in ARCH-GARCH models are (i) the time scaling of the probability density distribution of the returns, (ii) the known long-range power law correlations in the volatility (the latter being a measure of the local standard deviation of the returns), and (iii) the multifractality of the returns, i.e., volatility power law correlation exponents are non-unique and depend on the magnitude of the events considered. As an alternative model, Lévy mandel  and truncated (see ref. tlf ) Lévy distributions have been proposed to fit the observed fat tailed distributions. While Lévy processes can correctly reproduce the time scaling of distributions, they cannot explain one of the most relevant and characteristic phenomena: the multifractal long-range correlations in volatility which has been found to be responsible for clustering of volatility and the persistence of fat tails for long time lags gopi ; gandhi-condmat2002 .

Here we propose a new model that does not suffer from these drawbacks. In order to test the model, we compare the model with a typical dataset: the New York Exchange (NYSE) daily composite index price closes from January 1966 to September 2001 (a total of some 9000 data points).

The model is inspired in the recent finding that the probability density distribution of the local variance $`\sigma _t`$ is similar to a log-normal lrc1  . Other key findings include the long-range correlations found in the absolute value of the returns and the multifractal behavior of the returns lrc2 . Our proposal, based on these findings, is the following map for the evolution of the volatility:

$$\sigma _{t+1}=e^{[a+b\omega _{t+1}]}(\overline{\sigma _t})^d,$$
(1)
with returns defined as

$$r_t=\sigma _t\stackrel{~}{\omega }_t,$$
(2)
where $`\stackrel{~}{\omega _t}`$ and $`\omega _t`$ are independently distributed unitary Gaussian variables with zero mean that are also independent from each other. Note that the independence of these two variables contrasts with ARCH models, for which they coincide. We discuss why this is preferable below. The tunable constants $`a`$, $`b`$ and $`d`$ are related to the nearly log-normal form of $`\sigma _t.`$ The value of $`d`$ must be close but smaller than 1 to guaranty the stationarity of the returns; $`a`$ is related to the typical size of the daily volatility, and $`b`$ to the typical size of its fluctuations. The term $`\overline{\sigma _t}`$ represents the average with respect to last $`T`$ days, from $`tT`$ to $`t`$. (We have used below a simple unweighted average but, in principle, we could have used any type of moving average, for instance a power law or exponentially weighted moving average.)

Using this map, we generate a data set of $`9000`$ returns (Fig. 1(b)) with the following choice of constants: $`a=0.1`$, $`b=0.1`$, $`d=0.98`$ and $`T=10`$. We then compare the statistical and scaling properties of the time series generated with the real data set. Note that the above model can reproduce very well the clustering of volatility, as is evident comparing Fig. 1(a) with Fig. 1(b). Nevertheless, we proceed with a quantitative comparison of both the scaling of the distribution and the multifractal power law correlations.

We next study the probability density of returns for the proposed model. We estimate the probability density of returns using

$$P(r)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{2\pi }\mathrm{\Delta }}\mathrm{exp}\{\frac{(rr_i)^2}{2\mathrm{\Delta }^2}\}$$

with a smoothing window of $`\mathrm{\Delta }=0.001`$ (Fig 2a). The $`r_i`$ can be either the real or artificial data and $`N`$ is their number(about 9000 for both). The dashed and dotted lines in Fig. 2(a) show the distributions of daily (one business day) returns, for both real and artificial data. Typical of financial times series is their invariance under rescaling of time. Therefore, we estimate with the same smoothing window the distribution for monthly (25 business days) returns (Fig 2b)

$$\frac{_{i=1}^{25}r_{i+t}}{(25)^\delta }=\frac{1}{(25)^\delta }\mathrm{log}\frac{S_{t+25}}{S_t},$$

using the measured value $`\delta =0.53`$ of the scaling exponent both for real and artificial data. Note that the theoretical curve (solid line, see discussion below) is exactly the same in the two figures, therefore both real and artificial data exhibit almost perfect time scale invariance of the return distributions (see also gandhi-condmat2002 ). The lack of agreement for small monthly returns is due to short-range (1 day) correlations in the signs of the returns that are neglected in our model, but that could easily be incorporated gandhi-condmat2002 . Compared to other models of stock returns, this model shows remarkably good agreement with real data.

We next compare these probability distributions of real and artificial data with a “theoretical” prediction (solid line in Figs. 2). The model proposed leads to a distribution that is a Gaussian $`P_\text{t}`$ with log-normally distributed local variance. We can find $`P_\text{t}`$ by convolving Gaussians of varying widths:

$$P_\text{t}(r)=\rho (\sigma )\frac{e^{\frac{r^2}{2\sigma ^2}}}{\sqrt{2\pi }\sigma }𝑑\sigma $$
(3)
$$\rho (\sigma )=\frac{1}{\sqrt{2\pi }s\sigma }e^{\frac{1}{2}(\frac{\mathrm{log}\sigma m}{s})^2}$$
(4)
We have used the empirically found values $`s=0.41`$ and $`m=0.34`$. The solid lines in Figs. 2 leave no doubt that this distribution is one of the best candidates—if not the best—for describing the distribution of price variations. The agreement found is exceptionally strong.

It is known that daily returns have no auto-correlations for lags larger than a single day, consistent with the efficient market hypothesis. This fact can be also checked by using Detrended Fluctuation Analysis (DFA) pengdfa  and related methods. Consider the cumulative returns $`\varphi _t(L)`$, defined as the sum of $`L`$ successive returns divided by $`L`$:

$$\varphi _t(L)=\frac{1}{L}\underset{i=1}{\overset{L}{}}r_{t+i}$$
(5)
One can define $`N/L`$ non overlapping variables of this type, and compute the associated variance $`\sigma ^2(\varphi (L))`$, where $`N`$ 9000 is the number of data points. According to the central limit theorem, uncorrelated (or short-range correlated) $`r_t`$ would lead to power-law behavior:

$$\sigma ^2(\varphi (L))L^\alpha ,$$
(6)
with exponent $`\alpha =1`$ for large $`L`$. The exponent $`\alpha `$ both for the NYSE index and the model proposed here is near $`1`$ (see Fig. 3), confirming that returns are uncorrelated.

The lack of correlations does not hold true for quantities related to the absolute returns. In order to perform the appropriate scaling analysis, we introduce the generalized cumulative absolute returns defined as the sum of $`L`$ successive absolute return powers $`|r_t|^\gamma ,\mathrm{},|r_{t+L1}|^\gamma `$, divided by $`L`$

$$\varphi _t(L,\gamma )=\frac{1}{L}\underset{i=1}{\overset{L}{}}|r_{t+i}|^\gamma $$
(7)
where $`\gamma `$ is a real exponent (noting that these quantities do not overlap.) Using this method, if the autocorrelation for powers of absolute returns exhibits a power-law with exponent $`\alpha (\gamma )1`$ for large $`L`$, then we expect

$$\sigma ^2(\varphi (L,\gamma ))L^{\alpha (\gamma )}.$$
(8)
(Note, however, that if the $`|r_t|^\gamma `$ are short-range correlated or power-law correlated with an exponent $`\alpha (\gamma )>1`$, then we would not detect anomalous scaling in the analysis of variance, because it is not possible to detect $`\alpha >1`$ using the method discussed here).

We find that both the real data and the model show similar multifractal behavior with non-unique scaling exponents, i.e., $`\alpha (\gamma )`$ constant. Note that the function $`\alpha (\gamma )`$ is not universal lrc1 ; lrc2  but depends on the particular asset considered. Moreover, the values $`\alpha (\gamma )`$ we compute do not coincide for real and artificial data, nevertheless this is not a drawback if one takes into account the non-universality. For completeness, we note that in order to obtain full agreement for the multifractal exponents of any specific economic series one must generalize the volatility as follows:

$$x_{t+1}=e^{a+b\omega _{t+1}}(\overline{x_t^c})^d.$$
(9)
where the exponent $`c`$ in the volatility average would allow for a finer control of the effects caused by extremal events. Also the resulting returns could be written as

$$r_t=(\alpha +x_t)\stackrel{~}{\omega }_t\sigma _t\stackrel{~}{\omega }_t.$$
(10)
with $`\alpha `$ representing an intrinsic constant contribution to daily volatility.

We now comment on the motivation for using two independent Gaussian variables, one for returns and one for volatility. Indeed, this choice separates the temporal evolution of volatility and returns. Besides giving reasonable agreement with experimental data, there is a deeper underlying motivation for this approach. In ARCH-GARCH models the present-day volatility depends on the previous day’s absolute return. Any trader knows that today’s sentiment about volatility does not depend on the previous day’s variation of price but rather on other things like previous day intraday volatility, implicit volatility in derivative products and public (or less public) information. Therefore, the previous day’s absolute return (which can be small after a fearful market day with enormous variation of prices) does not directly influence the present day volatility and so the evolution of the latter should be kept separate, exactly as we have done in the model proposed.

In conclusion, the proposed stochastic multifractal model of long-range correlated absolute returns is able to reproduce features of real financial data that are not well accounted for by existing models. Specifically, our model is able to overcome the inability of ARCH-GARCH and Lévy models to adequately explain (i) the fat tailed distribution of the returns and its time scaling, (ii) the long-range power law volatility correlations, or (iii) the multifractality of the returns. The model presented here also represents an advance due to the exceptionally improved agreement with real data. We hope that this advance in modeling the kinematics of financial time series further contributes to the emerging study of econophysics and towards a better understanding of wider classes of complex systems.

We thank Michele Pasquini for useful discussions and for a critical reading of the manuscript. This work was supported by the Brazilian research agencies CAPES and CNPq and by the Alagoas State agency FAPEAL.


format: text

text:
# Microlensing of Strongly Interacting Binary Systems

## 1 Introduction

It is well known that about one half of the stars populating our Galaxy are part of gravitationally bound multiple systems. On the other hand, binary systems acting as microlenses produce drastic modifications into microlensing light curves (Mao & Paczynski 1991). At least $`10\%`$ of all microlensing events clearly show the typical signatures of binary lenses, including the presence of sharp caustic crossing peaks.

A conspicuous fraction of binary systems is in a state of tight interaction, where several physical phenomena show up and accompany the basic gravitational rotation. Tidal interactions, mass exchanges, the appearance of accretion disks and circumstellar matter around compact objects in their various shapes widely enrich the phenomenology of these systems. For a detailed review on strongly interacting binary systems (hereafter SIBS), we refer the reader to (Carroll & Ostlie 1996), where a full chapter is devoted to the subject.

In the present paper, we want to point out at all the possible signatures of SIBS in microlensing observations, with particular reference to chromatic absorption by circumstellar matter in microlensing light curves. These effects can be appreciated in careful observations and could provide an alternative tool to detect SIBS and determine their characteristics in a qualitative and quantitative way.

In the binary lens equation, two scales are present: the total Einstein radius $`R_\mathrm{E}`$, proportional to the square root of the total mass, and the separation between the binaries. In practice, the system will best show up its binary nature through microlensing when these two scales are of the same order. Of course, not all close binary systems will reveal their nature in a microlensing event. The separation between the two stars cannot be lower than $`0.1R_\mathrm{E}`$ (Gaudi & Gould 1997). If the two stars are too close, they will effectively behave as a single lens and their binary nature would be hidden to microlensing. Typical Einstein radii in galactic observations are of the order of AU, which are compatible with most classes of interacting binary systems, characterized by separations of the order of tenths of AU. However, if we also take into account astrometric observations of the center of light deformations during microlensing, the efficiency in the detection of close binaries may be decisively increased, so that even binaries with separations of the order $`0.01R_\mathrm{E}`$ can be detected (Chang & Han 1999).

The rotation periods of interacting binaries may range from few hours to several months. As regards those binaries which can yield distinctive features of binary microlensing (namely those which have separations larger than 0.1, 0.2 AU), for solar mass stars, the typical periods are of the order of some tens of days. This means that the rotation of the lens should be often included in a detailed study of the microlensing event. A comprehensive study of this problem has been already performed by Dominik (1998) and will not be discussed here.

## 2 Diffuse matter in binary systems and microlensing

If the size of one of the stars in the binary system exceeds its Roche lobe, a mass transfer process towards the other star initiates. If this happens, the star losing mass through its Roche lobe will be referred to as secondary star, while the star accreting mass will be the primary star. The matter leaving the secondary star spirals around the primary forming an accretion disk, whose size depends on the geometrical and dynamical properties of the system. The comprehension of the accretion disk around a compact object is important to obtain precious information about the nature of the collapsed object and its physical properties.

The matter in the accretion disk is hydrogen supplied by the external layers of the secondary star’s atmosphere. In typical situations, the temperature of the accretion disk allows the presence of neutral hydrogen. In degenerate primaries (neutron stars or black holes), however, the temperature of the inner regions of the accretion disk may rise to $`10^6`$ K, favouring a full ionization of matter.

The secondary star may also spontaneously expel mass as a natural stage in its evolutionary life or as a consequence of pulsations induced by tidal interactions. The mass lost by the the secondary star may set between the two stars, forming a circumstellar envelope of neutral hydrogen and dust (Dumm et al. 1999). If the ejection velocity is high enough, then a common envelope surrounds both stars for long periods of the life of the system.

From this discussion, we learn that the presence of gas and dust in binary systems is a typical product of strong gravitational interactions. Generically, in microlensing observations, the presence of diffuse matter shows up through chromatic absorption effects.

### 2.1 Chromatic Absorption

The light passing through the gas and the dust within the binary system is attenuated by scattering and absorption phenomena. The Lambert’s law

$$I_\lambda =I_{\lambda ,0}e^{\tau _\lambda },$$
(1)
describes the suppression experienced by light travelling through gas. $`I_{\lambda ,0}`$ is the original intensity of light, $`I_\lambda `$ is the transmitted light and $`\tau _\lambda `$ is the optical depth of the cloud, given by

$$\tau _\lambda =\underset{n}{}k_{\lambda ,n}\sigma _n(𝐱),$$
(2)
where $`k_{\lambda ,n}`$ is the absorption coefficient and $`\sigma _n(𝐱)`$ is the column density of the gas component labeled by $`n`$ at projected position $`𝐱`$.

Summing up, each image will have a total amplification given by the geometrical magnification due to gravitational light bending, weighted by the transmission coefficient $`e^{\tau _\lambda }`$.

The absorption coefficient can be expressed as the cross section over the mass of the elementary gas/dust particle. For instance, H2 Rayleigh scattering cross section is $`2.24\times 10^{27}`$ cm2 at $`\lambda =4400`$ Å. To appreciate a $`50\%`$ absorption, the column density should be about $`10^4`$ Kg/m2, which is the column density of one thousandth solar masses distributed on a uniform disk with radius equal to 2.7 AU. However, slight changes of the column density or the chemical nature of the gas may drive the absorption towards complete occultation or complete transparency. What is interesting for us is that the typical densities of circumstellar matter in interacting binary systems are compatible with partial or total absorption of background light and thus produce observable effects.

In binary microlensing, three images are normally present, while two more are created when the source is inside a caustic. Only the images which are formed inside the diffuse matter envelope are affected by absorption, while the others remain undisturbed. Of course, a complete classification of all possible light curve is beyond the purpose of this letter, since accretion disks and gaseous envelopes may appear in an extraordinary variety of shapes and can be oriented in whatever way with respect to the observer. However, in Fig. 1, we analyze some simple situations to give an idea of the anomalies that could appear in binary microlensing light curves, due to the presence of diffuse matter.

We consider a binary system composed by two half-solar mass stars, separated by 4.3 AU, at a distance of 8 kpc from the observer. With these values, rotational effects are negligible and will not be considered, so that we are allowed to concentrate just on chromatic absorption effects. The lensed source is 10 kpc distant and has a radius of 15.5 R. One of the lenses (the right one in the figures) is surrounded by a gaseous disk with projected density profile

$$\sigma (r)=1.27M_\mathrm{d}\frac{\left(r_\mathrm{d}^2r^2\right)^3}{r_\mathrm{d}^8}$$
(3)
where $`r_\mathrm{d}=2.87`$ AU is the radius of the disk, $`r`$ is the distance from the center of the right star, $`M_\mathrm{d}=0.01\mathrm{M}_{}`$ is the total mass of the disk. The choice of this density profile has been done just for computational reasons. In fact, the density profiles of circumstellar matter within interacting binaries are known only in some specific cases. This prevents us from operating more physical choices. However, since we are only interested in a qualitative description of the effects on microlensing light curves, we are allowed to pioneer the problem by Eq. (3), also neglecting tidal distortions, in a first extent. The gas in the disk will be taken to be a homogeneous mixture of H2-He with $`24\%`$ He by mass, so that the absorption coefficient is $`k_\lambda =5.16\times 10^5`$ m2/Kg at $`\lambda =4400`$ Å.

The distinctive characteristic of absorption by diffuse matter is given by its dependence on the wavelength of the light. Rayleigh scattering decreases as $`\lambda ^4`$, so that all the anomalies would be harder in blue and softer in red. This is very important for an unambiguous identification of the effect and we shall plot light curves both in blue ($`\lambda =0.44\mu \mathrm{m}`$) and in red band ($`\lambda =0.7\mu \mathrm{m}`$).

Gravitational and refractive bending induced by gaseous matter are generally negligible with respect to absorption for our parameters (Bozza et al. 2002) and will not be considered. We have chosen three examples of light curves showing the main anomalies that can be induced by the presence of diffuse matter in the system.

* Fig. 1 Col. (a) shows a microlensing light curve where the source crosses a caustic near a cusp. The right secondary image is destroyed and regenerated in two points which are completely inside the absorbing disk. This causes the two caustic crossing peaks to be highly suppressed in blue band, while very slight modifications appear in red. The secondary non-caustic peaks stay almost unmodified in both colors.
* In Fig. 1 Col. (b) the left secondary image is regenerated outside the disk and destroyed inside. As a result, only the second caustic crossing peak is depressed. Moreover, the right secondary peak, caused by the close approach between the source and the right top cusp, is suppressed as well. This is because the right secondary image, which is highly magnified by gravitational lensing during this peak, is heavily absorbed by the disk. These depressions almost disappear in red band. It is also possible to draw curves where no caustic-crossing peaks are suppressed but one secondary peak is.
* The scenario changes if also the principal image enters the disk. To show this case, in Figs 1 Col. (c) we have enlarged the disk of diffuse matter up to 7 AU, recentering it on the center of mass (again we are neglecting any ellipsoidal distortion). The curve is gradually suppressed in every region from the outer to the inner portions, as the principal image reaches the denser inner regions. The red band light curve still resembles an ordinary binary microlensing curve, but the blue band one is heavily suppressed.

By these three examples we do not mean to give an exhaustive picture, but at least we can learn that if the diffuse matter remains in the gravitational influence of one mass only, then we can expect slight distortions on the original microlensing light curve, which can be restricted to a limited number of peaks. If the diffuse matter involves the whole system, then a general suppression of the microlensing event in all its parts should be expected.

If the deviations from standard binary microlensing light curves are small, it should not be difficult to identify these events as binary microlensing and then examine a posteriori all possible anomalies comparing the light curves in different bands. If the modifications are too strong and, in particular, depress the event towards the background noise, then it may become very difficult to identify the event as microlensing. However, (as in Figs. 1 Col. (c)) it may happen that events strongly absorbed in blue are still recognizable as binary microlensing in red, so that they can be recovered anyway. Of course, to perform this kind of searches, it is necessary to relax the selection procedures based on microlensing achromaticity.

The chromaticity effects of SIBS are the key to identify them. Chromatic effects in caustic crossing peaks also arise when limb-darkening of finite size sources is taken into account (Han & Park 2001). However, in this case, the color shifts towards blue, while diffuse matter in binary systems provokes a reddening of the source light. The two effects are thus clearly distinguishable.

In exceptional cases, when one of the stars in the binary system is a low density red giant, occultations of one secondary image may occur (Bromley 1996; Agol 2002). The effects on the light curve are similar to those shown in Fig. 1 Cols. (a)-(b) for diffuse matter. The only difference is in the density profile that should be very close to a step function.

## 3 Discussion

The presence of diffuse matter within SIBS shows up through absorption of one or more microlensing images. If the diffuse matter is localized around one star, we expect a partial depression of caustic-crossing or secondary peaks, depending on the geometry of the event. If the matter englobes the whole system, then all microlensing images are simultaneously affected and a general suppression of the curve occurs, which may significantly alter the original microlensing curve. Since the modifications are more evident in blue than in red, a careful analysis in different color bands is mandatory to unambiguously detect and study these events in a serious way. Moreover, if it is possible to take a spectrum during the microlensing event, the appearance of characteristic absorption lines will be the clear signature of the presence of the diffuse medium. In this case, we would also have a way to analyze the chemical nature of this diffuse matter.

Microlensing may provide a very powerful tool for a deep physical investigation of SIBS. It may help to estimate the distribution of diffuse matter in these systems, also retrieving information on their density profiles. Eventual turbulent density fluctuations would be apparent as random noise on the light curve if their length scales are comparable to the extension of the absorbed images. In this case, a statistical study of this noise would provide information on these fluctuations. Moreover, microlensing comes with its standard information on the mass ratio and the separation of the stars along with other parameters which can be useful to give a complete characterization of the system.

An important observation is that interacting binaries are detected by traditional optical observations only when their orbital plane is edge-on with respect to the observer (unless they are in the solar neighbourhood, so that the rotation can be detected by astrometric measurements, or the two components are separated interferometrically). The binary nature of the object is then revealed by the mutual eclipses produced by rotation, or by spectroscopical measurements of radial velocities. If the orbital plane is disposed face-on with respect to the observer, it is necessary to resort to indirect proofs (e.g. X-ray emission from accretion disk, cataclysmic variables). In microlensing surveys, instead, both situations are practically detectable with the same chances. It is thus possible to enlarge the sample of known binary systems and study them from a different perspective.

###### Acknowledgements.

We wish to thank Philippe Jetzer and Gaetano Scarpetta for helpful comments.


text:
# Large Scattering Events and the Formation of Planetary Systems

## 1. The model

During the lifetime of the nebular gas, a protoplanetary disk with the profile of the standard Hayashi (1981) model, and a surface density several times the minimum, will tend to produce bodies large enough to serve as the solid cores of giant planets, i.e. around 10 M in mass, in an annulus roughly corresponding to the the Jupiter-Saturn region (Thommes, Duncan and Levison, hereafter TDL, 2002b). The inner bound results from the increase of the isolation mass (e.g., Lissauer 1987) with stellocentric distance; inside some radius, the final protoplanet masses are too small. The outer bound comes about from the increase of accretion timescale with stellocentric distance; beyond some radius, protoplanets cannot grow large enough during the lifetime of the nebular gas ($``$ 10 Myrs, e.g., Strom, Edwards and Skrutskie 1993).

In a system which forms gas giant planets by nucleated instability (e.g., Pollack et al. 1996), one of the giant protoplanets will eventually acquire a massive gas envelope. Using numerical simulations, TDL (1999, 2002a) demonstrated that this sudden mass increase of one of the protoplanets (by an order of magnitude or more over $`10^5`$ years) tends to destabilize the orbits of its neighbours, scattering them outward onto eccentric orbits. With most of its orbit now embedded in the still accretionally unevolved outer planetesimal disk, a scattered protoplanet then experiences dynamical friction (e.g., Stewart and Wetherill 1988), which damps its eccentricity down again. The end result is most often a protoplanet on a nearly circular orbit, albeit with a semimajor axis significantly larger than it originally had. In this way, a planetary system resembling the configuration of the giant planets of the Solar System is commonly produced, with the scattered and recircularized protoplanets standing in for Neptune, Uranus, and Saturn; the latter needs to acquire its own (less) massive gas envelope before the dispersal of the nebular gas, in order to match the Solar System.

In some of the simulations performed in the above work, protoplanets were so strongly scattered after “Jupiter”’s gas accretion phase that that they ended up with apocenter distances of hundreds of AU and ultimately became unbound from the Sun. However, in all these runs the planetesimal disk was only modeled out to a radius of 60 AU. Here, we investigate the possibility that, given a larger (and thus likely more realistic) planetesimal disk, such bodies might be retained on stable orbits. This would result in planetary systems with giant planets at stellocentric distances far beyond where significant accretion ought to have taken place.

## 2. Numerical simulations

The simulations performed here are very similar to those described in TDL (1999, 2002a). Four bodies, each of mass 10 M—thus large enough to serve as giant planet cores—are initially placed on stable, circular orbits between about 5 and 10 AU. Beyond this is a disk of small planetesimals, each of mass 0.1 M. Though such (Mars-mass) planetesimals are unrealistically large, the eccentricity damping rate due to dynamical friction depends principally on the overall density of the disk, rather than on the mass of the individual bodies which make it up, as long as the protoplanet eccentrities are high (details are given in TDL2002a). The main difference with respect to the earlier simulations is that the planetesimal disk extends not to 60 AU, but to 200 AU. Once the simulation commences, the mass of the innermost protoplanet is increased to 310 M, approximately the mass of Jupiter, over $`10^5`$ years, simulating the accretion of a massive gas envelope to form the first gas giant.

Fig. 1 shows the end states, after 10 Myrs, of a set of ten simulations. Initial conditions differ only in the orbital phases of the four protoplanets, which are randomly generated. The planetesimal disk is initially given the surface density of the minimum-mass protosolar disk, $`30(r/1\mathrm{AU})^{3/2}\mathrm{g}/\mathrm{cm}^2`$ (Hayashi 1981). In four of the cases, the end result is a system qualitatively similar to the outer Solar System: three protoplanets (gasless proto-Saturn, Uranus and Neptune) are on well-spaced, low-eccentricity orbits between about 10 and 30 AU. In three other cases, one of the protoplanets has been lost, either by becoming unbound, or by merging with another body. However, there are two cases, Runs C and G, where one protoplanet has been scattered much further, ending up with a semimajor axis greater than 200 AU, and between 100 and 200 AU, respectively. In both instances, the protoplanet retains a high eccentricity (0.7 - 0.8) even after 10 Myrs. This is because the dynamical friction timescale is inversely proportional to the density of the planetesimal disk, which falls off with stellocentric distance. For this same reason, simulations with more massive planetesimal disk show a lower incidence of such large scattering events; the higher density of planetesimals makes for stronger dynamical friction, and thus a stronger tendency for planetary systems to be relatively compact, like our own.

## 3. Discussion

A number of spatially resolved debris disks show features which may result from the influence of unseen (giant) planetary companions. The problem is, these features are seen at large stellocentric distances, far beyond where giant planet formation is expected to have occurred. Examples are the dust disks around HD 141569 (Weinberger et al. 1999), HR 4796A (Schneider et al. 1999), and HD 163296 (Grady et al. 2000). These systems possess, respectively, a 40 AU wide gap centered at a radius of 250 AU, a dust ring of radius 70 AU, and a 50 AU wide gap at a radius of 325 AU. A seemingly less extreme case is Vega, which displays two dust emission peaks 80 AU from the star. These have been successfully modeled as concentrations resulting from the trapping of dust particles (which spiral inward due to radiation pressure and Poynting-Robertson drag) in exterior mean-motion resonances with an eccentric gas giant (Wilner et al. 2002). A model with a 3 MJupiter planet having a semimajor axis of 40 AU and an eccentricity of 0.6 has been demonstrated to give a good match to the observations.

Large scattering events like the ones described here constitute a mechanism by which giant planets might, after forming at orbital radii similar to those of Jupiter and Saturn, be transported out to distances comparable to those of the disk features described above. This offers one way to avoid the difficult problem of explaining the formation of giant planets at 100 AU or more, just as the original model avoids the formation timescale problem of Uranus and Neptune in our Solar System. However, extending the model in this way introduces two principal problems. First, the mechanism works for “ice giants”—Uranus/Neptune mass bodies—not for bodies the mass of gas giants. In the latter case, one would be dealing with a scenario more akin to that described by Rasio and Ford (1996), in which multiple gas giants scatter each other. Eccentricities of any strongly-scattered planets would likely remain high indefinitely, since at $`>`$ 100 AU the mass in planetesimals is very small compared to that of a gas giant.

The other problem is that, as stated above, the probability of a very large scattering event decreases with the mass of the exterior planetesimal disk. A minimum-mass disk allows scattering of ice giants to semimajor axes of order 100 AU to occur readily, but a minimum-mass disk is also unlikely to produce giant planets in the first place. This will be a problem if giant planets at large stellocentric distances turn out to be ubiquitous. The above-described features do seem to be common, given that they appear in a significant fraction of the (small number of) disks that have been directly imaged thus far. In order to simultaneously form such large bodies and allow them to be readily scattered to large distances, one must invoke a steep density profile, or a local density enhancement where the giant planets form.

The above considerations provide ways in which the scenario proposed here can, in the near future, be observationally tested. First, newly-formed gas giants at large stellocentric are hot enough to be detectable, in principal, by current searches (e.g., Macintosh et al. 2000). Positive detections would immediately rule out the above mechanism as the sole source of planets at large orbital radii. Conversely, the continued absence of such detections would suggest that whatever is producing the features we see is smaller than a gas giant. Or, indeed, the features may have nothing to do with planets at all. Secondly, future observations will give us a better idea of the typical masses and density profiles of dust disks and, by extension, of the characteristics of the underlying planetesimal disks. It will then be easier to assess the likelihood of large scattering events. In particular, it will be interesting to see what correlations, if any, exist between disk mass and the occurence rate of rings, gaps and asymmetries at large radii.

## References

Grady et al. 2000. ApJ 544, 895.

Hayashi 1981. Prog. Theor. Phys. 70, 35.

Lissauer, J. J. 1987. Icarus 69, 249.

Macintosh et al. 2000. AAS Meeting 197, 64.01.

Pollack et al. 1996. Icarus 124, 62.

Schneider et al. 1999. ApJ 513, L127.

Stewart, G. R. & Wetherill, G. W. 1988. Icarus 74, 542.

Strom, S. E., Edwards, S. & Edwards, M. F. 1993. Protostars & Planets III, 837.

Thommes, E. W., Duncan, M. J. & Levison, H. F. 1999. Nature 402, 6762.

Thommes, E. W., Duncan, M. J. & Levison, H. F. 2002a. AJ 123,2862

Thommes, E. W., Duncan, M. J. & Levison, H. F. 2002b. Submitted to Icarus.

Weinberger et al. 1999. ApJ 525, L53.

Wilner et al. 2002. ApJ 569, L115.


text:
# Grand-Canonical simulation of 4D simplicial quantum gravity presented by S.Horata

## 1 Introduction

Until now 4D SQG has been investigated from several points of view, such as the phase structure, nature of the transition, effects of additional matter fields, and the modified measure action . Unfortunately the numerical results could not resolve the problem for the existence of the continuum theory of 4D QG clearly. Gradually the original enthusiasm to explore 4D QG has been toned down. However, we remind that many problems about 4D SQG are yet to be studied.

As one of those problems, we pay attention to additional matter effects in 4D Euclidean SQG (Eucl.SQG). According to the recent numerical simulations, the string susceptibility exponent ($`\gamma ^{(4)}`$) takes a negative value with adding a few $`\mathrm{U}(1)`$ gauge fields. It makes drastic change of the phase diagram, and continuous phase transition with one gauge matter has been found. Moreover, the numerically obtained $`\gamma ^{(4)}`$ as a function of the number of matter fields turned out to be approximately equal to that of the prediction of 4D conformal field theory.

In order to see more definitely how the matter fields act on 4D Eucl.SQG, we repeat the study with employing the Grand-Canonical method used for 2D SQG, which has given the satisfactory agreement of two results. For discussing the continuum theory of 4D QG based on the numerical analysis, we pay attention to the functional form of the partition function, which we expect enough to gain information for selecting the correct theory of gravity, if it exists at all. For obtaining the functional form of the partition function, we need to perform a simulation with varying the space volume. In practice, we calculate the probability ($`P(N_4)`$) as a function of the volume ($`N_4`$) by the Grand-Canonical method of the dynamical triangulation (DT) and estimate $`\gamma ^{(4)}`$ from the scaling behavior of $`P(N_4)`$.

## 2 Model and Grand-Canonical method

Let us explain our model and the numerical method of the Grand-Canonical simulation. We consider a system of the simplicial gravity coupled to $`\mathrm{U}(1)`$ gauge fields and massless scalar fields. The action ($`S`$) is given on 4D simplicial manifold as

$$S=S_G+S_A+S_X,$$
(1)
where $`S_G`$, $`S_A`$ and $`S_X`$ denote the action for the gravity, $`\mathrm{U}(1)`$ gauge fields ($`A`$) and scalar fields ($`X`$), respectively.

For the gravity part, we use the discretized Einstein-Hilbert action in 4D,

$$S_G[\kappa _2,\kappa _4]=\kappa _4N_4\kappa _2N_2,$$
(2)
where $`N_i`$ denotes the number of $`i`$-simplex. The two parameters, $`\kappa _2`$ and $`\kappa _4`$, correspond to the inverse of the gravitational constant and the cosmological constant, respectively. The action for the $`N_X`$ scalar fields $`X_i^a`$, $`a=1\mathrm{}N_X`$, on the vertex $`i`$ is given by

$$S_X=\underset{a=1}{\overset{N_X}{}}\underset{ij}{}o(l_{ij})(X_i^aX_j^a)^2,$$
(3)
where $`o(l_{ij})`$ is the number of four-simplices sharing a link $`l_{ij}`$. And the action for the $`N_A`$ vector field $`A_{ij}^a`$( $`a=1\mathrm{}N_A`$) on link $`l_{ij}`$ reads

$$S_A=\underset{a=1}{\overset{N_A}{}}\underset{ijk}{}o(t_{ijk})(A^a(l_{ij})+A^a(l_{jk})+A^a(l_{ki}))^2,$$
(4)
where $`o(t_{ijk})`$ is the number of four-simplices sharing a triangle $`t_{ijk}`$. These actions give the correct thermo-dynamical limit for large $`N_4`$ because they are proportional to $`N_4`$.

In order to simplify the calculation, the Monte-Carlo simulations have been performed with the Canonical method previously. This method requires that $`N_4`$ fluctuates around the given target volume ($`V_4`$). Thus, the volume fixing term ($`\delta S`$) is added to the total action Eq.(1),

$$\delta S(N_4,V_4)=\delta \kappa _4\left(N_4V_4\right)^2,$$
(5)
where $`\delta \kappa _4`$ is the free parameter. However, the Canonical method has limitations about the finite size effect and the ergodic property. Furthermore, the functional form of the size dependence of the partition function cannot be given clearly. For the improvement of numerical accuracy we perform the simulations varying the space volume by the Grand-Canonical method. The Grand-Canonical partition function with the matter fields is given as

$`Z(\kappa _2,\kappa _4,N_A,N_X)={\displaystyle \underset{T}{}}e^{S_G(\kappa _2,\kappa _4)}`$
$`{\displaystyle \underset{N_A,N_X}{}}\left({\displaystyle \underset{lT}{}dA_le^{S_A(A_l)}\underset{iT}{}dX_ie^{S_X(X_i)}}\right),`$ (6)
where all geometries $`T`$ and matter fields $`A,X`$ are updated using the Metropolis method keeping the topology of discretized surface to be $`S^4`$.

From the analogy of the 2D case , we assume that the partition function obeys the KPZ like formula for the large-$`N_4`$ limit as

$$Z\underset{N_4}{}N_4^{\gamma ^{(4)}(\kappa _2)3}e^{\left(\kappa _4\kappa _4^c\left(\kappa _2\right)\right)N_4}.$$
(7)
Then the probability $`P(N_4)`$ for the large-$`N_4`$ behaves approximately as

$$P(N_4)N_4^{\gamma ^{(4)}3}e^{(\kappa _4\kappa _4^c)N_4}.$$
(8)
By tuning $`\kappa _4`$ close to $`\kappa _4^c`$, the exponential growth is suppressed and $`P(N_4)`$ is given as the power function of $`N_4`$. The parameter $`\kappa _4`$ is determined in the following manner: If $`\kappa _4`$ is chosen to be less (or greater) than $`\kappa _4^c`$, the system will expand (or shrink) exponentially in the Grand-Canonical simulation. By requiring an exponential increase (or decrease) to disappear in the simulation, we can find $`\kappa _4`$ precisely with the following form,

$$\frac{\mathrm{ln}P(N_4)}{\mathrm{ln}N_4}=(\gamma ^{(4)}3)(\kappa _4\kappa _4^c)\frac{N_4}{\mathrm{ln}N_4}+O(1/\mathrm{ln}N_4).$$
(9)
After tuning $`\kappa _4`$, we generate the configurations constraining $`N_4`$ within an upper- and a lower-bound ($`N_4^{min}N_4N_4^{max}`$).

## 3 Numerical result

In this section, we report numerical results with the Grand-Canonical method.

In order to obtain $`P(N_4)`$, we allow $`N_4`$ to fluctuate between $`N_4^{min}=1000`$ and $`N_4^{max}=10000`$. In Fig.1, we show $`P(N_4)`$ versus $`N_4`$ at $`\kappa _4`$ close to $`\kappa _4^c`$. The numerical result shows that $`P(N_4)`$ is given as a power function of $`N_4`$. It is suggested that the partition function in 4D obeys the KPZ type partition function as the 2D case.

We can also compute $`\gamma ^{(4)}`$ from the slope of $`P(N_4)`$. Then, we show the parameter ($`b`$), which is calculated from $`\gamma ^4`$,

$$\gamma ^{(4)}=\frac{b}{2}\left(1\frac{4}{b}+\sqrt{1\frac{4}{b}}\right).$$
(10)
This relation is given from the analogy of 2D QG and the parameter $`b`$ is similar to the central charge in 2D QG. In Fig.2, we plot $`b`$ versus the number of matter fields, $`N_X+62N_A`$, on the critical point ($`\kappa _2^c`$) between the crumpled phase and the smooth phase. We compare $`b`$ obtained by two different numerical methods, i.e. the Canonical method with the MINBU algorithm and the Grand-Canonical method. Error bars include both the errors in determining $`\kappa _2^c`$ and the statistical errors, and we found that the former method has much larger errors than the latter method. We made the parameter fit of the data of the Grand-Canonical method by a linear function, $`b=0.0030(3)\left(N_X+62N_A\right)+3.98(3)`$. The slope $`b`$ is approximately equal to the analytical result.

## 4 Summary and Discussion

Let us summarize and discuss the main points made in the previous section. From the numerical results of the Grand-Canonical method, we find the KPZ type partition function also appears in 4D QG. The string susceptibility exponent $`\gamma ^{(4)}`$ obtained from the size dependence of the partition function can be obtained by the Grand-Canonical method rather accurately. The matter contribution to $`b`$, which is computed from $`\gamma ^{(4)}`$ with the analogy of 2D QG, turns out to be approximately equal to the conformal part. Our numerical results suggest that it may be possible to understand the continuum theory of 4D QG from the analogy of 2D QG. In order to look more details on the relation to the continuum theory of 4D QG, these results encourage us to continue investigation by the Grand-Canonical method and push more efforts on the analytical calculation of the conformal gravity at the same time.


text:
# Inventory of Black Hole Binaries

## 1. Definition and Source Selection

Black holes represent an extreme manifestation of Einstein’s theory of general relativity. As an observational astronomer, I will not consider any of the detailed theory of black holes but will instead consider only the practical question of how to find them. The usual route is by the process of elimination. In the current universe, black holes must form via the gravitational collapse of a dying star. Three outcomes are possible, depending on the mass of the degenerate core: a white dwarf, a neutron star, or a black hole. The mass of a white dwarf cannot exceed the well-known Chandrasekhar limit, and a neutron star has a somewhat analogous upper mass limit, generally thought to be on the order of $`3M_{}`$ (Rhoades & Ruffini 1974; Kalogera & Baym 1996). Once the mass of the degenerate core exceeds about $`3M_{}`$, no known force can halt the gravitational collapse, and a black hole must be formed. Given this, we have a relatively straightforward observational definition of a black hole: A black hole is a compact object with a mass greater than three solar masses, where “compact” in this context means an object with a radius smaller than about 100 km (i.e. much less than the radius of normal stars).

Since black holes are dark, the only way one could hope to observe them is through their gravitational influence on surrounding matter. Early attempts to search catalogs of spectroscopic binaries to look for single-lined binaries with massive and undetected companions were not successful (Zel’dovich & Guseynov 1966; Trimble & Thorne 1969). A far more efficient approach for source selection has its roots in the 1960s when Zel’dovich and others realized that if a black hole accreted material (either from a nearby companion star or the interstellar medium), then it might shine brightly in X-rays and $`\gamma `$ rays. Today we know of a few hundred bright X-ray sources which must be powered by accretion onto a compact object, either a neutron star or a black hole (i.e. rapid variability in many cases indicates a size scale of a few hundred km or less). All of the black hole candidates discussed below were selected on the basis of their X-ray activity, although two of the sources, V404 Cyg and V4641 Sgr, were previously known (optically) variable stars that were essentially ignored until they were associated with bright X-ray sources. Thus, one should keep in mind that owing to special circumstances needed to observe them (i.e. interacting binaries), the list of presently known stellar-mass black holes is probably a very biased sample and may not be representative of the general population of black holes.

## 2. Mass Limits

By my definition above, to establish the presence of a black hole in a given system one must show that there is a compact object with a mass greater than $`3M_{}`$. Establishing a lower limit to the mass is relatively easy since in a single-lined spectroscopic binary, the observed radial velocity curve of the visible star can be used to set a firm lower limit on the mass of the unseen companion:

$$f(M)\frac{PK_2^3}{2\pi G}=\frac{M_1^3\mathrm{sin}^3i}{(M_1+M_2)^2},$$
(1)
where $`P`$ is the binary orbital period, $`K_2`$ is the semiamplitude of the companion’s radial velocity curve, $`G`$ is the gravitation constant, $`i`$ is the inclination of the orbit to the line of sight ($`i=90^{}`$ for an orbit seen edge-on), and where $`M_1`$ and $`M_2`$ are the masses of the invisible object and the companion star, respectively. One can easily see that $`f(M)`$ (generally called the “optical mass function” or “mass function”) is always less than $`M_1`$ since $`M_2>0`$ and $`\mathrm{sin}i1`$ always. Thus a firm lower limit on the mass of the compact object can be obtained from two relatively simple measurements.

In real life, of course, things are always a bit more complicated. For starters, in many cases one cannot even see spectral lines from the companion star in many of the persistently X-ray bright sources. Needless to say this precludes any attempts at dynamical mass estimates. In other cases, the companion star may not be uniformly bright over its surface owing to heating by the X-ray source. If so, then the “centre of light” will not be the same as the “centre of mass” and the observed radial velocity curve will be different than the true one. Since the value of the mass function depends on $`K_2^3`$, small changes in $`K_2`$ can lead to relatively large changes in $`f(M)`$. Fortunately, many of the potential difficulties related to X-ray heating are largely absent in the subset of the X-ray binaries known as “X-ray transients”. As the name suggests, the X-ray luminosity in these sources is highly variable (dynamic ranges of $`10^6`$ to $`10^7`$ are not uncommon). Dynamical observations of these systems are generally reliable, provided one observes the system after the X-ray activity has completely ceased (see Phillips, Shahbaz, & Podsiadlowski 1999; Shahbaz et al. 1999a; and Shahbaz et al. 2000 for discussions regarding X-ray heating in GRO J1655-40).

Table 1 lists the mass functions for 17 systems that are widely believed to contain black holes. Fourteen of these systems are transient X-ray sources with relatively low-mass companions ($`M_2<1M_{}`$ generally), and the remaining three systems have massive O/B companions whose bolometric luminosity is comparable to or exceeds the X-ray luminosity of the accreting object. Six of the systems have mass functions which are above $`5M_{}`$, and six more have mass functions more than $`3M_{}`$. Based on the mass functions alone, we can safely conclude these systems contain black holes since the compact objects in them are too massive to be neutron stars. Although the other five systems have relatively small mass functions, we can be reasonably sure these systems also contain black holes since the currently best estimates for the compact object masses (see the next section) are all above $`3.5M_{}`$. In addition to the sources listed in Table 1, there are a few dozen sources which are believed to contain a black hole on the basis of their X-ray spectra and temporal properties (e.g. Tanaka & Lewin 1995; Liu, van Paradijs, & van den Heuvel, 2001).

## 3. Towards Actual Masses

If possible, we want to do better than simply establishing a lower limit to the mass of a black hole. If we think of the companion star as a test particle orbiting the black hole, then clearly in order to measure the actual mass of the black hole after its mass function has been measured, one must be able to convert projected radial velocities into actual radial velocities (i.e. the inclination $`i`$ is needed) and to find out how far the star is from the binary centre-of-mass (i.e. one needs to know the mass ratio or the mass of the companion star).

The procedure to derive astrophysically interesting parameters of a binary system is conceptually simple. One can make a model of a close binary system that is specified by several parameters such as the orbital period, the separation between the stars, the mass ratio, the relative radii of the stars, the temperatures of the stars, the inclination to the line of sight, etc. Then, given this model, one can compute various observable properties of this binary such as the radial velocity curves, the light curves in various bandpasses, eclipse duration(s) if any, the rotational velocities of the stars, etc. The goal then is to find a set of input model parameters that can best match the set of observed properties of the binary of interest. Given the optimal input parameters, one can then compute astrophysically interesting quantities such as the masses of the components.

Models of close binary stars are relatively easy to construct, and have been use for over thirty years (e.g. Wilson & Devinney 1971; Lyutyi, Syunyaev, & Cherepashchuk 1973; Hutchings 1974; Avni & Bahcall 1975). The difficult part has been related to the inverse problem of parameter optimisation since the parameter space to search is usually vast and often parameters can be tightly correlated. Common ways to tackle the inverse problem include “differential corrections” and related techniques (e.g. Wilson & Devinney 1971) and systematic “grid searches” (e.g. Kuiper, van Paradijs, & van der Klis 1988). More recently, I have been using optimisation techniques based on genetic algorithms (Orosz et al. 2002a). Genetic algorithms are very robust, and can allow for accurate estimation of uncertainties on the fitted and derived parameters.

In the case of the black hole binaries, the main observables usually are a radial velocity curve for the companion star (this basically sets the scale of the binary), a light curve for the companion in one or more bandpasses (these mainly constrain the inclination), and in most cases a measurement of the mean projected rotational velocity of the companion star (this sets the mass ratio). In nearly cases (with the possible exception of SAX J1819.3-2525), the X-ray source is not eclipsed, and one can set an upper limit to the inclination. Given these observations, it is not an unduly difficult problem (in principle!) to find a good binary model for most of the black hole systems listed in Table 1.

However, as before, in real life things are more complicated. In many cases, the observed light curve is not entirely due to the companion star owing to the presence of extra sources of light (usually from an accretion disk). In the optical, the mean light curve can change from one year to the next (e.g. Orosz et al. 1996; Leibowitz, Hemar, & Orio 1998; Pavlenko et al. 2001). Presumably these long-term changes are related to the accretion disk. From spectroscopy, one can estimate the light attributable to the accretion disk (e.g. Marsh, Robinson, & Wood 1994) and correct for it in the models. Alternatively, it is generally assumed that the contamination from the accretion disk is small in the $`J`$, $`H`$, and $`K`$ infrared bands (e.g. Shahbaz et al. 1994; Gelino, Harrison, & McNamara 2001). There have a few attempts to directly measure the disk contamination in the infrared with spectroscopy, but the results so far are not terribly constraining (Shahbaz et al. 1996; Shahbaz, Bandyopadhyay, & Charles 1999b).

I have summarised in Table 1 the best (in my opinion) measurements of the rotational velocities and inclinations for the 17 black hole binaries. Relatively uncertain values are denoted by ’?’. In some cases there is reasonably good agreement between inclination measurements made by different groups (e.g. GS 1124-683), while in other cases there have been a wide range of measurements (e.g. Cyg X-1, GRO J0422+32). I have computed $`1\sigma `$ ranges for the black hole mass and mass ratio using a simple Monte Carlo code that makes use of numerical Roche lobe integrations to map a rotational velocity into a mass ratio. In some cases, rotational velocities are not available and instead I have guessed the mass of the companion in order to arrive at a mass range for the black hole. Figure 1 shows graphically the black hole mass ranges and, for comparison, masses for neutron stars. All things considered, many of the masses are known reasonably well, but one should always keep in mind potential systematic errors such as those related to unstable light curves, etc.

## References

Avni, Y., & Bahcall, J. N. 1975, ApJ 197, 675

Brocksopp, C., Tarasov, A. E., Lyutyi, V. M., Roche, P. 1999, A&A 343, 86

Casares, J., & Charles, P. A. 1994, MNRAS 271, L5

Casares, J., Martin, E. L., Charles, P. A., Molaro, P., & Rebolo, R 1997, New Astronomy 1, 299

Chevalier, C., & Ilovaisky, S. A. 1993, A&A 269, 301

Cowley, A. P., Crampton, D., Hutchings, J. B., Remillard, R., & Penfold, J. E. 1983, ApJ 272, 118

Fender, R. P., et al. 1999, MNRAS 304, 865

Filippenko, A. V., Matheson, T., Leonard, D. C., Barth, A. J., & van Dyk, S. D. 1997, PASP 109, 461

Filippenko, A. V., Leonard, D. C., Matheson, T., Li, W., Moran, E. C., & Riess, A. G. 1999, PASP 111, 969

Filippenko, A. V., & Chornock, R. 2001, IAU Circ. #7644

Gelino, D. M., Harrison, T. E., & McNamara, B. J. 2001, AJ 122, 971

Gelino, D. M., Harrison, T. E., & Orosz, J. A. 2001, AJ 122, 2668

Gelino, D. M. 2001, Ph.D. Thesis, New Mexico State University

Gelino, D. M. 2002, BAAS 200, #08.11

Gelino, D. M., & Harrison, T. E. 2002, ApJ, in press

Gies, D. R., & Bolton, C. T. 1986, ApJ 304, 371

Greene, J., Bailyn, C. D., & Orosz, J. A. 2001, ApJ 554, 1290

Greiner, J., Cuby, J. G., & McCaughrean, M. J. 2001, Nature 414, 522

Harlaftis, E. T., Horne, K., & Filippenko, A. V. 1996, PASP 108, 762

Harlaftis, E., Steeghs, D., Horne, K., & Filippenko, A. V. 1997, AJ 114, 1170

Harlaftis, E., Collier, S., Horne, K., & Filippenko, A. V. 1999, A&A 341, 491

Herraro, A., Kudritzki, R. P., Gabler, R., Vilchez, J. M., & Gabler, A. 1995, A&A 297, 556

Hutchings, J. B. 1974, ApJ 193, L61

Hutchings, J. B., Crampton, D., Cowley, A. P., Bianchi, L., & Thompson, I. B. 1987, AJ 94, 340

Israelian, G., Rebolo, R., Basri, G., Casares, J., & Martin, E. L. 1999, Nature 401, 142

Kalogera, V., & Baym, G. 1996, ApJ 470, L61

Kuiper, L., van Paradijs, J., & van der Klis, M. 1988, A&A 203, 79

Leibowitz, E. M., Hemar, S., & Orio, M. 1998, MNRAS 300, 463

Liu, Q. Z., van Paradijs, J., & van den Heuvel, E. P. J 2001, A&A 368, 1021

Lyutyi, V. M., Syunyaev, R. A., & Cherepashchuk, A. M. 1973, Soviet Astron. 17, 1

Marsh, T. R., Robinson, E. L., & Wood, J. H. 1994, MNRAS 266, 137

McClintock, J. E., Garcia, M. R., Caldwell, M., Falco, E. E., Garnavich, P. M., & Zhao, P. 2001, ApJ 551, L147

Mirabel, I. F., & Rodriguez, L. F. 1994, Nature 371, 46

Ninkov, Z., Walker, G. A. H., & Yang, S. 1987, ApJ 321, 425

Orosz, J. A. 2001, ATel #67

Orosz, J. A., Bailyn, C. D., McClintock, J. E., & Remillard, R. A. 1996, ApJ 468, 380

Orosz, J. A., Kuulkers, E., van der Klis, M., McClintock, J. E., Garcia, M. R., Callanan, P. J., Bailyn, C. D., Jain, R. K., & Remillard, R. A. 2001, ApJ 555, 489

Orosz, J. A., Groot, P. J., van der Klis, M., McClintock, J. E., Garcia, M. R., Zhao, P., Jain, R. K., Bailyn, C. D., & Remillard, R. A. 2002a, ApJ 568, 845

Orosz, J. A., et al. 2002b, ApJ in preparation

Orosz, J. A., et al. 2002c, ApJ in preparation

Pavlenko, E. P., Kuznetsova, Y. G., Shugarov, S. Yu., & Petrov, V. S. 2001, ApSSS 276, 65

Phillips, S. N., Shahbaz, T., & Podsiadlowski, Ph. 1999, MNRAS 304, 839

Remillard, R. A., Orosz, J. A., McClintock, J. E., & Bailyn, C. D. 1996, ApJ 459, 226

Rhoades, C. E., Jr., & Ruffini, R. 1974, Phys. Rev. Lett 32, 324

Shahbaz, T., Ringwald, F. A., Bunn, J. C., Naylor, T., Charles, P. A., & Casares, J. 1994, MNRAS 271, L10

Shahbaz, T., Bandyopadhyay, R., Charles, P. A., & Naylor, T. 1996, MNRAS 282, 977

Shahbaz, T., Naylor, T., & Charles, P. A. 1997, MNRAS 285, 607

Shahbaz, T., van der Hooft, F., Casares, J., Charles, P. A., & van Paradijs, J. 1999a, MNRAS 306, 89

Shahbaz, T., Bandyopadhyay, R. M., & Charles, P. A. 1999b, A&A 346, 82

Shahbaz, T., Groot, P., Phillips, S. N., Casares, J., Charles, P. A., & van Paradijs, J. 2000, MNRAS 314, 747

Tanaka, Y., & Lewin, W. H. G. 1995, in: W. H. G. Lewin, J. van Paradijs, & E. P. J. van den Heuvel (eds.), X-ray Binaries (Cambridge: Cambridge University Press), Cambridge Astrophysical Series 26, 126

Thorsett, S. E., & Chakrabarty, D. 1999, ApJ 512, 288

Trimble, V. L., & Thorne, K. S. 1969, ApJ 156, 1013

van der Klis, M., Clausen, M., Jensen, K., Tjemkes, S., & van Paradijs, J. 1985, A&A 151, 322

Wagner, R. M., Foltz, C. B., Shahbaz, T., Casares, J., Charles, P. A., Starrfield, S. G., & Hewett, P. 2001, ApJ 556, 42

Webb, N. A., Naylor, T., Ioannou, Z., Charles, P. A., & Shahbaz, T. 2000, MNRAS 317, 528

Wilson, R. E., & Devinney, E. J. 1971, ApJ 166, 605

Zel’dovich, Ya. B., & Guseynov, O. H. 1965, ApJ 144, 840


text:
# Soliton back-action evading measurement using spectral filtering

## I Introduction

Quantum theory allows for the measurement of an observable with arbitrarily high precision. In a back-action evading (BAE) measurement the back action inherent in the process of a quantum measurement is confined to the observable which is conjugate to the measured one . If the initial state of a system is an eigenstate of the measurement operator, it is conserved and the BAE interaction is called quantum-nondemolition (QND). Initially designed for the detection of gravitational waves, the concept of QND measurements now is finding applications in the field of quantum information processing and communication .

Most BAE measurements have been carried out in quantum optics, where the quantum fluctuations are readily measurable. Quantitative criteria for QND and BAE experiments have been developed . There are three physical systems in which BAE detection has been pursued experimentally: The $`\chi ^{(3)}`$ interactions in optical fibers  and in atoms , and the $`\chi ^{(2)}`$ interaction in crystals . The QND criteria were fulfilled in BAE detection with atomic systems ($`\chi ^{(3)}`$)  and $`\chi ^{(2)}`$ interactions in crystals . Repeated BAE measurements were reported using OPAs  and a demonstration of principle using fiber optical pulses. Despite their large potential for QND measurements  fiber optical experiments have not previously been performed in the BAE regime .

In this Paper we report experimental BAE results following a recently proposed fiber-optical back-action evading measurement scheme of the soliton photon number using spectral filtering of solitons .

Optical solitons are a natural physical system for implementing quantum communication concepts in fibers. Recently, they were used for EPR-pair beam generation . Solitons are fiber optical pulses maintaining a stable shape due to a balance of group-velocity dispersion and nonlinearity. Moreover, after the collision of two solitons, the pulses recover energy, velocity and shape, and thus behave like particles. The nonstationary evolution of the soliton quantum fluctuations leads to remarkable effects. Quadrature squeezing was the first effect in a series of quantum optics experiments with solitons in fibers . In recent years amplitude squeezing of solitons using spectral filtering was discovered and explained as a multi mode quantum effect .

The technique of spectral filtering has since been generalized to a system of two colliding solitons . A new theoretical proposal for BAE measurements was derived, in which a ”signal” soliton collides with a second ”probe” soliton. During the collision, the frequency of the probe becomes coupled to the photon number of the signal. A measurement of the probe frequency therefore gives a BAE readout of the signal photon number. Using this scheme we have obtained experimental BAE operation in the quantum regime. For the first time in fiber QND detection, the standard QND criteria are clearly fulfilled.

## II Principle of the BAE measurement

The classical field evolution in a single-mode polarization preserving fiber is governed by the nonlinear Schrödinger equation (NLSE):

$$i\frac{u}{\xi }+\frac{1}{2}\frac{^2u}{\tau ^2}+|u|^2u=0,$$
(1)
where $`u`$ is the slowly varying envelope of the electric field and $`\xi `$ and $`\tau `$ are time and space variables in a reference frame moving with the optical field . $`u`$, $`\xi `$ and $`\tau `$ are given in soliton units . The single pulse soliton solutions, i.e. the fundamental soliton of amplitude $`A`$ is given by :

$`u(\tau ,\xi )=A\mathrm{sech}(A(\tau +\mathrm{\Omega }\xi ))\mathrm{exp}(i(A^2\mathrm{\Omega }^2)\xi /2i\mathrm{\Omega }\tau ).`$ (2)
$`\mathrm{\Omega }`$ is a dimensionless velocity of the soliton in units of inverse pulse length relative to the moving frame. At the same time $`\mathrm{\Omega }`$ is a dimensionless center frequency of the soliton, because of the group-velocity dispersion. The collision of two fundamental solitons of equal amplitudes $`A`$, moving before the collision with $`\pm \mathrm{\Omega }_0`$, is depicted in Figure 1.

As the solitons begin to overlap the enhanced intensity in the overlap region causes an attractive force and acceleration between the pulses. Therefore, the initial spectral separation and relative velocity of the solitons is transiently enhanced above the initial value of $`2\mathrm{\Omega }_0`$. The spectral separation of the pulses, ”signal” and ”probe”, changes to $`2\mathrm{\Omega }_0+d\mathrm{\Omega }_S+d\mathrm{\Omega }_P`$ and depends on the photon number of both solitons. Therefore, if the signal amplitude $`A_S`$ fluctuates around $`A=A_S`$, the fluctuations $`\mathrm{\Delta }A_S=A_SA`$ induce fluctuations in the frequency $`\mathrm{\Delta }\mathrm{\Omega }_P`$ of the probe soliton at the collision center. If the soliton spectra do not overlap ($`\mathrm{\Omega }_01`$), the probe shift is given by :

$$d\mathrm{\Omega }_P+\mathrm{\Delta }\mathrm{\Omega }_P=\frac{2A^2}{3\mathrm{\Omega }_0}(1+\frac{\mathrm{\Delta }A_S}{2A}+𝒪(\mathrm{\Delta }A_S^2))$$
(3)
Towards the end of the collision the solitons recover their initial velocities (Fig.1).

The main idea of the novel BAE scheme is to measure the signal soliton amplitude fluctuations via the probe soliton frequency. The coupling is described classically by Eq. 3. Signal and probe solitons can be spectrally separated at the collision center (Fig.1). Then the frequency fluctuations $`\mathrm{\Delta }\mathrm{\Omega }_P`$ of the probe are read out with a spectral edge filter. The signal, which is the photon number, is preserved in the interaction as well as in the free fiber propagation, because of negligible fiber loss. The back action of the measurement perturbs the frequency and the phase of the signal pulse. For a quantum readout, the probe readout has to determine the signal intensity to better than the signal shot noise uncertainty. This measurement scheme was assessed by a detailed theoretical investigation accounting for the quantization of the field $`u`$. The result indicates that the signal is measured more precisely than within the shot noise and the scheme fulfills the QND criteria .

## III Experimental realization

The experimental realization of the measurement requires three main steps: the pulse pair for the soliton collision has to be prepared and launched into the fiber. Then the detector is set up for the measurement of the conditional variance. Finally, the transfer coefficients are measured using an RF modulation on the signal input soliton. An overview of the experimental apparatus is presented in Fig. 2. The laser source is a $`\mathrm{Cr}^{4+}`$:YAG laser delivering 150-fs $`\mathrm{sech}`$-pulses at $`\lambda =1.5\mu `$m . In the pulse former every pulse is spectrally split using a wavelength dependent beam splitter and retroreflected, double passing the beam splitter such that two pulses of different center wavelengths and with an adjustable temporal separation are obtained. They exhibit no relative phase noise or timing jitter in contrast to pulses from different laser sources. Figure 3(a) displays the two pulses after the pulse former. They have a spectral width of 12 nm and are separated by 14 nm.

The temporal width of the autocorrelation traces indicates pulse lengths of 200 fs. The time-bandwidth product is slightly increased above the Fourier limit (sech-shape assumed). In consequence, a weak dispersive wave evolves in the fiber in addition to the two solitons. The timing of the input pulses was adjusted to the output wave autocorrelation trace, taken for every measurement, such that the fiber end is exactly in the middle of the collision ($`\xi =0`$). As the signal to be measured in the BAE apparatus we consider the amplitude of the signal pulse after it is coupled into the fiber; therefore the coupling losses do not degrade the signal transfer through the BAE apparatus. Before the fiber the input pulse pair was shown to be shot-noise limited around 20 MHz to within $`\pm 0.2`$ dB. We used 6.3 m of polarization preserving single mode fiber (3M, FS-PM 7811) corresponding to 3.2 and 2.9 soliton periods  for the signal and probe solitons, respectively. The collision length, i.e. the length of the fiber where the solitons overlap within their (power) half width, is $`L_{\mathrm{coll}}=3.2`$ m. The solitons experience a nonlinear cross-phase shift of $`\delta \varphi _{NL}1.0`$ rad. An output spectrum is shown in Fig. 3(c) and can be compared to the spectra of pulses which have propagated individually without collision (3(b)). The strong nonlinear coupling is clearly observable. The respective spectra shift about half a spectral width, $`d\lambda 5`$ nm. The undulations in the spectra are caused by the weak dispersive wave. Due to the partial spectral overlap ($`\mathrm{\Omega }_0=1.2`$) the soliton collision is phase dependent. This modifies the undulations, but the spectral shift changes by less than $`5\%`$.

The output pulses are spectrally dispersed using a grating (600 l/mm). Spectral filtering of the probe pulse is accomplished with a knife edge filter in the Fourier plane of the grating. Signal and probe pulses are separately directed onto two balanced two-port detectors (Epitaxx ETX500 photodiodes). The RF photocurrents are amplified after suppression of the repetition rate of the laser (163MHz). The sum current $`I_P`$ of the probe detector is attenuated with an attenuation $`a`$ and added to and subtracted from the signal sum current $`I_S`$. The shot-noise levels were determined by reading the DC photocurrents which display the average detected powers. They were carefully and repeatedly calibrated to the shot-noise level using the difference photocurrents of the balanced detectors which carry fluctuations of the signal and probe shot-noise levels.

The overall quantum efficiency for signal and probe detectors, including photodiode efficiencies, was measured to be 74.5% and 69.8%, respectively. The fiber loss ($`\alpha <0.002\mathrm{dB}/\mathrm{m}`$) was neglected. Coupling losses from the fiber were minimized by means of an antireflection coated gradient-index lens.

The conditional variance and the transfer coefficients quantify the performance of a BAE experiment . The conditional variance is a measure of the quantum state preparation ability of the apparatus. It can be written in terms of the photocurrents $`I`$ as:

$$V_{S/P}=\underset{a}{\mathrm{min}}\left(\frac{(I_S+aI_P)^2}{I_{S,SN}^2}\right)$$
(4)
assuming $`I=0`$. $`I_S`$ ($`I_P`$) measures the signal (probe) amplitude fluctuations, since the signal (probe) is directly detected. $`I_{S,SN}`$ denotes the corresponding shot noise.

## IV Results

A measurement of the conditional variance is shown in Fig. 4. The three noise power levels of signal, sum and difference currents are recorded simultaneously on three RF spectrum analyzers. The photocurrent fluctuations are measured at $`20`$ MHz with a $`300`$ kHz resolution bandwidth. Thermal noise in the detector electronics is $`6`$ dB below the lowest detected noise level and is subtracted off in Fig. 4. The noise powers are quadratic functions of the attenuation $`a`$ as indicated by the parabolic fitting curve of Eq. 4. The displayed noise levels are normalized to the signal output shot noise. Signal and probe pulse are separated at $`\lambda _{\mathrm{sep}}=1506`$ nm after the fiber; the frequency low pass filter for the probe was applied at $`\lambda _{\mathrm{filt}}=1490`$ nm. At this position roughly $`18\%`$ of the probe pulse energy is filtered off and absorbed. The sum noise level in Fig. 4 decreases below the signal noise. This noise suppression corresponds to a strong negative photon number correlation of $`C_{S/P}=0.62\pm 0.03`$, with the correlation coefficient defined as in .

An anti-correlation is anticipated since an increase in the photon number of the signal causes an enhancement in the mutual spectral shift, in turn causing more losses in the probe pulse at the spectral filter. This correlation is not corrected for the finite quantum efficiencies of the detection apparatus. The noise reduction reaches $`1.37\pm 0.3`$ dB below the shot noise level of the signal, corresponding to a conditional variance of $`V_{S/P}=0.73\pm 0.04`$. The necessary condition for a QND experiment is $`V<1`$. This is clearly fulfilled.

The second important criterion for a QND measurement is the sum of the transfer coefficients. They describe the nondemolition property of the measurement with the transfer of the optical signal-to-noise ratio from the input to the output:

$`\mathrm{T}_S={\displaystyle \frac{\mathrm{SNR}_S^{\mathrm{out}}}{\mathrm{SNR}_S^{\mathrm{in}}}}={\displaystyle \frac{(\mathrm{\Delta }n_S^{\mathrm{out}})^2n_S^{\mathrm{in}}}{(\mathrm{\Delta }n_S^{\mathrm{in}})^2n_S^{\mathrm{out}}}}`$ (5)
$`\mathrm{T}_P={\displaystyle \frac{\mathrm{SNR}_P^{\mathrm{out}}}{\mathrm{SNR}_S^{\mathrm{in}}}}={\displaystyle \frac{(\mathrm{\Delta }n_P^{\mathrm{out}})^2n_S^{\mathrm{in}}}{(\mathrm{\Delta }n_S^{\mathrm{in}})^2n_P^{\mathrm{out}}}},`$ (6)
where $`n`$ denotes the photon numbers and $`\mathrm{\Delta }n`$ is the modulation strength. Therefore we modulated the intensity of the signal beam with a $`\text{LiTaO}_3`$-electro-optic modulator at $`20`$ MHz before coupling into the fiber. The noise levels on and near to this frequency ($`50`$ kHz off the 20 MHz signal, recorded with a resolution bandwidth of 10 kHz.) were taken and the signal-to-noise ratio ($`\mathrm{SNR}`$) for the signal input, the signal output and the probe output was determined . If the probe soliton is absent, the unperturbed signal input exits the fiber (fiber damping and output coupling losses neglected). This way the signal input $`\mathrm{SNR}`$ is measured on the signal output detector. Therefore, the signal in- and output $`\mathrm{SNR}`$s can be measured with the same two-port detector. The signal transfer coefficient $`\mathrm{T}_S`$ is thus determined independently of detector efficiencies as the optical $`\mathrm{SNR}`$ transfer coefficient. The probe $`\mathrm{SNR}`$ is measured with the probe detector. As a consequence, the measured probe $`\mathrm{SNR}`$ has to be corrected for the small difference in signal and probe detector efficiencies of $`6.3\%`$. In a test, a weak phase modulation of the signal showed no $`\mathrm{SNR}`$ transfer to the probe. For an optimized filter, mainly the amplitude fluctuations of the signal couple to the detected part of the probe, despite the slight phase sensitivity of the collision. This also shows that the experiment is unaffected by thermal or quantum phase noise.

The best transfer coefficients measured were $`\mathrm{T}_S+\mathrm{T}_P=0.99+0.63=1.62(\pm 0.15)`$ at the filter settings of $`\lambda _{\mathrm{sep}}=1504`$ nm; $`\lambda _{\mathrm{filt}}=1492`$ nm. This was measured with an $`\mathrm{SNR}`$ of the incoming signal of 103, almost identical signal-out $`\mathrm{SNR}`$ and a probe-output $`\mathrm{SNR}`$ of 63. This high value of $`\mathrm{T}_S+\mathrm{T}_P`$ demonstrates the optical tapping property of the device.

Simultaneously optimizing the apparatus for both, a low conditional variance and high transfer coefficients, we observed $`V_{S/P}=0.92\pm 0.05`$ and $`\mathrm{T}_S+\mathrm{T}_P=1.37\pm 0.13`$ ($`\lambda _{\mathrm{sep}}=1503`$ nm; $`\lambda _{\mathrm{filt}}=1490`$ nm). This result clearly fulfills both the QND criteria. Depending on the relative phase of signal and probe as well as the filter positions ($`\lambda _{\mathrm{sep}}`$ and $`\lambda _{\mathrm{filt}}`$) the achieved transfer coefficients and conditional variances differ. The data obtained are presented in Fig.5. Stable results were obtained in all of the four quadrants. To achieve performance in the QND domain, we separated signal and probe at $`1500`$ nm$`<\lambda _{\mathrm{sep}}<1506.5`$ nm and low-pass filtered the probe at $`1488`$ nm$`<\lambda _{\mathrm{filt}}<1490`$ nm. The experimental results for $`V_{S/P}`$ and $`\mathrm{T}_S+\mathrm{T}_P`$ are not correlated. This suggests that additional noise is introduced to the solitons. Likely candidates are inter- and intra-soliton stimulated Raman scattering . The Raman effect pumps photons towards the red and depletes the probe. This assumption is supported by an observed imbalance in the power of the output pulses of $`\mathrm{\Delta }P/P\pm 8\%`$ and by the observation of excess noise in the signal output (Fig.4). The effect decreases the observed signal transfer coefficient. A second perturbative effect is the evolving dispersive wave which modifies the quantum fluctuations significantly, because spectral correlations not only arise between the solitons, but also between them and the continuum . The use of longer pulses would reduce the Raman effect; the dispersive wave can be avoided with improved input pulse shaping.

A full QND experiment can be performed with repeated BAE measurements. When the signal is freely propagated after a first BAE interaction, the soliton will reshape with only little losses  and can be used in a second BAE measurement. A way to nearly eliminate residual radiation loss for cascaded BAE detection is the insertion of a second probe pulse in the center of the collision, prepared to complete the collision.

## V Conclusion

In conclusion, we demonstrated a back-action evading measurement of the photon number based on the novel concept of spectral filtering. This new technique is robust against phase noise limiting previous experiments . It requires merely two pulses since no phase reference pulse is needed and only direct detection is used for the BAE readout. The experiment showed large capacity for optical tapping as well as the capability for quantum state reduction. The idea of coupling to a completely different degree of freedom in the probe, neither conjugate nor identical to the signal observable may be utilized to improve BAE or QND measurements in fibers as well as in other $`\chi ^{(3)}`$ and $`\chi ^{(2)}`$ systems. Applications can be found, e. g. in quantum information  or weak absorption measurements . Investigating the quantum properties of soliton collisions explores ultimate bounds in wavelength-division multiplexed optical transmission systems .

## Acknowledgement

We extend special thanks to Ch. Silberhorn, B. Pfeiffer and M. Meissner for their help in the data acquisition. We gratefully acknowledge fruitful discussions with P. K. Lam, S. Spälter and N. Korolkova. This project was supported by the Deutsche Forschungsgesellschaft.


text:
# Coupling parameter in synchronization of diluted neural networksPublished in Phys. Rev. E65, 041916(2002)

## Abstract

We study the critical features of coupling parameter in the synchronization of neural networks with diluted synapses. Based on simulations, the exponential decay form is observed in the extreme case of global coupling among subsystems and fully linking in each network: there exists maximum and minimum of the critical coupling intensity for synchronization in this spatially extended system. For the partial coupling, we present the primary result about the critical coupling fraction for various linking degrees of networks.

Synchronization of coupled complex systems has been an intensively studied subject since the pioneering work of Fujisaka and others. This phenomenon of synchronization is observed in many other fields, such as in neural networks, in biological populations and in chemical reactions. Recently, spatially extended system has inspired great interest.

Following the series of work contributed by Zanette, we consider a simple modified version of the neural network model described in . As a very important ubiquitous factor corresponding to real biotic neural systems, the neural activity and morphology of synaptic connectivity i.e. the dilution of neural networks must be introduced. That is, in natural neural systems, not all the neurons are linked together. So, there is a chance to investigate the critical features of coupling parameter and the function of structural topology in synchronization of extended systems.

We consider a neural network model that consists of $`N`$ analog neurons $`x_i\left(t\right)[0,1],i=1,\mathrm{},N`$. Each neuron $`x_i`$ is connected with other neurons $`x_j`$ by a random weighted coupling $`J_{ij}`$. Obviously, the connecting matrix $`𝐉`$ is asymmetric and the neural systems exhibit complex spatial oscillations. As a simple version of coupled neural systems designed by Zanette, we use the parallel dynamics for the updating neurons:

$$\begin{array}{c}x_i^1\left(t+1\right)=\left(1\epsilon \right)\mathrm{\Theta }\left(h_i^1\left(t\right)\right)+\epsilon \mathrm{\Theta }\left(h_i^1\left(t\right)+h_i^2\left(t\right)\right)\hfill \\ x_i^2\left(t+1\right)=\left(1\epsilon \right)\mathrm{\Theta }\left(h_i^2\left(t\right)\right)+\epsilon \mathrm{\Theta }\left(h_i^1\left(t\right)+h_i^2\left(t\right)\right)\hfill \end{array}$$
(1)
Here $`h_i^k(t)`$ is the local field of the $`i`$-th neuron and is expressed by

$$h_i^k(t)=\underset{j=1}{\overset{N}{}}C_{ij}J_{ij}x_j\left(t\right)$$
(2)
where $`C_{ij}\{0,1\}`$ is used to denote the linking status between the $`i`$-th neuron and the $`j`$-th neuron. The dilution factor $`C_{ij}`$ is independent identically distributed random variable. It is selected by

$`\mathrm{if}zd,\mathrm{then}C_{ij}=1;\mathrm{else}C_{ij}=0`$

where $`z[0,1]`$ is a random number and $`d[0,1]`$ denotes the linking degree of networks. The activation function $`\mathrm{\Theta }\left(r\right)`$ is defined as $`\mathrm{\Theta }\left(r\right)=\left[1+\mathrm{tanh}\left(\beta r\right)\right]/2`$. In there, $`\beta 1/T`$ characterizes a measure of the inverse magnitude of the amount of noise affecting this neuron, acting as the role of reciprocal of temperature in analogy to thermodynamic systems. For convenience, we set $`\beta =10`$ through all simulations.

Obviously, the first term on the right-hand side of Eq. (1) pictures the total response from its own internal units. The second term expresses the interaction of the summation of the received signals from the neurons with the same position in two networks. The parameter $`\epsilon [0,1]`$ named as coupling intensity, describes the interaction degree between coupling subsystems. When the intensity $`\epsilon 1`$, it is easy to see that the coupling sub-systems evolute independently. On the other limit case, $`\epsilon 1`$, the coupling subsystems are governed by the same dynamical law and will be synchronized very easily.

For measuring the coherence in the collective activity of the neural systems, a time-dependent important feature $`u_i\left(t\right)=_{k=1}^Nx_k^i\left(t\right)`$ for each network $`i=1,2`$ is introduced. When the global coupling of two systems is absent, $`\epsilon 0`$, the $`u_i(t)`$ will update independently and impossibly get to synchronization since the initial conditions in subsystems are different. On the other hand, the activity signals of two subsystems will be identical if the coupling systems come to be synchronous. Figure $`1\left(a\right)`$ shows that the synchronization for this extend system successfully takes place in $`t=275`$ for systems with $`N=100`$, while the coupling intensity $`\epsilon =0.34`$ and the linking degree in each subsystem $`d=0.2`$.

To show the degree of synchronization in this coupling system, the dispersion of activity patterns is defined as

$$D\left(t\right)=\frac{1}{2}\underset{i=1}{\overset{2}{}}\underset{k=1}{\overset{N}{}}\left[x_k^i\left(t\right)\overline{x}_k\left(t\right)\right]^2$$
(3)
where $`\overline{x}_k\left(t\right)=2^1_{i=1}^2x_k^i\left(t\right)`$ denotes the average activity of neurons occupying the $`k`$-th position in both subsystems at time $`t`$. Figure $`1\left(b\right)`$ shows that the dispersion with a logarithmic scale evolves in time with the same synchronous conditions in Fig. $`1\left(a\right)`$.

It is obvious that the larger coupling intensity $`\epsilon `$ makes the more easily synchronization arises for the system with the same other parameters. Concomitantly, the case is whether there exists a critical coupling parameter $`\epsilon _c`$, and furthermore, whether there exists a dependent relationship between $`\epsilon _c`$ and the topological structure in subsystems. In fact, since the evolution of networks is sensitive to the varied initial status and the different random connecting weight matrix in system with the same linking degree $`d`$, it is impossible to find an identical value of $`\epsilon _c`$. However, the fact that the dispersion for these $`\epsilon _c`$ corresponding to varied initial status and connecting matrices from our simulations is small brings our notice to investigate the qualitative curve of $`\epsilon _c`$ vs. $`d`$.

Figure $`2\left(a\right)`$ shows a plot of the critical coupling intensity $`\epsilon _c`$ versus the linking degree $`d`$ in subsystems with the size $`N=200`$. One can see that the qualitative relation between $`\epsilon _c`$ and $`d`$ is close to a sigmoidal curve. The larger $`d`$ arise, the larger $`\epsilon _c`$ become. This can be explained that the evolution of subsystems with larger $`d`$ is more stable and it needs more powerful coupling parameter to drive their evolutions into synchronization. In Fig. $`2\left(b\right)`$, we present the plot of simulation with the stepsize of linking degree $`\mathrm{}d=0.001`$ in the same conditions of Fig. $`2\left(a\right)`$. Comparing both plots of simulations, the agreement is excellent for the global tendency of the qualitative behavior of $`\epsilon _c`$ versus $`d`$.

From Fig. $`3`$, it follows that it is more difficult to come into synchronization with the increase of the size of subsystems. It is clear, however, that there exists a homologous asymptotic behavior in the area of larger linking degree. Now, the case is how the limit of coupling intensity depends on varied size of networks. In Fig. $`4`$, we present the plot of $`\epsilon _c`$ versus $`1/N`$ for the linking degree $`d=1`$ which is identical to the limit case. The form of the limit coupling intensity as a function of the inverse of size of networks calls for a fitting of these data with a exponential decay function

$$\epsilon _c=A+Be^{1/(CN)}$$
(4)
where the constant $`A=0.44\pm 0.024`$, $`B=0.44\pm 0.020`$ and $`C=0.0066\pm 0.00071`$. It follows that the maximal critical coupling intensity corresponding to $`N\mathrm{}`$ is set as $`0.88\pm 0.044`$, and vice versa, the minimal $`\epsilon _c`$ is $`0.44\pm 0.024`$ if both global connecting subsystems designed by Eq. (1) can be come to synchronization.

Another important topic is the fraction of coupling neurons between two subsystems. The considered coupling system can be viewed as a structure made of two horizontal layers of networks. Apparently, from the definition Eq. (1) of the above investigated systems, the neurons are involved in global vertical coupling interactions between two layers, or the dimensionality of coupling parameter is identical to the size of subsystems. Considering the real physical systems or the potential applications, the coupling interactions must be diluted and modified with time. As a result, the systems defined by Eq. (1) can be redefined as

$$\begin{array}{c}x_i^1\left(t+1\right)=\left(1\epsilon \xi _i\left(t\right)\right)\mathrm{\Theta }\left(h_i^1\left(t\right)\right)+\epsilon \xi _i\left(t\right)\mathrm{\Theta }\left(h_i^1\left(t\right)+h_i^2\left(t\right)\right)\hfill \\ x_i^2\left(t+1\right)=\left(1\epsilon \xi _i\left(t\right)\right)\mathrm{\Theta }\left(h_i^2\left(t\right)\right)+\epsilon \xi _i\left(t\right)\mathrm{\Theta }\left(h_i^1\left(t\right)+h_i^2\left(t\right)\right)\hfill \end{array}$$
(5)
where $`\xi _i\left(t\right)\{0,1\}`$ is a random number with probability $`1p`$ and $`p`$, respectively.

For revealing the association between the critical coupling fraction $`p_c`$ and the coupling intensity $`\epsilon `$, the qualitative diagram of $`p_c`$ versus $`\epsilon `$, is shown based on numerical simulations for systems with $`N=100`$, $`d=0.5`$ (see figure $`5\left(a\right)`$). It is easy to get an acceptable conclusion that $`p_c`$ decreases with the increase of $`\epsilon `$ . Note that the series of turning points corresponding to $`p_c=1`$ in Fig. $`5\left(b\right)`$ are equivalent to the points in plot of $`\epsilon _c`$ versus $`d`$ (cf. Fig. $`3`$).

In addition, another valuable informations about the minimal critical fraction for synchronization of this extend system denoted by Eq. (5) can be revealed from another critical point at $`\epsilon =1`$ in Fig. $`5\left(a\right)`$ and Fig. $`5\left(b\right)`$. In Fig. $`(6)`$, the minimal critical coupling probability for various linking degree of subsystem with $`N=400`$ is presented. During the evolution of networks, the neuron of each site in both networks updates due to the competing effect of the local rules and the coupling mechanism. It is clear that the region above the curve is the synchronization part, while the lower part is desynchronizatious, thus the curve embodies a competing relationship between local correlation and stochastic coupling. It is possible to give the minimal critical coupling fraction for various linking degree in synchronization of this extended system with $`N\mathrm{}`$ , which is more analogous to the case of real biotic systems, through analyzing curves of $`p_c`$ vs. $`d`$ for various size of subsystem. However, considering our computational device, the more intensive and detail work is left out in there.

In this paper, we have studied the critical features of coupling parameter in the synchronization of neural networks for various structural topology. We obtain the exponential decay form in the case of global coupling among subsystems and fully linking in each network. We find that it exists the maximal and minimal critical coupling intensity for synchronization in this extend systems. For the case of partial coupling, a primary result about the critical coupling fraction for various linking degrees of networks is shown. Considering the definition of our model is analogous to coupled map lattice, it is easy to generalize the present work to other extended systems, such as coupled ordinary differential equations and partial differential equations.

This work was supported by the Doctoral Research Foundation awarded by Lanzhou University and Innovation Project of CAS with Grant No. kzcx1-09. We wish to thank Prof. H. Zhao for the most constructive and fruitful discussions.

FIG. 1. The synchronization of diluted networks with $`N=100`$, $`\epsilon =0.34`$ and $`d=0.2`$. (a) The evolution of time-dependent activity of both subsystems comes into synchronization at $`t=275`$. (b) logarithm of dispersion of both networks.

FIG. 2. The simulations of relationship between the critical coupling parameter $`\epsilon _c`$ and the linking degree $`d`$ with $`N=200`$, corresponding to stepsize of linking degree (a) $`\mathrm{}d=0.05`$ and (b) $`\mathrm{}d=0.001`$.

FIG. 3. The qualitative relationship of $`\epsilon _c`$ versus $`d`$ for varied size of networks $`N`$.

FIG. 4. The plot of $`\epsilon _c`$ versus $`1/N`$ for varied size of networks in the limit case $`d=1`$. The exponential decay fitted curve of this relationship is shown by a dotted line.

FIG. 5. The synchronization diagram of minimal coupling probability $`p_c`$ versus the corresponding coupling intensity $`\epsilon `$, (a)for systems with size $`N=100`$, $`d=0.5`$ (b) for systems with varied linking degree $`d`$ and $`N=100`$ in the forms of qualitative curves.

FIG.6. The qualitative relationship of $`p_c`$ versus $`d`$ with the subsystem size $`N=400`$.