#  Exploring the potentiality of standard sirens to probe cosmic opacity at
high redshifts

Xiangyun Fu1111corresponding author<EMAIL_ADDRESS>Jianfei Yang1, Zhaoxia
Chen1, Lu Zhou1, and Jun Chen2 1Institute of Physics, Hunan University of
Science and Technology, Xiangtan, Hunan 411201, China
2School of Science, Kaili University, Kaili, Guizhou 556011, China

###### Abstract

In this work, using the Gaussian process, we explore the potentiality of
future gravitational wave (GW) measurements to probe cosmic opacity at high
redshifts through comparing its opacity-free luminosity distance (LD) with the
opacity-dependent one from the combination of Type Ia supernovae (SNIa) and
gamma-ray bursts (GRBs). The GW data, SNIa and GRB data are simulated from the
measurements of the future Einstein Telescope, the actual Pantheon compilation
and the latest observation of GRBs compiled by L. Amati et al, respectively. A
nonparametric method is proposed to probe the spatial homogeneity of cosmic
transparency at high redshift by comparing the LD reconstructed from the GW
data with that reconstructed from the Pantheon and GRB data. In addition, the
cosmic opacity is tested by using the parametrization for the optical depth,
and the results show that the constraints on cosmic opacity are more stringent
than the previous ones. It shows that the future GW measurements may be used
as an important tool to probe the cosmic opacity in the high redshift region.

$\mathbf{Keywords:}$ Cosmic opacity, gravitation wave, Type Ia supernovae,
gamma-ray bursts

###### pacs:

98.80.Es, 95.36.+x, 98.80.-k

##  I Introduction

In 1998, the evidence of the accelerating expansion of the Universe was first
revealed by the unexpected dimming of the type Ia supernovae (SNIa) Riess1998
; Perlmutter1999 . In the frame of the General Relativity, a cosmic
distribution of an exotic component with negative pressure, dubbed as dark
energy, has been suggested to explain the present acceleration. On the other
hand, a cosmological distribution of dust has been proposed to be an
alternative explanation for this dimming phenomenon Aguirre1999 . Indeed, the
photons may be absorbed or scattered by dust in the Milky Way, intervening
galaxies, the intergalactic medium, and their host galaxies Lima2011 .
However, it has always been a controversial topic that whether the extinction
effect of SNIa has a great impact on the conclusion of cosmic accelerating
expansion Lima2011 ; Kang2020 ; Rose2020 ; Tutusaus2017 ; Huillier2019 . In
addition, some other plausible mechanisms attempt to explain cosmic opacity,
such as scalar fields coupled nonminimally with the electromagnetic (EM)
Lagrangian Hees2014 ; Holanda2017 ; Holanda20161 ; Aguirre1999  or oscillation
of photons propagating in extragalactic magnetic fields into light axions
Csaki2002 ; Avgoustidis2010 ; Avgoustidis2009 ; Jaeckel2010 . Exotic
mechanisms for cosmic opacity are not fully understood yet. Any changes in the
photon flux during propagation towards the Earth will affect the luminosity
distance (LD) measurements of light sources. As the cosmic acceleration rate
and the cosmological parameters determined by the LD measurements are highly
dependent on the dimming effect, cosmic opacity still needs to be investigated
accurately in the era of precision cosmology.

The general tests on cosmic opacity have been performed by using the cosmic
distance-duality relation (CDDR), which connects the LD $D_{\rm L}$ with the
angular diameter distance (ADD) $D_{\rm A}$ at the same redshift $z$ through
the identity eth1933 : ${D_{\rm L}}{(1+z)}^{-2}/{D_{\rm A}}=1$. Provided that
light travels always along null geodesics in a Riemannian geometry and the
number of photons is conserved, this reciprocal relation holds true in
whatever cosmology ellis1971 ; ellis2007 . If Riemannian geometry is used as
the tool to describe the spacetime of the Universe and a photon traveling
along null geodesic is more fundamental than the conservation of the photon
number uzan2004 ; Santana2017 , any violations of the CDDR likely indicate the
opacity of the Universe. Many works have been devoted to probe cosmic opacity
with the SNIa, baryon acoustic oscillation (BAO), the galaxy cluster samples,
the Hubble data, the old passive galaxies, and the gas mass fraction of galaxy
clusters More20092 ; Avgoustidis2010 ; Avgoustidis2009 ; Liao2016 ; Jesus2016
; Holanda20172 ; Liao2011 ; Liao2015 ; Li2013 . No significant cosmic opacity
is obtained in these studies.

It should be noted that most of the measurements are using EM radiation,
although the cosmological parameters can be measured precisely from various
astronomical observations. Schuts first proposed that Schutz , using the fact
that the waveform signals of GWs from inspiraling and merging compact binaries
encode distance information, the Hubble constant can be determined from GW
observations. So, this type of GW sources can be considered as standard sirens
in astronomy, analogous to SNIa standard candles. Compared with the distance
estimation from SNIa measurements, the LDs can be obtained directly from the
GW signals without the need of a cosmic distance ladder since stand sirens are
self-calibrating. It is more important that GWs are able to propagate freely
through a perfect fluid without any absorption and dissipation in the
Friedman-Lemaître-Robertson-Walker metric. Therefore, most objects (e.g. the
Earth) are nearly transparent to GWs Weihao2019  and the standard siren can be
considered as an opacity-free distance indicator. Generally speaking, it is
hard to obtain the redshift of the GW source. But the source redshift might be
obtained from the observation of EM counterparts that occur coincidentally
with the GW events Zhao2011 ; ET ; Nissanke2010 ; Cai2017 , if compact
binaries are neutron star (NS)-NS or black hole (BH)-NS binaries. Recently,
the joint detection of GW event GW170817 with EM counterpart (GRB 170817A)
from the merger of binary NSs Abbott ; Abbott2 ; Daz ; Cowperthwaite  opened a
new era of multi-messenger cosmology. Thus, the LD-redshift relation can be
constructed by a cosmological model-independent way, and it can be used to
make constraints on the parameters of cosmology.

The Einstein Telescope (ET), which is the third generation of the ground-based
detector, will be able to detect GW signals ten times more sensitive in
amplitude than the advanced ground-based detectors, and it can cover a wide
range of frequency $1\sim 10^{4}$ Hz up to redshift $z\sim 2$ for the NS-NS
and $z>2$ for the BH-NS mergers systems. So, the measurements from ET will
provide us with an opportunity to make constraint on the parameters of
cosmology at relative higher redshifts. Up to now, the simulated GW data have
been used to measure the cosmological parameters Pozzo201217 ; Liao2017 ;
Cai2016 ; Wei2017 ; Belgacem2019 ; Jin2000 ; Agostino ; ZhangJ2019 ; YangW2019
, determine the total mass of neutrino Wang2018 , investigate the anisotropy
of the Universe Cai2018 ; Wei2018 , constrain the evolving Newton’s constant G
Zhao2018 , discuss the estimation of the Hubble parameter with the actual
background expansion model of the universe Shafieloo2020 ; Keeley2020 , and
test the CDDR Yang2017 ; fuxiangyun2019 ; Liaokai .

More recently, confronting the LD from Joint Light Analysis (JLA) SNIa or
Pantheon compilation with the opacity-free LD from future observational GW
data, Wei JUN-JIE WEI2019 , Qi et al. Jing-Zhao Qi2019 , and Zhou et al.
fuxiangyun20192  performed successively unbiased tests on the cosmic opacity
with different methods. The results showed that future GW measurements would
be at least competitive with current tests on cosmic opacity. However,
regarding the tests referred above, it is important to note that due to the
limitation of the redshift distribution observed by SNIa, the tests are only
limited to the redshift range $0<z<2.0$. The cosmic opacity in the redshift
region between the SNIa and cosmic microwave background (CMB) observations
remains to be explored. On the other hand, since current astronomical
observation involving SNIa and other available data cannot yet confirm whether
the density of dark energy is constant or time-varying, the LDs of gamma-ray
bursts (GRBs) have been used to constrain cosmological parameters to determine
the behaviors of dark energy at higher redshifts $2<z<8$. The measurements of
GRBs can be considered as a new means to make constraint on the property of
dark energy, once phenomenological relations are adopted to make the sources
of bursts to be standardized candles Amati ; Schaefer ; Samushia2010 ; Wei2013
. Thus, to have confidence in the analyses with the GRBs observations, it is
important to test cosmic opacity in the high redshift region as well. Holanda
et al. tried to probe the cosmic opacity with the LDs of SNIa, Hubble
parameter $H(z)$ , and GRBs data at high redshifts ($z>1.4$) for a flat
cosmological constant and cold dark matter ($\Lambda$CDM) model Holanda2014 ,
and found that the observational data were compatible with a transparent
universe at $1\sigma$ confidence level (CL). However, it should be noted that
the analyses in Ref. Holanda2014  are model-dependent, thus some bias might
result from some particular cosmological models. To draw firm and robust
conclusion on the behaves of dark energy, one needs to employ new
observational data and methods to probe the cosmic opacity by creasing the
redshift depth. Therefore, future GW measurement and GRB observation at high
redshifts will be reliable candidates for this task, which is also the
motivation of this paper.

In this work, we explore the potentiality of future GW measurements to probe
the cosmic opacity at high redshifts through comparing its opacity-free
luminosity distance (LD) with the opacity-dependent one from the combination
of SNIa and GRBs data. GW data points are simulated from the ET, the SNIa data
are simulated from the Pantheon compilation Scolnic2018 , and GRB data are
simulated from the latest catalogue which are obtained with a cosmological
model-independent technique to overcome the circularity problemAmati2019 . To
match the GW data with the combination of SNIa and GRB data at the same
redhsift, we employ the Gaussian process Seikel2012 ; Rasmussen2006  to
reconstruct the continuous LD function from the mock GW measurements or the
combination of SNIa and GRB data. In our analysis, we first probe the spatial
homogeneity of the cosmic transparency by using non-parametric method. Then,
we test the cosmic opacity by adopting two types of parameterizations for the
optical depth. Compared to previous results, our analyses show that future
measurements of GW will be an important tool to probe cosmic opacity at high
redshifts.

##  II data and the Gaussian process

To dodge the impacts of any systematics from the combination of SNIa and GRB
datasets, we compare the opacity-free LD from the mocked measurements of
future GW events with the observed one obtained from the combination of
‘Pantheon-like’ SNIa and ‘GRBs-like’ compilations.

###  II.1 Pantheon SNIa and GRB data

The Pantheon SNIa compilation is released by the Pan-STARRS1 Medium Deep
Survey, which consists of 1048 data points up to redshift $z\sim 2.26$. The
LDs of the Pantheon compilation are calibrated from the SALT2 light-curve
fitter through applying the Bayesian Estimation Applied to Multiple Species
with Bias Corrections method to determine the nuisance parameters and taking
account of the distance bias corrections. The distance modulus and its
uncertainty can be obtained from $\mu=m_{\rm B}-M_{\rm B}$ and
$\sigma_{\mu}=\sigma_{m_{\rm B}}$, where $m_{\rm B}$ and $M_{\rm B}$ is the
rest-frame peak magnitude in the B band and the absolute magnitude of SNIa,
respectively. Then the LD ($D_{\rm L}$) can be derived from the relation

$\displaystyle\mu(z)=5\log_{10}[D_{\rm L}(z)]+25,$ (1)

and its uncertainty $\sigma_{D_{\rm L}}={\ln{10}}~{}\sigma_{\mu}D_{\rm L}$/5.

