# Reaching high accuracy for energetic properties at second-order perturbation
cost by merging self-consistency and spin-opposite scaling

Nhan Tri Tran  University of Science, Vietnam National University, Ho Chi Minh
City, Vietnam  Hoang Thanh Nguyen  Institute of Applied Mechanics and
Informatics, Vietnam Academy of Science and Technology, Ho Chi Minh City,
Vietnam  Lan Nguyen Tran<EMAIL_ADDRESS>Department of Physics,
International University, Ho Chi Minh City, Vietnam  Vietnam National
University, Ho Chi Minh City, Vietnam

###### Abstract

Quantum chemical methods dealing with challenging systems while retaining low
computational costs have attracted attention. In particular, many efforts have
been devoted to developing new methods based on the second-order perturbation
that may be the simplest correlated method beyond Hartree-Fock. We have
recently developed a self-consistent perturbation theory named one-body
Møller-Plesset second-order perturbation theory (OBMP2) and shown that it can
resolve issues caused by the non-iterative nature of standard perturbation
theory. In the present work, we extend the method by introducing the spin-
opposite scaling to the double-excitation amplitudes, resulting in the O2BMP2
method. We assess the O2BMP2 performance on the triple-bond N2 dissociation,
singlet-triplet gaps, and ionization potentials. O2BMP2 performs much better
than standard MP2 and reaches the accuracy of coupled-cluster methods in all
cases considered in this work.

Second-order Møller-Plesset perturbation theory (MP2) on Hartree-Fock (HF)
orbitals1 is the simplest correlated wave-function method. Its accuracy
depends on the quality of reference wave functions, in particular, for open-
shell systems2, 3. To bypass the issue of poor references, many research
groups have actively developed orbital-optimized MP2 (OOMP2) and its spin-
scaled variants4, 5, 6, 7, 8, 9. In these methods, orbitals are optimized by
minimizing the Hylleraas functional. OOMP2 and its variants have outperformed
standard MP2 calculations for numerous properties. Apart from wave-function
methods, double-hybrid density functional (DHF) theory, in which a scaled
perturbative correction is performed on top of hybrid density functional
calculations, has attracted significant attention. These functionals are
considered the fifth rung of the DFT Jacob’s ladder and have been shown to
outperform conventional functionals in many cases10, 11.

It is well-known that perturbation theory is inadequate for multi-reference
systems, and the perturbative correlation energy diverges due to small gaps of
orbital energies. To eliminate these issues, several regularization schemes
that modify the MP2 amplitude with a function damping any divergent or
excessively large correlations have been recently developed 12, 13, 14, 15,
16. It has been shown that regularized (orbital-optimized) MP217, 18, 19 can
outperform standard MP2 across relevant chemical problems. In the meantime,
numerous efforts are devoted to developing alternative approaches to resolving
the abovementioned issues. These methods include Brillouin-Wigner perturbation
theory (BWPT) and its size-consistent variant 20, 21, retaining the excitation
degree MP2 (REMP2) and its orbital-optimized variant 22, 23. Empirical spin-
scaled methods, such as spin-component scaling (SCS) and spin-opposite scaling
(SOS), have also been widely used to improve the performance of perturbation
theory 24. Noticeably, SOS-MP2 does not only often improve the accuracy of
MP2, but it is also less costly ($N^{4}$) than standard MP2 ($N^{5}$ ).

In general, developing new methods based on low-cost perturbation theory able
to deal with challenging systems is still highly desirable. Recently, we have
developed a new self-consistent perturbation theory named one-body MP2
(OBMP2)25, 26, 27, 28. The key idea of OBMP2 is the use of canonical
transformation29, 30, 31, 32, 33, 34 followed by the cumulant approximation35,
36, 37, 38 to derive an effective one-body Hamiltonian. The resulting OBMP2
Hamiltonian is a sum of the standard Fock operator and a one-body correlation
MP2 potential. Molecular orbitals and orbital energies are relaxed in the
presence of correlation by diagonalizing correlated Fock matrix. The double-
excitation MP2 amplitudes are then updated using those new molecular orbitals
and orbital energies, resulting in a self-consistency. We have shown that the
self-consistency of OBMP2 can resolve issues caused by the non-iterative
nature of standard MP2 calculations for open-shell systems 26, 27. It is also
surprising that OBMP2 does not suddenly break down in bond stretching 28.

In this work, we present the extension of OBMP2 by introducing SOS into the
double-excitation amplitudes, denoted as the O2BMP2 method. We assess the
O2BMP2 performance on the triple-bond N2 dissociation curve, singlet-triplet
(ST) gaps of various sets of molecules, and ionization potentials (IPs)
obtained from the Koopmans’ approximation. We found that O2BMP2 can
dramatically outperform standard MP2 and reach the accuracy of coupled-cluster
methods in all cases considered in this work. Also, O2BMP2 performs better
than OBMP2 in most cases.

Details of OBMP2 theory are presented in Refs. 26, 27, 28, and it is
implemented in a local version of PySCF39. The OBMP2 Hamiltonian is derived
through the canonical transformation 29, 30, 31, 32, 33, 34:

$\displaystyle\hat{\bar{H}}=e^{\hat{A}^{\dagger}}\hat{H}e^{\hat{A}},$ (1)

with the molecular Hamiltonian as

$\displaystyle\hat{H}=\sum_{pq}h^{p}_{q}\hat{a}_{p}^{q}+\tfrac{1}{2}\sum_{pqrs}g^{pr}_{qs}\hat{a}_{pr}^{qs}.$
(2)

Here, $\left\\{p,q,r,\ldots\right\\}$ indices referring to general ($all$)
spin orbitals. One- and two-body second-quantized operators $\hat{a}_{p}^{q}$
and $\hat{a}_{pq}^{rs}$ are given by
$\hat{a}_{p}^{q}=\hat{a}^{\dagger}_{p}\hat{a}_{q}$ and
$\hat{a}_{pq}^{rs}=\hat{a}^{\dagger}_{p}\hat{a}^{\dagger}_{q}\hat{a}_{s}\hat{a}_{r}$.
$h_{pq}$ and $v_{pq}^{rs}$ are one- and two-electron integrals, respectively.
In OBMP2, the anti-Hermitian excited operator $\hat{A}$ includes only double
excitations.

$\displaystyle\hat{A}=\hat{A}_{\text{D}}=\tfrac{1}{2}\sum_{ij}^{occ}\sum_{ab}^{vir}T_{ij}^{ab}(\hat{a}_{ab}^{ij}-\hat{a}_{ij}^{ab})\,,$
(3)

with the MP2 amplitude

$\displaystyle
T_{ij}^{ab}=\frac{g_{ij}^{ab}}{\epsilon_{i}+\epsilon_{j}-\epsilon_{a}-\epsilon_{b}}\,,$
(4)

where $\left\\{i,j,k,\ldots\right\\}$ indices refer to occupied ($occ$) spin
orbitals and $\left\\{a,b,c,\ldots\right\\}$ indices refer to virtual ($vir$)
spin orbitals; $\epsilon_{i}$ is the orbital energy of the spin-orbital $i$.
The OBMP2 Hamiltonian is defined as

$\displaystyle\hat{H}_{\text{OBMP2}}=\hat{H}_{\text{HF}}+\left[\hat{H},\hat{A}_{\text{D}}\right]_{1}+\tfrac{1}{2}\left[\left[\hat{F},\hat{A}_{\text{D}}\right],\hat{A}_{\text{D}}\right]_{1}=\,\,\hat{H}_{\text{HF}}+\hat{v}_{\text{OBMP2}}.$
(5)

In Eq. 5, commutators with the subscription 1, $[\ldots]_{1}$, involve one-
body operators and constants that are reduced from many-body operators using
the cumulant approximation35, 36, 37, 38. $\hat{H}_{\text{HF}}$ is standard HF
Hamiltonian and $\hat{v}_{\text{OBMP2}}$ is a correlated potential composing
of one-body operators with the working expression given by

$\displaystyle\hat{v}_{\text{OBMP2}}=$
$\displaystyle\overline{T}_{ij}^{ab}\left[f_{a}^{i}\,\hat{\Omega}\left(\hat{a}_{j}^{b}\right)+g_{ab}^{ip}\,\hat{\Omega}\left(\hat{a}_{j}^{p}\right)-g^{aq}_{ij}\,\hat{\Omega}\left(\hat{a}^{b}_{q}\right)\right]$
$\displaystyle-2\overline{T}_{ij}^{ab}g^{ij}_{ab}+\,f_{a}^{i}\overline{T}_{ij}^{ab}\overline{T}_{jk}^{bc}\,\hat{\Omega}\left(\hat{a}_{c}^{k}\right)$
$\displaystyle+f_{c}^{a}T_{ij}^{ab}\overline{T}_{il}^{cb}\,\hat{\Omega}\left(\hat{a}^{l}_{j}\right)+f_{c}^{a}T_{ij}^{ab}\overline{T}_{kj}^{cb}\,\hat{\Omega}\left(\hat{a}^{k}_{i}\right)$
$\displaystyle-f^{k}_{i}T_{ij}^{ab}\overline{T}_{kl}^{ab}\,\hat{\Omega}\left(\hat{a}_{l}^{j}\right)-f^{p}_{i}T_{ij}^{ab}\overline{T}_{kj}^{ab}\,\hat{\Omega}\left(\hat{a}^{p}_{k}\right)$
$\displaystyle+f^{k}_{i}T_{ij}^{ab}\overline{T}_{kj}^{ad}\,\hat{\Omega}\left(\hat{a}_{b}^{d}\right)+f_{k}^{i}T_{ij}^{ab}\overline{T}_{kj}^{cb}\,\hat{\Omega}\left(\hat{a}_{a}^{c}\right)$
$\displaystyle-
f_{c}^{a}T_{ij}^{ab}\overline{T}_{ij}^{cd}\,\hat{\Omega}\left(\hat{a}^{b}_{d}\right)\,-f_{p}^{a}T_{ij}^{ab}\overline{T}_{ij}^{cb}\,\hat{\Omega}\left(\hat{a}^{p}_{c}\right)$
$\displaystyle-2f_{a}^{c}{T}_{ij}^{ab}\overline{T}_{ij}^{cb}+2f_{i}^{k}{T}_{ij}^{ab}\overline{T}_{kj}^{ab}.$
(6)

with $\overline{T}_{ij}^{ab}={T}_{ij}^{ab}-{T}_{ji}^{ab}$, the symmetrization
operator
$\hat{\Omega}\left(\hat{a}^{p}_{q}\right)=\hat{a}^{p}_{q}+\hat{a}^{q}_{p}$,
and the Fock matrix

$\displaystyle
f_{p}^{q}=h_{p}^{q}+\sum_{i}^{occ}\left(g^{pi}_{qi}-g^{pi}_{iq}\right).$ (7)

We rewrite $\hat{H}_{\text{OBMP2}}$ (Eq.5) in a similar form to standard HF as
follows:

$\displaystyle\hat{H}_{\text{OBMP2}}=$ $\displaystyle\hat{\bar{F}}+\bar{C},$
(8)

where the constant $\bar{C}$ is a sum of terms without excitation operators.
$\hat{\bar{F}}$ is the correlated Fock operator,
$\hat{\bar{F}}=\bar{f}^{p}_{q}\hat{a}_{p}^{q}$, with correlated Fock matrix
$\bar{f}^{p}_{q}$ written as

$\displaystyle\bar{f}^{p}_{q}$ $\displaystyle=f^{p}_{q}+v^{p}_{q}.$ (9)

$v^{p}_{q}$ serves as the correlation potential altering the uncorrelated HF
picture. The MO coefficients and energies then correspond to eigenvectors and
eigenvalues of $\bar{f}^{p}_{q}$.

