Probing Gravitational Waves from Pulsars in Brans–Dicke Theory

: This paper comprises the theoretical background for the data analysis of gravitational waves (GWs) from spinning neutron stars in Brans–Dicke (BD) theory. Einstein’s general theory of relativity (GR) predicts only two tensor polarization states, but a generic metric theory of gravity can also possess scalar and vector polarization states. The BD theory attempts to modify the GR by varying gravitational constant G, and it has three polarization states. The ﬁrst two states are the same as in GR, and the third one is scalar polarization. We derive the response of a laser interferometric detector to the GW signal from a spinning neutron star in BD theory. We obtain a statistic based on the maximum likelihood principle to identify the signal in BD theory in the detector’s noise. This statistic generalizes the well known F -statistic used in the case of GR. We perform Monte Carlo simulations in Gaussian noise to test the detectability of the signal and the accuracy of estimation of its parameters.


Introduction
Isaac Newton believed that inertial forces such as centrifugal forces must arise from the acceleration with respect to "absolute space". On the other hand, Ernst Mach argued that inertial forces are more likely caused by acceleration with respect to the mass of the celestial bodies, and this idea is known as Mach's principle. This principle made a substantial contribution to the development of a new scalar-tensor theory of gravitation which is known as Brans-Dicke (BD) theory, proposed by Carl Brans and Robert Dicke [1]. The foundations of BD theory are built on the previous work of Pascual Jordan [2], as well as Markus Fierz [3], and sometimes this theory is also referred to as Jordan-Fierz-Brans-Dicke theory.
In the last 100 years, GR has successfully passed various experimental tests, proving itself to be the most promising theory of gravity. But then a natural question arises, why is there a quest to find an alternative theory of gravity? There are plenty of reasons that justify this question, and one of them is related to understanding the accelerated expansion of the universe. It is more reasonable to substitute GR by an alternative theory of gravity instead of introducing the concept of dark energy in Einstein's field equations to expound this expansion. This is because dark energy does not fit the standard model, and also its nature is not perceived by the laws of modern physics [4]. However, while exploring new models of gravity, one must discern that any new theory must obey some constraints imposed by field theory. GR is also a field theory and the following rules determine the credibility of any new theory of gravity: 1.
There should exist vibrational modes which are freely excitable in the absence of any source. They are also called degrees of freedom (DOF). So, a theory is classified by the number of vibrational modes and the spin of these modes.

2.
The next characteristic is the propagation of these modes in an empty space. If the mode is massless, we have Coulomb potential, which results in the long-range force.
On the other hand, a massive mode generates Yukawa potential, which is responsible for short-range force. 3.
Finally, one needs to consider the interaction of these modes with themselves as well as other fields.
We can verify the efficacy of these requirements in BD theory. This theory has three DOF. The first two of them are the same as in GR. They are known as tensor polarization or spin-2 massless graviton. The third polarization is known as scalar polarization or spin-0 graviton. All three polarizations are massless because gravity is a long-range force and gravitational waves travel with the speed of light. The interaction of modes can be understood by considering the Einstein equivalence principle (EEP). The BD theory satisfies the EEP, and hence it is a metric theory of gravity. Furthermore, in metric theories of gravity, only the metric (g) interacts with the matter. So, the scalar field mediates in such a way that it generates a gravitational field along with the matter. Once the gravitational field is set up, it produces a metric which in turn acts back on the matter as described by EEP, but the scalar field cannot respond to the matter directly [5]. For example, let us consider the Schwarzchild metric generated by a compact object. The other astronomical objects in the universe (planets, galaxies, etc.) will source this scalar field and hence deviate the resultant metric from the Schwarzchild metric. In this theory, Newton's constant of gravitation G is not constant anymore, but G is determined by the totality of the matter in the universe through an auxiliary field equation. This theory still has general coordinate invariance, but it has an additional degree of freedom φ(x) that determines the strength of gravity. In the presence of the scalar field, Einstein's equation modifies to (1) where f (φ) is a coupling constant analogous to the Newton's constant in the GR, and it depends on φ, which in turn depends on position (x). T (matter) µν is the energy-momentum tensor produced by the matter, and T (φ) µν is an energy-momentum tensor produced by the scalar field φ.
Although BD theory satisfies the EEP, it does not meet the strong equivalence principle (SEP). This is true because the value of coupling constant depends on φ, which itself depends on the position. On the other hand, a form of velocity dependence in local physics can also enter indirectly if the value of the scalar field changes with time. Then, the rate of variation of the coupling constant could depend on the velocity of the frame. It is worth noting that GR is the only theory which satisfies SEP, otherwise all metric theories violate this principle at some level [5]. Furthermore, for the sake of the complete introduction of the BD theory, we should mention that it is a purely dynamical theory which means that coupled partial differential equations govern its structure and evolution [5].
The plan of the paper is as follows. In Section 2, we present the current status of the experimental verification of Barns-Dicke theory. In Section 3, we derive the response of a laser interferometric detector to the gravitational-wave signal from an asymmetric rotating neutron star. In Section 4, we obtain the statistic to detect the signal in the noise of the detector and estimate its parameters. We present results of the Monte Carlo simulations to verify our detection and parameter estimation method.

