text:
# The radio flaring behaviour of GRO J1655-40: an analogy with extragalactic radio sources?

## 1 Introduction

Microquasars, here defined as X-ray binaries containing a black hole accreting material from a companion star, are generally regarded as the Galactic counterparts of more powerful extragalactic radio sources, such as blazars and radio galaxies. 111The blazar class comprises the brightest and most variable of the extragalactic radio sources. Accordingly, for the purposes of this study, most observational data for extragalactic sources that allow comparison with the microquasars pertain to blazars. However, the ‘unified scheme’ for radio sources postulates that, to first order, the only difference between radio galaxies, quasars and blazars is one of orientation (e.g. Barthel 1989) so the same physical model should be applicable to all. This implies that the underlying physical processes that operate in these systems are the same, and thus are manifested in their radiative properties. At radio frequencies, the main argument for microquasars being Galactic miniatures of quasars is based on imaging observations of relativistic outflows in the former. This is further strengthened by the fact that some microquasars also display superluminal motion analogous to their extragalactic counterparts. For recent reviews of the radio properties of microquasars, see e.g. Mirabel & Rodriguez (1998, 1999) and Fender (2001, 2002).

The two unambiguous Galactic superluminal sources to date are GRS 1915+105 and GRO J1655$``$40. Their radio jets were discovered in the mid-1990s (Hunstead et al. 1994; Reynolds et al. 1994; Mirabel & Rodríguez 1994; Tingay et al. 1995; Hjellming & Rupen 1995). They often show radio flares during the X-ray active states, and their radio emission is polarized. The nature of GRS 1915+105 has been controversial but recent observations show that it is likely to be a low mass X-ray binary with a black-hole candidate and K-giant donor star (Greiner et al. 2001). It is certain that GRO J1655$``$40 is a binary system containing an evolved F-type star and a black hole (Orosz & Bailyn 1997; Soria et al. 1998).

Compact extragalactic radio sources usually show power law spectra at high-frequencies ($`>90`$ GHz). The radio emission is of synchrotron origin, is often highly polarized, and the power law spectrum implies a non-thermal energy distribution for the relativistic electrons and an optically thin emission region. The spectra often have a turn-over at low frequencies that is generally explained as the consequence of optical-depth effects; the emission is absorbed and the source becomes opaque below a certain critical frequency. The time evolution of the low-frequency (typically less than a few GHz) radio spectra was initially explained by means of a homogeneous cloud composed of relativistic electrons and magnetic field that expands adiabatically and emits synchrotron radiation (Shklovsky 1965; van der Laan 1966; Pauliny-Toth & Kellerman 1966). However, it was later found that a more satisfactory explanation is provided by models invoking shocks propagating along a relativistic jet, and hence accelerating electrons to relativistic energies (e.g. Marscher & Gear 1985; Hughes, Aller & Aller 1985, 1989a, b).

Like extragalactic radio sources, the Galactic superluminal source GRO J1655$``$40 also shows a power law spectrum at radio frequencies that becomes inverted during flaring (Hannikainen et al. 2000). To date the time-dependent radio properties of X-ray binaries have mostly been ascribed to the expanding-cloud (or synchrotron-bubble) model or variations thereof (Hjellming & Johnston 1988; Han & Hjellming 1992; Ball & Vlassis 1993) although an internal shock model has been applied successfully to a radio outburst of GRS 1915+105 (Kaiser, Sunyaev & Spruit 2000). Since the 1994 radio outburst of GRO J1655$``$40 is not readily accommodated by the expanding-cloud model (see Hannikainen et al. 2000 and Section 2 below) we propose alternatively that the processes that generate the radio emission from GRO J1655$``$40 and from extragalactic radio sources are the same. We make use of the shock model of Marscher & Gear (1985) to explain the spectral evolution of GRO J1655$``$40 during the 1994 outburst. The multifrequency variability is discussed in Section 2, basic features of the model are described in Section 3, and the quasar-microquasar connection is discussed in Section 4.

Throughout this paper, spectral index, $`\alpha `$, is defined as $`S_\nu \nu ^\alpha `$ where $`S_\nu `$ is the flux density at frequency $`\nu `$.

## 2 Time-dependent properties of GRO J1655$``$40

GRO J1655$``$40 flared dramatically in 1994 August, the flux density at 843 MHz increasing from a few hundred mJy to almost 8 Jy in 12 d (Hannikainen et al. 2000). Very Long Baseline Array (VLBA) images at 1.6 GHz taken at several epochs during the outburst show well-collimated relativistic jets. The jets are resolved into several knot-like features which expand outwards from a stationary compact core (Hjellming & Rupen 1995).

The multifrequency flux density light curves and spectra (from Hannikainen et al. 2000) are shown in Figs 1 and 2 respectively. The flux densities include the emission from the entire structure shown in the VLBA images, with the possible exception of the two highest frequencies at which the Australia Telescope Compact Array (ATCA) synthesized beam size is approximately 1 arcsec compared to a source size of 1–2 arcsec at 1.6 GHz. Our flux density measurements involve fitting a point-source model to the data which may lead to a slight systematic underestimate of the flux at 8.6 and 9.2 GHz. However, the source may well be more compact at higher frequencies due to the shorter synchrotron lifetimes of the emitting electrons (see Section 4), so it is unlikely that much/any flux density is missing from our measurements. Indeed, in Section 4, we argue that most of the flux density variability is confined to a region $`1`$ arcsec in extent. Since the available observations do not allow separation of emission from the compact variable component(s) and the extended regions we are forced to assume that the radio emission is dominated by a single variable component. The light curves show complex flaring behaviour with several minor events following the initial increase in flux density. The key features of the data are that the initial rise in flux density peaks simultaneously at all frequencies and that the amplitude of the flare, defined as maximum minus minimum flux density, increases towards lower frequency (see Section 4 for further discussion). We quantified the time delays between the data streams with the Discrete Correlation Function (DCF; Edelson & Krolik 1988) and the Interpolated Correlation Function (ICF; Gaskell & Peterson 1987). Both methods give essentially the same result (Fig. 3); the correlation functions are approximately symmetrical around zero lag in all cases. The flat peak to the curves suggest formal time delays of $`0\pm 1`$ d.

The first panel of Fig. 2 shows a spectrum from two days before flux density maximum, which is consistent with synchrotron emission that is optically thin at the highest frequencies but is partially self-absorbed at around 1 GHz. As the flux density increases, the entire spectrum becomes optically thin and remains so throughout the subsequent decline in flux density. In Fig. 2 we can see that the spectrum was relatively flat (with a spectral index of $`0.41`$ on TJD 9580). It steepened at the time of the maximum (with the spectral indices reaching $`0.65`$ on TJD 9582) and then flattened as the flux density declined to a level of $`0.4`$ around TJD 9586 and 9587. This behaviour can be explained by models incorporating emission from multiple regions. In the simple core-lobe model suggested by Hannikainen et al. (2000), the observed spectrum is the sum of the emission from a compact region (core), which has a flat spectrum, and the emission from extended regions (lobes/ejecta) with steeper spectra. The outburst can be seen as a consequence of the brightening of the compact region and then the rapid rise in the brightness of the ejecta. The emission from the ejecta eventually dominates as the radio flux densities reach their peak amplitude. At the same time the spectrum steepens. As the flare subsides, the emission from the ejecta drops to a level comparable to, or lower than, the emission from the compact region, and the spectra become flatter, similar to those seen at the onset of the flare.

The time-dependent properties of the radio outbursts of GRO J1655$``$40 and other black-hole X-ray binaries have been attributed to synchrotron emission from expanding clouds of relativistic electrons. Under the usual assumption that all radio frequencies are initially optically thick, the two major characteristics of these models are that the low-frequency emission has a time delay with respect to that at high frequencies and that the flare amplitudes fall off towards lower frequency (see e.g. Han & Hjellming 1992). The observations clearly do not conform to these predictions of the expanding-cloud model (see also Hannikainen et al. 2000).

If electrons are injected into an expanding cloud, and if expansion losses dominate over radiative losses, then it is possible to get a simultaneous peak at optically thin frequencies once the injection stops or once the expansion cooling dominates over the injection. Note that if radiative losses dominate under this scenario we would observe frequency-dependent peak times, assuming the data streams are sufficiently well sampled. However, the first panel of Fig. 3 provides evidence that the synchrotron self-absorption turnover was close to 1 GHz near the peak of the flare, implying that at least some of the radio frequencies were optically thick at earlier epochs, and thus arguing against this possibility. Similar scenarios can be produced by relaxing some of the initial assumptions of the model. For example, Ball & Vlassis (1993) considered cases where electrons are injected into the cloud with a constant energy spectrum. While both of these options can in principle reproduce the observed behaviour (optical depth arguments aside) they require a specific set of conditions that are unlikely to be met in practice. Indeed, Marscher & Gear (1985) found that they could fit the variability of 3C 273 with a uniform expanding source model but only if ‘the injection of relativistic electrons was allowed to vary with radius in a rather ad hoc fashion’. They attributed the success of this model to its large number of free parameters. Arguing along similar lines, we do not consider that such models provide a convincing physical explanation for the phenomena observed in GRO J1655$``$40. Instead, given the obvious analogy between Galactic and extragalactic jet sources, we argue alternatively that the variable radio emission arises in a small region behind a relativistic shock wave, and that the rise of the flare is driven by inverse-Compton cooling (cf. Marscher & Gear 1985). We note that the shock model developed by Kaiser et al. (2000) incorporates adiabatic and synchrotron cooling only and as such is mostly applicable to the decline of radio flares rather than to the rise phase discussed here.

## 3 Generalised shock model for extragalactic sources

Spectral evolution during the flares of compact extragalactic radio sources can be explained by the generalised shock model of Marscher & Gear (1985). (For specific applications see e.g. Valtaoja et al. 1988 and Stevens et al. 1994, 1996, 1998). In the model, temporal evolution of the radio emission is divided into growth, plateau and decay stages (see the schematic illustration in Fig. 4). The emission region consists of relativistic jets along which transverse shocks propagate and accelerate electrons to relativistic energies. The shock waves form in response to changes in the conditions in the jet (such as pressure or velocities of the bulk flow). The evolution of the shock is described in terms of the flux density ($`S_\nu `$) at the peak frequency ($`\nu _\mathrm{m}`$ where $`S_\nu [\nu _\mathrm{m}]=S_\mathrm{m}`$) of the synchrotron spectrum where the opacity is close to unity. This point is fixed at any one time by the dominant energy-loss mechanism. When the emitting region is compact, inverse-Compton losses predominate (Compton or growth stage) but these fall off rapidly with radius as the shock expands and are superseded by synchrotron losses (synchrotron or plateau stage). As the shock expands further, the radiative lifetime of the electrons becomes large with respect to the time needed to traverse the emitting region and losses due to adiabatic expansion (adiabatic or decay stage) become more important. All three stages are approximated by power laws on the logarithmic $`(S_\mathrm{m},\nu _\mathrm{m}`$) plane. As the emitting region expands, $`\nu _\mathrm{m}`$ is predicted to move to lower frequencies with time. $`S_\mathrm{m}`$ increases rapidly during the Compton stage, remains approximately constant during the synchrotron stage and decreases during the adiabatic stage.

The characteristics of spectral evolution predicted by the generalised shock model are summarised as follows (see Valtaoja et al. 1992 for more details). The maximum flare amplitude for any frequency on the Compton stage occurs when the spectrum transits onto the synchrotron stage. Light curves at such frequencies are thus predicted to peak simultaneously and because of the spectral shape, the flare amplitudes are expected to increase towards lower frequency – specifically, they should have the same power law form as the optically thin portion of the flare spectrum. Flare amplitudes at frequencies commensurate with the later stages are determined by the details of the spectral evolution, and are thus expected to be approximately constant during the synchrotron stage and to decrease during the adiabatic stage. The light curves will display time-lagged behaviour with a delay between any two frequencies equal to the time taken for the spectrum to evolve between them and become optically thin. Note that the development of the flare during the adiabatic stage is very similar to that predicted by the expanding cloud model, although, because the expansion is constrained to occur in a jet, the predicted time delays are shorter and the variation of flare amplitude with frequency is less pronounced.

## 4 Microquasars vs quasars

