2012 0

530.145:535.12:681.7:53.082.5

# Measurement-device-independent continuous variable quantum key distribution
protocol operation in optical transport networks

1I. Vorontsova  1R. Goncharov  1,2S. Kynev  1,2F. Kiselev  1,2V. Egorov  1
ITMO University,
Kronverkskiy, 49, St. Petersburg, 197101, Russia  2 SMARTS-Quanttelecom LLC,
6th Vasilyevskogo Ostrova Line, 59, St. Petersburg, 199178, Russia
<EMAIL_ADDRESS>
###### Abstract

Numerically, a theoretical analysis of the noise impact caused by spontaneous
Raman scattering, four-wave mixing, and linear channel crosstalk on the
measurement-device-independent continuous variable quantum key distribution
systems is conducted. The analysis considers symmetry and asymmetry of system
paths, as well as possible channel allocation schemes, for a quantum channel
located in C- and O-bans. Mathematical models for MDI CV-QKD system and the
contributing noises’ description are provided. The secure key generation rate
is estimated to state features of protocol operation when integrated with
existing DWDM systems in the context of its implementation into
telecommunication networks.

###### pacs:

03.67.-a, 42.50.-p

###### keywords:

device-independence, quantum key distribution, continuous variables

††volume: 0

##  1 Introduction

Quantum key distribution (QKD) is one of the most rapidly advancing fields of
quantum technologies [1, 2]. Its main idea is an opportunity to transport a
cryptographically secure key between two and more authenticated users
connected to each other through quantum and information channels. Guarantied
by the principles of quantum mechanics [3], the security of QKD to attacks
from an eavesdropper ensures safety of the transmitted data from all kinds of
hacking and known attacks.

One option to classify QKD protocols is based on them being discrete-variable
(DV) or continuous-variable (CV) [4]. Additionally, among many other QKD
protocol classifications, there is the one distinguishing between protocols in
terms of their device-dependence or (semi-)device-independence [5]. The
intersection of these two criteria gives rise to a new class of protocols,
that is measurement-device-independent (MDI) CV-QKD protocol. Device-
independence features particular practical importance, for it eliminates many
side channel attacks, though implies accurate theoretical analysis.

Not only does this work discuss the latter, but also for the first time
combines this analysis with the task of simultaneous propagation of
information and quantum signals in a single optical fiber [6]. The effects
considered as channel loss sources are spontaneous Raman scattering, four-wave
mixing, and linear channel crosstalk. A possible realization scheme is
discussed, as well as the allocation of classical channels on the standard
DWDM grid. The security analysis is carried out numerically, employing the
known theoretical security bounds to estimate the performance of the addressed
QKD system. The results are quite important in practice to be considered when
integrating QKD with the existing telecommunication networks.

##  2 Measurement-device-independent QKD

Let us discuss main principles underlying MDI QKD protocol operation [5, 7].
The essence of the approach lies in the fact that no assumptions are made
about the detectors involved in the protocol, such that they can even be
controlled by an eavesdropper (Eve). In a typical single-photon MDI QKD
protocol, two legitimate users (Alice and Bob) send quantum signals to an
untrusted central relay, often addressed as Charlie. A Bell state measurement
is performed then; both signals interfere at a 50:50 beam splitter (BS). Next,
output signals go through a polarizing beam splitter (PBS) to be projected
into either horizontal (H) or vertical (V) polarization state. The measurement
is pronounced successful if two of the four involved detectors click.

###  2.1 Continuous-variable MDI QKD protocol

Similarly to the conventional MDI QKD, the continuous variable (CV) version of
the protocol [8, 9] again implies there are the two senders and an untrusted
relay performing the measurements to be used during legitimate users’ post-
processing to generate the secure key.

The two known approaches to a general protocol description, namely, “prepare-
and-measure” (PM) and “entanglement-based” (EB) scenarios, apply to the case
of MDI CV-QKD as well. For these scenarios are equivalent in terms of their
mathematical description and effectiveness, we will consider a more
practically convenient PM version of the protocol. Gaussian modulation [4, 10]
(GG02 protocol) is considered, so Alice and Bob are operating with coherent
states with a two-dimensional Gaussian distribution. They first generate
coherent states $|x_{\mathrm{A}}+ip_{\mathrm{A}}\rangle$ and
$|x_{\mathrm{B}}+ip_{\mathrm{B}}\rangle$ with the quadratures $x$ and $p$
featuring variance $V_{\mathrm{A}(\mathrm{B})}-1$ (in shot noise units (SNU))
and then send their states to Charlie via quantum channels. Next, Alice’s and
Bob’s modes interfere at the beam splitter, while Charlie measures the
$\mathrm{C}$ and $\mathrm{D}$ modes’ quadratures on homodyne detectors and
announces the resulting state $\\{X_{\mathrm{C}},\>P_{\mathrm{D}}\\}$
publicly. It is Bob only who changes his state according to
$X_{\mathrm{B}}=x_{\mathrm{B}}+kX_{C},P_{\mathrm{B}}=p_{\mathrm{B}}-kP_{\mathrm{D}}$
(with $k$ standing for the gain associated with channel losses), whereas
Alice’s state remains unchanged. Finally, standard procedures are utilized for
parameter estimation, information reconciliation, and privacy amplification.

Since there is equivalency between the CV-QKD EB and PM scenarios’ security
proof against collective attacks [11, 12], we shall now switch to the well-
known covariance matrix formalism.

Implying that Eve has access to the relays, quantum channels, and even Bob’s
state displacement operation in the EB scheme, further security analysis of
the MDI CV-QKD protocol can be seen as a special case of a typical one-way CV-
QKD protocol [4, 10].

Then the secure key fraction can be estimated in accordance with the Devetak-
Winter bound [13, 14]:

$\displaystyle r=\beta
I(X_{\mathrm{A}},\>P_{\mathrm{A}}:X_{\mathrm{B}},\>P_{\mathrm{B}})-\chi(X_{\mathrm{B}},\>P_{\mathrm{B}}:E),$
(1)

where $0\leqslant\beta\leqslant 1$ is the reconciliation efficiency (assumed
to be ideal in the further numerical simulations), $I$ is the mutual
information between Alice and Bob,
$\chi(X_{\mathrm{B}},P_{\mathrm{B}}:E)=S(\widehat{\rho}_{E})-S(\widehat{\rho}_{E}|X_{\mathrm{B}},P_{\mathrm{B}})$
is the Holevo bound, and $S(\widehat{\rho}_{E})$ denotes the von Neumann
entropy of quantum state $\widehat{\rho}_{E}$.

The upper bound
$\chi(X_{\mathrm{B}},\>P_{\mathrm{B}}:\mathrm{A}_{1},\mathrm{B}^{\prime}_{1})$
is determined only using the corresponding covariance matrix. Supposing the
system is under two independent entangling cloner attacks [4], the covariance
matrix takes the form of:

$\displaystyle\Xi$
$\displaystyle=\left(\begin{array}[]{cc}{V_{\mathrm{A}}I_{2}}&{\sqrt{(T(V_{\mathrm{A}}^{2}-1)\sigma_{z})}}\\\
{\sqrt{(T(V_{\mathrm{A}}^{2}-1)\sigma_{z})}}&{[(V_{\mathrm{A}}-1)+1+T\xi^{\prime}]I_{2}}\end{array}\right),$
(4) $\displaystyle T$ $\displaystyle=\frac{\eta_{\mathrm{A}}}{2}g^{2},$ (5)
$\displaystyle\xi^{\prime}$
$\displaystyle=1+\frac{1}{\eta_{\mathrm{A}}}[\eta_{\mathrm{B}}(\Xi_{\mathrm{B}}-1)+\eta_{\mathrm{A}}\Xi_{\mathrm{A}}]+\frac{1}{\eta_{\mathrm{A}}}(\frac{\sqrt{2}}{g}\sqrt{V_{\mathrm{B}}-1}-\sqrt{\eta_{\mathrm{B}}}\sqrt{V_{\mathrm{B}}+1})^{2},$
(6) $\displaystyle\Xi_{\mathrm{A}}$
$\displaystyle=\frac{1-\eta_{\mathrm{A}}}{\eta_{\mathrm{A}}}+\xi_{\mathrm{A}},\>\Xi_{\mathrm{B}}=\frac{1-\eta_{\mathrm{B}}}{\eta_{\mathrm{B}}}+\xi_{\mathrm{B}},$
(7) $\displaystyle\eta_{\mathrm{A}}$ $\displaystyle=10^{-\alpha
L_{\mathrm{AC}}/10},\>\eta_{\mathrm{B}}=10^{-\alpha L_{\mathrm{BC}}/10},$ (8)

where $\eta_{\mathrm{A}}\>(\eta_{\mathrm{B}})$ is a channel (Alice-Charlie or
Bob-Charlie) transmittance, $\xi_{\mathrm{A}}(\xi_{\mathrm{B}})$ is the excess
noise, $g$ is the offset factor, $I_{2}$ is the identity matrix, and $\sigma$
is the Pauli $z$-matrix.

To minimize the excess noise, the offset factor is set as

$\displaystyle
g=\sqrt{\frac{2}{\eta_{\mathrm{B}}}}\sqrt{\frac{V_{\mathrm{B}}-1}{V_{\mathrm{B}}+1}}.$
(9)

Then, the excess noise is expressed as:

$\displaystyle\xi^{\prime}=\xi_{\mathrm{A}}+\frac{1}{\eta_{\mathrm{A}}}[\eta_{\mathrm{B}}(\xi_{\mathrm{B}}-2)+2],$
(10)

where the excess noise on Alice’s side (Bob’s side)

$\xi_{\mathrm{A(B)}}$

contain the corresponding channel noise converted to SNU. Alice’s and Bob’s
variances are considered equal $V_{\mathrm{A}}=V_{\mathrm{B}}=40$ in the
simulations.

##  3 Channel Noise Sources and their Mathematical Description

Naturally, losses are inevitable when it comes to signal propagation in any
medium, be it fiber-optical communication lines or free space. Regarding QKD,
three effects are addressed in terms of noise primarily, which are the
spontaneous Raman scattering (SpRS), the four-wave mixing (FWM) nonlinearity,
and the linear channel crosstalk (LCXT). Now, we will briefly summarize their
physical nature and corresponding mathematical description to then estimate
the negative contribution they make to the performance of MDI CV-QKD system
under consideration.

###  3.1 Spontaneous Raman Scattering

Firstly, the main contributor to the overall channel losses in case QKD
integrated with DWDM systems is the SpRS noise [15, 16]. Its impact is minor
for classical networks, though the contribution becomes substantial for joint
QKD and DWDM systems [15, 16]. The origin of the SpRS is different for the
cases of co- and counter-propagation of signals. Thus, two sub-types are
distinguished usually, that are forward (for co-propagating signals) and
backward (for counter-propagating signals) SpRS noise. When speaking in the
context of simultaneous QKD session and information transmission within single
oprical fiber, the mathematical representation of their contribution is given
by [6, 17]:

$P_{\mathrm{ram,f}}=P_{\mathrm{out}}L\sum_{c=1}^{N_{\mathrm{ch}}}\rho(\lambda_{\mathrm{c}},\lambda_{\mathrm{q}})\Delta\lambda,$
(11)

and

$P_{ram,b}=P_{out}\frac{\sinh(\xi
L)}{\xi}\sum_{c=1}^{N_{ch}}\rho(\lambda_{c},\lambda_{q})\Delta\lambda,$ (12)

where $P_{\mathrm{out}}$ denotes the output power for a single channel, $L$ is
the length of the optical fiber, $N_{ch}$ is the number of classical channels
present in a DWDM system, $\rho(\lambda_{\mathrm{c}},\lambda_{\mathrm{q}})$
describes the normalized scattering cross-section for the wavelengths of
classical ($\lambda_{\mathrm{c}}$) and quantum ($\lambda_{\mathrm{q}}$)
channels, and $\Delta\lambda$ is the bandwidth of the quantum channel
filtering system.

For the MDI CV-QKD realization considered here, both forward and backward SpRS
occur for different system paths. A detailed description of the configuration
will be provided in the following section.

The output power values appear in formulas instead of the input ones, so that
to meet the the bit error rate (BER) requirements of a DWDM system directly.
Thus, the value of $P_{\mathrm{out}}$ can be obtained as follows:

$P_{\mathrm{out}}\text{ (dBm)}=R_{\mathrm{x}}\text{ (dBm)}+IL\text{(dBm)},$
(13)

where $R_{\mathrm{x}}$ is the sensitivity of a receiver and $IL$ denotes the
insertion losses of the system.

###  3.2 Four-wave Mixing

Next channel noise source to consider is FWM. It is a third-order nonlinear
process, which sequence is creation of photons at new frequencies as a result
of the interaction between the initial ones [18]. These new frequencies might
coincide with the one of the quantum channel [19], thus contributing to the
overall noise in the band of the quantum channel.

To come up with the mathematical model for the FWM noise contribution
consideration, let us consider three pump channels with the frequencies
$f_{i},f_{j},\text{ and }f_{k}$. Then, the value of the resulting FWM noise
peak power $P_{ijk}$ featuring the frequency $f_{i}+f_{j}-f_{k}$ can be
expressed as [6]:

$P_{ijk}=\eta\gamma^{2}D^{2}p^{2}e^{-\xi L}\frac{(1-e^{-\xi
L})^{2}}{9\xi^{2}}P_{s}P_{l}P_{h},$ (14)

