## Abstract

In classical mechanics, the two-body problem has been well studied. The governing equations form a system of two-coupled second-order nonlinear differential equations for the radial and angular coordinates. The perturbation induced by the astronomical disturbance like ‘dust’ is normally not considered in the orbit dynamics. Distributed dust produces an additional random force on the orbiting particle, which can be modelled as a random force having ‘Gaussian statistics’. The estimation of accurate positioning of the orbiting particle is not possible without accounting for the stochastic perturbation felt by the orbiting particle. The objective of this paper is to use the stochastic differential equation (SDE) formalism to study the effect of such disturbances on the orbiting body. Specifically, in this paper, we linearize SDEs about the mean of the state vector. The linearization operation performed above, transforms the system of SDEs into another system of SDEs that resembles a bilinear system, as described in signal processing and control literature. However, the mean trajectory of the resulting bilinear stochastic differential model does not preserve the perturbation effect felt by the orbiting particle; only the variance trajectory includes the perturbation effect. For this reason, the effectiveness of the dust-perturbed model is examined on the basis of the bilinear and second-order approximations of the *system nonlinearity*. The bilinear and second-order approximations of the system nonlinearity allow substantial simplifications for the numerical implementation and preserve some of the properties of the original stochastically perturbed model. Most notably, this paper reveals that the Brownian motion process is accurate to model and study the effect of dust perturbation on the orbiting particle. In addition, analytical findings are supported with finite difference method-based numerical simulations.

## 1. Introduction

Linear models for physical processes come under two broad categories: time varying and time-invariant. In application, where forces depend both on the input and the output processes, it is more reasonable to employ a bilinear approximation rather than a linear approximation because the linear approximation will neglect the input–output coupling terms. Bilinear models appear in areas like channel equalization, echo cancellation, nonlinear tracking, multiplicative disturbance tracking as well as in many real processes in the fields of socio-economics, ecology, agriculture, biology, population dynamics and industry, etc. (Bruni *et al*. 1974; Mohler & Kolodziej 1980; Niedzwiecki & Cisowski 1996). Many other papers have been published on the control and filtering problems for the bilinear systems (Mohler 1973; Terdik1990; Elliott 1998; Carravetta *et al*. 2000).

Scheeres *et al*. (2001*a*,*b*) have published an interesting paper, which talks about solutions to the dynamics of an orbiter, which moves in the gravitational field of a planet and a satellite of the planet. In their paper, dynamics of the planetary satellite orbiter were cast in the form of a restricted three-body problem. Another problem involving perturbation to the two-body problem (planet–Sun system) is given in the classic book by Goldstein *et al*. (2002, pp.129–130), where they consider an additional deterministic force on a planet produced by a uniformly distributed dust between a planet and the Sun. A discussion that it is unrealistic to assume a uniform distribution of dust between a planet and the Sun is provided by Mann *et al*. (2003). Bierbaum *et al*. (2002) have described the two-body problem with stochastic perturbations using the Fokker–Planck approach. Their whole analysis is based upon the conditional probability density rather than the actual random trajectory as in the ‘Ito stochastic calculus theory’. The idea of using stochastic calculus and Brownian motion into the problem of mechanics was initiated by Chandrasekhar in his paper ‘stochastic problems in physics and astronomy’ (Wax 1954). Chandrasekhar had, in this paper, also considered the problem of a harmonically bound particle perturbed by white noise forces, modelled in the form of stochastic differential equations (SDEs). He also derived a Fokker–Planck equation from the equation of the motion and for field-free cases, he derived the solution to the Fokker–Planck equation. Sheeres *et al*. (2001*a*,*b*, p. 578) have examined the effect of the *stochastic acceleration* on the dynamical model of a spacecraft trajectory. In Scheeres' analysis, the integration of the Ricatti equation was carried out for orbit uncertainty of a one-dimensional force-free motion incorporating a stochastic analysis with correlation time *τ*. The evolutions of variances were derived using the Ricatti equation. For the simplified analysis, they further considered the uncertainties initially uncorrelated. All these examples illustrate the importance of stochastic models in physics and astronomy. However, previous investigations have not considered the two-body dynamics in stochastic differential systems framework. Mann *et al*. (2003) have described the comets and asteroids as the major sources of the dust production. Dust collisional fragmentation, sublimation, radiation pressure acceleration and rotational bursting are the major causes of the dust loss processes. In their paper, they talk about the recent observations of the Sun-grazing comets. The dust particles are small and randomly distributed. The randomly distributed dust produces additional random force on the moving particle. It is always strived that the estimated trajectory be close to the actual trajectory. The accurate modelling procedure increases the effectiveness of the estimated trajectory. The accurate modelling procedure is accomplished by introducing small perturbations felt by the orbiting particle. For these reasons, it is important to introduce the dust perturbation on the orbiting particle. However, the problem of determining the stochastically perturbed two-body model becomes quite difficult and several difficulties arise, since it involves investigations on ‘multi-dimensional diffusion equations’. This might be one of the reasons that system theorists have not attempted this problem previously.

