Removing the trend of drift induced from acceleration noise for LISA

# Removing the trend of drift induced from acceleration noise for LISA

## Abstract

In this paper we demonstrate a methodology to remove the power of the drift induced from random acceleration on LISA proof mass in the frequency domain. The drift must be cleaned from LISA time series data in advance of any further analysis. The cleaning is usually performed in the time domain by using a quadratic function to fit the time series data, and then removing the fitted part from the data. Having Fourier transformed the residuals, and then convolved with LISA transfer function, LISA sensitivity curve can be obtained. However, cosmic gravitational-wave background cannot be retrieved with this approach due to its random nature. Here we provide a new representation of power spectrum given by discrete Fourier transform, which is applied to find the function of the drift power for the cleaning in the frequency domain. We also give the probability distribution used to analyze the data in the frequency domain. We combine several techniques, including Markov Chain Monte Carlo method, simulated annealing, and Gelman & Rubin’s method, with Baye’s theorem to build the algorithm. The algorithm is utilized to analyze 24 simulations of LISA instrumental noise. We prove that the LISA sensitivity can be recovered through this approach. It can help us to build algorithms for some tasks which are must accomplished in the frequency domain for LISA data analysis. This method can be applied to other space-borne interferometers if charges on their proof masses cannot be perfectly cancelled.

## I Introduction

In the LISA data stream, intrinsic instrumental noise falls into two categories: shot noise and acceleration noise Bender et al. (1998). Acceleration noise, caused by the residual Coulomb force induced from the imperfect cancellation of charges on proof-masses, is dominant in the low-frequency range, resulting the sensitivity proportional to roughly below 2 mHz. Optical-path noise, including mainly shot noise and beam-pointing error, is dominant in the high frequency range, leading the sensitivity declining proportional to the frequency above 10 mHz due to the falloff of the antenna transfer function.

The theoretical LISA sensitivity can be obtained from the various types of noise spectral densities directly Bender et al. (1998), whereas when we deal with LISA time series raw data, the trend of the drift of proof mass in the time series due to the random acceleration shall be removed at first. In general, the removal is performed in the time domain. By Fourier transforming the residuals, the LISA sensitivity is recovered. However, if stochastic gravitational-wave background exists in the data, bias might be induced if the trend is firstly removed in the time domain and then the background is extracted in the frequency domain. The analysis given by such sequential subtraction might be inaccurate, particularly if signals have overlaps. For instance, if a data set containing two overlapped signals is fitted by a linear filter where the signal is parametrized by a rectangular function of amplitude and location, the filter may extract a stronger output around the overlapped region instead of one of the exact signals. This is because what fitting does is to minimize the difference between data and model. For that reason, the parameter estimation and the removal of the trend in the LISA data analysis should all be performed either in the time domain or in the frequency domain.

The cosmological sources in the very early universe are randomly distributed across the sky, emitting gravitational waves with various amplitude and frequency. If they are not strong enough to be located by LISA, their incoherent signals will form a continuum and be entangled with instrumental noise. In order to gain cosmological information, the functional form of the trend and cosmic gravitational-wave background (CGB) are both needed to be understood.

The waveform of CGB in the time domain cannot be obtained due to its stochastic nature, making the separation of CGB from the time series data very difficult. It is natural to extract the CGBs in the frequency domain since the function of the power spectra can be written analytically. Nevertheless, the method to remove the drift in the frequency domain is unknown. In this paper, the function of the power spectrum of the drift will be derived, and a method to remove the trend in the frequency domain for LISA data analysis will be developed.

In section II, we will review the time series of LISA instrumental noise, and demonstrate the approach conducted in the time domain to obtain the LISA sensitivity. In section III, we will find a new representation of Fourier power spectrum, and use the representation to derive the power spectrum of the drift trend. In section IV, the derivation of probability distribution of noise power will be provided. In section V, the algorithm of parameter estimation will be introduced, followed by the result of data analysis. Finally, in section VI, conclusion will be given.

## Ii Remove Displacement Caused by LISA Random Acceleration in the Time Domain