where the phase-matching efficiency for the FWM $\eta$ and parameter
$\Delta\beta$ are defined as:

$\eta=\frac{\xi^{2}}{\xi^{2}+\Delta\beta^{2}}\left[1+\frac{4e^{-\xi
L}\sin^{2}{(\Delta\beta L/2)}}{(1-e^{-\xi L})^{2}}\right],$ (15)

and

$\Delta\beta=\frac{2\pi\lambda^{2}}{c}\lvert f_{i}-f_{k}\rvert\lvert
f_{j}-f_{k}\rvert\cdot\left[D_{c}+\frac{dD_{c}}{d\lambda}\left(\frac{\lambda^{2}}{c}\right)\left(\lvert
f_{i}-f_{k}\rvert+\lvert f_{j}-f_{k}\rvert\right)\right],$ (16)

correspondingly. In the above equations, $L$ is the transmission distance of
the interacting light fields in the optical fiber, $D$ denotes the degeneracy
factor ($D=6$, $D=3$), $P_{i(j,k)}$ and $f_{i(j,k)}$ are the input power and
optical frequency of the interacting fields correspondingly, $\gamma$ stands
for the third-order nonlinear coefficient, $\xi$ is the loss coefficient,
$D_{c}$ and $dD_{c}/{d\lambda}$ are the dispersion coefficient of an optical
fiber and its slope respectively with $\lambda$ standing for the wavelength of
the FWM radiation.

Finally, performing the summation of all the powers of the resulting FWM
signals with frequencies coinciding with the frequency of the quantum channel,
one obtains:

$P_{\mathrm{FWM}}=\sum P_{ijk},f_{i}+f_{j}-f_{k}=f_{q}.$ (17)

###  3.3 Linear Channel Crosstalk

It is due to the imperfections of the demultiplexers that any practically
implemented DWDM system suffers LCXT losses [20].

Since information signals are orders of magnitude more powerful than quantum
ones, the insufficient isolation might cause considerable LCXT noise, that can
be estimated in the following way:

$P_{\mathrm{LCXT}}=P_{\mathrm{out}}\text{ (dBm)}-\text{ISOL (dB)}.$ (18)

Once the power values of the contributing to the overall channel noise effects
are calculated, it is needed to recalculate them to a photon detection
probability. To do so, the formula can be used:

$p_{\mathrm{ram,f(b)/FWM/LCXT}}=\frac{P_{\mathrm{ram,f(b)/FWM/LCXT}}}{hc/\lambda_{q}}\Delta
t\eta_{D}\eta_{B},$ (19)

where $\eta_{D}$ denotes the detector efficiency, $\eta_{B}=10^{-0.1IL}$ is
the transmission coefficient associated with the insertion losses of the
detection system, $h$ is the Planck constant, and $c$ is the light speed.

##  4 MDI CV-QKD scheme and channel allocation

Here, a possible realization scheme of MDI CV-QKD protocol will be addressed
to analyze its potential for creating telecommunication optical transport
networks integrated with DWDM systems. The notion of maximal achievable
distances of MDI CV-QKD systems employed to characterize the latter denotes
fiber length corresponding to the case, when the secure key generation rate is
non-zero.

The realization of MDI CV-QKD addressed in the work is shown in Figure 1).

Figure 1: Schematic illustration for the MDI CV-QKD protocol realization with
DWDM: information is transmitted from Alice to Bob via DWDM-information
channels; quantum channel is co-propagating with information channels in
Alice’s path, whereas counter-propagating with them in Bob’s path and, thus,
which means forward SpRS noise is induced in the Alice-Charlie path and
backward SpRS - in the Bob-Charlie path

Here, quantum signals are sent to the untrusted central relay to be homodyned
there, whereas the information is transferred from Alice to Bob directly. It
means, quantum and information signals are unidirectional in Alice’s path
(i.e., there is forward SpRS noise in the path) and counter-propagating
through the Bob’s path, so the SpRS noise features the backward type.

The performance of the MDI CV-QKD realization was then numerically analyzed in
terms of the possible channel allocation schemes and asymmetry coefficients
between Alice’s $L_{a}$ and Bob’s $L_{b}$ paths
$R_{\mathrm{asym}}=L_{a}/L_{b}$.

Similarly to the work [21], four allocation schemes for the quantum channel
located in C-band and O-band of the communication window were considered. The
criterion and complete explanation for such a choice is provided in detail in
the works [22, 23]. The final choice of the configurations considered in the
further numerical simulations is given in Table 1.

Table 1: Description of the optimal configurations chosen for numerical simulations Configuration |  | Number of
---
channels
| Quantum channel
---
wavelength, nm
Configuration 1 | 10 | 1536.61
Configuration 2 | 10 | 1310
Configuration 3 | 40 | 1537.40
Configuration 4 | 40 | 1310

The parameters of the DWDM system are the following: $\xi=0.18$ dB/km,
$\Delta\lambda=15$ GHz, $N_{ch}=10\text{ or }40$, $R_{x}=-32$ dBm and IL = 8
dB.

As for the asymmetry coefficient $R_{\mathrm{asym}}$, three different
relations are address here: a symmetric (i.e., $L_{a}/L_{b}=1$) and two
asymmetric realizations ($L_{a}/L_{b}=3/2$ and $L_{a}/L_{b}=2/1$).

##  5 Results and discussion

Using the mathematical models for the MDI CV-QKD secure key generation rate,
SpRS, FWM, and LCXT noises, the realization of the system depicted in Figure 1
was numerically simulated. The results are presented in Figure 2.

Figure 2: The dependence of the secure key generation rate on the optical
fiber length for MDI CV-QKD corresponding to the scheme in Figure 1.

It is a known fact that for MDI CV-QKD systems the secure key generation rate
decreases dramatically as the system approaches its symmetry ($L_{a}=L_{b}$),
with the best result corresponding to the situation when one of the paths
equals zero [24].

The results obtained confirm the stated conclusion, as the maximal achievable
distance increases together with the system’s paths asymmetry.

Interestingly, for the case of MDI CV-QKD, configurations for which the
quantum channel is located in the C-band appeared to be more efficient in
terms of their maximal achievable distances. For most cases, configurations
with O-band-located quantum channel are more beneficial (MDI QKD as well, see
[22]), as FWM noise does not contribute to the overall losses. Though the
overall contribution of channel noises here is less for such cases too, fiber
attenuation for the quantum channel wavelength of 1310 nm surpasses them
significantly. It can be seen also that the maximum achievable distances do
not exceed 6 km, thus, the secure key distribution over long distances is not
possible here. Still, such realizations can be utilized for short-distance
communication.

Regarding the number of information channels (10 and 40 in the work), for a
larger number of information channels, a decrease in the maximal achievable
distances takes place. This is a natural observation, as the more information
channels are there in the system, the greater the value of the overall channel
losses. The decrease is substantial for C-band configurations, whereas is
quite small in case of O-band.

##  6 Conclusion

In this work, the MDI CV-QKD protocol was addressed. A theoretical research
and numerical simulation of the noise influence caused by SpRS, FWM, and LCXT
effects on the performance of the MDI CV-QKD system performance was analyzed
for a proposed practical realization scheme, in terms of channel allocation
and the paths’ asymmetry coefficient. Increasing the number of information
channels naturally leads to a decrease in the maximal achievable distances. In
addition, the allocation of a wavelength of 1310 nm for a quantum channel
results in the shortening of maximal distance values for MDI CV-QKD,
regardless of the fact the overall channel noise is less for such
configurations. The superior contribution comes from the fiber attenuation,
which is larger for the wavelength of 1310 nm. It was confirmed that the more
asymmetric the paths for the MDI CV-QKD scheme are, the more efficient the
systems is. Moreover, MDI CV-QKD realizations feature significantly shorter
maximal achievable distances, that do not exceed several kilometers and can be
utilized for short-distance information transmission only. The results
obtained can be used in terms of practical implementation of MDI CV-QKD
systems, so that to obtain optimal results.

## Acknowledgements

The work was done by Leading Research Center ”National Center for Quantum
Internet” of ITMO University by the order of JSCo Russian Railways.

## References

  * [1] Valerio Scarani, Helle Bechmann-Pasquinucci, Nicolas J Cerf, Miloslav Dušek, Norbert Lütkenhaus, and Momtchil Peev.  The security of practical quantum key distribution.  Reviews of Modern Physics, 81(3):1301–1350, 9 2009.
  * [2] Stefano Pirandola, U. L. Andersen, Leonardo Banchi, Mario Berta, Darius Bunandar, Roger Colbeck, Dirk Englund, Tobias Gehring, Cosmo Lupo, Carlo Ottaviani, J. L. Pereira, Mohsen Razavi, Jesni Shamsul Shaari, Marco Tomamichel, Vladyslav C. Usenko, Giuseppe Vallone, Paolo Villoresi, and Petros Wallden.  Advances in quantum cryptography.  Advances in Optics and Photonics, 12(4):1012, 12 2020.
  * [3] Nicolas Gisin, Gré Goire Ribordy, Wolfgang Tittel, and Hugo Zbinden.  Quantum cryptography.  Reviews of modern physics, 74(1):145, 2002.
  * [4] Frédéric Grosshans, Gilles Van Assche, Jérôme Wenger, Rosa Brouri, Nicolas J. Cerf, and Philippe Grangier.  Quantum key distribution using gaussian-modulated coherent states.  Nature, 421(6920):238–241, 1 2003.
  * [5] Hoi Kwong Lo, Marcos Curty, and Bing Qi.  Measurement-device-independent quantum key distribution.  Physical Review Letters, 108(13), 3 2012.
  * [6] M. Mlejnek, Nikolay Kaliteevskiy, and D.A. Nolan.  Reducing spontaneous Raman scattering noise in high quantum bit rate QKD systems over optical fiber.  arXiv preprint, arXiv:1712.05891, 2017.
  * [7] Xiongfeng Ma and Mohsen Razavi.  Alternative schemes for measurement-device-independent quantum key distribution.  Physical Review A - Atomic, Molecular, and Optical Physics, 86(6), 12 2012.
  * [8] Panagiotis Papanastasiou, Carlo Ottaviani, and Stefano Pirandola.  Finite-size analysis of measurement-device-independent quantum cryptography with continuous variables.  Physical Review A, 96(4):042332, 10 2017.
  * [9] Cosmo Lupo, Carlo Ottaviani, Panagiotis Papanastasiou, and Stefano Pirandola.  Continuous-variable measurement-device-independent quantum key distribution: Composable security against coherent attacks.  Physical Review A, 97(5):1–10, 2018.
  * [10] Christian Weedbrook, Andrew M. Lance, Warwick P. Bowen, Thomas Symul, Timothy C. Ralph, and Ping Koy Lam.  Quantum Cryptography Without Switching.  Physical Review Letters, 93(17):170504, 10 2004.
  * [11] Frédéric Grosshans, Nicolas J Cerf, Jérôme Wenger, Rosa Tualle-Brouri, and Philippe Grangier.  Virtual Entanglement and Reconciliation Protocols for Quantum Cryptography with Continuous Variables.  Quantum Information and Computation, 3(7):535–552, 2003.
  * [12] Fabian Laudenbach, Christoph Pacher, Chi-Hang Fred Fung, Andreas Poppe, Momtchil Peev, Bernhard Schrenk, Michael Hentschel, Philip Walther, and Hannes Hübel.  Continuous-Variable Quantum Key Distribution with Gaussian Modulation-The Theory of Practical Implementations.  Advanced Quantum Technologies, 1(1):1800011, 8 2018.
  * [13] Igor Devetak and Andreas Winter.  Distillation of secret key and entanglement from quantum states.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 461(2053):207–235, 1 2005.
  * [14] Stefano Pirandola, Stefano Mancini, Seth Lloyd, and Samuel L. Braunstein.  Continuous-variable quantum cryptography using two-way quantum communication.  Nature Physics, 4(9):726–730, 9 2008.
  * [15] Rui Lin and Jiajia Chen.  Minimizing spontaneous raman scattering noise for quantum key distribution in wdm networks.  In Optical Fiber Communication Conference, pages F4E–6. Optica Publishing Group, 2021.
  * [16] Chun Cai, Yongmei Sun, and Yuefeng Ji.  Intercore spontaneous raman scattering impact on quantum key distribution in multicore fiber.  New Journal of Physics, 22(8):083020, 2020.
  * [17] P. Eraerds, N. Walenta, M. Legré, N. Gisin, and H. Zbinden.  Quantum key distribution and 1 Gbps data encryption over a single fibre.  New Journal of Physics, 12, 6 2010.
  * [18] Robert W. Boyd.  Nonlinear Optics.  Academic Press, 4th edition edition, 2020.
  * [19] Q. Lin, F. Yaman, and Govind P. Agrawal.  Photon-pair generation in optical fibers through four-wave mixing: Role of Raman scattering and pump polarization.  Physical Review A - Atomic, Molecular, and Optical Physics, 75(2), 1 2007.
  * [20] Alan M. Hill and David B. Payne.  LINEAR CROSSTALK IN WAVELENGTH-DIVISION-MULTIPLEXED OPTICAL-FIBER TRANSMISSION SYSTEMS.  Journal of Lightwave Technology, LT-3(3):643–651, 1985.
  * [21] Irina Vorontsova, Roman Goncharov, Angelina Tarabrina, Fedor Kiselev, and Vladimir Egorov.  Theoretical analysis of quantum key distribution systems when integrated with a dwdm optical transport network.  JOSA B, 40(1):63–71, 2023.
  * [22] IO Vorontsova, RK Goncharov, AD Tarabrina, DV Tupyakov, EA Bolychev, SV Smirnov, FD Kiselev, and VI Egorov.  Theoretical analysis of measurement-device-independent quantum key distribution systems integrated into fiber-optic communication lines using dense wavelength division multiplexing.  Journal of Optical Technology, 89(7):424–429, 2022.
  * [23] AD Tarabrina, DV Tupyakov, IO Vorontsova, RK Goncharov, AV Zinov’ev, SV Smirnov, FD Kiselev, and VI Egorov.  Application of optimization methods to minimize noise in quantum key distribution systems integrated into fiber-optic communication lines using dense wavelength division multiplexing.  Journal of Optical Technology, 89(9):549–554, 2022.
  * [24] Zhengyu Li, Yi-Chen Zhang, Feihu Xu, Xiang Peng, and Hong Guo.  Continuous-variable measurement-device-independent quantum key distribution.  Physical Review A, 89(5):052301, 5 2014.