Post-Newtonian and Parameterized Post-Newtonian Formalism
To understand the experimental tests of the BD theory in the weak-field limit, we need to digress a little bit and get an idea of the post-Newtonian (PN) expansion and parameterized post-Newtonian (PPN) framework. In the PN expansion, we expand the metric in the powers of v c . This approach is valid in the weak-field limit and slow-motion approximation (v << c). When the speed of gravity approaches infinity, we get back to the Newtonian potential. In the PPN formalism, we add parameters (coefficients) in front of the potentials obtained by PN expansion and add a few new potentials with their own parameters. This allows us to obtain a framework that encompasses a broad spectrum of alternative theories, and that can be used to calculate a wide range of testable phenomena. The only aspect that changes from one theory of gravity to the other is the numerical value of the various coefficients that appear in front of the potentials [5,6].
The value of γ is 1 in GR and 1+ω BD 2+ω BD in BD, whereas β is united in both theories. Furthermore, parameters α i and ζ i vanish in GR as well as BD, and hence both theories are conservative with no preferred frame effects. To test the BD theory, we need to measure the value of γ with the help of experiments. Once γ is obtained, we can infer the value of BD parameter ω BD . For a detailed analysis on PN and PPN formalism, refer to [5][6][7].

Experimental Tests of BD Theory in the Weak-Field Limit
The weak-field limit parameter is given by where G is the Newtonian gravitational constant, M is the characteristic mass scale of the phenomenon, R is the characteristic distance scale, and c is the speed of light. In the weak-field limit << 1. This is the case when we study a phenomenon in the solar system where < 10 −5 .
Measurement of γ: In the BD theory, the value of the parameter ω BD is constant. As ω BD tends to infinity, γ tends to unity, and the BD theory is no different from GR. Moreover, a higher value of ω BD diminishes the effect of the scalar field. Below we shall discuss the bounds obtained on the value of γ from different experimental tests. These methods are:

2.
The time delay of light: A radar signal sent across the solar system past the Sun to a planet or satellite and returned to the Earth suffers an additional non-Newtonian delay in its round-trip travel time. A significant improvement in the value of γ was reported in 2003 from Doppler tracking of Cassini spacecraft, while it was on its way to Saturn, with a result γ − 1 = (2.1 ± 2.3) × 10 −5 [13]. Using this result, one maintains that massless scalar-tensor theories must have ω > 40,000 to be compatible with this constraint [5].
Nordtvedt effect: Nordtvedt showed that many metric theories of gravity predict that massive bodies violate the weak equivalence principle-that is, fall with different accelerations depending on their gravitational self-energy [18]. This violation of the equivalence principle by massive bodies is known as the "Nordtvedt effect" and measured by the parameter η N . In GR, η N is zero, and a non-zero value of η N gives the deviation from GR. Various experiments have been carried out to measure this effect, but the one conducted by Eöt-Wash group is the most enhanced test. In this experiment, WEP was examined for laboratory bodies whose chemical compositions mimic that of the Earth and Moon, and they obtained the result |η N | = (4.4 ± 4.5) × 10 −4 [5].