The generalised shock model prediction of simultaneously peaking emission at all frequencies coupled with flare amplitudes that increase towards lower frequencies (i.e. the Compton/growth stage) is consistent with the observations of GRO J1655$``$40 in 1994. According to the shock model, we would also expect the flare amplitudes to fall off as a power law since they should have the same frequency dependence as the spectrum of the flare when it reaches the transition point between the growth and plateau stages. In Fig. 5 we plot the flare amplitudes normalised to 4.8 GHz. A power law is clearly a good fit with a slope of $`0.70\pm 0.05`$, consistent with that expected for the optically thin portion of the flare synchrotron spectrum (cf. Stevens et al. 1996).

One apparent difference between the flaring behaviour of GRO J1655$``$40 and its extragalactic analogues, the blazars, is that the flare remains in the growth stage to much lower frequencies. For example, many of the sources in the sample considered by Stevens et al. (1994) show time delays between the 37 GHz, or occasionally the 90 GHz, emission and that at higher frequencies, although those authors noted that some BL Lacertae objects remained in the growth stage out to 4.8 GHz. An obvious contributor to this difference is the Doppler effect which will shift the emitted frequency of radiation in the jet frame, $`\nu ^{}`$, to a higher observed frequency, $`\nu =\delta \nu ^{}`$. The relativistic Doppler factor, $`\delta =\mathrm{\Gamma }^1(1\beta \mathrm{cos}\theta )^1`$ where $`\mathrm{\Gamma }=(1\beta ^2)^{1/2}`$ is the bulk Lorentz factor, $`\theta `$ is the viewing angle and $`\beta `$ is the jet speed in units of $`c`$. It is generally accepted that blazars have their jet axes oriented quite closely towards the observer; for example Teräsranta & Valtaoja (1994) find that blazars on average are aligned to within about 20 deg of the line-of-sight (see also Barthel et al. 1989). Estimates of Doppler factors for blazars range from extreme values of 0.005 to 33 with a mean value of around 5 (Guijosa & Daly 1996). The blazars with measured time delays from the sample of Stevens et al. (1994) have ‘inverse-Compton’ Doppler factors in the range 3.4–16 (see Guijosa & Daly 1996). The microquasar GRO J1655$``$40, however, is angled further towards the plane of the sky ($`\theta =70`$ deg; Orosz & Bailyn 1997) and its estimated jet speed $`\beta =0.92`$ (Hjellming & Rupen 1995) leads to $`\delta 0.6`$. Higher values of the Doppler factor for blazars compared to that of GRO J1655$``$40 could thus easily lead to an order of magnitude difference in the frequency at which the flare exits the growth stage.

A second contributory factor might arise from physical differences between the jets of Galactic and extragalactic sources. Since the flare evolution is dependent on the photon energy density, the magnetic energy density and the lifetime of the emitting electrons compared to the time they take to cross the shock structure, two important parameters are the size of the emitting region and its magnetic field strength.

We can estimate the magnetic field from the flare decay time by assuming that the decay is radiative rather than expansion-loss driven. For the clear-cut case of a single flare we would expect the 1/e decay times, $`t`$, of the light curves to vary as $`\nu ^{1/2}`$ if radiative losses are the dominant energy-loss mechanism, or if adiabatic expansion losses predominate, then $`t`$ will be the same at all frequencies. Unfortunately, the flaring behaviour that we observe is complex. At the monitoring frequencies presented in Fig. 1 the decay of the flares is interrupted by subsequent events, making it impossible to measure an accurate e-folding timescale. However, the multifrequency light curves presented by Hjellming & Rupen (1995) include data up to 22.5 GHz. For the flare in question, the decay of the 22.5 GHz lightcurve is more rapid than those at lower frequencies, most probably because the monitoring frequency is now sufficiently high that the e-fold timescale is similar to the interval between the bursts. This result provides evidence that synchrotron losses are driving the flux decay. If this is not the case then the implication would be that the synchrotron-loss stage is short lived or non-existent, as observed for 3C 273 (Stevens et al. 1998). The rest of this section assumes that the decay of the flare is driven by synchrotron losses. If this is not the case, and adiabatic expansion losses predominate then we will overestimate the magnetic field strength and underestimate the electron Lorentz factor and source extent.

For an e-fold decay timescale $`t`$ and a peak frequency $`\nu _\mathrm{m}`$, the magnetic field is given by

$$B23\delta ^{1/3}\left(\frac{\nu _\mathrm{m}}{1\mathrm{GHz}}\right)^{1/3}\left(\frac{t}{5\mathrm{day}}\right)^{2/3}\mathrm{gauss}.$$
(1)
We estimated the flare decay time from the best sampled light curve (843 MHz). For the reasons discussed above this estimate is likely to be an upper limit but should not be out by more than a factor of two; using 5 d the estimated magnetic field strength is $``$25 gauss. This value is approximately two orders of magnitude higher than typical estimates of 0.1–1 Gauss in extragalactic sources (e.g. Brown et al. 1989).

The Lorentz factor of the electrons emitting at 843 MHz is given by $`\gamma (2\pi m_\mathrm{e}c\nu /eB)^{1/2}3.5`$, where $`e`$ and $`m_\mathrm{e}`$ are the electron charge and mass respectively. This value is approximately that required to explain the observed circular polarization as being intrinsic to the synchrotron radiation (Macquart et al. 2002) although the alternative mechanism, Faraday conversion, could also apply. In any case, we note that the presence of circularly polarized radiation is another characteristic that GRO J1655$``$40 shares with the blazars (see e.g. Wardle et al. 1998). Furthermore, strong evidence for our proposed model comes from the observed position angle of linear polarization at radio wavelengths which is very close to that of the jet direction (Hannikainen et al. 2000). For optically thin synchrotron radiation the implied magnetic field direction is perpendicular to the jet, as observed for many blazars (e.g. Cawthorne et al. 1993a,b), and as predicted by models in which a tangled component of magnetic field in the jet is compressed parallel to a transverse shock front (e.g. Hughes et al. 1989a).

Finally, we can estimate the source extent at the time when the spectrum transits from the growth stage onto the plateau stage. At this transition, the energy densities in particles ($`U_{\mathrm{ph}}`$) and magnetic field are equal and thus so are the synchrotron and inverse-Compton lifetimes $`t_{\mathrm{syn}}=t_{\mathrm{IC}}`$. Using

$$t_{\mathrm{IC}}\frac{3\times 10^7}{\gamma U_{\mathrm{ph}}}\frac{3\times 10^7cD^2}{\gamma \varphi L}\mathrm{seconds},$$
(2)
and assuming that, (1) the transition occurs when the 843 MHz light curve peaks (i.e. the flare amplitudes peak at this frequency), allowing us to use the values $`t_{\mathrm{IC}}=t_{\mathrm{syn}}5`$ d and $`\gamma 3.5`$ deduced above, (2) the radiation field is produced by a central source of luminosity, $`L10^{38}`$ erg s-1 and (3) that $`\varphi `$, a geometrical constant, is approximately unity, the calculated source size $`D10^{13}`$ cm. At a source distance of 3.2 kpc (Hjellming & Rupen 1995) this dimension translates to an angular size of about 0.2 mas. A 1.6 GHz VLBA map taken one day after the peak of the 843 MHz light curve shows a source of about 1 arcsec with a bright core and extended jets (Hjellming & Rupen 1995). Our proposed model requires that the radio flares originate in a small region close to the base of the jet as is observed in AGN.

## 5 Final remarks

We have shown that a model incorporating a propagating relativistic shock in a jet can qualitatively reproduce the multifrequency variability of GRO J1655$``$40. This model, which also fits the available polarization observations, is the canonical model used to explain the variability characteristics of extragalactic radio sources such as blazars. A detailed discussion of the applicability of this model to Galactic superluminal sources is beyond the scope of the present work, and in any case is limited by the availability of suitable data sets. It is important that other flares similar to the one discussed here are observed over a broad range of radio frequencies, and that they are first observed on the rising portion of the light curves. In this respect, we point out that the two sources discussed by Ball & Vlassis (1993), Nova Muscae 1991 and the Galactic Centre Transient, for which data were taken during the rise of the flares, both display the behaviour predicted by the Compton stage of the shock model. Both sources were observed to have an optically thin spectrum during the rising phase of the outburst which can be accommodated by the Compton stage if the self-absorption turnover had already passed the lowest observing frequency when the flux monitoring began.

## Acknowledgments

J.A.S. and D.C.H. acknowledge support from PPARC, and K.W. acknowledges support from the Australian Research Council through an Australian Research Fellowship. D.C.H also acknowledges travel support from the Academy of Finland, and thanks the University of Sydney and Jodrell Bank for their hospitality and financial support during her visits.


format: text

text:
# 1 The nucleon structure, as seen by a high energy, large 𝑄², probe. The nucleon breaks into unobserved final particles 𝑋; the ejected parton may fragment into an observed hadron.

What do we know about the proton spin structure? 111Talk delivered at the 2nd International Symposium on the Gerasimov-Drell-Hearn sum rule and the spin structure of the nucleon, GDH 2002, July 3-6 2002, Genova, Italy

Mauro Anselmino
Dipartimento di Fisica Teorica, Università di Torino and
INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy

## Abstract

A brief summary of the theoretical and experimental knowledge of the spin structure of the proton is presented. The helicity distributions of quark and gluons are discussed, together with their related sum rules. The transversity distribution is also introduced with possible strategies for its measurement. Novel spin dependent and $`𝒌_{}`$ unintegrated distribution and fragmentation functions are discussed, in connection with a new and rich phenomenology of transverse single spin asymmetries.

Introduction

The spin nucleon structure – as observed in high energy, short distance interactions – is schematically described in Fig. 1.

The large $`Q^2`$ probe – typically, a virtual photon – “sees” QCD partons, carrying a longitudinal momentum fraction $`x`$, and their interactions, with gluon and $`q\overline{q}`$ pair creation; the information about such a complicated structure is usually collected via measurements of the Deep Inelastic Scattering (DIS) cross-section and storaged in the structure functions which appear in the most general expression of the cross-section. When neglecting weak, parity violating contributions, there are 2 unpolarized ($`F_1`$, $`F_2`$) and 2 polarized ($`g_1`$, $`g_2`$) structure functions: perturbative QCD allows a simple partonic interpretation of $`F_2`$ and $`g_1`$ ($`F_1`$ is related to $`F_2`$ while $`g_2`$ does not have a partonic interpretation). We only consider here the polarized proton structure trying to summarize in a short time and space the main ideas, the most recent results and the open problems; many detailed and comprehensive reviews on the subject can be found in the literature .

The main issues and questions we are going to discuss here are:

* our knowledge about the polarized structure functions $`g_1(x,`$ $`Q^2)`$ and $`g_2(x,Q^2)`$ and about quark and gluon helicity distributions, $`\mathrm{\Delta }q(x,Q^2)`$ and $`\mathrm{\Delta }g(x,Q^2)`$; how well do we know them?
* are fundamental sum rules satisfied and what do we know about quark and gluon orbital angular momentum, $`L_q`$ and $`L_g`$?
* $`\mathrm{\Delta }q(x,Q^2)`$, $`\mathrm{\Delta }g(x,Q^2)`$, $`L_q`$ and $`L_g`$ are not the whole story: how and where do we learn about the transversity distribution $`h_1(x,Q^2)`$?
* could we learn more and understand more from intrinsic $`𝒌_{}`$ unintegrated distribution and fragmentation functions?

The longitudinally polarized proton and the helicity distributions

At NLO in the QCD parton model the structure function $`g_1`$ is given by

$$g_1(x,Q^2)=\frac{1}{2}\underset{q}{}e_q^2\left\{\mathrm{\Delta }C_q\left[\mathrm{\Delta }q+\mathrm{\Delta }\overline{q}\right]+\frac{1}{N_f}\mathrm{\Delta }C_g\mathrm{\Delta }g\right\}$$
(1)
where $`\mathrm{\Delta }q(x,Q^2)`$ and $`\mathrm{\Delta }g(x,Q^2)`$ are respectively the quark (of flavour $`q`$) and gluon helicity distributions; we have, as usual, defined the convolution

$$\mathrm{\Delta }C\mathrm{\Delta }q_x^1\frac{dy}{y}\mathrm{\Delta }C(\frac{x}{y},\alpha _s)\mathrm{\Delta }q(y,Q^2)$$
(2)
and the coefficients functions $`\mathrm{\Delta }C_i`$ have a perturbative expansion

$$\mathrm{\Delta }C_i(x,\alpha _s)=\mathrm{\Delta }C_i^0(x)+\frac{\alpha _s(Q^2)}{2\pi }\mathrm{\Delta }C_i^{(1)}(x)+\mathrm{}$$
(3)
The LO terms are simply