# Shadows and quasinormal modes of a charged non-commutative black hole by
different methods

Zening Yan1111Corresponding author<EMAIL_ADDRESS>İD , Xiaoji
Zhang2, Maoyuan Wan2 and Chen Wu3 1\. Department of Physics and Center for
Field Theory and Particle Physics, Fudan University, Shanghai 200433, China
2\. College of Science, University of Shanghai for Science and Technology,
Shanghai 200093, China
3\. Shanghai Advanced Research Institute, Chinese Academy of Sciences,
Shanghai 201210, China

(July 9, 2021)

###### Abstract

In this paper, we calculated the quasinormal modes (QNMs) of a charged non-
commutative black hole in scalar, electromagnetic and gravitational fields by
three methods. We gave the influence of non-commutative parameter $\theta$ and
charge $Q$ on QNMs in different fields. Thereafter, we calculated the shadow
radius of the black hole and provided the valid range of $\theta$ and $Q$
using the constraints on the shadow radius of $\text{M87}^{\ast}$ and
$\text{Sgr\;A}^{\ast}$ from the Event Horizon Telescope (EHT). In addition, we
estimated the “relative deviation” of the shadow radius ($\delta_{R_{s}}$)
between non-commutative spacetime and commutative spacetime. We found that the
maximum values of $\delta_{R_{s}}$ decreases with the increase of charge $Q$.
In other words, the non-commutativity of spacetime becomes harder to
distinguish as the charge of the black hole increases.

###### pacs:

04.70.Bw, 04.30.-w

##  I Introduction and motivation

In 1947, Snyder published the first paper on the non-commutativity of
spacetime Snyder (1947). Since then, non-commutative geometric black holes
have been studied more deeply. At present, there are two different methods to
study non-commutative quantum field theory: the Weyl-Wigner-Moyal
$\star$-product Groenewold (1946); Moyal (1949) and the other on coordinate
coherent state formalism Smailagic and Spallucci (2003a, b); Smailagic and
Spallucci (2004). The so-called “ $\star$-product ” method is that the
ordinary product between functions is replaced by the $\star$-product

$f(x)\star
g(x)\equiv\left.e^{\frac{i}{2}\theta^{\mu\nu}\frac{\partial}{\partial
x^{\mu}}\frac{\partial}{\partial\nu^{\nu}}}f(x)g(y)\right|_{y\rightarrow x},$
(1)

where $f(x)$ and $g(x)$ are two arbitrary infinitely differentiable functions.
It can be seen that the non-commutativity of spacetime is contained in
$\star$-product. The coordinate coherent state method means that the Dirac
$\delta$ distribution in ordinary spacetime is replaced by Gaussian
distribution. In the framework of this method, Nicolini believes that the non-
commutative effect only acts on the material source term and that there is no
need to change the Einstein tensor part of the field equation Nicolini et al.
(2006). Specifically, the mass density of the point-like function on the right
hand side of the Einstein field equation is replaced with the Gaussian smeared
matter distribution, while the left hand side of the equation is unchanged.

Nicolini first proposed the four-dimensional non-commutative - geometry -
inspired Schwarzschild black hole Nicolini et al. (2006). After that, this
black hole solution was extended to the case of charged Ansoldi et al. (2007).
In 2010, Modesto and Nicolini extended it to the general case of charged
rotating non-commutative black holes Modesto and Nicolini (2010). On this
basis, many researchers have explored the practical application of this black
hole, such as Douglas M. Gingrich, who explores the possibility of the
production and decay of non-commutative-geometry-inspired black holes on the
Large Hadron Collider (LHC) at the phenomenological level Gingrich (2010).
Chikun Ding and Jiliang Jing studied the influence of non-commutative
parameter on strong field gravitational lensing in the non-commutative
Reissner-Nordstr$\ddot{\text{o}}$m black-hole Ding and Jing (2011). In
addition, some authors use the observed data of gravitational waves and binary
pulsars to discuss the constraint of quantum fuzziness scale of non-
commutative spacetime and non-commutative gravity Kobakhidze et al. (2016);
Jenks et al. (2020).

On the other hand, the QNMs has always been an attractive topic in the
research field of black holes Konoplya and Zhidenko (2011); Berti et al.
(2009). The QNMs is the characteristic oscillations of the matter field under
the background of spacetime. The numerical result of the perturbation
frequency is a complex number, which has been found by the LIGO Scientific
Collaboration when detecting the gravitational wave generated by the merger of
two black holes Abbott et al. (2016). The reason why people are interested in
the QNMs stage is that it can be seen as a “characteristic sound” of a black
hole. This characteristic oscillation does not come from the matter inside the
black hole, but mostly depends on the spacetime metric outside the event
horizon of the black hole. This shows that spacetime itself is also a direct
participant in the oscillation.

The shadow of a black hole is a two-dimensional dark area formed on the
celestial sphere under the influence of a strong gravitational effect. In
2000, Falcke and his collaborators first proposed that the shadow of a black
hole could be observed Falcke et al. (2000). In 2019, the EHT first showed the
image of the supermassive black hole $\text{M87}^{\ast}$ Akiyama et al.
(2019a, b, c), from which shadows can be clearly observed. In 2022, the EHT
released an image of a supermassive black hole $\text{Sgr\;A}^{\ast}$ at the
center of the Milky way galaxy Akiyama et al. (2022a, b, c, d, e, f). In
general, the shadow cast by a non-rotating spherically symmetric black hole is
a circle, but for a rotating black hole, its shadow is similar to the shape of
the letter “D” Perlick and Tsupko (2022); Chael et al. (2021); Gralla et al.
(2019); Falcke et al. (2000).

The motivation of this paper is listed as follows: 1) In this paper, we will
test whether the corresponding relationship between eikonal perturbation and
shadow radius is valid in the case of non-vacuum Einstein’s equation solution;
2) In some papers, the WKB method fails to obtain correct numerical results
when calculating the QNMs of the non-commutative black hole. Therefore, we
will use other methods to calculate and give the correct numerical results of
QNMs in scalar, electromagnetic and gravitational fields; 3) We will use the
constraint range of shadow radius from EHT to limit the range of $\theta$ and
$Q$.

We organize the paper as follows: In Section 2, we introduce the basic
equations used to calculate the perturbation. And the valid range of $\theta$
corresponding to different $Q$ values is calculated in the case of ensuring
the existence of the event horizon of the black hole. In Section 3, the QNMs
of charged non-commutative black holes in different fields are calculated by
three methods. In Section 4, We test the relationship between shadow radius
and eikonal perturbation in charged non-commutative black hole spacetime, and
calculate the shadow radius of black hole corresponding to different
parameters. In Section 5, a brief summary of the whole paper is given. We use
natural units $\left(G=c=\hbar=1\right)$ throughout the paper.

##  II The basic equations

###  II.1 The metric of charged non-commutative black hole spacetime

When we define the spacetime coordinate $X^{\mu}$, non-commutative spacetime
can be expressed as

$\left[X^{\mu},X^{\nu}\right]=i\theta^{\mu\nu}=i\theta\epsilon^{\mu\nu}.$ (2)

The non-commutative parameter $\theta$ is a positive number and its dimension
is $length^{2}$, $\epsilon^{\mu\nu}$ is the anti-symmetric tensor Nicolini
(2009). We use the non-commutative spacetime metric constructed by coordinate
coherent state, that is, the “point particle” mass in spacetime is no longer
described by Dirac $\delta$ function, but by Gaussian function instead of
matter distribution. Because there is no point object, no physical distance
can be less than the minimum position uncertainty of the order of
$\sqrt{\theta}$ Spallucci et al. (2009).

The contours of the Gaussian distribution corresponding to matter and charge
density are

$\rho_{\text{matt}}(r)=\frac{\mathcal{M}}{(4\pi\theta)^{3/2}}\text{e}^{-\frac{r^{2}}{4\theta}},\quad\rho_{\text{el}}(r)=\frac{\mathcal{Q}}{(4\pi\theta)^{3/2}}\text{e}^{-\frac{r^{2}}{4\theta}},$
(3)

where $\mathcal{M}$ is the “bare mass” and $\mathcal{Q}$ is the total electric
charge. Here we need to consider a quasi-classical system, analogous to the
traditional Einstein-Maxwell system, which describes both electromagnetic and
gravitational fields Nicolini (2009)

$\mathcal{R}^{\mu}{}_{\nu}-\frac{1}{2}\delta^{\mu}{}_{\nu}\mathcal{R}=8\pi\left(\left.T^{\mu}_{\nu}\right|_{\text{matt}}+\left.T^{\mu}_{\nu}\right|_{\text{el}}\right),$
(4)
$\frac{1}{\sqrt{-g}}\partial_{\nu}\Big{(}\sqrt{-g}\delta^{0[\mu|}\delta^{r|\nu]}E(r)\Big{)}=J^{\mu},\quad
J^{\mu}(x)=4\pi\rho_{\text{el}}(r)\delta_{0}^{\mu}.$ (5)

The tensors $\left.T^{\mu}_{\nu}\right|_{\text{matt}}$ and
$\left.T^{\mu}_{\nu}\right|_{\text{el}}$ are the energy momentum tensors,
describing the matter and the electromagnetic content. $E(r)$ is the the
electric field and $J^{\mu}$ is the corresponding current density.

The line element for a static spacetime with spherical symmetry can be written
as follows:

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-A(r)dt^{2}+B(r)dr^{2}+D(r)\left(d\vartheta^{2}+\sin^{2}\vartheta
d\varphi^{2}\right),$ (6)

where the spherical coordinate $(t,r,\vartheta,\varphi)$ is used. This metric
satisfies $A(r)=B(r)^{-1}=f(r)$, $D(r)=r^{2}$, for the charged non-commutative
spacetime, and the lapse function Ansoldi et al. (2007); Nicolini (2009); Ding
and Jing (2011) is given by

$f(r)=1-\frac{4M}{r\sqrt{\pi}}\gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right)+\frac{Q^{2}}{r^{2}\pi}\left[\gamma^{2}\left(\frac{1}{2},\frac{r^{2}}{4\theta}\right)-\frac{r}{\sqrt{2\theta}}\gamma\left(\frac{1}{2},\frac{r^{2}}{2\theta}\right)+r\sqrt{\frac{2}{\theta}}\gamma\left(\frac{3}{2},\frac{r^{2}}{4\theta}\right)\right]\text{,}$
(7)

where the lower incomplete gamma function is written as

$\gamma(s,x)=\int_{0}^{x}t^{s-1}\mathrm{e}^{-t}\mathrm{~{}d}t\text{,}$ (8)

and notice that

$M=\frac{2\mathcal{M}}{\sqrt{\pi}}\Gamma\left(\frac{3}{2}\right)+\sqrt{\frac{2}{\theta}}\mathcal{Q}^{2}\Gamma\left(\frac{1}{2}\right),\quad
Q=\mathcal{Q}.$ (9)

When $r\rightarrow 0$ or $r\rightarrow\infty$, the limit of the lapse function
is

$\lim_{r\rightarrow 0\ \text{or}\ r\rightarrow\infty}f(r)\equiv 1\text{,}$
(10)

In order to explore the singularity of the charged non-commutative black hole
spacetime more rigorously, we calculated the Kretschmann scalar

$\mathcal{K}=\mathcal{R}_{\mu\nu\rho\sigma}\mathcal{R}^{\mu\nu\rho\sigma}.$
(11)

Then we show the variation of the Kretschmann scalar with parameters $\theta$
and $Q$ in Fig. 1. We can see from the figure that the Kretschmann scalar does
not diverge when $r\geq 0$. However, the Kretschmann scalar diverges when the
negative value $r$ approaches zero.

Figure 1: The function of the Kretschmann scalar with $r$ as the independent
variable. The parameter $M=1$ is selected.

When $M=1$ and $r$ is an infinitesimal positive number $(r\rightarrow 0)$, the
limit of the Kreichman scalar is

$\lim_{r\rightarrow
0}\mathcal{K}(r)=\frac{Q^{4}-4\sqrt{2\pi\theta}Q^{2}+8\pi\theta}{3\pi^{2}\theta^{4}}.$
(12)

As a result, we discover that there are no singularities in the $[0,\infty)$
range of the charged non-commutative black hole spacetime. Based on this, this
study only discusses the situation within the range of non-singularity.