Grimme40 found that the MP2 performance can be dramatically improved by
separating and scaling same-spin (SS) and opposite-spin (OS) contributions to
the correlation energy. Later, Jung et al. 41 extended Grimme’s method by only
considering the opposite-scaling component, SOS-MP2. Lochan and Head-Gordon 42
further developed the optimized second-order opposite-spin (O2) method by
optimizing orbitals with the SOS-MP2 energy. Kossmann and Neese 43 introduced
spin-component scaling to the OO-MP2 method by scaling the SS and OS
contributions to the MP2 amplitude. All these studies showed that SOS can
significantly improve the performance of conventional counterparts. In the
present work, we extend the OBMP2 method by incorporating the spin-opposite
scaling $c_{\text{os}}$ into the double-excitation amplitude (Eq. 4)

$\displaystyle
T_{ij}^{ab}=c_{\text{os}}\frac{g_{ij}^{ab}}{\epsilon_{i}+\epsilon_{j}-\epsilon_{a}-\epsilon_{b}}\,.$
(10)

The optimal value of $c_{\text{os}}$ for SOS-MP2 were found to be 1.3 42. In
the current work, we will assess three values $c_{\text{os}}=1.1,1.2,\text{and
}1.3$ to find the best scaling for O2BMP2.

In Ref. 28, we showed that the self-consistency of OBMP2 helps it avoid the
divergence in energy curves present in standard MP2 for H2 and LiH. We now
consider a more challenging system, N2. We use NEVPT2 with an active space of
(8e,8o) as the reference. Energies relative to the equilibrium energies of
each method are presented in Figure 1. Unsurprisingly, standard MP2 quickly
breaks down, whereas OBMP2 yields a better energy curve. However, beyond the
equilibrium bond length, the OBMP2 curve is far below the NEVPT2 reference.
O2BMP2, with all the scaling factors considered here, can improve the energy
curve upon OBMP2 and make curves close to NEVPT2. Among these factors, 1.2 may
perform best, particularly at long distances.

Figure 1: Potential energy curves N2 molecule in cc-pVDZ from different
methods. Energies are relative to the energy at the equilibrium geometry.

Let us now assess the performance of our methods on the prediction of singlet-
triplet (ST) gaps. We start with a test set including 38 small molecules. We
first examine the spin contamination presented in Figure S1 in Supporting
Information (SI). In the upper panel, we present some molecules for that HF
severely suffers from spin contamination. We can see that while MP2 cannot
eliminate the spin contamination in these cases, OBMP2 yields negligible spin
contamination. In the lower panel, we plot the change in spin contamination
with respect to OBMP2 iterations for two molecules CO and CO2. The spin
contamination at the first iteration is large and significantly reduced when
the loop converges, implying the importance of self-consistency in eliminating
the spin contamination. In Figure 2, we plot mean absolute deviations (MADs)
relative to CCSD(T) reference of ST gaps from different methods, including
MP2, SOS-MP2 with $c_{\text{OS}}$ = 1.2, OBMP2, and O2BMP2 with varying values
of $c_{\text{OS}}$. ST gaps of different methods are given in SI. We can see
that MP2 and SOS-MP2 give MADs larger than 0.3 eV, whereas OBMP2 and O2BMP2
with three scaling factors yield MADs smaller than 0.15 eV and comparable to
CCSD. For this set of small molecules, O2BMP2 is only marginally better than
OBMP2.

Figure 2: Mean absolute deviation (MAD) relative to CCSD(T) of singlet-triplet
(ST) gaps of 38 small molecules calculated using different methods. The basis
set cc-pVQZ was used.

We now consider some medium-size organic radicals adopted from Ref.44. All
results are shown in Table 1. We compare our results to experimental values
and other calculated results, including MP2, SOS-MP2($c_{\text{OS}}$ = 1.2),
CC2, and CCSD. MP2 and SOS-MP2 yield significant errors relative to the
experiment. While OBMP2 can dramatically improve MP2 ST gaps, its errors are
still quite large. Interestingly, O2BMP2 with $c_{\text{os}}=1.2$ performs
better than OBMP2 with a smaller MAD (0.19 eV).

Table 1:  Singlet-triplet gaps (in eV) of biradicals. CCSD, CC2, and experimental ST gaps are taken from Ref. 44. The basis set cc-pVTZ was employed. Molecules | exp | CCSD | CC2 | MP2 | SOS-MP2/1.2 | OBMP2 | O2BMP2/1.2
---|---|---|---|---|---|---|---
ethene | 4.36 | 4.42 | 4.52 | 4.59 | 4.55 | 4.60 | 4.35
butadiene | 3.22 | 3.25 | 3.34 | 3.52 | 3.55 | 3.44 | 3.37
hexatriene | 2.61 | 2.62 | 2.78 | 3.54 | 3.51 | 2.86 | 2.82
octatetraene | 2.1 | 2.23 | 2.4 | 3.07 | 3.06 | 2.46 | 2.45
cyclopropene | 4.16 | 4.3 | 4.44 | 4.52 | 4.49 | 4.45 | 4.24
cyclopentadiene | 3.1 | 3.18 | 3.36 | 3.51 | 3.46 | 3.41 | 3.26
furan | 4.02 | 4.17 | 4.30 | 4.51 | 4.33 | 4.42 | 4.11
pyrrole | 4.21 | 4.52 | 4.68 | 4.88 | 4.66 | 4.76 | 4.44
tetrazine | 1.69 | 1.99 | 1.86 | 2.10 | 2.63 | 1.52 | 2.10
MAD |  | 0.11 | 0.25 | 0.53 | 0.53 | 0.31 | 0.19
MAX |  | 0.31 | 0.47 | 0.97 | 0.96 | 0.55 | 0.41

The next set consists of 10 aryl carbenes adopted from Ref. 45. Determining
the ST gap of carbenes is a difficult task for both experiment and theory.
Among classes of carbenes, aryl carbenes have attracted extensive attention
due to the accessibility of the triplet state. It has been evident that HF
theory fails to accurately reproduce ST gaps of carbenes, whereas DFT cannot
guarantee consistently accurate predictions. One of the reasons for the
failure of HF and DFT in the ST gap prediction of carbenes may be the large
spin contamination. As shown in Figure S2, both HF and MP2 severely suffer
from spin contamination. Thanks to the self-consistency, OBMP2 can
significantly reduce the spin contamination. To further see the importance of
self-consistency, we plot in Figure S3 spin densities of three aryl carbenes
1, 4, and 9. We can see that MP2 predicts spin densities spreading over whole
molecules, which may lead to large spin contamination. On the other hand,
OBMP2 predicts spin densities localizing on the aryl group, which is
consistent with CCSD prediction. The ST gaps of aryl carbenes predicted by
MP2, OBMP2, CCSD, and CCSD(T) are presented in Figure 3. We use the scaling
$c_{\text{os}=1.2}$ for O2BMP2. Unsurprisingly, CCSD results are close to the
CCSD(T) reference. On the other hand, while HF underestimates the ST gaps, MP2
significantly overestimates them. Our methods yield results very close to
higher-cost methods, CCSD and CCSD(T).

Figure 3: ST gaps of 10 aryl carbenes from different methods. The basis set
cc-pVDZ was used.

The last set we used to test the OBMP2 and O2BMP2 prediction of ST gap is
polyaromatic hydrocarbons (PAHs). The prediction of accurate ST gaps of
polyaromatic hydrocarbons has been challenging for computational methods46,
47, 48, 49, 50. While the ST gaps of linear PAHs have shown an exponential
decay with system size, those of the non-linear PAHs are marginally sensitive
to system size50. Unfortunately, the latter has not been observed by single-
reference methods like DFT 51. Dey and Gosh have attributed the failure of DFT
to the multi-reference nature of each state of non-linear PAHs 50. In the
current work, we consider polyacenes (linear PAHs) and helicene (non-linear
PAHs) with geometries taken from Ref. 50. ST gaps from OBMP2 and MP2 compared
to the density matrix renormalization group (DMRG) are shown in Figure 4. For
both cases, while MP2 errors relative to DMRG are significant, OBMP2 and
O2BMP2/1.2 can dramatically improve ST gaps of PAHs. Our methods predict ST
gaps close to DMRG for polyacenes, whereas their errors are still quite
significant for helicene. It could be because of the stronger multi-reference
nature present in helicene. It is worth stressing that OBMP2 and O2BMP2 can
reproduce DMRG prediction on the less dependence of ST gaps on the system size
for helicene, which has not been observed by single-reference methods like
DFT51, 50. In Figure 5, we plot spin densities of helicene[3] and helicene[4].
While MP2 spin densities are delocalized over the structures, OBMP2 ones are
localized along the preferentially stable double bonds, entirely consistent
with DMRG prediction 50.

Figure 4: ST gaps of linear polyacenes (left) and helicene (right) from
different methods. The DMRG references are taken from Ref. 50. The basis set
cc-pVDZ was used. Figure 5: Spin densities of the triplet state of helicene[3]
and helicene[4] from MP2 and OBMP2.

We now move to assess the performance of our methods on the prediction of
molecular IPs. Other previous studies showed that Koopmans’ approximation with
MP2 cannot give satisfactory accuracy in the prediction of IPs52, 53. It is
interesting to check whether our methods can achieve accurate IPs in the
framework of Koopmans’ approximation. We previously derived the formula of
OBMP2 IPs within Koopmans’ approximation in Ref. 25. In the current work, we
implement it in the spin-unrestricted OBMP2 version, removing only one
electron instead of a pair of electrons in the restricted version.

We first consider a test set of 21 small molecules with 58 valence IPs. We
report IP-EOM-CCSD and G0W0 with HF and DFT (PBE) references for comparison.
For O2BMP2, we have tested three scaling factors $c_{\text{os}}=1.1,1.2,$ and
$1.3$. Calculated and experimental IPs are given in SI. We show in Figure 6
mean absolute deviations (MAD) and maximum absolute deviations (MAX) relative
to experimental values. We can see that both G0W0 yield large MAD and MAX,
whereas IP-EOM-CCSD can significantly reduce MAD to 0.23 eV. Although OBMP2
performs better than G0W0, its errors are still large. Regarding O2BMP2,
unlike ST gaps, IPs are sensitive to the scaling factor $c_{\text{os}}$. Among
the three values, 1.1 gives the smallest errors with MDA comparable to that of
IP-EOM-CCSD and MAX even smaller than IP-EOM-CCSD.

Figure 6: Mean absolute deviation (MAD) and maximum absolute deviation (MAX)
relative to experimental values of valence IPs of small molecules from
different methods. The basis set cc-pVQZ was used.

We finally evaluate the IPs of 10 organic acceptor molecules with medium size
adopted from Ref. 54. The above assessment for small molecules shows that
O2BMP2/1.1 performs best. We thus report its results in comparison with IP-
EOM-CCSD and G0W0@HF. All results are summarized in Table 2. G0W0@HF vastly
overestimates IPs of acceptor molecules with MAD up to 0.46 eV, consistent
with the error found in Ref. 54. IP-EOM-CCSD yields results close to
experimental values with MAD of 0.17 eV. Surprisingly, O2BMP2/1.1 can reach an
accuracy similar to EOM-CCSD with MAD of 0.16 eV. The maximum error of
O2BMP2/1.1 is 0.3 eV for mDCNB and benzonitrile molecules that may have a
strong multi-reference nature.

Table 2:  First ionization potential of 10 organic acceptor molecules. Experimental values are adopted from Ref. 54. The basis set aug-cc-pVDZ was employed. Molecules | exp | IP-EOM-CCSD | G0W0@HF | O2BMP2/1.1
---|---|---|---|---
Bezonquinone (BQ) | 9.95 | 10.04 | 10.30 | 10.10
tetrafluoro-BQ | 10.70 | 11.05 | 11.40 | 10.84
tetrachloro-BQ | 9.74 | 10.10 | 10.40 | 9.88
fumaronitrie | 11.15 | 11.30 | 11.38 | 11.10
maleic anhydride | 11.07 | 11.02 | 12.15 | 10.96
mDCNB | 10.20 | 10.26 | 10.77 | 9.90
nitrobenzene | 9.86 | 10.06 | 10.10 | 9.68
phtalic anhydride | 10.10 | 10.40 | 10.50 | 10.03
TCNE | 11.77 | 11.91 | 12.03 | 11.65
benzonitrile | 9.70 | 9.73 | 9.78 | 9.40
MAD |  | 0.17 | 0.46 | 0.16
MAX |  | 0.36 | 1.08 | 0.30

