long-distance high-fidelity continuous-variable quantum key distribution with non-Gaussian operations: An exact closed form solution

Khatereh Jafari ^1*, Mojtaba Golshani ^2, Alireza Bahrampour ^1

$^1$Department of Physics, Sharif University of Technology, Tehran, Iran

$^2$Faculty of Physics, Shahid Bahonar University of Kerman, Kerman, Iran

In this paper, we derive a closed-form expression for the output state of a CV-QKD protocol in the presence of zero-photon catalysis (ZPC) and a quantum scissor (QS). Then, based on this closed-form solution, we use a direct search algorithm to find the appropriate values of input state and QS parameters, which considerably enhance the range and the fidelity of a CV-QKD protocol. In the special case of pure loss channel, the largest range of the protocol is only $6.5\%$ less than the fundamental limit of repeaterless quantum communication. In addition, examination of the protocol for different values of excess noise, reveals that there is a trade-off between range and fidelity, and a high value of fidelity can be obtained at the cost of a slight reduction in protocol range.

Keywords: Continuous-Variable Quantum Key Distribution, Secure Key Rate, Zero Photon Catalysis, Quantum Scissor, Long-range, High-Fidelity

§ INTRODUCTION

Unlike classical cryptography, whose security is due to the high complexity of the mathematical problem, quantum key distribution (QKD) is a method that shares an unconditionally secure key between two users (traditionally called Alice and Bob) based on fundamental laws of quantum physics<cit.>. In the continuous variable QKD (CV-QKD), the information is encoded on the continuous characteristics of the light, e.g. its quadratures<cit.>.
Compared to its discrete counterpart<cit.>, CV-QKD has some advantages, namely compatibility with conventional technology, low cost, and high rate<cit.>. Nevertheless, this protocol has drawbacks such as short range and sensitivity to the excess noise<cit.>. To dispel these problems, different setups including quantum scissor (QS), photon addition, photon subtraction, and photon catalysis are proposed<cit.>.

QS is a non-deterministic operation that reduces the dimension of Hilbert space into two dimensions, and provides an arbitrary linear combination of vacuum and single photon states<cit.>. Based on the transmittance coefficients of scissor’s beam splitters, QS can act as a heralded noiseless linear amplifier for small amplitude input signals, and consequently can boost the performance of CV-QKD protocols at long distances<cit.>.

Photon Catalysis is a feasible approach to improve the performance of CV-QKD protocols<cit.>.
Zero photon catalysis (ZPC) is a special type of quantum catalysis that reduces the weight of higher-order Fock states and can be treated as a heralded noiseless attenuator<cit.>. This Gaussian operation, by reducing the modulation variance of the input signal, is able to increase the secure transmission distance of the CV-QKD protocols.

Quantum states become disturbed when passing through a quantum channel. In quantum communication, the main goal is to increase the secrete key rate and the secure transmission distance of the protocol, but it is also important to reduce errors caused by the channel. In fact, reducing the errors and its effect on the corruption of the transmitted state is important in the error correction process (post-processing stage) [1, 2, 3]. The fidelity is a measure of the closeness of Alice and Bob states, and hence can be used to quantify the corruption of states during the quantum channel <cit.>. In addition, some QKD protocols such as measurement-device-independent CV-QKD (MDI-CV-QKD) protocols are based on the teleportation of quantum states in which the fidelity plays a vital role [4, 5, 6, 7]. Therefore, it will be important to find the optimal parameters that can enhance the fidelity without significant reduction of the protocol range.

Recently, we offered a CV-QKD protocol equipped with both ZPC and QS, and demonstrated that this suggested scheme outstandingly enhances the fidelity  even at long distances <cit.>.
However, in that work, due to the structure of the protocol, it was not possible to obtain a closed-form solution for the output state of the quantum channel. Therefore,
due to the time-consuming numerical calculations, the optimal values of transmittance coefficients of QS's beam splitters, to achieve the maximum secure range, have not been obtained.
In this paper, we are looking to optimize the parameters to obtain the highest achievable values of the protocol range and fidelity.
In order to realize this, some parts of the previous protocol have been changed. In Ref. <cit.>, we used a modified QS that truncates the Hilbert space to third-order Fock states<cit.>. However, in this research, the standard QS is adopted.
This modification considerably simplifies both the experimental setup (because of the simpler sources and detectors of the standard QS) and computational complexity (optimization of one beam splitter instead of two).
Moreover, here, we assume Alice uses Gaussian modulation in contrast to the previous work that discrete (quadrature phase-shift keying (QPSK)) modulation was adopted<cit.>. In the entanglement-assisted protocol, this means Alice uses a two-mode squeezed vacuum (TMSV) state and performs heterodyne measurement on her mode<cit.>. In this case, exact closed-form expression can be obtained for the output state of the protocol, which significantly reduces the required computational time of the optimization process. Moreover, As we will see in the next section, performing of the ZPC operation on one mode of the TMSV state (in contrast to the QPSK state) is equivalent to reducing the modulation variance of that state. Therefore, it is possible to further simplify the protocol and eliminate the ZPC operation, and as an alternative optimize the modulation variance of the TMSV state.

In this paper, we derive an exact expression for the output state of a protocol in the presence of standard QS at the Bob side. Then, based on this result, the optimal modulation variance and transmission coefficient of QS are numerically calculated to achieve the maximum secure transmission distance. As will be seen, this optimization allows us to find the parameters that considerably  boost the secure range of the protocol, as well as, leads to a high fidelity between input and output states.
These findings
highlight the importance of the optimization in CV-QKD protocols.

This research is organized in four sections. The protocol description and theoretical consideration is described in Section 2. In section 3, we present the optimization of the protocol and discuss about its numerical results. Finally, conclusion is presented in section 4.

§ PROTOCOL DESCRIPTION AND THEORETICAL CONSIDERATION