In addition, Eq. (7) can be reduced to several well-known black hole solutions
when the non-commutative parameter $\theta$ and the parameter $Q$ take the
following limit cases:

$\left\\{\begin{aligned} \lim_{\theta\rightarrow
0}f(r)&=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\;&\text{Reissner-
Nordstr$\ddot{\text{o}}$m (Singularity)}\\\
f(r)|_{Q=0}&=1-\frac{4M}{r\sqrt{\pi}}\gamma\big{(}\frac{3}{2},\frac{r^{2}}{4\theta}\big{)}\;&\text{Non-
commutative-geometry-inspired (No singularity)}\\\ \lim_{\theta\rightarrow
0}f(r)|_{Q=0}&=1-\frac{2M}{r}\;&\text{Schwarzschild
(Singularity)}\end{aligned}\right.\text{.}$

It is well known that the charged non-commutative metric is a non-vacuum
static solution of the Einstein-Maxwell field equations.

Fig. 2 shows the change of the lapse functions of the charged non-commutative
black hole with parameters $\theta$ and $Q$.

Figure 2: The lapse functions of charged non-commutative black holes
corresponds to different parameters. The parameters $Q=0.2$ (left panel),
$\theta=0.12$ (right panel) and $M=1$ are selected.

Obviously, it can be seen that there is a critical value $\theta_{max}$ for
the non-commutative parameter $\theta$ to ensure that the black hole has at
least one event horizon on the left panel of Fig. 2.

As shown in Fig. 3, when a black hole has only one event horizon, the
corresponding critical value $\theta_{max}$ is calculated, implying that the
effective range of $\theta$ is $(0,\theta_{max})$. It can also be seen that
the critical value $\theta_{max}$ decreases as $Q$ increases.

Figure 3: The valid range of $\theta$ for various $Q$ values (cyan area) and
the critical value $\theta_{max}$ (red line). The parameter $M=1$ is selected.

###  II.2 Effective potential in the scalar, electromagnetic, and
gravitational fields

We consider the Klein-Gordon equation to solve the equation of massless
scalar:

$\qed\Psi=\frac{1}{\sqrt{-g}}\partial_{\mu}\big{[}\sqrt{-g}g^{\mu\nu}\partial_{\nu}\Psi\big{]}=0\text{,}$
(13)

and assume that the scalar field has the following form:

$\Psi(t,r,\vartheta,\varphi)=\frac{1}{r}\sum_{l,m}e^{-i\omega
t}\Phi(r)Y_{lm}(\vartheta,\varphi)\text{.}$ (14)

Substituting Eq. (14) into Eq. (13), we obtain the
Schr$\ddot{\text{o}}$dinger-like wave equation:

$\frac{\partial^{2}\Phi(r)}{\partial
r_{*}^{2}}+\big{[}\omega^{2}-V(r)_{scalar}\big{]}\Phi(r)=0\text{,}$ (15)

where tortoise coordinate is $r_{\ast}=\int dr/f(r)$. Therefore, the effective
potential of scalar field is

$V(r)_{scalar}=f(r)\left[\frac{l(l+1)}{r^{2}}+\frac{1}{r}f^{\prime}(r)\right]\text{,}$
(16)

where $l$ represents the multipole number.

We consider the evolution of a Maxwell field Cardoso and Lemos (2001). The
covariant equation is given by

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left[\left(A_{\sigma,\delta}-A_{\delta,\sigma}\right)g^{\delta\nu}g^{\sigma\mu}\sqrt{-g}\right]=0\text{.}$
(17)

The vector potential $A_{\mu}$ can be expanded into a four-dimensional vector
spherical harmonics

$A_{\mu}(t,r,\vartheta,\varphi)=\sum_{l,m}\begin{pmatrix}\begin{bmatrix}0\\\
0\\\
\frac{a^{lm}(t,r)}{\sin\vartheta}\partial_{\varphi}Y_{lm}(\vartheta,\varphi)\\\
-a^{lm}(t,r)\sin\vartheta\partial_{\vartheta}Y_{lm}(\vartheta,\varphi)\end{bmatrix}+\begin{bmatrix}f^{lm}(t,r)Y_{lm}(\vartheta,\varphi)\\\
h^{lm}(t,r)Y_{lm}(\vartheta,\varphi)\\\
k^{lm}(t,r)\partial_{\vartheta}Y_{lm}(\vartheta,\varphi)\\\
k^{lm}(t,r)\partial_{\varphi}Y_{lm}(\vartheta,\varphi)\end{bmatrix}\end{pmatrix}\text{,}$
(18)

the first term has parity $(-1)^{l+1}$ and the second ${(-1)}^{l}$.
Substituting Eq. (18) into Eq. (17), the wave equation can be obtained as
follows

$\frac{\partial^{2}\psi(r)}{\partial
r^{2}_{*}}+[\omega^{2}-V(r)_{ele}]\psi(r)=0\text{,}$ (19)

where $\psi(r)$ is a linear combination of functions $f^{lm}$, $h^{lm}$,
$k^{lm}$ and $a^{lm}$. The potential $V(r)_{ele}$ appearing in Eq. (19) is

$V(r)_{ele}=f(r)\left[\frac{l(l+1)}{r^{2}}\right].$ (20)

We consider the effective potential equations of gravitational perturbation
according to papers Kodama and Ishibashi (2003) and Lopez-Ortega (2006). Since
it is difficult to calculate the scalar type of gravitational perturbation,
only tensor and vector types of gravitational perturbation are calculated in
this paper. Their potential functions in four-dimensional spacetime are as
follows:

$V(r)_{grav(T)}=V(r)_{scalar}\text{,}$ (21)
$V(r)_{grav(V)}=f(r)\left[\frac{l(l+1)}{r^{2}}-\frac{r}{2}f^{\prime\prime\prime}(r)\right]\text{.}$
(22)

##  III The QNMs in different fields

###  III.1 Three calculation methods

The Mashhoon method is the P$\ddot{\text{o}}$schl-Teller potential
approximation method Poschl and Teller (1933), which uses the
P$\ddot{\text{o}}$schl-Teller potential $V_{PT}$ to approximate the effective
potential $V$ in the tortoise coordinate system.

The “asymptotic iteration method” (AIM) was applied to solve second order
differential equations for the first time in Ciftci et al. (2005). This method
was then used to obtain the QNMs of field perturbation in Schwarzschild black
hole Cho et al. (2010).

The Gundlach-Price-Pullin method, also known as the “time-domain integration
method” or “finite difference method” Gundlach et al. (1994). The time-domain
profile can be obtained by this method. Furthermore, we use the “least square
analysis” (LSA) to extract the QNMs in the time-domain profile. Specifically,
the LSA is used to fit the attenuated linear regression equation, and then the
approximate $\omega$ in time-domain profile are obtained.

There are specific explanations of this three methods in Yan et al. (2020,
2021), which will not be repeated here.

###  III.2 Analysis of numerical results

We use three different methods to calculate the QNMs of the charged non-
commutative black hole and give the numerical results corresponding to
different parameters in the scalar, electromagnetic and gravitational fields.
The WKB method is mostly used when calculating QNMs. However, due to the
complexity of the metric of charged non-commutative black holes, it is very
difficult to calculate with WKB. In addition, the WKB method will make serious
errors when calculating non-commutative black holes Yan et al. (2020, 2021).
Therefore, this paper chooses to use other methods for calculation, which can
not only ensure the accuracy of the results, but also greatly shorten the
calculation time.

Table 1: The QNMs of charged non-commutative black holes in scalar, electromagnetic and gravitational fields are calculated by the Msahhoon method, AIM and time-domain integration method. The parameters $M=1$, $\theta=0.1$, $l=3$ and $n=0$ are selected. $l$ is the multipole quantum number, $n$ is the overtone number. $Q$ | $\omega(\text{Mashhoon})$ | $\omega(\text{AIM})$ | $\omega(\text{Time-domain})$
---|---|---|---
| $V(r)_{scalar,grav(T)}$
0.0 | 0.678098 $-$ 0.0970906 $i$ | 0.675370 $-$ 0.0965014 $i$ | 0.675395 $-$ 0.0964935 $i$
0.1 | 0.679227 $-$ 0.0971423 $i$ | 0.676502 $-$ 0.0965541 $i$ | 0.676528 $-$ 0.0965400 $i$
0.2 | 0.682664 $-$ 0.0972957 $i$ | 0.679949 $-$ 0.0967101 $i$ | 0.679975 $-$ 0.0967086 $i$
0.3 | 0.688566 $-$ 0.0975444 $i$ | 0.685867 $-$ 0.0969634 $i$ | 0.685894 $-$ 0.0969463 $i$
0.4 | 0.697221 $-$ 0.0978746 $i$ | 0.694542 $-$ 0.0973012 $i$ | 0.694575 $-$ 0.0972955 $i$
| $V(r)_{ele}$
0.0 | 0.659686 $-$ 0.0962250 $i$ | 0.656900 $-$ 0.0956190 $i$ | 0.656924 $-$ 0.0956027 $i$
0.1 | 0.660805 $-$ 0.0962783 $i$ | 0.658023 $-$ 0.0956732 $i$ | 0.658046 $-$ 0.0956660 $i$
0.2 | 0.664212 $-$ 0.0964366 $i$ | 0.661440 $-$ 0.0958341 $i$ | 0.661464 $-$ 0.0958337 $i$
0.3 | 0.670065 $-$ 0.0966937 $i$ | 0.667310 $-$ 0.0960958 $i$ | 0.667334 $-$ 0.0960802 $i$
0.4 | 0.678654 $-$ 0.0970363 $i$ | 0.675921 $-$ 0.0964457 $i$ | 0.675948 $-$ 0.0964483 $i$
| $V(r)_{grav(V)}$
0.0 | 0.602419 $-$ 0.0933619 $i$ | 0.598781 $-$ 0.0927648 $i$ | 0.599109 $-$ 0.0923915 $i$
0.1 | 0.603694 $-$ 0.0934237 $i$ | 0.600072 $-$ 0.0929017 $i$ | 0.600379 $-$ 0.0924282 $i$
0.2 | 0.607583 $-$ 0.0936078 $i$ | 0.604036 $-$ 0.0933296 $i$ | 0.604256 $-$ 0.0926003 $i$
0.3 | 0.614294 $-$ 0.0939080 $i$ | 0.610970 $-$ 0.0940573 $i$ | 0.610945 $-$ 0.0928004 $i$
0.4 | 0.624215 $-$ 0.0943088 $i$ | 0.621410 $-$ 0.0948619 $i$ | 0.620833 $-$ 0.0930502 $i$
Table 2: Similar to Table 1, but the parameter $Q=0.4$ is selected. $\theta$ | $\omega(\text{Mashhoon})$ | $\omega(\text{AIM})$ | $\omega(\text{Time-domain})$
---|---|---|---
| $V(r)_{scalar,grav(T)}$
0.05 | 0.697221 $-$ 0.0978747 $i$ | 0.694550 $-$ 0.0973022 $i$ | 0.694587 $-$ 0.0972957 $i$
0.10 | 0.697221 $-$ 0.0978746 $i$ | 0.694542 $-$ 0.0973012 $i$ | 0.694575 $-$ 0.0973032 $i$
0.15 | 0.697226 $-$ 0.0978462 $i$ | 0.694316 $-$ 0.0971677 $i$ | 0.694245 $-$ 0.0971885 $i$
0.20 | 0.697348 $-$ 0.0972941 $i$ | 0.693019 $-$ 0.0960055 $i$ | 0.693263 $-$ 0.0959335 $i$
| $V(r)_{ele}$
0.05 | 0.678654 $-$ 0.0970364 $i$ | 0.675928 $-$ 0.0964458 $i$ | 0.675963 $-$ 0.0964321 $i$
0.10 | 0.678654 $-$ 0.0970363 $i$ | 0.675921 $-$ 0.0964457 $i$ | 0.675949 $-$ 0.0964380 $i$
0.15 | 0.678664 $-$ 0.0969893 $i$ | 0.675685 $-$ 0.0962744 $i$ | 0.675610 $-$ 0.0962781 $i$
0.20 | 0.678840 $-$ 0.0962754 $i$ | 0.674417 $-$ 0.0948607 $i$ | 0.674701 $-$ 0.0948216 $i$
| $V(r)_{grav(V)}$
0.05 | 0.624208 $-$ 0.0943465 $i$ | 0.621324 $-$ 0.0937073 $i$ | 0.621351 $-$ 0.0936956 $i$
0.10 | 0.624215 $-$ 0.0943088 $i$ | 0.621410 $-$ 0.0948619 $i$ | 0.620834 $-$ 0.0930569 $i$
0.15 | 0.625440 $-$ 0.0888988 $i$ | 0.623528 $-$ 0.0977381 $i$ | 0.624553 $-$ 0.0846251 $i$
0.20 | 0.641806 $-$ 0.0712106 $i$ | 0.631588 $-$ 0.0838247 $i$ | 0.647061 $-$ 0.0775725 $i$

It can be seen from Table 1 and Table 2 that these three methods are
effective, but the accuracy of QNMs calculated by the msahhoon method is low,
while the numerical results obtained by AIM and the time-domain integration
method are more accurate. Moreover, the numerical results of QNMs change
slightly with the increase of charge $Q$ and non-commutative parameter
$\theta$ in scalar, electromagnetic and gravitational fields. The change is so
small that it can almost be ignored. On the other hand, we also give the
profiles of perturbation evolution with time in scalar, electromagnetic and
gravitational fields corresponding to different parameters. As shown in the
Figs. 4 and 5. The case of a critical black hole is shown in Fig. 6. It can be
clearly seen that the damping of charged non-commutative black holes with time
has the stability of dynamical evolution under any parameters.

Figure 4: Time-domain profiles of $V(r)_{scalar,grav(T)}$ (left), $V(r)_{ele}$
(middle) and $V(r)_{grav(V)}$ (right). The parameters $M=1$, $Q=0.4$, $l=3$,
$n=0$, $v_{c}=10$ and $\sigma=3$ are selected. $v_{c}$ and $\sigma$ are the
parameters of the input initial Gaussian wave. Figure 5: Similar to Fig. 4,
but the parameter $\theta=0.1$ is selected. Figure 6: Similar to Fig. 4, but
the parameter $\theta=\theta_{max}$ is selected.

##  IV The shadow of a black hole

###  IV.1 The photon orbits follow the null geodesics in black hole spacetime

For the metric of Eq. (6), the Lagrangian of null geodesics is described by

$\mathcal{L}(x,\dot{x})=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}.$ (23)