The objective of this paper is to introduce the perturbation into the two-body problem using Brownian motion processes to model the effect of the interstellar dust on the orbiting body. In this paper, we use the concept of ‘central limit theorem’ in the theory of stochastic processes to model the dust population. This paper examines the effectiveness of the dust-perturbed model using *bilinear approximation*, i.e. bilinear stochastic differential systems framework and *second-order approximation* of the system nonlinearity. The original problem of estimating the trajectory of a particle goes back to Gauss. He suggested the least square method for estimating the trajectory of the particle. After Gauss, many more refined techniques were developed, but almost all of them revolve around modelling the motion using ‘deterministic dynamics’. Athans *et al*. (1968) examined the performance of the second-order filter estimating the states of nonlinear plant from discrete noisy observations. The second-order filter was realized by introducing the second-order partials of the system and measurement nonlinearities. In Athan's analysis, deterministic dynamics of the nonlinear plant were the subject of investigation. Julier *et al*. (1995) developed and analysed a linear estimation algorithm for filtering systems with nonlinear process and observation model. The linear estimation algorithm appears to be an alternative to the extended Kalman filter. However, their analysis is useful for the estimation procedure of the ‘stochastic dynamics’, accounting for ‘state-independent perturbations’ in contrast to ‘state-dependent perturbations’. For these reasons, a stochastic model, especially accounting for ‘state-dependent perturbation’ on the orbiting particle produced by the randomly distributed dust, is imperative, if we wish to confirm reality. Accurate modelling procedure prevents inaccurate estimation of the particle trajectory (Kessler & Cicci 1997; Vergez *et al*. 2004). Accurate modelling procedure is accomplished by introducing small perturbation effects felt by the orbiting particle, which leads to the stochastic framework. Thus, our model is different and better than that of the deterministic dynamics. The approach of this paper is a direct demonstration of the application of the theory of stochastic differential systems to the engineering problem of modelling the motion of a satellite in orbit, and to the physical problem of modelling the motion of a planet orbiting a star or the motion of a planetary satellite around a planet. The paper is of theoretical importance in astrophysics, where one wishes to infer something about the nature of the interstellar dust from the measurements of one orbiting body around another. One of the difficult problems is to account for the effect of the randomly distributed dust using the theory of random processes. This paper can be regarded as a starting point in this direction.

This paper is organized as follows: §2 begins by writing down SDEs for the perturbed two-body problem. In §3, we derive the approximate expressions for the propagations of the means and covariances of the state variables occurring in the model. In §4, simulation results are described, leading to the plots of the mean trajectory and the variance of the trajectory about the mean. In §5, we have tried to draw some conclusions regarding the simulation plots.

## 2. A stochastically perturbed model with Brownian motion inputs

