[1]Matthew Hayman [1]Earth Observing Lab, National Center for Atmospheric
Research, PO Box 3000, Boulder, 80307, CO, USA 2]Space Science and Engineering
Center, University of Wisconsin at Madison, 1225 West Dayton St, Madison,
53726, WI, USA

# Global Estimation of Range Resolved Thermodynamic Profiles from MicroPulse
Differential Absorption Lidar

<EMAIL_ADDRESS>Robert A. Stillwell<EMAIL_ADDRESS>Adam Karboski
<EMAIL_ADDRESS>Willem J. Marais<EMAIL_ADDRESS>Scott M.
Spuler<EMAIL_ADDRESS>*  [

###### Abstract

We demonstrate thermodynamic profile estimation with data obtained using the
MicroPulse DIAL such that the retrieval is entirely self contained. The only
external input is surface meteorological variables obtained from a weather
station installed on the instrument. The estimator provides products of
temperature, absolute humidity and backscatter ratio such that cross
dependencies between the lidar data products and raw observations are
accounted for and the final products are self consistent. The method described
here is applied to a combined oxygen DIAL, potassium HSRL, water vapor DIAL
system operating at two pairs of wavelengths (nominally centered at 770 and
828 nm). We perform regularized maximum likelihood estimation through the
Poisson Total Variation technique to suppress noise and improve the range of
the observations. A comparison to 119 radiosondes indicates that this new
processing method produces improved temperature retrievals, reducing total
errors to less than 2 K below 3 km altitude and extending the maximum altitude
of temperature retrievals to 5 km with less than 3 K error. The results of
this work definitively demonstrates the potential for measuring temperature
through the oxygen DIAL technique and furthermore that this can be
accomplished with low-power semiconductor-based lidar sensors.

##  1 Introduction

The MicroPulse DIAL (MPD) is designed to perform thermodynamic profiling of
the atmosphere using water vapor differential absorption lidar (DIAL), high
spectral resolution lidar (HSRL), and oxygen DIAL techniques. A testbed
network of five such MPD units has been developed at NSF NCAR. While water
vapor DIAL and HSRL have long legacies in lidar technology, oxygen DIAL is
relatively unexplored, with the first demonstrations of the approach occurring
on the MPD starting in 2020 [1]. The oxygen DIAL technique was definitively
described in [2] in the late 1990s, but the approach was dismissed as a viable
solution for temperature profiling due differences in oxygen absorption
between light scattered by aerosol and molecules. After scattering, aerosol
returns undergo the same narrow band absorption on the return trip, but
Doppler broadening from molecular scattering will result in a change in the
absorbing effect. This is commonly referred to as Rayleigh-Doppler effect in
DIAL and is typically a small error term that is either negligible, if the
absorption feature is sufficiently wide, or can be mitigated by operating at
an inflection point on the side of the absorption feature[3]. In the case of
oxygen DIAL, however, the effect of the error becomes significant in part due
to the narrow width of the oxygen lines under consideration (roughly the same
width as the broadening effect) and because it is necessary to operate near
the peak of the absorption feature to obtain sufficient sensitivity.
Uncorrected, this can result in temperature errors on the order of 10s of
degree C which [2] noted renders the technique relatively little value in the
field of atmospheric science.

In spite of these challenges, there are some important practical benefits to
performing temperature profiling using the oxygen DIAL technique. The DIAL
technique is well suited for low power instrument operation like that employed
for water vapor in the MPD. The technique is conducive to employing photon
counting with narrow band filters that enable daytime operation with eye-safe
(class 1M) semiconductor lasers. In addition, it has the potential to provide
observations which are internally calibrated. Such a capability, makes the MPD
one of the few instruments providing thermodynamic observations independent of
radiosondes. Sondes are often used as the gold standard thermodynamic
profiling, but it is important to recognize that they have biases that have to
be corrected and calibrated. Catching errors in sondes is extremely difficult
in part because they are used as the basis of calibration or priors for many
retrievals and ingested into models and reanalysis [4].

To enable temperature measurements through oxygen DIAL, the MPD team proposed
mitigating the Rayleigh-Doppler errors through an integrated HSRL, which
directly measures the composition of aerosol to molecular scattering in the
detection channels. This information then supplies the knowledge needed to
correct the oxygen spectroscopy and enable accurate temperature retrievals.
Thus when combined with existing HSRL and water vapor DIAL, (the oxygen DIAL’s
relationship to temperature is also dependent on water vapor concentration),
it should be possible to measure atmospheric temperature through the oxygen
DIAL technique.

While the combination of water vapor DIAL, HSRL and oxygen DIAL should carry
sufficient information to estimate atmospheric water vapor, backscatter ratio
and temperature, the means of estimation has remained a significant challenge.
The current method we employ relies on a serialized process, where water vapor
DIAL and HSRL data are initially processed using conventional techniques to
obtain absolute humidity and backscatter ratio estimates. Those products are
then inputs to a perturbative processing method [5, 6] to provide corrections
in the oxygen number density and absorption spectroscopy. While this approach
has shown promise, the following represent some key drawbacks:

  * •

Noisy data products are being used to provide corrections to the highly noise
sensitive oxygen DIAL retrieval, thereby compounding noise.

  * •

Water vapor and HSRL retrievals have cross dependencies and dependence on
atmospheric temperature, requiring assumptions to process.

  * •

Some MPD channels are under utilized because there is no seamless way of
integrating their overlapping information content into the standard processing
approach.

In this work we address these issues by performing a full estimate of absolute
humidity, temperature and backscatter ratio across all observation channels.
This is achieved by forward modeling the retrieved variables onto observations
of all six MPD channels and optimizing through a maximum likelihood estimation
framework while suppressing noise through total variation regularization. The
optimization approach leveraged in this work utilizes the Poisson Total
Variation (PTV) technique, which has been successfully applied to HSRL [7] and
MPD water vapor retrievals [8]. In this case, the temperature retrieval
requires significantly more retrieved products which has resulted in a need
for increased PTV processing speed and methods for reducing the problem size
so it can be solved with a manageable amount of time and resources.

This work describes the method we have developed for estimating range resolved
atmospheric water vapor, temperature and backscatter ratio using the MPD
instrument such that no external data sources regarding atmospheric state are
employed other than surface data obtained from a sensor installed on the
instrument. In this way, the retrievals are completely independent of all
other observations and models.

Because prior works have thoroughly addressed estimation of water vapor and
aerosol backscatter, the primary focus here is on the temperature retrievals,
which represent the most difficult retrieval and the technique with the
shortest legacy.

##  2 MicroPulse DIAL Architecture

The MPD is a semiconductor-based photon counting lidar designed around the
concept of quantitative atmospheric profiling in a package that accommodates
network deployment. The instrument is designed to be low cost, low
maintenance, portable and operate with an eye-safe (class 1M) laser. The
original MPD concept started with a water vapor only design at Montana State
University [9] which was advanced and demonstrated at NSF NCAR [10, 11, 12].
The instrument concept was later extended to include HSRL capability using a
rubidium vapor cell as an aerosol filter [13]. The HSRL filter was changed to
potassium, to allow for reduced hardware when further merging an oxygen DIAL
into the instrument for temperature measurements [1].

A schematic description of the MPD, useful for understanding the formulations
covered in this work, is shown in Fig. 1. As shown in the diagram, the MPD
operates at two main wavelength bands, each alternating between closely spaced
on and offline wavelengths while employing a total of three detectors. This
results in a total of six observations from the entire instrument.

The 828 nm wavelengths correspond to the water vapor (WV) DIAL system. As is
typical of DIAL, the laser output alternates between frequencies tuned on a
water vapor absorption line (online) and immediately off the line (offline). A
single detector is dedicated to the 828 wavelengths. The second set of
wavelengths are near 770 nm where temperature sensitive oxygen lines and the
potassium D1 transition are in close wavelength proximity to enable a combined
oxygen DIAL and HSRL. Like the water vapor DIAL, the laser frequencies are
alternated between a frequency tuned to the selected oxygen absorption line
and off the line. In this case the offline frequency is tuned to the potassium
D1 line to enable the HSRL functionality also integrated in the instrument. In
the receiver, the 770 nm light is split using a 70:30 beam splitter, where the
larger fraction is directed through the potassium cell to block aerosol
returns and provide a molecular only observation with the molecular detector
(offline wavelength only). The other 30 percent of the light is directed onto
the combined detector which operates similar to the water vapor detector.

While the oxygen offline frequency is used for the HSRL, the online frequency
produces nearly redundant measurements between the molecular and combined
detectors. These two online measurements are identical within a scalar
constant (due to differences in path efficiencies) and a differential overlap
term caused by slight differences in detector alignment and capture of the
return light. As such, these two channels carry very similar information
content with respect to the temperature profile but it is difficult to
leverage both channels effectively when estimating temperature through the
standard processing approach.

The reader should note that only the concept of the MPD architecture is
presented here. A number of components have been omitted or consolidated for
simplicity. For a complete description of the MPD see [12, 1, 10].

Figure 1: Conceptual diagram of the MPD operating at two main wavelengths
(each with alternating on- and offline laser frequency) and three detector
channels resulting in a total of six unique measurements.

##  3 Forward model

PTV employs a forward modeling approach to estimate derived products from raw
photon count observations. The derived products represent the estimates that
best explain the Poisson distributed, noisy observations in each channel.
Regularization is then used to prevent over fitting. In this way, cross
dependencies between channels can be accounted for and integrated into a
single retrieval which can suppress noise. The forward model employed here is
based on the lidar equation, which includes both time, range and frequency
dimensions to map the atmospheric derived variables absolute humidity,
backscatter ratio and temperature along with a series of instrument terms onto
the noisy observations.

Each MPD observation can be described through a general model where specific
terms in a channel may differ and are denoted by superscripts. These
differences may depend on laser wavelength (w) (e.g. o2on for oxygen DIAL
online wavelength), detector (d) (wv, comb or mol for water vapor, combined or
molecular detectors respectively – each detector in Figure 1) or both (c) such
that the term is channel specific. For this work, the general description for
the profile of a received signal as a function of range bin $k$ and relative
frequency $\nu$ where we account for the Rayleigh-Doppler effect in which the
narrow band laser light is broadened by molecular scattering [3]. In this work
$\nu$ is referenced to the laser transmit frequency and thus may assume
negative values. As such, the frequency resolved lidar equation separates the
outgoing and return narrow band absorption terms (i.e. water vapor and oxygen
resonances) where the absorption experienced by molecular backscattering is
different from that of aerosol backscattering. With the assumption that the
laser pulse is represented by a delta function in range and the laser spectrum
is a delta function in frequency, we employ the following forward model for
any MPD channel

$\displaystyle\hat{u}_{k}^{(c)}(\hat{\bm{n}},\hat{\bm{T}},\hat{\bm{B}},\hat{\bm{\phi}})=$
$\displaystyle\overbrace{\,G^{(w)}C_{k}^{(d)}\hat{\phi}_{k}e^{-\omega_{0,k}^{(w)}(\hat{\bm{n}},\hat{\bm{T}})}}^{\text{Scene
dependent terms and outgoing absorption}}$
$\displaystyle\times\biggl{[}\underbrace{\eta_{0}^{(c)}e^{-\omega_{0,k}^{(w)}(\hat{\bm{n}},\hat{\bm{T}})}\gamma^{(w)}\left(\hat{B}_{k}-1\right)}_{\text{Aerosol
scattering and return
absorption}}+\overbrace{\sum_{\nu}\eta_{\nu}^{(c)}e^{-\omega_{\nu,k}^{(w)}(\hat{\bm{n}},\hat{\bm{T}})}\tilde{\beta}_{\nu,k}^{(w)}(\hat{T}_{k})}^{\text{Molecular
scattering and return absorption}}\biggr{]}$ (1)

where $\hat{u}_{k}^{(c)}$ is the expected photon counts in range bin $k$ on
channel $(c)$ per laser shot as a function of estimated variables (where
estimated terms are indicated by $\,\hat{}$ ). The differential geometric
overlap function of the particular detector (relative to the combined
detector) is $C_{k}^{(d)}$, $G^{(w)}$ is a constant scalar accounting for
absorption below the first range bin and differences in laser power,
$\hat{\phi}_{k}$ represents the common terms across all the channels and
largely accounts for broadband atmospheric extinction (e.g. from aerosols and
clouds) and, for our parameterization, also captures the combined channel
geometric overlap function. The channel’s frequency dependent receiver
transmission is $\eta_{\nu}^{(c)}$ where the frequency index is indicated by
subscript $\nu$ and the transmission at the laser wavelength is
$\eta_{0}^{(c)}$ indicated by $\nu=0$, assuming a laser spectrum well
represented by a delta function. The frequency dependent optical depth
resulting from the estimated water vapor and oxygen absorption for the
specified wavelength is $\omega_{\nu,k}^{(c)}(\hat{\bm{n}},\hat{\bm{T}})$
(where when $\nu=0$ this is only evaluated at the laser center frequency).
This optical depth depends on the column of estimated water vapor number
density $\hat{\bm{n}}$ and temperature $\bm{\bm{T}}$, $\gamma^{(w)}$ is a
wavelength scaling term for when this channel operates at a different
wavelength than the HSRL system to account for differences in aerosol
scattering resulting in a different backscatter ratio from the estimate
$\hat{B}_{k}$ (which must be greater than 1). The frequency normalized
Rayleigh-Brillouin shape for the estimated temperature is
$\tilde{\beta}_{\nu,k}(\hat{T}_{k})$. Note that the summation term collapses
the laser frequency axis $\nu$ so that the forward model of the signal has
only range dimensions (where the time dimension is implicitly the result of
multiple observations).

The contributing terms in Eq. (3) may be physically understood as a series of
instrument and scene dependent scaling terms with an outgoing absorption term
(outside the brackets), an aerosol scattered term with return absorption
(first term in the brackets) and a molecular scattered term with return
absorption (second term in the brackets).

In this work retrievals are performed on 2D images of time and range to enable
PTV to capture correlated structure across both dimensions. However the
formulations presented here will omit the time component for simplicity. Terms
that are range dependent are represented as vectors in bold or are
individually indexed as in Eq. (3) in range $k$ or frequency $\nu$. Variables
without indices or bold are scalars.

Under the above instrument model, we assume that the laser wavelength can be
accurately represented as a delta function in frequency which creates the need
to ensure that the instrument operates within the constraints that this
description is valid.

The backscatter signal from Eq (3) is convolved in range with the laser pulse
shape $\bm{l}^{(w)}$, multiplied by the number of laser shots in the
accumulation interval $m^{(c)}$ and the background photon counts are added to
obtain the full estimate of the observed photon counts. In addition we include
an afterpulse baseline $\bm{a}^{(c)}$ induced by the outgoing laser pulse in
the transceiver and obtained through calibration.

$\hat{\bm{y}}^{(c)}=m^{(c)}\left(\bm{l}^{(w)}*\hat{\bm{u}}^{(c)}+\bm{a}^{(c)}\right)+b^{(c)}$
(2)

The optical depth $\omega$ in Eq. (3) depends on both oxygen and water vapor
absorption and is calculated using trapezoidal integration so

$\omega_{\nu,k}(\bm{n},\bm{T})=\omega_{\nu,k-1}(\bm{n},\bm{T})+\frac{1}{2}\left[\kappa_{k-1}(n_{k-1},T_{k-1})+\kappa_{k}(n_{k},T_{k})\right]\Delta
r_{k}$ (3)

where $\Delta r_{k}$ is the range spacing between points and $\kappa_{k}$ is
the combined extinction coefficient from both species

$\displaystyle\kappa_{k}(n_{k},T_{k})$
$\displaystyle=o_{\nu,k}(n_{k},T_{k})+q_{\nu,k}(n_{k},T_{k})$ (4)
$\displaystyle=f_{o2}\left(\frac{P_{k}}{k_{B}T_{k}}-n_{k}\right)\sigma_{\nu,k}(T_{k})+n_{k}\tau_{\nu,k}(T_{k})$
(5)

where $o_{\nu,k}(n_{k},T_{k})$ and $q_{\nu,k}(n_{k},T_{k})$ are the absorption
coefficients of oxygen and water vapor respectively, $f_{o2}$ is the molar
fraction of oxygen in dry air (0.211032), $P_{k}$ is pressure, $k_{B}$ is the
Boltzmann constant, $\sigma_{\nu,k}(T_{k})$ is the temperature dependent
absorption cross section of oxygen, $\tau_{\nu,k}(T_{k})$ is the absorption
cross section of water vapor (with weak temperature dependence). The pressure
is determined through a hydrostatic relationship which is dependent on the
surface pressure (measured by a weather station on the MPD unit) and the mean
(linear) temperature profile [14].

$P_{k}=P_{0}\left(\dfrac{\bar{T}_{k}}{\bar{T}_{0}}\right)^{C_{g}},$ (6)

where

$C_{g}=\frac{g_{0}M_{air}}{R_{0}L_{0}},$ (7)

where $g_{0}$ is the acceleration due to gravity (9.81 m/s), $M_{air}$ is the
molar mass of air (28.97 g/mol), $R_{0}$ is the universal gas constant
(8.314459 J/mol/K) and $L_{0}$ is the linear fit temperature lapse rate in
K/km. To obtain the this lapse rate, the mean temperature profile
$\bar{\bm{T}}$ is parameterized as a linear fit to the PTV estimated
temperature $\hat{\bm{T}}$ using the lapse rate $L_{0}$ and a zero range
intercept $\bar{T}_{0}$ such that

$\bar{T}_{k}=r_{k}L_{0}+\bar{T}_{0}$ (8)

where $r_{k}$ is the range at bin $k$. The lapse rate and zero intercept are
obtained from a linear fit to the temperature profile estimate

$\left[\begin{array}[]{c}\bar{T}_{0}\\\
L_{0}\end{array}\right]=\left[\begin{array}[]{cc}\bm{1}&\bm{r}\end{array}\right]^{-P}\hat{\bm{T}}.$
(9)

where the superscript -P indicates a pseudo-inverse operation. Both of these
linear fit terms are retrieved at the same temporal resolution as the
temperature profile estimate.

###  3.1 Absorption and Rayleigh-Brillouin spectrum

Frequency resolved oxygen and water vapor absorption cross sections
($\sigma_{\nu,k}$ and $\tau_{\nu,k}$ respectively) are key parameters for
performing thermodynamic estimation using the DIAL technique. In particular,
$\sigma_{\nu,k}$ provides the dominant sensitivity to temperature in the
oxygen DIAL technique, making it imperative that this quantity accurately and
rapidly update during optimization. It also must be differentiable with
respect to temperature. In this work we employ a principal component based
method for calculating the oxygen absorption, water vapor absorption and
Rayleigh-Brillouin spectra [15, 16]. In this technique, a series of principal
components are ”trained” on first principle calculated spectra across a range
of temperature and pressure combinations and the weights of these components
are fit to polynomials that are a function of temperature and pressure. This
approach is significantly faster than performing the calculations based on
first principles (which are prohibitively slow for an optimization routine
like PTV) and easily generalized to arbitrary spectra. Unlike lookup table
approaches, it is also analytically differentiable which is also a key feature
when speed is needed in numerical optimization.

This principal component approach to spectroscopy calculations runs in
contrast to the parameterized inversion method proposed in other papers on
temperature estimation with oxygen DIAL [2, 5, 6]. In these approaches, the
following approximations are made: 1) Rayleigh-Brillouin scattering is
approximated as Gaussian functions, 2) the laser wavelength is assumed to be
centered on the absorption line (which shifts in pressure/altitude so is not
strictly possible over a large altitude range) and 3) the absorption is
unaffected by additional interfering lines (there are a number of possible
interfering lines in the oxygen A-band due primarily to overtone absorption
and less abundant isotopes of oxygen).