The schematic of the protocol is illustrated in Fig. <ref>, which uses a ZPC at Alice’s side and a QS at Bob’s side. In general, the Alice state can be written in the number state basis as,
\ket{\psi}_{A A'}=\sum_{n_A,n_A'}{\Gamma(n_A,n_A')\ket{n_A}_{A}\ket{n_A'}_{A'}},
where $\Gamma(n_A,n_A')$ satisfies the normalization condition,
\sum_{n_A,n_A'}{\left|\Gamma(n_A,n_A')\right|^2}=1.
Alice performs a ZPC operation on mode $A'$, which converts $\Gamma(n_A,n_A')$ to
\frac{t_z^{n_A'}}{\sqrt{P_z}}\Gamma(n_A,n_A'),
where $t_z$ denotes the transmission coefficient of the beam splitter and the normalization factor
is the success probability of ZPC operation<cit.>. Since $\left| t_z\right|<1$, ZPC operation acts as an attenuator, reduces the contribution of higher-order number states by the factor $\left|t_z\right|^{2n_A'}$ and therefore degrades the impact of channel loss on the input state. Conversely, we use a QS at the end of the channel, which acts as a noiseless amplifier, compensates this weight reduction, and upgrades the performance of the protocol.
In this paper, we assume Alice uses a TMSV state
\ket{\psi}_{A A'}=\sqrt{1-\lambda_A^2}\sum_{n_A=0}^{\infty}{\lambda_A^{n_A}\ket{n_A}_{A}\ket{n_A}_{A'}},
which is a special case with
$\Gamma(n_A,n_A')=\sqrt{1-\lambda_A^2}\,\lambda_A^{n_A} \delta_{n_A,n_A'}$
($\delta_{n_A,n_A'}$ is the Kronecker delta).
ZPC operation on mode $A'$ changes this state to the state
\ket{\psi}_{A A''}=\sqrt{\frac{1-\lambda_A^2}{P_z}}\sum_{n_A=0}^{\infty}{t_z^{n_A}\lambda_A^{n_A}\ket{n_A}_{A}\ket{n_A}_{A''}},
with $P_z=\frac{1-\lambda_A^2}{1-t_z^2\lambda_A^2}$. For TMSV state, ZPC operation exclusively reduces $\lambda_A$ by the factor $t_z$, i.e., diminishes the quadrature variance from $\frac{1+\lambda_A^2}{1-\lambda_A^2}$ to $\frac{1+t_z^2\lambda_A^2}{1-t_z^2\lambda_A^2}$ <cit.>.
Accordingly, for this special input state, one can remove the ZPC operation (i.e., set $t_z=1$) and apply the optimization process to $\lambda_A$ (see Fig. <ref>).
It is important to note that this conclusion is not true for a general input state (for example, QPSK modulation scheme of Ref. <cit.>). Therefore, in general, both ZPC operation and QS are required to noiselessly diminish the loss of quantum channel and improve the performance of the protocol <cit.>.
Alice sends mode $A''$ through the quantum channel and performs heterodyne measurement on her mode $A$. In this paper, we assume that two parties use a thermal-loss channel with transitivity $t_c^2$ and an excess noise $\varepsilon$. This kind of channel can be modeled with a beam splitter with the transmission coefficient $t_c=10^{-0.01 L}$, where $L$ is the channel length (in km unit) <cit.>. This assumption corresponds to the case of a Gaussian attack enabled by an entangling cloner, in which one input of the beam splitter is Alice's mode $A''$ and the other is mode $E'_2$ of Eve's TMSV state
\ket{\psi}_{E_1 E'_2}=\sqrt{1-\lambda_E^2}\sum_{n_E=0}^{\infty}{\lambda_E^{n_E}\ket{n_E}_{E_1}\ket{n_E}_{E'_2}}.
Here,  $\lambda_E$ characterizes the excess noise of the channel, $\varepsilon=\frac{2\lambda_E^2}{1-\lambda_E^2}$<cit.>.
In this case, the equivalent excess noise at the transmitter side (input of the channel) is given by
$\varepsilon_{tm}=(1-t_c^2)\varepsilon/t_c^2$, while
$\varepsilon_{rec}=(1-t_c^2)\varepsilon$ represents the amount of excess noise at
the receiver side (output of the channel) <cit.>. Here, same as Refs. <cit.> we use
$\varepsilon$, and investigate this proposed protocol for different values of  $\varepsilon$.
It is important to note that  $\varepsilon_{rec}$ is always less than $\varepsilon$ and approaches it in the limit of large values of propagation distances $L$, but, $\varepsilon_{tm}$ is larger than $\varepsilon$ for $L\gtrsim 15\mathrm{km}$, and gets larger by increasing the propagation distances $L$. As we will see later, in the range of values considered in this study (see Table. <ref>), $\varepsilon$ is very close to $\varepsilon_{rec}$, and therefore, the excess noise $\varepsilon$ can be considered as the noise at the output of the channel.

Our protocol uses a QS that is a non-Gaussian operation. Therefore, Gaussian attack is not necessarily an optimal attack for this protocol. Finding optimal eavesdropper strategy may be very challenging and is not the intention of current research. Therefore, we restrict the eavesdropper to a Gaussian attack.  We will talk about this issue again later in the paper.

Entanglement-based scheme of the protocol.

Simplified version of the protocol.

The states
$\ket{\psi}_{A A''}$ and $\ket{\psi}_{E_1 E'_2} $ after mixing in the beam splitter are converted to<cit.>

\begin{eqnarray}
\label{eqOutputstate}
\ket{\psi}_{\scaleto{A B E_1 E_2}{5.5pt}}=\mathlarger{\sum}_{\substack{n_A,n_B \\ n_E}}{\Theta_{n_A,n_B}^{n_E} \scaleto{\ket{n_A}_{A}\ket{n_B}_{B}\ket{n_E}_{E_1}\ket{n_A+n_E-n_B}_{E_2}}{10.5pt}}.
\end{eqnarray}

\begin{eqnarray*}
\label{eqTheta1}
\begin{split}
&\Theta_{n_A,n_B}^{n_E}=\sqrt{(1-\lambda_A^2)(1-\lambda_E^2)}\lambda_A^{n_A}\lambda_E^{n_E}(-t_c)^{n_E-n_B} r_c^{n_A+n_B}\binom{n_E}{n_B}
\\
&\times \sqrt{\frac{n_B!(n_A+n_E-n_B)!}{n_A! n_E!}} {}_2\mathlarger{F}_1\left(-n_A,-n_B;1+n_E-n_B;-\frac{t_c^2}{r_c^2} \right),
\end{split}
\end{eqnarray*}

for $0 \le n_B \le n_E$, and

\begin{eqnarray*}
\label{eqTheta2}
\begin{split}
&\Theta_{n_A,n_B}^{n_E}=\sqrt{(1-\lambda_A^2)(1-\lambda_E^2)}\lambda_A^{n_A}\lambda_E^{n_E} t_c^{n_B-n_E} r_c^{n_A-n_B+2n_E}\binom{n_A}{n_B-n_E}
\\
&\times \sqrt{\frac{n_B!(n_A+n_E-n_B)!}{n_A! n_E!}} {}_2\mathlarger{F}_1\left(n_B-n_A-n_E,-n_E;1+n_B-n_E;-\frac{t_c^2}{r_c^2} \right),
\end{split}
\end{eqnarray*}

for $n_E <n_B \le n_A+n_E$.
Here, $r_c=\sqrt{1-t_c^2}$ is the reflection coefficient of the beam splitter and ${}_2F_1(a,b;c;z)$ is the Gaussian or ordinary hypergeometric function <cit.>.

The noise and loss of the channel limit the secure range of the protocol. To overcome this limitation, we assume Bob is equipped with a QS. QS is a heralded noiseless linear amplifier that can compensate the loss of quantum channel and boosts the secure transmission distance of the protocol<cit.>. Standard QS cuts off all Fock states with two or more photons and prepares a superposition of zero and one photon states <cit.>. Setting $n_B=0$ (zero photon) and $n_B=1$ (one photon) in $\Theta_{n_A,n_B}^{n_E}$, after straightforward calculations, one can find,

\begin{eqnarray}
\label{eqOutputstateScissor}
\begin{split}
\ket{\psi}_{\scaleto{A B_s E_1 E_2}{5.5pt}}=&\mathlarger{\sum}_{n_A,n_E}{V_{n_A}^{n_E} \scaleto{\ket{n_A}_{A}\ket{0}_{B_s}\ket{n_E}_{E_1}\ket{n_A+n_E}_{E_2}}{10.5pt}}
\\
+&\mathlarger{\sum}_{n_A,n_E}{W_{n_A}^{n_E} \scaleto{\ket{n_A}_{A}\ket{1}_{B_s}\ket{n_E}_{E_1}\ket{n_A+n_E-1}_{E_2}}{10.5pt}}.
\end{split}
\end{eqnarray}

\begin{eqnarray*}
\label{eqW}
V_{n_A}^{n_E}=t_s(-t_c)^{n_E}r_c^{n_A}\lambda_A^{n_A}\lambda_E^{n_E}\sqrt{\frac{(1-\lambda_A^2)(1-\lambda_E^2)(n_A+n_E)!}{2 n_A! n_E! P_Q}},
\end{eqnarray*}

\begin{eqnarray*}
\label{eqW}
\begin{split}
W_{n_A}^{n_E}=&gt_s(-t_c)^{n_E-1}r_c^{n_A-1}(n_E r_c^2-n_At_c^2) \lambda_A^{n_A}\lambda_E^{n_E}
\\
&\times \sqrt{\frac{(1-\lambda_A^2)(1-\lambda_E^2)(n_A+n_E-1)!}{2 n_A! n_E! P_Q}}.
\end{split}
\end{eqnarray*}

Here, $g=\sqrt{\frac{1-t_s^2}{t_s^2}}$ and $t_s$ are the amplification gain and transmission coefficient of QS, respectively <cit.>.
Moreover, the normalization factor

\begin{eqnarray}
\label{eqSuccessScissor}
\end{eqnarray}

is the success probability of QS, i.e., the probability that its detectors count zero and one photon at horizontal and vertical output arms, respectively (see the right bottom panel of Fig. <ref>).

The Alice-Bob state $\rho_{AB_s}$ is the partial trance of Eq. (<ref>) over eavesdropper modes $E_1$ and $E_2$, and is given by

\begin{eqnarray}
\label{eqRhoABs}
\begin{split}
\rho_{AB_s}=&\mathlarger{\sum}_{n_A=0}^{\infty}\left(\sigma^{00}_{n_A}\ket{n_A,0}\bra{n_A,0}+\sigma^{11}_{n_A}\ket{n_A,1}\bra{n_A,1} \right. \\
& \left. \sigma^{01}_{n_A}\ket{n_A,0}\bra{n_A+1,1}+\sigma^{10}_{n_A}\ket{n_A+1,1}\bra{n_A,0} \right),
\end{split}
\end{eqnarray}

\begin{eqnarray*}
\label{eqRhoABs00}
\sigma^{00}_{n_A}=\frac{t_s^2(1-\lambda_A^2)(1-\lambda_E^2)(r_c\lambda_A)^{2n_A}}{2P_Q(1-t_c^2\lambda_E^2)^{n_A+1}},
\end{eqnarray*}

\begin{eqnarray*}
\label{eqRhoABs11}
\sigma^{11}_{n_A}=\frac{g^2t_s^2(1-\lambda_A^2)(1-\lambda_E^2)(r_c\lambda_A)^{2n_A}\left[ n_At_c^2(1-\lambda_E^2)^2+r_c^4\lambda_E^2 \right]}{2P_Q r_c^2(1-t_c^2\lambda_E^2)^{n_A+2}},
\end{eqnarray*}

\begin{eqnarray*}
\label{eqRhoABs01}
\sigma^{01}_{n_A}=\sigma^{10}_{n_A}=\frac{gt_s^2t_c\sqrt{n_A+1}(1-\lambda_A^2)(1-\lambda_E^2)^2 r_c^{2n_A}\lambda_A^{2n_A+1}}{2P_Q (1-t_c^2\lambda_E^2)^{n_A+2}}.
\end{eqnarray*}

Eq. (<ref>) gives the analytical expression for output state of the protocol. Using this equation, the Alice $\rho_A$ and Bob $\rho_{B_s}$ states are given by

\begin{eqnarray}
\label{eqRhoA}
\rho_{A}=\mathlarger{\sum}_{n_A=0}^{\infty}{\gamma^{(A)}_{n_A}\ket{n_A}\bra{n_A} },
\end{eqnarray}

\begin{eqnarray}
\label{eqRhoBs}
\rho_{B_s}=\gamma^{(B)}_0\ket{0}\bra{0}+\gamma^{(B)}_1\ket{1}\bra{1},
\end{eqnarray}

\begin{eqnarray*}
\label{eqRhoAnA}
\begin{split}
\gamma^{(A)}_{n_A}&=\frac{(1-\lambda_A^2)(1-\lambda_E^2)(r_c\lambda_A)^{2n_A}}{2P_Q r_c^2(1-t_c^2\lambda_E^2)^{n_A+2}} \\
&\times \left[ n_A g^2 t_s^2 t_c^2(1-\lambda_E^2)^2+r_c^4\lambda_E^2 +r_c^2 t_s^2(1-\lambda_E^2) \right],
\end{split}
\end{eqnarray*}

\begin{eqnarray*}
\label{eqRhoB0}
\gamma^{(B)}_{0}=\frac{t_s^2(1-\lambda_A^2)(1-\lambda_E^2)}{2P_Q(1-r_c^2\lambda_A^2-t_c^2\lambda_E^2)},
\end{eqnarray*}

\begin{eqnarray*}
\label{eqRhoB1}
\gamma^{(B)}_{1}=\frac{g^2 t_s^2 (1-\lambda_A^2)(1-\lambda_E^2)[ t_c^2\lambda_A^2+\lambda_E^2(r_c^2-\lambda_A^2)]}{2P_Q(1-r_c^2\lambda_A^2-t_c^2\lambda_E^2)^2}.
\end{eqnarray*}

Here, same as <cit.>, we can investigate the non-Gaussian behavior of our protocol.
To this end, we focus on the distribution of homodyne measurement results on one of Bob’s quadratures, for example

The probability distribution of obtaining $x$ after measuring $\hat{x}_{B_s}$ is given by <cit.>

\begin{eqnarray}
\label{eqfxB}
f(x)=\bra{x}\rho_{B_s}\ket{x}=\sqrt{\frac{2}{\pi}}e^{-2x^2}\left( \gamma^{(B)}_{0} +4 \gamma^{(B)}_{1} x^2 \right),
\end{eqnarray}

where $\hat{x}_{B_s}\ket{x}=x\ket{x}$.
Accordingly, the Gaussian and non-Gaussian terms of $f(x)$ are proportional to $\gamma^{(B)}_{0}$ and $\gamma^{(B)}_{1}$, respectively. Therefore, the output state is approximately Gaussian as long as
$\gamma^{(B)}_{1}/\gamma^{(B)}_{0} \ll 1$. Considering the fact that $\gamma^{(B)}_{1}$ decreases by decreasing of $t_c$ (i.e., increasing of propagation length $L$) and takes its largest value at $t_c=1$, one can simply show that the output state is close to the Gaussian state if $g^2\lambda_A^2 \ll 1$. If this condition is met, the Gaussian attack can be considered as the optimal attack.

To demonstrate the benefit of optimized input state and QS, we consider the reverse reconciliation scheme, in which Alice corrects her raw key in order to match Bob's key. Therefore, under collective attack, the secret key rate is given by
$K=P_Q\left( \beta S[\rho_A] - S[\rho_{AB_s}] \right)$, where $\beta$ is the reconciliation efficiency and $S[\rho]=-tr\left( \rho \log_2 \rho \right)$ is the Von Neumann entropy <cit.>.
In this equation, $P_Q$ reduces the value of the key rate because the information is calculated only when QS is successful.
This is why QS diminishes the secure key rate at lower values of the propagation distance<cit.>. However, at larger values of the propagation distance where the signal's strength gets weaker, QS acts as a noiseless amplifier, distills the entanglement, and improves the secure range of the protocol.
Although the secret key rate and the secure propagation distance are the key parameters, however,  fidelity is also an important parameter and high fidelity can simplify the post-processing stage of the QKD protocol <cit.>. Therefore, in addition to the key rate, we will also discuss about the fidelity,
F=\left( tr \left( \sqrt{\sqrt{\rho_{B_s}} \rho_{A} \sqrt{\rho_{B_s}} } \right) \right)^2,
which is a popular measure of the closeness of the Alice $A$ and Bob $B_s$ modes <cit.>. Using equations (<ref>) and (<ref>), the fidelity can be written as

\begin{eqnarray}
\label{eqFidelity}
F=\left( \sqrt{\gamma^{(A)}_{0}\gamma^{(B)}_{0}}+\sqrt{\gamma^{(A)}_{1}\gamma^{(B)}_{1}} \right)^2.
\end{eqnarray}

As we will see in the next section, optimization of input state (i.e., $\lambda_A$) and QS (i.e., $t_s$) results in long-distance high-fidelity QKD.

§ NUMERICAL RESULTS AND DISCUSSIONS

Our goal is to find the optimal value of $\lambda_A$ and $t_s$ that maximize the protocol range.
TMSV state with $\lambda_A=1$ corresponds to maximally entangled (Einstein-Podolsky-Rosen) state, which is unnormalizable and unphysical<cit.>. In fact, achieving values of $\lambda_A$ close to unity is practically difficult.

In addition, QS acts as an amplifier ($g>1$) as long as $t_s^2<0.5$.
Therefore, we restrict the search domain to $\lambda_A \in (0,0.9)$ and $t_s \in (0,0.7)$.
Exact equation for the output state of the channel allows us to use a direct search algorithm for this two-parameter optimization problem. To this end, we divide the intervals $0<\lambda_A \le 0.9$ and $0<t_s \le 0.7$ into 900 and 700 partitions, respectively, i.e. the step size $\Delta \lambda_A=\Delta t_s=0.001$. Then, for each value of $\lambda_A$ and $t_s$, we use the bisection method to find the longest propagation distance, i.e., protocol range $R$, that has a key rate greater than $K_{min}$. In this paper, we set the minimum acceptable secrete key rate (bits per pulse) as $K_{min}=10^{-6}$, while in literature smaller values of $10^{-8}$ and $10^{-10}$ are also considered <cit.>.
Finally, by comparing the protocol range for all $\lambda_A$ and $t_s$ values, the optimal value of the protocol range can be achieved.

(Color online) Secret key rate (per pulse) versus transmission distance for $\lambda_A=0.815$, $t_s=0.041$ and $\epsilon \simeq 0$. For comparison, the PLOB bound is also shown. The inset shows fidelity versus transmission distance.

(Color online) The green area shows the regions with the range larger than $265km$. The yellow area shows the regions with the fidelity greater than $0.94$ ($\epsilon \simeq 0$).

We first consider a pure loss (noiseless) channel. Therefore, we set $\varepsilon=0$ ($\lambda_E \to 0$).
Our results show that, in this case, the maximum achievable protocol range is about $R=288 km$, which takes place at $\lambda_A=0.815$, at $t_s=0.041$.
Fig. <ref> shows the secret key rate versus transmission distance. For comparison, the PLOB bound (named after. Pirandola, Laurenza, Ottaviani, and Banchi) is also shown <cit.>. The PLOB bound $K_{PLOB}=-log_2(1-t_c^2)$ is the fundamental limit of repeaterless quantum communication and reaches the value of $K_{min}=10^{-6}$ at $L \simeq 308 km$. This indicates that the maximum achievable range of our optimized protocol is only $6.5\%$ less than the PLOB range.
It is important to note that, QS reduces the value of the key rate at short distances because the information is calculated for the post-selected data only when QS is successful. However, at large distances, QS acts as an amplifier and boosts the protocol range <cit.>.

The inset of Fig. <ref> shows the fidelity $F$ of Eq. <ref> for $\lambda_A=0.815$ and $t_s=0.041$. According to this figure, the fidelity is about $0.36$  at the protocol range ($R=288km$), which is a small value. Therefore, we search for other values of $\lambda_A$ and $t_s$ with high range $R$ that have considerable fidelity $F$ at that range. Hereafter, the term 'fidelity' ($F$) would be used for the value of fidelity at the protocol range, i.e. at the propagation distance $L=R$.
The green and yellow areas of Fig. <ref> show the regions in $\lambda_A-t_s$ plane with $R>265km$ and $F>0.94$, respectively.
From this figure, one can see that the protocol has high range for larger values of $\lambda_A$ and smaller values of $t_s$. This means, when the input state has a large number of photons (large value of $\lambda_A$), or QS has a high amplification gain (small value of $t_s$), in the absence of noise, the output state is less altered by the loss, and therefore, the range of protocol increases. In addition, at smaller values of $\lambda_A$, in which the weight of higher-order number states of Alice mode are small, the fidelity between Alice (Eq. <ref>) and Bob (Eq. <ref>) states increases.

(Color online) The $(\lambda_A,t_s)$ points with a range larger than $265km$ and a fidelity greater than $0.94$ ($\epsilon \simeq 0$).

(Color online) Secret key rate (per pulse) versus transmission distance for $\lambda_A=0.815$, $t_s=0.041$ and $\lambda_A=0.359$, $t_s=0.004$, for the case with $\epsilon\simeq 0$. The inset compares fidelities.

(Color online) The green area shows the regions with the range larger than $110km$. The yellow area shows the regions with the fidelity greater than $0.93$ ($\epsilon = 0.001$).

(Color online) The $(\lambda_A,t_s)$ points with a range larger than $110km$ and a fidelity greater than $0.93$ ($\epsilon = 0.001$).

(Color online) Secret key rate (per pulse) versus transmission distance for $\lambda_A=0.787$, $t_s=0.038$ and $\lambda_A=0.413$, $t_s=0.087$, for the case with $\epsilon= 0.001$. The inset compares fidelities.

Range and fidelity of the protocol for different values of the excess noise.

$\;\;\;\;$     $\epsilon\simeq 0\;\;\;\;$     $\epsilon=0.001\;\;\;\;$     $\epsilon=0.005\;\;\;\;$     $\epsilon=0.01\;\;\;\;$     $\epsilon=0.05\;\;\;\;$

$R (km)$     288     127     92     78     47

$F$     0.36     0.27     0.26     0.23     0.25

$\lambda_A$     0.815     0.787     0.798     0.819     0.849

$t_s$     0.041     0.038     0.085     0.121     0.246

$\varepsilon_{tm}$     0     0.346     0.341     0.353     0.385

$\varepsilon_{rec}$     0     0.001     0.005     0.010     0.044

$R (km)$     265     110     76     62     35

$F$     0.94     0.93     0.93     0.93     0.92

$\lambda_A$     0.359     0.413     0.413     0.423     0.473

$t_s$     0.004     0.087     0.182     0.257     0.464

$\varepsilon_{tm}$     0     0.157     0.161     0.164     0.201

$\varepsilon_{rec}$     0     0.001     0.005     0.009     0.040

$\Delta R (\%)$     -8     -13     -17     -21     -26

$\Delta F (\%)$     161     244     258     304     268

In the intersection of green and yellow regions of Fig. <ref>, the protocol has proper range and fidelity. The $(\lambda_A,t_s)$ points with $R>265km$ and $F>0.94$ are show in Fig. <ref>. In this figure, at a constant value of $t_s$, the fidelity decreases by increasing $\lambda_A$, while the range becomes greater or remains constant by increment of $\lambda_A$.
These results show that, in the case of $\epsilon \simeq 0$, by fine tuning of the modulation variance of the input state, and the gain of the QS, one can achieve a protocol with a range around $265 km$ and a fidelity of about $0.94$.

Fig. <ref> and its inset compare, respectively, the secret key rate and fidelity, for $\lambda_A=0.359$, $t_s=0.004$ (the first case with range $R\simeq265km$ and fidelity $F\simeq0.94$) and $\lambda_A=0.815$, $t_s=0.041$ (the second case with maximum range $R=288km$ and $F=0.36$).
As can be seen, there is a trade-off between protocol range (or key rate) and fidelity. However, it is obvious that although the first case has the range which is about $8\%$ less than the second case, but the fidelity of the first case is approximately $161\%$ larger than the second one. Therefore, a much better fidelity can be obtained at the cost of a slight reduction in protocol range.

Now, we repeat our analysis fo a nonzero value of the excess noise $\epsilon=0.001$. We find that, in this situation, when $\lambda_A=0.787$ and $t_s=0.038$, the maximum range of $R=127km$ is obtained, which is about $56\%$ less than the case with $\epsilon \simeq 0$ (note that the excess noise $\epsilon=0.001$ at distance $L=127km$ is equivalent to $\varepsilon_{tm}\simeq 0.346$ and $\varepsilon_{rec} \simeq 0.001$). In fact, the excess noise reduces the correlation between Alice and Bob modes and brings down the secure transmission distance of the protocol <cit.>. In addition, the fidelity is small and about $F=0.27$. As a result, we search for other values of $(\lambda_A,t_s)$ with considerable range and fidelity.
In Fig. <ref>, the green and yellow areas correspond to the regions with $R>110km$ and $F>0.93$, respectively.
The general results of this figure are similar to the results of Fig. <ref>, except that the fidelity is not large at low values of $t_s$. A small value of $t_s$ corresponds to a large value of gain $g$ for QS. Therefore, the reduction of the fidelity for low values of $t_s$ and $\lambda_A$ can be related to the amplification of the noise (instead of the low amplitude signal) at the Bob side, which diminishes the similarity between Alice and Bob states. The $(\lambda_A,t_s)$ points which lie in the intersection of green and yellow regions of Fig. <ref> are shown in Fig. <ref>, and have acceptable range, $R \simeq  110 km$, and fidelity, $F \simeq 0.93$ (the excess noise $\epsilon=0.001$ at distance $L=110km$ is equivalent to $\varepsilon_{tm}\simeq 0.157$ and $\varepsilon_{rec} \simeq 0.001$). In this case, the range is about $13\%$ less than the maximum achievable range (i.e., 127km) while the fidelity is approximately $244\%$ larger.
In Fig. <ref> and its inset, we plot the secrete key rate and fidelity versus propagation distance for $\lambda_A=0.413$, $t_s=0.087$ and compare them with the case of optimal range, i.e., $\lambda_A=0.787$, $t_s=0.038$. This figure indicates that both range (not the key rate) and fidelity can not be optimal at the same time, and by decreasing one, the other can be increased.

Finally, Table <ref> summarizes the range and fidelity of the protocol, for five different values of the excess noise $\varepsilon$. In this table, the first six rows correspond to the case with maximum achievable range, while the second six rows belong to a case with acceptable range and fidelity. In addition, the last two rows show the variation (in percent) of protocol range and fidelity. Please note that for completeness we also report the equivalent excess noise at the transmitter $\varepsilon_{tm}$ and receiver $\varepsilon_{rec}$ sides (for maximum achievable length $R$).

This table again confirms that we can obtain much higher fidelity (from $161\%$ to $304\%$), if we give up a smaller amount of the protocol range (from $8\%$ to $26\%$). For all values of the noise, this can be done by reduction of $\lambda_A$, i.e., the modulation variance of the input state. The higher-order number states have
a higher probability of being affected by channel loss. Accordingly, if, by decreasing $\lambda_A$, the weight of higher-order Fock states is reduced, and, at the end of the channel, is properly repaid by a QS (as an amplifier), the fidelity of the protocol can be effectively improved <cit.>.
Moreover, according to Table <ref>, at nonzero values of the excess noise, the amplification gain should be decreased (i.e., $t_s$ should be increased) to obtain higher fidelity. This action can be attributed to the prevention of noise amplification in the system.

§ CONCLUSIONS

In summary, by restricting eavesdropper to a Gaussian collective attack, we derived an exact expression for Alice and Bob's states and their fidelity, in a CV-QKD protocol which consists of a ZPC at the Alice side and a QS at the Bob side. For the case of Gaussian modulation, we realized that one can simply remove the ZPC operation and as an alternative reduce the modulation variance (or $\lambda_A$) of the input state. Based on these equations, a direct search algorithm was used to obtain the optimal value of $\lambda_A$ and $t_s$ that maximize the protocol range.
For the case of a pure loss channel ($\epsilon \simeq 0$), the results show that the maximum achievable range is only $6.5\%$ less than the fundamental limit of repeaterless quantum communication. On the other hand, our investigation for both noiseless and noisy channels revealed that, although at the optimal values of $\lambda_A$ and $t_s$, the protocol range is maximum, however, the fidelity is small and less than $0.36$. Therefore, other values of $\lambda_A$ and $t_s$ that give acceptable range and fidelity were obtained. We find that there is a trade-off between range and fidelity, and much higher fidelity can be obtained if someone ignores a smaller amount of the protocol range. This modification can be done by reduction of the modulation variance of the input state, and in the case of a noisy channel, by a decrease in the gain of QS.

Declaration of Competing Interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Data availability

Data will be made available on request.

R. Wolf, “Quantum key distribution protocols,” in Quantum Key Distribution, (Springer, 2021).

S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani et al., “Advances in quantum
cryptography,” Adv. optics photonics 12, 1012–1236 (2020).

N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, “Quantum cryptography,” Rev. modern physics 74, 145 (2002).

T. C. Ralph, “Continuous variable quantum cryptography,” Phys. Rev. A 61, 010303 (1999).

F. Grosshans and P. Grangier, “Continuous variable quantum cryptography using coherent states,” Phys. Rev. Lett. 88, 057902 (2002).

F. Grosshans, G. Van Assche, J. Wenger, R. Brouri, N. J. Cerf, and P. Grangier, “Quantum key distribution using gaussian-modulated coherent states,”
Nature 421, 238–241 (2003).

C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84,
621 (2012).

S. L. Braunstein and P. Van Loock, “Quantum information with continuous variables,” Rev. modern physics 77, 513 (2005).

G. Bennett, Charles H.; Brassard, “Quantum cryptography: Public key distribution and coin tossing,” Theor. Comput. Sci. 560, 7–11 (2014).

A. K. Ekert, “Quantum cryptography based on bell’s theorem,” Phys. Rev. Lett. 67, 661–663 (1991).

C. H. Bennett, F. Bessette, G. Brassard, L. Salvail, and J. Smolin, “Experimental quantum cryptography,” J. cryptology 5, 3–28 (1992).

T. Hirano, H. Yamanaka, M. Ashikaga, T. Konishi, and R. Namiki, “Quantum cryptography using pulsed homodyne detection,” Phys. review A 68, 042331

H. Yonezawa, S. L. Braunstein, and A. Furusawa, “Experimental demonstration of quantum teleportation of broadband squeezing,” Phys. review letters 99,
110503 (2007).

H.-X. Ma, P. Huang, D.-Y. Bai, T. Wang, S.-Y. Wang, W.-S. Bao, and G.-H. Zeng, “Long-distance continuous-variable measurement-device-independent
quantum key distribution with discrete modulation,” Phys. Rev. A 99, 022322 (2019).

L. Hu, M. Al-amri, Z. Liao, and M. S. Zubairy, “Continuous-variable quantum key distribution with non-gaussian operations,” Phys. Rev. A 102, 012608

Y. Guo, W. Ye, H. Zhong, and Q. Liao, “Continuous-variable quantum key distribution with non-gaussian quantum catalysis,” Phys. Rev. A 99, 032327

W. Ye, H. Zhong, Q. Liao, D. Huang, L. Hu, and Y. Guo, “Improvement of self-referenced continuous-variable quantum key distribution with quantum
photon catalysis,” Opt. express 27, 17186–17198 (2019).

R. Blandino, A. Leverrier, M. Barbieri, J. Etesse, P. Grangier, and R. Tualle-Brouri, “Improving the maximum transmission distance of continuous-variable
quantum key distribution using a noiseless amplifier,” Phys. Rev. A 86, 012327 (2012).

D. T. Pegg, L. S. Phillips, and S. M. Barnett, “Optical state truncation by projection synthesis,” Phys. Rev. Lett. 81, 1604 (1998).

T. C. Ralph and A. P. Lund, “Nondeterministic noiseless linear amplification of quantum systems,” arXiv: Quantum Phys. 1110, 155–160 (2009).

M. Ghalaii, C. Ottaviani, R. Kumar, S. Pirandola, and M. Razavi, “Long-distance continuous-variable quantum key distribution with quantum scissors,”
IEEE J. Sel. Top. Quantum Electron. 26, 1–12 (2020).

M. Ghalaii, C. Ottaviani, R. Kumar, S. Pirandola, and M. Razavi, “Discrete-modulation continuous-variable quantum key distribution enhanced by quantum
scissors,” IEEE J. on Sel. Areas Commun. 38, 506–516 (2020).

T. J. Bartley and I. A. Walmsley, “Directly comparing entanglement-enhancing non-gaussian operations,” New J. Phys. 17, 023038 (2015).

[1]
R. Matsumoto and M. Hagiwara, “A survey of quantum error correction,” IEICE Trans. on Fundam. Electron. Commun. Comput. Sci. 104, 1654–1664

[2]
S. I. Salim, A. Quaium, S. Chellappan, and A. A. Al Islam, “Enhancing fidelity of quantum cryptography using maximally entangled qubits,” in GLOBECOM
2020-2020 IEEE Global Communications Conference, (IEEE, 2020), pp. 1–6.

[3]
F. Bouchard, D. England, P. J. Bustard, K. Heshami, and B. Sussman, “Quantum communication with ultrafast time-bin qubits,” PRX Quantum 3, 010332

[4]
H.-K. Lo, M. Curty, and B. Qi, “Measurement-device-independent quantum key distribution,” Phys. Rev. Lett. 108(13), 130503 (2012).

[5]
Z. Li, Y.-C. Zhang, F. Xu, X. Peng, and H. Guo, “Continuous-variable measurement-device-independent quantum key distribution,” Phys. Rev. A 89(5), 052301 (2014).

[6]
C. Weedbrook, “Continuous-variable quantum key distribution with entanglement in the middle,” Phys. Rev. A 87(2),
022308 (2013).

[7]
S. Pirandola, C. Ottaviani, G. Spedalieri, C. Weedbrook, S. Braunstein, S. Lloyd, T. Gehring, C. Jacobsen, and U. Andersen, “High-rate measurement-device-independent quantum cryptography,” Nat. Photonics 9(6), 397–402 (2015).

K. Jafari, M. Golshani, and A. Bahrampour, “High-fidelity continuous-variable quantum key distribution via the proper combination of zero-photon catalysis
and quantum scissors,” JOSA B 40, 661–666 (2023).

K. Jafari, M. Golshani, and A. Bahrampour, “Discrete-modulation measurement-device-independent continuous-variable quantum key distribution with a
quantum scissor: exact non-gaussian calculation.” Opt. Express 30, 11400–11423 (2022).

A. Leverrier and P. Grangier, “Unconditional security proof of long-distance continuous-variable quantum key distribution with discrete modulation,” Phys.
review letters 102, 180504 (2009).

I. B. Djordjevic, Physical-layer security and quantum key distribution (Springer, 2019).

M. Mičuda, I. Straka, M.Miková, M. Dušek, N. J. Cerf, J. Fiurášek, and M. Ježek, “Noiseless loss suppression in quantum optical communication,” Phys. Rev. Lett. 109, 180503 (2012).

G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical methods for physicists: a comprehensive guide (Academic press, 2011).

M. S. Winnel, J. J. Guanzon, N. Hosseinidehaj, and T. C. Ralph, “Overcoming the repeaterless bound in continuous-variable quantum communication
without quantum memories,” arXiv preprint arXiv:2105.03586 (2021).

S. Pirandola, R. Laurenza, C. Ottaviani, and L. Banchi, “Fundamental limits of repeaterless quantum communications,” Nat. communications 8, 15043


iint TXFiint

# Gravitational Waves from the Pulsar Magnetosphere

Ioannis Contopoulos1, Demosthenes Kazanas2, Demetrios B. Papadopoulos3
1 Research Center for Astronomy and Applied Mathematics, Academy of Athens, 4
Soranou Efessiou Str., Athens 11527, Greece
2 Astrophysics Science Division, NASA/ Goddard Space Flight Center, Greenbelt
Maryland, 20771, USA
3 Aristotle University of Thessaloniki, Department of Physics, Section of
Astrophysics, Astronomy and Theoretical Mechanics, Thessaloniki 54124, Greece
Email<EMAIL_ADDRESS>demos.kazanas<EMAIL_ADDRESS>
(Accepted -. Received -; in original form -)

###### Abstract

We investigate the generation of gravitational waves from the rotation of an
orthogonal pulsar magnetosphere in flat space time. We calculate the first
order metric perturbation due to the rotation of the non-axisymmetric
distribution of electromagnetic energy density around the central star. We
show that gravitational waves from a strong magnetic field pulsar right after
its formation within a distance of 1 kpc may be detectable with the new
generation of gravitational wave detectors.

###### keywords:

gravitational waves; MHD; pulsars

††pagerange: Gravitational Waves from the Pulsar
Magnetosphere–LABEL:lastpage††pubyear: 2023

##  1 Introduction

Gravitational waves are ripples in the fabric of spacetime caused by the
acceleration of massive objects [Rainich 1925, Misner et al. 1973, Ivanov
2008]. In particular, rotating neutron stars, the remnants of massive stars
that have undergone supernova explosions, are potential sources of continuous
gravitational wave emission when irregularities or deformations of the star’s
spherical form are present. As the star’s rotation rate remains steady, the
emitted gravitational waves will possess consistent frequency and amplitude,
akin to a singer holding a sustained note. The study of these emissions, if
detected, can provide valuable insights into the properties of neutron stars.
Rotating neutron stars are also strongly magnetized, and are thus surrounded
by a magnetosphere that is known to emit electromagnetic (Poynting) radiation,
and is often observed to also emit periodic electromagnetic pulses at the
period of their rotation (milliseconds to tens of seconds).

In the present work we will show that a non-axisymmetric rotating
magnetosphere can also emit gravitational waves. We will ignore non-spherical
deformations of the neutron star, and we will investigate continuous
gravitational waves that are due solely to the imprint of the anisotropic
rotating distribution of electromagnetic energy on the metric. As a first
approximation, we will ignore the effect of the star’s own gravity and will
consider the gravitational perturbation of a rotating electromagnetic field on
flat space time. Notice that already in 1925, Rainich gave the necessary and
sufficient conditions for a gravitational field to originate from a no-null
electromagnetic field

We will assume that the metric $g_{\mu\nu}$ can be expanded in powers of a
small parameter $\varepsilon\ll 1$

$g_{\mu\nu}^{(exact)}=g_{\mu\nu}+\varepsilon h_{\mu\nu}+O(\varepsilon^{2})$
(1)

(see [Thorne et al. 1982, Wald 1984, Maggiore 2008, Bambi et al. 2022]) and
the stress-energy tensor $T_{\mu\nu}$ can be similarly expanded as

$T_{\mu\nu}=\varepsilon T_{\mu\nu}^{(1)}+O(\varepsilon^{2})$ (2)

(Greek indices take values $0,1,2,3$, Latin indices $1,2,3$). To expand the
Einstein equations $G_{\mu\nu}=8\pi T_{\mu\nu}$ in powers of $\varepsilon$, we
expand the Einstein tensor in powers of the perturbation $h$

$G_{\mu\nu}(g+h)=G_{\mu\nu}(g)+G_{\mu\nu}^{(1)}(h)+O(h^{2})$ (3)

where, $G_{\mu\nu}^{(1)}$ is obtained from the exact Riemann tensor [Wald
1984] (see Ch.4, page 74). For a vacuum background keeping only the first
order terms in our approximations we obtain

$G_{\mu\nu}^{(1)}=(g_{\mu}^{\alpha}g_{\nu}^{\beta}-\frac{1}{2}g_{\mu\nu}g^{\alpha\beta})R_{\alpha\beta}^{(1)}$
(4)

where

$R_{\alpha\beta}^{(1)}=-\frac{1}{2}(\Box
h_{\alpha\beta}+2g^{\mu\lambda}g^{\nu\kappa}R_{\alpha\lambda\beta\kappa}h_{\mu\nu}-2g^{\mu\sigma}\bar{h}_{\mu(\alpha;\sigma\beta)})$
(5)

Here, $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the d’Alembertian in
any background and

$\bar{h}_{\mu\nu}\equiv
h_{\mu\nu}-\frac{1}{2}g_{\mu\nu}g^{\alpha\beta}h_{\alpha\beta}$ (6)

A semicolon and $\nabla$ both denote the covariant derivative, compatible with
the flat space time $g_{\mu\nu}$. Noting that
$\bar{h}=g^{\mu\nu}\bar{h}_{\mu\nu}=h-2h=-h$, eq. (6) can be inverted to give

$h_{\mu\nu}=\bar{h}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\bar{h}$ (7)

Substitution of eqs. (1) and (2) into Einstein equations and equating only
powers of $\varepsilon$, we end up with the equation

$G_{\mu\nu}^{(1)}=8\pi T_{\mu\nu}^{(1)}$ (8)

This perturbative expansion comes with the freedom to perform gauge
transformations

$h_{\mu\nu}\rightarrow h_{\mu\nu}+{\cal L}_{\xi}g_{\mu\nu}$ (9)

where ${\cal L}_{\xi}$ is a Lie derivative, and $\xi^{\alpha}$ are the
components of a freely chosen vector field.

In self-force theory, this freedom is commonly used to impose the Lorenz gauge
condition

$\nabla_{\alpha}\bar{h}^{\alpha\beta}=0$ (10)

in which case the linearized Einstein tensor simplifies to

$G_{\mu\nu}^{(1)}=-\frac{1}{2}(\Box\bar{h}_{\mu\nu}+2g^{\alpha
k}g^{\lambda\beta}R_{\mu\kappa\nu\lambda}\bar{h}_{\alpha\beta})$ (11)

eq. (8) with the aid of eqs. (6) and (11) gives

$-\frac{1}{2}(\Box\bar{h}_{\mu\nu}+2g^{\alpha
k}g^{\lambda\beta}R_{\mu\kappa\nu\lambda}\bar{h}_{\alpha\beta})=8\pi
T_{\mu\nu}^{(1)}$ (12)

We will solve eq. (12) in a $(t,r,\theta,\phi)$ coordinate system in a
Minkowski background where $g_{\mu\nu}=\eta_{\mu\nu}$. To do so, we choose
$\xi$ such that the trace $\bar{h}=0$ and thus $h_{\mu\nu}=\bar{h}_{\mu\nu}$
[Maggiore 2008]. Next, we will consider the metric perturbations $h_{\mu\nu}$
that are due to a rotating non-axisymmetric strong electromagnetic field.

We recently came across the most interesting paper Nazari & Roshan (2020).
These authors investigate the same problem, their method is, however,
different. They implement a post-Newtonian (PN) approximation to the
energy–momentum tensor of the magnetized fluid, and a PN expansion of the
gravitational potential in the wave zone. We instead work directly with the
perturbation of the stress-energy tensor and study its effect on the
background metric. We will compare our respective results in the Discussion
section.

##  2 A Rotating Magnetosphere

The Minkowski space time in spherical coordinate system $(t,r,\theta,\phi)$ in
geometric units where $G=c=1$ is

$g_{\mu\nu}=\eta_{\mu\nu}=diag(1,-1,-r^{2},-r^{2}\sin^{2}{\theta})$ (13)

and eq. (12) reads [Wald 1984, Maggiore 2008, Bambi et al. 2022]

$\Box h_{\mu\nu}=-16\pi T_{\mu\nu}^{(1)}\ .$ (14)

We will assume that the energy-momentum tensor is due entirely to an
electromagneti field and is thus given by

$T^{(1)}_{\mu\nu}=-\frac{1}{4\pi}(F_{\mu\sigma}F_{\nu\rho}g^{\rho\sigma}-\frac{1}{4}F_{\rho\sigma}F^{\rho\sigma}g_{\mu\nu})\
,$ (15)

where the Faraday tensor $F^{\mu\nu}$ can be constructed from the spatial
vectors ${\bf E}$, ${\bf B}$ of the electric and magnetic field.

The electromagnetic field that we will investigate is the central magnetic
dipole of an orthogonally rotating neutron star, otherwise known as an
orthogonal pulsar. The rotating spherical source at the origin has size
$r_{0}\sim 10\ {\rm km}$ and angular velocity of rotation $\omega_{0}\equiv
2\pi/P$, where $P$ is the period of rotation. We will only consider the
stellar exterior, assuming that the mass and magnetic field distribution in
the stellar interior is more or less spherically symmetric and does not
generate gravitational waves. We will consider a rotating dipole field of the
form

${\bf B}=\frac{3\hat{\bf r}(\hat{\bf r}\cdot{\bf m})-{\bf m}}{r^{3}}$ (16)

where ${\bf m}=m(\cos({\omega_{0}t}),\sin{(\omega_{0}t)},0)$ in cartesian
coordinates $(x,y,z)$, $z$ is the axis of rotation, and $\hat{\bf r}$ is the
unit vector along the radial direction that connects a point in space with the
origin. Here, $m\equiv B_{0}r_{0}^{3}$ is the stellar magnetic dipole, and
$B_{0}$ is the equatorial value of the stellar magnetic field. In this
configuration, the covariant components of the magnetic field are

$\displaystyle B^{r}$ $\displaystyle=$ $\displaystyle
2B_{0}\frac{r_{0}^{3}}{r^{3}}\sin\theta\ \cos(\omega_{0}t-\phi)$ (17)
$\displaystyle B^{\theta}$ $\displaystyle=$ $\displaystyle-
B_{0}\frac{r_{0}^{3}}{r^{3}}\cos\theta\ \cos(\omega_{0}t-\phi)$ (18)
$\displaystyle B^{\phi}$ $\displaystyle=$ $\displaystyle-
B_{0}\frac{r_{0}^{3}}{r^{3}}\sin(\omega_{0}t-\phi)\ ,$ (19)

Assuming ideal conditions in the magnetosphere, it is straightforward to show
(see Appendix A) that, in steady state magnetospheric rotation,

${\bf E}=-\omega_{0}r\sin\theta\ \hat{\phi}\times{\bf B}\ ,$ (20)

which yields

$\displaystyle E^{r}$ $\displaystyle=$
$\displaystyle-\frac{B_{0}r_{0}^{3}\omega_{0}}{r^{2}}\sin\theta\cos\theta\
\cos(\omega_{0}t-\phi)$ (21) $\displaystyle E^{\theta}$ $\displaystyle=$
$\displaystyle-\frac{2B_{0}r_{0}^{3}\omega_{0}}{r^{2}}\sin^{2}\theta\
\cos(\omega_{0}t-\phi)$ (22) $\displaystyle E^{\phi}$ $\displaystyle=$
$\displaystyle 0\ .$ (23)

Here, $\hat{\phi}$ is the unit vector in the azimuthal direction. The reader
can directly check that the above sets of eqs. (17-19) and (21-23) satisfy the
ideal Maxwell’s equations which in geometric units read

$\displaystyle\nabla\cdot{\bf E}=4\pi\rho_{e}\ \ \nabla\cdot{\bf B}=0$
$\displaystyle\frac{\partial{\bf E}}{\partial t}=\nabla\times{\bf B}-4\pi{\bf
J}\approx 0\ \ \frac{\partial{\bf B}}{\partial t}=\nabla\times{\bf E}\approx
0\ .$ (24)

Here $\rho_{e}$ is the electric charge density, and ${\bf J}$ is the
magnetospheric electric current.

From the study of pulsar magnetospheres it is known that the above dipole
structure of the magnetic field changes significantly beyond the cylindrical
distance of the so-called light cylinder, namely $R_{\rm lc}\equiv
c/\omega_{0}$. It is interesting that, beyond that distance, the magnetic
field attains a configuration that may be described as a ‘wound monopole’,
with

$\displaystyle|B^{r}|(r\sin\theta\mathrel{\raise
2.15277pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}R_{\rm lc})$
$\displaystyle\propto$ $\displaystyle\frac{1}{r^{2}}$ (25) $\displaystyle
B^{\theta}(r\sin\theta\mathrel{\raise 2.15277pt\hbox{$>$}\mkern-14.0mu\lower
2.58334pt\hbox{$\sim$}}R_{\rm lc})$ $\displaystyle\approx$ $\displaystyle 0$
(26) $\displaystyle|B^{\phi}|(r\sin\theta\mathrel{\raise
2.15277pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}R_{\rm lc})$
$\displaystyle\propto$ $\displaystyle|B^{r}|\frac{r\sin\theta}{R_{\rm
lc}}\propto\frac{\sin\theta}{r}$ (27)

[Bogovalov 1999]. It is obvious from the above eqs. (25-27) that beyond the
light cylinder, the distribution of $B^{2}$ is axisymmetric to a good
approximation (i.e. it does not depend on the azimuthal angle $\phi$), and is,
therefore, incapable of generating gravitation radiation. In contrast, the
distribution of the magnetic and electric field well inside the light cylinder
as expressed via eqs. (17-19) and eqs. (21-23) respectively is non-
axisymmetric, and therefore, it generates gravitation radiation as we will see
below. In what follows, we will further assume that the source rotates slowly
enough so that $r_{0}\omega_{0}\ll 1$. This is equivalent to say that the size
$r_{0}$ of the source is much smaller than the light cylinder radius $R_{\rm
lc}\equiv c/\omega_{0}$, or that the period of rotation $P$ is much larger
than $2\pi r_{0}/c=0.2$ ms.

The Faraday tensor $F_{\mu\nu}$ can be constructed directly from the
expressions in eqs. (17)-(23), and using the equation [Thorne et al. 1982]

$F^{ab}=u^{a}E^{b}-u^{b}E^{a}=\varepsilon^{abcd}u_{c}B_{d}$ (28)

where $u^{a}=(1,0,0,0)$ and $\varepsilon^{abcd}$ is the Levi-Civita symbol we
obtain

$F_{\mu\nu}=\left(\begin{array}[]{cccc}0&-E^{r}&-rE^{\theta}&0\\\
E^{r}&0&rB^{\phi}&-r\sin{\theta}B^{\theta}\\\
rE^{\theta}&-rB^{\phi}&0&r^{2}\sin{\theta}B^{r}\\\
0&r\sin{\theta}B^{\theta}&-r^{2}\sin{\theta}B^{r}&0\end{array}\right)$ (29)

From eq. (15) and (29) we compute $T_{\mu\nu}$. We emphasize once again that
we work in the slow rotation approximation where second and higher order
powers of $r_{0}\omega_{0}$ are ignored. Therefore, according to eq. (15),

$\displaystyle T^{(1)}_{tt}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{6}r_{0}^{6}}{8\pi
r^{6}}(1+3\sin^{2}{\theta}\cos^{2}{w})$ $\displaystyle T^{(1)}_{tr}$
$\displaystyle=$ $\displaystyle-\frac{B_{0}^{2}r_{0}^{6}\omega_{0}}{4\pi
r^{5}}\sin{\theta}\sin{(2w)}$ $\displaystyle T^{(1)}_{t\theta}$
$\displaystyle=$ $\displaystyle\frac{B_{0}^{2}r_{0}^{6}\omega_{0}}{16\pi
r^{4}}\sin{(2\theta)}\sin{(2w)}$ $\displaystyle T^{(1)}_{t\phi}$
$\displaystyle=$ $\displaystyle-\frac{B_{0}^{2}r_{0}^{6}\omega_{0}}{8\pi
r^{4}}\sin{(2\theta)}(1+3\sin^{2}{\theta}\cos^{2}{w})$ (30) $\displaystyle
T^{(1)}_{rr}$ $\displaystyle=$ $\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{8\pi
r^{6}}(1-5\sin^{2}{\theta}\cos^{2}{w})$ $\displaystyle T^{(1)}_{r\theta}$
$\displaystyle=$ $\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{4\pi
r^{5}}\sin{(2\theta)}\cos^{2}{w}$ $\displaystyle T^{(1)}_{r\phi}$
$\displaystyle=$ $\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{4\pi
r^{5}}\sin^{2}{\theta}\sin{(2w)}$ (31) $\displaystyle T^{(1)}_{\theta\theta}$
$\displaystyle=$ $\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{8\pi
r^{4}}(\cos{(2w)}-5\sin^{2}{\theta}\cos^{2}{w})$ $\displaystyle
T^{(1)}_{\theta\phi}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{16\pi r^{4}}\sin{(2\theta)}\cos{(2w)}$
$\displaystyle T^{(1)}_{\phi\phi}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{8\pi
r^{4}}\sin^{2}{\theta}(\cos{(2w)}+3\sin^{2}{\theta}\cos^{2}{w})$ (32)

where $w\equiv\omega_{0}t-\phi$. The reader may check that the above energy
momentum tensor has zero divergence in the limit
$r_{0}^{2}\omega_{0}^{2}\simeq 0$ and zero trace , i.e. $Q^{\mu}\equiv
T^{(1)}\ _{;\nu}{}^{\mu\nu}=0$ and $T^{(1)}\ _{\mu}{}^{\mu}=0$. Because of
this result and since in flat space-time the d’Alembertian $\Box$ commutes
with $\partial_{\mu}$, the Lorentz gauge condition reads

$\partial^{\nu}\Box
h_{\mu\nu}=\partial^{\nu}T^{(1)}_{\mu\nu}=0\Rightarrow\Box\partial^{\nu}h_{\mu\nu}=0\Rightarrow\partial^{\nu}h_{\mu\nu}=0$
(33)

Following the standard procedure appearing in the literature, [Wald 1984,
Maggiore 2008, Bambi et al. 2022], the so-called TT-gauge may be applied, and
the only independent components of the gravitational wave are
$h_{\theta\theta}$, and $h_{\theta\phi}$. Thus, eqs. (14) yield two
independent inhomogeneous wave equations

$\Box h_{\theta\phi}=-16\pi T_{\theta\phi}^{(1)}$ (34) $\Box
h_{\theta\theta}=-16\pi T_{\theta\theta}^{(1)}$ (35)

where $T^{(1)}_{\theta\theta},T^{(1)}_{\theta\phi}$ are given by eqs. (2). The
source $T^{(1)}_{\phi\phi}$ yields one more component of the gravitational
wave that depends on $h_{\phi\phi}$ and satisfies the inhomogeneous wave
equation $\Box h_{\phi\phi}=-16\pi T_{\phi\phi}^{(1)}$, namely
$h_{\phi\phi}=-\sin^{2}{\theta}\ h_{\theta\theta}$. We will next solve the
above equations to determine $h_{\theta\phi}$ and $h_{\theta\theta}$.

##  3 The Solution

The inhomogeneous wave eqs. (34)-(35) are of the form

$\Box\Psi=4\pi f\>,$ (36)

where $\Box$ is the d’Alembertian in any background, and $f\equiv
f(t,r,\theta,\phi)$ is a source due to the electromagnetic field. In flat
space time, eq. (36) reads

$\displaystyle\frac{\partial^{2}\Psi}{\partial
t^{2}}-\frac{\partial^{2}\Psi}{\partial
r^{2}}-\frac{2}{r}\frac{\partial\Psi}{\partial r}$
$\displaystyle-\frac{1}{r^{2}\sin{\theta}}\left(\frac{\partial}{\partial\theta}(\sin{\theta}\frac{\partial\Psi}{\partial\theta})+\frac{1}{\sin{\theta}}\frac{\partial^{2}\Psi}{\partial\phi^{2}}\right)=4\pi
f$ (37)

Subsequently, we decompose $\Psi$ and the source term $f$ of eq. (3) in
spherical harmonics as

$\Psi=\sum_{lm}Y_{lm}(\theta,\phi)\int_{-\infty}^{\infty}d\omega e^{-i\omega
t}R_{lm}(r_{0},r;\omega)$ (38)
$f=\sum_{lm}Y_{lm}(\theta,\phi)\int_{-\infty}^{\infty}d\omega e^{-i\omega
t}f_{lm}(r_{0},r;\omega)\ ,$ (39)

where

$R_{lm}(r_{0},r;\omega)=-i\omega H_{l}^{(1)}(\omega
r)\int_{r_{0}}^{\infty}j_{l}(\omega\rho)f_{lm}(\rho;\omega)d\rho$ (40)
$f_{lm}(r;\omega)=-2\int_{-\infty}^{\infty}dte^{i\omega t}\int
f(t,r,\theta,\phi)Y_{lm}^{*}(\theta,\phi)d\Omega\ .$ (41)

Here $H_{l}^{(1)}$ and $j_{l}$ are the Hankel and Bessel functions of the
first kind respectively. Following the usual Fourier expansion method we end
up with the inhomogeneous radial equation

$\frac{d}{dr}\left(r^{2}\frac{dR_{lm}}{dr}\right)+(r^{2}\omega^{2}-l(l+1))R_{lm}=f_{lm}$
(42)

A star on the spherical harmonics and other quantities on our equations means
complex conjugate. Upon the consideration of eq. (40), eq. (38) becomes

$\Psi(t,r,\theta,\phi)=-i\sum_{l,m}Y_{lm}(\theta,\phi)\int_{-\infty}^{\infty}\omega\
d\omega e^{-i\omega t}H_{l}^{(1)}(\omega
r)\int_{r_{0}}^{\infty}j_{l}(\omega\rho)f_{lm}(\rho;\omega)d\rho$ (43)

Following the above decompositions, we are now ready to obtain the solutions
of the inhomogeneous eqs. (34) and (35).

###  3.1 The solution of eq. (34)

We start with eq. (34), which in the limit $r_{0}^{2}\omega_{0}^{2}\simeq 0$
and the aid of eqs. (2) yields

$\Box
h_{\theta\phi}=-\frac{B_{0}^{2}r_{0}^{6}}{r^{4}}\sin{(2\theta)}\cos{(2w)},~{}~{}\mbox{~{}with~{}}~{}~{}w=\omega_{0}t-\phi$
(44)

As it is apparent, its source term is

$f=-\frac{B_{0}^{2}r_{0}^{6}}{4\pi r^{4}}\sin{(2\theta)}\cos{(2w)}$ (45)

Eqs. (41) and (45) yield (see Appendix C)

$\displaystyle f_{lm}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}\sqrt{\frac{\pi}{5}}}{3r^{4}}\\{\delta(\omega-2\omega_{0})\int
d\Omega Y_{lm}(\theta,\phi)\frac{d}{d\theta}(Y_{22}(\theta,\phi))$ (46)
$\displaystyle+\ \delta(\omega+2\omega_{0})\int d\Omega
Y_{lm}(\theta,\phi)\frac{d}{d\theta}(Y_{22}^{*}(\theta,\phi))\\}$

Furthermore, from eq. (43) and eq. (94), we find the solution of eq. (44)

$\displaystyle
h_{\theta\phi}(t,r,\theta,\phi;r_{0},\omega_{0})=i\frac{4\sqrt{5\pi}B_{0}^{2}r_{0}^{6}}{3r}\cdot$
$\displaystyle\sum_{lm}\frac{2^{l}\omega_{0}^{l}}{(2l+1)!!}\frac{R_{\rm
lc}^{(l-3)}-r_{0}^{(l-3)}}{l-3}Y_{lm}(\theta,\phi)S_{lm}\delta_{m2}e^{\frac{i}{2}[4t\omega_{0}-4\omega_{0}r+l\pi]}$

where $S_{lm}$ is given by eq. (Appendix C) in Appendix C.

We found that the double sum in eq. (3.1) for $l=0,1$, and $m=-l$ to $m=+l\ $
becomes zero. In Appendix C, we calculated the term for $l=2$ and found that
it is also zero, since $S_{22}=0$. For $l=3$ and $m=2$, we have
$S_{32}=\frac{3}{128}\sqrt{210\pi}$, and using the de l’Hospital rule, we find
that $(R_{\rm lc}^{(l-3)}-r_{0}^{(l-3)})/(l-3)$ approaches $\ln(r_{0}/R_{\rm
lc})$ and thus

$h_{\theta\phi}=3.1\times 10^{-3}\ \ln(\frac{r_{0}}{R_{\rm
lc}})\frac{\sqrt{\pi}B_{0}^{2}r_{0}^{6}\omega_{0}^{3}}{r}\cos(2[\omega_{0}(t-r)+\phi])\cos\theta\sin^{2}\theta$
(48)

which is of order $(r_{0}\omega_{0})^{3}$.

###  3.2 The solution of eq. (35)

This equation

$\Box h_{\theta\theta}=-16\pi T_{\theta\theta}^{(1)}\equiv 4\pi
f=\frac{2B_{0}^{2}r_{0}^{6}}{r^{4}}(\cos{(2w)}-5\sin^{2}{\theta}\cos^{2}{w})\
.$ (49)

Eq (41) yields

$f_{lm}=-\frac{B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int(\cos{(2w)}-5\sin^{2}{\theta}\cos^{2}{w})Y_{lm}^{*}(\theta,\phi)d\Omega\\\
$

which becomes (see Appendix D)

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle-\frac{2\pi
B_{0}^{2}r_{0}^{6}}{r^{4}}\delta(\omega-2\omega_{0})M_{lm}\delta_{m2}$ (50)
$\displaystyle+$
$\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{r^{4}}\\{\frac{1}{3}\sqrt{\frac{\pi}{5}}\delta(\omega-2\omega_{0})\delta_{l2}\delta_{m2}$
$\displaystyle+$
$\displaystyle[\frac{4\sqrt{\pi}}{3}\delta_{l0}\delta_{m0}-\frac{1}{\sqrt{5}}\delta_{l2}\delta_{m0}]\delta(\omega)\\}$

The solution of the inhomogeneous equation (49) now reads

$\displaystyle h_{\theta\theta}(t,r,\theta,\phi;r_{0},\omega_{0})$
$\displaystyle=$ $\displaystyle i\frac{B_{0}^{2}r_{0}^{6}}{6\pi
r}\sum_{lm}\frac{2^{l}\omega_{0}^{l}}{(2l+1)!!}\frac{R_{\rm
lc}^{l-3}-r_{0}^{l-3}}{l-3}Y_{lm}(\theta,\phi)$ $\displaystyle\cdot\
e^{i[2\omega_{0}t-2\omega_{0}r+\frac{l\pi}{2}]}[6\pi\delta_{m2}M_{lm}+\sqrt{5\pi}\delta_{l2}\delta_{m2}]$

where $M_{lm}$ is given in Appendix D. Because of the form of eq. (3.2), the
first non-zero term of the double sum is obtained for $l=2$ and $m=2$. In this
case $M_{22}=4$.

Inserting all values in eq. (3.2) we find

$\displaystyle h_{\theta\theta}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}\omega_{0}^{2}\sin^{2}{\theta}}{15r}\left(\frac{1}{R_{\rm
lc}}-\frac{1}{r_{0}}\right)\sqrt{\frac{5}{\pi}}(24+\frac{5}{\sqrt{\pi}})$ (52)
$\displaystyle\cdot\sin(2[\omega_{0}(t-r)-\phi])$