In summary, we have extended our recently developed method, OBMP2, by
introducing the spin-opposite scaling to the double-excitation amplitudes,
termed O2BMP2. We assess the O2BMP2 performance on the triple-bond N2
dissociation, ST gaps, and IPs of medium-size organic compounds. O2BMP2
performs much better than standard MP2 and reaches the accuracy of coupled-
cluster methods in all cases considered in this work. Our method is then
expected to help tackle realistic, challenging systems with large sizes.
Working on further reducing computational costs of OBMP2 and O2BMP2 is in
progress.

## References

  * Møller and Plesset 1934 Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. _Phys. Rev._ 1934, _46_ , 618
  * Byrd et al. 2001 Byrd, E. F.; Sherrill, C. D.; Head-Gordon, M. The theoretical prediction of molecular radical species: a systematic study of equilibrium geometries and harmonic vibrational frequencies. _J. Phys. Chem. A_ 2001, _105_ , 9736
  * Stück et al. 2011 Stück, D.; Baker, T. A.; Zimmerman, P.; Kurlancheek, W.; Head-Gordon, M. On the nature of electron correlation in C60. _J. Chem. Phys._ 2011, _135_ , 11B608
  * Lochan and Head-Gordon 2007 Lochan, R. C.; Head-Gordon, M. Orbital-optimized opposite-spin scaled second-order correlation: An economical method to improve the description of open-shell molecules. _J. Chem. Phys._ 2007, _126_ , 164101
  * Neese et al. 2009 Neese, F.; Schwabe, T.; Kossmann, S.; Schirmer, B.; Grimme, S. Assessment of orbital-optimized, spin-component scaled second-order many-body perturbation theory for thermochemistry and kinetics. _J. Chem. Theory Comput._ 2009, _5_ , 3060
  * Bozkaya et al. 2011 Bozkaya, U.; Turney, J. M.; Yamaguchi, Y.; Schaefer III, H. F.; Sherrill, C. D. Quadratically convergent algorithm for orbital optimization in the orbital-optimized coupled-cluster doubles method and in orbital-optimized second-order Møller-Plesset perturbation theory. _J. Chem. Phys._ 2011, _135_ , 104103
  * Bozkaya and Sherrill 2013 Bozkaya, U.; Sherrill, C. D. Analytic energy gradients for the orbital-optimized second-order Møller–Plesset perturbation theory. _J. Chem. Phys._ 2013, _138_ , 184103
  * Bozkaya 2014 Bozkaya, U. Orbital-optimized second-order perturbation theory with density-fitting and cholesky decomposition approximations: An efficient implementation. _J. Chem. Theory Comput._ 2014, _10_ , 2371
  * Bozkaya 2014 Bozkaya, U. Analytic energy gradients and spin multiplicities for orbital-optimized second-order perturbation theory with density-fitting approximation: an efficient implementation. _J. Chem. Theory Comput._ 2014, _10_ , 4389–4399
  * Goerigk and Grimme 2014 Goerigk, L.; Grimme, S. Double-hybrid density functionals. _WIREs: Comput. Mol. Sci._ 2014, _4_ , 576
  * Martin and Santra 2020 Martin, J. M.; Santra, G. Empirical double-hybrid density functional theory: A “third way”in between WFT and DFT. _Isr. J. Chem_ 2020, _60_ , 787
  * Stück and Head-Gordon 2013 Stück, D.; Head-Gordon, M. Regularized orbital-optimized second-order perturbation theory. _J. Chem. Phys._ 2013, _139_ , 244109
  * Razban et al. 2017 Razban, R. M.; Stück, D.; Head-Gordon, M. Addressing first derivative discontinuities in orbital-optimised opposite-spin scaled second-order perturbation theory with regularisation. _Mol. Phys._ 2017, _115_ , 2102
  * Lee and Head-Gordon 2018 Lee, J.; Head-Gordon, M. Regularized orbital-optimized second-order Møller–Plesset perturbation theory: A reliable fifth-order-scaling electron correlation model with orbital energy dependent regularizers. _J. Chem. Theory Comput._ 2018, _14_ , 5203
  * Shee et al. 2021 Shee, J.; Loipersberger, M.; Rettig, A.; Lee, J.; Head-Gordon, M. Regularized second-order Møller–Plesset theory: A more accurate alternative to conventional MP2 for noncovalent interactions and transition metal thermochemistry for the same computational cost. _J. Phys. Chem. Lett._ 2021, _12_ , 12084
  * Rettig et al. 2022 Rettig, A.; Shee, J.; Lee, J.; Head-Gordon, M. Revisiting the orbital energy-dependent regularization of orbital-optimized second-order Møller–Plesset theory. _J. Chem. Theory Comput._ 2022, _18_ , 5382–5392
  * Keller et al. 2022 Keller, E.; Tsatsoulis, T.; Reuter, K.; Margraf, J. T. Regularized second-order correlation methods for extended systems. _J. Chem. Phys._ 2022, _156_
  * Santra and Martin 2022 Santra, G.; Martin, J. M. Do double-hybrid functionals benefit from regularization in the PT2 term? Observations from an extensive benchmark. _J. Phys. Chem. Lett._ 2022, _13_ , 3499–3506
  * Daas et al. 2023 Daas, K. J.; Kooi, D. P.; Peters, N. C.; Fabiano, E.; Della Sala, F.; Gori-Giorgi, P.; Vuckovic, S. Regularized and Opposite Spin-Scaled Functionals from Møller–Plesset Adiabatic Connection-Higher Accuracy at Lower Cost. _J. Phys. Chem. Lett._ 2023, _14_ , 8448–8459
  * Carter-Fenk and Head-Gordon 2023 Carter-Fenk, K.; Head-Gordon, M. Repartitioned Brillouin-Wigner perturbation theory with a size-consistent second-order correlation energy. _J. Chem. Phys._ 2023, _158_
  * Carter-Fenk et al. 2023 Carter-Fenk, K.; Shee, J.; Head-Gordon, M. Optimizing the Regularization in Size-Consistent Second-Order Brillouin-Wigner Perturbation Theory. _arXiv preprint arXiv:2309.01376_ 2023,
  * Behnle and Fink 2019 Behnle, S.; Fink, R. F. REMP: A hybrid perturbation theory providing improved electronic wavefunctions and properties. _J. Chem. Phys._ 2019, _150_
  * Behnle and Fink 2021 Behnle, S.; Fink, R. F. OO-REMP: Approaching chemical accuracy with second-order perturbation theory. _J. Chem. Theory Comput._ 2021, _17_ , 3259–3266
  * Grimme et al. 2012 Grimme, S.; Goerigk, L.; Fink, R. F. Spin-component-scaled electron correlation methods. _WIREs: Comput. Mol. Sci._ 2012, _2_ , 886
  * Tran and Yanai 2013 Tran, L. N.; Yanai, T. Correlated one-body potential from second-order Møller-Plesset perturbation theory: Alternative to orbital-optimized MP2 method. _J. Chem. Phys._ 2013, _138_ , 224108
  * Tran 2021 Tran, L. N. Improving perturbation theory for open-shell molecules via self-consistency. _J. Phys. Chem. A_ 2021, _125_ , 9242
  * Tran 2022 Tran, L. N. Can second-order perturbation theory accurately predict electron density of open-shell molecules? The importance of self-consistency. _Phys. Chem. Chem. Phys._ 2022, _24_ , 19393
  * Le and Tran 2023 Le, N. T.; Tran, L. N. Correlated reference-assisted variational quantum eigensolver. _J. Phys. Chem. A_ 2023, _127_ , 5222
  * Yanai and Chan 2006 Yanai, T.; Chan, G. K.-L. Canonical transformation theory for multireference problems. _J. Chem. Phys._ 2006, _124_ , 194106
  * Yanai and Chan 2007 Yanai, T.; Chan, G. K.-L. Canonical transformation theory from extended normal ordering. _J. Chem. Phys._ 2007, _127_ , 104107
  * Chan and Yanai 2007 Chan, G. K.-L.; Yanai, T. In _Adv. Chem. Phys., Volume 134_ , reduced-de ed.; Mazziotti, D. A., Ed.; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2007; Chapter 13, p 343
  * Neuscamman et al. 2009 Neuscamman, E.; Yanai, T.; Chan, G. K.-L. Quadratic canonical transformation theory and higher order density matrices. _J. Chem. Phys._ 2009, _130_ , 124102
  * Neuscamman et al. 2010 Neuscamman, E.; Yanai, T.; Chan, G. K.-L. Strongly contracted canonical transformation theory. _J. Chem. Phys._ 2010, _132_ , 024106
  * Neuscamman et al. 2010 Neuscamman, E.; Yanai, T.; Chan, G. K.-L. A review of canonical transformation theory. _Int. Rev. Phys. Chem._ 2010, _29_ , 231
  * Kutzelnigg and Mukherjee 1997 Kutzelnigg, W.; Mukherjee, D. Normal order and extended Wick theorem for a multiconfiguration reference wave function. _J. Chem. Phys._ 1997, _107_ , 432
  * Mazziotti 1998 Mazziotti, D. A. Contracted Schrödinger equation: Determining quantum energies and two-particle density matrices without wave functions. _Phys. Rev. A_ 1998, _57_ , 4219
  * Mazziotti 1998 Mazziotti, D. A. Approximate solution for electron correlation through the use of Schwinger probes. _Chem. Phys. Lett._ 1998, _289_ , 419
  * Kutzelnigg and Mukherjee 1999 Kutzelnigg, W.; Mukherjee, D. Cumulant expansion of the reduced density matrices. _J. Chem. Phys._ 1999, _110_ , 2800
  * Sun et al. 2018 Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S. et al. PySCF: the Python-based simulations of chemistry framework. _WIREs: Comput. Mol. Sci._ 2018, _8_ , e1340
  * Grimme 2003 Grimme, S. Improved second-order Møller–Plesset perturbation theory by separate scaling of parallel-and antiparallel-spin pair correlation energies. _J. Chem. Phys._ 2003, _118_ , 9095
  * Jung et al. 2004 Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled opposite-spin second order Møller–Plesset correlation energy: An economical electronic structure method. _J. Chem. Phys._ 2004, _121_ , 9793
  * Lochan and Head-Gordon 2007 Lochan, R. C.; Head-Gordon, M. Orbital-optimized opposite-spin scaled second-order correlation: An economical method to improve the description of open-shell molecules. _J. Chem. Phys._ 2007, _126_
  * Kossmann and Neese 2010 Kossmann, S.; Neese, F. Correlated ab initio spin densities for larger molecules: orbital-optimized spin-component-scaled MP2 method. _J. Phys. Chem. A_ 2010, _114_ , 11768
  * Schreiber et al. 2008 Schreiber, M.; Silva-Junior, M. R.; Sauer, S.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. _J. Chem. Phys._ 2008, _128_
  * Ghafarian Shirazi et al. 2018 Ghafarian Shirazi, R.; Neese, F.; Pantazis, D. A. Accurate Spin-State Energetics for Aryl Carbenes. _J. Chem. Theory Comput._ 2018, _14_ , 4733
  * Hachmann et al. 2007 Hachmann, J.; Dorando, J. J.; Avilés, M.; Chan, G. K. The radical character of the acenes: A density matrix renormalization group study. _J. Chem. Phys._ 2007, _127_
  * Ibeji and Ghosh 2015 Ibeji, C. U.; Ghosh, D. Singlet–triplet gaps in polyacenes: a delicate balance between dynamic and static correlations investigated by spin–flip methods. _Phys. Chem. Chem. Phys._ 2015, _17_ , 9849
  * Sharma et al. 2019 Sharma, P.; Bernales, V.; Knecht, S.; Truhlar, D. G.; Gagliardi, L. Density matrix renormalization group pair-density functional theory (DMRG-PDFT): singlet–triplet gaps in polyacenes and polyacetylenes. _Chem. Sci._ 2019, _10_ , 1716
  * Shee et al. 2019 Shee, J.; Arthur, E. J.; Zhang, S.; Reichman, D. R.; Friesner, R. A. Singlet–triplet energy gaps of organic biradicals and polyacenes with auxiliary-field quantum Monte Carlo. _J. Chem. Theory Comput._ 2019, _15_ , 4924
  * Dey and Ghosh 2022 Dey, M.; Ghosh, D. Curious Case of Singlet Triplet Gaps in Nonlinear Polyaromatic Hydrocarbons. _J. Phys. Chem. Lett._ 2022, _13_ , 11795
  * Rulíšek et al. 2007 Rulíšek, L.; Exner, O.; Cwiklik, L.; Jungwirth, P.; Starỳ, I.; Pospíšil, L.; Havlas, Z. On the convergence of the physicochemical properties of [n] helicenes. _J. Phys. Chem. C_ 2007, _111_ , 14948
  * Ayala et al. 2001 Ayala, P. Y.; Kudin, K. N.; Scuseria, G. E. Atomic orbital Laplace-transformed second-order Møller–Plesset theory for periodic systems. _J. Chem. Phys._ 2001, _115_ , 9698
  * Maksić and Vianello 2002 Maksić, Z. B.; Vianello, R. How good is Koopmans’ approximation? G2 (MP2) study of the vertical and adiabatic ionization potentials of some small molecules. _J. Phys. Chem. A_ 2002, _106_ , 6515
  * Knight et al. 2016 Knight, J. W.; Wang, X.; Gallandi, L.; Dolgounitcheva, O.; Ren, X.; Ortiz, J. V.; Rinke, P.; Körzdörfer, T.; Marom, N. Accurate ionization potentials and electron affinities of acceptor molecules III: a benchmark of GW methods. _J. Chem. Theory Comput._ 2016, _12_ , 615