The GRBs can be observed up to redshift $z\sim 10$ due to the intense
explosions in the Universe. An important observational aspect of long GRB is
the several correlation between the spectral and intensity properties, which
suggests that the GRB measurement can be used as a complementary cosmic probe
to the standard candles Schaefer ; Bromm ; Amati ; Norris ; Fenimore ;
Schaefer2003 ; Ghirlanda ; Liang2005 ; Firmani ; Yu2009 ; Demianski2016 ,
although the mechanism behind GRBs explosions is not completely known yet. One
of the most successful proposals is the Amati relation ($E_{\rm p}$-$E_{\rm
iso}$ relation) Amati , which relates the rest frame spectral peak energy
$E_{\rm p}$ and the bolometric isotropic-equivalent radiated energy $E_{\rm
iso}$. However, the application of GRB observations for cosmology is also
affected by the so-called circularity problem Amati ; Ghirlanda ;
Demianski2016 , which arises from the fact that, given the lack of low
redshift measurements of GRBs, the $E_{\rm p}$-$E_{\rm iso}$ correlation is
obtained through assuming a background cosmology. For example, using the
standard $\Lambda$CDM model to calibrate the GRB measurements, it will
inevitably return an overall agreement with the concordance model to make
constraints on cosmological parameters of any dark energy framework. Thus,
calibrating the Amati relation is a challenge if one dodges the problem of
circularity. Recently, Amati et al. proposed a model-independent technique to
overcome the circularity problem with the Amati relation Amati2019 , and built
up a compilation consisting of 193 GRB data points. To calibrate this
relation, they fitted the most recent Hubble data with a B$\acute{\rm e}$zier
parametric curve, and employed the Hubble data to approximate the cosmic
evolution. Then, with the combination of this GRB data set and SNIa JLA
compilation, they made observational constraints on the flat $\Lambda$CDM
model, and obtained that $(\Omega_{\rm m},\alpha,\beta,{\rm M}_{\rm
B})=(0.397\pm_{0.039}^{0.040},0.126\pm_{0.012}^{0.011},2.610\pm
0.130,19.090\pm 0.037)$ at $2\sigma$ CL. Here, $\Omega_{\rm m}$ denotes the
present dark matter density parameter. The constraints on the cosmological
parameters are obviously different from the results of current observations
Scolnic2018 ; SDSS2014 .

It should be noted that, due to the little knowledge of the physical processes
driving the explosion, GRB observations, as the cosmological probing tools,
are not considered as important as the other astronomical measurements, such
as SNIa, BAO, and CMB. There is a strong debate in Ref. Heussaff2013
concerning whether the Amati relation is an intrinsic property of GRBs or
merely a combination of selection effects. Moreover, Ref. Salvaterra  studied
the strong evolution in the LD function of GRB, and the results showed that
GRBs might be intrinsically more luminous at high redshifts. If the Amati
relation is used to calibrate the high-redshift GRB with the low-redshift
ones, a smaller distance to GRB would be obtained.

Since the main purpose of this work is to probe the potentiality of future GW
measurements, it can be assumed that the problems of the GRB explosion will be
solved in the next few decades, and the GRBs observation can be used as a
complementary tool for SNIa measurements. In this work, in order to doge any
systematic impact while combining the SNIa data with the GRB data, SNIa and
GRB data are simulated from a flat $\Lambda$CDM. 137 deta points from the GRB
compilation in the redshift range $1<z<5$ are added to Pantheon SNIa data to
test cosmic opacity. The redshifts and uncertainties are taken from the actual
Pantheon Scolnic2018  and GRBs compilation Amati2019 . The fiducial LDs
($D_{\rm L}^{\rm fid}$) of the mock data are obtained from a known flat
$\Lambda$CDM with the most recent Planck results Plank2015 ,

$\displaystyle h_{0}=0.678,~{}~{}~{}\Omega_{\rm m}=0.308,$ (2)

where $H_{0}=100h_{0}{\rm kms^{-1}Mpc^{-1}}$ denotes the present Hubble
constant. The measurements of LD ($D_{\rm L}$) can be derived from the random
normal distribution through the equation $D_{\rm L}=\mathcal{N}(D_{\rm L}^{\rm
fid},\sigma_{D_{\rm L}})$. The mock SNIa and GRB data refer to ‘SNIa-like’ and
‘GRB-like’ data, respectively. The simulated method has been applied in
discussing the estimation of the Hubble parameter with different background
expansion model Keeley2020 . Our simulated results are presented in the left
panel of Fig. 1.

###  II.2 Simulated GW data

We simulate the GW data based on the ET which is the third-generation of
ground-based GW detector. The ET will be able to detect GW signals to be ten
times more sensitive in amplitude than the advanced ground-based detectors. It
can cover a wide range of frequency $1\sim 10^{4}$ Hz up to redshift $z\sim
5$. The strategy implemented with the future GW detectors has been discussed
in Refs. Regimbau2012 ; Regimbau2014 ; Regimbau2015 ; Regimbau20152 ;
Regimbau2016 ; Regimbau2017 . In the following of this subsection, we will
summarize the process of mocked GW data for simplification.

Given that the effect, which is from the spin of the binary system and the
change of orbital frequency over a single period, is negligible, for the
waveform of GW, one can apply the stationary phase approximation to compute
the Fourier transform $\mathcal{H}(f)$ of the time domain waveform $h(t)$,

$\displaystyle\mathcal{H}(f)=\mathcal{A}f^{-7/6}\exp[i(2\pi
ft_{0}-\pi/4+2\psi(f/2)-\varphi_{(2.0)})],$ (3)

where $\psi$ denotes the polarization angle, $\varphi_{(2.0)}$ denotes the
phase parameter, and the Fourier amplitude $\mathcal{A}$ is given by

$\displaystyle\mathcal{A}=$ $\displaystyle~{}~{}\frac{1}{D_{\rm
L}}\sqrt{F_{+}^{2}(1+\cos^{2}(\iota))^{2}+4F_{\times}^{2}\cos^{2}(\iota)}$
$\displaystyle~{}~{}\times\sqrt{5\pi/96}\pi^{-7/6}\mathcal{M}_{\rm
c}^{5/6}\,.$ (4)

Here, $\mathcal{M}_{\rm c}$ presents the chirp mass, $F_{+,\times}$ denotes
the beam pattern functions, ${D_{\rm L}}$ is the LD from GW signals, $t_{0}$
presents the epoch of the merger, and $\iota$ responds to the angle of
inclination. For the simulation estimation of LD, we use the flat $\Lambda$CDM
as the fiducial cosmological model with the model parameters in Equ. 2.

The signal-to-noise ratio (SNR) for the network from the three independent ET
interferometers can be written as

$\displaystyle\rho=\sqrt{\sum_{i=1}^{3}\langle\mathcal{H}^{(i)},\mathcal{H}^{(i)}\rangle}\,.$
(5)

Here, the inner product indicates the following equation

$\displaystyle\langle a,b\rangle=4\int_{f_{\rm lower}}^{f_{\rm
upper}}{\widetilde{a}(f)\widetilde{b}^{\ast}(f)+\widetilde{a}^{\ast}(f)\widetilde{b}(f)\over{2}}{df\over{S_{h}(f)}}\,,$
(6)

here ${S_{h}(f)}$ represents the one-side noise power spectral density
characterizing the performance of a GW detector. $f_{\rm lower}$ represents
the lower cutoff frequency which is fixed at 1 Hz, and $f_{\rm upper}$ does
the upper cutoff one decided by the last stable orbit (LSO), and $f_{\rm
upper}=2f_{\rm LSO}$, where $f_{\rm LSO}=1/(6^{3/2}2\pi M_{\rm obs})$ is the
orbit frequency at the last stable orbit. In this work, the masses of NS and
BH are simulated with uniform distribution in the intervals $[1,2]M_{\odot}$
and $[3,10]M_{\odot}$, respectfully, where $M_{\odot}$ denotes the solar mass.
The ratio between BH-NS and NS-NS binary systems is taken to be nearly 0.03
Cai2018 ; Abadie2010a . As presented in the Advanced LIGO-Virgo network
Abadie2010a ; Zhao2011 , the observational signal might be identified as a GW
events, only if the interferometers have a network SNR of $\rho>8.0$.

Using the Fisher information matrix, one can obtain the instrumental
uncertainty of the measurement of GW LD. It is assumed that the LD ($D_{\rm
L}$) is uncorrelated with the errors on the remaining GW parameters (the
inclination angle $\iota=0$), $\mathcal{H}\propto D_{\rm L}^{-1}$, and the
corresponding instrumental uncertainty can be written as, $\sigma_{D_{\rm
L}}^{\rm inst}\simeq{2D_{\rm L}/{\rho}}\,$. Furthermore, the LD is also
affected by the effect of the weak lensing, and the lensing uncertainty can be
assumed as $\sigma_{D_{\rm L}}^{\rm lens}=0.05zD_{\rm L}$ Cai2017 ; Zhao2011 ;
Belgacem2019 . Thus, the total uncertainty of GW LD measurements is taken to
be,

$\displaystyle\sigma_{D_{\rm L}}$ $\displaystyle=\sqrt{(\sigma_{D_{\rm
L}}^{\rm inst})^{2}+(\sigma_{D_{\rm L}}^{\rm lens})^{2}}$
$\displaystyle=\sqrt{\bigg{(}{2D_{\rm L}\over{\rho}}\bigg{)}^{2}+(0.05zD_{\rm
L})^{2}}\,.$ (7)

It is expected that the source redshift can be obtained from the identified EM
counterparts from the NS-NS and BH-NS binary systems. The redshift
distribution of GW sources is taken from the form of Zhao2011 ;
Sathyaprakash2010

$\displaystyle P(z)\propto{4\pi d^{2}_{\rm C}(z)R(z)\over{H(z)(1+z)}}\,,$ (8)

where $H(z)$ represents the Hubble parameter from the fiducial cosmological
model, $d_{\rm C}$ is comoving distance, and $R(z)$ is the merger rate of
binary systems with following expression Cai2017 ; Zhao2011 ,

$R(z)=\begin{cases}1+2z,&(z\leq 1)\\\ \frac{3}{4}(5-z),&(1<z<5)\\\ 0,&(z\geq
5).\end{cases}$ (9)

It is expected for the ET that the rates of NS-NS and BH-NS binary detections
per year are about the order $10^{3}-10^{7}$ ET . However, this prediction
about the rate is very uncertain. Then, Cai and Yang Cai2017 ; Yang2017
predicted that about $10^{2}$ GW measurements with EM counterparts will be
observed per year, if the middle detection rate around $10^{5}$ is viable and
the efficiency of EM counterparts are around $10^{-3}$ of the total number of
binary coalescence. In this paper, following the process from Ref. Zhao2011 ;
Yang2017 ; Cai2017 , we simulate 1000 data points in the redshift range
$0<z<5$, and show the results in the right panel of Fig. 1.

###  II.3 Gaussian process

In this work, we reconstruct the continuous function of LD from the mock GW or
EM observational data with the Gaussian process. The advantages of the
Gaussian process method are nonparametric and cosmological model-independent
smoothing technique used to reconstruct a continuous function from the
observed data. The reconstructed function, at each redshift point $z$, is
related to a Gaussian distribution with a mean value and Gaussian error bands.
The results from observational data points at any two redshifts $z_{i}$ and
$z_{j}$ are correlated through a covariance function $\kappa(z_{i},z_{j})$ due
to their nearness to each other. There is a wide range of possible candidates
for this function, and it depends on a set of hyperparameters. The different
choices of the covariance function may affect the reconstruction to some
extent. The detailed description and analysis of the Gaussian process can be
found in Refs. Seikel2012 ; Seikel20121 . This method has been widely used to
reconstruct the equation of state of dark energy Seikel2012 ; Cai2017 , and
test cosmography Shafieloo  and the CDDR Yang2017 ; fuxiangyun2019 . The
simplest covariance function usually takes the squared exponential as