###  3.2 Estimate parameterization

In this work, the variables we estimate using PTV represent parameterizations
of the physical variables described in Eq. (3). These parameterizations are
used to impose soft limits on the physical variable range, impose some amount
of normalization to the retrieved variable domain or better match the
retrieval basis to the expected form of the physical variable. The notation
employed here is that the total estimated state is described by a single
variable consisting of several specific state variables
$\bm{X}=\left\\{\bm{x}_{\phi},\bm{x}_{T},...\right\\}$. In the notation of
this work, these estimated state variables have a direct relationship to the
physical variables in Eq. (3) denoted by their subscript. For example the
common terms $\bm{\phi}$ are obtained from the estimated variable
$\bm{x}_{\phi}$ through a parameterizing function (defined in the third column
of Table 1) $\bm{\phi}=f_{\phi}(\bm{x}_{\phi})$. Table 1 also includes the
name used to refer to the terms and the dimensions (time $t$ and/or range $r$)
in which the variable is represented. The formulations are such that
constraints do not have to be set to avoid non-physical values in most
channels. Thus the range of variables are bounded only by the parameterization
function.

Table 1: Estimated Variables Name | Symbol | Parameterization | Dims
---|---|---|---
Common Terms | $\bm{\phi}$ | $\phi_{k}=\exp\left(x_{\phi,k}\right)$ | $t,r$
Backscatter Ratio | $\bm{B}$ | $B_{k}=1+\exp\left(x_{B,k}\right)$ | $t,r$
WV Differential Overlap | $\bm{C}^{(wv)}$ | $C^{(wv)}_{k}=\exp\left(x_{C^{(wv)},k}\right)$ | $t,r$
WV Number Density | $\bm{n}$ | $n_{k}=x_{n,k}N_{A}/M_{H_{2}O}$ | $t,r$
Mol. Differential Overlap | $\bm{C}^{(mol)}$ | $C^{(mol)}_{k}=\exp\left(x_{C^{(mol)},k-1}+x_{C^{(mol)},k}\right)$ | $r$
Temperature | $\bm{T}$ | $T_{k}=T_{k-1}+\Delta r_{k}x_{T,k}$ | $t,r$
Online Gain | $G^{(o2on)}$ | $G^{(o2on)}=\exp\left(x_{G^{(o2on)}}\right)$ | $t$

While most of these parameterizations are straightforward we should note that
the Molecular Differential Overlap is a cumulative sum operation in range
where the initial value $x_{C^{(mol)},0}=0$. Also, the retrieval of
temperature is actually obtained by estimating a localized lapse rate for each
point. This is because temperature is typically better represented as a sparse
set of piecewise linear functions in range than piecewise constant functions
typically obtained through PTV processing. The variable $\Delta r_{k}$ is the
spacing in range between each retrieved point. The surface temperature profile
$T_{0}$ (not to be confused with $\bar{T}_{0}$ which is the zero intercept of
the linear fit to the temperature profile) is obtained directly from the
surface station measurement on the MPD.

In the definition for WV Number Density, $M_{H_{2}O}$ is the molar mass of
water (18.018 g/mol) and $N_{A}$ is Avogadro’s number. The forward model
requires number density to calculate the lidar signals, but the retrieved
variable is absolute humidity, with a domain better conditioned for
optimization because typical values are not large floating point numbers.