The position of the orbiting particle at time *t* in spherical coordinate is given by (*r*(*t*), *θ*(*t*), *ϕ*(*t*)). From the Euler–Lagrange equation, a system of equations of motion iswhere *m* is the mass of the orbiting particle, *V*(*r*)=−*GM*/*r* is the central gravitation potential and *V*′(*r*)=d*V*(*r*)/d*r*. If we wish to take other forces, i.e. the gravitational forces caused by the dust perturbation on the orbiting particle into the account, we assume that their generalized components along (*r*(*t*), *θ*(*t*), *ϕ*(*t*)) are (*F*_{r}(*t*), *F*_{θ}(*t*), *F*_{ϕ}(*t*)). The above equations can be stated asHere, we explain the structure of the forces acting on the orbiting particle caused by dust perturbation. If we assume that apart from the central force, there is ‘dust’ present having uniform dust density *ρ*(*t*) in space, then the contribution to the radial force from this dust at a distance *r* comes entirely from that present inside the radius *r*. The Gauss' law in integral form implies that the flux of the gravitational field outside a sphere of radius *r* comes from the amount of matter inside the sphere. From this, we can deduce that the dust perturbation on the orbiting particle accounts for the dust matter inside the sphere. The mass of dust inside the sphere is (4/3)*πr*^{3}*ρ*(*t*) and hence from inverse attractive square law, magnitude of this radial force is given by 4*πGrρ*(*t*)/3. This model has been described by Goldstein *et al*. (2002, pp. 129–130). If dust is of constant density, then the deterministic model of dust population would be permissible, but in all probability, dust density exhibits highly fluctuating behaviour and white noise model for its effects on the orbiting particle is imperative. To account for the random fluctuations in dust density with time, we take the dust density *ρ*(*t*) proportional to d*B*_{r}/d*t* where d*B*_{r} is the Brownian motion process. We now explain the modelling of dust effects as Brownian motion processes. One natural question is that why should we use the Brownian motion model for the stochastic perturbation of the particle trajectory? The answer is pretty straightforward. The cumulative effect of the gravitational field produced by distributed dust is modelled as the sum of a large number of small effects; each small effect is attributed to one dust particle. When large numbers of small effects are summed up, *central limit theorem* shows that limiting force would be normally distributed. Further, we show that over a small duration [*t*, *t*+Δ*t*], one dust particle acts on the moving particle, then in a finite duration *t*, the net force would be the sum of *t*/Δ*t* small effects and the variance of the cumulative effect will be proportional to the number, i.e. *t*. The Brownian motion process grows proportionally to time *t* and this justifies to model the sum of the small random effects. The Poisson process will need to be introduced into the stochastic model, if one wishes to model the effect of spikes, i.e. dust acts on the particle at the discrete time instants in the form of spurts. But most of the dust in the Solar System is small and distributed. Goldstein has shown that if dust is distributed with a constant density within a sphere of radius *r*, then the radial gravitational force by dust particle acting on the moving particle has the form -*Cρ*(*t*)*r* where *C* is a constant. In our model, radial gravitational force is given by −*σ*_{r}(d*B*_{r}/d*t*)*r* where *σ*_{r}=*Cα*, and *α* is a constant. (d*B*_{r}/d*t*) can assume negative values with equal probability as the positive values. This illustrates that if dust is distributed outside the sphere of motion, sometimes the radial gravitational pull will be attractive and sometimes repulsive. Here, −*σ*_{r}(d*B*_{r}/d*t*)*r* is replaced with *σ*_{r}(d*B*_{r}/d*t*)*r*, since both lead to same statistical characteristics. We also consider adding a random forcing term in the angular equation of the motion. In this way, we are led to the following system of random differential equations:After simplification and setting *θ*=*π*/2, which corresponds to planar motion of the orbiting particle, we haveand *r*, *ϕ* are the radial and angular coordinates, respectively. The left side of the above equation is the radial component of the particle acceleration and the latter denotes the angular component of the particle acceleration. The radial and angular components of the stochastic acceleration are *σ*_{r}*r*(d*B*_{r}/d*t*) and *σ*_{ϕ}(d*B*_{ϕ}/d*t*), respectively. The proofs of the stochastic equation of motion and the conservation of energy and angular momentum are given separately in appendix A. The system of equations can be cast in the SDE format as(2.1)where *v*_{r} and *ω* are the radial and angular velocities of the orbiting particle, respectively. The radial and angular components of the stochastic velocity are *σ*_{r}*r*d*B*_{r} and (*σ*_{ϕ}/*r*)d*B*_{ϕ}, respectively, that confirm the structure of the stochastic term *G*(*x*(*t*), *t*) d*B*(*t*) of the standard differential equation formalism. The special case, i.e. *θ*=*π*/2, corresponds to the planar motion and is valid, provided the dust perturbations act along the radial and angular components of the motion.

The system of SDEs in equation (2.1) can be stated in the standard format aswhere *x*(*t*) are the states, *f*(*x*(*t*), *t*) is the system nonlinearity, *G*(*x*(*t*), *t*), *σ*_{r}, *σ*_{ϕ}) is the noise coefficient matrix, and *σ*_{r} and *σ*_{ϕ} are the unknown parameters. The online estimation of dust density, i.e. the diffusion coefficients *σ*_{r} and *σ*_{ϕ} can be accomplished from experiments by taking measurements on the particle trajectory at discrete time instants. The method used should be maximum likelihood estimate (MLE). This involves discretization of the stochastic equations of motion and writing down the joint density of the discrete observations in terms of the parameters *σ*_{r} and *σ*_{ϕ} and then maximizing the *joint probability density* over these parameters *σ*_{r} and *σ*_{ϕ}. After determining the diffusion coefficients on the basis of MLE, the diffusion coefficients are plugged in the above diffusion equation, i.e. SDEs. As a result of this, *standard diffusion equation* can be taken as(2.2)whereand the diffusion coefficient *Gφ*_{ω}(*t*)*G*^{T}(*x*_{t}, *t*) is given aswhere denote the intensity of the Brownian motion process d*B*(*t*). For the convenience of notations, we consider that the terms and are absorbed into *σ*_{r} and *σ*_{ϕ}, respectively. As a result of this, we have