and also

$\displaystyle h_{\phi\phi}$ $\displaystyle=$ $\displaystyle-\sin^{2}\theta\
h_{\theta\theta}$ (53)

Figure 1: Thick red lines: approximate Advanced LIGO and Einstein Telescope
sensitivity limits. Dashed line: amplitude of gravitational wave radiation on
earth due to a $10^{16}$ G orthogonal magnetic rotator with initial rotation
period of 1 ms (at neutron star formation) at a distance of 1 kpc. Dotted
line: same for a $10^{15}$ G stellar magnetic field. The arrows show the
direction of evolution of the signal as the pulsar spins down. The magnetic
field contribution to the gravitational wave radiation of newly formed
strongly magnetized pulsars is borderline observable with advanced
gravitational wave detectors.

##  4 Discussion and Conclusions

We found that, in the TT-gauge, the $h_{\theta\phi}$ and $h_{\theta\theta}$
components of the gravitational wave due to a rotating orthogonal dipole
magnetic field behave as

$h_{\theta\theta}\sim
h_{\phi\phi}\sim(r_{0}\omega_{0})^{2}~{}~{}\mbox{~{}and~{}}~{}~{}h_{\theta\phi}\sim(r_{0}\omega_{0})^{3}$
(54)

Inserting real units in eqs. (52) and (53) we obtain