Implicit in the parameterization step is also a resampling function such that
any estimated state variable may be up-sampled to the resolution of the
observations. That is, a state variable may be estimated at lower resolution
than the observations. In this work, this is only employed for the temperature
($\bm{T}$) retrieval. This effectively constrains the resolution of the
parameter which tends to be slowly varying while the retrieval is highly noise
sensitive.

###  3.3 Calibrations

In order to evaluate the forward model for each channel, the MPD requires a
series of internal calibrations to obtain accurate estimates of terms in Eq.
(3). It’s important to recognize the distinction between internally and
externally calibrated instrumentation. Under internal calibration, we aim to
directly measure the terms needed in the instrument model. As such, we employ
a series of hardware, software and signal processing solutions that enable
direct measurement of those specific terms.

MPD data products and calibrations are never modified to obtain agreement with
an external reference. Sonde data is highly useful in the development of the
instrument but only used for validation.

External calibration methods often similarly aim to characterize instrument
terms in the physical model, but are obtained by forcing agreement between the
instrument retrieved product and an external reference (typically a sonde). In
practice those physical terms represent ”catch-all” corrections for the
instrument data product while appearing to be physically justified. That is,
if the physical term were directly measured, it would not necessarily correct
the retrieval by the same amount or in the same way. Instrument models can be
oversimplified or contain assumptions that are not fully justified which
inherently represent ”unknown unknowns”. External calibration allows one to
avoid these complexities, enabling more immediate pursuit of scientific goals
at the cost of rigorous understanding of the instrument, inheriting errors
from the reference and the risk of inadvertently imposing biases due to
incomplete sampling of factors contributing to the calibration terms.

The principle calibrations performed on MPD consist of a receiver scan to
characterize the frequency dependent transmission of each channel
$\eta_{\nu}^{(c)}$ and an afterpulse calibration to obtain the baseline signal
induced noise from the outgoing laser pulse $\bm{a}^{(c)}$.

The receiver scan is performed by sweeping the wavelength of a seed laser over
the passband for each channel, as described in [12]. For most channels, this
captures the etalon transmission as a function of frequency while accounting
for angular spread in the light passing through the receiver. A physical model
of the etalon, and potassium absorption (in the case of the oxygen offline
molecular channel), is used to obtain a denoised estimate of the receiver’s
characteristics. In addition to the receiver spectrum shape, this also
provides the relative scaling efficiency of each detector which largely
mitigates the need for any additional scaling of the receiver channels
($G^{(w)}$, the channel gain accounting for differences in collected power, is
strictly laser dependent and even then, 1 for most instances).

The afterpulse baselines are obtained by covering the instrument ports to the
atmosphere so that any returns can be attributed to afterpulsing effects from
the outgoing pulse. These profiles are accumulated over approximately a half
hour and processed using a physical model based on [17, 18] where the
afterpulse is a weighted sum of decaying exponentials. Thus the forward model
used to fit to noisy afterpulse data is

$\displaystyle\rho_{k}=$ $\displaystyle b+a_{k}$ (10) $\displaystyle\rho_{k}=$
$\displaystyle b+\sum_{i}c_{i}\exp\left(-\frac{r_{k}}{p_{i}}\right)$ (11)

where $\rho_{k}$ is the expected photon counts per shot at the $kth$ range
bin, $b$ is a constant background term (not used in actual lidar processing),
$a_{k}$ is the afterpulse function employed in Eq. (2), $r_{k}$ is the range,
$c_{i}$ is the exponential coefficient that is used for fitting to the noisy
data and $p_{i}$ is the decay constant in range of the afterpulse mode.

In addition to directly calibrated quantities, we can obtain the differential
overlap between the oxygen molecular and combined channels $C^{(mol)}$ through
the retrieval itself. This is possible because overlap is the only major
distinction between the two channels when the laser is transmitting the online
wavelength. To simplify and accelerate the estimation process, we only
estimate this quantity as needed (e.g. when other calibrations are updated)
and it is generally stored as a calibration parameter.

##  4 Optimization

We employ the Poisson Total Variation (PTV) technique in order to obtain
denoised estimates of atmospheric and instrument variables. This technique is
able to efficiently perform total variation regularized maximum likelihood
estimation. The technique described in the seminal work [7], was applied to
water vapor retrievals in [8] and its implementation is described in further
detail in for high resolution estimation in [19]. Here we provide a high-level
summary of the technique.

We aim to obtain some set of estimated state variables $\bm{X}$ (consisting of
all terms summarized in Table 1) which may consist of both time and range
dimensions. The state variables are related to noisy photon count observations
$\bm{y}^{(c)}$ (consisting of both time and range dimensions) through a
forward model $\bm{\hat{y}}^{(c)}(\bm{X})$ as described in Eq. (2).

Our goal is to obtain an estimate of $\bm{X}$ which minimizes the objective
function

$\hat{\bm{X}}=\operatorname*{argmin}\mathcal{O}(\bm{X})$ (12)

where the objective function consists of a negative log-likelihood term
$\mathcal{L}$ and a total variation (TV) regularization term with scalar
$\alpha_{i}$ for each estimated variable $\bm{x}_{i}$ in $\bm{X}$ such that

$\mathcal{O}(\bm{X})=\mathcal{L}(\bm{X})+\sum_{i}\alpha_{i}||\bm{x}_{i}||_{TV}.$
(13)

The negative log-likelihood of Poisson distributed observations is used to
evaluate the fit between the forward model and the noisy observations

$\mathcal{L}(\bm{X})=\sum_{c}\bm{w}^{(c)}\left[\bm{\hat{y}}^{(c)}(\bm{X})-\bm{y}^{(c)}\ln\bm{\hat{y}}^{(c)}(\bm{X})\right],$
(14)

where $\bm{w}^{(c)}$ is the channel weight which includes a mask for non-
Poisson data and is described in further detail in the next section. Note that
in the validation steps described later in the regularizer search process,
this same equation is used to evaluate the estimate where $\bm{y}^{(c)}$ is
replaced with independent validation data.

In the estimation process, the objective function is minimized through the
Spiral-TV [20] framework with gradients calculated using fast iterative
shrinkage/thresholding algorithm (FISTA) [21].

The PTV software used in previous works was based on a Cython implementation
designed for CPU hardware. We have developed an independent PTV package for
this new application. This new package is built on the PyTorch machine
learning package [22] and is designed to leverage GPU hardware to accelerate
the gradient calculations and enable faster prototyping through the package’s
”autograd” feature. In addition to moving most of the data and model onto a
GPU, a CUDA kernel was developed to further accelerate the implementation of
FISTA [21] used for gradient calculations. This effectively removed FISTA as a
bottleneck in the processing and enabled retrievals in realistic time frames.
With this new PTV processing capability, fully processing MPD data consisting
of a 24 hour period typically takes between 3 and 5 hours (depending on
cluster computing queue wait times) and consumes approximately 50 GPU hours on
the NSF Derecho cluster computer [23]. Processing speed typically depends on
the size of data being processed, so shorter time periods and lower
resolutions are expected to complete in shorter amounts of time and fewer GPU
hours.

###  4.1 Masks and weights

Similar to the methods described in [8], masking must be applied to photon
count data that is poorly approximated as Poisson distributed. The primary
example of this is observations of clouds, where DIAL data is notorious for
producing non-physical data products. We suspect that these regions cannot be
corrected using standard deadtime corrections because the photon flux is both
high and rapidly varying. As a result, the captured photon counts tend to have
high error that is poorly captured in the Poisson noise model and can corrupt
the retrieval of the full scene.

In previous works, we employed a binary masking scheme, however in this work
we apply a continuous weighting function based on a logistic function

$w_{k}^{\prime(c)}(z_{k})=\frac{1}{1+\exp(z_{k})}$ (15)

where the argument of the function $z_{k}$ is determined based on the maximum
observed count rate (determined at a capture resolution of 2 seconds) relative
to a configuration set maximum count rate $\rho_{max}$ (2 MHz for this work)
which is scaled by the observed standard deviation in photon counts divided by
that expected for a Poisson random number (square root of the mean counts+1).
Each processed time bin is an accumulation of two second raw time bins $\Delta
t$ (the typical capture resolution of the instrument) which are used to
directly calculate the max, mean ($\mu_{y_{k}}$) and standard deviation
($\sigma_{y_{k}}$) of the photon counts accumulated in $y_{k}$ and the
argument to the logistic function in Eq. (15) is

$z_{k}=\frac{\max(y_{k})/\Delta
t}{\rho_{max}}\frac{\sigma_{y_{k}}}{\sqrt{\mu_{y_{k}}+1}}.$ (16)

Through this masking scheme, data points with high backscatter and high
variability (in excess of what is expected from a Poisson distribution) have
reduced weight on the optimization loss function and therefore have lower
priority in the optimization process.

In addition to weighting each bin, we also include a scalar weighting for the
entire channel based on the available photon counts for fitting. This was also
done in [8]. This scalar channel weight is given by

$w^{(c)}=\frac{1}{\sqrt{\sum_{k}\left(w^{\prime(c)}_{k}y^{(c)}_{k}\right)^{2}}}$
(17)

where the scalar weight is dependent on the non-scalar weight
$\bm{w}^{\prime(c)}$ in Eq. (15). The total weight in Eq. (14) is then

$\bm{w}^{(c)}=w^{(c)}\bm{w}^{\prime(c)}$ (18)

###  4.2 Regularizer search space

PTV is built on a linear piecewise basis set, meaning the retrieved parameters
are assumed to be well represented by an image composed of patches of constant
value. In the noisy observations, changes from pixel to pixel for each
estimated parameter can either caused by the actual parameter fluctuations or
noise. Regularization (forcing the retrieval to prefer simpler solutions)
allows the retrieval to balance adding complexity to the retrieved parameter
while not over fitting noise. As stated in the objective function definition
in Eq. (13), each estimated variable $\bm{x}_{i}$ in $\bm{X}$ (described for
this specific case in Section 3.2) has an associated regularizer $\alpha_{i}$
which determines the amount of total variation penalty.

The standard practice in PTV is to optimize each variable regularizer
$\alpha_{i}$, by fully mapping the associated validation error as a function
of regularizer values. This is done by solving PTV at each regularizer (or
combination of regularizer) value and projecting the forward model solution
onto a statistically independent set of observations of the same scene using
the loss function in Eq. (14). Obtaining this independent validation data is
possible because Poisson thinning [24, 25, 26] enables us to split Poisson
distributed observations of one scene into two independent and identically
distributed (i.i.d.) sets of observations. As such, fitting to noise is
penalized when evaluating the result against the independent validation data
(using Poisson thinning to optimize tuning parameters in lidar processing is
covered in more detail in [26, 25, 19]).

Most previous PTV instances optimize the regularization parameters by
thoroughly mapping out the regularizer space across all retrieved variables.
However, in the case described here, there are six independent variables.
Fully evaluating a six dimensional regularizer space is not practical with our
current software speed and computing capabilities (the computational expense
increases by the power of the search space). In order to address this we have
taken a more heuristic approach that allows us to break the search process
into a series of smaller problems. This is accomplished by performing
estimation in stages to obtain accurate initial conditions (which accelerate
processing) and significantly localize the regularizer search space in the
final global estimation step. At each stage, the solutions obtained from
previous steps are employed either as initial conditions (if the variable is
estimated at the new step) or constants in the forward model. The breakdown of
the PTV processing steps are as follows:

  1. 1.

Estimate the Molecular Differential Overlap $\bm{C}^{o2_{mol}}$ and Common
Terms $\bm{\phi}$ by fitting to the Oxygen Molecular Online and Oxygen
Combined Online channel observations.

  2. 2.

Estimate Common Terms $\bm{\phi}$ by fitting to only the Oxygen Offline
Molecular channel observations.

  3. 3.