In the literature of stochastic processes, *standard diffusion equation* is used (Wax 1954; Jazwinski 1970) for *known parameters*. This is the reason why the paper is of theoretical importance in astrophysics, if one wishes to infer something about the nature of the interstellar dust from the measurements on the orbiting particle trajectory. This is in the Ito formalism of SDEs. One can equivalently express this in *Stratonovich framework*. Conversion from ‘Ito’ to Stratonovich formalism is present in standard textbooks (Jazwinski 1970; Pugachev & Sinitsyn 1977) and we do not delve into that here. It is the Ito formalism that has been standard, which is used by all applied mathematicians and engineers, since a process, which satisfies SDE in Stratonovich sense, will not be the Markovian, see Sage & Melsa (1971, p. 102). We now consider 〈*r*〉, 〈*ϕ*〉, 〈*v*_{r}〉 and 〈*ω*〉 for the mean values of the quantities *r*, *ϕ*, *v*_{r}, *ω* and define the deviationsThen, we linearize the system of stochastic differentials in equation (2.1), about the means to get(2.3)(2.4)(2.5)(2.6)Note that a given SDE in equation (2.2) or equivalently by its component version, can be linearized about its mean value by the equationwhich is a bilinear model, owing to the presence of the product term *δx*_{m}d*B*_{j}(*t*). What we have for the stochastic two-body problem is a special case of this. We now cast the system of equations (2.3)–(2.6) in the matrix-vector format. Definethen the bilinear model is given by(2.7)Noting that the matrix *H* can be expressed as *H*_{0}+*H*_{1}, with *H*_{0} the function of time alone and independent of the state vector, and the matrix *H*_{1} is a linear function of the state vector only, it follows that equation (2.7) is a bilinear approximate model for the original perturbed two-body. Let *p*=*p*(*t*, *δr*, *δϕ*, *δv*_{r}, *δω*) denote the probability density of the state vector at time *t*. From the theory of SDEs, it satisfies a Fokker–Planck equationUsing the solution to the Fokker–Planck equation, we can calculate the conditional probability density of the state vector at an arbitrary set of the discrete time points. Specially, we know that the conditional probability density *p*(*t*, *δr*, *δϕ*, *δv*_{r}, *δω*|*δr*′, *δϕ*′, , *δω*′) of the state variables at the time *t*, given those at time *s*<*t* satisfies the Fokker–Planck equation with respect to the variables at time *t* and the initial condition.

Using shorthand notationthen, the Markov property of the solution implies that the joint density of the state at time *t*_{1}<*t*_{2}<*t*_{3}< … <*t*_{n}, is given byThe exact evolutions of the mean and variance using the Fokker–Planck equation (Sage & Melsa 1971, p.100) and the definition of the differential of the expectation of the scalar function of the state vector *x*_{t}, arewhere *f*(*x*_{t}, *t*) is the system nonlinearity and *G*(*x*_{t}, *t*) are the noise coefficient matrix. In matrix-vector format, the evolutions are

The evolutions, using the second-order approximation, areThe evolutions, using the bilinear approximation, become a special case of the evolutions using the second-order by vanishing the second-order partial terms in the mean and variance evolutions.

## 3. Mean and variance evolutions of the dust-perturbed model