Substituting Eq. (6) into Eq. (23), and for convenience, we only consider the
geodesic on the equatorial plane, that is, $\vartheta=\pi/2$. Therefore, we
obtain the following

$\begin{aligned}
\mathcal{L}(x,\dot{x})&=\frac{1}{2}\left(g_{00}\dot{t}^{2}+g_{11}\dot{r}^{2}+g_{22}\dot{\vartheta}^{2}+g_{33}\dot{\varphi}^{2}\right)\\\
&=\frac{1}{2}\left[-A(r)\dot{t}^{2}+B(r)\dot{r}^{2}+D(r)\dot{\vartheta}^{2}+D(r)\sin^{2}\vartheta\dot{\varphi}^{2}\right]\\\
&=\frac{1}{2}\left[-A(r)\dot{t}^{2}+B(r)\dot{r}^{2}+D(r)\dot{\varphi}^{2}\right]\\\
\end{aligned}.$ (24)

The motion of photons is described by the Euler-Lagrange equation

$\frac{d}{d\lambda}\left(\frac{\partial\mathcal{L}}{\partial\dot{x}^{\mu}}\right)-\frac{\partial\mathcal{L}}{\partial
x^{\mu}}=0,$ (25)

where $\lambda$ is the affine parameter, $\dot{x}^{\mu}$ represents the four-
velocity components of light ray. Substituting Eq. (24) into Eq. (25),

$\begin{aligned} &\frac{d}{d\lambda}\left(-A(r)\dot{t}\right)=0\\\
&\frac{d}{d\lambda}\left(B(r)\dot{r}\right)=0\\\
&\frac{d}{d\lambda}\left(D(r)\dot{\varphi}\right)=0\end{aligned}.$ (26)

Obviously, there exist two conserved quantities (constant of motion):

$E=-\frac{\partial\mathcal{L}}{\partial\dot{t}}=A(r)\dot{t},\quad
L=\frac{\partial\mathcal{L}}{\partial\dot{\varphi}}=D(r)\dot{\varphi},$ (27)

which are defined as energy and angular momentum of photons, respectively.

We use a first integral of the geodesic equation, namely
$\mathcal{L}(x,\dot{x})=0$ (for light). Therefore, we obtain the equation as
follows

$-A(r)\dot{t}^{2}+B(r)\dot{r}^{2}+D(r)\dot{\varphi}^{2}=0,$ (28)

Using $\dot{r}^{2}/\dot{\varphi}^{2}=(dr/d\varphi)^{2}$ and Eq. (27), the
orbit equation of the lightlike geodesic is given as follows:

$\left(\frac{dr}{d\varphi}\right)^{2}=\frac{D(r)}{B(r)}\left(\frac{D(r)}{A(r)}\cdot\frac{1}{b^{2}}-1\right),$
(29)

where $L/E=b$ is defined as the impact parameter, it is interpreted as the
vertical distance between the geodesic line and the parallel line passing
through the origin.

When a photon with a certain impact parameter $b$ is in the critical state of
being captured or escaping by a black hole, the photon is in an unstable
circular orbit near the black hole. Many circular orbits form a photon sphere
shell around the black hole. This critical photon sphere radius is defined as
$r_{c}$, and its corresponding critical impact parameter is defined as
$b_{c}$. The motion of light on the critical photon sphere must have a certain
value, so

$\dot{r}=0,\quad\ddot{r}=0,$ (30)

should be satisfied. According to Eq. (28), we give the expressions of
$\dot{r}$ and $\ddot{r}$ as follows:

$\dot{r}=\left(\frac{A(r)\dot{t}^{2}-D(r)\dot{\varphi}^{2}}{B(r)}\right)^{\frac{1}{2}},$
(31) $\begin{aligned} \ddot{r}&=\frac{\partial
r}{\partial\lambda}\left[\frac{\partial}{\partial
r}\left(\dot{r}\right)\right]=\dot{r}\left(\frac{\partial\dot{r}}{\partial
r}\right)\\\
&=\frac{\left(\dot{t}^{2}A^{\prime}(r)-\dot{\varphi}^{2}D^{\prime}(r)\right)B(r)-\left(\dot{t}^{2}A(r)-\dot{\varphi}^{2}D(r)\right)B^{\prime}(r)}{2B^{2}(r)}\\\
\end{aligned},$ (32)

where “ $\prime$ ” is the derivative of $r$. On the other hand, another
condition of Eq. (26) can be written as

$\frac{d}{d\lambda}\left(B(r)\dot{r}\right)=\frac{dr}{d\lambda}\left[\frac{d}{dr}\left(B(r)\dot{r}\right)\right]=B^{\prime}(r)\dot{r}^{2}+B(r)\ddot{r}=0,$
(33)

namely,

$\left(\dot{t}^{2}A(r)-\dot{\varphi}^{2}D(r)\right)B^{\prime}(r)+\left(\dot{t}^{2}A^{\prime}(r)-\dot{\varphi}^{2}D^{\prime}(r)\right)B(r)=0.$
(34)

Furthermore, Eq. (32) can be simplified to

$\ddot{r}=\frac{\dot{t}^{2}A^{\prime}(r)-\dot{\varphi}^{2}D^{\prime}(r)}{B(r)}.$
(35)

Therefore, Eq. (30) can be written as