## Supporting Information to:
Reaching high accuracy for energetic properties at second-order perturbation
cost by merging self-consistency and spin-opposite scaling

Nhan Tri Tran
University of Science, Vietnam National University, Ho Chi Minh City, Vietnam
Hoang Thanh Nguyen
Institute of Applied Mechanics and Informatics, Vietnam Academy of Science and
Technology, Ho Chi Minh City, Vietnam
Lan Nguyen Tran*
Email<EMAIL_ADDRESS>
Department of Physics, International University, Ho Chi Minh City, Vietnam
Vietnam National University, Ho Chi Minh City, Vietnam

File supporting-information.xls includes:

  * •

Singlet-triplet gaps (in eV) of 39 small molecules

  * •

Valence IPs (in eV) of 21 small molecules.

Figure S1: Upper panel: $\left<S^{2}\right>$ errors for the triplet state of
small radicals from HF, MP2, and OBMP2. Lower panel: $\left<S^{2}\right>$
errors for CO and CO2 at different OBMP2 iterations. Figure S2:
$\left<S^{2}\right>$ errors for the triplet state of 10 aryl carbenes from HF,
MP2, and OBMP2 Figure S3: Spin densities of the triplet state of three aryl
carbenes 1, 4, and 9 from MP2, OBMP2, and CCSD.


# Unveiling extra dimensions through the shadow of a dark star

Sagnik Roy<EMAIL_ADDRESS>Soham Chatterjee<EMAIL_ADDRESS>Ratna
Koley<EMAIL_ADDRESS>Department of Physics, Presidency
University, 86/1 College Street, Kolkata 700073, India

###### Abstract

The shadow of a black hole or a collapsing star is of great importance as we
can extract important properties of the object and of the surrounding
spacetime from the shadow profile. It can also be used to distinguish
different types of black holes and ultra compact objects. In this work, we
have analytically calculated the shadow of a higher dimensional collapsing
dark star, described by higher dimensional Vaidya metric, by choosing a
slightly generalized version of Misner–-Sharp mass function. We have also
numerically investigated the properties of the shadows of the black holes and
the collapsing stars for a slightly more general mass function. Examining the
potential influence of extra spatial dimensions on the shadow, we have
explored the possibility of distinguishing higher dimensions from the standard
four-dimensional spacetime.

shadow, collapsing star, Vaidya spacetime, higher dimensions, photon sphere

###### pacs:

04.50.Gh, 04.40.Dg, 95.30.Sf, 98.80.Jk, 04.70.Bw

##  I Introduction

_Photon sphere_ plays a very important role in the study of light trajectories
around a blackhole. Light passing near a black hole gets deflected by the
strong gravitational field in such a way that it moves along a circular orbit
around the black hole and for the case of a spherically symmetric black hole,
these circular light trajectories form a sphere called the _photon sphere_.
Observation of a black hole against the backdrop of light sources reveals a
black disc – known as the _shadow_ of a black hole. For a Schwarzschild black
hole, the shadow is circular due to the spherical symmetry of the space time
and the light rays spiral in a orbit with radius $r=3m^{bh}$, where $m^{bh}$
is the mass parameter of the black hole. On the other hand, for a rotating
Kerr black hole the shadow gets deformed and it is flattened on one side. In
the minimal set up, an observer will see the shadow if the space around the
black hole is filled with light rays and there lies no light source between
the observer and the black hole. All past-oriented light rays that starts at
observer’s position can be categorised into two classes: (i) light rays
getting deflected by the black hole and meeting with a source on their way and
(ii) light rays going directly to the horizon. We thus assign brightness to
the first class of light rays and darkness to the second group. There is a
third category of light rays in between these two aforementioned classes and
they asymptotically spiral in a photon sphere.

In 2000, Falcke et. al. [1] proposed that a shadow can be observed practically
and in 2019, Event Horizon Telescope (EHT) [2, 3, 4, 5, 6] successfully
observed the gravitational light deflection by the supermassive black hole of
M87. This is a major breakthrough in this field, followed by the observation
of the shadow of supermassive black hole Sagittarius A* in the center of our
Milky Way in 2022 [7]. Testing of compact gravitational objects, arising out
of different theoretical propositions is being done through this new window of
observation.

Analytical calculations of shadow is generally done by considering a black
hole against the backdrop of light sources with the assumption that there is
no source between the black hole and the observer. The eternal black holes
have been the popular choice for these calculations so far. In 1966, Synge [8]
analytically calculated the angular radius of the shadow of a Schwarzschild
black hole and Bardeen [9] found the shape of the shadow of a Kerr black hole
for an observer at infinity by choosing eternal black holes. According to our
present knowledge, the black holes are formed by gravitational collapse and
evolve by accretion of surrounding matter or mergers with other black holes.
Thus we find it quite interesting to explore the shadow of evolving systems.
Visual appearance of a star collapsing through its gravitational radius was
first studied by Ames and Thorne [10], after that Jaffe [11], Lake and Roeder
[12] and Frolov et. al. [13] studied the frequency shift of light coming from
the surface of a collapsing non-transparent star. On the other hand, work by
Kong et. al. [14, 15] and Ortiz et. al. [16, 17] dealt with the frequency
shift of light passing through a collapsing transparent star, thereby
contrasting the collapse to a black hole with the collapse to a naked
singularity. General properties of photon sphere of a static, spherically
symmetric spacetime is discussed in [18]. A detailed review of analytical
calculations of black shadow is presented in [19]. Recently Schneider and
Perlick [20] studied the shadow of a dark and non-transparent collapsing star.
Also Solanki and Perlick [21] studied the shadow of a time-dependent black
hole described by Vaidya metric. The above-mentioned works are all in
$4$-dimensions. However, the shadows of different kinds of black holes in
higher dimensions have also been studied with great interest by several
physicists. Singh and Ghosh [22] obtained the shadow of a
Schwarzschild–Tangherlini black hole and Papnoi et. al. [23] studied the
shadow of $5$-dimensional Myres–Perry black hole. The shadows of other types
of higher dimensional black holes and the effect of extra dimensions on the
shadow have been explored in [24, 25, 26, 27, 28].

In this work, we consider a _higher dimensional collapsing dark star_ emitting
radiation and non-relativistic matter particles. In the situation when the
amount of radiation is too feeble to be detected, the star will appear as a
non-transparent, black disc against the background of light sources. In order
to explore the evolution of the black disc, we choose the generalized Vaidya
metric which is a non-vacuum solution of the Einstein’s field equation around
a spherically symmetric body and generalized it to the higher dimensions. The
photon sphere and shadow profile have been studied thereafter. In a recent
work by Banerjee et. al. [29] it has been argued that the observed shadow of
the M87 and Sgr A* may carry a signature of extra dimension through a tidal
charge as the hair. It was shown that the observed shadow and the image
diameter of the Sgr A* always predict a non-zero value of the tidal charge
parameter within the allowed range of angular diameter of the shadow. However,
a more complete study is due where one should consider the effect of the
environment on the shadow profile in the higher dimensional scenario. The
degeneracy in the existing explanations of the shadow profiles can only be
lifted by means of independent experimental results. That is why we have
adopted a different way to suggest an observational tool for the extra
dimensions. Obtaining the shadow of a collapsing star both analytically and
numerically, we have shown how the mass of the star can be used to speculate a
signature of the extra spatial dimension(s).

The work is organised in the following way. In Section II, we have discussed
the generalized $D$-dimensional Vaidya metric in Eddington–Finkelstein like
coordinates and the generalized mass function. A new conformal symmetry has
been identified through a coordinate transformation. The null geodesic has
been studied in this new coordinate for a special type of mass parameter.
Section III deals with the angular radius and the shape of a higher
dimensional evaporating black hole in outgoing Vaidya metric. We have used a
slightly generalized version of Misner–Sharp mass function. The escape angle
of a higher dimensional Vaidya black hole has been calculated in the frame of
a static observer. In Section IV, we first obtained the shadow of a collapsing
star for a static observer in the transformed frame and then reverted back to
the original coordinates (Eddington–Finkelstein-like coordinates). Numerical
analysis of the photon sphere and the shadow profile is shown in Section V.
Next we explored the detection possibilities of extra dimensions using the
shadow of a collapsing star. Finally concluded in Section VII.

##  II Null Geodesics in Higher dimensional Vaidya Spacetime

The matter fields present in the nature can be classified in general into two
broad categories – (i) with energy-momentum tensor having one timelike and
three spacelike eigenvectors (this includes dust, perfect fluid, called type-I
matter field) and (ii) with energy-momentum tensor having double null
eigenvectors (which includes radiation, null dust, called type-II matter
field). The most general spherically symmetric metric for any arbitrary
combination of these two types of matter fields in outgoing
Eddington–Finkelstein-like coordinates is given by [30, 31],

$ds^{2}=-e^{2\psi_{4}(u,r)}\left[1-\dfrac{2m_{4}(u,r)}{r}\right]du^{2}-2e^{\psi_{4}(u,r)}dudr+r^{2}d\Omega_{2}^{2}$
(1)

where $m_{4}(u,r)$ is the Misner–Sharp mass function in $4$-dimensions [31,
32]. This tells about the amount of energy inside radial distance $r$ at
retarded time $u$ and $d\Omega_{2}^{2}$ is the metric on a $2$-dimensional
unit sphere. Wang and Wu [33] proved that a specific combination of the two
types of matter fields mentioned above can give rise to a generalized Vaidya
metric as follows

$ds_{4}^{2}=-\left(1-\dfrac{2m_{4}(u,r)}{r}\right)du^{2}-2dudr+r^{2}d\Omega_{2}^{2}$
(2)

This is a straight forward generalisation of the original Vaidya metric [34].
The same metric was generalized for $D$-dimensions by Iyer and Vishveshwara in
[35] (further generalisation for charged blackholes can be found in [36]),

$ds_{D}^{2}=-\left(1-\dfrac{2m_{D}(u,r)}{r^{D-3}}\right)du^{2}-2dudr+r^{2}d\Omega_{D-2}^{2}$
(3)

where $d\Omega_{D-2}^{2}$ denotes the metric on a $(D-2)$-dimensional unit
sphere, defined by,

$d\Omega_{D-2}^{2}=d\theta_{1}^{2}+\displaystyle\sum_{i=2}^{D-2}\left(\displaystyle\prod_{j=1}^{i-1}\sin{{}^{2}\theta_{j}}\right)d\theta_{i}^{2}$
(4)

and $m_{D}(u,r)$ is the $D$-dimensional Misner–Sharp mass function [37]
defined as