The analytical and numerical solutions of the mean and variance evolutions for the nonlinear stochastic systems are not possible, since their evolutions are infinite dimensional and require the knowledge of higher-order moment evolutions (Maybeck 1982). To address these, we introduce nearly ‘Gaussian assumption’ and approximations to the system nonlinearity and noise coefficient matrix around the mean of the state. As a result of these, we derive approximate mean and variance evolutions. Approximate evolutions preserve some of the qualitative characteristics of the exact evolutions. In this paper, the bilinear and second-order approximations are the subject of investigation.

### (a) Bilinear approximation

Given the stochastic model of the kind of equation (2.2), the mean propagation up to the bilinear approximation can be obtained by replacing the state vector *x*(*t*) on the right-hand side by taking its mean value. Specifically, taking its expectation and using 〈*G*(*x*(*t*), *t*) d*B*(*t*)〉=0 givesThe bilinear approximation of this equation yieldsWhen this equation is specialized to the perturbed two-body problem, described by equation (2.7), then the equations obtained for 〈*r*〉, 〈*ϕ*〉, 〈*v*_{r}〉 and 〈*ω*〉 are those corresponding to the unperturbed orbit.These are equivalent to the pairwhich yield the energy and angular momentum integrals for the mean trajectory. These are to be understood as only approximations to the true mean evolution. The stochastic term contributes to a *bias* in the trajectory, and to calculate the bias, we must make a second-order approximation to the driving functions. In other words, a second-order approximation of the function *f* giveswhere *P*=(*P*_{ij})=(〈*δx*_{i}*δx*_{j}〉) is the covariance matrix of the state vector *x*(*t*).

Here, we shall be content with bilinear approximation. The covariance matrix elements of the state variables areThe evolution equations for these covariance matrix elements with time have been derived in Appendix B. These variance evolutions can be also derived using the results given in Jazwinski (1970, p.363). We list the final equations here.

### (b) Second-order approximation

The mean trajectory for the dust-perturbed particle using bilinear approximation does not include *variance term* in mean evolution. The term 〈*Gφ*_{w}(*t*)*G*^{T}(*x*(*t*), *t*)〉 in variance evolution accounts (Jazwinski 1970, p. 363) for the stochastic perturbation felt by the orbiting particle. For this reason, the bilinear approximation does not preserve perturbation effects in mean evolution. On the other hand, the variance evolution using bilinear approximation for the dust-perturbed model includes perturbation effects, i.e. *Gφ*_{ω}(*t*)*G*^{T}(〈*x*(*t*)〉, *t*). The main interest in bilinear approximate model is owing to computational simplicity as well as it representing a ‘nearly’ linear yet nonlinear system (Bruni *et al*. 1974; Elliott 1998).

In order to account for the stochastic perturbation in the mean evolution, we use the second-order approximation in the mean evolution. The second-order approximation includes ‘the second-order partials’ of the system nonlinearity *f*(*x*(*t*), *t*) and variance terms in the mean trajectory, which leads to better *estimation* of the trajectory (Athans *et al*. 1968). Making use of second-order approximation of the mean evolution (Jazwinski 1970, p. 363), and equation (2.2) for the nonlinear problem of concern here, we haveMaking use of the second-order approximation of the variance evolution (Jazwinski 1970, p. 363) and equation (2.2) for the nonlinear problem of this paper, we haveThe variance evolution using the second-order approximation is different from the bilinear approximation owing to the second-order partial of the diffusion coefficient, i.e. the term Owing to this term, the expression for the variance evolution of the radial velocity of this paper, i.e. , using the second-order approximation, involves an additional term in contrast to the bilinear approximation. The expression for the variance evolution of the angular velocity, i.e. d*P*_{ω}, using the second-order approximation, accounts for a correction term , in contrast to the bilinear approximation as well. The expressions for the variance evolution, i.e. d*P*_{r} using both the approximations will be the same, since the contribution from the second-order partial term will be zero; see the expression for corresponding to equation (2.1). The same argument holds for the variance evolution of the phase angle, i.e. d*P*_{ϕ}. The variance evolutions using the second-order and bilinear approximations would be different from the Ricatti equation, since this paper discusses the nonlinear stochastic differential system, see equation (2.1), on the other, the Ricatti equation is useful for the linear stochastic differential system.

## 4. Simulation results