Energetic particles keep hitting the proof mass, producing random accelerations on it continuously. The perturbations from random accelerations accumulate, and gradually depart the proof mass from its free-falling trajectory Cornish et al. (). The LISA sensitivity curve cannot be obtained just by directly Fourier transforming the drift induced from such accelerations. The drift must be fitted by a quadratic function and the best fit must be removed. Then the Fourier transformed residuals can represent the LISA noise level. In this section we will demonstrate how LISA sensitivity curve is obtained from fitting in the time domain. Firstly we will simulate the drift induced from the acceleration, and then fit the simulated data. We utilized Gaussian distribution to simulate the shifts Cornish et al. (). The acceleration noise power spectral density is suggested as Bender et al. (1998). Dividing it by sampling rate and then taking square root, we will obtain the average of acceleration. Using this average as the standard deviations, we can draw a time series of random accelerations. Having double integrated the random accelerations, we will have the time series drift. Here the sampling rate is set as 1.5 second, and 8192 data are simulated. The simulation of drift is shown by the black curve in Fig. 1.

Dividing the drift by arm-length the dimensionless drift can be obtained. Having Fourier transformed the dimensionless drift, we will get the strain amplitude, as shown by the solid line in Fig. 2. As shown in the Fig. 2, the frequency dependence of the strain is approximately rather than as indicated by the sensitivity curve in the LISA Prephase A study. In order to recover the LISA sensitivity the trend of drift must be removed.

Since the drift is induced from random acceleration, it is natural to model its trend by a quadratic function where a, b, c are unknown parameters. Using the function to fit the drift, and then remove the trend, which is identified as green dash-dot curve in the Fig. 1, from the drift, the residual can be obtained. By Fourier transforming the residual, the true noise level can be recovered. The strain amplitude of the residual is shown by the green dot curve in Fig. 2. Its frequency dependance is proportional to as shown in the LISA Prephase A study Bender et al. (1998).

In addition to the quadratic equation, the linear equation and the cubic equation are tested to remove the trend as well. The best fits given by the linear and cubic equation are indicated by red and blue line, respectively. From Fig. 2 it is noticed that the amplitude of the residual given by the fitting with linear function is lower than the uncleaned displacement by one order of magnitude, but it is inversely proportional to the frequency , not to . The amplitude of the residual given by the fitting with cubic function is still proportional to , and is as the same level of the amplitude corresponding to quadratic fitting. One more parameter used in the cubic fitting does not provide extra benefit.

Here we demonstrated that the LISA sensitivity curve is obtained by removing the trend of the drift. There is no problem if only point sources involve in data since point sources can be analyzed simultaneously with removing the trend in the time domain. However, if data contain the waves produced by stochastic sources, such as astrophysical foreground or cosmological background, which can only be analyzed in the frequency domain, removing the trend beforehand in the time domain would induce a bias in the analysis of the sources. A technique to deal with the trend in the frequency domain is necessary.

## Iii Expected Power of Instrumental Noise

In last section the way to remove the trend of the drift in the time domain was described. In this section, a method to remove the trend directly in the frequency domain will be shown.

Suppose is the time series of data, and is the Fourier transform of the data, the power is given by

 |Hn|2 = Hn×H∗n (1) = N−1∑k=0N−1∑k′=0hkhk′exp{2πi(k−k′)nN}, (2)

where is index, and is total number of data. The summation is usually calculated firstly over one index and then the other, or vice versa. This is implied by the functionality of summation. However, what does matter is summing the term over all k and k’ on the grid. The implied procedure is not the only approach to carry out the calculation. We calculate the summation along the diagonal arrays as shown in the Figure 3 rather than the regular procedure. The terms along red lines have the property that the difference of and is fixed. Thus we introduce an index to indicate their difference . equals on the diagonal line, so their phase is cancelled. The term turns to be . Considering the arrays corresponding to and , the terms along those diagonal arrays have symmetric mathematical expression. Their phases, and , have same magnitude but opposite sign. Because of that, the sum of those two is where . Similarly the other symmetric terms with can be combined to where . Therefore, we can rewrite Eq. (2) as the following form

 |Hn|2=N−1∑k=0h2k+2N−1∑t=1N−1−t∑k=0hkhk+tcos2πntN. (3)

The benefit of new representation is that the time series data can be separated into autocorrelation terms and cross-correlation terms . When we deal with the Fourier component of random noise, the autocorrelation term will remain and cross-correlation term will vanish as their ensemble average is taken. The ensemble average is associated with the standard deviation of the time series random noise. As a conclusion, with this representation, the power spectrum of random noise can be expressed by some statistical properties of the time series data.

### iii.1 Shot Noise