$\displaystyle h_{\theta\theta}$ $\displaystyle=$ $\displaystyle
2.26\times\frac{GB_{0}^{2}r_{0}^{3}}{c^{4}}\frac{1}{r}\frac{r_{0}^{2}\omega_{0}^{2}\sin^{2}{\theta}}{c^{2}}\sin(2[\omega_{0}(t-r)-\phi])$
(55) $\displaystyle=$ $\displaystyle 7\times
10^{-26}\left(\frac{B_{0}}{10^{15}\ {\rm G}}\right)^{2}\left(\frac{P}{10{\rm
ms}}\right)^{-2}\left(\frac{r}{{\rm kpc}}\right)^{-1}$ $\displaystyle\ \ \ \ \
\ \ \ \ \ \cdot\sin^{2}{\theta}\cdot\sin(2[\omega_{0}(t-r)-\phi])$

and

$\displaystyle h_{\phi\phi}$ $\displaystyle=$ $\displaystyle-\sin^{2}\theta\
h_{\theta\theta}$ (56) $\displaystyle=$ $\displaystyle-7\times 10^{-26}\
\left(\frac{B_{0}}{10^{15}\ {\rm G}}\right)^{2}\left(\frac{P}{10{\rm
ms}}\right)^{-2}\left(\frac{r}{{\rm kpc}}\right)^{-1}$ $\displaystyle\ \ \ \ \
\ \ \ \ \ \cdot\sin^{4}\theta\cdot\sin(2[\omega_{0}(t-r)-\phi])$