Experimental Tests of BD Theory in Strong-Field Limit
The value of ≈ 1 defines the strong-field regime. It corresponds to a region in the vicinity of a neutron star or black hole. To test the theories of gravitation in a strong gravity regime, one can consider a binary system. It is known that the orbit of inspiralling binary will decay because of the emission of gravitational radiation. In GR, there is only quadrupolar emission, but scalar-tensor theories also predict dipolar emission. One can test the dipole radiation by observing the rate of change of the orbital period. Before we consider known binary systems, it is worth mentioning two propositions: • In a binary system with identical objects, dipole emission is suppressed. For example, a system with two neutron stars will result in weak dipole radiation, if there is any [5]. • Roger Penrose proposed that black holes in BD theory are identical to their GR counterparts. Motivated by this remark, Thorne and Dykla showed that during gravitational collapse to form a black hole, the BD scalar field is radiated away, in accordance with Price's theorem, leaving only its constant asymptotic value, and a GR black hole [5,19].
Although Hulse-Taylor binary (PSR 1913+16) confirmed the existence of gravitational radiation in GR, it could not provide a reliable test for dipole radiation because of the mass ratio of the system. However, the discovery of binary pulsar systems with a white dwarf companion, such as J1738+0333, J1141-6545 and J0348+0432, has made it possible to perform robust tests of the existence of dipole radiation. This is because such systems are necessarily asymmetrical. Already, significant bounds have been placed on dipole radiation using J1738+0333 and J1141-6545 [20,21]. Furthermore, constraints have been put on the dipole radiation using the event GW170817 [22].

Search for Nontensorial Gravitational Waves
Another approach to test the BD theory in the strong-field regime is to search for gravitational radiation from isolated rotating neutron stars [23]. One can search for scalar, vector or tensor polarizations in LIGO, Virgo, and KAGRA gravitational detector data [24][25][26]. In [27], the first search has been carried out in LIGO detector data from its first observational run for around 200 known pulsars without relying on any particular alternative theory of gravity. This investigation could not discover any GW signal, but it imposed upper limits for scalar, vector and tensor amplitudes [28].

Generation of Gravitational Waves in the Brans-Dicke Theory
The action in the BD theory is given by where and In the above expressions, g µν is the metric tensor, φ is the scalar field, ω BD is a parameter, ψ m are the matter fields and L m is the Lagrangian of matter fields.
The equations obtained after varying the action are Our goal is to solve Equations (6) and (7) simultaneously and retrieve the wave solution in the vacuum. This derivation is motivated by Section 13.4 in [6]. To do so, first we linearize the field equations as and then make an infinitesimal gauge transformation, along with introducing a tracereversed metric. Finally, we impose the Lorentz gauge to extract the independent degrees of freedom. After some algebra, the perturbation metric can be written as where h + corresponds to the plus polarization, h × corresponds to the cross polarization and h S denotes the scalar polarization. This form of the metric assumes that the wave travels in the z-direction, whereas distortion occurs in the xy-plane. In general, a metric theory can have six polarization states, but there are only three in BD theory, two tensor polarizations h + , h × and one scalar polarization h S . In the following, it will be convenient to use the parameter ζ instead of the parameter ω BD : To obtain approximate expressions for the two polarizations h + and h × , we use the standard quadrupole formalism (see Chapter 1.7 of [29]). Assuming that the wave propagates in the +z direction, we have the following formulae for h + and h × polarizations where G is the gravitational constant, c is the speed of light, t is the retarded time, r is the distance of the source,Q ij W is the second time derivative of the mass quadrupole moment in the wave frame. In the limit ω BD → ∞, i.e., ζ → 0 the above expressions for the two polarizations reduce to the expression in classical general relativity given by Equation (1.114) of [29].
The scalar polarization is derived in Chapter 13.5 of [6] and it is given by the following expression (see Equations (13.164) and (13.168a) in [6]) whereḊ i W is the first time derivative of the mass dipole moment in the wave frame and M is the mass monopole moment. In GR, there is no contribution from the mass monopole moment and dipole moment because of mass and linear momentum conservation. However, monopole and dipole radiation appear in BD theory because the scalar field φ does not satisfy a conservation law.