$\kappa(z_{i},z_{j})=\sigma^{2}_{f}\exp\bigg{[}-{(z_{i}-z_{j})^{2}\over{2l^{2}}}\bigg{]},$
(10)

to implement the reconstructing process. Here, $\sigma_{f}$ and $l$ are two
hyperparameters. $\sigma_{f}$ denotes the output variance and fixes the
overall amplitude of the correction in the $y$ direction, and $l$ gives the
measure of the coherence length scale of the correlation in the $x$ direction.
As discussed in Ref. Keeley2020 , the hyperparameters play important role in
determining the error bars of observation data. Both of the hyperparameters
are optimized by the GP with the observed data set. In this work, we use the
simplest squared exponential covariance function. In order to avoid bias
caused by different initial values, we change the initial value of the super
parameter. After being optimized by the program, the returned values are about
$\sigma_{f}\sim 34.48$ and $l\sim 2.49$. The continuous functions of the LD
from the mock data are shown in Fig. 1. It can be seen from this figure that,
in the redshift region $2<z<5$, due to the small amount of GRB observation
data (only $71$ data points) and its large error bar, the uncertainty of LDs
reconstructed from the GRB data is greater than that of LDs reconstructed from
the GW measurements.

##  III Method

The LDs from future detectable GW sources are expected to be opacity-free. The
LDs from the electromagnetic observations, i.e. SNIa and GRB, are assumed to
be systematically influenced by the cosmic opacity, as pointed out by Ref.
Avgoustidis2009 . The straightforward method to probe the cosmic opacity is to
compare the observed LDs with the corresponding ones of GW measurements at the
same redshifts. The observed LD can be expressed as More2009 ,

${D_{\rm L,obs}}(z)={D_{\rm L,true}}(z)e^{\tau(z)/2},$ (11)

if the photon flux received by the observers will be reduced by the factor
$e^{-\tau(z)}$. Here $\tau(z)$ represents the optical depth related to the
cosmic absorption.

Figure 1:  The sample catalogs of the ‘pantheon-like’ data and ‘GRBs-like’
data (left panel) , and 1000 mock GW events of redshifts (right panel). The
corresponding smoothed function of LD obtained from the Gaussian process
method is also presented.

Then the function of cosmic optical depth at any redshifts can be obtained
through the following equation

${\tau(z)}=2\ln{D_{\rm L,GW}(z)\over{D_{\rm L,\,EM}(z)}}$ (12)

by comparing the continuous function of LD $D_{\rm L}(z)$ reconstructed from
the mock GW data with that from the EM observational data. Here and latter,
the subscript “EM” represents the observed variable obtained from the
electromagnetic wave measurements, such as the SNIa and GRB observations.
Deviations from a transparent universe at any redshifts will be encoded in the
function $\tau{(z)}$. The uncertainties of the optical depth from the
observations, $\sigma_{\tau_{\rm obs}}$, can be obtained from following
equation,

$\sigma_{\tau_{\rm obs}}=2\sqrt{\left[\left({\sigma_{D_{\rm
L,\,EM}(z)}\over{D_{\rm L,\,EM}(z)}}\right)^{2}+\left(\sigma_{D_{\rm
L,GW}(z)}\over{D_{\rm L,GW}(z)}\right)^{2}\right]}\,.$ (13)

The results are shown in the left panel of Fig. 2. It should be noted that the
constraint on cosmic opacity is independent of any parametrization.

Figure 2:  The cosmic optical depth function $\tau(z)$ (left panel) and the
likelihood distribution functions (right panel) obtained from the combination
of ‘Pantheon-like’ and ‘GRB-like’ data with $\Lambda$CDM.

In addtion, we employ two typical parameterizations for the optical depth
$\tau{(z)}$, i.e. $\tau(z)=2\epsilon z$ (P1) and $\tau(z)=(1+z)^{2\epsilon}-1$
(P2), to detect cosmic opacity in the redshift region $0<z<5$. While
$\epsilon=0$, the Universe is transparent. Here, the observational data of
opacity-free LD is from the simulated GW observation. To match the GW
measurement with the EM observation at the same redshift and employ all GW
data to test the cosmic opacity, the opacity-dependent one is from the
function of LD reconstructed from EM measurements with Gaussian process. The
probability density of $\epsilon$ can be obtained with the expression
$P(\epsilon)=A\,{\rm exp}[-\chi^{2}(\epsilon)/2]$, where, $A$ is a normalized
coefficient and $\chi^{2}(\epsilon)$ has the form

$\chi^{2}(\epsilon)=\sum\frac{{\left[\tau(z,\epsilon)-\tau_{\rm
obs}\right]}^{2}}{\sigma^{2}_{\tau_{\rm obs}}}\,.$ (14)

The results are shown in the right panel of Fig. 2.

Furthermore, in order to verify the effectiveness of this method, we also
simulated SNIa data and GRB data by assuming that the universe is opaque in
the case of optical depth function $\tau(z)=2\epsilon z$. The values of
$\epsilon$ are set to $\epsilon=0.01$ and $0.02$, respectively. Then, the
cosmic opacity are tested, and the results from the tests are shown in Fig. 3.

Figure 3:  The cosmic optical depth function $\tau(z)$ (left panel) and the
likelihood distribution functions (right panel) obtained from the combination
of ‘Pantheon-like’ and ‘GRB-like’ data with a certain model of opacity in
parametrization $\tau(z)=2\epsilon z$ with $\epsilon=0.01$ (top panel) and
$\epsilon=0.02$ (bottom panel), respectively.

##  IV Results and analyses

For the test on the spatial homogeneity of cosmic transparency in Fig. 2, it
can be seen that the best-fit value of cosmic optical depth $\tau(z)$ evolves
with the increase of the redshift. No deviation from a transparent universe is
found at $1\sigma$ CL. Our results are similar to those obtained in Refs.
Chen2012 ; fuxiangyun20192  in the range of SNIa measurement redshifts, in
which the best fit value of cosmic opacity varies with the increase of the
redshift. In addition, the value of cosmic optical depth $\tau(z)$ is
constrained between $-0.09$ and $0.05$ at $1\sigma$ CL in the redshift region
$0.1<z<2.0$, and it is constrained between $-0.14$ and $0.60$ in the redshift
region $2.0<z<5.0$. It suggests that the potentiality of observational data is
relatively weak to probe the spatial homogeneity of cosmic opacity at the high
redshift. As presented in Section III, due to a small number and large
uncertainty of GRB observations in this redshift region, the uncertainty of
LDs reconstructed from the GRB data is greater than that of LDs reconstructed
from the GW measurements. Therefore, the improvements of electromagnetic
observations in the high redshift region, including the improvements of data
quantity and quality, will impose more stringent constraints on the test of
cosmic opacity. The figure also shows that the cosmic optical depth $\tau(z)$
is divergent in the redshift range $z<0.1$, and this bias might result from
the absence of the mock GW data in this redshift region.

For the case of tests on cosmic opacity with parameterizations, as shown in
Tab. 1, we also present a comparison between our forecast results and the
previous ones. Compared the error bars of the constraint on cosmic opacity
with the previous ones from the SNIa data in the lower redshift range $0<z<2$,
our results are at least 50% less than those obtained from the simulated ET GW
data and the JLA SNIa data in Refs. JUN-JIE WEI2019 ; Jing-Zhao Qi2019 ;
fuxiangyun20192 , and they are at least 30% less than those from the simulated
ET GW data and Pantheon SNIa compilation in Ref. fuxiangyun20192 . Therefore,
the future GW measurements can be considered as an important tool to probe
cosmic opacity. For the results obtained from parameterizations for the
optical depth $\tau(z)$, it is easy to see that parametrization P1 may offer
more stringent constraint on cosmic opacity than P2.

Data |  $\epsilon\,\,\,(P_{1})$ | $\epsilon\,\,\,(P_{2})$
---|---|---
ET $\times$ 1048 Pantheon + 137 GRB∗ | ${-0.001{\pm{0.004}\pm{0.007}}}$ | ${-0.002{\pm{0.007}\pm{0.013}}}$
ET $\times$ 1048 Pantheon + 137 GRB${}^{\ast}\lozenge$ | ${0.010{\pm{0.004}\pm{0.007}}}$ | ${0.009{\pm{0.006}\pm{0.012}}}$
ET $\times$ 1048 Pantheon + 137 GRB${}^{\ast}\lozenge$ | ${0.019{\pm{0.004}\pm{0.007}}}$ | ${0.018{\pm{0.006}\pm{0.012}}}$
ET $\times$ 740 JLA (A) fuxiangyun20192  | ${0.004{\pm{0.008}}}$ | ${0.002{\pm{0.010}}}$
ET $\times$ 740 JLA (B) fuxiangyun20192  | ${0.009{\pm{0.011}}}$ | ${0.008{\pm^{0.014}_{0.012}}}$
ET $\times$ 1048 Pantheon (B) fuxiangyun20192  | ${0.005{\pm{0.006}}}$ | ${0.006{\pm{0.008}}}$
ET $\times$ 1048 Panthoen (B) JUN-JIE WEI2019  | ${0.004\pm 0.026}$ | $\Box$
ET $\times$ 740 JLA (B) Jing-Zhao Qi2019  | ${0.002\pm 0.035}$ | ${-0.006\pm 0.053}$
ET $\times$ 1048 Panthoen (B) Jing-Zhao Qi2019  | ${0.009\pm 0.016}$ | ${0.015\pm 0.025}$
581 SNIa + 19 H(z) $\times$ $\Lambda$CDM (A) Holanda2014  |  ${0.02\pm 0.055}$ $\triangle$ | $\Box$
59 GRB + 19 H(z) $\times$ $\Lambda$CDM (A) Holanda2014  |  ${0.06\pm^{0.18}_{0.18}}$ $\triangle$ | $\Box$
Table 1: Constraints on $\epsilon$ represented by the best fit value at
$1\sigma$ or $2\sigma$ CL for each data set. The superscripts ‘$\ast$’
represent the results obtained from the Gaussian process in this work, the
symbol ‘$\lozenge$’ refers to the results obtained from a certain model of
opacity, A and B represent the cosmological model-dependent method and the
model-independent method, respectively, and the triangle symbols ‘$\triangle$’
denote the parametrization is adopted in the form $\tau(z)=\epsilon z$.

For the test on the spatial homogeneity of the optical depth $\tau(z)$ from a
certain model of cosmic opacity while $\epsilon=0.01$, an opaque universe can
be found at $2\sigma$ CL in the redshift region $0.46<z<0.82$, as shown in the
left-top panel of Fig. 3. It is worth mentioning that the test on cosmic
opacity obtained directly from the observational data is independent of any
parametrization. For the test from the parametrization of the optical depth
$\tau(z)$, an opaque universe can be marginally found at $3\sigma$ CL for the
parametrization of P1, as shown in the right-top panel of this figure. While
$\epsilon=0.02$, it can be seen from bottom panel of Fig. 3 that an opaque
universe can be found at $3\sigma$ CL. It suggests that future gravitational
wave data can be used as an effective tool to identify an opaque universe.

##  V conclusion

In the past few years, some astronomical observations, such as the luminosity
distance (LD) of SNIa, baryon acoustic oscillations, gas mass fraction,
angular diameter distances from galaxy clusters, and Hubble parameter
measurements have been used to probe cosmic opacity. As a gravitational wave
propagates freely in a perfect fluid without any absorption and dissipation,
the measurement of the luminosity distance to a GW source provides us with the
opacity-free distance to probe cosmic opacity. More recently, confronting the
mock GW measurements to the SNIa compilations, some tests on cosmic opacity
have been performed to explore the potentiality of future GW measurements.
However, due to the limitation from the redshift distributions of
observational data, most of the tests are limited in the redshift region
$0<z<2$, and the cosmic opacity in the redshift region beyond the SNIa
observations remains unexplored.