For galactic pulsars that are born with periods $P\sim 1\ {\rm ms}$, radius
$r_{0}=12$ km and extra strong magnetic fields $B_{0}>10^{15}\ {\rm G}$ within
a distance of about 1 kpc, this gravitational signal is borderline observable
by the new generation of gravitational wave detectors (Advanced LIGO and
Einstein Telescope; see figure 1) around the time of the neutron star
formation when their frequency of rotation is higher than about 100 Hz. If
observed, the shape of the gravitational wave signal would yield details about
the distribution of the magnetic field in the inner neutron star
magnetosphere. We have assumed here an orthogonal dipole field magnetosphere,
but our results may be generalized for higher order multipole magnetic field
configurations (quadrupole, etc.). It may be difficult to differentiate
between gravitational wave signals coming from the magnetosphere and from a
potential non-sphericity of the neutron star matter distribution. Notice that
a result similar to eqs. (55) and (56) was obtained in [Contopoulos et al.
2022] with an order of magnitude calculation of the magnetic field asymmetry
in the gravitational collapse of supermassive star cores. In that case, we
predicted a particular waveform for the resulting gravitation wave (figure 5
in that paper). In the present case of a young strongly magnetized neutron
star, the waveform will possess consistent frequency and amplitude, akin to a
singer holding a sustained note.