Rotating Neutron Star
If the star is perfectly symmetric, both dipole and quadrupole moments vanish, and there is no gravitational radiation. However, a real neutron star is not perfectly symmetrical because there may be some elastic deformations of its solid crust, which may be due to strong magnetic fields present in the neutron stars. We model the deformation as a mountain on a neutron star. This asymmetry leads both to quadrupolar and dipolar gravitational-wave emission.
The multipole expansion I <L> in terms of symmetric trace-free (STF) tensors is given by [6] where x = (x, y, z) are Cartesian coordinates, ρ is the density of the star and d 3 x is the volume element. < L > denotes the STF tensor and L represents a collection of l individual indices. When l = 1 and < L >=< i >, we get and when l = 2 and < L >=< ij > we have where r 2 is equal to x 2 + y 2 + z 2 and δ ij is the Kronecker delta function. Using the above relations, the i component of the dipole moment D i s in the star's frame is given by Similarly, the ij component of the quadrupole tensor Q ij s in the star's frame is calculated as Q ij For simplicity, we consider the star to be a sphere of radius a, and we assume that the size of the mountain is minimal compared to the radius of the star. Consequently, we can model the mountain as a point mass. Moreover, we assume that the Cartesian coordinates of the mountain are (a, 0, 0). Thus, the mass density ρ of the mountain is given by where δ is the Dirac delta function and m is the mass of the mountain. Using the equations above, we can easily determine the components of the dipole and quadrupole moments for this model. and Since our system is nearly spherical, it can also be expanded in terms of multipole moments I lm given by is the complex conjugate of spherical harmonics functions and d 3 x = r 2 sin θdrdθdφ. The density of the mountain in spherical coordinates can be written as For the quadrupolar emission, the most dominant moment is I 22 , and for the case of our model, it can be simplified as

Gravitational Wave Signal from a Rotating Neutron Star
To obtain the explicit equations for the polarization function h + , h × , and h S , we first transform the dipole vector D s and quadrupole matrix Q s from the source to wave frame. Following the construction in Chapter 2.5 of [29], we have and where S is the transformation matrix from the source frame to an inertial frame, and R is the transformation matrix from the inertial frame to the wave frame. The matrix R is given by where φ s (t) is the instantaneous rotational phase of the star. The matrix S has the form where ι is the angle between the angular momentum vector of the rotating neutron star and the direction along which the wave travels. We assume that the rotational phase φ s (t) is slowly varying, and it can be modelled by Taylor expansion as where φ o is a constant phase offset, f (k) 0 is the kth time derivative of the instantaneous frequency evaluated at t = 0 at and s is the number of spin-down parameters.
Using Equation (13) and neglecting the second order time derivatives of the phase φ(t), we obtain the following expressions for the two tensor polarizations where h o is the constant amplitude given by where quadrupole parameter Q is given by and f 0 is the spin frequency of the star. We notice that when ζ → 0, the expressions for h + and h × reduce to the polarizations in classical general relativity.
To calculate the scalar polarization, we make a simplifying assumption that the only non-vanishing component of the dipole moment in the star's frame is in the x-direction, like in the simple mountain model in Section 3.3, i.e., we have With this assumption, we obtain the following formula for the scalar polarization, neglecting the mass monopole's contribution (see [6]). The scalar polarization has contributions both from the dipole and quadrupole parts.
where constant amplitudes h d o and h q o are given by We see that the quadrupole part of the radiation gives in a signal at twice the spin frequency of the star, whereas dipole radiation results in a signal at once in spin frequency.