The standard discretization method for simulating an ordinary differential equation by means of a difference equation is employed to simulate the mean and covariance evolution equations. The simulations of the mean and variance evolutions are accomplished using a simple, but effective, finite difference method-based numerical scheme. Thus, the simulating SDEwe would replace d*x*(*t*) by *x*(*t*_{k+1})−*x*(*t*_{k}), d*t* by *t*_{k+1}−*t*_{k} and d*B*_{t} by , is a standard normal variable (Jazwinski 1970, p. 177). For the nonlinear problem considered here, the component of *x*(*t*) has four means 〈*r*〉, 〈*ϕ*〉, 〈*v*_{r}〉 and 〈*ω*〉, and sixteen covariances *P*_{rr}, *P*_{rϕ}, … *P*_{ωω}. Note that the covariance matrix is a positive semi-definite matrix at each time instant, but by symmetry, it has only 1+2+3+4=10 independent components, which can be taken as the diagonal entries, above and below the diagonal. Thus, the state vector *x*(*t*) has 14 independent components for the simulations. The initial conditions are chosen asfor the states *x* and *y*.The initial state covariance *P*_{xy}(0) is considered zero. The initial conditions considered here are in canonical system of units. Astronomers adopt a normalized system of units, i.e. ‘canonical units’, for the simplification purposes (Bate *et al*. 1971, pp. 40–41). In canonical units, the physical quantities are expressed in terms of time unit (TU) and astronomical unit (AU), where 1 TU=5.0226757×10^{6} s and 1 AU=1.4959965×10^{8} km (Bate *et al*. 1971, p. 429). The diffusion coefficients are chosen as *σ*_{r}=0.0121(TU)^{−3/2} and *σ*_{ϕ}=2.2×10^{−4}(AU/(TU)^{3/2}). The initial covariance matrix being chosen as identically zero, which corresponds to the physical situation that, initially, the particle is located at a definite point in the space and imparted a definite radial and angular velocities. This paper considers that there are no uncertainties in the initial conditions, i.e. the initial states of the orbiting particle are known accurately, which imply that the system is deterministic at *t*=*t*_{0} and becomes stochastic at *t*≻*t*_{0} owing to the perturbation. This allows the initial covariance matrix to be zero. This makes the contribution to the variance evolution from the ‘system nonlinearity’ coupled with ‘initial error terms’ will be zero at *t*=*t*_{1}. Only the contribution to the variance evolution at *t*=*t*_{1} comes from the perturbation term (*Gψ*_{w}*G*^{T})(*x*_{t}, *t*). For *t*≻*t*_{1} the contribution to the variance evolution comes from the system nonlinearity as well as the perturbation term. This assumption allows to study the effect of the dust perturbation explicitly on the orbiting particle, the intention of this paper. The problem of analysing stochastic systems, involving uncertainties in initial conditions and the stochastic perturbation, can be accomplished via numerical simulations. However, the problem of separating the contributions from uncertainties in initial conditions and the stochastic perturbation becomes quite difficult, owing to the coupling terms, e.g. The initial condition chosen corresponds to those values that ensure that the unperturbed motion corresponds to the orbiting of the particle about the centre at a given radial distance. Thus, these could be applied to both planetary motion about the Sun and the Earth satellites. When applied to the latter situation, we can make predictions about how smoothly the satellite would revolve around the Earth and to what extent would the true perturbed trajectory deviate from the desired one. Other initial conditions can also be chosen, but they invariably lead to either the satellite collapsing in spirals towards the centre or its eventual escape from the gravitational field of the centre particle. The problem of simulating nonlinear filters involves the parameter *φ*_{η}, intensity of the noise d*η*_{t} introduced into observations. The standard approach of simulating the observation is that the parameter *φ*_{η} in variance of the noise is chosen using the criterion that the statistical contribution is a smaller part of the deterministic contribution. The same approach we have adopted for the simulation of the stochastic equation of motion, since both equations, i.e. observation and motion equations use the SDE formalism. The values of diffusion coefficients are selected so that the contribution to the force coming from the random part is smaller than the force coming from the deterministic part. It has been chosen for simulational convenience only. However, the online estimation of dust density, i.e. the diffusion coefficients *σ*_{r} and *σ*_{ϕ} can be accomplished from experiments by taking measurements on the particle trajectory at discrete time instants. The MLE should be used for the estimation of the diffusion parameters. This paper is not intended to study the estimation of the diffusion parameters *σ*_{r} and *σ*_{ϕ} from the measurements, but only to study the effectiveness of the dust-perturbed particle model using numerical simulations.