Estimate the Common Terms $\bm{\phi}$ and Backscatter Ratio $\bm{B}$ by
fitting to the Oxygen Offline Molecular and Combined channel observations.

  4. 4.

Estimate the Water Vapor Differential Overlap $\bm{C}^{wv}$ by fitting to only
the Water Vapor Offline channel observations.

  5. 5.

Estimate the Water Vapor Differential Overlap $\bm{C}^{wv}$ and Water Vapor
Number Density $\bm{n}$ by fitting to only the Water Vapor Offline and Water
Vapor Online channel observations.

  6. 6.

Estimate the scalar Online Gain $G^{(o2_{on})}$ by applying an exponentially
decaying fit weight (in range) and fitting to the Oxygen Online Combined
channel observations.

  7. 7.

Perform global estimation of all state variables fitting to all channel
observations but freeze all regularizers except temperature to their
previously determined solutions.

  8. 8.

Perform global estimation of all state variables fitting to all channel
observations where all regularizers are allowed to randomly vary around the
minimum obtained in the previous step.

Note that in all steps before (7.) no more than two variables are ever
estimated. This allows us to perform a full regularizer search and obtain the
optimal regularizers for those preprocess steps.

Step (1.) represents an exception case for a few reasons. First, the common
terms estimate is used purely as a proxy variable for everything other than
the differential overlap term and so that specific solution is not used in any
subsequent steps. In addition this step is not run every time because it adds
significant time to the processing (it turns out to be one of the slowest
optimization steps). Since the key term obtained from this step is an
instrument parameter that should be stable over relatively long time periods,
the step is only occasionally run and the result is stored as an instrument
calibration. Most processing runs use this calibration instead of running step
(1.).

The regularizer search method applied to each of these steps also differs from
the typical PTV approach. Here, in a processing run, the regularizer value is
randomly generated from a Gaussian random variable (with a mean updated based
on each search step). This approach is conducive to employing a Gaussian
process search (often used for hyperparameter optimization in machine
learning), however we have found that centering the search on the regularizer
point with the lowest validation loss is typically more robust (but slower to
converge) because it is less affected by bad optimization runs (e.g.
optimization fails due to significantly excessive regularization). For the
purpose of routinely processing data, this robustness is more appealing than
any potential time savings (which are lost as soon as the process fails).

The initial conditions at the start of this optimization routine are
summarized in Table 2. Most terms are initialized as constant values. The
temperature profile is initialized as a linear lapse from the surface station
temperature. Also, the water vapor differential overlap term is typically
initialized based on the relative scale of signals between water vapor and
oxygen channels.

Table 2: Initial Conditions of Variables Name | Initial Condition
---|---
Common Terms | 1.0
Backscatter Ratio | 1.01
WV Differential Overlap | system dependent constant
WV Number Density | 0
Mol. Differential Overlap | 1.0 or prior calibration
Temperature | 9 K/km lapse from surface
Online Gain | 1.0

##  5 Uncertainty Estimation

Uncertainty estimation represents an important element for observational data
products. For the performance assessment presented here, we require the metric
as a means of providing quality control (QC) to the data products. Data known
to be invalid based on the uncertainty obtained from the instrument
observations is not included in the subsequent performance analysis.

###  5.1 Bootstrapping

In the standard MPD processing routines, we leverage a combination of Poisson
thinning and bootstrapping to obtain estimates of derived parameter
uncertainties [12]. Under this method, we split each channel’s photon counts
into two i.i.d. scenes, process each, and calculate the square difference of
the data products. This process is repeated several times and averaged to
obtain a numerical estimate of the error due to shot noise. We have found this
method consistently captures the derived product errors better than linear
propagation of error.

The challenge in implementing bootstrapping with this PTV processing routine
is that repeating the calculations many times becomes prohibitively expensive
in processing time. To fully encompass the impact of uncertainty, the full
processing chain must be repeated along with variations in initial conditions
at each bootstrapping step. Because of the number of parameters and
sensitivities, it also implies a rigorous approach would likely require more
bootstrapping iterations than are typical of our water vapor and HSRL
retrievals. To provide some amount of uncertainty information, we have
implemented a “bootstrap-lite” approach in which the final global processing
step is repeated 12 times with independently thinned photon counts and random
Gaussian scalars are added to the initial conditions of each state variable.
The resulting standard deviation of each data product is stored from this step
to help establish errors that are essential for data QC in the data products.
While this method works relatively well with water vapor products, it is very
difficult to interpret in the temperature product, sometimes reporting low
errors in regions we know have high uncertainty (i.e. masked regions and
regions blocked by clouds).

###  5.2 Evidential Neural Net

The bootstrap-lite approach provides some information about the potential
variances of the estimated variables. However, this result is difficult to
interpret, particularly for the final temperature estimate which is obtained
by range integrating the estimated lapse rate. In order to help with the
uncertainty estimation process, we turned to using an evidential machine
learning technique originally proposed for performing inference on machine
learning output while providing uncertainty [27, 28].

The formulation used in [27] leverages a parameterization of the epistemic
uncertainty (uncertainty in the mean) which allows the marginal distribution
to analytically evaluate. In their approach, the epistemic uncertainty is a
scaled version of the aleatoric uncertainty (the variance and its associated
uncertainty) where the scale term is predicted by the neural net. However in
this work there is no profile-by-profile capability to measure a mean error
(effectively a bias) based on instrument observed parameters. This is only
possible as an ensemble analysis when comparing to a reference over many
profiles. As a result, this parameterization of epistemic uncertainty is
poorly suited to our application.

In this work we modify the approach to explain all uncertainties in the
temperature using the epistimic (Gaussian variance) term alone. This avoids
asserting an estimate of the mean error which cannot be predicted from the MPD
data on a profile-by-profile basis and would require a composite analysis (and
ultimately be an external calibration). The result of this approach is that
the error data is not truly Gaussian distributed because a persistent 1.5 K
bias (seen later in Section 6) is not removed in the uncertainty estimator.
These results are nevertheless quite useful in performing data QC on the
temperature estimates.

Where the original evidential technique estimates a mean and variance
conditioned on a Normal Inverse-Gamma prior, our application drops the mean
estimation component. We aim to estimate the variance of the temperature
output by assuming that the error is Normally distributed with zero mean

$f(\Delta;\sigma^{2})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{\Delta^{2}}{2\sigma^{2}}\right)$
(19)

where $\Delta$ is the difference between the MPD temperature retrieval and a
temperature reference (ECMWF ERA 5 reanalysis in this case [29]) and
$\sigma^{2}$ is the distribution’s variance.

In accordance with [27] we then impose an inverse-Gamma prior on the variance
estimate