Response of the Interferometric Detector to a Gravitational-Wave Signal Forming a Rotating Neutron Star
To derive the detector's response function to gravitational-wave signal from a rotating neutron star, we follow the formalism presented in Section II of [30]. We obtain the response separately for the tensor and the scalar polarization. Thus, we divide the three-dimensional matrix H(t) of the spatial metric perturbation produced in the wave frame into tensor part H T (t) and the scalar part H S (t) defined by and Then, the response function h BD T to the tensor polarization and the response function h BD S to the scalar polarization are given by and In the above expressions, n 1 and n 2 denote the unit vectors parallel to the arm number 1 and 2 respectively. n 1 × n 2 points outwards from the surface of the Earth. The responses are defined as the difference between the wave-induced relative length changes of the two interferometer arms. We assume that the detector arms are orthogonal and choose the x-axis to be along the first interferometer arm and the y-axis to be along the second arm. Hence, we have n 1 = (1, 0, 0), n 2 = (0, 1, 0).
The matricesH T (t) andH S (t) are tensorial and scalar matrices respectively of the spatial metric perturbation produced by the wave in the proper reference frame of the detector. They are related to matrices in the wave frame by the following transformations The matrix M is a three-dimensional orthogonal matrix that transforms from the wave Cartesian coordinates in the wave frame to the Cartesian coordinates in the detector's proper reference frame. We have assumed that the wave travels in the +z-direction. The matrix M can be expressed as M 1 is the matrix of transformation from wave to celestial sphere frame coordinates given by where δ is the declination of the gravitational-wave source, α is its right ascension and ψ is the polarization angle. M 2 is the matrix of transformation from celestial coordinates to cardinal coordinates where λ is the latitude of the detector's site, Ω r is the rotational angular velocity of the Earth, φ r is the deterministic phase which defines the position of the Earth in its diurnal motion at time t = 0. (φ r + Ω r t) coincides with the local sidereal time of the detector's site. M 3 is the matrix of transformation from cardinal coordinates to detector of proper reference frame coordinates.
where the angle γ determines the orientation of the detector's arms with respect to local geographical directions. Tensor response. The tensor response has exactly the same form as the response of the detector in classical general relativity, except for the factor 1 − ζ in the constant amplitude h o (see Equation (33)). The detailed derivation is given in Section II of [30]. Here, we only summarize the basic formulae. The tensor response can be expressed as a linear combination of four time-dependent functions h i (t).
The phase φ(t) of the signal in the detector's frame is approximately given by where n 0 is the constant unit vector in the direction of the star in the SSB reference frame and r d is the position vector of the detector in that frame. The approximation that led to Equation (61) above is discussed in detail in Section IIB and Appendix A of [30]. where c(t) is the amplitude modulation function given by c(t) = 1 8 (−1 + 3 cos 2δ) cos 2 λ sin 2γ + 2 cos(α − φ r − tΩ r ) sin 2γ sin 2δ sin 2λ The scalar polarization can be written as a linear combination of four time-dependent functions h iS (t) where the four constant amplitudes A iS are given by We see that the scalar part of the response has contributions both from the dipole and quadrupole radiation. The functions h 1S (t) and h 2S (t) determine the response corresponding to the dipole contribution, whereas h 3S (t) and h 4S (t) define quadrupole contribution.
It is helpful to rewrite the above equations as a correction to the GR part. To do so, first, we rewrite constant amplitudes for tensor polarizations as where The combined tensor and scalar response can be written as which can be simplified to We notice that the last two terms in Equation (80) are corrections to the GR.
It is also advantageous to express the constant amplitude functions in terms of numerical astrophysical values.