The third-generation GW detector, i.e., the Einstein Telescope (ET), is
expected to detect GW signals up to redshift range $z\sim 5$. Thus, future GW
measurements will provide us with an opportunity to probe cosmic opacity at
high redshifts. In this work, using the Gaussian process, we explore the
potentiality of future gravitational wave measurements to probe the cosmic
opacity at high redshift through comparing LD from GW with that from the
combination of the SNIa and GRBs compilations. One thousand GW data points are
simulated from the future ET. To dodge any impacts of systematics while
combining the SNIa and GRB data sets, the ‘SNIa-like’ and ‘GRBs-like’ data are
simulated from the actual Pantheon SNIa compilation and the latest GRB
catalogue of Ref. Amati2019  in the $\Lambda$CDM cosmological model. 137 GRB
data points in the redshift region $1<z<5$ are added to the Pantheon
compilation. The purpose of using the Gaussian process is to reconstruct the
continuous function of LD from the future GW data or the combination of SNIa
and GRB data sets. A non-parametric method are proposed to probe the spatial
homogeneity of the cosmic transparency at any redshifts through comparing the
LD reconstructed from GW measurements with that reconstructed from SNIa and
GRB data. Furthermore, using the simulated GW data and the continuous function
of LDs reconstructed from the combination of SNIa data and GRB data, the
cosmic opacity is tested with two types of parameterizations for the optical
depth. All of the available GW data points can be used to test cosmic opacity.
To test the effectiveness of the method in this work, SNIa data and GRB data
are also simulated with a certain cosmic opacity model to probe the cosmic
spatial homogeneity and the cosmic opacity.

For the test on the spatial homogeneity of cosmic opacity, the results show
that the best fit value of cosmic opacity varies as the redshift increases,
and no deviation from a transparent universe is found at $1\sigma$ CL. Due to
the large uncertainty of the GRB observation in the high redshift region and
the small amount of data, the potentiality of astronomical observation to test
cosmic opacity in high redshift region is relatively weaker than that in the
low redshift region dominated by the SNIa observations. It suggests that
improvements from the electromagnetic observations at high reshifts, including
the improvements of data quantity and quality, would lead to more stringent
constraints on the cosmic opacity. Results, obtained from the test with
parameterizations for the optical depth, show that the error bar of constraint
on cosmic opacity can be reduced to $\sigma_{\epsilon}\sim 0.004$ at $1\sigma$
CL. The uncertainty of the constraint is less than the previous results
obtained from the simulated GW data, JLA SNIa data, and Pantheon SNIa data.
The tests from a certain model of opacity show that the GW measurements in the
future may be used as an effective tool to probe cosmic opacity. Therefore, it
can be concluded that, given the GW detector can be carried out as expected by
the programme of the ET, future measurements of GW may not only provide us
with an opportunity to investigate the spatial homogeneity of the cosmic
transparency at high redshifts, but also play an important role in probing
cosmic opacity.

###### Acknowledgements.

We very much appreciate helpful comments and suggestions from anonymous
referees. We are also thankful for helpful discussion from Puxun Wu and
Zhengxiang Li. This work was supported by the National Natural Science
Foundation of China under Grants No. 11147011, No. 11405052, No. 11465011, and
No. 11865011; the Hunan Provincial Natural Science Foundation of China under
Grant No. 12JJA001; the Foundation of Department of science and technology of
Guizhou Province of China under Grant No. J [2014] 2150; and the Foundation of
the Guizhou Provincial Education Department of China under Grant No. KY
[2016]104.

## References

  * (1) A. G. Riess, et al., Astron. J. 116, 1009 (1998).
  * (2) S. Perlmutter et al., Astrophys. J. 517, 565 (1999).
  * (3) A. Aguirre, Astrophys. J. 525, 583 (1999).
  * (4) J. A. S. Lima, J. V. Cunha, and V. T. Zanchin, Astrophys. J. 742, L26 (2011).
  * (5) Y. Kang, et al., Astrophys. J. 889, 8 (2020).
  * (6) B. M. Rose, et al. The Astrophysical Journal, Astrophys. J. 896, L4 (2020).
  * (7) I. Tutusaus, B., Lamine, A. Dupays, et al. Astron. Astrophys. 602, A73 (2017).
  * (8) B. L’Huillier, A. Shafieloo, E. V. Linder, and A. G. Kim, Mon. Not. R. Astron. Soc. 485(2), 2783 (2019).
  * (9) A. Hees, O. Minazzoli, and J. Larena, Phys. Rev. D 90, 124064 (2014).
  * (10) R. F. L. Holanda and S. H. Pereira, Phys. Rev. D 94, 104037 (2016).
  * (11) R. F. L. Holanda, S. H. Pereira, and S. Santos da Costa, Phys. Rev. D 95, 084006 (2017).
  * (12) A. Avgoustidis, L. Verde, and R. Jimenez, J. Cosmol. Astropart. Phys. 06 (2009) 012 .
  * (13) A. Avgoustidis, C. Burrage, J. Redondo, L. Verde, and R. Jimenez, J. Cosmol. Astropart. Phys. 10 (2010) 024 .
  * (14) J. Jaeckel and A. Ringwald, Annu. Rev. Nucl. Part. Sci. 60, 405 (2010).
  * (15) C. Csaki, N. Kaloper, and J. Terning, Phys. Rev. L 88, 161302 (2002).
  * (16) I. M. H. Etherington, Phil. Mag. 15, 761 (1933); reprinted in Gen. Relativ. Gravit. 39, 1055 (2007).
  * (17) G. F. R. Ellis, Proc. Int. School of Physics “Enrico Fermi” (Varenna), Course XLVII, edited by R. K. Sachs (Academic Press, New York, 1971), pp. 104 C182; reprinted in Gen. Relativ. Gravit. 41, 581 (2009).
  * (18) G. F. R. Ellis, Gen. Relativ. Gravit. 39, 1047 (2007).
  * (19) J. P. Uzan, N. Aghanim, and Y. Mellier, Phys. Rev. D 70, 083533 (2004).
  * (20) L. T. Santana, M. O. Calvao, R. R. R. Reis, and B. B. Siffert, Phys. Rev. D 95, 061501(R) (2017).
  * (21) S. More, J. Bovy, and D. W. Hogg, Astrophys. J. 696, 1727 (2009).
  * (22) K. Liao, et al., Astrophys. J. 822, 74 (2016).
  * (23) J. F. Jesus and R. F. L. Holanda, Gen. Relativ. and Gravit. 49, 150 (2017).
  * (24) R. F. L. Holanda, K. V. R. A. Silva, and V. C. Busti, J. Cosmol. Astropart. Phys. 03 (2018) 031.
  * (25) K. Liao, Z. Li, J. Ming, and Z. Zhu, Phys. Lett. B 718, 1166 (2013).
  * (26) K. Liao, A. Avgoustidis, and Z. Li, Phys. Rev. D 92, 123539 (2015).
  * (27) Z. Li, P. Wu, H. Yu, and Z. H. Zhu, Phys. Rev. D 87, 103013 (2013).
  * (28) B. P. Abbott et al., (LIGO Scientific and Virgo Collab- orations), Phys. Rev. Lett. 119, 161101 (2017).
  * (29) B. P. Abbott et al., (LIGO Scientific and Virgo and Fermi- GBM and INTEGRAL Collaborations), Astrophys. J. 848, L13 (2017).
  * (30) M. C. Daz, et al., (TOROS Collaboration), Astrophys. J. 848, L29 (2017).
  * (31) P. S. Cowperthwaite et al., Astrophys. J. 848, L17 (2017).
  * (32) B. F. Schutz, Nature (Longdon) 323, 310 (1986).
  * (33) H. Wei, D. C. Qiang, Z. X. Yu, and H. K. Deng, Communications in Theoretical Physics 72, 065401.
  * (34) W. Zhao, C. VanDenBroeck, D. Baskaran, and T. G. F. Li, Phys. Rev. D 83, 023005 ( 2011).
  * (35) The Einstein Telescope Project, https://www.et-gw.eu/.
  * (36) S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, and J. L. Sievers, Astrophys. J. 725, 496 (2010).
  * (37) R. G. Cai and T. Yang, Phys. Rev. D 95, 044024 (2017).
  * (38) W. Del Pozzo, Phys. Rev. D 86, 043011 (2012); W. Del Pozzo, T. G. F. Li, and C. Messenger, Phys. Rev. D 95, 043502 (2017).
  * (39) R. G. Cai, Z. K. Guo, and T. Yang, Phys. Rev. D 93, 043517 (2016).
  * (40) K. Liao, X. Fan, X. Ding, M. Biesiada, and Z. Zhu, Nature Commun. 8, 1148 (2017).
  * (41) J. Wei and X. Wu, Mon. Not. R. Astron. Soc. 472, 2906 (2017); J. Wei and X. Wu, and H. Gao, Astrophys. J. 860, L7 (2018).
  * (42) E. Belgacem et al., J. Cosmol. Astropart. Phys. 08 (2019) 015.
  * (43) S.-J. Jin, et al. J. Cosmol. Astropart. Phys. 03 (2020) 051.
  * (44) R. D’Agostino and R. C. Nunes, Phys. Rev. D 100, 044041 (2019).
  * (45) J.-F. Zhang, et al., J. Cosmol. Astropart. Phys. 09 (2019) 068.
  * (46) W. Yang, et al., J. Cosmol. Astropart. Phys. 07 (2019) 037.
  * (47) L. Wang, X. Zhang, J. Zhang, and X. Zhang, Phys. Lett. B 782, 87 (2018).
  * (48) J. Wei, Astrophys. J. 868, 29 (2018).
  * (49) R. G. Cai, T. B. Liu, X. W. Liu, S. J. Wang, and T. Yang, Phys. Rev. D 97, 103005 (2018).
  * (50) W. Zhao, B. S. Wright, and B. Li, J. Cosmol. Astropart. Phys. 10 (2018) 052.
  * (51) A. Shafieloo, R. E. Keeley, and E. V. Linder, J. Cosmol. Astropart. Phys. 03 ( 2020) 019.
  * (52) R. E. Keeley, et al., Mon. Not. R. Astron. Soc. 491, 3983 (2020).
  * (53) T. Yang, R. F. L. Holanda, and B. Hu, Astropart. Phys. 108, 57 (2019).
  * (54) X. Fu, L. Zhou, and J. Chen, Phys. Rev. D 99, 083523 (2019).
  * (55) K. Liao, Astrophys. J. 885, 70 (2019).
  * (56) J. Wei, Astrophys. J. 876, 66 (2019).
  * (57) J. Qi, S. Cao, Y. Pan, and J. Li, Phys. Dark Universe 26, 100338 (2019).
  * (58) L. Zhou, X. Fu, and J. Chen, Phys. Rev. D 100, 123539 (2019).
  * (59) B. E. Schaefer, Astrophys. J. 660, 16 (2007).
  * (60) L. Amati, et al., Astron. Astrophys. 390, 81 (2002); L. Amati, Mon. Not. R. Astron. Soc. 372, 233 (2006).
  * (61) L. Samushia and B. Ratra, Astrophys. J. 714, 1347 (2010).
  * (62) J. J. Wei, X. Wu, and F. Melia, Astrophys. J. 772, 43 (2013).
  * (63) R. F. L. Holanda and V. C. Busti, Phys. Rev. D 89, 103517 (2014).
  * (64) D. Scolnic, D. Jones, A. Rest, Y. Pan, R. Chornock, R. Foley, M. Huber, R. Kessler, G. Narayan, A. Riess, et al., Astrophys. J. 859, 101 (2018).
  * (65) L. Amati, et al., Mon. Not. R. Astron. Soc., 486, L46 (2019).
  * (66) C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (MIT Press, Cambridge, MA 2006), ISBN 026218253X.c.
  * (67) M. Seikel, C. Clarkson, and M. Smith, J. Cosmol. Astropart. Phys. 06 (2012) 036.
  * (68) M. Betoule et al. (SDSS collaboration), Astron. Astrophys. 568, A22 (2014).
  * (69) M. Sako, et al., Publ. Astron. Soc. Pac. 130, 064002 (2018).
  * (70) J. Guy et al. (SNLS collaboration), Astron. Astrophys. 523, A7 (2010).
  * (71) A. G. Riess et al., Astrophys. J. 659, 98 (2007).
  * (72) V. Bromm and A. Loeb, Astrophys. J. 575, 111 (2002); J. R. Lin, S. N. Zhang, and T. P. Li, Astrophys. J. 605, 819 (2004).
  * (73) J. P. Norris, G. F. Marani , and J. T. Bonnell, Astrophys. J. 534, 248 (2000); T. N. Ukwatta, et al., Astrophys. J. 711, 1073 (2010).
  * (74) D. E. Riechart, et al., Astrophys. J. 552, 57 (2001).
  * (75) B. E. Schaefer, Astrophys. J. 583, 71 (2003); D. Yonetoku, et al., Astrophys. J. 609, 935 (2004).
  * (76) G. Ghirlanda, G. Ghisellini, and D. Lazzati, Astrophys. J. 616, 331 (2004).
  * (77) E. W. Liang and B. Zhang, Astrophys. J. 633, 611 (2005).
  * (78) C. Firmani, et al., Mon. Not. R. Astron. Soc. 370, 185 (2006); F. Rossi, et al., Mon. Not. R. Astron. Soc. 388, 1284 (2008).
  * (79) B. Yu, et al., Astrophys. J. 705, 15 (2009).
  * (80) M. Demianski, E. Piedipalumbo, D. Sawant, and L. Amati, Astron. Astrophys. 598, A112 (2017); M. Demianski, E. Piedipalumbo, D. Sawant, and L. Amati, Astron. Astrophys. 598, A113 (2017).
  * (81) V. Heussaff, J.-L. Atteia, and Y. Zolnierowski, Astron. Astrophys. 557, A1100 (2013).
  * (82) R. Salvaterra, et al. Astrophys. J. 749, 68 (2012).
  * (83) P. A. R. Ade et al. [Planck Collaboration], Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594, A13 (2016).
  * (84) T. Regimbau, et al., Phys. Rev. D 86, 122001 (2012).
  * (85) T. Regimbau, D. Meacher, and M. Coughlin, Phys. Rev. D 89, 084046 (2014).
  * (86) T. Regimbau, K. Siellez, D. Meacher, B. Gendre, and M. Boer, Astrophys. J. 799, 69 (2015).
  * (87) D. Meacher, M. Coughlin, S. Morris, T. Regimbau, N. Christensen, S. Kandhasamy, V. Mandic, J. D. Romano, and E. Thrane, Phys. Rev. D 92, 063002 (2015).
  * (88) D. Meacher, K. Cannon, C. Hanna, T. Regimbau, and B. S. Sathyaprakash, Phys. Rev. D 93, 024018 (2016).
  * (89) T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis, B. Sathyaprakash, and S. Vitale, Phys. Rev. Lett. 118, 151105 (2017).
  * (90) J. Abadie et al. (LIGO Scientific Collabora-tion), Nucl. Instrum. Methods A624, 223 (2010).
  * (91) B. Sathyaprakash, B. Schutz, and C. Van Den Broeck, Class. Quant. Gravit. 27, 215006 (2010).
  * (92) M. Seikel, S. Yahya, R. Maartens, and C. Clarkson, Phys. Rev. D 86, 083001 (2012); S. Yahya, M. Seikel, C. Clarkson, R. Maartens, and M. Smith, Phys. Rev. D 89, 023503 (2014); V. C. Busti, C. Clarkson, and M. Seikel, Mon. Not. R. Astron. Soc. 441, L11 (2014).
  * (93) A. Shafieloo, A. G. Kim, and E. V. Linder, Phys. Rev. D 85, 123530 (2012); A. Shafieloo, A. G. Kim, and E. V. Linder, Phys. Rev. D 87, 023520 (2013).
  * (94) B. Chen and R. Kantowski, Phys. Rev. D 79 (2009) 104007; B. Chen and R. Kantowski, Phys. Rev. D 80, 044019 (2009).
  * (95) J. Chen, P. Wu, H. Yu, and Z. Li, J. Cosmol. Astropart. Phys. 10 (2012) 029.