$g(\sigma^{2};\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\left(\frac{1}{\sigma^{2}}\right)e^{-\beta/\sigma^{2}}$
(20)

Applying Bayesian probability theory, we can obtain the PDF of the temperature
error conditioned by the inverse-Gamma prior

$\displaystyle p(\Delta;\alpha,\beta)$
$\displaystyle=\int_{0}^{\infty}f(\Delta;\sigma^{2})g(\sigma^{2};\alpha,\beta)d\sigma^{2}$
(21)
$\displaystyle=\frac{1}{2\pi}\frac{\beta^{\alpha}\Gamma\left(\alpha+\frac{1}{2}\right)}{\Gamma(\alpha)}\left(\frac{2}{\Delta^{2}+2\beta}\right)^{\alpha+1/2}$
(22)

With this PDF analytically defined, the resulting loss function for the neural
net is defined as the negative log-likelihood of the PDF

$\mathcal{L}^{E}(\alpha,\beta)=\frac{1}{2}\ln
2\pi-\alpha\ln\beta+\left(\alpha+\frac{1}{2}\right)\ln\left(\frac{\Delta^{2}+2\beta}{2}\right)+\ln\left(\frac{\Gamma(\alpha)}{\Gamma(\alpha+1/2)}\right).$
(23)

In addition to the negative log-likelihood, [27] employs a regularization term
which we adapt to the zero mean application as

$\mathcal{L}^{R}(\alpha)=|\Delta|(2+\alpha)$ (24)

so that the total loss function in the training process is

$\mathcal{L}(\alpha,\beta)=\mathcal{L}^{E}(\alpha,\beta)+\lambda\mathcal{L}^{R}(\alpha).$
(25)

By leveraging this loss function, the neural network is trained with inputs of
vertical temperature profile, the standard deviation in temperature across
bootstrap-lite estimates and the standard deviation of lapse rate across
bootstrap-lite estimates to predict values for $\alpha$ and $\beta$ which can
statistically describe the error in the temperature measurements. Once
trained, the neural network provides $\alpha$ and $\beta$ for a given profile
and the variance can be calculated directly from these parameters

$\sigma^{2}_{e}=E\left[\Delta^{2}\right]=\frac{\beta\Gamma(\alpha-1)}{\Gamma(\alpha)}.$
(26)

We note that this uncertainty estimate is somewhat different from that
provided in [27]. Here we obtain the variance by calculating the expected
value of the second moment of $\Delta$ using $p(\Delta;\alpha,\beta)$ which
accounts for spread in the prior. In [27], they computed the expected value of
$\sigma^{2}$ directly from the Inverse Gamma PDF $g(\sigma^{2};\alpha,\beta)$.
In most instances we find the differences are small, however in some cases,
accounting for the spread imposed by the Inverse Gamma function will result in
larger error estimates.

The neural net used to perform the inference consists of two 512 dense layers
with ReLu activation functions followed by a linear output layer. The neural
net was built and trained using PyTorch. The inputs to the neural net consist
of the vertical temperature profile, the standard deviation in temperature
across bootstrap-lite estimates and the standard deviation of lapse rate
across bootstrap-lite estimates. The temperature error ($\Delta$) is
calculated by comparing the MPD estimate to ECMWF ERA 5 temperature [29, 30].
Training, validation and test data for the process were grouped by days to
minimize temporal correlation between training, validation and test datasets
which require statistical independence to fulfill their respective purposes.
The selection of days for each of these sets however, was random.

While we would argue reanalysis data is not a reasonable truth standard for
state variable estimation, we claim it makes a reasonable standard when
performing uncertainty estimates, where the accuracy requirements are more
relaxed. For the purposes of this work, the uncertainty is strictly used as a
quality control tool to determine where data should be masked. As a practical
consideration, the use of reanalysis enables us to obtain a much larger
training dataset than would be possible with the relatively few sondes. This
large dataset is essential for training neural nets.

An analysis of the uncertainty estimation was performed by comparing the
predicted uncertainties to the difference between MPD temperature estimates
and 119 sonde measurements from the M2HATS field project [31]. Note that none
of these sondes were used in the training, validation, or test process of the
neural net. The histogrammed results of that analysis for the bootstrap-lite
process (described in the previous section) and the evidential neural net
(which uses the bootstrap-lite process as an input) are shown in Fig. 2.
Notably the bootstrap-lite output alone has poor correlation to the MPD/sonde
differences. The evidential neural net results, while far from perfect,
represent a considerable improvement in interpreting bootstrap results with a
clear correlation to MPD/sonde differences. The performance of the evidential
neural net is sufficient to perform data QC analysis critical for evaluating
the performance of the retrieval method in the next section.

Figure 2: Comparison of bootstrap standard deviation (left) and evidential
neural net uncertainty (right) to temperature difference between PTV processed
MPD data and sondes. The dashed line is the 1:1 line. Note that the bootstrap
estimates are used as inputs to the neural net. The color scales are
logarithmic.

We would note that while this method works adequately well for temperature
uncertainty estimates, using reanalysis for water vapor error estimates is not
likely to be as productive because reanalysis water vapor products tend to be
much less accurate than the temperature products. For that purpose we are
forced to rely directly on the bootstrap-lite results which, fortunately, are
a more direct predictor of uncertainty in the water vapor product.

##  6 Results

NSF NCAR has recently developed a testbed network of five MPD units that all
measure temperature, absolute humidity and calibrated aerosol backscatter.
From July 23 to Sept 24, 2023, MPD 3 was deployed to Tonopah, NV, USA as part
of the M2HATS [32] field campaign. During that time, radiosondes [31] were
launched at the MPD site twice a day typically near 1700 and 2200 UTC. Data
from the unit was processed for the full time period to obtain temperature,
absolute humidity and backscatter ratio. We applied both the standard and PTV
processing approaches to the collected data and used the radiosonde products
to evaluate the two approaches.

The PTV data processing employed the global thermodynamic estimation approach
described in this work. It was conducted at a base resolution of five minute
time bin intervals with 37.5 m range bins. All products except temperature
were obtained at these resolutions. Temperature, due to the high noise
sensitivity, was estimated at a resolution of 40 minutes and 150 m.

The standard thermodynamic processing utilized standard water vapor DIAL [12,
10, 33, 2] and HSRL [13, 34] processing to obtain absolute humidity and
backscatter ratio estimates. The results of those processing outputs were then
supplied to the perturbative processing method described in [5, 6] which
treats these products as constants in the temperature estimation process.
Temperature retrievals were processed at a base resolution of 1 minute and
37.5 m then smoothed in each perturbative step to 20 minutes and 300 m.
Finally a 40 minute median filter was applied to the data. Uncertainty was
assessed using the bootstrapping method, but importantly only applied to the
temperature processing step, so uncertainties in the retrieved backscatter
ratio and absolute humidity that compounded in the temperature estimate were
not accounted for.

An example of temperature estimates from both methods is shown for a three day
period spanning August 28 through August 30 in Fig. 3. The data is masked
based on uncertainty analysis employed by the two techniques (bootstrapping
for the standard method and evidential NN for PTV) and times of radiosonde
launches are shown with vertical dashed lines. This time period was selected
because it contains a radiosonde launched near 0440 UTC, the only sonde
launched at night. Comparisons of the MPD temperature retrievals to those
sondes are shown below the temperature curtain plots where the shaded regions
indicate uncertainty for the respective methods determined in accordance with
the aforementioned techniques. In addition, those sonde comparisons also show
the initial condition used in the PTV temperature retrieval as black dotted
line (obtained by assuming a constant lapse rate of 9 K/km from the surface
station). This allows us to see that the retrieval is deviating from the
initial condition and not just a product of a lucky initial guess.

Note that in this example, the PTV retrieval of temperature appears to be
responsive to inversion layers. This is important because temperature
retrievals often lack features that significantly distinguish them from an
assumed linear lapse rate (such as the initial condition in the optimization
routine). While the retrieval lacks sufficient resolution to track sharp
changes such as that seen in 2023-08-29 22:33:38 near 3 km, it does follow the
overall trend of the change. The PTV retrieval is less responsive to noise
than the standard method, generally making these changes easier to recognize
and enabling the retrieval to extend high enough to capture them.

Figure 3: Curtain plots of temperature estimates from August 28-30, 2023 using
PTV processing (top) and standard processing (middle) with sonde times
indicated by dashed lines. The bottom row shows the temperature comparisons of
the two methods (blue: PTV, orange: standard method) with the sondes (heavy
gray) as well as the initial conditions (dotted black) used for the PTV
retrieval.

We conducted an analysis of both temperature retrieval methods based on their
agreement with 119 radiosondes from the M2HATS project [31]. The results of
that analysis are summarized in Fig. 4. For each processing approach we only
considered data that was not flagged due to high uncertainty. In order to
penalize overzealous masking, we report the data availability as a function of
altitude for all sonde comparisons. We calculate RMS difference with sonde
temperature to provide a full analysis of the retrieval errors, mean
difference to estimate bias in the retrievals and standard deviation to
estimate spread in the retrieval errors. We should note that all sondes,
except one, were launched during the day, which represents the highest noise
conditions for the lidar.

We can see that PTV tends to have slightly lower data availability in regions
where the standard method is consistently valid. This is consistent with
results we reported in [8] and is the result of masking due to high photon
rates (particularly near clouds). For PTV to perform accurate retrievals it
requires that the Poisson noise model accurately describes the photon
detection process. However in regions of high photon flux, this noise model is
no longer valid. In order to avoid such cases corrupting the retrievals (as it
impacts the entire scene, not just local to the problem area), we are forced
to aggressively mask clouds and other high flux regions. This adversely
impacts the data availability. This issue highlights the need for noise models
that are accurate across larger ranges of photon flux and observational
conditions. If a noise model existed that were able to encapsulate nonlinear
detection processes with heterogeneous targets, no such masking would be
required to process these scenes and PTV would have higher data availability.

The temperature error in the PTV processing method is generally lower than the
standard method across all overlapping altitudes (300 m to 3.5 km). We note
that the mean error tends to agree relatively well between the two methods
(where there are significant data points in the standard method). This bias of
approximately 1.5 K appears to be a result of a bias in the instrument itself.
It could be attributed to errors in the measured instrument wavelength, laser
spectral purity, or other aspects of the hardware that are either inaccurately
captured or modeled in the processing routines. It represents an area of
ongoing development and research for the MPD team. We note that this bias is
the limiting factor in the performance of PTV retrievals below 2.4 km (where
the standard deviation drops below 1 K). This indicates that if the bias in
the temperature estimate can be eliminated, the MPD PTV temperature error
would significantly improve. By contrast, the standard method error still
tends to be dominated by the standard deviation for most of the valid altitude
range, so reduction in biases would have a less significant impact on the data
product performance.

Overall, this analysis suggests that the global estimation of thermodynamic
variables using PTV produces lower data product error compared to the standard
method at the cost of slightly lower data availability below 2.2 km. However
PTV also delivers much higher data product availability above 2.2 km.
Importantly PTV seems better suited to capture temperature inversions in the
atmosphere, where the standard method struggles due to (typically) lower
maximum range and noise which is often difficult to distinguish from
atmospheric structure.

Figure 4: Evaluation of MPD PTV (blue) and Standard (orange) Temperature
retrievals based on data availability and agreement with sondes from M2HATS
campaign.

##  7 Conclusion

We have developed and demonstrated an advanced statistical signal processing
approach to retrieving thermodynamic profiles from MicroPulse DIAL
observations without the use of external data for processing or calibration.
The processing technique is an adaptation of PTV which has previously been
demonstrated for HSRL and water vapor DIAL processing. In this work we have
introduced a temperature processing component leveraging the oxygen DIAL
technique. By using GPU-accelerated processing on a cluster computer, we are
able to process a full day of data in approximately four hours at a cost of 50
GPU hours.

A comparison to sondes shows the PTV temperature estimates are an improvement
over the standard method. The approach not only extends the range of the
retrieval beyond the standard method but also reduces the errors –– to a limit
imposed from a bias in the hardware itself.

While observing temperature with oxygen DIAL is a novel concept in lidar
remote sensing, this work demonstrates its potential viability. Through
advances in signal processing and continued improvement in hardware
performance we expect continued improvement in the temperature data products.
This is an important step toward network-capable temperature monitoring using
lidar, as oxygen DIAL is well suited to low-power (eye-safe) semiconductor
laser-based instrument architectures that offer high reliability and low
maintenance.

## Acknowledgements

This material is based upon work supported by NSF AGS-1930907 and the NSF
National Center for Atmospheric Research, which is a major facility sponsored
by the National Science Foundation under Cooperative Agreement No. 1852977. We
would like to acknowledge high-performance computing support from Derecho [23]
provided by NCAR’s Computational and Information Systems Laboratory (CISL),
sponsored by the National Science Foundation.

## Disclosures

The authors declare no conflicts of interest.

## References

  *   * Stillwell et al. [2020] Stillwell, R.A., Spuler, S.M., Hayman, M., Repasky, K.S., Bunn, C.E.: Demonstration of a combined differential absorption and high spectral resolution lidar for profiling atmospheric temperature. Optics Express 28(1), 71–93 (2020) https://doi.org/10.1364/OE.379804
  * Bösenberg [1998] Bösenberg, J.: Ground-based differential absorption lidar for water-vapor and temperature profiling: methodology. Applied optics 37(18), 3845–3860 (1998)
  * Späth et al. [2020] Späth, F., Behrendt, A., Wulfmeyer, V.: Minimization of the Rayleigh-Doppler error of differential absorption lidar by frequency tuning: a simulation study. Opt. Express 28(20), 30324 (2020) https://doi.org/10.1364/OE.396568
  * Cady-Pereira et al. [2008] Cady-Pereira, K.E., Shephard, M.W., Turner, D.D., Mlawer, E.J., Clough, S.A., Wagner, T.J.: Improved daytime column-integrated precipitable water vapor from vaisala radiosonde humidity sensors. J. Atmos. Oceanic Technol. 25, 873–883 (2008) https://doi.org/10.1175/2007JTECHA1027.1
  * Bunn et al. [2018] Bunn, C.E., Repasky, K.S., Hayman, M., Stillwell, R.A., Spuler, S.M.: Perturbative solution to the two-component atmosphere DIAL equation for improving the accuracy of the retrieved absorption coefficient. Applied Optics 57(16), 4440 (2018) https://doi.org/10.1364/AO.57.004440
  * Repasky et al. [2019] Repasky, K.S., Bunn, C.E., Hayman, M., Stillwell, R.A., Spuler, S.M.: Modeling the performance of a diode laser-based (DLB) Micro-Pulse Differential Absorption Lidar (MPD) for temperature profiling in the lower troposphere. Opt. Express 27(23), 33543–33563 (2019) https://doi.org/10.1364/OE.27.033543
  * Marais et al. [2016] Marais, W.J., Holz, R.E., Hu, Y.H., Kuehn, R.E., Eloranta, E.E., Willett, R.M.: Approach to simultaneously denoise and invert backscatter and extinction from photon-limited atmospheric lidar observations. Applied Opt. 55, 8316–8334 (2016) https://doi.org/10.1364/AO.55.008316
  * Marais and Hayman [2022] Marais, W.J., Hayman, M.: Extending water vapor measurement capability of photon limited differential absorption lidars through simultaneous denoising and inversion. Atmos. Meas. Tech. 15, 5159–5180 (2022) https://doi.org/10.5194/amt-15-5159-2022
  * Nehrir et al. [2011] Nehrir, A., Repasky, K., Carlsten, J.: Eye-safe diode-laser-based micropulse differential absorption lidar (dial) for water vapor profiling in the lower troposphere. J. Atmos. Oceanic Technol. 28, 131–147 (2011) https://doi.org/10.1175/2010JTECHA1452.1
  * Spuler et al. [2015] Spuler, S.M., Repasky, K.S., Morley, B., Moen, D., Hayman, M., Nehrir, A.R.: Field-deployable diode-laser-based differential absorption lidar (dial) for profiling water vapor. Atmospheric Measurement Techniques 8(3), 1073–1087 (2015) https://doi.org/10.5194/amt-8-1073-2015
  * Weckwerth et al. [2016] Weckwerth, T.M., Weber, K.J., Turner, D.D., Spuler, S.M.: Validation of a Water Vapor Micropulse Differential Absorption Lidar (DIAL). Journal of Atmospheric and Oceanic Technology 33(11), 2353–2372 (2016) https://doi.org/10.1175/JTECH-D-16-0119.1
  * Spuler et al. [2021] Spuler, S.M., Hayman, M., Stillwell, R.A., Carnes, J., Bernatsky, T., Repasky, K.S.: Micropulse dial (mpd) – a diode-laser-based lidar architecture for quantitative atmospheric profiling. Atmos. Meas. Tech. 14, 4593–4616 (2021) https://doi.org/10.5194/amt-14-4593-2021
  * Hayman and Spuler [2017] Hayman, M., Spuler, S.: Demonstration of a diode-laser-based high spectral resolution lidar (HSRL) for quantitative profiling of clouds and aerosols. Optics Express 25(24), 1096 (2017) https://doi.org/10.1364/OE.25.0A1096
  * Lente and Ösz [2020] Lente, G., Ösz, K.: Barometric formulas: various derivations and comparisons to environmentally relevant observations. ChemTexts 6 (2020) https://doi.org/10.1007/s40828-020-0111-6
  * Hayman et al. [2019] Hayman, M., Stillwell, R.A., Spuler, S.M.: Fast computation of absorption spectra for lidar data processing using principal component analysis. Opt. Lett. 44(8), 1900 (2019) https://doi.org/10.1364/OL.44.001900
  * Binietoglou et al. [2016] Binietoglou, I., Giampouras, P., Belegante, L.: Linear approximation of Rayleigh–Brillouin scattering spectra. Applied Optics 55(27), 7707–7711 (2016) https://doi.org/10.1364/AO.55.007707
  * Horoshko et al. [2017] Horoshko, D.B., Chizhevsky, V.N., Kilin, S.Y.: Afterpulsing model based on the quasi-continuous distribution of deep levels in single-photon avalanche diodes. J. of Modern Optics 64(2), 191–195 (2017) https://doi.org/10.1080/09500340.2016.1220643
  * Ziarkash et al. [2018] Ziarkash, A.W., Joshi, S.K., Stipčevi/’c, M., Ursin, R.: Comparative study of afterpulsing behavior and models in single photon counting avalanche photo diode detectors. Scientific Reports 8(1), 1–8 (2018) https://doi.org/%****␣main.bbl␣Line␣325␣****10.1038/s41598-018-23398-z
  * Hayman et al. [2023] Hayman, M., Stillwell, R.A., Carnes, J., Kirchhoff, G.J., Spuler, S.M., Thayer, J.P.: 2D signal estimation for sparse distributed target photon counting data. arXiv:2311.18037 (2023) https://doi.org/10.48550/arXiv.2311.18037 arXiv:2311.18037 [eess.SP]
  * Harmany et al. [2012] Harmany, Z.T., Marcia, R.F., Willett, R.M.: This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms—theory and practice. IEEE Transactions on Image Processing 21, 1084–1096 (2012) https://doi.org/10.1109/TIP.2011.2168410
  * Beck and Teboulle [2009] Beck, A., Teboulle, M.: Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Transactions on Image Processing 18, 2419–2434 (2009) https://doi.org/10.1109/TIP.2009.2028250
  * Paszke et al. [2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc., ??? (2019). http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * Computational and Information Systems Laboratory, CISL [2023] Computational and Information Systems Laboratory, CISL: Derecho: HPE Cray EX System (NCAR Community Computing) (2023). https://doi.org/10.5065/qx9a-pg09
  * Oh et al. [2013] Oh, A.K., Harmany, Z.T., Willett, R.M.: Logarithmic total variation regularization for cross-validation in photon-limited imaging. In: 2013 IEEE International Conference on Image Processing, pp. 484–488 (2013)
  * Marais [2017] Marais, W.J.: Poisson Inverse and Denoising Problems in Atmospheric Lidar Imaging. PhD thesis, The University of Wisconsin - Madison (January 2017)
  * Hayman et al. [2020] Hayman, M., Stillwell, R.A., Spuler, S.M.: Optimization of linear signal processing in photon counting lidar using poisson thinning. Opt. Lett. 45, 5213–5216 (2020) https://doi.org/10.1364/OL.396498
  * Amini et al. [2020] Amini, A., Schwarting, W., Soleimany, A., Rus, D.: Deep evidential regression. Advances in Neural Information Processing Systems 33, 20004–20016 (2020)
  * Schreck et al. [2023] Schreck, J.S., II, D.J.G., Becker, C., Chapman, W.E., Elmore, K., Gantos, G., Kim, E., Kimpara, D., Martin, T., Molina, M.J., Pryzbylo, V.M., Radford, J., Saavedra, B., Willson, J., Wirz, C.: Evidential deep learning: Enhancing predictive uncertainty estimation for earth system science applications. arXiv:2309.13207 (2023) https://doi.org/10.48550/arXiv.2309.13207
  * [29] Copernicus Climate Change Service (C3S): ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS). Accessed Nov. 2023. https://doi.org/10.26023/WKM7-HNCF-FX0B . https://cds.climate.copernicus.eu/cdsapp#!/home
  * Hersbach et al. [2020] Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková́, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut, J.-N.: The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146(730), 1999–2049 (2020) https://doi.org/10.1002/qj.3803
  * [31] NCAR/EOL In-situ Sensing Facility and Clemson University and California State University at Chico: M2HATS: ISS Radiosonde Data. Version 1.0. UCAR/NCAR - Earth Observing Laboratory. Accessed 2023. https://doi.org/10.26023/WKM7-HNCF-FX0B
  * [32] NCAR/EOL Remote Sensing Facility: M2HATS: NCAR MPD data. UCAR/NCAR - Earth Observing Laboratory. Accessed 2023. https://doi.org/10.26023/2286-1446-JV0E
  * Fix et al. [2010] Fix, A., Ehret, G., Löhring, J., Hoffmann, D., Alpers, M.: Water vapor differential absorption lidar measurements using a diode-pumped all-solid-state laser at 935 nm. Applied Physics B (2010) https://doi.org/10.1007/s00340-010-4310-5
  * Piironen and Eloranta [1994] Piironen, P., Eloranta, E.W.: Demonstration of a high-spectral-resolution lidar based on an iodine absorption filter. Optics letters 19(3), 234 (1994)


# Can electromagnetic charge inhabit in Rastall gravity?

Bobir Toshmatov<EMAIL_ADDRESS>Zdeněk Stuchlík
<EMAIL_ADDRESS>Bobomurat Ahmedov<EMAIL_ADDRESS>New
Uzbekistan University, Mustaqillik Ave. 54, Tashkent 100007, Uzbekistan
Research Centre of Theoretical Physics and Astrophysics, Institute of Physics,
Silesian University in Opava, Bezručovo nám. 13, CZ-74601 Opava, Czech
Republic  Ulugh Beg Astronomical Institute, Astronomy 33, Tashkent 100052,
Uzbekistan  IFAR, “TIIAME” NRU, Kori Niyoziy 39, Tashkent 100000, Uzbekistan
National University of Uzbekistan, Tashkent 100174, Uzbekistan

###### Abstract

One of the eminent generalizations of theory of general relativity is the
Rastall gravity which was constructed based on the assumption of the non-
conserved energy-momentum tensor of the matter field. Despite in the
literature several solutions of black holes in the Rastall gravity coupled to
the electromagnetic field have been presented, in the current paper we argue
that the Rastall gravity with non-conserved energy-momentum tensor (with
$\lambda\neq 0$ and $R\neq 0$) cannot couple to the electrodynamics, i.e., the
electromagnetically charged black hole solution cannot be obtained in this
case. This statement is adequate for both linear and nonlinear electrodynamics
with the electric, magnetic, or dyonic charges coupled to the Rastall gravity.

###### keywords:

General relativity , Rastall gravity , electrodynamics

††journal: Journal of LaTeX Templates

##  1 Introduction

General relativity is the prevailing theory of gravity, which describes the
force as a result of massive objects warping spacetime. It is well-known that
the theory of general relativity is one of the simplest and most elegant
theories of gravitation that has explained several hidden mysteries of the
nature and universe. The correctness of general relativity has been verified
by several experimental tests in the weak field regimes [1]. However, its
correctness in the strong gravitational fields is still under debate, as
recent observational breakthrough events such as detection of the
gravitational waves from the coalescence of two black holes or neutron stars
in binaries [2, 3, 4, 5, 6] and obtained images of the supermassive black
holes in the centers of Milky Way [7] and M87 [8] galaxies have still left
some room for the alternative and modified theories of gravity. Furthermore,
there are still big numbers of shortcomings in theoretical and observational
aspects of the science, such as accelerated universe, the existence of dark
matter and energy, on account of which modification of general relativity is
inevitable.

There are possible modified and alternative theories of gravity that have been
constructed in mainly two ways. In the first method, the fundamental
assumptions of general relativity are still kept but new additional terms are
added to the Lagrangian density of the general relativity on account of which
all field equations are modified with additional terms. In the second method,
the fundamental assumptions of general relativity are changed. The Rastall
gravity belongs to the latter group. It is clear that one of the fundamental
assumptions of general relativity is the null covariant divergence of the
energy-momentum tensor. However, Rastall followed the second method and
attempted to modify the theory of general relativity by assuming that the
covariant derivative of the energy-momentum tensor varies directly as the
derivative of the Ricci scalar [9, 10]. Clearly, it explicitly represents a
nonminimal coupling between matter and geometry in curved spacetimes. In flat
spacetimes, the energy-momentum tensor satisfies a conservation law as the
curvature of spacetime vanishes. As a result, the field equations become
equivalent to the vacuum Einstein equations. Even in the curved spacetime with
matter field, the field equations of the Rastall gravity can be as the
Einstein equations in terms of the effective energy-momentum tensor which we
will present later in the next section. From this point of view, The Rastall
gravity can be considered as the equivalent to general relativity as stated in
[11]. Moreover, these theories are equivalent in vacuum and in non-vacuum case
with conserved energy-momentum tensors. However, in the matter field with non-
conserved energy-momentum tensor, these two theories are different. In spite
of that, there are still various issues that these two theories result from
different consequences [12]. Over the last decades, Rastall gravity has gained
attention of scientists and as a result of it, a vast number of papers related
to the new solutions, or properties of the solutions of black holes in Rastall
gravity, have been published [13, 14, 15, 16, 17]. Among them there are the
ones related to the charged solutions [18, 19, 20, 21, 22, 23, 24].
Unsurprisingly, in this paper we present results that question possibility of
obtaining charged black hole solutions in Rastall gravity.

We study the possibility of the construction of the black hole solution in
Rastall gravity coupled to the electrodynamics. The paper is organized as
follows: In section 2 we revisit Rastall gravity and equations of motion.
Section 3 is devoted to a system of Rastall gravity coupled to the Maxwell
electrodynamics. In section 4 we repeat the calculations presented in the
previous section for the nonlinear electrodynamics. Finally, in sections 5 and
6 we discuss and summarize the main results obtained in the paper. Throughout
the paper, we use the geometrized units in which the Newtonian gravitational
constant $G_{N}$, and speed of light $c$ are set as $G_{N}=c=1$.

##  2 Rastall gravity revisited

According to Rastall’s theory of gravity, the conservation of the energy-
momentum tensor, $T^{\mu\nu}$, is assumed to be given by [9]

$\displaystyle\nabla_{\nu}T^{\mu\nu}=\lambda\nabla^{\mu}R\ ,$ (1)

where $\lambda$ is a free parameter, so-called Rastall coupling parameter and
$R$ is a Ricci scalar. The nondivergence-free energy-momentum assumption that
varies directly as the derivative of Ricci scalar, represents a nonminimal
coupling between matter and geometry. In other words, if there is no matter
field (vacuum), the theory of general relativity is recovered. Thus, Rastall
gravity being generalization of general relativity, it includes general
relativity in which the energy-momentum tensor conserved that corresponds to
$\lambda=0$ or $R={\rm const.}$ (inner region in Fig. 1) and the case in which
the energy-momentum tensor is not conserved that corresponds to $\lambda\neq
0$ and $R\neq{\rm const.}$ (outer region in Fig. 1).

Figure 1: Schematic picture of the Rastall gravity that generalizes general
relativity (GR) is presented.

According to the assumption of a non-divergence-free energy-momentum tensor in
the curved spacetime, Rastall obtained a consistent set of field equations

$\displaystyle\tilde{G}_{\mu\nu}\equiv G_{\mu\nu}+\lambda\kappa
g_{\mu\nu}R=\kappa T_{\mu\nu}\ ,$ (2)

with the Einstein tensor defined by

$\displaystyle G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\ ,$ (3)

where $\kappa=8\pi$. In the case of general relativity which corresponds to
$\lambda=0$, equation (2) reduces to the well-known Einstein’s field equations
$G_{\mu\nu}=\kappa T_{\mu\nu}$. The trace of the field equation (2) results
the following equation:

$\displaystyle R=\frac{\kappa}{4\kappa\lambda-1}T\
,\qquad\kappa\lambda\neq\frac{1}{4}\ .$ (4)

In terms of the expression (4), one can rewrite the field equation (2) in the
alternative form as

$\displaystyle G_{\mu\nu}=\kappa\tilde{T}_{\mu\nu}\ ,$ (5)

where

$\displaystyle\tilde{T}_{\mu\nu}=T_{\mu\nu}-\frac{\kappa\lambda}{4\kappa\lambda-1}g_{\mu\nu}T\
.$ (6)

One can see from (5) that the Rastall term in the field equation can be
considered as a new additional matter field to the right-hand side of the
Einstein equations. Therefore, in this sense, Rastall gravity can be
considered formally equivalent to general relativity [11]. Now by choosing the
appropriate matter field and energy-momentum tensor and fixing the background
ansatz, one can construct the solution in Rastall gravity. For simplicity, we
choose the static, spherically symmetric spacetime as

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\ .$ (7)

For the spacetime line element (7) non-zero components of the left hand side
of the field equations (2) are given as

$\displaystyle\tilde{G}^{t}_{\ t}=\frac{rf^{\prime}+f-1}{r^{2}}+\lambda\kappa
R=\tilde{G}^{r}_{\ r}\ ,$ $\displaystyle\tilde{G}^{\theta}_{\
\theta}=\frac{f^{\prime\prime}}{2}+\frac{f^{\prime}}{r}+\lambda\kappa
R=\tilde{G}^{\phi}_{\ \phi}\ ,$ (8)

where the Ricci scalar equals

$\displaystyle
R=-f^{\prime\prime}-\frac{2\left(2rf^{\prime}+f-1\right)}{r^{2}}\ .$ (9)

In the following section, we obtain such solutions in Rastall gravity for
electrically charged matter fields. The solutions will be obtained through
mathematical calculations and analysis in the framework of the newly developed
formalism.

##  3 Maxwell electrodynamics coupled to Rastall gravity

Since in the previous section all equations related to the geometry of the
spacetime are given, in this section we solve these equations by setting the
energy-momentum tensor of the matter field by the Maxwell electrodynamics:

$\displaystyle T_{\mu\nu}=\frac{2}{\kappa}\left(F_{\mu}^{\
\alpha}F_{\nu\alpha}-\frac{1}{4}g_{\mu\nu}F_{\delta\gamma}F^{\delta\gamma}\right)\
,$ (10)

where the electromagnetic field tensor is determined by the vector potential
$A_{\mu}$ as $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. Their
explicit forms are identical for both electrically and magnetically charged
cases as

$\displaystyle T^{t}_{\ t}=T^{r}_{\ r}=-\frac{Q^{2}}{\kappa r^{4}}\ ,\qquad
T^{\theta}_{\ \theta}=T^{\phi}_{\ \phi}=\frac{Q^{2}}{\kappa r^{4}}\ .$ (11)

By simplifying the field equations (5) with (2) and (11), we obtain the
following two independent differential equations:

$\displaystyle\frac{1-rf^{\prime}-f}{r^{2}}-\lambda\kappa
R=\frac{Q^{2}}{r^{4}},$ (12)
$\displaystyle\frac{f^{\prime\prime}}{2}+\frac{f^{\prime}}{r}+\lambda\kappa
R=\frac{Q^{2}}{r^{4}}\ .$ (13)

If the above differential equations are solved, separately, in general
relativity limit ($\lambda=0$), both equations (12) and (13) would result in
the well-known Reissner-Nordström solution

$\displaystyle f(r)=1-\frac{2M}{r}+\frac{Q_{e}^{2}}{r^{2}}\ .$ (14)

From the fact that the Reissner-Nordström solution (14) is obtained by solving
both differential equations (12) and (13), we are convinced that the linear
electrodynamics can couple to the general relativity.

In the next step, we solve the differential equations (12) and (13) separately
in Rastall gravity ($\lambda\neq 0$). By solving the first differential
equation (12), we obtain the metric function as following:

$\displaystyle
f(r)=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}+cr^{-2+\frac{1}{\lambda\kappa}}\ ,$
(15)