$m_{D}(u,r)=\dfrac{4\pi^{-(D-3)/2}}{(D-2)}\Gamma\left(\dfrac{D-1}{2}\right)G_{D}M_{D}(u,r)~{}\propto~{}G_{D}M_{D}(u,r)$
(5)

where $G_{D}$ is Newton’s gravitational constant in $D$-dimensions and
$M_{D}(u,r)$ is the mass of the star for its radius $r$ at time $u$. Let us
consider that the star is collapsing by radiating shell of null dusts and
other matter particles. We take _outgoing_ Eddington–Finkelstein time in the
metric. Also throughout the paper, we consider a _soft_ version of the
Oppenheimer–Snyder model [38], i.e. a very small amount of outward pressure
resides at the surface of the dark star, so that a radial timelike inward
geodesic of a free particle can represent the collapse of the star surface. A
careful examination reveals there is no other conformal Killing vector field
except $\frac{\partial}{\partial\phi}$, where we have set $\theta_{D-2}=\phi$.
Therefore, it is _almost impossible_ (barring some exceptions) to do any
analytical calculation of the shadow profile due to lack of enough constants
of motion. In that case, one has to import some kind of numerical techniques
(which we will do later).

In the following we will show that analytically one can proceed with a special
form for $m_{D}(u,r)$. Ojako et. al. stated in [32] that there exists a
homothetic Killing vector field, for which the spacetime is self-similar and
the mass function for the same is as follows.

$m_{4}(u,r)=\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}u^{n}}{r^{n-1}}$ (6)

where $\mathbb{Z}$ represents the set of all integers on the real line. In the
very similar fashion the mass function will take the following form in
$D$-dimensions

$m_{D}(u,r)=\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}u^{n}}{r^{n+3-D}}$
(7)

Obviously, the values of $a_{n}$’s should be so chosen that the convergence of
the series is confirmed. It is trivial to notice if all $a_{n}$’s are zero
except for $a_{1}$, the mass function will correspond to constant mass
accretion or expel in $4$-dimensions. We now perform a coordinate
transformation to non–angular coordinates, $(u,r)\longrightarrow(T,R)$,
following [21], to identify a new conformal symmetry,

$u=r_{0}~{}e^{T/r_{0}}~{}~{}\text{and}~{}~{}r=R~{}e^{T/r_{0}}$ (8)

where $r_{0}$ is an arbitrary constant with the dimension of length. In the
transformed coordinates, the line element (3) becomes,

$ds_{D}^{2}=e^{2T/r_{0}}\left\\{-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)dT^{2}-2dTdR+R^{2}d\Omega_{D-2}^{2}\right\\}$
(9)

The Lagrangian density corresponding to a massless particle (on the equatorial
plane i.e. $\theta_{i}=\pi/2$ for $1\leq i\leq D-3$) in this set of new
coordinates is given by

$\mathcal{L}_{D}=e^{2T/r_{0}}\left\\{-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}+\dfrac{2R}{r_{0}}\right)\dot{T}^{2}-2\dot{T}\dot{R}+R^{2}\dot{\phi}^{2}\right\\}=0$
(10)

It is trivial to identify the conformally conserved angular momentum from the
above equation. (10):

$L=R^{2}e^{2T/r_{0}}\dot{\phi}$ (11)

In addition to this we also get $\frac{\partial}{\partial T}$ from Eq. (10),
as a conformal Killing vector field leading to a constant of motion
$\mathcal{E}$, along every null geodesic.

$\mathcal{E}=e^{2T/r_{0}}\left\\{-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\dot{T}-\dot{R}\right\\}$
(12)

It is also apparent from Eq. (10) that for radial light rays
$\left(\dot{\phi}=0\right)$ one can have following relations for $T$.

$T=T_{0}~{}~{}\text{(constant) ~{}
or}~{}~{}T=T_{0}-2\displaystyle\int\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{-1}dR$
(13)

where the first equation is valid for outgoing and the second stands for
incoming radial light rays. We now explore the geodesic motion of massless
particles in this background. For this purpose, we have to solve two coupled
equations, one appearing from the Lagrangian density in Eq. (10) and the other
from the combined version of the conservation laws in Eqs. (11) and (12):

$\displaystyle\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\left(\dfrac{dT}{d\phi}\right)^{2}+2~{}\dfrac{dT}{d\phi}\dfrac{dR}{d\phi}=R^{2}$
(14)
$\displaystyle\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\dfrac{dT}{d\phi}+\dfrac{dR}{d\phi}=\dfrac{\mathcal{E}R^{2}}{L}$
(15)

Solution of the above set of equations (14) and (15), will provide us with the
motion of a test photon on the equatorial plane in this spacetime. Simplifying
the above equations we obtain

$\displaystyle\dfrac{dT}{d\phi}$ $\displaystyle=$
$\displaystyle\dfrac{-\dfrac{\mathcal{E}R^{2}}{L}\mp\sqrt{\dfrac{\mathcal{E}^{2}R^{4}}{L^{2}}+\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n-2}}+\dfrac{2R^{3}}{r_{0}}-R^{2}}}{1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}}$
(16) $\displaystyle\dfrac{dR}{d\phi}$ $\displaystyle=$
$\displaystyle\pm\sqrt{\dfrac{\mathcal{E}^{2}R^{4}}{L^{2}}+\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n-2}}+\dfrac{2R^{3}}{r_{0}}-R^{2}}$
(17)

Note that the above equations will give $T$ and $R$ as function of $\phi$. We
do not find the geodesics explicitly. The expression in Eq. (17) has been used
in Section III to calculate the shadow of an evaporating black hole as well as
the shadow of a collapsing dark star in Section IV.1.

##  III Shadow Profile of a Higher-dimensional Black Hole described by Vaidya
Metric

Let us consider a higher dimensional evaporating black hole. In four
dimensions an evaporating black hole can be – very idealistically – modelled
by a $4D$ Vaidya metric with decreasing mass function [39, 40, 41, 42]. We
assume this model remains valid in higher dimensions also. The shadow profile
of a $4$D Vaidya spacetime was studied by Solanki and Perlick [21] for a
linearly increasing and decreasing Misner–Sharp mass function. They have shown
for a black hole with linearly increasing mass there exist two horizons. This
is in contrast to the linearly decreasing mass, which admits only one horizon.
The angular radius of an evaporating Vaidya black hole shadow is shown to be
time-independent for a conformally static observer though the area of the
photon sphere is decreasing. On the contrary to this work, we have generalized
the functional dependence of the mass function slightly as given in Eq. (7).

In this section, we find the expression for the photon sphere and angular
radius of a Vaidya black hole as seen by a static observer in $(T,R)$ frame.
Let us identify the tetrad basis for the metric in Eq. (3),

$\displaystyle e_{T}$ $\displaystyle=$ $\displaystyle
e^{-T/r_{0}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{-1/2}\dfrac{\partial}{\partial
T}$ (18) $\displaystyle e_{R}$ $\displaystyle=$ $\displaystyle
e^{-T/r_{0}}\left[\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{-1/2}\dfrac{\partial}{\partial
T}~{}-~{}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{1/2}\dfrac{\partial}{\partial
R}\right]$ (19) $\displaystyle e_{\phi}$ $\displaystyle=$
$\displaystyle\dfrac{e^{-T/r_{0}}}{R}\dfrac{\partial}{\partial\phi}$ (20)

We have set $\theta_{D-2}=\phi$, the azimuthal angle of the coordinate system
and have restricted ourselves on the equatorial plane (i.e. $\theta_{i}=\pi/2$
for $1\leq i\leq D-3$) of the collapsing star. This can always be done in case
of spherical symmetry without any loss of generality.

We now consider a null geodesic $(T(s),R(s),\phi(s))$ on the equatorial plane,
where $s$ is the affine parameter along the geodesic of a massless particle.
Let us expand the tangent vector with respect to the tetrad (18). Since the
tangent vector is light-like, the expansion can be written in terms of the
angle $(\alpha)$ between the photon geodesic and the radial direction in the
rest system of the observer in the following manner.

$\dot{T}~{}\dfrac{\partial}{\partial T}+\dot{R}~{}\dfrac{\partial}{\partial
R}+\dot{\phi}~{}\dfrac{\partial}{\partial\phi}=\kappa\left(e_{T}+e_{R}\cos{\alpha}+e_{\phi}\sin{\alpha}\right)$
(21)

where $\kappa$ is a scale factor $(\kappa>0)$, $\alpha$ denotes the celestial
coordinate and $\dot{T}$ represents $\frac{dT}{ds}$ and so on. Comparison of
the coefficients of $\frac{\partial}{\partial T}$, $\frac{\partial}{\partial
R}$ and $\frac{\partial}{\partial\phi}$ on both sides of Eq. (21) results in

$\displaystyle\dot{T}$ $\displaystyle=$
$\displaystyle\phantom{+}\kappa~{}e^{-T/r_{0}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{-1/2}\left(1+\cos{\alpha}\right)~{},$
(22) $\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle-\kappa~{}e^{-T/r_{0}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{1/2}\cos{\alpha}~{}~{}\text{and}$
(23) $\displaystyle\dot{\phi}$ $\displaystyle=$
$\displaystyle\phantom{+}\dfrac{\kappa~{}e^{-T/r_{0}}}{R}~{}\sin{\alpha}$ (24)

From the last two equations above we can write

$\left(\dfrac{dR}{d\phi}\right)^{2}=R^{2}\cot{{}^{2}\alpha}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)$
(25)

Further comparing Eqs. (17) and (25) we obtain

$\sin{{}^{2}\alpha}=\dfrac{L^{2}}{\mathcal{E}^{2}R^{2}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)$
(26)

where $\mathcal{E}/L$ represents some proportionality of the impact parameter
of a photon within the black hole spacetime. Here it is a constant of motion.
However, in practice it is not that much reliable while observing a shadow.
Instead a photon sphere or something like that may be a more useful parameter.
That is why we express $L^{2}/\mathcal{E}^{2}$ in terms of the position of the
photon sphere (or any equivalent quantity) of this spacetime. Evaluating the
condition $\frac{dR}{d\phi}=0$ one can get the extremum points of light paths,
$R=R_{m}$ [21] from the following equation

$\dfrac{L^{2}}{\mathcal{E}^{2}}=R_{m}^{2}\left(1-\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R_{m}^{n}}-\dfrac{2R_{m}}{r_{0}}\right)^{-1}$
(27)

The radial coordinate $R_{m}$ will be the position of photon sphere if it
satisfies the condition $\left.\frac{d^{2}R}{d\phi^{2}}\right|_{R=R_{m}}=0$
[21]. Imposing this we get

$\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}(\frac{n}{2}+1)r_{0}^{n}}{R_{m}^{n}}+\dfrac{R_{m}}{r_{0}}=1$
(28)

The solution of this equation will give the position of photon sphere. Now by
using Eqs. (26) and (27), we have

$\sin{\alpha}=\dfrac{R_{m}}{R}\sqrt{\dfrac{1-\sum\frac{a_{n}r_{0}^{n}}{R^{n}}-\frac{2R}{r_{0}}}{1-\sum\frac{a_{n}r_{0}^{n}}{R_{m}^{n}}-\frac{2R_{m}}{r_{0}}}}$
(29)

An observer situated at $R=R_{O}$ will measure the angular radius of the
shadow of a $D$-dimensional evaporating Vaidya black hole as

$\sin{\alpha_{sh}^{bh}}=\dfrac{R_{m}}{R_{O}}\sqrt{\dfrac{1-\sum\frac{a_{n}r_{0}^{n}}{R_{O}^{n}}-\frac{2R_{O}}{r_{0}}}{1-\sum\frac{a_{n}r_{0}^{n}}{R_{m}^{n}}-\frac{2R_{m}}{r_{0}}}}$
(30)

where the subscript $sh$ and the superscript $bh$ in $\alpha$ will mean the
shadow of the black hole for the corresponding observer. This formula clearly
depicts that the angular radius of this shadow is time independent for an
observer on the constant–$R$ line i.e. the boundary of the shadow profile,
which will be circular in shape due to the spherical symmetry of the
spacetime, will not evolve over time $T$ and this is quite obvious from the
fact that the spacetime is conformally static in $(T,R)$-frame, as shown
earlier.