Detection Statistic and Parameter Estimation
To obtain the detection statistic for the signal derived in the previous section, we use the maximum likelihood method described in Chapter 6 of [29]. We assume that the noise n(t) in the detector is an additive and zero-mean Gaussian process. Thus, when the signal s(t) is present, the data can be written as For the Gaussian case, the log-likelihood function is given by where, for the case of white noise with one-sided spectral density S o , the scalar product (·|·) is defined as Our signal is represented by a linear combination of time-dependent functions h l (t).
where the eight amplitudes depend on the six parameters h 0 , h d 0 , h q 0 , ψ, ι and φ 0 . The first six terms in the signal (90) originate from the quadrupolar radiation, whereas the last two are the dipolar contribution.
The signal s(t) can be expressed in the following compact form where A T is the transpose of A and · denotes matrix multiplication. A is a vector of constant amplitudes and h represents a vector of time-dependent functions h l (t) given.
By substituting Equation (91) into the likelihood function (88), we get where N is a vector with components and M is a matrix with components The maximum likelihood estimators are found by maximizing the Equation (93) with respect to amplitudes and this results inÂ By substituting these estimators back in Equation (93), we obtain the L-statistic. The final result is Two fundamental quantities for the data analysis method presented in this paper are the signal-to-noise ratio ρ and the Fisher matrix Γ. The signal-to-noise ratio is defined as and for the signal (91), it is given by The Fisher matrix for the case of Gaussian noise is defined by and for the case of signal (91), we simply have For the case of Gaussian noise, the signal-to-noise ratio determines the probability of detection of the signal. For a large signal-to-noise ratio, the inverse of the Fisher matrix approximates the covariance matrix of the estimators of the parameters. For the case of signal (91), which is a linear function of the amplitude parameters, one can show that the maximum likelihood estimators of the amplitudes given by Equation (97) are unbiased, and their covariance matrix is precisely equal to the inverse of the Fisher matrix.
We shall now obtain an explicit form of the L-statistic for the case of the gravitationalwave signal from a rotating neutron star in BD theory. First, we need the matrix M given by Equation (95). It can approximately be computed as where and A ≡ (a|a) , B ≡ (b|b) , C ≡ (a|b) , E ≡ (a|c), G ≡ (b|c) and H ≡ (c|c). To obtain the above approximation, we assume that the phase φ(t) has many oscillations over the observation time T 0 . Consequently, we have Similarly, we can split the column vector N into a quadrupolar part as well as a dipolar part and write it as where and Consequently, the likelihood ratio statistic L can be expressed as a sum of the two statistics F BD Q and F BD D for the quadrupole and dipole part of the signal, respectively. and The explicit expressions for F BD Q and F BD D are given by where and From Equation (97) the amplitude estimators take the following explicit form It should be noted that amplitude estimators (Â 1 ,Â 2 ,Â 3 ,Â 4 ) are related to the quadrupolar part in the tensor polarization, the estimators (Â 3S ,Â 4S ) correspond to the quadrupolar part in scalar polarization and (Â 1S ,Â 2S ) represent the dipole part in the scalar polarization. The quadrupole part has six amplitudes, four from tensor polarization and two from scalar polarization. On the other hand, the dipole part has only two amplitudes from the scalar polarization. The observations indicate that the ζ coefficient in the Brans-Dicke theory is small. Consequently, the quadrupole part from the scalar polarization is small compared to the quadrupole part from the tensor polarizations.
Furthermore, the amplitude h 0 is deviated from its GR counterpart by a factor of (1 − ζ), and to first approximation, we can neglect parameter ζ in h 0 .
Thus, we can safely approximate the quadrupole radiation in BD theory by quadrupole radiation in ordinary GR.
Thus, the statistic L in Equation (109) can be approximated by where F is the statistic for the quadrupole signal in Einstein's GR and D is the statistic for the dipole signal in BD theory. The F -statistic is given by and corresponding amplitude estimators arê The D-statistic is equal to F BD D given by Equation (114). F-statistic has been used extensively in searching GWs from known pulsars [31][32][33].
With the approximations above, we are left only with the six amplitude param- ψ, ι, φ 0 ) and the GW signal in Equation (90) simplifies to

Data Analysis Method
In this section, we shall present the data analysis method to detect the signal s(t) given by Equation (130) and estimate its parameters. The signal s(t) is a linear combination of six amplitudes, and these amplitudes are functions of five astrophysical parameters. To detect the signal, we adopt the likelihood ratio test which is equivalent to comparing the statistic L given by Equation (109) to a threshold. Once the signal is detected, we obtain the maximum likelihood estimatorsÂ l of the six amplitude parameters A l using the Equations (121), (122), (125)-(128). It is important to estimate the five astrophysical parameters h o , h d o , ψ, ι, and φ 0 . We estimate these five parameters by a least-squares (LS) fit determined by six amplitude parameters.
The estimators are obtained by minimizing the following function with respect to the astrophysical parameters: where Γ A l A l are components of the Fisher matrix given by Equation (102). The least-squares fit involves a non-linear minimization procedure for which we need the initial values for the five parameters (h o , h d o , ψ, ι, φ 0 ) with respect to which the LS function is minimized. We use an analytic solution for the six parameters in terms of the amplitude estimators for the initial values. Several such solutions exist, but we use the one given by the equations below. Firstly, we introduce the following auxiliary quantities.
The quadrupole amplitude h Q 0 is then given by The expressions for parameters ψ and ι are and We see that the quadrupole amplitude and angles ψ and ι can be calculated from the four amplitudes of the quadrupole signal. With the value of inclination angle ι obtained from Equation (138) above, we easily obtain expressions for dipole amplitude h d o and phase angle φ o from the two dipole amplitudes.