Time series of shot noise can be characterised as independent Gaussian noise. This implies that their ensemble average is zero, and they are not correlated, which can be expressed by for any . The standard deviation of Gaussian distribution should equal to the ensemble average of square of data if data length is infinite long. Knowing this we can apply Eq. (3) to find a formula for the Fourier component of shot noise. Substituting time series shot noise into Eq. (3) we will obtain

 ⟨Pn⟩ = ⟨|Hn|2⟩=1NN−1∑k=0⟨h2k⟩ (4) + 2NN−1∑t=1N−1−t∑k=0⟨hkhk+t⟩cos2πntN = σ2. (5)

It is not surprised that the formula for describing the power spectrum is a constant since shot noise is white noise.

### iii.2 Acceleration Noise

To find the expression for the power spectrum of acceleration noise, we review the characteristics of time series drift induced from random acceleration at first. The gross feature of can be recognised as independent Gaussian noise as well

 P(ai|I)=1√2πaexp{−a2i2a2}. (6)

First of all, the ensemble average is zero. Secondly, the random acceleration noise is independent, so we have . Thirdly, the ensemble average of square of equals to where is the standard deviation of the distribution.

Now we construct the Fourier component of the displacement noise step by step. We begin this work by establishing the Fourier component of velocity noise. Suppose is the time series of velocity noise associated with random acceleration and its initial value is zero. We assume that it obeys the equation

 vk=vk−1+akΔ. (7)

With the initial condition we can derive

 vk=(a1+a2+...+ak)Δ, (8)

and

 vk+t=(a1+a2+...+ak+t)Δ. (9)

Then, we know

 v2k = Δ2(k∑i=1ai)(k∑j=1aj) (10) = Δ2(k∑i=1a2i+2k∑i=1k∑j=i+1aiaj), (11)

and

 vkvk+t = Δ2(k∑i=1ai)(k+t∑j=1aj) (12) = Δ2(k∑i=1ai)(k∑j=1aj+k+t∑j=k+1aj) (13) = Δ2[k∑i=1a2i+2k∑i=1k∑j=i+1aiaj + k∑i=1k+t∑j=k+1aiaj]. (14)

Substituting Eq. (11) and (14) into Eq. (3), we can expand the power of velocity noise as

 |Vn|2 = Δ2NN−1∑k=1k∑i=1a2i+2Δ2NN−1∑k=1k∑i=1k∑j=i+1aiaj (15) + 2Δ2NN−2∑t=1N−1−t∑k=1cos2πntN×[k∑i=1a2i + 2k∑i=1k∑j=i+1aiaj+k∑i=1k+t∑j=k+1aiaj].

Then we can calculate the expected power of velocity noise by taking ensemble average of

 ⟨|Vn|2⟩ = Δ2NN−1∑k=1k∑i=1a2+2Δ2NN−2∑t=1cos2πntNN−1−t∑k=1k∑i=1a2 (16) = Δ2a2NN(N−1)2+Δ2a2N × N−2∑t=1[t2−(2N−1)t+N(N−1)]cos2πntN

and all other terms are zero. The term and are given as

 N∑t=1costx=sinNx2sinx2cos(N+12x), (17)

and

 N−1∑t=1tcostx=Nsin2N−12x2sinx2−1−cosNx4sin2x2. (18)

Substituting into in Eq. (17), the right hand side turns to zero. The left hand side of Eq. (17) is

 N∑t=1cos2πntN=N−2∑t=1cos2πntN+cos2πn+cos2πnN. (19)

With the information on the both sides of Eq. (17), we obtain

 N−2∑t=1cos2πntN=−1−cos2πnN. (20)

Substituting into in Eq. (18), it is then simplified as

 N−1∑t=1tcos2πntN = (N−1)cos2πnN+N−2∑t=1tcos2πntN (21) = −NsinπnN2sinπnN=−N2.

Thus, we acquire

 N−2∑t=1tcos2πntN=−N2−(N−1)cos2πnN. (22)

For the expression of , it can be derived by differentiating with . The close form of is

 N−1∑t=1tsintx=sinNx4sin2x2−Ncos2N−12x2sinx2. (23)

Differentiating it with respect to on the both sides, we get

 N−1∑t=1t2costx = N(2N−1)sin2N−12x4sinx2 (24) + N[32cosNx+12cos(N−1)x]4sin2x2 − sinNx4sin3x2cosx2.