##  IV Shadow of a Collapsing dark Star for a Static Observer

We start by considering the free fall of a radially ingoing massive particle
(residing just above the star surface) on the equatorial plane of this
spacetime. This will mimic the star collapse. The Lagrangian and the
conformally conserved energy of such a vertically infalling massive test
particle are given by Eqs. (10) and (12) for $\dot{\phi}=0$. The Lagrangian of
a vertically infalling massive particle following Eq. (10) implies

$\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\dot{T}^{2}+2\dot{T}\dot{R}=e^{-2T/r_{0}}$
(31)

and the conformal constant of motion (from Eq. (12)) corresponding to the
Killing vector $\frac{\partial}{\partial T}$ takes the form

$-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\dot{T}-\dot{R}=e^{-2T/r_{0}}\mathcal{E}.$
(32)

In the following, we restrict ourselves to the condition $\dot{u}>0$, which
depicts that time is future oriented and $\dot{r}<0$, i.e. the particle is
infalling. Recalling the inverse coordinate transformations:
$T=r_{0}\ln{\dfrac{u}{r_{0}}}~{}~{}\text{and}~{}~{}R=r_{0}\dfrac{r}{u}$ we get

$\displaystyle\dot{T}$ $\displaystyle=$ $\displaystyle
r_{0}\dfrac{\dot{u}}{u}~{}~{}\text{and}~{}~{}\dot{R}=\dfrac{r_{0}}{u^{2}}\left(u\dot{r}-r\dot{u}\right).$
(33)

The conditions of infall of a massive particle mentioned above automatically
imply $\dot{T}>0$ and $\dot{R}<0$. One does not find any solution set for
$\dot{T}$ and $\dot{R}$ admitting this infalling condition while solving Eqs.
(31) and (32). However, there is no physical reason behind the non-existence
of an infalling timelike geodesic in this spacetime as we know Vaidya
spacetime can describe the infall of massive particles. The spacetime under
consideration is just a coordinate-transformed version of the Vaidya metric.
Therefore, it should possess an infalling timelike geodesic which can
represent the radial infall of the star surface and hence star collapse. One
good way is to identify the conformally conserved energy, $E$ as
$-\mathcal{E}=-\partial\mathcal{L}_{D}/\partial\dot{T}$. Accordingly, we get

$\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)\dot{T}+\dot{R}=e^{-2T/r_{0}}E$
(34)

Solving Eq. (31) and Eq. (34), we obtain (for $\dot{T}>0$ and $\dot{R}<0$)

$\displaystyle\dot{T}$ $\displaystyle=$ $\displaystyle
e^{-2T/r_{0}}\left\\{E+\sqrt{E^{2}-e^{2T/r_{0}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)}\right\\}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)^{-1}$
(35) $\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle-e^{-2T/r_{0}}\sqrt{E^{2}-e^{2T/r_{0}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)}$
(36)

Further simplification of Eqs. (35) and (36) yields

$\dfrac{dR}{dT}=e^{-2T/r_{0}}E^{2}-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)-e^{-2T/r_{0}}E^{2}\sqrt{1-\dfrac{e^{2T/r_{0}}}{E^{2}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)}$
(37)

The relative velocity between a radially infalling observer and a static
observer outside the star in $(T,R)$ frame is found to be 111This has been
discussed in [43] for a general static $4$-dimensional case and the formula
for the relative velocity
$\left(\vec{v}\equiv\left(v_{R},0,0,...\right)\right)$ between a radially
infalling observer and a static observer has been derived in [44, page-250];
and it turns out that the proof is easily extendable to any higher dimensional
static case.

$v_{R}=\sqrt{1-\dfrac{e^{2T/r_{0}}}{E^{2}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)}$
(38)

We now calculate the orthogonal tetrads for a radially infalling observer in
the $(T,R)$ frame using the above results. For this purpose it is necessary to
generalize the family of orthogonal tetrads (denoted by $\tilde{e}_{\mu}$) for
any arbitrary observer in $D$-dimensions, moving with velocity
$\vec{v}\equiv\left(v_{1},v_{2},v_{3},...\right)$ with respect to a static
frame 222The method of generating the family of orthogonal tetrads in
$4$-dimensions with respect to a static frame has been discussed in [45, 46]..
Careful observations reveal that the extra dimensional coordinates being
angular in nature, the number of non-zero components of basis tetrads reduces
due to spherical symmetry and we are left with only $T$, $R$ and $\phi$
components. Here the $\phi$ component is identical with its counterpart in
$4$-dimensions for obvious reasons.

The orthogonal tetrads for a freely infalling observer with respect to a
static frame thus become

$\displaystyle\tilde{e}_{T}$ $\displaystyle=$
$\displaystyle\dfrac{e_{T}+v_{R}~{}e_{R}}{\sqrt{1-v_{R}^{2}}}=e^{-2T/r_{0}}\left[\dfrac{E-\sqrt{E^{2}-m_{*}e^{2T/r_{0}}}}{m_{*}}\dfrac{\partial}{\partial
T}-\sqrt{E^{2}-m_{*}e^{2T/r_{0}}}~{}\dfrac{\partial}{\partial R}\right]$ (39)
$\displaystyle\tilde{e}_{R}$ $\displaystyle=$
$\displaystyle\dfrac{v_{R}~{}e_{T}+e_{R}}{\sqrt{1-v_{R}^{2}}}=e^{-2T/r_{0}}\left[\dfrac{E-\sqrt{E^{2}-m_{*}e^{2T/r_{0}}}}{m_{*}}\dfrac{\partial}{\partial
T}-E~{}\dfrac{\partial}{\partial R}\right]$ (40)
$\displaystyle\tilde{e}_{\phi}$ $\displaystyle=$ $\displaystyle
e_{\phi}=\dfrac{e^{-T/r_{0}}}{R}\dfrac{\partial}{\partial\phi}~{}~{}\text{(on
the equatorial plane)}$ (41)

where
$m_{*}=1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}$
has been defined for the sake of brevity.

###  IV.1 Shadow profile of a collapsing star from a static frame in (T,R)
coordinates

In order to determine the escape angle for the collapsing star, we will follow
exactly the same method, as we did earlier in the case of a static observer
(in $(T,R)$ frame) in section (III). Let us denote this angle (celestial
angle) between the photon geodesic and the radial direction in the rest frame
of the observer by $\tilde{\alpha}$ and therefore (the scale factor in this
case is $\tilde{\kappa}>0$),

$\dot{T}~{}\dfrac{\partial}{\partial T}+\dot{R}~{}\dfrac{\partial}{\partial
R}+\dot{\phi}~{}\dfrac{\partial}{\partial\phi}=\tilde{\kappa}\left(\tilde{e}_{T}+\tilde{e}_{R}\cos{\tilde{\alpha}}+\tilde{e}_{\phi}\sin{\tilde{\alpha}}\right)$
(42)

Comparing the coefficients of $\dfrac{\partial}{\partial T}$,
$\dfrac{\partial}{\partial R}$ and $\dfrac{\partial}{\partial\phi}$ on both
sides of Eq. (42), we have

$\left(\dfrac{dR}{d\phi}\right)^{2}=\dfrac{R^{2}e^{-2T/r_{0}}}{\sin{{}^{2}\tilde{\alpha}}}\left(\sqrt{E^{2}-m_{*}e^{2T/r_{0}}}-E~{}\cos{\tilde{\alpha}}\right)^{2}$
(43)

Now using the expression of $\dfrac{dR}{d\phi}$ from Eq. (17), we have,

$\displaystyle\dfrac{E^{2}R^{2}}{L^{2}}-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)=\dfrac{e^{-2T/r_{0}}}{\sin{{}^{2}\tilde{\alpha}}}$
$\displaystyle\left\\{\sqrt{E^{2}-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)e^{2T/r_{0}}}+E~{}\cos{\tilde{\alpha}}\right\\}^{2}$
(44)

and finally by substituting $E^{2}/L^{2}=\mathcal{E}^{2}/L^{2}$ from Eq. (27),
we now get,

$\displaystyle\dfrac{R^{2}}{R_{m}^{2}}\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R_{m}^{n}}-\dfrac{2R_{m}}{r_{0}}\right)-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)$
$\displaystyle=\dfrac{e^{-2T/r_{0}}}{\sin{{}^{2}\tilde{\alpha}}}\left\\{\sqrt{E^{2}-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R^{n}}-\dfrac{2R}{r_{0}}\right)e^{2T/r_{0}}}+E~{}\cos{\tilde{\alpha}}\right\\}^{2}$
(45)

One can obtain the expression for the celestial angle $\tilde{\alpha}$ from
the above equation. However, that will not be useful in calculating the shadow
of a collapsing star for a static observer. So, we skip that part. The shadow
is obtained by considering the null geodesics that are grazing the surface of
the star. One can (in principle) determine $R_{m}$ by equating $R=R_{s}$
(where $R_{s}$ is the radius of the star) and $\tilde{\alpha}=\pi/2$, if such
a null geodesic passes through a minimum radius value $R_{m}$. An observer on
the star surface will see its entire sky to be covered by the shadow.
Therefore we get from the above equation (45)

$\displaystyle\dfrac{R_{s}^{2}}{R_{m}^{2}}\left(1-\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R_{m}^{n}}-\dfrac{2R_{m}}{r_{0}}\right)-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R_{s}^{n}}-\dfrac{2R_{s}}{r_{0}}\right)$
$\displaystyle=e^{-2T/r_{0}}\left\\{E^{2}-\left(1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}r_{0}^{n}}{R_{s}^{n}}-\dfrac{2R_{s}}{r_{0}}\right)e^{2T/r_{0}}\right\\}$
(46)

Substituting the above expression in the equation of the escape angle of a
black hole in (30), we obtain

$\sin{\alpha_{sh}^{cs}}=\dfrac{R_{s}}{R_{O}}\frac{e^{T/r_{0}}}{\left|E\right|}\sqrt{1-\displaystyle\sum_{n\in\mathbb{Z}}a_{n}\dfrac{r_{0}^{n}}{R_{O}^{n}}-\dfrac{2R_{O}}{r_{0}}}$
(47)

where the superscript $cs$ stands for the shadow of the collapsing star. Using
the expressions of $\dot{T}$ in (35) and $\dot{R}$ in (36) one can in
principle obtain $\left|E\right|$ from the energy of the infalling observer
when the process of collapse starts, let’s say at radius $R_{m}=R_{i}$.
Substituting this in Eq. (47) we get the final expression of the escape angle
of the collapsing star as seen by a static observer in $(T,R)$ coordinate
system.

###  IV.2 Analysis in the original Eddington–Finkelstein-like Coordinates

In this section, our aim is to find the outcomes in the Eddington–Finkelstein
coordinates. We perform a coordinate transformation back to the original
coordinate system $(u,r,\theta,\phi)$, which will eventually enable us to
examine the characteristics of the collapsing star shadow from the viewpoint
of an observer sitting on a line with constant ($r,\theta,\phi$). Using the
inverse coordinate transformation and considering the radius of the star at
time $u(T)$ is $r_{s}=R_{s}e^{T/r_{0}}$ we obtain the escape angle of the
collapsing star as

$\sin{\alpha_{sh}^{cs}}=\dfrac{r_{s}u}{r_{0}r_{O}\left|E\right|}\sqrt{1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}u^{n}}{r_{O}^{n}}-\dfrac{2r_{O}}{u}}$
(48)

where the photon orbit is situated at $r_{m}$ and the radial position of the
observer is at $r_{O}$. Rescaling $u/r_{0}$, $r_{s}/r_{0}$ and $r_{O}/r_{0}$
(i.e. setting $r_{0}$ to unity) we get the angular radius of the shadow of a
collapsing star from the following relation

$\sin{\alpha_{sh}^{cs}}=\dfrac{r_{s}}{r_{O}}\dfrac{u}{\left|E\right|}\sqrt{1-\displaystyle\sum_{n\in\mathbb{Z}}\dfrac{a_{n}u^{n}}{r_{O}^{n}}-\dfrac{2r_{O}}{u}}.$
(49)