while from the second differential equation (13) we arrive at the metric
function

$\displaystyle
f(r)=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}+cr^{\frac{4\lambda\kappa}{1-2\lambda\kappa}}\
,$ (16)

where $c$ is the integration constant. For $f(r)$ to be the solution, both
metric functions (15) and (16) must be identical. This can happen only in the
following two cases:

  * (i)

the integration constant is zero ($c=0$). In this case both solutions are
reduced to the Reissner-Nordström black hole solution. We should note that the
Reissner-Nordström solution is also the solution of the Rastall gravity only
because of the fact that it is a solution of general relativity coupled to the
Maxwell electrodynamics. We strengthen our this statement with the fact that
for the Reissner-Nordström solution the Ricci scalar vanishes ($R=0$) and
consequently, energy-momentum tensor becomes conserved and all field equations
(including (1), (2) and (5)) take the form of the ones of pure general
relativity;

  * (ii)

exponents of $r$ in the last terms of both metric functions (15) and (16) are
equal when $\lambda\kappa=1/4$. In this case the metric function takes the de-
Sitter-like form as (it was obtained in [21])

$\displaystyle f(r)=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}+cr^{2}\ .$ (17)

However, even if the metric function (17) satisfies both differential
equations (12) and (13), it does not satisfy the trace of the field equation
(4), as $\lambda\kappa=1/4$ is the forbidden value (see also [9]). This means
that even though the metric function satisfies both differential equations, it
fails to satisfy a crucial condition that needs to be met in order to be
considered a valid solution. The restricted value of $\lambda\kappa=1/4$ is an
essential requirement that the metric function must fulfill in order to be
considered a solution to the field equations. Without satisfying this
condition, it can be concluded that the metric function is not a valid
solution. Therefore, it is not a solution.