# Probing NV and SiV charge state dynamics using high-voltage nanosecond pulse
and photoluminescence spectral analysis

Artur Pambukhchyan  Department of Chemistry, University of Southern
California, Los Angeles CA 90089, USA  Sizhe Weng  Department of Electrical
Engineering, University of Southern California, Los Angeles CA 90089, USA
Indu Aravind  Department of Physics & Astronomy, University of Southern
California, Los Angeles CA 90089, USA  Stephen B. Cronin  Department of
Electrical Engineering, University of Southern California, Los Angeles CA
90089, USA  Susumu Takahashi<EMAIL_ADDRESS>Department of
Chemistry, University of Southern California, Los Angeles CA 90089, USA
Department of Physics & Astronomy, University of Southern California, Los
Angeles CA 90089, USA

###### Abstract

Nitrogen-vacancy (NV) and silicon-vacancy (SiV) color defects in diamond are
promising systems for applications in quantum technology. The NV and SiV
centers have multiple charge states, and their charge states have different
electronic, optical and spin properties. For the NV centers, most
investigations for quantum sensing applications are targeted on the negatively
charged NV (NV-), and it is important for the NV centers to be in the NV-
state. However, it is known that the NV centers are converted to the neutrally
charged state (NV0) under laser excitation. An energetically favorable charge
state for the NV and SiV centers depends on their local environments. It is
essential to understand and control the charge state dynamics for their
quantum applications. In this work, we discuss the charge state dynamics of NV
and SiV centers under high-voltage nanosecond pulse discharges. The NV and SiV
centers coexist in the diamond crystal. The high-voltage pulses enable
manipulating the charge states efficiently. These voltage-induced changes in
charge states are probed by their photoluminescence spectral analysis. The
analysis result from the present experiment shows that the high-voltage
nanosecond pulses cause shifts of the chemical potential and can convert the
charge states of NV and SiV centers with the transition rates of $\sim$ MHz.
This result also indicates that the major population of the SiV centers in the
sample is the doubly negatively charged state (SiV2-), which is often
overlooked because of its non-fluorescent and non-magnetic nature. This
demonstration paves a path for a method of rapid manipulation of the NV and
SiV charge states in the future.

††: Materials for Quantum Technology

  * July 2023

##  1 Introduction

Diamond is a wide bandgap semiconductor and contains many impurities and
defects in the diamond lattice [1, 2]. The types and concentrations of these
defect centers are responsible for the color of diamond crystals. These color
centers were first studied by the gem industry, and have also become research
interests for the development of high power electronics devices [3]. More
recently, color defect centers in diamond have been investigated for
applications of quantum technology, including quantum sensing and quantum
communications.

The nitrogen-vacancy (NV) center in diamond is a fluorescent spin defect in
the diamond lattice. The NV center is promising for applications in quantum
information science because of its unique properties, including the capability
to polarize the spin state with optical excitation, stable photoluminescence
(PL) and its long spin coherence times [4, 5, 6, 7]. NV-based quantum sensing
techniques, nuclear magnetic resonance and electron spin resonance
spectroscopy have been demonstrated with nanoscale sample volumes [8, 9, 10,
11, 12, 13, 14]. The NV center is also a great candidate for single photon
sources. Electroluminescence (EL) using a single NV center has been
demonstrated [15, 16]. When the NV center is formed in the diamond lattice,
the NV center is in one of several charge states, including negatively charged
NV (NV-), neutrally charged state (NV0) and positively charged state NV (NV+))
[17]. While most NV research is based on the NV- center, it is important to
understand and control the dynamics of the NV charge state for the formation
and stabilization of the NV- state as well as applications [18, 19, 20, 21,
22]. The NV charge state can be converted under optical excitation [23, 19,
24], electrochemical control [25, 26] and gate-operation of the NV centers
[27, 28, 29].

The silicon-vacancy (SiV) center is another important optically active defect
in diamond. The PL of the SiV center has been found to be particularly
sensitive to the charge state of this defect. In particular, the negatively
charged SiV- emits PL at 738 nm [30, 31, 32], whereas the charge neutral
(SiV0) emits at 946 nm [33, 34]. Unlike NV emission, which is quite broad
($\sim$ 200 nm) at room temperature due to a wide array of phonon sidebands,
the SiV feature is spectrally narrow (4 nm at room temperature). These
spectrally narrow photons have a high degree of indistinguishability and are
better able to couple to photonic crystal cavities and waveguides making them
promising candidates for solid-state-based quantum technologies. Single SiV
defects have been detected in single crystal diamond [35]. The SiV defect
provides a unique molecular structure with the interstitial silicon atom being
placed between two unoccupied carbon sites [36]. This results in an inversion
symmetry ($D_{3d}$) of the split-vacancy configuration, which is subsequently
responsible for the high spectral stability as it protects the optical
transition from local electric field fluctuations [37]. Studies of the
electronic structure of the defect have revealed an $S=1$ ground state for the
neutral SiV [32] and a $S=1/2$ ground state for negative configuration [33,
38, 30]. In addition, a double negatively charged state (SiV2-) is also
expected to be stable in diamond, however, SiV2- has no optical transition and
accessible spin levels, and it may be difficult to be detected [39, 40]. It is
known that the charge state of SiV- centers can be converted under optical
illumination [38, 41]. However, the nature of the charge state and the
mechanisms of the conversion are not known well.

Charge state control of the NV and SiV centers by applied electric fields may
be useful for some applications in quantum devices. Doi et al. reported
electrical charge state manipulation of NV centers in diamond using a p-i-n
heterostructure device [42]. Here, the charge state was modulated via a
current driven effect (rather than a voltage driven effect), in which NV0
states were induced by current injection. In contrast, the more recent work of
Weng et al. uses high voltage pulses to overcome the large Schottky barrier at
the copper/diamond interface, injecting charge, and resulting in SiV- charge
state conversion [43]. In addition, this high voltage approach does not
require heterostructure fabrication and is particularly advantageous
considering the relatively large contact resistances associated with this wide
band gap semiconductor.

In this work, we probe the charge state dynamics of NV and SiV centers in
diamond with the application of high voltage nanosecond pulses and the
analysis of their PL spectra. We employ a type-IIa diamond, which contains
both NV and SiV centers. Recent work by Weng et al. shows the efficient
increase of the SiV population under high voltage nanosecond pulses [43].
However, the efficiency of the charge conversion and the mechanism of the
charge state dynamics have yet to be fully understood. Here, we discuss the
charge conversion of NV centers under high voltage nanosecond pulses and
characterize their conversion rates quantitatively. The high voltage pulses
are typically 10-20 ns long and are $\sim$ 5 kV in magnitude. The high voltage
pulses are repetitively applied to convert the NV charge states. The
repetition rate of the high-voltage pulse is varied from 1 Hz to 1 kHz to
modulate the charge conversion dynamics. The PL measurement of both NV and SiV
centers is performed by acquiring the PL signal with application of a
continuous-wave (CW) laser excitation and repetitive high-voltage pulses. We
analyze the PL spectra to extract the populations of the negatively charged
state of NV centers ([NV-]) and neutrally charged states of NV centers
([NV0]). We also analyze the change in the population of the SiV- states. The
result shows that $\sim$13 % of the NV- population is converted to NV0 with
the application of the high-voltage pulses with a repetition rate of 1 kHz.
The population of the SiV- is also increased by $\sim$20 times. This study
paves a way for the method development of fast manipulation of the NV and SiV
charge states.