The only issue that can arise regarding the evaluation of $r_{s}$ at any
particular instant of time, can be managed by studying the collapse of the
star (or alternatively by studying a radial timelike inward geodesic in this
spacetime) with appropriate physical conditions (which has been studied by
Schneider and Perlick [20] in $4$-dimensions for the Schwarzschild metric at
the exterior).

##  V Simulating a collapsing star with a mass function asymptotically
approaching Schwarzschild–Tangherlini Black Hole Mass

The profile of a collapsing dark star surface for an evaporating black hole is
studied in this section following the method presented in [47]. We have
already emphasized that a radially infalling timelike geodesic can mimic the
collapse of the star surface. It will be sufficient to study the ingoing
geodesics of a massive free particle to visualize the collapse. The Lagrangian
for a radially infalling massive particle in retarded Eddington–Finkelstein
like coordinates is given by

$-\left(1-\dfrac{2m_{D}(u)}{r^{D-3}}\right)\dot{u}^{2}-2\dot{u}\dot{r}=-1$
(50)

Here we choose the Misner–Sharp mass function dependent only on $u$. The
equations of motion reduce to the following form with the above choice.

$\displaystyle\ddot{r}$ $\displaystyle+$
$\displaystyle\left(1-\dfrac{2m_{D}}{r^{D-3}}\right)\ddot{u}-\dfrac{1}{r^{D-3}}~{}\dfrac{dm_{D}}{du}~{}\dot{u}^{2}+\dfrac{2(D-3)m_{D}}{r^{D-2}}~{}\dot{u}\dot{r}=0$
(51) $\displaystyle\ddot{u}$ $\displaystyle=$
$\displaystyle(D-3)\dfrac{m_{D}}{r^{D-2}}~{}\dot{u}^{2}$ (52)

Substituting in Eq. (51) the expression of $\ddot{u}$ from Eq. (52) and
$\left(1-\dfrac{2m_{D}}{r^{D-3}}\right)\dot{u}^{2}$ from Eq. (50) we get a
simplified version of the same equation

$\ddot{r}+(D-3)\dfrac{m_{D}}{r^{D-2}}-\dfrac{1}{r^{D-3}}~{}\dfrac{dm_{D}}{du}~{}\dot{u}^{2}=0$
(53)

Let us now consider a simple model as suggested in [48] to proceed with the
solutions. According to the choice, the mass function approaches a constant
value $m_{D}^{bh}$ at a sufficiently large time, i.e. $m_{D}(u\to-\infty)\to
m_{D}^{bh}=G_{D}M_{D}^{bh}$, where $M_{D}^{bh}$ is interpreted as the
(resulting) black hole mass in $D$-dimensions. A detailed discussion on this
condition can be found in [48, 49]. Let us choose $M_{D}^{bh}=1$ to achieve
the following form of the mass function.

$m_{D}(u)=\dfrac{1}{2}m_{D}^{bh}\left(2+\operatorname{sech}{u}\right)=\dfrac{1}{2}G_{D}\left(2+\operatorname{sech}{u}\right)$
(54)

The evolution of the mass function and the expel profile (represented by the
derivative of the mass function) over the course of time for a collapsing star
obeying Vaidya metric at the exterior has been shown in the Fig. (1) below.
The mass function is peaked at $u=0$ and gradually decreases as one probes
deep into negative $u$ direction.

Figure 1: Evolution of Miser–Sharp mass parameter, $m_{D}(u)$ and its time
derivative

Now we find the solution of Eqs. (52) and (53) subject to the initial
conditions – (i) at $s=u=0$, $r=r_{i}$ ($r_{i}$ is the initial radius of the
star) and (ii) $\dot{u}=\dot{r}=0$. The solution has been shown below for
$D=4,5~{}\text{and}~{}6$-dimensions. The typical size of extra dimension has
been chosen as $\mathcal{R}$. Due to lack of our knowledge about the exact
size of the higher dimension, we have plotted $r(s)/\mathcal{R}$ versus
$s/\mathcal{R}$ instead of plotting only $s$ or $r(s)$.

Figure 2: Radius vs. proper time diagram of a collapsing star. The star starts
collapsing from outgoing Eddington–Finkelstein time $u=0$ with radius
$r_{i}=50$.

It is apparent from the graph that the rate of collapse is faster in four
dimensions compared to higher dimensions. The figure (2) shows a star with
same initial radius $\left(r_{i}=50\right)$ will collapse more rapidly in
four-dimensions than in the higher dimension(s). A more massive star in higher
dimension will exhibit same rate of collapse compared to a less heavy star in
four dimensions. This is, undoubtedly, a very interesting as well as important
feature. We now explore the shadow and photon sphere cast by this collapsing
star. For this we refer to the general Vaidya metric (3):

$ds_{D}^{2}=-f_{D}(u,r)du^{2}-2dudr+r^{2}d\Omega_{D-2}^{2}$ (55)

where $f_{D}(u,r)\equiv 1-\dfrac{2m_{D}(u)}{r^{D-3}}$. The position of the
photon sphere and the corresponding shadow in $4D$ have been obtained in [48].
We generalize the results in $D$-dimensions and extend it to the case of a
collapsing star possessing the same metric. Consider $\alpha$ and $\beta$ as
the celestial coordinates that define a two-dimensional celestial plane
perpendicular to the observer’s line of sight and located at spatial infinity
[50]. The resulting shadow configuration turns out to be

$\alpha(u)^{2}+\beta(u)^{2}=\dfrac{r_{m}^{2}}{f(u,r_{m})}\left[1-\left\\{\dfrac{\frac{dr_{m}}{du}}{f(u,r_{m})+\frac{dr_{m}}{du}}\right\\}^{2}\right]$
(56)

The photon sphere at $r_{m}\equiv r_{m}(u)$ in this spacetime is obtained by
solving the following differential equation [48, Eq. 2.9]

$\displaystyle\dfrac{d^{2}r_{m}}{du^{2}}-\left\\{\dfrac{3}{r_{m}}f(u,r_{m})-\dfrac{3}{2}\left.\dfrac{\partial
f}{\partial
r}\right|_{u,r_{m}}\right\\}\dfrac{dr_{m}}{du}-\dfrac{2}{r_{m}}\left(\dfrac{dr_{m}}{du}\right)^{2}$
$\displaystyle+\dfrac{1}{2}\left\\{f(u,r_{m})\left.\dfrac{\partial f}{\partial
r}\right|_{u,r_{m}}+\left.\dfrac{\partial f}{\partial
u}\right|_{u,r_{m}}\right\\}-\dfrac{1}{r_{m}}f(u,r_{m})^{2}$ $\displaystyle=0$
(57)

Substituting $f(u,r)$ for $D$ dimensions and doing some manipulations we get
the following form of the above equation

$\displaystyle\dfrac{d^{2}r_{m}}{du^{2}}$ $\displaystyle-$ $\displaystyle
3\left\\{\dfrac{1}{r_{m}}-(D-1)\dfrac{m_{D}(u)}{r_{m}^{D-2}}\right\\}\dfrac{dr_{m}}{du}-\dfrac{2}{r_{m}}\left(\dfrac{dr_{m}}{du}\right)^{2}-\dfrac{1}{r_{m}}$
(58) $\displaystyle-$ $\displaystyle
2(D-1)\dfrac{m_{D}(u)^{2}}{r_{m}^{2D-5}}+(D+1)\dfrac{m_{D}(u)}{r_{m}^{D-2}}-\dfrac{1}{r_{m}^{D-3}}\dfrac{dm_{D}}{du}=0$

The time evolution of the photon sphere of a Vaidya spacetime, possessing mass
function (54), has been studied by solving Eq. (58). We impose the asymptotic
condition – that after a sufficiently large time, the star will collapse to
form a Schwarzschild–Tangherlini black hole to solve the above equation. The
mass function is chosen accordingly such that it follows the above condition.

$\displaystyle\displaystyle\lim_{u\to-\infty}r_{m}(u)$ $\displaystyle=$
$\displaystyle\left[\dfrac{4\Gamma{\left(\frac{D-1}{2}\right)}G_{D}}{\pi^{\frac{D-3}{2}}}\left(\dfrac{D-1}{D-2}\right)\right]^{\frac{1}{D-3}}$
(59)
$\displaystyle\text{and}~{}~{}\displaystyle\lim_{u\to-\infty}\dfrac{dr_{m}(u)}{du}$
$\displaystyle=$ $\displaystyle 0$ (60)

From the above conditions we find that a Schwarzschild–Tangherlini black hole
has its photon sphere at a position given by the first condition [51],[37] and
the second condition indicates that the position is not changing over the
course of time. We have shown time variation of the position of photon sphere
of a Vaidya black hole with mass function (54) in Figure (3) where
$\mathcal{R}$ represents the size of each extra dimension. As evident from the
graph, the photon sphere will shrink more slowly as we go to higher
dimension(s). To combat with our lack of knowledge about the exact size of the
extra dimension, we have plotted $u/\mathcal{R}$ in the $x$-axis and
$r_{m}/\mathcal{R}$ in the $y$-axis instead of plotting only $u$ or $r_{m}$.

Figure 3: Evolution of the photon sphere over the course of time

The relative time evolution of the photon sphere can be visualized more
clearly from the evolution of the corresponding collapsing star shadow with
$u$. Here we have substituted numerical values of $r_{m}$ and its first
derivatives in Eq. (56) to plot the shadow at nine different snapshots of
time.

(a)

(b)

(c)

Figure 4: Profile of collapsing star shadow for the same set of physical
conditions, but in different dimensions

These plots depict the fact that in spite of the _similar_ physical conditions
(i.e. the boundary conditions) in all the cases, the collapsing star shadow
shrinks more faster in $4$-dimensions than that in the higher dimensions and
this is quite expected and obvious from our previous analyses.

##  VI Signature of Extra Spatial Dimensions

Let us critically examine the analytical results obtained in the preceding
sections. The analytical expression of the escape angle of a collapsing star
in Eq. (49) apparently does not depend on the dimension of spacetime. This is
not the actual picture. A close look at the expression reveals that the
apparent independence of the dimension of spacetime is due to the specific
form of the mass parameter $m_{D}(u,r)$. It has been chosen in such a way in
Eq. (7), that the metric on the equatorial plane has become independent of the
dimension of spacetime. Clearly, a degenerate shadow formation is possible if
one changes the mass parameter in the higher dimension in an appropriate way.
_Shadow can be entirely identical for two suitably chosen different mass
parameters in four and higher dimensions._

This can also be realized from the numerical analysis done in the previous
section. Figures (2), (3) and (4) clearly show that the stars with same
initial radius and final mass will collapse at different rates in different
dimensions and the shadow will also follow the same. One may end up predicting
higher or lower value of the mass of a star if the analysis is based only on
the observation of the shadow. If extra dimensions exist, a degeneracy in the
shadow profile will be inevitable. The possibility of existence of such a
shadow degeneracy and “isospectral” spacetimes have been discussed in [52] but
there all the “isospectral” spacetimes are in four dimensions and the approach
is slightly different.

We now ask the obvious question – whether the extra dimensions can be detected
from the isospectrality of spacetime? Or in other words, if different masses
cast the same shadow in different dimensions and by some means we get to know
which of those masses are actually from the star – can we detect the extra
dimensions? We present a proposal to give a positive answer to the above
question. To demonstrate this, we find the relationship between the mass of
the star in $D$-dimensional spacetime and the mass in four dimension, that
will cast the identical shadow. Note that we are proposing just a model. So,
instead of going into the complications of establishing the relationship
between a general mass parameter and the dimensionality of spacetime, we only
talk about the analytical form of the mass parameter, that has been considered
earlier (7).

$\dfrac{m_{D}}{m_{4}}=\dfrac{r_{O}^{D-3}\displaystyle\sum
a_{n}\left(u/r_{O}\right)^{n}}{r_{O}\displaystyle\sum
a_{n}\left(u/r_{O}\right)^{n}}=r_{O}^{D-4}$ (61)

Therefore,

$M_{D}\sim M_{4}\left(\dfrac{r_{O}}{\mathcal{R}}\right)^{D-4}$ (62)