As we acknowledged in the Introduction, Nazari & Roshan (2020) solved
essentially the same problem, namely the direct contribution of a rotating
dipolar magnetic field to the generation of gravitational waves. Their results
are more general than ours since they consider all inclinations of the dipole
magnetic field with respect to the axis of rotation, but their method is
different. They implement a PN approximation to the energy–momentum tensor of
the magnetized fluid and to the gravitational potential in the wave zone. We
instead work directly with the perturbation of the stress-energy tensor
$T^{(1)}_{\mu\nu}$ (eq. 2), and study its effect on the background metric
$h_{\mu\nu}$ (eq. 1) by solving the wave equation that connects the two (eq.
14) in the TT-gauge. In principle, the two methods are expected to yield the
same result in the far zone. However, our final estimate of the strain of the
gravitational wave observed on earth (eqs. 55 & 56) is two orders of magnitude
higher than their corresponding numerical estimate (their eq. 48). We conclude
that the gravitational wave due to the rotating electromagnetic field is
borderline observable with the new generation of gravitational wave detectors,
while they conclude that it is unobservable and in order to obtain an
observable signal they resort to signal-to-noise (SN) amplification with
integration times of several days (Moore, Cole & Berry 2014).

Their main result in their eq. (47) looks very similar to ours in eqs. (55)
and (56), only the orders of magnitude differ. After checking their
calculations more carefully we realized that their magnetic field value is the
polar value of the stellar dipole magnetic field which is double our $B_{0}$
value. When we recalculated ourselves their result we obtained a different
numerical factor, namely

$\displaystyle h_{0F}$ $\displaystyle=$
$\displaystyle\frac{32}{105}\frac{\pi^{3}G}{\mu_{0}c^{6}}\frac{(2B_{0})^{2}r_{0}^{5}}{P^{2}r}$
(57) $\displaystyle=$ $\displaystyle 3\times 10^{-26}\
\left(\frac{B_{0}}{10^{15}\ {\rm G}}\right)^{2}\left(\frac{P}{10{\rm
ms}}\right)^{-2}\left(\frac{r}{{\rm kpc}}\right)^{-1}$

which is now of the same order of magnitude as our result (notice that the
magnetic permeability of vacuum $\mu_{0}=1$ in the CGS units that we are both
using). These extra two orders of magnitude that we obtained in this
recalculation allow us to claim that our results are in agreement with the
corrected estimates of Nazari & Roshan (2020).

We decided to further check the order of magnitude of our result through an
analogy with the calculation of the strain of the gravitational wave $h_{M}$
generated by a rotating deformed neutron star. In that case,

$\displaystyle h_{M}$ $\displaystyle=$
$\displaystyle\frac{4GI_{zzM}\Omega^{2}}{c^{4}r}\epsilon_{M}=\frac{8GM_{*}r_{0}^{2}\Omega^{2}}{5c^{4}r}\epsilon_{M}$
(58) $\displaystyle=$ $\displaystyle 7\times
10^{-20}\epsilon_{M}\left(\frac{P}{10{\rm ms}}\right)^{-2}\left(\frac{r}{{\rm
kpc}}\right)^{-1}$

(eq. 31a from Nazari & Roshan 2020), where,

$\displaystyle I_{zzM}$ $\displaystyle=$
$\displaystyle\int_{r=0}^{r_{0}}\int_{\theta=0}^{\pi}\int_{\phi=0}^{2\pi}r^{2}\sin\theta\
dr\ d\theta\ d\phi\ \rho_{M}(x^{2}+y^{2})\approx\frac{2}{5}M_{*}r_{0}^{2}$
(59) $\displaystyle I_{xxM}$ $\displaystyle=$
$\displaystyle\int_{r=0}^{r_{0}}\int_{\theta=0}^{\pi}\int_{\phi=0}^{2\pi}r^{2}\sin\theta\
dr\ d\theta\ d\phi\ \rho_{M}(y^{2}+z^{2})$ (60) $\displaystyle\epsilon_{M}$
$\displaystyle=$ $\displaystyle 1-\frac{I_{xxM}}{I_{zzM}}$ (61)