Thus, as we have discussed in above two cases that the charged black hole
solution can exist only in general relativity that is a subset of Rastall
gravity.

##  4 Nonlinear electrodynamics coupled to Rastall gravity

We have shown in the previous section that Maxwell electrodynamics cannot
couple to Rastall gravity to produce charged black hole solution. As a
continuation, we generalize the scenario to the nonlinear electrodynamics
coupled to Rastall gravity. We again adopt the spacetime in a static,
spherically symmetric form. The energy-momentum tensor of the nonlinear
electrodynamics is given by

$\displaystyle T_{\mu\nu}=\frac{2}{\kappa}\left(\mathscr{L}_{F}F_{\mu}^{\
\alpha}F_{\nu\alpha}-\frac{1}{4}g_{\mu\nu}\mathscr{L}\right)\ ,$ (18)

where $\mathscr{L}_{F}$ denotes a derivative of the Lagrangian density of the
nonlinear electrodynamics with respect to invariant, $F\equiv
F_{\mu\nu}F^{\mu\nu}$, $\mathscr{L}_{F}=\partial_{F}\mathscr{L}$. In general,
the electromagnetic four-potential, $A_{\mu}$, is written for the spherically
symmetric spacetimes in terms of the 1-forms of Carter tetrad for the
background line element (7) as [25]

$\displaystyle
A=-\frac{\varphi}{\sqrt{f}}\omega^{t}-\frac{Q_{m}\cot\theta}{r}\omega^{\phi}\
.$ (19)

where the 1-forms of the Carter tetrad is defined as

$\displaystyle\omega^{t}=\sqrt{f}dt\ ,\qquad\omega^{r}=\frac{1}{\sqrt{f}}dr\
,$ $\displaystyle\omega^{\theta}=rd\theta\ ,\qquad\omega^{\phi}=r\sin\theta
d\phi\ .$ (20)

Thus, the electromagnetic field four-potential (19) can be written explicitly
as

$\displaystyle A=\varphi(r)dt-Q_{m}\cos\theta d\phi\ ,$ (21)

where the time and azimuthal components of the potential (21) correspond to
the electrically and magnetically charged black holes, respectively, as
$\varphi$ and $Q_{m}$ are electric scalar potential and magnetic charge,
respectively. Unlike the case of linear electrodynamics in a previous section,
in the case of nonlinear electrodynamics electrically and magnetically charged
electromagnetic fields do not generate the identical solutions [26, 27, 28,
29, 30] and their characteristic oscillations are not isospectral [31, 32, 33,
34, 35]. Therefore, below we present a construction of the electrically and
magnetically charged black hole solutions in Rastall gravity coupled to the
nonlinear electrodynamics, separately.

###  4.1 Electrically charged black holes

In the case of the electrically charged black hole, the time component of the
four-potential (21) of the electromagnetic field survives as

$\displaystyle A_{\mu}=\\{\varphi(r),\ 0,\ 0,\ 0\\}\ .$ (22)

After calculating the components of the electromagnetic field tensor, we find
the electromagnetic field invariant as

$\displaystyle F_{e}=-2\varphi^{\prime 2}\ .$ (23)

The nonvanishing components of the energy-momentum tensor of the nonlinear
electrodynamics are found to be

$\displaystyle T^{t}_{\ t}$ $\displaystyle=$
$\displaystyle-\frac{\mathscr{L}-2\mathscr{L}_{F}F_{e}}{\kappa}=T^{r}_{\ r}\
,$ $\displaystyle T^{\theta}_{\ \theta}$ $\displaystyle=$
$\displaystyle-\frac{\mathscr{L}}{\kappa}=T^{\phi}_{\ \phi}\ .$ (24)

By solving the field equations (2) with (2) and (4.1), we obtain the following
expressions for the Lagrangian density of the nonlinear electrodynamics:

$\displaystyle\mathscr{L}$ $\displaystyle=$
$\displaystyle-f^{\prime\prime}-\frac{2f^{\prime}}{r}-2\lambda\kappa R,$ (25)
$\displaystyle\mathscr{L}_{F}$ $\displaystyle=$
$\displaystyle-\frac{r^{2}f^{\prime\prime}-2f+2}{2r^{2}F_{e}}\ .$ (26)

It can be easily noticed from the expressions (25) and (26) that the
Lagrangian density of the nonlinear electrodynamics is the function of Rastall
coupling parameter ($\partial\mathscr{L}/\partial\lambda\neq 0$), while its
derivative with respect to the electromagnetic invariant, $\mathscr{L}_{F}$,
is independent of the Rastall coupling parameter
($\partial\mathscr{L}_{F}/\partial\lambda=0$). At this point, we come across a
problem with the equations, as expressions of the Lagrangian density (25) and
(26) contradict each other. Mathematically, if the Lagrangian density is the
function of the parameter $\lambda$, its derivative with respect to the
electromagnetic field invariant, $F_{e}$, which is not a function of the
parameter $\lambda$ ($\partial F/\partial\lambda=0$) must also be the function
of $\lambda$, as $\mathscr{L}_{F}=\mathscr{L}^{\prime}/F^{\prime}$ with the
prime being the derivative with respect to the radial coordinate, $r$.
However, it is not the case here. In short, if the differential equations (25)
and (26) are solved for any Lagrangian density of the nonlinear
electrodynamics separately, one would obtain the metric function with the
Rastall term in the first solution (solution of (25)), while in the second one
(solution of (26)) there is no term with the Rastall parameter. In other
words, if one solves the differential equations for the nonlinear
electrodynamics Lagrangian density separately, the first solution (the
solution of (25)) would include the Rastall term, while the second solution
(the solution of (26)) would not include it. Therefore, there are two cases to
address this contradiction:

  * (i)