where Newton’s constant in $D$ and four dimensions are related by the relation
$G_{D}\sim G_{4}\mathcal{R}^{D-4}$ for $\mathcal{R}$ being the typical size of
each of the extra dimensions. Using the current bound from LHC we present a
crude estimate on the mass $M_{D}$. Substituting the highest achievable energy
scale of LHC, i.e. $\sim 13.6$ TeV we derive an upper bound on the size
$\mathcal{R}$ as

$\mathcal{R}<\dfrac{hc}{13.6~{}\text{TeV}}=9.12\times 10^{-23}~{}\text{km}$
(63)

From Eq. 62 we find a lower bound on mass in $D>4$ dimensions

$M_{D}>\left(3.56\times 10^{5}~{}r_{O}\right)^{D-4}M_{4}$ (64)

Note that the typical size of a standard galaxy lies around a few kpc. That is
why we have chosen $r_{O}$ in kpc. Now to determine $M_{D}$ or $M_{4}$ which
one is actually the mass of that collapsing star we have to rely on
astrophysical methods that are generally used to determine the properties of
various dark objects in the sky.

It is worth noting that the particular shadow given by Eq. (49) can also be
cast by some other completely different spacetime. The degeneracy is a serious
issue here. We plan to explore this in detail in a future study. The focus of
present work was investigation of the $4$-dimensional isospectral counterpart
of the higher dimensional spacetime (3) and presenting the method of how the
higher dimensional feature of the spacetime can be extracted. The treatment is
quite general and can be used in specific case studies.

##  VII Concluding Remarks

In this work, higher dimensional Vaidya metric has been considered to model a
collapsing dark star. The analytical part of the work is restricted to a
special type of mass function for which we have calculated the escape angle
and the position of the photon sphere for black hole and collapsing dark star.
We have also explored numerically shadow of the black hole and the collapsing
star for an explicit choice of mass function, which is extendable for any mass
function under certain physical conditions. This has been elaborated in the
text. We have further shown how this shadow can be used to study the possible
signatures and distinctive features of higher dimensions. The analysis is done
for a spherically symmetric collapse, however, this method can also be applied
for a (slowly) rotating dark star. We plan to elaborate on this in a future
work. The major issues along this direction are two fold - (i) lack of a
suitable metric around a collapsing rotating star and (ii) of course the
mathematical complications involved. The second one can still be handled
somehow, but the first one is a major issue.

The procedure we discussed here can also be applied to search the higher
spatial dimensions from a black hole shadow [24, 37, 23, 25, 26, 27, 28].
However, in the case of a black hole shadow the difference between the masses
in different dimensions are generally much smaller and the effects are
observationally less prominent than that in our case. The mass of the star has
to be much different in four-dimensions than that in the higher dimensions for
the formation of same shadow by the collapsing dark star. This makes it much
easier from the astrophysical point of view to identify the correct one using
other physical parameters as mentioned in Section VI. Overall the process of
studying the shadow of a collapsing dark star is very much complex compared to
a blackhole shadow. As the Vaidya metric is very general and our work is based
on this geometry, the methodology presented here is widely applicable.

## References

  * Falcke  _et al._ [1999a] H. Falcke, F. Melia, and E. Agol, Astrophys. J. 528, L13 (1999a).
  * EHT1 [2019] EHT1, Astrophys. J. Lett. 875, L1 (2019).
  * Akiyama  _et al._ [2019a] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L3 (2019a).
  * 3 [2019] T. E. H. T. C. 3, Astrophys. J. Lett. 875, L4 (2019).
  * Akiyama  _et al._ [2019b] K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L5 (2019b).
  * 5 [2019] T. E. H. T. C. 5, Astrophys. J. Lett. 875, L6 (2019).
  * Falcke  _et al._ [1999b] H. Falcke, F. Melia, and E. Agol, Astrophys. J. 528, L13 (1999b).
  * Synge [1966] J. L. Synge, MNRAS 131, 463 (1966).
  * J. Bardeen and DeWitt [1973] B. D. J. Bardeen and C. DeWitt,  _Black Holes_ (New York : Gordon and Breach, 1973).
  * Ames and Thorne [1968] W. L. Ames and K. S. Thorne, Astrophys. J. 151, 659 (1968).
  * Jaffe [1969] J. Jaffe, Ann. Phys. 55, 10.1016/0003-4916(69)90184-5 (1969).
  * Lake and Roeder [1979] K. Lake and R. C. Roeder, Astrophys. J. 232, 10.1086/157286 (1979).
  * Frolov  _et al._ [2007] V. P. Frolov, K. Kim, and H. K. Lee, Phys. Rev. D 75, 10.1103/physrevd.75.087501 (2007).
  * Kong  _et al._ [2014] L. Kong, D. Malafarina, and C. Bambi, Eur. Phys. J. C 74, 10.1140/epjc/s10052-014-2983-3 (2014).
  * Kong  _et al._ [2015] L. Kong, D. Malafarina, and C. Bambi, Phys. Lett. B 741, 82 (2015).
  * Ortiz  _et al._ [2015a] N. Ortiz, O. Sarbach, and T. Zannias, Phys. Rev. D 92, 10.1103/physrevd.92.044035 (2015a).
  * Ortiz  _et al._ [2015b] N. Ortiz, O. Sarbach, and T. Zannias, Observational distinction between black holes and naked singularities: the role of the redshift function (2015b), arXiv:1401.4227 [gr-qc] .
  * Claudel  _et al._ [2001] C.-M. Claudel, K. S. Virbhadra, and G. F. R. Ellis, J. Math. Phys. 42, 818 (2001).
  * Perlick and Tsupko [2022] V. Perlick and O. Y. Tsupko, Phys. Rep. 947, 1 (2022).
  * Schneider and Perlick [2018] S. Schneider and V. Perlick, Gen. Relativ. Gravit. 50, 10.1007/s10714-018-2379-z (2018).
  * Solanki and Perlick [2022] J. Solanki and V. Perlick, Phys. Rev. D 105, 10.1103/PhysRevD.105.064056 (2022).
  * Singh and Ghosh [2017] B. P. Singh and S. G. Ghosh, Shadow of schwarzschild–tangherlini black holes (2017), arXiv:1707.07125 [gr-qc] .
  * Papnoi  _et al._ [2014] U. Papnoi, F. Atamurotov, S. G. Ghosh, and B. Ahmedov, Phys. Rev. D 90, 10.1103/physrevd.90.024073 (2014).
  * Singh  _et al._ [2023] B. P. Singh, R. Kumar, and S. G. Ghosh, New Astron. 99, 10.1016/j.newast.2022.101945 (2023).
  * Ahmed  _et al._ [2020] F. Ahmed, D. V. Singh, and S. G. Ghosh, 5d rotating regular myers–perry black holes and their shadow (2020), arXiv:2008.10241 [gr-qc] .
  * Amir  _et al._ [2018] M. Amir, B. P. Singh, and S. G. Ghosh, Eur. Phys. J. C 78, 10.1140/epjc/s10052-018-5872-3 (2018).
  * Belhaj  _et al._ [2022] A. Belhaj, M. Benali, and Y. Hassouni, Eur. Phys. J. C 82, 10.1140/epjc/s10052-022-10564-x (2022).
  * Nozari and Saghafi [2023] K. Nozari and S. Saghafi, Asymptotically locally flat and ads higher-dimensional black holes of einstein–horndeski–maxwell gravity in the light of eht observations: shadow behavior and deflection angle (2023), arXiv:2305.17237 [gr-qc] .
  * Banerjee  _et al._ [2022] I. Banerjee, S. Chakraborty, and S. Sengupta, Phys. Rev. D 106, 084051 (2022).
  * Barrabès and Israel [1991] C. Barrabès and W. Israel, Phys. Rev. D 43, 10.1103/PhysRevD.43.1129 (1991).
  * Lake and Zannias [1991] K. Lake and T. Zannias, Phys. Rev. D 43, 10.1103/PhysRevD.43.1798 (1991).
  * Ojako  _et al._ [2020] S. Ojako, R. Goswami, S. D. Maharaj, and R. Narain, Class. Quantum Grav. 37, 10.1088/1361-6382/ab5e2d (2020).
  * Wang and Wu [1999] A. Wang and Y. Wu, Gen. Relativ. Gravit. 31, 10.1023/a:1018819521971 (1999).
  * Vaidya [1951] P. C. Vaidya, Phys. Rev. 83, 10.1103/PhysRev.83.10 (1951).
  * Iyer and Vishveshwara [1989] B. Iyer and C. Vishveshwara, Pramana 32, 10.1007/BF02845995 (1989).
  * Chatterjee  _et al._ [1990] S. Chatterjee, B. Bhui, and A. Banerjee, J. Math. Phys. 31, 10.1063/1.528626 (1990).
  * Singh and Ghosh [2018] B. P. Singh and S. G. Ghosh, Ann. Phys. 395, 10.1016/j.aop.2018.05.010 (2018).
  * Oppenheimer and Snyder [1939] J. R. Oppenheimer and H. Snyder, Phys. Rev. 56, 10.1103/PhysRev.56.455 (1939).
  * Hiscock [1981a] W. A. Hiscock, Phys. Rev. D 23, 10.1103/PhysRevD.23.2813 (1981a).
  * Hiscock [1981b] W. A. Hiscock, Phys. Rev. D 23, 10.1103/PhysRevD.23.2823 (1981b).
  * Kuroda [1984] Y. Kuroda, Prog. Theor. Phys. 71, 10.1143/PTP.71.1422 (1984).
  * Beciu [1984] M. Beciu, Phys. Lett. A 100, 10.1016/0375-9601(84)90665-0 (1984).
  * Crawford and Tereno [2002] P. Crawford and I. Tereno, Gen. Relativ. Gravit. 34, 10.1023/a:1021131401034 (2002).
  * Landau and Lifshitz [1975] L. D. Landau and E. M. Lifshitz,  _The Classical Theory of Fields_ , 3rd ed., Course of Theoretical Physics, Vol. 2 (Pergamon Press, Oxford, 1975).
  * Grenzebach [2015] A. Grenzebach, in  _Fundamental Theories of Physics_ (Springer International Publishing, 2015) pp. 823–832.
  * Grenzebach [2016] A. Grenzebach, The shadow of black holes, in  _The Shadow of Black Holes: An Analytic Description_ (Springer International Publishing, 2016) pp. 55–79.
  * Piesnack and Kassner [2022] J. Piesnack and K. Kassner, Am. J. Phys. 90, 10.1119/10.0006367 (2022).
  * Mishra  _et al._ [2019] A. K. Mishra, S. Chakraborty, and S. Sarkar, Phys. Rev. D 99, 10.1103/PhysRevD.99.104080 (2019).
  * Coudray and Nicolas [2021] A. Coudray and J.-P. Nicolas, Gen. Relativ. Gravit. 53, 10.1007/s10714-021-02839-7 (2021).
  * Vázquez and Esteban [2004] S. Vázquez and E. Esteban, Il Nuovo Cim. B 119, 10.1393/ncb/i2004-10121-y (2004).
  * Décanini  _et al._ [2011] Y. Décanini, A. Folacci, and B. Raffaelli, Class. Quantum Grav. 28, 10.1088/0264-9381/28/17/175021 (2011).
  * Junior  _et al._ [2021] H. C. L. Junior, L. C. Crispino, P. V. Cunha, and C. A. Herdeiro, Phys. Rev. D 103, 10.1103/physrevd.103.084040 (2021).
  * Vertogradov and Kudryavcev [2022] V. Vertogradov and D. Kudryavcev, Generalized vaidya spacetime: horizons, conformal symmetries, surface gravity and diagonalization (2022), arXiv:2212.07130 [gr-qc] .
  * Frolov and Lee [2005] V. P. Frolov and H. K. Lee, Phys. Rev. D 71, 10.1103/PhysRevD.71.044002 (2005).
  * Yoshino  _et al._ [2019] H. Yoshino, K. Takahashi, and K. ichi Nakao, Phys. Rev. D 100, 10.1103/physrevd.100.084062 (2019).