Here, $I_{zzM}$ is the stellar matter moment of inertia around the $z$ axis of
rotation, $\epsilon_{M}$ is the ellipticity of the stellar matter distribution
(eq. 4 of [Mastrano et al. 2015]), and $M_{*}\approx 1.4M_{\odot}$ is the mass
of the neutron star. The main difference between our calculation and eqs.
(58-61) is that, in our case, the source of the gravitational wave is the
energy density distribution outside the surface of the star, therefore, we
must calculate the above volume integrals for $r>r_{0}$. If we replace
$\rho_{M}$ in eqs. (60-61) by $B^{2}/(8\pi
c^{2})=[(B^{r})^{2}+(B^{\theta})^{2}+(B^{\phi})^{2}]/(8\pi c^{2})$ where the
magnetic field components are given by eqs. (17-19) for $t=0$, and we also
replace the radial limits of integration from $(0,r_{0})$ to $(r_{0},\infty)$
we obtain

$\displaystyle I_{zzF}$ $\displaystyle=$
$\displaystyle\int_{r=r_{0}}^{\infty}\int_{\theta=0}^{\pi}\int_{\phi=0}^{2\pi}r^{2}\sin\theta\
dr\ d\theta\ d\phi\ B^{2}(x^{2}+y^{2})/(8\pi c^{2})$ (62) $\displaystyle=$
$\displaystyle\frac{11}{15}\frac{B_{0}^{2}r_{0}^{5}}{c^{2}}$ $\displaystyle
I_{xxF}$ $\displaystyle=$
$\displaystyle\int_{r=r_{0}}^{\infty}\int_{\theta=0}^{\pi}\int_{\phi=0}^{2\pi}r^{2}\sin\theta\
dr\ d\theta\ d\phi\ B^{2}(y^{2}+z^{2})/(8\pi c^{2})$ (63) $\displaystyle=$
$\displaystyle\frac{8}{15}\frac{B_{0}^{2}r_{0}^{5}}{c^{2}}$
$\displaystyle\epsilon_{F}$ $\displaystyle=$ $\displaystyle
1-\frac{I_{xxF}}{I_{zzF}}=\frac{3}{11}=0.27\ .$ (64)

Using the above expressions, we can obtain an order of magnitude estimate of
the electromagnetic contribution to the gravitational wave signal $h_{F}$
through an analogy with the calculation that led to eq. (58), namely

$\displaystyle h_{F}$ $\displaystyle\sim$ $\displaystyle
h_{M}\frac{I_{zzF}}{I_{zzM}}\frac{\epsilon_{F}}{\epsilon_{M}}=\frac{1}{2}\left(\frac{B_{0}^{2}r_{0}^{3}}{M_{*}c^{2}}\right)\frac{8GM_{*}r_{0}^{2}\Omega^{2}}{5c^{4}r}$
(65) $\displaystyle=$ $\displaystyle 2.4\times 10^{-26}\
\left(\frac{B_{0}}{10^{15}\ {\rm G}}\right)^{2}\left(\frac{P}{10{\rm
ms}}\right)^{-2}\left(\frac{r}{{\rm kpc}}\right)^{-1}\ .$

This estimate is of the same order of magnitude as our results in eqs. (55),
(56) and (57).

In the present work we considered the perturbation of the electromagnetic
field energy on the Minkowski (empty) space time. In a forthcoming publication
we will investigate the generation of gravitational waves from rotating
neutron star magnetospheres in the Schwarzschild and the Kerr space time.

## Data availability statement

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * [Arfken et al. 2001] Arfken, G. B., Weber, H. J. 2001, Mathematical Methods for Physicists, Academic Press, San Diego, CA (USA)
  * [Bambi et al. 2022] Bambi,C., Katsanevas, S., Kokkotas, K. D. 2022, Handbook of Gravitational Wave Astronomy, Springer
  * [Bogovalov 1999] Bogovalov, S. 1999, Astron. Astrophys., 349, 1017
  * [Contopoulos et al. 2022] Contopoulos, I., Strantzalis, A., Papadopoulos, D. B., Kazanas, D. 2022, Mon. Not. R. Astr. Soc., 509, 174
  * [Ivanov 2008] Ivanov, B. V. 2008, Phys.Rev.D., 77, 044007
  * [Jackson 1999] Jackson, J. D. 1999, Classical Electrodynamics, John Wiley and Sons, Inc., third Ed.
  * [1] Korn, G. A., Korn, T. M. 1968, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems and Formulas for Reference and Review (MacGraw-Hill Book Co, (NY)
  * [Maggiore 2008] Maggiore, M. 2008, Gravitational Waves, Oxford University Press
  * [Mastrano et al. 2015] Mastrano, A., Suvorov, A. G., Melatos, A. 2015, Mon. Not. Roy. Astron. Soc. , 447, 3475
  * [Misner et al. 1973] Misner, C. W., Thorne, K., Wheller, J. 1973, Gravitation, W.H.Freeman, San Francisco
  * [Mpore et al. 2014] Moore C. J., Cole R. H., Berry C. P. 2014, CQG, 32, 015014
  * [Muslimov & Harding 2009] Muslimov, A. G., Harding, A. K. 2009, Astrophys. J., 692, 140
  * [Nazari & ROshan 2020] Nazari, E., Roshan, M. 2020, Mon. Not. Roy. Astron. Soc. , 498, 110
  * [Persides 1974] Persides, S. 1974, J. Math. Phys., 15, 885
  * [Rainich 1925] Rainich, G. Y. 1925, Trans. Amer. Math. Soc., 27, 106
  * [Thorne et al. 1982] Thorne, K., Macdonald, D. 1982, Mon. Not. Roy. Astron. Soc. , 198, 339
  * [Turyshev et al. 2018] Turyshev, S. G., Toth, V. T. 2018, arXiv:1801.06253 [physics.optics]
  * [Wald 1984] Wald, R. M. 1984, General relativity, The University of Chicago, Ltd., London

## Appendix A

We will be searching for the steady-state solution to eqs. (24) in spherical
coordinates $(r,\theta,\phi)$ centered onto the central star. We will follow
the approach of [Muslimov & Harding 2009], and we will define the
transformation of partial time derivatives between our lab frame, and a
reference frame that rotates with the pulsar angular velocity $\omega_{0}$
around the axis $\theta=0$, namely

$\displaystyle\frac{\partial{\bf B}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial{\bf B}}{\partial t}\right|_{\rm
corot}-\nabla\times(\omega_{0}r\sin\theta\ \hat{\phi}\times{\bf B})\ ,$ (66)
$\displaystyle\frac{\partial{\bf E}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial{\bf E}}{\partial t}\right|_{\rm
corot}-\nabla\times(\omega_{0}r\sin\theta\ \hat{\phi}\times{\bf
E})+\omega_{0}r\sin\theta\ \hat{\phi}\ \nabla\cdot{\bf E}\ ,$ (67)
$\displaystyle\frac{\partial\rho_{\rm e}}{\partial t}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial\rho_{\rm e}}{\partial t}\right|_{\rm
corot}+\omega_{0}r\sin\theta\ \hat{\phi}\cdot\nabla\rho_{\rm e}\ .$ (68)

Time derivatives in that frame (denoted by the subscript “corot”) vanish in
steady state, and therefore, Maxwell’s equations that describe the steady
state of the force-free pulsar magnetosphere become

$\displaystyle\nabla\times\left({\bf E}+\frac{\omega_{0}r\sin\theta\
\hat{\phi}}{c}\times{\bf B}\right)$ $\displaystyle=$ $\displaystyle 0\ ,$ (69)
$\displaystyle\nabla\times\left({\bf B}-\frac{\omega_{0}r\sin\theta\
\hat{\phi}}{c}\times{\bf E}\right)$ $\displaystyle=$
$\displaystyle\frac{4\pi}{c}{\bf J}-\frac{\omega_{0}r\sin\theta\
\hat{\phi}}{c}\nabla\cdot{\bf E}\ .$ (70)

From eq. (69) above we obtain that

$\displaystyle{\bf E}=-\frac{\omega_{0}r\sin\theta}{c}\hat{\phi}\times{\bf B}\
.$ (71)

## Appendix B: Usefull relations

We use the following relations from [Arfken et al. 2001] for Legendre
polynomials $P_{n}(x)$ and spherical harmonics $Y_{lm}(\theta,\phi)$:

$P_{n}(x)=\frac{1}{2^{n}n!}(\frac{d}{dx})^{n}(x^{2}-1)^{n}$ (72)
$P_{n}^{m}(x)=(1-x^{2})^{m/2}(\frac{d}{dx})^{m}P_{n}(x)$ (73)
$Y_{lm}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos{\theta})e^{im\phi}$
(74)

and $Y_{l,-m}=(-1)^{m}Y^{*}_{lm}$ where $Y^{*}_{lm}$ is the complex conjugate
of $Y^{*}_{lm}$. From eq. (74) we compute the

$\displaystyle Y_{00}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{4\pi}}$ $\displaystyle Y_{20}(\theta,\phi)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\frac{5}{\pi}}P_{2}(\cos{\theta})$
$\displaystyle Y_{22}(\theta,\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\frac{5}{\pi}}P_{2}^{2}(\cos{\theta})e^{2i\phi}=\frac{3}{2}\sqrt{\frac{5}{\pi}}\sin^{2}{\theta}e^{2i\phi}$
$\displaystyle\sin^{2}{\theta}$ $\displaystyle=$
$\displaystyle\frac{4\sqrt{\pi}}{3}[Y_{00}(\theta,\phi)-\frac{Y_{20}(\theta,\phi)}{\sqrt{5}}]$
$\displaystyle\sin^{2}{\theta}e^{2i\phi}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\sqrt{\frac{\pi}{5}}Y_{22}(\theta,\phi)$

In our calculations we also use the spherical Bessel functions
$j_{l}(x),n_{l}(x)$ and the Hankel functions $H_{l}^{(1)},H_{l}^{(2)}$ defined
as

$H_{l}^{(1)}(x)\equiv j_{l}(x)+in_{l}(x),~{}~{}H_{l}^{(2)}(x)\equiv
j_{l}(x)-in_{l}(x)$ (77)

For approximate calculations, we may use limiting values of the above
fucntions to calculate the integrals in eq. (43). For $x\rightarrow\infty$, we
have

$\displaystyle j_{l}(x)=\frac{1}{x}\sin\left(x-\frac{l\pi}{2}\right)$
$\displaystyle
H_{l}^{(1)}(x)=-\frac{i}{x}\exp\left(-i(x-\frac{l\pi}{2})\right)$

Also, for $x<<1$, we have