Substituting into , it gives

 N−1∑t=1t2cos2πntN = (N−1)2cos2πnN (25) + N−2∑t=1t2cos2πntN = −N(2N−1)4+N[32+12cos2πnN]4sin2πnN = N2sin2πnN−N22. (26)

Moving the first term in Eq. (25) to Eq. (26), we obtain the expression

 N−2∑t=1t2cos2πntN=−(N−1)2cos2πnN+N2sin2πnN−N22. (27)

Substituting Eq. (20), Eq. (22), and Eq. (27) into Eq. (16), we have a compact form for the power of velocity noise

 ⟨|Vn|2⟩=Δ2a22sin2πnN. (28)

Next, we consider the displacement noise of proof mass due to random acceleration noise. Suppose is the time series of position noise with the initial condition of . From Newtonian dynamics, the time series can be described by the following regression relationship

 xi=xi−1+viΔ ∀i:1∼N−1. (29)

From Eq. (29) and the initial condition we know

 xi = xi−1+viΔ=xi−2+(vi−1+vi)Δ (30) = (v1+⋯+vi)Δ.

By Eq. (8) the velocity can be calculated from the random acceleration noise

 xi = [a1+(a1+a2)+⋯+i∑k=1ak]Δ2 (31) = Δ2i∑k=1(i−k+1)ak.

From Eq. (31) we can compute and

 x2i = [Δ2i∑j=1(i−j+1)aj][Δ2i∑k=1(i−k+1)aj] (32) = Δ4[i∑j=1(i−j+1)2a2j + 2i∑j=2j−1∑k=1(i−j+1)(i−k+1)ajak],
 xixi+t = [Δ2i∑j=1(i−j+1)aj][Δ2i+t∑k=1(i+t−k+1)aj] (33) = Δ4i∑j=1i+t∑k=1(i−j+1)(i+t−k+1)ajak.

Substituting Eq. (32) and Eq. (33) into Eq. (3), then we can obtain the expected power of the drift

 ⟨|Xn|2⟩ = 1NN−1∑i=0⟨x2i⟩+2NN−1∑t=1N−1−t∑i=0⟨xixi+t⟩cos2πntN (34) = Δ4a2NN−1∑i=1i∑j=1(i−j+1)2+2Δ4a2NN−2∑t=1cos2πntN ×N−1−t∑i=1i∑j=1i+t∑k=1(i−j+1)(i+t−k+1)δjk = Δ4a2N(N2−1)12+Δ4a26NN−2∑t=1cos2πntN ×[−t4+2Nt3+t2−2N3t+N2(N2−1)].

We now have the expressions of and but we still need those of and to complete Eq. (34). They can be obtained by differentiating and with respect to x once and twice, respectively, and then substituting for . Then we will get

 N−1∑t=1t3cos2πntN=−N32+3N24sin2πnN, (35)

and

 N−1∑t=1t4cos2πntN=−N42+N(N2+1)sin2πnN−3N2sin4πnN. (36)

Substituting Eq. (18), Eq. (27), Eq. (35), and Eq. (36) into Eq. (34), we finally obtain

 ⟨|Xn|2⟩=NΔ4a24⎡⎣1sin4πnN+N2−13sin2πnN⎤⎦. (37)

There are two terms in Eq. (37) involved in the power of drifts. One is inversely proportional to , the other is . If n is small, we can approximate and by and , respectively. The strain amplitude corresponding to these two terms can be obtained by taking square root. There is not only term, but also term. Moreover, the strain amplitude of drifts is dominant by term especially in the high frequency range. As we have done in the analysis in the time domain, in order to reveal the true noise level, the power of drift must be removed in the frequency domain. Eq. (37) can be employed to estimate the drift power.

## Iv Probability Distribution of Noise Power

The cosmic gravitational-wave background entangles with noise in time domain, forbidding the analysis performed. Hence, we intend to analyze a data set in the frequency domain through its power spectrum. The measurement is the sum of signal and Gaussian random noise where the Gaussian noise is drawn from the distribution

 P(n)=1√2πσexp{−n22σ2}. (38)

#### Fourier Amplitude of Noise

Giving a time series of noise drawn from Eq. (38), its Fourier amplitude is given by

 ~nk = N−1∑j=0Wkjnj (39) = Re(n0+Wn1+⋯+Wk(N−1)nN−1) + i Im(n0+Wn1+⋯+Wk(N−1)nN−1) (40) = N−1∑j=0Rkjnj+iN−1∑j=0Ikjnj (41) ≡ ~nrk+i~nik (42)