$\left\\{\begin{array}[]{l}A(r)\dot{t}^{2}-D(r)\dot{\varphi}^{2}=0\\\
\dot{t}^{2}A^{\prime}(r)-\dot{\varphi}^{2}D^{\prime}(r)=0\\\
\end{array},\right.$ (36)

and according to Eq. (27), the equations expressed by variable $b$ can be
given as

$\left\\{\begin{array}[]{l}D(r)-A(r)b^{2}=0\\\
D(r)^{2}A^{\prime}(r)-D^{\prime}(r)A(r)^{2}b^{2}=0\\\ \end{array}.\right.$
(37)

By solving the above equation, the critical photon sphere radius $r_{c}$ and
critical impact parameter $b_{c}$ of any spherically symmetric black hole can
be obtained.

It can be seen that the values of $r_{c}$ and $b_{c}$ are independent of
$B(r)$, but the photon orbits are affected by $B(r)$ at other positions (not
in the critical photon orbit), as can be seen from Eq. (29) and Eq. (37).

###  IV.2 The relationship between eikonal perturbation and shadow radius

There is a deep relationship between the QNMs of a black hole and its shadow.
The relationship between the eikonal (short wavelengths or large multipole
number) perturbation and the shadow radius of black holes in static
spherically symmetric asymptotically flat spacetime is given in the paper
Cardoso et al. (2009). Recently, the relationship between the QNMs and the
shadow of a rotating black hole has been given in Jusufi (2020) and Yang
(2021), respectively.

Since scalar, electromagnetic and gravitational perturbations of high-
dimensional static black holes have the same behavior in the eikonal limit
Kodama and Ishibashi (2003); Ishibashi and Kodama (2003); Kodama and Ishibashi
(2004). Therefore, when the multipole number $l\rightarrow+\infty$, the
general form of the potential function Cardoso et al. (2009); Konoplya and
Stuchlík (2017); Churilova (2019); Moura and Rodrigues (2021) is expressed as

$V_{eik}(r)=\frac{f(r)}{r^{2}}l^{2}\text{.}$ (38)

In Einstein gravity, this limit can be applied to scalar, electromagnetic,
Dirac and all types (tensor, vector, scalar) of gravitational perturbations.
However, when considering the coupled with nonlinear electromagnetic field
Toshmatov et al. (2019); Toshmatov et al. (2018a, b) and the modified
gravitational theory perturbations Konoplya and Stuchlík (2017), there will be
different situations.

The radius of the orbit is determined by

$\left.\frac{\partial V_{eik}(r)}{\partial r}\right|_{r=r_{c}}=0\text{,}$ (39)

therefore, it needs to satisfy

$2f(r_{c})=r_{c}f^{\prime}(r_{c})\text{.}$ (40)

The solution $r_{c}$ of the above equation is defined as the radius of the
circular null geodesic, that is, the critical photon sphere orbit.

Under eikonal approximation, the QNMs of static spherically symmetric
asymptotically flat spacetime in any dimension can be expressed as

$\omega=l\Omega-i\left(n+\frac{1}{2}\right)\Lambda\text{.}$ (41)

The real part is determined by the angular velocity of the unstable null
geodesic $\Omega$, and the imaginary part is determined by the Lyapunov
exponent $\Lambda$. The expression Konoplya and Stuchlík (2017); Moura and
Rodrigues (2021) is as follows:

$\Omega=\frac{\sqrt{f\left(r_{c}\right)}}{r_{c}}\text{,}$ (42)
$\Lambda=\sqrt{-\frac{r_{c}^{2}}{2f\left(r_{c}\right)}\left[\frac{d^{2}}{dr_{\ast}^{2}}\left(\frac{f}{r^{2}}\right)\right]_{r=r_{c}}}=\sqrt{-\frac{r_{c}^{2}}{2}\left[f^{\prime}\left(\frac{f}{r^{2}}\right)^{\prime}+f\left(\frac{f}{r^{2}}\right)^{\prime\prime}\right]_{r=r_{c}}}\text{.}$
(43)

The shadow radius of a spherically symmetric black hole is given by

$R_{s}=\frac{1}{\Omega}=\frac{r_{c}}{\sqrt{f\left(r_{c}\right)}}\text{.}$ (44)

We can see that (44) and (40) combine to form Eq. (37).

In order to show this two-dimensional dark area, we use celestial coordinates
$(\text{X},\text{Y})$ Shaikh (2019); Lee et al. (2021) to describe it, namely

$R_{s}\equiv\sqrt{\text{X}^{2}+\text{Y}^{2}}\text{,}$ (45)

where

$\text{X}=\lim_{r_{o}\rightarrow\infty}\left[-\left.r_{o}^{2}\sin\vartheta_{o}\frac{d\varphi}{dr}\right|_{\left(r_{o},\vartheta_{o}\right)}\right]\text{,}\quad\text{Y}=\lim_{r_{o}\rightarrow\infty}\left[\left.r_{o}^{2}\frac{d\vartheta}{dr}\right|_{\left(r_{o},\vartheta_{o}\right)}\right]\text{,}$
(46)

and $(r_{o},\vartheta_{o})$ denotes the position coordinates of the observer
at infinity, $\vartheta_{o}$ is the angle between the rotation axis of the
black hole and the observer’s line of sight. The shadow of non-rotating
spherically symmetric black holes is not affected by $\vartheta_{o}$.

###  IV.3 Discussion and evaluation of results

First of all, Fig. 7 is obtained by calculating Eq. (29), which clearly shows
the schematic diagram of black hole shadow, that is, the critical impact
parameter $b_{c}$ is the shadow radius $R_{s}$.

Figure 7: The left panel shows the photon orbits around the Schwarzschild
black hole on the equatorial plane ($\vartheta_{o}=0$). Some photon
trajectories do not circle the equatorial plane once (blue and gray), some
trajectories circle the equatorial plane once (orange and magenta), and some
trajectories circle the equatorial plane twice (red). The right panel shows
the event horizon $R_{EH}=2$ (the boundary of the black disc), the critical
photon sphere orbit $r_{c}=3$ (the yellow dashed line) and the silhouette of
shadow $b_{c}=R_{s}=3\sqrt{3}$ Perlick and Tsupko (2022); Gralla et al. (2019)
(the boundary of the gray disc) from the perspective of the vertical
equatorial plane ($\vartheta_{o}=\pi/2$). The parameter $M=1$ is selected.

By solving the solution of Eq. (37) with metric Eq. (7), we can obtain $r_{c}$
and $b_{c}$ corresponding to different parameters $\theta$ and $Q$ in Table 3.
It can be seen that $r_{c}$ and $b_{c}$ decrease with the increase of $\theta$
or $Q$.

Table 3: Numerical results for $r_{c}$ and $b_{c}$ of charged non-commutative black holes. The parameter $M=1$ is selected. | $Q=0.0$ | $Q=0.4$
---|---|---
$\theta$ | $r_{c}$ | $b_{c}$ | $r_{c}$ | $b_{c}$
0.10 | 2.999999956 | 5.196152418 | 2.889244213 | 5.052977266
0.15 | 2.999955737 | 5.196145250 | 2.889131851 | 5.052957923
0.20 | 2.998727753 | 5.195884333 | 2.886728313 | 5.052416616
0.275811 | 2.980325376 | 5.190775292 | $\ast\ast\ast$ | $\ast\ast\ast$

Next, we use the critical photon orbit method mentioned in the previous
section to calculate the eikonal perturbation and compare the results with
other methods. Considering the time and accuracy of the calculation, we choose
the AIM to calculate and compare with it. In addition, the “relative
deviation” of the two methods is defined as

$\delta_{method}=\left|\frac{\omega(\text{AIM})-\omega(\Omega,\Lambda)}{\omega(\Omega,\Lambda)}\right|\text{.}$
(47)

The numerical results are shown in Table 4. It can be seen that when the value
of multipole number $l$ tends to be very large, the relative deviation
$\delta_{method}$ is very small, which shows that the numerical results
obtained by the two methods are very close, and it is also verified that the
eikonal perturbation calculated by the critical photon orbit method is
correct. On the other hand, it is well known that when $l\rightarrow+\infty$,
the imaginary part of QNMs tends to a fixed value, as also can be seen in Fig.
8.

Table 4: The eikonal perturbation of charged non-commutative black hole is calculated by AIM and critical photon orbit method, and the relative deviations of the two numerical results are compared. The parameters $M=1$, $\theta=0.1$, $Q=0.2$ and $n=0$ are selected. | QNMs | $\delta_{method}$
---|---|---
$l$ | $\omega(\text{AIM})$ | $\omega(\Omega,\Lambda)$ | $\text{Re}(\omega)$ | $\text{Im}(\omega)$
0 | 0.0803340 $-$ 0.254736 $i$ | 0.000000 $-$ 0.0964366 $i$ | $\infty$ | 1.64149
1 | 0.158243 $-$ 0.0881689 $i$ | 0.193752 $-$ 0.0964366 $i$ | 0.224395 | 0.0937712
2 | 0.370671 $-$ 0.0947315 $i$ | 0.387504 $-$ 0.0964366 $i$ | 0.0454122 | 0.0179993
5 | 0.962009 $-$ 0.0961490 $i$ | 0.968760 $-$ 0.0964366 $i$ | 0.00701761 | 0.00299119
10 | 1.93414 $-$ 0.0963649 $i$ | 1.93752 $-$ 0.0964366 $i$ | 0.00174755 | 0.000744047
50 | 9.68692 $-$ 0.0964337 $i$ | 9.68760 $-$ 0.0964366 $i$ | $7.01978\times 10^{-5}$ | $3.00725\times 10^{-5}$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
100 | 19.3749 $-$ 0.0964367 $i$ | 19.3752 $-$ 0.0964366 $i$ | $1.54840\times 10^{-5}$ | $1.03695\times 10^{-6}$
Figure 8: Time-domain profile of perturbation approaching the eikonal limit
($l\gg 1$). The parameters $M=1$, $\theta=0.1$, $Q=0.2$, $n=0$, $v_{c}=10$ and
$\sigma=3$ are selected.

And then, we obtain the shadow radius of some charged non-commutative black
holes, as shown in Figs. 9 and 10.

Figure 9: The radius of shadow corresponding to different non-commutative
parameter $\theta$. The parameters $Q=0$ (left panel), $Q=0.9$ (right panel),
$M=1$ and $\vartheta_{o}=\pi/2$ are selected. Figure 10: The radius of shadow
corresponding to different charge $Q$. The parameters $\theta=0.02$ (left
panel), $\theta=0.25$ (right panel), $M=1$ and $\vartheta_{o}=\pi/2$ are
selected.

It can be seen from Fig. 9 that when the black hole is uncharged and the
charge is large, the influence of the non-commutativity on the shadow radius
is not easy to distinguish. Therefore, the $R_{s}$ decreases slightly when the
$\theta$ increases. In Fig. 10, in the case of small non-commutativity, the
change of charge has an obvious effect on the shadow radius, which is due to
the large change range of charge in the case of small non-commutativity. As
shown in Fig. 3.

For spherically symmetric ($\tilde{r}_{\mathrm{sh}}$) and axisymmetric
($r_{\mathrm{sh},A}$) black holes, the EHT gives the constraint range of the
shadow radius of $\text{M87}^{\ast}$ Kocherlakota et al. (2021); Psaltis et
al. (2020) and the constraint range of the shadow radius of
$\text{Sgr\;A}^{\ast}$ from (Keck) and (VLTI) Akiyama et al. (2022f) as
follows:

$4.31M\approx r_{\text{sh,EHT-
min}}\leq\tilde{r}_{\mathrm{sh}},r_{\mathrm{sh},A}\leq r_{\mathrm{sh},\text{
EHT-max}}\approx 6.08M,$ (48) $\begin{aligned}
&4.5M\lesssim\tilde{r}_{\mathrm{sh}}\lesssim 5.5M,\quad\text{(Keck)}\\\
&4.3M\lesssim\tilde{r}_{\mathrm{sh}}\lesssim 5.3M,\quad\text{(VLTI)}\\\
\end{aligned},$ (49)

which has a confidence levels of $68\%$ and $G=c=1$ has been set. Therefore,
we can constrain the parameters of the charged non-commutative black hole by
using equations Eq. (48) and Eq. (LABEL:rangeofshadow2).

Figure 11: The values of $R_{s}$ with the change of $\theta$ when $Q$ takes
different values. The constraint range of $R_{s}$ by $\text{M87}^{\ast}$ (blue
area), Keck (green area) and VLTI (red area). The parameter $M=1$ is selected.

In Fig. 11, we can see that when $Q$ is close to the critical maximum, the
valid range of $\theta_{critical}$ will be less than $\theta_{max}$. In other
words, the valid range of $\theta$ will become smaller when $Q\rightarrow
Q_{max}$,

We define the “relative deviation” between the shadow radius corresponding to
different $\theta^{\prime}$ and the shadow radius in the commutative spacetime
as

$\delta_{R_{s}}=\left|\frac{R_{s}(\theta\rightarrow
0)-R_{s}^{\prime}(\theta=\theta^{\prime})}{R_{s}^{\prime}(\theta=\theta^{\prime})}\right|\text{,}$
(50)

where $\theta^{\prime}$ satisfies the range
$0<{\theta^{\prime}\leq\theta_{max}}$ corresponding to different $Q$ values.
Obviously, $\delta_{R_{s}}$ is a dimensionless quantity. Therefore, it can be
used to evaluate the deviation between non-commutative spacetime and
commutative spacetime at the level of shadow radius. The range of
$\delta_{R_{s}}$ is shown in Fig. 12.

Figure 12: The valid range of $\delta_{R_{s}}$ for various $Q$ values. The
parameter $M=1$ is selected.

Therefore, we can find that the range of $\delta_{R_{s}}$ decreases with the
increase of charge $Q$. In other words, the non-commutativity of spacetime is
more difficult to distinguish when the black hole carries more charge.

##  V Conclusions

In this paper, different methods are used to calculate the perturbation of
charged non-commutative black holes in scalar, electromagnetic and
gravitational fields, and these numerical results are analyzed in detail. In
addition, we have verified the relationship between the eikonal perturbation
and the shadow radius in the non-vacuum Einstein’s equation solution and
calculated its shadow radius. And then, we evaluated these results. We give
five conclusions as follows:

  1. 1.

In scalar, electromagnetic and gravitational fields, we obtain the accurate
numerical results of QNMs in the charged non-commutative black hole spacetime.
The results obtained by AIM and time-domain integration method are more
accurate than Msahhoon method, and the change of QNMs is almost be ignored
with the increase of charge $Q$ and non-commutative parameter $\theta$. In
addition, the change of damping with time can ensure the stability of
dynamical evolution under any parameter.

  2. 2.

We verified that the relationship between eikonal perturbation and shadow
radius is valid in the case of non-vacuum Einstein’s equation solution.

  3. 3.

We determined that the maximum value of $Q$, which comes from the constraint
conditions (48) and (LABEL:rangeofshadow2). And when $Q$ approaches the
maximum value, the effective range of $\theta$ will becomes smaller, that is,
$\theta_{critical}<\theta_{max}$, as shown in Fig. 11.

  4. 4.

The shadow radius $R_{s}$ of a black hole decreases slightly as the non-
commutative parameter $\theta$ increases, and the $R_{s}$ decreases as $Q$
increases. When the black hole carries more charge, the range of
$\delta_{R_{s}}$ value is smaller, which means that it is more difficult to
distinguish the non-commutativity of spacetime by shadow radius.

In addition, some papers have also discussed the shadow of non-commutative
black holes. The conclusions of these papers are listed: 1) The matter
accretion rate increases rapidly with the increase of non-commutative
parameter $\theta$ Gangopadhyay et al. (2018); 2) In the low-frequency limit,
the value of scattering/absorption cross section decreases with the increase
of non-commutative parameter $\theta$ Anacleto et al. (2020); 3) In paper
Sharif and Iftikhar (2016), the shadow of a rotating charged non-commutative
black hole is discussed, and it is found that the shadow is affected not only
by non-commutative parameter $\theta$ and charge $Q$, but also by spin $a$ and
the angle $\vartheta_{o}$. When the charge $Q$ increases, the deformation
parameter $\delta_{s}$ of the silhouette of shadow decreases, that is, the
shadow of the rotating black hole will maintain a circular shape with the
increase of the charge $Q$. The particle orbit is also affected with the
increase of $Q$.

## Acknowledgements:

The authors of this article thank developers of the AIM for their opening
codes.

## Data Availability Statement:

No data associated in the manuscript.