$$\mathrm{\Delta }C_q^0=\delta (1x)\mathrm{\Delta }C_g^0=0,$$
(4)
and the NLO corrections are scheme dependent; typical choices differ in the amount of gluon contribution to the quark singlet distributions, while quark non-singlet distributions are scheme independent . Finally, the $`Q^2`$ evolution of the parton densities obeys the DGLAP evolution equations , and, if known at an initial scale $`\mu ^2`$, the r.h.s. of Eq. (1) can be computed at any perturbative $`Q^2`$ value.

By comparing data on $`g_1(x,Q^2)`$ with Eq. (1) one obtains information on the quark and gluon helicity distributions; the more data one has and the wider the $`x`$ and $`Q^2`$ range is, the more stringent the comparison is. The normal procedure is that of using a simple ansatz for the unknown distribution functions at the initial scale $`\mu ^2`$, with some assumptions regarding the sea quark densities (for example, whether $`SU(3)_F`$ symmetric or not) and some constraints from $`SU(3)_F`$ hyperon decay sum rules on the first moments $`\mathrm{\Delta }q(1,Q^2)_0^1\mathrm{\Delta }q(x,Q^2)𝑑x`$.

In Fig. 2 a most recent analysis of the world data on $`xg_1(x)`$ is shown together with a fit from Ref. , where the resulting helicity distributions can also be seen. Several similar analyses can be found in the literature; a complete list of references is given in Ref. .

Some data on $`xg_2(x)`$ are also available and the most recent ones  are shown in Fig. 3.

Let us shortly comment on these experimental results and the information which they offer.

* We have now good data on $`g_1`$ and $`g_2`$, although not yet comparable with the amount and quality of similar data obtained on the unpolarized structure functions.
* $`g_1`$ \[see Eq. (1)\] allows to obtain information on linear combinations of ($`\mathrm{\Delta }q+\mathrm{\Delta }\overline{q}`$). We still need a better flavour resolution; this might come from semi-inclusive DIS which gives information on $`_q\mathrm{\Delta }qD_q^h`$ rather than $`_qe_q^2(\mathrm{\Delta }q+\mathrm{\Delta }\overline{q})`$, where $`D_q^h`$ is the quark $`q`$ fragmentation function into the observed hadron $`h`$. Flavour separation in inclusive DIS would naturally be possible in neutrino iniziated charged current processes .
* Eq. (1) also offers indirect (via QCD evolution) information, on $`\mathrm{\Delta }g`$. This is not stringent enough and a more direct measurement of the gluon helicity distribution is needed. This might come from the study of spin dependences in processes like $`\mathrm{}p\mathrm{}+2jets`$, $`\mathrm{}p\mathrm{}+c+\overline{c}+X`$, $`pN\gamma +X`$, etc. which could be performed at HERMES, COMPASS, RHIC.

Sum rules and orbital angular momentum

In extracting information from experimental data – or in testing theories – a special role is plaid by sum rules. Let us mention a few of them.

The Bjorken sum rule ($`g_A/g_V=1.2670\pm 0.0035`$),

$`{\displaystyle _0^1}[g_1^p(x,Q^2)g_1^n(x,Q^2)]𝑑x`$
$`=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle \frac{g_A}{g_V}}\left\{1{\displaystyle \frac{\alpha _s}{\pi }}{\displaystyle \frac{43}{12}}{\displaystyle \frac{\alpha _s^2}{\pi ^2}}20.2{\displaystyle \frac{\alpha _s^3}{\pi ^3}}+\mathrm{}\right\},`$
is used in many ways. One can simply assume the validity of Eq. (S0.Ex1) and deduct from data on $`g_1`$ the value of $`\alpha _s`$ ; or one can use in it a value of $`\alpha _s`$ otherwise obtained, to check whether data on $`g_1`$ obey the sum rule or not (the answer is yes); or, also , one can assume to know the r.h.s. of Eq. (S0.Ex1), and see – among the poorly known behaviours of $`g_1(x)`$ at small $`x`$ – which one best satisfies the sum rule.

Another, more debated, sum rule is the so called Burkhardt-Cottingham sum rule, according to which $`_0^1g_2(x,Q^2)𝑑x=0`$, provided the integral exists. The recent E155 data  of Fig. 3 seem to indicate $`_{0.02}^{0.8}g_2^p(x,Q^2)𝑑x=0.042\pm 0.008`$, which, taking into account uncertainties in the extrapolation to $`x=0`$ and $`x=1`$, might be the first indication of a violation of the sum rule, assuming that no $`\delta `$-function contributes at the origin .

The last, fundamental sum rule which we mention is the spin sum rule:

$$\frac{1}{2}=\frac{1}{2}\mathrm{\Delta }\mathrm{\Sigma }(1)+\mathrm{\Delta }g(1)+L_q+L_g$$
(6)
where $`\mathrm{\Delta }\mathrm{\Sigma }(1)`$ is the first moment of $`_q[\mathrm{\Delta }q(x)+\mathrm{\Delta }\overline{q}(x)]`$ and $`L_{q,g}`$ is the third component of the orbital angular momentum carried by quarks, gluons. This last quantity is unavoidable in a picture of the proton like that of Fig. 1: a spin 1/2 massless quark can emit a spin 1 massless gluon, via a helicity conserving coupling, only if some orbital angular momentum restores the total angular momentum conservation. However, there is little agreement at the moment both about the proper formal definition of a $`\widehat{L}_{q,g}`$ operator and about a possible measurement of its expectation value between proton states .

Transversity distribution

The trasverse polarization of quarks inside a trasversely polarized nucleon, denoted by $`h_1`$, $`\delta q`$ or $`\mathrm{\Delta }_Tq`$, is a fundamental twist-2 quantity, as important as the unpolarized distributions $`q`$ and the helicity distributions $`\mathrm{\Delta }q`$. It is given by

$$h_1(x,Q^2)=q_{}^{}(x,Q^2)q_{}^{}(x,Q^2),$$
(7)
that is the difference between the number density of quarks with transverse spin parallel and antiparallel to the nucleon spin. It is the same as the helicity distribution only in a non relativistic approximation, but it is expected to differ from it for a relativistic nucleon.

When represented in the helicity basis (see Fig. 4) $`h_1`$ relates quarks with different helicities, revealing its chiral-odd nature. This is the reason why this important quantity has never been measured in DIS: the electromagnetic or QCD interactions are helicity conserving, there is no perturbative way of flipping helicities and $`h_1`$ decouples from inclusive DIS dynamics, as shown in Fig. 4a.

However, it can be accessed in semi-inclusive DIS, where some non perturbative chiral-odd effects may take place in the non perturbative fragmentation process, Fig. 4b. Indeed, a serious program to measure $`h_1`$ in semi-inclusive DIS at HERMES, where a transversely polarized proton target is now available, is in progress. A similar program, in different, complementary, kinematical regions, is planned at COMPASS.

The transversity distribution is also accessible at RHIC, where transversely polarized proton beams are available; by measuring double transverse spin asymmetries in Drell-Yan processes one obtains an observable which depends on the convolution of two transversity distributions, which might make the overall effect rather tiny . In general, $`h_1`$ must appear in a physical observable coupled to another chiral-odd quantity, which is either the transversity itself or a new unknown function.

Unintegrated distribution and fragmentation functions

We conclude by mentioning a new phenomenological approach to the description of many single transverse spin asymmetries which have been measured and keep being measured, with unexpected and interesting results . The apparent problem with these asymmetries is related to the fact that, within perturbative QCD and the collinear factorization scheme, they should be vanishing, which is not true experimentally.

Recently, a series of papers  have shown how single spin asymmetries may occurr at the level of parton distributions and fragmentations, provided one takes into account the intrinsic motion of partons inside hadrons and of hadrons relatively to the fragmenting parton. For example, there might be a correlation between the transverse spin of a quark and the $`𝒌_{}`$ of a resulting hadron, say a pion. This is the so-called Collins effect , pictorially shown in Fig. 5.

Similar spin-$`𝒌_{}`$ correlations may occurr also in the fragmentation of an unpolarized quark into a polarized hadron (the so-called polarizing fragmentation functions ), in the distribution of unpolarized quarks inside polarized nucleons (the Sivers effect ) and in the distribution of polarized quarks inside an unpolarized hadron .

When generalizing the factorization scheme with the inclusion of intrinsic $`𝒌_{}`$, both in the distribution/fragmentation functions and in the elementary interactions, single transverse spin asymmetries appear immediately as possible and even sizeable. A phenomenological approach can be developed in which experimental information on the new functions is obtained from some processes and then used to make predictions in other cases.

The spin structure of the nucleon is subtle and challenging. Enormous progress has been achieved in the last years; yet, new surprising experimental results keep beeing obtained and fresh, interesting ideas keep being suggested. A lot more good work, both experimental and theoretical, is in progress.

Acknowledgments I would like to thank the organizers of the Symposium for the invitation and for the beautiful and stimulating organization.


text:
# Sur la “solution analytique générale” d’une équation différentielle chaotique du troisième ordre 11footnote 1 IRMA Lectures in Mathematics and Theoretical Physics (de Gruyter, Berlin, 2003). Journées de calcul formel, IRMA, Strasbourg, 21–22 mars 2002. Correspondance : RC. S2002/070. nlin.PS/0302056

## 1 Situation du problème

D’après un résultat classique , un système différentiel ne peut présenter un comportement chaotique que si son ordre égale au moins trois. D’autre part, malgré l’opposition intuitive entre chaos et analyticité, il n’est pas interdit à un système génériquement chaotique de posséder des solutions particulières globalement analytiques, la seule impossibilité étant que le nombre de constantes d’intégration ne peut égaler l’ordre. Le but de cet article est de passer en revue les diverses méthodes connues permettant de caractériser (idéalement d’exhiber) cette solution analytique, forcément particulière, qui dépend du nombre maximal possible de constantes d’intégration. Convenons de l’appeler solution analytique générale, par opposition à la solution générale dont l’expression analytique n’existe pas.

L’exemple physique optimal pour illustrer ces méthodes est à notre avis l’équation de Ginzburg et Landau cubique unidimensionnelle (CGL3)

$`i\psi _t+p\psi _{xx}+q|\psi |^2\psi i\gamma A=0,pq\gamma 0,(\psi ,p,q)𝒞,\gamma ,`$ (1)
pour les raisons suivantes,

1. elle est génériquement chaotique ,
2. elle admet une limite intégrable $`\mathrm{}(p)=\mathrm{}(q)=\gamma =0`$ (équation de Schrödinger nonlinéaire (NLS)) dont tous les éléments d’íntégrabilité sont connus analytiquement,
3. les phénomènes physiques qu’elle modélise sont très importants : propagation du signal dans une fibre optique , intermittence spatio-temporelle en turbulence, supraconductivité, etc.

Nous nous contenterons d’étudier ici une équation plus simple, celle de Kuramoto et Sivashinsky (KS), déduite de celle pour la phase $`\phi =\mathrm{arg}\psi `$ du champ $`\psi `$ de CGL3 en prenant une certaine limite ,

$`\phi _t+\nu \phi _{xxxx}+b\phi _{xxx}+\mu \phi _{xx}+\phi \phi _x=0,\phi 𝒞,(\nu ,b,\mu ),\nu 0.`$ (2)
Sa réduction en onde propagative

$`\phi (x,t)=c+u(\xi ),\xi =xct,`$ (3)
$`E(u,\xi )\nu u^{\prime \prime \prime }+bu^{\prime \prime }+\mu u^{}+{\displaystyle \frac{u^2}{2}}+A=0,(\nu ,b,\mu ),\nu 0,`$ (4)
$`A`$ est une constante d’intégration, a également un comportement chaotique , et elle dépend de deux paramètres sans dimension, $`b^2/(\mu \nu )`$ et $`\nu A/\mu ^3`$. Nous noterons désormais $`x`$ pour $`\xi `$. Le problème à résoudre est le suivant.

Problème. Trouver l’expression explicite de la “solution analytique générale” (définition supra) de l’EDO (4).

La motivation vient de ce que cette solution (pour CGL3 et donc pour KS) est “observée” sous forme de texture bien définie, tant dans des simulations numériques que dans de véritables expériences physiques .

## 2 Suppression locale de la contribution chaotique