##  2 Experiment

Figure 1:  Overview of the PL measurement. (a) Schematic diagram of the
experimental setup showing a typical lineshape of the high-voltage pulse. (b)
Diagram representing the experimental sequence, in which the high-voltage
pulse has the duration of $t_{hv}$, the repetition period is $T_{rep}$, and
the acquisition time ($T_{acq}$) is 10 s in the present case. (c) PL spectra
plotted as a function of wavelength. The spectra were collected with various
repetition rates. The full-widths at half-maximum of NV PL spectra are
indicated by orange dashed lines. The peak position of the SiV PL emission is
indicated by the purple dashed line.

In this experiment, polycrystalline type IIa CVD diamond (Element Six, Ltd.,
Thermal Management Grade, TM220) was studied. The size of the diamond crystal
is 10$\times$10$\times$0.5 mm3. No processing was applied on the diamond
sample. The experimental setup, as depicted in Fig. 1(a), involved attaching
strips of copper tape that were manually cut to create sharp tips with a
radius of 20 $\mu$m on both the top and bottom of the diamond sample. These
sharp electrode tips served to enhance the local electric field across the
diamond sample when the high voltage is applied. The separation between the
two electrodes is 0.5 mm in the lateral direction [43]. A continuous wave (CW)
laser with a wavelength of 532 nm was focused on the top surface of the
diamond using a microscope objective lens with a numerical aperture of 0.6 and
a working distance of 4 mm (Olympus ULCPlanFLN). The bottom electrode was
connected to a pulse generator (SSPG-20X, Transient Plasma Systems, Inc.)
capable of producing high-voltage pulses with a rise time of 5-10 ns
($dV/dt=10^{12}$ V/s) at a pulse repetition rate ranging from 1 to 1000 Hz.
The top electrode was grounded. Figure 1(a) displays the typical output
characteristics of the nanosecond pulses applied to the bottom electrode. PL
spectra of the diamond sample were captured from the diamond sample’s top
surface, near the tip of the top electrode, using an inVia micro-spectrometer
(Renishaw, Inc.).

As shown in Fig. 1(b), the PL measurement was performed under CW laser
excitation, PL acquisition and repetitive excitation of high voltage
nanosecond pulses. The power of the excitation laser was 37 $\mu$W. The PL
acquisition time ($T_{acq}$) was set to 10 s for the PL measurement. The PL
measurements were performed with and without the application of high voltage
nanosecond pulses. The repetition rate of high-voltage nanosecond pulses was
varied between 1 - 1000 Hz (the corresponding repetition period ($T_{rep}$) is
from 1 ms to 1 s). As shown in Fig. 1(a), the high-voltage pulses typically
have a duration of 10-20 ns and a peak voltage of up to $\sim$ 5 kV.

Figure 1(c) shows the PL spectra of NV centers with and without the
application of the high-voltage pulses. The repetition rate of the high-
voltage pulses was varied in 1-1000 Hz. As shown in Fig. 1(c), the NV spectra
become broader by applying the high-voltage pulses at higher repetition rates.
This broadening effect is much larger than that due to Stark shifts [44]. In
addition, as reported previously [43], the high-voltage pulses also increases
the PL intensity of the SiV- emission observed at 738 nm.

##  3 Analysis methods

Figure 2:  Population determination method using diamond PL spectra. (a)
Depiction of the PL spectral decomposition method to extract the population of
NV- and NV0 centers. The PL spectrum shown here was taken without the high-
voltage pulse excitation. (b) Background-subtracted PL spectrum used for the
determination of the SiV- population. The PL spectrum was taken with a
repetition rate of 1000 Hz.

###  3.1 Determination of [NV-] and [NV0] using PL decomposition analysis

Here, we discuss a method to estimate the populations of NV- and NV0 centers
using the PL decomposition analysis [45]. The population of the NV+ state is
not considered here. In the decomposition analysis, we consider an
experimentally observed PL spectrum (I($\lambda$)) to be the sum of the NV-
and NV0 PL spectra. Namely, I($\lambda$) is expressed by the following
equation,

$I(\lambda)=a_{1}I^{-}(\lambda)+a_{2}I^{0}(\lambda)$ (1)

where $I^{-}(\lambda)$ and $I^{0}(\lambda)$ are the normalized PL spectrum of
NV- and NV0 centers, respectively. $I^{-}(\lambda)$ and $I^{0}(\lambda)$ may
be the data reported previously. In the present case, we used $I^{-}(\lambda)$
and $I^{0}(\lambda)$ extracted from Ref. [46]. In addition, $a_{1}$ and
$a_{2}$ in Eq. 1, representing the contributions from NV- and NV0 centers on
the PL spectrum, are fitting parameters. $a_{1}$ and $a_{2}$ are proportional
to the populations of NV- and NV0 centers as well as their PL emission
efficiency. Therefore, by taking into account the emission efficiency ratio
between NV- and NV0 centers, the populations of NV- and NV0 can be expressed
by,

$\eqalign{[NV^{-}]=\frac{a_{1}}{a_{1}+\kappa_{532}a_{2}},\cr[NV^{0}]=1-[NV^{-}]=\frac{\kappa_{532}a_{2}}{a_{1}+\kappa_{532}a_{2}},}$
(2)

where $\kappa_{532}$ is the PL emission efficiency ratio between NV- and NV0
centers under 532 nm wavelength excitation, and was estimated to be 2.5
previously [45]. Fig. 2(a) shows the analysis result on the PL data without
high voltage pulse. As can be seen, the fit agrees well with the experimental
data. With the fit and Eqs. 1 and 2, we obtained [NV-] $=46\pm 1\%$ and [NV0]
$=54\pm 1\%$.

###  3.2 Estimate of SiV- population change

The population of charged SiV- centers under high-voltage pulses is estimated
by analyzing the intensity of the PL spectrum. As shown in Fig. 1(c), the SiV-
emission was observed at a wavelength of 738 nm. Since the PL peak area is
proportional to the SiV- population ([SiV-]), the population change due to the
application of high voltage pulses can be estimated by analyzing the peak
areas. Figure 2(b) shows the SiV PL spectrum with background correction. In
the correction, the background from the NV PL signal is removed by dividing
the background from the signal data. Then, the peak area ($I_{SiV-}$) is
obtained by fitting the data with a Lorentzian function, as shown in Fig.
2(b). In the present case, we analyzed the PL data taken with a repetition
rate of 1000 Hz and obtained $I_{SiV}=6.2\pm 0.1$ Using the value of
$I_{SiV}=6.2$ for 1000 Hz, the normalized peak area ($\tilde{I}_{SiV-}$) at
each repetition rate was calculated, namely, $\tilde{I}_{SiV-}=I_{SiV}/6.2$.

Figure 3:  Depiction of the simulation method for NV and SiV charge state
dynamics. (a) Two-state model to understand diamond PL spectra with charge
conversion dynamics. (b) Time evolution of the $p_{1}$ population. The
calculation was performed with the repetition rate of 1 Hz ($T_{rep}=1$ s),
$k^{0}_{12}=0.3$ $s^{-1}$, $k^{0}_{21}=(46/54)k^{0}_{12}$,
$k^{hv}_{12}=18.5\times 10^{6}$ $s^{-1}$ and $k^{hv}_{21}=9.2\times 10^{6}$
$s^{-1}$. $\bar{p_{1}}=42$% in the present case.

###  3.3 Simulation of two-charge state dynamics

In order to understand the repetition rate dependence of the PL spectra of NV
and SiV centers, we model their charge state dynamics as a two-state system.
Figure 3(b) shows the model, consisting of the population proportions of State
1 and 2 (denoted as $p_{1}$ and $p_{2}$). $k_{12}$ and $k_{21}$ are the
transition rates from State 1 to 2 and from State 2 to 1, respectively. The
time evolution of $p_{1}$ and $p_{2}$ are given by the following system rate
equations,

$\eqalign{\frac{dp_{1}(t)}{dt}=-k_{12}p_{1}(t)+k_{21}p_{2}(t),\cr\frac{dp_{2}(t)}{dt}=k_{12}p_{1}(t)-k_{21}p_{2}(t),}$
(3)

where

$k_{12(21)}=\cases{k^{0}_{12(21)}+k^{hv}_{12(21)}&during $t_{hv}$-period\\\
k^{0}_{12(21)}&during $T_{rep}$-period\\\ }$

$k^{0}_{12}$ and $k^{0}_{21}$ are the transition rates of $p_{1}$
$\rightarrow$ $p_{2}$ and $p_{2}$ $\rightarrow$ $p_{1}$ processes without the
application of high-voltage pulse. $k^{hv}_{12}$ and $k^{hv}_{21}$ are the
transition rates induced by the high-voltage pulse. Using the measurement
sequence shown in Fig. 1(b) and Eq. 3, the time evolution of $p_{1}(t)$ and
$p_{2}(t)=1-p_{1}(t)$ is calculated. In the calculation, we consider that the
high-voltage pulse is 10 ns long (i.e., $t_{hv}=10$ ns) and $k^{hv}_{12}$ and
$k^{hv}_{21}$ are constant during $t_{hv}$. Figure 3(b) shows the result of
$p_{1}(t)$. As shown in Fig. 3(b), $p_{1}(t)$ decreases rapidly by the
nanosecond voltage pulse, and then $p_{1}(t)$ relaxes back after the
nanosecond pulse. However, $p_{1}(t)$ does not reach the initial state before
the next nanosecond pulse because both $k^{0}_{12}$ and $k^{0}_{21}$ rates are
slower than the repetition rate of $1/T_{rep}=1$ Hz. Then, when the second
nanosecond pulse is applied, $p_{1}(t)$ decreases further. In a longer time
scale, $p_{1}(t)$ is oscillatory and the center of the oscillations decreases
asymptotically. From the time evolution result with $T_{acq}=10$ s, the mean
values of $p_{1}(t)$ and $p_{2}(t)$ ($\overline{p_{1}}$ and
$\overline{p_{2}}$) are calculated to compare with [NV-] and [NV0] as well as
the change of the normalized PL area of SiV- ($\tilde{I}_{SiV-}$). In the next
section, we discuss the population analyses for the NV and SiV centers. The NV
population is analyzed by considering the NV- population as $p_{1}(t)$ and the
NV0 population as $p_{2}(t)$. In the SiV analysis, SiV- population is
considered as $p_{2}(t)$.

##  4 Results and Discussion