$\displaystyle j_{l}(x)\approx\frac{x^{l}}{(2l+1)!!}$
$\displaystyle\eta_{l}(x)\approx-\frac{((2l-1)!!}{x^{l+1}}$ $\displaystyle
H_{l}^{(1),(2)}\approx\pm i\eta_{l}(x)$ (79)

[Arfken et al. 2001]. In our further calculations we use eqs. (Appendix B:
Usefull relations) in the limit $x=\omega_{0}r\ll 1$.

## Appendix C

eqs. (41) and (45) yield

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle-2\int_{-\infty}^{\infty}dte^{i\omega
t}\int(-\frac{B_{0}^{2}r_{0}^{6}}{4\pi
r^{4}}\sin{(2\theta)}\cos{(2w)})Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle=$ $\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{2\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int\sin{(2\theta)}[\frac{e^{2iw}+e^{-2iw}}{2}]Y_{lm}^{*}(\theta,\phi)d\Omega$

We define

$\displaystyle W_{lm}^{(1)}$ $\displaystyle\equiv$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{4\pi
r^{4}}\int_{-\infty}^{\infty}dte^{it(\omega-2\omega_{0})}\int\sin{(2\theta)}e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
(81) $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{2r^{4}})\delta(\omega-2\omega_{0})\int\sin{(2\theta)}e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$

and

$\displaystyle W_{lm}^{(2)}$ $\displaystyle\equiv$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{4\pi
r^{4}}\int_{-\infty}^{\infty}dte^{it(\omega+2\omega_{0})}\int\sin{(2\theta)}e^{-2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
(82) $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{2r^{4}}\delta(\omega+2\omega_{0})\int\sin{(2\theta)}e^{-2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$

Using the relations

$\displaystyle
2\sin{\theta}cos{\theta}e^{2i\phi}=\frac{2}{3}\sqrt{\frac{\pi}{5}}\frac{d}{d\theta}Y_{22}(\theta,\phi)\
,$ $\displaystyle
2\sin{\theta}cos{\theta}e^{-2i\phi}=\frac{2}{3}\sqrt{\frac{\pi}{5}}\frac{d}{d\theta}Y_{22}^{*}(\theta,\phi)\
.$

eq. (Appendix C) becomes

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle W_{lm}^{(1)}+W_{lm}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{2r^{4}}\\{\delta(\omega-2\omega_{0})\int\sin{(2\theta)}e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle+\delta(\omega+2\omega_{0})\int\sin{(2\theta)}e^{-2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega\\}$
$\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{3r^{4}}\sqrt{\frac{\pi}{5}}$
$\displaystyle\cdot\ \\{\delta(\omega-2\omega_{0})\int d\Omega
Y_{lm}^{*}(\theta,\phi)\frac{d}{d\theta}Y_{22}(\theta,\phi)$
$\displaystyle+\delta(\omega+2\omega_{0})\int d\Omega
Y_{lm}^{*}(\theta,\phi)\frac{d}{d\theta}Y_{22}^{*}(\theta,\phi)\\}$

We define

$A_{lm}\equiv\int d\Omega
Y_{lm}^{*}(\theta,\phi)\frac{d}{d\theta}Y_{22}(\theta,\phi)$ (85)

and

$\tilde{A}_{lm}\equiv\int d\Omega
Y_{lm}^{*}(\theta,\phi)\frac{d}{d\theta}Y_{22}^{*}(\theta,\phi)$ (86)

To calculate the $A_{lm}$ and $\tilde{A}_{lm}$, we use the orthogonality
condition

$\int_{0}^{2\pi}e^{-im_{1}\phi}e^{im_{2}\phi}d\phi=2\pi\delta_{m_{1}m_{2}}$
(87)

[Arfken et al. 2001] and we find

$\displaystyle A_{lm}$ $\displaystyle=$ $\displaystyle
2\int\sin{\theta}\cos{\theta}e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$ (88)
$\displaystyle=$ $\displaystyle
2\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}$ $\displaystyle\cdot\
\int_{0}^{2\pi}e^{2i\phi}e^{-im\phi}d\phi\int_{0}^{\pi}\sin{\theta}\cos{\theta}P_{l}^{m}(\cos{\theta})\sin{\theta}d\theta$
$\displaystyle=$ $\displaystyle
4\pi\delta_{m2}\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}$
$\displaystyle\cdot\
\int_{-1}^{1}\frac{1}{3}P_{2}^{1}(\cos{\theta})P_{l}^{m}(\cos{\theta})d\cos{\theta}$

For simplicity we set

$\displaystyle S_{lm}$ $\displaystyle\equiv$
$\displaystyle\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\int_{-1}^{1}\frac{1}{3}P_{2}^{1}(\cos{\theta})P_{l}^{m}(\cos{\theta})d\cos{\theta}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\int_{-1}^{1}\cos{\theta}\sqrt{1-\cos^{2}{\theta}}P_{l}^{m}(\cos{\theta})d\cos{\theta}$

and eq. (88) reads

$A_{lm}=4\pi\delta_{m2}S_{lm}$ (90)

Because of the relation

$\int_{0}^{2\pi}e^{-im_{1}\phi}e^{-im_{2}\phi}d\phi=0$ (91)

$\tilde{A}_{lm}=0$. Eventually, eq. (Appendix C) becomes

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{2r^{4}}\delta(\omega-2\omega_{0})\int\sin{(2\theta)}e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
(92) $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{3r^{4}}\sqrt{\frac{\pi}{5}}\delta(\omega-2\omega_{0})4\pi\delta_{m2}S_{lm}$

As we mention in Sec.3, the solution of the inhomogeneous eq. (34) will be of
the form of eq. (43). Therfore, we need to compute the following integrals:

$\displaystyle J_{lm}(r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle\int_{r_{0}}^{R_{\rm lc}}j_{l}(\omega
r)f_{lm}(r;\omega;\omega_{0})dr$ $\displaystyle=$
$\displaystyle\frac{B_{0}^{2}r_{0}^{6}}{3}\sqrt{\frac{\pi}{5}}\int_{r_{0}}^{R_{\rm
lc}}\frac{\omega^{l}r^{l}}{r^{4}(2l+1)!!}[A_{lm}\delta(\omega-2\omega_{0})]dr$
$\displaystyle=$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{3}\sqrt{\frac{\pi}{5}}\frac{\omega^{l}}{(2l+1)!!}\frac{r_{0}^{(l-3)}-R_{\rm
lc}^{(l-3)}}{l-3}4\pi\delta_{lm}S_{lm}\delta(\omega-2\omega_{0})$

and

$\displaystyle I_{lm}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\omega d\omega e^{-i\omega
t}H_{l}^{(1)}(\omega r)J_{lm}(r_{0};\omega,\omega_{0})$ (94) $\displaystyle=$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{3r}\sqrt{\frac{5}{\pi}}\frac{2^{l}\omega_{0}^{l}}{(2l+1)!!}\frac{r_{0}^{(l-3)}-R_{\rm
lc}^{(l-3)}}{l-3}$ $\displaystyle\cdot\
4\pi\delta_{lm}S_{lm}e^{\frac{i}{2}[4t\omega_{0}-4\omega_{0}r+l\pi]}$

Now, as it is apparent, with the aid of eqs. (Appendix C) and (94) we can find
the solution of the inhomogeneous eq. (44).

## Appendix D

We re-write

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle-2\int_{-\infty}^{\infty}dte^{i\omega t}$ $\displaystyle\cdot\
\int\frac{B_{0}^{2}r_{0}^{6}}{2\pi
r^{4}}[\cos{(2w)}-5\sin^{2}{\theta}\cos^{2}{w}]Y_{lm}^{*}(\theta,\phi)d\Omega$

We define

$\displaystyle Q_{lm}^{(1)}$ $\displaystyle\equiv$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int[\frac{e^{2iw}+e^{-2iw}}{2}]Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle Q_{lm}^{(2)}$ $\displaystyle\equiv$
$\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int\sin^{2}{\theta}\cos^{2}{w}Y_{lm}^{*}(\theta,\phi)d\Omega$ (96)

Below, we calculate $Q_{lm}^{(1)}$ and $Q_{lm}^{(2)}$ and upon the
consideration of eqs. (74), (87) and (91), we find

$\displaystyle Q_{lm}^{(1)}$ $\displaystyle\equiv$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int[\frac{e^{2iw}+e^{-2iw}}{2}]Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle=$ $\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{2\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int[e^{-2i\omega_{0}t}e^{2i\phi}+e^{2i\omega_{0}t}e^{-2i\phi}]Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle=$ $\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{2\pi
r^{4}}\\{2\pi\delta(\omega-2\omega_{0})\int
e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle+2\pi\delta(\omega+2\omega_{0})\int
e^{-2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega\\}$ $\displaystyle=$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{r^{4}}\delta(\omega-2\omega_{0})\int
e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$

Working as previously and taking into account eqs. (74), (87) and (91) we find

$\displaystyle Q_{lm}^{(2)}$ $\displaystyle\equiv$
$\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int[\sin^{2}{\theta}\cos^{2}{w}]Y_{lm}^{*}(\theta,\phi)d\Omega$ (98)
$\displaystyle=$ $\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int\sin^{2}{\theta}Y_{lm}^{*}(\theta,\phi)\frac{e^{2iw}+e^{-2iw}+2}{4}$
$\displaystyle=$ $\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{4\pi
r^{4}}\int_{-\infty}^{\infty}dte^{i\omega
t}\int\sin^{2}{\theta}Y_{lm}^{*}(\theta,\phi)[e^{-2i\omega_{0}t}e^{2i\phi}+e^{2i\omega_{0}t}e^{-2i\phi}+2]$
$\displaystyle=$ $\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{4\pi
r^{4}}\\{2\pi\delta(\omega-2\omega_{0})\frac{2}{3}\sqrt{\frac{\pi}{5}}\int
d\Omega Y_{lm}^{*}(\theta,\phi)Y_{22}(\theta,\phi)$
$\displaystyle+2\pi\delta(\omega+2\omega_{0})\frac{2}{3}\sqrt{\frac{\pi}{5}}\int
d\Omega Y_{lm}^{*}(\theta,\phi)Y_{22}^{*}(\theta,\phi)$
$\displaystyle+4\pi\delta(\omega)\int d\Omega
Y_{lm}^{*}(\theta,\phi)\sin^{2}{\theta}\\}$ $\displaystyle=$
$\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{2r^{4}}\\{\frac{2}{3}\sqrt{\frac{\pi}{5}}\delta(\omega-2\omega_{0})\delta_{l2}\delta_{m2}$
$\displaystyle+2[\frac{4\sqrt{\pi}}{3}\delta_{l0}\delta_{m0}-\frac{1}{\sqrt{5}}\delta_{l2}\delta_{m0}]\delta(\omega)\\}$

Thus, the source term (LABEL:d7xx) reads

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle Q_{lm}^{(1)}+Q_{lm}^{(2)}$ (99) $\displaystyle=$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}}{r^{4}}\delta(\omega-2\omega_{0})\int
e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$
$\displaystyle+\frac{5B_{0}^{2}r_{0}^{6}}{r^{4}}\\{\frac{1}{3}\sqrt{\frac{\pi}{5}}\delta(\omega-2\omega_{0})\delta_{l2}\delta_{m2}$
$\displaystyle+[\frac{4\sqrt{\pi}}{3}\delta_{l0}\delta_{m0}-\frac{1}{\sqrt{5}}\delta_{l2}\delta_{m0}]\delta(\omega)\\}$

Below, we calculate the expression

$\displaystyle B_{lm}$ $\displaystyle\equiv$ $\displaystyle\int
e^{2i\phi}Y_{lm}^{*}(\theta,\phi)d\Omega$ (100) $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}e^{-im\phi}e^{2i\phi}d\phi\int_{0}^{\pi}\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos{\theta})\sin{\theta}d\theta$
$\displaystyle=$ $\displaystyle
2\pi\delta_{m2}\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\int_{-1}^{1}P_{l}^{m}(\cos{\theta})d\cos{\theta}$

For simplicity we call

$M_{lm}\equiv\sqrt{\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\int_{-1}^{1}P_{l}^{m}(\cos{\theta})d\cos{\theta}$
(101)

and thus, eq. (100) reads

$B_{lm}=2\pi M_{lm}\delta_{m2}$ (102)

Eqs. (99) and (102) then yield

$\displaystyle f_{lm}(r,r_{0};\omega;\omega_{0})$ $\displaystyle=$
$\displaystyle-\frac{2\pi
B_{0}^{2}r_{0}^{6}}{r^{4}}\delta(\omega-2\omega_{0})M_{lm}\delta_{m2}$ (103)
$\displaystyle+$
$\displaystyle\frac{5B_{0}^{2}r_{0}^{6}}{r^{4}}\\{\frac{1}{3}\sqrt{\frac{\pi}{5}}\delta(\omega-2\omega_{0})\delta_{l2}\delta_{m2}$
$\displaystyle+$
$\displaystyle[\frac{4\sqrt{\pi}}{3}\delta_{l0}\delta_{m0}-\frac{1}{\sqrt{5}}\delta_{l2}\delta_{m0}]\delta(\omega)\\}$

In Sec.3, we point out that eq. (49) is of the form of eq. (43) and thus, we
need to compute the following integrals: First, we compute

$J_{lm}=\int_{r_{0}}^{R_{\rm lc}}j_{l}(\omega\rho)f_{lm}(\rho;\omega)d\rho$
(104)

Using the approximate eqs. (Appendix B: Usefull relations), we find

$\displaystyle J_{lm}$ $\displaystyle=$ $\displaystyle\int_{r_{0}}^{R_{\rm
lc}}j_{l}(\omega r)f_{lm}(r,r_{0};\omega,\omega_{0})dr$ $\displaystyle=$
$\displaystyle\int_{r_{0}}^{R_{\rm
lc}}\frac{\omega^{l}r^{l}}{(2l+1)!!}f_{lm}(r,r_{0};\omega,\omega_{0})dr$
$\displaystyle=$
$\displaystyle-\frac{B_{0}^{2}r_{0}^{6}\omega^{l}}{(2l+1)!!}\frac{R_{\rm
lc}^{l-3}-r_{0}^{l-3}}{l-3}\\{2\pi\delta_{m2}M_{lm}\delta(\omega-2\omega_{0})$
$\displaystyle+\frac{\sqrt{5\pi}}{3}\delta(\omega-2\omega_{0})\delta_{l2}\delta_{m2}+\frac{20\sqrt{\pi}}{3}\delta(\omega)\delta_{l0}\delta_{m0}-\sqrt{5}\delta(\omega)\delta_{l2}\delta_{m0}\\}$

Second, using the approximate expressions (Appendix B: Usefull relations), we
calculate

$\displaystyle I_{lm}(r,r_{0};\omega,\omega_{0})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\omega d\omega e^{-i\omega
t}H_{l}^{(1)}(\omega r)J_{lm}(r,r_{0};\omega,\omega_{0})$ $\displaystyle=$
$\displaystyle\frac{iB_{0}^{2}r_{0}^{6}}{6\pi
r}\frac{2^{l}\omega_{0}^{l}}{(2l+1)!!}\frac{r_{0}^{l-3}-R_{\rm
lc}^{l-3}}{l-3}$ $\displaystyle\cdot\
e^{i[2\omega_{0}t-2\omega_{0}r+\frac{l\pi}{2}]}[6\pi\delta_{m2}M_{lm}+\sqrt{5\pi}\delta_{l2}\delta_{m2}]$

With the aid of eq. (43) and the eqs. (Appendix D) and (Appendix D) we find
the solution of the inhomogeneous eq. (49).