Monte Carlo Simulations
We have performed Monte Carlo simulations to test the performance of the data analysis method proposed in the previous Section. Each simulation consists of generating a signal given by Equation (130) and adding it to white, Gaussian noise. Then, for each simulation, we calculate the F and D statistics given by Equations (124) and (114) and we estimate the four quadrupole amplitudes and the two dipole amplitudes from Equations (121), (122), (125)-(128) respectively. We perform the simulations for a range of signal-to-noise ratios. For each value of the signal-to-noise ratio, we repeat our simulations 1000 times for different realizations of noise. Figure 1 shows the plots of the means and standard deviations of the D-statistic and the F -statistic against the signal-to-noise ratio of the signals added, whereas Figure 2 shows the sum of D and F statistics. We compare the simulated values with the theoretical predictions. In Gaussian noise, the 2 × F -statistic has the non-central χ 2 distribution with 4 degrees of freedom and non-centrality parameter equal to ρ 2 , whereas 2 × D-statistic has the non-central χ 2 distribution with 2 degrees of freedom. For any non-central χ 2 distribution with k independent degrees of freedom, the mean (µ) and standard deviations (σ) are given by and where λ is the non-centrality parameter. In our case, λ = ρ 2 , where ρ is the signal-to-noise ratio defined by Equation (99).  For the case of amplitude parameters, we calculate biases and variances of the parameter estimators and compare them with the theoretical values. For a parameter, the numerical bias is calculated as numerical value−theoretical value theoretical value × 100. The theoretical biases of parameter estimators are zero, whereas the variances are square roots of the diagonal elements of the inverse of the Fisher matrix given by Equation (102). Figures 3 and 4 show the biases and standard deviations of the estimators of the four quadrupole amplitudes.  In Figure 5, we present the bias and standard deviation of the estimators of two amplitudes from the dipole radiation. From our simulations, we also estimate the five astrophysical parameters θ = (h o , h d o , ψ, ι, φ o ) using the least squares procedure presented in the previous section. The theoretical estimates of the standard deviation can be obtained from the Fisher matrix Γ(θ) for the astrophysical parameters, which is given by the following formula where Γ(A) is the Fisher matrix for the amplitude parameters given by Equation (102) and J is the 6 × 5 Jacobi matrix with elements J lm = ∂A l ∂θ m (l ∈ [1, 6] and m ∈ [1, 5]). Figures 6-8 present the biases and standard deviations of the five astrophysical parameters. We have plotted the theoretical mean and standard deviation against the initial estimators and estimators obtained from the LS procedure.
The astrophysical parameters are related to amplitude parameters by non-linear transformation; consequently, we can expect that their estimators will be unbiased and minimum variance only asymptotically, for a large signal-to-noise ratio. Moreover, the estimators obtained from the LS procedure agree well with the theoretical predictions.

Conclusions
In this paper, we have derived the response of a laser interferometric detector to a GW signal from a non-axisymmetric rotating neutron star in the BD theory. We have also obtained the maximum likelihood ratio statistic to detect such a signal in the detector data and estimate its parameters. We find that in the Brans-Dicke theory, a non-negligible contribution to gravitational radiation from a rotating neutron star can come from the time-varying dipole moment. We have proposed a scheme whereby detecting both the quadrupole and the dipole radiation, we can estimate the GW amplitude of dipole radiation in the BD theory. We have performed Monte Carlo simulations in white Gaussian noise, which indicate that our scheme detects the signal and estimates its parameters with the theoretically predicted efficiency and accuracy.

Funding:
The work was supported by the Polish National Science Centre Grant No. 2017/26/M/ST9/00978.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.