Comptons d’abord le nombre (nécessairement inférieur à trois) de constantes arbitraires dans la solution analytique générale. La recherche d’un comportement local algébrique au voisinage d’une singularité mobile $`x_0`$ (mobile signifie : qui dépend des conditions initiales),

$`u_{xx_0}u_0\chi ^p,u_00,\chi =xx_0,`$ (5)
conduit à la série de Laurent

$`u^{(0)}=120\nu \chi ^315b\chi ^2+\left({\displaystyle \frac{60\mu }{19}}{\displaystyle \frac{15b^2}{76\nu }}\right)\chi ^1+O(\chi ^0),`$ (6)
d’où sont absentes deux des trois constantes arbitraires. Celles-ci apparaissent en perturbation ,

$`u=u^{(0)}+\epsilon u^{(1)}+\epsilon ^2u^{(2)}+\mathrm{},`$ (7)
où le petit paramètre $`\epsilon `$ ne figure pas dans l’EDO (4). L’équation linéarisée en $`u^{(0)}`$

$`\left(\nu {\displaystyle \frac{\text{d}^3}{\text{d}x^3}}+b{\displaystyle \frac{\text{d}^2}{\text{d}x^2}}+\mu {\displaystyle \frac{\text{d}}{\text{d}x}}+u^{(0)}\right)u^{(1)}=0,`$ (8)
est alors du type de Fuchs au voisinage de $`x=x_0`$, avec pour équation indicielle ($`q=6`$ désigne le degré de singularité du premier membre $`E`$ de (4))

$`\underset{\chi 0}{lim}\chi ^{jq}(\nu _x^3+u_0\chi ^p)\chi ^{j+p}`$ (9)
$`=\nu (j3)(j4)(j5)+120\nu =\nu (j+1)(j^213j+60)`$ (10)
$`=\nu (j+1)\left(j{\displaystyle \frac{13+i\sqrt{71}}{2}}\right)\left(j{\displaystyle \frac{13i\sqrt{71}}{2}}\right)=0.`$ (11)
La représentation locale de la solution générale,

$`u(x_0,\epsilon c_+,\epsilon c_{})`$ $`=`$ $`120\nu \chi ^3\{\text{Taylor}(\chi )`$
$`+\epsilon [c_+\chi ^{(13+i\sqrt{71})/2}\text{Taylor}(\chi )`$
$`+c_{}\chi ^{(13i\sqrt{71})/2}\text{Taylor}(\chi )]+𝒪(\epsilon ^2)\},`$
où “Taylor” désigne des séries convergentes de $`\chi `$, dépend bien de trois constantes arbitraires $`(x_0,\epsilon c_+,\epsilon c_{})`$ (l’indice de Fuchs $`1`$ ne représente qu’une translation de $`x_0`$). Le branchement mobile dense provenant des deux indices irrationnels caractérise  le comportement chaotique, et l’unique moyen de l’éliminer est d’exiger $`\epsilon c_+=\epsilon c_{}=0`$, c’est-à-dire $`\epsilon =0`$, restreignant ainsi à un seul arbitraire la partie analytique de la solution. Le problème est donc de trouver une expression compacte pour la série de Laurent (6).

## 3 Sur la suppression globale de la contribution chaotique

Par élimination de la constante mobile $`x_0`$ (donc de $`xx_0`$) entre la série de Laurent (6) et sa dérivée, la solution analytique générale satisfait à une EDO autonome d’ordre un,

$`F(u^{},u)=0,`$ (12)
que nous appellerons sous-équation de (4), puisque cette dernière en est une conséquence différentielle.

Trouver la solution analytique générale équivaut alors à trouver une sous-équation autonome d’ordre un. En effet, si l’on admet que $`F`$ est polynomiale et que sa solution générale est uniforme, alors

1. (Hermite \[32, Vol. II, §139\]) le genre de la courbe algébrique $`F=0`$ est zéro ou un,
2. (Briot et Bouquet \[24, pages 58–59\]) sa forme nécessaire est

$`F(u,u^{}){\displaystyle \underset{k=0}{\overset{m}{}}}{\displaystyle \underset{j=0}{\overset{2m2k}{}}}a_{j,k}u^ju_{}^{}{}_{}{}^{k}=0,a_{0,m}=1,`$ (13)
$`m`$ est un entier positif et les $`a_{j,k}`$ sont des constantes,
3. (Briot et Bouquet) la solution est, pour le genre zéro, une fraction rationnelle de $`e^{ax}`$ ou de $`x`$ ($`a`$ désignant une constante); pour le genre un, une fonction elliptique.

Dans les deux cas, il existe un algorithme, implanté en Maple par Mark van Hoeij , qui donne l’expression explicite de la solution.

Malheureusement, à notre connaissance, on ne sait pas effectuer l’élimination de $`xx_0`$ entre la série de Laurent et sa dérivée qui, si elle était possible, fournirait (12).

## 4 Rappel sur la classe d’équations $`P(u^{(n)},u)=0`$

$`P`$ désignant un polynôme de deux variables à coefficients complexes, la classe d’équations autonomes

$`P(u^{(n)},u)=0,`$ (14)
indépendantes des dérivées d’ordre $`2`$ à $`n1`$ inclus, a été classiquement beaucoup étudiée. Soit $`W`$ la classe de fonctions d’une variable complexe $`x`$

$`W:=\{\text{rationnel}(x),\text{rationnel}(e^{ax}),\text{elliptique}(x)\},`$ (15)
$`a`$ est une constante complexe. C’est en fait la classe des fonctions elliptiques et de leurs dégénérescences, car

$`\underset{a0}{lim}\text{rationnel}(e^{ax})=\underset{a0}{lim}\text{rationnel}((e^{ax}1)/a)=\text{rationnel}(x).`$ (16)
Les résultats connus sont les suivants.

1. (Briot et Bouquet) Pour $`n=1`$, si la solution est uniforme, alors elle est dans $`W`$.
2. (Picard) Pour $`n=2`$, toute solution méromorphe est dans $`W`$.
3. (Eremenko ) Si $`n`$ est impair, si le genre de la courbe algébrique $`P=0`$ est nul, et s’il existe une solution particulière méromorphe qui possède au moins un pôle, alors cette solution est dans $`W`$.

Ce troisième cas (qui se démontre par la théorie de Nevanlinna) s’applique à (4) pour $`b=\mu =0`$, la solution est alors connue

$`b=0,\mu =0:u=60\nu \mathrm{}^{}(xx_0,0,g_3),g_2=0,g_3={\displaystyle \frac{A}{1080\nu ^2}},`$ (17)
$`\mathrm{}`$ représentant la fonction elliptique de Weierstrass, et il en existe même une extrapolation de codimension un

$`b^216\mu \nu =0:u=60\nu \mathrm{}^{}15b\mathrm{}{\displaystyle \frac{b\mu }{4\nu }},g_2={\displaystyle \frac{\mu ^2}{12\nu ^2}},g_3={\displaystyle \frac{13\mu ^3+\nu A}{1080\nu ^3}}.`$ (18)
Notre problème revient à extrapoler cette solution pour une valeur arbitraire de $`b^2/(\mu \nu )`$.

## 5 Solutions de codimension non nulle

Elles appartiennent toutes (du moins celles qui sont connues) à la classe $`W`$, Eq. (15). L’unique solution connue de codimension trois est rationnelle,

$`b=0,\mu =0,A=0u=120\nu (xx_0)^3.`$ (19)
On connaît six solutions de codimension deux , toutes rationnelles en $`e^{kx}`$,

$`u`$ $`=`$ $`{\displaystyle \frac{5}{2}}bk^2{\displaystyle \frac{13b^3}{32\times 19\nu ^2}}+{\displaystyle \frac{7\mu b}{4\times 19\nu }}+\left({\displaystyle \frac{60}{19}}\mu 30\nu k^2{\displaystyle \frac{15b^2}{4\times 19\nu }}\right){\displaystyle \frac{k}{2}}\mathrm{tanh}{\displaystyle \frac{k}{2}}(xx_0)`$ (20)
$`15b\left({\displaystyle \frac{k}{2}}\mathrm{tanh}{\displaystyle \frac{k}{2}}(xx_0)\right)^2+120\nu \left({\displaystyle \frac{k}{2}}\mathrm{tanh}{\displaystyle \frac{k}{2}}(xx_0)\right)^3,`$
les valeurs permises figurant dans la Table 1.

Enfin, l’unique solution connue de codimension un est elliptique, cf. (18), et c’est une extrapolation de la dernière ligne de la Table 1.

Une belle propriété commune à toutes ces solutions est d’admettre la représentation

$`u`$ $`=`$ $`𝒟Log\sigma ,𝒟=60\nu {\displaystyle \frac{\text{d}}{\text{d}x^3}}+15b{\displaystyle \frac{\text{d}}{\text{d}x^2}}+{\displaystyle \frac{15(16\mu \nu b^2)}{76\nu }}{\displaystyle \frac{\text{d}}{\text{d}x}},`$ (21)
$`\sigma `$ est une fonction entière dont l’EDO est facile à construire.

Avant de rechercher plus avant la solution inconnue de codimension zéro, il convient de se poser la question si elle est ou non dans la classe $`W`$, c.à.d. si elle est elliptique.

On pourrait craindre en effet que cette classe $`W`$ soit insuffisante, au vu d’une EDO qui ne diffère de (4) que par le terme $`uu^{}`$ ,

$`Eu^{\prime \prime \prime }+9ku^{\prime \prime }+(a_1u26k^2)u^{}2a_1ku^2+24k^2u,`$ (22)
et dont la solution sort de la classe $`W`$,

$`u={\displaystyle \frac{12}{a_1}}e^{2kx}\mathrm{}({\displaystyle \frac{e^{kx}1}{k}}X_0,g_2,g_3),(X_0,g_2,g_3)\text{ arbitraires}.`$ (23)
Une solution dans cette nouvelle classe est exclue pour (4), car l’EDO d’ordre un résultant de l’élimination de $`X_0`$ entre (23) et sa dérivée,

$`(u^{}2ku)^24u^3+e^{4kx}g_2u+e^{6kx}g_3=0,`$ (24)
n’est autonome que pour la dégénérescence rationnelle $`g_2=g_3=0`$.

## 6 Vers une approche numérique de la solution elliptique

À partir de l’unique élément d’information que constitue la série de Laurent (6), nous allons définir plusieurs approches numériques pour tester la validité de l’hypothèse elliptique (appartenance à la classe $`W`$). L’outil naturel est celui des approximants de Padé (AP) .

Étant donné les $`N`$ premiers termes d’une série de Taylor en $`x=0`$,

$`S_N={\displaystyle \underset{j=0}{\overset{N}{}}}c_jx^j,`$ (25)
l’approximant de Padé $`[L,M]`$ de la série est défini comme l’unique fraction rationnelle

$`[L,M]={\displaystyle \frac{_{l=0}^La_lx^l}{_{m=0}^Mb_mx^m}},b_0=1,`$ (26)
satisfaisant à la condition

$`S_N[L,M]=𝒪(x^{N+1}),L+M=N.`$ (27)
L’extension aux séries de Laurent ne présente pas de difficulté. En particulier, les AP sont exacts sur les fractions rationnelles dès que $`L`$ et $`M`$ sont suffisamment grands.

Le calcul préalable d’un grand nombre (environ $`200`$) de termes de la série de Laurent se fait par un algorithme linéaire, la seule difficulté, aisée à surmonter, étant de ne jamais engendrer de termes qui s’avéreront inutiles. Il importe également de calculer les $`c_j`$, puis les $`(a_l,b_m)`$, sur $`𝒬`$, jamais sur $``$, à cause d’importantes compensations. C’est pour cela que nous avons dû écarter le programme pade de Numerical recipes , qui engendre des doublets de Froissart  indésirables (couples d’un zéro et d’un pôle très proches, ces couples se répartissant sur un cercle sans définir de coupure comme ce serait le cas s’il ne s’agissait pas d’erreurs numériques). La fonction pade du langage formel Maple  est en revanche tout à fait adaptée.

Remarque. Il est possible d’éviter l’étape de calcul des AP et d’obtenir directement et rapidement  les valeurs numériques des zéros et des pôles du Padé $`[L,M]`$.

### 6.1 Premier test d’ellipticité

La structure des pôles et des zéros des AP de (6) caractérise la classe de fonctions à laquelle appartient sa somme inconnue.