where is the number of noise data, and . and are the real part and imaginary part of , and and denote the real part and imaginary part of , respectively.

#### Probability of ~nrk and ~nik

From the probability distribution of time series instrumental noise given by Eq. (38), we can derive the probability distribution of by marginalizing the probability distribution over :

 P(~nrk) = ∫⋯∫∞−∞dn0dn1⋯dnN−1 (43) × P(~nrk,n0,n1,⋯,nN−1| I)

where is the probability distribution of , being true given the background information . From Baye’s theorem, can be decomposed into the product of likelihood function which is the probability distribution of for a given , and prior function

 P(~nrk) = ∫⋯∫∞−∞dn0⋯dnN−1P(n0,⋯,nN−1) (44) × P(~nrk| n0,⋯,nN−1)

Since is completely determined by Eq. (42) if , , , are known, must satisfy

 P(~nrk| n0,n1,⋯,nN−1)=δ(~nrk−N−1∑j=0Rkjnj) (45)

where is the delta function of x. Moreover, if , , , are independent random variables, can be written by

 P(n0,n1,⋯,nN−1)=P(n0)P(n1)⋯P(nN−1). (46)

Therefore, Eq. (44) can be rewritten as

 P(~nrk) = ∫⋯∫∞−∞δ(~nrk−N−1∑j=1Rkjnj−n0)P(n0)dn0 (47) × P(n1)⋯P(nN−1) dn1⋯dnN−1.

Since the delta function is involved, the integration over produces directly where equals . Since is a Gaussian distribution, we will have

 P(~nrk)=1√2πσ∫⋯∫∞−∞N−1∏l=1P(nl)exp{−(~nrk−∑N−1j=1Rkjnj)22σ2}dn1⋯dnN−1. (48)

Combining with the exponential term, it becomes

 P(~nrk)=12πσ2∫⋯∫∞−∞dn2⋯dnN−1N−1∏l=2P(nl)∫∞−∞dn1exp{−n21+(~nrk−∑N−1j=2Rkjnj−Rkn1)22σ2}. (49)

where is given by Eq. (38) for any integer j. Denoting K as the integral term over in Eq. (49), and as , the integral can be simplified as

 K = ∫∞−∞dn1exp{−(1+R2k)n21−2αRkn1+α22σ2} (50) = √2πσ21+R2kexp{−α22(1+R2k)σ2}.

Substituting Eq. (50) into Eq. (49) we have

 P(~nrk) =1√2π(1+R2k)σ2∫⋯∫∞−∞dn2⋯dnN−1⋯ (51) N−1∏l=2P(nl)exp{−(~nrk−∑N−1j=2Rkjnj)22(1+R2k)σ2}.

We can integrate over from to with the same process used in integrating . Then can be written as follows

 P(~nrk)=1√2πσ2~nrkexp{−~nr 2k2σ2~nrk} (52)

where denotes . With the same steps can be found as

 P(~nik)=1√2πσ2~nikexp{−~ni 2k2σ2~nik} (53)

where denotes .

#### Probability distribution of P~nrk and P~nik

The total power contained in the kth Fourier component is where and . Now our goal is to derive the probability distribution of from Eq. (52) and Eq. (53). This can be done by a series of changing variable. At the first place the variables of the probability distributions were changed from and to and , respectively, and then changed from and to . Although the variables were changed, the integral of probability over entire region should be the same (and equal to one). Therefore, it is known that

 ∫∞0P(P~nrk) dP~nrk = ∫∞−∞P(~nrk) d~nrk (54) = 2∫∞0P(~nrk) d~nrk.

The range of is from zero to infinity since is positive. The second line is obtained from the symmetry of about zero. From Eq. (54) we know

 P(P~nrk)=2P(~nrk) ∣∣d~nrkdP~nrk∣∣ (55)

where is Jacobian. Since , the Jacobian is . Substituting Eq. (52) and the Jacobian into Eq. (55) we have

 P(P~nrk)=1√2πσ~nrk1√P~nrkexp{−P~nrk2σ2~nrk}. (56)

With the same steps we can derive

 P(P~nik)=1√2πσ~nik1√P~nikexp{−P~nik2σ2~n