Figure 4:  Analysis of the NV and SiV charge state modulation. (a) Populations of NV- and NV0 plotted as a function of the repetition rate of the high-voltage pulses. The y-axis shows the population proportions of NV- and NV0. The x-axis represents the repetition rate in Hz (1/s). The data without the high-voltage pulse (indicated as OFF) is shown by dashed lines. The error bars represent 95 % confidence intervals. The error bars are calculated using the obtained uncertainties of $a_{1}$ and $a_{2}$ and the error propagation. (b) The normalized peak areas of SiV- PL spectra plotted as a function of the repetition rates. (c) Schematic diagram depicting the charge dynamics with the application of the high-voltage pulses. (d) The formation energy of the NV and SiV charge states plotted as a function of the chemical potential ($E_{F}$). See Ref. [39, 47] for details about the charge formation energy levels.  Rep. rate (Hz) | [NV-] | [NV0] | $2\sigma_{NV-}$ | $2\sigma_{NV0}$ | $\tilde{I}_{SiV-}$ | $2\sigma_{SiV-}$
---|---|---|---|---|---|---
OFF | 46 | 54 | 1 | 1 | 5 | 1
1 | 41 | 59 | 1 | 1 | 21 | 1
5 | 38 | 62 | 1 | 1 | 36 | 1
10 | 35 | 65 | 1 | 1 | 58 | 1
50 | 34 | 66 | 1 | 1 | 83 | 1
500 | 33 | 67 | 1 | 1 | 100 | 1
1000 | 33 | 67 | 1 | 1 | 100 | 1
Table 1: Obtained [NV-], [NV0] and the normalized SiV PL intensity
($\tilde{I}_{SiV-}=I_{SiV-}/6.2$) values and their uncertainties. 2$\sigma$
represents 95 % confidence intervals.

###  4.1 NV Centers

Figure 4(a) and Table 1 show the NV- and NV0 populations proportions ([NV-]
and [NV0], respectively) as a function of the repetition rate of the applied
high-voltage pulses. The populations were obtained using the PL decomposition
analysis with Eqs. 1 and 2. As can be seen in Fig. 4, when the repetition rate
increases, [NV-] decreases while [NV0] increases. In the present experiment,
[NV-] ([NV0]) was changed from 46 % (54 %) to 33 % (67 %) when the repetition
rate is varied from 0 Hz (corresponding to no voltage pulse applied) to 1000
Hz. This population changes caused the PL broadening shown in Fig. 2. As shown
in Fig. 4, the changes of the populations is more pronounced in the repetition
range of 1-100 Hz, and the populations are nearly constant above the
repetition of rate of 100 Hz. To understand the charge state dynamics, the
result was fitted with the simulation described in the previous section. We
consider $p_{1}$ and $p_{2}$ in Fig. 3 to be NV- and NV0 states, respectively.
In the calculation of Eq. 3, we set the initial values of [NV-] and [NV0] to
be the values from the OFF result, namely, 46 % and 54 %, respectively. In
addition, the OFF high-voltage pulse result sets
$k^{0}_{21}/k^{0}_{12}=46/54$. Thus, the model has three parameters, namely,
$k^{0}_{12}$, $k^{hv}_{12}$ and $k^{hv}_{21}$. Figure 4(a) shows the fit
result of the experiment and model. As one can see, the model explains the
experimental data very well. From the fit, we obtained $k^{0}_{12}=0.3\pm 0.1$
Hz and $k^{hv}_{12}=(18.5\pm 0.1)\times 10^{6}$ Hz and $k^{hv}_{21}=(9.2\pm
0.3)\times 10^{6}$ Hz. Therefore, the result showed that the NV charge
conversion rates by the high-voltage nanosecond pulse is in the range of MHz.

###  4.2 SiV Centers

Figure 4(b) shows the normalized PL intensity of the SiV- centers plotted as a
function of the repetition rate. The PL intensity without the application of
high voltage pulses is also indicated by the purple dashed line. The PL
intensity is proportional to the SiV- population. As can be seen, the PL
intensity increased from OFF to 1000 Hz by $\sim$20 times, indicating a
significant increase of the SiV- population. Moreover, we perform the PL
measurement of SiV0 at 946 nm (not shown), and no PL signal was observed.
However, it is possible that phonon broadening effects on SiV0 PL hamper the
observation at room temperature since most of the previous PL measurements of
SiV0 have been done at low temperatures [33, 34, 48]. Therefore, in Fig. 4(b),
we showed the SiV- population normalized by the value measured at the
repetition rate of 1000 Hz ($\tilde{I}_{SiV-}$) instead of the SiV- portion in
all SiV charge state populations since we have no information of the
populations of other SiV charge states such as SiV0 and SiV2-. Similarly to
the NV centers (see Fig. 4(a)), the PL intensity changes are pronounced in the
range of 1-100 Hz, indicating the time scale of the SiV charge state dynamics
is similar to that of the NV centers. Finally, we determine the charge
conversion rates using the model shown in Sect. 3.3. In the analysis, we set
$k^{0}_{12}/k^{0}_{21}=5/95$ and $k^{0}_{21}=0$ because of the initial and
final states of $\tilde{I}_{SiV-}$. By fitting with Eq. 3, we obtained
$k^{0}_{12}=0.033\pm 0.007$ Hz and $k^{hv}_{12}=(8.58\pm 0.01)\times 10^{6}$
Hz.

###  4.3 Charge dynamics with high-voltage nanosecond pulse

As can be seen in many semiconductors, band bending occurs when an electric
field created by a high voltage pulse ($\epsilon^{hv}$) is applied to the
diamond crystal. Then, the band bending causes electrons moves to the high
voltage side and creates the unbalance of the charge density over the sample
volume, with which the anode size has higher electron density and the cathode
side has lower electron density (see Fig. 4(c)). In addition, the charge
unbalance induces an electric field oppose to the high-voltage pulse
($\epsilon^{ind}$ in Fig. 4(c)). Once the high-voltage pulse is off, electrons
diffuse over the sample volume to cancel the charge separation. With the use
of the high-voltage nanosecond pulse, electrons get a high kinetic energy and
the charge movement process happens in short time while the charge diffusion
is a slow process. The modulation of the electron density also changes the
chemical potential. Since the present experiment measures the PL signals from
the small detection volume on the cathode side, the NV and SiV centers in the
detection volume experience a lower chemical potential than that with no high-
voltage pulse. Figure 4 shows the formation energies of NV and SiV charge
states obtained from Ref. [39, 47]. As can be seen, NV- charge state is stable
when the chemical potential is $\sim$2.6 eV above the valence band energy (the
chemical potential ($E_{F}$) is defined from the valence band energy here).
When $E_{F}$ is below $\sim$2.6 eV, NV0 is more stable. For SiV centers, the
SiV2- charge state is energetically favorable when $E_{F}$ is above $\sim$2.2
eV. When $E_{F}$ is between 2.2 and 1.5 eV, the SiV- is stable. Then, SiV0 is
favorable, when $E_{F}$ is below 1.5 eV.

In the present experiment, we observed that the [NV-] and [NV0] values are
similar when no high-voltage pulse is applied. This indicates the chemical
potential with no high-voltage pulse ($E_{F}^{0}$) is $\sim$2.6 eV. Then, if
$E_{F}$ decreases with the application of the high-voltage pulses
($E_{F}^{hv}$), [NV-] will decrease and [NV0] will increase. This is
consistent with the observation as shown in Fig. 4(a). For the SiV center,
since $E_{F}^{hv}$ gets closer to the SiV- and SiV2- crossing energy of 2.2
eV, SiV- will also increase. Again, this is consistent with the observation on
the SiV- population change. Therefore, the SiV result strongly supports that
the SiV- centers are converted from the SiV2- centers.

##  5 Summary

In summary, we investigate the charge state dynamics of NV and SiV centers
under the application of high voltage nanosecond pulses and analyze their PL
spectra. The results show that the application of the high-voltage pulses
converts the NV populations from NV- to NV0. Moreover, a significant increase
of the SiV- populations was observed with the application of the high-voltage
pulses. We show that the observed charge state dynamics can be explained by a
simple two state model and transition rates of MHz were extracted. The changes
of the populations of both NV and SiV centers are consistent with the lowering
of the chemical potential of the diamond by the high voltage pulses. The
result also supports that the charge state of SiV centers in the diamond
sample is non-fluorescent and non-magnetic SiV2- centers. The experimental
method can be further developed by combining with time-domain measurements to
study the charge dynamics at faster time scales and multi-charge state
processes as well as to demonstrate the manipulation of NV and SiV charge
states with a fast electric field pulse.

This research was supported by the National Science Foundation (NSF) Award No.
ECCS-2204667 (S.W. and S.B.C.) and the Army Research Office (ARO) Award No.
W911NF2210284 (I.A. and S.B.C.). A.P. and S.T. also thanks the supports from
NSF (ECCS-2204667 and CHE-2004252 with partial co-funding from the Quantum
Information Science program in the Division of Physics), the USC Anton B. Burg
Foundation, and the Searle scholars program.