Une somme elliptique est caractérisée par une double infinité de pôles et de zéros situés aux nœuds d’un réseau doublement périodique, une somme rationnelle en $`e^{ax}`$, avec $`a`$ complexe, par une simple infinité de pôles et de zéros, une somme rationnelle en $`x`$ par un nombre fini de pôles et de zéros, etc. C’est pour tester les structures doublement périodiques que l’on a besoin d’un grand nombre de termes dans la série de Laurent.

Pour tous les cas de solution connue ((20), (18)), en choisissant pour les constantes fixes $`(\nu ,b,\mu ,A)`$ des valeurs rationnelles, le résultat est conforme aux espérances, voir Fig. 1. Pour plusieurs cas génériques de valeurs de $`b^2/(\mu \nu )`$ (par exemple $`b^2/(\mu \nu )=4,8`$, Fig. 2), nous trouvons une solide indication numérique que les textures observées sont toutes doublement périodiques, ce qui est la signature des fonctions elliptiques. Dans tous les cas, connus et inconnus, il est aisé de compter le nombre (commun) de zéros et de pôles par période, et les calculs indiquent que ce nombre serait égal à trois.

Remarque. La connaissance du nombre et des positions des zéros et des pôles d’une période permet la détermination numérique des invariants $`g_2`$ et $`g_3`$.

### 6.2 Deuxième test d’ellipticité

Si la solution $`u`$ est elliptique, alors, d’après un résultat classique, il existe deux fractions rationnelles $`R_1,R_2`$ telles que $`u=R_1(\mathrm{})+R_2(\mathrm{})\mathrm{}^{}`$, plus précisément

$$u(x,g_2,g_3)=\underset{\text{partie paire}}{\underset{}{\frac{\text{Poly}_{N_1}(\mathrm{})}{\text{Poly}_D(\mathrm{})}}}+\frac{\text{Poly}_{N_2}(\mathrm{})}{\text{Poly}_D(\mathrm{})}\sqrt{4\mathrm{}^3g_2\mathrm{}g_3},$$
(28)
$`\text{Poly}_D`$ désigne un polynôme de degré $`D`$.

En vue d’une recherche directe de la solution sous la forme (28), il est en principe possible de calculer numériquement les degrés $`N_1,N_2,D`$ des fractions rationnelles de $`\mathrm{}`$, en remplaçant le plan complexe de $`x`$ par celui de $`\mathrm{}(x)`$, de la manière suivante.

Tout d’abord, la nécessité d’une seule variable complexe dans les AP oblige à éliminer $`\mathrm{}^{}`$, donc à introduire le branchement algébrique indésirable $`\sqrt{4\mathrm{}^3g_2\mathrm{}g_3}`$. Un moyen simple de le supprimer est alors, puisque $`\mathrm{}`$ est une fonction paire, de se restreindre à la partie paire de la série de Laurent $`u(\chi =xx_0)`$,

$`u^{\mathrm{paire}}=15b\chi ^2+{\displaystyle \frac{b(56\mu \nu 13b^2)}{19\times 32\nu ^2}}{\displaystyle \frac{b\left(10\mu \nu 3b^2\right)^2}{19^2\times 64\nu ^4}}\chi ^2+𝒪(\chi ^4),`$ (29)
d’inverser \[14, vol. II, p. 527\] la série de Laurent $`\mathrm{}`$ de $`x`$,

$`\mathrm{}(x,g_2,g_3)=x^2+{\displaystyle \frac{g_2}{20}}x^2+{\displaystyle \frac{g_3}{28}}x^4+𝒪(x^6)`$ (30)
en la série de Laurent $`x^2`$ de $`\mathrm{}`$,

$`x^2`$ $`=`$ $`\mathrm{}^1[1+{\displaystyle \frac{g_2}{20}}\mathrm{}^2+{\displaystyle \frac{g_3}{28}}\mathrm{}^3+{\displaystyle \frac{7g_2^2}{1200}}\mathrm{}^4+{\displaystyle \frac{29g_2g_3}{3080}}\mathrm{}^5`$ (31)
$`+\left({\displaystyle \frac{11g_2^3}{12480}}+{\displaystyle \frac{5g_3^2}{1274}}\right)\mathrm{}^6+{\displaystyle \frac{167g_2^2g_3}{73920}}\mathrm{}^7`$
$`+({\displaystyle \frac{77g_2^4}{509184}}+{\displaystyle \frac{669g_2g_3^2}{340340}})\mathrm{}^8+𝒪(\mathrm{}^9)],`$
puis de composer les deux séries (29) et (31) pour obtenir

$`u^{\mathrm{paire}}`$ $`=`$ $`15b\mathrm{}+{\displaystyle \frac{b(56\mu \nu 13b^2)}{19\times 32\nu ^2}}+\left({\displaystyle \frac{3bg_2}{4}}{\displaystyle \frac{b\left(10\mu \nu 3b^2\right)^2}{19^2\times 64\nu ^4}}\right)\mathrm{}^1`$ (32)
$`+𝒪(\mathrm{}^2).`$
La question est alors de tester si cette série (32), dont les coefficients dépendent de $`(\nu ,b,\mu ,A,g_2,g_3)`$, a bien pour somme une fraction rationnelle

$$u^{\mathrm{paire}}(x,g_2,g_3)=\frac{\text{Poly}_{N_1}(\mathrm{})}{\text{Poly}_D(\mathrm{})}.$$
(33)
La difficulté provient de la nécessité de donner, dans la série (32), des valeurs numériques, de préférence rationnelles, à $`(g_2,g_3)`$ avant de pouvoir calculer les zéros et les pôles des AP. Puisque le parallélogramme période du plan $`x`$ est déjà déterminé, il permet de calculer des valeurs approchées (complexes) de $`(g_2,g_3)`$. Il importe donc d’étudier la sensibilité du résultat (valeurs de $`N_1,D`$) aux valeurs de $`(g_2,g_3)`$, et cette sensibilité ne peut être testée que sur l’unique solution elliptique connue (18). La Figure 3 montre la structure des AP de (32) dans le plan complexe $`\mathrm{}`$ pour $`b^2=16\mu \nu `$ et pour un choix de $`(g_2,g_3)`$ voisin de celui de ce cas.

### 6.3 Détermination des degrés de la sous-équation

Suivant une suggestion de Belabas faite lors de la conférence, les degrés de la courbe algébrique (13) pourraient se calculer par une technique de Padé, comme suit. À partir de la série de Laurent (6), on construit d’abord une autre série de Laurent à pôle simple, par exemple $`u^{}/u=\text{Laurent}(xx_0)`$, on l’inverse en $`xx_0=\text{Laurent}(u^{}/u)`$, et on la substitue dans la série $`u=\text{Laurent}(xx_0)`$. La série résultante $`u=\text{Laurent}(u^{}/u)`$, qui réalise sous forme locale l’élimination de $`xx_0`$ entre $`u`$ et $`u^{}`$, est alors analysée par Padé. Du branchement algébrique ainsi trouvé, il devrait être possible de déduire le degré en $`U`$ de la courbe algébrique $`G(U,V)=0`$, avec $`U=u,V=u^{}/u`$, donc les degrés de (13).

## 7 Approches analytiques

L’étude numérique précédente semble conclure à la nature elliptique de la solution cherchée. Nous présentons maintenant brièvement trois méthodes analytiques susceptibles de fournir les coefficients de la sous-équation (12). Si de plus celle-ci est supposée algébrique et à solution générale uniforme, alors elle a la forme nécessaire (13).

### 7.1 Équation de Monge-Ampère équivalente

La fonction inconnue de deux variables $`F(u,u^{})`$ qui définit la sous-équation d’ordre un (12) satisfait à une équation aux dérivées partielles très simple, obtenue par élimination de $`u^{\prime \prime }`$ et $`u^{\prime \prime \prime }`$ entre (4) et les deux équations

$`{\displaystyle \frac{\text{d}}{\text{d}x}}F(u(x),u^{}(x))=0,`$ (34)
$`{\displaystyle \frac{\text{d}^2}{\text{d}x^2}}F(u(x),u^{}(x))=0.`$ (35)
Dans la notation classique

$`X=u,Y=u^{},P={\displaystyle \frac{F}{X}},Q={\displaystyle \frac{F}{Y}},`$
$`R={\displaystyle \frac{^2F}{X^2}},S={\displaystyle \frac{^2F}{XY}},T={\displaystyle \frac{^2F}{Y^2}},`$ (36)
c’est une équation de Monge-Ampère qui s’écrit

$`\nu {\displaystyle \frac{Y^2(Q^2R+2PQSP^2T)+YP^2Q}{Q^3}}bY{\displaystyle \frac{P}{Q}}+\mu Y+{\displaystyle \frac{X^2}{2}}+A=0.`$ (37)
Notre problème équivaut alors à en trouver une solution particulière à zéro paramètre qui soit polynomiale. Ce nouveau problème pourrait n’être pas plus simple que l’original, car la méthode des caractéristiques est inapplicable puisque l’équation (37) hérite de la nature chaotique de KS. En particulier nous n’avons pas réussi à obtenir les degrés de cette solution polynomiale, si elle existe.

### 7.2 Construction de la sous-équation par l’algorithme LLL

Il existe une méthode , héritée de la cryptographie, qui fournit en nombres entiers les coefficients de (13) pour $`m`$ donné (et $`a_{0,m}`$ non normalisé à un). Connue sous le nom d’algorithme de réduction du réseau (lattice reduction algorithm) ou d’algorithme LLL, elle reçoit en entrée l’entier $`m`$, la suite des inconnues $`a_{j,k}`$ de (13), et $`L`$ valeurs numériques $`(u_i,u_i^{}),i=1,\mathrm{},L`$ (donc dans $``$ et non pas dans $`𝒬`$) de la somme des séries de Laurent $`u`$ et $`u^{}`$ aux points $`x_ix_0,i=1,\mathrm{},L`$. En sortie, si $`L`$ est suffisamment grand devant le nombre $`(m+1)^2`$ d’inconnues $`a_{j,k}`$, elle exhibe en un temps polynomial une suite d’entiers $`a_{j,k}`$, par résolution du système linéaire homogène surdéterminé

$`l=1,\mathrm{},L:{\displaystyle \underset{k=0}{\overset{m}{}}}{\displaystyle \underset{j=0}{\overset{2m2k}{}}}a_{j,k}u_l^ju_{l}^{}{}_{}{}^{k}=0.`$ (38)
Si les entiers ainsi trouvés sont tous nuls, il faut augmenter $`m`$ et recommencer. Nous n’avons pas encore mis cette méthode en œuvre car celle de la section suivante semble beaucoup plus prometteuse.

### 7.3 Construction de la sous-équation par la série de Laurent

Deux d’entre nous ont récemment proposé  une méthode analytique permettant de construire la sous-équation (13) à partir de la série de Laurent (6). L’algorithme est le suivant.

1. Choisir un entier positif $`m`$ et définir l’EDO de Briot and Bouquet (13). Elle contient $`(m+1)^21`$ constantes inconnues $`a_{j,k}`$.
2. Calculer $`J`$ termes de la série de Laurent, où $`J`$ est légèrement supérieur à $`(m+1)^21`$,

$`u=\chi ^p\left({\displaystyle \underset{j=0}{\overset{J}{}}}u_j\chi ^j+𝒪(\chi ^{j+1})\right),`$ (39)
avec dans le cas de KS $`p=3`$.
3. Exiger que la série de Laurent satisfasse à l’EDO de Briot and Bouquet, c.à.d. demander l’annulation identique de la série de Laurent de $`F(u,u^{})`$ à l’ordre $`J`$

$`F\chi ^D\left({\displaystyle \underset{j=0}{\overset{J}{}}}F_j\chi ^j+𝒪(\chi ^{J+1})\right),D=m(p1),j:F_j=0.`$ (40)
S’il n’existe pas de solution pour $`a_{j,k}`$, augmenter $`m`$ et revenir à l’étape 1.
4. Pour toute solution, intégrer l’EDO autonome d’ordre un (13).

Le cœur de la méthode est la troisième étape, où le système de $`J+1`$ équations $`F_j=0`$ aux $`(m+1)^21`$ inconnues $`a_{j,k}`$ est linéaire and surdéterminé, donc très facile à résoudre.

Quant à la quatrième étape, elle est résolue par le progiciel algcurves , qui met en œuvre des algorithmes de Poincaré.