## References

  * Snyder (1947) H. S. Snyder, Phys. Rev. 71, 38 (1947).
  * Groenewold (1946) H. J. Groenewold, Physica 12, 405 (1946).
  * Moyal (1949) J. E. Moyal, Proc. Cambridge Phil. Soc. 45, 99 (1949).
  * Smailagic and Spallucci (2003a) A. Smailagic and E. Spallucci, J. Phys. A 36, L467 (2003a), eprint hep-th/0307217.
  * Smailagic and Spallucci (2003b) A. Smailagic and E. Spallucci, J. Phys. A 36, L517 (2003b), eprint hep-th/0308193.
  * Smailagic and Spallucci (2004) A. Smailagic and E. Spallucci, J. Phys. A 37, 7169 (2004), eprint hep-th/0406174.
  * Nicolini et al. (2006) P. Nicolini, A. Smailagic, and E. Spallucci, Phys. Lett. B 632, 547 (2006), eprint gr-qc/0510112.
  * Ansoldi et al. (2007) S. Ansoldi, P. Nicolini, A. Smailagic, and E. Spallucci, Phys. Lett. B 645, 261 (2007), eprint gr-qc/0612035.
  * Modesto and Nicolini (2010) L. Modesto and P. Nicolini, Phys. Rev. D 82, 104035 (2010), eprint 1005.5605.
  * Gingrich (2010) D. M. Gingrich, JHEP 05, 022 (2010), eprint 1003.1798.
  * Ding and Jing (2011) C. Ding and J. Jing, JHEP 10, 052 (2011), eprint 1106.1974.
  * Kobakhidze et al. (2016) A. Kobakhidze, C. Lagger, and A. Manning, Phys. Rev. D 94, 064033 (2016), eprint 1607.03776.
  * Jenks et al. (2020) L. Jenks, K. Yagi, and S. Alexander, Phys. Rev. D 102, 084022 (2020), eprint 2007.09714.
  * Konoplya and Zhidenko (2011) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011), eprint 1102.4014.
  * Berti et al. (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), eprint 0905.2975.
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), eprint 1602.03837.
  * Falcke et al. (2000) H. Falcke, F. Melia, and E. Agol, Astrophys. J. Lett. 528, L13 (2000), eprint astro-ph/9912263.
  * Akiyama et al. (2019a) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019a), eprint 1906.11238.
  * Akiyama et al. (2019b) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 875, L6 (2019b), eprint 1906.11243.
  * Akiyama et al. (2019c) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 875, L5 (2019c), eprint 1906.11242.
  * Akiyama et al. (2022a) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022a).
  * Akiyama et al. (2022b) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L13 (2022b).
  * Akiyama et al. (2022c) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L14 (2022c).
  * Akiyama et al. (2022d) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L15 (2022d).
  * Akiyama et al. (2022e) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L16 (2022e).
  * Akiyama et al. (2022f) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L17 (2022f).
  * Perlick and Tsupko (2022) V. Perlick and O. Y. Tsupko, Phys. Rept. 947, 1 (2022), eprint 2105.07101.
  * Chael et al. (2021) A. Chael, M. D. Johnson, and A. Lupsasca, Astrophys. J. 918, 6 (2021), eprint 2106.00683.
  * Gralla et al. (2019) S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100, 024018 (2019), eprint 1906.00873.
  * Nicolini (2009) P. Nicolini, Int. J. Mod. Phys. A 24, 1229 (2009), eprint 0807.1939.
  * Spallucci et al. (2009) E. Spallucci, A. Smailagic, and P. Nicolini, Phys. Lett. B 670, 449 (2009), eprint 0801.3519.
  * Cardoso and Lemos (2001) V. Cardoso and J. P. S. Lemos, Phys. Rev. D 64, 084017 (2001), eprint gr-qc/0105103.
  * Kodama and Ishibashi (2003) H. Kodama and A. Ishibashi, Prog. Theor. Phys. 110, 701 (2003), eprint hep-th/0305147.
  * Lopez-Ortega (2006) A. Lopez-Ortega, Gen. Rel. Grav. 38, 1565 (2006), eprint gr-qc/0605027.
  * Poschl and Teller (1933) G. Poschl and E. Teller, Z. Phys. 83, 143 (1933).
  * Ciftci et al. (2005) H. Ciftci, R. L. Hall, and N. Saad, Phys. Lett. A 340, 388 (2005), eprint math-ph/0504056.
  * Cho et al. (2010) H. T. Cho, A. S. Cornell, J. Doukas, and W. Naylor, Class. Quant. Grav. 27, 155004 (2010), eprint 0912.2740.
  * Gundlach et al. (1994) C. Gundlach, R. H. Price, and J. Pullin, Phys. Rev. D 49, 883 (1994), eprint gr-qc/9307009.
  * Yan et al. (2020) Z. Yan, C. Wu, and W. Guo, Nucl. Phys. B 961, 115217 (2020), eprint 2012.00320.
  * Yan et al. (2021) Z. Yan, C. Wu, and W. Guo, Nucl. Phys. B 973, 115595 (2021), eprint 2012.03004.
  * Cardoso et al. (2009) V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79, 064016 (2009), eprint 0812.1806.
  * Jusufi (2020) K. Jusufi, Phys. Rev. D 101, 124063 (2020), eprint 2004.04664.
  * Yang (2021) H. Yang, Phys. Rev. D 103, 084010 (2021), eprint 2101.11129.
  * Ishibashi and Kodama (2003) A. Ishibashi and H. Kodama, Prog. Theor. Phys. 110, 901 (2003), eprint hep-th/0305185.
  * Kodama and Ishibashi (2004) H. Kodama and A. Ishibashi, Prog. Theor. Phys. 111, 29 (2004), eprint hep-th/0308128.
  * Konoplya and Stuchlík (2017) R. A. Konoplya and Z. Stuchlík, Phys. Lett. B 771, 597 (2017), eprint 1705.05928.
  * Churilova (2019) M. S. Churilova, Eur. Phys. J. C 79, 629 (2019), eprint 1905.04536.
  * Moura and Rodrigues (2021) F. Moura and J. a. Rodrigues, Phys. Lett. B 819, 136407 (2021), eprint 2103.09302.
  * Toshmatov et al. (2019) B. Toshmatov, Z. Stuchlík, B. Ahmedov, and D. Malafarina, Phys. Rev. D 99, 064043 (2019), eprint 1903.03778.
  * Toshmatov et al. (2018a) B. Toshmatov, Z. Stuchlík, J. Schee, and B. Ahmedov, Phys. Rev. D 97, 084058 (2018a), eprint 1805.00240.
  * Toshmatov et al. (2018b) B. Toshmatov, Z. Stuchlík, and B. Ahmedov, Phys. Rev. D 98, 085021 (2018b), eprint 1810.06383.
  * Shaikh (2019) R. Shaikh, Phys. Rev. D 100, 024028 (2019), eprint 1904.08322.
  * Lee et al. (2021) B.-H. Lee, W. Lee, and Y. S. Myung, Phys. Rev. D 103, 064026 (2021), eprint 2101.04862.
  * Kocherlakota et al. (2021) P. Kocherlakota et al. (Event Horizon Telescope), Phys. Rev. D 103, 104047 (2021), eprint 2105.09343.
  * Psaltis et al. (2020) D. Psaltis et al. (Event Horizon Telescope), Phys. Rev. Lett. 125, 141104 (2020), eprint 2010.01055.
  * Gangopadhyay et al. (2018) S. Gangopadhyay, R. Mandal, and B. Paik, Int. J. Mod. Phys. A 33, 1850084 (2018), eprint 1703.10057.
  * Anacleto et al. (2020) M. A. Anacleto, F. A. Brito, J. A. V. Campos, and E. Passos, Phys. Lett. B 803, 135334 (2020), eprint 1907.13107.
  * Sharif and Iftikhar (2016) M. Sharif and S. Iftikhar, Eur. Phys. J. C 76, 630 (2016), eprint 1611.00611.


# On the Use of Dimension Reduction or Signal Separation Methods for Nitrogen
River Pollution Source Identification

Hatipoğlu, Güray

(April 27, 2022)

## Abstract

Identification of the current and expected future pollution sources to rivers
is crucial for sound environmental management. For this purpose numerous
approaches were proposed that can be clustered under physical based models,
stable isotope analysis and mixing methods, mass balance methods, time series
analysis, land cover analysis, and spatial statistics. Another extremely
common method is Principal Component Analysis, as well as its modifications,
such as Absolute Principal Component Score. they have been applied to the
source identification problems for nitrogen entry to rivers. This manuscript
is checking whether PCA can really be a powerful method to uncover nitrogen
pollution sources considering its theoretical background and assumptions.
Moreover, slightly similar techniques, Independent Component Analysis and
Factor Analysis will also be considered.

##  1 Introduction

The excess nitrogen load into aquatic environment result in unwanted
ecological conditions, as well as nitrogenous greenhouse gas emission and
valuable reactive nitrogen source. Among the mediums, rivers are very
significant. They usually do not have very high capacity to accumulate
nitrogen to excessive concentrations by themselves, yet they can carry it to
other environments, coastal or inland, to their predicament (Galloway, 2005).
Governing this issue requires working knowledge over how nitrogen behaves in
the environment and what are the loads of nitrogen to a specific region.

Nitrogen can reside in all mediums. A summary of river relevant inorganic
nitrogen species is given in the following table.

tableRelevant inorganic nitrogen species for river environments  Name Formula
Molecular Mass (g/mol) Ammonium NH4+ 17 Dinitrogen N2 28 Nitrate NO3- 62
Nitrite NO2- 46 Nitrogen Dioxide NO2 78 Nitrogen Monoxide NO 62 Nitrous Oxide
N2O 44

Other relevant environmental measurements for nitrogen are Total Nitrogen (N
(from ammonia, N-NH4) + N-NO3 + N-NO2 (from nitrite) + organic N) and Total
Kjeldahl Nitrogen (Organic N + N-NH4). They will be referred at the later
parts of this manuscript again.

Atmosphere mostly has inert dinitrogen, where biological or technological
fixation can make it available to plants and others, in addition to the
lightning bolts. Atmosphere also transports and deposits land-generated
combustion byproduct nitrous oxides, evaporated ammonia, denitrified nitrous
oxide and dinitrogen as well as their counterparts from aquatic environments.
Very high majority of the plants can only utilize nitrate and partially
nitrite and ammonium, while only a fraction of photosynthetic cyanobacteria
and several other bacteria type can fix dinitrogen directly. Consequently,
ammonium and nitrate concentrations are monitored to estimate state of the
aquatic ecosystem and possible eutrophication or harmful algal bloom
conditions. More complex transformation of these chemical species can be
simulated in biogeochemical models to understand nitrogen fate and transport
better.

Figure 1: Main ways of entry to rivers by nitrogen species. From left top to
right top; atmospheric wet deposition, sewage discharge, urban runoff,
atmospheric dry deposition, upstream river flow. From bottom right to bottom
left; soil decomposition, industrial wastewater discharge, agricultural runoff
(manure), agricultural runoff (fertilizer)

The part above introduced the basic movements of nitrogenous species through
the environment. In the following section, a close-up view to riverine
dynamics of nitrate will be examined to permit the selection of a better
algorithm to make sense of collected in situ data.

##  2 What is Happening While Nitrogen is Polluting Rivers?

The nitrogen dynamics in rivers is simply too complex to satisfactorily
include here as a whole. However, the critical parts will be tried to be
included to be able to select the necessary statistical tool for later
analyses of source apportionment. For a full investigation of global nitrogen
cycle, the readers are kindly referred to Ussiri and Lal (2013).

One complication of nitrogen dynamics in river is that many other seemingly
unrelated components can directly or indirectly impact its fate. Phytoplankton
and cyanobacteria both consumes and generates nitrogen species while they are
affected by many other environmental conditions, including but not limited to
salinity, temperature, phosphorus. Their numbers, species and spatial
distributions vary with time even without considering their predators, such as
zooplanktons, which further complicates the situation. Changes in river
dynamics, flows, surface run-off and groundwater inflows to rivers also
impacts both concentration and amount of nitrogen species. As a result, many
water quality variables and physical characteristics are interdependent, hence
if not impossible, not causally related variable set generation is very
difficult, let alone simple covariance/correlation. Henceforth, the first
remark is, the approaches assuming that water quality and physical
characteristics are independent from each other are not suitable for nitrogen
source apportionment. Another issue is spatial and temporal correlation
between nitrate measurement values of samples collected throughout the study
area in rivers. Many statistical analyses and machine learning routines work
much better with larger data volumes. However, the samples fed to the
algorithms are in some degree not independent from each other. In line with
intuition, upstream nearby sampling point is similar to the downstream ones.
From the temporal perspective, there are daily, seasonally, and yearly cycles
in many nitrogen emitting processes. There are several ways to check and act
according to these intersample correlations. For instance, Moran’s I value
related calculations have been used to determine minimum distance between
river sampling network that will result in spatially uncorrelated locations.
However, the fact that many point and non-point pollution sources enter the
river from different points, and consequent differences in different set of
nitrogen sources in different locations should be kept in mind. For temporal
analyses, there are various decomposition techniques to get rid of temporal
autocorrelations. Yet, one critical point is the complexity of source
apportionment will be higher if the proportions of those sources on their
nitrogen contribution to the rivers changes significantly with time. Another
possible way is taking a ”temporal snapshot” while considering spatial
autocorrelation as well, and make separate analyses, with obvious trade-offs
from data volume. In short, the individual measurement points, or simply
sample measurements’ independency from each other should be thoroughly checked
prior to the analyses, especially the analyses assuming that samples are
independent.
Yet another point is related to the sources and sinks. Moreover, there are
several factors that affect different sources in varying degrees. The
following examples will try to clarify these statements:
Urban runoff and sewage discharge are both related to the size of city and
consequently to population casually. Agricultural runoff from manure and
fertilizer can be correlated since manure is also being used as a type of
organic fertilizer. The nitrogen from the atmospheric deposition is, in
addition to the combustion and lightning bolts, sourced by nearby agricultural
nitrogen emissions as well.
For factors affecting many of these sources, the first one is changes in
precipitation. Surface runoff from urban or agricultural environments are
either creathed or enhanced by precipitation. Wastewater discharges can be
diluted by precipitation, or after exceeding the capacity of treatment
facilities it might result in untreated wastewater discharge. Wet deposition
is directly coming with precipitation. The second factor is temperature. If it
results in a decline in relative humidity, crops will draw more water, hence
fertilizer loss might be lower. On the other hand, overall increase in aerobic
bacterial activity will decompose soil organic matter faster, hence more
nitrogen loss from this source. Additionally, increase in water temperature
will also increase denitrification, which will reduce nitrate level overall in
water. Last factor is related to environmental management in a municipal or
regional level. Urban runoffs, sewerage, wastewater treatment plants and other
relevant infrastructures may or may not be robust enough to constrain the
adverse impacts of anthropogenic pollution sources. Henceforth, pollution
sources (as well as sinks) can be directly or indirectly related to each
other, they are not independent.
Nevertheless, a number of techniques has been developed to date with non-
trivial successes that can still be able to predict the nitrogen pollution in
the midst of all this complexity. The following three sections will explain
three relatively similar methods, Principal Component Analysis, Independent
Component Analysis, and Factor Analysis one by one with now their assumptions.
In this way, the suitability conditions of these methods for the case of
nitrogen pollution source apportionment will be discussed.

##  3 What Does Principal Component Analysis (PCA) Do?

PCA is generally known as a dimensional reduction method or a way to rotate
axises to maximize variance. To comprehend the theory behind the PCA and its
advantages/shortcomings, readers are referred to the tutorial of Shlens (2014)
and book of Everitt and Hothorn (2011). Several points of PCA are considered
below:
Covariance matrix construction
Assume that we have a dataset with N samples and M variables. Firstly, mean of
each variable will be subtracted from its corresponding sample value for every
sample. After this, how the variables covary with each other will be examined
via covariance matrix.

$\begin{bmatrix}\sigma_{1}\sigma_{1}&\sigma_{1}\sigma_{2}&\sigma_{1}\sigma_{3}&\sigma_{1}\sigma_{4}\\\
\sigma_{2}\sigma_{1}&\sigma_{2}\sigma_{2}&\sigma_{2}\sigma_{3}&\sigma_{2}\sigma_{4}\\\
\sigma_{3}\sigma_{1}&\sigma_{3}\sigma_{2}&\sigma_{3}\sigma_{3}&\sigma_{3}\sigma_{4}\\\
\sigma_{4}\sigma_{1}&\sigma_{4}\sigma_{2}&\sigma_{4}\sigma_{3}&\sigma_{4}\sigma_{4}\end{bmatrix}$