## Data availability statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * [1] J Walker.  Optical absorption and luminescence in diamond.  Rep. Prog. Phys., 42(10):1605, 1979.
  * [2] Michael N. R. Ashfold, Jonathan P. Goss, Ben L. Green, Paul W. May, Mark E. Newton, and Chloe V. Peaker.  Nitrogen in diamond.  Chem. Rev., 120(12):5745–5794, 2020.
  * [3] Satoshi Koizumi, Hitoshi Umezawa, Julien Pernot, and Mariko Suzuki, editors.  Power Electronics Device Applications of Diamond Semiconductors.  Woodhead Publishing Series in Electronic and Optical Materials. Woodhead Publishing, 2018.
  * [4] A. Gruber, A. Dräbenstedt, C. Tietz, L. Fleury, J. Wrachtrup, and C. von Borczyskowski.  Scanning confocal optical microscopy and magnetic resonance on single defect centers.  Science, 276(5321):2012, 1997.
  * [5] F. Jelezko, T. Gaebel, I. Popa, A. Gruber, and J. Wrachtrup.  Observation of coherent oscillations in a single electron spin.  Phys. Rev. Lett., 92:076401, 2004.
  * [6] R. J. Epstein, F. M. Mendoza, Y. K. Kato, and D. D. Awschalom.  Anisotropic interactions of a single spin and dark-spin spectroscopy in diamond.  Nat. Phys., 1(2):94–98, 2005.
  * [7] Susumu Takahashi, Ronald Hanson, Johan van Tol, Mark S. Sherwin, and David D. Awschalom.  Quenching spin decoherence in diamond through spin bath polarization.  Phys. Rev. Lett., 101:047601, 2008.
  * [8] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, and J. Wrachtrup.  Ultralong spin coherence time in isotopically engineered diamond.  Nat. Mater., 8:383, 2008.
  * [9] J R Maze, P L Stanwix, J S Hodges, S Hong, J M Taylor, P Cappellaro, L Jiang, M V Gurudev Dutt, E Togan, A S Zibrov, A Yacoby, R L Walsworth, and M D Lukin.  Nanoscale magnetic sensing with an individual electronic spin in diamond.  Nature, 455:644, 2008.
  * [10] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, R. Walsworth A. Yacoby, and M. D. Lukin.  High-sensitivity diamond magnetometer with nanoscale resolution.  Nat. Phys., 4:810, 2008.
  * [11] Chathuranga Abeywardana, Viktor Stepanov, Franklin H. Cho, and Susumu Takahashi.  Electron spin resonance spectroscopy of small ensemble paramagnetic spins using a single nitrogen-vacancy center in diamond.  J. Appl. Phys., 120(12):123907, 2016.
  * [12] Benjamin Fortman, Junior Pena, Karoly Holczer, and Susumu Takahashi.  Demonstration of NV-detected ESR spectroscopy at 115 GHz and 4.2 T.  Appl. Phys. Lett., 116(17):174004, 2020.
  * [13] Benjamin Fortman, Laura Mugica-Sanchez, Noah Tischler, Cooper Selco, Yuxiao Hang, Karoly Holczer, and Susumu Takahashi.  Electron–electron double resonance detected NMR spectroscopy using ensemble NV centers at 230 GHz and 8.3 T.  J. Appl. Phys., 130(8):083901, 2021.
  * [14] Shang Li, Huijie Zheng, Zaili Peng, Mizuki Kamiya, Tomoyuki Niki, Viktor Stepanov, Andrey Jarmola, Yasuhiro Shimizu, Susumu Takahashi, Arne Wickenbrock, and Dmitry Budker.  Determination of local defect density in diamond by double electron-electron resonance.  Phys. Rev. B, 104:094307, 2021.
  * [15] M. Haruyama, H. Kato, M. Ogura, Y. Kato, D. Takeuchi, S. Yamasaki, T. Iwasaki, H. Morishita, M. Fujiwara, N. Mizuochi, and T. Makino.  Electroluminescence of negatively charged single NV centers in diamond.  Appl. Phys. Lett., 122(7), 2023.
  * [16] Yuhang Guo, Wei Zhu, Jiaxin Zhao, Shengran Lin, Yuanjie Yang, Liren Lou, and Guanzhong Wang.  Electroluminescence of NV by impact excitation and Stark shift in a MIM diamond structure.  Appl. Phys. Lett., 119(25):252102, 2021.
  * [17] Peter Deák, Bálint Aradi, Moloud Kaviani, Thomas Frauenheim, and Adam Gali.  Formation of nv centers in diamond: A theoretical study based on calculated transitions and migration of nitrogen and vacancy related defects.  Phys. Rev. B, 89:075203, 2014.
  * [18] K.-M. C. Fu, C. Santori, P. E. Barclay, and R. G. Beausoleil.  Conversion of neutral nitrogen-vacancy centers to negatively charged nitrogen-vacancy centers through selective oxidation.  Appl. Phys. Lett., 96(12):121907, 2010.
  * [19] B. J. Shields, Q. P. Unterreithmeier, N. P. de Leon, H. Park, and M. D. Lukin.  Efficient readout of a single spin state in diamond via spin-to-charge conversion.  Phys. Rev. Lett., 114:136402, 2015.
  * [20] Harishankar Jayakumar, Siddharth Dhomkar, Jacob Henshaw, and Carlos A. Meriles.  Spin readout via spin-to-charge conversion in bulk diamond nitrogen-vacancy ensembles.  Appl. Phys. Lett., 113(12):122404, 2018.
  * [21] Amila Ariyaratne, Dolev Bluvstein, Bryan A. Myers, and Ania C. Bleszynski Jayich.  Nanoscale electrical conductivity imaging using a nitrogen-vacancy center in diamond.  Nat. Commun., 9:2406, 2018.
  * [22] S. Baier, C. E. Bradley, T. Middelburg, V. V. Dobrovitski, T. H. Taminiau, and R. Hanson.  Orbital and spin dynamics of single neutrally-charged nitrogen-vacancy centers in diamond.  Phys. Rev. Lett., 125:193601, 2020.
  * [23] N Aslam, G Waldherr, P Neumann, F Jelezko, and J Wrachtrup.  Photo-induced ionization dynamics of the nitrogen vacancy defect in diamond investigated by single-shot charge state detection.  N. J. Phys., 15(1):013064, 2013.
  * [24] Peng Ji and M. V. Gurudev Dutt.  Charge state dynamics of the nitrogen vacancy center in diamond under 1064-nm laser excitation.  Phys. Rev. B, 94:024101, Jul 2016.
  * [25] Sinan Karaveli, Ophir Gaathon, Abraham Wolcott, Reyu Sakakibara, Or A. Shemesh, Darcy S. Peterka, Edward S. Boyden, Jonathan S. Owen, Rafael Yuste, and Dirk Englund.  Modulation of nitrogen vacancy charge state and fluorescence in nanodiamonds using electrochemical potential.  Proc. Natl. Acad. Sci. U.S.A., 113(15):3938–3943, 2016.
  * [26] M. V. Hauf, B. Grotz, B. Naydenov, M. Dankerl, S. Pezzagna, J. Meijer, F. Jelezko, J. Wrachtrup, M. Stutzmann, F. Reinhard, and J. A. Garrido.  Chemical control of the charge state of nitrogen-vacancy centers in diamond.  Phys. Rev. B, 83:081304, 2011.
  * [27] Bernhard Grotz, Moritz V. Hauf, Markus Dankerl, Boris Naydenov, Sébastien Pezzagna, Jan Meijer, Fedor Jelezko, Jörg Wrachtrup, Martin Stutzmann, Friedemann Reinhard, and Jose A. Garrido.  Charge state manipulation of qubits in diamond.  Nat. Commun., 3:729, 2012.
  * [28] Christoph Schreyvogel, Marco Wolfer, Hiromitsu Kato, Matthias Schreck, and Christoph E. Nebel.  Tuned nv emission by in-plane al-schottky junctions on hydrogen terminated diamond.  Sci. Rep., 4:3634, 2014.
  * [29] Moritz V. Hauf, Patrick Simon, Nabeel Aslam, Matthias Pfender, Philipp Neumann, Sébastien Pezzagna, Jan Meijer, Jörg Wrachtrup, Martin Stutzmann, Friedemann Reinhard, and José A. Garrido.  Addressing single nitrogen-vacancy centers in diamond with transparent in-plane gate structures.  Nano Lett., 14(5):2359–2364, 2014.
  * [30] B. L. Green, M. W. Doherty, E. Nako, N. B. Manson, U. F. S. D’Haenens-Johansson, S. D. Williams, D. J. Twitchen, and M. E. Newton.  Electronic structure of the neutral silicon-vacancy center in diamond.  Phys. Rev. B, 99:161112, 2019.
  * [31] Tina Müller, Christian Hepp, Benjamin Pingault, Elke Neu, Stefan Gsell, Matthias Schreck, Hadwig Sternschulte, Christoph Becher Doris Steinmüller-Nethl, and Mete Atatüre.  Optical signatures of silicon-vacancy spins in diamond.  Nat. Commun., 5:3328, 2014.
  * [32] Lachlan J. Rogers, Kay D. Jahnke, Marcus W. Doherty, Andreas Dietrich, Liam P. McGuinness, Christoph Müller, Tokuyuki Teraji, Hitoshi Sumiya, Junichi Isoya, Neil B. Manson, and Fedor Jelezko.  Electronic structure of the negatively charged silicon-vacancy center in diamond.  Phys. Rev. B, 89:235101, 2014.
  * [33] U. F. S. D’Haenens-Johansson, A. M. Edmonds, B. L. Green, M. E. Newton, G. Davies, P. M. Martineau, R. U. A. Khan, and D. J. Twitchen.  Optical properties of the neutral silicon split-vacancy center in diamond.  Phys. Rev. B, 84:245208, 2011.
  * [34] B. L. Green, S. Mottishaw, B. G. Breeze, A. M. Edmonds, U. F. S. D’Haenens-Johansson, M. W. Doherty, S. D. Williams, D. J. Twitchen, and M. E. Newton.  Neutral silicon-vacancy center in diamond: Spin polarization and lifetimes.  Phys. Rev. Lett., 119:096402, 2017.
  * [35] Elke Neu, David Steinmetz, Janine Riedrich-Möller, Stefan Gsell, Martin Fischer, Matthias Schreck, and Christoph Becher.  Single photon emission from silicon-vacancy colour centres in chemical vapour deposition nano-diamonds on iridium.  N. J. Phys., 13(2):025012, 2011.
  * [36] J. P. Goss, R. Jones, S. J. Breuer, P. R. Briddon, and S. Öberg.  The twelve-line 1.682 ev luminescence center in diamond and the vacancy-silicon complex.  Phys. Rev. Lett., 77:3041–3044, 1996.
  * [37] Christian Hepp, Tina Müller, Victor Waselowski, Jonas N. Becker, Benjamin Pingault, Hadwig Sternschulte, Doris Steinmüller-Nethl, Adam Gali, Jeronimo R. Maze, Mete Atatüre, and Christoph Becher.  Electronic structure of the silicon vacancy color center in diamond.  Phys. Rev. Lett., 112:036405, 2014.
  * [38] Siddharth Dhomkar, Pablo R. Zangara, Jacob Henshaw, and Carlos A. Meriles.  On-demand generation of neutral and negatively charged silicon-vacancy centers in diamond.  Phys. Rev. Lett., 120:117401, 2018.
  * [39] Adam Gali and Jeronimo R. Maze.  Ab initio study of the split silicon-vacancy defect in diamond: Electronic structure and related properties.  Phys. Rev. B, 88:235205, 2013.
  * [40] B. G. Breeze, C. J. Meara, X. X. Wu, C. P. Michaels, R. Gupta, P. L. Diggle, M. W. Dale, B. L. Cann, T. Ardon, U. F. S. D’Haenens-Johansson, I. Friel, M. J. Rayson, P. R. Briddon, J. P. Goss, M. E. Newton, and B. L. Green.  Doubly charged silicon vacancy center, si-n complexes, and photochromism in n and si codoped diamond.  Phys. Rev. B, 101:184115, 2020.
  * [41] Louis Nicolas, Tom Delord, Paul Huillery, Clément Pellet-Mary, and Gabriel Hétet.  Sub-ghz linewidth ensembles of siv centers in a diamond nanopyramid revealed by charge state conversion.  ACS Photonics, 6(10):2413–2420, 2019.
  * [42] Y. Doi, T. Makino, H. Kato, D. Takeuchi, M. Ogura, H. Okushi, H. Morishita, T. Tashima, S. Miwa, S. Yamasaki, P. Neumann, J. Wrachtrup, Y. Suzuki, and N. Mizuochi.  Deterministic electrical charge-state initialization of single nitrogen-vacancy center in diamond.  Phys. Rev. X, 4:011057, 2014.
  * [43] Sizhe Weng, Christopher Coleman, Indu Aravind, Yu Wang, Bofan Zhao, and Stephen B. Cronin.  Voltage-induced modulation in the charge state of Si-vacancy defects in diamond using high voltage nanosecond pulses.  Appl. Phys. Lett., 119(17):171101, 2021.
  * [44] Ph. Tamarat, T. Gaebel, J. R. Rabeau, M. Khan, A. D. Greentree, H. Wilson, L. C. L. Hollenberg, S. Prawer, P. Hemmer, F. Jelezko, and J. Wrachtrup.  Stark shift control of single optical centers in diamond.  Phys. Rev. Lett., 97:083002, 2006.
  * [45] Scott T. Alsid, John F. Barry, Linh M. Pham, Jennifer M. Schloss, Michael F. O’Keeffe, Paola Cappellaro, and Danielle A. Braje.  Photoluminescence decomposition analysis: A technique to characterize $\mathrm{N}$-$v$ creation in diamond.  Phys. Rev. Appl., 12:044003, 2019.
  * [46] Jan Jeske, Desmond W. M. Lau, Xavier Vidal, Liam P. McGuinness, Philipp Reineck, Brett C. Johnson, Marcus W. Doherty, Jeffrey C. McCallum, Shinobu Onoda, Fedor Jelezko, Takeshi Ohshima, Thomas Volz, Jared H. Cole, Brant C. Gibson, and Andrew D. Greentree.  Stimulated emission from nitrogen-vacancy centres in diamond.  Nat. Commun., 8:14000, 2017.
  * [47] Peter Deák, Bálint Aradi, Moloud Kaviani, Thomas Frauenheim, and Adam Gali.  Calculation ofthe transitions and migration of nitrogen and vacancy related defects,with implications on the formation of nv centers in bulk diamond.  arXiv, page 1311.6598, 2013.
  * [48] Brendon C. Rose, Ding Huang, Zi-Huai Zhang, Paul Stevenson, Alexei M. Tyryshkin, Sorawis Sangtawesin, Srikanth Srinivasan, Lorne Loudin, Matthew L. Markham, Andrew M. Edmonds, Daniel J. Twitchen, Stephen A. Lyon, and Nathalie P. de Leon.  Observation of an environmentally insensitive solid-state spin defect in diamond.  Science, 361(6397):60–63, 2018.