Dans le cas de l’équation (4) étudiée ici, du fait de la valeur $`p=3`$, la valeur minimale de $`m`$ est $`m=3`$, ce qui correspond à $`15`$ inconnues $`a_{j,k}`$ et même à $`10`$ inconnues seulement, à cause de la règle de sélection supplémentaire $`m(p1)jp+k(p1)`$, qui provient du degré de singularité $`p=3`$ de $`u`$. Ces dix coefficients sont déterminés par le système de Cramer défini par les dix équations $`F_j=0,j=0:6,8,9,12`$. Le système résiduel, qui est non-linéaire et surdéterminé aux inconnues fixes $`(\nu ,b,\mu ,A)`$, a pour PGCD $`b^216\mu \nu `$ (cf. (18)), ce qui définit la solution de codimension un,

$`{\displaystyle \frac{b^2}{\mu \nu }}=16,\left(u^{}+{\displaystyle \frac{b}{2\nu }}u_s\right)^2\left(u^{}{\displaystyle \frac{b}{4\nu }}u_s\right)+{\displaystyle \frac{9}{40\nu }}\left(u_s^2+{\displaystyle \frac{60\mu ^3}{\nu }}+{\displaystyle \frac{10A}{3}}\right)^2=0,`$
$`u_s=u+{\displaystyle \frac{3b\mu }{2\nu }}`$ (41)
Après division par ce PGCD, le système restant aux inconnues $`(\nu ,b,\mu ,A)`$ admet exactement quatre solutions (il est suffisant d’arrêter la série à $`j=16`$ pour obtenir ce résultat),

$`b=0,`$
$`\left(u^{}+{\displaystyle \frac{180\mu ^2}{19^2\nu }}\right)^2\left(u^{}{\displaystyle \frac{360\mu ^2}{19^2\nu }}\right)+{\displaystyle \frac{9}{40\nu }}\left(u^2+{\displaystyle \frac{30\mu }{19}}u^{}{\displaystyle \frac{30^2\mu ^3}{19^2\nu }}\right)^2=0,`$ (42)
$`b=0,u_{}^{}{}_{}{}^{3}+{\displaystyle \frac{9}{40\nu }}\left(u^2+{\displaystyle \frac{30\mu }{19}}u^{}+{\displaystyle \frac{30^2\mu ^3}{19^3\nu }}\right)^2=0,`$ (43)
$`{\displaystyle \frac{b^2}{\mu \nu }}={\displaystyle \frac{144}{47}},u_s=u{\displaystyle \frac{5b\mu }{47\nu }},\left(u^{}+{\displaystyle \frac{b}{4\nu }}u_s\right)^3+{\displaystyle \frac{9}{40\nu }}u_s^4=0,`$ (44)
$`{\displaystyle \frac{b^2}{\mu \nu }}={\displaystyle \frac{256}{73}},u_s=u{\displaystyle \frac{45b\mu }{584\nu }},`$
$`\left(u^{}+{\displaystyle \frac{b}{8\nu }}u_s\right)^2\left(u^{}+{\displaystyle \frac{b}{2\nu }}u_s\right)+{\displaystyle \frac{9}{40\nu }}\left(u_s^2+{\displaystyle \frac{5b^3}{1024\nu ^2}}u_s+{\displaystyle \frac{5b^2}{128\nu }}u^{}\right)^2=0.`$ (45)
À l’étape 4, on trouve alors la valeur un pour le genre de (41), et la valeur zéro pour celui de (42)–(45), puis les expressions (18) et (20) pour les solutions.

Avec la valeur minimale $`m=3`$, on retrouve donc tous les résultats connus. Le calcul pour $`m4`$ est actuellement en cours.

## 8 Conclusion

L’étude par approximants de Padé semble conforter la nature probablement elliptique de l’unique solution analytique de codimension nulle de l’EDO chaotique (4). Parmi les autres pistes, mi-analytiques, mi-numériques, qu’il serait possible de mettre en œuvre, il en existe une particulièrement intéressante car constructive, c’est celle d’obtention de la sous-équation par la série de Laurent. Son seul inconvénient est l’absence de borne supérieure pour le degré $`m`$, ce qui contraint à de volumineux calculs qui sont encore en cours.

## Remerciements

Nous remercions les organisateurs et l’IRMA pour l’invitation à donner cette conférence, Karim Belabas pour plusieurs suggestions fructueuses, ainsi que Jean-Robert Burgan, Chiang Yik-man, Pierre Moussa, Jean-Pierre Ramis, Simon Ruijnesaars et Jean Thomann pour des discussions enrichissantes. La facilité de programmation du langage de calcul formel AMP  a constitué une aide précieuse.

T.-L. Yee remercie la Fondation Croucher de Hong-Kong, et MM le CEA pour leur soutien financier.


text:
# Continuous variable entanglement by radiation pressure

## 1 Introduction

The radiation pressure acting on a movable mirror realizes an optomechanical coupling between the incident optical modes and the various vibrational modes of the mirror. The use of this coupling has been proposed many years ago for the implementation of quantum limited measurements of mechanical forces , as the interferometric detection of gravitational-waves , or atomic force microscopy . Then optomechanical coupling has been proposed for quantum state engineering: for example, it has been shown that it may lead to nonclassical states of both the radiation field , and the motional degree of freedom of the mirror .

The appearance of quantum effects in ponderomotive systems, paves the way for using them also for quantum information purposes . In particular, quantum information can be encoded in the continuous quadratures of electromagnetic modes , and many quantum communication protocols can be implemented using entangled states of optical fields . A typical example of continuous variable entanglement is represented by the two-mode squeezed states at the output of a nondegenerate parametric amplifier . However, it has been showed in Refs.  that radiation pressure can also be used to entangle two or more modes of a cavity with a movable and perfectly reflecting mirror. In the present paper, we get rid of the cavity and we show that radiation pressure provides a new way to entangle travelling electromagnetic modes. We shall consider the very simple case of an intense optical mode, incident on a single, vibrating, and perfectly reflecting mirror. A vibrational mode of the mirror will induce sideband modes on the reflected field, which will be shown to be entangled. Entanglement proves to be very robust with respect to the thermal noise acting on the mirror. In particular, for a well defined interaction time, i.e., a given time duration of the incident field, the two first sideband modes become an EPR-correlated pair, identical to the two-mode squeezed state, with the ratio between the optical and the mechanical frequency playing the role of the squeezing parameter.

The paper is organized as follows. In Section II we shall derive the effective Hamiltonian of the system. In Section III we shall exactly solve the dynamics and provide a complete characterization of the reduced state of the system composed by the two reflected sideband modes. Section IV is for concluding remarks.

## 2 Effective Hamiltonian of the system

Let us consider a perfectly reflecting mirror and an intense laser beam impinging on its surface (see Fig. 1). For simplicity we consider only the motion and the elastic deformations of the mirror taking place along the spatial direction $`x`$, orthogonal to its reflecting surface.

The electromagnetic field exerts a force on the mirror proportional to its intensity and, at the same time, it is phase-shifted by the mirror displacement from the equilibrium position . In the limit of small mirror displacements, and in the interaction picture with respect to the free Hamiltonian of the electromagnetic field and the mirror displacement field $`\widehat{x}(𝐫,t)`$ ($`𝐫`$ is the coordinate on the mirror surface), one has the following Hamiltonian

$`\widehat{H}`$ $`=`$ $`{\displaystyle d^2𝐫\widehat{P}(𝐫,t)\widehat{x}(𝐫,t)},`$ (1)
where $`\widehat{P}(𝐫,t)`$ is the radiation pressure force . All the continuum of electromagnetic modes with positive longitudinal wave vector $`q`$ and transverse wave vector $`𝐤`$ contributes to the radiation pressure force. Following Ref., and considering linearly polarized radiation with the electric field parallel to the mirror surface, we have

$`\widehat{P}(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{8\pi ^3}}{\displaystyle 𝑑𝐤𝑑q𝑑𝐤^{}𝑑q^{}\frac{c^2qq^{}}{\sqrt{\omega \omega ^{}}}(𝐮_k𝐮_k^{})𝐮_q}`$ (2)
$`\times `$ $`\{\widehat{a}(𝐤,q)\widehat{a}(𝐤^{},q^{})\mathrm{exp}[i(\omega +\omega ^{})t+i(𝐤+𝐤^{})𝐫]`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}^{}(𝐤^{},q^{})\mathrm{exp}[i(\omega +\omega ^{})ti(𝐤+𝐤^{})𝐫]`$
$`+\widehat{a}(𝐤,q)\widehat{a}^{}(𝐤^{},q^{})\mathrm{exp}[i(\omega \omega ^{})t+i(𝐤𝐤^{})𝐫]`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}(𝐤^{},q^{})\mathrm{exp}[i(\omega \omega ^{})ti(𝐤𝐤^{})𝐫]\},`$
where $`\widehat{a}(𝐤,q)`$ are continuous mode destruction operators obeying the commutation relations

$$[\widehat{a}(𝐤,q),\widehat{a}^{}(𝐤^{},q^{})]=\delta (𝐤𝐤^{})\delta (qq^{}).$$
(3)
Furthermore, the electromagnetic wave frequencies $`\omega `$ and $`\omega ^{}`$ are given by $`\omega =\sqrt{c^2(k^2+q^2)}`$ and $`\omega ^{}=\sqrt{c^2(k^2+q^2)}`$ ($`c`$ is the light speed in vacuum), and $`𝐮_k`$, $`𝐮_q`$ denote dimensionless unit vectors parallel to $`𝐤`$, $`q`$ respectively.

The mirror displacement $`\widehat{x}(𝐫,t)`$ is generally given by a superposition of many acoustic modes ; however, a single vibrational mode description can be adopted whenever detection is limited to a frequency bandwidth including a single mechanical resonance. In particular, focused light beams are able to excite Gaussian acoustic modes, in which only a small portion of the mirror, localized at its center, vibrates. These modes have a small waist $`w`$, a large mechanical quality factor $`Q`$, a small effective mass $`M`$ , and the simplest choice is to choose the fundamental Gaussian mode with frequency $`\mathrm{\Omega }`$ and annihilation operator $`\widehat{b}`$, with $`[\widehat{b},\widehat{b}^{}]=1`$,

$$\widehat{x}(𝐫,t)=\sqrt{\frac{\mathrm{}}{2M\mathrm{\Omega }}}\left[\widehat{b}e^{i\mathrm{\Omega }t}+\widehat{b}^{}e^{i\mathrm{\Omega }t}\right]\mathrm{exp}(r^2/w^2).$$
(4)
By inserting Eqs. (2) and (4) in Eq. (1) and integrating over the variable $`𝐫`$, one obtains

$`\widehat{H}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}w^2}{8\pi ^2}}\sqrt{{\displaystyle \frac{\mathrm{}}{2M\mathrm{\Omega }}}}{\displaystyle 𝑑𝐤𝑑q𝑑𝐤^{}𝑑q^{}\frac{c^2qq^{}}{\sqrt{\omega \omega ^{}}}(𝐮_k𝐮_k^{})𝐮_q}`$ (5)
$`\times `$ $`\{\widehat{a}(𝐤,q)\widehat{a}(𝐤^{},q^{})\mathrm{exp}[i(\omega +\omega ^{})t(𝐤+𝐤^{})^2w^2/4]`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}^{}(𝐤^{},q^{})\mathrm{exp}[i(\omega +\omega ^{})t(𝐤+𝐤^{})^2w^2/4]`$
$`+\widehat{a}(𝐤,q)\widehat{a}^{}(𝐤^{},q^{})\mathrm{exp}[i(\omega \omega ^{})t(𝐤𝐤^{})^2w^2/4]`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}(𝐤^{},q^{})\mathrm{exp}[i(\omega \omega ^{})t(𝐤𝐤^{})^2w^2/4]\}`$
$`\times \left\{\widehat{b}e^{i\mathrm{\Omega }t}+\widehat{b}^{}e^{i\mathrm{\Omega }t}\right\}.`$
In common situations, the acoustical waist $`w`$ is much larger than typical optical wavelengths , and therefore we can approximate $`\mathrm{exp}\left\{(𝐤\pm 𝐤^{})^2w^2/4\right\}w^2/4\pi \delta (𝐤\pm 𝐤^{})`$ and then integrate Eq. (5) over $`𝐤^{}`$, obtaining