These are standard deviations of i and y variables, and their multiplication
results in covariance of i and y variables. If two variables vary together the
multiplication of their standard deviations will be nearer to 1, and nearer to
0 in the opposite case. As this is a very powerful way to check the overall
situation between variables and that many statistical methods have this in
their processes, R Statistical software, MATLAB and Python languages have a
single function to create it. After this step, creation the correlation matrix
of the same dataset is easy, permitting to check the traces of
multicollinearity among variables. To the covariance matrix (or correlation
matrix) a matrix decomposition method, eigenvalue decomposition can be
applied. The entire dataset will be rotated from the original variables to a
new space with principal components ordered from maximum explanation of
variance in data to the lowest one, which is in general noise in the data.
Each principal component is a linear combinations of variables with different
weights
The actual pollution source apportionment part comes now. Numerous research
has been done via directly interpreting these different principal components
as sources that vary each variable in a certain degree, with high weights
correspond to higher dominance by a certain source. The main intuition behind
this approach is that if a group of variables are strongly and significantly
covary according to a single pirncipal component, then they are ought to be
impacted by a certain source that simultaneously introduce them or
directly/indirectly affect them together. Everitt and Hothorn found this as an
overinterpretation, as PCA is not and does not do what is directly expected
from a source apportionment algorithm. Albeit many researchers reported quite
high variance explanation and well-sounding principal component to source
elaborations, PCA is not itself appropriate for non-linear cases, and sources
and changes in nitrogen dynamics are frequently non-linear in nature, with
mentioned spatio-temporal autocorrelations. To encounter these,
1-) many other modified versions of PCA (kernel PCA etc.) might be used. An
introduction to a comprehensive discussion over PCA itself and comments on
this discussion can be found in Joliffe’s book (2002), or
2-) variables might be aggregated (see Total Nitrogen above) or discarded, or
3-) other methods that are insensitive against mentioned limitations might be
selected, or
4-) samples might be reduced to get rid of some autocorrelations, or
decomposed, or
5-) in case temporal data is used, for instance on a single station, rather
than the use of actual nitrate concentration values, difference of nitrate
concentration between two consecutive time-steps might be used.
These considerations (except 3) will be further discussed in the following
application seciton, after considering ICA and FA approaches.

##  4 Can Independent Component Analysis (ICA) be More Appropriate?

Independent Component Analysis, as its name implies, looks for independent
components and it is directly being used for source apportionment, blind
source separation. Readers are kindly referred to Shlens (2014b) tutorial over
how it works, theoretical basis and comparisons with PCA.
The base under the ICA is Singular Value Decomposition.

$A=UAV^{T}$ (1)

First half of the ICA methdology is as if doing PCA and normalizing the
variances to have dimensions in standard units. This removes linear
dependencies. Later, statistical independence will be ensured via putting a
constraint on data from information theory. Here, the aim is making the
statistical dependence near-zero:

$I(y)=\int{P(y)log_{2}\frac{P(y)}{\prod_{i=i}P(y_{i})}dy}$ (2)

This very strict criterion also result in removal of higher-order dependencies
in the dataset compared to the PCA. The critical point is because of the
”requirement” of non-gaussian distributions in the dataset, there are local
minima in the independent component searching second part. Infact, ICA
produces similar results to PCA when the distribution of data is Gaussian. The
readers are referred to the book of Independent Component Analysis by
Hyvarinen, Karhunen and Oja (2001)’s page 161 for the caveat of this non-
Gaussian criterion.

##  5 Is This Topic More Suited to be Studied by Exploratory Factor Analysis?

FA simply analyzes the factors, i.e. latent variables, over the dataset. The
latent variable is hard to directly measure, while its effects can be measured
and later related to latent variable quantitatively and linearly. For the
theoretical basis and examples on its use, the Exploratory Factor Analysis and
Confirmatory Factor Analysis and Structural Equation Models chapters of
Everitt and Hothorn (2011), and Adachi (2016) are recommended. Several caveats
in the factor analysis as follow: Factor analysis can be made via ”principal
factor analysis”, similar to PCA but focusing on common factors and tries the
find the variance that is related to other variables as well. Another way to
achieve it is maximum likelihood method.
Overall, what FA assumes in the actual situartion is below: Let A, B, C, D, E
are variables, and Factor1 and Factor2 are latent variables, or factors, ”c”s
are intercepts of the equations and ”e”s are error terms.

$\displaystyle A=\textit{a}_{11}\times Factor_{1}+\textit{a}_{12}\times
Factor_{2}+c_{1}+e_{1}$ $\displaystyle B=\textit{a}_{21}\times
Factor_{1}+\textit{a}_{22}\times Factor_{2}+c_{2}+e_{2}$ $\displaystyle
C=\textit{a}_{31}\times Factor_{1}+\textit{a}_{32}\times
Factor_{2}+c_{3}+e_{3}$ $\displaystyle D=\textit{a}_{41}\times
Factor_{1}+\textit{a}_{42}\times Factor_{2}+c_{4}+e_{4}$ $\displaystyle
E=\textit{a}_{51}\times Factor_{1}+\textit{a}_{52}\times
Factor_{2}+c_{5}+e_{5}$

The error term e has its own variance for each variable, while the Factors are
common with all variables with only different weights. This separation is
quite important as the FA look for common factors/variances and not the unique
variance of each variable

Determining the number of factors to be retrieved is a priori, and the same
dataset will generate different factors and variable loadings with changing
number of factors. Slightly similar to the k-means clustering, one way is
starting form a low number of factor and keep increasing the number until the
minimum factor’s significance is low enough.
As it will be discussed in the application part, factors hypthetically may
have common variances owing to not only from estimated pollution sources, but
also from meterological and socioeconomic factors, too.

##  6 A Case Study for Applying PCA, ICA, FA to Identify Nitrogen Pollution
Sources: Sacramento River

After describing overall characteristics of source apportionment of nitrogen
pollution and several statistical methods used for source separation, in this
chapter I will apply these algorithms to a real dataset collected by USGS from
Sacramento River’s site SACRAMENTO R A FREEPORT CA coded with USGS-11447650
(USGS, 2022). This site is chosen as it has relatively high number of
historical data in terms of different variables, and for future works there
are ancillary data regarding point source pollution, pesticide use nearby,
land use, additional water quality analyses campaigns, atmospheric deposition
data.
For now, only one location is chosen to practically eliminate spatial
autocorrelation impact. 32 parameters are filtered having more than 400 data
points between year 1958 and 2022, collected from surface water (WS code). The
data obtained via ”dataRetrieval” package in R (R Core Team, 2022) further
filtered with ”not having NA for 00618” (nitrate related parameter) and later
removing all rows with NA values. This resulted in 124 rows and 30 variables
spanning from October 22,1985 to September 9, 2000. One last and significant
step before the analyses is temporal preprocessing. The aim is getting rid of
irregularly spaced data problem and removing the impact of temporal
correlation. Annual means are taken to mask the impacts of monthly and
seasonal dynamics. The variables with ”NA” values after annual mean operation
were dropped, resulted in remaining 17 variables. Deleting the variables that
are directly related to each other (such as ”Nitrate nitrite” in the presence
of Nitrate and Nitrite) the list further trimmed down to 11 variables.
Autocorrelation function for the variables in the temporal dimension did not
point out major relations. Hence, this dataframe will be entered to PCA,ICA,FA
methods.
In summary, the sources of changes in yearly differences on nitrogen water
quality variables are sought after. The reason of the inclusion of non-
nitrogen species and measurements are that they can assist interpretation.

###  6.1 PCA

To the difference data.frame with 38 rows and 11 variables, PCA was applied
with prcomp function of R. The input data is centered, as well as scaled since
there are many different scales within the variable set. Principal Components
(PCs) with standard deviations higher than 1 were retained, resulted in 4 PCs.

tablePC loadings on selected variables (Un stands for unfiltered samples
including solid fraction in the analysis, capital Ns at the and and Po4 means
values given in ”as nitrogen” and ”as phosphate”) Variables PC1 PC2 PC3 PC4
Water Temperature -0.1755410 -0.38984238 -0.01303532 0.2804046 Dissolved
Oxygen 0.1187128 0.27379533 0.31632429 -0.4363757 pHUn -0.4619997 -0.26795873
0.11140634 -0.2611894 CarbonDioxide 0.2833118 0.38008688 -0.14124896 0.4847708
Organic N -0.2860378 0.02734127 -0.53273642 -0.3712082 Ammonia -0.3408886
-0.09170747 0.02860327 0.3678343 Nitrite -0.1451199 0.13136421 0.30351669
-0.2083277 Nitrate N -0.1543979 0.59988497 -0.06531265 -0.1205199 Phosphorus
Po4 -0.4482934 0.14167495 0.22421733 0.1691096 PhosphorusUn -0.2950863
0.23006018 -0.56526880 0.0667281 Chloride -0.3611361 0.31227442 0.34145247
0.2522957 Total stdev 1.7091970 1.4538088 1.2109395 1.0650152

Reiterating briefly, PCA aims to find out if there is any interesting patterns
in the in the loadings of the variables (annual average concentrations) under
different principal components. The four PCs retained according to the Kaiser
criterion together explains 64 percent of the variation in this dataset, which
looks relatively low but there is no expected spatial and temporal correlation
processed dataset. And considering the last row of Table. 2, already low
standard deviations levels off after PC3, indicating that it would be an elbow
in a scree plot. When the results are examined while keeping aforementioned
situation in mind, PC1 has, in decreasing explanatory power had high loadings
for Phosphorus as phosphate, unfiltered pH, chloride and ammonia all
negatively. PC2 has distinct positive loading for Nitrate N followed by
comparatively insignificant loadings. PC3 has relatively high loading of
nitrite and negative loading in organic N and unfiltered phosphorus. From this
much of information, it may be inferred that PC1 might reflect the variation
in domestic wastewater quality or entry to the system and/or phytoplankton
filtering of the water , while PC2 is a process introducing nitrate, which
might be agricultural fertilizer runoff, and PC3 might be a phytoplankton
existence as well with decrease in total phosphorus and increase in dissolved
oxygen

###  6.2 ICA

For conducting ICA, ”fastICA” package in R was used with changing number of
sources 2 to 5, among them 3 did not converge to a result. The parameters used
were: maximum iteration of 200 and tolerance of 0.0001, no standardization,
method ”C”, simultaneous extraction (alg.typ = ”parallel”), and logcosh
functional form. The following Figure 2 shows the independent components
together. The ICA algorithm illustrated a way to generate the variances in
this study’s dataset with 3 independent components. Curiously enough all these
components have spikes and dips in different places. Let’s first check the
nitrate and ammonia annual average values and find out if there is any change
in phase.

Figure 2: The Independent Components uncovered by ICA with 3 components Figure
3: Annual average time series of nitrate and ammonia

Direct estimation of nitrate, ammonia or their year-to-year differences
linearly with these components fails. Another approach is trying to go from
known contaminators to these components. Nitrogen atmospheric deposition in
National Trend Network CA88 station has similar iin-phase trends to the second
Independent Component. Schlegel and Domagalski’s peer reviewed study (2015)
has also total N fertilizer sold in Sacramento Basin between 1987-2010, which
again has some spikes reminiscent of Second ICA here. Nevertheless, there are
other even harder to quantify nitrogen sources, such as soil decomposition,
urban runoff; and we have quite sharp uncompensateed drops in ICs for many
years, which require further data/analysis, and interpretation to infer
something concrete.

###  6.3 FA

”factanal” function of R is used to conduct exploratory factor analysis on
this difference data. Factor analysis was done with 1,2,3, 4 and 5 factors
first, after getting the result from PCA with no rotation and maximum
likelihood estimation as the method to find factors. The ”p-value” for the
resulting FAs follow as in order from FA for 1 factor to FA for 5 factor:
1.17e-06, 0.0321, and 0.113. Although 2-factor analysis was chosen to further
examine the outputs, the element sof the residual matrix were too distant to
being zero so the entire FA undertaking was deemed unsuccessful to explain the
underlying factors impacting the changes in annual average values of the
selected variables.

##  7 Concluding Remarks

In this paper, it was sought to approach to the source apportionment of river
nitrogen pollution from many perspectives of PCA and similar analyses. As
suggested by several authors when PCA does not seem to be enough, ICA looks
promising to solve this problem as well. The next steps are the inclusion of
spatial and temporal dimensions, time-series analysis and the so-called
temporal-PCA method: Singular Spectrum Analysis. This and following studies
aim to include more rigorous theretical basis on the use of these dimension
reduction and source separation technics for nitrogen river pollution source
apportionment task.

##  8 References

B. Everitt, T. Hothorn. An Introduction to Applied Multivariate Analysis with
R. Springer. (2011)
B. Schlegel, J.L. Domagalski. San Francisco Estaury and Watershed Science, 13
(4). (2015)
D. Ussiri, R. Lal. Global Nitrogen Cycle, chapter in ”Soil emission of Nitrous
Oxide and Its Mitigation 29-62. (2012)
I. T. Joliffe. Principal Component Analysis. Springer. 2nd Edition. (2002).
J. N. Galloway, Science in China Series C: Life Sciences, 48, 669-678. (2005)
J. Shlens. A Tutorial on Principal Component Analysis. (2014). arXiv:1404.1100
J. Shlens. A Tutorial on Independent Component Analysis. (2014).
arXiv:1404.2986
K. Adachi. Matrix Based Introduction to Multivariate Data Analysis. Springer.
(2016)
R Core Team. R: Language and Environment for Statistical Computing. R
Foundation for Statistical Computing. (2022).
USGS. USGS 11447650 Sacramento R A Freeport CA.
url:https://waterdata.usgs.gov/ (2022).