$\lambda=0$ that corresponds to general relativity coupled to the nonlinear
electrodynamics, as the energy-momentum tensor becomes conserved. The obtained
solution can be considered the solution of the Rastall gravity only because of
the fact that general relativity is a subset of the Rastall gravity;

  * (ii)

$R=0$ that corresponds to general relativity coupled to the nonlinear
electrodynamics. However, since $R=0$ implies the solution to be Reissner-
Nordström-like and it is the solution of general relativity coupled to the
linear electrodynamics, this system does not provide any solution not only in
the general relativity coupled to the nonlinear electrodynamics, but also in
Rastall gravity coupled to the nonlinear electrodynamics.

Thus, in both cases the final solution would be the one in general relativity
and not in Rastall gravity with nonconserved energy-momentum tensor. Thus, we
are convinced that one cannot obtain the electrically charged black hole
solution in Rastall gravity with $\lambda\neq 0$ and $R\neq 0$ coupled to the
nonlinear electrodynamics.

###  4.2 Magnetically charged black holes

In the magnetically charged black hole solution, the four potential of the
electromagnetic field (21) is given by

$\displaystyle A_{\mu}=\\{0,0,0,-Q_{m}\cos\theta\\}\ .$ (27)

In this case, the electromagnetic field invariant takes the following form:

$\displaystyle F_{m}=\frac{2Q_{m}^{2}}{r^{4}}\ .$ (28)

The independent components of the energy-momentum tensor of the nonlinear
electrodynamics are given as

$\displaystyle T^{t}_{\ t}$ $\displaystyle=$
$\displaystyle-\frac{\mathscr{L}}{\kappa}=T^{r}_{\ r}\ ,$ $\displaystyle
T^{\theta}_{\ \theta}$ $\displaystyle=$
$\displaystyle-\frac{\mathscr{L}-\mathscr{L}_{F}F_{m}}{\kappa}=T^{\phi}_{\
\phi}\ .$ (29)

By solving the field equations (2) with (2) and (4.2), we obtain the following
expressions for the Lagrangian density:

$\displaystyle\mathscr{L}$ $\displaystyle=$
$\displaystyle\frac{2\left(1-rf^{\prime}-f\right)}{r^{2}}-2\lambda\kappa R,$
(30) $\displaystyle\mathscr{L}_{F}$ $\displaystyle=$
$\displaystyle\frac{r^{2}f^{\prime\prime}-2f+2}{2r^{2}F_{m}}\ .$ (31)

Thus, the situation is the same as in the case of the electrically charged
black hole solution in a previous subsection, and consequently the conclusion
is also similar. Since the results in both electrically and magnetically
charged cases are the same, it is obvious that it remains unchanged in the
dyonically charged case too. Therefore, we do not here report this case.

##  5 Discussion

In previous sections, we have shown by solving the field equations (2) that
the electromagnetic field cannot couple to Rastall gravity in black hole
solutions. However, in that process, we have not used the equation on non-
conserved energy-momentum tensor in Rastall gravity (1) and the Maxwell
equations for the nonlinear electrodynamics which are given by

$\displaystyle\nabla_{\mu}\left(\mathscr{L}_{F}F^{\mu\nu}\right)=0\ .$ (32)

These two equations contradict each other, as equation (32) requires the
covariant derivative of the right hand side of the field equation (2) to
vanish, while the covariant derivative of the left hand side of (2) is nonzero
via (1). To overcome this problem, we must consider either $\lambda=0$ or
$R=0$ cases. The former case represents the absence of the Rastall term in the
equations of motion, i.e., it corresponds to general relativity with Einstein
equations coupled to the nonlinear electrodynamics. This type of solutions are
obtained and well studied in the vast literature such as [26, 27, 28]. In the
latter case when the Ricci scalar (or trace of the energy-momentum tensor)
vanishes, there is only one solution that corresponds to the Reissner-
Nordström-like one as

$\displaystyle R=0\quad\Rightarrow\quad f(r)=1-\frac{2M}{r}+\frac{c}{r^{2}}\
.$ (33)

with $c$ being the integration constant. From this solution it turns out that
to have all field equations satisfied, the following equation also must give
the same metric function as (33):

$\displaystyle R_{\mu\nu}=\kappa T_{\mu\nu}\ .$ (34)

Despite the field equations have changed, they still do not represent the
general Rastall gravity containing nonzero $\lambda$. The solution of the
field equation (34) for the Maxwell electrodynamics represents the Reissner-
Nordström black hole solution which is again independent from the Rastall
parameter.

By solving equations (34) for the nonlinear electrodynamics (18) with metric
function (33), we obtain that the Lagrangian density of the nonlinear
electrodynamics is linearly related to the electromagnetic field invariant,
$F$, as

$\displaystyle\mathscr{L}\propto F,\qquad\mathscr{L}_{F}\propto{\rm constant}\
,$ (35)

whose solution is independent from the Rastall parameter. From this approach
also we have confirmed that Rastall gravity (with $\lambda\neq 0$ and $R\neq
0$) cannot coexist with the electrodynamics.

##  6 Conclusion

Rastall gravity is a simple mathematical generalization of general relativity
on account of the non-minimal coupling between the geometry and energy-
momentum tensor of the matter field. The original physical motivation of this
theory was that the local energy-momentum conservation law in the flat
spacetime does not necessarily imply its conservation in the curved spacetime
with matter field. This on turn implies several new interesting properties in
the gravitational theory. One of such properties has been addressed in this
paper. Despite in many literature new black hole solutions have been obtain
within the framework of Rastall gravity coupled to the electrodynamics
(sometimes some additional second matter field is also considered), we have
shown here that in this system one cannot obtain the electromagnetically
(electrically, magnetically, or dyonically) charged black hole solutions in
Rastall gravity with non-conserved energy-momentum tensor ($\lambda\neq 0$ and
$R\neq 0$) coupled not only to the linear but also to the nonlinear
electrodynamics. In other words, the electromagnetically charged solutions of
black hole cannot be obtained in Rastall gravity corresponding to the outer
region of Fig. 1.

## Acknowledgement

BT acknowledges the support of program ”CZ.02.2.69/0.0/0.0/18-053/0017871:
Podpora mezinárodní mobility na Slezské univerzitě v Opavě” at the Institute
of Physics, Silesian University in Opava. Authors acknowledge the support of
Ministry of Innovative Development of the Republic of Uzbekistan Grants No.
F-FA-2021-432 and MRB-2021-527.

## References

  * [1] Will C M 2014 Living Rev. Rel. 17 4 (Preprint 1403.7377)
  * [2] The LIGO Scientific Collaboration and the Virgo Collaboration 2016 Phys. Rev. Lett. 116 241103 (Preprint 1606.04855)
  * [3] The LIGO Scientific Collaboration and the Virgo Collaboration 2017 Phys. Rev. Lett. 118 221101 (Preprint 1706.01812)
  * [4] The LIGO Scientific Collaboration and the Virgo Collaboration 2017 Phys. Rev. Lett. 119 141101 (Preprint 1709.09660)
  * [5] The LIGO Scientific Collaboration and the Virgo Collaboration 2017 Astrophys. J. Lett. 851 L35 (Preprint 1711.05578)
  * [6] The LIGO Scientific Collaboration and the Virgo Collaboration 2017 Phys. Rev. Lett. 119 161101 (Preprint 1710.05832)
  * [7] The Event Horizon Telescope Collaboration 2022 Astrophys. J. Lett. 930 L12
  * [8] The Event Horizon Telescope Collaboration 2019 Astrophys. J. 875 L1 (Preprint 1906.11238)
  * [9] Rastall P 1972 Phys. Rev. D 6 3357–3359
  * [10] Rastall P 1976 Can. J. Phys. 54 66–75
  * [11] Visser M 2018 Phys. Lett. B 782 83–86 (Preprint %****␣main.bbl␣Line␣50␣****1711.11500)
  * [12] Darabi F, Moradpour H, Licata I, Heydarzade Y and Corda C 2018 Eur. Phys. J. C 78 25 (Preprint 1712.09307)
  * [13] Batista C E M, Daouda M H, Fabris J C, Piattella O F and Rodrigues D C 2012 Phys. Rev. D 85 084008 (Preprint 1112.4141)
  * [14] Oliveira A M, Velten H E S, Fabris J C and Casarini L 2015 Phys. Rev. D 92 044020 (Preprint 1506.00567)
  * [15] Ghosh S, Dey S, Das A, Chanda A and Paul B C 2021 JCAP 07 004 (Preprint 2102.01524)
  * [16] Maulana H and Sulaksono A 2019 Phys. Rev. D 100 124014
  * [17] El Hanafy W 2022 Astrophys. J. 940 51 (Preprint 2209.10287)
  * [18] Sakti M F A R, Suroso A and Zen F P 2020 Annals Phys. 413 168062 (Preprint 1901.09163)
  * [19] Sakti M F A R, Suroso A, Sulaksono A and Zen F P 2022 Phys. Dark Univ. 35 100974 (Preprint 2110.03525)
  * [20] Guo S, He K J, Li G R and Li G P 2021 Class. Quant. Grav. 38 165013 (Preprint 2205.07242)
  * [21] Nashed G G L 2022 Universe 8 510 (Preprint 2210.05376)
  * [22] Gogoi D J and Goswami U D 2021 Phys. Dark Univ. 33 100860 (Preprint 2104.13115)
  * [23] Shao C Y, Tan Y J, Shao C G, Lin K and Qian W L 2022 Chin. Phys. C 46 105103 (Preprint 2208.09988)
  * [24] Narzilloev B, Hussain I, Abdujabbarov A and Ahmedov B 2022 Eur. Phys. J. Plus 137 645 (Preprint 2205.11760)
  * [25] Znajek R L 1977 Mon. Not. Roy. Astron. Soc. 179 457–472
  * [26] Bronnikov K 2000 Phys. Rev. Lett. 85 4641
  * [27] Bronnikov K A 2001 Phys. Rev. D 63 044005 (Preprint gr-qc/0006014)
  * [28] Fan Z Y and Wang X 2016 Phys. Rev. D 94 124027 (Preprint 1610.02636)
  * [29] Bronnikov K A 2017 Phys. Rev. D 96 128501 (Preprint 1712.04342)
  * [30] Toshmatov B, Stuchlík Z and Ahmedov B 2018 Phys. Rev. D 98 028501 (Preprint 1807.09502)
  * [31] Toshmatov B, Stuchlík Z, Schee J and Ahmedov B 2018 Phys. Rev. D 97 084058 (Preprint 1805.00240)
  * [32] Toshmatov B, Ahmedov B and Malafarina D 2021 Phys. Rev. D 103 024026 (Preprint 2101.05496)
  * [33] Moreno C and Sarbach O 2003 Phys. Rev. D 67 024028 (Preprint gr-qc/0208090)
  * [34] Nomura K, Yoshida D and Soda J 2020 Phys. Rev. D 101 124026 (Preprint 2004.07560)
  * [35] Nomura K and Yoshida D 2022 Phys. Rev. D 105 044006 (Preprint 2111.06273)