The potential of the central force has been taken as *V*(*r*)=−(1/*r*) The following basic idea regarding the order of the magnitudes is used. Given an SDE of the form , if we assume that the unperturbed orbit defined by the differential equation d*x*(*t*)=*f*(*x*(*t*), *t*)d*t* is periodic with time *T*, then we compute the average of the function *f* over one period of the unperturbed motion and denote this nominal value of *f* by *f*_{nom}. We also compute the average of *g*(*x*(*t*), *t*) over one period and denote this nominal value of *g* by *g*_{nom}. Then, using the Ito's formula (d*B*(*t*))^{2}=d*t*, we find that the nominal magnitude of the random perturbing force equals , and the force that yields the unperturbed motion is |*f*_{nom}|d*t*. We let d*t*=Δ, the sampling interval and choose the parameter *σ* so that the perturbing force is around *α*=(1/10)th part of the deterministic component of the force. This gives or equivalently, . The time period of the unperturbed motion can also be estimated using elementary mechanics. The nominal radial distance can be taken with its initial condition as *r*(0) and with the potential as −(1/*r*). The centripetal acceleration satisfies the relation *rω*^{2}=1/*r*^{2}, where *ω* is the angular velocity. Thus, *ω*=1/*r*^{3/2} and the time period *T*=2*π*(*r*(0))^{3/2} Note that the mass has been chosen as the unit.

The natural question to ask after performing the simulation of noisy trajectory in perturbed two-body problem is that how our graphs relate to the actual statistics of the processes considered here. To understand this, we consider two separate cases. The first is the analysis of the perturbations in the mean trajectory due to noise and the second is that variance–covariance evolutions of the state variables in the perturbed model. The differential equation for the mean differs from the ideal unperturbed case in a sense that the covariances enter into the picture and the evolutions of covariances have no counterpart for the unperturbed problem. The approximate mean and covariance evolutions can be solved by numerical techniques. We can draw some useful conclusions about the nature of the perturbation introduced by the random forces by looking at the graph of the mean and covariance evolutions of the state variables. The graph of the mean radial distance nearly coincides with the graph of the unperturbed trajectory for a considerably long duration and thereafter slightly deviates above the latter, see figure 1. Noting that , it follows that for a long time 〈*v*_{r}(*t*)〉 coincides with the velocity of the unperturbed trajectory and thereafter becomes larger. From the differential equation for 〈*v*_{r}(*t*)〉, it can be inferred that the correction term becomes significant only after a long time and the integral of this quantity is positive after a long time. The graph of the mean phase trajectory almost coincides with unperturbed phase trajectory, see figure 2. This means that 〈*ω*〉 almost coincides with the angular velocity of the unperturbed trajectory, see figure 4. This further means that the perturbation term in the differential equation for 〈*ω*〉, i.e. is not significant. We infer that *P*_{rr}, , *P*_{rω}, remain negligible over a long duration. The graph of the mean radial velocity 〈*v*_{r}〉 appears to show a small oscillatory behaviour about the unperturbed graph after a long duration, see figure 3. This means that the perturbation term changes sign frequently. This can be attributed to the fact that *P*_{rω} is oscillatory. The integrated version of this term, which appears in the equation for 〈*r*〉, as we have seen, is nearly zero for a long time and then positive. Then we conclude that the first two terms give a result after a long time, the integral of which keeps accumulating to give an increasing function, causing 〈*r*〉 to deviate above the unperturbed radial distance trajectory after a long time, while oscillatory term 2〈*ω*〉*P*_{rω} averages out to be zero, not affecting the graph of the mean radial distance 〈*r*〉. The variance analysis of nonlinear problem of this paper shows the bounded variance characteristics, i.e. an increase in the variance is followed by a decrease. This confirms the efficacy of the bilinear and second-order approximations for analysing the problem, see figures 5–8. This variance property continues for larger durations. The smaller size of sampling interval *t*_{k+1}−*t*_{k} and larger number of sampling intervals, i.e. *n*=*t*/(*t*_{k+1}−*t*_{k}), for a given duration *t* must be chosen for accurate numerical simulations of the mean and variance trajectories. Owing to inadequate choice of the sampling interval, the numerical simulation may show instability in the trajectories (Julier *et al*. 1995). In this paper, Matlab programs have been written on the basis of analytical results derived using sampling interval 0.01 TU for a duration 15 TU, which is believed to be an accurate choice of the interval. Numerical simulations demonstrated in figures 5–8 indicate that the variance evolution using the second-order approximation results reduced variances of the state variables rather than bilinear. This illustrates that the second-order approximation of the mean evolution produces less random fluctuation in the mean trajectory. The second-order approximation of the mean evolution involves the second-order partials of the system nonlinearity *f*(*x*(*t*), *t*), i.e. which accounts for the biasing correction term in mean evolution (Athans *et al*. 1968). As a result of this, the second-order approximation produces less biased estimate rather than the bilinear. The third-order approximate estimation procedure may be the alternative to the approximations, i.e. the bilinear and second-order approximations. The third-order approximate estimation procedure can be realized by introducing the third-order partials of the system nonlinearity *f*(*x*_{t}, *t*) and diffusion coefficient (*Gφ*_{ω}*G*^{T})_{ij}(*x*_{t}, *t*) into the mean and variance evolutions (Sharma *et al*. 2006). However, the popularity of the higher-order approximate estimation procedure may be hampered owing to formidable complexity introduced by the higher-order partials of the system nonlinearity and diffusion coefficient. The mean evolution, using the bilinear approximation, does not preserve the perturbation effect, since it does not involve the variance term in the mean evolution. For these reasons, the second-order approximation is ‘somewhat’ more accurate for estimating the states rather than the bilinear approximation and more attractive in contrast to the third-order approximation. The variance evolution for unperturbed particle motion at every discrete instant will be zero, since initial covariances are zero. The variance evolution for unperturbed trajectory does not include the perturbation term 〈*Gφ*_{ω}(*t*)*G*^{T}(*x*_{t}, *t*)〉 as well. For these reasons, comparison of the variance evolutions of the dust-perturbed model is not given with *unperturbed two-body dynamics*.

The Brownian motion process is accurate to describe the effect of dust perturbation on the orbiting particle, since it justifies the fact that the dust perturbation introduces slight modifications into the dynamics of the orbiting particle, see figures 1–4. Such analysis confirms the physical situation, i.e. the effect of dust perturbation on the orbiting particle, see Goldstein *et al*. (2002, pp. 129–130). However, this paper discusses the dust perturbation in stochastic sense rather than the deterministic (Goldestein *et al*. 2002, pp. 129–130).

## 5. Conclusion

In this paper, we have started with a stochastic model for describing the effect of the gravitational field produced by the randomly distributed dust on the motion of an orbiting body in a central potential. The resulting stochastic model was linearized about the mean trajectory, leading to a bilinear model for the stochastic motion. From the bilinear model, we arrived at the Fokker–Planck equation for the conditional probability density of the perturbed trajectory. We then derived a set of approximate differential equations, which describe the evolutions of the mean and variance random trajectories with time. However, bilinear approximation for the mean trajectory produces biased estimates, since it does not include the second-order partials of the system nonlinearity *f*(*x*(*t*), *t*). The effectiveness of the dust-perturbed two-body model has been examined on the basis of the bilinear and the second-order approximations. It is shown that the second-order approximation for the mean trajectory gives better estimates rather than the bilinear approximation because second-order partials account for biasing correction terms in mean trajectory. As a result of this, the second-order approximation increases the accuracy of the estimated trajectory. Numerical simulations indicate that the second-order approximation for the variance evolution results in reduced variances rather than bilinear approximation as well.

The main contribution of this paper is to develop the dust-perturbed two-body model, if we wish to confirm the reality. The motivation for analysing the dust-perturbed model has come from the fact that *accurate modelling procedure* coupled with an effective filtering algorithm increases the accuracy of the estimation procedure. In this paper, accurate modelling procedure for the orbiting particle has been accomplished by introducing the dust perturbation on the orbiting body. Another objective achieved by this paper is that the theory of ‘Brownian motion and stochastic calculus’ can be used to study and model the effect of dust on the orbiting planet. This paper suggests that the ‘gap’ between classical mechanics, statistical mechanics and the theory of stochastic processes is not as great as generally believed.

A good source that describes the nature of the dust population in the planet–Sun region is given in Mann *et al*. (2003).

## Footnotes

- Received October 4, 2006.
- Accepted November 28, 2006.

- © 2007 The Royal Society