$`\widehat{H}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2\pi }}\sqrt{{\displaystyle \frac{\mathrm{}}{2M\mathrm{\Omega }}}}{\displaystyle 𝑑𝐤𝑑q𝑑q^{}\frac{c^2qq^{}}{\sqrt{\omega \omega ^{}}}}`$ (6)
$`\times `$ $`\{\widehat{a}(𝐤,q)\widehat{a}(𝐤,q^{})\mathrm{exp}[i(\omega +\omega ^{})t]`$
$`\widehat{a}^{}(𝐤,q)\widehat{a}^{}(𝐤,q^{})\mathrm{exp}[i(\omega +\omega ^{})t]`$
$`+\widehat{a}(𝐤,q)\widehat{a}^{}(𝐤,q^{})\mathrm{exp}[i(\omega \omega ^{})t]`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}(𝐤,q^{})\mathrm{exp}[i(\omega \omega ^{})t]\}\times \{\widehat{b}e^{i\mathrm{\Omega }t}+\widehat{b}^{}e^{i\mathrm{\Omega }t}\}.`$
We now make the Rotating Wave Approximation (RWA), that is, we neglect all the terms oscillating in time faster than the mechanical frequency $`\mathrm{\Omega }`$. This means averaging the Hamiltonian over a time $`\tau `$ such that $`\mathrm{\Omega }\tau 1`$, yielding the following replacements in Eq. (6)

$$\mathrm{exp}\left\{\pm i(\omega ^{}\pm \omega \pm \mathrm{\Omega })t\right\}\frac{2\pi }{\tau }\delta (\omega ^{}\pm \omega \pm \mathrm{\Omega }).$$
(7)
The parameter $`\tau `$ is not arbitrary, but its inverse, $`1/\tau =\mathrm{\Delta }\nu _{det}`$, is the detection bandwidth, that is, the spectral resolution of the detection apparatus employed.

Since $`\omega `$ and $`\omega ^{}`$ are positive and $`\mathrm{\Omega }`$ is much smaller than typical optical frequencies, the two terms $`\delta (\omega ^{}+\omega \pm \mathrm{\Omega })`$ give no contribution, while the other two terms can be rewritten as

$$\frac{2\pi }{\tau }\delta (\omega ^{}\omega \pm \mathrm{\Omega })=2\pi \mathrm{\Delta }\nu _{det}\delta (q^{}\overline{q}_\pm )\frac{\omega ^{}(\overline{q}_\pm )}{c^2\overline{q}_\pm },$$
(8)
where $`\overline{q}_\pm =\sqrt{(\omega \pm \mathrm{\Omega })^2/c^2k^2}`$. Integrating over $`q^{}`$ we get

$`\widehat{H}`$ $`=`$ $`\mathrm{}\mathrm{\Delta }\nu _{det}\sqrt{{\displaystyle \frac{\mathrm{}}{2M\mathrm{\Omega }}}}{\displaystyle }d𝐤{\displaystyle }dq{\displaystyle \frac{q}{\sqrt{\omega }}}\{\widehat{a}(𝐤,q)\widehat{a}^{}(𝐤,\overline{q}_+)\widehat{b}\sqrt{\omega +\mathrm{\Omega }}`$ (9)
$`+\widehat{a}(𝐤,q)\widehat{a}^{}(𝐤,\overline{q}_{})\widehat{b}^{}\sqrt{\omega \mathrm{\Omega }}`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}(𝐤,\overline{q}_+)\widehat{b}^{}\sqrt{\omega +\mathrm{\Omega }}`$
$`+\widehat{a}^{}(𝐤,q)\widehat{a}(𝐤,\overline{q}_{})\widehat{b}\sqrt{\omega \mathrm{\Omega }}\},`$
where we have used the fact that $`\omega ^{}(\overline{q}_\pm )=\omega \pm \mathrm{\Omega }`$.

We now consider the situation where the radiation field incident on the mirror is characterized by an intense, quasi-monochromatic, laser field with transversal wave vector $`𝐤_\mathrm{𝟎}`$, longitudinal wave vector $`q_0`$, cross-sectional area $`A`$, and power $`\mathrm{}`$. Since this component is very intense, it can be treated as classical and one can approximate $`\widehat{a}(𝐤,q)\alpha (𝐤,q)`$ in Eq. (9), where (with an appropriate choice of phases)

$$\alpha (𝐤,q)=i\sqrt{\frac{(2\pi )^3\mathrm{}}{\mathrm{}\omega _0cA}}\delta (𝐤𝐤_\mathrm{𝟎})\delta (qq_0),$$
(10)
with $`\omega _0=c\sqrt{𝐤_{\mathrm{𝟎}}^{}{}_{}{}^{2}+q_0^2}`$.

Due to the Dirac delta, the only nonvanishing terms in the optomechanical interaction driven by the intense laser beam involve only two back-scattered waves, that is, the sidebands of the driving beam at frequencies $`\omega _0\pm \mathrm{\Omega }`$, as described by

$`\widehat{H}`$ $`=`$ $`i\mathrm{}\mathrm{\Delta }\nu _{det}\sqrt{{\displaystyle \frac{\mathrm{}}{2M\mathrm{\Omega }}}}q_0\sqrt{{\displaystyle \frac{\mathrm{}}{\mathrm{}\omega _0cA}}}`$ (11)
$`\times `$ $`\{\sqrt{{\displaystyle \frac{\omega _0+\mathrm{\Omega }}{\omega _0}}}\widehat{a}^{}(𝐤_\mathrm{𝟎},\overline{q}_+)\widehat{b}+\sqrt{{\displaystyle \frac{\omega _0\mathrm{\Omega }}{\omega _0}}}\widehat{a}^{}(𝐤_\mathrm{𝟎},\overline{q}_{})\widehat{b}^{}`$
$`\sqrt{{\displaystyle \frac{\omega _0+\mathrm{\Omega }}{\omega _0}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_+)\widehat{b}^{}\sqrt{{\displaystyle \frac{\omega _0\mathrm{\Omega }}{\omega _0}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_{})\widehat{b}\},`$
where now $`\overline{q}_\pm =\sqrt{(\omega _0\pm \mathrm{\Omega })^2/c^2k_0^2}`$. The physical process described by this interaction Hamiltonian is very similar to a stimulated Brillouin scattering , even though in this case the Stokes and anti-Stokes component are back-scattered by the acoustic waves at reflection, and the optomechanical coupling is provided by the radiation pressure and not by the dielectric properties of the mirror.

In practice, either the driving laser beam and the back-scattered modes are never monochromatic, but have a nonzero bandwidth. In general the bandwidth of the back-scattered modes is determined by the bandwidth of the driving laser beam and that of the acoustic mode. However, due to its high mechanical quality factor, the spectral width of the mechanical resonance is negligible (about $`1`$ Hz) and, in practice, the bandwidth of the two sideband modes $`\mathrm{\Delta }\nu _{mode}`$ coincides with that of the incident laser beam. It is then convenient to consider this nonzero bandwidth to redefine the bosonic operators of the Stokes and anti-Stokes modes to make them dimensionless,

$`\widehat{a}_1`$ $`=`$ $`2\pi \sqrt{{\displaystyle \frac{2\pi \mathrm{\Delta }\nu _{mode}}{cA}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_{})=2\pi \sqrt{{\displaystyle \frac{\mathrm{\Delta }q}{A}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_{})`$ (12)
$`\widehat{a}_2`$ $`=`$ $`2\pi \sqrt{{\displaystyle \frac{2\pi \mathrm{\Delta }\nu _{mode}}{cA}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_+)=2\pi \sqrt{{\displaystyle \frac{\mathrm{\Delta }q}{A}}}\widehat{a}(𝐤_\mathrm{𝟎},\overline{q}_+),`$ (13)
so that Eq.(11) reduces to an effective Hamiltonian

$$\widehat{H}_{eff}=i\mathrm{}\chi (\widehat{a}_1\widehat{b}\widehat{a}_1^{}\widehat{b}^{})i\mathrm{}\theta (\widehat{a}_2\widehat{b}^{}\widehat{a}_2^{}\widehat{b}),$$
(14)
where the couplings $`\chi `$ and $`\theta `$ are given by

$`\chi `$ $`=`$ $`q_0\mathrm{\Delta }\nu _{det}\sqrt{{\displaystyle \frac{\mathrm{}}{2M\mathrm{\Omega }}}}\sqrt{{\displaystyle \frac{\mathrm{}}{\mathrm{\Delta }\nu _{mode}\mathrm{}\omega _0}}}\sqrt{{\displaystyle \frac{\omega _0\mathrm{\Omega }}{\omega _0}}}`$ (15)
$`=`$ $`\mathrm{cos}\varphi _0\sqrt{{\displaystyle \frac{\mathrm{}\mathrm{\Delta }\nu _{det}^2(\omega _0\mathrm{\Omega })}{2M\mathrm{\Omega }c^2\mathrm{\Delta }\nu _{mode}}}}`$
$`\theta `$ $`=`$ $`\chi \sqrt{{\displaystyle \frac{\omega _0+\mathrm{\Omega }}{\omega _0\mathrm{\Omega }}}},`$ (16)
with $`\varphi _0=\mathrm{arccos}(cq_0/\omega _0)`$, is the angle of incidence of the driving beam. It is possible to verify that with the above definitions, the Stokes and anti-Stokes annihilation operators $`a_1`$ and $`a_2`$ satisfy the usual commutation relations $`[a_i,a_j^{}]=\delta _{i,j}`$.

## 3 System dynamics

Eq. (14) contains two interaction terms: the first one, between modes $`\widehat{a}_1`$ and $`\widehat{b}`$, is a parametric-type interaction leading to squeezing in phase space , and it is able to generate the EPR-like entangled state which has been used in the continuous variable teleportation experiment of Ref. . The second interaction term, between modes $`\widehat{a}_2`$ and $`\widehat{b}`$, is a beam-splitter-type interaction , which may degrade the entanglement between modes $`\widehat{a}_1`$ and $`\widehat{b}`$ generated by the first term. In general, one has a system of three bosonic modes, coupled by a bilinear interaction. The corresponding Hamiltonian evolution of the system can be straightforwardly obtained. This Hamiltonian description satisfactorily reproduces the dynamics as long as the dissipative coupling of the mirror vibrational mode with its environment is negligible. This happens in the case of modes with a high-Q mechanical quality factor. In this case, the mechanical frequency $`\mathrm{\Omega }`$ is sufficiently high (some MHz) so that the RWA of Eq. (7) can be made, and at the same time we can consider an interaction time, i.e., a time duration of the incident laser pulse, much smaller than the relaxation time of the vibrational mode (which can be of order of one second).

The system dynamics can be easily studied through the (normally ordered) characteristic function $`\mathrm{\Phi }(\mu ,\nu ,\zeta )`$, where $`\mu ,\nu ,\zeta `$ are the complex variables corresponding to the operators $`\widehat{a}_1,\widehat{b},\widehat{a}_2`$ respectively. From the Hamiltonian (14) the dynamical equation for $`\mathrm{\Phi }`$ results

$`\dot{\mathrm{\Phi }}`$ $`=`$ $`\chi \left(\mu \nu +\mu ^{}\nu ^{}\mu ^{}{\displaystyle \frac{}{\nu }}\mu {\displaystyle \frac{}{\nu ^{}}}\nu ^{}{\displaystyle \frac{}{\mu }}\nu {\displaystyle \frac{}{\mu ^{}}}\right)\mathrm{\Phi }`$ (17)
$`+\theta \left(\zeta ^{}{\displaystyle \frac{}{\nu ^{}}}+\zeta {\displaystyle \frac{}{\nu }}\nu ^{}{\displaystyle \frac{}{\zeta ^{}}}\nu {\displaystyle \frac{}{\zeta }}\right)\mathrm{\Phi },`$
with the initial condition

$$\mathrm{\Phi }(t=0)=\mathrm{exp}\left[\overline{n}|\nu |^2\right],$$
(18)
corresponding to the vacuum for modes $`\widehat{a}_1`$, $`\widehat{a}_2`$ and to a thermal state for mode $`\widehat{b}`$. The latter is characterized by an average number of excitations $`\overline{n}=[\mathrm{exp}(\mathrm{}\mathrm{\Omega }/k_BT)1]^1`$, $`T`$ being the equilibrium temperature and $`k_B`$ the Boltzmann constant. Since the initial condition is a Gaussian state of the three-mode system, and the system is linear, the joint state of the whole system at time $`t`$ is still Gaussian, with characteristic function

$`\mathrm{\Phi }`$ $`=`$ $`\mathrm{exp}[𝒜|\mu |^2|\nu |^2|\zeta |^2`$ (19)
$`+𝒞\mu \nu +𝒞\mu ^{}\nu ^{}+\mu \zeta +\mu ^{}\zeta ^{}+𝒟\nu \zeta ^{}+𝒟\nu ^{}\zeta ],`$
where

$`𝒜(t)`$ $`=`$ $`{\displaystyle \frac{1}{2(r^21)^2}}\{[1\mathrm{cos}\left(2\mathrm{\Theta }t\right)][\overline{n}(r^21)1]`$ (20)
$`+4r^2[1\mathrm{cos}\left(\mathrm{\Theta }t\right)]\},`$
$`(t)`$ $`=`$ $`{\displaystyle \frac{1}{2(r^21)}}\left\{\left[1+\mathrm{cos}\left(2\mathrm{\Theta }t\right)\right]\overline{n}(r^21)+1\mathrm{cos}\left(2\mathrm{\Theta }t\right)\right\},`$ (21)
$`𝒞(t)`$ $`=`$ $`{\displaystyle \frac{1}{2(r^21)^{3/2}}}\{2r^2\mathrm{sin}\left(\mathrm{\Theta }t\right)`$ (22)
$`+[\overline{n}(r^21)1]\mathrm{sin}\left(2\mathrm{\Theta }t\right)\},`$
$`𝒟(t)`$ $`=`$ $`{\displaystyle \frac{r}{2(r^21)^{3/2}}}\left\{\mathrm{sin}\left(\mathrm{\Theta }t\right)+\left[\overline{n}(r^21)1\right]\mathrm{sin}\left(2\mathrm{\Theta }t\right)\right\},`$ (23)
$`(t)`$ $`=`$ $`{\displaystyle \frac{r^2}{2(r^21)^2}}\{[1\mathrm{cos}\left(2\mathrm{\Theta }t\right)][\overline{n}(r^21)1]`$ (24)
$`+4[1\mathrm{cos}\left(\mathrm{\Theta }t\right)]\},`$
$`(t)`$ $`=`$ $`{\displaystyle \frac{r}{2(r^21)^2}}\{[1\mathrm{cos}\left(2\mathrm{\Theta }t\right)][\overline{n}(r^21)1]`$ (25)
$`+2(1+r^2)[1\mathrm{cos}\left(\mathrm{\Theta }t\right)]\},`$
$`r=\theta /\chi =\left[(\omega _0+\mathrm{\Omega })/(\omega _0\mathrm{\Omega })\right]^{1/2}`$ and $`\mathrm{\Theta }=\chi (r^21)^{1/2}`$. The corresponding density operator can be expressed as

$`\widehat{\rho }_{1b2}={\displaystyle \frac{d^2\mu }{\pi }\frac{d^2\nu }{\pi }\frac{d^2\zeta }{\pi }}`$ $`\mathrm{\Phi }(\mu ,\nu ,\zeta ,t)e^{|\mu |^2|\nu |^2|\zeta |^2}`$ (26)
$`\times \widehat{D}_1(\mu )\widehat{D}_b(\nu )\widehat{D}_2(\zeta ),`$
where $`\widehat{D}_i`$ indicates the corresponding normally ordered displacement operator .

Here we are interested in the reduced state of the system composed by the two reflected optical sideband modes. They do not interact directly but their interaction is mediated by the optomechanical coupling of each mode with the moving mirror. This reduced state can be immediately obtained by tracing over the mirror mode $`\widehat{b}`$, obtaining

$`\mathrm{\Phi }_{12}(\mu ,\zeta )=\mathrm{exp}\left[𝒜|\mu |^2|\zeta |^2+\left(\mu \zeta +\mu ^{}\zeta ^{}\right)\right].`$ (27)
Introducing the vector of field quadratures $`\widehat{𝐯}=(\widehat{X}_1,\widehat{P}_1,\widehat{X}_2,\widehat{P}_2)`$, where

$$\widehat{X}_j=\frac{\widehat{a}_j+\widehat{a}_j^{}}{\sqrt{2}},\widehat{P}_j=\frac{\widehat{a}_j\widehat{a}_j^{}}{i\sqrt{2}},j=1,2,$$
(28)
it is possible to connect the characteristic function (27) with the correlation matrix $`\mathrm{\Gamma }`$ of the two sideband modes, defined as $`\mathrm{\Gamma }_{i,j}=\widehat{𝐯}_i\widehat{𝐯}_j+\widehat{𝐯}_j\widehat{𝐯}_i/2\widehat{𝐯}_i\widehat{𝐯}_j`$, and which, in this case, is equal to

$$\mathrm{\Gamma }=\left(\begin{array}{cccc}𝒜+\frac{1}{2}& 0& & 0\\ 0& 𝒜+\frac{1}{2}& 0& \\ & 0& +\frac{1}{2}& 0\\ 0& & 0& +\frac{1}{2}\end{array}\right).$$
(29)
The Gaussian state of the two optical modes is completely characterized by this correlation matrix. In particular we can employ the necessary and sufficient criterion for entanglement in the case of Gaussian states derived by Simon , which allows us to determine the parameter region where the two modes are entangled by the optomechanical interaction with the mirror vibrational mode. Simon’s necessary and sufficient criterion for entanglement becomes, for the two-sideband modes state considered here,

$`𝒯`$ $``$ $`\left(𝒜+{\displaystyle \frac{1}{4}}+{\displaystyle \frac{𝒜+}{2}}^2\right)^2+{\displaystyle \frac{1}{16}}{\displaystyle \frac{^2}{2}}`$ (30)
$`\left({\displaystyle \frac{𝒜}{2}}+{\displaystyle \frac{1}{4}}\right)^2\left({\displaystyle \frac{}{2}}+{\displaystyle \frac{1}{4}}\right)^2<0`$
As it can be immediately seen from Eqs. (20)-(25), the dynamics of the system is determined only by three dimensionless parameters, $`\mathrm{\Theta }t`$, $`r`$, and $`\overline{n}`$, and it is periodic in $`\mathrm{\Theta }t`$ with period $`2\pi `$. The parameter $`r`$ is determined by the ratio between the mechanical and the driving laser frequency, $`\mathrm{\Omega }/\omega _0`$ and it is always very close to $`1`$. For this reason in Fig. 2 we show the marker of entanglement $`𝒯`$ as a function of the scaled time $`\mathrm{\Theta }t`$ for $`0\mathrm{\Theta }t2\pi `$, for different values of the initial mean thermal vibrational number $`\overline{n}`$ and at a fixed value of $`r`$, i.e., $`r1=2.5\times 10^7`$. Fig. 2 clearly shows that the two sideband modes are almost always entangled (obviously except for $`t=2k\pi `$, for integer $`k`$, due to the factorized initial condition) and that this entanglement is extremely robust with respect to the thermal noise acting on the mirror. In fact, the marker of entanglement $`𝒯`$ remains always negative, even at extremely large values of $`\overline{n}`$. Fig. 2 shows $`𝒯`$ for $`\overline{n}=0`$ (full line), $`\overline{n}=10^5`$ (dashed line, almost coinciding with that for $`\overline{n}=0`$), $`\overline{n}=5\times 10^6`$ (dotted-dashed line), and $`\overline{n}=10^7`$ (dotted line). This result is extremely interesting because it shows that radiation pressure proves to be an efficient source of entangled travelling optical fields, not particularly affected by thermal noise. In this treatment we have only considered this kind of noise and we have neglected other more technical noise sources, such as the fluctuations of the intensity and the frequency of the incident driving field. In fact, the latter, differently from thermal noise, are negligible in today experiments involving optomechanical systems .

From Fig. 2, it is evident that $`𝒯`$ at exactly half period, $`\mathrm{\Theta }t=\pi `$, is unaffected by thermal noise, i.e., it is independent of $`\overline{n}`$. In fact, using Eqs (20)-(25) and (30), it is easy to see that $`𝒯(\mathrm{\Theta }t=\pi )=4(r+r^3)^2/(r^21)^4`$ (actually $`𝒯`$ has been scaled exactly by the modulus of this value in Fig. 2). This independence of $`\overline{n}`$ is a manifestation of quantum interference between the dynamical effects of the two interaction terms in Eq. (14). It is therefore interesting to see the reduced state of the two reflected sideband modes just at this interaction time, $`\mathrm{\Theta }t=\pi `$. Using Eqs. (20)-(25), one has the following correlation matrix

$`\mathrm{\Gamma }(\mathrm{\Theta }t=\pi )=`$ (31)
$`=\left(\begin{array}{cccc}\frac{4r^2}{(r^21)^2}+\frac{1}{2}& 0& \frac{2r(1+r^2)}{(r^21)^2}& 0\\ 0& \frac{4r^2}{(r^21)^2}+\frac{1}{2}& 0& \frac{2r(1+r^2)}{(r^21)^2}\\ \frac{2r(1+r^2)}{(r^21)^2}& 0& \frac{4r^2}{(r^21)^2}+\frac{1}{2}& 0\\ 0& \frac{2r(1+r^2)}{(r^21)^2}& 0& \frac{4r^2}{(r^21)^2}+\frac{1}{2}\end{array}\right).`$ (36)
It is easy to check that this is just the correlation matrix of an EPR-like entangled state, i.e., identical to the two-mode squeezed state generated by a parametric amplifier . The correspondence between the squeezing parameter $`\xi `$ and our parameters is $`\mathrm{sinh}\xi 2r(1+r^2)/(r^21)^2`$, so that, in practice one can reach very high two-mode squeezing by choosing $`r10`$, i.e., by decreasing the ratio $`\mathrm{\Omega }/\omega _0`$. This fact can be equivalently seen by evaluating the variances of the linear combination of field quadratures, typical of the two-mode squeezed state, i.e.,

$`\mathrm{\Delta }_{}`$ $`=`$ $`(X_1X_2)^2=(P_1+P_2)^2,`$ (37)
$`\mathrm{\Delta }_+`$ $`=`$ $`(X_1+X_2)^2=(P_1P_2)^2.`$ (38)
For $`\mathrm{\Theta }t=\pi `$, one has

$`\mathrm{\Delta }_{}(\mathrm{\Theta }t=\pi )`$ $`=`$ $`\left({\displaystyle \frac{r1}{r+1}}\right)^2,`$ (39)
$`\mathrm{\Delta }_+(\mathrm{\Theta }t=\pi )`$ $`=`$ $`\left({\displaystyle \frac{r+1}{r1}}\right)^2.`$ (40)
The time behavior of the variances (37) and (38) is depicted in Fig. 3. It is worth noticing that the typical EPR correlations are available for a wide time range around $`\mathrm{\Theta }t=\pi `$. Furthermore, they are robust against thermal noise since the curves are shown for $`\overline{n}=10^5`$, and they practically coincide with those for $`\overline{n}=0`$. That is in agreement with the results shown in Fig.2.

## 4 Results and Conclusions

In conclusion, we have seen that, under appropriate conditions, the radiation pressure of an intense optical pulse impinging on a perfectly reflecting mirror couples efficiently the first two sideband modes induced by a vibrational mode of the mirror. This optomechanical interaction actually entangles the two modes, and quite surprisingly, the resulting entanglement is extremely robust with respect to the thermal noise acting on the mirror. In particular, due to quantum interference, if the interaction time, i.e., the duration of the incident driving pulse, is appropriately tailored, the reduced state of the reflected sideband modes is a two-mode squeezed state. Therefore, rather unexpectedly, radiation pressure could be employed as a new source of entangled two-mode squeezed state, in which the role of the squeezing parameter is played by the ratio between the optical and the mechanical frequency $`\omega _0/\mathrm{\Omega }`$. These interesting results hold if the time duration of the pulse satisfies two conditions: i) it has to be much longer than the period of mechanical oscillations, because we have to average over these oscillations; ii) it has to be shorter than the vibrational mode relaxation time (mechanical damping has been neglected in our treatment). These conditions are not extremely easy but they could be realized using micro-opto-electro-mechanical-systems (MOEMS) for example . Possible values are: damping rates $`\gamma _m1`$ Hz, $`\mathrm{}=10`$W, $`\omega _02\times 10^{15}`$ Hz, $`\mathrm{\Omega }5\times 10^8`$ Hz, $`\mathrm{\Delta }\nu _{det}10^7`$ Hz, $`\mathrm{\Delta }\nu _{mode}t^110^3`$ Hz, and $`M10^{10}`$ Kg, yielding $`\chi \theta 5\times 10^5`$ Hz, and $`\mathrm{\Theta }10^3`$ Hz. Therefore, using well tailored and long pulses (order of $`1`$ ms), one could realize this new source of two mode squeezing.

## References